Changeset 12200

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

Files:

• Makefile.am (modified) (3 diffs)
• modules/Krigingx/Krigingx.cpp (modified) (5 diffs)
• modules/Krigingx/Krigingx.h (modified) (1 diff)
• objects/Constraints/Constraint.h (modified) (1 diff)
• objects/Kriging (added)
• objects/Kriging/ExponentialVariogram.cpp (added)
• objects/Kriging/ExponentialVariogram.h (added)
• objects/Kriging/GaussianVariogram.cpp (added)
• objects/Kriging/GaussianVariogram.h (added)
• objects/Kriging/SphericalVariogram.cpp (added)
• objects/Kriging/SphericalVariogram.h (added)
• objects/Kriging/Variogram.h (added)
• objects/objects.h (modified) (1 diff)

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

@@ -609 +609 @@
                                 ./modules/BamgTriangulatex/BamgTriangulatex.h
 #}}}
+#Kriging sources  {{{1
+kriging_sources = ./objects/Kriging/Variogram.h \
+                                         ./objects/Kriging/GaussianVariogram.h\
+                                         ./objects/Kriging/GaussianVariogram.cpp\
+                                         ./objects/Kriging/ExponentialVariogram.h\
+                                         ./objects/Kriging/ExponentialVariogram.cpp\
+                                         ./objects/Kriging/SphericalVariogram.h\
+                                         ./objects/Kriging/SphericalVariogram.cpp\
+                                         ./modules/Krigingx/Krigingx.cpp\
+                                         ./modules/Krigingx/Krigingx.h
+
+#}}}
 #Kml sources  {{{1
 kml_sources = ./modules/Exp2Kmlx/Exp2Kmlx.h\
@@ -820 +832 @@
                         ./modules/HoleFillerx/HoleFillerx.cpp\
                         ./modules/HoleFillerx/HoleFillerx.h\
-                        ./modules/Krigingx/Krigingx.cpp\
-                        ./modules/Krigingx/Krigingx.h\
                         ./modules/AverageFilterx/AverageFilterx.cpp\
                         ./modules/AverageFilterx/AverageFilterx.h\
@@ -946 +956 @@
 libISSMModules_a_SOURCES = $(module_sources)
 libISSMModules_a_SOURCES += $(bamg_sources)
+libISSMModules_a_SOURCES += $(kriging_sources)
 libISSMModules_a_SOURCES += $(kml_sources)
 libISSMModules_a_CXXFLAGS = $(ALLCXXFLAGS)

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);
 

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

@@ -11 +11 @@
 /*Headers:*/
 /*{{{1*/
+class Nodes;
 #include "../Object.h"
-
-class Nodes;
-
 #include "../../toolkits/toolkits.h"
 /*}}}*/

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

@@ -173 +173 @@
 #include "./Bamg/SetOfE4.h"
 
+/*Kriging*/
+#include "./Kriging/Variogram.h"
+#include "./Kriging/GaussianVariogram.h"
+#include "./Kriging/ExponentialVariogram.h"
+#include "./Kriging/SphericalVariogram.h"
+
 #endif