refactor

1 files changed, 83 insertions, 0 deletions

diff --git a/quantity.hpp b/quantity.hpp
new file mode 100644
index 0000000..55d5a28
--- /dev/null
+++ b/quantity.hpp
@@ -0,0 +1,83 @@
+
+#pragma once
+
+#include <vector>
+#include <list>
+#include <cmath>
+#include <array>
+
+class Quantity {
+private:
+  double total;
+  double total2;
+  std::vector<double> C;
+  unsigned wait;
+
+public:
+  unsigned n;
+  std::list<double> hist;
+
+  Quantity(unsigned lag, unsigned wait) : C(lag), wait(wait) {
+    n = 0;
+    total = 0;
+    total2 = 0;
+  }
+
+  void add(double x) {
+    if (n > wait) {
+      hist.push_front(x);
+      if (hist.size() > C.size()) {
+        hist.pop_back();
+        unsigned t = 0;
+        for (double a : hist) {
+          C[t] += a * x;
+          t++;
+        }
+        total += x;
+        total2 += pow(x, 2);
+      }
+    }
+    n++;
+  }
+
+  double avg() const { return total / (n - wait); }
+
+  double avg2() const { return total2 / (n - wait); }
+
+  std::vector<double> ρ() const {
+    double C0 = C.front() / (n - wait);
+    double avg2 = pow(total / (n - wait), 2);
+
+    std::vector<double> ρtmp;
+
+    for (double Ct : C) {
+      ρtmp.push_back((Ct / (n - wait) - avg2) / (C0 - avg2));
+    }
+
+    return ρtmp;
+  }
+
+  std::array<double, 2> τ() const {
+    double τtmp = 0.5;
+    unsigned M = 1;
+    double c = 8.0;
+
+    std::vector<double> ρ_tmp = this->ρ();
+
+    while (c * τtmp > M && M < C.size()) {
+      τtmp += ρ_tmp[M];
+      M++;
+    }
+
+    return {τtmp, 2.0 * (2.0 * M + 1) * pow(τtmp, 2) / (n - wait)};
+  }
+
+  double σ() const {
+    return 2.0 / (n - wait) * this->τ()[0] * (C[0] / (n - wait) - pow(this->avg(), 2));
+  }
+
+  double serr() const { return sqrt(this->σ()); }
+
+  unsigned num_added() const { return n - wait; }
+};
+