Changeset 12230

Timestamp:

05/10/12 14:19:39 (13 years ago)

Author:

Mathieu Morlighem

Message:

Added power variogram

Location:

issm/trunk-jpl/src/c

Files:

• Container/Observations.cpp (modified) (1 diff)
• Container/Observations.h (modified) (1 diff)
• Makefile.am (modified) (1 diff)
• modules/Krigingx/Krigingx.cpp (modified) (2 diffs)
• objects/Kriging/PowerVariogram.cpp (added)
• objects/Kriging/PowerVariogram.h (added)
• objects/Kriging/Quadtree.cpp (modified) (1 diff)
• objects/objects.h (modified) (1 diff)

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

@@ -150 +150 @@
 
 }/*}}}*/
+/*FUNCTION Observations::Variomap{{{*/
+void Observations::Variomap(double* gamma,double *x,int n){
+
+        /*Output and Intermediaries*/
+        int          nobs,i,j,k;
+        double       range;
+        Observation *observation1 = NULL;
+        Observation *observation2 = NULL;
+        int         *indices      = NULL;
+
+        int *counter= (int*)xmalloc(n*sizeof(int));
+        for(j=0;j<n;j++) counter[j] = 0;
+        for(j=0;j<n;j++) gamma[j]   = 0.0;
+
+        for(i=0;i<this->Size();i++){
+                observation1=(Observation*)this->GetObjectByOffset(i);
+
+                for(j=0;j<n;j++){
+                        range=x[j]; _assert_(range>=0.);
+
+                        /*Find all observations that are in range*/
+                        this->quadtree->RangeSearch(&indices,&nobs,observation1->x,observation1->y,range);
+
+                        for (k=0;k<nobs;k++){
+                                observation2 = (Observation*)this->GetObjectByOffset(indices[k]);
+                                gamma[j]    += 1./2.*pow(observation1->value - observation2->value,2.);
+                        }
+
+                        counter[j] += nobs;
+                        xfree((void**)&indices);
+                }
+        }
+
+        /*Normalize semivariogram*/
+        for(j=0;j<n;j++){
+                if(counter[j]) gamma[j] = gamma[j]/double(counter[j]);
+        }
+
+        /*Assign output pointer*/
+        xfree((void**)&counter);
+}/*}}}*/

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

@@ -25 +25 @@
                 void ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double range);
                 void QuadtreeColoring(double* A,double *x,double *y,int n);
+                void Variomap(double* gamma,double *x,int n);
 
 };

issm/trunk-jpl/src/c/Makefile.am

@@ -619 +619 @@
                                                 ./objects/Kriging/SphericalVariogram.h\
                                                 ./objects/Kriging/SphericalVariogram.cpp\
+                                                ./objects/Kriging/PowerVariogram.h\
+                                                ./objects/Kriging/PowerVariogram.cpp\
                                                 ./objects/Kriging/Quadtree.h\
                                                 ./objects/Kriging/Quadtree.cpp\

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

@@ -54 +54 @@
         if(strcmp(output,"quadtree")==0){
                 observations->QuadtreeColoring(predictions,x_interp,y_interp,n_interp);
+                *ppredictions=predictions;
+                return 1;
+        }
+        if(strcmp(output,"variomap")==0){
+                observations->Variomap(predictions,x_interp,n_interp);
                 *ppredictions=predictions;
                 return 1;
@@ -129 +134 @@
                 else if(strcmp(model,"exponential")==0) variogram = new ExponentialVariogram(options);
                 else if(strcmp(model,"spherical")==0)   variogram = new SphericalVariogram(options);
-                else _error_("variogram %s not supported yet (list of supported variogram: gaussian, exponential and spherical)",model);
+                else if(strcmp(model,"power")==0)       variogram = new PowerVariogram(options);
+                else _error_("variogram %s not supported yet (list of supported variogram: gaussian, exponential, spherical and power)",model);
         }
         else variogram = new GaussianVariogram(options);

issm/trunk-jpl/src/c/objects/Kriging/Quadtree.cpp

@@ -433 +433 @@
 
         /*Return 2 if the this box is included in the range*/
-        if(
-                                this->xcenter+this->length/2 <= x+range &&
-                                this->ycenter+this->length/2 <= y+range &&
-                                this->xcenter-this->length/2 >= x-range &&
-                                this->ycenter-this->length/2 >= y-range
-          ) return 2;
+        if(this->xcenter+this->length/2 <= x+range &&
+                this->ycenter+this->length/2 <= y+range &&
+                this->xcenter-this->length/2 >= x-range &&
+                this->ycenter-this->length/2 >= y-range) return 2;
 
         /*This is a simple overlap*/

issm/trunk-jpl/src/c/objects/objects.h

@@ -178 +178 @@
 #include "./Kriging/ExponentialVariogram.h"
 #include "./Kriging/SphericalVariogram.h"
+#include "./Kriging/PowerVariogram.h"
 #include "./Kriging/Quadtree.h"
 #include "./Kriging/Observation.h"