EmpiricMeasures.hpp 5.57 KB
Newer Older
niklas.baumgarten's avatar
niklas.baumgarten committed
1
2
#ifndef EMPIRICMEASURES_HPP
#define EMPIRICMEASURES_HPP
3

4
#include "basics/Sample.hpp"
niklas.baumgarten's avatar
niklas.baumgarten committed
5
#include "Parallel.hpp"
6

7
8

struct Sums {
9
10
11
12
13
14
15
16
17
  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;
18
19
20
21
22

//    Vector U; Todo; Also be careful with instabilities of difference
//
//    Vector V;

23
  Sums() {};
24

25
26
27
  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) {}
28

29
  void Update(SampleSolution &fineSolution, SampleSolution &coarseSolution);
30

31
32
33
34
35
36
  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;
  }
37
38
39
};

struct Averages {
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
  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;
  }
61
62
63
};

struct Variances {
64
65
  double Q = 0.0;
  double Y = 0.0;
66

67
  Variances() {};
68

69
  Variances(double Q, double Y) : Q(Q), Y(Y) {}
70

71
  void Update(Averages avgs);
72

73
74
75
  friend Logging &operator<<(Logging &s, const Variances &vars) {
    return s << "V[Q]=" << vars.Q << " V[Y]=" << vars.Y << endl;
  }
76
77
78
};

struct Kurtosis {
79
80
  double Q = 0.0;
  double Y = 0.0;
81

82
  void Update(Averages avgs, Variances vars);
83

84
85
86
  friend Logging &operator<<(Logging &s, const Kurtosis &kurtosis) {
    return s << "K[Q]=" << kurtosis.Q << " K[Y]=" << kurtosis.Y << endl;
  }
87
88
};

niklas.baumgarten's avatar
niklas.baumgarten committed
89
struct WelfordAggregate { // todo make class
90
91
  int Mcomm = 0;
  int dMcomm = 0;
92

93
94
95
  double MeanQcomm = 0.0;
  double Q2comm = 0.0;
  double SVarQcomm = 0.0;
96

97
98
99
  double MeanYcomm = 0.0;
  double Y2comm = 0.0;
  double SVarYcomm = 0.0;
niklas.baumgarten's avatar
niklas.baumgarten committed
100

101
102
103
  double MeanQ = 0.0;
  double Q2 = 0.0;
  double SVarQ = 0.0;
104

105
106
107
  double MeanY = 0.0;
  double Y2 = 0.0;
  double SVarY = 0.0;
108

109
110
  int dM = 0;
  int M = 0;
niklas.baumgarten's avatar
niklas.baumgarten committed
111

112
113
  int commSplit = 0;
  double Cost = 0.0;
114

115
  bool parallel = true;
niklas.baumgarten's avatar
niklas.baumgarten committed
116

117
public:
118
  void Update(const SampleSolution &fineSolution, const SampleSolution &coarseSolution) {
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
    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() {
niklas.baumgarten's avatar
niklas.baumgarten committed
146
147
    M = PPM->SumOnCommSplit(Mcomm, 0) / PPM->Size(commSplit);
    dM = PPM->SumOnCommSplit(dMcomm, 0) / PPM->Size(commSplit);
niklas.baumgarten's avatar
niklas.baumgarten committed
148

149
150
    MeanQ = (PPM->SumOnCommSplit(Mcomm * MeanQcomm, 0) / PPM->Size(commSplit)) / M;
    MeanY = (PPM->SumOnCommSplit(Mcomm * MeanYcomm, 0) / PPM->Size(commSplit)) / M;
151

152
153
    Q2 = PPM->SumOnCommSplit(Q2, 0) / PPM->Size(commSplit);
    Y2 = PPM->SumOnCommSplit(Y2, 0) / PPM->Size(commSplit);
154

niklas.baumgarten's avatar
niklas.baumgarten committed
155
156
157
158
159
160
//        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
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
  }

  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;
niklas.baumgarten's avatar
niklas.baumgarten committed
176
        }
177
178
179
180
      }
    } else {
      dMcomm = dM;
      commSplit = 0;
niklas.baumgarten's avatar
niklas.baumgarten committed
181
    }
182
183
184
185
186
187
188
189
  }

  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;
  }
niklas.baumgarten's avatar
niklas.baumgarten committed
190
191
192
};

struct SampleCounter {
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
  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);
  }

niklas.baumgarten's avatar
niklas.baumgarten committed
213
  void UpdateSampleCounter(int dM);
niklas.baumgarten's avatar
niklas.baumgarten committed
214

215
216
  SampleCounter(int M, int dM) :
    M(M), dM(dM) {}
niklas.baumgarten's avatar
niklas.baumgarten committed
217

218
219
220
221
222
  void Update() {
    M += dM;
    dMComm = 0;
    dM = 0;
  }
niklas.baumgarten's avatar
niklas.baumgarten committed
223

224
225
226
  friend Logging &operator<<(Logging &s, const SampleCounter &ctr) {
    return s << "M=" << ctr.M << " dM=" << ctr.dM << " dMComm=" << ctr.dMComm << endl;
  }
227
228
};

niklas.baumgarten's avatar
niklas.baumgarten committed
229
#endif //EMPIRICMEASURES_HPP