Changeset 12421


Ignore:
Timestamp:
06/15/12 11:48:36 (13 years ago)
Author:
Mathieu Morlighem
Message:

Added other interpolation for Kriging

Location:
issm/trunk-jpl/src/c/modules/Krigingx
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r12377 r12421  
    5656                observations->Variomap(predictions,x_interp,n_interp);
    5757        }
    58         else if(strcmp(output,"prediction")==0){
    59 
     58        else if(strcmp(output,"delaunay")==0){
     59                int nobs,nel;
     60                double *x     = NULL;
     61                double *y     = NULL;
     62                double *data  = NULL;
     63                int    *index = NULL;
     64                double *indexd= NULL;
     65
     66                observations->ObservationList(&x,&y,&data,&nobs);
     67
     68                printf("Generation Delaunay Triangulation\n");
     69                BamgTriangulatex(&index,&nel,x,y,nobs);
     70                indexd =(double*)xcalloc(nel*3,sizeof(double));
     71                for(int i=0;i<nel*3;i++) indexd[i]=(double)index[i];
     72                xfree((void**)&index);
     73
     74                printf("Interpolating\n");
     75                xfree((void**)&predictions);
     76                InterpFromMeshToMesh2dx(&predictions,indexd,x,y,nobs,nel,data,nobs,1,x_interp,y_interp,n_interp,NULL,0,new DataSet());
     77                xfree((void**)&x);
     78                xfree((void**)&y);
     79                xfree((void**)&index);
     80                xfree((void**)&data);
     81        }
     82        else if(strcmp(output,"nearestneighbor")==0){
    6083                /*initialize thread parameters: */
    6184                gate.n_interp     = n_interp;
    … …  
    7295
    7396                /*launch the thread manager with Krigingxt as a core: */
     97                LaunchThread(NearestNeighbort,(void*)&gate,num);
     98                printf("\r      interpolation progress:  100.00%%\n");
     99                xfree((void**)&gate.percent);
     100        }
     101        else if(strcmp(output,"idw")==0){ //Inverse distance weighting
     102                /*initialize thread parameters: */
     103                gate.n_interp     = n_interp;
     104                gate.x_interp     = x_interp;
     105                gate.y_interp     = y_interp;
     106                gate.radius       = radius;
     107                gate.mindata      = mindata;
     108                gate.maxdata      = maxdata;
     109                gate.variogram    = variogram;
     110                gate.observations = observations;
     111                gate.predictions  = predictions;
     112                gate.error        = error;
     113                gate.percent      = (double*)xcalloc(num,sizeof(double));
     114
     115                /*launch the thread manager with Krigingxt as a core: */
     116                LaunchThread(idwt,(void*)&gate,num);
     117                printf("\r      interpolation progress:  100.00%%\n");
     118                xfree((void**)&gate.percent);
     119        }
     120        else if(strcmp(output,"prediction")==0){
     121
     122                /*initialize thread parameters: */
     123                gate.n_interp     = n_interp;
     124                gate.x_interp     = x_interp;
     125                gate.y_interp     = y_interp;
     126                gate.radius       = radius;
     127                gate.mindata      = mindata;
     128                gate.maxdata      = maxdata;
     129                gate.variogram    = variogram;
     130                gate.observations = observations;
     131                gate.predictions  = predictions;
     132                gate.error        = error;
     133                gate.percent      = (double*)xcalloc(num,sizeof(double));
     134
     135                /*launch the thread manager with Krigingxt as a core: */
    74136                LaunchThread(Krigingxt,(void*)&gate,num);
    75137                printf("\r      interpolation progress:  100.00%%\n");
    … …  
    118180
    119181        /*Intermediaries*/
     182        double        localpercent;
     183
     184        /*partition loop across threads: */
     185        PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
     186        for(int idx=i0;idx<i1;idx++){
     187
     188                /*Print info*/
     189                percent[my_thread]=double(idx-i0)/double(i1-i0)*100.;
     190                localpercent=percent[0];
     191                for(int i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]);
     192                printf("\r      interpolation progress: %5.2lf%%",localpercent);
     193
     194                /*Kriging interpolation*/
     195                observations->InterpolationKriging(&predictions[idx],&error[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,variogram);
     196        }
     197
     198        return NULL;
     199}/*}}}*/
     200/*FUNCTION NearestNeighbort{{{*/
     201void* NearestNeighbort(void* vpthread_handle){
     202
     203        /*gate variables :*/
     204        KrigingxThreadStruct *gate        = NULL;
     205        pthread_handle       *handle      = NULL;
     206        int my_thread;
     207        int num_threads;
     208        int i0,i1;
     209
     210        /*recover handle and gate: */
     211        handle      = (pthread_handle*)vpthread_handle;
     212        gate        = (KrigingxThreadStruct*)handle->gate;
     213        my_thread   = handle->id;
     214        num_threads = handle->num;
     215
     216        /*recover parameters :*/
     217        int           n_interp     = gate->n_interp;
     218        double       *x_interp     = gate->x_interp;
     219        double       *y_interp     = gate->y_interp;
     220        double        radius       = gate->radius;
     221        int           mindata      = gate->mindata;
     222        int           maxdata      = gate->maxdata;
     223        Variogram    *variogram    = gate->variogram;
     224        Observations *observations = gate->observations;
     225        double       *predictions  = gate->predictions;
     226        double       *error        = gate->error;
     227        double       *percent      = gate->percent;
     228
     229        /*Intermediaries*/
    120230        int           i,j,n_obs;
    121         double        numerator,denominator,ratio,localpercent;
    122         double       *x            = NULL;
    123         double       *y            = NULL;
    124         double       *obs          = NULL;
    125         double       *Gamma        = NULL;
    126         double       *GinvG0       = NULL;
    127         double       *Ginv1        = NULL;
    128         double       *GinvZ        = NULL;
    129         double       *gamma0       = NULL;
    130         double       *ones         = NULL;
     231        double        x,y,obs;
     232        double        localpercent;
    131233
    132234        /*partition loop across threads: */
    … …  
    140242                printf("\r      interpolation progress: %5.2lf%%",localpercent);
    141243
    142                 /*Get list of observations for current point*/
     244                /*Get closest observation to current point*/
     245                observations->ClosestObservation(&x,&y,&obs,x_interp[idx],y_interp[idx],radius);
     246                predictions[idx] = obs;
     247        }
     248
     249        return NULL;
     250}/*}}}*/
     251/*FUNCTION idwt{{{*/
     252void* idwt(void* vpthread_handle){
     253
     254        /*gate variables :*/
     255        KrigingxThreadStruct *gate        = NULL;
     256        pthread_handle       *handle      = NULL;
     257        int my_thread;
     258        int num_threads;
     259        int i0,i1;
     260
     261        /*recover handle and gate: */
     262        handle      = (pthread_handle*)vpthread_handle;
     263        gate        = (KrigingxThreadStruct*)handle->gate;
     264        my_thread   = handle->id;
     265        num_threads = handle->num;
     266
     267        /*recover parameters :*/
     268        int           n_interp     = gate->n_interp;
     269        double       *x_interp     = gate->x_interp;
     270        double       *y_interp     = gate->y_interp;
     271        double        radius       = gate->radius;
     272        int           mindata      = gate->mindata;
     273        int           maxdata      = gate->maxdata;
     274        Variogram    *variogram    = gate->variogram;
     275        Observations *observations = gate->observations;
     276        double       *predictions  = gate->predictions;
     277        double       *error        = gate->error;
     278        double       *percent      = gate->percent;
     279
     280        /*Intermediaries*/
     281        int           i,j,n_obs;
     282        double        localpercent;
     283        double        numerator,denominator,h,weight;
     284        double *x     = NULL;
     285        double *y     = NULL;
     286        double *obs   = NULL;
     287        double  power = 2.;
     288
     289        /*partition loop across threads: */
     290        PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
     291        for(int idx=i0;idx<i1;idx++){
     292
     293                /*Print info*/
     294                percent[my_thread]=double(idx-i0)/double(i1-i0)*100.;
     295                localpercent=percent[0];
     296                for(i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]);
     297                printf("\r      interpolation progress: %5.2lf%%",localpercent);
     298
     299                /*Get closest observation to current point*/
    143300                observations->ObservationList(&x,&y,&obs,&n_obs,x_interp[idx],y_interp[idx],radius,maxdata);
    144                 if(n_obs<mindata){
    145                         predictions[idx] = -999.0;
    146                         error[idx]       = -999.0;
    147                         continue;
     301                if(n_obs){
     302                        numerator   = 0;
     303                        denominator = 0;
     304                        for(j=0;j<n_obs;j++){
     305                                h      = sqrt( (x[j]-x_interp[idx])*(x[j]-x_interp[idx]) + (y[j]-y_interp[idx])*(y[j]-y_interp[idx]));
     306                                weight = 1./pow(h,power);
     307                                numerator   += weight*obs[j];
     308                                denominator += weight;
     309                        }
     310                        predictions[idx] = numerator/denominator;
    148311                }
    149 
    150                 /*Allocate intermediary matrix and vectors*/
    151 
    152                 Gamma  = (double*)xmalloc(n_obs*n_obs*sizeof(double));
    153                 gamma0 = (double*)xmalloc(n_obs*sizeof(double));
    154                 ones   = (double*)xmalloc(n_obs*sizeof(double));
    155 
    156                 /*First: Create semivariogram matrix for observations*/
    157                 for(i=0;i<n_obs;i++){
    158                         for(j=0;j<=i;j++){
    159                                 //Gamma[i*n_obs+j] = variogram->SemiVariogram(x[i]-x[j],y[i]-y[j]);
    160                                 Gamma[i*n_obs+j] = variogram->Covariance(x[i]-x[j],y[i]-y[j]);
    161                                 Gamma[j*n_obs+i] = Gamma[i*n_obs+j];
    162                         }
     312                else{
     313                        predictions[idx] = UNDEF;
    163314                }
    164                 for(i=0;i<n_obs;i++) ones[i]=1;
    165 
    166                 /*Get semivariogram vector associated to this location*/
    167                 //for(i=0;i<n_obs;i++) gamma0[i] = variogram->SemiVariogram(x[i]-x_interp[idx],y[i]-y_interp[idx]);
    168                 for(i=0;i<n_obs;i++) gamma0[i] = variogram->Covariance(x[i]-x_interp[idx],y[i]-y_interp[idx]);
    169 
    170                 /*Solve the three linear systems*/
    171                 GslSolve(&GinvG0,Gamma,gamma0,n_obs); // Gamma^-1 gamma0
    172                 GslSolve(&Ginv1, Gamma,ones,n_obs);   // Gamma^-1 ones
    173                 GslSolve(&GinvZ, Gamma,obs,n_obs);    // Gamma^-1 Z
    174 
    175                 /*Prepare predictor*/
    176                 numerator=-1.; denominator=0.;
    177                 for(i=0;i<n_obs;i++) numerator  +=GinvG0[i];
    178                 for(i=0;i<n_obs;i++) denominator+=Ginv1[i];
    179                 ratio=numerator/denominator;
    180 
    181                 predictions[idx] = 0.;
    182                 error[idx]       = - numerator*numerator/denominator;
    183                 for(i=0;i<n_obs;i++) predictions[idx] += (gamma0[i]-ratio)*GinvZ[i];
    184                 for(i=0;i<n_obs;i++) error[idx] += gamma0[i]*GinvG0[i];
    185 
    186                 /*clean-up*/
    187                 xfree((void**)&x);
    188                 xfree((void**)&y);
    189                 xfree((void**)&obs);
    190                 xfree((void**)&Gamma);
    191                 xfree((void**)&gamma0);
    192                 xfree((void**)&ones);
    193                 xfree((void**)&GinvG0);
    194                 xfree((void**)&Ginv1);
    195                 xfree((void**)&GinvZ);
    196315        }
    197316
    … …  
    219338        *pvariogram = variogram;
    220339}/*}}}*/
    221 void GslSolve(double** pX,double* A,double* B,int n){/*{{{*/
    222 #ifdef _HAVE_GSL_
    223 
    224                 /*GSL Matrices and vectors: */
    225                 int              s;
    226                 gsl_matrix_view  a;
    227                 gsl_vector_view  b;
    228                 gsl_vector      *x = NULL;
    229                 gsl_permutation *p = NULL;
    230 
    231                 /*A will be modified by LU decomposition. Use copy*/
    232                 double* Acopy = (double*)xmalloc(n*n*sizeof(double));
    233                 memcpy(Acopy,A,n*n*sizeof(double));
    234 
    235                 /*Initialize gsl matrices and vectors: */
    236                 a = gsl_matrix_view_array (Acopy,n,n);
    237                 b = gsl_vector_view_array (B,n);
    238                 x = gsl_vector_alloc (n);
    239 
    240                 /*Run LU and solve: */
    241                 p = gsl_permutation_alloc (n);
    242                 gsl_linalg_LU_decomp (&a.matrix, p, &s);
    243                 gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x);
    244 
    245                 //printf ("x = \n");
    246                 //gsl_vector_fprintf (stdout, x, "%g");
    247 
    248                 /*Copy result*/
    249                 double* X = (double*)xmalloc(n*sizeof(double));
    250                 memcpy(X,gsl_vector_ptr(x,0),n*sizeof(double));
    251 
    252                 /*Clean up and assign output pointer*/
    253                 xfree((void**)&Acopy);
    254                 gsl_permutation_free(p);
    255                 gsl_vector_free(x);
    256                 *pX=X;
    257 #else
    258                 _error_("GSL support required");
    259 #endif
    260         }/*}}}*/

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

    r12377 r12421  
    1515int  pKrigingx(double** ppredictions,double **perror,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp,Options* options);
    1616void ProcessVariogram(Variogram **pvariogram,Options* options);
    17 void GslSolve(double** pX,double* A,double* B,int n);
    1817
    1918/*threading: */
    … …  
    3332
    3433void* Krigingxt(void*);
     34void* NearestNeighbort(void*);
     35void* idwt(void*);
    3536#endif /* _KRIGINGX_H */

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

    r12409 r12421  
    5656        else if(strcmp(output,"prediction")==0){
    5757
    58                 /*Intermediaries*/
    59                 int           i,j,n_obs;
    60                 double        numerator,denominator,ratio,localpercent;
    61                 double       *x            = NULL;
    62                 double       *y            = NULL;
    63                 double       *obs          = NULL;
    64                 double       *Gamma        = NULL;
    65                 double       *GinvG0       = NULL;
    66                 double       *Ginv1        = NULL;
    67                 double       *GinvZ        = NULL;
    68                 double       *gamma0       = NULL;
    69                 double       *ones         = NULL;
    70 
    7158                /*partition loop across threads: */
    7259                for(int idx=my_rank;idx<n_interp;idx+=num_procs){
    73 
    74                         /*Print info*/
    7560                        _printf_(true,"\r      interpolation progress: %5.2lf%%",double(idx)/double(n_interp)*100);
    76 
    77                         /*Get list of observations for current point*/
    78                         observations->ObservationList(&x,&y,&obs,&n_obs,x_interp[idx],y_interp[idx],radius,maxdata);
    79                         if(n_obs<mindata){
    80                                 predictions[idx] = -999.0;
    81                                 error[idx]       = -999.0;
    82                                 continue;
    83                         }
    84 
    85                         /*Allocate intermediary matrix and vectors*/
    86 
    87                         Gamma  = (double*)xmalloc(n_obs*n_obs*sizeof(double));
    88                         gamma0 = (double*)xmalloc(n_obs*sizeof(double));
    89                         ones   = (double*)xmalloc(n_obs*sizeof(double));
    90 
    91                         /*First: Create semivariogram matrix for observations*/
    92                         for(i=0;i<n_obs;i++){
    93                                 for(j=0;j<=i;j++){
    94                                         //Gamma[i*n_obs+j] = variogram->SemiVariogram(x[i]-x[j],y[i]-y[j]);
    95                                         Gamma[i*n_obs+j] = variogram->Covariance(x[i]-x[j],y[i]-y[j]);
    96                                         Gamma[j*n_obs+i] = Gamma[i*n_obs+j];
    97                                 }
    98                         }
    99                         for(i=0;i<n_obs;i++) ones[i]=1;
    100 
    101                         /*Get semivariogram vector associated to this location*/
    102                         //for(i=0;i<n_obs;i++) gamma0[i] = variogram->SemiVariogram(x[i]-x_interp[idx],y[i]-y_interp[idx]);
    103                         for(i=0;i<n_obs;i++) gamma0[i] = variogram->Covariance(x[i]-x_interp[idx],y[i]-y_interp[idx]);
    104 
    105                         /*Solve the three linear systems*/
    106                         GslSolve(&GinvG0,Gamma,gamma0,n_obs); // Gamma^-1 gamma0
    107                         GslSolve(&Ginv1, Gamma,ones,n_obs);   // Gamma^-1 ones
    108                         GslSolve(&GinvZ, Gamma,obs,n_obs);    // Gamma^-1 Z
    109 
    110                         /*Prepare predictor*/
    111                         numerator=-1.; denominator=0.;
    112                         for(i=0;i<n_obs;i++) numerator  +=GinvG0[i];
    113                         for(i=0;i<n_obs;i++) denominator+=Ginv1[i];
    114                         ratio=numerator/denominator;
    115 
    116                         predictions[idx] = 0.;
    117                         error[idx]       = - numerator*numerator/denominator;
    118                         for(i=0;i<n_obs;i++) predictions[idx] += (gamma0[i]-ratio)*GinvZ[i];
    119                         for(i=0;i<n_obs;i++) error[idx] += gamma0[i]*GinvG0[i];
    120 
    121                         /*clean-up*/
    122                         xfree((void**)&x);
    123                         xfree((void**)&y);
    124                         xfree((void**)&obs);
    125                         xfree((void**)&Gamma);
    126                         xfree((void**)&gamma0);
    127                         xfree((void**)&ones);
    128                         xfree((void**)&GinvG0);
    129                         xfree((void**)&Ginv1);
    130                         xfree((void**)&GinvZ);
     61                        observations->InterpolationKriging(&predictions[idx],&error[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,variogram);
    13162                }
     63                _printf_(true,"\r      interpolation progress: %5.2lf%%\n",100.);
    13264
    13365#ifdef _HAVE_MPI_
    … …  
    13971                xfree((void**)&predictions); predictions=sumpredictions;
    14072#endif
    141 
    14273        }
    14374        else{
Note: See TracChangeset for help on using the changeset viewer.