#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) {
	// t = const
	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);
	t1	.add_columns({"DAi", "DSi", "lnViV0"})
		.apply([&t1](vecarg a) -> f64 {
			return prak::discrete_integral_trapezoid(
				t1.begin("V"), std::find(t1.begin("V"), t1.end("V"), a[1])+1,
				t1.begin("P"), std::find(t1.begin("P"), t1.end("P"), a[0])+1
			);
		}, {"P", "V"}, "DAi")
		.apply([&t1](vecarg a){return a[0] / a[1];}, {"DAi", "T"}, "DSi")
		.apply([&t1](vecarg a){ return std::log(a[0] / t1["V", 0]); }, "V", "lnViV0")
		.apply([](vecarg a) {return std::abs(a[0]); }, "lnViV0", "lnViV0");
	f64p DS_teor = nu * prak::R<f64> * std::abs(std::log(*(t1.end("V")-1) / *t1.begin("V")));
	f64p A = std::abs(t1.least_squares_prop("lnViV0", "DSi", std::nullopt, 0.1));
	f64p R = A / nu;
	std::cout << "Упражнение 1\n" var(DA) var(DS) var(nu) var(DS_teor) << "доп данные (упражнение 4):\n" 
			var(A) var(R) << std::endl;
	t1	.write_plot(std::string("PV") + name + ".plot", "V", "P", std::nullopt)
		.write_plot(std::string("PVV") + name + ".plot", "V", "x", std::nullopt)
		.write_plot(std::string("DSilnViV0_Tconst") + name + ".plot", "lnViV0", "DSi", std::nullopt);
	std::ofstream f(std::string("ex1_") + name + ".print");
	f << "Данные упражнения 1 (T = const)\n";
	t1.print(f);
	f.close();
}

void ex2(table &t) {
	// p = const
	std::cout << "Упражнение 2\n";
	std::ofstream f("ex2.print");
	f << "Данные упражнения 2 (P = const)\n";
	f64p nu = add_nu(t);
	t	.add_columns({"lnTiT0", "DSi"})
		.apply([&t](vecarg a) -> f64 {return std::log(a[0] / t["T", 0]); }, "T", "lnTiT0")
		.apply([&t, &nu](vecarg a) -> f64 { 
			return nu.val * (prak::R<f64> * std::log(a[1] / t["V", 0]) + 5.0/2 * prak::R<f64> * a[2]); 
		}, {"T", "V", "lnTiT0"}, "DSi")
		.write_plot("DSilnTiT0_Pconst.plot", "lnTiT0", "DSi", std::nullopt)
		.write_plot("VTV.plot", "V", "x", std::nullopt)
		.print(f);
	f64p A = t.least_squares_prop("lnTiT0", "DSi", std::nullopt, 0.1);
	f64p Cp = A / nu;
	std::cout var(nu) var(A) var(Cp);
}

void ex3(table &t) {
	// v = const
	std::cout << "Упражнение 3\n";
	std::ofstream f("ex3.print");
	f << "Данные упражнения 3 (V = const)\n";
	f64p nu = add_nu(t);
	t	.add_columns({"lnTiT0", "DSi"})
		.apply([&t](vecarg a) -> f64 {return std::log(a[0] / t["T", 0]); }, "T", "lnTiT0")
		.apply([&t, &nu](vecarg a) -> f64 { 
			return nu.val * (-prak::R<f64> * std::log(a[0] / t["P", 0]) + 7.0/2 * prak::R<f64> * a[1]); 
		}, {"P", "lnTiT0"}, "DSi")
		.write_plot("DSilnTiT0_Vconst.plot", "lnTiT0", "DSi", std::nullopt)
		.write_plot("PTP.plot", "P", "x", std::nullopt)
		.print(f);
	auto [A, B] = t.least_squares_linear("lnTiT0", "DSi", std::nullopt, 0.1);
	f64p Cv = A / nu;
	std::cout var(nu) var(A) var(Cv);

}

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]);
}

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