Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 12206) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 12207) @@ -617,4 +617,8 @@ ./objects/Kriging/SphericalVariogram.h\ ./objects/Kriging/SphericalVariogram.cpp\ + ./objects/Kriging/Quadtree.h\ + ./objects/Kriging/Quadtree.cpp\ + ./objects/Kriging/Observation.h\ + ./objects/Kriging/Observation.cpp\ ./modules/Krigingx/Krigingx.cpp\ ./modules/Krigingx/Krigingx.h Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp (revision 12206) +++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp (revision 12207) @@ -15,5 +15,5 @@ #include "../../objects/Kriging/GaussianVariogram.h" -int Krigingx(double** ppredictions,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp,Options* options){ +int Krigingx(double** ppredictions,double* obs_x, double* obs_y, double* obs_list, int obs_length,double* x_interp,double* y_interp,int n_interp,Options* options){ /*output*/ @@ -21,22 +21,91 @@ /*Intermediaries*/ - int i,j; + int i,j,n_obs; double numerator,denominator,ratio; + double *x = NULL; + double *y = NULL; + double *obs = NULL; + double *Gamma = NULL; + double *GinvG0 = NULL; + double *Ginv1 = NULL; + double *GinvZ = NULL; + double *gamma0 = NULL; + double *ones = NULL; + Variogram *variogram = NULL; + DataSet *observations = NULL; + + /*Get Variogram from Options*/ + ProcessVariogram(&variogram,options); + + /*Process observation dataset*/ + ProcessObservations(&observations,obs_list,obs_x,obs_y,obs_length); + + /*Allocation output*/ + predictions =(double*)xmalloc(n_interp*sizeof(double)); + + /*Loop over all interpolations*/ + printf(" interpolation progress: %5.2lf %%",0.0); + for(int idx=0;idxSemiVariogram(x[i]-x[j],y[i]-y[j]); + Gamma[j*n_obs+i] = Gamma[i*n_obs+j]; + } + } + for(i=0;iSemiVariogram(x[i]-x_interp[idx],y[i]-y_interp[idx]); + + /*Solve the three linear systems*/ + GslSolve(&GinvG0,Gamma,gamma0,n_obs); // Gamma^-1 gamma0 + GslSolve(&Ginv1, Gamma,ones,n_obs); // Gamma^-1 ones + GslSolve(&GinvZ, Gamma,obs,n_obs); // Gamma^-1 Z + + /*Prepare predictor*/ + numerator=-1.; denominator=0.; + for(i=0;iGetOption("model")){ options->Get(&model,"model"); @@ -44,94 +113,95 @@ else if(strcmp(model,"exponential")==0) variogram = new ExponentialVariogram(options); else if(strcmp(model,"spherical")==0) variogram = new SphericalVariogram(options); - else _error_("only gaussian semivariogram models supported yet"); + else _error_("variogram %s not supported yet (list of supported variogram: gaussian, exponential and spherical)",model); } else variogram = new GaussianVariogram(options); - variogram->Echo(); - - /*First: Create semivariogram matrix for observations*/ - for(i=0;iSemiVariogram(x[i]-x[j],y[i]-y[j]); - Gamma[j*n_obs+i] = Gamma[i*n_obs+j]; - } - } - for(i=0;iSemiVariogram(x[i]-x_interp[idx],y[i]-y_interp[idx]); - - /*Solve the three linear systems*/ - GslSolve(&GinvG0,Gamma,gamma0,n_obs); // Gamma^-1 gamma0 - GslSolve(&Ginv1, Gamma,ones,n_obs); // Gamma^-1 ones - GslSolve(&GinvZ, Gamma,observations,n_obs); // Gamma^-1 Z - - /*Prepare predictor*/ - numerator=-1.; denominator=0.; - for(i=0;iAddObject(new Observation(x[i],y[i],0,0,observations_list[i])); + } + + /*Assign output pointer*/ + *pobservations = observations; +}/*}}}*/ +void ObservationList(double **px,double **py,double **pobs,int* pnobs,DataSet* observations,double x_interp,double y_interp){/*{{{*/ + + /*Output and Intermediaries*/ + int nobs,i; + double *x = NULL; + double *y = NULL; + double *obs = NULL; + Observation *observation = NULL; + + /*Get number of observations*/ + nobs = observations->Size(); + + /*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;iSize();i++){ + observation=(Observation*)observations->GetObjectByOffset(i); + observation->WriteXYObs(&x[i],&y[i],&obs[i]); + } + + /*Assign output pointer*/ + *px=x; + *py=y; + *pobs=obs; + *pnobs=nobs; +}/*}}}*/ +void GslSolve(double** pX,double* A,double* B,int n){/*{{{*/ #ifdef _HAVE_GSL_ - /*GSL Matrices and vectors: */ - int s; - gsl_matrix_view a; - gsl_vector_view b; - gsl_vector *x = NULL; - gsl_permutation *p = NULL; - - /*A will be modified by LU decomposition. Use copy*/ - double* Acopy = (double*)xmalloc(n*n*sizeof(double)); - memcpy(Acopy,A,n*n*sizeof(double)); - - /*Initialize gsl matrices and vectors: */ - a = gsl_matrix_view_array (Acopy,n,n); - b = gsl_vector_view_array (B,n); - x = gsl_vector_alloc (n); - - /*Run LU and solve: */ - p = gsl_permutation_alloc (n); - gsl_linalg_LU_decomp (&a.matrix, p, &s); - gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x); - - //printf ("x = \n"); - //gsl_vector_fprintf (stdout, x, "%g"); - - /*Copy result*/ - double* X = (double*)xmalloc(n*sizeof(double)); - memcpy(X,gsl_vector_ptr(x,0),n*sizeof(double)); - - /*Clean up and assign output pointer*/ - xfree((void**)&Acopy); - gsl_permutation_free(p); - gsl_vector_free(x); - *pX=X; + /*GSL Matrices and vectors: */ + int s; + gsl_matrix_view a; + gsl_vector_view b; + gsl_vector *x = NULL; + gsl_permutation *p = NULL; + + /*A will be modified by LU decomposition. Use copy*/ + double* Acopy = (double*)xmalloc(n*n*sizeof(double)); + memcpy(Acopy,A,n*n*sizeof(double)); + + /*Initialize gsl matrices and vectors: */ + a = gsl_matrix_view_array (Acopy,n,n); + b = gsl_vector_view_array (B,n); + x = gsl_vector_alloc (n); + + /*Run LU and solve: */ + p = gsl_permutation_alloc (n); + gsl_linalg_LU_decomp (&a.matrix, p, &s); + gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x); + + //printf ("x = \n"); + //gsl_vector_fprintf (stdout, x, "%g"); + + /*Copy result*/ + double* X = (double*)xmalloc(n*sizeof(double)); + memcpy(X,gsl_vector_ptr(x,0),n*sizeof(double)); + + /*Clean up and assign output pointer*/ + xfree((void**)&Acopy); + gsl_permutation_free(p); + gsl_vector_free(x); + *pX=X; #else - _error_("GSL support required"); + _error_("GSL support required"); #endif -} + }/*}}}*/ Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h (revision 12206) +++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h (revision 12207) @@ -10,4 +10,7 @@ int Krigingx(double** ppredictions,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp,Options* options); +void ProcessVariogram(Variogram **pvariogram,Options* options); +void ProcessObservations(DataSet **pobservations,double* observations_list,double* x,double* y,int n); +void ObservationList(double **px,double **py,double **pobs,int* pnobs,DataSet* observations,double x_interp,double y_interp); void GslSolve(double** pX,double* A,double* B,int n); Index: /issm/trunk-jpl/src/c/objects/Kriging/Observation.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/Observation.cpp (revision 12207) +++ /issm/trunk-jpl/src/c/objects/Kriging/Observation.cpp (revision 12207) @@ -0,0 +1,49 @@ +/*!\file Observation.c + * \brief: implementation of the Observation object + */ + +#include "../objects.h" + +/*Observation constructors and destructor*/ +/*FUNCTION Observation::Observation(){{{1*/ +Observation::Observation(){ + return; +} +/*}}}*/ +/*FUNCTION Observation::Observation(double x,double y,int xi,int yi,double value){{{1*/ +Observation::Observation(double x_in,double y_in,int xi_in,int yi_in,double value_in){ + + this->x = x_in; + this->y = y_in; + this->xi = xi_in; + this->yi = yi_in; + this->value = value_in; + +} +/*}}}*/ +/*FUNCTION Observation::~Observation(){{{1*/ +Observation::~Observation(){ + return; +} +/*}}}*/ + +/*Object virtual functions definitions:*/ +/*FUNCTION Observation::Echo {{{1*/ +void Observation::Echo(void){ + printf("Observation\n"); + printf(" x : %g\n",this->x); + printf(" y : %g\n",this->y); + printf(" xi : %i\n",this->xi); + printf(" yi : %i\n",this->yi); + printf(" value : %g\n",this->value); +} +/*}}}*/ + +/*Observations functions*/ +/*FUNCTION Observation::WriteXYObs(double* px,double* py,double* pobs){{{1*/ +void Observation::WriteXYObs(double* px,double* py,double* pobs){ + *px = this->x; + *py = this->y; + *pobs = this->value; +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Kriging/Observation.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/Observation.h (revision 12207) +++ /issm/trunk-jpl/src/c/objects/Kriging/Observation.h (revision 12207) @@ -0,0 +1,33 @@ +/*! \file Observation.h + * \brief: header file for Observation object + */ + +#ifndef _OBSERVATION_H_ +#define _OBSERVATION_H_ + +#include "../Object.h" + +class Observation: public Object{ + + public: + double x,y; + int xi,yi; + double value; + + /*Observation constructors, destructors*/ + Observation(); + Observation(double x_in,double y_in,int xi_in,int yi_in,double value_in); + ~Observation(); + + /*Object virtual functions definitions*/ + void Echo(); + void DeepEcho() {_error_("Not implemented yet"); }; + int Id() {_error_("Not implemented yet"); }; + int MyRank() {_error_("Not implemented yet"); }; + int ObjectEnum(){_error_("Not implemented yet"); }; + Object *copy() {_error_("Not implemented yet"); }; + + /*Management*/ + void WriteXYObs(double* px,double* py,double* pobs); +}; +#endif /* _EXPONENTIALVARIOGRAM_H */ Index: /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp (revision 12207) +++ /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp (revision 12207) @@ -0,0 +1,98 @@ +#include "../objects.h" +/*DOCUMENTATION What is a Quadtree? {{{1 + * A Quadtree is a very simple way to group vertices according + * to their locations. A square that holds all the points of the mesh + * (or the geometry) is divided into 4 boxes. As soon as one box + * hold more than 4 vertices, it is divided into 4 new boxes, etc... + * There cannot be more than MAXDEEP (=30) subdivision. + * This process is like a Dichotomy in dimension 2 + * + * + - - - - - - - - + - - + - + - + - - - - + + * | | | | X | | + * + - + - + + * | | | | | | + * + - - + - + - + + + * | | | | | + * + * | | | | | + * + - - - - - - - - + - - + - - + - - - - + + * | | | | + * + * | | | | + * + * | | | | + * | | | | + * + - - - - - - - - + - - - - + - - - - + + * | | | + * + * | | | + * + * | | | + * + * | | | + * | | | + * | | | + * | | | + * | | | + * + - - - - - - - - + - - - - - - - - + + * + * The coordinate system used in a quadtree are integers to avoid + * round-off errors. The vertex in the lower left box has the coordinates + * (0 0) + * The upper right vertex has the follwing coordinates: + * 2^30 -1 2^30 -1 in decimal + * 0 1 1 1 .... 1 0 1 1 1 .... 1 in binary + * \-- 29 --/ \-- 29 --/ + * Using binaries is therefore very easy to locate a vertex in a box: + * we just need to look at the bits from the left to the right (See ::Add) + }}}1*/ + + /*Constructors/Destructors*/ + /*FUNCTION Quadtree::Quadtree(){{{1*/ + Quadtree::Quadtree(){ + + /*Number of boxes and vertices*/ + NbQuadtreeBox=0; + NbPoints=0; + + /*Create container*/ + boxcontainer=new DataSet(); + + /*Create Root, pointer toward the main box*/ + root=NewQuadtreeBox(); + + } + /*}}}1*/ + /*FUNCTION Quadtree::~Quadtree(){{{1*/ + Quadtree::~Quadtree() { + + delete boxcontainer; + root=NULL; + + } + /*}}}1*/ + + /*Methods*/ + /*FUNCTION Quadtree::NewQuadtreeBox {{{1*/ + Quadtree::QuadtreeBox* Quadtree::NewQuadtreeBox(void){ + + /*Output*/ + QuadtreeBox* newbox=NULL; + + /*Create and initialize a new box*/ + newbox=new QuadtreeBox(); + newbox->nbitems=0; + newbox->box[0]=NULL; + newbox->box[1]=NULL; + newbox->box[2]=NULL; + newbox->box[3]=NULL; + + /*Add root to the container*/ + boxcontainer->AddObject(newbox); + + /*Increase counter*/ + NbQuadtreeBox++; + + /*currentbox now points toward next quadtree box*/ + return newbox; + }/*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.h (revision 12207) +++ /issm/trunk-jpl/src/c/objects/Kriging/Quadtree.h (revision 12207) @@ -0,0 +1,43 @@ + +#ifndef _QUADTREEK_H +#define _QUADTREEK_H + +class Quadtree{ + + private: + /* A quadtree box contains up to 4 points (locations). 4 new quadtree boxes are + * created if a fifth point is added to the same box. A Quadtree box is therefore + * composed of EITHER: + * - up to 4 vertices + * - 4 "sub" quadtree boxes*/ + + class QuadtreeBox: public Object{ + public: + int nbitems; // number of current vertices in the box + union{ + QuadtreeBox *box[4]; + int index[4]; + }; + + /*Object functions (Needed because the Quadtree uses a Container*/ + void Echo() {_error_("not implemented yet"); }; + void DeepEcho() {_error_("not implemented yet"); }; + int Id() {_error_("not implemented yet"); }; + int MyRank() {_error_("not implemented yet"); }; + int ObjectEnum(){_error_("not implemented yet"); }; + Object *copy() {_error_("not implemented yet"); }; + }; + + /*Quadtree private Fields*/ + DataSet* boxcontainer; + + public: + QuadtreeBox* root; // main box + long NbQuadtreeBox; // total number of boxes + long NbPoints; // number of points + + Quadtree(); + ~Quadtree(); + QuadtreeBox* NewQuadtreeBox(void); +}; +#endif //_QUADTREEK_H Index: /issm/trunk-jpl/src/c/objects/Kriging/Variogram.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/Variogram.h (revision 12206) +++ /issm/trunk-jpl/src/c/objects/Kriging/Variogram.h (revision 12207) @@ -2,5 +2,4 @@ * \brief abstract class for Variogram object */ - #ifndef _VARIOGRAM_H_ Index: /issm/trunk-jpl/src/c/objects/objects.h =================================================================== --- /issm/trunk-jpl/src/c/objects/objects.h (revision 12206) +++ /issm/trunk-jpl/src/c/objects/objects.h (revision 12207) @@ -178,4 +178,6 @@ #include "./Kriging/ExponentialVariogram.h" #include "./Kriging/SphericalVariogram.h" +#include "./Kriging/Quadtree.h" +#include "./Kriging/Observation.h" #endif