Changeset 12207


Ignore:
Timestamp:
05/04/12 11:27:33 (13 years ago)
Author:
Mathieu Morlighem
Message:

Added from objects for Kriging

Location:
issm/trunk-jpl/src/c
Files:
4 added
5 edited

Legend:

Unmodified
Added
Removed

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

    r12200 r12207  
    617617                                         ./objects/Kriging/SphericalVariogram.h\
    618618                                         ./objects/Kriging/SphericalVariogram.cpp\
     619                                         ./objects/Kriging/Quadtree.h\
     620                                         ./objects/Kriging/Quadtree.cpp\
     621                                         ./objects/Kriging/Observation.h\
     622                                         ./objects/Kriging/Observation.cpp\
    619623                                         ./modules/Krigingx/Krigingx.cpp\
    620624                                         ./modules/Krigingx/Krigingx.h

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

    r12204 r12207  
    1515
    1616#include "../../objects/Kriging/GaussianVariogram.h"
    17 int Krigingx(double** ppredictions,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp,Options* options){
     17int Krigingx(double** ppredictions,double* obs_x, double* obs_y, double* obs_list, int obs_length,double* x_interp,double* y_interp,int n_interp,Options* options){
    1818
    1919        /*output*/
    … …  
    2121
    2222        /*Intermediaries*/
    23         int        i,j;
     23        int        i,j,n_obs;
    2424        double     numerator,denominator,ratio;
     25        double    *x            = NULL;
     26        double    *y            = NULL;
     27        double    *obs          = NULL;
     28        double    *Gamma        = NULL;
     29        double    *GinvG0       = NULL;
     30        double    *Ginv1        = NULL;
     31        double    *GinvZ        = NULL;
     32        double    *gamma0       = NULL;
     33        double    *ones         = NULL;
     34        Variogram *variogram    = NULL;
     35        DataSet   *observations = NULL;
     36
     37        /*Get Variogram from Options*/
     38        ProcessVariogram(&variogram,options);
     39
     40        /*Process observation dataset*/
     41        ProcessObservations(&observations,obs_list,obs_x,obs_y,obs_length);
     42
     43        /*Allocation output*/
     44        predictions =(double*)xmalloc(n_interp*sizeof(double));
     45
     46        /*Loop over all interpolations*/
     47        printf("      interpolation progress:  %5.2lf %%",0.0);
     48        for(int idx=0;idx<n_interp;idx++){
     49                if(idx%10==0) printf("\b\b\b\b\b\b\b%5.2lf %%",(double)idx/n_interp*100);
     50
     51                /*Get list of observations for current point*/
     52                ObservationList(&x,&y,&obs,&n_obs,observations,x_interp[idx],y_interp[idx]);
     53
     54                /*Allocate intermediary matrix and vectors*/
     55                Gamma       =(double*)xmalloc(n_obs*n_obs*sizeof(double));
     56                gamma0      =(double*)xmalloc(n_obs*sizeof(double));
     57                ones        =(double*)xmalloc(n_obs*sizeof(double));
     58
     59                /*First: Create semivariogram matrix for observations*/
     60                for(i=0;i<n_obs;i++){
     61                        for(j=0;j<=i;j++){
     62                                Gamma[i*n_obs+j] = variogram->SemiVariogram(x[i]-x[j],y[i]-y[j]);
     63                                Gamma[j*n_obs+i] = Gamma[i*n_obs+j];
     64                        }
     65                }
     66                for(i=0;i<n_obs;i++) ones[i]=1;
     67
     68                /*Get semivariogram vector associated to this location*/
     69                for(i=0;i<n_obs;i++) gamma0[i] = variogram->SemiVariogram(x[i]-x_interp[idx],y[i]-y_interp[idx]);
     70
     71                /*Solve the three linear systems*/
     72                GslSolve(&GinvG0,Gamma,gamma0,n_obs);       // Gamma^-1 gamma0
     73                GslSolve(&Ginv1, Gamma,ones,n_obs);         // Gamma^-1 ones
     74                GslSolve(&GinvZ, Gamma,obs,n_obs); // Gamma^-1 Z
     75
     76                /*Prepare predictor*/
     77                numerator=-1.; denominator=0.;
     78                for(i=0;i<n_obs;i++) numerator  +=GinvG0[i];
     79                for(i=0;i<n_obs;i++) denominator+=Ginv1[i];
     80                ratio=numerator/denominator;
     81
     82                predictions[idx]=0;
     83                for(i=0;i<n_obs;i++) predictions[idx] += (gamma0[i]-ratio)*GinvZ[i];
     84
     85                /*clean-up*/
     86                xfree((void**)&x);
     87                xfree((void**)&y);
     88                xfree((void**)&obs);
     89                xfree((void**)&Gamma);
     90                xfree((void**)&gamma0);
     91                xfree((void**)&ones);
     92                xfree((void**)&GinvG0);
     93                xfree((void**)&Ginv1);
     94                xfree((void**)&GinvZ);
     95        }
     96        printf("\b\b\b\b\b\b\b\b%5.2lf %%\n",100.0);
     97
     98        /*clean-up and Assign output pointer*/
     99        delete variogram;
     100        delete observations;
     101        *ppredictions=predictions;
     102}
     103
     104void ProcessVariogram(Variogram **pvariogram,Options* options){/*{{{*/
     105
     106        /*Intermediaries*/
     107        Variogram* variogram = NULL;
    25108        char      *model     = NULL;
    26         double    *Gamma     = NULL;
    27         double    *GinvG0    = NULL;
    28         double    *Ginv1     = NULL;
    29         double    *GinvZ     = NULL;
    30         double    *gamma0    = NULL;
    31         double    *ones      = NULL;
    32         Variogram *variogram = NULL;
    33 
    34         /*Memory allocation*/
    35         predictions =(double*)xmalloc(n_interp*sizeof(double));
    36         Gamma       =(double*)xmalloc(n_obs*n_obs*sizeof(double));
    37         gamma0      =(double*)xmalloc(n_obs*sizeof(double));
    38         ones        =(double*)xmalloc(n_obs*sizeof(double));
    39 
    40         /*Create Semi-Variogram object*/
     109
    41110        if(options->GetOption("model")){
    42111                options->Get(&model,"model");
    … …  
    44113                else if(strcmp(model,"exponential")==0) variogram = new ExponentialVariogram(options);
    45114                else if(strcmp(model,"spherical")==0)   variogram = new SphericalVariogram(options);
    46                 else _error_("only gaussian semivariogram models supported yet");
     115                else _error_("variogram %s not supported yet (list of supported variogram: gaussian, exponential and spherical)",model);
    47116        }
    48117        else variogram = new GaussianVariogram(options);
    49         variogram->Echo();
    50 
    51         /*First: Create semivariogram matrix for observations*/
    52         for(i=0;i<n_obs;i++){
    53                 for(j=0;j<=i;j++){
    54                         Gamma[i*n_obs+j] = variogram->SemiVariogram(x[i]-x[j],y[i]-y[j]);
    55                         Gamma[j*n_obs+i] = Gamma[i*n_obs+j];
    56                 }
    57         }
    58         for(i=0;i<n_obs;i++) ones[i]=1;
    59 
    60         /*Loop over all interpolations*/
    61         printf("      interpolation progress:  %5.2lf %%",0.0);
    62         for(int idx=0;idx<n_interp;idx++){
    63                 if(idx%100==0) printf("\b\b\b\b\b\b\b%5.2lf %%",(double)idx/n_interp*100);
    64 
    65                 /*Get semivariogram vector associated to this location*/
    66                 for(i=0;i<n_obs;i++) gamma0[i] = variogram->SemiVariogram(x[i]-x_interp[idx],y[i]-y_interp[idx]);
    67 
    68                 /*Solve the three linear systems*/
    69                 GslSolve(&GinvG0,Gamma,gamma0,n_obs);       // Gamma^-1 gamma0
    70                 GslSolve(&Ginv1, Gamma,ones,n_obs);         // Gamma^-1 ones
    71                 GslSolve(&GinvZ, Gamma,observations,n_obs); // Gamma^-1 Z
    72 
    73                 /*Prepare predictor*/
    74                 numerator=-1.; denominator=0.;
    75                 for(i=0;i<n_obs;i++) numerator  +=GinvG0[i];
    76                 for(i=0;i<n_obs;i++) denominator+=Ginv1[i];
    77                 ratio=numerator/denominator;
    78 
    79                 predictions[idx]=0;
    80                 for(i=0;i<n_obs;i++) predictions[idx] += (gamma0[i]-ratio)*GinvZ[i];
    81 
    82                 /*clean-up*/
    83                 xfree((void**)&GinvG0);
    84                 xfree((void**)&Ginv1);
    85                 xfree((void**)&GinvZ);
    86         }
    87         printf("\b\b\b\b\b\b\b\b%5.2lf %%\n",100.0);
    88 
    89         /*clean-up and Assign output pointer*/
    90         delete variogram;
    91         xfree((void**)&Gamma);
    92         xfree((void**)&gamma0);
    93         xfree((void**)&ones);
     118
     119        /*Assign output pointer*/
    94120        xfree((void**)&model);
    95         *ppredictions=predictions;
    96 }
    97 
    98 void GslSolve(double** pX,double* A,double* B,int n){
     121        *pvariogram = variogram;
     122}/*}}}*/
     123void ProcessObservations(DataSet **pobservations,double* observations_list,double* x,double* y,int n){/*{{{*/
     124
     125        int i;
     126        DataSet* observations = NULL;
     127
     128        /*Initialize Observation Dataset*/
     129        observations = new DataSet();
     130
     131        /*Add observations one by one*/
     132        for(i=0;i<n;i++){
     133                observations->AddObject(new Observation(x[i],y[i],0,0,observations_list[i]));
     134        }
     135
     136        /*Assign output pointer*/
     137        *pobservations = observations;
     138}/*}}}*/
     139void ObservationList(double **px,double **py,double **pobs,int* pnobs,DataSet* observations,double x_interp,double y_interp){/*{{{*/
     140
     141        /*Output and Intermediaries*/
     142        int          nobs,i;
     143        double      *x           = NULL;
     144        double      *y           = NULL;
     145        double      *obs         = NULL;
     146        Observation *observation = NULL;
     147
     148        /*Get number of observations*/
     149        nobs = observations->Size();
     150
     151        /*Allocate vectors*/
     152        x   = (double*)xmalloc(nobs*sizeof(double));
     153        y   = (double*)xmalloc(nobs*sizeof(double));
     154        obs = (double*)xmalloc(nobs*sizeof(double));
     155
     156        /*Loop over all observations and fill in x, y and obs*/
     157        for (i=0;i<observations->Size();i++){
     158                observation=(Observation*)observations->GetObjectByOffset(i);
     159                observation->WriteXYObs(&x[i],&y[i],&obs[i]);
     160        }
     161
     162        /*Assign output pointer*/
     163        *px=x;
     164        *py=y;
     165        *pobs=obs;
     166        *pnobs=nobs;
     167}/*}}}*/
     168void GslSolve(double** pX,double* A,double* B,int n){/*{{{*/
    99169#ifdef _HAVE_GSL_
    100170
    101         /*GSL Matrices and vectors: */
    102         int              s;
    103         gsl_matrix_view  a;
    104         gsl_vector_view  b;
    105         gsl_vector      *x = NULL;
    106         gsl_permutation *p = NULL;
    107 
    108         /*A will be modified by LU decomposition. Use copy*/
    109         double* Acopy = (double*)xmalloc(n*n*sizeof(double));
    110         memcpy(Acopy,A,n*n*sizeof(double));
    111 
    112         /*Initialize gsl matrices and vectors: */
    113         a = gsl_matrix_view_array (Acopy,n,n);
    114         b = gsl_vector_view_array (B,n);
    115         x = gsl_vector_alloc (n);
    116 
    117         /*Run LU and solve: */
    118         p = gsl_permutation_alloc (n);
    119         gsl_linalg_LU_decomp (&a.matrix, p, &s);
    120         gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x);
    121 
    122         //printf ("x = \n");
    123         //gsl_vector_fprintf (stdout, x, "%g");
    124 
    125         /*Copy result*/
    126         double* X = (double*)xmalloc(n*sizeof(double));
    127         memcpy(X,gsl_vector_ptr(x,0),n*sizeof(double));
    128 
    129         /*Clean up and assign output pointer*/
    130         xfree((void**)&Acopy);
    131         gsl_permutation_free(p);
    132         gsl_vector_free(x);
    133         *pX=X;
     171                /*GSL Matrices and vectors: */
     172                int              s;
     173                gsl_matrix_view  a;
     174                gsl_vector_view  b;
     175                gsl_vector      *x = NULL;
     176                gsl_permutation *p = NULL;
     177
     178                /*A will be modified by LU decomposition. Use copy*/
     179                double* Acopy = (double*)xmalloc(n*n*sizeof(double));
     180                memcpy(Acopy,A,n*n*sizeof(double));
     181
     182                /*Initialize gsl matrices and vectors: */
     183                a = gsl_matrix_view_array (Acopy,n,n);
     184                b = gsl_vector_view_array (B,n);
     185                x = gsl_vector_alloc (n);
     186
     187                /*Run LU and solve: */
     188                p = gsl_permutation_alloc (n);
     189                gsl_linalg_LU_decomp (&a.matrix, p, &s);
     190                gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x);
     191
     192                //printf ("x = \n");
     193                //gsl_vector_fprintf (stdout, x, "%g");
     194
     195                /*Copy result*/
     196                double* X = (double*)xmalloc(n*sizeof(double));
     197                memcpy(X,gsl_vector_ptr(x,0),n*sizeof(double));
     198
     199                /*Clean up and assign output pointer*/
     200                xfree((void**)&Acopy);
     201                gsl_permutation_free(p);
     202                gsl_vector_free(x);
     203                *pX=X;
    134204#else
    135         _error_("GSL support required");
     205                _error_("GSL support required");
    136206#endif
    137 }
     207        }/*}}}*/

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

    r12200 r12207  
    1010
    1111int  Krigingx(double** ppredictions,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp,Options* options);
     12void ProcessVariogram(Variogram **pvariogram,Options* options);
     13void ProcessObservations(DataSet **pobservations,double* observations_list,double* x,double* y,int n);
     14void ObservationList(double **px,double **py,double **pobs,int* pnobs,DataSet* observations,double x_interp,double y_interp);
    1215void GslSolve(double** pX,double* A,double* B,int n);
    1316

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

    r12200 r12207  
    22 * \brief abstract class for Variogram object
    33 */
    4 
    54
    65#ifndef _VARIOGRAM_H_

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

    r12200 r12207  
    178178#include "./Kriging/ExponentialVariogram.h"
    179179#include "./Kriging/SphericalVariogram.h"
     180#include "./Kriging/Quadtree.h"
     181#include "./Kriging/Observation.h"
    180182
    181183#endif
Note: See TracChangeset for help on using the changeset viewer.