/*!\file: CreateParametersAutodiff.cpp
 * \brief driver for creating parameters dataset, for autodiff analysis.
 */ 

#include "../../../toolkits/toolkits.h"
#include "../../../classes/classes.h"
#include "../../../shared/shared.h"
#include "../ModelProcessorx.h"

void CreateParametersAutodiff(Parameters* parameters,IoModel* iomodel){

	int         i;
	bool        autodiff_analysis;
	int         num_dependent_objects;
	int         num_dep=0;
	char**      names=NULL;
	int*        types=NULL;
	int         dummy;
	char*       autodiff_driver=NULL;
	int*        indices=NULL;
	int         num_indices;
	char* options=NULL;

	IssmDouble* xp=NULL;
	IssmDouble* xp_backup=NULL;
	int         num_ind,local_num_ind;
	DataSet*    dependent_objects=NULL;

	/*retrieve some parameters: */
	iomodel->Constant(&autodiff_analysis,AutodiffIsautodiffEnum);

	if(autodiff_analysis){

		#ifdef _HAVE_ADOLC_

		/*Copy some parameters from IoModel to parameters dataset: {{{*/
		parameters->AddObject(iomodel->CopyConstantObject(AutodiffKeepEnum));
		parameters->AddObject(iomodel->CopyConstantObject(AutodiffObufsizeEnum));
		parameters->AddObject(iomodel->CopyConstantObject(AutodiffCbufsizeEnum));
		parameters->AddObject(iomodel->CopyConstantObject(AutodiffLbufsizeEnum));
		parameters->AddObject(iomodel->CopyConstantObject(AutodiffTbufsizeEnum));
		parameters->AddObject(iomodel->CopyConstantObject(AutodiffGcTriggerRatioEnum));
		parameters->AddObject(iomodel->CopyConstantObject(AutodiffGcTriggerMaxSizeEnum));
		/*}}}*/
		/*retrieve driver: {{{*/
		iomodel->Constant(&autodiff_driver,AutodiffDriverEnum);
		parameters->AddObject(iomodel->CopyConstantObject(AutodiffDriverEnum));

		if(strcmp(autodiff_driver,"fos_forward")==0){
			parameters->AddObject(iomodel->CopyConstantObject(AutodiffFosForwardIndexEnum));
		}
		else if(strcmp(autodiff_driver,"fos_reverse")==0){
			parameters->AddObject(iomodel->CopyConstantObject(AutodiffFosReverseIndexEnum));
		}
		else if(strcmp(autodiff_driver,"fov_forward")==0){
			/*Retrieve list of indices: */
			iomodel->FetchData(&indices,&num_indices,&dummy,AutodiffFovForwardIndicesEnum);
			parameters->AddObject(new IntMatParam(AutodiffFovForwardIndicesEnum,indices,num_indices,1));
			xDelete<int>(indices);
		}
		xDelete<char>(autodiff_driver);
		/*}}}*/
		/*Deal with dependents first: {{{*/
		iomodel->Constant(&num_dependent_objects,AutodiffNumDependentObjectsEnum);
		dependent_objects=new DataSet();
		num_dep=0;

		if(num_dependent_objects){
			iomodel->FetchData(&names,&dummy,AutodiffDependentObjectNamesEnum);
			iomodel->FetchData(&types,&dummy,&dummy,AutodiffDependentObjectTypesEnum);
			iomodel->FetchData(&indices,&dummy,&dummy,AutodiffDependentObjectIndicesEnum);

			for(i=0;i<num_dependent_objects;i++){
				DependentObject* dep=new DependentObject(names[i],types[i],indices[i]);
				dependent_objects->AddObject(dep);
				num_dep+=dep->NumDependents();
			}

			/*Free ressources:*/
			for(i=0;i<num_dependent_objects;i++){
				char* string=names[i]; xDelete<char>(string);
			}
			xDelete<char*>(names);
			xDelete<int>(types);
			xDelete<int>(indices);
		}
		parameters->AddObject(new DataSetParam(AutodiffDependentObjectsEnum,dependent_objects));
		parameters->AddObject(new IntParam(AutodiffNumDependentsEnum,num_dep));

		delete dependent_objects;
		/*}}}*/
		/*Deal with independents: {{{*/

		/*Independents have already been recovered in iomodel->DeclareIndependents. Just do some more processing. 
		 *In particular, figure out num_independents, and create the state vector xp, or size num_independents x 1 :*/
		num_ind=0;
		for(i=0;i<iomodel->independent_objects->Size();i++){
			IndependentObject* ind=(IndependentObject*)iomodel->independent_objects->GetObjectByOffset(i);
			num_ind+=ind->NumIndependents();
		}
		if(num_ind){
			xp=xNew<IssmDouble>(num_ind);
			xp_backup=xp;
			for(i=0;i<iomodel->independent_objects->Size();i++){
				IndependentObject* ind=(IndependentObject*)iomodel->independent_objects->GetObjectByOffset(i);
				ind->FillIndependents(iomodel->data,xp);
				local_num_ind=ind->NumIndependents(); xp=xp+local_num_ind;
			}
			xp=xp_backup; parameters->AddObject(new DoubleVecParam(AutodiffXpEnum,xp,num_ind));
			xDelete<IssmDouble>(xp);
		}
		parameters->AddObject(new IntParam(AutodiffNumIndependentsEnum,num_ind));

		/*Don't forget to copy  iomodel->independent_objects to parameters: */
		parameters->AddObject(new DataSetParam(AutodiffIndependentObjectsEnum,iomodel->independent_objects));
		/*}}}*/

		#endif
	}

	#ifdef _HAVE_ADOLC_
		/*initialize a placeholder to store solver pointers: {{{*/
		GenericParam<Adolc_edf> *theAdolcEDF_p=new GenericParam<Adolc_edf>(AdolcParamEnum);

		/*Solver pointers depend on what type of solver we are implementing: */
		options=OptionsFromAnalysis(parameters,DefaultAnalysisEnum); //options database is not filled in yet, use default.
		ToolkitOptions::Init(options);

		switch(IssmSolverTypeFromToolkitOptions()){
			case MumpsEnum:{
				#ifdef _HAVE_MUMPS_
				theAdolcEDF_p->GetParameterValue().myEDF_for_solverx_p=reg_ext_fct(mumpsSolveEDF);
				#else
				_error_("requesting mumps solver without MUMPS being compiled in!");
				#endif
				break;
							}
			case GslEnum: {
				#ifdef _HAVE_GSL_
				theAdolcEDF_p->GetParameterValue().myEDF_for_solverx_p=reg_ext_fct(EDF_for_solverx);
				#endif
			    break;
						  }
			default:
				_error_("solver type not supported yet!");
		}

		// to save some space:
		// we know we won't use adolc inside of  the solver:
		theAdolcEDF_p->GetParameterValue().myEDF_for_solverx_p->nestedAdolc=false;
		// the solution vector is just allocated and doesn't have a meaningful prior value
		theAdolcEDF_p->GetParameterValue().myEDF_for_solverx_p->dp_y_priorRequired=false;
		// the solver wrapper makes sure the matrix and the right hand side don't change
		theAdolcEDF_p->GetParameterValue().myEDF_for_solverx_p->dp_x_changes=false;
		parameters->AddObject(theAdolcEDF_p);

		/*Free ressources: */
		xDelete<char>(options);
		/*}}}*/
	#endif

}
