/*!\file TransientInput2.c
 * \brief: implementation of the TransientInput2 object
 */
/*Headers*/
#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include "./TransientInput2.h"
#include "./TriaInput2.h"
#include "./PentaInput2.h"
#include "../../shared/shared.h"
#include "../Params/Parameters.h"

/*TransientInput2 constructors and destructor*/
TransientInput2::TransientInput2(){/*{{{*/

	enum_type=UNDEF;
	inputs=NULL;
	this->numtimesteps=0;
	this->parameters=NULL;
	this->timesteps=NULL;

	this->current_input=NULL;
	this->current_step=-1;

}
/*}}}*/
TransientInput2::TransientInput2(int in_enum_type,int nbe,int nbv,IssmDouble* timesin,int N){/*{{{*/

	/*Set Enum*/
	this->enum_type=in_enum_type;
	this->numberofelements_local = nbe;
	this->numberofvertices_local = nbv;

	/*Allocate values and timesteps, and copy: */
	_assert_(N>=0 && N<1e6);
	this->numtimesteps=N;
	if(N>0){
		this->timesteps=xNew<IssmDouble>(N);
		xMemCpy(this->timesteps,timesin,N);

		this->inputs     = xNew<Input2*>(N);
		for(int i=0;i<N;i++) this->inputs[i] = NULL;
	}
	else{
		this->timesteps=0;
		this->inputs   =0;
	}
	this->parameters = NULL;
	this->current_input=NULL;
	this->current_step=-1;
}
/*}}}*/
TransientInput2::~TransientInput2(){/*{{{*/

	for(int i=0;i<this->numtimesteps;i++){
		delete this->inputs[i];
	}
	xDelete<Input2*>(this->inputs);
	xDelete<IssmDouble>(this->timesteps);

	if(this->current_input) delete this->current_input;
}
/*}}}*/

/*Object virtual functions definitions:*/
Input2* TransientInput2::copy() {/*{{{*/

	TransientInput2* output=NULL;

	output = new TransientInput2();
	output->enum_type=this->enum_type;
	output->numtimesteps=this->numtimesteps;
	if(this->numtimesteps>0){
		output->timesteps=xNew<IssmDouble>(this->numtimesteps);
		xMemCpy(output->timesteps,this->timesteps,this->numtimesteps);
		output->inputs = xNew<Input2*>(this->numtimesteps);
		for(int i=0;i<this->numtimesteps;i++){
			output->inputs[i] = this->inputs[i]->copy();
		}
	}
	output->parameters=this->parameters;

	return output;
}/*}}}*/
void TransientInput2::DeepEcho(void){/*{{{*/

	int i;

	_printf_("TransientInput2:\n");
	_printf_("   enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")\n");
	_printf_("   numtimesteps: " << this->numtimesteps << "\n");
	_printf_("---inputs: \n");
	for(i=0;i<this->numtimesteps;i++){
		_printf_("   time: " << this->timesteps[i]<<"  ");
		if(this->inputs[i]) this->inputs[i]->Echo();
		else                _printf_(" NOT SET! \n");
	}
}
/*}}}*/
void TransientInput2::Configure(Parameters* params){/*{{{*/
	this->parameters=params;
}
/*}}}*/
void TransientInput2::Echo(void){/*{{{*/
	this->DeepEcho();
}
/*}}}*/
int  TransientInput2::Id(void){ return -1; }/*{{{*/
/*}}}*/
void TransientInput2::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/

	if (marshall_direction == MARSHALLING_BACKWARD){
		_error_("not implmented");
		//inputs = new Inputs();
	}

	MARSHALLING_ENUM(TransientInput2Enum);

	MARSHALLING(enum_type);
	MARSHALLING(numtimesteps);
	MARSHALLING_DYNAMIC(this->timesteps,IssmDouble,numtimesteps);
	//inputs->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
	_error_("not implemented");
}
/*}}}*/
int  TransientInput2::ObjectEnum(void){/*{{{*/

	return TransientInput2Enum;

}
/*}}}*/

/*Intermediary*/
void TransientInput2::AddTriaTimeInput(IssmDouble time,int numindices,int* indices,IssmDouble* values_in,int interp_in){/*{{{*/


	/*Check whether this is the last time step that we have*/
	if(this->numtimesteps){
		if(fabs(this->timesteps[this->numtimesteps-1]-time)<1.0e-5){
			this->AddTriaTimeInput(this->numtimesteps-1,numindices,indices,values_in,interp_in);
			return;
		}
	}

	/*This is a new time step! we need to add it to the list*/
	if(this->numtimesteps>0 && time<this->timesteps[this->numtimesteps-1]) _error_("timestep values must increase sequentially");

	IssmDouble *old_timesteps = NULL;
	Input2    **old_inputs    = NULL;
	if (this->numtimesteps > 0){
		old_timesteps=xNew<IssmDouble>(this->numtimesteps);
		xMemCpy(old_timesteps,this->timesteps,this->numtimesteps);
		xDelete<IssmDouble>(this->timesteps);
		old_inputs=xNew<Input2*>(this->numtimesteps);
		xMemCpy(old_inputs,this->inputs,this->numtimesteps);
		xDelete<Input2*>(this->inputs);
	}

	this->numtimesteps=this->numtimesteps+1;
	this->timesteps=xNew<IssmDouble>(this->numtimesteps);
	this->inputs   = xNew<Input2*>(this->numtimesteps);

	if (this->numtimesteps > 1){
		xMemCpy(this->inputs,old_inputs,this->numtimesteps-1);
		xMemCpy(this->timesteps,old_timesteps,this->numtimesteps-1);
		xDelete(old_timesteps);
		xDelete<Input2*>(old_inputs);
	}

	/*go ahead and plug: */
	this->timesteps[this->numtimesteps-1] = time;
	this->inputs[this->numtimesteps-1]    = NULL;
	this->AddTriaTimeInput(this->numtimesteps-1,numindices,indices,values_in,interp_in);

}
/*}}}*/
void TransientInput2::AddPentaTimeInput(IssmDouble time,int numindices,int* indices,IssmDouble* values_in,int interp_in){/*{{{*/

	_error_("not implemented yet, look at TransientInput2::AddTriaTimeInput");

}
/*}}}*/
void TransientInput2::AddTriaTimeInput(int step,int numindices,int* indices,IssmDouble* values_in,int interp_in){/*{{{*/

	_assert_(step>=0 && step<this->numtimesteps);

	/*Create it if necessary*/
	if(this->inputs[step]){
		if(this->inputs[step]->ObjectEnum()!=TriaInput2Enum) _error_("cannot add Element values to a "<<EnumToStringx(this->inputs[step]->ObjectEnum()));
	}
	else{
		this->inputs[step] = new TriaInput2(this->numberofelements_local,this->numberofvertices_local,interp_in);
	}

	/*Set input*/
	TriaInput2* input = xDynamicCast<TriaInput2*>(this->inputs[step]);
	input->SetInput(interp_in,numindices,indices,values_in);

}
/*}}}*/
void TransientInput2::AddPentaTimeInput(int step,int numindices,int* indices,IssmDouble* values_in,int interp_in){/*{{{*/

	_assert_(step>=0 && step<this->numtimesteps);

	/*Create it if necessary*/
	if(this->inputs[step]){
		if(this->inputs[step]->ObjectEnum()!=PentaInput2Enum) _error_("cannot add Element values to a "<<EnumToStringx(this->inputs[step]->ObjectEnum()));
	}
	else{
		this->inputs[step] = new PentaInput2(this->numberofelements_local,this->numberofvertices_local,interp_in);
	}

	/*Set input*/
	PentaInput2* input = xDynamicCast<PentaInput2*>(this->inputs[step]);
	input->SetInput(interp_in,numindices,indices,values_in);

}
/*}}}*/
void TransientInput2::GetAllTimes(IssmDouble** ptimesteps,int* pnumtimesteps){/*{{{*/

	if(ptimesteps){
		*ptimesteps=xNew<IssmDouble>(this->numtimesteps);
		xMemCpy(*ptimesteps,this->timesteps,this->numtimesteps);
	}
	if(pnumtimesteps){
		*pnumtimesteps = this->numtimesteps;
	}

}
/*}}}*/
TriaInput2* TransientInput2::GetTriaInput(){/*{{{*/

	IssmDouble time;
	this->parameters->FindParam(&time,TimeEnum);
	return this->GetTriaInput(time);

}
/*}}}*/
TriaInput2* TransientInput2::GetTriaInput(IssmDouble time){/*{{{*/

	/*Set current time input*/
	this->SetCurrentTimeInput(time);
	_assert_(this->current_input);

	/*Cast and return*/
	if(this->current_input->ObjectEnum()!=TriaInput2Enum){
		_error_("Cannot return a TriaInput2");
	}
	return xDynamicCast<TriaInput2*>(this->current_input);

}
/*}}}*/
TriaInput2* TransientInput2::GetTriaInput(IssmDouble start_time, IssmDouble end_time){/*{{{*/

	/*Set current time input*/
	this->SetAverageAsCurrentTimeInput(start_time,end_time);
	_assert_(this->current_input);

	/*Cast and return*/
	if(this->current_input->ObjectEnum()!=TriaInput2Enum){
		_error_("Cannot return a TriaInput2");
	}
	return xDynamicCast<TriaInput2*>(this->current_input);

}
/*}}}*/
TriaInput2* TransientInput2::GetTriaInput(int offset){/*{{{*/

	/*Check offset*/
	if(offset<0 || offset>this->numtimesteps-1){
		_error_("Cannot return input for offset "<<offset);
	}
	Input2* input = this->inputs[offset];

	/*Cast and return*/
	_assert_(input);
	if(input->ObjectEnum()!=TriaInput2Enum) _error_("Cannot return a TriaInput2");
	return xDynamicCast<TriaInput2*>(input);

}
/*}}}*/
PentaInput2* TransientInput2::GetPentaInput(){/*{{{*/

	IssmDouble time;
	this->parameters->FindParam(&time,TimeEnum);
	return this->GetPentaInput(time);
}
/*}}}*/
PentaInput2* TransientInput2::GetPentaInput(IssmDouble time){/*{{{*/

	/*Set current time input*/
	this->SetCurrentTimeInput(time);
	_assert_(this->current_input);

	/*Cast and return*/
	if(this->current_input->ObjectEnum()!=PentaInput2Enum){
		_error_("Cannot return a PentaInput2");
	}
	return xDynamicCast<PentaInput2*>(this->current_input);

}
/*}}}*/
PentaInput2* TransientInput2::GetPentaInput(int offset){/*{{{*/


	/*Check offset*/
	if(offset<0 || offset>this->numtimesteps-1){
		_error_("Cannot return input for offset "<<offset);
	}
	Input2* input = this->inputs[offset];

	/*Cast and return*/
	if(input->ObjectEnum()!=PentaInput2Enum) _error_("Cannot return a PentaInput2");
	return xDynamicCast<PentaInput2*>(input);

}
/*}}}*/
void TransientInput2::SetCurrentTimeInput(IssmDouble time){/*{{{*/

	/*First, recover current time from parameters: */
	bool linear_interp;
	this->parameters->FindParam(&linear_interp,TimesteppingInterpForcingsEnum);

	/*Figure step out*/
	int offset;
	if(!binary_search(&offset,time,this->timesteps,this->numtimesteps)){
		_error_("Input not found (is TransientInput sorted ?)");
	}

	if (offset==-1){

		/*get values for the first time: */
		_assert_(time<this->timesteps[0]);

		/*If already processed return*/
		if(this->current_step==0.) return;

		/*Prepare input*/
		if(this->current_input) delete this->current_input;
		this->current_step = 0.;
		this->current_input = this->inputs[0]->copy();

	}
	else if(offset==(this->numtimesteps-1) || !linear_interp){

		/*get values for the last time: */
		_assert_(time>=this->timesteps[offset]);

		/*If already processed return*/
		if(this->current_step==reCast<IssmDouble>(offset)) return;

		/*Prepare input*/
		if(this->current_input) delete this->current_input;
		this->current_step  = reCast<IssmDouble>(offset);
		this->current_input = this->inputs[offset]->copy();
	}
	else {

		/*Interpolate */
		_assert_(time>=this->timesteps[offset] && time<this->timesteps[offset+1]);

		/*get values between two times [offset:offset+1[, Interpolate linearly*/
		IssmDouble deltat=this->timesteps[offset+1]-this->timesteps[offset];
		IssmDouble this_step = reCast<IssmDouble>(offset) + (time - this->timesteps[offset])/deltat;

		/*If already processed return*/
		if(this->current_step>this_step-1.e-5 && this->current_step<this_step+1.e-5) return;

		/*Prepare input*/
		if(this->current_input) delete this->current_input;
		this->current_step = this_step;
		IssmDouble alpha2=(time-this->timesteps[offset])/deltat;
		IssmDouble alpha1=(1.0-alpha2);

		Input2* input1=this->inputs[offset];
		Input2* input2=this->inputs[offset+1];

		this->current_input = input1->copy();
		this->current_input->Scale(alpha1);
		this->current_input->AXPY(input2,alpha2);
	}

}/*}}}*/
void TransientInput2::SetAverageAsCurrentTimeInput(IssmDouble start_time,IssmDouble end_time){/*{{{*/

	IssmDouble  dt;
	IssmPDouble eps=1.0e-6;
	IssmDouble  dtsum=0;
	int         found,start_offset,end_offset;
	int         averaging_method = 0;

	/*go through the timesteps, and grab offset for start and end*/
	IssmDouble temp = start_time-eps;
	found=binary_search(&start_offset,temp,this->timesteps,this->numtimesteps);
	if(!found) _error_("Input not found (is TransientInput sorted ?)");
	temp = end_time+eps;
	found=binary_search(&end_offset,temp,this->timesteps,this->numtimesteps);
	if(!found) _error_("Input not found (is TransientInput sorted ?)");

	int offset=start_offset;
	if(this->current_input) delete this->current_input;
	while(offset < end_offset){
		if (offset==-1){
			dt=this->timesteps[0]-start_time;
			_assert_(start_time<this->timesteps[0]);
		}
		else if(offset==this->numtimesteps-1)dt=end_time-this->timesteps[offset];
		else if(offset==start_offset && this->timesteps[offset]<start_time) dt=this->timesteps[offset+1]-start_time;
		else if(offset==end_offset && this->timesteps[offset]>end_time) dt=end_time-this->timesteps[offset];
		else dt=this->timesteps[offset+1]-this->timesteps[offset];
		_assert_(dt>0.);

		Input2* stepinput=this->inputs[offset+1];

		switch(averaging_method){
			case 0: /*Arithmetic mean*/
				if(offset==start_offset){
					this->current_input = stepinput->copy();
					this->current_input->Scale(dt);
				}
				else{
					this->current_input->AXPY(stepinput,dt);
				}
				break;
			case 1: /*Geometric mean*/
				_error_("Geometric not implemented yet");
			case 2: /*Harmonic mean*/
				_error_("Harmonic not implemented yet");
			default:
				_error_("averaging method is not recognised");
		}
		dtsum+=dt;
		offset+=1;
	}
		/*Integration done, now normalize*/
	switch(averaging_method){
		case 0: //Arithmetic mean
			this->current_input->Scale(1/(dtsum));
			break;
		case 1: /*Geometric mean*/
			_error_("Geometric not implemented yet");
		case 2: /*Harmonic mean*/
			_error_("Harmonic not implemented yet");
		default:
			_error_("averaging method is not recognised");
	}
}/*}}}*/
IssmDouble  TransientInput2::GetTimeByOffset(int offset){/*{{{*/
	if(offset<0) offset=0;
	_assert_(offset<this->numtimesteps);
	return this->timesteps[offset];
}
/*}}}*/
int  TransientInput2::GetTimeInputOffset(IssmDouble time){/*{{{*/

	int offset;

	/*go through the timesteps, and figure out which interval we
	 *     *fall within. Then interpolate the values on this interval: */
	int found=binary_search(&offset,time,this->timesteps,this->numtimesteps);
	if(!found) _error_("Input not found (is TransientInput sorted ?)");

	return offset;
}
/*}}}*/
