/*!\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** pparameters,IoModel* iomodel,int solution_type,int analysis_type){

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

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

	/*Get parameters: */
	parameters=*pparameters;

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

		/*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);
		}
		/*}}}*/
		/*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,&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:*/
			xDelete<int>(names);
			xDelete<int>(types);
		}
		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));
		}
		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);
		#ifdef _HAVE_GSL_
		theAdolcEDF_p->GetParameterValue().myEDF_for_solverx_p=reg_ext_fct(EDF_for_solverx);
		#endif
		#ifdef _HAVE_MUMPS_
		theAdolcEDF_p->GetParameterValue().myEDF_for_solverx_p=reg_ext_fct(mumpsSolveEDF);
		#endif
		// 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);
		/*}}}*/
	#endif

	/*Assign output pointer: */
	*pparameters=parameters;

}
