added Welford files

mlmc/src/montecarlo/CMakeLists.txt

@@ -4,6 +4,7 @@ add_library(MONTECARLO STATIC
         datastructure/Errors.cpp
         datastructure/Exponents.cpp
         datastructure/EmpiricMeasures.cpp
+        datastructure/WelfordAggregate.cpp
         datastructure/EmpiricMeasureLevelMaps.cpp
         )
 target_link_libraries(MONTECARLO PDESOLVER MPP_LIBRARIES)

mlmc/src/montecarlo/datastructure/WelfordAggregate.cpp

+#include "WelfordAggregate.hpp"
+

mlmc/src/montecarlo/datastructure/WelfordAggregate.hpp

+#ifndef WELFORDAGGREGATE_HPP
+#define WELFORDAGGREGATE_HPP
+
+#include "basics/Sample.hpp"
+#include "Parallel.hpp"
+
+
+struct SampleCounter {
+  int M = 0;
+
+  int dM = 0;
+
+  int Mcomm = 0;
+
+  int dMcomm = 0;
+
+  void Update() {
+    Mcomm++;
+    dMcomm--;
+  }
+
+  void UpdateParallel(int commSplit) {
+    pout << "test" << endl;
+    M = PPM->SumOnCommSplit(Mcomm, 0) / PPM->Size(commSplit);
+    dM = PPM->SumOnCommSplit(dMcomm, 0) / PPM->Size(commSplit);
+  }
+
+  friend Logging &operator<<(Logging &s, const SampleCounter &ctr) {
+    return s << "M=" << ctr.M << " dM=" << ctr.dM
+             << " MComm=" << ctr.Mcomm << " dMComm=" << ctr.dMcomm << endl;
+  }
+};
+
+class WelfordAggregate {
+public:
+  SampleCounter ctr;
+
+
+  double MeanQcomm = 0.0;
+  double Q2comm = 0.0;
+  double SVarQcomm = 0.0;
+
+  double MeanYcomm = 0.0;
+  double Y2comm = 0.0;
+  double SVarYcomm = 0.0;
+
+  double MeanQ = 0.0;
+  double Q2 = 0.0;
+  double SVarQ = 0.0;
+
+  double MeanY = 0.0;
+  double Y2 = 0.0;
+  double SVarY = 0.0;
+
+  int commSplit = 0;
+  double Cost = 0.0;
+
+  bool parallel = true;
+
+public:
+  void Update(const SampleSolution &fineSolution, const SampleSolution &coarseSolution) {
+    double newC = fineSolution.Cost;
+    double newQ = fineSolution.Q;
+    double newY = fineSolution.Q - coarseSolution.Q;
+    Update(newC, newQ, newY);
+  }
+
+  void Update(double newC, double newQ, double newY) {
+    ctr.Update();
+
+    double dCost = newC - Cost;
+    Cost += dCost / ctr.Mcomm;
+
+    double dQ = newQ - MeanQcomm;
+    MeanQcomm += dQ / ctr.Mcomm;
+    double dQ2 = newQ - MeanQcomm;
+    Q2comm += dQ * dQ2;
+    SVarQcomm = Q2comm / (ctr.Mcomm - 1);
+
+    double dY = newY - MeanYcomm;
+    MeanYcomm += dY / ctr.Mcomm;
+    double dY2 = newY - MeanYcomm;
+    Y2comm += dY * dY2;
+    SVarYcomm = Y2comm / (ctr.Mcomm - 1);
+  }
+
+  void UpdateParallel() {
+    ctr.UpdateParallel(commSplit);
+
+    MeanQ =
+      (PPM->SumOnCommSplit(ctr.Mcomm * MeanQcomm, 0) / PPM->Size(commSplit)) / ctr.M;
+    MeanY =
+      (PPM->SumOnCommSplit(ctr.Mcomm * MeanYcomm, 0) / PPM->Size(commSplit)) / ctr.M;
+
+    Q2 = PPM->SumOnCommSplit(Q2, 0) / PPM->Size(commSplit);
+    Y2 = PPM->SumOnCommSplit(Y2, 0) / PPM->Size(commSplit);
+
+//        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
+  }
+
+  int index() {
+    return ctr.M + ctr.Mcomm + ctr.dMcomm * PPM->Proc(0);
+//    return M  + dMcomm * PPM->Proc(0);
+  }
+
+  void UpdateSampleCounter(int dM) {
+    this->ctr.dM = dM;
+    if (parallel) {
+      for (int i = 0; i < ceil(log2(PPM->Size(0))) + 1; i++) {
+        if (dM >= PPM->Size(i)) {
+          ctr.dMcomm = ceil((double) dM / PPM->Size(i));
+          commSplit = PPM->MaxCommSplit() - i;
+          break;
+        }
+      }
+    } else {
+      ctr.dMcomm = ctr.dM;
+      commSplit = 0;
+    }
+  }
+
+  friend Logging &operator<<(Logging &s, const WelfordAggregate &aggregate) {
+    return s << "M=" << aggregate.ctr.M
+             << " dM=" << aggregate.ctr.dM
+             << " Mcomm=" << aggregate.ctr.Mcomm
+             << " dMcomm=" << aggregate.ctr.dMcomm << endl
+             << "MeanQ=" << aggregate.MeanQ << " MeanY=" << aggregate.MeanY
+             << " SVarQ=" << aggregate.SVarQ << " SVarY=" << aggregate.SVarY << endl;
+  }
+};
+
+#endif //WELFORDAGGREGATE_HPP