source: issm/oecreview/Archive/12221-12240/ISSM-12230-12231.diff@ 12325

Last change on this file since 12325 was 12325, checked in by Eric.Larour, 13 years ago

11990 to 12321 oec compliance
File size: 8.3 KB

  • proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/Container/Observations.cpp

     
    112112
    113113        if(nobs==0){
    114114                /*No observation found, double range*/
    115                 printf("No observation found within range, doubling range\n");
     115                //printf("No observation found within range, doubling range\n");
    116116                xfree((void**)&indices);
    117117                this->ObservationList(&x,&y,&obs,&nobs,x_interp,y_interp,range*2);
    118118        }
    119119        else{
    120                 if(nobs>1000) printf("Taking more than 1000 observations\n");
     120                //if(nobs>1000) printf("Taking more than 1000 observations\n");
    121121                /*Allocate vectors*/
    122122                x   = (double*)xmalloc(nobs*sizeof(double));
    123123                y   = (double*)xmalloc(nobs*sizeof(double));

  • proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp

     
    1515#endif
    1616
    1717#include "../../objects/Kriging/GaussianVariogram.h"
     18/*FUNCTION Krigingx{{{*/
    1819int 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){
    1920
    2021        /*output*/
    2122        double *predictions = NULL;
    2223
    2324        /*Intermediaries*/
    24         int           i,j,n_obs;
    2525        double        range;
    26         double        numerator,denominator,ratio;
    27         double       *x            = NULL;
    28         double       *y            = NULL;
    29         double       *obs          = NULL;
    30         double       *Gamma        = NULL;
    31         double       *GinvG0       = NULL;
    32         double       *Ginv1        = NULL;
    33         double       *GinvZ        = NULL;
    34         double       *gamma0       = NULL;
    35         double       *ones         = NULL;
    3626        char         *output       = NULL;
    3727        Variogram    *variogram    = NULL;
    3828        Observations *observations = NULL;
    3929
     30        /*threading: */
     31        KrigingxThreadStruct gate;
     32        int num=1;
     33#ifdef _MULTITHREADING_
     34        num=_NUMTHREADS_;
     35#endif
     36
    4037        /*Get Variogram from Options*/
    4138        ProcessVariogram(&variogram,options);
    4239        options->Get(&range,"searchrange",0.);
    … …  
    4441        /*Process observation dataset*/
    4542        observations=new Observations(obs_list,obs_x,obs_y,obs_length,options);
    4643
    47         /*Allocation output*/
     44        /*Allocate output*/
    4845        predictions =(double*)xmalloc(n_interp*sizeof(double));
    49         for(i=0;i<n_interp;i++) predictions[i]=0;
     46        for(int i=0;i<n_interp;i++) predictions[i]=0;
    5047
    5148        /*Get output*/
    52         options->Get(&output,"output","quadtree");
     49        options->Get(&output,"output","prediction");
    5350
    5451        if(strcmp(output,"quadtree")==0){
    5552                observations->QuadtreeColoring(predictions,x_interp,y_interp,n_interp);
    56                 *ppredictions=predictions;
    57                 return 1;
    5853        }
    59         if(strcmp(output,"variomap")==0){
     54        else if(strcmp(output,"variomap")==0){
    6055                observations->Variomap(predictions,x_interp,n_interp);
    61                 *ppredictions=predictions;
    62                 return 1;
    6356        }
     57        else if(strcmp(output,"prediction")==0){
    6458
    65         /*Loop over all interpolations*/
    66         printf("      interpolation progress:  %5.2lf %%",0.0);
    67         for(int idx=0;idx<n_interp;idx++){
    68                 if(idx%10==0) printf("\b\b\b\b\b\b\b%5.2lf %%",(double)idx/n_interp*100);
     59                /*initialize thread parameters: */
     60                gate.n_interp     = n_interp;
     61                gate.x_interp     = x_interp;
     62                gate.y_interp     = y_interp;
     63                gate.range        = range;
     64                gate.variogram    = variogram;
     65                gate.observations = observations;
     66                gate.predictions     = predictions;
    6967
     68                /*launch the thread manager with Krigingxt as a core: */
     69                LaunchThread(Krigingxt,(void*)&gate,num);
     70        }
     71        else{
     72                _error_("output '%s' not supported yet",output);
     73        }
     74
     75        /*clean-up and Assign output pointer*/
     76        delete variogram;
     77        delete observations;
     78        xfree((void**)&output);
     79        *ppredictions=predictions;
     80        return 1;
     81}/*}}}*/
     82/*FUNCTION Krigingxt{{{*/
     83void* Krigingxt(void* vpthread_handle){
     84
     85        /*gate variables :*/
     86        KrigingxThreadStruct *gate        = NULL;
     87        pthread_handle       *handle      = NULL;
     88        int my_thread;
     89        int num_threads;
     90        int i0,i1;
     91
     92        /*recover handle and gate: */
     93        handle      = (pthread_handle*)vpthread_handle;
     94        gate        = (KrigingxThreadStruct*)handle->gate;
     95        my_thread   = handle->id;
     96        num_threads = handle->num;
     97
     98        /*recover parameters :*/
     99        int           n_interp     = gate->n_interp;
     100        double       *x_interp     = gate->x_interp;
     101        double       *y_interp     = gate->y_interp;
     102        double        range        = gate->range;
     103        Variogram    *variogram    = gate->variogram;
     104        Observations *observations = gate->observations;
     105        double       *predictions  = gate->predictions;
     106
     107        /*Intermediaries*/
     108        int           i,j,n_obs;
     109        double        numerator,denominator,ratio;
     110        double       *x            = NULL;
     111        double       *y            = NULL;
     112        double       *obs          = NULL;
     113        double       *Gamma        = NULL;
     114        double       *GinvG0       = NULL;
     115        double       *Ginv1        = NULL;
     116        double       *GinvZ        = NULL;
     117        double       *gamma0       = NULL;
     118        double       *ones         = NULL;
     119
     120        /*partition loop across threads: */
     121        if(my_thread==0) printf("      interpolation progress:  %5.2lf %%",0.0);
     122        PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
     123        for(int idx=i0;idx<i1;idx++){
     124                if(my_thread==0 && idx%10==0) printf("\b\b\b\b\b\b\b%5.2lf %%",double(idx-i0)/double(i1-i0)*100);
     125
    70126                /*Get list of observations for current point*/
    71127                observations->ObservationList(&x,&y,&obs,&n_obs,x_interp[idx],y_interp[idx],range);
    72128
    … …  
    88144                for(i=0;i<n_obs;i++) gamma0[i] = variogram->SemiVariogram(x[i]-x_interp[idx],y[i]-y_interp[idx]);
    89145
    90146                /*Solve the three linear systems*/
    91                 GslSolve(&GinvG0,Gamma,gamma0,n_obs);       // Gamma^-1 gamma0
    92                 GslSolve(&Ginv1, Gamma,ones,n_obs);         // Gamma^-1 ones
    93                 GslSolve(&GinvZ, Gamma,obs,n_obs); // Gamma^-1 Z
     147                GslSolve(&GinvG0,Gamma,gamma0,n_obs); // Gamma^-1 gamma0
     148                GslSolve(&Ginv1, Gamma,ones,n_obs);   // Gamma^-1 ones
     149                GslSolve(&GinvZ, Gamma,obs,n_obs);    // Gamma^-1 Z
    94150
    95151                /*Prepare predictor*/
    96152                numerator=-1.; denominator=0.;
    … …  
    112168                xfree((void**)&Ginv1);
    113169                xfree((void**)&GinvZ);
    114170        }
    115         printf("\b\b\b\b\b\b\b\b%5.2lf %%\n",100.0);
     171        if(my_thread==0) printf("\b\b\b\b\b\b\b\b%5.2lf %%\n",100.0);
    116172
    117         /*clean-up and Assign output pointer*/
    118         delete variogram;
    119         delete observations;
    120         xfree((void**)&output);
    121         *ppredictions=predictions;
    122         return 1;
    123 }
     173        return NULL;
     174}/*}}}*/
    124175
    125176void ProcessVariogram(Variogram **pvariogram,Options* options){/*{{{*/
    126177

  • proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/modules/Krigingx/Krigingx.h

     
    88#include "../../objects/objects.h"
    99#include "../../toolkits/toolkits.h"
    1010
     11class Observations;
     12class Variogram;
     13
    1114int  Krigingx(double** ppredictions,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp,Options* options);
    1215void ProcessVariogram(Variogram **pvariogram,Options* options);
    1316void GslSolve(double** pX,double* A,double* B,int n);
    1417
     18/*threading: */
     19typedef struct{
     20        int           n_interp;
     21        double       *x_interp;
     22        double       *y_interp;
     23        double        range;
     24        Variogram    *variogram;
     25        Observations *observations;
     26        double       *predictions;
     27}KrigingxThreadStruct;
     28
     29void* Krigingxt(void*);
    1530#endif /* _KRIGINGX_H */

  • proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/modules/Shp2Kmlx/Shp2Kmlx.h

     
    2121#include "../../objects/objects.h"
    2222
    2323/* local prototypes: */
    24 int Shp2Kmlx(char* filshp,char* filkml,
    25                          int sgn);
    26 int Shp2Kmlx(char* filshp,char* filkml,
    27                          int sgn,double cm,double sp);
     24int Shp2Kmlx(char* filshp,char* filkml, int sgn);
     25int Shp2Kmlx(char* filshp,char* filkml, int sgn,double cm,double sp);
    2826
    2927#endif  /* _SHP2KMLX_H */
Note: See TracBrowser for help on using the repository browser.