MonteCarlo.h 2.92 KB
Newer Older
1
2
3
4
5
6
7
8
#ifndef MLMC_MC_HPP
#define MLMC_MC_HPP

#include "Elliptic.h"

/*
 * MonteCarlo
 *
9
 * l:          Level, should be set as the index in the maps
10
11
12
13
14
15
 * baseLevel:  Bool to indicate if on base level -> Only fine grid calculations
 * dM:         Samples to go
 * M:          Samples so far
 * functional: Type of functional used
 * sum_cost:   Summed cost
 * avg_cost:   Average cost
16
 *
17
18
19
20
 * Qf:         Functional on fine level
 * Qc:         Functional on coarse level
 * Q = Qf
 * Y = Qf - Qc
21
 *
22
23
24
25
 * sum_Y, sum_Q:       Sum of differences, functionals
 * avg_Y, avg_Q:       Average of differences, functionals
 * sum_Y2, sum_Q2:     Sum of squared differences, functionals
 * E(Y^2), E(Q^2):     Average of squared differences, functionals
26
27
28
 *
 * var_Y = E(Y^2) - (E(Y))^2, V(Q)
 *       = E(Q^2) - (E(Q))^2
29
30
31
32
33
34
35
36
37
38
39
40
 *
 * MCPlotting = 0     for no plots at all
 * MCPlotting = 1     for average plots
 * MCPlotting = 2     for sample and average plots
 * MCPlotting = 3     for sample, average and interpolated fields
 *
 */

class MonteCarlo {
public:
    int verbose = 0, plotting = 0;

41
    int l = 0, plevel = 0;
42
43
44
    bool baseLevel = true;
    string functional = "L2";

45
    Meshes &meshes;
46
47
48
    StochasticFields &stochFields;
    MatrixGraphs &graphs;

49
50
51
    Vector u_sum, u_avg, u_f, u_c;

    Plot plot_c, plot_f;
52
53
54
55
56
57
58
59
60
61
62
63

    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;

64
65
66
67
68
69
70
71
72
73
    MonteCarlo(int l, int dM, bool baseLevel, Meshes &meshes,
               StochasticFields &stochFields, MatrixGraphs &graphs) :

            l(l), dM(dM), baseLevel(baseLevel), meshes(meshes), plevel(meshes.pLevel()),
            stochFields(stochFields), graphs(graphs),
            u_sum(graphs[l - plevel]), u_avg(u_sum),
            u_f(u_sum), u_c(graphs[l - plevel - 1]),
            plot_c(meshes[l - 1], meshes.dim(), 2),
            plot_f(meshes[l], meshes.dim(), 2) {

74
75
76
        ReadConfig(Settings, "MCVerbose", verbose);
        ReadConfig(Settings, "MCPlotting", plotting);
        ReadConfig(Settings, "functional", functional);
77

78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
    }

    virtual ~MonteCarlo() = default;;

    void updateSum(const double *sums);

    void updateAvg();

    void updateStatistic();

    void upscaleU();

    void updateUSum();

    void updateUAvg();

    virtual void method() = 0;

96
    virtual void solvePDE(double &Q, double &cost, Vector &u) {};
97
98
};

99
100
101
102
MonteCarlo *
getMonteCarlo(int l, int dM, bool baseLevel, Meshes &meshes,
              StochasticFields &stochFields, Assemble *assemble,
              MatrixGraphs &graphs);
103
104
105
106
107
108

void estimateExponents(const map<int, MonteCarlo *> &df_mc,
                       double &alpha, double &beta, double &gamma,
                       bool excludeBaseLevel = true);

#endif //MLMC_MC_HPP