MonteCarlo.h 3.95 KB
Newer Older
1
2
3
#ifndef MLMC_MC_HPP
#define MLMC_MC_HPP

4
#include "m++.h"
5
#include "stochastics/StochasticField.h"
6
#include "MonteCarloLogger.h"
niklas.baumgarten's avatar
niklas.baumgarten committed
7
8
9
#include <sys/stat.h>

#include <utility>
10

11

12
13
class MonteCarlo {
public:
niklas.baumgarten's avatar
niklas.baumgarten committed
14
    int plotting = 0;
niklas.baumgarten's avatar
niklas.baumgarten committed
15
16
17
    Plot *coarsePlot = nullptr;
    Plot *finePlot = nullptr;
    MonteCarloLogger *logger = nullptr;
18

niklas.baumgarten's avatar
niklas.baumgarten committed
19
20
    int l = 0;
    int pLevel = 0;
21
    bool onlyFineLevel = true;
22
23
    string functional = "L2";

niklas.baumgarten's avatar
niklas.baumgarten committed
24
25
    Meshes *meshes = nullptr;
    StochasticField *stochasticField = nullptr;
26
27
28

    int M = 0, dM = 0;

29
    double sumCost = 0.0, avgCost = 0.0;
30
    double sumY1 = 0.0, avgY = 0.0;
31
32
33
    double sumY2 = 0.0, avgY2 = 0.0, varY = 0.0;
    double sumY3 = 0.0, avgY3 = 0.0;
    double sumY4 = 0.0, avgY4 = 0.0;
34
    double sumQ1 = 0.0, avgQ = 0.0;
35
    double sumQ2 = 0.0, avgQ2 = 0.0, varQ = 0.0;
36
37
    double kurtosis = 0.0;

38
39
40
41
42
    MonteCarlo(int l, int dM, bool onlyFineLevel,
               Meshes *meshes, StochasticField *stochFields) :
        l(l), dM(dM), onlyFineLevel(onlyFineLevel),
        meshes(meshes), pLevel(meshes->pLevel()),
        stochasticField(stochFields) {
43

44
45
        config.get("MCPlotting", plotting);
        config.get("functional", functional);
46

niklas.baumgarten's avatar
niklas.baumgarten committed
47
48
        coarsePlot = new Plot((*meshes)[l - 1], meshes->dim());
        finePlot = new Plot((*meshes)[l], meshes->dim());
49
        logger = new MonteCarloLogger();
niklas.baumgarten's avatar
niklas.baumgarten committed
50
    }
51

52
    ~MonteCarlo() {
niklas.baumgarten's avatar
niklas.baumgarten committed
53
54
        delete coarsePlot;
        delete finePlot;
55
56
        delete logger;
    }
57

niklas.baumgarten's avatar
niklas.baumgarten committed
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
    void Method() {
        logger->StartMethodMsg();
        if (!onlyFineLevel)
            logger->LogMsgv1("l: " + to_string(l)
                                 + " dM: " + to_string(dM));
        else
            logger->LogMsgv1("l: " + to_string(l)
                                 + " dM: " + to_string(dM)
                                 + " onlyFineLevel: True");

        method();
        updateAvg();
        updateStatistic();

        logger->LogMsgv1("|E[Y]|: " + to_string(avgY)
                             + " V[Y]: " + to_string(varY));
        logger->EndMethodMsg();
    }
76
77

protected:
niklas.baumgarten's avatar
niklas.baumgarten committed
78
79
80
81
82
83
84
85
    virtual void initialize() = 0;

    virtual void method() = 0;

    virtual void solvePDE(int l, bool fineLevel,
                          double &valueQ, double &cost,
                          Vector &solution) = 0;

86
87
    void updateSums(double fineCost, double diffQ, double fineQ) {
        sumCost += fineCost;
88
89
90
91
92
93
        sumQ1 += fineQ;
        sumQ2 += fineQ * fineQ;
        sumY1 += diffQ;
        sumY2 += diffQ * diffQ;
        sumY3 += diffQ * diffQ * diffQ;
        sumY4 += diffQ * diffQ * diffQ * diffQ;
94
    }
95

96
97
98
    void updateAvg() {
        M += dM;
        dM = 0;
99
        avgCost = sumCost / M;
100
        avgY = abs(sumY1 / M);
101
102
103
        avgY2 = sumY2 / M;
        avgY3 = sumY3 / M;
        avgY4 = sumY4 / M;
104
        avgQ = abs(sumQ1 / M);
105
        avgQ2 = sumQ2 / M;
106
    }
107

108
    void updateStatistic() {
109
110
111
112
        varY = fmax(avgY2 - avgY * avgY, 1e-10);
        varQ = fmax(avgQ2 - avgQ * avgQ, 1e-10);
        kurtosis = (avgY4 - 4.0 * avgY3 * avgY + 6.0 * avgY2 * avgY * avgY -
            3.0 * avgY * avgY * avgY * avgY) / varY / varY;
113
    }
niklas.baumgarten's avatar
niklas.baumgarten committed
114

niklas.baumgarten's avatar
niklas.baumgarten committed
115
116
117
118
119
    string mkSampleDirName(int m, bool fineLevel) {
        if (fineLevel)
            return buildName(m, "sample");
        else
            return buildName(m, "sample", "coarse");
niklas.baumgarten's avatar
niklas.baumgarten committed
120
121
    }

niklas.baumgarten's avatar
niklas.baumgarten committed
122
    static void mkSampleDir(int m, const string &sampleDirName) {
niklas.baumgarten's avatar
niklas.baumgarten committed
123
        string pathToData = "data/vtk/";
niklas.baumgarten's avatar
niklas.baumgarten committed
124
125
126
        char *sampleDirNameCharArr =
            const_cast<char *> ((pathToData.append(sampleDirName)).c_str());
        int status = mkdir(sampleDirNameCharArr, 0777);
niklas.baumgarten's avatar
niklas.baumgarten committed
127
128
    }

niklas.baumgarten's avatar
niklas.baumgarten committed
129
130
    string buildName(int m, const string &name, const string &suffix = "") {
        string returnName = name;
niklas.baumgarten's avatar
niklas.baumgarten committed
131
        if (!suffix.empty())
niklas.baumgarten's avatar
niklas.baumgarten committed
132
            returnName = returnName.append("_").append(suffix);
niklas.baumgarten's avatar
niklas.baumgarten committed
133
        if (l != -1)
niklas.baumgarten's avatar
niklas.baumgarten committed
134
            returnName = returnName.append("_").append(to_string(l));
niklas.baumgarten's avatar
niklas.baumgarten committed
135
        if (m != -1)
niklas.baumgarten's avatar
niklas.baumgarten committed
136
137
            returnName = returnName.append("_").append(to_string(m));
        return returnName;
niklas.baumgarten's avatar
niklas.baumgarten committed
138
    }
139
140
141
};

#endif //MLMC_MC_HPP