MultilevelMonteCarlo.C 9.94 KB
Newer Older
1
#include "MultilevelMonteCarlo.h"
niklas.baumgarten's avatar
niklas.baumgarten committed
2

3

niklas.baumgarten's avatar
niklas.baumgarten committed
4
5
using namespace std;

6
MonteCarlo *MultilevelMonteCarlo::getMonteCarlo(int l, int dM, bool onlyFineLevel) {
7
    auto ellipticAssemble = dynamic_cast<EllipticAssemble *>(assemble);
8
    auto transportAssemble = dynamic_cast<DGTransportAssemble *>(assemble);
9
    if (ellipticAssemble != nullptr)
10
        return new MonteCarloElliptic(l, dM, onlyFineLevel, meshes,
11
                                      stochasticField, ellipticAssemble, plots);
12
    if (transportAssemble != nullptr)
13
        return new MonteCarloTransport(l, dM, onlyFineLevel, meshes, plots,
14
                                       stochasticField, transportAssemble);
niklas.baumgarten's avatar
niklas.baumgarten committed
15
    Exit("Wrong Assembly Class")
16
17
}

18
19
20
21
22
23
void MultilevelMonteCarlo::initializeMapMonteCarlo() {
    clearMapMonteCarlo();
    for (unsigned long i = 0; i < initLevels.size(); i++) {
        bool onlyFineLevel = (i == 0) || (mcOnly);
        int l = initLevels[i], M = initSampleAmount[i];
        mapMonteCarlo[l] = getMonteCarlo(l, M, onlyFineLevel);
niklas.baumgarten's avatar
niklas.baumgarten committed
24
25
26
    }
}

27
void MultilevelMonteCarlo::initializeValues() {
28
29
30
    levels = {}, numsamples = {};
    value = 0.0, cost = 0.0;
    alpha = 0.0, beta = 0.0, gamma = 0.0;
niklas.baumgarten's avatar
niklas.baumgarten committed
31
    numErr = 0.0, statErr = 0.0, totalErr = 0.0;
32
33
}

34
35
void MultilevelMonteCarlo::clearMapMonteCarlo() {
    for (auto &iter: mapMonteCarlo) {
niklas.baumgarten's avatar
niklas.baumgarten committed
36
37
        delete iter.second;
    }
38
    mapMonteCarlo.clear();
niklas.baumgarten's avatar
niklas.baumgarten committed
39
40
}

41
void MultilevelMonteCarlo::Method() {
42
    logger->StartMethod("Start MLMC Method", verbose);
43
44
    for (auto &mc : mapMonteCarlo)
        mc.second->Method();
45
    logger->EndMethod(verbose);
46
47
}

48
void MultilevelMonteCarlo::Method(const double eps) {
49
50
    logger->StartMethod("Start MLMC Method", verbose);
    logger->InnerMsg("eps: " + to_string(eps), verbose);
niklas.baumgarten's avatar
niklas.baumgarten committed
51
52
    bool converged = false;
    while (!converged) {
53
        for (auto &mc : mapMonteCarlo)
54
            if (mc.second->dM > 0)
55
                mc.second->Method();
niklas.baumgarten's avatar
niklas.baumgarten committed
56
        estimateExponents();
57
        updateSampleAmount(eps);
58
        if (checkForNewLevel()) {
niklas.baumgarten's avatar
niklas.baumgarten committed
59
60
            estimateNumericalError();
            if (numErr > eps) {
61
                appendLevel(eps);
niklas.baumgarten's avatar
niklas.baumgarten committed
62
63
64
            } else {
                estimateStatisticalError();
                totalErr = statErr + numErr;
niklas.baumgarten's avatar
niklas.baumgarten committed
65
                converged = true;
niklas.baumgarten's avatar
niklas.baumgarten committed
66
            }
niklas.baumgarten's avatar
niklas.baumgarten committed
67
68
        }
    }
69
    wrapUpResults(value, cost, levels, numsamples);
70
    logger->EndMethod(verbose);
niklas.baumgarten's avatar
niklas.baumgarten committed
71
72
}

73
void MultilevelMonteCarlo::Method(const vector<double> &epsLst) {
74
    logger->StartMethod("Start MLMC Method over epsilon", verbose);
75
    for (auto &eps : epsLst) {
76
77
78
79
80
81
82
        initializeMapMonteCarlo();
        initializeValues();
        Method(eps);
        levelsOverEpsilon.push_back(levels);
        sampleAmountOverEpsilon.push_back(numsamples);
        valueOverEpsilon.push_back(value);
        costOverEpsilon.push_back(cost);
niklas.baumgarten's avatar
niklas.baumgarten committed
83
    }
84
    logger->EndMethod(verbose);
85
86
}

niklas.baumgarten's avatar
niklas.baumgarten committed
87
void MultilevelMonteCarlo::estimateExponents(bool excludeBaseLevel) {
88
    vector<double> levelList, avgsY, varsY, costs;
89
    auto iter = mapMonteCarlo.begin();
niklas.baumgarten's avatar
niklas.baumgarten committed
90
91
92
    if (excludeBaseLevel)
        iter++;

93
    for (; iter != mapMonteCarlo.end(); iter++) {
94
95
96
        levelList.push_back((double) iter->first);
        avgsY.push_back(-log2(iter->second->avgY));
        varsY.push_back(-log2(iter->second->varY));
97
        costs.push_back(log2(iter->second->avgCost));
niklas.baumgarten's avatar
niklas.baumgarten committed
98
    }
99

niklas.baumgarten's avatar
niklas.baumgarten committed
100
101
102
103
104
105
    pair<double, double> res = linear_fit(levelList, avgsY);
    alpha = res.first;
    res = linear_fit(levelList, varsY);
    beta = res.first;
    res = linear_fit(levelList, costs);
    gamma = res.first;
106

107
    logger->InnerMsg("alpha: " + to_string(alpha) +
108
        " beta: " + to_string(beta) +
109
        " gamma: " + to_string(gamma), verbose);
niklas.baumgarten's avatar
niklas.baumgarten committed
110
111
}

112
void MultilevelMonteCarlo::updateSampleAmount(const double &eps) {
113
    int optimalM;
niklas.baumgarten's avatar
niklas.baumgarten committed
114
    double factor = 0.0;
115
116
117
    for (auto &mc : mapMonteCarlo)
        factor += sqrt(mc.second->varY * mc.second->avgCost);
    for (auto &mc : mapMonteCarlo) {
118
        optimalM = (int) (ceil(2 * pow(eps, -2) * factor *
119
            sqrt(mc.second->varY / mc.second->avgCost)));
120
121
        if (optimalM == 1) optimalM++; // Hack
        mapMonteCarlo[mc.first]->dM = optimalM - mapMonteCarlo[mc.first]->M;
niklas.baumgarten's avatar
niklas.baumgarten committed
122
    }
123
124

    string msg = "dM:";
125
    for (auto &mc:mapMonteCarlo) msg += " " + to_string(mc.second->dM);
126
    logger->InnerMsg(msg, verbose);
127
128
129
130
}

bool MultilevelMonteCarlo::checkForNewLevel() {
    bool newLevel = true;
131
    for (auto &mc : mapMonteCarlo) {
132
133
134
135
136
137
        if (mc.second->dM > 0) {
            newLevel = false;
            break;
        }
    }
    return newLevel;
niklas.baumgarten's avatar
niklas.baumgarten committed
138
139
}

niklas.baumgarten's avatar
niklas.baumgarten committed
140
141
void MultilevelMonteCarlo::estimateNumericalError() {
    numErr = 0.0;
142
143
144
    unsigned long exponent = mapMonteCarlo.size() - 2;
    for (auto &mc : mapMonteCarlo) {
        if (mc.first == mapMonteCarlo.begin()->first) continue;
niklas.baumgarten's avatar
niklas.baumgarten committed
145
146
        numErr = max(mc.second->avgY / (pow(2.0, alpha) - 1) / (pow(2.0, exponent)),
                     numErr);
niklas.baumgarten's avatar
niklas.baumgarten committed
147
148
        exponent--;
    }
149

niklas.baumgarten's avatar
niklas.baumgarten committed
150
151
152
153
154
155
156
157
158
    logger->InnerMsg("estimated numerical error: " + to_string(numErr), verbose);
}

void MultilevelMonteCarlo::estimateStatisticalError() {
    statErr = 0.0;
    for (auto &mc : mapMonteCarlo) {
        statErr += mc.second->varY / mc.second->M;
    }
    statErr = sqrt(statErr);
159
160
}

161
void MultilevelMonteCarlo::appendLevel(double eps) {
162
163
    int L = mapMonteCarlo.rbegin()->first + 1;
    if (L > maxLevel) {
164
165
        Exit ("Maximum level has been reached.")
    } else {
166
        logger->InnerMsg("new level: " + to_string(L), verbose);
167
168
169
170
        mapMonteCarlo[L] = getMonteCarlo(L, 0, false);
        mapMonteCarlo[L]->varY = mapMonteCarlo[L - 1]->varY / pow(2.0, beta);
        mapMonteCarlo[L]->avgCost = mapMonteCarlo[L - 1]->avgCost * pow(2.0, gamma);
        updateSampleAmount(eps);
171
172
173
    }
}

174
175
176
177
void MultilevelMonteCarlo::wrapUpResults(double &val,
                                         double &mlmcCost,
                                         vector<int> &listLevels,
                                         vector<int> &samples) {
178
179
180
181
182
    int L = mapMonteCarlo.rbegin()->first;
    for (auto &mc : mapMonteCarlo) {
        val += mc.second->avgY;
        mlmcCost += mc.second->sumCost;
        listLevels.push_back(mc.first);
niklas.baumgarten's avatar
niklas.baumgarten committed
183
184
185
186
        samples.push_back(mc.second->M);
    }
}

niklas.baumgarten's avatar
niklas.baumgarten committed
187
188
189
190
191
192
193
194
195
void MultilevelMonteCarlo::PrintInfo() {
    Vout(1) << endl
            << "MLMC Info:" << endl
            << "  initLevels:       " << vec2str(initLevels) << endl
            << "  initSample:       " << vec2str(initSampleAmount) << endl
            << "  maxLevel:         " << maxLevel << endl
            << endl;
}

196
void MultilevelMonteCarlo::ShowMCTableHead() {
197
    mout << endl
198
199
200
201
202
         << " l     M    E[Qf-Qc]       E[Qf]    V[Qf-Qc]"
            "       V[Qf]    kurtosis        cost"
         << endl << line() << endl;
}

203
void MultilevelMonteCarlo::ShowResultsMLMCRunHead() {
204
205
206
207
208
209
    mout << endl
         << "   eps       value   MLMC cost             l"
            "                                   M"
         << endl << line() << endl;
}

210
void MultilevelMonteCarlo::ShowTableBottom() {
211
212
213
214
215
216
217
    mout << line() << endl;
}

string MultilevelMonteCarlo::line() {
    return "----------------------------------------"
           "----------------------------------------";
}
218

219
220
221
void MultilevelMonteCarlo::ShowMCTable() {
    ShowMCTableHead();
    for (auto &mc : mapMonteCarlo) {
222
223
        int l = mc.first;
        mout << setw(2) << l << setw(6) << mc.second->M
224
225
226
227
             << setw(12) << mc.second->avgY
             << setw(12) << mc.second->avgQ << setw(12) << mc.second->varY
             << setw(12) << mc.second->varQ << setw(12) << mc.second->kurtosis
             << setw(12) << mc.second->avgCost << endl;
228
    }
229
    ShowTableBottom();
230
231
}

232
233
void MultilevelMonteCarlo::ShowKurtosisWarning() {
    for (auto &mc : mapMonteCarlo) {
234
235
236
237
238
239
240
241
242
        if (mc.second->kurtosis > 100.0) {
            mout << endl << "WARNING: kurtosis on level " << setw(2) << mc.first
                 << " is " << setw(6) << mc.second->kurtosis << endl
                 << "indicates MLMC correction is dominated by a few rare paths!"
                 << endl;
        }
    }
}

243
void MultilevelMonteCarlo::ShowExponentResults() {
niklas.baumgarten's avatar
niklas.baumgarten committed
244
    estimateExponents();
niklas.baumgarten's avatar
niklas.baumgarten committed
245
246
247
248
249
250
251
252
253
254
255
256
    mout << endl << "alpha   = " << setw(10) << alpha
         << "   (exponent for the weak convergence of the FEM)"
         << endl << " beta   = " << setw(10) << beta
         << "   (exponent for the decay of the variance)"
         << endl << "gamma   = " << setw(10) << gamma
         << "   (exponent for the growth of the cost)"
         << endl << "numErr  = " << setw(10) << numErr
         << "   (estimated numerical error)"
         << endl << "statErr = " << setw(10) << statErr
         << "   (estimated statistical error)"
         << endl << "error   = " << setw(10) << totalErr
         << "   (estimated total error)"
257
258
259
         << endl;
}

260
261
void MultilevelMonteCarlo::ShowResultsMLMCRun(double eps) {
    ShowResultsMLMCRunHead();
262
263
264
265
    string level_str = vec2str(levels);
    string sample_str = vec2str(numsamples);
    mout << setw(6) << eps << setw(12) << value << setw(12) << cost
         << setw(14) << level_str << setw(36) << sample_str << endl;
266
    ShowTableBottom();
267
268
}

269
void MultilevelMonteCarlo::ShowResultsOverEpsilon(const vector<double> &epsLst) {
270
    ShowResultsMLMCRunHead();
271
    int ctr = 0;
272
    for (auto &eps: epsLst) {
273
        string level_str, sample_str;
274
275
276
277
278
        mout << setw(6) << eps << setw(12) << valueOverEpsilon[ctr]
             << setw(12) << costOverEpsilon[ctr];
        for (unsigned long i = 0; i < levelsOverEpsilon[ctr].size(); i++) {
            if (i != levelsOverEpsilon[ctr].size() - 1)
                level_str += to_string(levelsOverEpsilon[ctr][i]) + " ";
279
            else
280
                level_str += to_string(levelsOverEpsilon[ctr][i]);
281
        }
niklas.baumgarten's avatar
niklas.baumgarten committed
282

283
        mout << setw(14) << level_str;
284
285
286
        for (unsigned long i = 0; i < sampleAmountOverEpsilon[ctr].size(); i++) {
            if (i != sampleAmountOverEpsilon[ctr].size() - 1)
                sample_str += to_string(sampleAmountOverEpsilon[ctr][i]) + " ";
287
            else
288
                sample_str += to_string(sampleAmountOverEpsilon[ctr][i]);
289
290
291
        }
        mout << setw(36) << sample_str << endl;
        ctr++;
niklas.baumgarten's avatar
niklas.baumgarten committed
292
    }
293
    ShowTableBottom();
niklas.baumgarten's avatar
niklas.baumgarten committed
294
}