Index: /issm/trunk-jpl/src/c/Container/Observations.cpp =================================================================== --- /issm/trunk-jpl/src/c/Container/Observations.cpp (revision 12291) +++ /issm/trunk-jpl/src/c/Container/Observations.cpp (revision 12292) @@ -94,47 +94,71 @@ /*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; - } - - /*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; - } - } - /*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) _error_("Unexpected error: no observation within current range and more than 50 observations if doubling range"); - - /*Allocate vectors*/ - x = (double*)xmalloc(nobs*sizeof(double)); - y = (double*)xmalloc(nobs*sizeof(double)); - obs = (double*)xmalloc(nobs*sizeof(double)); - - /*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]); + + /*Compute radius square*/ + if(radius==0) radius=this->quadtree->root->length; + radius2 = radius*radius; + + /*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++; + } + if(nobs==1) continue; + + /*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]); + } } Index: /issm/trunk-jpl/src/c/Container/Observations.h =================================================================== --- /issm/trunk-jpl/src/c/Container/Observations.h (revision 12291) +++ /issm/trunk-jpl/src/c/Container/Observations.h (revision 12292) @@ -23,5 +23,5 @@ /*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: /issm/trunk-jpl/src/c/Container/Options.cpp =================================================================== --- /issm/trunk-jpl/src/c/Container/Options.cpp (revision 12291) +++ /issm/trunk-jpl/src/c/Container/Options.cpp (revision 12292) @@ -65,4 +65,42 @@ 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; + } } /*}}}*/ Index: /issm/trunk-jpl/src/c/Container/Options.h =================================================================== --- /issm/trunk-jpl/src/c/Container/Options.h (revision 12291) +++ /issm/trunk-jpl/src/c/Container/Options.h (revision 12292) @@ -22,4 +22,6 @@ 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); Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp (revision 12291) +++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp (revision 12292) @@ -24,5 +24,6 @@ /*Intermediaries*/ - double range; + int mindata,maxdata; + double radius; char *output = NULL; Variogram *variogram = NULL; @@ -38,5 +39,7 @@ /*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*/ @@ -62,5 +65,7 @@ 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; @@ -106,5 +111,7 @@ 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; @@ -137,8 +144,8 @@ /*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: /issm/trunk-jpl/src/c/objects/Options/OptionDouble.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Options/OptionDouble.h (revision 12291) +++ /issm/trunk-jpl/src/c/objects/Options/OptionDouble.h (revision 12292) @@ -39,4 +39,5 @@ int NDims(); int* Size(); + void Get(int* pvalue); void Get(double* pvalue); void Get(bool* pvalue){ _error_("An OptionDouble object cannot return a bool");}; Index: /issm/trunk-jpl/src/c/objects/Options/OptionLogical.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Options/OptionLogical.h (revision 12291) +++ /issm/trunk-jpl/src/c/objects/Options/OptionLogical.h (revision 12292) @@ -39,4 +39,5 @@ 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); Index: /issm/trunk-jpl/src/c/objects/Options/OptionStruct.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Options/OptionStruct.h (revision 12291) +++ /issm/trunk-jpl/src/c/objects/Options/OptionStruct.h (revision 12292) @@ -39,4 +39,5 @@ int NDims(); int* Size(); + void Get(int* pvalue){_error_("An OptionStruct object cannot return a int");}; void Get(double* pvalue){_error_("An OptionStruct object cannot return a double");}; void Get(bool* pvalue){ _error_("An OptionStruct object cannot return a bool");}; Index: /issm/trunk-jpl/src/c/shared/Threads/LaunchThread.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Threads/LaunchThread.cpp (revision 12291) +++ /issm/trunk-jpl/src/c/shared/Threads/LaunchThread.cpp (revision 12292) @@ -67,4 +67,3 @@ function((void*)&handle); #endif - }