kriging.cpp

Go to the documentation of this file.

 
#include "./issm.h"
 
/*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);
 
int main(int argc,char **argv){
 
   /*I/O: */
   FILE *output_fid = NULL;
   FILE *input_fid  = NULL;
 
   /*File names*/
   char *lockfilename   = NULL;
   char *binfilename    = NULL;
   char *outbinfilename = NULL;
   char *rootpath       = NULL;
 
   /*Input*/
   int         ninterp,nobs;
   IssmDouble *x        = NULL;
   IssmDouble *y        = NULL;
   IssmDouble *data     = NULL;
   IssmDouble *x_interp = NULL;
   IssmDouble *y_interp = NULL;
   Options    *options  = NULL;
 
   /*Output*/
   IssmDouble *predictions = NULL;
   IssmDouble *error       = NULL;
 
   /*Initialize exception trapping: */
   ExceptionTrapBegin();
 
   /*Initialize environment (MPI, PETSC, MUMPS, etc ...)*/
   ISSM_MPI_Comm comm=EnvironmentInit(argc,argv);
   IssmComm::SetComm(comm);
 
   ProcessArguments2(&binfilename,&outbinfilename,&lockfilename,&rootpath,argc,argv);
 
   /*Process input files*/
   input_fid=pfopen(binfilename,"rb");
   ProcessInputfile(&x,&y,&data,&nobs,&x_interp,&y_interp,&ninterp,&options,input_fid);
   pfclose(input_fid,binfilename);
 
   _printf0_("call computational core:\n");
   pKrigingx(&predictions,&error,x,y,data,nobs,x_interp,y_interp,ninterp,options);
 
   _printf0_("write results to disk:\n");
   Results *results = new Results();
   if(IssmComm::GetRank()==0){
      output_fid=pfopen0(outbinfilename,"wb");
      results->AddObject(new GenericExternalResult<double*>(results->Size()+1,"predictions",predictions,ninterp,1,1,0));
      results->AddObject(new GenericExternalResult<double*>(results->Size()+1,"error",error,ninterp,1,1,0));
      for(int i=0;i<results->Size();i++){
         ExternalResult* result=xDynamicCast<ExternalResult*>(results->GetObjectByOffset(i));
         result->WriteData(output_fid,1);
      }
      pfclose(output_fid,outbinfilename);
   }
 
   /*Close output and toolkits options file and write lock file if requested*/
   _printf0_("write lock file:\n");
   WriteLockFile(lockfilename);
 
   /*Free ressources */
   xDelete<char>(lockfilename);
   xDelete<char>(binfilename);
   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);
   delete options;
   delete results;
 
   /*Finalize environment:*/
   EnvironmentFinalize();
 
   /*Finalize exception trapping: */
   ExceptionTrapEnd();
 
   return 0; //unix success return;
}
 
void ProcessArguments2(char** pbinfilename,char** poutbinfilename,char** plockfilename,char** prootpath,int argc,char **argv){
   char *modelname      = NULL;
   char *binfilename    = NULL;
   char *outbinfilename = NULL;
   char *lockfilename   = NULL;
   char *rootpatharg    = NULL;
   char *rootpath       = NULL;
 
   if(argc<1)_error_("Usage error: no execution path provided");
   if(argc<2)_error_("Usage error: missing model name");
 
   rootpatharg=argv[1];
   if(strcmp(strstr(rootpatharg,"/"),"/")!=0){ 
      rootpath       = xNew<char>(strlen(rootpatharg)+2); sprintf(rootpath,"%s/",rootpatharg);
   }  
   else{
      rootpath       = xNew<char>(strlen(rootpatharg)+1); sprintf(rootpath,"%s",rootpatharg);
   } 
 
   modelname=argv[2];
   if(strstr(modelname,rootpath)==NULL){
      binfilename    = xNew<char>(strlen(rootpath)+strlen(modelname)+strlen(".bin")   +1); sprintf(binfilename,   "%s%s%s",rootpath,modelname,".bin");
      outbinfilename = xNew<char>(strlen(rootpath)+strlen(modelname)+strlen(".outbin")+1); sprintf(outbinfilename,"%s%s%s",rootpath,modelname,".outbin");
      lockfilename   = xNew<char>(strlen(rootpath)+strlen(modelname)+strlen(".lock")  +1); sprintf(lockfilename,  "%s%s%s",rootpath,modelname,".lock");
   }
   else{
      binfilename    = xNew<char>(strlen(modelname)+strlen(".bin")   +1); sprintf(binfilename,   "%s%s",modelname,".bin");
      outbinfilename = xNew<char>(strlen(modelname)+strlen(".outbin")+1); sprintf(outbinfilename,"%s%s",modelname,".outbin");
      lockfilename   = xNew<char>(strlen(modelname)+strlen(".lock")  +1); sprintf(lockfilename,  "%s%s",modelname,".lock");
   }
 
   /*Clean up and assign output pointer*/
   *pbinfilename=binfilename;
   *poutbinfilename=outbinfilename;
   *plockfilename=lockfilename;
   *prootpath=rootpath;
}
 
void ProcessInputfile(IssmDouble **px,IssmDouble **py,IssmDouble **pdata,int *pnobs,IssmDouble **px_interp,IssmDouble **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;
   Options    *options  = NULL;
 
   int      M,N;
   IoModel* iomodel = new IoModel();
   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_interp,&M,&N,"md.x_interp"); ninterp=M*N;
   iomodel->FetchData(&y_interp,&M,&N,"md.y_interp"); _assert_(M*N==ninterp);
 
   /*Read options*/
   options = new Options();
   iomodel->FetchData(options,"md.y_interp");
 
   /*Assign output pointer*/
   *px        = x;
   *py        = y;
   *pdata     = data;
   *pnobs     = nobs;
   *px_interp = x_interp;
   *py_interp = y_interp;
   *pninterp  = ninterp;
   *poptions  = options;
   delete iomodel;
}