#include "reverseTask_solver.h"



#include <iomanip>

#include <fstream>



void swap(double* link1, double* link2)

{

	double *temp;

	temp = link1;

	link1 = link2;

	link2 = temp;

}



void vcopy(double* _output, double* _input, int size)

{

	for (int i = 0; i < size; ++i)

	{

		_output[i] = _input[i];

	}

}





CReverseTask_solver::CReverseTask_solver()

{

	iter_num = 0;

	cur_q = NULL;

	prev_q = NULL;

	curGr = NULL;

	prevGr = NULL;

	curS = NULL;

	prevS = NULL;

	accuracy = 0.01; 

}



CReverseTask_solver::~CReverseTask_solver()

{

	if (cur_q != NULL)

		delete [] cur_q;

}



void CReverseTask_solver::solve_task()

{

	stepT = maxT / K;



	//??? phi1 ??????????? ???????? q ? ?? ?????? ?????!! (q - ??????? ???????? ???????)

	CParabolicSolver* direct_solver = new CParabolicSolver();

	direct_solver->setGrid(N, K, xb, xe, maxT);

	direct_solver->setTask(a, b, c, f, alpha1, betta1, alpha2, betta2, phi1, phi2, psi);

	direct_solver->solve_implicit();

	direct_solver->CheckSolve("direct_solve.txt");

	directTable = direct_solver->GetTable();

	analytic = direct_solver->getLeftBorder();



	CParabolicVector_solver* directVector_solver = new CParabolicVector_solver();

	directVAq_solver = new CParabolicVector_solver();

	directVAs_solver = new CParabolicVector_solver();

	CBackTimeSolve* backTime_solver = new CBackTimeSolve();



	cur_q = new double[K+1];

	prev_q = new double[K+1];

	curGr = new double[K+1];

	prevGr = new double[K+1];

	curS = new double[K+1];

	prevS = new double[K+1];



	// ????????? ??????????? ??????? q

	for (int i = 0; i < K + 1; ++i)

	{

		cur_q[i] = 1.0;

	}



	// ?????? ??? ?????

		//????????? cur_q 



	// ?????? ??? ???



	//?????????? ?????????



	//curGr = ... ???????? ????????



	directVector_solver->setGrid(N, K, xb, xe, maxT);

	//????????? ????????? (????? ?????...)

	directVector_solver->setTask(a, b, c, f, task_lambda, 0, 0, 1, cur_q, phi2, psi); // border - ??????, ??????? ???? ???????? ? ?????? ???????? ?????????

//	directVector_solver->setTask(a, b, c, f, task_lambda, 0, 1, 0, cur_q, phi2, psi); // border - ??????, ??????? ???? ???????? ? ?????? ???????? ?????????

	directVector_solver->solve_implicit();



	print_table(directVector_solver->getTable(), "directVectorGradient.txt");



	double* right_border = directVector_solver->getRightBorder();



	//????????? ????????? (????? ?????...) ???? ??????????? ??? ??????? ? ?????? ???????.

	backTime_solver->setTask(a, b, c, task_lambda, 0, task_lambda, 0, xb, xe, maxT, N, K, right_border, analytic);

	backTime_solver->getSolve();

	//print_table(backTime_solver->getTable(), "backTimeGradient.txt");

	backTime_solver->getGradient(curGr);



	iter_gamma = 0.0;



	vcopy(curS, curGr, K + 1);	//curS = curGr;	// ????????? ???????????? //??????????? ???????????! ??? ???????????????? ??????!



	iter_betta = get_betta(curGr, curS);



	vcopy(prev_q, cur_q, K + 1);		// ???????? ??????????

	get_q(curS, prev_q, iter_betta, cur_q);

	//????? ??? k = 0: curS_0, cur_q_1, prev_q_0, curGr_0



	double eps;

	double prev_eps = 1.0;

	eps = get_discrepancy(cur_q, analytic);



	while ( eps > accuracy )

	{

		iter_num++;

	

		vcopy(prevS, curS, K + 1);		//prevS = curS;			//???????? ?????????? S		// ??????

		vcopy(prevGr, curGr, K + 1);	//prevGr = curGr;		//???????? ?????????? ????????



		// ?????? ??? ?????

		//????????? cur_q 



		// ?????? ??? ???



		//?????????? ?????????



		//curGr = ... ???????? ????????

		directVector_solver->setGrid(N, K, xb, xe, maxT);

		//????????? ????????? (????? ?????...)

//		directVector_solver->setTask(a, b, c, f, task_lambda, 0, 1, 0, cur_q, phi2, psi); // border - ??????, ??????? ???? ???????? ? ?????? ???????? ?????????

		directVector_solver->setTask(a, b, c, f, task_lambda, 0, 0, 1, cur_q, phi2, psi); // border - ??????, ??????? ???? ???????? ? ?????? ???????? ?????????

		directVector_solver->solve_implicit();

		double* right_border = directVector_solver->getRightBorder();



		//????????? ????????? (????? ?????...) ???? ??????????? ??? ??????? ? ?????? ???????.

		backTime_solver->setTask(a, b, c, task_lambda, 0, task_lambda, 0, xb, xe, maxT, N, K, right_border, analytic);

		backTime_solver->getSolve();

		backTime_solver->getGradient(curGr);

		delete [] right_border;



		iter_gamma = get_gamma(prevGr, curGr);

		get_S(curGr, prevS, iter_gamma, curS);	//curS - output

		iter_betta = get_betta(curGr, curS);

		vcopy(prev_q, cur_q, K + 1);

		//swap (prev_q, cur_q);		// ???????? ??????????

		get_q(curS, prev_q, iter_betta, cur_q);		//cur_q - output



		prev_eps = eps;

		eps = get_discrepancy(cur_q, analytic);



		if (prev_eps < eps)

		{

			printf("Bad luck, your process is foolish\n");

			break;

		}

	}



	printf("eps = %lf\n", eps);



	printf("q:\n");

	for (int k = 0; k < K + 1; ++k)

		printf("q[%d] = %lf\n", k, cur_q[k]);



	delete [] prev_q;

	delete [] curGr;

	delete [] prevGr;

	delete [] curS;

	delete [] prevS;

}



double CReverseTask_solver::get_discrepancy(double* _cur_q, double* _analytic)

{

	directVAq_solver->setGrid(N, K, xb, xe, maxT);

	directVAq_solver->setTask(a, b, c, f, alpha1, betta1, alpha2, betta2, _cur_q, phi2, psi);

	directVAq_solver->solve_implicit();

	//print_table(directVAq_solver->getTable(), "q_k+1.txt");



	double* Aq;

	Aq = directVAq_solver->getLeftBorder();

	double sum = 0.0;



	for (int i = 0; i < K + 1; ++i)

	{

		sum += pow(Aq[i] - _analytic[i], 2.0);

	}



	sum *= stepT;



	return sum;

}



double CReverseTask_solver::get_betta(double* _curGr, double* _curS)

{

	double answer = 0.0;

	double* VectorAs;



	directVAs_solver->setGrid(N, K, xb, xe, maxT);

	directVAs_solver->setTask(a, b, c, f, alpha1, betta1, alpha2, betta2, _curS, phi2, psi);

	directVAs_solver->solve_implicit();

	VectorAs = directVAs_solver->getLeftBorder();



	double sum_up = 0.0;

	double sum_down = 0.0;



	for(int i = 0; i < K + 1; i++)

	{

		sum_up += _curGr[i] * _curS[i];

		sum_down += VectorAs[i] * VectorAs[i];

	}



	answer = sum_up / (2.0 * sum_down);



	return answer;

}



double CReverseTask_solver::get_gamma(double* _prevGr, double* _curGr)

{

	double ans = 0.0;

	double sum_cur = 0.0;

	double sum_prev = 0.0;



	for(int i = 0; i < K + 1; i++)

	{

		sum_cur += _curGr[i] * _curGr[i];

		sum_prev += _prevGr[i] * _prevGr[i];

	}



	ans = sum_cur / sum_prev;

	return ans;

}



void CReverseTask_solver::get_S(double* _curGr, double* _prevS, double _iter_gamma, double* _curS)	// ???????? ?????? ? curS

{

	for (int i = 0; i < K + 1; ++i)

	{

		_curS[i] = _curGr[i] + _iter_gamma * _prevS[i];

	}

}



void CReverseTask_solver::get_q(double* _curS, double* _prev_q, double _iter_betta, double* _cur_q)	// ???????? ?????? ? cur_q

{

	for (int i = 0; i < K + 1; ++i)

	{

		_cur_q[i] = _prev_q[i] - _iter_betta * _curS[i];

	}

}



void CReverseTask_solver::print_table(double** a, string fname)

{

	ofstream out(fname.c_str());

	for (int i = 0; i < K+1; ++i)

	{

		for (int j = 0; j < N + 1; ++j)

		{

			out<< fixed;

		out <<  setprecision (5)  <<"t = " << i * stepT <<"   x = " << j * (xe - xb) / N << "   u(x,t) = " << a[i][j] << endl;

		}

			out << endl;

	}

}



//////////////////////////////////////////////////////////////////////////