Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/Container/Options.h =================================================================== --- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/Container/Options.h (revision 12291) +++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/Container/Options.h (revision 12292) @@ -21,6 +21,8 @@ Option* GetOption(const char* name); void Get(double* pvalue,const char* name); void Get(double* pvalue,const char* name,double default_value); + void Get(int* pvalue,const char* name); + void Get(int* pvalue,const char* name,int default_value); void Get(bool* pvalue,const char* name); void Get(bool* pvalue,const char* name,bool default_value); void Get(char** pvalue,const char* name); Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/Container/Observations.cpp =================================================================== --- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/Container/Observations.cpp (revision 12291) +++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/Container/Observations.cpp (revision 12292) @@ -93,49 +93,73 @@ /*Methods*/ /*FUNCTION Observations::ObservationList{{{*/ -void Observations::ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double range){ +void Observations::ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double radius,int maxdata){ /*Output and Intermediaries*/ - int nobs,i,index; + bool stop; + int nobs,i,j,k,n,counter; + double h2,radius2; + int *indices = NULL; + double *dists = NULL; double *x = NULL; double *y = NULL; double *obs = NULL; Observation *observation = NULL; - int *indices = NULL; - /*Treat range*/ - if(range==0){ - range=this->quadtree->root->length; - } + /*Compute radius square*/ + if(radius==0) radius=this->quadtree->root->length; + radius2 = radius*radius; - /*Find all observations that are in range*/ - nobs=0; - while(nobs==0){ - xfree((void**)&indices); - this->quadtree->RangeSearch(&indices,&nobs,x_interp,y_interp,range); - if(!nobs){ - /*No observation found, double range*/ - range=2*range; + /*Find all observations that are in radius*/ + indices = (int*)xmalloc(this->Size()*sizeof(int)); + dists = (double*)xmalloc(this->Size()*sizeof(double)); + nobs = 0; + + for (i=0;iSize();i++){ + observation=(Observation*)this->GetObjectByOffset(i); + h2 = (observation->x-x_interp)*(observation->x-x_interp) + (observation->y-y_interp)*(observation->y-y_interp); + + if(nobs==maxdata && h2>radius2) continue; + if(nobs<=maxdata){ + indices[nobs] = i; + dists[nobs] = h2; + nobs++; } - } - /*Check that we do not have more than 50 observations*/ - //while(nobs>50){ - // range=range/2; - // xfree((void**)&indices); - // this->quadtree->RangeSearch(&indices,&nobs,x_interp,y_interp,range); - //} + if(nobs==1) continue; - if(!nobs) _error_("Unexpected error: no observation within current range and more than 50 observations if doubling range"); + /*Sort all dists up to now*/ + n=nobs-1; + stop = false; + for(k=0;kGetObjectByOffset(indices[i]); - observation->WriteXYObs(&x[i],&y[i],&obs[i]); + /*Loop over all observations and fill in x, y and obs*/ + for (i=0;iGetObjectByOffset(indices[i]); + observation->WriteXYObs(&x[i],&y[i],&obs[i]); + } } /*Assign output pointer*/ Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/Container/Observations.h =================================================================== --- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/Container/Observations.h (revision 12291) +++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/Container/Observations.h (revision 12292) @@ -22,7 +22,7 @@ ~Observations(); /*Methods*/ - void ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double range); + void ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double radius,int maxdata); void QuadtreeColoring(double* A,double *x,double *y,int n); void Variomap(double* gamma,double *x,int n); Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/Container/Options.cpp =================================================================== --- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/Container/Options.cpp (revision 12291) +++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/Container/Options.cpp (revision 12292) @@ -66,6 +66,44 @@ return 1; } /*}}}*/ +/*FUNCTION Options::Get(int* pvalue, char* name){{{1*/ +void Options::Get(int* pvalue,const char* name){ + + vector::iterator object; + Option* option=NULL; + + /*Get option*/ + option=GetOption(name); + + /*If the pointer is not NULL, the option has been found*/ + if(option){ + option->Get(pvalue); + } + /*Else, the Option does not exist, no default provided*/ + else{ + _error_("option of name \"%s\" not found, and no default value has been provided",name); + } +} +/*}}}*/ +/*FUNCTION Options::Get(int* pvalue, char* name,int default_value){{{1*/ +void Options::Get(int* pvalue,const char* name,int default_value){ + + vector::iterator object; + Option* option=NULL; + + /*Get option*/ + option=GetOption(name); + + /*If the pointer is not NULL, the option has been found*/ + if(option){ + option->Get(pvalue); + } + /*Else, the Option does not exist, a default is provided here*/ + else{ + *pvalue=default_value; + } +} +/*}}}*/ /*FUNCTION Options::Get(double* pvalue, char* name){{{1*/ void Options::Get(double* pvalue,const char* name){ Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/shared/Threads/LaunchThread.cpp =================================================================== --- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/shared/Threads/LaunchThread.cpp (revision 12291) +++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/shared/Threads/LaunchThread.cpp (revision 12292) @@ -66,5 +66,4 @@ function((void*)&handle); #endif - } Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp =================================================================== --- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp (revision 12291) +++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp (revision 12292) @@ -23,7 +23,8 @@ double *error = NULL; /*Intermediaries*/ - double range; + int mindata,maxdata; + double radius; char *output = NULL; Variogram *variogram = NULL; Observations *observations = NULL; @@ -37,7 +38,9 @@ /*Get Variogram from Options*/ ProcessVariogram(&variogram,options); - options->Get(&range,"searchrange",0.); + 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); @@ -61,7 +64,9 @@ gate.n_interp = n_interp; gate.x_interp = x_interp; gate.y_interp = y_interp; - gate.range = range; + gate.radius = radius; + gate.mindata = mindata; + gate.maxdata = maxdata; gate.variogram = variogram; gate.observations = observations; gate.predictions = predictions; @@ -105,7 +110,9 @@ int n_interp = gate->n_interp; double *x_interp = gate->x_interp; double *y_interp = gate->y_interp; - double range = gate->range; + double radius = gate->radius; + int mindata = gate->mindata; + int maxdata = gate->maxdata; Variogram *variogram = gate->variogram; Observations *observations = gate->observations; double *predictions = gate->predictions; @@ -136,10 +143,10 @@ printf("\r interpolation progress: %5.2lf%%",localpercent); /*Get list of observations for current point*/ - observations->ObservationList(&x,&y,&obs,&n_obs,x_interp[idx],y_interp[idx],range); - if(n_obs==0){ - predictions[idx] = -999.; - error[idx] = -999.; + observations->ObservationList(&x,&y,&obs,&n_obs,x_interp[idx],y_interp[idx],radius,maxdata); + if(n_obsNumEl()!=1){ + _error_("option \"%s\" has %i elements and cannot return a single int",this->name,this->NumEl()); + } + + /*Assign output pointer*/ + *pvalue=(int)this->values[0]; +} +/*}}}*/ /*FUNCTION OptionDouble::Get(double* pvalue) {{{1*/ void OptionDouble::Get(double* pvalue){ Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/OptionDouble.h =================================================================== --- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/OptionDouble.h (revision 12291) +++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/OptionDouble.h (revision 12292) @@ -38,6 +38,7 @@ int NumEl(); int NDims(); int* Size(); + void Get(int* pvalue); void Get(double* pvalue); void Get(bool* pvalue){ _error_("An OptionDouble object cannot return a bool");}; void Get(char** pvalue){ _error_("An OptionDouble object cannot return a string");}; Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/OptionLogical.h =================================================================== --- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/OptionLogical.h (revision 12291) +++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/OptionLogical.h (revision 12292) @@ -38,6 +38,7 @@ int NumEl(); int NDims(); int* Size(); + void Get(int* pvalue){_error_("An OptionLogical object cannot return a int");}; void Get(double* pvalue){_error_("An OptionLogical object cannot return a double");}; void Get(bool* pvalue); void Get(char** pvalue){ _error_("An OptionLogical object cannot return a string");}; Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/OptionChar.h =================================================================== --- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/OptionChar.h (revision 12291) +++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/OptionChar.h (revision 12292) @@ -38,6 +38,7 @@ int NumEl(); int NDims(); int* Size(); + void Get(int* pvalue){_error_("An OptionChar object cannot return a int");}; void Get(double* pvalue){_error_("An OptionChar object cannot return a double");}; void Get(bool* pvalue){ _error_("An OptionChar object cannot return a bool");}; void Get(char** pvalue); Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/Option.h =================================================================== --- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/Option.h (revision 12291) +++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/Option.h (revision 12292) @@ -41,6 +41,7 @@ virtual int NumEl()=0; virtual int NDims()=0; virtual int* Size()=0; + virtual void Get(int* pvalue)=0; virtual void Get(double* pvalue)=0; virtual void Get(bool* pvalue)=0; virtual void Get(char** pvalue)=0; Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/OptionCell.h =================================================================== --- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/OptionCell.h (revision 12291) +++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/OptionCell.h (revision 12292) @@ -38,6 +38,7 @@ int NumEl(); int NDims(); int* Size(); + void Get(int* pvalue){_error_("An OptionCell object cannot return a int");}; void Get(double* pvalue){_error_("An OptionCell object cannot return a double");}; void Get(bool* pvalue){ _error_("An OptionCell object cannot return a bool");}; void Get(char** pvalue){ _error_("An OptionCell object cannot return a string");}; Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Kriging/GaussianVariogram.cpp =================================================================== --- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Kriging/GaussianVariogram.cpp (revision 12291) +++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Kriging/GaussianVariogram.cpp (revision 12292) @@ -74,7 +74,7 @@ a = 1./3.; cova = (sill-nugget)*exp(-h2/(a*range*range)); - //if(h2<0.0000001) cova = sill; + if(h2<0.0000001) cova = sill; return cova; }