Merge branch '34-derive-sprasegridgenerator-from-samplegenerator' into 'feature'

added weight solution - other possible approach could be to add the weight as...

Closes #34

See merge request !47

mluq/src/estimators/Estimator.cpp

@@ -7,7 +7,10 @@
 
 Estimator *EstimatorCreator::Create() {
   if (_estimatorName == "StochasticCollocation")
-    return new StochasticCollocation();
+    return new StochasticCollocation(
+      _epsilon, _level, _stochLevel, _onlyFine, _parallel,
+      _meshesCreator, _pdeSolverCreator
+    );
 
   if (_estimatorName == "MultilevelMonteCarlo")
     return new MultilevelEstimator(

mluq/src/estimators/Estimator.hpp

@@ -71,6 +71,8 @@ private:
 
   int _initSamples;
 
+  int _stochLevel;
+
   double _epsilon;
 
   Levels _levelVec;
@@ -120,6 +122,11 @@ public:
     return *this;
   }
 
+  EstimatorCreator WithStochLevel(int stochLevel) {
+    _stochLevel = stochLevel;
+    return *this;
+  }
+
   EstimatorCreator WithOnlyFine(bool onlyFine) {
     _onlyFine = onlyFine;
     return *this;

mluq/src/estimators/MonteCarlo.hpp

@@ -7,6 +7,8 @@
 
 #include "MeshesCreator.hpp"
 
+#include <utility>
+
 
 class MonteCarlo : public Estimator {
 protected:
@@ -34,7 +36,7 @@ public:
              PDESolverCreator pdeCreator = PDESolverCreator()) :
     Estimator(epsilon, level, initSamples), onlyFine(onlyFine), parallel(parallel),
     fineId(SampleID(level, 0, false)), coarseId(SampleID(level, 0, true)),
-    meshesCreator(meshesCreator), pdeSolverCreator(pdeCreator) {
+    meshesCreator(std::move(meshesCreator)), pdeSolverCreator(std::move(pdeCreator)) {
 
     config.get("MCVerbose", verbose);
     config.get("MCParallel", parallel);

mluq/src/estimators/StochasticCollocation.cpp

-//
-// Created by niklas on 11.05.21.
-//
-
 #include "StochasticCollocation.hpp"
+
+
+void StochasticCollocation::Method() {
+  mout.StartBlock("Stochastic Collocation l=" + to_string(level));
+  vout(1) << "Start with: " << aggregate;
+  method();
+  aggregate.UpdateParallel();
+  vout(1) << "End with: " << aggregate;
+  mout.EndBlock(verbose == 0);
+}
+
+void StochasticCollocation::method() {
+  std::unique_ptr<Meshes> meshes = meshesCreator.
+    WithMeshName(pdeSolverCreator.GetMeshName()).
+    WithCommSplit(aggregate.commSplit).
+    WithPLevel(level - 1).
+    WithLevel(level).
+    CreateUnique();
+
+  std::unique_ptr<PDESolver> pdeSolver = pdeSolverCreator.
+    CreateUnique(*meshes);
+
+  SampleSolution fineSolution(pdeSolver->GetDisc(), fineId);
+  SampleSolution coarseSolution(pdeSolver->GetDisc(), coarseId);
+
+  int numGridPoints = pdeSolver->GetProblem()->NumOfSamples();
+  aggregate.UpdateSampleCounter(numGridPoints);
+  while (aggregate.ctr.dMcomm != 0) {
+    updateSampleIds(fineSolution, coarseSolution);
+
+    pdeSolver->DrawSample(fineId);
+    pdeSolver->Run(fineSolution);
+
+    aggregate.Update(fineSolution, coarseSolution);
+  }
+}
+
+void StochasticCollocation::updateSampleIds(SampleSolution &fSol, SampleSolution &cSol) {
+  fineId.number = aggregate.index();
+  coarseId.number = aggregate.index();
+  fSol.id = fineId;
+  cSol.id = coarseId;
+}

mluq/src/estimators/StochasticCollocation.hpp

@@ -3,20 +3,41 @@
 
 #include "Estimator.hpp"
 
+#include <utility>
+
 
 class StochasticCollocation : public Estimator {
-private:
+protected:
+  int verbose = 1;
+
+  bool onlyFine;
+
+  bool parallel;
+
+  SampleID fineId;
+
+  SampleID coarseId;
+
+  MeshesCreator meshesCreator;
+
+  PDESolverCreator pdeSolverCreator;
+
+  void updateSampleIds(SampleSolution &fSol, SampleSolution &cSol);
+
+  void method();
 
 public:
-  StochasticCollocation() : Estimator(0.0, 0, 0) {}
+  StochasticCollocation(double epsilon, int level, int stochLevel,
+                        bool onlyFine, bool parallel,
+                        MeshesCreator meshesCreator = MeshesCreator(),
+                        PDESolverCreator pdeCreator = PDESolverCreator()) :
+    Estimator(epsilon, level, 0),
+    fineId(SampleID(level, 0, false)), coarseId(SampleID(level, 0, true)),
+    meshesCreator(std::move(meshesCreator)), pdeSolverCreator(std::move(pdeCreator)) {}
 
-  std::string Name() const override {
-    return "StochasticCollocation";
-  }
+  std::string Name() const override { return "StochasticCollocation"; }
 
-  void Method() override {
-    // Todo here is some stuff todo
-  }
+  void Method() override;
 };
 
 #endif //STOCHASTICCOLLOCATION_HPP

mluq/src/estimators/datastructure/WelfordAggregate.hpp

@@ -86,7 +86,10 @@ struct SVar {
     Ycomm = Y2comm / (ctr.Mcomm - 1);
   }
 
-  void UpdateParallel(int commSplit) {
+  void UpdateParallel(double C2, double Q2, double Y2, SampleCounter ctr) {
+    C = C2 / (ctr.M - 1);
+    Q = Q2 / (ctr.M - 1);
+    Y = Y2 / (ctr.M - 1);
   }
 
   friend Logging &operator<<(Logging &s, const SVar &sVar) {
@@ -136,7 +139,8 @@ struct Kurtosis {
 
   double Ycomm = 0.0;
 
-  void Update() {//Averages avgs, Variances vars) {
+  void Update() {
+//Averages avgs, Variances vars) {
 //    Q = (avgs.Q4 - 4.0 * avgs.Q3 * avgs.Q + 6.0 * avgs.Q2 * avgs.Q * avgs.Q -
 //      3.0 * avgs.Q * avgs.Q * avgs.Q * avgs.Q) / vars.Q / vars.Q;
 //    Y = (avgs.Y4 - 4.0 * avgs.Y3 * avgs.Y + 6.0 * avgs.Y2 * avgs.Y * avgs.Y -
@@ -194,7 +198,7 @@ public:
   bool parallel = true;
 
 public:
-  WelfordAggregate(int dM) {
+  explicit WelfordAggregate(int dM) {
     UpdateSampleCounter(dM);
   }
 
@@ -233,24 +237,11 @@ public:
     Q2 = PPM->SumAcrossComm(Q2comm, commSplit);
     Y2 = PPM->SumAcrossComm(Y2comm, commSplit);
 
-//        delta = avg_b - avg_a
-//        M2 = M2_a + M2_b + delta ** 2 * n_a * n_b / n
-
-    sVar.C = C2 / (ctr.M - 1);
-    sVar.Q = Q2 / (ctr.M - 1);
-    sVar.Y = Y2 / (ctr.M - 1);
-
-//    def parallel_variance(n_a, avg_a, M2_a, n_b, avg_b, M2_b):
-//        n = n_a + n_b
-//        delta = avg_b - avg_a
-//        M2 = M2_a + M2_b + delta ** 2 * n_a * n_b / n
-//        var_ab = M2 / (n - 1)
-//        return var_ab
+    sVar.UpdateParallel(C2, Q2, Y2, ctr);
   }
 
-  int index() {
+  int index() const {
     return ctr.M + ctr.Mcomm + ctr.dMcomm * PPM->Proc(0);
-//    return M  + dMcomm * PPM->Proc(0);
   }
 
   void UpdateSampleCounter(int dM) {

mluq/src/generators/SparseGridGenerator.hpp

@@ -2,11 +2,11 @@
 #define SPARSEGRIDGENERATOR_HPP
 
 #include "TasmanianSparseGrid.hpp"
+#include "SampleGenerator.hpp"
 
 
-// Todo Use RVector and such
-
-class SparseGridGenerator {
+//template<typename T>
+class SparseGridGenerator : public SampleGenerator<RVector> {
 protected:
   TasGrid::TasmanianSparseGrid grid;
 
@@ -18,32 +18,68 @@ protected:
 
   int outputs;
 
+  double weight = 0.0;
+
+  RVector sample {};
+
+  std::vector<double> weights{};
+
+  std::vector<RVector> samples{};
+
+  void fillSampleVector() {
+    std::vector<double> points = grid.getPoints();
+    for (auto it_start = points.begin(); it_start != points.end(); it_start += dimension)
+      samples.emplace_back(std::vector<double>(it_start, it_start + dimension));
+  }
+
+  void drawSample(const SampleID &id) override {
+    sample = samples[id.number];
+    weight = weights[id.number];
+  };
+
 public:
-  SparseGridGenerator(int dimension, int outputs,
+  SparseGridGenerator(const Meshes &meshes, int dimension, int outputs, int level,
                       TasGrid::TypeDepth depth = TasGrid::type_level,
                       TasGrid::TypeOneDRule rule = TasGrid::rule_clenshawcurtis) :
-    dimension(dimension), outputs(outputs), depth(depth), rule(rule) {
+    SampleGenerator(meshes), dimension(dimension), outputs(outputs),
+    depth(depth), rule(rule) {
+    grid.makeGlobalGrid(dimension, outputs, level, depth, rule);
+    weights = grid.getQuadratureWeights();
+    fillSampleVector();
   }
 
   void CreateGlobalGrid(int level) {
     grid.makeGlobalGrid(dimension, outputs, level, depth, rule);
   }
 
-  std::vector<double> GetPoints() {
-    return grid.getPoints();
+  std::vector<RVector> GetSamples() {
+    return samples;
   }
 
   std::vector<double> GetWeights() {
-    return grid.getQuadratureWeights();
+    return weights;
+  }
+
+  double SampleWeight(const SampleID &id) {
+    return weights[id.number];
   }
 
   int GetNumPoints() {
     return grid.getNumPoints();
   }
 
-  double Quadrature(double func (double, double)) {
+  int GetStochDimension() {
+    return grid.getNumDimensions();
+  }
+
+  RVector EvalSample() const override {
+    return sample;
+  }
+
+  // todo move to stochastic collocation class
+  double Quadrature(double func(double, double)) {
     double I = 0.0;
-    std::vector<double> points = GetPoints();
+    std::vector<double> points = grid.getPoints();
     std::vector<double> weights = GetWeights();
     int num_points = GetNumPoints();
     for (int i = 0; i < num_points; i++) {
@@ -53,6 +89,8 @@ public:
     }
     return I;
   }
+
+  string Name() const override { return "SparseGridGenerator"; };
 };
 
 #endif //SPARSEGRIDGENERATOR_HPP

mluq/src/pdesolvers/EllipticPDESolver.cpp

@@ -5,6 +5,7 @@ void EllipticPDESolver::run(SampleSolution &solution) {
   newton->operator()(*assemble, solution.U);
 }
 
+// Todo maybe multiply weights to Q
 void EllipticPDESolver::computeQ(SampleSolution &solution) {
   if (quantity == "L2") solution.Q = assemble->L2(solution.U);
   else if (quantity == "H1") solution.Q = assemble->H1(solution.U);

mluq/src/pdesolvers/PDESolver.hpp

@@ -32,6 +32,10 @@ protected:
 
   virtual void plotSolution(SampleSolution &solution) = 0;
 
+  void weightSolution(SampleSolution &solution) const {
+    solution.Q = GetProblem()->SampleWeight(solution.id) * solution.Q;
+  }
+
 public:
   PDESolver(const Meshes &meshes, const std::string &quantity,
             const std::string &costMeasure) :
@@ -56,6 +60,8 @@ public:
     computeQ(solution);
     computeCost(solution);
     plotSolution(solution);
+    weightSolution(solution);
+    // Todo other idea: Add weight as class member of solution
     vout(2) << "Q=" << solution.Q << " cost=" << solution.Cost << endl;
     mout.EndBlock(verbose <= 1);
   }
@@ -79,7 +85,7 @@ protected:
   void run(SampleSolution &solution) override {}
 
   void computeQ(SampleSolution &solution) override {
-    if (quantity == "GeneratorValue") solution.Q = assemble->GeneratorValue();
+    if (quantity == "FunctionEvaluation") solution.Q = assemble->FunctionEvaluation();
     else Exit("Quantity of interest not implemented")
   }
 
@@ -98,12 +104,7 @@ public:
                  const Meshes &meshes,
                  const std::string &quantity = "L2",
                  const std::string &costMeasure = "size") :
-    PDESolver(meshes, quantity, costMeasure), assemble(assemble) {
-    if (verbose > 0) {
-      PrintInfo();
-      assemble->PrintInfo();
-    }
-  }
+    PDESolver(meshes, quantity, costMeasure), assemble(assemble) {}
 
   IAssemble *GetAssemble() const override { return assemble; }
 

mluq/src/pdesolvers/PDESolverCreator.hpp

@@ -10,7 +10,7 @@ private:
   std::string _pc;
   std::string _model = "DummyPDESolver";
   std::string _problem = "StochasticDummyProblem";
-  std::string _quantity = "GeneratorValue";
+  std::string _quantity = "FunctionEvaluation";
   std::string _costMeasure = "size";
   std::string _linearSolver = "GMRES";
 

mluq/src/pdesolvers/assembling/IStochasticDummyAssemble.hpp

@@ -55,8 +55,8 @@ public:
         return 0.0;
     };
 
-    double GeneratorValue() const {
-        return problem->GeneratorValue();
+    double FunctionEvaluation() const {
+        return problem->FunctionEvaluation();
     }
 };
 

mluq/src/problems/IStochasticProblem.cpp

 #include "IStochasticProblem.hpp"
 
+
 StochasticDummyProblem *
 CreateStochasticDummyProblem(const std::string &problemName, const Meshes &meshes) {
-    if (problemName == "StochasticDummyScalarGeneratorProblem")
-        return new StochasticDummyScalarGeneratorProblem(meshes);
-    if (problemName == "StochasticDummyComplexGeneratorProblem")
-        return new StochasticDummyComplexGeneratorProblem(meshes);
-    if (problemName == "StochasticDummyVectorFieldGeneratorProblem")
-        return new StochasticDummyVectorFieldGeneratorProblem(meshes);
-    if (problemName == "StochasticDummyTensorGeneratorProblem")
-        return new StochasticDummyTensorGeneratorProblem(meshes);
-    if (problemName == "StochasticDummyRVectorGeneratorProblem")
-        return new StochasticDummyRVectorGeneratorProblem(meshes);
-    if (problemName == "StochasticDummyCVectorGeneratorProblem")
-        return new StochasticDummyCVectorGeneratorProblem(meshes);
-    if (problemName == "StochasticDummyRMatrixGeneratorProblem")
-        return new StochasticDummyRMatrixGeneratorProblem(meshes);
-    if (problemName == "StochasticDummyCMatrixGeneratorProblem")
-        return new StochasticDummyCMatrixGeneratorProblem(meshes);
-    Exit(problemName + " not found")
+  if (problemName == "ScalarGeneratorProblem")
+    return new ScalarGeneratorProblem(meshes);
+  if (problemName == "ComplexGeneratorProblem")
+    return new ComplexGeneratorProblem(meshes);
+  if (problemName == "VectorFieldGeneratorProblem")
+    return new VectorFieldGeneratorProblem(meshes);
+  if (problemName == "TensorGeneratorProblem")
+    return new TensorGeneratorProblem(meshes);
+  if (problemName == "RVectorGeneratorProblem")
+    return new RVectorGeneratorProblem(meshes);
+  if (problemName == "CVectorGeneratorProblem")
+    return new CVectorGeneratorProblem(meshes);
+  if (problemName == "RMatrixGeneratorProblem")
+    return new RMatrixGeneratorProblem(meshes);
+  if (problemName == "CMatrixGeneratorProblem")
+    return new CMatrixGeneratorProblem(meshes);
+  if (problemName == "SparseGrid1DGeneratorProblem")
+    return new SparseGrid1DGeneratorProblem(meshes);
+  if (problemName == "SparseGrid2DGeneratorProblem")
+    return new SparseGrid2DGeneratorProblem(meshes);
+  Exit(problemName + " not found")
 }
\ No newline at end of file

mluq/src/problems/IStochasticProblem.hpp

@@ -21,6 +21,10 @@ public:
 
   virtual void DrawSample(const SampleID &id) = 0;
 
+  virtual double SampleWeight(const SampleID &id) { return 1.0; }
+
+  virtual int NumOfSamples() { return 0; };
+
   virtual string Name() const = 0;
 };
 
@@ -29,13 +33,13 @@ public:
   explicit StochasticDummyProblem(const Meshes &meshes) :
     IStochasticProblem(meshes) {}
 
-  virtual double GeneratorValue() = 0;
+  virtual double FunctionEvaluation() = 0;
 };
 
-class StochasticDummyScalarGeneratorProblem : public StochasticDummyProblem {
+class ScalarGeneratorProblem : public StochasticDummyProblem {
   NormalDistributionReal scalarGenerator;
 public:
-  explicit StochasticDummyScalarGeneratorProblem(const Meshes &meshes) :
+  explicit ScalarGeneratorProblem(const Meshes &meshes) :
     StochasticDummyProblem(meshes),
     scalarGenerator(NormalDistributionReal(meshes)) {}
 
@@ -43,19 +47,19 @@ public:
     scalarGenerator.DrawSample(id);
   }
 
-  double GeneratorValue() override {
+  double FunctionEvaluation() override {
     return this->scalarGenerator.EvalSample();
   }
 
   string Name() const override {
-    return "StochasticDummyScalarGeneratorProblem";
+    return "ScalarGeneratorProblem";
   }
 };
 
-class StochasticDummyComplexGeneratorProblem : public StochasticDummyProblem {
+class ComplexGeneratorProblem : public StochasticDummyProblem {
   NormalDistributionComplex complexGenerator;
 public:
-  explicit StochasticDummyComplexGeneratorProblem(const Meshes &meshes) :
+  explicit ComplexGeneratorProblem(const Meshes &meshes) :
     StochasticDummyProblem(meshes),
     complexGenerator(NormalDistributionComplex(meshes)) {}
 
@@ -63,53 +67,53 @@ public:
     complexGenerator.DrawSample(id);
   }
 
-  double GeneratorValue() override {
+  double FunctionEvaluation() override {
     return abs_compl(this->complexGenerator.EvalSample());
   }
 
   string Name() const override {
-    return "StochasticDummyComplexGeneratorProblem";
+    return "ComplexGeneratorProblem";
   }
 };
 
 // TODO
-class StochasticDummyVectorFieldGeneratorProblem : public StochasticDummyProblem {
+class VectorFieldGeneratorProblem : public StochasticDummyProblem {
 public:
-  explicit StochasticDummyVectorFieldGeneratorProblem(const Meshes &meshes)
+  explicit VectorFieldGeneratorProblem(const Meshes &meshes)
     : StochasticDummyProblem(meshes) {}
 
   void DrawSample(const SampleID &id) override {}
 
-  double GeneratorValue() override {
+  double FunctionEvaluation() override {
     return 0.0;
   }
 
   string Name() const override {
-    return "StochasticDummyVectorFieldGeneratorProblem";
+    return "VectorFieldGeneratorProblem";
   }
 };
 
 // TODO
-class StochasticDummyTensorGeneratorProblem : public StochasticDummyProblem {
+class TensorGeneratorProblem : public StochasticDummyProblem {
 public:
-  explicit StochasticDummyTensorGeneratorProblem(const Meshes &meshes) :
+  explicit TensorGeneratorProblem(const Meshes &meshes) :
     StochasticDummyProblem(meshes) {}
 
   void DrawSample(const SampleID &id) override {}
 
-  double GeneratorValue() override {
+  double FunctionEvaluation() override {
     return 0.0;
   }
 
   string Name() const override {
-    return "StochasticDummyTensorGeneratorProblem";
+    return "TensorGeneratorProblem";
   }
 };
 
-class StochasticDummyRVectorGeneratorProblem : public StochasticDummyProblem {
+class RVectorGeneratorProblem : public StochasticDummyProblem {
   NormalDistributionRVector rVectorGenerator;
 public:
-  explicit StochasticDummyRVectorGeneratorProblem(const Meshes &meshes) :
+  explicit RVectorGeneratorProblem(const Meshes &meshes) :
     StochasticDummyProblem(meshes),
     rVectorGenerator(NormalDistributionRVector(meshes)) {}
 
@@ -117,19 +121,19 @@ public:
     rVectorGenerator.DrawSample(id);
   }
 
-  double GeneratorValue() override {
+  double FunctionEvaluation() override {
     return this->rVectorGenerator.EvalSample().Mean();
   }
 
   string Name() const override {
-    return "StochasticDummyRVectorGeneratorProblem";
+    return "RVectorGeneratorProblem";
   }
 };
 
-class StochasticDummyCVectorGeneratorProblem : public StochasticDummyProblem {
+class CVectorGeneratorProblem : public StochasticDummyProblem {
   NormalDistributionCVector cVectorGenerator;
 public:
-  explicit StochasticDummyCVectorGeneratorProblem(const Meshes &meshes) :
+  explicit CVectorGeneratorProblem(const Meshes &meshes) :
     StochasticDummyProblem(meshes),
     cVectorGenerator(NormalDistributionCVector(meshes)) {}
 
@@ -137,19 +141,19 @@ public:
     cVectorGenerator.DrawSample(id);
   }
 
-  double GeneratorValue() override {
+  double FunctionEvaluation() override {
     return abs_compl(this->cVectorGenerator.EvalSample().Mean());
   }
 
   string Name() const override {
-    return "StochasticDummyCVectorGeneratorProblem";
+    return "CVectorGeneratorProblem";
   }
 };
 
-class StochasticDummyRMatrixGeneratorProblem : public StochasticDummyProblem {
+class RMatrixGeneratorProblem : public StochasticDummyProblem {
   NormalDistributionRMatrix rMatrixGenerator;
 public:
-  explicit StochasticDummyRMatrixGeneratorProblem(const Meshes &meshes) :
+  explicit RMatrixGeneratorProblem(const Meshes &meshes) :
     StochasticDummyProblem(meshes),
     rMatrixGenerator(NormalDistributionRMatrix(meshes)) {}
 
@@ -157,19 +161,19 @@ public:
     rMatrixGenerator.DrawSample(id);
   }
 
-  double GeneratorValue() override {
+  double FunctionEvaluation() override {
     return abs_compl(this->rMatrixGenerator.EvalSample().Mean());
   }