Changeset 12422

Timestamp:

06/15/12 12:13:38 (13 years ago)

Author:

Mathieu Morlighem

Message:

Simpler kriging

Location:

issm/trunk-jpl/src/c

Files:

• Container/Observations.cpp (modified) (2 diffs)
• Container/Observations.h (modified) (1 diff)
• modules/Krigingx/Krigingx.cpp (modified) (4 diffs)
• modules/Krigingx/pKrigingx.cpp (modified) (5 diffs)

issm/trunk-jpl/src/c/Container/Observations.cpp

@@ -276 +276 @@
         *pnobs=nobs;
 }/*}}}*/
+/*FUNCTION Observations::InterpolationIDW{{{*/
+void Observations::InterpolationIDW(double *pprediction,double x_interp,double y_interp,double radius,int mindata,int maxdata,double power){
+
+        /*Intermediaries*/
+        int    i,n_obs;
+        double prediction;
+        double numerator,denominator,h,weight;
+        double *x   = NULL;
+        double *y   = NULL;
+        double *obs = NULL;
+
+        /*Some checks*/
+        _assert_(maxdata>0);
+        _assert_(pprediction);
+        _assert_(power>0);
+
+        /*If radius is not provided or is 0, return all observations*/
+        if(radius==0) radius=this->quadtree->root->length;
+
+        /*Get list of observations for current point*/
+        this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata);
+
+        /*If we have less observations than mindata, return UNDEF*/
+        if(n_obs<mindata){
+                prediction = UNDEF; 
+        }
+        else{
+                numerator   = 0.;
+                denominator = 0.;
+                for(i=0;i<n_obs;i++){
+                        h = sqrt( (x[i]-x_interp)*(x[i]-x_interp) + (y[i]-y_interp)*(y[i]-y_interp));
+                        if (h<0.0000001){
+                                numerator   = obs[i];
+                                denominator = 1.;
+                                break;
+                        }
+                        weight = 1./pow(h,power);
+                        numerator   += weight*obs[i];
+                        denominator += weight;
+                }
+                prediction = numerator/denominator; 
+        }
+
+        /*clean-up*/
+        *pprediction = prediction;
+        xfree((void**)&x);
+        xfree((void**)&y);
+        xfree((void**)&obs);
+}/*}}}*/
 /*FUNCTION Observations::InterpolationKriging{{{*/
 void Observations::InterpolationKriging(double *pprediction,double *perror,double x_interp,double y_interp,double radius,int mindata,int maxdata,Variogram* variogram){
@@ -363 +412 @@
 
 }/*}}}*/
+/*FUNCTION Observations::InterpolationNearestNeighbor{{{*/
+void Observations::InterpolationNearestNeighbor(double *pprediction,double x_interp,double y_interp,double radius){
+
+        /*Intermediaries*/
+        double        x,y,obs;
+
+        /*Get clostest observation*/
+        this->ClosestObservation(&x,&y,&obs,x_interp,y_interp,radius);
+
+        /*Assign output pointer*/
+        *pprediction = obs;
+}/*}}}*/
 /*FUNCTION Observations::QuadtreeColoring{{{*/
 void Observations::QuadtreeColoring(double* A,double *x,double *y,int n){

issm/trunk-jpl/src/c/Container/Observations.h

@@ -25 +25 @@
                 /*Methods*/
                 void ClosestObservation(double *px,double *py,double *pobs,double x_interp,double y_interp,double radius);
+                void InterpolationIDW(double *pprediction,double x_interp,double y_interp,double radius,int mindata,int maxdata,double power);
                 void InterpolationKriging(double *pprediction,double *perror,double x_interp,double y_interp,double radius,int mindata,int maxdata,Variogram* variogram);
+                void InterpolationNearestNeighbor(double *pprediction,double x_interp,double y_interp,double radius);
                 void ObservationList(double **px,double **py,double **pobs,int* pnobs);
                 void ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double radius,int maxdata);

issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp

@@ -228 +228 @@
 
         /*Intermediaries*/
-        int           i,j,n_obs;
-        double        x,y,obs;
+        int           i;
         double        localpercent;
 
@@ -242 +241 @@
                 printf("\r      interpolation progress: %5.2lf%%",localpercent);
 
-                /*Get closest observation to current point*/
-                observations->ClosestObservation(&x,&y,&obs,x_interp[idx],y_interp[idx],radius);
-                predictions[idx] = obs; 
+                observations->InterpolationNearestNeighbor(&predictions[idx],x_interp[idx],y_interp[idx],radius);
         }
 
@@ -279 +276 @@
 
         /*Intermediaries*/
-        int           i,j,n_obs;
-        double        localpercent;
-        double        numerator,denominator,h,weight;
-        double *x     = NULL;
-        double *y     = NULL;
-        double *obs   = NULL;
+        double localpercent;
         double  power = 2.;
 
@@ -294 +286 @@
                 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]);
+                for(int i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]);
                 printf("\r      interpolation progress: %5.2lf%%",localpercent);
 
-                /*Get closest observation to current point*/
-                observations->ObservationList(&x,&y,&obs,&n_obs,x_interp[idx],y_interp[idx],radius,maxdata);
-                if(n_obs){
-                        numerator   = 0;
-                        denominator = 0;
-                        for(j=0;j<n_obs;j++){
-                                h      = sqrt( (x[j]-x_interp[idx])*(x[j]-x_interp[idx]) + (y[j]-y_interp[idx])*(y[j]-y_interp[idx]));
-                                weight = 1./pow(h,power);
-                                numerator   += weight*obs[j];
-                                denominator += weight;
-                        }
-                        predictions[idx] = numerator/denominator; 
-                }
-                else{
-                        predictions[idx] = UNDEF; 
-                }
-        }
-
+                observations->InterpolationIDW(&predictions[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,power);
+        }
         return NULL;
 }/*}}}*/

issm/trunk-jpl/src/c/modules/Krigingx/pKrigingx.cpp

@@ -9 +9 @@
 #include "../../objects/objects.h"
 #include "../../Container/Container.h"
-#include "../modules.h"
-
-#ifdef _HAVE_GSL_
-#include <gsl/gsl_linalg.h>
-#endif
+#include "../../io/io.h"
 
 /*FUNCTION pKrigingx{{{*/
@@ -32 +28 @@
         Observations *observations = NULL;
 
-        /*Get Variogram from Options*/
-        ProcessVariogram(&variogram,options);
+        /*Get some Options*/
         options->Get(&radius,"searchradius",0.);
         options->Get(&mindata,"mindata",1);
@@ -56 +51 @@
         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){
@@ -72 +70 @@
 #endif
         }
+        else if(strcmp(output,"nearestneighbor")==0){
+
+                /*partition loop across threads: */
+                for(int idx=my_rank;idx<n_interp;idx+=num_procs){
+                        _printf_(true,"\r      interpolation progress: %5.2lf%%",double(idx)/double(n_interp)*100);
+                        observations->InterpolationNearestNeighbor(&predictions[idx],x_interp[idx],y_interp[idx],radius);
+                }
+                _printf_(true,"\r      interpolation progress: %5.2lf%%\n",100.);
+
+#ifdef _HAVE_MPI_
+                double *sumpredictions =(double*)xmalloc(n_interp*sizeof(double));
+                MPI_Allreduce(predictions,sumpredictions,n_interp,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+                xfree((void**)&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){
+                        _printf_(true,"\r      interpolation progress: %5.2lf%%",double(idx)/double(n_interp)*100);
+                        observations->InterpolationIDW(&predictions[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,power);
+                }
+                _printf_(true,"\r      interpolation progress: %5.2lf%%\n",100.);
+
+#ifdef _HAVE_MPI_
+                double *sumpredictions =(double*)xmalloc(n_interp*sizeof(double));
+                MPI_Allreduce(predictions,sumpredictions,n_interp,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
+                xfree((void**)&predictions); predictions=sumpredictions;
+#endif
+        }
         else{
                 _error_("output '%s' not supported yet",output);
         }
-        _printf_(true,"\r      interpolation progress: %5.2lf%%\n",100.);
 
         /*clean-up and Assign output pointer*/
@@ -105 +134 @@
         *pvariogram = variogram;
 }/*}}}*/
-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;
-#else
-                _error_("GSL support required");
-#endif
-        }/*}}}*/