Changeset 25425


Ignore:
Timestamp:
08/18/20 08:18:44 (5 years ago)
Author:
Mathieu Morlighem
Message:

CHG: enable compilation of Kriging with AD on

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

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/c/classes/IoModel.cpp

    r25379 r25425  
    12631263}
    12641264/*}}}*/
     1265#ifdef _HAVE_AD_
     1266void  IoModel::FetchData(IssmPDouble** pmatrix,int* pM,int* pN,const char* data_name){/*{{{*/
     1267
     1268        /*First, look if has already been loaded (might be an independent variable)*/
     1269        vector<IoData*>::iterator iter;
     1270        for(iter=data.begin();iter<data.end();iter++){
     1271                IoData* iodata=*iter;
     1272                if(strcmp(iodata->name,data_name)==0){
     1273                        _error_(data_name <<" has already been loaded as in IssmDouble and cannot be converted to IssmPDouble");
     1274                }
     1275        }
     1276
     1277        /*output: */
     1278        int          M,N;
     1279        IssmPDouble *matrix = NULL;
     1280        int          code   = 0;
     1281
     1282        /*recover my_rank:*/
     1283        int my_rank=IssmComm::GetRank();
     1284
     1285        /*Set file pointer to beginning of the data: */
     1286        fid=this->SetFilePointerToData(&code,NULL,data_name);
     1287        if(code!=5 && code!=6 && code!=7 && code!=10)_error_("expecting a IssmDouble, integer or boolean matrix for \""<<data_name<<"\""<<" (Code is "<<code<<")");
     1288
     1289        /*Now fetch: */
     1290
     1291        /*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */
     1292        /*numberofelements: */
     1293        if(my_rank==0){
     1294                if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix ");
     1295        }
     1296        ISSM_MPI_Bcast(&M,1,ISSM_MPI_INT,0,IssmComm::GetComm());
     1297
     1298        if(my_rank==0){
     1299                if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix ");
     1300        }
     1301        ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm());
     1302
     1303        /*Now allocate matrix: */
     1304        if(M*N){
     1305                if(code==10){
     1306                        /*Special case for Compressed mat*/
     1307                        IssmPDouble offset,range;
     1308                        if(my_rank==0) if(fread(&offset,sizeof(IssmPDouble),1,fid)!=1) _error_("could not read offset");
     1309                        ISSM_MPI_Bcast(&offset,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
     1310
     1311                        if(my_rank==0) if(fread(&range,sizeof(IssmPDouble),1,fid)!=1) _error_("could not read range");
     1312                        ISSM_MPI_Bcast(&range,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
     1313
     1314                        *pmatrix=xNew<IssmPDouble>(M*N);
     1315
     1316                        /*Read matrix*/
     1317                        uint8_t* rawmatrix=xNew<uint8_t>((M-1)*N);
     1318                        if(my_rank==0) if(fread(rawmatrix,(M-1)*N*sizeof(char),1,fid)!=1) _error_("could not read matrix ");
     1319                        ISSM_MPI_Bcast(rawmatrix,(M-1)*N,ISSM_MPI_CHAR,0,IssmComm::GetComm());
     1320
     1321                        for(int i=0;i<(M-1)*N;++i) (*pmatrix)[i]=offset+range*reCast<IssmPDouble>(rawmatrix[i])/255.;
     1322                        xDelete<uint8_t>(rawmatrix);
     1323
     1324                        /*read time now*/
     1325                        IssmPDouble* timematrix=xNew<IssmPDouble>(N);
     1326                        if(my_rank==0) if(fread(timematrix,N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read time in compressed matrix");
     1327                        ISSM_MPI_Bcast(timematrix,N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
     1328
     1329                        for(int i=0;i<N;++i) (*pmatrix)[(M-1)*N+i]=timematrix[i];
     1330                        xDelete<IssmPDouble>(timematrix);
     1331
     1332                }
     1333                else{
     1334                        /*Read matrix on node 0, then broadcast: */
     1335                        matrix=xNew<IssmPDouble>(M*N);
     1336                        if(my_rank==0) if(fread(matrix,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix ");
     1337                        ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
     1338                        *pmatrix=matrix;
     1339                }
     1340        }
     1341        else{
     1342                *pmatrix=NULL;
     1343        }
     1344        /*Assign output pointers: */
     1345        if(pM) *pM=M;
     1346        if(pN) *pN=N;
     1347}
     1348/*}}}*/
     1349#endif
    12651350void  IoModel::FetchData(IssmDouble*** pmatrices,int** pmdims,int** pndims, int* pnumrecords,const char* data_name){/*{{{*/
    12661351
  • issm/trunk-jpl/src/c/classes/IoModel.h

    r25379 r25425  
    136136                void        FetchData(int** pmatrix,int* pM,int* pN,const char* data_name);
    137137                void        FetchData(IssmDouble**  pscalarmatrix,int* pM,int* pN,const char* data_name);
     138#if _HAVE_AD_  && !defined(_WRAPPERS_)
     139                void        FetchData(IssmPDouble**  pscalarmatrix,int* pM,int* pN,const char* data_name);
     140#endif
    138141                void        FetchData(IssmDouble*** pmatrixarray,int** pmdims,int** pndims, int* pnumrecords,const char* data_name);
    139142                void        FetchData(Options *options,const char* data_name);
  • issm/trunk-jpl/src/c/classes/kriging/Observations.cpp

    r24506 r25425  
    4444
    4545        /*Get tree type (FIXME)*/
    46         IssmDouble dtree = 0.;
     46        double dtree = 0.;
    4747        options->Get(&dtree,"treetype",1.);
    4848        this->treetype = reCast<int>(dtree);
     
    552552        IssmPDouble* B = xNew<IssmPDouble>(n_obs+1);
    553553
    554         IssmDouble unbias = variogram->Covariance(0.,0.);
     554        IssmPDouble unbias = variogram->Covariance(0.,0.);
    555555        /*First: Create semivariogram matrix for observations*/
    556556        for(i=0;i<n_obs;i++){
  • issm/trunk-jpl/src/c/classes/kriging/Observations.h

    r20827 r25425  
    2525                /*constructors, destructors*/
    2626                Observations();
    27                 Observations(IssmDouble* observations_list,IssmDouble* x,IssmDouble* y,int n,Options* options);
     27                Observations(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options);
    2828                ~Observations();
    2929
    3030                /*Initialize data structures*/
    31                 void InitCovertree(IssmDouble* observations_list,IssmDouble* x,IssmDouble* y,int n,Options* options);
    32                 void InitQuadtree(IssmDouble* observations_list,IssmDouble* x,IssmDouble* y,int n,Options* options);
     31                void InitCovertree(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options);
     32                void InitQuadtree(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options);
    3333
    3434                /*Methods*/
    35                 void ClosestObservation(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);
    36                 void ClosestObservationCovertree(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);
    37                 void ClosestObservationQuadtree(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);
     35                void ClosestObservation(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius);
     36                void ClosestObservationCovertree(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius);
     37                void ClosestObservationQuadtree(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius);
    3838                void Distances(IssmPDouble* distances,IssmPDouble *x,IssmPDouble *y,int n,IssmPDouble radius);
    39                 void InterpolationIDW(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,IssmDouble power);
    40                 void InterpolationV4(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata);
    41                 void InterpolationKriging(IssmDouble *pprediction,IssmDouble *perror,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,Variogram* variogram);
    42                 void InterpolationNearestNeighbor(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);
    43                 void ObservationList(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs);
    44                 void ObservationList(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata);
    45                 void ObservationListCovertree(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata);
    46                 void ObservationListQuadtree(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata);
    47                 void QuadtreeColoring(IssmDouble* A,IssmDouble *x,IssmDouble *y,int n);
    48                 void Variomap(IssmDouble* gamma,IssmDouble *x,int n);
     39                void InterpolationIDW(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,IssmPDouble power);
     40                void InterpolationV4(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata);
     41                void InterpolationKriging(IssmPDouble *pprediction,IssmPDouble *perror,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,Variogram* variogram);
     42                void InterpolationNearestNeighbor(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius);
     43                void ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs);
     44                void ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata);
     45                void ObservationListCovertree(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata);
     46                void ObservationListQuadtree(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata);
     47                void QuadtreeColoring(IssmPDouble* A,IssmPDouble *x,IssmPDouble *y,int n);
     48                void Variomap(IssmPDouble* gamma,IssmPDouble *x,int n);
    4949
    5050};
  • issm/trunk-jpl/src/c/main/kriging.cpp

    r21068 r25425  
    77/*Local prototypes*/
    88void ProcessArguments2(char** pbinfilename,char** poutbinfilename,char** plockfilename,char** prootpath,int argc,char **argv);
    9 void ProcessInputfile(IssmDouble **px,IssmDouble **py,IssmDouble **pdata,int *pnobs,IssmDouble **px_interp,IssmDouble **py_interp,int *pninterp,Options **poptions,FILE* fid);
     9void ProcessInputfile(double **px,double **py,double **pdata,int *pnobs,double **px_interp,double **py_interp,int *pninterp,Options **poptions,FILE* fid);
    1010
    1111int main(int argc,char **argv){
     
    2323        /*Input*/
    2424        int         ninterp,nobs;
    25         IssmDouble *x        = NULL;
    26         IssmDouble *y        = NULL;
    27         IssmDouble *data     = NULL;
    28         IssmDouble *x_interp = NULL;
    29         IssmDouble *y_interp = NULL;
     25        double *x        = NULL;
     26        double *y        = NULL;
     27        double *data     = NULL;
     28        double *x_interp = NULL;
     29        double *y_interp = NULL;
    3030        Options    *options  = NULL;
    3131
    3232        /*Output*/
    33         IssmDouble *predictions = NULL;
    34         IssmDouble *error       = NULL;
     33        double *predictions = NULL;
     34        double *error       = NULL;
    3535
    3636        /*Initialize exception trapping: */
     
    7373        xDelete<char>(outbinfilename);
    7474        xDelete<char>(rootpath);
    75         xDelete<IssmDouble>(x);
    76         xDelete<IssmDouble>(y);
    77         xDelete<IssmDouble>(data);
    78         xDelete<IssmDouble>(x_interp);
    79         xDelete<IssmDouble>(y_interp);
    80         xDelete<IssmDouble>(predictions);
    81         xDelete<IssmDouble>(error);
     75        xDelete<double>(x);
     76        xDelete<double>(y);
     77        xDelete<double>(data);
     78        xDelete<double>(x_interp);
     79        xDelete<double>(y_interp);
     80        xDelete<double>(predictions);
     81        xDelete<double>(error);
    8282        delete options;
    8383        delete results;
     
    130130}
    131131
    132 void ProcessInputfile(IssmDouble **px,IssmDouble **py,IssmDouble **pdata,int *pnobs,IssmDouble **px_interp,IssmDouble **py_interp,int *pninterp,Options **poptions,FILE* fid){
     132void ProcessInputfile(double **px,double **py,double **pdata,int *pnobs,double **px_interp,double **py_interp,int *pninterp,Options **poptions,FILE* fid){
    133133
    134         int         ninterp    ,nobs;
    135         IssmDouble *x        = NULL;
    136         IssmDouble *y        = NULL;
    137         IssmDouble *data     = NULL;
    138         IssmDouble *x_interp = NULL;
    139         IssmDouble *y_interp = NULL;
     134        int     ninterp,nobs;
     135        double *x        = NULL;
     136        double *y        = NULL;
     137        double *data     = NULL;
     138        double *x_interp = NULL;
     139        double *y_interp = NULL;
    140140        Options    *options  = NULL;
    141141
     
    144144        iomodel->fid=fid;
    145145        iomodel->CheckFile();
    146         iomodel->FetchData(&x,&M,&N,"md.x");        nobs=M*N;
    147         iomodel->FetchData(&y,&M,&N,"md.y");        _assert_(M*N==nobs);
    148         iomodel->FetchData(&data,&M,&N,"md.data");     _assert_(M*N==nobs);
     146        iomodel->FetchData(&x,&M,&N,"md.x");               nobs=M*N;
     147        iomodel->FetchData(&y,&M,&N,"md.y");               _assert_(M*N==nobs);
     148        iomodel->FetchData(&data,&M,&N,"md.data");         _assert_(M*N==nobs);
    149149        iomodel->FetchData(&x_interp,&M,&N,"md.x_interp"); ninterp=M*N;
    150150        iomodel->FetchData(&y_interp,&M,&N,"md.y_interp"); _assert_(M*N==ninterp);
Note: See TracChangeset for help on using the changeset viewer.