CirculantEmbedding.cpp 14.3 KB
Newer Older
1
2
#include "CirculantEmbedding.hpp"
#include "main/MultilevelPlotter.hpp"
3
4
5
6
7
#include "dof/BasicDoFs.hpp"

void CirculantEmbedding1D::generateFineSample(SampleID id,
                                              Vector *&fineSample,
                                              Vector *&coarseSample) {
niklas.baumgarten's avatar
niklas.baumgarten committed
8
    fineSample = new Vector(cellMGraphs[id.level - this->meshes.pLevel()]);
9

niklas.baumgarten's avatar
naming    
niklas.baumgarten committed
10
    if (internalCounter == 0)
11
12
        fineComplexField = generateField();

13
14
    for (cell c = fineSample->cells(); c != fineSample->cells_end(); c++) {
        int i = floor(c()[0] * sqrt(fineSample->size()));
niklas.baumgarten's avatar
naming    
niklas.baumgarten committed
15
        if (fieldType == "Gaussian") {
16
            if (internalCounter == 0)
17
                fineSample->operator()(c(), 0) = fineComplexField[i].real() + mean;
18
            else
19
                fineSample->operator()(c(), 0) = fineComplexField[i].imag() + mean;
niklas.baumgarten's avatar
naming    
niklas.baumgarten committed
20
        } else if (fieldType == "LogNormal") {
21
            if (internalCounter == 0)
22
                fineSample->operator()(c(), 0) = exp(fineComplexField[i].real() + mean);
23
            else
24
                fineSample->operator()(c(), 0) = exp(fineComplexField[i].imag() + mean);
niklas.baumgarten's avatar
niklas.baumgarten committed
25
26
        } else Exit("Has to be Gaussian or LogNormal")
    }
niklas.baumgarten's avatar
naming    
niklas.baumgarten committed
27
28
    internalCounter += 1;
    internalCounter %= 2;
niklas.baumgarten's avatar
niklas.baumgarten committed
29

niklas.baumgarten's avatar
niklas.baumgarten committed
30
31
    if (coarseSample) delete coarseSample;

niklas.baumgarten's avatar
niklas.baumgarten committed
32
    if (plotting)
33
        plotter->PlotVector("kappa", *fineSample,
niklas.baumgarten's avatar
niklas.baumgarten committed
34
                            1, l, "CellData");
35
36
}

37
38
39
void CirculantEmbedding1D::generateCoarseSample(SampleID id,
                                                Vector *&fineSample,
                                                Vector *&coarseSample) {
niklas.baumgarten's avatar
niklas.baumgarten committed
40
    coarseSample = new Vector(cellMGraphs[l - this->meshes.pLevel() - 1]);
41
42
43
44

    (*coarseSample) = 0;
    for (cell c = coarseSample->cells(); c != coarseSample->cells_end(); c++) {
        row coarseRow = coarseSample->find_row(c());
45
        for (int k = 0; k < c.Children(); ++k) {
46
            row fineRow = fineSample->find_row(c.Child(k));
47
            for (int i = 0; i < coarseRow.n(); ++i)
48
49
                coarseSample->operator()(coarseRow, i) +=
                    fineSample->operator()(fineRow, i);
50
51
        }
        for (int i = 0; i < coarseRow.n(); ++i)
52
            coarseSample->operator()(coarseRow, i) /= c.Children();
53
    }
niklas.baumgarten's avatar
niklas.baumgarten committed
54

niklas.baumgarten's avatar
niklas.baumgarten committed
55
56
    if (fineSample) delete fineSample;

niklas.baumgarten's avatar
niklas.baumgarten committed
57
    if (plotting)
58
        plotter->PlotVector("kappa", *coarseSample, 1,
niklas.baumgarten's avatar
niklas.baumgarten committed
59
                            l - 1, "CellData");
niklas.baumgarten's avatar
niklas.baumgarten committed
60
61
}

62
63
64
ComplexField1D CirculantEmbedding1D::generateField() {
    ComplexField1D Xf;
    ComplexField1D normallyDistributed;
65
66
    normallyDistributed.resize(numberOfCellsEmbedded);
    fineComplexField.clear();
67
    Xf.resize(numberOfCellsEmbedded);
niklas.baumgarten's avatar
niklas.baumgarten committed
68
    double z1, z2;
69
    for (int i = 0; i < numberOfCellsEmbedded; i++) {
niklas.baumgarten's avatar
niklas.baumgarten committed
70
        rnmg->get_normal_rn(z1, z2);
71
        normallyDistributed[i] = {z1, z2};
niklas.baumgarten's avatar
niklas.baumgarten committed
72
73
    }

74
    ComplexField1D product;
75
76
77
    product.resize(numberOfCellsEmbedded);
    for (int i = 0; i < numberOfCellsEmbedded; i++) {
        product[i] = sqrtEigenvalues[i] * normallyDistributed[i];
niklas.baumgarten's avatar
niklas.baumgarten committed
78
    }
79
    fft(product, Xf);
80
81
    double divisor = sqrt(numberOfCellsEmbedded);
    for (int i = 0; i < numberOfCellsEmbedded; i++)
82
83
        Xf[i] /= divisor;
    return Xf;
niklas.baumgarten's avatar
niklas.baumgarten committed
84
85
}

86
87
void CirculantEmbedding1D::createCovarianceMatrix(
    vector<double> &ar, vector<double> &ac) {
88
89
90
    const vector<double> &coordinate1 = linspace(domain_start, domain_end, numberOfCells);
    ar.resize(numberOfCells);
    ac.resize(numberOfCells);
niklas.baumgarten's avatar
niklas.baumgarten committed
91
92
    double y_rows = coordinate1.front();

93
    for (int index_c = 0; index_c < numberOfCells; index_c++) {
niklas.baumgarten's avatar
niklas.baumgarten committed
94
95
96
97
98
99
100
101
102
103
104
105
106
        double c1 = coordinate1[index_c];
        double x_rows = c1;
        double tau_rows = x_rows - y_rows;

        double x_columns = coordinate1.front();
        double y_columns = c1;
        double tau_columns = x_columns - y_columns;

        ar[index_c] = covariance_fct(&tau_rows);
        ac[index_c] = covariance_fct(&tau_columns);
    }
}

107
108
void CirculantEmbedding1D::embedCovarianceMatrix(vector<double> &ar, vector<double> &ac,
                                                 vector<double> &br) {
niklas.baumgarten's avatar
niklas.baumgarten committed
109
110
111
112
113
    br.reserve(2 * ar.size() - 1);
    br.insert(br.end(), ar.begin(), ar.end());
    br.insert(br.end(), ac.rbegin(), ac.rend() - 1);
}

niklas.baumgarten's avatar
naming    
niklas.baumgarten committed
114
void CirculantEmbedding1D::computeEV(vector<double> &br, vector<double> &ev) {
niklas.baumgarten's avatar
niklas.baumgarten committed
115
116
117
118
    fft(br, ev);
}

void CirculantEmbedding1D::padding(vector<double> &ar, vector<double> &ac) {
119
    double delta = (domain_end - domain_start) / (numberOfCells - 1);
niklas.baumgarten's avatar
niklas.baumgarten committed
120
121

    // Add 50% more points
122
    int N_padded = numberOfCells + int(numberOfCells / 2);
niklas.baumgarten's avatar
niklas.baumgarten committed
123
124
125
126
127
128
129
130
131
132

    double extended_domain_start = domain_start;
    double extended_domain_end = delta * (N_padded - 1);

    ar.resize(N_padded);
    ac.resize(N_padded);

    const vector<double> &coordinate1 = linspace(extended_domain_start,
                                                 extended_domain_end, N_padded);
    double y_rows = coordinate1.front();
133
    for (int i = numberOfCells; i < N_padded; i++) {
niklas.baumgarten's avatar
niklas.baumgarten committed
134
135
136
137
138
139
140
141
142
143
144
145
        double x_rows = coordinate1[i];
        double tau_rows = x_rows - y_rows;

        double x_columns = coordinate1.front();
        double y_columns = coordinate1[i];
        double tau_columns = x_columns - y_columns;

        ar[i] = covariance_fct(&tau_rows);
        ac[i] = covariance_fct(&tau_columns);
    }
}

niklas.baumgarten's avatar
naming    
niklas.baumgarten committed
146
vector<double> CirculantEmbedding1D::getSqrtEV() {
niklas.baumgarten's avatar
niklas.baumgarten committed
147
148
    vector<double> sqrt_ev_return, ar, ac, br, ev;

niklas.baumgarten's avatar
naming    
niklas.baumgarten committed
149
150
151
    createCovarianceMatrix(ar, ac);
    embedCovarianceMatrix(ar, ac, br);
    computeEV(br, ev);
niklas.baumgarten's avatar
niklas.baumgarten committed
152

153
    for (double &i : ev) {
154
155
        if (i > 0) continue;
        else if (i < 0 && i > evtol) i = 0.0;
niklas.baumgarten's avatar
niklas.baumgarten committed
156
157
158
159
160
161
162
163
164
        else padding(ar, ac);
    }

    sqrt_ev_return.resize(ev.size());
    for (int i = 0; i < ev.size(); i++)
        sqrt_ev_return[i] = sqrt(ev[i]);
    return sqrt_ev_return;
}

165
166
167
void CirculantEmbedding2D::generateFineSample(SampleID id,
                                              Vector *&fineSample,
                                              Vector *&coarseSample) {
niklas.baumgarten's avatar
naming    
niklas.baumgarten committed
168
    if (internalCounter == 0)
169
170
        fineComplexField = generateField();

171
172
173
    for (cell c = fineSample->cells(); c != fineSample->cells_end(); c++) {
        int i = floor(c()[0] * sqrt(fineSample->size()));
        int j = floor(c()[1] * sqrt(fineSample->size()));
174
175
        if (fieldType == "Gaussian") {
            if (internalCounter == 0)
176
                fineSample->operator()(c(), 0) = fineComplexField[i][j].real() + mean;
177
            else
178
                fineSample->operator()(c(), 0) = fineComplexField[i][j].imag() + mean;
179
180
        } else if (fieldType == "LogNormal") {
            if (internalCounter == 0)
181
                fineSample->operator()(c(), 0) = exp(fineComplexField[i][j].real() + mean);
182
            else
183
                fineSample->operator()(c(), 0) = exp(fineComplexField[i][j].imag() + mean);
184
        } else Exit("Has to be Gaussian or LogNormal")
niklas.baumgarten's avatar
niklas.baumgarten committed
185
    }
niklas.baumgarten's avatar
naming    
niklas.baumgarten committed
186
187
    internalCounter += 1;
    internalCounter %= 2;
188

niklas.baumgarten's avatar
niklas.baumgarten committed
189
    if (plotting)
190
        plotter->PlotVector("kappa", *fineSample,
niklas.baumgarten's avatar
niklas.baumgarten committed
191
                            1, l, "CellData");
192
193
}

194
195
196
197
198
199
void CirculantEmbedding2D::generateCoarseSample(SampleID id,
                                                Vector *&fineSample,
                                                Vector *&coarseSample) {
    *coarseSample = 0;
    for (cell c = coarseSample->cells(); c != coarseSample->cells_end(); ++c) {
        row coarseRow = coarseSample->find_row(c());
200
        for (int k = 0; k < c.Children(); ++k) {
201
            row fineRow = fineSample->find_row(c.Child(k));
202
            for (int i = 0; i < coarseRow.n(); ++i)
203
204
                coarseSample->operator()(coarseRow, i) +=
                    fineSample->operator()(fineRow, i);
205
206
        }
        for (int i = 0; i < coarseRow.n(); ++i)
207
            coarseSample->operator()(coarseRow, i) /= c.Children();
208
    }
209

niklas.baumgarten's avatar
niklas.baumgarten committed
210
    if (plotting)
211
        plotter->PlotVector("kappa", *coarseSample,
niklas.baumgarten's avatar
niklas.baumgarten committed
212
                            1, l - 1, "CellData");
niklas.baumgarten's avatar
niklas.baumgarten committed
213
214
}

215
216
217
ComplexField2D CirculantEmbedding2D::generateField() {
    ComplexField2D Xf;
    ComplexField2D Z;
218
219
    Z.resize(numberOfCellsEmbedded[1]);
    fineComplexField.clear();
220
    Xf.resize(numberOfCellsEmbedded[1]);
niklas.baumgarten's avatar
niklas.baumgarten committed
221
    double z1, z2;
222
223
    for (int i = 0; i < numberOfCellsEmbedded[1]; i++) {
        Z[i].resize(numberOfCellsEmbedded[0]);
224
        Xf[i].resize(numberOfCellsEmbedded[0]);
225
        for (int j = 0; j < numberOfCellsEmbedded[0]; j++) {
niklas.baumgarten's avatar
niklas.baumgarten committed
226
227
228
229
230
            rnmg->get_normal_rn(z1, z2);
            Z[i][j] = {z1, z2};
        }
    }

231
    ComplexField2D product;
232
233
234
235
236
    product.resize(numberOfCellsEmbedded[1]);
    for (int i = 0; i < numberOfCellsEmbedded[1]; i++) {
        product[i].resize(numberOfCellsEmbedded[0]);
        for (int j = 0; j < numberOfCellsEmbedded[0]; j++) {
            product[i][j] = sqrtEigenvalues[i][j] * Z[i][j];
niklas.baumgarten's avatar
niklas.baumgarten committed
237
238
        }
    }
239
    fft2(product, Xf);
240
241
242
    double divisor = sqrt(numberOfCellsEmbedded[0] * numberOfCellsEmbedded[1]);
    for (int i = 0; i < numberOfCellsEmbedded[0]; i++)
        for (int j = 0; j < numberOfCellsEmbedded[1]; j++)
243
244
            Xf[i][j] /= divisor;
    return Xf;
niklas.baumgarten's avatar
niklas.baumgarten committed
245
246
}

247
248
void CirculantEmbedding2D::createCovarianceMatrix(vector<vector<double>> &ar,
                                                  vector<vector<double>> &ac) {
249
250
251
252
    const vector<double>
        &coordinate1 = linspace(domain1_start, domain1_end, numberOfCells[0]);
    const vector<double>
        &coordinate2 = linspace(domain2_start, domain2_end, numberOfCells[1]);
niklas.baumgarten's avatar
niklas.baumgarten committed
253

254
255
    ar.resize(numberOfCells[1]);
    ac.resize(numberOfCells[1]);
niklas.baumgarten's avatar
niklas.baumgarten committed
256
257

    double y_rows[2] = {coordinate1.front(), coordinate2.front()};
258
259
260
261
    for (int index_c2 = 0; index_c2 < numberOfCells[1]; index_c2++) {
        ar[index_c2].resize(numberOfCells[0]);
        ac[index_c2].resize(numberOfCells[0]);
        for (int index_c1 = 0; index_c1 < numberOfCells[0]; index_c1++) {
niklas.baumgarten's avatar
niklas.baumgarten committed
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
            double c1 = coordinate1[index_c1];
            double c2 = coordinate2[index_c2];
            double x_rows[2] = {c1, c2};
            double tau_rows[2] = {x_rows[0] - y_rows[0],
                                  x_rows[1] - y_rows[1]};

            double x_columns[2] = {coordinate1.front(), c2};
            double y_columns[2] = {c1, coordinate2.front()};
            double tau_columns[2] = {x_columns[0] - y_columns[0],
                                     x_columns[1] - y_columns[1]};

            ac[index_c2][index_c1] = covariance_fct(tau_columns);
            ar[index_c2][index_c1] = covariance_fct(tau_rows);
        }
    }
}

279
280
void CirculantEmbedding2D::embedCovarianceMatrix(vector<vector<double>> &ar,
                                                 vector<vector<double>> &ac,
niklas.baumgarten's avatar
naming    
niklas.baumgarten committed
281
                                                 vector<vector<double>> &br) {
282
    int br_size = (int) (2 * ar.size() - 1);
niklas.baumgarten's avatar
niklas.baumgarten committed
283
284
285
286
287
288
289
290
291
292
293
294
    br.resize(br_size);
    for (int i = 0; i < br.size(); i++) {
        br[i].reserve(br_size);
        if (i < ar.size()) {
            br[i].insert(br[i].end(), ar[i].begin(), ar[i].end());
            br[i].insert(br[i].end(), ac[i].rbegin(), ac[i].rend() - 1);
        } else {
            br[i] = br[br_size - i];
        }
    }
}

295
296
void CirculantEmbedding2D::computeEV(vector<vector<double>> &br,
                                     vector<vector<double>> &ev) {
niklas.baumgarten's avatar
niklas.baumgarten committed
297
298
299
    fft2(br, ev);
}

300
301
void
CirculantEmbedding2D::padding(vector<vector<double>> &ar, vector<vector<double>> &ac) {
302
303
    double delta1 = (domain1_end - domain1_start) / (numberOfCells[0] - 1);
    double delta2 = (domain2_end - domain2_start) / (numberOfCells[1] - 1);
niklas.baumgarten's avatar
niklas.baumgarten committed
304
305

    // Add 50% more points
306
307
    int N_padded1 = numberOfCells[0] + int(numberOfCells[0] / 2);
    int N_padded2 = numberOfCells[1] + int(numberOfCells[1] / 2);
niklas.baumgarten's avatar
niklas.baumgarten committed
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324

    double extended_domain1_start = domain1_start;
    double extended_domain1_end = delta1 * (N_padded1 - 1);
    double extended_domain2_start = domain2_start;
    double extended_domain2_end = delta2 * (N_padded2 - 1);

    ar.resize(N_padded2);
    ac.resize(N_padded2);

    const vector<double> &coordinate1 = linspace(extended_domain1_start,
                                                 extended_domain1_end, N_padded1);
    const vector<double> &coordinate2 = linspace(extended_domain2_start,
                                                 extended_domain2_end, N_padded2);
    double y_rows[2] = {coordinate1.front(), coordinate2.front()};
    for (int index_c2 = 0; index_c2 < N_padded2; index_c2++) {
        ar[index_c2].resize(N_padded1);
        ac[index_c2].resize(N_padded1);
325
        for (int index_c1 = numberOfCells[0]; index_c1 < N_padded1; index_c1++) {
niklas.baumgarten's avatar
niklas.baumgarten committed
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
            double c1 = coordinate1[index_c1];
            double c2 = coordinate2[index_c2];
            double x_rows[2] = {c1, c2};
            double tau_rows[2] = {x_rows[0] - y_rows[0],
                                  x_rows[1] - y_rows[1]};

            double x_columns[2] = {coordinate1.front(), c2};
            double y_columns[2] = {c1, coordinate2.front()};
            double tau_columns[2] = {x_columns[0] - y_columns[0],
                                     x_columns[1] - y_columns[1]};

            ac[index_c2][index_c1] = covariance_fct(tau_columns);
            ar[index_c2][index_c1] = covariance_fct(tau_rows);
        }
    }
341
342
    for (int index_c2 = numberOfCells[1]; index_c2 < N_padded2; index_c2++) {
        for (int index_c1 = 0; index_c1 < numberOfCells[0]; index_c1++) {
niklas.baumgarten's avatar
niklas.baumgarten committed
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
            double c1 = coordinate1[index_c1];
            double c2 = coordinate2[index_c2];
            double x_rows[2] = {c1, c2};
            double tau_rows[2] = {x_rows[0] - y_rows[0],
                                  x_rows[1] - y_rows[1]};

            double x_columns[2] = {coordinate1.front(), c2};
            double y_columns[2] = {c1, coordinate2.front()};
            double tau_columns[2] = {x_columns[0] - y_columns[0],
                                     x_columns[1] - y_columns[1]};

            ac[index_c2][index_c1] = covariance_fct(tau_columns);
            ar[index_c2][index_c1] = covariance_fct(tau_rows);
        }
    }
}

niklas.baumgarten's avatar
naming    
niklas.baumgarten committed
360
vector<vector<double>> CirculantEmbedding2D::getSqrtEV() {
niklas.baumgarten's avatar
niklas.baumgarten committed
361
362
    vector<vector<double>> sqrt_ev_return, ar, ac, br, ev;

niklas.baumgarten's avatar
naming    
niklas.baumgarten committed
363
364
365
    createCovarianceMatrix(ar, ac);
    embedCovarianceMatrix(ar, ac, br);
    computeEV(br, ev);
niklas.baumgarten's avatar
niklas.baumgarten committed
366
367
368
369

    for (int i = 0; i < ev[0].size(); i++) {
        for (int j = 0; j < ev[1].size(); j++) {
            if (ev[i][j] > 0) continue;
niklas.baumgarten's avatar
naming    
niklas.baumgarten committed
370
            else if (ev[i][j] < 0 && ev[i][j] > evtol) ev[i][j] = 0.0;
niklas.baumgarten's avatar
niklas.baumgarten committed
371
372
373
374
375
376
377
378
379
380
381
382
383
            else padding(ar, ac);
        }
    }

    sqrt_ev_return.resize(ev[0].size());
    for (int i = 0; i < ev[0].size(); i++) {
        sqrt_ev_return[i].resize(ev[1].size());
        for (int j = 0; j < ev[1].size(); j++)
            sqrt_ev_return[i][j] = sqrt(ev[i][j]);
    }
    return sqrt_ev_return;
}