Changeset 12289

Timestamp:

05/25/12 15:23:23 (13 years ago)

Author:

Mathieu Morlighem

Message:

Added Covariance for Kriging calculation (similar to gslib)

Location:

issm/trunk-jpl/src/c

Files:

• Container/Observations.cpp (modified) (2 diffs)
• modules/Krigingx/Krigingx.cpp (modified) (3 diffs)
• objects/Kriging/ExponentialVariogram.cpp (modified) (1 diff)
• objects/Kriging/ExponentialVariogram.h (modified) (1 diff)
• objects/Kriging/GaussianVariogram.cpp (modified) (2 diffs)
• objects/Kriging/GaussianVariogram.h (modified) (1 diff)
• objects/Kriging/PowerVariogram.cpp (modified) (2 diffs)
• objects/Kriging/PowerVariogram.h (modified) (1 diff)
• objects/Kriging/SphericalVariogram.cpp (modified) (1 diff)
• objects/Kriging/SphericalVariogram.h (modified) (1 diff)
• objects/Kriging/Variogram.h (modified) (1 diff)
• shared/Sorting/Quicksort.c (deleted)
• shared/Sorting/binary_search.cpp (modified) (2 diffs)
• shared/Sorting/sorting.h (modified) (1 diff)

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

@@ -162 +162 @@
 
         /*Output and Intermediaries*/
-        int          nobs,i,j,k;
-        double       range;
+        int          i,j,k;
+        double       distance;
         Observation *observation1 = NULL;
         Observation *observation2 = NULL;
-        int         *indices      = NULL;
 
         int *counter= (int*)xmalloc(n*sizeof(int));
@@ -175 +174 @@
                 observation1=(Observation*)this->GetObjectByOffset(i);
 
-                for(j=0;j<n;j++){
-                        range=x[j]; _assert_(range>=0.);
+                for(j=i+1;j<this->Size();j++){
+                        observation2=(Observation*)this->GetObjectByOffset(j);
 
-                        /*Find all observations that are in range*/
-                        this->quadtree->RangeSearch(&indices,&nobs,observation1->x,observation1->y,range);
+                        distance=sqrt(pow(observation1->x - observation2->x,2.) + pow(observation1->y - observation2->y,2.));
+                        if(distance>x[n-1]) continue;
 
-                        for (k=0;k<nobs;k++){
-                                observation2 = (Observation*)this->GetObjectByOffset(indices[k]);
-                                gamma[j]    += 1./2.*pow(observation1->value - observation2->value,2.);
-                        }
+                        int index = int(distance/(x[1]-x[0]));
+                        if(index>n-1) index = n-1;
+                        if(index<0)   index = 0;
 
-                        counter[j] += nobs;
-                        xfree((void**)&indices);
+                        gamma[index]   += 1./2.*pow(observation1->value - observation2->value,2.);
+                        counter[index] += 1;
                 }
         }
 
         /*Normalize semivariogram*/
-        for(j=0;j<n;j++){
-                if(counter[j]) gamma[j] = gamma[j]/double(counter[j]);
+        gamma[0]=0;
+        for(k=0;k<n;k++){
+                if(counter[k]) gamma[k] = gamma[k]/double(counter[k]);
         }
 

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

@@ -138 +138 @@
                 /*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.; 
+                        continue;
+                }
 
                 /*Allocate intermediary matrix and vectors*/
+
                 Gamma  = (double*)xmalloc(n_obs*n_obs*sizeof(double));
                 gamma0 = (double*)xmalloc(n_obs*sizeof(double));
@@ -147 +153 @@
                 for(i=0;i<n_obs;i++){
                         for(j=0;j<=i;j++){
-                                Gamma[i*n_obs+j] = variogram->SemiVariogram(x[i]-x[j],y[i]-y[j]);
+                                //Gamma[i*n_obs+j] = variogram->SemiVariogram(x[i]-x[j],y[i]-y[j]);
+                                Gamma[i*n_obs+j] = variogram->Covariance(x[i]-x[j],y[i]-y[j]);
                                 Gamma[j*n_obs+i] = Gamma[i*n_obs+j];
                         }
@@ -154 +161 @@
 
                 /*Get semivariogram vector associated to this location*/
-                for(i=0;i<n_obs;i++) gamma0[i] = variogram->SemiVariogram(x[i]-x_interp[idx],y[i]-y_interp[idx]);
+                //for(i=0;i<n_obs;i++) gamma0[i] = variogram->SemiVariogram(x[i]-x_interp[idx],y[i]-y_interp[idx]);
+                for(i=0;i<n_obs;i++) gamma0[i] = variogram->Covariance(x[i]-x_interp[idx],y[i]-y_interp[idx]);
 
                 /*Solve the three linear systems*/

issm/trunk-jpl/src/c/objects/Kriging/ExponentialVariogram.cpp ¶

@@ -59 +59 @@
 
 /*Variogram function*/
+/*FUNCTION ExponentialVariogram::Covariance{{{1*/
+double ExponentialVariogram::Covariance(double deltax,double deltay){
+        /*The covariance can be deduced from the variogram from the following
+         * relationship:
+         *    2 gamma = C(x,x) + C(y,y) -2 C(x,y)
+         * so
+         *    C(h) = sill - gamma                                            */
+        double h,a,cova;
+
+        /*Calculate length*/
+        h=sqrt(pow(deltax,2.)+pow(deltay,2.));
+
+        /*return covariance*/
+        a     = 1./3.;
+        cova = (sill-nugget)*exp(-h/(a*range));
+        return cova;
+}
+/*}}}*/
 /*FUNCTION ExponentialVariogram::SemiVariogram{{{1*/
 double ExponentialVariogram::SemiVariogram(double deltax,double deltay){

issm/trunk-jpl/src/c/objects/Kriging/ExponentialVariogram.h ¶

@@ -31 +31 @@
                 /*Variogram functions*/
                 double SemiVariogram(double deltax,double deltay);
+                double Covariance(double deltax,double deltay);
 };
 #endif  /* _EXPONENTIALVARIOGRAM_H */

issm/trunk-jpl/src/c/objects/Kriging/GaussianVariogram.cpp ¶

@@ -59 +59 @@
 
 /*Variogram function*/
+/*FUNCTION GaussianVariogram::Covariance{{{1*/
+double GaussianVariogram::Covariance(double deltax,double deltay){
+        /*The covariance can be deduced from the variogram from the following
+         * relationship:
+         *    2 gamma = C(x,x) + C(y,y) -2 C(x,y)
+         * so
+         *    C(h) = sill - gamma                                            */
+        double h2,a,cova;
+
+        /*Calculate length square*/
+        h2=pow(deltax,2.)+pow(deltay,2.);
+
+        /*return covariance*/
+        a     = 1./3.;
+        cova = (sill-nugget)*exp(-h2/(a*range*range));
+
+        //if(h2<0.0000001) cova = sill;
+
+        return cova;
+}
+/*}}}*/
 /*FUNCTION GaussianVariogram::SemiVariogram{{{1*/
 double GaussianVariogram::SemiVariogram(double deltax,double deltay){
@@ -70 +91 @@
         a     = 1./3.;
         gamma = (sill-nugget)*(1.-exp(-h2/(a*range*range))) + nugget;
+
+        //if(h2>1000*1000) printf("gamma = %g h= %g\n",gamma,sqrt(h2));
+        printf("h = %g gamma = %g\n",sqrt(h2),gamma);
         return gamma;
 }

issm/trunk-jpl/src/c/objects/Kriging/GaussianVariogram.h ¶

@@ -31 +31 @@
                 /*Variogram functions*/
                 double SemiVariogram(double deltax,double deltay);
+                double Covariance(double deltax,double deltay);
 };
 #endif  /* _GAUSSIANVARIOGRAM_H */

issm/trunk-jpl/src/c/objects/Kriging/PowerVariogram.cpp ¶

@@ -60 +60 @@
 
 /*Variogram function*/
+/*FUNCTION PowerVariogram::Covariance{{{1*/
+double PowerVariogram::Covariance(double deltax,double deltay){
+        /*The covariance can be deduced from the variogram from the following
+         * relationship:
+         *    2 gamma = C(x,x) + C(y,y) -2 C(x,y)
+         * so
+         *    C(h) = sill - gamma                                            */
+        double h,cova;
+
+        /*Calculate length square*/
+        h=sqrt(pow(deltax,2.)+pow(deltay,2.));
+
+        /*return covariance*/
+        cova = 9999. - this->slope*pow(h,this->power);
+
+        return cova;
+}
+/*}}}*/
 /*FUNCTION PowerVariogram::SemiVariogram{{{1*/
 double PowerVariogram::SemiVariogram(double deltax,double deltay){
@@ -71 +89 @@
         gamma = this->nugget + this->slope*pow(h,this->power);
 
+        //if(h>1000) printf("gamma = %g h=%g\n",gamma,h);
         return gamma;
 }

issm/trunk-jpl/src/c/objects/Kriging/PowerVariogram.h ¶

@@ -31 +31 @@
                 /*Variogram functions*/
                 double SemiVariogram(double deltax,double deltay);
+                double Covariance(double deltax,double deltay);
 };
 #endif  /* _POWERVARIOGRAM_H */

issm/trunk-jpl/src/c/objects/Kriging/SphericalVariogram.cpp ¶

@@ -59 +59 @@
 
 /*Variogram function*/
+/*FUNCTION SphericalVariogram::Covariance{{{1*/
+double SphericalVariogram::Covariance(double deltax,double deltay){
+        /*The covariance can be deduced from the variogram from the following
+         * relationship:
+         *    2 gamma = C(x,x) + C(y,y) -2 C(x,y)
+         * so
+         *    C(h) = sill - gamma                                            */
+        double h,cova;
+
+        /*Calculate length square*/
+        h=sqrt(pow(deltax,2.)+pow(deltay,2.));
+
+        /*return covariance*/
+        if(h<=range)
+         cova = (sill-nugget)*(1 - (3*h)/(2*range) + pow(h,3.)/(2*pow(range,3.)) );
+        else
+         cova = 0.;
+
+        return cova;
+}
+/*}}}*/
 /*FUNCTION SphericalVariogram::SemiVariogram{{{1*/
 double SphericalVariogram::SemiVariogram(double deltax,double deltay){

issm/trunk-jpl/src/c/objects/Kriging/SphericalVariogram.h ¶

@@ -31 +31 @@
                 /*Variogram functions*/
                 double SemiVariogram(double deltax,double deltay);
+                double Covariance(double deltax,double deltay);
 };
 #endif  /* _SPHERICALVARIOGRAM_H */

issm/trunk-jpl/src/c/objects/Kriging/Variogram.h ¶

@@ -13 +13 @@
                 virtual ~Variogram(){};
                 virtual double SemiVariogram(double deltax,double deltay)=0;
+                virtual double Covariance(double deltax,double deltay)=0;
 
 };

issm/trunk-jpl/src/c/shared/Sorting/binary_search.cpp ¶

@@ -11 +11 @@
 #include <stdio.h>
 
-int binary_search(int* poffset,int target, int* sorted_integers,int num_integers){
+int binary_search(int* poffset,int target,int* sorted_integers,int num_integers){
 
         /*output: */
         int offset;  //offset, if found
-        int found=0;  //found=0 if target is not found, 1 otherwise.
+        int found=0; //found=0 if target is not found, 1 otherwise.
 
         /*intermediary: */
@@ -64 +64 @@
         return found;
 }
-

issm/trunk-jpl/src/c/shared/Sorting/sorting.h ¶

@@ -6 +6 @@
 #define  _SORTING_H_
 
-int binary_search(int* poffset,int target, int* sorted_integers,int num_integers);
-
+int binary_search(int* poffset,int target,int* sorted_integers,int num_integers);
 
 #endif //ifndef _SORTING_H_
-