/*!\file:  Kriging.cpp
 * \brief  "c" core code for Kriging
 */ 

#include "./Krigingx.h"
#include "../../shared/shared.h"
#include "../../include/include.h"
#include "../../toolkits/toolkits.h"
#include "../../classes/objects/objects.h"
#include "../../Container/Container.h"
#include "../../io/io.h"

/*FUNCTION pKrigingx{{{*/
int pKrigingx(double** ppredictions,double **perror,double* obs_x, double* obs_y, double* obs_list, int obs_length,double* x_interp,double* y_interp,int n_interp,Options* options){

	int num_procs;
	int my_rank;

	/*output*/
	double *predictions = NULL;
	double *error       = NULL;

	/*Intermediaries*/
	int           mindata,maxdata;
	double        radius;
	char         *output       = NULL;
	Variogram    *variogram    = NULL;
	Observations *observations = NULL;

	/*Get my_rank: */
	my_rank=IssmComm::GetRank();
	num_procs=IssmComm::GetSize();

	/*Get some Options*/
	options->Get(&radius,"searchradius",0.);
	options->Get(&mindata,"mindata",1);
	options->Get(&maxdata,"maxdata",50);

	/*Process observation dataset*/
	observations=new Observations(obs_list,obs_x,obs_y,obs_length,options);

	/*Allocate output*/
	predictions =xNewZeroInit<double>(n_interp);
	error       =xNewZeroInit<double>(n_interp);

	/*Get output*/
	options->Get(&output,"output",(char*)"prediction");

	if(strcmp(output,"quadtree")==0){
		observations->QuadtreeColoring(predictions,x_interp,y_interp,n_interp);
	}
	else if(strcmp(output,"variomap")==0){
		observations->Variomap(predictions,x_interp,n_interp);
	}
	else if(strcmp(output,"prediction")==0){

		/*Process Variogram*/
		ProcessVariogram(&variogram,options);

		/*partition loop across threads: */
		for(int idx=my_rank;idx<n_interp;idx+=num_procs){
			_pprintLine_("      interpolation progress: "<<setw(6)<<setprecision(2)<<double(idx)/double(n_interp)*100<<"%");
			observations->InterpolationKriging(&predictions[idx],&error[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,variogram);
		}
		_pprintLine_("      interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%");

#ifdef _HAVE_MPI_
		double *sumpredictions =xNew<double>(n_interp);
		double *sumerror       =xNew<double>(n_interp);
		MPI_Allreduce(predictions,sumpredictions,n_interp,MPI_DOUBLE,MPI_SUM,IssmComm::GetComm());
		MPI_Allreduce(error,sumerror,n_interp,MPI_DOUBLE,MPI_SUM,IssmComm::GetComm());
		xDelete<double>(error); error=sumerror;
		xDelete<double>(predictions); predictions=sumpredictions;
#endif
	}
	else if(strcmp(output,"nearestneighbor")==0){

		/*partition loop across threads: */
		for(int idx=my_rank;idx<n_interp;idx+=num_procs){
			_pprintLine_("      interpolation progress: "<<setw(6)<<setprecision(2)<<double(idx)/double(n_interp)*100<<"%");
			observations->InterpolationNearestNeighbor(&predictions[idx],x_interp[idx],y_interp[idx],radius);
		}
		_pprintLine_("      interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%");

#ifdef _HAVE_MPI_
		double *sumpredictions =xNew<double>(n_interp);
		MPI_Allreduce(predictions,sumpredictions,n_interp,MPI_DOUBLE,MPI_SUM,IssmComm::GetComm());
		xDelete<double>(predictions); predictions=sumpredictions;
#endif
	}
	else if(strcmp(output,"idw")==0){
		double power;
		options->Get(&power,"power",2.);

		/*partition loop across threads: */
		for(int idx=my_rank;idx<n_interp;idx+=num_procs){
			_pprintLine_("      interpolation progress: "<<setw(6)<<setprecision(2)<<double(idx)/double(n_interp)*100<<"%");
			observations->InterpolationIDW(&predictions[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,power);
		}
		_pprintLine_("      interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%");

#ifdef _HAVE_MPI_
		double *sumpredictions =xNew<double>(n_interp);
		MPI_Allreduce(predictions,sumpredictions,n_interp,MPI_DOUBLE,MPI_SUM,IssmComm::GetComm());
		xDelete<double>(predictions); predictions=sumpredictions;
#endif
	}
	else{
		_error_("output '" << output << "' not supported yet");
	}

	/*clean-up and Assign output pointer*/
	delete variogram;
	delete observations;
	xDelete<char>(output);
	*ppredictions = predictions;
	*perror       = error;
	return 1;
}/*}}}*/
void ProcessVariogram(Variogram **pvariogram,Options* options){/*{{{*/

	/*Intermediaries*/
	Variogram* variogram = NULL;
	char      *model     = NULL;

	if(options->GetOption("model")){
		options->Get(&model,"model");
		if     (strcmp(model,"gaussian")==0)    variogram = new GaussianVariogram(options);
		else if(strcmp(model,"exponential")==0) variogram = new ExponentialVariogram(options);
		else if(strcmp(model,"spherical")==0)   variogram = new SphericalVariogram(options);
		else if(strcmp(model,"power")==0)       variogram = new PowerVariogram(options);
		else _error_("variogram " << model << " not supported yet (list of supported variogram: gaussian, exponential, spherical and power)");
	}
	else variogram = new GaussianVariogram(options);

	/*Assign output pointer*/
	xDelete<char>(model);
	*pvariogram = variogram;
}/*}}}*/
