source: issm/branches/trunk-jpl-damage/src/c/solutions/kriging.cpp@ 13101

/*!\file:  kriging.cpp
 * \brief: kriging main parallel program
 */ 

#include "../issm.h"
#include "../include/globals.h"

/*Local prototypes*/
void ProcessArguments2(char** pbinfilename,char** poutbinfilename,char** plockfilename,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;
        bool  waitonlock = false;

        /*File names*/
        char *lockfilename   = NULL;
        char *binfilename    = NULL;
        char *outbinfilename = 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;

        ISSMBOOT();

        /*Initialize environments: Petsc, MPI, etc...: */
#ifdef _HAVE_PETSC_
        int ierr=PetscInitialize(&argc,&argv,(char*)0,"");  
        if(ierr) _error_("Could not initialize Petsc");
#else
#ifdef _HAVE_MPI_
        MPI_Init(&argc,&argv);
#endif
#endif

        /*Size and rank: */
#ifdef _HAVE_MPI_
        MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);  
        MPI_Comm_size(MPI_COMM_WORLD,&num_procs); 
#endif

        /*First process inputs*/
        _pprintLine_("");
        _pprintLine_("Ice Sheet System Model (" << PACKAGE_NAME << ") version " << PACKAGE_VERSION);
        _pprintLine_("(website: " << PACKAGE_URL << " contact: " << PACKAGE_BUGREPORT << ")");
        _pprintLine_("");
        ProcessArguments2(&binfilename,&outbinfilename,&lockfilename,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);

        _pprintLine_("call computational core:");
        pKrigingx(&predictions,&error,x,y,data,nobs,x_interp,y_interp,ninterp,options);

        _pprintLine_("write results to disk:");
        Results *results = new Results();
        if(my_rank==0){
                output_fid=pfopen(outbinfilename,"wb");
                results->AddObject(new DoubleVecExternalResult(results->Size()+1,0,predictions,ninterp,1,0));
                results->AddObject(new DoubleVecExternalResult(results->Size()+1,1,error,ninterp,1,0));
                for(int i=0;i<results->Size();i++){
                        ExternalResult* result=(ExternalResult*)results->GetObjectByOffset(i);
                        result->WriteData(output_fid,1);
                }
                pfclose(output_fid,outbinfilename);
        }

        /*Close output and petsc options file and write lock file if requested*/
        _pprintLine_("write lock file:");
        WriteLockFile(lockfilename);

        /*Free ressources */
        xDelete<char>(lockfilename);
        xDelete<char>(binfilename);
        xDelete<char>(outbinfilename);
        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;

#ifdef _HAVE_PETSC_
        _pprintLine_("closing MPI and Petsc");
        PetscFinalize(); 
#else
#ifdef _HAVE_MPI_
        _pprintLine_("closing MPI and Petsc");
        MPI_Finalize();
#endif
#endif

        /*end module: */
        ISSMEND();

        return 0; //unix success return;
}

void ProcessArguments2(char** pbinfilename,char** poutbinfilename,char** plockfilename,int argc,char **argv){

        char *modelname      = NULL;
        char *binfilename    = NULL;
        char *outbinfilename = NULL;
        char *lockfilename   = NULL;

        if(argc<2)_error_("Usage error: missing model name");
        modelname=argv[2];
        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;
}

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,numoptions;
        IssmDouble  *x        = NULL;
        IssmDouble  *y        = NULL;
        IssmDouble  *data     = NULL;
        IssmDouble  *x_interp = NULL;
        IssmDouble  *y_interp = NULL;
        Options *options  = NULL;
        Option  *option   = NULL;

        int dummy,M,N;
        IoModel* iomodel = new IoModel();
        iomodel->fid=fid;
        iomodel->CheckEnumSync();
        iomodel->FetchData(&x,&M,&N,0);        nobs=M*N;
        iomodel->FetchData(&y,&M,&N,1);        _assert_(M*N==nobs);
        iomodel->FetchData(&data,&M,&N,2);     _assert_(M*N==nobs);
        iomodel->FetchData(&x_interp,&M,&N,3); ninterp=M*N;
        iomodel->FetchData(&y_interp,&M,&N,4); _assert_(M*N==ninterp);

        /*Read options*/
        options = new Options();
        iomodel->LastIndex(&N);
        numoptions=(int)((N-4)/2);
        for(int i=0;i<numoptions;i++){
                iomodel->FetchData(&option,5+2*i);
                options->AddOption(option);
                option=NULL;
        }

        /*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;
}