From 89d26d5d0e8ff2b026cc99606868699f6fcc08db Mon Sep 17 00:00:00 2001 From: justanothercatgirl Date: Mon, 18 Nov 2024 23:40:21 +0300 Subject: added tempalte --- libprakpp/include/prakcommon.hpp | 98 +++++++++++ libprakpp/include/prakmath.hpp | 319 ++++++++++++++++++++++++++++++++++++ libprakpp/include/prakmatrix.hpp | 185 +++++++++++++++++++++ libprakpp/include/praktable.hpp | 345 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 947 insertions(+) create mode 100644 libprakpp/include/prakcommon.hpp create mode 100644 libprakpp/include/prakmath.hpp create mode 100644 libprakpp/include/prakmatrix.hpp create mode 100644 libprakpp/include/praktable.hpp (limited to 'libprakpp/include') diff --git a/libprakpp/include/prakcommon.hpp b/libprakpp/include/prakcommon.hpp new file mode 100644 index 0000000..749d5ef --- /dev/null +++ b/libprakpp/include/prakcommon.hpp @@ -0,0 +1,98 @@ +#pragma once + +#include +#include +#include +#include + +#if defined(_MSC_VER) || !defined(__cpp_multidimensional_subscript) || __cplusplus < 202110L +#warning "can not use multidimentional subscript operator: falling back to `operator()`" +#undef MDSUBSCRIPT +#else +#define MDSUBSCRIPT +#endif + +#ifdef NO_MDSUBSCRIPT +#undef MDSUBSCRIPT +#endif + +#ifdef MDSUBSCRIPT + #define SUBSCR_OPN [ + #define SUBSCR_CLS ] + #define SUBSCR_OPRTR operator[] +#else + #define SUBSCR_OPN ( + #define SUBSCR_CLS ) + #define SUBSCR_OPRTR operator() +#endif + +template +using function_t = std::function &)>; + +namespace prak { + +/// stolen from [cppreference.com](https://en.cppreference.com/w/cpp/types/numeric_limits/epsilon) +/// Compares 2 floating-point values up to `ulps` ULPS (units in the last place) +template +std::enable_if_t::is_integer, bool> +fequal(T x, T y, std::size_t ulps = 1) +{ + // Since `epsilon()` is the gap size (ULP, unit in the last place) + // of floating-point numbers in interval [1, 2), we can scale it to + // the gap size in interval [2^e, 2^{e+1}), where `e` is the exponent + // of `x` and `y`. + + // If `x` and `y` have different gap sizes (which means they have + // different exponents), we take the smaller one. Taking the bigger + // one is also reasonable, I guess. + const T m = std::min(std::fabs(x), std::fabs(y)); + + // Subnormal numbers have fixed exponent, which is `min_exponent - 1`. + const int exp = m < std::numeric_limits::min() + ? std::numeric_limits::min_exponent - 1 + : std::ilogb(m); + + // We consider `x` and `y` equal if the difference between them is + // within `n` ULPs. + return std::fabs(x - y) <= ulps * std::ldexp(std::numeric_limits::epsilon(), exp); +} + + +/// prints a vector +template +void printv(const T &vec) { + std::cout << "std::vector { "; + for (const auto &x : vec) { + std::cout << x << ' '; + } + std::cout << '}' << std::endl; +} + +/// An allocator that aligns memory to 64 bytes. Needed for AVX instructions. +/// C++17 required for std::align_val_t +template +struct align_alloc : public std::allocator { + constexpr T* allocate( std::size_t n ) { + return static_cast(::operator new(sizeof (T) * n, std::align_val_t{64})); + } + + constexpr void deallocate( T* p, std::size_t n ) { + ::operator delete(p, std::align_val_t{64}); + } +};; + +/// alias prak::vector that is the same as std::vector, but uses aligned allocator +template +using vector = std::vector>; + +/// prak value / pair value: a value with an error +template +struct pvalue { T val, err; }; + +template +std::ostream &operator<<(std::ostream &os, const struct pvalue &p) { + /* return os << "value {" << p.val << "±" << p.err << "}"; */ + return os << p.val << "±" << p.err; +} + +} // namespace prak diff --git a/libprakpp/include/prakmath.hpp b/libprakpp/include/prakmath.hpp new file mode 100644 index 0000000..3a4e4db --- /dev/null +++ b/libprakpp/include/prakmath.hpp @@ -0,0 +1,319 @@ +#pragma once + +#include +#include + +#ifdef __x86_64__ +#include +#endif + +#include "prakcommon.hpp" +#include "prakmatrix.hpp" + +namespace prak { + +// :) +constexpr const double PI = 3.141592657; +constexpr const float fPI = 3.141592657f; + +/// defines a type that supports arithmetic operation +template +concept arithmetic = requires(T a) { + a + a; + a * a; + a - a; + a / a; +}; + +/// TODO: remove +enum struct operation { mul, div, add, sub }; + +/// vector multiply: fallback for non-floating-point types +template +void vmul(const T *op1, const T *op2, T *dest, size_t s) { + for (size_t i = 0; i < s; ++i) { + dest[i] = op1[i] * op2[i]; break; + } +} + +/// vector multiply: float implementation +template <> +inline void vmul(const float *op1, const float *op2, float *dest, size_t s) { + if (s < 8) goto scalar; + __m256 b1; + __m256 b2; + + for (size_t i = 0; i < s / 8; ++i) { + b1 = _mm256_load_ps(op1 + 8*i); + b2 = _mm256_load_ps(op2 + 8*i); + b1 = _mm256_mul_ps(b1, b2); + _mm256_store_ps(dest + 8*i, b1); + } +scalar: + for (size_t i = s - s % 8; i < s; ++i) { + dest[i] = op1[i] * op2[i]; + } +} + +/// vector multiply: double implementation +template <> +inline void vmul(const double *op1, const double *op2, double *dest, size_t s) { + if (s < 4) goto scalar; + __m256d b1; + __m256d b2; + + for (size_t i = 0; i < s / 4; ++i) { + b1 = _mm256_load_pd(op1 + 4*i); + b2 = _mm256_load_pd(op2 + 4*i); + b1 = _mm256_mul_pd(b1, b2); + _mm256_store_pd(dest + 4*i, b1); + } +scalar: + for (size_t i = s - s % 4; i < s; ++i) { + dest[i] = op1[i] * op2[i]; + } +} + +/// vector division: non-floating-point types fallback +template +void vdiv(const T *op1, const T *op2, T *dest, size_t s) { + for (size_t i = 0; i < s; ++i) { + dest[i] = op1[i] / op2[i]; break; + } +} + +/// vector division: floating point: single precision +template <> +inline void vdiv(const float *op1, const float *op2, float *dest, size_t s) { + if (s < 8) goto scalar; + __m256 b1; + __m256 b2; + + for (size_t i = 0; i < s / 8; ++i) { + b1 = _mm256_load_ps(op1 + 8*i); + b2 = _mm256_load_ps(op2 + 8*i); + b1 = _mm256_div_ps(b1, b2); + _mm256_store_ps(dest + 8*i, b1); + } +scalar: + for (size_t i = s - s % 8; i < s; ++i) { + dest[i] = op1[i] / op2[i]; + } +} + +/// vector division: floating point: double precision +template <> +inline void vdiv(const double *op1, const double *op2, double *dest, size_t s) { + if (s < 4) goto scalar; + __m256d b1; + __m256d b2; + + for (size_t i = 0; i < s / 4; ++i) { + b1 = _mm256_load_pd(op1 + 4*i); + b2 = _mm256_load_pd(op2 + 4*i); + b1 = _mm256_div_pd(b1, b2); + _mm256_store_pd(dest + 4*i, b1); + } +scalar: + for (size_t i = s - s % 4; i < s; ++i) { + dest[i] = op1[i] / op2[i]; + } +} + +inline float finalize(__m256 reg) { + // reg: x0 | x1 | x2 | x3 | x4 | x5 | x6 | x7 + __m128 extr = _mm256_extractf128_ps(reg, 1); + // extr: x0 | x1 | x2 | x3 + // __mm256_castps256_ps128: x4 | x5 | x6 | x7 + extr = _mm_add_ps(extr, _mm256_castps256_ps128(reg)); + // extr: x0+x4 | x1 + x5 | x2 + x6 | x3 + x7 + extr = _mm_hadd_ps(extr, extr); + // extr: x0 + x1 + x2 + x5 | x2 + x3 + x6 x7 | ... | ... + return extr[0] + extr[1]; +} + +inline double finalize(__m256d reg) { + // reg: x0 | x1 | x2 | x3 + reg = _mm256_hadd_pd(reg, reg); + // reg: x0 + x1 | x2 + x3 | ... | ... + return reg[0] + reg[1]; +} + +/// Fallback generic method to sum a vector +template +T vsum(const T *op, size_t size) { + T res; + for (size_t i = 0; i < size; ++i) + res += op[i]; + return res; +} + +/// Sum a vector, float implementation +template <> +inline float vsum(const float *op, size_t size) { + __m256 buf; + buf = _mm256_setzero_ps(); + if (size < 8) goto scalar; + for (size_t i = 0; i < size / 8; ++i) { + __m256 next = _mm256_load_ps(op + 8*i); + buf = _mm256_add_ps(buf, next); + } +scalar: + float linear = 0; + for (size_t i = size - size%8; i < size; ++i) { + linear += op[i]; + } + return linear + finalize(buf); +} + +/// Sum a vector, double implementation +template <> +inline double vsum(const double *op, size_t size) { + __m256d buf; + buf = _mm256_setzero_pd(); + if (size < 4) goto scalar; + for (size_t i = 0; i < size / 4; ++i) { + __m256d next = _mm256_load_pd(op + 4*i); + buf = _mm256_add_pd(buf, next); + } +scalar: + double linear = 0; + for (size_t i = size - size%8; i < size; ++i) { + linear += op[i]; + } + return linear + finalize(buf); +} + +// take a derivative of function `f` with arguments `at` and with differentiable variable being at index `varidx` +template +T deriv4(function_t f, std::vector /*I actually do want a copy here*/ at, size_t varidx) { + T h = 1e-5; + T sum = 0; + + at[varidx] += 2*h; + sum -= f(at); + at[varidx] -= h; + sum += 8 * f(at); + at[varidx] -= 2*h; + sum -= 8*f(at); + at[varidx] -= h; + sum += f(at); + + return sum / (12*h); +} + +/// get error of the calculated value +template +T sigma(function_t f, const std::vector &args, const std::vector sigmas) noexcept(false) { + if (args.size() != sigmas.size()) throw dimension_error("function prak::sigma : Args.size() does not match sigmas.size()"); + T sum = 0; + for (size_t i = 0; i < args.size(); ++i) { + T tmp = deriv4(f, args, i) * sigmas[i]; + sum += tmp*tmp; + } + return std::sqrt(sum); +} + +/// calculate least-squares linear approximation to fit data +/// ax+b = y (ss is an error of Y value) +template +void least_squares_linear( + const prak::vector &xs, + const prak::vector &ys, + const prak::vector &ss, + struct pvalue &a, + struct pvalue &b) +{ + if (xs.size() != ys.size() || xs.size() != ss.size()) { + std::cout << "x.size() = " << xs.size() + << ", y.size() = " << ys.size() + << ", s.size() = " << ss.size() << std::endl; + return; + } + [[assume(xs.size() == ys.size() && ys.size() == ss.size())]]; + size_t sz = xs.size(); + + prak::vector ones(sz); + prak::vector ssq(sz); + prak::vector ssq_m1(sz); + prak::vector buf(sz); + prak::vector buf_xs_ssq_m1(sz); + std::fill(ones.begin(), ones.end(), (T)1); // ones: [1] + vmul(ss.data(), ss.data(), ssq.data(), sz); // ssq: [ss*ss] + vdiv(ones.data(), ssq.data(), ssq_m1.data(), sz); // ssq_m1: [1/ss^2] + vmul(xs.data(), ssq_m1.data(), buf_xs_ssq_m1.data(), sz); // [xs / ss^2] + + const T *ysd = ys.data(), + *xsd = xs.data(); + T *ssqd = ssq.data(), + *ssq_m1d = ssq_m1.data(), + *onesd = ones.data(), + *bufd = buf.data(), + *buf_xs_ssq_m1d = buf_xs_ssq_m1.data(); + + vmul(buf_xs_ssq_m1d, xsd, bufd, sz); // buf: [xs^2 / ss^2] + T d1 = vsum(bufd, sz); // sum([xs^2 / ss^2]) + T ssq_m1sum = vsum(ssq_m1d, sz); // sum(1 / ss^2) + vmul(xsd, ssq_m1d, bufd, sz); // buf: [xs / ss^2] + T d2 = vsum(buf_xs_ssq_m1d, sz); // sum([xs / ss^2]) + T D = d1 * ssq_m1sum - d2*d2; // sum((xs/ss)^2) * sum(1/ss^2) - sum(xs/ss^2)^2 + + vmul(ysd, buf_xs_ssq_m1d, bufd, sz); // buf: [ys*xs/ss^2] + T da1 = vsum(bufd, sz); // sum([ys*xs/ss^2]) + vmul(ysd, ssq_m1d, bufd, sz); // buf: [ys/ss^2] + T da2 = vsum(bufd, sz); // sum([ys/ss^2]) + T DA = da1 * ssq_m1sum - da2 * d2; // sum([ys*xs/ss^2]) * sum([1/ss^2]) - sum([ys/ss^2]) * sum(xs/ss^2) + + T DB = d1 * da2 - d2 * da1; // sum([xs^2/ss^2]) * sum([ys/ss^2]) - sum([xs/ss^2]) * sum(ys*xs/ss^2) + + a.val = DA/D; + a.err = sqrt(ssq_m1sum / D); + + b.val = DB/D; + b.err = sqrt(d1 / D); +} + + +/// May throw std::bad_optional_access +template +std::enable_if, std::vector>> +polynomial_regression( + size_t degree, + std::vector data_x, + std::vector data_y, + std::optional> data_errors = std::nullopt) +{ + ++degree; // hack) + size_t data_size = data_x.size(); + if (data_size != data_y.size() || (data_errors.has_value() && data_errors->size() != data_size)) + throw dimension_error("Xs, Ys or Sigmas do not match sizes"); + struct matrix X(data_size, degree), + Y(data_size, 1), + B(degree, 1), + W = matrix::identity(data_size); + // initialize X + for (size_t row = 0; row < X.rows; ++row) { + X.SUBSCR_OPRTR(row, 0) = 1; + for (size_t col = 1; col < X.cols; ++col) { + X.SUBSCR_OPRTR(row, col) = X.SUBSCR_OPRTR(row, col-1) * data_x[row]; + } + } + // initialize Y + std::memcpy(Y.data(), data_y.data(), sizeof (T) * data_size); + // initialize W + if (data_errors.has_value()) { + std::vector &err_value = *data_errors; + for (size_t i = 0; i < err_value.size(); ++i) + W.data()[i * (data_size + 1)] = 1 / err_value[i] / err_value[i]; + } + std::cerr << X << '\n' << Y << '\n' << W << '\n'; + matrix X_T_W = X.tr() * W; + matrix tmp1 = (X_T_W * X).inv().value(); + B = tmp1 * X_T_W * Y; + B.print(); + // TODO: FINISH (with covariation matrix) + return {}; +} + + +} // namespace prak diff --git a/libprakpp/include/prakmatrix.hpp b/libprakpp/include/prakmatrix.hpp new file mode 100644 index 0000000..87b38b4 --- /dev/null +++ b/libprakpp/include/prakmatrix.hpp @@ -0,0 +1,185 @@ +#pragma once + +#include "prakcommon.hpp" + +#include +#include +#include +#include + +namespace prak { + +struct dimension_error : public std::logic_error { + explicit dimension_error(const char *what): std::logic_error(what) {} +}; + +/// A class to represent a matrix. Works as an extension to std::vector, so all methods of std::vector are available on this class +/// BIG TODO: replave `throw` with `std::optional` return type +template +struct matrix : public std::vector { + size_t rows, cols; + size_t feq_precision = 10; + const T &SUBSCR_OPRTR(size_t row, size_t col) const { + return this->at(row * cols + col); + } + T &SUBSCR_OPRTR(size_t row, size_t col) { + return this->at(row * cols + col); + } + matrix(void) = default; + matrix(struct matrix &&) = default; + matrix(const struct matrix &) = default; + struct matrix &operator=(const struct matrix &) = default; + explicit matrix(size_t rows, size_t cols, T def) + : rows{rows}, cols{cols}, std::vector(rows * cols, def) {} + explicit matrix(size_t rows, size_t cols, T *data = nullptr) + : rows{rows}, cols{cols}, std::vector(rows * cols) { + if (data != nullptr) std::memcpy(this->data(), data, rows * cols); + } + explicit matrix (size_t rows, size_t cols, std::vector &&data) noexcept(false) + : rows{rows}, cols{cols}, std::vector(std::move(data)) { + if (this->size() != rows * cols) throw std::out_of_range("data.size() != rows*cols"); + } + /// add 2 matrixes + struct matrix operator+(const struct matrix &other) noexcept(false) { + if (other.rows != rows || other.cols != cols) + throw dimension_error("matrix add: dimensions of 'other' do not match dimensions of 'this'"); + struct matrix ret(rows, cols); + for (size_t i = 0; i < this->size(); ++i) { + ret.data()[i] = this->data()[i] + other.data()[i]; + } + return ret; + } + /// perform matrix multiplication + struct matrix operator*(const struct matrix &other) noexcept(false) { + if (cols != other.rows) + throw dimension_error("Matrix mul: 'other.rows' does not match 'this.cols'"); + struct matrix ret(rows, other.cols, T{}); + for (size_t i = 0; i < rows; ++i) { + for (size_t j = 0; j < other.cols; ++j) { + for (size_t k = 0; k < cols; ++k) + ret.SUBSCR_OPRTR(i, j) += SUBSCR_OPRTR(i, k) * other.SUBSCR_OPRTR(k, j); + } + } + return ret; + } + /// multiply matrix by a vector [other.size() x 1] + std::enable_if, std::vector> + operator*(const std::vector &other) noexcept(false) { + if (cols != other.size()) throw dimension_error("Matrix mul: 'other.rows' does not match 'this.cols'"); + std::vector ret(rows, T{}); + for (size_t row = 0; row < rows; ++rows) { + for (size_t k = 0; k < cols; ++k) + ret[row] += SUBSCR_OPRTR(row, k) * other[k]; + } + return ret; + } + /// multiply matrix by a scalar + template + std::enable_if_t, struct matrix> + operator*(const Arg other) const { + struct matrix ret(rows, cols); + for (size_t i = 0; i < this->size(); ++i) + ret.data()[i] = this->data()[i] * other; + return ret; + } + /// Tries to find an inverse of the matrix by using gaussian elimination + /// TODO: handle the case where the inverse does not exist + std::enable_if_t, std::optional>> + inv(void) const noexcept(false) { + if (rows != cols) throw dimension_error("Matrix must be square to be inverted"); + struct matrix gaussian(rows, cols * 2, T{}); + for (size_t row = 0; row < rows; ++row) { + std::memcpy(gaussian.data() + (row * cols * 2 ), this->data() + row * cols, cols * sizeof (T)); + gaussian.data()[row * cols * 2 + cols + row] = 1; + } + for (size_t row = 0; row < rows; ++row) { + T &fst_row_v = gaussian.SUBSCR_OPRTR(row, row); + bool row_found = true; + if (fequal(fst_row_v, 0, feq_precision)) row_found = gaussian.find_suitable_row(row); + if (!row_found) return std::nullopt; + gaussian.mul_row(row, 1/gaussian.SUBSCR_OPRTR(row, row)); + for (size_t pr_row = 0; pr_row < rows; ++pr_row) { + if (row == pr_row) continue; + T& row_n_v = gaussian.SUBSCR_OPRTR(pr_row, row); + if (fequal(row_n_v, 0, feq_precision)) continue; + gaussian.muladd_row(row, pr_row, -row_n_v/fst_row_v); + } + } + struct matrix ret(rows, cols); + for (size_t row = 0; row < rows; ++row) { + std::memcpy(ret.data() + row * cols, gaussian.data() + row * 2*cols + cols, cols * sizeof (T)); + } + return ret; + + } + struct matrix tr(void) const { + struct matrix ret(cols, rows); + for (size_t row = 0; row < rows; ++row) + for (size_t col = 0; col < cols; ++col) + ret.SUBSCR_OPRTR(col, row) = SUBSCR_OPRTR(row, col); + return ret; + } + /// swaps rows `r1`, `r2` + void swap_rows(size_t r1, size_t r2) { + T* tmpbuf = new T[cols]; + size_t rowsize = cols * sizeof (T); + std::memcpy(tmpbuf, this->data() + cols * r1, rowsize); + std::memcpy(this->data() + cols * r1, this->data() + cols * r2, rowsize); + std::memcpy(this->data() + cols * r2, tmpbuf, rowsize); + delete[] tmpbuf; + } + /// multiplies `r_src` row by number `mul` and adds it ro row `r_dest` + void muladd_row(size_t r_src, size_t r_dest, T mul) { + for (size_t col = 0; col < cols; ++col) { + size_t r_src_ofs = cols * r_src + col; + size_t r_dest_ofs = cols * r_dest + col; + this->data()[r_dest_ofs] += this->data()[r_src_ofs] * mul; + } + } + /// multiplies a row by a number + void mul_row(size_t r1, T mul) { + for (size_t ind = r1 * cols; ind < (r1 + 1) * cols; ++ind) + this->data()[ind] *= mul; + } + /// finds a row that contains non-zero number at index [row = i > `row`, column=`row`]. + /// will only search for the rows below `row` and swap them + /// returns whether the swap occured (whether the row was found) + bool find_suitable_row(size_t row) { + for (size_t i = row + 1; i < rows; ++i) { + if (SUBSCR_OPRTR(i, row) != 0) { + swap_rows(i, row); + return true; + } + } + return false; + } + + static struct matrix identity(size_t rank) { + struct matrix ret(rank, rank); + for (size_t i = 0; i < rank * rank; i += (rank + 1)) { + ret.data()[i] = 1; + } + return ret; + } + + std::ostream &print(std::ostream &out = std::cout, size_t width = 10) const { + for (size_t row = 0; row < rows; ++row) { + for (size_t col = 0; col < cols; ++col) + out << "|" << std::setw(width) << SUBSCR_OPRTR(row, col); + std::cout << "|\n"; + } + return out; + } +}; + +template +void operator/(const int x, const struct matrix &m) { + std::cout << __PRETTY_FUNCTION__ << std::endl; +} + +template +std::ostream &operator<<(std::ostream &out, const struct matrix &m) { + return m.print(out); +} + +} // namespace prak diff --git a/libprakpp/include/praktable.hpp b/libprakpp/include/praktable.hpp new file mode 100644 index 0000000..f32d00d --- /dev/null +++ b/libprakpp/include/praktable.hpp @@ -0,0 +1,345 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +#include "prakcommon.hpp" +#include "prakmath.hpp" + +namespace prak { + +/// truncate a string to n characters, or return in unmodified +inline std::string truncate_or(const std::string &str, size_t n) { + if (str.size() >= n) return str.substr(0, n); + return str; +} + + +/// Unused for now. +/// TODO: Remove +template +struct opt_value { + enum struct t: unsigned char { + val = 0, + up = 1 << 0, + left = 1 << 1, + down = 1 << 2, + right = 1 << 3, + } tag; + T val; + + explicit opt_value(const T &v): tag{t::val}, val{v} {} + explicit opt_value(const opt_value::t t): tag{t} {} + opt_value &operator=(const T &v) { + val = v; + tag = t::val; + return *this; + } + opt_value &operator=(const opt_value::t t) { + tag = t; + return *this; + } + +}; + + +/// Class that can store, print, and apply operations to a table used by lab works +template +class table { + std::vector data; + + using function_type = function_t; + using stringvec = std::vector; + + size_t index(const std::string &str) const { + return std::distance(names.cbegin(), std::find(names.cbegin(), names.cend(), str)); + } + + FILE* open_gnuplot() { +#if defined(_POSIX_VERSION) || defined(__unix) + return popen("gnuplot", "w"); +#else +#warning Not implemented for non-unix operating systems + return NULL; +#endif + } + +public: + size_t rows = 0, columns = 0; + friend class iterator; + struct iterator { + table *parent; + size_t columns; + size_t col_index = 0, data_index = 0; + using iterator_category = std::forward_iterator_tag; + using value_type = dtype; + using difference_type = ptrdiff_t; + using pointer = dtype *; + using reference = dtype &; + iterator(table *new_parent, const std::string &column, size_t row_idx = 0) + : parent{new_parent}, columns(new_parent->columns) + { + col_index = parent->index(column); + data_index = col_index + row_idx * new_parent->columns; + } + iterator &operator++() { data_index += columns; return *this; } + iterator operator++(int) { iterator ret = *this; ++(*this); return ret; } + bool operator==(iterator other) { return data_index == other.data_index && parent == other.parent && col_index == other.col_index; } + bool operator!=(iterator other) { return data_index != other.data_index || parent != other.parent || col_index != other.col_index; } + value_type operator*() { return parent->data[data_index]; }; + }; + /// Optional rownames: names of rows + std::vector opt_rownames; + /// Mandatory columnnames: names of columns + std::vector names; + /// width used for printing, defaults to 8 + size_t column_width = 12; + + /// default constructor + table() = default; + + /// Data: array of ararys, not freed automatically. + explicit table(const std::vector &strings, dtype **data, size_t rows /*size_t columns = strings.size()*/) + : rows{rows}, columns{strings.size()} + { + names = strings; + + data = std::vector(rows * strings.size()); + for (size_t i = 0; i < rows; ++i) + for (size_t j = 0; j < strings.size(); ++j) + data[i * strings.size() + j] = data[i][j]; + } + + /// Strings: names for columns + /// `new_data` format: { { entry1_a, entry2_a, ...} { entry1_b, entry2_b, ... }, ... } + /// where `a`, `b`, ... are columns + explicit table(std::vector &&strings, std::vector> &&new_data) + : rows{new_data.size() ? new_data[0].size() : 0}, columns{strings.size()} + { + names = strings; + data = std::vector(new_data.size() * rows); + auto dit = data.begin(); + for (size_t j = 0; j < rows; ++j) { + for (size_t i = 0; i < new_data.size(); ++i) { + *dit = new_data[i][j]; + ++dit; + } + } + } + + iterator begin(std::string column) { return iterator(this, column); } + iterator end(std::string column) { return iterator(this, column, rows); } + + dtype & SUBSCR_OPRTR (const std::string &column, size_t row) { + return data.at(names.size() * row + index(column)); + } + + dtype & SUBSCR_OPRTR (size_t column, size_t row) { + return data.at(names.size() * row + column); + } + + /// prints a table. defaults to using std::cout, but any std::ostream can be passed in it. + void print(std::ostream &stream = std::cout) const { + stream << "columns: " << columns << ", rows: " << rows << std::endl; + bool print_names = opt_rownames.size() == data.size() / columns; + size_t topchars = (column_width + 1) * (columns + print_names) + 1; + + std::string rowsep(column_width, '_'); + for (size_t i = 0; i < columns + print_names - 1; ++i) { + rowsep.append(1, '|'); + rowsep.append(column_width, '_'); + } + + stream << std::left << std::string(topchars, '_') << std::endl; + if (print_names) stream << '|' << std::setw(column_width) << ' '; + + for (const auto &s : names) + stream << '|' << std::setw(column_width) << truncate_or(s, column_width-1); + + for (size_t i = 0; auto x : data) { + if (i % columns == 0) { + stream << '|' << std::endl; + stream << '|' << rowsep << '|' << std::endl; + if (print_names) + stream << '|' << std::setw(column_width) << truncate_or(opt_rownames[i / columns], column_width); + } + stream << '|' << std::setw(column_width) << x; + ++i; + } + stream << '|' << std::endl << '|' << rowsep << '|' << std::endl; + } + + /// Returns whether the amount of names is correct + /// If it is incorrect, names won't be displayed during printing + bool set_rownames(std::vector &&names) { + opt_rownames = names; + return opt_rownames.size() == data.size() / names.size(); + } + + /// apply a function to several columns and store result in another column + /// function must accept std::vector or arguments + table &apply(function_type function, stringvec args, std::optional result) { + size_t result_index = result.has_value() ? index(*result) : 0; + for (size_t i = 0; i < rows; ++i) { + std::vector v(args.size()); + for (size_t j = 0; j < args.size(); ++j) + v[j] = SUBSCR_OPRTR(args[j], i); + if (result.has_value()) data[columns * i + result_index] = function(v); + else (void)function(v); + } + return *this; + } + + /// adds a column with name `name` and data `column_data` + void add_column(std::string name, std::vector column_data) { + std::vector data_new(rows * (++columns)); + + for (size_t row = 0; row < rows; ++row) { + for (size_t column = 0; column < names.size(); ++column) // columns variable is incremented already + data_new[row * columns + column] = data[row * names.size() + column]; + data_new[(row + 1) * columns - 1] = column_data[row]; + } + data = std::move(data_new); + names.push_back(name); + } + + /// Appends a column to the table. if name is set, appends it to `opt_rownames` + void add_row(std::vector values, std::optional name = std::nullopt) { + data.resize(columns * (++rows)); + std::copy_n(values.cbegin(), columns, data.end() - columns); + if (name.has_value()) opt_rownames.push_back(*name); + } + + friend std::ostream& operator<<(std::ostream &os, table &t) { + t.print(os); + return os; + } + + /// Reads a table from a file in a format: + /// ``` + /// col1 col2 col3 ... + /// val1 val2 val3 ... + /// val4 val5 val6 ... + /// ... + /// ``` + /// Note tha `val` may either be a real number or a question mark, denoting that the value is unknown + /// `col` may be any string without whitespaeces. + /// if the first column is named "__name__" (as in python), first val in each row is a string used as + /// a row name. + void read(std::ifstream& f) { + std::string header; + std::getline(f >> std::ws, header); + std::istringstream h_stream(header); + std::string buffer; + + bool read_names = false; + h_stream >> buffer; + if (buffer == "__name__") read_names = true; + else names.push_back(buffer); + + for (size_t i = read_names ? 0 : 1; h_stream >> buffer; ++i) + names.push_back(buffer); + + std::vector tmp_row(names.size()); + while (!(f >> std::ws).eof()) { + if (read_names) { + f >> buffer; + opt_rownames.push_back(buffer); + } + for (auto it = tmp_row.begin(); it != tmp_row.end(); ++it) { + if ((f >> std::ws).peek() == '?') { + *it = NAN; + f >> buffer; + } + else f >> *it; + } + data.resize(data.size() + names.size()); + std::copy_n(tmp_row.cbegin(), names.size(), data.end() - names.size()); + ++rows; + } + columns = names.size(); + } + + /// returns an std::pair with coefficients A and B in that order + std::pair, prak::pvalue> least_squares_linear(std::string x, std::string y, std::optional sigma = std::nullopt) { + prak::vector _x(rows); + prak::vector _y(rows); + prak::vector _s(rows); + + std::copy(begin(x), end(x), _x.begin()); + std::copy(begin(y), end(y), _y.begin()); + if (sigma.has_value()) std::copy(begin(*sigma), end(*sigma), _s.begin()); + else _s = prak::vector(rows, static_cast(1)); + + std::pair, prak::pvalue> ret; + prak::least_squares_linear(_x, _y, _s, ret.first, ret.second); + return ret; + } + + /// Serialize data in format `data[args[0]][i] data[args[1]][i] data[args[2]][i]...` + void print_plot(const stringvec &args, std::ostream &out = std::cout) const { + std::vector offsets(args.size()); + for (size_t i = 0; i < args.size(); ++i) { + offsets[i] = index(args[i]); + } + for (size_t row = 0; row < rows; ++row) { + size_t row_offset = columns * row; + for (const auto column_offset : offsets) + std::cout << data[row_offset + column_offset] << ' '; + std::cout << std::endl; + } + } + + /// Serialize data in format `data[xs][i] data[ys][i] `, readable by gnuplot with yerrorbars + void write_plot(const std::string &xs, const std::string &ys, std::optional yss = std::nullopt, std::ostream &out = std::cout) const { + size_t nosigma = std::numeric_limits::max(); + size_t xsi = index(xs), ysi = index(ys), ssi = nosigma; + if (yss.has_value()) ssi = index(*yss); + for (size_t row = 0; row < rows; ++row) { + size_t offset = columns * row; + out << data[offset + xsi] << ' ' << data[offset + ysi]; + if (ssi != nosigma) out << ' ' << data[offset+ssi]; + out << std::endl; + } + } + + void plot_png( + const std::string output_filename, + const std::string &x, + const std::string &y, + const std::optional &xlabel = std::nullopt, + const std::optional &ylabel = std::nullopt, + const std::optional &sigma = std::nullopt) { + // TODO: Finish + } + + struct plot { + std::string x; + std::string y; + std::optional sigma; + bool plot_points = true; + std::optional label; + }; + + void plot_png(const std::string output, const std::vector &plots) { + // TODO: Finish later + /* FILE *gnuplot = open_gnuplot(); */ + /* fprintf(gnuplot, */ + /* "set terminal pngcairo enhanced size 800,600 dpi 300\n" */ + /* "set output '%s'\n" */ + /* , output.data() */ + /* ); */ + } +}; + +} // namespace prak -- cgit v1.2.3-70-g09d2