Changeset 12200 for issm/trunk-jpl/src/c/modules

Timestamp:

05/03/12 16:14:45 (13 years ago)

Author:

Mathieu Morlighem

Message:

Added Kriging object for semivariograms

Location:

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

Files:

• Krigingx.cpp (modified) (5 diffs)
• Krigingx.h (modified) (1 diff)

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

@@ -11 +11 @@
 #include <gsl/gsl_linalg.h>
 
+#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){
 
@@ -17 +18 @@
 
         /*Intermediaries*/
-        int     i,j;
-        double  numerator,denominator,ratio;
-        double *Gamma     = NULL;
-        double *GinvG0    = NULL;
-        double *Ginv1     = NULL;
-        double *GinvZ     = NULL;
-        double *gamma0    = NULL;
-        double *ones      = NULL;
+        int        i,j;
+        double     numerator,denominator,ratio;
+        char      *model     = NULL;
+        double    *Gamma     = NULL;
+        double    *GinvG0    = NULL;
+        double    *Ginv1     = NULL;
+        double    *GinvZ     = NULL;
+        double    *gamma0    = NULL;
+        double    *ones      = NULL;
+        Variogram *variogram = NULL;
 
         /*Memory allocation*/
@@ -32 +35 @@
         ones        =(double*)xmalloc(n_obs*sizeof(double));
 
+        /*Create Semi-Variogram object*/
+        if(options->GetOption("model")){
+                options->Get(&model,"model");
+                if     (strcmp(model,"gaussian")==0)    variogram = new GaussianVariogram(options);
+                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 variogram = new GaussianVariogram(options);
+        variogram->Echo();
+
         /*First: Create semivariogram matrix for observations*/
         for(i=0;i<n_obs;i++){
                 for(j=0;j<=i;j++){
-                        SemiVariogram(&Gamma[i*n_obs+j],x[i],y[i],x[j],y[j]);
+                        Gamma[i*n_obs+j] = variogram->SemiVariogram(x[i]-x[j],y[i]-y[j]);
                         Gamma[j*n_obs+i] = Gamma[i*n_obs+j];
                 }
@@ -47 +61 @@
 
                 /*Get semivariogram vector associated to this location*/
-                for(i=0;i<n_obs;i++) SemiVariogram(&gamma0[i],x[i],y[i],x_interp[idx],y_interp[idx]);
+                for(i=0;i<n_obs;i++) gamma0[i] = variogram->SemiVariogram(x[i]-x_interp[idx],y[i]-y_interp[idx]);
 
                 /*Solve the three linear systems*/
@@ -71 +85 @@
 
         /*clean-up and Assign output pointer*/
+        delete variogram;
         xfree((void**)&Gamma);
         xfree((void**)&gamma0);
         xfree((void**)&ones);
+        xfree((void**)&model);
         *ppredictions=predictions;
-}
-
-void SemiVariogram(double* gamma,double x1,double y1,double x2,double y2){
-
-        /*Intermediaries*/
-        double a,c0,c1;
-
-        /*Calculate distance*/
-        double r=sqrt(pow(x1-x2,2.)+pow(y1-y2,2.));
-
-        /*Switch between variogram models*/
-        switch(2){
-                case 1:/*Exponential*/
-                        c0=0.2;
-                        c1=0.8;
-                        a =1;
-                        *gamma = c0 + c1*(1-exp(-r/a));
-                        return;
-                case 2:/*Gaussian*/
-                        c0=0.2;
-                        c1=0.8;
-                        a =1;
-                        *gamma = c0 + c1*(1-exp(-pow(r,2.)/pow(a,2.)));
-                        return;
-                default:
-                        _error_("Not implemented yet");
-        }
 }
 

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

@@ -10 +10 @@
 
 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 SemiVariogram(double* gamma,double x1,double y1,double x2,double y2);
 void GslSolve(double** pX,double* A,double* B,int n);