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

4
5
#include "m++.h"
#include "stochfields/StochasticFields.h"
6
#include "MonteCarloLogger.h"
7
8
9
10

class MonteCarlo {
public:
    int verbose = 0, plotting = 0;
11
    MonteCarloLogger *logger;
12

13
    int l = 0, plevel = 0;
14
15
16
    bool baseLevel = true;
    string functional = "L2";

17
18
19
    Meshes *meshes;
    StochasticFields *stochFields;
    MatrixGraphs *graphs;
20

niklas.baumgarten's avatar
niklas.baumgarten committed
21
    Vector u_sum, u_avg, u_sum_diff, u_avg_diff;
22

niklas.baumgarten's avatar
niklas.baumgarten committed
23
    Plot plot;
24
25
26
27
28
29
30
31
32
33
34
35

    int M = 0, dM = 0;

    double sum_cost = 0.0, avg_cost = 0.0;
    double sum_Y = 0.0, avg_Y = 0.0;
    double sum_Y2 = 0.0, avg_Y2 = 0.0, var_Y = 0.0;
    double sum_Y3 = 0.0, avg_Y3 = 0.0;
    double sum_Y4 = 0.0, avg_Y4 = 0.0;
    double sum_Q = 0.0, avg_Q = 0.0;
    double sum_Q2 = 0.0, avg_Q2 = 0.0, var_Q = 0.0;
    double kurtosis = 0.0;

36
37
    MonteCarlo(int l, int dM, bool baseLevel,
               Meshes *meshes, StochasticFields *stochFields, MatrixGraphs *graphs) :
38

niklas.baumgarten's avatar
niklas.baumgarten committed
39
            l(l), dM(dM), baseLevel(baseLevel),
40
            meshes(meshes), plevel(meshes->pLevel()),
41
            stochFields(stochFields), graphs(graphs),
42
            u_sum((*graphs)[l - plevel]), u_avg(u_sum),
niklas.baumgarten's avatar
niklas.baumgarten committed
43
            u_sum_diff(u_sum), u_avg_diff(u_sum),
44
            plot((*meshes)[l], meshes->dim(), 2) {
45

46
47
        ReadConfig(Settings, "MCPlotting", plotting);
        ReadConfig(Settings, "functional", functional);
48

49
50
        logger = new MonteCarloLogger();

niklas.baumgarten's avatar
niklas.baumgarten committed
51
        u_sum = 0.0, u_avg = 0.0, u_sum_diff = 0.0, u_avg_diff = 0.0;
52

niklas.baumgarten's avatar
niklas.baumgarten committed
53
    }
54

55
56
57
    ~MonteCarlo() {
        delete logger;
    }
58

59
60
61
62
63
64
65
66
67
    void updateSum(const double *sums) {
        sum_cost += sums[0];
        sum_Y += sums[1];
        sum_Y2 += sums[2];
        sum_Y3 += sums[3];
        sum_Y4 += sums[4];
        sum_Q += sums[5];
        sum_Q2 += sums[6];
    }
68

69
70
71
72
73
74
75
76
77
78
79
    void updateAvg() {
        M += dM;
        dM = 0;
        avg_cost = sum_cost / M;
        avg_Y = abs(sum_Y / M);
        avg_Y2 = sum_Y2 / M;
        avg_Y3 = sum_Y3 / M;
        avg_Y4 = sum_Y4 / M;
        avg_Q = abs(sum_Q / M);
        avg_Q2 = sum_Q2 / M;
    }
80

81
82
83
84
85
86
    void updateStatistic() {
        var_Y = fmax(avg_Y2 - avg_Y * avg_Y, 1e-10);
        var_Q = fmax(avg_Q2 - avg_Q * avg_Q, 1e-10);
        kurtosis = (avg_Y4 - 4.0 * avg_Y3 * avg_Y + 6.0 * avg_Y2 * avg_Y * avg_Y -
                    3.0 * avg_Y * avg_Y * avg_Y * avg_Y) / var_Y / var_Y;
    }
87

88
    virtual void upscale(const Vector &uc, Vector &uc_up) = 0;
89

90
91
92
93
    virtual void method() = 0;
};

#endif //MLMC_MC_HPP