Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 12199) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 12200) @@ -609,4 +609,16 @@ ./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,6 +832,4 @@ ./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,4 +956,5 @@ libISSMModules_a_SOURCES = $(module_sources) libISSMModules_a_SOURCES += $(bamg_sources) +libISSMModules_a_SOURCES += $(kriging_sources) libISSMModules_a_SOURCES += $(kml_sources) libISSMModules_a_CXXFLAGS = $(ALLCXXFLAGS) Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp (revision 12199) +++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp (revision 12200) @@ -11,4 +11,5 @@ #include +#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,12 +18,14 @@ /*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,8 +35,19 @@ 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;iSemiVariogram(x[i]-x[j],y[i]-y[j]); Gamma[j*n_obs+i] = Gamma[i*n_obs+j]; } @@ -47,5 +61,5 @@ /*Get semivariogram vector associated to this location*/ - for(i=0;iSemiVariogram(x[i]-x_interp[idx],y[i]-y_interp[idx]); /*Solve the three linear systems*/ @@ -71,35 +85,10 @@ /*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"); - } } Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h (revision 12199) +++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h (revision 12200) @@ -10,5 +10,4 @@ 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); Index: /issm/trunk-jpl/src/c/objects/Constraints/Constraint.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Constraints/Constraint.h (revision 12199) +++ /issm/trunk-jpl/src/c/objects/Constraints/Constraint.h (revision 12200) @@ -11,8 +11,6 @@ /*Headers:*/ /*{{{1*/ +class Nodes; #include "../Object.h" - -class Nodes; - #include "../../toolkits/toolkits.h" /*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Kriging/ExponentialVariogram.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/ExponentialVariogram.cpp (revision 12200) +++ /issm/trunk-jpl/src/c/objects/Kriging/ExponentialVariogram.cpp (revision 12200) @@ -0,0 +1,74 @@ +/*!\file ExponentialVariogram.c + * \brief: implementation of the ExponentialVariogram object + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include +#include +#include "../objects.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../shared/shared.h" +#include "../../include/include.h" + +/*ExponentialVariogram constructors and destructor*/ +/*FUNCTION ExponentialVariogram::ExponentialVariogram(){{{1*/ +ExponentialVariogram::ExponentialVariogram(){ + this->nugget = 0.2; + this->sill = 1; + this->range = SQRT3; + return; +} +/*}}}*/ +/*FUNCTION ExponentialVariogram::ExponentialVariogram(Options* options){{{1*/ +ExponentialVariogram::ExponentialVariogram(Options* options){ + + /*Defaults*/ + this->nugget = 0.2; + this->sill = 1; + this->range = SQRT3; + + /*Overwrite from options*/ + if(options->GetOption("nugget")) options->Get(&this->nugget,"nugget"); + if(options->GetOption("sill")) options->Get(&this->sill,"sill"); + if(options->GetOption("range")) options->Get(&this->range,"range"); + + /*Checks*/ + if(nugget==sill) _error_("nugget and sill cannot be equal (constant semivariogram not allowed)"); +} +/*}}}*/ +/*FUNCTION ExponentialVariogram::~ExponentialVariogram(){{{1*/ +ExponentialVariogram::~ExponentialVariogram(){ + return; +} +/*}}}*/ + +/*Object virtual functions definitions:*/ +/*FUNCTION ExponentialVariogram::Echo {{{1*/ +void ExponentialVariogram::Echo(void){ + printf("ExponentialVariogram\n"); + printf(" nugget: %g\n",this->nugget); + printf(" sill : %g\n",this->sill); + printf(" range : %g\n",this->range); +} +/*}}}*/ + +/*Variogram function*/ +/*FUNCTION ExponentialVariogram::SemiVariogram{{{1*/ +double ExponentialVariogram::SemiVariogram(double deltax,double deltay){ + /*http://en.wikipedia.org/wiki/Variogram*/ + double h,a,gamma; + + /*Calculate length*/ + h=sqrt(pow(deltax,2.)+pow(deltay,2.)); + + /*return semi-variogram*/ + a = 1./3.; + gamma = (sill-nugget)*(1-exp(-h/(a*range))) + nugget; + return gamma; +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Kriging/ExponentialVariogram.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/ExponentialVariogram.h (revision 12200) +++ /issm/trunk-jpl/src/c/objects/Kriging/ExponentialVariogram.h (revision 12200) @@ -0,0 +1,34 @@ +/*! \file ExponentialVariogram.h + * \brief: header file for triavertexinput object + */ + +#ifndef _EXPONENTIALVARIOGRAM_H_ +#define _EXPONENTIALVARIOGRAM_H_ + +/*Headers:*/ +#include "./Variogram.h" + +class ExponentialVariogram: public Variogram{ + + public: + double nugget; //The height of the jump of the semivariogram at the discontinuity at the origin + double sill; //Limit of the variogram tending to infinity lag distances + double range; //The distance in which the difference of the variogram from the sill becomes negligible + + /*ExponentialVariogram constructors, destructors*/ + ExponentialVariogram(); + ExponentialVariogram(Options* options); + ~ExponentialVariogram(); + + /*Object virtual functions definitions*/ + void Echo(); + void DeepEcho(){_error_("Not implemented yet");}; + int Id(){_error_("Not implemented yet");}; + int MyRank(){_error_("Not implemented yet");}; + int ObjectEnum(){_error_("Not implemented yet");}; + Object* copy(){_error_("Not implemented yet");}; + + /*Variogram functions*/ + double SemiVariogram(double deltax,double deltay); +}; +#endif /* _EXPONENTIALVARIOGRAM_H */ Index: /issm/trunk-jpl/src/c/objects/Kriging/GaussianVariogram.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/GaussianVariogram.cpp (revision 12200) +++ /issm/trunk-jpl/src/c/objects/Kriging/GaussianVariogram.cpp (revision 12200) @@ -0,0 +1,74 @@ +/*!\file GaussianVariogram.c + * \brief: implementation of the GaussianVariogram object + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include +#include +#include "../objects.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../shared/shared.h" +#include "../../include/include.h" + +/*GaussianVariogram constructors and destructor*/ +/*FUNCTION GaussianVariogram::GaussianVariogram(){{{1*/ +GaussianVariogram::GaussianVariogram(){ + this->nugget = 0.2; + this->sill = 1; + this->range = SQRT3; + return; +} +/*}}}*/ +/*FUNCTION GaussianVariogram::GaussianVariogram(Options* options){{{1*/ +GaussianVariogram::GaussianVariogram(Options* options){ + + /*Defaults*/ + this->nugget = 0.2; + this->sill = 1; + this->range = SQRT3; + + /*Overwrite from options*/ + if(options->GetOption("nugget")) options->Get(&this->nugget,"nugget"); + if(options->GetOption("sill")) options->Get(&this->sill,"sill"); + if(options->GetOption("range")) options->Get(&this->range,"range"); + + /*Checks*/ + if(nugget==sill) _error_("nugget and sill cannot be equal (constant semivariogram not allowed)"); +} +/*}}}*/ +/*FUNCTION GaussianVariogram::~GaussianVariogram(){{{1*/ +GaussianVariogram::~GaussianVariogram(){ + return; +} +/*}}}*/ + +/*Object virtual functions definitions:*/ +/*FUNCTION GaussianVariogram::Echo {{{1*/ +void GaussianVariogram::Echo(void){ + printf("GaussianVariogram\n"); + printf(" nugget: %g\n",this->nugget); + printf(" sill : %g\n",this->sill); + printf(" range : %g\n",this->range); +} +/*}}}*/ + +/*Variogram function*/ +/*FUNCTION GaussianVariogram::SemiVariogram{{{1*/ +double GaussianVariogram::SemiVariogram(double deltax,double deltay){ + /*http://en.wikipedia.org/wiki/Variogram*/ + double h2,a,gamma; + + /*Calculate length square*/ + h2=pow(deltax,2.)+pow(deltay,2.); + + /*return semi-variogram*/ + a = 1./3.; + gamma = (sill-nugget)*(1-exp(-h2/(a*pow(range,2.)))) + nugget; + return gamma; +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Kriging/GaussianVariogram.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/GaussianVariogram.h (revision 12200) +++ /issm/trunk-jpl/src/c/objects/Kriging/GaussianVariogram.h (revision 12200) @@ -0,0 +1,34 @@ +/*! \file GaussianVariogram.h + * \brief: header file for triavertexinput object + */ + +#ifndef _GAUSSIANVARIOGRAM_H_ +#define _GAUSSIANVARIOGRAM_H_ + +/*Headers:*/ +#include "./Variogram.h" + +class GaussianVariogram: public Variogram{ + + public: + double nugget; //The height of the jump of the semivariogram at the discontinuity at the origin + double sill; //Limit of the variogram tending to infinity lag distances + double range; //The distance in which the difference of the variogram from the sill becomes negligible + + /*GaussianVariogram constructors, destructors*/ + GaussianVariogram(); + GaussianVariogram(Options* options); + ~GaussianVariogram(); + + /*Object virtual functions definitions*/ + void Echo(); + void DeepEcho(){_error_("Not implemented yet");}; + int Id(){_error_("Not implemented yet");}; + int MyRank(){_error_("Not implemented yet");}; + int ObjectEnum(){_error_("Not implemented yet");}; + Object* copy(){_error_("Not implemented yet");}; + + /*Variogram functions*/ + double SemiVariogram(double deltax,double deltay); +}; +#endif /* _GAUSSIANVARIOGRAM_H */ Index: /issm/trunk-jpl/src/c/objects/Kriging/SphericalVariogram.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/SphericalVariogram.cpp (revision 12200) +++ /issm/trunk-jpl/src/c/objects/Kriging/SphericalVariogram.cpp (revision 12200) @@ -0,0 +1,77 @@ +/*!\file SphericalVariogram.c + * \brief: implementation of the SphericalVariogram object + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include +#include +#include "../objects.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../shared/shared.h" +#include "../../include/include.h" + +/*SphericalVariogram constructors and destructor*/ +/*FUNCTION SphericalVariogram::SphericalVariogram(){{{1*/ +SphericalVariogram::SphericalVariogram(){ + this->nugget = 0.2; + this->sill = 1; + this->range = SQRT3; + return; +} +/*}}}*/ +/*FUNCTION SphericalVariogram::SphericalVariogram(Options* options){{{1*/ +SphericalVariogram::SphericalVariogram(Options* options){ + + /*Defaults*/ + this->nugget = 0.2; + this->sill = 1; + this->range = SQRT3; + + /*Overwrite from options*/ + if(options->GetOption("nugget")) options->Get(&this->nugget,"nugget"); + if(options->GetOption("sill")) options->Get(&this->sill,"sill"); + if(options->GetOption("range")) options->Get(&this->range,"range"); + + /*Checks*/ + if(nugget==sill) _error_("nugget and sill cannot be equal (constant semivariogram not allowed)"); +} +/*}}}*/ +/*FUNCTION SphericalVariogram::~SphericalVariogram(){{{1*/ +SphericalVariogram::~SphericalVariogram(){ + return; +} +/*}}}*/ + +/*Object virtual functions definitions:*/ +/*FUNCTION SphericalVariogram::Echo {{{1*/ +void SphericalVariogram::Echo(void){ + printf("SphericalVariogram\n"); + printf(" nugget: %g\n",this->nugget); + printf(" sill : %g\n",this->sill); + printf(" range : %g\n",this->range); +} +/*}}}*/ + +/*Variogram function*/ +/*FUNCTION SphericalVariogram::SemiVariogram{{{1*/ +double SphericalVariogram::SemiVariogram(double deltax,double deltay){ + /*http://en.wikipedia.org/wiki/Variogram*/ + double h,gamma; + + /*Calculate length square*/ + h=sqrt(pow(deltax,2.)+pow(deltay,2.)); + + /*return semi-variogram*/ + if(h<=range) + gamma = (sill-nugget)*( (3*h)/(2*range) - pow(h,3.)/(2*pow(range,3.)) ) + nugget; + else + gamma = sill; + + return gamma; +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/objects/Kriging/SphericalVariogram.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/SphericalVariogram.h (revision 12200) +++ /issm/trunk-jpl/src/c/objects/Kriging/SphericalVariogram.h (revision 12200) @@ -0,0 +1,34 @@ +/*! \file SphericalVariogram.h + * \brief: header file for triavertexinput object + */ + +#ifndef _SPHERICALVARIOGRAM_H_ +#define _SPHERICALVARIOGRAM_H_ + +/*Headers:*/ +#include "./Variogram.h" + +class SphericalVariogram: public Variogram{ + + public: + double nugget; //The height of the jump of the semivariogram at the discontinuity at the origin + double sill; //Limit of the variogram tending to infinity lag distances + double range; //The distance in which the difference of the variogram from the sill becomes negligible + + /*SphericalVariogram constructors, destructors*/ + SphericalVariogram(); + SphericalVariogram(Options* options); + ~SphericalVariogram(); + + /*Object virtual functions definitions*/ + void Echo(); + void DeepEcho(){_error_("Not implemented yet");}; + int Id(){_error_("Not implemented yet");}; + int MyRank(){_error_("Not implemented yet");}; + int ObjectEnum(){_error_("Not implemented yet");}; + Object* copy(){_error_("Not implemented yet");}; + + /*Variogram functions*/ + double SemiVariogram(double deltax,double deltay); +}; +#endif /* _SPHERICALVARIOGRAM_H */ Index: /issm/trunk-jpl/src/c/objects/Kriging/Variogram.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Kriging/Variogram.h (revision 12200) +++ /issm/trunk-jpl/src/c/objects/Kriging/Variogram.h (revision 12200) @@ -0,0 +1,18 @@ +/*!\file: Variogram.h + * \brief abstract class for Variogram object + */ + + +#ifndef _VARIOGRAM_H_ +#define _VARIOGRAM_H_ + +#include "../Object.h" + +class Variogram: public Object{ + + public: + virtual ~Variogram(){}; + virtual double SemiVariogram(double deltax,double deltay)=0; + +}; +#endif Index: /issm/trunk-jpl/src/c/objects/objects.h =================================================================== --- /issm/trunk-jpl/src/c/objects/objects.h (revision 12199) +++ /issm/trunk-jpl/src/c/objects/objects.h (revision 12200) @@ -173,3 +173,9 @@ #include "./Bamg/SetOfE4.h" +/*Kriging*/ +#include "./Kriging/Variogram.h" +#include "./Kriging/GaussianVariogram.h" +#include "./Kriging/ExponentialVariogram.h" +#include "./Kriging/SphericalVariogram.h" + #endif