Commit e3ef4f7b authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

refactored MonteCarlo

parent 5b1d3bfa
#include "MonteCarlo.hpp"
// Todo use here
//void MultilevelMonteCarlo::checkKurtosis(int l, MonteCarlo *mc) {
// if (mc->kurtosis > 100.0) {
// vout(1) << endl << "WARNING: kurtosis on level " << setw(2) << l
// << " is " << setw(6) << mc->kurtosis << endl
// << "indicates MLMC correction is dominated by a few rare paths!"
// << endl;
// }
//}
void MonteCarlo::Method() {
mout.StartBlock("MC Method");
vout(1) << "l=" << l << " dM=" << dM << " fine=" << onlyFineLevel << endl;
vout(1) << "l=" << level << " dM=" << ctr.dM << " fine=" << onlyFine << endl;
method();
updateAvg();
updateStatistic();
vout(1) << "|E[Y]|=" << avgY << " V[Y]=" << varY << endl;
mout.EndBlock(verbose == 0);
}
void MonteCarlo::SolvePDE(SampleID id, double &valueQ, double &cost, Vector &solution) {
mout.StartBlock("Solve PDE");
vout(1) << id.Str() << endl;
solvePDE(id, valueQ, cost, solution);
vout(1) << "Q=" << valueQ << " cost=" << cost << endl;
avgs.Update(sums, ctr.M);
vars.Update(avgs);
kurtosis.Update(avgs, vars);
vout(1) << "|E[Y]|=" << avgs.Y << " V[Y]=" << vars.Y << endl;
mout.EndBlock(verbose == 0);
}
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;
void MonteCarlo::method() {
for (int m = ctr.M; m < ctr.M + ctr.dM; m++) {
fineId.number = m;
fineSolution.id = fineId;
plotter->SetDirectory(fineId);
assemble->DrawSample(fineId);
pdeSolver->Run(assemble, fineSolution);
if (!onlyFine) {
coarseId.number = m;
coarseSolution.id = coarseId;
plotter->SetDirectory(coarseId);
assemble->DrawSample(coarseId);
pdeSolver->Run(assemble, coarseSolution);
} else {
coarseSolution.Init();
}
sums.Update(fineSolution.Cost,
fineSolution.Q - coarseSolution.Q,
fineSolution.Q);
}
}
......@@ -2,61 +2,65 @@
#define MLMC_MC_HPP
#include "main/MultilevelPlotter.hpp"
#include "StatisticMeasures.hpp"
#include "PDESolver.hpp"
#include "utility/SampleID.hpp"
#include "assemble/IStochasticAssemble.hpp"
class MonteCarlo {
public:
protected:
int plotting = 0;
int verbose = 0;
int l = 0;
int pLevel = 0;
bool onlyFineLevel = true;
string functional = "L2";
Meshes *meshes = 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;
int verbose = 0;
double fineQ = 0.0, coarseQ = 0.0;
double fineCost = 0.0, coarseCost = 0.0;
virtual void method();
MonteCarlo(int l, int dM, bool onlyFineLevel, Meshes *meshes) :
l(l), dM(dM), onlyFineLevel(onlyFineLevel),
meshes(meshes), pLevel(meshes->pLevel()) {
public:
int level = 0;
int pLevel = 0;
bool onlyFine = true;
SampleCounter ctr;
Sums sums;
Averages avgs;
Variances vars;
Kurtosis kurtosis;
IStochasticAssemble *assemble;
PDESolver *pdeSolver;
MatrixGraphs solMGraphs;
SampleID coarseId;
SampleID fineId;
SampleSolution coarseSolution;
SampleSolution fineSolution;
MonteCarlo(int level,
int dM,
bool onlyFine,
Meshes &meshes,
IStochasticAssemble *assemble,
PDESolver *pdeSolver) :
level(level), onlyFine(onlyFine),
assemble(assemble), pdeSolver(pdeSolver),
solMGraphs(MatrixGraphs(meshes, *assemble->GetDiscretization())),
coarseSolution(SampleSolution(solMGraphs, level - 1)),
fineSolution(SampleSolution(solMGraphs, level)) {
config.get("MCPlotting", plotting);
config.get("MCVerbose", verbose);
config.get("Functional", functional);
ctr.dM = dM;
fineId.level = level;
fineId.coarse = false;
coarseId.level = level - 1;
coarseId.coarse = true;
}
virtual ~MonteCarlo() {};
void Method();
void SolvePDE(SampleID id, double &valueQ, double &cost, Vector &solution);
virtual void Initialize() = 0;
protected:
virtual void method() = 0;
virtual void solvePDE(SampleID id, double &valueQ, double &cost, Vector &solution) = 0;
void updateSums(double fineCost, double diffQ, double fineQ);
void updateAvg();
void updateStatistic();
};
#endif //MLMC_MC_HPP
\ No newline at end of file
#include "MonteCarloElliptic.hpp"
using namespace std;
void MonteCarloElliptic::Initialize() {
fineSolution = 0.0, coarseSolution = 0.0;
fineQ = 0.0, coarseQ = 0.0, fineCost = 0.0, coarseCost = 0.0;
}
void MonteCarloElliptic::method() {
Initialize();
SampleID id;
for (int m = M; m < M + dM; m++) {
id = SampleID(l, m, false);
plotter->SetDirectory(id);
assemble->problem->DrawSample(id);
SolvePDE(id, fineQ, fineCost, fineSolution);
if (!onlyFineLevel) {
id = SampleID(l, m, true);
plotter->SetDirectory(id);
assemble->problem->DrawSample(id);
SolvePDE(SampleID(), coarseQ, coarseCost, coarseSolution);
} else {
coarseQ = 0.0, coarseCost = 0.0, coarseSolution = 0.0;
}
updateSums(fineCost, fineQ - coarseQ, fineQ);
if (plotting) {
plotter->PlotVector("u", fineSolution,
1, l, "VertexData");
plotter->PlotVector("u", coarseSolution,
1, l - 1, "VertexData");
}
}
}
void MonteCarloElliptic::solvePDE(SampleID id,
double &valueQ,
double &cost,
Vector &solution) {
newton->operator()(*assemble, solution);
cost = solution.pSize();
if (functional == "L2")
valueQ = assemble->L2(solution);
if (functional == "Energy")
valueQ = assemble->Energy(solution);
if (functional == "Outflow")
valueQ = assemble->InflowOutflow(solution).second;
if (functional == "Goal")
valueQ = assemble->GoalFunctional(solution);
}
\ No newline at end of file
#ifndef MLMC_MONTECARLOELLIPTIC_H
#define MLMC_MONTECARLOELLIPTIC_H
#include "MonteCarlo.hpp"
#include "assemble/LagrangeEllipticAssemble.hpp"
#include "solver/Solver.h"
#include "solver/Newton.h"
#include "dof/BasicDoFs.hpp"
class MonteCarloElliptic : public MonteCarlo {
protected:
void Initialize() override;
void method() override;
void solvePDE(SampleID id, double &valueQ, double &cost, Vector &solution) override;
public:
LagrangeEllipticAssemble *assemble;
MatrixGraphs solMatrixGraphs;
Vector fineSolution;
Vector coarseSolution;
Preconditioner *pc;
Solver *solver;
Newton *newton;
MonteCarloElliptic(int l, int dM, bool onlyFineLevel,
Meshes *meshes, LagrangeEllipticAssemble *assemble) :
MonteCarlo(l, dM, onlyFineLevel, meshes),
assemble(assemble),
solMatrixGraphs(MatrixGraphs(*meshes, *assemble->disc)),
fineSolution(Vector(solMatrixGraphs[l - pLevel])),
coarseSolution(Vector(solMatrixGraphs[l - pLevel - 1])) {
solver = new Solver(GetPC("SuperLU"), "GMRES");
newton = new Newton(*solver);
}
virtual ~MonteCarloElliptic() {
delete solver;
delete newton;
}
};
#endif //MLMC_MONTECARLOELLIPTIC_H
#include "MonteCarloTransport.hpp"
using namespace std;
void MonteCarloTransport::Initialize() {
fineInput = 0.0, coarseInput = 0.0;
fineSolution = 0.0, coarseSolution = 0.0;
fineQ = 0.0, coarseQ = 0.0, fineCost = 0.0, coarseCost = 0.0;
}
void MonteCarloTransport::method() {
Initialize();
SampleID id;
for (int m = M; m < M + dM; m++) {
id = SampleID{l, m, false};
plotter->SetDirectory(id);
assemble->problem->DrawSample(id);
SolvePDE(id, fineQ, fineCost, fineSolution);
if (onlyFineLevel)
coarseQ = 0.0, coarseCost = 0.0;
else {
id = SampleID{l, m, true};
plotter->SetDirectory(id);
assemble->problem->DrawSample(id);
SolvePDE(SampleID(),
coarseQ, coarseCost, coarseSolution);
}
updateSums(fineCost, fineQ - coarseQ, fineQ);
}
}
void MonteCarloTransport::solvePDE(SampleID id,
double &valueQ,
double &cost,
Vector &solution) {
Preconditioner *pc = GetPC("PointBlockGaussSeidel");
Solver solver(pc, "GMRES");
TimeSeries timeSeries = getTimeSeries(!id.coarse);
assemble->resetTAssemble();
timeInt->operator()(assemble, timeSeries, solution);
cost = solution.pSize() * assemble->getStep();
if (functional == "Mass")
valueQ = assemble->Mass(solution);
if (functional == "Energy")
valueQ = assemble->Energy(solution);
if (functional == "Outflow")
valueQ = assemble->InFlowOutFlowRate(solution).second;
}
TimeSeries MonteCarloTransport::getTimeSeries(bool fineLevel) {
int scaling = 1;
config.get("scaling", scaling);
if (fineLevel)
return TimeSeries(startTime, endTime,
scaling * (int) (pow(2, l) * (endTime - startTime)));
else
return TimeSeries(startTime, endTime,
scaling * (int) (pow(2, l - 1) * (endTime - startTime)));
}
#ifndef MLMC_MONTECARLOTRANSPORT_HPP
#define MLMC_MONTECARLOTRANSPORT_HPP
#include "MonteCarlo.hpp"
#include "assemble/DGTransportAssemble.hpp"
#include "timestepping/TimeIntegrator.hpp"
#include "matrixgraph/CellMatrixGraph.hpp"
#include "TimeSeries.h"
#include "dof/BasicDoFs.hpp"
class MonteCarloTransport : public MonteCarlo {
protected:
void Initialize() override;
void method() override;
void solvePDE(SampleID id, double &valueQ, double &cost, Vector &solution) override;
private:
double startTime = 0.0;
double endTime = 0.0;
TimeSeries getTimeSeries(bool fineLevel);
public:
DGTransportAssemble *assemble;
Preconditioner *pc;
Solver *solver;
TimeIntegrator *timeInt;
MatrixGraphs cellMatrixGraphs;
MatrixGraphs faceMatrixGraphs;
CellMatrixGraphs solMatrixGraphs;
Vector fineInput;
Vector coarseInput;
Vector fineSolution;
Vector coarseSolution;
MonteCarloTransport(int l, int dM, bool baseLevel,
Meshes *meshes, DGTransportAssemble *assemble) :
MonteCarlo(l, dM, baseLevel, meshes),
assemble(assemble),
cellMatrixGraphs(MatrixGraphs(*meshes, dof(new CellDoF(3)))),
faceMatrixGraphs(MatrixGraphs(*meshes, dof(new FaceDoF(3)))),
solMatrixGraphs(*meshes, *assemble->disc),
fineInput(Vector(faceMatrixGraphs[l - pLevel])),
coarseInput(Vector(faceMatrixGraphs[l - pLevel - 1])),
fineSolution(Vector(solMatrixGraphs[l - pLevel])),
coarseSolution(Vector(solMatrixGraphs[l - pLevel - 1])),
pc(GetPC("PointBlockGaussSeidel")),
solver(new Solver(pc, "GMRES")),
timeInt(new TimeIntegrator(*solver)) {
config.get("startTime", startTime);
config.get("endTime", endTime);
}
};
#endif //MLMC_MONTECARLOTRANSPORT_HPP
Markdown is supported
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