Commit 9f617e8c authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

Adapted MultilevelMonteCarlo

parent e3ef4f7b
......@@ -6,20 +6,15 @@ using namespace std;
void MultilevelMonteCarlo::initializeMapMonteCarlo() {
clearMapMonteCarlo();
for (unsigned long i = 0; i < initLevels.size(); i++) {
bool onlyFineLevel = (i == 0) || (mcOnly);
bool onlyFine = (i == 0) || (mcOnly);
int l = initLevels[i], M = initSampleAmount[i];
mapMonteCarlo[l] = createMonteCarlo(l, M, onlyFineLevel);
mapMonteCarlo[l]->Initialize();
mapMonteCarlo[l] = createMonteCarlo(l, M, onlyFine);
}
}
MonteCarlo *MultilevelMonteCarlo::createMonteCarlo(int l, int dM, bool onlyFineLevel) {
auto ellipticAssemble = dynamic_cast<LagrangeEllipticAssemble *>(assemble);
auto transportAssemble = dynamic_cast<DGTransportAssemble *>(assemble);
if (ellipticAssemble != nullptr)
return new MonteCarloElliptic(l, dM, onlyFineLevel, meshes, ellipticAssemble);
if (transportAssemble != nullptr)
return new MonteCarloTransport(l, dM, onlyFineLevel, meshes, transportAssemble);
MonteCarlo *MultilevelMonteCarlo::createMonteCarlo(int l, int dM, bool onlyFine) {
return new MonteCarlo(l, dM, onlyFine, *meshes, assemble, new EllipticPDESolver());
// return new MonteCarlo(l, onlyFine, *meshes, assemble, new TransportPDESolver());
Exit("Wrong Assembly Class")
}
......@@ -46,7 +41,7 @@ void MultilevelMonteCarlo::Method(const double eps) {
bool converged = false;
while (!converged) {
for (auto &mc : mapMonteCarlo)
if (mc.second->dM > 0)
if (mc.second->ctr.dM > 0)
mc.second->Method();
estimateExponents();
updateSampleAmount(eps);
......@@ -73,9 +68,9 @@ void MultilevelMonteCarlo::estimateExponents(bool excludeBaseLevel) {
for (; iter != mapMonteCarlo.end(); iter++) {
levelList.push_back((double) iter->first);
avgsY.push_back(-log2(iter->second->avgY));
varsY.push_back(-log2(iter->second->varY));
costs.push_back(log2(iter->second->avgCost));
avgsY.push_back(-log2(iter->second->avgs.Y));
varsY.push_back(-log2(iter->second->vars.Y));
costs.push_back(log2(iter->second->avgs.Cost));
}
pair<double, double> res = linear_fit(levelList, avgsY);
......@@ -92,23 +87,23 @@ void MultilevelMonteCarlo::updateSampleAmount(const double &eps) {
int optimalM;
double factor = 0.0;
for (auto &mc : mapMonteCarlo)
factor += sqrt(mc.second->varY * mc.second->avgCost);
factor += sqrt(mc.second->vars.Y * mc.second->avgs.Cost);
for (auto &mc : mapMonteCarlo) {
optimalM = (int) (ceil(2 * pow(eps, -2) * factor *
sqrt(mc.second->varY / mc.second->avgCost)));
sqrt(mc.second->vars.Y / mc.second->avgs.Cost)));
if (optimalM == 1) optimalM++; // Hack
mapMonteCarlo[mc.first]->dM = optimalM - mapMonteCarlo[mc.first]->M;
mapMonteCarlo[mc.first]->ctr.dM = optimalM - mapMonteCarlo[mc.first]->ctr.M;
}
vout(1) << "dM=[";
for (auto &mc:mapMonteCarlo) vout(1) << " " << mc.second->dM;
for (auto &mc:mapMonteCarlo) vout(1) << " " << mc.second->ctr.dM;
vout(1) << "]" << endl;
}
bool MultilevelMonteCarlo::checkForNewLevel() {
bool newLevel = true;
for (auto &mc : mapMonteCarlo) {
if (mc.second->dM > 0) {
if (mc.second->ctr.dM > 0) {
newLevel = false;
break;
}
......@@ -121,7 +116,7 @@ void MultilevelMonteCarlo::estimateNumericalError() {
unsigned long exponent = mapMonteCarlo.size() - 2;
for (auto &mc : mapMonteCarlo) {
if (mc.first == mapMonteCarlo.begin()->first) continue;
numErr = max(mc.second->avgY / (pow(2.0, alpha) - 1) / (pow(2.0, exponent)),
numErr = max(mc.second->avgs.Y / (pow(2.0, alpha) - 1) / (pow(2.0, exponent)),
numErr);
exponent--;
}
......@@ -131,7 +126,7 @@ void MultilevelMonteCarlo::estimateNumericalError() {
void MultilevelMonteCarlo::estimateStatisticalError() {
statErr = 0.0;
for (auto &mc : mapMonteCarlo) {
statErr += mc.second->varY / mc.second->M;
statErr += mc.second->vars.Y / mc.second->ctr.M;
}
statErr = sqrt(statErr);
vout(1) << "statErr=" << to_string(numErr) << endl;
......@@ -144,25 +139,25 @@ void MultilevelMonteCarlo::appendLevel(double eps) {
} else {
vout(1) << "L_new=" << to_string(L) << endl;
mapMonteCarlo[L] = createMonteCarlo(L, 0, false);
mapMonteCarlo[L]->varY = mapMonteCarlo[L - 1]->varY / pow(2.0, beta);
mapMonteCarlo[L]->avgCost = mapMonteCarlo[L - 1]->avgCost * pow(2.0, gamma);
mapMonteCarlo[L]->vars.Y = mapMonteCarlo[L - 1]->vars.Y / pow(2.0, beta);
mapMonteCarlo[L]->avgs.Cost = mapMonteCarlo[L - 1]->avgs.Cost * pow(2.0, gamma);
updateSampleAmount(eps);
}
}
void MultilevelMonteCarlo::wrapUpResults() {
for (auto &mc : mapMonteCarlo) {
value += mc.second->avgY;
cost += mc.second->sumCost;
value += mc.second->avgs.Y;
cost += mc.second->sums.Cost;
levels.push_back(mc.first);
samples.push_back(mc.second->M);
avgY.push_back(mc.second->avgY);
avgQ.push_back(mc.second->avgQ);
varY.push_back(mc.second->varY);
varQ.push_back(mc.second->varQ);
kurtosis.push_back(mc.second->kurtosis);
avgCost.push_back(mc.second->avgCost);
samples.push_back(mc.second->ctr.M);
avgY.push_back(mc.second->avgs.Y);
avgQ.push_back(mc.second->avgs.Q);
varY.push_back(mc.second->vars.Y);
varQ.push_back(mc.second->vars.Q);
kurtosis.push_back(mc.second->kurtosis.Y);
avgCost.push_back(mc.second->avgs.Cost);
}
}
......
#ifndef MULTILEVELMONTECARLO_HPP
#define MULTILEVELMONTECARLO_HPP
#include "montecarlo/MonteCarloElliptic.hpp"
#include "montecarlo/MonteCarloTransport.hpp"
#include "montecarlo/MonteCarlo.hpp"
#include "main/Utils.hpp"
typedef std::map<int, MonteCarlo *> MapMonteCarlo;
class MultilevelMonteCarlo {
private:
int verbose = 1;
public:
int verbose = -1;
int plotting = -1;
int maxLevel, pLevel;
vector<int> &initLevels, &initSampleAmount;
MapMonteCarlo mapMonteCarlo;
Meshes *meshes;
IAssemble *assemble;
IStochasticAssemble *assemble;
vector<int> levels = {};
vector<int> samples = {};
......@@ -41,13 +41,12 @@ public:
MultilevelMonteCarlo(vector<int> &initLevels,
vector<int> &initSampleAmount,
Meshes *meshes,
IAssemble *assemble) :
IStochasticAssemble *assemble) :
initLevels(initLevels),
initSampleAmount(initSampleAmount),
meshes(meshes),
assemble(assemble) {
config.get("MLMCPlotting", plotting);
config.get("mcOnly", mcOnly);
config.get("MLMCVerbose", verbose);
config.get("epsilon", epsilon);
......@@ -79,7 +78,7 @@ private:
void initializeValues();
MonteCarlo *createMonteCarlo(int l, int dM, bool onlyFineLevel);
MonteCarlo *createMonteCarlo(int l, int dM, bool onlyFine);
void clearMapMonteCarlo();
......
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