#include <iostream>

#include "include/praktable.hpp"

using table = prak::table<f64>;
using f64p = prak::pvalue<f64>;
using f64v = std::vector<f64>;
using vecarg = const std::vector<f64> &;
using argvec = const std::vector<f64> &;

f64 factorial(u64 times) {
	f64 res = 1;
	for (; times; res *= times--);
	return res;
}

f64 poisson(u64 at, f64 avg) {
	return std::pow(avg, at) * std::exp(-avg) / factorial(at);
}

template <typename lambda, typename...Args>
f64 distr_sum(lambda f, u64 from, u64 to, Args...args) {
	f64 res = 0;
	for(; from <= to; res += f(from++, args...));
	return res;
}

void ex1() {
	f64 n = 1.053, s = 0.582553;
	std::cout << "ex 1:\n"
		var(distr_sum(poisson, n, n, n)) 
		var(distr_sum(poisson, n, n+1, n))
		var(distr_sum(poisson, std::round(n-s), std::round(n+s), n))
		var(distr_sum(poisson, std::round(n-2*s), std::round(n+2*s), n))
	<< std::endl;
}
table ex3(const char* file) {
	table t(file);
	t	.apply([](argvec a){return a[0]*a[0];}, "s", "s2")
		.apply([](argvec a){return a[0]/a[1];}, std::vector<std::string>{"s", "n"}, "sn")
		.write_plot("s(n).plot" , "n", "s", std::nullopt)
		.write_plot("s2(n).plot", "n", "s2", std::nullopt)
		.write_plot("sn(n).plot", "n", "sn", std::nullopt)
		.print();
	return t;
}
void ex4(table t0) {
	table t1({"count", "N", "n1", "s1", "n2", "s2", "s", "s/sn1", "s/sn2", "sqrtn1", "sqrtn2"}, 6, 0);
	f64 T1 = 45e-3;
	f64 T2 = 0.625;
	for (int i = 0; i < 6; ++i) {
		t1["count", i] = i+1;
		t1["N", i] = t0["n", i] / t0["t", i] * 1000;
		t1["s", i] = t0["s", i] / std::sqrt(t0["t", i] / 1000);
	}
	t1	.apply([T1](vecarg a){return a[0]*T1;}, "N", "n1")
		.apply([T1](vecarg a){return a[0]*std::sqrt(T1);}, "s", "s1")
		.apply([T2](vecarg a){return a[0]*T2;}, "N", "n2")
		.apply([T2](vecarg a){return a[0]*std::sqrt(T2);}, "s", "s2");
	f64	avgn1 = t1.col_avg("n1"),
		avgn2 = t1.col_avg("n2"),
		savgn1 = t1.col_stddev("n1"),
		savgn2 = t1.col_stddev("n2");
	t1	.apply([savgn1](vecarg a){return a[0]/savgn1;}, "s1", "s/sn1")
		.apply([savgn2](vecarg a){return a[0]/savgn2;}, "s2", "s/sn2")
		.apply([](vecarg a){return std::sqrt(a[0]);}, "n1", "sqrtn1")
		.apply([](vecarg a){return std::sqrt(a[0]);}, "n2", "sqrtn2")
		.print();
	std::cout var(avgn1)var(avgn2)var(savgn1)var(savgn2);
}

int main() {
	ex1();
	ex4(ex3("data3"));
	return 0;
}