Commit 85341af1 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

started to incorporate estimator interfaces

parent d7b0084f
Pipeline #149112 failed with stages
in 10 minutes and 36 seconds
#ifndef MAIN_HPP
#define MAIN_HPP
#include "MultilevelStochasticCollocation.hpp"
#include "StochasticCollocation.hpp"
#include "MultilevelMonteCarlo.hpp"
#include "MonteCarlo.hpp"
#include "IEstimator.hpp"
bool PowerOfTwo(int n) {
if(n==0) return false;
if (n == 0) return false;
return (ceil(log2(n)) == floor(log2(n)));
}
std::unique_ptr<IEstimator> 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;
std::string estimatorName = "MultilevelMonteCarlo";
public:
MultilevelMonteCarlo mlmc;
std::unique_ptr<IEstimator> estimator;
MainProgram() : mlmc(MultilevelMonteCarlo()) {
MainProgram() {
config.get("MainVerbose", verbose);
config.get("Estimator", estimatorName);
estimator = CreateEstimator(estimatorName);
}
int Run() {
config.PrintInfo();
mout.StartBlock("MLMC Experiment");
mout.StartBlock(estimator->Name() + " Experiment");
mout << "Run" << endl;
mlmc.Method();
estimator->Method();
mout.EndBlock();
mout << endl;
mlmc.PrintMCResults();
mlmc.PrintExponentResults();
mlmc.PrintMLMCResults();
estimator->MultilevelResults();
estimator->ExponentResults();
estimator->EstimatorResults();
return 0;
}
......
......@@ -15,95 +15,95 @@
template<typename T>
struct LevelMap {
protected:
std::map<int, T> _levelMap;
std::map<int, T> _levelMap;
public:
LevelMap() {};
LevelMap() {};
LevelMap(std::initializer_list<std::pair<int, T>> levelMap) {
for (auto const &pair : levelMap) {
_levelMap.emplace(pair.first, T(pair.second));
}
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<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> 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> 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<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> 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;
}
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 size() const { return _levelMap.size(); }
auto clear() { _levelMap.clear(); }
auto clear() { _levelMap.clear(); }
auto insert(std::pair<int, T> pair) { _levelMap.insert(pair); }
auto insert(std::pair<int, T> pair) { _levelMap.insert(pair); }
auto begin() { return _levelMap.begin(); }
auto begin() { return _levelMap.begin(); }
auto rbegin() { return _levelMap.rbegin(); }
auto rbegin() { return _levelMap.rbegin(); }
auto end() { return _levelMap.end(); }
auto end() { return _levelMap.end(); }
auto rend() { return _levelMap.end(); }
auto rend() { return _levelMap.end(); }
auto begin() const { return _levelMap.begin(); }
auto begin() const { return _levelMap.begin(); }
auto rbegin() const { return _levelMap.rbegin(); }
auto rbegin() const { return _levelMap.rbegin(); }
auto end() const { return _levelMap.end(); }
auto end() const { return _levelMap.end(); }
auto rend() const { return _levelMap.end(); }
auto rend() const { return _levelMap.end(); }
auto operator[](int level) { return _levelMap[level]; }
auto operator[](int level) { return _levelMap[level]; }
auto at(int level) { return _levelMap.at(level); }
auto at(int level) { return _levelMap.at(level); }
auto find(int level) { return _levelMap.find(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 << "]";
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
#ifndef IESTIMATOR_HPP
#define IESTIMATOR_HPP
//#include "Errors.hpp"
struct Exponents;
class IEstimator {
protected:
double epsilon = 0.0;
// Errors errors;
public:
virtual double TotalError() const {
return 0.0;
}
double Epsilon() const { return epsilon; }
virtual void Method() = 0;
virtual std::string Name() const = 0;
virtual void EstimatorResults() const = 0;
virtual void MultilevelResults() const {};
virtual void ExponentResults() const {};
};
class IMultilevelEstimator : public IEstimator {
public:
virtual void Method() = 0;
virtual void EstimatorResults() const = 0;
virtual void MultilevelResults() const {};
virtual void ExponentResults() const {};
};
#endif //IESTIMATOR_HPP
......@@ -3,11 +3,12 @@
#include "WelfordAggregate.hpp"
#include "PDESolverCreator.hpp"
#include "IEstimator.hpp"
#include "MeshesCreator.hpp"
class MonteCarlo {
class MonteCarlo : public IEstimator {
protected:
int verbose = 1;
......@@ -47,6 +48,27 @@ public:
Init(dM);
}
MonteCarlo() :
pdeSolverCreator(PDESolverCreator()),
meshesCreator(MeshesCreator(pdeSolverCreator.GetMeshName())) {
config.get("maxLevel", level);
}
std::string Name() const override {
return "MonteCarlo";
}
void EstimatorResults() const override {
mout.PrintInfo("MLMC Results", verbose,
PrintInfoEntry("Value", aggregate.mean.Q),
PrintInfoEntry("Cost", aggregate.mean.C * aggregate.ctr.M),
PrintInfoEntry("Epsilon", epsilon));
// PrintInfoEntry("Total Error", errors.total),
// PrintInfoEntry("Stochastical Error", errors.stochastic),
// PrintInfoEntry("Numerical Error", errors.numeric));
}
void Init(int dM) {
config.get("MCVerbose", verbose);
config.get("MCParallel", parallel);
......@@ -74,7 +96,7 @@ public:
if (!meshes) delete meshes;
}
void Method();
void Method() override;
private:
void method();
......
......@@ -67,7 +67,7 @@ void MultilevelMonteCarlo::PrintInfo() const {
PrintInfoEntry("initSample", vec2str(initSampleAmount)));
}
void MultilevelMonteCarlo::PrintMCResults() const {
void MultilevelMonteCarlo::MultilevelResults() const {
mout.PrintInfo("MC Results", verbose,
PrintInfoEntry("E(Qf-Qc)", vec2str(data.avgs.GetYVector())),
PrintInfoEntry("E(Qf)", vec2str(data.avgs.GetQVector())),
......@@ -79,7 +79,7 @@ void MultilevelMonteCarlo::PrintMCResults() const {
PrintInfoEntry("Used Samples", vec2str(data.ctrs.GetMVector())));
}
void MultilevelMonteCarlo::PrintMLMCResults() const {
void MultilevelMonteCarlo::EstimatorResults() const {
mout.PrintInfo("MLMC Results", verbose,
PrintInfoEntry("Value", data.avgs.Q),
PrintInfoEntry("Cost", data.avgs.Cost),
......@@ -89,7 +89,7 @@ void MultilevelMonteCarlo::PrintMLMCResults() const {
PrintInfoEntry("Numerical Error", errors.numeric));
}
void MultilevelMonteCarlo::PrintExponentResults() const {
void MultilevelMonteCarlo::ExponentResults() const {
mout.PrintInfo("Exponents", verbose,
PrintInfoEntry("Final alpha", exponents.alpha),
PrintInfoEntry("Final beta", exponents.beta),
......
#ifndef MULTILEVELMONTECARLO_HPP
#define MULTILEVELMONTECARLO_HPP
#include "MonteCarlo.hpp"
#include "EmpiricMeasureLevelMaps.hpp"
#include "MonteCarlo.hpp"
#include "Exponents.hpp"
#include "Errors.hpp"
#include "Utilities.hpp"
#include "IEstimator.hpp"
// This could be MultilevelMethod using either MC MCForMLMC Collocation : Estimator
class MultilevelMonteCarlo {
class MultilevelMonteCarlo : public IMultilevelEstimator {
private:
int verbose = 1;
......@@ -48,15 +48,23 @@ public:
mcMap.clear();
}
void Method();
std::string Name() const override {
return "MultilevelMonteCarlo";
}
void Method() override;
void PrintInfo() const;
void PrintMCResults() const;
double TotalError() const override {
return errors.total;
}
void EstimatorResults() const override;
void PrintMLMCResults() const;
void MultilevelResults() const override;
void PrintExponentResults() const;
void ExponentResults() const override;
private:
void method();
......
#ifndef MULTILEVELSTOCHASTICCOLLOCATION_HPP
#define MULTILEVELSTOCHASTICCOLLOCATION_HPP
class MultilevelStochasticCollocation {
#include "IEstimator.hpp"
class MultilevelStochasticCollocation : public IMultilevelEstimator {
private:
public:
std::string Name() const override {
return "MultilevelStochasticCollocation";
}
void Method() override {
// todo here is more stuff to be done
}
void EstimatorResults() const override {
// todo do stuffs
}
};
#endif //MULTILEVELSTOCHASTICCOLLOCATION_HPP
#ifndef STOCHASTICCOLLOCATION_HPP
#define STOCHASTICCOLLOCATION_HPP
class StochasticCollocation {
#include "IEstimator.hpp"
class StochasticCollocation : public IEstimator {
private:
public:
std::string Name() const override {
return "StochasticCollocation";
}
void Method() override {
// Todo here is some stuff todo
}
void EstimatorResults() const override {
// todo do stuffs
}
};
#endif //STOCHASTICCOLLOCATION_HPP
......@@ -2,50 +2,49 @@
void Errors::Estimate(const MonteCarloMap &mcMap) {
Estimate(MultiLevelMonteCarloData(mcMap));
Estimate(MultiLevelMonteCarloData(mcMap));
}
void Errors::Estimate(const MultiLevelMonteCarloData &data) {
exponents.ComputeExponents(data);
numeric = EstimateNumeric(data);
stochastic = EstimateStochastic(data);
total = numeric + stochastic;
exponents.ComputeExponents(data);
numeric = EstimateNumeric(data);
stochastic = EstimateStochastic(data);
total = numeric + stochastic;
}
double Errors::EstimateNumeric(const MultiLevelMonteCarloData &data) {
exponents.ComputeExponents(data);
return EstimateNumeric(data.avgs);
exponents.ComputeExponents(data);
return EstimateNumeric(data.avgs);
}
double Errors::EstimateNumeric(const MeanMap& avgs) {
auto _levels = avgs.GetLevelVector();
auto _avgs = avgs.GetYVector();
double alpha = exponents.alpha;
double Errors::EstimateNumeric(const MeanMap &avgs) {
auto _levels = avgs.GetLevelVector();
auto _avgs = avgs.GetYVector();
double alpha = exponents.alpha;
if (alpha == 0.0)
alpha = exponents.ComputeAlpha(avgs);
if (alpha == 0.0) alpha = exponents.ComputeAlpha(avgs);
double err = 0.0;
int integer = avgs.size() - 2;
for (int i = 1; i < _levels.size(); i++) {
err = max(_avgs[i] / (pow(2.0, alpha) - 1) / (pow(2.0, integer)), err);
integer--;
}
double err = 0.0;
int integer = avgs.size() - 2;
for (int i = 1; i < _levels.size(); i++) {
err = max(_avgs[i] / (pow(2.0, alpha) - 1) / (pow(2.0, integer)), err);
integer--;
}
return err;
return err;
}
double Errors::EstimateStochastic(const MultiLevelMonteCarloData &data) {
return EstimateStochastic(data.ctrs, data.vars);
return EstimateStochastic(data.ctrs, data.vars);
}
double Errors::EstimateStochastic(const SampleCounterMap &ctrs, const SVarMap &vars) {
auto _ctrs = ctrs.GetMVector();
auto _vars = vars.GetYVector();
auto _ctrs = ctrs.GetMVector();
auto _vars = vars.GetYVector();
double err = 0.0;
for (int i = 0; i < _ctrs.size(); i++)
err += _vars[i] / double(_ctrs[i]);
double err = 0.0;
for (int i = 0; i < _ctrs.size(); i++)
err += _vars[i] / double(_ctrs[i]);
return sqrt(err);
return sqrt(err);
}
......@@ -49,10 +49,10 @@ INSTANTIATE_TEST_SUITE_P(TestMainProgram, TestMainProgramElliptic, Values(
TEST_P(TestMainProgramElliptic, TestRun) {
EXPECT_EQ(mainProgram->Run(), 0);
EXPECT_TRUE(mainProgram->mlmc.errors.total < mainProgram->mlmc.epsilon);
EXPECT_TRUE(mainProgram->mlmc.exponents.alpha > 0);
EXPECT_TRUE(mainProgram->mlmc.exponents.beta > 0);
EXPECT_TRUE(mainProgram->mlmc.exponents.gamma > 0);
EXPECT_TRUE(mainProgram->estimator->TotalError() < mainProgram->estimator->Epsilon());
// EXPECT_TRUE(mainProgram->estimator.exponents.alpha > 0);
// EXPECT_TRUE(mainProgram->estimator.exponents.beta > 0);
// EXPECT_TRUE(mainProgram->estimator.exponents.gamma > 0);
}
//class TestMainProgramTransport : public TestMainProgram {
......
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