/*!\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 "../ParseToolkitsOptionsx/ParseToolkitsOptionsx.h"
#include "./ModelProcessorx.h"

void CreateParameters(Parameters* parameters,IoModel* iomodel,char* rootpath,FILE* toolkitsoptionsfid,const int solution_type){

	int         i,j,m,k;
	int         numoutputs,materialtype,smb_model,basalforcing_model;
	char**      requestedoutputs = NULL;
	IssmDouble  time;
	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;

	/*Make sure current dataset is empty*/
	_assert_(parameters->Size()==0); 

	/*Copy some constants from iomodel */
	parameters->AddObject(iomodel->CopyConstantObject(DomainTypeEnum));
	parameters->AddObject(iomodel->CopyConstantObject(DomainDimensionEnum));
	parameters->AddObject(iomodel->CopyConstantObject(MeshElementtypeEnum));
	parameters->AddObject(iomodel->CopyConstantObject(SettingsOutputFrequencyEnum));
	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(TimesteppingInterpForcingsEnum));
	parameters->AddObject(iomodel->CopyConstantObject(SettingsLowmemEnum));
	parameters->AddObject(iomodel->CopyConstantObject(DebugProfilingEnum));
	parameters->AddObject(iomodel->CopyConstantObject(MeshAverageVertexConnectivityEnum));
	parameters->AddObject(iomodel->CopyConstantObject(ConstantsReferencetemperatureEnum));
	parameters->AddObject(iomodel->CopyConstantObject(SettingsWaitonlockEnum));
	parameters->AddObject(iomodel->CopyConstantObject(MeshNumberofelementsEnum));
	parameters->AddObject(iomodel->CopyConstantObject(MeshNumberofverticesEnum));
	parameters->AddObject(iomodel->CopyConstantObject(SettingsResultsOnNodesEnum));
	parameters->AddObject(iomodel->CopyConstantObject(SettingsIoGatherEnum));
	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(TransientIslevelsetEnum));
	parameters->AddObject(iomodel->CopyConstantObject(TransientIsdamageevolutionEnum));
	parameters->AddObject(iomodel->CopyConstantObject(TransientIshydrologyEnum));
	parameters->AddObject(iomodel->CopyConstantObject(MaterialsRheologyLawEnum));
	parameters->AddObject(iomodel->CopyConstantObject(AutodiffIsautodiffEnum));
	parameters->AddObject(iomodel->CopyConstantObject(QmuIsdakotaEnum));
	parameters->AddObject(iomodel->CopyConstantObject(InversionIscontrolEnum));
	parameters->AddObject(iomodel->CopyConstantObject(InversionTypeEnum));
	parameters->AddObject(iomodel->CopyConstantObject(GiaCrossSectionShapeEnum));

	/*For stress balance only*/
	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsFSEnum));
	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRiftPenaltyThresholdEnum));
	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceMaxiterEnum));
	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRestolEnum));
	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceReltolEnum));
	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceAbstolEnum));
	if(iomodel->domaintype==Domain3DEnum)
	 parameters->AddObject(iomodel->CopyConstantObject(MeshNumberoflayersEnum));

	/*Surface mass balance parameters*/
	parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsEnum));
	iomodel->Constant(&smb_model,SurfaceforcingsEnum);
	switch(smb_model){
		case SMBEnum:
			/*Nothing to add to parameters*/
			break;
		case SMBpddEnum:
			parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIsdelta18oEnum));
			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);
			}
			break;
		case SMBgradientsEnum:
			/*Nothing to add to parameters*/
			break;
		case SMBhenningEnum:
			/*Nothing to add to parameters*/
			break;
		case SMBcomponentsEnum:
			/*Nothing to add to parameters*/
			break;
		case SMBmeltcomponentsEnum:
			/*Nothing to add to parameters*/
			break;
		default:
			_error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
	}

	/*Basal forcing parameters*/
	parameters->AddObject(iomodel->CopyConstantObject(BasalforcingsEnum));
	iomodel->Constant(&basalforcing_model,BasalforcingsEnum);
	switch(basalforcing_model){
		case FloatingMeltRateEnum:
			/*Nothing to add to parameters*/
			break;
		case LinearFloatingMeltRateEnum:
			parameters->AddObject(iomodel->CopyConstantObject(BasalforcingsDeepwaterMeltingRateEnum));
			parameters->AddObject(iomodel->CopyConstantObject(BasalforcingsDeepwaterElevationEnum));
			parameters->AddObject(iomodel->CopyConstantObject(BasalforcingsUpperwaterElevationEnum));
			break;
		default:
			_error_("Basal forcing model "<<EnumToStringx(smb_model)<<" not supported yet");
	}

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

	iomodel->Constant(&time,TimesteppingStartTimeEnum);
	parameters->AddObject(new DoubleParam(TimeEnum,time));  
	parameters->AddObject(new IntParam(StepEnum,1));  

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

	/*Requested outputs */
	iomodel->FetchData(&requestedoutputs,&numoutputs,TransientRequestedOutputsEnum);
	parameters->AddObject(new IntParam(TransientNumRequestedOutputsEnum,numoutputs));
	if(numoutputs)parameters->AddObject(new StringArrayParam(TransientRequestedOutputsEnum,requestedoutputs,numoutputs));
	iomodel->DeleteData(&requestedoutputs,numoutputs,TransientRequestedOutputsEnum);

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

	iomodel->Constant(&materialtype,MaterialsEnum);
	if(materialtype==MatdamageiceEnum){
		iomodel->FetchData(&requestedoutputs,&numoutputs,DamageEvolutionRequestedOutputsEnum);
		parameters->AddObject(new IntParam(DamageEvolutionNumRequestedOutputsEnum,numoutputs));
		if(numoutputs)parameters->AddObject(new StringArrayParam(DamageEvolutionRequestedOutputsEnum,requestedoutputs,numoutputs));
		iomodel->DeleteData(&requestedoutputs,numoutputs,DamageEvolutionRequestedOutputsEnum);
	}

	/*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];
			_assert_(temp_n==5);

			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: */
	CreateParametersControl(parameters,iomodel,solution_type);

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

	/*Now, deal with toolkits options, which need to be put into the parameters dataset: */
	ParseToolkitsOptionsx(parameters,toolkitsoptionsfid);

 	#ifdef _HAVE_ADOLC_
	CreateParametersAutodiff(parameters,iomodel);
	#endif
}
