Commit 1bbb666e authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

Merge branch '42-create-all-level-estimators-in-estimator-map' into 'feature'

Revert "trying onlyFine fix"

Closes #42

See merge request !54
parents f086e2fd 0705c9e6
Pipeline #162778 passed with stages
in 21 minutes and 29 seconds
......@@ -14,11 +14,11 @@ variables:
HOST: 'horeka'
DOCKER_TLS_CERTDIR: "/certs"
services:
- docker:dind
# ------------------------ Start of CI ------------------------
build-mluq:
stage: build
......@@ -42,6 +42,7 @@ build-mluq:
dependencies: [ ]
tags: [ docker ]
ctest-mluq:
stage: test
variables:
......@@ -55,7 +56,36 @@ ctest-mluq:
dependencies: [ "build-mluq" ]
tags: [ docker ]
mpitest-mluq:
mpitest-1-mluq:
stage: test
variables:
GIT_STRATEGY: none
only:
variables: [ $RUN_EXPERIMENTS == 'none' ]
image: ${REGISTRY}/${IMAGE_NAME_MLUQ}
script:
- cd /mpp/build/
- python3 mppyrun.py -n 1 --mpi_tests=1 --mute=1
dependencies: [ "build-mluq" ]
tags: [ docker ]
mpitest-2-mluq:
stage: test
variables:
GIT_STRATEGY: none
only:
variables: [ $RUN_EXPERIMENTS == 'none' ]
image: ${REGISTRY}/${IMAGE_NAME_MLUQ}
script:
- cd /mpp/build/
- python3 mppyrun.py -n 2 --mpi_tests=1 --mute=1
dependencies: [ "build-mluq" ]
tags: [ docker ]
mpitest-4-mluq:
stage: test
variables:
GIT_STRATEGY: none
......@@ -64,10 +94,25 @@ mpitest-mluq:
image: ${REGISTRY}/${IMAGE_NAME_MLUQ}
script:
- cd /mpp/build/
- python3 mppyrun.py --mpi_tests=1 --mute=0
- python3 mppyrun.py -n 4 --mpi_tests=1 --mute=1
dependencies: [ "build-mluq" ]
tags: [ docker ]
mpitest-8-mluq:
stage: test
variables:
GIT_STRATEGY: none
only:
variables: [ $RUN_EXPERIMENTS == 'none' ]
image: ${REGISTRY}/${IMAGE_NAME_MLUQ}
script:
- cd /mpp/build/
- python3 mppyrun.py -n 8 --mpi_tests=1 --mute=1
dependencies: [ "build-mluq" ]
tags: [ docker ]
deploy-mluq:
stage: deploy
variables:
......@@ -79,10 +124,12 @@ deploy-mluq:
${REGISTRY}/mluq-${CI_COMMIT_TAG}-${OS}${OS_VERSION_MLUQ}
- docker push ${REGISTRY}/mluq-${CI_COMMIT_TAG}-${OS}${OS_VERSION_MLUQ}
only: [ tags ]
dependencies: [ "ctest-mluq", "mpitest-mluq" ]
dependencies: [ "ctest-mluq", "mpitest-1-mluq", "mpitest-2-mluq",
"mpitest-4-mluq", "mpitest-8-mluq"]
tags: [ docker ]
# ------------------------- End of CI -------------------------
# ------------------------ Start of CD ------------------------
# See https://www.nhr.kit.edu/userdocs/ci/ci-level2
......
......@@ -26,7 +26,7 @@ public:
config.get("MainVerbose", verbose);
config.get("Estimator", estimatorName);
estimator = EstimatorCreator(estimatorName).
estimator = EstimatorCreator(EstimatorEnum(estimatorName)).
CreateUnique();
}
......
......@@ -6,19 +6,25 @@
Estimator *EstimatorCreator::Create() {
if (_estimatorName == "StochasticCollocation")
if (_estimator == MULTILEVEL_STOCHASTIC_COLLOCATION)
return new MultilevelEstimator(
_epsilon, _levelVec, _samplesVec, STOCHASTIC_COLLOCATION, _parallel,
_meshesCreator, _pdeSolverCreator
);
if (_estimator == STOCHASTIC_COLLOCATION)
return new StochasticCollocation(
_epsilon, _level, _stochLevel, _onlyFine, _parallel,
_meshesCreator, _pdeSolverCreator
);
if (_estimatorName == "MultilevelMonteCarlo")
if (_estimator == MULTILEVEL_MONTE_CARLO)
return new MultilevelEstimator(
_epsilon, _levelVec, _samplesVec, _subEstName, _parallel,
_epsilon, _levelVec, _samplesVec, MONTE_CARLO, _parallel,
_meshesCreator, _pdeSolverCreator
);
if (_estimatorName == "MonteCarlo")
if (_estimator == MONTE_CARLO)
return new MonteCarlo(
_epsilon, _level, _initSamples, _onlyFine, _parallel,
_meshesCreator, _pdeSolverCreator
......@@ -31,3 +37,15 @@ std::unique_ptr<Estimator> EstimatorCreator::CreateUnique() {
return std::unique_ptr<Estimator>(Create());
}
ESTIMATOR EstimatorEnum(const string &estimatorName) {
if (estimatorName == "MC" || estimatorName == "MonteCarlo")
return MONTE_CARLO;
else if (estimatorName == "MLMC" || estimatorName == "MultilevelMonteCarlo")
return MULTILEVEL_MONTE_CARLO;
else if (estimatorName == "SC" || estimatorName == "StochasticCollocation")
return STOCHASTIC_COLLOCATION;
else if (estimatorName == "MLSC" || estimatorName == "MultilevelStochasticCollocation")
return MULTILEVEL_STOCHASTIC_COLLOCATION;
else return NO_ESTIMATOR;
}
#ifndef ESTIMATOR_HPP
#define ESTIMATOR_HPP
#include <utility>
#include "WelfordAggregate.hpp"
#include "PDESolverCreator.hpp"
#include "Exponents.hpp"
......@@ -10,28 +12,35 @@
#include "MeshesCreator.hpp"
enum ESTIMATOR {
MONTE_CARLO = 1,
STOCHASTIC_COLLOCATION = 2,
MULTILEVEL_MONTE_CARLO = 3,
MULTILEVEL_STOCHASTIC_COLLOCATION = 4,
NO_ESTIMATOR = 5,
};
ESTIMATOR EstimatorEnum(const std::string &estimatorName);
class Estimator {
protected:
int verbose = 1;
int level;
double epsilon;
double initSamples;
double cost = 0.0;
double value;
double value = 0.0;
double cost;
double epsilon = 0.0;
Errors errors;
Exponents exponents;
public:
WelfordAggregate aggregate;
Estimator(double epsilon, int level, int initSamples) :
epsilon(epsilon), level(level), initSamples(initSamples),
aggregate(WelfordAggregate(initSamples)) {}
Estimator(double epsilon, WelfordAggregate aggregate) :
epsilon(epsilon), aggregate(aggregate) {}
double TotalError() const { return errors.total; }
......@@ -59,40 +68,47 @@ public:
virtual void MultilevelResults() const {};
virtual void ExponentResults() const {};
virtual Exponents GetExponents() const { return Exponents(); };
};
/*
* Todo use initMap
*/
typedef std::vector<int> Levels;
typedef std::vector<int> Samples;
class EstimatorCreator {
private:
int _level;
int _level{};
int _initSamples;
int _stochLevel{};
int _stochLevel;
double _epsilon;
int _initSamples{};
Levels _levelVec;
Samples _samplesVec;
double _epsilon = 0.0;
bool _parallel = true;
bool _onlyFine = false;
std::string _subEstName = "MonteCarlo";
std::string _estimatorName = "MonteCarlo";
MeshesCreator _meshesCreator;
PDESolverCreator _pdeSolverCreator;
ESTIMATOR _estimator = MONTE_CARLO;
ESTIMATOR _subEstimator = MONTE_CARLO;
public:
EstimatorCreator(const std::string &estimatorName) : _estimatorName(estimatorName) {
explicit EstimatorCreator(ESTIMATOR estimator) :
_estimator(estimator) {
config.get("Level", _level);
config.get("epsilon", _epsilon);
config.get("OnlyFine", _onlyFine);
......@@ -143,12 +159,12 @@ public:
}
EstimatorCreator WithMeshesCreator(MeshesCreator meshesCreator) {
_meshesCreator = meshesCreator;
_meshesCreator = std::move(meshesCreator);
return *this;
}
EstimatorCreator WithPDESolverCreator(PDESolverCreator pdeSolverCreator) {
_pdeSolverCreator = pdeSolverCreator;
_pdeSolverCreator = std::move(pdeSolverCreator);
return *this;
}
......@@ -165,6 +181,10 @@ struct EstimatorMap : public LevelMap<Estimator *> {
config.get("maxLevel", maxLevel);
};
EstimatorMap(int minLevel, int maxLevel) {
config.get("maxLevel", maxLevel);
};
EstimatorMap(std::initializer_list<std::pair<int, Estimator *>> estimatorMap) :
LevelMap<Estimator *>(estimatorMap) {
config.get("maxLevel", maxLevel);
......@@ -172,31 +192,41 @@ struct EstimatorMap : public LevelMap<Estimator *> {
void UpdateSampleCounter(double epsilon);
void AppendLevel(double epsilon, Exponents exponents, int newLevel);
void UpdateNewLevel(double epsilon, Exponents exponents, int newLevel);
void SetEstimatorMapData(const MultilevelData &_data) {
for (auto &[level, estimator] : _levelMap) {
estimator->aggregate.ctr = _data.ctrs.at(level);
estimator->aggregate.mean = _data.means.at(level);
estimator->aggregate.sVar = _data.sVars.at(level);
estimator->aggregate.kurtosis = _data.kurtosis.at(level);
for (auto &[level, ctr] : _data.ctrs) {
this->at(level)->aggregate.ctr = _data.ctrs.at(level);
this->at(level)->aggregate.mean = _data.means.at(level);
this->at(level)->aggregate.sVar = _data.sVars.at(level);
this->at(level)->aggregate.kurtosis = _data.kurtosis.at(level);
}
}
void Fill(Levels initLevels, Samples initSamples,
const std::string &estimatorName, bool parallel,
MeshesCreator meshesCreator = MeshesCreator(),
PDESolverCreator pdeCreator = PDESolverCreator()) {
void Fill(Levels initLevels, Samples initSamples, ESTIMATOR estimator, bool parallel,
const MeshesCreator &meshesCreator = MeshesCreator(),
const PDESolverCreator &pdeCreator = PDESolverCreator()) {
for (auto level_index = 0; level_index < initLevels.size(); level_index++)
this->insert({initLevels[level_index],
EstimatorCreator(estimator).
WithInitSamples(initSamples[level_index]).
WithInitLevel(initLevels[level_index]).
WithPDESolverCreator(pdeCreator).
WithMeshesCreator(meshesCreator).
WithOnlyFine((level_index == 0)).
WithParallel(parallel).
Create()});
for (auto i = 0; i < initLevels.size(); i++)
this->insert({initLevels[i],
EstimatorCreator(estimatorName).
for (auto level = initLevels.back(); level <= maxLevel; level++)
this->insert({level,
EstimatorCreator(estimator).
WithPDESolverCreator(pdeCreator).
WithMeshesCreator(meshesCreator).
WithInitSamples(initSamples[i]).
WithInitLevel(initLevels[i]).
WithOnlyFine((i == 0)).
WithParallel(parallel).
WithInitLevel(level).
WithOnlyFine(false).
WithInitSamples(0).
Create()});
}
};
......
......@@ -2,22 +2,38 @@
void MonteCarlo::Method() {
mout.StartBlock("MonteCarlo l=" + to_string(level));
vout(1) << "Start with: " << aggregate;
method();
aggregate.UpdateParallel();
vout(1) << "End with: " << aggregate;
std::string parallelStr = aggregate.parallel ? "Parallel" : "Seriell";
mout.StartBlock(parallelStr + " MC l=" + to_string(level));
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;
}
value = aggregate.mean.Q;
cost = aggregate.mean.C;
// errors = 0.0;
mout.EndBlock(verbose == 0);
}
void MonteCarlo::updateSampleIds(SampleSolution &fSol, SampleSolution &cSol) {
fineId.number = aggregate.index();
coarseId.number = aggregate.index();
fineId.number = aggregate.Index();
coarseId.number = aggregate.Index();
fSol.id = fineId;
cSol.id = coarseId;
}
void MonteCarlo::method() {
// WithPLevel(onlyFine ? level : level - 1)
std::unique_ptr<Meshes> meshes = meshesCreator.
WithMeshName(pdeSolverCreator.GetMeshName()).
WithCommSplit(aggregate.commSplit).
......@@ -25,6 +41,8 @@ void MonteCarlo::method() {
WithLevel(level).
CreateUnique();
meshes->PrintInfo();
std::unique_ptr<PDESolver> pdeSolver = pdeSolverCreator.
CreateUnique(*meshes);
......@@ -43,10 +61,8 @@ void MonteCarlo::method() {
}
aggregate.Update(fineSolution, coarseSolution);
// errors.Estimate(aggregate);
//
// if (errors.stochastic > epsilon && epsilon != 0.0)
// mout << "TODO" << endl;
// aggregate.
}
}
void MonteCarlo::adaptiveMethod() {
}
......@@ -14,9 +14,9 @@ class MonteCarlo : public Estimator {
protected:
int verbose = 1;
bool onlyFine;
int level;
bool parallel;
bool onlyFine;
SampleID fineId;
......@@ -30,19 +30,20 @@ protected:
void method();
void adaptiveMethod();
public:
MonteCarlo(double epsilon, int level, int initSamples, bool onlyFine, bool parallel,
MeshesCreator meshesCreator = MeshesCreator(),
PDESolverCreator pdeCreator = PDESolverCreator()) :
Estimator(epsilon, level, initSamples), onlyFine(onlyFine), parallel(parallel),
fineId(SampleID(level, 0, false)), coarseId(SampleID(level, 0, true)),
meshesCreator(std::move(meshesCreator)), pdeSolverCreator(std::move(pdeCreator)) {
Estimator(epsilon, WelfordAggregate(initSamples, parallel)),
onlyFine(onlyFine), level(level),
fineId(SampleID(level, 0, false)),
coarseId(SampleID(level, 0, true)),
meshesCreator(std::move(meshesCreator)),
pdeSolverCreator(std::move(pdeCreator)) {
config.get("MCVerbose", verbose);
config.get("MCParallel", parallel);
aggregate.parallel = parallel;
aggregate.UpdateSampleCounter(initSamples);
}
std::string Name() const override { return "MonteCarlo"; }
......
......@@ -2,6 +2,7 @@
void MultilevelEstimator::Method() {
std::string parallelStr = aggregate.parallel ? "Parallel" : "Seriell";
mout.StartBlock("Multilevel Method");
if (epsilon == 0.0) {
......@@ -24,7 +25,8 @@ void MultilevelEstimator::Method() {
void MultilevelEstimator::method() {
for (auto &[level, estimator] : estMap)
estimator->Method();
if (estimator->aggregate.ctr.dM > 0)
estimator->Method();
}
void MultilevelEstimator::adaptiveMethod() {
......@@ -46,7 +48,7 @@ void MultilevelEstimator::adaptiveMethod() {
vout(1) << "numErr=" << mlErrors.numeric << endl;
if (mlErrors.numeric > epsilon / 2) {
estMap.AppendLevel(epsilon, exponents, estMap.rbegin()->first + 1);
estMap.UpdateNewLevel(epsilon, exponents, estMap.rbegin()->first + 1);
} else {
mlErrors.Estimate(data);
vout(1) << "statErr=" << mlErrors.stochastic << endl;
......@@ -66,3 +68,7 @@ void MultilevelEstimator::ExponentResults() const {
PrintInfoEntry("Final beta", exponents.beta),
PrintInfoEntry("Final gamma", exponents.gamma));
}
Exponents MultilevelEstimator::GetExponents() const {
return exponents;
};
......@@ -30,13 +30,15 @@ public:
MultilevelErrors mlErrors;
MultilevelEstimator(double epsilon, Levels initLevels, Samples initSamples,
const std::string &estName, bool parallel,
ESTIMATOR subEstimator, bool parallel,
MeshesCreator meshesCreator = MeshesCreator(),
PDESolverCreator pdeCreator = PDESolverCreator()) :
Estimator(epsilon, initLevels.back(), initSamples.front()) {
Estimator(epsilon, WelfordAggregate(0, parallel)) {
config.get("MLMCVerbose", verbose);
estMap.Fill(initLevels, initSamples, estName, parallel, meshesCreator, pdeCreator);
estMap.Fill(initLevels, initSamples,
subEstimator, parallel,
meshesCreator, pdeCreator);
}
~MultilevelEstimator() {
......@@ -52,6 +54,8 @@ public:
void MultilevelResults() const override;
void ExponentResults() const override;
Exponents GetExponents() const override;
};
#endif //MULTILEVELESTIMATOR_HPP
......@@ -37,8 +37,8 @@ void StochasticCollocation::method() {
}
void StochasticCollocation::updateSampleIds(SampleSolution &fSol, SampleSolution &cSol) {
fineId.number = aggregate.index();
coarseId.number = aggregate.index();
fineId.number = aggregate.Index();
coarseId.number = aggregate.Index();
fSol.id = fineId;
cSol.id = coarseId;
}
......@@ -10,9 +10,9 @@ class StochasticCollocation : public Estimator {
protected:
int verbose = 1;
bool onlyFine;
int level;
bool parallel;
bool onlyFine;
SampleID fineId;
......@@ -27,11 +27,12 @@ protected:
void method();
public:
StochasticCollocation(double epsilon, int level, int stochLevel,
bool onlyFine, bool parallel,
StochasticCollocation(double epsilon, int level, int stochLevel, bool onlyFine,
bool parallel,
MeshesCreator meshesCreator = MeshesCreator(),
PDESolverCreator pdeCreator = PDESolverCreator()) :
Estimator(epsilon, level, 0),
Estimator(epsilon, WelfordAggregate(0, parallel)),
onlyFine(onlyFine), level(level),
fineId(SampleID(level, 0, false)), coarseId(SampleID(level, 0, true)),
meshesCreator(std::move(meshesCreator)), pdeSolverCreator(std::move(pdeCreator)) {}
......
#include "Exponents.hpp"
void Exponents::ComputeExponents(const EstimatorMap &mcMap) {
ComputeExponents(MultilevelData(mcMap));
ComputeExponents(MultilevelData(mcMap));
}
void Exponents::ComputeExponents(const MultilevelData &data) {
ComputeExponents(data.means, data.sVars);
ComputeExponents(data.means, data.sVars);
}
void Exponents::ComputeExponents(const MeanMap &mean, const SVarMap &sVar) {
alpha = ComputeAlpha(mean);
beta = ComputeBeta(sVar);
gamma = ComputeGamma(mean);
alpha = ComputeAlpha(mean);
beta = ComputeBeta(sVar);
gamma = ComputeGamma(mean);
}
double Exponents::ComputeAlpha(const MeanMap &avgs) {
auto _levels = avgs.GetLevelVector();
auto _avgsY = avgs.GetYVector();
auto _levels = avgs.GetLevelVector();
auto _avgsY = avgs.GetYVector();
vector<double> levels, avgsY;
for (int i = 1; i < _levels.size(); i++) {
levels.push_back(double(_levels[i]));
avgsY.push_back(-log2(_avgsY[i]));
}
vector<double> levels, avgsY;
for (int i = 1; i < _levels.size(); i++) {
if (_avgsY[i] == 0.0) continue;
levels.push_back(double(_levels[i]));
avgsY.push_back(-log2(_avgsY[i]));
}
return linearFit(levels, avgsY).first;
return linearFit(levels, avgsY).first;
}
double Exponents::ComputeBeta(const SVarMap &vars) {
auto _levels = vars.GetLevelVector();
auto _varsY = vars.GetYVector();