Skip to content

GitLab

Commit 6b552232 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

Merge branch '18-remove-level-and-cellcount-in-ce' into 'feature' 

Resolve "Remove Level and CellCount in CE"

Closes #18

See merge request mpp/mlmc!22
parents 505c0a23 59c25a34
Pipeline #148080 passed with stages
in 60 minutes and 59 seconds
Hide whitespace changes
Inline Side-by-side
mlmc/src/basics/Level.cpp
View file @ 6b552232
#include "Level.hpp"
void Level::Update(int level) {
fine = level;
coarse = level - 1;
if (coarse < 0) coarse = 0;
pLevel = coarse;
}
bool Level::operator==(const Level &level) const {
if (this->fine != level.fine)
return false;
if (this->coarse != level.coarse)
return false;
if (this->pLevel != level.pLevel)
return false;
return true;
}
bool Level::operator==(int level) const {
if (this->fine != level)
return false;
return true;
}
bool Level::operator!=(const Level &level) const {
return !(*this == level);
}
bool Level::operator!=(int level) const {
return !(*this == level);
}
bool Level::operator<=(const Level &level) const {
if (this->fine < level.fine) return true;
else return false;
}
bool Level::operator>=(const Level &level) const {
if (this->fine > level.fine) return true;
else return false;
}
bool Level::operator<(const Level &level) const {
if (this->fine < level.fine) return true;
else return false;
}
bool Level::operator>(const Level &level) const {
if (this->fine > level.fine) return true;
else return false;
}
bool Level::operator<(int level) const {
if (this->fine < level) return true;
else return false;
}
bool Level::operator>(int level) const {
if (this->fine > level) return true;
else return false;
}
bool Level::operator<=(int level) const {
if (this->fine <= level) return true;
else return false;
}
bool Level::operator>=(int level) const {
if (this->fine >= level) return true;
else return false;
}
mlmc/src/basics/Level.hpp
View file @ 6b552232
......@@ -6,75 +6,23 @@
#include <map>
/*
* Todo
* Only Fine true: pLevel = fine = coarse
* Only Fine false: plevel = coarse, fine = coarse + 1
*
* Add commSplit to Level
*/
struct Level {
int fine;
int coarse;
int pLevel;
Level() {}
explicit Level(int level) {
Update(level);
}
void Update(int level);
int ActualLevel(bool getCoarse) const {
if (!getCoarse) return fine;
else return coarse;
};
bool operator==(const Level &level) const;
bool operator==(int level) const;
bool operator!=(const Level &level) const;
bool operator!=(int level) const;
bool operator<=(const Level &level) const;
bool operator>=(const Level &level) const;
bool operator<(const Level &level) const;
bool operator>(const Level &level) const;
bool operator<(int level) const;
bool operator>(int level) const;
bool operator<=(int level) const;
bool operator>=(int level) const;
};
/*
* Todo:
* * make second key template for degree
* * idea: make second key template for degree
* to realize MultidegreeMonteCarlo
* * also key (thus Level) should include proc information
* && ?
*/
template<typename T>
struct LevelMap {
protected:
std::map<Level, T> _levelMap;
std::map<int, T> _levelMap;
public:
LevelMap() {};
LevelMap(std::initializer_list<std::pair<Level, T>> levelMap) {
for (std::pair<Level, T> pair : levelMap) {
LevelMap(std::initializer_list<std::pair<int, T>> levelMap) {
for (std::pair<int, T> pair : levelMap) {
_levelMap.emplace(pair.first, T(pair.second));
}
}
......@@ -82,7 +30,7 @@ public:
std::vector<int> GetLevelVector() const {
std::vector<int> levelVector;
for (auto &&[level, entry]: _levelMap)
levelVector.push_back(level.fine);
levelVector.push_back(level);
return levelVector;
}
......@@ -125,7 +73,7 @@ public:
auto clear() { _levelMap.clear(); }
auto insert(std::pair<Level, T> pair) { _levelMap.insert(pair); }
auto insert(std::pair<int, T> pair) { _levelMap.insert(pair); }
auto begin() { return _levelMap.begin(); }
......@@ -143,11 +91,11 @@ public:
auto rend() const { return _levelMap.end(); }
auto operator[](Level level) { return _levelMap[level]; }
auto operator[](int level) { return _levelMap[level]; }
auto at(Level level) { return _levelMap.at(level); }
auto at(int level) { return _levelMap.at(level); }
auto find(Level level) { return _levelMap.find(level); }
auto find(int level) { return _levelMap.find(level); }
friend Logging &operator<<(Logging &s, const LevelMap<T> &levelMap) {
s << "[";
......
mlmc/src/basics/Sample.hpp
View file @ 6b552232
......@@ -9,76 +9,80 @@
struct SampleID {
Level level;
public:
int fLevel;
std::string name;
int cLevel;
int number;
std::string name;
bool coarse; // Todo should be; SampleID makes more sense then
int number;
SampleID() {}
bool coarse; // Todo should be; SampleID makes more sense then
SampleID(Level level, int number, bool coarse, const std::string &name = "U") :
level(level), number(number), coarse(coarse), name(name) {}
SampleID() {}
SampleID(int level, int number, bool coarse, const std::string &name = "U") :
level(Level(level)), number(number), coarse(coarse), name(name) {}
SampleID(int level, int number, bool coarse, const std::string &name = "U") :
fLevel(level), cLevel( (level > 0) ? (level - 1) : 0),
number(number), coarse(coarse), name(name) {}
int l() const { return level.ActualLevel(coarse); };
int level() const {
if (coarse) return cLevel;
else return fLevel;
};
std::string IdString(const std::string &name_) const {
return name_ + "." + std::to_string(level.fine)
+ "." + std::to_string((int) coarse)
+ "." + std::to_string(number);
}
std::string IdString(const std::string &name_) const {
return name_ + "." + std::to_string(fLevel)
+ "." + std::to_string((int) coarse)
+ "." + std::to_string(number);
}
std::string IdString() const {
return IdString(name);
}
std::string IdString() const {
return IdString(name);
}
};
struct SampleSolution {
SampleID id;
SampleID id;
public:
double Q;
double Cost;
Vector U;
SampleSolution(IDiscretization *disc, const std::string &name = "U") :
U(Vector((*disc))) {
this->id.name = name;
Init();
}
SampleSolution(IDiscretization *disc, const SampleID &id,
const std::string &name = "U") :
id(id), U(Vector(*disc, id.l())) {
this->id.name = name;
Init();
}
SampleSolution(const Vector &U, const SampleID &id, const std::string &name = "U") :
id(id), U(U) {
this->id.name = name;
}
std::string IdString() const {
return id.IdString();
}
const SampleID& Id() const {
return id;
}
void Init() {
Q = 0.0;
Cost = 0.0;
U = 0.0;
}
double Q;
double Cost;
Vector U;
SampleSolution(IDiscretization *disc, const std::string &name = "U") :
U(Vector((*disc))) {
this->id.name = name;
Init();
}
SampleSolution(IDiscretization *disc, const SampleID &id,
const std::string &name = "U") :
id(id), U(Vector(*disc, id.level())) {
this->id.name = name;
Init();
}
SampleSolution(const Vector &U, const SampleID &id, const std::string &name = "U") :
id(id), U(U) {
this->id.name = name;
}
std::string IdString() const {
return id.IdString();
}
const SampleID &Id() const {
return id;
}
void Init() {
Q = 0.0;
Cost = 0.0;
U = 0.0;
}
};
#endif //SAMPLE_HPP
mlmc/src/generators/algorithms/CirculantEmbedding.cpp
View file @ 6b552232
......@@ -26,8 +26,8 @@ RMatrix Downsample(const RMatrix &fSample) {
*/
RVector CirculantEmbedding1D::ComputeSqrtEV() {
RVector toepRow(meshes.fine().CellCountGeometry());
RVector toepCol(meshes.fine().CellCountGeometry());
RVector toepRow(cellCount);
RVector toepCol(cellCount);
RVector circRow = covariance.EmbeddedToeplitzMatrix(toepRow, toepCol);
RVector eigenVal(circRow);
......@@ -59,12 +59,14 @@ RVector CirculantEmbedding1D::GenerateLogNormalField(const SampleID &id) {
if (internalCounter == 0)
fineComplexField = GenerateField(id);
RVector logNormalField(0.0, meshes.fine().CellCountGeometry());
RVector logNormalField(0.0, cellCount);
for (int i = 0; i < logNormalField.size(); i++)
if (internalCounter == 0)
logNormalField[i] = exp(fineComplexField[i].real() + mean);
else
logNormalField[i] = exp(fineComplexField[i].imag() + mean);
// internalCounter++;
// internalCounter = internalCounter % 2;
return logNormalField;
}
......@@ -72,12 +74,14 @@ RVector CirculantEmbedding1D::GenerateGaussianField(const SampleID &id) {
if (internalCounter == 0)
fineComplexField = GenerateField(id);
RVector gaussianField(0.0, meshes.fine().CellCountGeometry());
RVector gaussianField(0.0, cellCount);
for (int i = 0; i < gaussianField.size(); i++)
if (internalCounter == 0)
gaussianField[i] = fineComplexField[i].real() + mean;
else
gaussianField[i] = fineComplexField[i].imag() + mean;
// internalCounter++;
// internalCounter = internalCounter % 2;
return gaussianField;
}
......@@ -87,12 +91,15 @@ Tensor CirculantEmbedding1D::EvalSample(const cell &c) {
void CirculantEmbedding1D::initVec(const SampleID &id) {
if (sampleVec) delete sampleVec;
sampleVec = new Vector(disc, id.l());
sampleVec = new Vector(disc, id.level());
*sampleVec = 0.0;
}
void CirculantEmbedding1D::plotSample(const SampleID &id) {
mpp::plot(id.IdString("Kappa")) << *sampleVec << mpp::endp;
if (sampleVec != nullptr)
mpp::plot(id.IdString("Kappa")) << *sampleVec << mpp::endp;
else
Exit("Vector is nullptr")
}
void CirculantEmbedding1D::drawSample(const SampleID &id) {
......@@ -112,8 +119,8 @@ void CirculantEmbedding1D::drawSample(const SampleID &id) {
*/
RMatrix CirculantEmbedding2D::ComputeSqrtEV() {
RMatrix toepRows((int) sqrt(meshes.fine().CellCountGeometry()));
RMatrix toepCols((int) sqrt(meshes.fine().CellCountGeometry()));
RMatrix toepRows((int) sqrt(cellCount));
RMatrix toepCols((int) sqrt(cellCount));
RMatrix circRows = covariance.EmbeddedToeplitzMatrix(toepRows, toepCols);
RMatrix eigenVal(circRows);
......@@ -146,13 +153,15 @@ RMatrix CirculantEmbedding2D::GenerateLogNormalField(const SampleID &id) {
if (internalCounter == 0)
fineComplexField = GenerateField(id);
RMatrix logNormalField(0.0, (int) sqrt(meshes.fine().CellCountGeometry()));
RMatrix logNormalField(0.0, (int) sqrt(cellCount));
for (int i = 0; i < logNormalField.rows(); i++)
for (int j = 0; j < logNormalField.cols(); j++)
if (internalCounter == 0)
logNormalField[i][j] = exp(fineComplexField[i][j].real() + mean);
else
logNormalField[i][j] = exp(fineComplexField[i][j].imag() + mean);
// internalCounter++;
// internalCounter = internalCounter % 2;
return logNormalField;
}
......@@ -160,13 +169,15 @@ RMatrix CirculantEmbedding2D::GenerateGaussianField(const SampleID &id) {
if (internalCounter == 0)
fineComplexField = GenerateField(id);
RMatrix gaussianField(0.0, (int) sqrt(meshes.fine().CellCountGeometry()));
RMatrix gaussianField(0.0, (int) sqrt(cellCount));
for (int i = 0; i < gaussianField.rows(); i++)
for (int j = 0; j < gaussianField.cols(); j++)
if (internalCounter == 0)
gaussianField[i][j] = fineComplexField[i][j].real() + mean;
else
gaussianField[i][j] = fineComplexField[i][j].imag() + mean;
// internalCounter++;
// internalCounter = internalCounter % 2;
return gaussianField;
}
......@@ -176,7 +187,7 @@ Tensor CirculantEmbedding2D::EvalSample(const cell &c) {
void CirculantEmbedding2D::initVec(const SampleID &id) {
if (sampleVec) delete sampleVec;
sampleVec = new Vector(disc, id.l());
sampleVec = new Vector(disc, id.level());
*sampleVec = 0.0;
}
......
mlmc/src/generators/algorithms/CirculantEmbedding.hpp
View file @ 6b552232
......@@ -36,6 +36,8 @@ private:
RVector sample;
int cellCount;
protected:
void initVec(const SampleID &id);
......@@ -51,6 +53,7 @@ public:
config.get("StochasticField", fieldType);
config.get("Mean", mean);
cellCount = meshes.fine().CellCountGeometry();
sqrtEigenvalues = ComputeSqrtEV();
}
......@@ -91,6 +94,8 @@ private:
RMatrix sample;
int cellCount;
protected:
void initVec(const SampleID &id);
......@@ -107,6 +112,7 @@ public:
config.get("StochasticField", fieldType);
config.get("Mean", mean);
cellCount = meshes.fine().CellCountGeometry();
sqrtEigenvalues = ComputeSqrtEV();
}
......
mlmc/src/montecarlo/MonteCarlo.cpp
View file @ 6b552232
......@@ -2,7 +2,7 @@
void MonteCarlo::Method() {
mout.StartBlock("MonteCarlo l=" + to_string(level.fine));
mout.StartBlock("MonteCarlo l=" + to_string(level));
vout(1) << "Start with: " << aggregate;
method();
aggregate.UpdateParallel();
......
mlmc/src/montecarlo/MonteCarlo.hpp
View file @ 6b552232
......@@ -9,90 +9,92 @@
class MonteCarlo {
protected:
int verbose = 1;
int verbose = 1;
int plotting = 0;
int plotting = 0;
bool onlyFine = false;
bool onlyFine = false;
bool parallel = true;
bool parallel = true;
Level level;
int level;
Meshes *meshes;
Meshes *meshes;
PDESolver *pdeSolver;
PDESolver *pdeSolver;
public:
WelfordAggregate aggregate;
WelfordAggregate aggregate;
SampleCounter ctr; // todo remove, is still here for test
SampleCounter ctr; // todo remove, is still here for test
Sums sums;
Sums sums;
Averages avgs;
Averages avgs;
Variances vars;
Variances vars;
Kurtosis kurtosis;
Kurtosis kurtosis;
PDESolverCreator pdeSolverCreator;
PDESolverCreator pdeSolverCreator;
MeshesCreator meshesCreator;
MeshesCreator meshesCreator;
SampleID coarseId;
SampleID coarseId;
SampleID fineId;
SampleID fineId;
MonteCarlo(Level level, int dM, bool onlyFine) :
level(level),
onlyFine(onlyFine),
pdeSolverCreator(PDESolverCreator()),
meshesCreator(MeshesCreator(pdeSolverCreator.GetMeshName()).
WithCommSplit(aggregate.commSplit). // Todo: Is this correct?
WithPLevel(level.coarse).
WithLevel(level.fine)) {
Init(dM);
}
MonteCarlo(int level, int dM, bool onlyFine) :
level(level),
onlyFine(onlyFine),
pdeSolverCreator(PDESolverCreator()),
meshesCreator(MeshesCreator(pdeSolverCreator.GetMeshName()).
WithCommSplit(aggregate.commSplit). // Todo: Is this correct?
WithPLevel(level - 1).
WithLevel(level)) {
Init(dM);
}
MonteCarlo(Level level, int dM, bool onlyFine,
MeshesCreator meshesCreator, PDESolverCreator pdeCreator) :
level(level),
onlyFine(onlyFine),
pdeSolverCreator(pdeCreator),
meshesCreator(meshesCreator) {
Init(dM);
}
MonteCarlo(int level, int dM, bool onlyFine,
MeshesCreator meshesCreator, PDESolverCreator pdeCreator) :
level(level),
onlyFine(onlyFine),
pdeSolverCreator(pdeCreator),
meshesCreator(meshesCreator) {
Init(dM);
}
void Init(int dM) {
config.get("MCPlotting", plotting);
config.get("MCVerbose", verbose);
config.get("MCParallel", parallel);
void Init(int dM) {
config.get("MCPlotting", plotting);
config.get("MCVerbose", verbose);
config.get("MCParallel", parallel);
meshes = meshesCreator.Create();
pdeSolver = pdeSolverCreator.Create(*meshes);
meshes = meshesCreator.Create();
pdeSolver = pdeSolverCreator.Create(*meshes);
fineId.level = level;
fineId.coarse = false;
coarseId.level = level;
coarseId.coarse = true;
ctr.parallel = parallel;
ctr.UpdateSampleCounter(dM);
aggregate.parallel = parallel;
aggregate.UpdateSampleCounter(dM);
}
fineId.fLevel = level;
fineId.cLevel = level - 1;
fineId.coarse = false;
coarseId.fLevel = level;
coarseId.cLevel = level - 1;