#ifndef EMPIRICMEASURES_HPP #define EMPIRICMEASURES_HPP #include "basics/Sample.hpp" #include "Parallel.hpp" struct Sums { double Q = 0.0; double Q2 = 0.0; double Q3 = 0.0; double Q4 = 0.0; double Y = 0.0; double Y2 = 0.0; double Y3 = 0.0; double Y4 = 0.0; double Cost = 0.0; // Vector U; Todo; Also be careful with instabilities of difference // // Vector V; Sums() {}; Sums(double Q, double Q2, double Q3, double Q4, double Y, double Y2, double Y3, double Y4, double Cost) : Q(Q), Q2(Q2), Q3(Q3), Q4(Q4), Y(Y), Y2(Y2), Y3(Y3), Y4(Y4), Cost(Cost) {} void Update(SampleSolution &fineSolution, SampleSolution &coarseSolution); 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 << "Sum(Y)=" << sums.Y << " Sum(Y^2)=" << sums.Y2 << " Sum(Y^3)=" << sums.Y3 << " Sum(Y^4)=" << sums.Y4 << endl; } }; struct Averages { double Q = 0.0; double Q2 = 0.0; double Q3 = 0.0; double Q4 = 0.0; double Y = 0.0; double Y2 = 0.0; double Y3 = 0.0; double Y4 = 0.0; double Cost = 0.0; Averages() {}; Averages(double Q, double Q2, double Q3, double Q4, double Y, double Y2, double Y3, double Y4, double Cost) : Q(Q), Q2(Q2), Q3(Q3), Q4(Q4), Y(Y), Y2(Y2), Y3(Y3), Y4(Y4), Cost(Cost) {} void Update(Sums sums, int M); friend Logging &operator<<(Logging &s, const Averages &avgs) { return s << "|E[Q]|=" << avgs.Q << " |E[Y]|=" << avgs.Y << endl; } }; struct Variances { double Q = 0.0; double Y = 0.0; Variances() {}; Variances(double Q, double Y) : Q(Q), Y(Y) {} void Update(Averages avgs); friend Logging &operator<<(Logging &s, const Variances &vars) { return s << "V[Q]=" << vars.Q << " V[Y]=" << vars.Y << endl; } }; struct Kurtosis { double Q = 0.0; double Y = 0.0; void Update(Averages avgs, Variances vars); friend Logging &operator<<(Logging &s, const Kurtosis &kurtosis) { return s << "K[Q]=" << kurtosis.Q << " K[Y]=" << kurtosis.Y << endl; } }; struct WelfordAggregate { // todo make class int Mcomm = 0; int dMcomm = 0; 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 dM = 0; int M = 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) { Mcomm++; dMcomm--; double dCost = newC - Cost; Cost += dCost / Mcomm; double dQ = newQ - MeanQcomm; MeanQcomm += dQ / Mcomm; double dQ2 = newQ - MeanQcomm; Q2comm += dQ * dQ2; SVarQcomm = Q2comm / (Mcomm - 1); double dY = newY - MeanYcomm; MeanYcomm += dY / Mcomm; double dY2 = newY - MeanYcomm; Y2comm += dY * dY2; SVarYcomm = Y2comm / (Mcomm - 1); } void UpdateParallel() { M = PPM->SumOnCommSplit(Mcomm, 0) / PPM->Size(commSplit); dM = PPM->SumOnCommSplit(dMcomm, 0) / PPM->Size(commSplit); MeanQ = (PPM->SumOnCommSplit(Mcomm * MeanQcomm, 0) / PPM->Size(commSplit)) / M; MeanY = (PPM->SumOnCommSplit(Mcomm * MeanYcomm, 0) / PPM->Size(commSplit)) / 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 M + Mcomm + dMcomm * PPM->Proc(0); // 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 << " Mcomm=" << aggregate.Mcomm << " dMcomm=" << aggregate.dMcomm << endl << "MeanQ=" << aggregate.MeanQ << " MeanY=" << aggregate.MeanY << " 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); 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