/*!\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 "../../objects/objects.h"
#include "../../Container/Observations.h"
#include "../modules.h"
#ifdef _HAVE_GSL_
#include <gsl/gsl_linalg.h>
#endif
/*FUNCTION Krigingx{{{*/
int Krigingx(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){

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

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

	/*threading: */
	KrigingxThreadStruct gate;
	int num=1;
#ifdef _MULTITHREADING_
	num=_NUMTHREADS_;
#endif

	/*Get Variogram from Options*/
	ProcessVariogram(&variogram,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","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,"delaunay")==0){
		int nobs,nel;
		double *x     = NULL;
		double *y     = NULL;
		double *data  = NULL;
		int    *index = NULL;

		observations->ObservationList(&x,&y,&data,&nobs);

		_printLine_("Generation Delaunay Triangulation");
		BamgTriangulatex(&index,&nel,x,y,nobs);

		_printLine_("Interpolating");
		xDelete<double>(predictions);
		InterpFromMeshToMesh2dx(&predictions,index,x,y,nobs,nel,data,nobs,1,x_interp,y_interp,n_interp,options);
		xDelete<double>(x);
		xDelete<double>(y);
		xDelete<double>(data);
		xDelete<int>(index);
	}
	else if(strcmp(output,"nearestneighbor")==0){
		/*initialize thread parameters: */
		gate.n_interp     = n_interp;
		gate.x_interp     = x_interp;
		gate.y_interp     = y_interp;
		gate.radius       = radius;
		gate.mindata      = mindata;
		gate.maxdata      = maxdata;
		gate.variogram    = variogram;
		gate.observations = observations;
		gate.predictions  = predictions;
		gate.error        = error;
		gate.percent      = xNewZeroInit<double>(num);

		/*launch the thread manager with Krigingxt as a core: */
		LaunchThread(NearestNeighbort,(void*)&gate,num);
		_printLine_("\r      interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%");
		xDelete<double>(gate.percent);
	}
	else if(strcmp(output,"idw")==0){ //Inverse distance weighting
		/*initialize thread parameters: */
		gate.n_interp     = n_interp;
		gate.x_interp     = x_interp;
		gate.y_interp     = y_interp;
		gate.radius       = radius;
		gate.mindata      = mindata;
		gate.maxdata      = maxdata;
		gate.variogram    = variogram;
		gate.observations = observations;
		gate.predictions  = predictions;
		gate.error        = error;
		gate.percent      = xNewZeroInit<double>(num);

		/*launch the thread manager with Krigingxt as a core: */
		LaunchThread(idwt,(void*)&gate,num);
		_printLine_("\r      interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%");
		xDelete<double>(gate.percent);
	}
	else if(strcmp(output,"prediction")==0){

		/*initialize thread parameters: */
		gate.n_interp     = n_interp;
		gate.x_interp     = x_interp;
		gate.y_interp     = y_interp;
		gate.radius       = radius;
		gate.mindata      = mindata;
		gate.maxdata      = maxdata;
		gate.variogram    = variogram;
		gate.observations = observations;
		gate.predictions  = predictions;
		gate.error        = error;
		gate.percent      = xNewZeroInit<double>(num);

		/*launch the thread manager with Krigingxt as a core: */
		LaunchThread(Krigingxt,(void*)&gate,num);
		_printLine_("\r      interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%");
		xDelete<double>(gate.percent);
	}
	else{
		_error2_("output '" << output << "' not supported yet");
	}

	/*clean-up and Assign output pointer*/
	delete variogram;
	delete observations;
	xDelete<char>(output);
	*ppredictions = predictions;
	*perror       = error;
	return 1;
}/*}}}*/
/*FUNCTION Krigingxt{{{*/
void* Krigingxt(void* vpthread_handle){

	/*gate variables :*/
	KrigingxThreadStruct *gate        = NULL;
	pthread_handle       *handle      = NULL;
	int my_thread;
	int num_threads;
	int i0,i1;

	/*recover handle and gate: */
	handle      = (pthread_handle*)vpthread_handle;
	gate        = (KrigingxThreadStruct*)handle->gate;
	my_thread   = handle->id;
	num_threads = handle->num;

	/*recover parameters :*/
	int           n_interp     = gate->n_interp;
	double       *x_interp     = gate->x_interp;
	double       *y_interp     = gate->y_interp;
	double        radius       = gate->radius;
	int           mindata      = gate->mindata;
	int           maxdata      = gate->maxdata;
	Variogram    *variogram    = gate->variogram;
	Observations *observations = gate->observations;
	double       *predictions  = gate->predictions;
	double       *error        = gate->error;
	double       *percent      = gate->percent;

	/*Intermediaries*/
	double        localpercent;

	/*partition loop across threads: */
	PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
	for(int idx=i0;idx<i1;idx++){

		/*Print info*/
		percent[my_thread]=double(idx-i0)/double(i1-i0)*100.;
		localpercent=percent[0];
		for(int i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]);
		if(my_thread==0) _printString_("\r      interpolation progress: "<<setw(6)<<setprecision(2)<<localpercent<<"%");

		/*Kriging interpolation*/
		observations->InterpolationKriging(&predictions[idx],&error[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,variogram);
	}

	return NULL;
}/*}}}*/
/*FUNCTION NearestNeighbort{{{*/
void* NearestNeighbort(void* vpthread_handle){

	/*gate variables :*/
	KrigingxThreadStruct *gate        = NULL;
	pthread_handle       *handle      = NULL;
	int my_thread;
	int num_threads;
	int i0,i1;

	/*recover handle and gate: */
	handle      = (pthread_handle*)vpthread_handle;
	gate        = (KrigingxThreadStruct*)handle->gate;
	my_thread   = handle->id;
	num_threads = handle->num;

	/*recover parameters :*/
	int           n_interp     = gate->n_interp;
	double       *x_interp     = gate->x_interp;
	double       *y_interp     = gate->y_interp;
	double        radius       = gate->radius;
	int           mindata      = gate->mindata;
	int           maxdata      = gate->maxdata;
	Variogram    *variogram    = gate->variogram;
	Observations *observations = gate->observations;
	double       *predictions  = gate->predictions;
	double       *error        = gate->error;
	double       *percent      = gate->percent;

	/*Intermediaries*/
	int           i;
	double        localpercent;

	/*partition loop across threads: */
	PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
	for(int idx=i0;idx<i1;idx++){

		/*Print info*/
		percent[my_thread]=double(idx-i0)/double(i1-i0)*100.;
		localpercent=percent[0];
		for(i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]);
		if(my_thread==0) _printString_("\r      interpolation progress: "<<setw(6)<<setprecision(2)<<localpercent<<"%");

		observations->InterpolationNearestNeighbor(&predictions[idx],x_interp[idx],y_interp[idx],radius);
	}

	return NULL;
}/*}}}*/
/*FUNCTION idwt{{{*/
void* idwt(void* vpthread_handle){

	/*gate variables :*/
	KrigingxThreadStruct *gate        = NULL;
	pthread_handle       *handle      = NULL;
	int my_thread;
	int num_threads;
	int i0,i1;

	/*recover handle and gate: */
	handle      = (pthread_handle*)vpthread_handle;
	gate        = (KrigingxThreadStruct*)handle->gate;
	my_thread   = handle->id;
	num_threads = handle->num;

	/*recover parameters :*/
	int           n_interp     = gate->n_interp;
	double       *x_interp     = gate->x_interp;
	double       *y_interp     = gate->y_interp;
	double        radius       = gate->radius;
	int           mindata      = gate->mindata;
	int           maxdata      = gate->maxdata;
	Variogram    *variogram    = gate->variogram;
	Observations *observations = gate->observations;
	double       *predictions  = gate->predictions;
	double       *error        = gate->error;
	double       *percent      = gate->percent;

	/*Intermediaries*/
	double localpercent;
	double  power = 2.;

	/*partition loop across threads: */
	PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
	for(int idx=i0;idx<i1;idx++){

		/*Print info*/
		percent[my_thread]=double(idx-i0)/double(i1-i0)*100.;
		localpercent=percent[0];
		for(int i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]);
		if(my_thread==0) _printString_("\r      interpolation progress: "<<setw(6)<<setprecision(2)<<localpercent<<"%");

		observations->InterpolationIDW(&predictions[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,power);
	}
	return NULL;
}/*}}}*/

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 _error2_("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;
}/*}}}*/
