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

Merge branch '46-epsilon-for-montecarlo' into 'feature'

Resolve "Epsilon for MonteCarlo"

Closes #46

See merge request !56
parents be4251e0 33346672
Pipeline #168574 passed with stages
in 23 minutes and 4 seconds
......@@ -216,6 +216,7 @@ struct EstimatorMap : public LevelMap<Estimator *> {
WithMeshesCreator(meshesCreator).
WithOnlyFine((level_index == 0)).
WithParallel(parallel).
WithEpsilon(0.0).
Create()});
for (auto level = initLevels.back(); level <= maxLevel; level++)
......@@ -227,6 +228,7 @@ struct EstimatorMap : public LevelMap<Estimator *> {
WithInitLevel(level).
WithOnlyFine(false).
WithInitSamples(0).
WithEpsilon(0.0).
Create()});
}
};
......
......@@ -4,23 +4,24 @@
void MonteCarlo::Method() {
std::string parallelStr = aggregate.parallel ? "Parallel" : "Seriell";
mout.StartBlock(parallelStr + " MC l=" + to_string(level));
vout(1) << "epsilon=" << to_string(epsilon) << endl;
if (epsilon == 0.0) {
vout(1) << "Start with: " << aggregate;
method();
aggregate.UpdateParallel();
vout(1) << "End with: " << aggregate;
} else {
vout(1) << "epsilon=" << to_string(epsilon) << endl;
vout(1) << "Start with: " << aggregate;
adaptiveMethod();
aggregate.UpdateParallel();
vout(1) << "End with: " << aggregate;
}
std::unique_ptr<Meshes> meshes = meshesCreator.
WithMeshName(pdeSolverCreator.GetMeshName()).
WithCommSplit(aggregate.commSplit).
WithPLevel(level - 1).
WithLevel(level).
CreateUnique();
meshes->PrintInfo();
std::unique_ptr<PDESolver> pdeSolver = pdeSolverCreator.
CreateUnique(*meshes);
method(*pdeSolver);
value = aggregate.mean.Q;
cost = aggregate.mean.C;
// errors = 0.0;
mout.EndBlock(verbose == 0);
}
......@@ -32,37 +33,39 @@ void MonteCarlo::updateSampleIds(SampleSolution &fSol, SampleSolution &cSol) {
cSol.id = coarseId;
}
void MonteCarlo::method() {
// WithPLevel(onlyFine ? level : level - 1)
std::unique_ptr<Meshes> meshes = meshesCreator.
WithMeshName(pdeSolverCreator.GetMeshName()).
WithCommSplit(aggregate.commSplit).
WithPLevel(level - 1).
WithLevel(level).
CreateUnique();
void MonteCarlo::updateSampleCounter() {
int optimalM = (int) (ceil(aggregate.sVar.Q / (epsilon * epsilon)));
aggregate.UpdateSampleCounter(optimalM - aggregate.ctr.M);
}
meshes->PrintInfo();
void MonteCarlo::method(PDESolver &pdeSolver) {
std::unique_ptr<PDESolver> pdeSolver = pdeSolverCreator.
CreateUnique(*meshes);
SampleSolution fineSolution(pdeSolver.GetDisc(), fineId);
SampleSolution coarseSolution(pdeSolver.GetDisc(), coarseId);
SampleSolution fineSolution(pdeSolver->GetDisc(), fineId);
SampleSolution coarseSolution(pdeSolver->GetDisc(), coarseId);
while (aggregate.ctr.dM != 0) {
while (aggregate.ctr.dMcomm != 0) {
updateSampleIds(fineSolution, coarseSolution);
pdeSolver->DrawSample(fineId);
pdeSolver->Run(fineSolution);
vout(1) << "Start with: " << aggregate;
if (!onlyFine) {
pdeSolver->DrawSample(coarseId);
pdeSolver->Run(coarseSolution);
while (aggregate.ctr.dMcomm != 0) {
updateSampleIds(fineSolution, coarseSolution);
pdeSolver.Run(fineSolution);
if (!onlyFine) pdeSolver.Run(coarseSolution);
aggregate.Update(fineSolution, coarseSolution);
}
aggregate.Update(fineSolution, coarseSolution);
}
}
aggregate.UpdateParallel();
errors.Estimate(aggregate);
void MonteCarlo::adaptiveMethod() {
vout(1) << "End with: " << aggregate;
if (epsilon != 0.0) {
if (errors.numeric < epsilon) {
if (errors.total > epsilon)
updateSampleCounter();
} else {
Exit("Numeric Error is too big too reach error bound")
}
}
}
}
......@@ -28,9 +28,9 @@ protected:
void updateSampleIds(SampleSolution &fSol, SampleSolution &cSol);
void method();
void method(PDESolver &pdeSolver);
void adaptiveMethod();
void updateSampleCounter();
public:
MonteCarlo(double epsilon, int level, int initSamples, bool onlyFine, bool parallel,
......
......@@ -31,7 +31,6 @@ void StochasticCollocation::method() {
// Todo implement coarse solution solving
pdeSolver->DrawSample(fineId);
pdeSolver->Run(fineSolution);
aggregate.Update(fineSolution, coarseSolution);
......
......@@ -15,14 +15,16 @@ struct Errors {
void Estimate(const WelfordAggregate &aggregate) {
numeric = EstimateNumeric(aggregate.mean);
stochastic = EstimateStochastic(aggregate.ctr, aggregate.sVar);
total = numeric + stochastic;
}
double EstimateNumeric(const Mean &mean) {
return mean.Y; // todo
if (mean.Q == mean.Y) return 0.0;
return mean.Y;
}
double EstimateStochastic(const SampleCounter &ctr, const SVar &sVar) {
return sVar.Q / double(ctr.M);
return sqrt(sVar.Q / double(ctr.M));
}
};
......
......@@ -20,7 +20,6 @@ void HybridFaceNormalFluxGenerator::createPDESolver() {
// Todo delete sample solutions
void HybridFaceNormalFluxGenerator::drawSample(const SampleID &id) {
pdeSolver->DrawSample(id);
solutionFaceValues = new SampleSolution(pdeSolver->GetDisc(), id);
solutionFaceFlux = new SampleSolution(pdeSolver->GetDisc(), id);
pdeSolver->Run(*solutionFaceValues);
......
......@@ -55,7 +55,7 @@ public:
IStochasticProblem *GetProblem() const override { return assemble->GetProblem(); }
void DrawSample(const SampleID &id) override { assemble->DrawSample(id); }
void drawSample(const SampleID &id) override { assemble->DrawSample(id); }
std::string Name() const override { return "EllipticPDESolver"; }
};
......
......@@ -3,6 +3,7 @@
void PDESolver::Run(SampleSolution &solution) {
mout.StartBlock("PDE Solver");
drawSample(solution.id);
vout(2) << solution.id << endl;
run(solution);
computeQ(solution);
......
......@@ -20,6 +20,8 @@ protected:
virtual void run(SampleSolution &solution) = 0;
virtual void drawSample(const SampleID &id) = 0;
virtual void computeQ(SampleSolution &solution) = 0;
virtual void computeCost(SampleSolution &solution) = 0;
......@@ -46,8 +48,6 @@ public:
virtual IStochasticProblem *GetProblem() const = 0;
virtual void DrawSample(const SampleID &id) = 0;
virtual std::string Name() const = 0;
};
......@@ -78,7 +78,7 @@ public:
IStochasticProblem *GetProblem() const override { return assemble->GetProblem(); }
void DrawSample(const SampleID &id) override { assemble->DrawSample(id); }
void drawSample(const SampleID &id) override { assemble->DrawSample(id); }
std::string Name() const override { return "DummyPDESolver"; }
};
......
......@@ -61,7 +61,7 @@ public:
IStochasticProblem *GetProblem() const override { return nullptr; }
void DrawSample(const SampleID &id) override {}
void drawSample(const SampleID &id) override {}
std::string Name() const override { return "ParabolicPDESolver"; }
};
......
......@@ -56,7 +56,7 @@ public:
IStochasticProblem *GetProblem() const override { return assemble->GetProblem(); }
void DrawSample(const SampleID &id) override { assemble->DrawSample(id); }
void drawSample(const SampleID &id) override { assemble->DrawSample(id); }
std::string Name() const override { return "TransportPDESolver"; }
};
......
......@@ -2,7 +2,7 @@
INSTANTIATE_TEST_SUITE_P(TestMonteCarlo, TestMonteCarloWithoutEpsilon, Values(
TestParams{"ScalarGeneratorProblem", "FunctionEvaluation", "DummyPDESolver"},
TestParams{"ScalarGeneratorProblem", "FunctionEvaluation", "DummyPDESolver", 1, true},
TestParams{"StochasticLaplace1D", "L2", "LagrangeElliptic", 5},
TestParams{"StochasticLaplace2D", "L2", "LagrangeElliptic"},
TestParams{"StochasticLaplace2DTest", "L2", "LagrangeElliptic"},
......@@ -26,23 +26,37 @@ TEST_P(TestMonteCarloWithoutEpsilon, TestSeriellAgainstParallel) {
EXPECT_NEAR(mcParallel->aggregate.sVar.Y, mcSeriell->aggregate.sVar.Y, SVarTol());
}
//INSTANTIATE_TEST_SUITE_P(TestMonteCarlo, TestMonteCarloWithEpsilon, Values(
// TestParams{"ScalarGeneratorProblem", "FunctionEvaluation", "DummyPDESolver"},
// TestParams{"StochasticLaplace2DTest", "L2", "LagrangeElliptic"},
// TestParams{"StochasticLaplace2DTest", "Outflow", "HybridElliptic"}
//));
//
//TEST_P(TestMonteCarloWithEpsilon, TestWithEpsilon) {
// mout << GetParam() << endl;
//
// mout.StartBlock("Monte Carlo parallel");
// mout << "Start" << endl;
// mcParallel->Method();
// mout.EndBlock();
// mout << endl;
//
// mcParallel->EstimatorResults();
//}
INSTANTIATE_TEST_SUITE_P(TestMonteCarlo, TestMonteCarloWithEpsilon, Values(
TestParams{"ScalarGeneratorProblem", "FunctionEvaluation", "DummyPDESolver", 1, true},
TestParams{"StochasticLaplace1D", "L2", "LagrangeElliptic", 5},
TestParams{"StochasticLaplace2D", "L2", "LagrangeElliptic"},
TestParams{"StochasticLaplace2DTest", "L2", "LagrangeElliptic"},
TestParams{"StochasticLaplace2DTest", "Outflow", "HybridElliptic"}
));
TEST_P(TestMonteCarloWithEpsilon, TestSeriell) {
mout << GetParam() << endl;
mcSeriell->Method();
mout << endl;
mcSeriell->EstimatorResults();
EXPECT_LE(mcSeriell->StochasticError(), mcSeriell->Epsilon());
EXPECT_LE(mcSeriell->NumericError(), mcSeriell->Epsilon());
EXPECT_LE(mcSeriell->TotalError(), mcSeriell->Epsilon());
}
TEST_P(TestMonteCarloWithEpsilon, TestParallel) {
mout << GetParam() << endl;
mcParallel->Method();
mout << endl;
mcParallel->EstimatorResults();
EXPECT_LE(mcParallel->StochasticError(), mcParallel->Epsilon());
EXPECT_LE(mcParallel->NumericError(), mcParallel->Epsilon());
EXPECT_LE(mcParallel->TotalError(), mcParallel->Epsilon());
}
int main(int argc, char **argv) {
return MppTest(
......
......@@ -17,7 +17,7 @@ struct TestParams {
int level = 3;
bool onlyFine = true;
bool onlyFine = false;
};
Logging &operator<<(Logging &s, const TestParams &testParams) {
......@@ -40,12 +40,17 @@ protected:
std::unique_ptr<Estimator> mcParallel;
double MeanTol() const { return sqrt(1.0 / samples); }
double MeanTol() const {
return 2 * max(mcSeriell->StochasticError(), mcParallel->StochasticError());
}
double SVarTol() const { return sqrt(10.0 / samples); }
double SVarTol() const {
return 2 * max(mcSeriell->StochasticError(), mcParallel->StochasticError());
}
TestMonteCarlo(double epsilon, int samples) :
samples(samples), epsilon(epsilon),
meshesCreator(MeshesCreator().
WithDistribute("RCB").
WithoutOverlap()),
......@@ -87,7 +92,7 @@ class TestMonteCarloWithEpsilon : public TestMonteCarlo {
// 0.01
// 0.003
public:
TestMonteCarloWithEpsilon() : TestMonteCarlo(0.03, 10) {}
TestMonteCarloWithEpsilon() : TestMonteCarlo(0.03, 20) {}
};
#endif //TESTMONTECARLO_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