#include <iostream>
#include "least_sqrs_est.h"

int main();

int main() {
	std::cout.precision( 16 );
	std::clog.precision( 16 );

	// Use this to turn logging off
	//  - comment this out if you want to see the coefficients
	// std::clog.setstate( std::ios_base::failbit );

	// The following data was generated by taking
	//    the linear polynomial      0.7 n + 5.6,
	//                                    2
	//    the quadratic polynomial  -0.1 n + 0.6 n + 1.5
	// evaluated at n = 0, ..., 9
	// and added to this a normally distributed error with
	// standard deviation sigma = 0.25 and 1.0, respectively.

	double y[2][10]{
		{ 5.7344,  6.7585,  6.4353,  7.9155,  8.4797,
		  8.7731,  9.6916, 10.5857, 12.0946, 12.5924 },
		{ 2.1715,  0.7925,  3.0172,  4.0302,  2.7889,
		  3.0347,  2.2269,  0.4966,  0.1939, -1.9873}
	};


	// The current value should therefore be estimated as
	//     6.3 + 5.6 = 11.9
	// and the value one time step into the future should be
	//     7.0 + 5.6 = 12.6

	for ( unsigned int n{0}; n < 2; ++n ) {
		for ( unsigned int k{5}; k <= 10; k += 5 ) {
			double   y9_L{  lin_y0( y[n], 9, k ) };
			double  y10_L{  lin_y1( y[n], 9, k ) };
			double  dy9_L{ lin_dy0( y[n], 9, k ) };
			double  iy9_L{ lin_iy0( y[n], 9, k ) };

			double    r_L{    lin_root( y[n], 9, k ) };
			double  jy9_L{ lin_jit( y[n], 9, k, 0.1 ) };
   
			double   y9_Q{   quad_y0( y[n], 9, k ) };
			double  y10_Q{   quad_y1( y[n], 9, k ) };
			double  dy9_Q{  quad_dy0( y[n], 9, k ) };
			double ddy9_Q{ quad_ddy0( y[n], 9, k ) };
			double  iy9_Q{  quad_iy0( y[n], 9, k ) };

			double  ext_Q{    quad_ext( y[n], 9, k ) };
			auto      r_Q{  quad_roots( y[n], 9, k ) };
			double  jy9_Q{  quad_jit( y[n], 9, k, 0.1 ) };

		 	std::cout << "With " << k << " points:" << std::endl;

		 	for ( unsigned int i{10 - k}; i < 10 - 1; ++i ) {
				std::cout << y[n][i] << ", ";
		 	}

			std::cout << y[n][9] << " = y[n]" << std::endl;

			std::cout << "Linear" << std::endl;
			std::cout << "       Approximating y[ 9]: " <<  y9_L << "\t(11.9 or -1.2)"  << std::endl;
			std::cout << "       Approximating y[10]: " << y10_L << "\t(12.6 or -2.5)"  << std::endl;
			std::cout << "            The slope at 9: " << dy9_L << "\t(0.7 or -1.2)"   << std::endl;
			std::cout << "  The integral from 8 to 9: " << iy9_L << "\t(11.55 or -0.63333)"
			  << std::endl << std::endl;

			std::cout << "Quadratic" << std::endl;
			std::cout << "       Approximating y[ 9]: " <<   y9_Q << "\t(11.9 or -1.2)"  << std::endl;
			std::cout << "       Approximating y[10]: " <<  y10_Q << "\t(12.6 or -2.5)"  << std::endl;
			std::cout << "            The slope at 9: " <<  dy9_Q << "\t(0.7 or -1.2)"   << std::endl;
			std::cout << "        The concavity at 9: " << ddy9_Q << "\t(0.0 or -0.2)"   << std::endl;
			std::cout << "  The integral from 8 to 9: " <<  iy9_Q << "\t(11.55 or -0.63333)"
				  << std::endl << std::endl;
	
			std::cout << "      A maximum or minimum: " << ext_Q << "\t(undefined or -6)" << std::endl;
			std::cout << "           A root (linear): " <<   r_L << "\t(-17 or undefined)" << std::endl;
			std::cout << "          Root (quadratic): " <<   r_Q.first << ", "
				                                    << r_Q.second << "\t(-17, undefined or -10.89898, -1.10102)"
			          << std::endl << std::endl;
		}
	}

	return 0;
}