refactoring and fitting MC classes to new main structure

mlmc/src/MonteCarlo.C

@@ -32,7 +32,9 @@ void MonteCarlo::updateStatistic() {
 }
 
 void MonteCarlo::upscaleU() {
-
+//        TransferMatrix transferMatrix(TransferMatrixGraph(graphs[l - plevel],
+//                                                          graphs[l - plevel - 1]));
+//        assemble.AssembleTransfer(transferMatrix);
 }
 
 void MonteCarlo::updateUSum() {
@@ -45,9 +47,12 @@ void MonteCarlo::updateUAvg() {
 
 class MonteCarloElliptic : public MonteCarlo {
 public:
-    MonteCarloElliptic(int l, int dM, bool baseLevel, StochasticFields &stochFields,
-                       EllipticAssemble &assemble, MatrixGraphs &graphs) :
-            MonteCarlo(l, dM, baseLevel,stochFields, assemble, graphs) {}
+    EllipticAssemble &assemble;
+
+    MonteCarloElliptic(int l, int dM, bool baseLevel, Meshes &meshes,
+                       StochasticFields &stochFields, MatrixGraphs &graphs,
+                       EllipticAssemble &assemble) : assemble(assemble),
+            MonteCarlo(l, dM, baseLevel, meshes, stochFields, graphs) {}
 
     void method() override {
         Date Start;
@@ -56,47 +61,38 @@ public:
         vout(1) << " with " << dM << " samples" << endl;
         stochFields.setCurrentLevel(l);
 
-        Vector u_f(u_sum);
-        Vector u_c(graphs[l - 2]);
-//        TransferMatrix transferMatrix(TransferMatrixGraph(graphs[l - 2],
-//                                                          graphs[l - 1]));
-//        assemble.AssembleTransfer(transferMatrix);
-
         for (int m = M; m < M + dM; m++) {
             stochFields.generate();
 
             if (verbose >= 2) {
                 vout(2) << "    Generated field:" << endl;
                 pair<double, double> min_max_f, min_max_c;
-                min_max_f = stochFields.getMinMaxValues();
+                min_max_f = stochFields.getMinMaxValues(true);
                 vout(2) << "      min_f=" << min_max_f.first
                         << " max_f=" << min_max_f.second << endl;
                 if (!baseLevel) {
-                    stochFields.setCurrentLevel(l - 1);
-                    min_max_c = stochFields.getMinMaxValues();
+                    min_max_c = stochFields.getMinMaxValues(false);
                     vout(2) << "      min_c=" << min_max_c.first
                             << " max_c=" << min_max_c.second << endl;
-                    stochFields.setCurrentLevel(l);
                 }
             }
 
             double Qf, Qc;
             double costf, costc;
 
-            u_f = 0.0;
-            u_c = 0.0;
+            u_f = 0.0, u_c = 0.0;
 
             Date Start_solving;
             vout(2) << "    Solving PDE on level " << l;
-            solvePDE(Qf, costf, u_f, m, l);
+            solvePDE(Qf, costf, u_f);
             vout(2) << " took " << Start_solving - Date() << endl;
 
             if (!baseLevel) {
                 Date Start_solving_c;
                 vout(2) << "    Solving PDE on level " << l - 1;
-                graphs.SetLevel(l - 2);
-                solvePDE(Qc, costc, u_c, m, l - 1);
-                graphs.SetLevel(l - 1);
+                stochFields.setFineField(false);
+                solvePDE(Qc, costc, u_c);
+                stochFields.setFineField(true);
                 vout(2) << " took " << Start_solving_c - Date() << endl;
             } else {
                 Qc = 0.0;
@@ -111,20 +107,26 @@ public:
             updateSum(sums);
             updateUSum();
 
-            if (plotting >= 2) assemble.plotSolution(u_f, l);
-            if (plotting == 3) assemble.plotSolution(u_c, l);
+            if (plotting >= 2)
+                assemble.plotSolution(plot_f, u_f, meshes, l, m);
+            if (plotting == 3 && !baseLevel) {
+                stochFields.setFineField(false);
+                assemble.plotSolution(plot_c, u_c, meshes, l - 1, m, "ip");
+                stochFields.setFineField(true);
+            }
         }
         updateAvg();
         updateUAvg();
         updateStatistic();
 
-        if (plotting >= 1) assemble.plotSolution(u_avg, l);
+        if (plotting >= 1)
+//            assemble.plotSolution(plot_f, u_avg, meshes, l, -1, "avg");
 
         vout (1) << "  End after " << Date() - Start << ": |E(Y)|=" << avg_Y
                  << " V(Y)=" << var_Y << endl;
     }
 
-    void solvePDE(double &Q, double &cost, Vector &u, int m, int l) override {
+    void solvePDE(double &Q, double &cost, Vector &u) override {
 //        Preconditioner *preconditioner = GetPC(*graphs, assemble, "Multigrid");
 //        Solver solver(preconditioner, "GMRES");
         Solver solver(graphs, assemble);
@@ -134,20 +136,21 @@ public:
 
         if (functional == "L2") Q = assemble.L2(u);
         if (functional == "Energy") Q = assemble.Energy(u);
-//        if (functional == "Max") Q = assemble.Max(u);
-        if (functional == "Outflow") Q = assemble.FluxError(u);
+        if (functional == "Outflow") Q = assemble.getInflowOutflow(u).second;
     }
 };
 
-MonteCarlo *getMonteCarlo(int l, int dM, bool baseLevel, StochasticFields &stochFields,
-                          EllipticAssemble &assemble, MatrixGraphs &graphs,
-                          string MonteCarlo) {
-    if (MonteCarlo.empty())
-        ReadConfig(Settings, "MonteCarlo", MonteCarlo);
-    if (MonteCarlo == "MonteCarloElliptic")
-        return new MonteCarloElliptic(l, dM, baseLevel,
-                                      stochFields, assemble, graphs);
-    else Exit("This MonteCarlo is not implemented")
+MonteCarlo *getMonteCarlo(int l, int dM, bool baseLevel, Meshes &meshes,
+                          StochasticFields &stochFields, Assemble *assemble,
+                          MatrixGraphs &graphs) {
+
+    auto ellipticAssemble = dynamic_cast<EllipticAssemble *>(assemble);
+    auto transportAssemble = dynamic_cast<EllipticAssemble *>(assemble);
+    if (ellipticAssemble != nullptr)
+        return new MonteCarloElliptic(l, dM, baseLevel, meshes,
+                                      stochFields, graphs, (*ellipticAssemble));
+    if(transportAssemble != nullptr)
+        return nullptr;
     return nullptr;
 }
 

mlmc/src/MonteCarlo.h

@@ -6,23 +6,26 @@
 /*
  * MonteCarlo
  *
- * l:          Level, should be set as the index in the data frames
+ * l:          Level, should be set as the index in the maps
  * 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
+ *
  * Qf:         Functional on fine level
  * Qc:         Functional on coarse level
  * Q = Qf
  * Y = Qf - Qc
+ *
  * 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
- * var_Y = E(Y^2) - (E(Y))^2, V(Q) = E(Q^2) - (E(Q))^2:
- *                     Variance of differences, functionals
+ *
+ * var_Y = E(Y^2) - (E(Y))^2, V(Q)
+ *       = E(Q^2) - (E(Q))^2
  *
  * MCPlotting = 0     for no plots at all
  * MCPlotting = 1     for average plots
@@ -35,15 +38,17 @@ class MonteCarlo {
 public:
     int verbose = 0, plotting = 0;
 
-    int l = 0;
+    int l = 0, plevel = 0;
     bool baseLevel = true;
     string functional = "L2";
 
+    Meshes &meshes;
     StochasticFields &stochFields;
-    EllipticAssemble &assemble;
     MatrixGraphs &graphs;
 
-    Vector u_sum, u_avg;
+    Vector u_sum, u_avg, u_f, u_c;
+
+    Plot plot_c, plot_f;
 
     int M = 0, dM = 0;
 
@@ -56,14 +61,20 @@ public:
     double sum_Q2 = 0.0, avg_Q2 = 0.0, var_Q = 0.0;
     double kurtosis = 0.0;
 
-    MonteCarlo(int l, int dM, bool baseLevel, StochasticFields &stochFields,
-               EllipticAssemble &assemble, MatrixGraphs &graphs) :
-            l(l), dM(dM), baseLevel(baseLevel),
-            stochFields(stochFields), assemble(assemble),
-            graphs(graphs), u_sum(graphs[l - 1]), u_avg(graphs[l - 1]) {
+    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) {
+
         ReadConfig(Settings, "MCVerbose", verbose);
         ReadConfig(Settings, "MCPlotting", plotting);
         ReadConfig(Settings, "functional", functional);
+
     }
 
     virtual ~MonteCarlo() = default;;
@@ -82,12 +93,13 @@ public:
 
     virtual void method() = 0;
 
-    virtual void solvePDE(double &Q, double &cost, Vector &u, int m, int l) {};
+    virtual void solvePDE(double &Q, double &cost, Vector &u) {};
 };
 
-MonteCarlo *getMonteCarlo(int l, int dM, bool baseLevel, StochasticFields &stochFields,
-                          EllipticAssemble &assemble, MatrixGraphs &graphs,
-                          string MonteCarlo = "");
+MonteCarlo *
+getMonteCarlo(int l, int dM, bool baseLevel, Meshes &meshes,
+              StochasticFields &stochFields, Assemble *assemble,
+              MatrixGraphs &graphs);
 
 void estimateExponents(const map<int, MonteCarlo *> &df_mc,
                        double &alpha, double &beta, double &gamma,

mlmc/src/MultilevelMonteCarlo.C

@@ -2,16 +2,25 @@
 
 using namespace std;
 
+void MultilevelMonteCarlo::method() {
+    Date Start;
+    vout (1) << "Start multilevel Monte Carlo method only for initial values" << endl;
+
+    for (auto &mc : map_mc)
+        mc.second->method();
+
+    vout (1) << "End after " << Date() - Start << endl;
+}
+
 void MultilevelMonteCarlo::method(const double eps) {
     Date Start;
-    vout (1) << "Start multilevel Monte Carlo method for epsilon=" << eps << endl;
+    vout (1) << "Start multilevel Monte Carlo method for epsilon " << eps << endl;
 
     bool converged = false;
     while (!converged) {
-        for (auto &iter : map_mc)
-            if (iter.second->dM > 0) {
-                iter.second->method();
-            }
+        for (auto &mc : map_mc)
+            if (mc.second->dM > 0)
+                mc.second->method();
 
         estimateExponents(map_mc, alpha, beta, gamma);
         vout (2) << "  Estimated exponents: " << "alpha=" << alpha
@@ -20,8 +29,8 @@ void MultilevelMonteCarlo::method(const double eps) {
         updateNumSamples(eps);
 
         bool check_for_new_level = true;
-        for (auto &iter : map_mc) {
-            if (iter.second->dM > 0) {
+        for (auto &mc : map_mc) {
+            if (mc.second->dM > 0) {
                 check_for_new_level = false;
                 break;
             }
@@ -35,10 +44,9 @@ void MultilevelMonteCarlo::method(const double eps) {
                     Exit ("Maximum level has been reached.")
                 } else {
                     vout(2) << "  Append level " << L << endl;
-                    map_mc[L] = getMonteCarlo(L, 0, false,
+                    map_mc[L] = getMonteCarlo(L, 0, false, meshes,
                                               stochFields,
-                                              assemble, graphs,
-                                              "MonteCarloElliptic");
+                                              assemble, graphs);
                     map_mc[L]->var_Y = map_mc[L - 1]->var_Y / pow(2.0, beta);
                     map_mc[L]->avg_cost = map_mc[L - 1]->avg_cost * pow(2.0, gamma);
                     updateNumSamples(eps);
@@ -49,38 +57,117 @@ void MultilevelMonteCarlo::method(const double eps) {
         }
     }
 
-    for (auto &iter : map_mc) {
-        value += iter.second->avg_Y;
-        cost += iter.second->sum_cost;
-        levels.push_back(iter.first);
-        numsamples.push_back(iter.second->M);
+    for (auto &mc : map_mc) {
+        value += mc.second->avg_Y;
+        cost += mc.second->sum_cost;
+        levels.push_back(mc.first);
+        numsamples.push_back(mc.second->M);
     }
     vout (1) << "End after " << Date() - Start << endl;
 }
 
+void MultilevelMonteCarlo::method(const vector<double> &eps) {
+    Date Start;
+    vout (1) << "Start multilevel Monte Carlo method for several eps" << endl;
+
+    for (auto &e : eps)
+        method(e);
+
+    vout (1) << "End after " << Date() - Start << endl;
+}
+
 void MultilevelMonteCarlo::updateNumSamples(const double &eps) {
     int M_optimal;
     double factor = 0.0;
-    for (auto &iter : map_mc)
-        factor += sqrt(iter.second->var_Y * iter.second->avg_cost);
-    for (auto &iter : map_mc) {
+    for (auto &mc : map_mc)
+        factor += sqrt(mc.second->var_Y * mc.second->avg_cost);
+    for (auto &mc : map_mc) {
         M_optimal = (int) (ceil(2 * pow(eps, -2) * factor *
-                                sqrt(iter.second->var_Y / iter.second->avg_cost)));
-        map_mc[iter.first]->dM = M_optimal - map_mc[iter.first]->M;
+                                sqrt(mc.second->var_Y / mc.second->avg_cost)));
+        map_mc[mc.first]->dM = M_optimal - map_mc[mc.first]->M;
     }
     vout(2) << "  Missing samples:";
-    for (auto &iter:map_mc) vout(1) << " " << iter.second->dM;
+    for (auto &mc:map_mc) vout(1) << " " << mc.second->dM;
     vout(2) << endl;
 }
 
 double MultilevelMonteCarlo::estimateNumericalError() {
     double maxErr = 0.0;
     unsigned long exponent = map_mc.size() - 2;
-    for (auto &iter : map_mc) {
-        if (iter.first == map_mc.begin()->first) continue;
-        maxErr = max(iter.second->avg_Y / (pow(2.0, alpha) - 1) / (pow(2.0, exponent)),
+    for (auto &mc : map_mc) {
+        if (mc.first == map_mc.begin()->first) continue;
+        maxErr = max(mc.second->avg_Y / (pow(2.0, alpha) - 1) /
+                     (pow(2.0, exponent)),
                      maxErr);
         exponent--;
     }
     return maxErr;
-}
+}
\ No newline at end of file
+
+void MultilevelMonteCarlo::showMCTable() {
+    mout << endl
+         << " l     M    E[Qf-Qc]       E[Qf]    V[Qf-Qc]       V[Qf]    kurtosis        cost"
+         << endl
+         << "--------------------------------------------------------------------------------"
+         << endl;
+
+    for (auto &mc : map_mc) {
+        int l = mc.first;
+        mout << setw(2) << l << setw(6) << mc.second->M
+             << setw(12) << mc.second->avg_Y
+             << setw(12) << mc.second->avg_Q << setw(12) << mc.second->var_Y
+             << setw(12) << mc.second->var_Q << setw(12) << mc.second->kurtosis
+             << setw(12) << mc.second->avg_cost << endl;
+    }
+    mout
+            << "--------------------------------------------------------------------------------"
+            << endl;
+}
+
+void MultilevelMonteCarlo::showKurtosisWarning() {
+    for (auto &mc : map_mc) {
+        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;
+        }
+    }
+}
+
+void MultilevelMonteCarlo::showExponentResults() {
+    estimateExponents(map_mc, alpha, beta, gamma, true);
+    mout << endl << "alpha = " << setw(8) << alpha
+         << "  (exponent for the weak convergence of the FEM)"
+         << endl << " beta = " << setw(8) << beta
+         << "  (exponent for the decay of the variance)"
+         << endl << "gamma = " << setw(8) << gamma
+         << "  (exponent for the growth of the cost)"
+         << endl;
+}
+
+void MultilevelMonteCarlo::showResultsMLMCRun(double eps) {
+    string level_str, sample_str;
+
+    mout << endl
+         << "   eps       value   MLMC cost             l                                   M"
+         << endl
+         << "--------------------------------------------------------------------------------"
+         << endl;
+
+    mout << setw(6) << eps << setw(12) << value << setw(12) << cost;
+    for (unsigned long i = 0 ; i < levels.size(); i++) {
+        if (i != levels.size() - 1) level_str += to_string(levels[i]) + " ";
+        else level_str += to_string(levels[i]);
+    }
+
+    mout << setw(14) << level_str;
+    for (unsigned long i = 0; i < numsamples.size(); i++) {
+        if (i != numsamples.size() - 1) sample_str += to_string(numsamples[i]) + " ";
+        else sample_str += to_string(numsamples[i]);
+    }
+    mout << setw(36) << sample_str << endl;
+}
+
+
+

mlmc/src/MultilevelMonteCarlo.h

@@ -9,30 +9,32 @@ class MultilevelMonteCarlo {
 public:
     int verbose = 0, plotting = 0;
 
-    int Lmax = 9;
+    int Lmax;
     map<int, MonteCarlo *> map_mc;
 
+    Meshes &meshes;
     StochasticFields &stochFields;
-    EllipticAssemble &assemble;
+    Assemble *assemble;
     MatrixGraphs &graphs;
 
     vector<int> levels, numsamples;
     double value = 0.0, cost = 0.0;
     double alpha = 0.0, beta = 0.0, gamma = 0.0;
 
-    MultilevelMonteCarlo(vector<int> l_init, vector<int> M_init, int Lmax,
-                         StochasticFields &stochFields, EllipticAssemble &assemble,
+    MultilevelMonteCarlo(vector<int> l_init, vector<int> M_init, Meshes &meshes,
+                         StochasticFields &stochFields, Assemble *assemble,
                          MatrixGraphs &graphs) :
-            stochFields(stochFields), assemble(assemble), graphs(graphs) {
+            stochFields(stochFields), assemble(assemble), graphs(graphs),
+            meshes(meshes), Lmax(meshes.Level()) {
 
         ReadConfig(Settings, "MLMCVerbose", verbose);
         ReadConfig(Settings, "MLMCPlotting", plotting);
 
         for (unsigned long i = 0; i < l_init.size(); i++) {
             bool baseLevel = i == 0;
-            map_mc[l_init[i]] = getMonteCarlo(l_init[i], M_init[i], baseLevel,
-                                              stochFields, assemble, graphs,
-                                              "MonteCarloElliptic");
+            int l = l_init[i], M = M_init[i];
+            map_mc[l] = getMonteCarlo(l, M, baseLevel, meshes,
+