/*!\file ad_core
 * \brief: compute outputs from the AD mode,  using our dependents and independents, and drivers available in Adolc.
 */

/*Includes: {{{*/
#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include <set>
#include "./analyses.h"
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../shared/shared.h"
#include "../modules/modules.h"
#include "../solutionsequences/solutionsequences.h"
/*}}}*/

void ad_core(FemModel* femmodel){

	/*diverse: */
	int     i;
	int     dummy;
	int     num_dependents=0;
	int     num_independents=0;
	bool    isautodiff       = false;
	char   *driver           = NULL;
	size_t  tape_stats[11];

	/*state variables: */
	IssmDouble *axp = NULL;
	double     *xp  = NULL;
	int my_rank=IssmComm::GetRank();


	/*AD mode on?: */
	femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);

	if(isautodiff){

		#ifdef _HAVE_ADOLC_

			/*First, stop tracing: */
			trace_off();

			/*preliminary checks: */
			if (my_rank==0) {
				femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
				femmodel->parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum);
				if(!(num_dependents*num_independents)) return;
			}

			if(VerboseAutodiff())_printf0_("   start ad core\n");

			/*retrieve state variable: */
			femmodel->parameters->FindParam(&axp,&dummy,AutodiffXpEnum);

			/* driver argument */
			xp=xNew<double>(num_independents);
			for(i=0;i<num_independents;i++){
				xp[i]=reCast<double,IssmDouble>(axp[i]);
			}

			/*get the EDF pointer:*/
			ext_diff_fct *anEDF_for_solverx_p=dynamic_cast<GenericParam<Adolc_edf> * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p;

			/*Branch according to AD driver: */
			femmodel->parameters->FindParam(&driver,AutodiffDriverEnum);

			/* these are always needed regardless of the interpreter */
			anEDF_for_solverx_p->dp_x=xNew<double>(anEDF_for_solverx_p->max_n);
			anEDF_for_solverx_p->dp_y=xNew<double>(anEDF_for_solverx_p->max_m);

			if (strcmp(driver,"fos_forward")==0){

				int     anIndepIndex;
				double *tangentDir         = NULL;
				double *jacTimesTangentDir = NULL;
				double *theOutput          = NULL;

				/*retrieve direction index: */
				femmodel->parameters->FindParam(&anIndepIndex,AutodiffFosForwardIndexEnum);

				if (anIndepIndex<0 || anIndepIndex>=num_independents) _error_("index value for AutodiffFosForwardIndexEnum should be in [0,num_independents-1]");

				tangentDir=xNewZeroInit<double>(num_independents);
				tangentDir[anIndepIndex]=1.0;

				jacTimesTangentDir=xNew<double>(num_dependents);
				theOutput=xNew<double>(num_dependents);

				/*set the forward method function pointer: */
#ifdef _HAVE_GSL_
				anEDF_for_solverx_p->fos_forward=EDF_fos_forward_for_solverx;
#endif

				/*allocate the space for the parameters to invoke the EDF fos_forward:*/
				anEDF_for_solverx_p->dp_X=xNew<double>(anEDF_for_solverx_p->max_n);
				anEDF_for_solverx_p->dp_Y=xNew<double>(anEDF_for_solverx_p->max_m);

				/*call driver: */
				fos_forward(my_rank,num_dependents,num_independents, 0, xp, tangentDir, theOutput, jacTimesTangentDir );

				/*add to results*/
				femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,jacTimesTangentDir,num_dependents,1,1,0.0));

				/*free resources :*/
				xDelete(theOutput);
				xDelete(jacTimesTangentDir);
				xDelete(tangentDir);
			}
			else if ((strcmp(driver,"fov_forward")==0) || (strcmp(driver,"fov_forward_all")==0)){

				int      tangentDirNum;
				int      dummy;
				int     *indepIndices  = NULL;
				double **jacTimesSeed  = NULL;
				double **seed          = NULL;
				double  *theOutput     = NULL;
				std::set<unsigned int> anIndexSet;

				/*retrieve directions:*/
				if (strcmp(driver,"fov_forward_all")==0){
					tangentDirNum=num_independents;
					indepIndices=xNewZeroInit<int>(tangentDirNum);
					for(i=0;i<num_independents;i++)indepIndices[i]=1;
				}
				else{
					femmodel->parameters->FindParam(&indepIndices,&tangentDirNum,&dummy,AutodiffFovForwardIndicesEnum);
				}

				/*Some checks: */
				if (tangentDirNum<1 || tangentDirNum>num_independents) _error_("tangentDirNum should be in [1,num_independents]");

				/* full Jacobian or Jacobian projection:*/
				jacTimesSeed=xNew<double>(num_dependents,tangentDirNum);

				/*set the forward method function pointers: */
#ifdef _HAVE_GSL_
				anEDF_for_solverx_p->fov_forward=EDF_fov_forward_for_solverx;
#endif
				// anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx;

				/*allocate the space for the parameters to invoke EDF fov_forward:*/
				anEDF_for_solverx_p->dpp_X=xNew<double>(anEDF_for_solverx_p->max_n, tangentDirNum);
				anEDF_for_solverx_p->dpp_Y=xNew<double>(anEDF_for_solverx_p->max_m, tangentDirNum);

				/*seed matrix: */
				seed=xNewZeroInit<double>(num_independents,tangentDirNum);

				/*collect indices in a set to prevent accidental duplicates as long as we don't do compression:*/
				for (int i=0; i<tangentDirNum; ++i) {
					/* make sure the index is in range*/
					if (indepIndices[i]>num_independents) {
						_error_("indepIndices values must be in [0,num_independents-1]");
					}
					if (anIndexSet.find(indepIndices[i])!=anIndexSet.end()) {
						_error_("duplicate indepIndices values are not allowed until we implement Jacobian decompression");
					}
					anIndexSet.insert(indepIndices[i]);
					/* now populate the seed matrix from the set of independent indices;
					 * simple setup with a single 1.0 per column and at most a single 1.0 per row*/
					seed[indepIndices[i]][i]=1.0;
				}

				/*allocate output: */
				theOutput=xNew<double>(num_dependents);

				/*call driver: */
				fov_forward(my_rank,num_dependents,num_independents, tangentDirNum, xp, seed, theOutput, jacTimesSeed );
				/*Free resources: */
				xDelete(theOutput);
				xDelete(indepIndices);
				xDelete(seed);

				/*add to results: */
				femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,*jacTimesSeed,num_dependents*tangentDirNum,1,1,0.0));

				/*Free resources: */
				xDelete(jacTimesSeed);
				xDelete(indepIndices);
			}
			else if (strcmp(driver,"fos_reverse")==0) {

				int     aDepIndex=0;
				double *aWeightVector=NULL;
				double *weightVectorTimesJac=NULL;

				/*retrieve direction index: */
				femmodel->parameters->FindParam(&aDepIndex,AutodiffFosReverseIndexEnum);
				aWeightVector=xNewZeroInit<double>(num_dependents);
				if (my_rank==0) {
					if (aDepIndex<0 || aDepIndex>=num_dependents) _error_("index value for AutodiffFosReverseIndexEnum should be in [0,num_dependents-1]");
					aWeightVector[aDepIndex]=1.0;
				}
				weightVectorTimesJac=xNew<double>(num_independents);

				/*set the forward method function pointer: */
#ifdef _HAVE_GSL_
				anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx;
#endif
#ifdef _HAVE_MUMPS_
				anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF;
#endif

				/*allocate the space for the parameters to invoke the EDF fos_reverse :*/
				anEDF_for_solverx_p->dp_U=xNew<double>(anEDF_for_solverx_p->max_m);
				anEDF_for_solverx_p->dp_Z=xNew<double>(anEDF_for_solverx_p->max_n);

				/*call driver: */
				fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac );
				if(VerboseAutodiff())_printf_("   done with fos_reverse\n");

				/*add to results*/
				femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,weightVectorTimesJac,num_independents,1,1,0.0));

				/*free resources :*/
				xDelete(weightVectorTimesJac);
				xDelete(aWeightVector);
			}
			else if ((strcmp(driver,"fov_reverse")==0) || (strcmp(driver,"fov_reverse_all")==0)){

				int* depIndices=NULL;
				int weightNum;
				int dummy;
				double **weightsTimesJac=NULL;
				double **weights=NULL;
				std::set<unsigned int> anIndexSet;

				/*retrieve directions:*/
				if (strcmp(driver,"fov_reverse_all")==0){
					weightNum=num_dependents;
					depIndices=xNewZeroInit<int>(weightNum);
					for(i=0;i<num_dependents;i++)depIndices[i]=1;
				}
				else{
					femmodel->parameters->FindParam(&depIndices,&weightNum,&dummy,AutodiffFovForwardIndicesEnum);
				}

				/*Some checks: */
				if (weightNum<1 || weightNum>num_dependents) _error_("tangentDirNum should be in [1,num_dependents]");

				/* full Jacobian or Jacobian projection:*/
				weightsTimesJac=xNew<double>(weightNum,num_independents);

				/*set the forward method function pointers: */
#ifdef _HAVE_GSL_
				anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx;
#endif

				/*allocate the space for the parameters to invoke the EDF fos_reverse :*/
				anEDF_for_solverx_p->dpp_U=xNew<double>(weightNum,anEDF_for_solverx_p->max_m);
				anEDF_for_solverx_p->dpp_Z=xNew<double>(weightNum,anEDF_for_solverx_p->max_n);

				/*seed matrix: */
				weights=xNewZeroInit<double>(weightNum,num_dependents);

				/*collect indices in a set to prevent accidental duplicates as long as we don't do compression:*/
				for (int i=0; i<weightNum; ++i) {
					/* make sure the index is in range*/
					if (depIndices[i]>num_dependents) {
						_error_("depIndices values must be in [0,num_dependents-1]");
					}
					if (anIndexSet.find(depIndices[i])!=anIndexSet.end()) {
						_error_("duplicate depIndices values are not allowed until we implement Jacobian decompression");
					}
					anIndexSet.insert(depIndices[i]);
					/* now populate the seed matrix from the set of independent indices;
					 * simple setup with a single 1.0 per column and at most a single 1.0 per row*/
					weights[depIndices[i]][i]=1.0;
				}

				/*call driver: */
				fov_reverse(my_rank,num_dependents,num_independents, weightNum, weights, weightsTimesJac );

				/*add to results: */
				femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,*weightsTimesJac,weightNum*num_independents,1,1,0.0));

				/*Free resources: */
				xDelete(weights);
				xDelete(weightsTimesJac);
				xDelete(depIndices);
			}
			else _error_("driver: " << driver << " not yet supported!");

			/* delete the allocated space for the parameters:*/
			xDelete(anEDF_for_solverx_p->dp_x);
			xDelete(anEDF_for_solverx_p->dp_X);
			xDelete(anEDF_for_solverx_p->dpp_X);
			xDelete(anEDF_for_solverx_p->dp_y);
			xDelete(anEDF_for_solverx_p->dp_Y);
			xDelete(anEDF_for_solverx_p->dpp_Y);
			xDelete(anEDF_for_solverx_p->dp_U);
			xDelete(anEDF_for_solverx_p->dpp_U);
			xDelete(anEDF_for_solverx_p->dp_Z);
			xDelete(anEDF_for_solverx_p->dpp_Z);

			/*Print statistics:*/
			tapestats(my_rank,tape_stats); //reading of tape statistics
			if(VerboseAutodiff()){
				_printf0_("   ADOLC statistics: \n");
				_printf0_("   "<<setw(45)<<left<<"Number of independents: " <<tape_stats[0] << "\n");
				_printf0_("   "<<setw(45)<<left<<"Number of dependents: " <<tape_stats[1] << "\n");
				_printf0_("   "<<setw(45)<<left<<"Maximal number of live active variables: " <<tape_stats[2] << "\n");
				_printf0_("   "<<setw(45)<<left<<"Size of value stack (number of overwrites): " <<tape_stats[3] << "\n");
				_printf0_("   "<<setw(45)<<left<<"Buffer size (a multiple of eight): " <<tape_stats[4] << "\n");
				_printf0_("   "<<setw(45)<<left<<"Total number of operations recorded: " <<tape_stats[5] << "\n");
			}
			if(VerboseAutodiff())_printf0_("   end AD core\n");

			/*Free resources: */
			xDelete(xp);
			xDelete(axp); 
		#else
			_error_("Should not be requesting AD drivers when an AD library is not available!");
		#endif
	}
}
