CirculantEmbedding.cpp 14.1 KB
Newer Older
1
#include "CirculantEmbedding.hpp"
2
3


4
void CirculantEmbedding1D::generateFineSample(const SampleID &id,
5
6
                                              Vector *&fineSample,
                                              Vector *&coarseSample) {
7
    fineSample = new Vector(cellMGraphs[id.level.mGfine]);
8

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

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

niklas.baumgarten's avatar
niklas.baumgarten committed
29
    if (coarseSample) delete coarseSample;
30
    coarseSample = nullptr;
31
32
}

33
void CirculantEmbedding1D::generateCoarseSample(const SampleID &id,
34
35
                                                Vector *&fineSample,
                                                Vector *&coarseSample) {
36
    coarseSample = new Vector(cellMGraphs[id.level.mGcoarse]);
37
38
39
40

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

niklas.baumgarten's avatar
niklas.baumgarten committed
51
    if (fineSample) delete fineSample;
52
    fineSample = nullptr;
niklas.baumgarten's avatar
niklas.baumgarten committed
53
54
}

55
ComplexSequence1D CirculantEmbedding1D::generateField(const SampleID &id) {
niklas.baumgarten's avatar
niklas.baumgarten committed
56
    generator.DrawSample(id);
57
58
    ComplexSequence1D normalDist = generator.EvalSample();
    ComplexSequence1D Xf;
59
    fineComplexField.clear();
60
    Xf.resize(numberOfCellsEmbedded);
niklas.baumgarten's avatar
niklas.baumgarten committed
61

62
    ComplexSequence1D product;
63
64
    product.resize(numberOfCellsEmbedded);
    for (int i = 0; i < numberOfCellsEmbedded; i++) {
niklas.baumgarten's avatar
niklas.baumgarten committed
65
        product[i] = sqrtEigenvalues[i] * normalDist[i];
niklas.baumgarten's avatar
niklas.baumgarten committed
66
    }
67
    fft(product, Xf);
68
69
    double divisor = sqrt(numberOfCellsEmbedded);
    for (int i = 0; i < numberOfCellsEmbedded; i++)
70
71
        Xf[i] /= divisor;
    return Xf;
niklas.baumgarten's avatar
niklas.baumgarten committed
72
73
}

74
void CirculantEmbedding1D::createCovarianceMatrix(ToeplitzRow &ar, ToeplitzColumn &ac) {
75
76
77
    const vector<double> &coordinate1 = linspace(domain_start, domain_end, numberOfCells);
    ar.resize(numberOfCells);
    ac.resize(numberOfCells);
niklas.baumgarten's avatar
niklas.baumgarten committed
78
79
    double y_rows = coordinate1.front();

80
    for (int index_c = 0; index_c < numberOfCells; index_c++) {
niklas.baumgarten's avatar
niklas.baumgarten committed
81
82
83
84
85
86
87
88
89
90
91
92
93
        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);
    }
}

94
95
96
void CirculantEmbedding1D::embedCovarianceMatrix(const ToeplitzRow &ar,
                                                 const ToeplitzColumn &ac,
                                                 CirculantRow &br) {
niklas.baumgarten's avatar
niklas.baumgarten committed
97
98
99
100
101
    br.reserve(2 * ar.size() - 1);
    br.insert(br.end(), ar.begin(), ar.end());
    br.insert(br.end(), ac.rbegin(), ac.rend() - 1);
}

102
103
void CirculantEmbedding1D::computeEV(const CirculantRow &br,
                                     EigenValues1D &ev) {
niklas.baumgarten's avatar
niklas.baumgarten committed
104
105
106
    fft(br, ev);
}

107
108
void CirculantEmbedding1D::padding(ToeplitzRow &ar,
                                   ToeplitzColumn &ac) {
109
    double delta = (domain_end - domain_start) / (numberOfCells - 1);
niklas.baumgarten's avatar
niklas.baumgarten committed
110
111

    // Add 50% more points
112
    int N_padded = numberOfCells + int(numberOfCells / 2);
niklas.baumgarten's avatar
niklas.baumgarten committed
113
114
115
116
117
118
119
120
121
122

    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();
123
    for (int i = numberOfCells; i < N_padded; i++) {
niklas.baumgarten's avatar
niklas.baumgarten committed
124
125
126
127
128
129
130
131
132
133
134
135
        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);
    }
}

136
SqrtEigenValues1D CirculantEmbedding1D::computeSqrtEV() {
niklas.baumgarten's avatar
niklas.baumgarten committed
137
138
    vector<double> sqrt_ev_return, ar, ac, br, ev;

niklas.baumgarten's avatar
naming    
niklas.baumgarten committed
139
140
141
    createCovarianceMatrix(ar, ac);
    embedCovarianceMatrix(ar, ac, br);
    computeEV(br, ev);
niklas.baumgarten's avatar
niklas.baumgarten committed
142

143
    for (double &i : ev) {
144
145
        if (i > 0) continue;
        else if (i < 0 && i > evtol) i = 0.0;
niklas.baumgarten's avatar
niklas.baumgarten committed
146
147
148
149
150
151
152
153
154
        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;
}

155
void CirculantEmbedding2D::generateFineSample(const SampleID &id,
156
157
                                              Vector *&fineSample,
                                              Vector *&coarseSample) {
158
    fineSample = new Vector(cellMGraphs[id.level.mGfine]);
159

niklas.baumgarten's avatar
naming    
niklas.baumgarten committed
160
    if (internalCounter == 0)
161
        fineComplexField = generateField(id);
162

163
164
165
    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()));
166
167
        if (fieldType == "Gaussian") {
            if (internalCounter == 0)
168
                fineSample->operator()(c(), 0) = fineComplexField[i][j].real() + mean;
169
            else
170
                fineSample->operator()(c(), 0) = fineComplexField[i][j].imag() + mean;
171
172
        } else if (fieldType == "LogNormal") {
            if (internalCounter == 0)
173
174
                fineSample->operator()(c(), 0) =
                    exp(fineComplexField[i][j].real() + mean);
175
            else
176
177
                fineSample->operator()(c(), 0) =
                    exp(fineComplexField[i][j].imag() + mean);
178
        } else Exit("Has to be Gaussian or LogNormal")
niklas.baumgarten's avatar
niklas.baumgarten committed
179
    }
niklas.baumgarten's avatar
naming    
niklas.baumgarten committed
180
181
    internalCounter += 1;
    internalCounter %= 2;
182

183
184
    if (coarseSample) delete coarseSample;
    coarseSample = nullptr;
185
186
}

187
void CirculantEmbedding2D::generateCoarseSample(const SampleID &id,
188
189
                                                Vector *&fineSample,
                                                Vector *&coarseSample) {
190
    coarseSample = new Vector(cellMGraphs[id.level.mGcoarse]);
191

192
193
194
    *coarseSample = 0;
    for (cell c = coarseSample->cells(); c != coarseSample->cells_end(); ++c) {
        row coarseRow = coarseSample->find_row(c());
195
        for (int k = 0; k < c.Children(); ++k) {
196
            row fineRow = fineSample->find_row(c.Child(k));
197
            for (int i = 0; i < coarseRow.n(); ++i)
198
199
                coarseSample->operator()(coarseRow, i) +=
                    fineSample->operator()(fineRow, i);
200
201
        }
        for (int i = 0; i < coarseRow.n(); ++i)
202
            coarseSample->operator()(coarseRow, i) /= c.Children();
203
    }
204

205
206
    if (fineSample) delete fineSample;
    fineSample = nullptr;
niklas.baumgarten's avatar
niklas.baumgarten committed
207
208
}

209
210
211
ComplexSequence2D CirculantEmbedding2D::generateField(const SampleID &id) {
    ComplexSequence2D Xf;
    ComplexSequence2D Z;
212
213
    Z.resize(numberOfCellsEmbedded[1]);
    fineComplexField.clear();
214
    Xf.resize(numberOfCellsEmbedded[1]);
215
216
    for (int i = 0; i < numberOfCellsEmbedded[1]; i++) {
        Z[i].resize(numberOfCellsEmbedded[0]);
217
        Xf[i].resize(numberOfCellsEmbedded[0]);
218
        for (int j = 0; j < numberOfCellsEmbedded[0]; j++) {
219
            normalDist.DrawSample(id);
220
            Z[i][j] = normalDist.EvalSample();
niklas.baumgarten's avatar
niklas.baumgarten committed
221
222
223
        }
    }

224
    ComplexSequence2D product;
225
226
227
228
229
    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
230
231
        }
    }
232
    fft2(product, Xf);
233
234
235
    double divisor = sqrt(numberOfCellsEmbedded[0] * numberOfCellsEmbedded[1]);
    for (int i = 0; i < numberOfCellsEmbedded[0]; i++)
        for (int j = 0; j < numberOfCellsEmbedded[1]; j++)
236
237
            Xf[i][j] /= divisor;
    return Xf;
niklas.baumgarten's avatar
niklas.baumgarten committed
238
239
}

240
241
void CirculantEmbedding2D::createCovarianceMatrix(vector<vector<double>> &ar,
                                                  vector<vector<double>> &ac) {
242
243
244
245
    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
246

247
248
    ar.resize(numberOfCells[1]);
    ac.resize(numberOfCells[1]);
niklas.baumgarten's avatar
niklas.baumgarten committed
249
250

    double y_rows[2] = {coordinate1.front(), coordinate2.front()};
251
252
253
254
    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
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
            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);
        }
    }
}

272
273
274
void CirculantEmbedding2D::embedCovarianceMatrix(const BlockToeplitzRow &ar,
                                                 const BlockToeplitzColumn &ac,
                                                 BlockCirculantRow &br) {
275
    int br_size = (int) (2 * ar.size() - 1);
niklas.baumgarten's avatar
niklas.baumgarten committed
276
277
278
279
280
281
282
283
284
285
286
287
    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];
        }
    }
}

288
289
void CirculantEmbedding2D::computeEV(const BlockCirculantRow &br,
                                     EigenValues2D &ev) {
niklas.baumgarten's avatar
niklas.baumgarten committed
290
291
292
    fft2(br, ev);
}

293
294
void CirculantEmbedding2D::padding(BlockToeplitzRow &ar,
                                   BlockToeplitzColumn &ac) {
295
296
    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
297
298

    // Add 50% more points
299
300
    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
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317

    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);
318
        for (int index_c1 = numberOfCells[0]; index_c1 < N_padded1; index_c1++) {
niklas.baumgarten's avatar
niklas.baumgarten committed
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
            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);
        }
    }
334
335
    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
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
            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
353
vector<vector<double>> CirculantEmbedding2D::getSqrtEV() {
niklas.baumgarten's avatar
niklas.baumgarten committed
354
355
    vector<vector<double>> sqrt_ev_return, ar, ac, br, ev;

niklas.baumgarten's avatar
naming    
niklas.baumgarten committed
356
357
358
    createCovarianceMatrix(ar, ac);
    embedCovarianceMatrix(ar, ac, br);
    computeEV(br, ev);
niklas.baumgarten's avatar
niklas.baumgarten committed
359
360
361
362

    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
363
            else if (ev[i][j] < 0 && ev[i][j] > evtol) ev[i][j] = 0.0;
niklas.baumgarten's avatar
niklas.baumgarten committed
364
365
366
367
368
369
370
371
372
373
374
375
            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;
}