Merge branch 'use-notebooks-in-pipeline' into 'feature'

Use notebooks in pipeline

See merge request mpp/mlmc!19

.gitlab-ci.yml

@@ -60,17 +60,18 @@ elliptic-exercises-mlmc:
     RUN_EXPERIMENTS: ''
   image: ${REGISTRY}/${IMAGE_NAME_MLMC}
   timeout: '1h'
-  only:
-    variables: [ $RUN_EXPERIMENTS == 'all', $RUN_EXPERIMENTS == 'elliptic' ]
+  #only:
+  #  variables: [ $RUN_EXPERIMENTS == 'all', $RUN_EXPERIMENTS == 'elliptic' ]
   script:
     - cd /mpp/notebooks
     - ls -al
     - jupyter nbconvert --ExecutePreprocessor.timeout=2400
       --execute --to html Elliptic\ Experiments.ipynb
+    - cp Elliptic\ Experiments.html $CI_PROJECT_DIR/
   dependencies: [ "build-mlmc" ]
   tags: [ docker ]
   artifacts:
-    paths: [ notebooks/Elliptic\ Experiments.html ]
+    paths: [ Elliptic\ Experiments.html ]
 
 
 deploy-mlmc:

CMakeLists.txt

@@ -4,8 +4,8 @@ project(MLMC)
 set(SPACE_DIM 2 CACHE STRING "SPACE_DIM")
 set(USE_FFTW ON CACHE STRING "USE_FFTW")
 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/conf/mlmc_elliptic.conf

@@ -23,11 +23,10 @@ initSampleAmount = [12, 6, 3]
 
 # ----- Stochastic Field -----
 StochasticField = LogNormal
-mean = 0.0
+Mean = 0.0
 sigma = 1.0
 lambda = [0.15, 0.15]
 smoothing = 1.0
-evtol = 1e-10
 
 # ----- Plotting -----
 GeneratorPlotting = 0
@@ -38,8 +37,9 @@ MCVerbose = 1
 MLMCVerbose = 1
 MainVerbose = 1
 MeshVerbose = 1
-ConfigVerbose = 0
+ConfigVerbose = 1
 LinearVerbose = 0
 NewtonVerbose = 0
-PDESolverVerbose = 1
-GeneratorVerbose = 0
\ No newline at end of file
+AssembleVerbose = 0
+PDESolverVerbose = 0
+GeneratorVerbose = 0

mlmc/src/Main.hpp

@@ -6,28 +6,29 @@
 
 class MainProgram {
 private:
-    int verbose = 1;
+  int verbose = 1;
 
 public:
-    MultilevelMonteCarlo mlmc;
-
-    MainProgram() : mlmc(MultilevelMonteCarlo()){
-        config.get("MainVerbose", verbose);
-    }
-
-    int Run() {
-        mout.StartBlock("MLMC Experiment");
-        mout << "Run" << endl;
-        mlmc.Method();
-        mout.EndBlock();
-        mout << endl;
-
-        mlmc.PrintMCResults();
-        mlmc.PrintExponentResults();
-        mlmc.PrintMLMCResults();
-
-        return 0;
-    }
+  MultilevelMonteCarlo mlmc;
+
+  MainProgram() : mlmc(MultilevelMonteCarlo()) {
+    config.get("MainVerbose", verbose);
+  }
+
+  int Run() {
+    config.PrintInfo();
+    mout.StartBlock("MLMC Experiment");
+    mout << "Run" << endl;
+    mlmc.Method();
+    mout.EndBlock();
+    mout << endl;
+
+    mlmc.PrintMCResults();
+    mlmc.PrintExponentResults();
+    mlmc.PrintMLMCResults();
+
+    return 0;
+  }
 };
 
 #endif //MAIN_HPP

mlmc/src/basics/Level.hpp

@@ -62,6 +62,7 @@ struct Level {
  * * make second key template for degree
  *   to realize MultidegreeMonteCarlo
  * * also key (thus Level) should include proc information
+ * && ?
  */
 
 template<typename T>

mlmc/src/basics/Sample.hpp

@@ -20,17 +20,21 @@ struct SampleID {
     SampleID() {}
 
     SampleID(Level level, int number, bool coarse, const std::string &name = "U") :
-        level(level), number(number), coarse(coarse) {}
+        level(level), number(number), coarse(coarse), name(name) {}
 
     SampleID(int level, int number, bool coarse, const std::string &name = "U") :
-        level(Level(level)), number(number), coarse(coarse) {}
+        level(Level(level)), number(number), coarse(coarse), name(name) {}
 
     int l() const { return level.ActualLevel(coarse); };
 
+    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 {
-        return name + "." + std::to_string(level.fine)
-            + "." + std::to_string((int) coarse)
-            + "." + std::to_string(number);
+        return IdString(name);
     }
 };
 

mlmc/src/generators/SampleGenerator.hpp

@@ -14,154 +14,158 @@
 typedef std::complex<double> Complex;
 
 static double abs_compl(const Complex &complex) {
-    return sqrt(pow(complex.real(), 2) + pow(complex.imag(), 2));
+  return sqrt(pow(complex.real(), 2) + pow(complex.imag(), 2));
 }
 
 template<typename T>
 class SampleGenerator {
 protected:
-    int plotting = 0;
+  int plotting = 0;
 
-    int verbose = 0;
+  int verbose = 0;
 
-    virtual void drawSample(const SampleID &id) = 0;
+  virtual void drawSample(const SampleID &id) = 0;
+
+  virtual void plotSample(const SampleID &id) {};
 
 public:
-    Meshes &meshes;
+  Meshes &meshes;
 
-    explicit SampleGenerator(Meshes &meshes) : meshes(meshes) {
-        config.get("GeneratorPlotting", plotting);
-        config.get("GeneratorVerbose", verbose);
-    }
+  explicit SampleGenerator(Meshes &meshes) : meshes(meshes) {
+    // Todo add prefix
+    config.get("GeneratorPlotting", plotting);
+    config.get("GeneratorVerbose", verbose);
+  }
 
-    virtual ~SampleGenerator() {}
+  virtual ~SampleGenerator() {}
 
-    void DrawSample(const SampleID &id) {
-        mout.StartBlock(Name());
-        vout(1) << id.IdString() << endl;
-        drawSample(id);
-        mout.EndBlock(verbose == 0);
-    }
+  void DrawSample(const SampleID &id) {
+    mout.StartBlock(Name());
+    vout(1) << id.IdString() << endl;
+    drawSample(id);
+    if (plotting) plotSample(id);
+    mout.EndBlock(verbose == 0);
+  }
 
-    virtual string Name() const = 0;
+  virtual string Name() const = 0;
 
-    virtual size_t SampleSize() const { return 0; };
+  virtual size_t SampleSize() const { return 0; };
 
-    virtual T EvalSample() {
-        Exit("Not implemented")
-    }
+  virtual T EvalSample() {
+    Exit("Not implemented")
+  }
 
-    virtual T EvalSample(const Point &) {
-        Exit("Not implemented")
-    }
+  virtual T EvalSample(const Point &) {
+    Exit("Not implemented")
+  }
 
-    virtual T EvalSample(int, const Point &) {
-        Exit("Not implemented")
-    }
+  virtual T EvalSample(int, const Point &) {
+    Exit("Not implemented")
+  }
 
-    virtual T EvalSample(double, const Point &) {
-        Exit("Not implemented")
-    }
+  virtual T EvalSample(double, const Point &) {
+    Exit("Not implemented")
+  }
 
-    virtual T EvalSample(const cell &) {
-        Exit("Not implemented")
-    }
+  virtual T EvalSample(const cell &) {
+    Exit("Not implemented")
+  }
 
-    virtual T EvalSample(int, const cell &) {
-        Exit("Not implemented")
-    }
+  virtual T EvalSample(int, const cell &) {
+    Exit("Not implemented")
+  }
 
-    virtual T EvalSample(double, const cell &) {
-        Exit("Not implemented")
-    }
+  virtual T EvalSample(double, const cell &) {
+    Exit("Not implemented")
+  }
 };
 
 class ScalarDummy : public SampleGenerator<Scalar> {
 public:
-    ScalarDummy(Meshes &meshes) : SampleGenerator(meshes) { verbose = 0; }
+  ScalarDummy(Meshes &meshes) : SampleGenerator(meshes) { verbose = 0; }
 
-    void drawSample(const SampleID &id) override {
-        vout(1) << "Dummy Generator" << endl;
-    };
+  void drawSample(const SampleID &id) override {
+    vout(1) << "Dummy Generator" << endl;
+  };
 
-    virtual string Name() const override { return "DummyGenerator"; };
+  virtual string Name() const override { return "DummyGenerator"; };
 };
 
 class ComplexDummy : public SampleGenerator<Complex> {
 public:
-    ComplexDummy(Meshes &meshes) : SampleGenerator(meshes) { verbose = 0; }
+  ComplexDummy(Meshes &meshes) : SampleGenerator(meshes) { verbose = 0; }
 
-    void drawSample(const SampleID &id) override {
-        vout(1) << "Dummy Generator" << endl;
-    };
+  void drawSample(const SampleID &id) override {
+    vout(1) << "Dummy Generator" << endl;
+  };
 
-    virtual string Name() const override { return "DummyGenerator"; };
+  virtual string Name() const override { return "DummyGenerator"; };
 };
 
 class VectorFieldDummy : public SampleGenerator<VectorField> {
 public:
-    VectorFieldDummy(Meshes &meshes) : SampleGenerator(meshes) { verbose = 0; }
+  VectorFieldDummy(Meshes &meshes) : SampleGenerator(meshes) { verbose = 0; }
 
-    void drawSample(const SampleID &id) override {
-        vout(1) << "Dummy Generator" << endl;
-    };
+  void drawSample(const SampleID &id) override {
+    vout(1) << "Dummy Generator" << endl;
+  };
 
-    virtual string Name() const override { return "DummyGenerator"; };
+  virtual string Name() const override { return "DummyGenerator"; };
 };
 
 class TensorDummy : public SampleGenerator<Tensor> {
 public:
-    TensorDummy(Meshes &meshes) : SampleGenerator(meshes) { verbose = 0; }
+  TensorDummy(Meshes &meshes) : SampleGenerator(meshes) { verbose = 0; }
 
-    void drawSample(const SampleID &id) override {
-        vout(1) << "Dummy Generator" << endl;
-    };
+  void drawSample(const SampleID &id) override {
+    vout(1) << "Dummy Generator" << endl;
+  };
 
-    virtual string Name() const override { return "DummyGenerator"; };
+  virtual string Name() const override { return "DummyGenerator"; };
 };
 
 class RVectorDummy : public SampleGenerator<RVector> {
 public:
-    RVectorDummy(Meshes &meshes) : SampleGenerator(meshes) { verbose = 0; }
+  RVectorDummy(Meshes &meshes) : SampleGenerator(meshes) { verbose = 0; }
 
-    void drawSample(const SampleID &id) override {
-        vout(1) << "Dummy Generator" << endl;
-    };
+  void drawSample(const SampleID &id) override {
+    vout(1) << "Dummy Generator" << endl;
+  };
 
-    virtual string Name() const override { return "DummyGenerator"; };
+  virtual string Name() const override { return "DummyGenerator"; };
 };
 
 class CVectorDummy : public SampleGenerator<CVector> {
 public:
-    CVectorDummy(Meshes &meshes) : SampleGenerator(meshes) { verbose = 0; }
+  CVectorDummy(Meshes &meshes) : SampleGenerator(meshes) { verbose = 0; }
 
-    void drawSample(const SampleID &id) override {
-        vout(1) << "Dummy Generator" << endl;
-    };
+  void drawSample(const SampleID &id) override {
+    vout(1) << "Dummy Generator" << endl;
+  };
 
-    virtual string Name() const override { return "DummyGenerator"; };
+  virtual string Name() const override { return "DummyGenerator"; };
 };
 
 class RMatrixDummy : public SampleGenerator<RMatrix> {
 public:
-    RMatrixDummy(Meshes &meshes) : SampleGenerator(meshes) { verbose = 0; }
+  RMatrixDummy(Meshes &meshes) : SampleGenerator(meshes) { verbose = 0; }
 
-    void drawSample(const SampleID &id) override {
-        vout(1) << "Dummy Generator" << endl;
-    };
+  void drawSample(const SampleID &id) override {
+    vout(1) << "Dummy Generator" << endl;
+  };
 
-    virtual string Name() const override { return "DummyGenerator"; };
+  virtual string Name() const override { return "DummyGenerator"; };
 };
 
 class CMatrixDummy : public SampleGenerator<CMatrix> {
 public:
-    CMatrixDummy(Meshes &meshes) : SampleGenerator(meshes) { verbose = 0; }
+  CMatrixDummy(Meshes &meshes) : SampleGenerator(meshes) { verbose = 0; }
 
-    void drawSample(const SampleID &id) override {
-        vout(1) << "Dummy Generator" << endl;
-    };
+  void drawSample(const SampleID &id) override {
+    vout(1) << "Dummy Generator" << endl;
+  };
 
-    virtual string Name() const override { return "DummyGenerator"; };
+  virtual string Name() const override { return "DummyGenerator"; };
 };
 
 #endif //SAMPLEGENERATOR_HPP

mlmc/src/generators/SampleGeneratorContainer.cpp

@@ -2,45 +2,49 @@
 
 
 void SampleGeneratorContainer::init(Meshes &meshes) {
-    scalarGenerator = new ScalarDummy(meshes);
-    complexGenerator = new ComplexDummy(meshes);
-    vectorFieldGenerator = new VectorFieldDummy(meshes);
-    tensorGenerator = new TensorDummy(meshes);
-    rVectorGenerator = new RVectorDummy(meshes);
-    cVectorGenerator = new CVectorDummy(meshes);
-    rMatrixGenerator = new RMatrixDummy(meshes);
-    cMatrixGenerator = new CMatrixDummy(meshes);
+  scalarGenerator = new ScalarDummy(meshes);
+  complexGenerator = new ComplexDummy(meshes);
+  vectorFieldGenerator = new VectorFieldDummy(meshes);
+  tensorGenerator = new TensorDummy(meshes);
+  rVectorGenerator = new RVectorDummy(meshes);
+  cVectorGenerator = new CVectorDummy(meshes);
+  rMatrixGenerator = new RMatrixDummy(meshes);
+  cMatrixGenerator = new CMatrixDummy(meshes);
 }
 
 void SampleGeneratorContainer::initWithGenNames(Meshes &meshes, GeneratorNames names) {
-    for (auto const &genName: names) {
-        if (genName == "NormalDistributionReal")
-            scalarGenerator = new NormalDistributionReal(meshes);
-        if (genName == "NormalDistributionComplex")
-            complexGenerator = new NormalDistributionComplex(meshes);
-        if (genName == "NormalDistributionRVector")
-            rVectorGenerator = new NormalDistributionRVector(meshes);
-        if (genName == "NormalDistributionCVector")
-            cVectorGenerator = new NormalDistributionCVector(meshes);
-        if (genName == "NormalDistributionRMatrix")
-            rMatrixGenerator = new NormalDistributionRMatrix(meshes);
-        if (genName == "NormalDistributionCMatrix")
-            cMatrixGenerator = new NormalDistributionCMatrix(meshes);
+  for (auto const &genName: names) {
+    if (genName == "NormalDistributionReal")
+      scalarGenerator = new NormalDistributionReal(meshes);
+    if (genName == "NormalDistributionComplex")
+      complexGenerator = new NormalDistributionComplex(meshes);
+    if (genName == "NormalDistributionRVector")
+      rVectorGenerator = new NormalDistributionRVector(meshes);
+    if (genName == "NormalDistributionCVector")
+      cVectorGenerator = new NormalDistributionCVector(meshes);
+    if (genName == "NormalDistributionRMatrix")
+      rMatrixGenerator = new NormalDistributionRMatrix(meshes);
+    if (genName == "NormalDistributionCMatrix")
+      cMatrixGenerator = new NormalDistributionCMatrix(meshes);
 
-        if (genName == "CirculantEmbedding")
-            tensorGenerator = new MultilevelCirculantEmbedding(meshes);
-        if (genName == "HybridFluxGenerator")
-            new HybridFaceNormalFluxGenerator(meshes);
+    if (genName == "CirculantEmbedding") {
+      if (meshes.dim() == 1)
+        tensorGenerator = new CirculantEmbedding1D(meshes);
+      else
+        tensorGenerator = new CirculantEmbedding2D(meshes);
     }
+    if (genName == "HybridFluxGenerator")
+      new HybridFaceNormalFluxGenerator(meshes);
+  }
 }
 
 void SampleGeneratorContainer::DrawSample(const SampleID &id) {
-    scalarGenerator->DrawSample(id);
-    complexGenerator->DrawSample(id);
-    vectorFieldGenerator->DrawSample(id);
-    tensorGenerator->DrawSample(id);
-    rVectorGenerator->DrawSample(id);
-    cVectorGenerator->DrawSample(id);
-    rMatrixGenerator->DrawSample(id);
-    cMatrixGenerator->DrawSample(id);
+  scalarGenerator->DrawSample(id);
+  complexGenerator->DrawSample(id);
+  vectorFieldGenerator->DrawSample(id);
+  tensorGenerator->DrawSample(id);
+  rVectorGenerator->DrawSample(id);
+  cVectorGenerator->DrawSample(id);
+  rMatrixGenerator->DrawSample(id);
+  cMatrixGenerator->DrawSample(id);
 }

mlmc/src/generators/algorithms/CirculantEmbedding.cpp

 #include "CirculantEmbedding.hpp"
 
 
-void CirculantEmbedding1D::generateFineSample(const SampleID &id,
-                                              Vector *&fineSample,
-                                              Vector *&coarseSample) {
-  fineSample = new Vector(disc[id.level.fine]);
-  *fineSample = 0.0;
-  if (internalCounter == 0)
-    fineComplexField = GenerateField(id);
-
-  for (cell c = fineSample->cells(); c != fineSample->cells_end(); c++) {
-    int i = floor(c()[0] * fineSample->size());
-    if (fieldType == "Gaussian") {
-      if (internalCounter == 0)
-        fineSample->operator()(c(), 0) = fineComplexField[i].real() + mean;
-      else
-        fineSample->operator()(c(), 0) = fineComplexField[i].imag() + mean;
-    } else if (fieldType == "LogNormal") {
-      if (internalCounter == 0)
-        fineSample->operator()(c(), 0) = exp(fineComplexField[i].real() + mean);
-      else
-        fineSample->operator()(c(), 0) = exp(fineComplexField[i].imag() + mean);
-    } else Exit("Has to be Gaussian or LogNormal")
-  }
-  internalCounter += 1;
-  internalCounter %= 2;
-
-  if (coarseSample) delete coarseSample;
-  coarseSample = nullptr;
-}
-
-void CirculantEmbedding1D::generateCoarseSample(const SampleID &id,
-                                                Vector *&fineSample,
-                                                Vector *&coarseSample) {
-  coarseSample = new Vector(disc[id.level.coarse]);
-
-  (*coarseSample) = 0;
-//    for (cell r = coarseSample->rows(); r != coarseSample->rows_end(); r++) {
-  for (cell c = coarseSample->cells(); c != coarseSample->cells_end(); c++) {
-    row coarseRow = coarseSample->find_row(c());
-    for (int k = 0; k < c.Children(); ++k) {
-      row fineRow = fineSample->find_row(c.Child(k));
-      for (int i = 0; i < coarseRow.n(); ++i)
-        coarseSample->operator()(coarseRow, i) +=
-          fineSample->operator()(fineRow, i);
-    }
-    for (int i = 0; i < coarseRow.n(); ++i)
-      coarseSample->operator()(coarseRow, i) /= c.Children();
-  }
-
-  if (fineSample) delete fineSample;
-  fineSample = nullptr;
+RVector Downsample(const RVector &fSample) {
+  if (fSample.size() % 2 != 0) Exit(
+    "Vector should be of size 2 * size ")
+  RVector cSample(fSample.size() / 2);
+  for (int i = 0; i < cSample.size(); i++)
+    cSample[i] = (fSample[2 * i] + fSample[2 * i + 1]) / 2;
+  return cSample;
 }
 
-void CirculantEmbedding2D::generateFineSample(const SampleID &id,
-                                              Vector *&fineSample,
-                                              Vector *&coarseSample) {
-  fineSample = new Vector(disc[id.level.fine]);
-  *fineSample = 0.0;
-  if (internalCounter == 0)
-    fineComplexField = GenerateField(id);
-
-  for (cell c = fineSample->cells(); c != fineSample->cells_end(); c++) {
-    int i = floor(c()[0] * sqrt(fineSample->size()));
-    int j = floor(c()[1] * sqrt(fineSample->size()));
-    if (fieldType == "Gaussian") {
-      if (internalCounter == 0)
-        fineSample->operator()(c(), 0) = fineComplexField[i][j].real() + mean;
-      else
-        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);
-      else
-        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;
-}
-
-void CirculantEmbedding2D::generateCoarseSample(const SampleID &id,
-                                                Vector *&fineSample,
-                                                Vector *&coarseSample) {
-  coarseSample = new Vector(disc[id.level.coarse]);
-  *coarseSample = 0;
-  for (cell c = coarseSample->cells(); c != coarseSample->cells_end(); ++c) {
-    row coarseRow = coarseSample->find_row(c());
-    for (int k = 0; k < c.Children(); ++k) {
-      row fineRow = fineSample->find_row(c.Child(k));
-      for (int i = 0; i < coarseRow.n(); ++i)
-        coarseSample->operator()(coarseRow, i) +=
-          fineSample->operator()(fineRow, i);
-    }
-    for (int i = 0; i < coarseRow.n(); ++i)
-      coarseSample->operator()(coarseRow, i) /= c.Children();
-  }
-
-  if (fineSample) delete fineSample;
-  fineSample = nullptr;