#include <cmath>
#include <iostream>
#include <ranges>

#include "include/praktable.hpp"
#include "include/prakphys.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> &;

f64p add_nu(table& t) {
	t.add_column("nu").apply([](vecarg a) { return a[0] * a[1] / (prak::R<f64> * a[2]); }, {"P", "V", "T"}, "nu");
	return {t.col_avg("nu"), t.col_stddev("nu")};
}

void ex1(table &t1, const char *name) {
	f64 DA = prak::discrete_integral_trapezoid(t1.cbegin("V"), t1.cend("V"), t1.cbegin("P"), t1.cend("P"));
	f64 DS = DA / t1.col_avg("T");
	f64p nu = add_nu(t1);
	f64p DS_teor = nu * prak::R<f64> * std::abs(std::log(*(t1.end("V")-1) / *t1.begin("V")));
	std::cout << "Упражнение 1\n" var(DA) var(DS) var(nu) var(DS_teor) << std::endl;
	t1.write_plot(std::string("P(V)") + name + ".plot", "V", "P", std::nullopt);
	t1.write_plot(std::string("PV(V)") + name + ".plot", "V", "x", std::nullopt);
	std::ofstream f(std::string("ex1_") + name + ".print");
	f << "Данные упражнения 1";
	t1.print(f);
	f.close();
}

void ex2(const table &t) {
	
}

void ex3(const table &t) {

}

void ex4(void) {

}

void ex(void) {
	table tables[] = {table("data1"), table("data2"), table("data3"), table("data4")};

	// approximate how the volume should have changed in ex3 using parabola
	std::vector<f64> mrxdata; mrxdata.reserve(9);
	std::vector<f64> mrydata; mrydata.reserve(3);
	for (size_t i = 0; i < tables[2].rows; ++i) {
		if (f64 v; !std::isnan(v = tables[2]["V", i])) {
			f64 n = tables[2]["n", i];
			mrxdata.push_back(n*n);
			mrxdata.push_back(n);
			mrxdata.push_back(1);
			mrydata.push_back(v);
		}
	}
	prak::matrix<f64> xs(3, 3, std::move(mrxdata)), ys(3, 1, std::move(mrydata));
	prak::matrix<f64> poly = xs.inv().value() * ys;
	for (size_t i = 0; i < tables[2].rows; ++i)
		tables[2]["V", i] = (poly.tr() * prak::matrix<f64>(3, 1, {(double)(i+1)*(i+1), (double)(i+1), 1}))[0, 0];

	std::function<f64(vecarg)> fs[] = {
		[](vecarg a) { return a[0] * a[1]; },
		[](vecarg a) { return a[0] / a[1]; },
		[](vecarg a) { return a[0] / a[1]; },
		[](vecarg a) { return a[0] * a[1]; },
	};
	std::vector<std::string> vs[] = {
		{"P", "V"}, {"V", "T"}, {"P", "T"}, {"P", "V"},
	};
	for (auto elem : std::views::zip(tables, fs, vs))
		std::get<0>(elem).apply(std::get<1>(elem), std::get<2>(elem), "x");
	ex1(tables[0], "1");
	ex1(tables[3], "4");
	ex2(tables[1]);
	ex3(tables[2]);
	ex4();
}

int main() {
	ex();	
	return 0;
}