Commit 55db0b1d authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

Refactored MonteCarlo class

parent 3ad3ca7e
#include "MonteCarlo.h"
void MonteCarlo::Method() {
logger->StartMethod("Start MC Method", verbose);
if (!onlyFineLevel)
logger->InnerMsg("l: " + to_string(l) + " dM: " + to_string(dM), verbose);
else
logger->InnerMsg("l: " + to_string(l) + " dM: " + to_string(dM)
+ " onlyFineLevel: True", verbose);
method();
updateAvg();
updateStatistic();
logger->InnerMsg("|E[Y]|: " + to_string(avgY) + " V[Y]: " + to_string(varY),
verbose);
logger->EndMethod(verbose);
}
void MonteCarlo::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 MonteCarlo::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 MonteCarlo::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 MonteCarlo::mkSampleDirName(int m, bool fineLevel) {
if (fineLevel)
return buildName(m, "sample");
else
return buildName(m, "sample", "coarse");
}
void MonteCarlo::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 MonteCarlo::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;
}
......@@ -10,6 +10,8 @@
class MonteCarlo {
private:
public:
int plotting = 0;
Plot *coarsePlot = nullptr;
......@@ -61,85 +63,31 @@ public:
delete finePlot;
}
void Method() {
logger->StartMethod("Start MC Method", verbose);
if (!onlyFineLevel)
logger->InnerMsg("l: " + to_string(l) + " dM: " + to_string(dM), verbose);
else
logger->InnerMsg("l: " + to_string(l) + " dM: " + to_string(dM)
+ " onlyFineLevel: True", verbose);
method();
updateAvg();
updateStatistic();
logger->InnerMsg("|E[Y]|: " + to_string(avgY) + " V[Y]: " + to_string(varY),
verbose);
logger->EndMethod(verbose);
}
void Method();
protected:
virtual void initialize() = 0;
virtual void method() = 0;
virtual void solvePDE(int l, bool fineLevel,
double &valueQ, double &cost,
virtual void solvePDE(int l,
int m,
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 updateSums(double fineCost, double diffQ, double fineQ);
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 updateAvg();
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;
}
void updateStatistic();
string mkSampleDirName(int m, bool fineLevel) {
if (fineLevel)
return buildName(m, "sample");
else
return buildName(m, "sample", "coarse");
}
string mkSampleDirName(int m, bool fineLevel);
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);
}
static void mkSampleDir(int m, const string &sampleDirName);
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;
}
string buildName(int m, const string &name, const string &suffix = "");
};
#endif //MLMC_MC_HPP
\ No newline at end of file
......@@ -32,7 +32,7 @@ void MonteCarloElliptic::method() {
assemble->plot = finePlot;
assemble->problem->LoadNewSample(&finePerm);
solvePDE(l, true, fineQ, fineCost, fineSolution);
solvePDE(l, 0, true, fineQ, fineCost, fineSolution);
if (plotting) assemble->PlotResults(fineSolution);
......@@ -48,7 +48,7 @@ void MonteCarloElliptic::method() {
assemble->plot = coarsePlot;
assemble->problem->LoadNewSample(&coarsePerm);
solvePDE(l, false, coarseQ, coarseCost, coarseSolution);
solvePDE(l, 0, false, coarseQ, coarseCost, coarseSolution);
if (plotting) assemble->PlotResults(coarseSolution);
......@@ -59,7 +59,11 @@ void MonteCarloElliptic::method() {
}
}
void MonteCarloElliptic::solvePDE(int l, bool fineLevel, double &valueQ, double &cost,
void MonteCarloElliptic::solvePDE(int l,
int m,
bool fineLevel,
double &valueQ,
double &cost,
Vector &solution) {
Preconditioner *pc = GetPC("SuperLU");
Solver solver(pc, "GMRES");
......
......@@ -11,8 +11,11 @@ protected:
void method() override;
void solvePDE(int l, bool fineLevel,
double &valueQ, double &cost,
void solvePDE(int l,
int m,
bool fineLevel,
double &valueQ,
double &cost,
Vector &solution) override;
public:
......@@ -29,9 +32,6 @@ public:
solMatrixGraphs(MatrixGraphs(*meshes, *assemble->disc)) {
}
~MonteCarloElliptic() {
}
};
#endif //MLMC_MONTECARLOELLIPTIC_H
......@@ -22,7 +22,7 @@ void MonteCarloTransport::method() {
assemble->plot = finePlot;
assemble->problem->LoadNewSample(&fineNormalFlux);
solvePDE(l, true, fineQ, fineCost, fineSolution);
solvePDE(l, 0, true, fineQ, fineCost, fineSolution);
if (!onlyFineLevel) {
......@@ -37,7 +37,7 @@ void MonteCarloTransport::method() {
assemble->plot = coarsePlot;
assemble->problem->LoadNewSample(&coarseNormalFlux);
solvePDE(l, false, coarseQ, coarseCost, coarseSolution);
solvePDE(l, 0, false, coarseQ, coarseCost, coarseSolution);
} else {
coarseQ = 0.0, coarseCost = 0.0;
......@@ -47,11 +47,14 @@ void MonteCarloTransport::method() {
}
void MonteCarloTransport::solvePDE(int l,
int m,
bool fineLevel,
double &valueQ,
double &cost,
Vector &solution) {
assemble->logger->StartMethod("Start Solving PDE", verbose);
logger->StartMethod("Start Solving PDE on ", verbose);
logger->InnerMsg("l: " + to_string(l) +
" m: " + to_string(m), verbose);
Preconditioner *pc = GetPC("PointBlockGaussSeidel"); //TODO RMV after usage
Solver solver(pc, "GMRES");
......@@ -69,7 +72,7 @@ void MonteCarloTransport::solvePDE(int l,
if (functional == "Outflow")
valueQ = assemble->InFlowOutFlowRate(solution).second;
mout << endl;
assemble->logger->EndMethod(verbose);
logger->EndMethod(verbose);
}
TimeSeries MonteCarloTransport::getTimeSeries(bool fineLevel) {
......
......@@ -13,8 +13,11 @@ protected:
void method() override;
void solvePDE(int l, bool fineLevel,
double &valueQ, double &cost,
void solvePDE(int l,
int m,
bool fineLevel,
double &valueQ,
double &cost,
Vector &solution) override;
private:
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment