summaryrefslogtreecommitdiff
path: root/quantity.hpp
blob: 296a6e8cfe40760af7676ed6cb8c8831a2f02ce5 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110

#pragma once

#include <array>
#include <cmath>
#include <list>
#include <vector>
#include <fstream>
#include <iomanip>

class Quantity {
private:
  unsigned n;
  std::list<double> hist;
  double total;
  double total2;
  std::vector<double> C;

public:
  Quantity(unsigned lag) : C(lag) {
    n = 0;
    total = 0;
    total2 = 0;
  }

  void read(std::string filename) {
    std::ifstream inFile(filename);
    inFile >> n;
    inFile >> total;
    inFile >> total2;

    for (double& Ci : C) {
      inFile >> Ci;
    }
  }

  void write(std::string filename) const {
    std::ofstream outFile(filename);
    outFile << std::setprecision(15) << n << " " << total << " " << total2 << std::endl;
    for (double Ci : C) {
      outFile << Ci << " ";
    }
    outFile << std::endl;
    outFile.close();
  }

  void reset() {
    total = 0;
    total2 = 0;
    std::fill(C.begin(), C.end(), 0);
    n = 0;
    hist = {};
  }

  void add(double x) {
    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; }

  double avg2() const { return total2 / n; }

  std::vector<double> ρ() const {
    double C0 = C.front() / n;
    double avg2 = pow(total / n, 2);

    std::vector<double> ρtmp;

    for (double Ct : C) {
      ρtmp.push_back((Ct / n - 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};
  }

  double σ() const {
    return 2.0 / n * this->τ()[0] * (C[0] / n - pow(this->avg(), 2));
  }

  double serr() const { return sqrt(this->σ()); }

  unsigned num_added() const { return n; }
};