MonteCarlo.h 4.07 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
    Plot *coarsePlot = nullptr;
    Plot *finePlot = nullptr;
17
18
19

    int verbose = 0;
    IndentedLogger *logger = nullptr;
20

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

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

    int M = 0, dM = 0;

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

40
41
42
    double fineQ = 0.0, coarseQ = 0.0;
    double fineCost = 0.0, coarseCost = 0.0;

43
44
45
46
47
    MonteCarlo(int l, int dM, bool onlyFineLevel,
               Meshes *meshes, StochasticField *stochFields) :
        l(l), dM(dM), onlyFineLevel(onlyFineLevel),
        meshes(meshes), pLevel(meshes->pLevel()),
        stochasticField(stochFields) {
48

49
50
        config.get("MCPlotting", plotting);
        config.get("functional", functional);
51

niklas.baumgarten's avatar
niklas.baumgarten committed
52
53
        coarsePlot = new Plot((*meshes)[l - 1], meshes->dim());
        finePlot = new Plot((*meshes)[l], meshes->dim());
54
55
56

        config.get("MCVerbose", verbose);
        logger = IndentedLogger::GetInstance();
niklas.baumgarten's avatar
niklas.baumgarten committed
57
    }
58

niklas.baumgarten's avatar
niklas.baumgarten committed
59
    virtual ~MonteCarlo() {
niklas.baumgarten's avatar
niklas.baumgarten committed
60
61
        delete coarsePlot;
        delete finePlot;
62
    }
63

niklas.baumgarten's avatar
niklas.baumgarten committed
64
    void Method() {
65
        logger->StartMethod("Start MC Method", verbose);
niklas.baumgarten's avatar
niklas.baumgarten committed
66
        if (!onlyFineLevel)
67
            logger->InnerMsg("l: " + to_string(l) + " dM: " + to_string(dM), verbose);
niklas.baumgarten's avatar
niklas.baumgarten committed
68
        else
69
70
            logger->InnerMsg("l: " + to_string(l) + " dM: " + to_string(dM)
                                 + " onlyFineLevel: True", verbose);
niklas.baumgarten's avatar
niklas.baumgarten committed
71
72
73
74
75

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

76
77
78
        logger->InnerMsg("|E[Y]|: " + to_string(avgY) + " V[Y]: " + to_string(varY),
                       verbose);
        logger->EndMethod(verbose);
niklas.baumgarten's avatar
niklas.baumgarten committed
79
    }
80
81

protected:
niklas.baumgarten's avatar
niklas.baumgarten committed
82
83
84
85
86
87
88
89
    virtual void initialize() = 0;

    virtual void method() = 0;

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

90
91
    void updateSums(double fineCost, double diffQ, double fineQ) {
        sumCost += fineCost;
92
93
94
95
96
97
        sumQ1 += fineQ;
        sumQ2 += fineQ * fineQ;
        sumY1 += diffQ;
        sumY2 += diffQ * diffQ;
        sumY3 += diffQ * diffQ * diffQ;
        sumY4 += diffQ * diffQ * diffQ * diffQ;
98
    }
99

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

112
    void updateStatistic() {
113
114
115
116
        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;
117
    }
niklas.baumgarten's avatar
niklas.baumgarten committed
118

niklas.baumgarten's avatar
niklas.baumgarten committed
119
120
121
122
123
    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
124
125
    }

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

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

#endif //MLMC_MC_HPP