/*!\file: CreateParameters.cpp
 * \brief general driver for creating parameters dataset
 */ 

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

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

void CreateParameters(Parameters** pparameters,IoModel* iomodel,char* rootpath,const int solution_type,int analysis_type,int analysis_counter){

	int         i,j,m,k;
	int         numoutputs;
	Parameters *parameters       = NULL;
	IssmDouble *requestedoutputs = NULL;
	IssmDouble  time;
	bool        ispdd,isdelta18o;

	/*parameters for mass flux: {{{*/
	int          mass_flux_num_profiles     = 0;
	bool         qmu_mass_flux_present      = false;
	bool         autodiff_mass_flux_present = false;
	bool         mass_flux_present          = false;
	IssmDouble **array                      = NULL;
	int         *mdims_array                = NULL;
	int         *ndims_array                = NULL;
	IssmDouble  *temp_matrix                = NULL;
	int          temp_m,temp_n;
	IssmDouble  *matrix                     = NULL;
	int          count;
	/*}}}*/

	if(*pparameters)return; //do not create parameters twice!

	/*Initialize dataset: */
	parameters = new Parameters();

	/*Copy some constants from iomodel */
	parameters->AddObject(iomodel->CopyConstantObject(MeshDimensionEnum));
	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSIAEnum));
	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSSAEnum));
	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsL1L2Enum));
	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsHOEnum));
	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsFSEnum));
	parameters->AddObject(iomodel->CopyConstantObject(FlowequationFeFSEnum));
	parameters->AddObject(iomodel->CopyConstantObject(SettingsOutputFrequencyEnum));
	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRestolEnum));
	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceReltolEnum));
	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceAbstolEnum));
	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceIsnewtonEnum));
	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceMaxiterEnum));
	parameters->AddObject(iomodel->CopyConstantObject(StressbalancePenaltyFactorEnum));
	parameters->AddObject(iomodel->CopyConstantObject(SteadystateReltolEnum));
	parameters->AddObject(iomodel->CopyConstantObject(SteadystateMaxiterEnum));
	parameters->AddObject(iomodel->CopyConstantObject(ConstantsYtsEnum));
	parameters->AddObject(iomodel->CopyConstantObject(TimesteppingStartTimeEnum));
	parameters->AddObject(iomodel->CopyConstantObject(TimesteppingFinalTimeEnum));
	parameters->AddObject(iomodel->CopyConstantObject(TimesteppingTimeAdaptEnum));
	parameters->AddObject(iomodel->CopyConstantObject(TimesteppingTimeStepEnum));
	parameters->AddObject(iomodel->CopyConstantObject(TimesteppingCflCoefficientEnum));
	parameters->AddObject(iomodel->CopyConstantObject(MasstransportIsfreesurfaceEnum));
	parameters->AddObject(iomodel->CopyConstantObject(MasstransportHydrostaticAdjustmentEnum));
	parameters->AddObject(iomodel->CopyConstantObject(MasstransportStabilizationEnum));
	parameters->AddObject(iomodel->CopyConstantObject(MasstransportMinThicknessEnum));
	parameters->AddObject(iomodel->CopyConstantObject(MasstransportPenaltyFactorEnum));
	parameters->AddObject(iomodel->CopyConstantObject(ThermalPenaltyFactorEnum));
	parameters->AddObject(iomodel->CopyConstantObject(SettingsLowmemEnum));
	parameters->AddObject(iomodel->CopyConstantObject(DebugProfilingEnum));
	parameters->AddObject(iomodel->CopyConstantObject(MeshAverageVertexConnectivityEnum));
	parameters->AddObject(iomodel->CopyConstantObject(ConstantsReferencetemperatureEnum));
	parameters->AddObject(iomodel->CopyConstantObject(BalancethicknessStabilizationEnum));
	parameters->AddObject(iomodel->CopyConstantObject(GroundinglineMeltingRateEnum));
	parameters->AddObject(iomodel->CopyConstantObject(ThermalMaxiterEnum));
	parameters->AddObject(iomodel->CopyConstantObject(ThermalStabilizationEnum));
	parameters->AddObject(iomodel->CopyConstantObject(ThermalPenaltyThresholdEnum));
	parameters->AddObject(iomodel->CopyConstantObject(ThermalPenaltyLockEnum));
	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRiftPenaltyThresholdEnum));
	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceFSreconditioningEnum));
	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceShelfDampeningEnum));
	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceViscosityOvershootEnum));
	parameters->AddObject(iomodel->CopyConstantObject(SettingsWaitonlockEnum));
	parameters->AddObject(iomodel->CopyConstantObject(MeshNumberofelementsEnum));
	parameters->AddObject(iomodel->CopyConstantObject(MeshNumberofverticesEnum));
	parameters->AddObject(iomodel->CopyConstantObject(MeshNumberoflayersEnum));
	parameters->AddObject(iomodel->CopyConstantObject(SettingsIoGatherEnum));
	parameters->AddObject(iomodel->CopyConstantObject(SettingsResultsAsPatchesEnum));
	parameters->AddObject(iomodel->CopyConstantObject(GroundinglineMigrationEnum));
	parameters->AddObject(iomodel->CopyConstantObject(TransientIsstressbalanceEnum));
	parameters->AddObject(iomodel->CopyConstantObject(TransientIsmasstransportEnum));
	parameters->AddObject(iomodel->CopyConstantObject(TransientIsthermalEnum));
	parameters->AddObject(iomodel->CopyConstantObject(TransientIsgroundinglineEnum));
	parameters->AddObject(iomodel->CopyConstantObject(TransientIsgiaEnum));
	parameters->AddObject(iomodel->CopyConstantObject(ThermalIsenthalpyEnum));
	parameters->AddObject(iomodel->CopyConstantObject(MaterialsRheologyLawEnum));
	parameters->AddObject(iomodel->CopyConstantObject(AutodiffIsautodiffEnum));
	parameters->AddObject(iomodel->CopyConstantObject(QmuIsdakotaEnum));
	parameters->AddObject(iomodel->CopyConstantObject(InversionIscontrolEnum));
	parameters->AddObject(iomodel->CopyConstantObject(InversionTaoEnum));
	parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIspddEnum));
	parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsdelta18oEnum));
	parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIssmbgradientsEnum));
	parameters->AddObject(iomodel->CopyConstantObject(GiaCrossSectionShapeEnum));

	iomodel->Constant(&ispdd,SurfaceforcingsIspddEnum);
	if(ispdd){

		parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsDesfacEnum));
		parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsS0pEnum));

		iomodel->Constant(&isdelta18o,SurfaceforcingsIsdelta18oEnum);
		if(isdelta18o){
			IssmDouble *temp = NULL;
			IssmDouble  yts;
			int         N,M;

			iomodel->Constant(&yts,ConstantsYtsEnum);

			iomodel->FetchData(&temp,&N,&M,SurfaceforcingsDelta18oEnum); _assert_(N==2);
			for(i=0;i<M;i++) temp[M+i]=yts*temp[M+i];
			parameters->AddObject(new TransientParam(SurfaceforcingsDelta18oEnum,&temp[0],&temp[M],M));
			iomodel->DeleteData(temp,SurfaceforcingsDelta18oEnum);

			iomodel->FetchData(&temp,&N,&M,SurfaceforcingsDelta18oSurfaceEnum); _assert_(N==2);
			for(i=0;i<M;i++) temp[M+i]=yts*temp[M+i];
			parameters->AddObject(new TransientParam(SurfaceforcingsDelta18oSurfaceEnum,&temp[0],&temp[M],M));
			iomodel->DeleteData(temp,SurfaceforcingsDelta18oSurfaceEnum);
		}
	}

	/*some parameters that did not come with the iomodel: */
	parameters->AddObject(new IntParam(SolutionTypeEnum,solution_type));
	parameters->AddObject(new IntParam(AnalysisTypeEnum,analysis_type));
	parameters->AddObject(new IntParam(AnalysisCounterEnum,analysis_counter));

	iomodel->Constant(&time,TimesteppingStartTimeEnum);
	parameters->AddObject(new DoubleParam(TimeEnum,time));  //start at time 0 by default for all solutions FIXME: to be deleted
	parameters->AddObject(new IntParam(StepEnum,1));  //start at time 0 by default for all solutions FIXME: to be deleted

	/*By default, save all results*/
	parameters->AddObject(new BoolParam(SaveResultsEnum,true));

	/*Requested outputs*/
	iomodel->FetchData(&requestedoutputs,&numoutputs,NULL,StressbalanceRequestedOutputsEnum);
	parameters->AddObject(new IntParam(StressbalanceNumRequestedOutputsEnum,numoutputs));
	if(numoutputs)parameters->AddObject(new IntVecParam(StressbalanceRequestedOutputsEnum,requestedoutputs,numoutputs));
	iomodel->DeleteData(requestedoutputs,StressbalanceRequestedOutputsEnum);

	iomodel->FetchData(&requestedoutputs,&numoutputs,NULL,TransientRequestedOutputsEnum);
	parameters->AddObject(new IntParam(TransientNumRequestedOutputsEnum,numoutputs));
	if(numoutputs)parameters->AddObject(new IntVecParam(TransientRequestedOutputsEnum,requestedoutputs,numoutputs));
	iomodel->DeleteData(requestedoutputs,TransientRequestedOutputsEnum);

	iomodel->FetchData(&requestedoutputs,&numoutputs,NULL,SteadystateRequestedOutputsEnum);
	parameters->AddObject(new IntParam(SteadystateNumRequestedOutputsEnum,numoutputs));
	if(numoutputs)parameters->AddObject(new IntVecParam(SteadystateRequestedOutputsEnum,requestedoutputs,numoutputs));
	iomodel->DeleteData(requestedoutputs,SteadystateRequestedOutputsEnum);

	iomodel->FetchData(&requestedoutputs,&numoutputs,NULL,MasstransportRequestedOutputsEnum);
	parameters->AddObject(new IntParam(MasstransportNumRequestedOutputsEnum,numoutputs));
	if(numoutputs)parameters->AddObject(new IntVecParam(MasstransportRequestedOutputsEnum,requestedoutputs,numoutputs));
	iomodel->DeleteData(requestedoutputs,MasstransportRequestedOutputsEnum);

	/*Deal with mass flux segments: {{{*/
	iomodel->FetchData(&qmu_mass_flux_present,QmuMassFluxSegmentsPresentEnum);
	iomodel->FetchData(&autodiff_mass_flux_present,AutodiffMassFluxSegmentsPresentEnum);

	if(qmu_mass_flux_present || autodiff_mass_flux_present)mass_flux_present=true;
	else mass_flux_present=false;
	parameters->AddObject(new BoolParam(MassFluxSegmentsPresentEnum,mass_flux_present));

	if(mass_flux_present){

		/*Fetch the mass flux segments necessary to compute the mass fluxes.  Build a DoubleMatArrayParam object out of them: */ 
		iomodel->FetchData(&array,&mdims_array,&ndims_array,&mass_flux_num_profiles,MassFluxSegmentsEnum);
		if(mass_flux_num_profiles==0)_error_("mass_flux_num_profiles is 0, when MassFlux computations were requested!");

		/*Go through segments, and extract those that belong to this cpu: */
		for(i=0;i<mass_flux_num_profiles;i++){
			temp_matrix=array[i];
			temp_m=mdims_array[i];
			temp_n=ndims_array[i];

			m=0;
			for(j=0;j<temp_m;j++){
				if (  iomodel->my_elements[reCast<int>(*(temp_matrix+5*j+4))-1] )m++;
			}
			if(m){
				matrix=xNewZeroInit<IssmDouble>(5*m);
				count=0;
				for(j=0;j<temp_m;j++){
					if (iomodel->my_elements[reCast<int>(*(temp_matrix+5*j+4))-1]){
						for(k=0;k<5;k++)*(matrix+5*count+k)=*(temp_matrix+5*j+k);
						count++;
					}
				}
			}
			else{
				matrix=NULL;
			}

			/*Assign: */
			array[i]=matrix;
			mdims_array[i]=m;
			ndims_array[i]=5;

			/*Free temporary matrix: */
			xDelete<IssmDouble>(temp_matrix);
		}

		/*Ok, we have an array of segments, different on every cpu. Create a DoubleMatArrayParam object with it: */
		parameters->AddObject(new DoubleMatArrayParam(MassFluxSegmentsEnum,array,mass_flux_num_profiles,mdims_array,ndims_array));

		/*Free data: */
		for(i=0;i<mass_flux_num_profiles;i++){
			IssmDouble* matrix=array[i];
			xDelete<IssmDouble>(matrix);
		}
		xDelete<int>(mdims_array); 
		xDelete<int>(ndims_array);
		xDelete<IssmDouble*>(array);
	}
	/*}}}*/

	/*Solution specific parameters (FIXME: extend to other params)*/
	#ifdef _HAVE_HYDROLOGY_
	CreateParametersHydrologyShreve(&parameters,iomodel,solution_type,analysis_type);
	CreateParametersHydrologyDCInefficient(&parameters,iomodel,solution_type,analysis_type);
	CreateParametersHydrologyDCEfficient(&parameters,iomodel,solution_type,analysis_type);
	#endif

	/*Before returning, create parameters in case we are running Qmu or control types runs: */
	#ifdef _HAVE_CONTROL_
	CreateParametersControl(&parameters,iomodel,solution_type,analysis_type);
	#endif

	#ifdef _HAVE_DAKOTA_
	CreateParametersDakota(&parameters,iomodel,rootpath,solution_type,analysis_type);
	#endif

 	#ifdef _HAVE_ADOLC_
	CreateParametersAutodiff(&parameters,iomodel,solution_type,analysis_type);
	#endif

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