MonteCarlo.h

#ifndef MLMC_MC_HPP
#define MLMC_MC_HPP

#include "m++.h"
#include "stochastics/StochasticField.h"
#include "MonteCarloLogger.h"
#include <sys/stat.h>

#include <utility>


class MonteCarlo {
public:
    int plotting = 0;
    Plot *coarsePlot = nullptr;
    Plot *finePlot = nullptr;
    MonteCarloLogger *logger = nullptr;

    int l = 0;
    int pLevel = 0;
    bool onlyFineLevel = true;
    string functional = "L2";

    Meshes *meshes = nullptr;
    StochasticField *stochasticField = nullptr;

    int M = 0, dM = 0;

    double sumCost = 0.0, avgCost = 0.0;
    double sumY1 = 0.0, avgY = 0.0;
    double sumY2 = 0.0, avgY2 = 0.0, varY = 0.0;
    double sumY3 = 0.0, avgY3 = 0.0;
    double sumY4 = 0.0, avgY4 = 0.0;
    double sumQ1 = 0.0, avgQ = 0.0;
    double sumQ2 = 0.0, avgQ2 = 0.0, varQ = 0.0;
    double kurtosis = 0.0;

    MonteCarlo(int l, int dM, bool onlyFineLevel,
               Meshes *meshes, StochasticField *stochFields) :
        l(l), dM(dM), onlyFineLevel(onlyFineLevel),
        meshes(meshes), pLevel(meshes->pLevel()),
        stochasticField(stochFields) {

        config.get("MCPlotting", plotting);
        config.get("functional", functional);

        coarsePlot = new Plot((*meshes)[l - 1], meshes->dim());
        finePlot = new Plot((*meshes)[l], meshes->dim());
        logger = new MonteCarloLogger();
    }

    ~MonteCarlo() {
        delete coarsePlot;
        delete finePlot;
        delete logger;
    }

    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();
    }

protected:
    virtual void initialize() = 0;

    virtual void method() = 0;

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

    void updateSums(double fineCost, double diffQ, double fineQ) {
        sumCost += fineCost;
        sumQ1 += fineQ;
        sumQ2 += fineQ * fineQ;
        sumY1 += diffQ;
        sumY2 += diffQ * diffQ;
        sumY3 += diffQ * diffQ * diffQ;
        sumY4 += diffQ * diffQ * diffQ * diffQ;
    }

    void updateAvg() {
        M += dM;
        dM = 0;
        avgCost = sumCost / M;
        avgY = abs(sumY1 / M);
        avgY2 = sumY2 / M;
        avgY3 = sumY3 / M;
        avgY4 = sumY4 / M;
        avgQ = abs(sumQ1 / M);
        avgQ2 = sumQ2 / M;
    }

    void updateStatistic() {
        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;
    }

    string mkSampleDirName(int m, bool fineLevel) {
        if (fineLevel)
            return buildName(m, "sample");
        else
            return buildName(m, "sample", "coarse");
    }

    static void mkSampleDir(int m, const string &sampleDirName) {
        string pathToData = "data/vtk/";
        char *sampleDirNameCharArr =
            const_cast<char *> ((pathToData.append(sampleDirName)).c_str());
        int status = mkdir(sampleDirNameCharArr, 0777);
    }

    string buildName(int m, const string &name, const string &suffix = "") {
        string returnName = name;
        if (!suffix.empty())
            returnName = returnName.append("_").append(suffix);
        if (l != -1)
            returnName = returnName.append("_").append(to_string(l));
        if (m != -1)
            returnName = returnName.append("_").append(to_string(m));
        return returnName;
    }
};

#endif //MLMC_MC_HPP