/*!\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 "../../Container/Container.h"
#include "../../toolkits/toolkits.h"
#include "../../EnumDefinitions/EnumDefinitions.h"
#include "../../classes/objects/objects.h"
#include "../../shared/shared.h"
#include "../MeshPartitionx/MeshPartitionx.h"
#include "../../io/io.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;
	bool        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(FlowequationIshutterEnum));
	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsmacayealpattynEnum));
	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsl1l2Enum));
	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsstokesEnum));
	parameters->AddObject(iomodel->CopyConstantObject(SettingsOutputFrequencyEnum));
	parameters->AddObject(iomodel->CopyConstantObject(DiagnosticRestolEnum));
	parameters->AddObject(iomodel->CopyConstantObject(DiagnosticReltolEnum));
	parameters->AddObject(iomodel->CopyConstantObject(DiagnosticAbstolEnum));
	parameters->AddObject(iomodel->CopyConstantObject(DiagnosticIsnewtonEnum));
	parameters->AddObject(iomodel->CopyConstantObject(DiagnosticMaxiterEnum));
	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(PrognosticHydrostaticAdjustmentEnum));
	parameters->AddObject(iomodel->CopyConstantObject(PrognosticStabilizationEnum));
	parameters->AddObject(iomodel->CopyConstantObject(DiagnosticPenaltyFactorEnum));
	parameters->AddObject(iomodel->CopyConstantObject(PrognosticMinThicknessEnum));
	parameters->AddObject(iomodel->CopyConstantObject(PrognosticPenaltyFactorEnum));
	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(HydrologyStabilizationEnum));
	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(DiagnosticRiftPenaltyThresholdEnum));
	parameters->AddObject(iomodel->CopyConstantObject(DiagnosticStokesreconditioningEnum));
	parameters->AddObject(iomodel->CopyConstantObject(DiagnosticShelfDampeningEnum));
	parameters->AddObject(iomodel->CopyConstantObject(DiagnosticViscosityOvershootEnum));
	parameters->AddObject(iomodel->CopyConstantObject(SettingsWaitonlockEnum));
	parameters->AddObject(iomodel->CopyConstantObject(MeshNumberofelementsEnum));
	parameters->AddObject(iomodel->CopyConstantObject(MeshNumberofverticesEnum));
	parameters->AddObject(iomodel->CopyConstantObject(SettingsIoGatherEnum));
	parameters->AddObject(iomodel->CopyConstantObject(SettingsResultsAsPatchesEnum));
	parameters->AddObject(iomodel->CopyConstantObject(GroundinglineMigrationEnum));
	parameters->AddObject(iomodel->CopyConstantObject(TransientIsdiagnosticEnum));
	parameters->AddObject(iomodel->CopyConstantObject(TransientIsprognosticEnum));
	parameters->AddObject(iomodel->CopyConstantObject(TransientIsthermalEnum));
	parameters->AddObject(iomodel->CopyConstantObject(TransientIsgroundinglineEnum));
	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));

	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));
	parameters->AddObject(new DoubleParam(TimeEnum,0.0));  //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,DiagnosticRequestedOutputsEnum);
	parameters->AddObject(new IntParam(DiagnosticNumRequestedOutputsEnum,numoutputs));
	if(numoutputs)parameters->AddObject(new IntVecParam(DiagnosticRequestedOutputsEnum,requestedoutputs,numoutputs));
	iomodel->DeleteData(requestedoutputs,DiagnosticRequestedOutputsEnum);

	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,PrognosticRequestedOutputsEnum);
	parameters->AddObject(new IntParam(PrognosticNumRequestedOutputsEnum,numoutputs));
	if(numoutputs)parameters->AddObject(new IntVecParam(PrognosticRequestedOutputsEnum,requestedoutputs,numoutputs));
	iomodel->DeleteData(requestedoutputs,PrognosticRequestedOutputsEnum);

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

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

	/*Go through all parameters, and convert units to SI: */
	parameters->UnitConversion(ExtToIuEnum);

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