Skip to content
GitLab
Commit d914eb5f authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

refactoring CirculantEmbedding

parent a1f57880
Hide whitespace changes
Inline Side-by-side
mlmc/src/stochastics/CirculantEmbedding.cpp
View file @ d914eb5f
#include "CirculantEmbedding.hpp"
#include "main/MultilevelPlotter.hpp"
#include "dof/BasicDoFs.hpp"
void CirculantEmbedding1D::drawSample(SampleID id,
Vector *&fineSample,
Vector *&coarseSample) {
if (internalCounter == 0) fineComplexField = generateField();
}
void CirculantEmbedding1D::generateFineSample(SampleID id,
Vector *&fineSample,
......@@ -11,7 +17,7 @@ void CirculantEmbedding1D::generateFineSample(SampleID id,
fineComplexField = generateField();
for (cell c = fineSample->cells(); c != fineSample->cells_end(); c++) {
int i = floor(c()[0] * sqrt(fineSample->size()));
int i = floor(c()[0] * fineSample->size());
if (fieldType == "Gaussian") {
if (internalCounter == 0)
fineSample->operator()(c(), 0) = fineComplexField[i].real() + mean;
......@@ -28,6 +34,7 @@ void CirculantEmbedding1D::generateFineSample(SampleID id,
internalCounter %= 2;
if (coarseSample) delete coarseSample;
coarseSample = nullptr;
if (plotting)
plotter->PlotVector("kappa", *fineSample,
......@@ -53,6 +60,7 @@ void CirculantEmbedding1D::generateCoarseSample(SampleID id,
}
if (fineSample) delete fineSample;
fineSample = nullptr;
if (plotting)
plotter->PlotVector("kappa", *coarseSample, 1,
......@@ -83,13 +91,15 @@ ComplexField1D CirculantEmbedding1D::generateField() {
return Xf;
}
void CirculantEmbedding1D::createCovarianceMatrix(
vector<double> &ar, vector<double> &ac) {
void CirculantEmbedding1D::createCovarianceMatrix(ToeplitzRow &ar, ToeplitzColumn &ac) {
const vector<double> &coordinate1 = linspace(domain_start, domain_end, numberOfCells);
ar.resize(numberOfCells);
ac.resize(numberOfCells);
double y_rows = coordinate1.front();
for (cell c = meshes[l].cells(); c != meshes[l].cells_end(); ++c) {
}
for (int index_c = 0; index_c < numberOfCells; index_c++) {
double c1 = coordinate1[index_c];
double x_rows = c1;
......@@ -104,18 +114,21 @@ void CirculantEmbedding1D::createCovarianceMatrix(
}
}
void CirculantEmbedding1D::embedCovarianceMatrix(vector<double> &ar, vector<double> &ac,
vector<double> &br) {
void CirculantEmbedding1D::embedCovarianceMatrix(const ToeplitzRow &ar,
const ToeplitzColumn &ac,
CirculantRow &br) {
br.reserve(2 * ar.size() - 1);
br.insert(br.end(), ar.begin(), ar.end());
br.insert(br.end(), ac.rbegin(), ac.rend() - 1);
}
void CirculantEmbedding1D::computeEV(vector<double> &br, vector<double> &ev) {
void CirculantEmbedding1D::computeEV(const CirculantRow &br,
EigenValues1D &ev) {
fft(br, ev);
}
void CirculantEmbedding1D::padding(vector<double> &ar, vector<double> &ac) {
void CirculantEmbedding1D::padding(ToeplitzRow &ar,
ToeplitzColumn &ac) {
double delta = (domain_end - domain_start) / (numberOfCells - 1);
// Add 50% more points
......@@ -143,7 +156,7 @@ void CirculantEmbedding1D::padding(vector<double> &ar, vector<double> &ac) {
}
}
vector<double> CirculantEmbedding1D::getSqrtEV() {
SqrtEigenValues1D CirculantEmbedding1D::computeSqrtEV() {
vector<double> sqrt_ev_return, ar, ac, br, ev;
createCovarianceMatrix(ar, ac);
......@@ -165,6 +178,8 @@ vector<double> CirculantEmbedding1D::getSqrtEV() {
void CirculantEmbedding2D::generateFineSample(SampleID id,
Vector *&fineSample,
Vector *&coarseSample) {
fineSample = new Vector(cellMGraphs[id.level - this->meshes.pLevel()]);
if (internalCounter == 0)
fineComplexField = generateField();
......@@ -178,14 +193,19 @@ void CirculantEmbedding2D::generateFineSample(SampleID id,
fineSample->operator()(c(), 0) = fineComplexField[i][j].imag() + mean;
} else if (fieldType == "LogNormal") {
if (internalCounter == 0)
fineSample->operator()(c(), 0) = exp(fineComplexField[i][j].real() + mean);
fineSample->operator()(c(), 0) =
exp(fineComplexField[i][j].real() + mean);
else
fineSample->operator()(c(), 0) = exp(fineComplexField[i][j].imag() + mean);
fineSample->operator()(c(), 0) =
exp(fineComplexField[i][j].imag() + mean);
} else Exit("Has to be Gaussian or LogNormal")
}
internalCounter += 1;
internalCounter %= 2;
if (coarseSample) delete coarseSample;
coarseSample = nullptr;
if (plotting)
plotter->PlotVector("kappa", *fineSample,
1, l, "CellData");
......@@ -194,6 +214,8 @@ void CirculantEmbedding2D::generateFineSample(SampleID id,
void CirculantEmbedding2D::generateCoarseSample(SampleID id,
Vector *&fineSample,
Vector *&coarseSample) {
coarseSample = new Vector(cellMGraphs[l - this->meshes.pLevel() - 1]);
*coarseSample = 0;
for (cell c = coarseSample->cells(); c != coarseSample->cells_end(); ++c) {
row coarseRow = coarseSample->find_row(c());
......@@ -207,6 +229,9 @@ void CirculantEmbedding2D::generateCoarseSample(SampleID id,
coarseSample->operator()(coarseRow, i) /= c.Children();
}
if (fineSample) delete fineSample;
fineSample = nullptr;
if (plotting)
plotter->PlotVector("kappa", *coarseSample,
1, l - 1, "CellData");
......@@ -276,9 +301,9 @@ void CirculantEmbedding2D::createCovarianceMatrix(vector<vector<double>> &ar,
}
}
void CirculantEmbedding2D::embedCovarianceMatrix(vector<vector<double>> &ar,
vector<vector<double>> &ac,
vector<vector<double>> &br) {
void CirculantEmbedding2D::embedCovarianceMatrix(const BlockToeplitzRow &ar,
const BlockToeplitzColumn &ac,
BlockCirculantRow &br) {
int br_size = (int) (2 * ar.size() - 1);
br.resize(br_size);
for (int i = 0; i < br.size(); i++) {
......@@ -292,13 +317,13 @@ void CirculantEmbedding2D::embedCovarianceMatrix(vector<vector<double>> &ar,
}
}
void CirculantEmbedding2D::computeEV(vector<vector<double>> &br,
vector<vector<double>> &ev) {
void CirculantEmbedding2D::computeEV(const BlockCirculantRow &br,
EigenValues2D &ev) {
fft2(br, ev);
}
void
CirculantEmbedding2D::padding(vector<vector<double>> &ar, vector<vector<double>> &ac) {
void CirculantEmbedding2D::padding(BlockToeplitzRow &ar,
BlockToeplitzColumn &ac) {
double delta1 = (domain1_end - domain1_start) / (numberOfCells[0] - 1);
double delta2 = (domain2_end - domain2_start) / (numberOfCells[1] - 1);
......
mlmc/src/stochastics/CirculantEmbedding.hpp
View file @ d914eb5f
......@@ -7,27 +7,52 @@
#include "Algebra.h"
#include "SampleGenerator.hpp"
#include "dof/BasicDoFs.hpp"
#include "main/MultilevelPlotter.hpp"
typedef std::vector<std::complex<double>> ComplexField1D;
typedef std::vector<double> EigenValues1D;
typedef std::vector<double> SqrtEigenValues1D;
typedef std::vector<double> ToeplitzRow;
typedef std::vector<double> ToeplitzColumn;
typedef std::vector<double> CirculantRow;
typedef std::vector<std::vector<std::complex<double>>> ComplexField2D;
typedef std::vector<std::vector<double>> EigenValues2D;
typedef std::vector<std::vector<double>> SqrtEigenValues2D;
typedef std::vector<std::vector<double>> BlockToeplitzRow;
typedef std::vector<std::vector<double>> BlockToeplitzColumn;
typedef std::vector<std::vector<double>> BlockCirculantRow;
class CirculantEmbedding : public SampleGenerator<Vector> {
protected:
int internalCounter = 0;
RNManager *rnmg;
public:
double evtol = 1e-10;
string fieldType = "LogNormal";
double mean = 0.0;
MatrixGraphs cellMGraphs;
public:
explicit CirculantEmbedding(int l, Meshes &meshes) :
SampleGenerator(l, meshes),
cellMGraphs(meshes, dof(new CellDoF(1))) {
config.get("evtol", evtol);
config.get("StochasticField", fieldType);
config.get("mean", mean);
......@@ -36,13 +61,17 @@ public:
rnmg = new RNManager(streamnum, nstreams, gtype);
}
virtual ~CirculantEmbedding() { delete rnmg; };
virtual ~CirculantEmbedding() {
delete rnmg;
};
};
class CirculantEmbedding1D : public CirculantEmbedding, public CovarianceFunction1D {
public:
ComplexField1D fineComplexField;
vector<double> sqrtEigenvalues;
SqrtEigenValues1D sqrtEigenvalues;
int numberOfCellsEmbedded = 0;
int numberOfCells = 0;
......@@ -51,8 +80,8 @@ public:
explicit CirculantEmbedding1D(int l, Meshes &meshes) :
CirculantEmbedding(l, meshes) {
numberOfCells = (int) pow(2, l);
sqrtEigenvalues = getSqrtEV();
numberOfCells = meshes[l].CellCount();
sqrtEigenvalues = computeSqrtEV();
numberOfCellsEmbedded = sqrtEigenvalues.size();
}
......@@ -64,24 +93,43 @@ public:
Vector *&fineSample,
Vector *&coarseSample) override;
private:
void drawSample(SampleID id,
Vector *&fineSample,
Vector *&coarseSample) override;
//private:
ComplexField1D generateField();
void createCovarianceMatrix(vector<double> &ar, vector<double> &ac);
void createCovarianceMatrix(ToeplitzRow &ar,
ToeplitzColumn &ac);
static void embedCovarianceMatrix(vector<double> &ar,
vector<double> &ac,
vector<double> &br);
static void embedCovarianceMatrix(const ToeplitzRow &ar,
const ToeplitzColumn &ac,
CirculantRow &br);
static void computeEV(vector<double> &br, vector<double> &ev);
static void computeEV(const CirculantRow &br,
EigenValues1D &ev);
void padding(vector<double> &ar, vector<double> &ac);
void padding(ToeplitzRow &ar,
ToeplitzColumn &ac);
vector<double> getSqrtEV();
SqrtEigenValues1D computeSqrtEV();
};
class CirculantEmbedding2D : public CirculantEmbedding, public CovarianceFunction2D {
public:
ComplexField2D fineComplexField;
SqrtEigenValues2D sqrtEigenvalues;
int numberOfCellsEmbedded[2] = {0, 0};
int numberOfCells[2] = {0, 0};
double domain1_start = 0;
double domain1_end = 1;
double domain2_start = 0;
double domain2_end = 1;
explicit CirculantEmbedding2D(int l, Meshes &meshes) :
CirculantEmbedding(l, meshes) {
numberOfCells[0] = (int) pow(2, l);
......@@ -99,30 +147,23 @@ public:
Vector *&fineSample,
Vector *&coarseSample) override;
private:
ComplexField2D fineComplexField;
vector<vector<double>> sqrtEigenvalues;
int numberOfCellsEmbedded[2] = {0, 0};
int numberOfCells[2] = {0, 0};
double domain1_start = 0;
double domain1_end = 1;
double domain2_start = 0;
double domain2_end = 1;
//private:
ComplexField2D generateField();
void createCovarianceMatrix(vector<vector<double>> &ar, vector<vector<double>> &ac);
void createCovarianceMatrix(BlockToeplitzRow &ar,
BlockToeplitzColumn &ac);
static void embedCovarianceMatrix(vector<vector<double>> &ar,
vector<vector<double>> &ac,
vector<vector<double>> &br);
static void embedCovarianceMatrix(const BlockToeplitzRow &ar,
const BlockToeplitzColumn &ac,
BlockCirculantRow &br);
static void computeEV(vector<vector<double>> &br, vector<vector<double>> &ev);
static void computeEV(const BlockCirculantRow &br,
EigenValues2D &ev);
void padding(vector<vector<double>> &ar, vector<vector<double>> &ac);
void padding(BlockToeplitzRow &ar,
BlockToeplitzColumn &ac);
vector<vector<double>> getSqrtEV();
SqrtEigenValues2D getSqrtEV();
};
#endif //M_CIRCULANTEMBEDDING_H
\ No newline at end of file
Supports Markdown
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

Imprint