Changeset 25425

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:

• classes/IoModel.cpp (modified) (1 diff)
• classes/IoModel.h (modified) (1 diff)
• classes/kriging/Observations.cpp (modified) (2 diffs)
• classes/kriging/Observations.h (modified) (1 diff)
• main/kriging.cpp (modified) (5 diffs)

issm/trunk-jpl/src/c/classes/IoModel.cpp

@@ -1263 +1263 @@
 }
 /*}}}*/
+#ifdef _HAVE_AD_
+void  IoModel::FetchData(IssmPDouble** pmatrix,int* pM,int* pN,const char* data_name){/*{{{*/
+
+        /*First, look if has already been loaded (might be an independent variable)*/
+        vector<IoData*>::iterator iter;
+        for(iter=data.begin();iter<data.end();iter++){
+                IoData* iodata=*iter;
+                if(strcmp(iodata->name,data_name)==0){
+                        _error_(data_name <<" has already been loaded as in IssmDouble and cannot be converted to IssmPDouble");
+                }
+        }
+
+        /*output: */
+        int          M,N;
+        IssmPDouble *matrix = NULL;
+        int          code   = 0;
+
+        /*recover my_rank:*/
+        int my_rank=IssmComm::GetRank();
+
+        /*Set file pointer to beginning of the data: */
+        fid=this->SetFilePointerToData(&code,NULL,data_name);
+        if(code!=5 && code!=6 && code!=7 && code!=10)_error_("expecting a IssmDouble, integer or boolean matrix for \""<<data_name<<"\""<<" (Code is "<<code<<")");
+
+        /*Now fetch: */
+
+        /*We have to read a matrix from disk. First read the dimensions of the matrix, then the whole matrix: */
+        /*numberofelements: */
+        if(my_rank==0){
+                if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows for matrix ");
+        }
+        ISSM_MPI_Bcast(&M,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+
+        if(my_rank==0){
+                if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns for matrix ");
+        }
+        ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+
+        /*Now allocate matrix: */
+        if(M*N){
+                if(code==10){
+                        /*Special case for Compressed mat*/
+                        IssmPDouble offset,range;
+                        if(my_rank==0) if(fread(&offset,sizeof(IssmPDouble),1,fid)!=1) _error_("could not read offset");
+                        ISSM_MPI_Bcast(&offset,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
+
+                        if(my_rank==0) if(fread(&range,sizeof(IssmPDouble),1,fid)!=1) _error_("could not read range");
+                        ISSM_MPI_Bcast(&range,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
+
+                        *pmatrix=xNew<IssmPDouble>(M*N);
+
+                        /*Read matrix*/
+                        uint8_t* rawmatrix=xNew<uint8_t>((M-1)*N);
+                        if(my_rank==0) if(fread(rawmatrix,(M-1)*N*sizeof(char),1,fid)!=1) _error_("could not read matrix ");
+                        ISSM_MPI_Bcast(rawmatrix,(M-1)*N,ISSM_MPI_CHAR,0,IssmComm::GetComm());
+
+                        for(int i=0;i<(M-1)*N;++i) (*pmatrix)[i]=offset+range*reCast<IssmPDouble>(rawmatrix[i])/255.;
+                        xDelete<uint8_t>(rawmatrix);
+
+                        /*read time now*/
+                        IssmPDouble* timematrix=xNew<IssmPDouble>(N);
+                        if(my_rank==0) if(fread(timematrix,N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read time in compressed matrix");
+                        ISSM_MPI_Bcast(timematrix,N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
+
+                        for(int i=0;i<N;++i) (*pmatrix)[(M-1)*N+i]=timematrix[i];
+                        xDelete<IssmPDouble>(timematrix);
+
+                }
+                else{
+                        /*Read matrix on node 0, then broadcast: */
+                        matrix=xNew<IssmPDouble>(M*N);
+                        if(my_rank==0) if(fread(matrix,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix ");
+                        ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
+                        *pmatrix=matrix;
+                }
+        }
+        else{
+                *pmatrix=NULL;
+        }
+        /*Assign output pointers: */
+        if(pM) *pM=M;
+        if(pN) *pN=N;
+}
+/*}}}*/
+#endif
 void  IoModel::FetchData(IssmDouble*** pmatrices,int** pmdims,int** pndims, int* pnumrecords,const char* data_name){/*{{{*/
 

issm/trunk-jpl/src/c/classes/IoModel.h

@@ -136 +136 @@
                 void        FetchData(int** pmatrix,int* pM,int* pN,const char* data_name);
                 void        FetchData(IssmDouble**  pscalarmatrix,int* pM,int* pN,const char* data_name);
+#if _HAVE_AD_  && !defined(_WRAPPERS_)
+                void        FetchData(IssmPDouble**  pscalarmatrix,int* pM,int* pN,const char* data_name);
+#endif
                 void        FetchData(IssmDouble*** pmatrixarray,int** pmdims,int** pndims, int* pnumrecords,const char* data_name);
                 void        FetchData(Options *options,const char* data_name);

issm/trunk-jpl/src/c/classes/kriging/Observations.cpp

@@ -44 +44 @@
 
         /*Get tree type (FIXME)*/
-        IssmDouble dtree = 0.;
+        double dtree = 0.;
         options->Get(&dtree,"treetype",1.);
         this->treetype = reCast<int>(dtree);
@@ -552 +552 @@
         IssmPDouble* B = xNew<IssmPDouble>(n_obs+1);
 
-        IssmDouble unbias = variogram->Covariance(0.,0.);
+        IssmPDouble unbias = variogram->Covariance(0.,0.);
         /*First: Create semivariogram matrix for observations*/
         for(i=0;i<n_obs;i++){

issm/trunk-jpl/src/c/classes/kriging/Observations.h

@@ -25 +25 @@
                 /*constructors, destructors*/
                 Observations();
-                Observations(IssmDouble* observations_list,IssmDouble* x,IssmDouble* y,int n,Options* options);
+                Observations(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options);
                 ~Observations();
 
                 /*Initialize data structures*/
-                void InitCovertree(IssmDouble* observations_list,IssmDouble* x,IssmDouble* y,int n,Options* options);
-                void InitQuadtree(IssmDouble* observations_list,IssmDouble* x,IssmDouble* y,int n,Options* options);
+                void InitCovertree(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options);
+                void InitQuadtree(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options);
 
                 /*Methods*/
-                void ClosestObservation(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);
-                void ClosestObservationCovertree(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);
-                void ClosestObservationQuadtree(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);
+                void ClosestObservation(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius);
+                void ClosestObservationCovertree(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius);
+                void ClosestObservationQuadtree(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius);
                 void Distances(IssmPDouble* distances,IssmPDouble *x,IssmPDouble *y,int n,IssmPDouble radius);
-                void InterpolationIDW(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,IssmDouble power);
-                void InterpolationV4(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata);
-                void InterpolationKriging(IssmDouble *pprediction,IssmDouble *perror,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,Variogram* variogram);
-                void InterpolationNearestNeighbor(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);
-                void ObservationList(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs);
-                void ObservationList(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata);
-                void ObservationListCovertree(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata);
-                void ObservationListQuadtree(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata);
-                void QuadtreeColoring(IssmDouble* A,IssmDouble *x,IssmDouble *y,int n);
-                void Variomap(IssmDouble* gamma,IssmDouble *x,int n);
+                void InterpolationIDW(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,IssmPDouble power);
+                void InterpolationV4(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata);
+                void InterpolationKriging(IssmPDouble *pprediction,IssmPDouble *perror,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,Variogram* variogram);
+                void InterpolationNearestNeighbor(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius);
+                void ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs);
+                void ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata);
+                void ObservationListCovertree(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata);
+                void ObservationListQuadtree(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata);
+                void QuadtreeColoring(IssmPDouble* A,IssmPDouble *x,IssmPDouble *y,int n);
+                void Variomap(IssmPDouble* gamma,IssmPDouble *x,int n);
 
 };

issm/trunk-jpl/src/c/main/kriging.cpp

@@ -7 +7 @@
 /*Local prototypes*/
 void ProcessArguments2(char** pbinfilename,char** poutbinfilename,char** plockfilename,char** prootpath,int argc,char **argv);
-void ProcessInputfile(IssmDouble **px,IssmDouble **py,IssmDouble **pdata,int *pnobs,IssmDouble **px_interp,IssmDouble **py_interp,int *pninterp,Options **poptions,FILE* fid);
+void ProcessInputfile(double **px,double **py,double **pdata,int *pnobs,double **px_interp,double **py_interp,int *pninterp,Options **poptions,FILE* fid);
 
 int main(int argc,char **argv){
@@ -23 +23 @@
         /*Input*/
         int         ninterp,nobs;
-        IssmDouble *x        = NULL;
-        IssmDouble *y        = NULL;
-        IssmDouble *data     = NULL;
-        IssmDouble *x_interp = NULL;
-        IssmDouble *y_interp = NULL;
+        double *x        = NULL;
+        double *y        = NULL;
+        double *data     = NULL;
+        double *x_interp = NULL;
+        double *y_interp = NULL;
         Options    *options  = NULL;
 
         /*Output*/
-        IssmDouble *predictions = NULL;
-        IssmDouble *error       = NULL;
+        double *predictions = NULL;
+        double *error       = NULL;
 
         /*Initialize exception trapping: */
@@ -73 +73 @@
         xDelete<char>(outbinfilename);
         xDelete<char>(rootpath);
-        xDelete<IssmDouble>(x);
-        xDelete<IssmDouble>(y);
-        xDelete<IssmDouble>(data);
-        xDelete<IssmDouble>(x_interp);
-        xDelete<IssmDouble>(y_interp);
-        xDelete<IssmDouble>(predictions);
-        xDelete<IssmDouble>(error);
+        xDelete<double>(x);
+        xDelete<double>(y);
+        xDelete<double>(data);
+        xDelete<double>(x_interp);
+        xDelete<double>(y_interp);
+        xDelete<double>(predictions);
+        xDelete<double>(error);
         delete options;
         delete results;
@@ -130 +130 @@
 }
 
-void ProcessInputfile(IssmDouble **px,IssmDouble **py,IssmDouble **pdata,int *pnobs,IssmDouble **px_interp,IssmDouble **py_interp,int *pninterp,Options **poptions,FILE* fid){
+void ProcessInputfile(double **px,double **py,double **pdata,int *pnobs,double **px_interp,double **py_interp,int *pninterp,Options **poptions,FILE* fid){
 
-        int         ninterp    ,nobs;
-        IssmDouble *x        = NULL;
-        IssmDouble *y        = NULL;
-        IssmDouble *data     = NULL;
-        IssmDouble *x_interp = NULL;
-        IssmDouble *y_interp = NULL;
+        int     ninterp,nobs;
+        double *x        = NULL;
+        double *y        = NULL;
+        double *data     = NULL;
+        double *x_interp = NULL;
+        double *y_interp = NULL;
         Options    *options  = NULL;
 
@@ -144 +144 @@
         iomodel->fid=fid;
         iomodel->CheckFile();
-        iomodel->FetchData(&x,&M,&N,"md.x");        nobs=M*N;
-        iomodel->FetchData(&y,&M,&N,"md.y");        _assert_(M*N==nobs);
-        iomodel->FetchData(&data,&M,&N,"md.data");     _assert_(M*N==nobs);
+        iomodel->FetchData(&x,&M,&N,"md.x");               nobs=M*N;
+        iomodel->FetchData(&y,&M,&N,"md.y");               _assert_(M*N==nobs);
+        iomodel->FetchData(&data,&M,&N,"md.data");         _assert_(M*N==nobs);
         iomodel->FetchData(&x_interp,&M,&N,"md.x_interp"); ninterp=M*N;
         iomodel->FetchData(&y_interp,&M,&N,"md.y_interp"); _assert_(M*N==ninterp);