diff options
Diffstat (limited to 'quantity.hpp')
-rw-r--r-- | quantity.hpp | 119 |
1 files changed, 119 insertions, 0 deletions
diff --git a/quantity.hpp b/quantity.hpp new file mode 100644 index 0000000..400a0dd --- /dev/null +++ b/quantity.hpp @@ -0,0 +1,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; } +}; + |