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

#include "../toolkits/toolkits.h"
#include "../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;
	double  tol_cm;
	bool    cm_gradient;
	int     dim;
	int     solution_type;
	bool    isstokes;
	bool    dakota_analysis=false;

	int*    control_type = NULL;
	double* responses=NULL;
	int*    step_responses=NULL;
	double* maxiter=NULL;
	double* cm_jump=NULL;
		
	/*intermediary: */
	double  search_scalar=1;
	OptArgs optargs;
	OptPars optpars;

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

	/*output: */
	double* 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);

	/*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*/
	_printf_(VerboseControl(),"%s\n","   preparing initial solution");
	if (isstokes) solutioncore(femmodel);

	/*Initialize responses: */
	J=(double*)xmalloc(nsteps*sizeof(double));
	step_responses=(int*)xmalloc(num_responses*sizeof(int));
		
	/*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*/
		_printf_(VerboseControl(),"\n%s%i%s%i\n","   control method step ",n+1,"/",nsteps);
		for(i=0;i<num_responses;i++) step_responses[i]=(int)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);

		_printf_(VerboseControl(),"%s\n","   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;
		}

		_printf_(VerboseControl(),"%s\n","   optimizing along gradient direction");
		optpars.maxiter=(int)maxiter[n]; optpars.cm_jump=cm_jump[n];
		BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs);

		_printf_(VerboseControl(),"%s\n","   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;
	}

	_printf_(VerboseControl(),"%s\n","   preparing final solution");
	femmodel->parameters->SetParam(false,InversionIscontrolEnum); //needed to turn control result output in solutioncore
	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]);
		femmodel->results->AddObject(new DoubleVecExternalResult(femmodel->results->Size()+1,JEnum,J,nsteps,1,0));
	}

	cleanup_and_return:
	/*Free ressources: */
	xfree((void**)&control_type);
	xfree((void**)&responses);
	xfree((void**)&step_responses);
	xfree((void**)&maxiter);
	xfree((void**)&cm_jump);
	xfree((void**)&J);
	
	/*control_core might be used in Qmu, so leave everything similar to where it started: */
	femmodel->parameters->SetParam(true,InversionIscontrolEnum);
}
