blob: 400a0dddd0e79df95957f0a8d5945586d831bd35 (
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
111
112
113
114
115
116
117
118
119
|
#pragma once
#include <array>
#include <cmath>
#include <list>
#include <vector>
#include <fstream>
#include <iomanip>
template <class T>
class Quantity {
private:
uint64_t N;
uint64_t n;
unsigned skip;
std::list<T> hist;
double total;
double total2;
std::vector<double> C;
public:
Quantity(unsigned lag, unsigned s = 1) : C(lag) {
skip = s;
N = 0;
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(const T& x) {
if (N % skip == 0) {
hist.push_front(x);
if (hist.size() > C.size()) {
hist.pop_back();
unsigned t = 0;
for (T a : hist) {
C[t] += a * x;
t++;
}
double norm = x * x;
total += sqrt(norm);
total2 += norm;
n++;
}
}
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 {skip * τtmp, skip * 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; }
};
|
