Index: /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_dakota.cpp =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_dakota.cpp (revision 25614) +++ /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_dakota.cpp (revision 25615) @@ -4,5 +4,4 @@ #include "./issm.h" -#include /*Dakota includes: */ @@ -14,8 +13,4 @@ #include "DakotaInterface.hpp" #endif - -/*prototypes:*/ -bool dirstructure(int argc,char** argv); -int issm_dakota_statistics(int argc,char** argv); int main(int argc,char **argv){ /*{{{*/ @@ -46,5 +41,5 @@ /*Create directory structure for model outputs:*/ - statistics=dirstructure(argc,argv); + statistics=DakotaDirStructure(argc,argv); /* Parse input and construct Dakota LibraryEnvironment, performing input data checks*/ @@ -93,5 +88,5 @@ /* Run statistics if requested:*/ - if(statistics)issm_dakota_statistics(argc,argv); + if(statistics)DakotaStatistics(argc,argv); /*free allocations:*/ @@ -108,210 +103,2 @@ } /*}}}*/ -bool dirstructure(int argc,char** argv){ /*{{{*/ - - char* input_file; - FILE* fid; - IoModel* iomodel=NULL; - int check; - - //qmu statistics - bool statistics = false; - int numdirectories = 0; - - /*First things first, set the communicator as a global variable: */ - IssmComm::SetComm(MPI_COMM_WORLD); - - /*Barrier:*/ - ISSM_MPI_Barrier(IssmComm::GetComm()); - _printf0_("Preparing directory structure for model outputs:" << "\n"); - - //open model input file for reading - input_file=xNew((strlen(argv[2])+strlen(argv[3])+strlen(".bin")+2)); - sprintf(input_file,"%s/%s%s",argv[2],argv[3],".bin"); - fid=fopen(input_file,"rb"); - if (fid==NULL) Cerr << "dirstructure error message: could not open model " << input_file << " to retrieve qmu statistics parameters" << std::endl; - - //initialize IoModel, but light version, we just need it to fetch one constant: - iomodel=new IoModel(); - iomodel->fid=fid; - iomodel->FetchConstants(); - - //early return if statistics not requested: - iomodel->FindConstant(&statistics,"md.qmu.statistics"); - if(!statistics){ - delete iomodel; - fclose(fid); - return false; //important return value! - } - - iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories"); - - /*Ok, we have everything we need to create the directory structure:*/ - if(IssmComm::GetRank()==0){ - for (int i=0;i((strlen(argv[2])+strlen(argv[3])+strlen(".bin")+2)); - sprintf(input_file,"%s/%s%s",argv[2],argv[3],".bin"); - fid=fopen(input_file,"rb"); - if (fid==NULL) Cerr << "issm_dakota_statistics error message: could not open model " << input_file << " to retrieve qmu statistics parameters" << std::endl; - - //initialize IoModel, but light version, we'll need it to fetch constants: - iomodel=new IoModel(); - iomodel->fid=fid; - iomodel->FetchConstants(); - - //early return if statistics not requested: - iomodel->FindConstant(&statistics,"md.qmu.statistics"); - if(!statistics){ - delete iomodel; - fclose(fid); - return 0; - } - - //create parameters datasets with al the qmu statistics settings we need: - if(statistics){ - - /*Initialize parameters and results:*/ - results = new Results(); - parameters=new Parameters(); - - //solution type: - parameters->AddObject(new IntParam(SolutionTypeEnum,StatisticsSolutionEnum)); - - //root directory - directory=xNew(strlen(argv[2])+1); - xMemCpy(directory,argv[2],strlen(argv[2])+1); - parameters->AddObject(new StringParam(DirectoryNameEnum,directory)); - - //model name - model=xNew(strlen(argv[3])+1); - xMemCpy(model,argv[3],strlen(argv[3])+1); - parameters->AddObject(new StringParam(InputFileNameEnum,model)); - - //nsamples - iomodel->FindConstant(&nsamples,"md.qmu.method.params.samples"); - parameters->AddObject(new IntParam(QmuNsampleEnum,nsamples)); - - //ndirectories - iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories"); - parameters->AddObject(new IntParam(QmuNdirectoriesEnum,numdirectories)); - - //nfiles per directory - iomodel->FindConstant(&nfilesperdirectory,"md.qmu.statistics.nfiles_per_directory"); - parameters->AddObject(new IntParam(QmuNfilesPerDirectoryEnum,nfilesperdirectory)); - - //At this point, we don't want to go forward any longer, we want to create an MPI - //communicator on which to carry out the computations: - if ((my_rank+1)*nfilesperdirectory>nsamples)color=MPI_UNDEFINED; - else color=0; - ISSM_MPI_Comm_split(ISSM_MPI_COMM_WORLD,color, my_rank, &statcomm); - - - iomodel->FindConstant(&numstatistics,"md.qmu.statistics.numstatistics"); - for (int i=1;i<=numstatistics;i++){ - - char* directory=NULL; - char* model=NULL; - int nsamples; - _printf0_("Dealing with qmu statistical computation #" << i << "\n"); - - sprintf(string,"md.qmu.statistics.method(%i).name",i); - iomodel->FindConstant(&name,string); - - sprintf(string,"md.qmu.statistics.method(%i).fields",i); - iomodel->FindConstant(&fields,&nfields,string); - parameters->AddObject(new StringArrayParam(FieldsEnum,fields,nfields)); - - sprintf(string,"md.qmu.statistics.method(%i).steps",i); - iomodel->FetchData(&steps,&dummy,&nsteps,string); - parameters->AddObject(new IntVecParam(StepsEnum,steps,nsteps)); - - if (strcmp(name,"Histogram")==0){ - /*fetch nbins: */ - sprintf(string,"md.qmu.statistics.method(%i).nbins",i); - iomodel->FindConstant(&nbins,string); - parameters->AddObject(new IntParam(NbinsEnum,nbins)); - ComputeHistogram(parameters,results,color,statcomm); - } - else if (strcmp(name,"SampleSeries")==0){ - /*fetch indices: */ - sprintf(string,"md.qmu.statistics.method(%i).indices",i); - iomodel->FetchData(&indices,&dummy,&nindices,string); - parameters->AddObject(new IntVecParam(IndicesEnum,indices,nindices)); - - ComputeSampleSeries(parameters,results,color,statcomm); - } - else if (strcmp(name,"MeanVariance")==0){ - ComputeMeanVariance(parameters,results,color,statcomm); - } - else _error_(" error creating qmu statistics methods parameters: unsupported method " << name); - } - } - //close model file: - fclose(fid); - - /*output results:*/ - OutputStatistics(parameters,results,color,statcomm); - - /*all meet here: */ - ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); _printf0_("Output file.\n"); - - /*Delete ressources:*/ - delete parameters; - delete results; - - - -} /*}}}*/ Index: /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp (revision 25614) +++ /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp (revision 25615) @@ -1,138 +1,21 @@ /*!\file: issm_post.cpp - * \brief: ISSM POST processing main program + * \brief: ISSM DAKOTA post-processing of statistics */ -/*includes and prototypes:*/ + #include "./issm.h" +#include int main(int argc,char **argv){ /*{{{*/ - - char* method=NULL; - int nfields; - char* string=NULL; - char** fields=NULL; - char* field=NULL; - int offset; - char* stepstring=NULL; - int step1,step2; - char* pattern=NULL; - int nsteps; - int* steps=NULL; - int nindices; - int* indices=NULL; - int nbins; - - /*datasets*/ - Parameters *parameters = NULL; - Results *results = NULL; + char* dakota_input_file=NULL; + char* dakota_output_file = NULL; + char* dakota_error_file = NULL; + bool statistics=false; - /*Initialize environment (MPI, PETSC, MUMPS, etc ...)*/ - ISSM_MPI_Comm comm=EnvironmentInit(argc,argv); + /*Initialize MPI: */ + ISSM_MPI_Init(&argc, &argv); // initialize MPI - /*First things first, store the communicator, and set it as a global variable: */ - IssmComm::SetComm(comm); - - /*Initialize and retrieve parameters:{{{*/ - parameters = new Parameters(); - results = new Results(); - - /*Method and Solution:*/ - method=argv[1]; - parameters->AddObject(new IntParam(SolutionTypeEnum,StatisticsSolutionEnum)); - parameters->AddObject(new StringParam(AnalysisTypeEnum,method)); - - /*Directory:*/ - parameters->AddObject(new StringParam(DirectoryNameEnum,argv[2])); - - /*Model name:*/ - parameters->AddObject(new StringParam(InputFileNameEnum,argv[3])); - - /*Number of samples:*/ - parameters->AddObject(new IntParam(QmuNsampleEnum,atoi(argv[4]))); - - /*Retrieve fields:*/ - nfields=atoi(argv[5]); - fields=xNew(nfields); - for(int i=0;i(strlen(string)+1); - xMemCpy(field,string,(strlen(string)+1)); - fields[i]=field; - } - parameters->AddObject(new StringArrayParam(FieldsEnum,fields,nfields)); - offset=6+nfields; - - /*free some memory: */ - for(int i=0;i(field); - } - xDelete(fields); - - /*Retrieve time steps:*/ - stepstring=argv[offset]; - pattern=strstr(stepstring,":"); - if (pattern==NULL){ - step1=atoi(stepstring); - step2=step1; - } - else{ - step2=atoi(pattern+1); - stepstring[pattern-stepstring]='\0'; - step1=atoi(stepstring); - } - nsteps=step2-step1+1; - steps=xNew(nsteps); - for (int i=step1;i<=step2;i++)steps[i-step1]=i; - parameters->AddObject(new IntVecParam(StepsEnum,steps,nsteps)); - offset++; - - /*free some memory:*/ - xDelete(steps); - - /*}}}*/ - - /*Key off method:*/ - if(strcmp(method,"MeanVariance")==0){ - - ComputeMeanVariance(parameters,results,0,ISSM_MPI_COMM_WORLD); - } - else if(strcmp(method,"Histogram")==0){ - - /*Retrieve the size of the histogram (number of bins):*/ - nbins=atoi(argv[offset]); offset++; - parameters->AddObject(new IntParam(NbinsEnum,nbins)); - - ComputeHistogram(parameters,results,0,ISSM_MPI_COMM_WORLD); - - } - else if(strcmp(method,"SampleSeries")==0){ - - /*Retrieve the vertex indices where we'll retrieve our sample time series:*/ - nindices=atoi(argv[offset]); offset++; - indices=xNew(nindices); - for (int i=0;iAddObject(new IntVecParam(IndicesEnum,indices,nindices)); - - /*free some memory:*/ - xDelete(indices); - - ComputeSampleSeries(parameters,results,0,ISSM_MPI_COMM_WORLD); - } - - else _error_("unknown method: " << method << "\n"); - - /*output results:*/ - ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); _printf0_("Output file.\n"); - OutputStatistics(parameters,results,0,ISSM_MPI_COMM_WORLD); - - /*Delete ressources:*/ - delete parameters; - delete results; - - /*Finalize ISSM:*/ - ISSM_MPI_Finalize(); + /*Run statistics:*/ + DakotaStatistics(argc,argv); /*Return unix success: */ Index: /issm/branches/trunk-larour-SLPS2020/src/c/modules/QmuStatisticsx/QmuStatisticsx.cpp =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/modules/QmuStatisticsx/QmuStatisticsx.cpp (revision 25614) +++ /issm/branches/trunk-larour-SLPS2020/src/c/modules/QmuStatisticsx/QmuStatisticsx.cpp (revision 25615) @@ -2,4 +2,5 @@ */ /*includes and prototypes:*/ +#include #include "./QmuStatisticsx.h" #include "../OutputResultsx/OutputResultsx.h" @@ -433,4 +434,6 @@ IssmDouble mi=*minxs[counter]; int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index--; + if(ma==mi)index=0; + //_printf_( index << "|" << scalar << "|" << mi << "|" << ma << "|" << nbins << "\n"); localhistogram[index]++; histogram[counter]=localhistogram; @@ -443,4 +446,5 @@ IssmDouble scalar=doublemat[k]; int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index--; + if (mi[k]==ma[k])index=0; _assert_(scalar<=ma[k]); _assert_(scalar>=mi[k]); _assert_(index((strlen(argv[2])+strlen(argv[3])+strlen(".bin")+2)); + sprintf(input_file,"%s/%s%s",argv[2],argv[3],".bin"); + fid=fopen(input_file,"rb"); + if (fid==NULL) Cerr << "dirstructure error message: could not open model " << input_file << " to retrieve qmu statistics parameters" << std::endl; + + //initialize IoModel, but light version, we just need it to fetch one constant: + iomodel=new IoModel(); + iomodel->fid=fid; + iomodel->FetchConstants(); + + //early return if statistics not requested: + iomodel->FindConstant(&statistics,"md.qmu.statistics"); + if(!statistics){ + delete iomodel; + fclose(fid); + return false; //important return value! + } + + iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories"); + + /*Ok, we have everything we need to create the directory structure:*/ + if(IssmComm::GetRank()==0){ + for (int i=0;i((strlen(argv[2])+strlen(argv[3])+strlen(".bin")+2)); + sprintf(input_file,"%s/%s%s",argv[2],argv[3],".bin"); + fid=fopen(input_file,"rb"); + if (fid==NULL) Cerr << "issm_dakota_statistics error message: could not open model " << input_file << " to retrieve qmu statistics parameters" << std::endl; + + //initialize IoModel, but light version, we'll need it to fetch constants: + iomodel=new IoModel(); + iomodel->fid=fid; + iomodel->FetchConstants(); + + //early return if statistics not requested: + iomodel->FindConstant(&statistics,"md.qmu.statistics"); + if(!statistics){ + delete iomodel; + fclose(fid); + return 0; + } + + //create parameters datasets with al the qmu statistics settings we need: + if(statistics){ + + /*Initialize parameters and results:*/ + results = new Results(); + parameters=new Parameters(); + + //solution type: + parameters->AddObject(new IntParam(SolutionTypeEnum,StatisticsSolutionEnum)); + + //root directory + directory=xNew(strlen(argv[2])+1); + xMemCpy(directory,argv[2],strlen(argv[2])+1); + parameters->AddObject(new StringParam(DirectoryNameEnum,directory)); + + //model name + model=xNew(strlen(argv[3])+1); + xMemCpy(model,argv[3],strlen(argv[3])+1); + parameters->AddObject(new StringParam(InputFileNameEnum,model)); + + //nsamples + iomodel->FindConstant(&nsamples,"md.qmu.method.params.samples"); + parameters->AddObject(new IntParam(QmuNsampleEnum,nsamples)); + + //ndirectories + iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories"); + parameters->AddObject(new IntParam(QmuNdirectoriesEnum,numdirectories)); + + //nfiles per directory + iomodel->FindConstant(&nfilesperdirectory,"md.qmu.statistics.nfiles_per_directory"); + parameters->AddObject(new IntParam(QmuNfilesPerDirectoryEnum,nfilesperdirectory)); + + //At this point, we don't want to go forward any longer, we want to create an MPI + //communicator on which to carry out the computations: + if ((my_rank+1)*nfilesperdirectory>nsamples)color=MPI_UNDEFINED; + else color=0; + ISSM_MPI_Comm_split(ISSM_MPI_COMM_WORLD,color, my_rank, &statcomm); + + + iomodel->FindConstant(&numstatistics,"md.qmu.statistics.numstatistics"); + for (int i=1;i<=numstatistics;i++){ + + char* directory=NULL; + char* model=NULL; + int nsamples; + _printf0_("Dealing with qmu statistical computation #" << i << "\n"); + + sprintf(string,"md.qmu.statistics.method(%i).name",i); + iomodel->FindConstant(&name,string); + + sprintf(string,"md.qmu.statistics.method(%i).fields",i); + iomodel->FindConstant(&fields,&nfields,string); + parameters->AddObject(new StringArrayParam(FieldsEnum,fields,nfields)); + + sprintf(string,"md.qmu.statistics.method(%i).steps",i); + iomodel->FetchData(&steps,&dummy,&nsteps,string); + parameters->AddObject(new IntVecParam(StepsEnum,steps,nsteps)); + + if (strcmp(name,"Histogram")==0){ + /*fetch nbins: */ + sprintf(string,"md.qmu.statistics.method(%i).nbins",i); + iomodel->FindConstant(&nbins,string); + parameters->AddObject(new IntParam(NbinsEnum,nbins)); + ComputeHistogram(parameters,results,color,statcomm); + } + else if (strcmp(name,"SampleSeries")==0){ + /*fetch indices: */ + sprintf(string,"md.qmu.statistics.method(%i).indices",i); + iomodel->FetchData(&indices,&dummy,&nindices,string); + parameters->AddObject(new IntVecParam(IndicesEnum,indices,nindices)); + + ComputeSampleSeries(parameters,results,color,statcomm); + } + else if (strcmp(name,"MeanVariance")==0){ + ComputeMeanVariance(parameters,results,color,statcomm); + } + else _error_(" error creating qmu statistics methods parameters: unsupported method " << name); + } + } + //close model file: + fclose(fid); + + /*output results:*/ + OutputStatistics(parameters,results,color,statcomm); + + /*all meet here: */ + ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); _printf0_("Output file.\n"); + + /*Delete ressources:*/ + delete parameters; + delete results; + + + +} /*}}}*/ Index: /issm/branches/trunk-larour-SLPS2020/src/c/modules/QmuStatisticsx/QmuStatisticsx.h =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/modules/QmuStatisticsx/QmuStatisticsx.h (revision 25614) +++ /issm/branches/trunk-larour-SLPS2020/src/c/modules/QmuStatisticsx/QmuStatisticsx.h (revision 25615) @@ -12,4 +12,6 @@ int ComputeHistogram(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm); int readdata(IssmDouble** pdoublemat, int* pdoublematsize, IssmDouble* pdouble, FILE* fid,char* field,int step); +bool DakotaDirStructure(int argc,char** argv); +int DakotaStatistics(int argc,char** argv); /* local prototypes: */