source: issm/oecreview/Archive/12281-12300/ISSM-12291-12292.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: 15.2 KB

  • proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/Container/Options.h

     
    2121                Option* GetOption(const char* name);
    2222                void Get(double*  pvalue,const char* name);
    2323                void Get(double*  pvalue,const char* name,double default_value);
     24                void Get(int*  pvalue,const char* name);
     25                void Get(int*  pvalue,const char* name,int default_value);
    2426                void Get(bool*    pvalue,const char* name);
    2527                void Get(bool*    pvalue,const char* name,bool default_value);
    2628                void Get(char**   pvalue,const char* name);

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

     
    9393
    9494/*Methods*/
    9595/*FUNCTION Observations::ObservationList{{{*/
    96 void Observations::ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double range){
     96void Observations::ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double radius,int maxdata){
    9797
    9898        /*Output and Intermediaries*/
    99         int          nobs,i,index;
     99        bool         stop;
     100        int          nobs,i,j,k,n,counter;
     101        double       h2,radius2;
     102        int         *indices      = NULL;
     103        double      *dists        = NULL;
    100104        double      *x            = NULL;
    101105        double      *y            = NULL;
    102106        double      *obs          = NULL;
    103107        Observation *observation  = NULL;
    104         int         *indices      = NULL;
    105108
    106         /*Treat range*/
    107         if(range==0){
    108                 range=this->quadtree->root->length;
    109         }
     109        /*Compute radius square*/
     110        if(radius==0) radius=this->quadtree->root->length;
     111        radius2 = radius*radius;
    110112
    111         /*Find all observations that are in range*/
    112         nobs=0;
    113         while(nobs==0){
    114                 xfree((void**)&indices);
    115                 this->quadtree->RangeSearch(&indices,&nobs,x_interp,y_interp,range);
    116                 if(!nobs){
    117                         /*No observation found, double range*/
    118                         range=2*range;
     113        /*Find all observations that are in radius*/
     114        indices = (int*)xmalloc(this->Size()*sizeof(int));
     115        dists   = (double*)xmalloc(this->Size()*sizeof(double));
     116        nobs    = 0;
     117
     118        for (i=0;i<this->Size();i++){
     119                observation=(Observation*)this->GetObjectByOffset(i);
     120                h2 = (observation->x-x_interp)*(observation->x-x_interp) + (observation->y-y_interp)*(observation->y-y_interp);
     121
     122                if(nobs==maxdata && h2>radius2) continue;
     123                if(nobs<=maxdata){
     124                        indices[nobs]   = i;
     125                        dists[nobs]     = h2;
     126                        nobs++;
    119127                }
    120         }
    121         /*Check that we do not have more than 50 observations*/
    122         //while(nobs>50){
    123         //      range=range/2;
    124         //      xfree((void**)&indices);
    125         //      this->quadtree->RangeSearch(&indices,&nobs,x_interp,y_interp,range);
    126         //}
     128                if(nobs==1) continue;
    127129
    128         if(!nobs) _error_("Unexpected error: no observation within current range and more than 50 observations if doubling range");
     130                /*Sort all dists up to now*/
     131                n=nobs-1;
     132                stop = false;
     133                for(k=0;k<n-1;k++){
     134                        if(h2<dists[k]){
     135                                counter=1;
     136                                for(int jj=k;jj<n;jj++){
     137                                        j  = n-counter;
     138                                        dists[j+1]   = dists[j];
     139                                        indices[j+1] = indices[j];
     140                                        counter++;
     141                                }
     142                                dists[k]   = h2;
     143                                indices[k] = i;
     144                                stop = true;
     145                                break;
     146                        }
     147                        if(stop) break;
     148                }
     149        } 
     150        xfree((void**)&dists);
    129151
    130         /*Allocate vectors*/
    131         x   = (double*)xmalloc(nobs*sizeof(double));
    132         y   = (double*)xmalloc(nobs*sizeof(double));
    133         obs = (double*)xmalloc(nobs*sizeof(double));
     152        if(nobs){
     153                /*Allocate vectors*/
     154                x   = (double*)xmalloc(nobs*sizeof(double));
     155                y   = (double*)xmalloc(nobs*sizeof(double));
     156                obs = (double*)xmalloc(nobs*sizeof(double));
    134157
    135         /*Loop over all observations and fill in x, y and obs*/
    136         for (i=0;i<nobs;i++){
    137                 observation=(Observation*)this->GetObjectByOffset(indices[i]);
    138                 observation->WriteXYObs(&x[i],&y[i],&obs[i]);
     158                /*Loop over all observations and fill in x, y and obs*/
     159                for (i=0;i<nobs;i++){
     160                        observation=(Observation*)this->GetObjectByOffset(indices[i]);
     161                        observation->WriteXYObs(&x[i],&y[i],&obs[i]);
     162                }
    139163        }
    140164
    141165        /*Assign output pointer*/

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

     
    2222                ~Observations();
    2323
    2424                /*Methods*/
    25                 void ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double range);
     25                void ObservationList(double **px,double **py,double **pobs,int* pnobs,double x_interp,double y_interp,double radius,int maxdata);
    2626                void QuadtreeColoring(double* A,double *x,double *y,int n);
    2727                void Variomap(double* gamma,double *x,int n);
    2828

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

     
    6666        return 1;
    6767}
    6868/*}}}*/
     69/*FUNCTION Options::Get(int* pvalue, char* name){{{1*/
     70void Options::Get(int* pvalue,const char* name){
     71
     72        vector<Object*>::iterator object;
     73        Option* option=NULL;
     74
     75        /*Get option*/
     76        option=GetOption(name);
     77
     78        /*If the pointer is not NULL, the option has been found*/
     79        if(option){
     80                option->Get(pvalue);
     81        }
     82        /*Else, the Option does not exist, no default provided*/
     83        else{
     84                _error_("option of name \"%s\" not found, and no default value has been provided",name);
     85        }
     86}
     87/*}}}*/
     88/*FUNCTION Options::Get(int* pvalue, char* name,int default_value){{{1*/
     89void Options::Get(int* pvalue,const char* name,int default_value){
     90
     91        vector<Object*>::iterator object;
     92        Option* option=NULL;
     93
     94        /*Get option*/
     95        option=GetOption(name);
     96
     97        /*If the pointer is not NULL, the option has been found*/
     98        if(option){
     99                option->Get(pvalue);
     100        }
     101        /*Else, the Option does not exist, a default is provided here*/
     102        else{
     103                *pvalue=default_value;
     104        }
     105}
     106/*}}}*/
    69107/*FUNCTION Options::Get(double* pvalue, char* name){{{1*/
    70108void Options::Get(double* pvalue,const char* name){
    71109

  • proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/shared/Threads/LaunchThread.cpp

     
    6666       
    6767        function((void*)&handle);
    6868        #endif
    69 
    7069}

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

     
    2323        double *error       = NULL;
    2424
    2525        /*Intermediaries*/
    26         double        range;
     26        int           mindata,maxdata;
     27        double        radius;
    2728        char         *output       = NULL;
    2829        Variogram    *variogram    = NULL;
    2930        Observations *observations = NULL;
    … …  
    3738
    3839        /*Get Variogram from Options*/
    3940        ProcessVariogram(&variogram,options);
    40         options->Get(&range,"searchrange",0.);
     41        options->Get(&radius,"searchradius",0.);
     42        options->Get(&mindata,"mindata",1);
     43        options->Get(&maxdata,"maxdata",50);
    4144
    4245        /*Process observation dataset*/
    4346        observations=new Observations(obs_list,obs_x,obs_y,obs_length,options);
    … …  
    6164                gate.n_interp     = n_interp;
    6265                gate.x_interp     = x_interp;
    6366                gate.y_interp     = y_interp;
    64                 gate.range        = range;
     67                gate.radius       = radius;
     68                gate.mindata      = mindata;
     69                gate.maxdata      = maxdata;
    6570                gate.variogram    = variogram;
    6671                gate.observations = observations;
    6772                gate.predictions  = predictions;
    … …  
    105110        int           n_interp     = gate->n_interp;
    106111        double       *x_interp     = gate->x_interp;
    107112        double       *y_interp     = gate->y_interp;
    108         double        range        = gate->range;
     113        double        radius       = gate->radius;
     114        int           mindata      = gate->mindata;
     115        int           maxdata      = gate->maxdata;
    109116        Variogram    *variogram    = gate->variogram;
    110117        Observations *observations = gate->observations;
    111118        double       *predictions  = gate->predictions;
    … …  
    136143                printf("\r      interpolation progress: %5.2lf%%",localpercent);
    137144
    138145                /*Get list of observations for current point*/
    139                 observations->ObservationList(&x,&y,&obs,&n_obs,x_interp[idx],y_interp[idx],range);
    140                 if(n_obs==0){
    141                         predictions[idx] = -999.;
    142                         error[idx]       = -999.;
     146                observations->ObservationList(&x,&y,&obs,&n_obs,x_interp[idx],y_interp[idx],radius,maxdata);
     147                if(n_obs<mindata){
     148                        predictions[idx] = -999.0;
     149                        error[idx]       = -999.0;
    143150                        continue;
    144151                }
    145152

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

     
    2020        int           n_interp;
    2121        double       *x_interp;
    2222        double       *y_interp;
    23         double        range;
     23        double        radius;
     24        int           mindata;
     25        int           maxdata;
    2426        Variogram    *variogram;
    2527        Observations *observations;
    2628        double       *predictions;

  • proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/OptionStruct.h

     
    3838                int   NumEl();
    3939                int   NDims();
    4040                int*  Size();
     41                void  Get(int* pvalue){_error_("An OptionStruct object cannot return a int");};
    4142                void  Get(double* pvalue){_error_("An OptionStruct object cannot return a double");};
    4243                void  Get(bool* pvalue){  _error_("An OptionStruct object cannot return a bool");};
    4344                void  Get(char** pvalue){ _error_("An OptionStruct object cannot return a string");};

  • proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/OptionDouble.cpp

     
    119119        return(Option::Size());
    120120}
    121121/*}}}*/
     122/*FUNCTION OptionDouble::Get(int* pvalue) {{{1*/
     123void OptionDouble::Get(int* pvalue){
     124
     125        /*We should first check that the size is one*/
     126        if(this->NumEl()!=1){
     127                _error_("option \"%s\" has %i elements and cannot return a single int",this->name,this->NumEl());
     128        }
     129
     130        /*Assign output pointer*/
     131        *pvalue=(int)this->values[0];
     132}
     133/*}}}*/
    122134/*FUNCTION OptionDouble::Get(double* pvalue) {{{1*/
    123135void OptionDouble::Get(double* pvalue){
    124136

  • proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/OptionDouble.h

     
    3838                int   NumEl();
    3939                int   NDims();
    4040                int*  Size();
     41                void  Get(int* pvalue);
    4142                void  Get(double* pvalue);
    4243                void  Get(bool* pvalue){  _error_("An OptionDouble object cannot return a bool");};
    4344                void  Get(char** pvalue){ _error_("An OptionDouble object cannot return a string");};

  • proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/OptionLogical.h

     
    3838                int   NumEl();
    3939                int   NDims();
    4040                int*  Size();
     41                void  Get(int* pvalue){_error_("An OptionLogical object cannot return a int");};
    4142                void  Get(double* pvalue){_error_("An OptionLogical object cannot return a double");};
    4243                void  Get(bool* pvalue);
    4344                void  Get(char** pvalue){ _error_("An OptionLogical object cannot return a string");};

  • proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/OptionChar.h

     
    3838                int   NumEl();
    3939                int   NDims();
    4040                int*  Size();
     41                void  Get(int* pvalue){_error_("An OptionChar object cannot return a int");};
    4142                void  Get(double* pvalue){_error_("An OptionChar object cannot return a double");};
    4243                void  Get(bool* pvalue){  _error_("An OptionChar object cannot return a bool");};
    4344                void  Get(char** pvalue);

  • proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/Option.h

     
    4141                virtual int   NumEl()=0;
    4242                virtual int   NDims()=0;
    4343                virtual int*  Size()=0;
     44                virtual void  Get(int* pvalue)=0;
    4445                virtual void  Get(double* pvalue)=0;
    4546                virtual void  Get(bool* pvalue)=0;
    4647                virtual void  Get(char** pvalue)=0;

  • proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Options/OptionCell.h

     
    3838                int   NumEl();
    3939                int   NDims();
    4040                int*  Size();
     41                void  Get(int* pvalue){_error_("An OptionCell object cannot return a int");};
    4142                void  Get(double* pvalue){_error_("An OptionCell object cannot return a double");};
    4243                void  Get(bool* pvalue){  _error_("An OptionCell object cannot return a bool");};
    4344                void  Get(char** pvalue){ _error_("An OptionCell object cannot return a string");};

  • proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Kriging/GaussianVariogram.cpp

     
    7474        a     = 1./3.;
    7575        cova = (sill-nugget)*exp(-h2/(a*range*range));
    7676
    77         //if(h2<0.0000001) cova = sill;
     77        if(h2<0.0000001) cova = sill;
    7878
    7979        return cova;
    8080}
Note: See TracBrowser for help on using the repository browser.