Commit 813d2a1e authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

deploy WelfordAggregate

parent 6c874445
......@@ -3,16 +3,16 @@
void MonteCarlo::Method() {
mout.StartBlock("MC Method on l=" + to_string(level.fine));
vout(1) << ctr;
vout(1) << aggregate;
method();
ctr.Update(); // todo should move to mc method and use while
vout(1) << ctr;
aggregate.UpdateParallel();
// ctr.Update(); // todo should move to mc method and use while
vout(1) << aggregate;
// sums.UpdateParallel(ctr);
vout(2) << sums;
avgs.Update(sums, ctr.M);
avgs.Update(sums, aggregate.M);
vout(1) << avgs;
vars.Update(avgs);
......@@ -27,11 +27,12 @@ void MonteCarlo::Method() {
void MonteCarlo::method() {
SampleSolution coarseSolution(pdeSolver->GetDisc(), coarseId);
SampleSolution fineSolution(pdeSolver->GetDisc(), fineId);
int m = ctr.M + ctr.dMComm * PPM->Proc(0); // start index count
for (; m < ctr.M + ctr.dMComm * PPM->Proc(0) + ctr.dMComm; m++) { // change to while dMComm != 0
computeSampleSolution(m, fineId, fineSolution);
// int m = ctr.M + ctr.dMComm * PPM->Proc(0); // start index count
// for (; m < ctr.M + ctr.dMComm * PPM->Proc(0) + ctr.dMComm; m++) { // change to while dMComm != 0
while(aggregate.dMcomm != 0) {
computeSampleSolution(aggregate.index(), fineId, fineSolution);
if (onlyFine) coarseSolution.Init();
else computeSampleSolution(m, coarseId, coarseSolution);
else computeSampleSolution(aggregate.index(), coarseId, coarseSolution);
sums.Update(fineSolution, coarseSolution);
vout(3) << sums;
aggregate.Update(fineSolution, coarseSolution);
......
......@@ -26,7 +26,7 @@ protected:
public:
WelfordAggregate aggregate;
SampleCounter ctr;
SampleCounter ctr; // todo remove, is still here for test
Sums sums;
......@@ -49,7 +49,7 @@ public:
onlyFine(onlyFine),
pdeSolverCreator(PDESolverCreator()),
meshesCreator(MeshesCreator(pdeSolverCreator.GetMeshName()).
WithCommSplit(ctr.commSplit).
WithCommSplit(aggregate.commSplit). // Todo: Is this correct?
WithPLevel(level.coarse).
WithLevel(level.fine)) {
Init(dM);
......@@ -76,8 +76,8 @@ public:
fineId.coarse = false;
coarseId.level = level;
coarseId.coarse = true;
ctr.parallel = parallel;
ctr.UpdateSampleCounter(dM);
aggregate.parallel = parallel;
aggregate.UpdateSampleCounter(dM);
}
~MonteCarlo() {
......
......@@ -37,7 +37,7 @@ void MultilevelMonteCarlo::adaptiveMethod() {
bool converged = false;
while (!converged) {
for (auto &[level, mc] : mcMap)
if (mc.ctr.dM > 0)
if (mc.aggregate.dM > 0)
mc.Method();
mcMap.UpdateSampleCounter(epsilon);
......
......@@ -13,7 +13,7 @@ void MonteCarloMap::UpdateSampleCounter(double epsilon) {
optimalM = (int) (ceil(2 * pow(epsilon, -2) * factor *
sqrt(mc.vars.Y / mc.avgs.Cost)));
if (optimalM == 1) optimalM++; // Hack
mc.ctr.UpdateSampleCounter(optimalM - mc.ctr.M);
mc.aggregate.UpdateSampleCounter(optimalM - mc.aggregate.M);
}
}
......@@ -33,7 +33,7 @@ void MonteCarloMap::AppendLevel(double epsilon, Exponents exponents, Level newLe
void SampleCounterMap::Update(const MonteCarloMap &mcMap) {
for (auto &[level, mc]: mcMap)
_levelMap[level] = mc.ctr;
_levelMap[level] = SampleCounter(mc.aggregate);
}
bool SampleCounterMap::NoSamplesLeft() {
......
#ifndef EMPIRICMEASURES_HPP
#define EMPIRICMEASURES_HPP
#include "SampleCounter.hpp"
#include "basics/Sample.hpp"
#include "Parallel.hpp"
struct Sums {
......@@ -29,8 +28,6 @@ struct Sums {
void Update(SampleSolution &fineSolution, SampleSolution &coarseSolution);
void UpdateParallel(SampleCounter ctr);
friend Logging &operator<<(Logging &s, const Sums &sums) {
return s << "Sum(Q)=" << sums.Q << " Sum(Q^2)=" << sums.Q2
<< " Sum(Q^3)=" << sums.Q3 << " Sum(Q^4)=" << sums.Q4 << endl
......@@ -96,11 +93,11 @@ struct WelfordAggregate { // todo make class
double SVarQ = 0.0;
double MeanY = 0.0;
double MeanY2 = 0.0;
double Y2 = 0.0;
double SVarY = 0.0;
int commM = 0;
int dMComm = 0;
int Mcomm = 0;
int dMcomm = 0;
int dM = 0;
int M = 0;
......@@ -114,29 +111,33 @@ public:
double newC = fineSolution.Cost;
double newQ = fineSolution.Q;
double newY = fineSolution.Q - coarseSolution.Q;
Update(newC, newQ, newY);
}
commM += 1;
void Update(double newC, double newQ, double newY) {
Mcomm++;
dMcomm--;
double dCost = newC - Cost;
Cost += dCost / commM;
Cost += dCost / Mcomm;
double dQ = newQ - MeanQ;
MeanQ += dQ / commM;
MeanQ += dQ / Mcomm;
double dQ2 = newQ - MeanQ;
Q2 += dQ * dQ2;
SVarQ = Q2 / (commM - 1);
SVarQ = Q2 / (Mcomm - 1);
double dY = newY - MeanY;
MeanY += dY / commM;
MeanY += dY / Mcomm;
double dY2 = newY - MeanY;
MeanY2 += dY * dY2;
SVarY = MeanY2 / (commM - 1);
Y2 += dY * dY2;
SVarY = Y2 / (Mcomm - 1);
}
void UpdateParallel() {
M = PPM->SumOnCommSplit(commM, commSplit);
M = PPM->SumOnCommSplit(Mcomm, 0) / PPM->Size(commSplit);
dM = PPM->SumOnCommSplit(dMcomm, 0) / PPM->Size(commSplit);
pout << M << endl;
// 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
......@@ -145,6 +146,55 @@ public:
// return var_ab
}
int index() {
return M + dMcomm * PPM->Proc(0);
}
void UpdateSampleCounter(int dM) {
this->dM = dM;
if (parallel) {
for (int i = 0; i < ceil(log2(PPM->Size(0))) + 1; i++) {
if (dM >= PPM->Size(i)) {
dMcomm = ceil((double) dM / PPM->Size(i));
commSplit = PPM->MaxCommSplit() - i;
break;
}
}
} else {
dMcomm = dM;
commSplit = 0;
}
}
friend Logging &operator<<(Logging &s, const WelfordAggregate &aggregate) {
return s << "M=" << aggregate.M << " dM=" << aggregate.dM << endl
<< "Mcomm=" << aggregate.Mcomm << " dMcomm=" << aggregate.dMcomm << endl
<< "MeanQ=" << aggregate.MeanQ << " MeanY=" << aggregate.MeanY << endl
<< "SVarQ=" << aggregate.SVarQ << " SVarY=" << aggregate.SVarY << endl;
}
};
struct SampleCounter {
int M = 0;
int dM = 0;
int dMComm = 0;
int commSplit = 0;
int parallel = true;
SampleCounter() {};
SampleCounter(WelfordAggregate aggregate) {
M = aggregate.M;
dM = aggregate.dM;
dMComm = aggregate.dMcomm;
parallel = aggregate.parallel;
commSplit = aggregate.commSplit;
}
SampleCounter(int dM, bool parallel = true) : parallel(parallel) {
UpdateSampleCounter(dM);
}
void UpdateSampleCounter(int dM) {
this->dM = dM;
if (parallel) {
......@@ -160,6 +210,19 @@ public:
commSplit = 0;
}
}
SampleCounter(int M, int dM) :
M(M), dM(dM) {}
void Update() {
M += dM;
dMComm = 0;
dM = 0;
}
friend Logging &operator<<(Logging &s, const SampleCounter &ctr) {
return s << "M=" << ctr.M << " dM=" << ctr.dM << " dMComm=" << ctr.dMComm << endl;
}
};
#endif //EMPIRICMEASURES_HPP
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