/*!\file: control_core.cpp
 * \brief: core of the control solution 
 */ 

#include "../toolkits/toolkits.h"
#include "../classes/objects/objects.h"
#include "../shared/shared.h"
#include "../io/io.h"
#include "../EnumDefinitions/EnumDefinitions.h"
#include "./solutions.h"
#include "../modules/modules.h"
#include "../include/include.h"
#include "../solvers/solvers.h"

void control_core(FemModel* femmodel){

	int     i,n;

	/*parameters: */
	int     num_controls,num_responses;
	int     nsteps;
	IssmDouble  tol_cm;
	bool    cm_gradient;
	int     dim;
	int     solution_type;
	bool    isstokes;
	bool    dakota_analysis=false;

	int*    control_type = NULL;
	IssmDouble* responses=NULL;
	int*    step_responses=NULL;
	IssmDouble* maxiter=NULL;
	IssmDouble* cm_jump=NULL;

	/*intermediary: */
	IssmDouble  search_scalar=1;
	OptArgs optargs;
	OptPars optpars;

	/*Solution and Adjoint core pointer*/
	void (*solutioncore)(FemModel*)=NULL;
	void (*adjointcore)(FemModel*)=NULL;

	/*output: */
	IssmDouble* J=NULL;

	/*Recover parameters used throughout the solution*/
	femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
	femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
	femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
	femmodel->parameters->FindParam(&responses,NULL,NULL,InversionCostFunctionsEnum);
	femmodel->parameters->FindParam(&nsteps,InversionNstepsEnum);
	femmodel->parameters->FindParam(&maxiter,NULL,InversionMaxiterPerStepEnum);
	femmodel->parameters->FindParam(&cm_jump,NULL,InversionStepThresholdEnum);
	femmodel->parameters->FindParam(&tol_cm,InversionCostFunctionThresholdEnum);
	femmodel->parameters->FindParam(&cm_gradient,InversionGradientOnlyEnum);
	femmodel->parameters->FindParam(&dim,MeshDimensionEnum);
	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
	femmodel->parameters->FindParam(&isstokes,FlowequationIsstokesEnum);
	femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
	femmodel->parameters->SetParam(false,SaveResultsEnum);

	/*out of solution_type, figure out solution core and adjoint function pointer*/
	CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
	AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);

	/*Launch once a complete solution to set up all inputs*/
	if(VerboseControl()) _pprintLine_("   preparing initial solution");
	if(isstokes) solutioncore(femmodel);

	/*Initialize responses: */
	J=xNew<IssmDouble>(nsteps);
	step_responses=xNew<int>(num_responses);

	/*Initialize some of the BrentSearch arguments: */
	optargs.femmodel=femmodel;
	optpars.xmin=0; optpars.xmax=1;

	/*Start looping: */
	for(n=0;n<nsteps;n++){

		/*Display info*/
		if(VerboseControl()) _pprintLine_("\n" << "   control method step " << n+1 << "/" << nsteps);
		for(i=0;i<num_responses;i++) step_responses[i]=reCast<int,IssmDouble>(responses[n*num_responses+i]);
		femmodel->parameters->SetParam(step_responses,1,num_responses,StepResponsesEnum);

		/*In steady state inversion, compute new temperature field now*/
		if(solution_type==SteadystateSolutionEnum) solutioncore(femmodel);

		if(VerboseControl()) _pprintLine_("   compute adjoint state:");
		adjointcore(femmodel);
		gradient_core(femmodel,n,search_scalar==0);

		/*Return gradient if asked: */
		if (cm_gradient){
			InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,GradientEnum);
			goto cleanup_and_return;
		}

		if(VerboseControl()) _pprintLine_("   optimizing along gradient direction");
		optpars.maxiter=reCast<int,IssmDouble>(maxiter[n]); optpars.cm_jump=cm_jump[n];
		BrentSearch(&search_scalar,J+n,&optpars,&objectivefunction,&optargs);

		if(VerboseControl()) _pprintLine_("   updating parameter using optimized search scalar"); //true means update save controls
		InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,true);

		if(controlconvergence(J[n],tol_cm)) break;
	}

	if(VerboseControl()) _pprintLine_("   preparing final solution");
	femmodel->parameters->SetParam(true,SaveResultsEnum);
	solutioncore(femmodel);

	/*some results not computed by steadystate_core or diagnostic_core: */
	if(!dakota_analysis){ //do not save this if we are running the control core from a qmu run!
		for(i=0;i<num_controls;i++) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i]);

		#ifdef _HAVE_ADOLC_
		IssmPDouble* J_passive=xNew<IssmPDouble>(nsteps);
		for(int i=0;i<nsteps;i++)J_passive[i]=reCast<IssmPDouble>(J[i]);
		femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,J_passive,nsteps,1,1,0));
		xDelete<IssmPDouble>(J_passive);
		#else
		femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,J,nsteps,1,1,0));
		#endif
	}

	cleanup_and_return:
	/*Free ressources: */
	xDelete<int>(control_type);
	xDelete<int>(step_responses);
	xDelete<IssmDouble>(maxiter);
	xDelete<IssmDouble>(responses);
	xDelete<IssmDouble>(cm_jump);
	xDelete<IssmDouble>(J);
}
