multilevelcirculant embedding

mlmc/src/stochastics/CirculantEmbedding.cpp

 #include "CirculantEmbedding.hpp"
 
 
-void CirculantEmbedding1D::drawSample(SampleID id,
-                                      Vector *&fineSample,
-                                      Vector *&coarseSample) {
-    if (internalCounter == 0) fineComplexField = generateField();
-
-}
-
 void CirculantEmbedding1D::generateFineSample(SampleID id,
                                               Vector *&fineSample,
                                               Vector *&coarseSample) {
@@ -35,16 +28,12 @@ void CirculantEmbedding1D::generateFineSample(SampleID id,
 
     if (coarseSample) delete coarseSample;
     coarseSample = nullptr;
-
-    if (plotting)
-        plotter->PlotVector("kappa", *fineSample,
-                            1, l, "CellData");
 }
 
 void CirculantEmbedding1D::generateCoarseSample(SampleID id,
                                                 Vector *&fineSample,
                                                 Vector *&coarseSample) {
-    coarseSample = new Vector(cellMGraphs[l - this->meshes.pLevel() - 1]);
+    coarseSample = new Vector(cellMGraphs[id.level - this->meshes.pLevel() - 1]);
 
     (*coarseSample) = 0;
     for (cell c = coarseSample->cells(); c != coarseSample->cells_end(); c++) {
@@ -61,10 +50,6 @@ void CirculantEmbedding1D::generateCoarseSample(SampleID id,
 
     if (fineSample) delete fineSample;
     fineSample = nullptr;
-
-    if (plotting)
-        plotter->PlotVector("kappa", *coarseSample, 1,
-                            l - 1, "CellData");
 }
 
 ComplexField1D CirculantEmbedding1D::generateField() {
@@ -97,9 +82,6 @@ void CirculantEmbedding1D::createCovarianceMatrix(ToeplitzRow &ar, ToeplitzColum
     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;
@@ -205,16 +187,12 @@ void CirculantEmbedding2D::generateFineSample(SampleID id,
 
     if (coarseSample) delete coarseSample;
     coarseSample = nullptr;
-
-    if (plotting)
-        plotter->PlotVector("kappa", *fineSample,
-                            1, l, "CellData");
 }
 
 void CirculantEmbedding2D::generateCoarseSample(SampleID id,
                                                 Vector *&fineSample,
                                                 Vector *&coarseSample) {
-    coarseSample = new Vector(cellMGraphs[l - this->meshes.pLevel() - 1]);
+    coarseSample = new Vector(cellMGraphs[id.level - this->meshes.pLevel() - 1]);
 
     *coarseSample = 0;
     for (cell c = coarseSample->cells(); c != coarseSample->cells_end(); ++c) {
@@ -231,10 +209,6 @@ void CirculantEmbedding2D::generateCoarseSample(SampleID id,
 
     if (fineSample) delete fineSample;
     fineSample = nullptr;
-
-    if (plotting)
-        plotter->PlotVector("kappa", *coarseSample,
-                            1, l - 1, "CellData");
 }
 
 ComplexField2D CirculantEmbedding2D::generateField() {

mlmc/src/stochastics/CirculantEmbedding.hpp

@@ -34,7 +34,7 @@ typedef std::vector<std::vector<double>> BlockToeplitzColumn;
 
 typedef std::vector<std::vector<double>> BlockCirculantRow;
 
-class CirculantEmbedding : public SampleGenerator<Vector> {
+class CirculantEmbedding {
 protected:
     int internalCounter = 0;
 
@@ -48,10 +48,11 @@ protected:
 
     MatrixGraphs cellMGraphs;
 
+    Meshes &meshes;
+
 public:
-    explicit CirculantEmbedding(int l, Meshes &meshes) :
-        SampleGenerator(l, meshes),
-        cellMGraphs(meshes, dof(new CellDoF(1))) {
+    explicit CirculantEmbedding(Meshes &meshes) :
+        meshes(meshes), cellMGraphs(meshes, dof(new CellDoF(1))) {
 
         config.get("evtol", evtol);
         config.get("StochasticField", fieldType);
@@ -61,6 +62,14 @@ public:
         rnmg = new RNManager(streamnum, nstreams, gtype);
     }
 
+    virtual void generateFineSample(SampleID id,
+                                    Vector *&fineSample,
+                                    Vector *&coarseSample) = 0;
+
+    virtual void generateCoarseSample(SampleID id,
+                                      Vector *&fineSample,
+                                      Vector *&coarseSample) = 0;
+
     virtual ~CirculantEmbedding() {
         delete rnmg;
     };
@@ -79,7 +88,7 @@ public:
     double domain_end = 1;
 
     explicit CirculantEmbedding1D(int l, Meshes &meshes) :
-        CirculantEmbedding(l, meshes) {
+        CirculantEmbedding(meshes) {
         numberOfCells = meshes[l].CellCount();
         sqrtEigenvalues = computeSqrtEV();
         numberOfCellsEmbedded = sqrtEigenvalues.size();
@@ -93,10 +102,6 @@ public:
                               Vector *&fineSample,
                               Vector *&coarseSample) override;
 
-    void drawSample(SampleID id,
-                    Vector *&fineSample,
-                    Vector *&coarseSample) override;
-
 //private:
     ComplexField1D generateField();
 
@@ -131,7 +136,7 @@ public:
     double domain2_end = 1;
 
     explicit CirculantEmbedding2D(int l, Meshes &meshes) :
-        CirculantEmbedding(l, meshes) {
+        CirculantEmbedding(meshes) {
         numberOfCells[0] = (int) pow(2, l);
         numberOfCells[1] = (int) pow(2, l);
         sqrtEigenvalues = getSqrtEV();
@@ -166,4 +171,54 @@ public:
     SqrtEigenValues2D getSqrtEV();
 };
 
+class MultilevelCirculantEmbedding : public SampleGenerator {
+private:
+    std::map<int, CirculantEmbedding *> circulantEmbeddings;
+
+    Vector *fineSample = nullptr;
+
+    Vector *coarseSample = nullptr;
+
+protected:
+    void drawSample(SampleID id) override {
+        if (!id.coarse) {
+            circulantEmbeddings[id.level]->generateFineSample(id, fineSample,
+                                                              coarseSample);
+            if (plotting)
+                plotter->PlotVector("kappa", *fineSample, 1,
+                                    id.level - 1, "CellData");
+        } else {
+            circulantEmbeddings[id.level]->generateCoarseSample(id, fineSample,
+                                                                coarseSample);
+            if (plotting)
+                plotter->PlotVector("kappa", *coarseSample, 1,
+                                    id.level, "CellData");
+        }
+    }
+
+public:
+    MultilevelCirculantEmbedding(Meshes &meshes) : SampleGenerator(meshes) {
+        for (int l = meshes.pLevel(); l <= meshes.Level(); l++) {
+            if (meshes.dim() == 1)
+                circulantEmbeddings[l] = new CirculantEmbedding1D(l, meshes);
+            if (meshes.dim() == 2)
+                circulantEmbeddings[l] = new CirculantEmbedding2D(l, meshes);
+        }
+    }
+
+    ~MultilevelCirculantEmbedding() {
+        for (auto &circEmb : circulantEmbeddings)
+            delete circEmb.second;
+        circulantEmbeddings.clear();
+    }
+
+    Tensor EvalTensorSample(const cell &c) override {
+        if (fineSample && !coarseSample) return fineSample->operator()(c(), 0) * One;
+        if (coarseSample && !fineSample) return coarseSample->operator()(c(), 0) * One;
+        Exit("Error in Generator")
+    }
+
+    string Name() const override { return "Circulant Embedding"; }
+};
+
 #endif //M_CIRCULANTEMBEDDING_H
\ No newline at end of file