CirculantEmbedding.hpp 6.73 KB
Newer Older
1
2
3
#ifndef M_CIRCULANTEMBEDDING_H
#define M_CIRCULANTEMBEDDING_H

4
#include "CovarianceFunction.hpp"
5
#include "basics/Utilities.hpp"
niklas.baumgarten's avatar
niklas.baumgarten committed
6
#include "Algebra.hpp"
7
#include "SampleGenerator.hpp"
8
#include "NormalDistribution.hpp"
niklas.baumgarten's avatar
niklas.baumgarten committed
9
#include "basics/PlotMap.hpp"
10

11

12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
typedef std::vector<double> EigenValues1D;

typedef std::vector<double> SqrtEigenValues1D;

typedef std::vector<double> ToeplitzRow;

typedef std::vector<double> ToeplitzColumn;

typedef std::vector<double> CirculantRow;

typedef std::vector<std::vector<double>> EigenValues2D;

typedef std::vector<std::vector<double>> SqrtEigenValues2D;

typedef std::vector<std::vector<double>> BlockToeplitzRow;

typedef std::vector<std::vector<double>> BlockToeplitzColumn;

typedef std::vector<std::vector<double>> BlockCirculantRow;

32
class CirculantEmbedding {
33
34
protected:
    int internalCounter = 0;
35

36
    ComplexNormalDistribution normalDist;
37

38
    double evtol = 1e-10;
39

40
    string fieldType = "LogNormal";
41

42
43
    double mean = 0.0;

niklas.baumgarten's avatar
niklas.baumgarten committed
44
45
    MatrixGraphs cellMGraphs;

46
47
    Meshes &meshes;

48
public:
49
    explicit CirculantEmbedding(Meshes &meshes) :
50
51
        meshes(meshes),
        normalDist(meshes),
niklas.baumgarten's avatar
niklas.baumgarten committed
52
        cellMGraphs(meshes, dof(new LagrangeDoF(0))) {
53

54
        config.get("evtol", evtol);
niklas.baumgarten's avatar
niklas.baumgarten committed
55
        config.get("StochasticField", fieldType);
56
57
58
        config.get("mean", mean);
    }

59
    virtual void generateFineSample(const SampleID &id,
60
61
62
                                    Vector *&fineSample,
                                    Vector *&coarseSample) = 0;

63
    virtual void generateCoarseSample(const SampleID &id,
64
65
66
                                      Vector *&fineSample,
                                      Vector *&coarseSample) = 0;

67
    virtual ~CirculantEmbedding() {};
68
69
};

70
71
class CirculantEmbedding1D : public CirculantEmbedding, public CovarianceFunction1D {
public:
72
    NormalDistributionComplexSequence1D generator;
niklas.baumgarten's avatar
niklas.baumgarten committed
73

74
    ComplexSequence1D fineComplexField;
75
76
77

    SqrtEigenValues1D sqrtEigenvalues;

niklas.baumgarten's avatar
niklas.baumgarten committed
78

79
80
81
82
83
84
    int numberOfCellsEmbedded = 0;
    int numberOfCells = 0;

    double domain_start = 0;
    double domain_end = 1;

85
    explicit CirculantEmbedding1D(int l, Meshes &meshes) :
niklas.baumgarten's avatar
niklas.baumgarten committed
86
87
        CirculantEmbedding(meshes),
        generator(meshes) {
88
89
        numberOfCells = meshes[l].CellCount();
        sqrtEigenvalues = computeSqrtEV();
90
        numberOfCellsEmbedded = sqrtEigenvalues.size();
niklas.baumgarten's avatar
niklas.baumgarten committed
91
        generator.Resize(numberOfCellsEmbedded);
92
93
    }

94
    void generateFineSample(const SampleID &id,
95
96
                            Vector *&fineSample,
                            Vector *&coarseSample) override;
97

98
    void generateCoarseSample(const SampleID &id,
99
100
                              Vector *&fineSample,
                              Vector *&coarseSample) override;
101

102
//private:
103
    ComplexSequence1D generateField(const SampleID &id);
104

105
106
    void createCovarianceMatrix(ToeplitzRow &ar,
                                ToeplitzColumn &ac);
107

108
109
110
    static void embedCovarianceMatrix(const ToeplitzRow &ar,
                                      const ToeplitzColumn &ac,
                                      CirculantRow &br);
111

112
113
    static void computeEV(const CirculantRow &br,
                          EigenValues1D &ev);
114

115
116
    void padding(ToeplitzRow &ar,
                 ToeplitzColumn &ac);
117

118
    SqrtEigenValues1D computeSqrtEV();
119
120
121
122
};

class CirculantEmbedding2D : public CirculantEmbedding, public CovarianceFunction2D {
public:
123
    ComplexSequence2D fineComplexField;
124
125
126
127
128
129
130
131
132
133
134

    SqrtEigenValues2D sqrtEigenvalues;

    int numberOfCellsEmbedded[2] = {0, 0};
    int numberOfCells[2] = {0, 0};

    double domain1_start = 0;
    double domain1_end = 1;
    double domain2_start = 0;
    double domain2_end = 1;

135
    explicit CirculantEmbedding2D(int l, Meshes &meshes) :
136
        CirculantEmbedding(meshes) {
137
138
139
140
141
142
143
        numberOfCells[0] = (int) pow(2, l);
        numberOfCells[1] = (int) pow(2, l);
        sqrtEigenvalues = getSqrtEV();
        numberOfCellsEmbedded[0] = sqrtEigenvalues.size();
        numberOfCellsEmbedded[1] = sqrtEigenvalues[0].size();
    }

144
    void generateFineSample(const SampleID &id,
145
146
                            Vector *&fineSample,
                            Vector *&coarseSample) override;
147

148
    void generateCoarseSample(const SampleID &id,
149
150
                              Vector *&fineSample,
                              Vector *&coarseSample) override;
151

152
//private:
153
    ComplexSequence2D generateField(const SampleID &id);
154

155
156
    void createCovarianceMatrix(BlockToeplitzRow &ar,
                                BlockToeplitzColumn &ac);
157

158
159
160
    static void embedCovarianceMatrix(const BlockToeplitzRow &ar,
                                      const BlockToeplitzColumn &ac,
                                      BlockCirculantRow &br);
161

162
163
    static void computeEV(const BlockCirculantRow &br,
                          EigenValues2D &ev);
164

165
166
    void padding(BlockToeplitzRow &ar,
                 BlockToeplitzColumn &ac);
167

168
    SqrtEigenValues2D getSqrtEV();
169
170
};

171
class MultilevelCirculantEmbedding : public SampleGenerator<Tensor> {
172
private:
niklas.baumgarten's avatar
niklas.baumgarten committed
173
174
    PlotMap plotMap;

175
176
177
178
179
180
181
    std::map<int, CirculantEmbedding *> circulantEmbeddings;

    Vector *fineSample = nullptr;

    Vector *coarseSample = nullptr;

protected:
182
    void drawSample(const SampleID &id) override {
183
        if (!id.coarse) {
184
185
            circulantEmbeddings[id.level.fine]->generateFineSample(id, fineSample,
                                                                   coarseSample);
niklas.baumgarten's avatar
niklas.baumgarten committed
186
187
188
189
            SampleSolution _fineSample(*fineSample, id); // Todo
            _fineSample.U = 1.0;
            _fineSample.name = "Kappa";
            plotMap.VtkPlot(_fineSample);
190
        } else {
191
192
            circulantEmbeddings[id.level.fine]->generateCoarseSample(id, fineSample,
                                                                     coarseSample);
niklas.baumgarten's avatar
niklas.baumgarten committed
193
194
195
            SampleSolution _coarseSample(*coarseSample, id); // Todo
            _coarseSample.name = "Kappa";
            plotMap.VtkPlot(_coarseSample);
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
        }
    }

public:
    MultilevelCirculantEmbedding(Meshes &meshes) : SampleGenerator(meshes) {
        for (int l = meshes.pLevel(); l <= meshes.Level(); l++) {
            if (meshes.dim() == 1)
                circulantEmbeddings[l] = new CirculantEmbedding1D(l, meshes);
            if (meshes.dim() == 2)
                circulantEmbeddings[l] = new CirculantEmbedding2D(l, meshes);
        }
    }

    ~MultilevelCirculantEmbedding() {
        for (auto &circEmb : circulantEmbeddings)
            delete circEmb.second;
        circulantEmbeddings.clear();
    }

215
    Tensor EvalSample(const cell &c) override {
216
217
218
219
220
221
222
223
        if (fineSample && !coarseSample) return fineSample->operator()(c(), 0) * One;
        if (coarseSample && !fineSample) return coarseSample->operator()(c(), 0) * One;
        Exit("Error in Generator")
    }

    string Name() const override { return "Circulant Embedding"; }
};

224
#endif //M_CIRCULANTEMBEDDING_H