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

4
#include "CovarianceFunction.hpp"
niklas.baumgarten's avatar
niklas.baumgarten committed
5
#include "LagrangeDiscretization.hpp"
6
#include "basics/Utilities.hpp"
niklas.baumgarten's avatar
niklas.baumgarten committed
7
#include "Algebra.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
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;

22
class CirculantEmbedding {
23
24
protected:
    int internalCounter = 0;
25

26
    ComplexNormalDistribution normalDist;
27

28
    double evtol = 1e-10;
29

30
    string fieldType = "LogNormal";
31

32
33
    double mean = 0.0;

niklas.baumgarten's avatar
niklas.baumgarten committed
34
    LagrangeDiscretization disc;
niklas.baumgarten's avatar
niklas.baumgarten committed
35

36
37
    Meshes &meshes;

38
public:
39
    explicit CirculantEmbedding(Meshes &meshes) :
niklas.baumgarten's avatar
niklas.baumgarten committed
40
        meshes(meshes), normalDist(meshes), disc(meshes, 0) {
41

42
        config.get("evtol", evtol);
niklas.baumgarten's avatar
niklas.baumgarten committed
43
        config.get("StochasticField", fieldType);
44
45
46
        config.get("mean", mean);
    }

47
    virtual void generateFineSample(const SampleID &id,
48
49
50
                                    Vector *&fineSample,
                                    Vector *&coarseSample) = 0;

51
    virtual void generateCoarseSample(const SampleID &id,
52
53
54
                                      Vector *&fineSample,
                                      Vector *&coarseSample) = 0;

55
    virtual ~CirculantEmbedding() {};
56
57
};

58
59
class CirculantEmbedding1D : public CirculantEmbedding, public CovarianceFunction1D {
public:
60
    NormalDistributionComplexSequence1D generator;
niklas.baumgarten's avatar
niklas.baumgarten committed
61

niklas.baumgarten's avatar
niklas.baumgarten committed
62
    CVector fineComplexField;
63

niklas.baumgarten's avatar
niklas.baumgarten committed
64
    RVector sqrtEigenvalues;
65

66
67
68
69
70
71
    int numberOfCellsEmbedded = 0;
    int numberOfCells = 0;

    double domain_start = 0;
    double domain_end = 1;

72
    explicit CirculantEmbedding1D(int l, Meshes &meshes) :
niklas.baumgarten's avatar
niklas.baumgarten committed
73
74
        CirculantEmbedding(meshes),
        generator(meshes) {
75
76
        numberOfCells = meshes[l].CellCount();
        sqrtEigenvalues = computeSqrtEV();
77
        numberOfCellsEmbedded = sqrtEigenvalues.size();
niklas.baumgarten's avatar
niklas.baumgarten committed
78
        generator.Resize(numberOfCellsEmbedded);
79
80
    }

81
    void generateFineSample(const SampleID &id,
82
83
                            Vector *&fineSample,
                            Vector *&coarseSample) override;
84

85
    void generateCoarseSample(const SampleID &id,
86
87
                              Vector *&fineSample,
                              Vector *&coarseSample) override;
88

niklas.baumgarten's avatar
niklas.baumgarten committed
89
    CVector generateField(const SampleID &id);
90

niklas.baumgarten's avatar
niklas.baumgarten committed
91
92
    void createCovarianceMatrix(RVector &ar,
                                RVector &ac);
93

niklas.baumgarten's avatar
niklas.baumgarten committed
94
95
96
    static void embedCovarianceMatrix(const RVector &ar,
                                      const RVector &ac,
                                      RVector &br);
97

niklas.baumgarten's avatar
niklas.baumgarten committed
98
99
    static void computeEV(const RVector &br,
                          RVector &ev);
100

niklas.baumgarten's avatar
niklas.baumgarten committed
101
102
    void padding(RVector &ar,
                 RVector &ac);
103

niklas.baumgarten's avatar
niklas.baumgarten committed
104
    RVector computeSqrtEV();
105
106
107
108
};

class CirculantEmbedding2D : public CirculantEmbedding, public CovarianceFunction2D {
public:
109
    ComplexSequence2D fineComplexField;
110
111
112
113
114
115
116
117
118
119
120

    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;

121
    explicit CirculantEmbedding2D(int l, Meshes &meshes) :
122
        CirculantEmbedding(meshes) {
123
124
125
126
127
128
129
        numberOfCells[0] = (int) pow(2, l);
        numberOfCells[1] = (int) pow(2, l);
        sqrtEigenvalues = getSqrtEV();
        numberOfCellsEmbedded[0] = sqrtEigenvalues.size();
        numberOfCellsEmbedded[1] = sqrtEigenvalues[0].size();
    }

130
    void generateFineSample(const SampleID &id,
131
132
                            Vector *&fineSample,
                            Vector *&coarseSample) override;
133

134
    void generateCoarseSample(const SampleID &id,
135
136
                              Vector *&fineSample,
                              Vector *&coarseSample) override;
137

138
    ComplexSequence2D generateField(const SampleID &id);
139

140
141
    void createCovarianceMatrix(BlockToeplitzRow &ar,
                                BlockToeplitzColumn &ac);
142

143
144
145
    static void embedCovarianceMatrix(const BlockToeplitzRow &ar,
                                      const BlockToeplitzColumn &ac,
                                      BlockCirculantRow &br);
146

147
148
    static void computeEV(const BlockCirculantRow &br,
                          EigenValues2D &ev);
149

150
151
    void padding(BlockToeplitzRow &ar,
                 BlockToeplitzColumn &ac);
152

153
    SqrtEigenValues2D getSqrtEV();
154
155
};

156
class MultilevelCirculantEmbedding : public SampleGenerator<Tensor> {
157
private:
niklas.baumgarten's avatar
niklas.baumgarten committed
158
159
    PlotMap plotMap;

160
161
162
163
164
165
166
    std::map<int, CirculantEmbedding *> circulantEmbeddings;

    Vector *fineSample = nullptr;

    Vector *coarseSample = nullptr;

protected:
167
    void drawSample(const SampleID &id) override {
168
        if (!id.coarse) {
169
170
            circulantEmbeddings[id.level.fine]->generateFineSample(id, fineSample,
                                                                   coarseSample);
niklas.baumgarten's avatar
niklas.baumgarten committed
171
172
173
174
            SampleSolution _fineSample(*fineSample, id); // Todo
            _fineSample.U = 1.0;
            _fineSample.name = "Kappa";
            plotMap.VtkPlot(_fineSample);
175
        } else {
176
177
            circulantEmbeddings[id.level.fine]->generateCoarseSample(id, fineSample,
                                                                     coarseSample);
niklas.baumgarten's avatar
niklas.baumgarten committed
178
179
180
            SampleSolution _coarseSample(*coarseSample, id); // Todo
            _coarseSample.name = "Kappa";
            plotMap.VtkPlot(_coarseSample);
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
        }
    }

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();
    }

200
    Tensor EvalSample(const cell &c) override {
201
202
203
204
205
206
207
208
        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"; }
};

209
#endif //M_CIRCULANTEMBEDDING_H