Skip to content

GitLab

Commit 65a65deb authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

major restructuring

parent d52e044f
Pipeline #149253 failed with stages
in 21 minutes and 10 seconds
Hide whitespace changes
Inline Side-by-side
CMakeLists.txt
View file @ 65a65deb
......@@ -5,8 +5,8 @@ set(SPACE_DIM 2 CACHE STRING "SPACE_DIM")
set(USE_FFTW ON CACHE STRING "USE_FFTW")
set(USE_SPLITTED_COMMS ON CACHE STRING "USE_SPLITTED_COMMS")
set(USE_SPACETIME OFF CACHE STRING "USE_SPACETIME")
#set(COMPILER_OPTIMIZE -O0 CACHE STRING "COMPILER_OPTIMIZE")
set(COMPILER_OPTIMIZE -O3 CACHE STRING "COMPILER_OPTIMIZE")
set(COMPILER_OPTIMIZE -O0 CACHE STRING "COMPILER_OPTIMIZE")
#set(COMPILER_OPTIMIZE -O3 CACHE STRING "COMPILER_OPTIMIZE")
set(NO_DEPRECATED OFF CACHE STRING "NO_DEPRECATED")
set(AFFINE_LINEAR_TRAFO ON CACHE STRING "AFFINE_LINEAR_TRAFO")
set(BUILD_TESTS OFF CACHE STRING "BUILD_TESTS")
......
mlmc/src/Main.hpp
View file @ 65a65deb
#ifndef MAIN_HPP
#define MAIN_HPP
#include "MultilevelStochasticCollocation.hpp"
#include "StochasticCollocation.hpp"
#include "MultilevelMonteCarlo.hpp"
#include "MultilevelEstimator.hpp"
#include "MonteCarlo.hpp"
#include "Estimator.hpp"
......@@ -14,18 +13,6 @@ bool PowerOfTwo(int n) {
return (ceil(log2(n)) == floor(log2(n)));
}
std::unique_ptr<Estimator> CreateEstimator(const std::string &estimator) {
if (estimator == "MultilevelStochasticCollocation")
return std::make_unique<MultilevelStochasticCollocation>();
if (estimator == "StochasticCollocation")
return std::make_unique<StochasticCollocation>();
if (estimator == "MultilevelMonteCarlo")
return std::make_unique<MultilevelMonteCarlo>();
if (estimator == "MonteCarlo")
return std::make_unique<MonteCarlo>();
else Exit("Estimator is not implemented")
}
class MainProgram {
private:
int verbose = 1;
......@@ -39,7 +26,7 @@ public:
config.get("MainVerbose", verbose);
config.get("Estimator", estimatorName);
estimator = CreateEstimator(estimatorName);
estimator = CreateEstimatorUniquePtr(estimatorName);
}
int Run() {
......
mlmc/src/basics/Level.hpp deleted 100644 → 0
View file @ d52e044f
#ifndef LEVEL_HPP
#define LEVEL_HPP
#include "Logging.hpp"
#include <vector>
#include <map>
/*
* Todo: make second key template for degree
* to realize MultidegreeMonteCarlo
*/
template<typename T>
struct LevelMap {
protected:
std::map<int, T> _levelMap;
public:
LevelMap() {};
LevelMap(std::initializer_list<std::pair<int, T>> levelMap) {
for (auto const &pair : levelMap) {
_levelMap.emplace(pair.first, T(pair.second));
}
}
std::vector<int> GetLevelVector() const {
std::vector<int> levelVector;
for (auto const &[level, entry]: _levelMap)
levelVector.push_back(level);
return levelVector;
}
std::vector<double> GetQVector() const {
std::vector<double> qVector;
for (auto const &[level, entry]: _levelMap)
qVector.push_back(entry.Q);
return qVector;
}
std::vector<double> GetYVector() const {
std::vector<double> yVector;
for (auto const &[level, entry]: _levelMap)
yVector.push_back(entry.Y);
return yVector;
}
std::vector<double> GetCostVector() const {
std::vector<double> costVector;
for (auto const &[level, entry]: _levelMap)
costVector.push_back(entry.C);
return costVector;
}
std::vector<int> GetMVector() const {
std::vector<int> mVector;
for (auto const &[level, entry]: _levelMap)
mVector.push_back(entry.M);
return mVector;
}
std::vector<int> GetdMVector() const {
std::vector<int> dMVector;
for (auto const &[level, entry]: _levelMap)
dMVector.push_back(entry.dM);
return dMVector;
}
auto size() const { return _levelMap.size(); }
auto clear() { _levelMap.clear(); }
auto insert(std::pair<int, T> pair) { _levelMap.insert(pair); }
auto begin() { return _levelMap.begin(); }
auto rbegin() { return _levelMap.rbegin(); }
auto end() { return _levelMap.end(); }
auto rend() { return _levelMap.end(); }
auto begin() const { return _levelMap.begin(); }
auto rbegin() const { return _levelMap.rbegin(); }
auto end() const { return _levelMap.end(); }
auto rend() const { return _levelMap.end(); }
auto operator[](int level) { return _levelMap[level]; }
auto at(int level) { return _levelMap.at(level); }
auto find(int level) { return _levelMap.find(level); }
friend Logging &operator<<(Logging &s, const LevelMap<T> &levelMap) {
s << "[";
for (auto &[level, entry] : levelMap) {
if (level != levelMap.rbegin()->first) s << entry << ", ";
else s << entry;
}
return s << "]";
}
};
#endif //LEVEL_HPP
mlmc/src/basics/Sample.hpp
View file @ 65a65deb
#ifndef SAMPLE_HPP
#define SAMPLE_HPP
#include "Level.hpp"
#include "Algebra.hpp"
#include <string>
......
mlmc/src/estimators/CMakeLists.txt
View file @ 65a65deb
add_library(MONTECARLO STATIC
Estimator.cpp
MonteCarlo.cpp
MultilevelMonteCarlo.cpp
MultilevelEstimator.cpp
StochasticCollocation.cpp
datastructure/MultilevelErrors.cpp
datastructure/Exponents.cpp
datastructure/WelfordAggregate.cpp
datastructure/EmpiricMeasureLevelMaps.cpp
datastructure/LevelMaps.cpp
)
target_link_libraries(MONTECARLO PDESOLVER MPP_LIBRARIES)
mlmc/src/estimators/MultilevelMonteCarlo.cpp → mlmc/src/estimators/MultilevelEstimator.cpp
View file @ 65a65deb
#include "MultilevelMonteCarlo.hpp"
#include "MultilevelEstimator.hpp"
void MultilevelMonteCarlo::initializeMonteCarloMap() {
mcMap.clear();
void MultilevelEstimator::initializeCollocationMap() {
}
void MultilevelEstimator::initializeMonteCarloMap() {
estimatorMap.clear();
for (unsigned long i = 0; i < initLevels.size(); i++) {
bool onlyFine = (i == 0);
int M = initSampleAmount[i];
mcMap.insert({initLevels[i], MonteCarlo(initLevels[i], M, onlyFine, true)});
estimatorMap.insert(
{initLevels[i], new MonteCarlo(initLevels[i], M, onlyFine, true)}
);
}
}
void MultilevelMonteCarlo::Method() {
void MultilevelEstimator::Method() {
mout.StartBlock("MLMC Method");
if (epsilon == 0.0) {
......@@ -21,26 +26,26 @@ void MultilevelMonteCarlo::Method() {
adaptiveMethod();
}
data.ExtractDataFrom(mcMap);
data.ExtractDataFrom(estimatorMap);
exponents.ComputeExponents(data);
mout.EndBlock(verbose == 0);
}
void MultilevelMonteCarlo::method() {
for (auto &[level, mc] : mcMap)
mc.Method();
void MultilevelEstimator::method() {
for (auto &[level, estimator] : estimatorMap)
estimator->Method();
}
void MultilevelMonteCarlo::adaptiveMethod() {
void MultilevelEstimator::adaptiveMethod() {
bool converged = false;
while (!converged) {
for (auto &[level, mc] : mcMap)
if (mc.aggregate.ctr.dM > 0)
mc.Method();
for (auto &[level, estimator] : estimatorMap)
if (estimator->aggregate.ctr.dM > 0)
estimator->Method();
mcMap.UpdateSampleCounter(epsilon);
data.ExtractDataFrom(mcMap);
estimatorMap.UpdateSampleCounter(epsilon);
data.ExtractDataFrom(estimatorMap);
exponents.ComputeExponents(data);
vout(1) << exponents;
......@@ -51,7 +56,7 @@ void MultilevelMonteCarlo::adaptiveMethod() {
vout(1) << "numErr=" << errors.numeric << endl;
if (errors.numeric > epsilon / 2) {
mcMap.AppendLevel(epsilon, exponents, mcMap.rbegin()->first + 1);
estimatorMap.AppendLevel(epsilon, exponents, estimatorMap.rbegin()->first + 1);
} else {
errors.Estimate(data);
vout(1) << "statErr=" << errors.stochastic << endl;
......@@ -61,17 +66,17 @@ void MultilevelMonteCarlo::adaptiveMethod() {
}
}
void MultilevelMonteCarlo::PrintInfo() const {
mout.PrintInfo("MLMC", verbose,
void MultilevelEstimator::PrintInfo() const {
mout.PrintInfo("MultilevelEstimator", verbose,
PrintInfoEntry("initLevels", vec2str(initLevels)),
PrintInfoEntry("initSample", vec2str(initSampleAmount)));
}
void MultilevelMonteCarlo::MultilevelResults() const {
void MultilevelEstimator::MultilevelResults() const {
data.MultilevelResults();
}
void MultilevelMonteCarlo::EstimatorResults() const {
void MultilevelEstimator::EstimatorResults() const {
mout.PrintInfo("MLMC Results", verbose,
PrintInfoEntry("Value", data.means.Q),
PrintInfoEntry("Cost", data.means.Cost),
......@@ -81,7 +86,7 @@ void MultilevelMonteCarlo::EstimatorResults() const {
PrintInfoEntry("Numerical Error", errors.numeric));
}
void MultilevelMonteCarlo::ExponentResults() const {
void MultilevelEstimator::ExponentResults() const {
mout.PrintInfo("Exponents", verbose,
PrintInfoEntry("Final alpha", exponents.alpha),
PrintInfoEntry("Final beta", exponents.beta),
......
mlmc/src/estimators/MultilevelMonteCarlo.hpp → mlmc/src/estimators/MultilevelEstimator.hpp
View file @ 65a65deb
#ifndef MULTILEVELMONTECARLO_HPP
#define MULTILEVELMONTECARLO_HPP
#ifndef MULTILEVELESTIMATOR_HPP
#define MULTILEVELESTIMATOR_HPP
#include "EmpiricMeasureLevelMaps.hpp"
#include "LevelMaps.hpp"
#include "MonteCarlo.hpp"
#include "Exponents.hpp"
#include "MultilevelErrors.hpp"
......@@ -9,18 +9,42 @@
#include "Estimator.hpp"
class MultilevelMonteCarlo : public Estimator {
struct EstimatorMap : public LevelMap<Estimator *> {
int maxLevel = 10;
EstimatorMap() {
config.get("maxLevel", maxLevel);
};
EstimatorMap(std::initializer_list<std::pair<int, Estimator *>> estimatorMap) :
LevelMap<Estimator *>(estimatorMap) {
config.get("maxLevel", maxLevel);
};
void UpdateSampleCounter(double epsilon);
void AppendLevel(double epsilon, Exponents exponents, int newLevel);
};
class MultilevelEstimator : public Estimator {
private:
int verbose = 1;
// Levels init{3, 4, 5};
std::vector<int> initLevels{3, 4, 5};
// Samples init{12, 6, 3};
std::vector<int> initSampleAmount{12, 6, 3};
void initializeCollocationMap();
void initializeMonteCarloMap();
void adaptiveMethod();
void method();
public:
EstimatorMap mcMap;
EstimatorMap estimatorMap;
MultilevelData data;
......@@ -28,7 +52,7 @@ public:
MultilevelErrors errors;
MultilevelMonteCarlo() : Estimator() {
MultilevelEstimator() : Estimator() {
config.get("MLMCVerbose", verbose);
config.get("initLevels", initLevels);
config.get("initSampleAmount", initSampleAmount);
......@@ -36,8 +60,12 @@ public:
initializeMonteCarloMap();
}
~MultilevelMonteCarlo() {
mcMap.clear();
MultilevelEstimator(const EstimatorMap &mcMap) : Estimator(), estimatorMap(mcMap) {
config.get("MLMCVerbose", verbose);
}
~MultilevelEstimator() {
estimatorMap.clear();
}
std::string Name() const override {
......@@ -48,22 +76,11 @@ public:
void PrintInfo() const;
double TotalError() const override {
return errors.total;
}
void EstimatorResults() const override;
void MultilevelResults() const override;
void ExponentResults() const override;
private:
void method();
void adaptiveMethod();
void initializeMonteCarloMap();
};
#endif //MULTILEVELMONTECARLO_HPP
#endif //MULTILEVELESTIMATOR_HPP
mlmc/src/estimators/datastructure/Exponents.hpp
View file @ 65a65deb
#ifndef EXPONENTS_HPP
#define EXPONENTS_HPP
#include "EmpiricMeasureLevelMaps.hpp"
#include "LevelMaps.hpp"
struct Exponents {
......
mlmc/src/estimators/datastructure/EmpiricMeasureLevelMaps.cpp → mlmc/src/estimators/datastructure/LevelMaps.cpp
View file @ 65a65deb
#include "EmpiricMeasureLevelMaps.hpp"
#include "Exponents.hpp"
#include "LevelMaps.hpp"
#include "MultilevelEstimator.hpp"
void EstimatorMap::UpdateSampleCounter(double epsilon) {
......@@ -7,28 +7,28 @@ void EstimatorMap::UpdateSampleCounter(double epsilon) {
double factor = 0.0;
for (auto &[level, estimator] : _levelMap)
factor += sqrt(estimator.aggregate.sVar.Y * estimator.aggregate.mean.C);
factor += sqrt(estimator->aggregate.sVar.Y * estimator->aggregate.mean.C);
for (auto &[level, estimator] : _levelMap) {
optimalM = (int) (ceil(2 * pow(epsilon, -2) * factor *
sqrt(estimator.aggregate.sVar.Y / estimator.aggregate.mean.C)));
sqrt(estimator->aggregate.sVar.Y / estimator->aggregate.mean.C)));
if (optimalM == 1) optimalM++; // Hack
estimator.aggregate.UpdateSampleCounter(optimalM - estimator.aggregate.ctr.M);
estimator->aggregate.UpdateSampleCounter(optimalM - estimator->aggregate.ctr.M);
}
}
void EstimatorMap::AppendLevel(double epsilon, Exponents exponents, int newLevel) {
if (newLevel <= maxLevel) {
auto oldLevel = this->rbegin();
double sVarY = oldLevel->second.aggregate.sVar.Y / pow(2.0, exponents.beta);
double avgsCost = oldLevel->second.aggregate.mean.C * pow(2.0, exponents.gamma);
double sVarY = oldLevel->second->aggregate.sVar.Y / pow(2.0, exponents.beta);
double avgsCost = oldLevel->second->aggregate.mean.C * pow(2.0, exponents.gamma);
_levelMap.insert(
{newLevel, MonteCarlo(newLevel, 0, false, true)}
{newLevel, new MonteCarlo(newLevel, 0, false, true)}
);
_levelMap.find(newLevel)->second.aggregate.sVar.Y = sVarY;
_levelMap.find(newLevel)->second.aggregate.mean.C = avgsCost;
_levelMap.find(newLevel)->second->aggregate.sVar.Y = sVarY;
_levelMap.find(newLevel)->second->aggregate.mean.C = avgsCost;
this->UpdateSampleCounter(epsilon);
} else {
......@@ -38,7 +38,7 @@ void EstimatorMap::AppendLevel(double epsilon, Exponents exponents, int newLevel
void SampleCounterMap::Update(const EstimatorMap &estimatorMap) {
for (auto &[level, estimator]: estimatorMap)
_levelMap[level] = SampleCounter(estimator.aggregate.ctr);
_levelMap[level] = SampleCounter(estimator->aggregate.ctr);
}
bool SampleCounterMap::NoSamplesLeft() {
......@@ -53,18 +53,18 @@ bool SampleCounterMap::NoSamplesLeft() {
void MeanMap::Update(const EstimatorMap &estimatorMap) {
for (auto &&[level, estimator]: estimatorMap) {
_levelMap[level] = estimator.aggregate.mean;
Q += estimator.aggregate.mean.Y;
Cost += estimator.aggregate.mean.C; // Todo at the right place?
_levelMap[level] = estimator->aggregate.mean;
Q += estimator->aggregate.mean.Y;
Cost += estimator->aggregate.mean.C; // Todo at the right place?
}
}
void SVarMap::Update(const EstimatorMap &estimatorMap) {
for (auto &[level, estimator]: estimatorMap)
_levelMap[level] = estimator.aggregate.sVar;
_levelMap[level] = estimator->aggregate.sVar;
}
void KurtosisMap::Update(const EstimatorMap &estimatorMap) {
for (auto &[level, estimator]: estimatorMap)
_levelMap[level] = estimator.aggregate.kurtosis;
_levelMap[level] = estimator->aggregate.kurtosis;
}
mlmc/src/estimators/datastructure/EmpiricMeasureLevelMaps.hpp → mlmc/src/estimators/datastructure/LevelMaps.hpp
View file @ 65a65deb
#ifndef EMPIRICMEASURELEVELMAPS_HPP
#define EMPIRICMEASURELEVELMAPS_HPP
#ifndef LEVELMAPS_HPP
#define LEVELMAPS_HPP
#include "WelfordAggregate.hpp"
#include "Level.hpp"
#include "MonteCarlo.hpp"
#include <math.h>
struct Exponents;
#include "Logging.hpp"
struct EstimatorMap : public LevelMap<Estimator> {
#include <vector>
#include <map>
int maxLevel = 10;
EstimatorMap() {
config.get("maxLevel", maxLevel);
};
/*
* Todo: make second key template for degree
* to realize MultidegreeMonteCarlo
*/
EstimatorMap(std::initializer_list<std::pair<int, Estimator>> mcMap) :
LevelMap<Estimator>(mcMap) {
config.get("maxLevel", maxLevel);
};
template<typename T>
struct LevelMap {
protected:
std::map<int, T> _levelMap;
public:
LevelMap() {};
~LevelMap() {
clear();
}
LevelMap(std::initializer_list<std::pair<int, T>> levelMap) {
for (auto const &pair : levelMap) {
_levelMap.emplace(pair.first, T(pair.second));
}
}
std::vector<int> GetLevelVector() const {
std::vector<int> levelVector;
for (auto const &[level, entry]: _levelMap)
levelVector.push_back(level);
return levelVector;
}
std::vector<double> GetQVector() const {
std::vector<double> qVector;
for (auto const &[level, entry]: _levelMap)
qVector.push_back(entry.Q);
return qVector;
}
std::vector<double> GetYVector() const {
std::vector<double> yVector;
for (auto const &[level, entry]: _levelMap)
yVector.push_back(entry.Y);
return yVector;
}
std::vector<double> GetCostVector() const {
std::vector<double> costVector;
for (auto const &[level, entry]: _levelMap)
costVector.push_back(entry.C);
return costVector;
}
std::vector<int> GetMVector() const {
std::vector<int> mVector;
for (auto const &[level, entry]: _levelMap)
mVector.push_back(entry.M);
return mVector;
}
std::vector<int> GetdMVector() const {
std::vector<int> dMVector;
for (auto const &[level, entry]: _levelMap)
dMVector.push_back(entry.dM);
return dMVector;
}
void UpdateSampleCounter(double epsilon);
auto size() const { return _levelMap.size(); }
void AppendLevel(double epsilon, Exponents exponents, int newLevel);
auto clear() {
for (auto const &[level, entry] : _levelMap)
if constexpr (std::is_pointer<T>::value)
delete entry;
_levelMap.clear();
}
auto insert(std::pair<int, T> &&pair) { _levelMap.insert(pair); }
auto begin() { return _levelMap.begin(); }
auto rbegin() { return _levelMap.rbegin(); }
auto end() { return _levelMap.end(); }
auto rend() { return _levelMap.end(); }
auto begin() const { return _levelMap.begin(); }
auto rbegin() const { return _levelMap.rbegin(); }