Index: /issm/branches/trunk-larour-SLPS2020/src/c/Makefile.am =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/Makefile.am (revision 25595) +++ /issm/branches/trunk-larour-SLPS2020/src/c/Makefile.am (revision 25596) @@ -635,5 +635,6 @@ ./modules/NodeConnectivityx/NodeConnectivityx.cpp \ ./modules/ElementConnectivityx/ElementConnectivityx.cpp \ - ./modules/PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.cpp + ./modules/PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.cpp \ + ./modules/QmuStatisticsx/QmuStatisticsx.cpp if CHACO Index: /issm/branches/trunk-larour-SLPS2020/src/c/classes/FemModel.cpp =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/classes/FemModel.cpp (revision 25595) +++ /issm/branches/trunk-larour-SLPS2020/src/c/classes/FemModel.cpp (revision 25596) @@ -70,4 +70,5 @@ char *restartfilename = NULL; char *rootpath = NULL; + char *modelname = NULL; /*First things first, store the communicator, and set it as a global variable: */ @@ -85,9 +86,9 @@ /*From command line arguments, retrieve different filenames needed to create the FemModel: */ - ProcessArguments(&solution_type,&binfilename,&outbinfilename,&petscfilename,&lockfilename,&restartfilename,&rootpath,argc,argv); + ProcessArguments(&solution_type,&binfilename,&outbinfilename,&petscfilename,&lockfilename,&restartfilename,&rootpath,&modelname,argc,argv); /*Create femmodel from input files: */ profiler->Start(MPROCESSOR); - this->InitFromFiles(rootpath,binfilename,outbinfilename,petscfilename,lockfilename,restartfilename, solution_type,trace,NULL); + this->InitFromFiles(rootpath,binfilename,outbinfilename,petscfilename,lockfilename,restartfilename, modelname, solution_type,trace,NULL); profiler->Stop(MPROCESSOR); @@ -148,5 +149,5 @@ } /*}}}*/ -FemModel::FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X){ /*{{{*/ +FemModel::FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, char* modelname, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X){ /*{{{*/ bool traceon=true; @@ -159,5 +160,5 @@ /*Create femmodel from input files, with trace activated: */ profiler->Start(MPROCESSOR); - this->InitFromFiles(rootpath,inputfilename,outputfilename,toolkitsfilename,lockfilename,restartfilename, solution_type,traceon,X); + this->InitFromFiles(rootpath,inputfilename,outputfilename,toolkitsfilename,lockfilename,restartfilename, modelname,solution_type,traceon,X); profiler->Stop(MPROCESSOR); @@ -394,5 +395,5 @@ } /*}}}*/ -void FemModel::InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, const int in_solution_type,bool trace,IssmPDouble* X){/*{{{*/ +void FemModel::InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, char* modelname, const int in_solution_type,bool trace,IssmPDouble* X){/*{{{*/ /*intermediary*/ @@ -418,4 +419,5 @@ /*Now save all of these file names into parameters, you never know when you might need them: */ this->parameters->AddObject(new StringParam(ToolkitsFileNameEnum,toolkitsfilename)); + this->parameters->AddObject(new StringParam(ModelnameEnum,modelname)); this->parameters->AddObject(new StringParam(RootPathEnum,rootpath)); this->parameters->AddObject(new StringParam(InputFileNameEnum,inputfilename)); Index: /issm/branches/trunk-larour-SLPS2020/src/c/classes/FemModel.h =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/classes/FemModel.h (revision 25595) +++ /issm/branches/trunk-larour-SLPS2020/src/c/classes/FemModel.h (revision 25596) @@ -69,5 +69,5 @@ FemModel(void); FemModel(int argc,char** argv,ISSM_MPI_Comm comm_init,bool trace=false); - FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X); + FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, char* modelname, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X); ~FemModel(); @@ -79,5 +79,5 @@ void Echo(); int GetElementsWidth(){return 3;};//just tria elements in this first version - void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, char* restartfilename, const int solution_type,bool trace,IssmPDouble* X=NULL); + void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, char* restartfilename, char* modelname, const int solution_type,bool trace,IssmPDouble* X=NULL); void InitFromFids(char* rootpath, FILE* IOMODEL, FILE* toolkitsoptionsfid, int in_solution_type, bool trace, IssmPDouble* X=NULL); void Marshall(char** pmarshalled_data, int* pmarshalled_data_size, int marshall_direction); Index: /issm/branches/trunk-larour-SLPS2020/src/c/cores/ProcessArguments.cpp =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/cores/ProcessArguments.cpp (revision 25595) +++ /issm/branches/trunk-larour-SLPS2020/src/c/cores/ProcessArguments.cpp (revision 25596) @@ -8,5 +8,5 @@ #include "../shared/shared.h" -void ProcessArguments(int* solution_type,char** pbinfilename,char** poutbinfilename,char** ptoolkitsfilename,char** plockfilename,char** prestartfilename, char** prootpath, int argc,char **argv){ +void ProcessArguments(int* solution_type,char** pbinfilename,char** poutbinfilename,char** ptoolkitsfilename,char** plockfilename,char** prestartfilename, char** prootpath, char** pmodelname, int argc,char **argv){ /*Check input arguments*/ @@ -18,5 +18,6 @@ *solution_type = StringToEnumx(argv[1]); char* rootpatharg = argv[2]; - char* modelname = argv[3]; + char* modelname = xNew(strlen(argv[3])+1); + xMemCpy(modelname,argv[3],strlen(argv[3])+1); /*Recover myrank and length of string "my_rank" */ @@ -42,4 +43,5 @@ *prestartfilename=restartfilename; *prootpath=rootpath; + *pmodelname=modelname; } Index: /issm/branches/trunk-larour-SLPS2020/src/c/cores/cores.h =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/cores/cores.h (revision 25595) +++ /issm/branches/trunk-larour-SLPS2020/src/c/cores/cores.h (revision 25596) @@ -72,5 +72,5 @@ //diverse -void ProcessArguments(int* solution,char** pbinname,char** poutbinname,char** ptoolkitsname,char** plockname,char** prestartname, char** prootpath,int argc,char **argv); +void ProcessArguments(int* solution,char** pbinname,char** poutbinname,char** ptoolkitsname,char** plockname,char** prestartname, char** prootpath,char** pmodelname, int argc,char **argv); void WriteLockFile(char* filename); void ResetBoundaryConditions(FemModel* femmodel, int analysis_type); Index: /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_dakota.cpp =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_dakota.cpp (revision 25595) +++ /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_dakota.cpp (revision 25596) @@ -4,4 +4,5 @@ #include "./issm.h" +#include /*Dakota includes: */ @@ -14,5 +15,9 @@ #endif -int main(int argc,char **argv){ +/*prototypes:*/ +int dirstructure(int argc,char** argv); +int issm_dakota_statistics(int argc,char** argv); + +int main(int argc,char **argv){ /*{{{*/ #if defined(_HAVE_DAKOTA_) && _DAKOTA_MAJOR_ >= 6 @@ -38,4 +43,7 @@ dakota_error_file=xNew((strlen(argv[2])+strlen(argv[3])+strlen(".qmu.err")+2)); sprintf(dakota_error_file,"%s/%s%s",argv[2],argv[3],".qmu.err"); + + /*Create directory structure for model outputs:*/ + dirstructure(argc,argv); /* Parse input and construct Dakota LibraryEnvironment, performing input data checks*/ @@ -83,4 +91,8 @@ env.execute(); + /* Run statistics if requested:*/ + issm_dakota_statistics(argc,argv); + + /*free allocations:*/ xDelete(dakota_input_file); xDelete(dakota_output_file); @@ -94,3 +106,203 @@ #endif -} +} /*}}}*/ +int 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 0; + } + + 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(); + + //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,&nsteps,&dummy,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,&nindices,&dummy,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:*/ + ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); _printf0_("Output file.\n"); + OutputStatistics(parameters,results); + + /*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 25595) +++ /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp (revision 25596) @@ -4,10 +4,4 @@ /*includes and prototypes:*/ #include "./issm.h" -int ComputeMeanVariance(Parameters* parameters,Results* results); -int ComputeSampleSeries(Parameters* parameters,Results* results); -int readdata(IssmDouble** pdoublemat, int* pdoublematsize, - IssmDouble* pdouble, FILE* fid,char* field,int step); -int OutputStatistics(Parameters* parameters,Results* results); -int ComputeHistogram(Parameters* parameters,Results* results); int main(int argc,char **argv){ /*{{{*/ @@ -146,1159 +140,2 @@ } /*}}}*/ -int ComputeMeanVariance(Parameters* parameters,Results* results){ /*{{{*/ - - int nsamples; - char* directory=NULL; - char* model=NULL; - char** fields=NULL; - int* steps=NULL; - int nsteps; - int nfields; - int range,lower_row,upper_row; - - /*intermediary:*/ - IssmDouble* doublemat=NULL; - int doublematsize; - IssmDouble scalar; - - /*computation of average and variance itself:*/ - IssmDouble* x = NULL; - IssmDouble* x2 = NULL; - IssmDouble** xs = NULL; - IssmDouble** xs2 = NULL; - int* xtype=NULL; - int* xsize=NULL; - - IssmDouble** meanx=NULL; - IssmDouble** meanx2=NULL; - int* meantype=NULL; - int* meansize=NULL; - - /*Retrieve parameters:*/ - parameters->FindParam(&nsamples,QmuNsampleEnum); - parameters->FindParam(&directory,DirectoryNameEnum); - parameters->FindParam(&model,InputFileNameEnum); - parameters->FindParam(&fields,&nfields,FieldsEnum); - parameters->FindParam(&steps,&nsteps,StepsEnum); - - /*Get rank:*/ - int my_rank=IssmComm::GetRank(); - - /*Open files and read them complelety, in a distributed way:*/ - range=DetermineLocalSize(nsamples,IssmComm::GetComm()); - GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm()); - - /*Initialize arrays:*/ - xs=xNew(nfields*nsteps); - xs2=xNew(nfields*nsteps); - xtype=xNew(nfields*nsteps); - xsize=xNew(nfields*nsteps); - - meantype=xNew(nfields); - meansize=xNew(nfields); - meanx=xNew(nfields); - meanx2=xNew(nfields); - - /*Start opening files:*/ - for (int i=(lower_row+1);i<=upper_row;i++){ - _printf0_("reading file #: " << i << "\n"); - char file[1000]; - long int length; - char* buffer=NULL; - - /*string:*/ - sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i); - - /*open file: */ - _printf0_(" opening file: " << file << "\n"); - FILE* fid=fopen(file,"rb"); - if(fid==NULL) _error_(" could not open file: " << file << "\n"); - - /*figure out size of file, and read the whole thing:*/ - _printf0_(" reading file:\n"); - fseek (fid, 0, SEEK_END); - length = ftell (fid); - fseek (fid, 0, SEEK_SET); - buffer = xNew(length); - fread (buffer, sizeof(char), length, fid); - - /*close file:*/ - fclose (fid); - - /*create a memory stream with this buffer:*/ - _printf0_(" processing file:\n"); - fid=fmemopen(buffer, length, "rb"); - - /*start reading data from the buffer directly:*/ - for (int f=0;f(1); - xs2[counter]=xNew(1); - *xs[counter]=scalar; - *xs2[counter]=pow(scalar,2.0); - xsize[counter]=1; - } - else if (xtype[counter]==3){ - IssmDouble* doublemat2=xNew(doublematsize); - for(int k=0;k(1); - meanx2[f]=xNewZeroInit(1); - meansize[f]=1; - } - else{ - meanx[f]=xNewZeroInit(doublematsize); - meanx2[f]=xNewZeroInit(doublematsize); - meansize[f]=doublematsize; - } - } - fseek(fid,0,SEEK_SET); - if(meantype[f]==1){ - IssmDouble sc=0; - IssmDouble sc2=0; - for(int j=0;j(doublematsize); - IssmDouble* timemean2=xNewZeroInit(doublematsize); - - for(int j=0;j(buffer); - } - ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); - _printf0_("Done reading files, now computing mean and variance.\n"); - - /*We have agregated x and x^2 across the cluster, now gather across the cluster onto - *cpu0 and then compute statistics:*/ - for (int f=0;fAddResult(new GenericExternalResult(results->Size()+1,fieldname,mean,j+1,0)); - sprintf(fieldname,"%s%s",fields[f],"Stddev"); - results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,stddev,j+1,0)); - } - - } - } /*}}}*/ - else{ /*deal with arrays:{{{*/ - - int size=xsize[counter0]; - - IssmDouble* mean=xNew(size); - IssmDouble* stddev=xNew(size); - - for (int j=0;j(size); - IssmDouble* sumx2=xNew(size); - - ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm()); - ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm()); - - /*Build average and standard deviation. For standard deviation, use the - *following formula: sigma^2=E(x^2)-mu^2:*/ - for (int k=0;kAddResult(new GenericExternalResult(results->Size()+1,fieldname,mean,size,1,j+1,0)); - sprintf(fieldname,"%s%s",fields[f],"Stddev"); - results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,stddev,size,1,j+1,0)); - } - } - } /*}}}*/ - } - /*Do the same but for the time mean:*/ - for (int f=0;fAddResult(new GenericExternalResult(results->Size()+1,fieldname,mean,1,0)); - sprintf(fieldname,"%s%s",fields[f],"TimeStddev"); - results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,stddev,1,0)); - } - } /*}}}*/ - else{ /*deal with arrays:{{{*/ - - int size=meansize[f]; - IssmDouble* mean=xNew(size); - IssmDouble* stddev=xNew(size); - - /*we are broadcasting double arrays:*/ - x=meanx[f]; - x2=meanx2[f]; - - IssmDouble* sumx=xNew(size); - IssmDouble* sumx2=xNew(size); - - ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm()); - ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm()); - - /*Build average and standard deviation. For standard deviation, use the - *following formula: sigma^2=E(x^2)-mu^2:*/ - for (int k=0;kAddResult(new GenericExternalResult(results->Size()+1,fieldname,mean,size,1,1,0)); - sprintf(fieldname,"%s%s",fields[f],"TimeStddev"); - results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,stddev,size,1,1,0)); - } - } /*}}}*/ - } - - -} /*}}}*/ -int ComputeSampleSeries(Parameters* parameters,Results* results){ /*{{{*/ - - int nsamples; - char* directory=NULL; - char* model=NULL; - char** fields=NULL; - int* steps=NULL; - int nsteps; - int nfields; - int range,lower_row,upper_row; - int* indices=NULL; - int nindices; - - /*intermediary:*/ - IssmDouble* doublemat=NULL; - int doublematsize; - IssmDouble scalar; - - /*computation of average and variance itself:*/ - IssmDouble* x = NULL; - IssmDouble* allx=NULL; - IssmDouble** xs = NULL; - int* xtype=NULL; - int* xsize=NULL; - - /*Retrieve parameters:*/ - parameters->FindParam(&nsamples,QmuNsampleEnum); - parameters->FindParam(&directory,DirectoryNameEnum); - parameters->FindParam(&model,InputFileNameEnum); - parameters->FindParam(&fields,&nfields,FieldsEnum); - parameters->FindParam(&steps,&nsteps,StepsEnum); - parameters->FindParam(&indices,&nindices,IndicesEnum); - - /*Get rank:*/ - int my_rank=IssmComm::GetRank(); - - /*Open files and read them complelety, in a distributed way:*/ - range=DetermineLocalSize(nsamples,IssmComm::GetComm()); - GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm()); - - /*Initialize arrays:*/ - xs=xNew(nfields*nsteps); - xtype=xNew(nfields*nsteps); - xsize=xNew(nfields*nsteps); - - /*Start opening files:*/ - for (int i=(lower_row+1);i<=upper_row;i++){ - _printf0_("reading file #: " << i << "\n"); - char file[1000]; - long int length; - char* buffer=NULL; - - /*string:*/ - sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i); - - /*open file: */ - _printf0_(" opening file:\n"); - FILE* fid=fopen(file,"rb"); - - /*figure out size of file, and read the whole thing:*/ - _printf0_(" reading file:\n"); - fseek (fid, 0, SEEK_END); - length = ftell (fid); - fseek (fid, 0, SEEK_SET); - buffer = xNew(length); - fread (buffer, sizeof(char), length, fid); - - /*close file:*/ - fclose (fid); - - /*create a memory stream with this buffer:*/ - _printf0_(" processing file:\n"); - fid=fmemopen(buffer, length, "rb"); - - /*start reading data from the buffer directly:*/ - for (int f=0;f(range); - x[0]=scalar; - xs[counter]=x; - xsize[counter]=range; - } - else if (xtype[counter]==3){ - x=xNew(nindices*range); - for(int k=0;k(buffer); - } - ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); - _printf0_("Done reading files, now assembling time series.\n"); - - for (int f=0;f(nsamples); - MPI_Gather(x, range, ISSM_MPI_PDOUBLE,allx, range, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm()); - /*add to results:*/ - if(my_rank==0){ - char fieldname[1000]; - - sprintf(fieldname,"%s%s",fields[f],"Samples"); - results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,allx,nsamples,1,j+1,0)); - } - } - else{ - /*we are broadcasting double arrays:*/ - x=xs[counter]; - allx=xNew(nsamples*nindices); - - MPI_Gather(x, range*nindices, ISSM_MPI_PDOUBLE,allx, range*nindices, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm()); - - /*add to results:*/ - if(my_rank==0){ - char fieldname[1000]; - sprintf(fieldname,"%s%s",fields[f],"Samples"); - results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,allx,nsamples,nindices,j+1,0)); - } - } - } - } - - -} /*}}}*/ -int ComputeHistogram(Parameters* parameters,Results* results){ /*{{{*/ - - int nsamples; - char* directory=NULL; - char* model=NULL; - char** fields=NULL; - int* steps=NULL; - int nsteps; - int nfields; - int nbins; - int range,lower_row,upper_row; - - /*intermediary:*/ - IssmDouble* doublemat=NULL; - int doublematsize; - IssmDouble scalar; - - /*computation of average and variance itself:*/ - IssmDouble** maxxs = NULL; - IssmDouble** minxs = NULL; - int* xtype=NULL; - int* xsize=NULL; - - IssmDouble** maxmeans=NULL; - IssmDouble** minmeans=NULL; - int* meanxtype=NULL; - int* meanxsize=NULL; - - /*Retrieve parameters:*/ - parameters->FindParam(&nsamples,QmuNsampleEnum); - parameters->FindParam(&directory,DirectoryNameEnum); - parameters->FindParam(&model,InputFileNameEnum); - parameters->FindParam(&fields,&nfields,FieldsEnum); - parameters->FindParam(&steps,&nsteps,StepsEnum); - parameters->FindParam(&nbins,NbinsEnum); - - /*Get rank:*/ - int my_rank=IssmComm::GetRank(); - - /*Open files and read them complelety, in a distributed way:*/ - range=DetermineLocalSize(nsamples,IssmComm::GetComm()); - GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm()); - - /*Initialize arrays:*/ - - maxmeans=xNew(nfields); - minmeans=xNew(nfields); - meanxtype=xNew(nfields); - meanxsize=xNew(nfields); - - maxxs=xNew(nfields*nsteps); - minxs=xNew(nfields*nsteps); - xtype=xNew(nfields*nsteps); - xsize=xNew(nfields*nsteps); - - /*Start opening files:*/ - for (int i=(lower_row+1);i<=upper_row;i++){ - _printf0_("reading file #: " << i << "\n"); - char file[1000]; - long int length; - char* buffer=NULL; - - /*string:*/ - sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i); - - /*open file: */ - _printf0_(" opening file: " << file << "\n"); - FILE* fid=fopen(file,"rb"); - if(fid==NULL)_error_("cound not open file: " << file << "\n"); - - /*figure out size of file, and read the whole thing:*/ - _printf0_(" reading file:\n"); - fseek (fid, 0, SEEK_END); - length = ftell (fid); - fseek (fid, 0, SEEK_SET); - buffer = xNew(length); - fread (buffer, sizeof(char), length, fid); - - /*close file:*/ - fclose (fid); - - /*create a memory stream with this buffer:*/ - _printf0_(" processing file:\n"); - fid=fmemopen(buffer, length, "rb"); - - /*start reading data from the buffer directly:*/ - for (int f=0;f(1); - minxs[counter]=xNew(1); - *maxxs[counter]=scalar; - *minxs[counter]=scalar; - xsize[counter]=1; - } - else if (xtype[counter]==3){ - maxxs[counter]=xNew(doublematsize); - xMemCpy(maxxs[counter],doublemat,doublematsize); - minxs[counter]=xNew(doublematsize); - xMemCpy(minxs[counter],doublemat,doublematsize); - xsize[counter]=doublematsize; - xDelete(doublemat); - } - else _error_("cannot carry out statistics on type " << xtype[counter]); - } - else{ - if(xtype[counter]==1){ - *maxxs[counter]=max(*maxxs[counter],scalar); - *minxs[counter]=min(*minxs[counter],scalar); - } - else if (xtype[counter]==3){ - IssmDouble* newmax=maxxs[counter]; - IssmDouble* newmin=minxs[counter]; - for(int k=0;knewmax[k])newmax[k]=doublemat[k]; - if(doublemat[k](doublemat); - } - else _error_("cannot carry out statistics on type " << xtype[counter]); - } - } - } - _printf0_(" average in time:\n"); - - /*Deal with average in time: */ - for (int f=0;f(1); - minmeans[f]=xNewZeroInit(1); - *maxmeans[f]=timemean; - *minmeans[f]=timemean; - } - else{ - *maxmeans[f]=max(*maxmeans[f],timemean); - *minmeans[f]=min(*minmeans[f],timemean); - } - } - else{ - meanxsize[f]=doublematsize; - fseek(fid,0,SEEK_SET); - IssmDouble* timemean=xNewZeroInit(doublematsize); - for (int j=0;j(doublemat); - } - - if(i==(lower_row+1)){ - maxmeans[f]=xNew(doublematsize); - xMemCpy(maxmeans[f],timemean,doublematsize); - minmeans[f]=xNew(doublematsize); - xMemCpy(minmeans[f],timemean,doublematsize); - } - else{ - IssmDouble* maxx=maxmeans[f]; - IssmDouble* minx=minmeans[f]; - - for(int k=0;k(buffer); - } - ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); - _printf0_("Done reading files, now computing min and max.\n"); - - /*We have agregated minx and max across the cluster, now gather across the cluster onto - *cpu0 and then compute statistics:*/ - for (int f=0;f(size); - IssmDouble* allmin=xNew(size); - - ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm()); - ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm()); - - /*Store broadcasted value for later computation of histograms:*/ - maxxs[counter]=allmax; - minxs[counter]=allmin; - } - } /*}}}*/ - } - - /*Now do the same for the time mean fields:*/ - for (int f=0;f(size); - IssmDouble* allmin=xNew(size); - - ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm()); - ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm()); - - - /*Store for later use in histogram computation:*/ - maxmeans[f]=allmax; - minmeans[f]=allmin; - - } /*}}}*/ - } - - /*Now that we have the min and max, we can start binning. First allocate - * histograms, then start filling them:*/ - IssmDouble** histogram=xNew(nfields*nsteps); - IssmDouble** timehistogram=xNew(nfields); - - _printf0_("Start reading files again, this time binning values in the histogram:\n"); - /*Start opening files:*/ - for (int i=(lower_row+1);i<=upper_row;i++){ - _printf0_("reading file #: " << i << "\n"); - char file[1000]; - long int length; - char* buffer=NULL; - - /*string:*/ - sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i); - - /*open file: */ - _printf0_(" opening file:\n"); - FILE* fid=fopen(file,"rb"); - if(fid==NULL)_error_("cound not open file: " << file << "\n"); - - /*figure out size of file, and read the whole thing:*/ - _printf0_(" reading file:\n"); - fseek (fid, 0, SEEK_END); - length = ftell (fid); - fseek (fid, 0, SEEK_SET); - buffer = xNew(length); - fread (buffer, sizeof(char), length, fid); - - /*close file:*/ - fclose (fid); - - /*create a memory stream with this buffer:*/ - _printf0_(" processing file:\n"); - fid=fmemopen(buffer, length, "rb"); - - /*start reading data from the buffer directly:*/ - for (int f=0;f(nbins); - IssmDouble ma=*maxxs[counter]; - IssmDouble mi=*minxs[counter]; - int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index--; - localhistogram[index]++; - histogram[counter]=localhistogram; - } - else if (xtype[counter]==3){ - IssmDouble* localhistogram=xNewZeroInit(doublematsize*nbins); - IssmDouble* ma=maxxs[counter]; - IssmDouble* mi=minxs[counter]; - for (int k=0;k=mi[k]); _assert_(index(doublemat); - } - else _error_("cannot carry out statistics on type " << xtype[counter]); - } - else{ - if(xtype[counter]==1){ - IssmDouble* localhistogram=histogram[counter]; - IssmDouble ma=*maxxs[counter]; - IssmDouble mi=*minxs[counter]; - int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1; - localhistogram[index]++; - } - else if (xtype[counter]==3){ - IssmDouble* localhistogram=histogram[counter]; - IssmDouble* ma=maxxs[counter]; - IssmDouble* mi=minxs[counter]; - for (int k=0;k(doublemat); - } - else _error_("cannot carry out statistics on type " << xtype[counter]); - } - } - } - _printf0_(" average in time:\n"); - - /*Deal with average in time: */ - for (int f=0;f(nbins); - IssmDouble ma=*maxmeans[f]; - IssmDouble mi=*minmeans[f]; - int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1; - localhistogram[index]++; - timehistogram[f]=localhistogram; - } - else{ - IssmDouble* localhistogram=timehistogram[f]; - IssmDouble ma=*maxmeans[f]; - IssmDouble mi=*minmeans[f]; - int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1; - localhistogram[index]++; - } - } - else{ - fseek(fid,0,SEEK_SET); - IssmDouble* timemean=xNewZeroInit(doublematsize); - for (int j=0;j(doublemat); - } - - if(i==(lower_row+1)){ - IssmDouble* localhistogram=xNewZeroInit(doublematsize*nbins); - IssmDouble* ma=maxmeans[f]; - IssmDouble* mi=minmeans[f]; - - for (int k=0;k(buffer); - } - _printf0_("Start aggregating histogram:\n"); - - /*We have agregated histograms across the cluster, now gather them across the cluster onto - *cpu0: */ - for (int f=0;f(nbins); - - ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); - - /*add to results:*/ - if(my_rank==0){ - char fieldname[1000]; - - sprintf(fieldname,"%s%s",fields[f],"Histogram"); - results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,allhisto,1,nbins,j+1,0)); - - sprintf(fieldname,"%s%s",fields[f],"Max"); - results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,*maxxs[counter],j+1,0)); - sprintf(fieldname,"%s%s",fields[f],"Min"); - results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,*minxs[counter],j+1,0)); - } - } - } /*}}}*/ - else{ /*deal with arrays:{{{*/ - - int size=xsize[counter0]; - for (int j=0;j(size*nbins); - - ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); - xDelete(histo); - - /*add to results:*/ - if(my_rank==0){ - char fieldname[1000]; - - sprintf(fieldname,"%s%s",fields[f],"Histogram"); - results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,allhisto,size,nbins,j+1,0)); - - sprintf(fieldname,"%s%s",fields[f],"Max"); - results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,maxxs[counter],size,1,j+1,0)); - sprintf(fieldname,"%s%s",fields[f],"Min"); - results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,minxs[counter],size,1,j+1,0)); - } - } - } /*}}}*/ - } - _printf0_("Start aggregating time mean histogram:\n"); - - /*Now do the same for the time mean fields:*/ - for (int f=0;f(nbins); - - ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm()); - - /*add to results at time step 1:*/ - if(my_rank==0){ - char fieldname[1000]; - - sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram"); - results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,allhisto,1,nbins,1,0)); - - sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax"); - results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,*maxmeans[f],1,0)); - sprintf(fieldname,"%s%s",fields[f],"TimeMeaMin"); - results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,*minmeans[f],1,0)); - } - } /*}}}*/ - else{ /*deal with arrays:{{{*/ - - int size=meanxsize[f]; - - /*we are broadcasting double arrays:*/ - IssmDouble* histo=timehistogram[f]; - IssmDouble* allhisto=xNewZeroInit(size*nbins); - - ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); - xDelete(histo); - /*add to results at step 1:*/ - if(my_rank==0){ - char fieldname[1000]; - - sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram"); - results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,allhisto,size,nbins,1,0)); - sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax"); - results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,maxmeans[f],size,1,1,0)); - sprintf(fieldname,"%s%s",fields[f],"TimeMeanMin"); - results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,minmeans[f],size,1,1,0)); - } - } /*}}}*/ - } -} /*}}}*/ -int OutputStatistics(Parameters* parameters,Results* results){ /*{{{*/ - - char outputfilename[1000]; - char* directory=NULL; - char* model=NULL; - char* method=NULL; - int nsamples; - int* steps=NULL; - int nsteps; - - FemModel* femmodel=new FemModel(); - - /*Some parameters that will allow us to use the OutputResultsx module:*/ - parameters->AddObject(new BoolParam(QmuIsdakotaEnum,false)); - parameters->AddObject(new BoolParam(SettingsIoGatherEnum,true)); - - parameters->FindParam(&directory,DirectoryNameEnum); - parameters->FindParam(&model,InputFileNameEnum); - parameters->FindParam(&method,AnalysisTypeEnum); - parameters->FindParam(&nsamples,QmuNsampleEnum); - parameters->FindParam(&steps,&nsteps,StepsEnum); - - sprintf(outputfilename,"%s/%s-%s.stats-samp%i-time%i:%i",directory,model,method,nsamples,steps[0],steps[nsteps-1]); - parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename)); - - /*Call OutputResults module:*/ - femmodel->parameters=parameters; - femmodel->results=results; - - OutputResultsx(femmodel); -} /*}}}*/ -int readdata(IssmDouble** pdoublemat, int* pdoublematsize, IssmDouble* pdouble, FILE* fid,char* field,int step){ /*{{{*/ - - int length; - char fieldname[1000]; - int fieldname_size; - IssmDouble rtime; - int rstep; - int M,N; - - //fields that we retrive: - IssmDouble dfield; - char* sfield = NULL; - IssmDouble* dmatfield = NULL; - int* imatfield = NULL; - - //type of the returned field: - int type; - int found=0; - - while(1){ - - size_t ret_code = fread(&fieldname_size, sizeof(int), 1, fid); - if(ret_code != 1) break; //we are done. - - fread(fieldname, sizeof(char), fieldname_size, fid); - //_printf0_("fieldname: " << fieldname << "\n"); - - fread(&rtime, sizeof(IssmDouble), 1, fid); - fread(&rstep, sizeof(int), 1, fid); - - //check on field: - if ((step==rstep) && (strcmp(field,fieldname)==0)){ - - //ok, go read the result really: - fread(&type,sizeof(int),1,fid); - fread(&M,sizeof(int),1,fid); - if (type==1){ - fread(&dfield,sizeof(IssmDouble),1,fid); - } - else if (type==2){ - fread(&M,sizeof(int),1,fid); - sfield=xNew(M); - fread(sfield,sizeof(char),M,fid); - } - else if (type==3){ - fread(&N,sizeof(int),1,fid); - dmatfield=xNew(M*N); - fread(dmatfield,sizeof(IssmDouble),M*N,fid); - } - else if (type==4){ - fread(&N,sizeof(int),1,fid); - imatfield=xNew(M*N); - fread(imatfield,sizeof(int),M*N,fid); - } - else _error_("cannot read data of type " << type << "\n"); - found=1; - break; - } - else{ - //just skim to next results. - fread(&type,sizeof(int),1,fid); - fread(&M,sizeof(int),1,fid); - if (type==1){ - fseek(fid,sizeof(IssmDouble),SEEK_CUR); - } - else if(type==2){ - fseek(fid,M*sizeof(char),SEEK_CUR); - } - else if(type==3){ - fread(&N,sizeof(int),1,fid); - fseek(fid,M*N*sizeof(IssmDouble),SEEK_CUR); - } - else if(type==4){ - fread(&N,sizeof(int),1,fid); - fseek(fid,M*N*sizeof(int),SEEK_CUR); - } - else _error_("cannot read data of type " << type << "\n"); - } - } - if(found==0)_error_("cound not find " << field << " at step " << step << "\n"); - - /*assign output pointers:*/ - *pdoublemat=dmatfield; - *pdoublematsize=M*N; - *pdouble=dfield; - - /*return:*/ - return type; - -} -/*}}}*/ Index: /issm/branches/trunk-larour-SLPS2020/src/c/modules/ModelProcessorx/Dakota/CreateParametersDakota.cpp =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/modules/ModelProcessorx/Dakota/CreateParametersDakota.cpp (revision 25595) +++ /issm/branches/trunk-larour-SLPS2020/src/c/modules/ModelProcessorx/Dakota/CreateParametersDakota.cpp (revision 25596) @@ -41,4 +41,9 @@ int M,N; + //qmu statistics + bool statistics = false; + int numdirectories = 0; + int nfilesperdirectory = 0; + /*recover parameters: */ iomodel->FindConstant(&dakota_analysis,"md.qmu.isdakota"); @@ -72,4 +77,15 @@ /*Ok, we have all the response descriptors. Build a parameter with it: */ parameters->AddObject(new StringArrayParam(QmuResponsedescriptorsEnum,responsedescriptors,numresponsedescriptors)); + + /*Deal with statistics: */ + iomodel->FindConstant(&statistics,"md.qmu.statistics"); + parameters->AddObject(new BoolParam(QmuStatisticsEnum,statistics)); + if(statistics){ + iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories"); + parameters->AddObject(new IntParam(QmuNdirectoriesEnum,numdirectories)); + + iomodel->FindConstant(&nfilesperdirectory,"md.qmu.statistics.nfiles_per_directory"); + parameters->AddObject(new IntParam(QmuNfilesPerDirectoryEnum,nfilesperdirectory)); + } /*Load partitioning vectors specific to variables:*/ @@ -109,5 +125,4 @@ /*}}}*/ - /*Deal with data needed because of qmu variables*/ DataSet* dataset_variable_descriptors = new DataSet(QmuVariableDescriptorsEnum); @@ -143,5 +158,5 @@ delete dataset_variable_descriptors; - /*clean-up*/ + /*clean-up {{{*/ for(i=0;i(qmuerrname); xDelete(qmuoutname); - - - + xDelete(name); + /*}}}*/ } - /*Free data*/ - xDelete(name); } Index: /issm/branches/trunk-larour-SLPS2020/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp (revision 25595) +++ /issm/branches/trunk-larour-SLPS2020/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp (revision 25596) @@ -19,5 +19,5 @@ int i,j,n; - const int maxels=50; + const int maxels=100; const int width=maxels+1; Index: /issm/branches/trunk-larour-SLPS2020/src/c/modules/OutputResultsx/OutputResultsx.cpp =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/modules/OutputResultsx/OutputResultsx.cpp (revision 25595) +++ /issm/branches/trunk-larour-SLPS2020/src/c/modules/OutputResultsx/OutputResultsx.cpp (revision 25596) @@ -63,7 +63,25 @@ else{ if(my_rank==0){ + /*a little bit complicated. Either statistic computations are requested, which means we + * put our outbin files in subidirectories with numbers, or we don't, and we dump our + * outbins directly in the current directory:*/ int currEvalId ; + int nfilesperdirectory; + bool statistics=false; + char* root=NULL; + char* modelname=NULL; + femmodel->parameters->FindParam(&currEvalId,QmuCurrEvalIdEnum); - sprintf(outputfilename2,"%s.%i",outputfilename,currEvalId); + femmodel->parameters->FindParam(&statistics,QmuStatisticsEnum); + femmodel->parameters->FindParam(&nfilesperdirectory,QmuNfilesPerDirectoryEnum); + + if(statistics){ + femmodel->parameters->FindParam(&root,RootPathEnum); + femmodel->parameters->FindParam(&modelname,ModelnameEnum); + sprintf(outputfilename2,"%s/%i/%s.outbin.%i",root,(int)(floor((currEvalId-1)/nfilesperdirectory)+1),modelname,currEvalId); + } + else{ + sprintf(outputfilename2,"%s.%i",outputfilename,currEvalId); + } fid=pfopen0(outputfilename2,"ab+"); } Index: /issm/branches/trunk-larour-SLPS2020/src/c/modules/QmuStatisticsx/QmuStatisticsx.cpp =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/modules/QmuStatisticsx/QmuStatisticsx.cpp (revision 25596) +++ /issm/branches/trunk-larour-SLPS2020/src/c/modules/QmuStatisticsx/QmuStatisticsx.cpp (revision 25596) @@ -0,0 +1,1170 @@ +/*!\file: QmuStatisticsx routines + */ +/*includes and prototypes:*/ +#include "./QmuStatisticsx.h" +#include "../OutputResultsx/OutputResultsx.h" + +int readdata(IssmDouble** pdoublemat, int* pdoublematsize, IssmDouble* pdouble, FILE* fid,char* field,int step){ /*{{{*/ + + int length; + char fieldname[1000]; + int fieldname_size; + IssmDouble rtime; + int rstep; + int M,N; + + //fields that we retrive: + IssmDouble dfield; + char* sfield = NULL; + IssmDouble* dmatfield = NULL; + int* imatfield = NULL; + + //type of the returned field: + int type; + int found=0; + + while(1){ + + size_t ret_code = fread(&fieldname_size, sizeof(int), 1, fid); + if(ret_code != 1) break; //we are done. + + fread(fieldname, sizeof(char), fieldname_size, fid); + //_printf0_("fieldname: " << fieldname << "\n"); + + fread(&rtime, sizeof(IssmDouble), 1, fid); + fread(&rstep, sizeof(int), 1, fid); + + //check on field: + if ((step==rstep) && (strcmp(field,fieldname)==0)){ + + //ok, go read the result really: + fread(&type,sizeof(int),1,fid); + fread(&M,sizeof(int),1,fid); + if (type==1){ + fread(&dfield,sizeof(IssmDouble),1,fid); + } + else if (type==2){ + fread(&M,sizeof(int),1,fid); + sfield=xNew(M); + fread(sfield,sizeof(char),M,fid); + } + else if (type==3){ + fread(&N,sizeof(int),1,fid); + dmatfield=xNew(M*N); + fread(dmatfield,sizeof(IssmDouble),M*N,fid); + } + else if (type==4){ + fread(&N,sizeof(int),1,fid); + imatfield=xNew(M*N); + fread(imatfield,sizeof(int),M*N,fid); + } + else _error_("cannot read data of type " << type << "\n"); + found=1; + break; + } + else{ + //just skim to next results. + fread(&type,sizeof(int),1,fid); + fread(&M,sizeof(int),1,fid); + if (type==1){ + fseek(fid,sizeof(IssmDouble),SEEK_CUR); + } + else if(type==2){ + fseek(fid,M*sizeof(char),SEEK_CUR); + } + else if(type==3){ + fread(&N,sizeof(int),1,fid); + fseek(fid,M*N*sizeof(IssmDouble),SEEK_CUR); + } + else if(type==4){ + fread(&N,sizeof(int),1,fid); + fseek(fid,M*N*sizeof(int),SEEK_CUR); + } + else _error_("cannot read data of type " << type << "\n"); + } + } + if(found==0)_error_("cound not find " << field << " at step " << step << "\n"); + + /*assign output pointers:*/ + *pdoublemat=dmatfield; + *pdoublematsize=M*N; + *pdouble=dfield; + + /*return:*/ + return type; + +} +/*}}}*/ +int ComputeHistogram(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){ /*{{{*/ + + int nsamples; + char* directory=NULL; + char* model=NULL; + char** fields=NULL; + int* steps=NULL; + int nsteps; + int nfields; + int nbins; + int range,lower_row,upper_row; + int nfilesperdirectory; + + /*intermediary:*/ + IssmDouble* doublemat=NULL; + int doublematsize; + IssmDouble scalar; + + /*computation of average and variance itself:*/ + IssmDouble** maxxs = NULL; + IssmDouble** minxs = NULL; + int* xtype=NULL; + int* xsize=NULL; + + IssmDouble** maxmeans=NULL; + IssmDouble** minmeans=NULL; + int* meanxtype=NULL; + int* meanxsize=NULL; + + /*only work on the statistical communicator: */ + if (color==MPI_UNDEFINED)return 0; + + /*Retrieve parameters:*/ + parameters->FindParam(&nfilesperdirectory,QmuNfilesPerDirectoryEnum); + parameters->FindParam(&nsamples,QmuNsampleEnum); + parameters->FindParam(&directory,DirectoryNameEnum); + parameters->FindParam(&model,InputFileNameEnum); + parameters->FindParam(&fields,&nfields,FieldsEnum); + parameters->FindParam(&steps,&nsteps,StepsEnum); + parameters->FindParam(&nbins,NbinsEnum); + + /*Get rank from the stat comm communicator:*/ + IssmComm::SetComm(statcomm); + int my_rank=IssmComm::GetRank(); + + /*Open files and read them complelety, in a distributed way:*/ + range=DetermineLocalSize(nsamples,IssmComm::GetComm()); + GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm()); + + /*Initialize arrays:*/ + maxmeans=xNew(nfields); + minmeans=xNew(nfields); + meanxtype=xNew(nfields); + meanxsize=xNew(nfields); + + maxxs=xNew(nfields*nsteps); + minxs=xNew(nfields*nsteps); + xtype=xNew(nfields*nsteps); + xsize=xNew(nfields*nsteps); + + /*Start opening files:*/ + for(int i=(lower_row+1);i<=upper_row;i++){ + _printf0_("reading file #: " << i << "\n"); + char file[1000]; + long int length; + char* buffer=NULL; + + /*string:*/ + sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i); + + /*open file: */ + _printf0_(" opening file: " << file << "\n"); + FILE* fid=fopen(file,"rb"); + if(fid==NULL)_error_("cound not open file: " << file << "\n"); + + /*figure out size of file, and read the whole thing:*/ + _printf0_(" reading file:\n"); + fseek (fid, 0, SEEK_END); + length = ftell (fid); + fseek (fid, 0, SEEK_SET); + buffer = xNew(length); + fread (buffer, sizeof(char), length, fid); + + /*close file:*/ + fclose (fid); + + /*create a memory stream with this buffer:*/ + _printf0_(" processing file:\n"); + fid=fmemopen(buffer, length, "rb"); + + /*start reading data from the buffer directly:*/ + for (int f=0;f(1); + minxs[counter]=xNew(1); + *maxxs[counter]=scalar; + *minxs[counter]=scalar; + xsize[counter]=1; + } + else if (xtype[counter]==3){ + maxxs[counter]=xNew(doublematsize); + xMemCpy(maxxs[counter],doublemat,doublematsize); + minxs[counter]=xNew(doublematsize); + xMemCpy(minxs[counter],doublemat,doublematsize); + xsize[counter]=doublematsize; + xDelete(doublemat); + } + else _error_("cannot carry out statistics on type " << xtype[counter]); + } + else{ + if(xtype[counter]==1){ + *maxxs[counter]=max(*maxxs[counter],scalar); + *minxs[counter]=min(*minxs[counter],scalar); + } + else if (xtype[counter]==3){ + IssmDouble* newmax=maxxs[counter]; + IssmDouble* newmin=minxs[counter]; + for(int k=0;knewmax[k])newmax[k]=doublemat[k]; + if(doublemat[k](doublemat); + } + else _error_("cannot carry out statistics on type " << xtype[counter]); + } + } + } + _printf0_(" average in time:\n"); + + /*Deal with average in time: */ + for (int f=0;f(1); + minmeans[f]=xNewZeroInit(1); + *maxmeans[f]=timemean; + *minmeans[f]=timemean; + } + else{ + *maxmeans[f]=max(*maxmeans[f],timemean); + *minmeans[f]=min(*minmeans[f],timemean); + } + } + else{ + meanxsize[f]=doublematsize; + fseek(fid,0,SEEK_SET); + IssmDouble* timemean=xNewZeroInit(doublematsize); + for (int j=0;j(doublemat); + } + + if(i==(lower_row+1)){ + maxmeans[f]=xNew(doublematsize); + xMemCpy(maxmeans[f],timemean,doublematsize); + minmeans[f]=xNew(doublematsize); + xMemCpy(minmeans[f],timemean,doublematsize); + } + else{ + IssmDouble* maxx=maxmeans[f]; + IssmDouble* minx=minmeans[f]; + + for(int k=0;k(buffer); + } + ISSM_MPI_Barrier(IssmComm::GetComm()); + _printf0_("Done reading files, now computing min and max.\n"); + + /*We have agregated minx and max across the cluster, now gather across the cluster onto + *cpu0 and then compute statistics:*/ + for (int f=0;f(size); + IssmDouble* allmin=xNew(size); + + ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm()); + ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm()); + + /*Store broadcasted value for later computation of histograms:*/ + maxxs[counter]=allmax; + minxs[counter]=allmin; + } + } /*}}}*/ + } + + /*Now do the same for the time mean fields:*/ + for (int f=0;f(size); + IssmDouble* allmin=xNew(size); + + ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm()); + ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm()); + + + /*Store for later use in histogram computation:*/ + maxmeans[f]=allmax; + minmeans[f]=allmin; + + } /*}}}*/ + } + + /*Now that we have the min and max, we can start binning. First allocate + * histograms, then start filling them:*/ + IssmDouble** histogram=xNew(nfields*nsteps); + IssmDouble** timehistogram=xNew(nfields); + + _printf0_("Start reading files again, this time binning values in the histogram:\n"); + /*Start opening files:*/ + for (int i=(lower_row+1);i<=upper_row;i++){ + _printf0_("reading file #: " << i << "\n"); + char file[1000]; + long int length; + char* buffer=NULL; + + /*string:*/ + sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i); + + /*open file: */ + _printf0_(" opening file:\n"); + FILE* fid=fopen(file,"rb"); + if(fid==NULL)_error_("cound not open file: " << file << "\n"); + + /*figure out size of file, and read the whole thing:*/ + _printf0_(" reading file:\n"); + fseek (fid, 0, SEEK_END); + length = ftell (fid); + fseek (fid, 0, SEEK_SET); + buffer = xNew(length); + fread (buffer, sizeof(char), length, fid); + + /*close file:*/ + fclose (fid); + + /*create a memory stream with this buffer:*/ + _printf0_(" processing file:\n"); + fid=fmemopen(buffer, length, "rb"); + + /*start reading data from the buffer directly:*/ + for (int f=0;f(nbins); + IssmDouble ma=*maxxs[counter]; + IssmDouble mi=*minxs[counter]; + int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index--; + localhistogram[index]++; + histogram[counter]=localhistogram; + } + else if (xtype[counter]==3){ + IssmDouble* localhistogram=xNewZeroInit(doublematsize*nbins); + IssmDouble* ma=maxxs[counter]; + IssmDouble* mi=minxs[counter]; + for (int k=0;k=mi[k]); _assert_(index(doublemat); + } + else _error_("cannot carry out statistics on type " << xtype[counter]); + } + else{ + if(xtype[counter]==1){ + IssmDouble* localhistogram=histogram[counter]; + IssmDouble ma=*maxxs[counter]; + IssmDouble mi=*minxs[counter]; + int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1; + localhistogram[index]++; + } + else if (xtype[counter]==3){ + IssmDouble* localhistogram=histogram[counter]; + IssmDouble* ma=maxxs[counter]; + IssmDouble* mi=minxs[counter]; + for (int k=0;k(doublemat); + } + else _error_("cannot carry out statistics on type " << xtype[counter]); + } + } + } + _printf0_(" average in time:\n"); + + /*Deal with average in time: */ + for (int f=0;f(nbins); + IssmDouble ma=*maxmeans[f]; + IssmDouble mi=*minmeans[f]; + int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1; + localhistogram[index]++; + timehistogram[f]=localhistogram; + } + else{ + IssmDouble* localhistogram=timehistogram[f]; + IssmDouble ma=*maxmeans[f]; + IssmDouble mi=*minmeans[f]; + int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1; + localhistogram[index]++; + } + } + else{ + fseek(fid,0,SEEK_SET); + IssmDouble* timemean=xNewZeroInit(doublematsize); + for (int j=0;j(doublemat); + } + + if(i==(lower_row+1)){ + IssmDouble* localhistogram=xNewZeroInit(doublematsize*nbins); + IssmDouble* ma=maxmeans[f]; + IssmDouble* mi=minmeans[f]; + + for (int k=0;k(buffer); + } + _printf0_("Start aggregating histogram:\n"); + + /*We have agregated histograms across the cluster, now gather them across the cluster onto + *cpu0: */ + for (int f=0;f(nbins); + + ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); + + /*add to results:*/ + if(my_rank==0){ + char fieldname[1000]; + + sprintf(fieldname,"%s%s",fields[f],"Histogram"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,allhisto,1,nbins,j+1,0)); + + sprintf(fieldname,"%s%s",fields[f],"Max"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,*maxxs[counter],j+1,0)); + sprintf(fieldname,"%s%s",fields[f],"Min"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,*minxs[counter],j+1,0)); + } + } + } /*}}}*/ + else{ /*deal with arrays:{{{*/ + + int size=xsize[counter0]; + for (int j=0;j(size*nbins); + + ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); + xDelete(histo); + + /*add to results:*/ + if(my_rank==0){ + char fieldname[1000]; + + sprintf(fieldname,"%s%s",fields[f],"Histogram"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,allhisto,size,nbins,j+1,0)); + + sprintf(fieldname,"%s%s",fields[f],"Max"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,maxxs[counter],size,1,j+1,0)); + sprintf(fieldname,"%s%s",fields[f],"Min"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,minxs[counter],size,1,j+1,0)); + } + } + } /*}}}*/ + } + _printf0_("Start aggregating time mean histogram:\n"); + + /*Now do the same for the time mean fields:*/ + for (int f=0;f(nbins); + + ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm()); + + /*add to results at time step 1:*/ + if(my_rank==0){ + char fieldname[1000]; + + sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,allhisto,1,nbins,1,0)); + + sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,*maxmeans[f],1,0)); + sprintf(fieldname,"%s%s",fields[f],"TimeMeaMin"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,*minmeans[f],1,0)); + } + } /*}}}*/ + else{ /*deal with arrays:{{{*/ + + int size=meanxsize[f]; + + /*we are broadcasting double arrays:*/ + IssmDouble* histo=timehistogram[f]; + IssmDouble* allhisto=xNewZeroInit(size*nbins); + + ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); + xDelete(histo); + /*add to results at step 1:*/ + if(my_rank==0){ + char fieldname[1000]; + + sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,allhisto,size,nbins,1,0)); + sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,maxmeans[f],size,1,1,0)); + sprintf(fieldname,"%s%s",fields[f],"TimeMeanMin"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,minmeans[f],size,1,1,0)); + } + } /*}}}*/ + } + _printf0_("Done aggregating time mean histogram:\n"); + IssmComm::SetComm(ISSM_MPI_COMM_WORLD); +} +/*}}}*/ +int ComputeMeanVariance(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){ /*{{{*/ + + int nsamples; + char* directory=NULL; + char* model=NULL; + char** fields=NULL; + int* steps=NULL; + int nsteps; + int nfields; + int range,lower_row,upper_row; + + /*intermediary:*/ + IssmDouble* doublemat=NULL; + int doublematsize; + IssmDouble scalar; + + /*computation of average and variance itself:*/ + IssmDouble* x = NULL; + IssmDouble* x2 = NULL; + IssmDouble** xs = NULL; + IssmDouble** xs2 = NULL; + int* xtype=NULL; + int* xsize=NULL; + + IssmDouble** meanx=NULL; + IssmDouble** meanx2=NULL; + int* meantype=NULL; + int* meansize=NULL; + + /*Retrieve parameters:*/ + parameters->FindParam(&nsamples,QmuNsampleEnum); + parameters->FindParam(&directory,DirectoryNameEnum); + parameters->FindParam(&model,InputFileNameEnum); + parameters->FindParam(&fields,&nfields,FieldsEnum); + parameters->FindParam(&steps,&nsteps,StepsEnum); + + /*Get rank:*/ + int my_rank=IssmComm::GetRank(); + + /*Open files and read them complelety, in a distributed way:*/ + range=DetermineLocalSize(nsamples,IssmComm::GetComm()); + GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm()); + + /*Initialize arrays:*/ + xs=xNew(nfields*nsteps); + xs2=xNew(nfields*nsteps); + xtype=xNew(nfields*nsteps); + xsize=xNew(nfields*nsteps); + + meantype=xNew(nfields); + meansize=xNew(nfields); + meanx=xNew(nfields); + meanx2=xNew(nfields); + + /*Start opening files:*/ + for (int i=(lower_row+1);i<=upper_row;i++){ + _printf0_("reading file #: " << i << "\n"); + char file[1000]; + long int length; + char* buffer=NULL; + + /*string:*/ + sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i); + + /*open file: */ + _printf0_(" opening file: " << file << "\n"); + FILE* fid=fopen(file,"rb"); + if(fid==NULL) _error_(" could not open file: " << file << "\n"); + + /*figure out size of file, and read the whole thing:*/ + _printf0_(" reading file:\n"); + fseek (fid, 0, SEEK_END); + length = ftell (fid); + fseek (fid, 0, SEEK_SET); + buffer = xNew(length); + fread (buffer, sizeof(char), length, fid); + + /*close file:*/ + fclose (fid); + + /*create a memory stream with this buffer:*/ + _printf0_(" processing file:\n"); + fid=fmemopen(buffer, length, "rb"); + + /*start reading data from the buffer directly:*/ + for (int f=0;f(1); + xs2[counter]=xNew(1); + *xs[counter]=scalar; + *xs2[counter]=pow(scalar,2.0); + xsize[counter]=1; + } + else if (xtype[counter]==3){ + IssmDouble* doublemat2=xNew(doublematsize); + for(int k=0;k(1); + meanx2[f]=xNewZeroInit(1); + meansize[f]=1; + } + else{ + meanx[f]=xNewZeroInit(doublematsize); + meanx2[f]=xNewZeroInit(doublematsize); + meansize[f]=doublematsize; + } + } + fseek(fid,0,SEEK_SET); + if(meantype[f]==1){ + IssmDouble sc=0; + IssmDouble sc2=0; + for(int j=0;j(doublematsize); + IssmDouble* timemean2=xNewZeroInit(doublematsize); + + for(int j=0;j(buffer); + } + ISSM_MPI_Barrier(IssmComm::GetComm()); + _printf0_("Done reading files, now computing mean and variance.\n"); + + /*We have agregated x and x^2 across the cluster, now gather across the cluster onto + *cpu0 and then compute statistics:*/ + for (int f=0;fAddResult(new GenericExternalResult(results->Size()+1,fieldname,mean,j+1,0)); + sprintf(fieldname,"%s%s",fields[f],"Stddev"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,stddev,j+1,0)); + } + + } + } /*}}}*/ + else{ /*deal with arrays:{{{*/ + + int size=xsize[counter0]; + + IssmDouble* mean=xNew(size); + IssmDouble* stddev=xNew(size); + + for (int j=0;j(size); + IssmDouble* sumx2=xNew(size); + + ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm()); + ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm()); + + /*Build average and standard deviation. For standard deviation, use the + *following formula: sigma^2=E(x^2)-mu^2:*/ + for (int k=0;kAddResult(new GenericExternalResult(results->Size()+1,fieldname,mean,size,1,j+1,0)); + sprintf(fieldname,"%s%s",fields[f],"Stddev"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,stddev,size,1,j+1,0)); + } + } + } /*}}}*/ + } + /*Do the same but for the time mean:*/ + for (int f=0;fAddResult(new GenericExternalResult(results->Size()+1,fieldname,mean,1,0)); + sprintf(fieldname,"%s%s",fields[f],"TimeStddev"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,stddev,1,0)); + } + } /*}}}*/ + else{ /*deal with arrays:{{{*/ + + int size=meansize[f]; + IssmDouble* mean=xNew(size); + IssmDouble* stddev=xNew(size); + + /*we are broadcasting double arrays:*/ + x=meanx[f]; + x2=meanx2[f]; + + IssmDouble* sumx=xNew(size); + IssmDouble* sumx2=xNew(size); + + ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm()); + ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm()); + + /*Build average and standard deviation. For standard deviation, use the + *following formula: sigma^2=E(x^2)-mu^2:*/ + for (int k=0;kAddResult(new GenericExternalResult(results->Size()+1,fieldname,mean,size,1,1,0)); + sprintf(fieldname,"%s%s",fields[f],"TimeStddev"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,stddev,size,1,1,0)); + } + } /*}}}*/ + } + + +} /*}}}*/ +int ComputeSampleSeries(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){ /*{{{*/ + + int nsamples; + char* directory=NULL; + char* model=NULL; + char** fields=NULL; + int* steps=NULL; + int nsteps; + int nfields; + int range,lower_row,upper_row; + int* indices=NULL; + int nindices; + + /*intermediary:*/ + IssmDouble* doublemat=NULL; + int doublematsize; + IssmDouble scalar; + + /*computation of average and variance itself:*/ + IssmDouble* x = NULL; + IssmDouble* allx=NULL; + IssmDouble** xs = NULL; + int* xtype=NULL; + int* xsize=NULL; + + /*Retrieve parameters:*/ + parameters->FindParam(&nsamples,QmuNsampleEnum); + parameters->FindParam(&directory,DirectoryNameEnum); + parameters->FindParam(&model,InputFileNameEnum); + parameters->FindParam(&fields,&nfields,FieldsEnum); + parameters->FindParam(&steps,&nsteps,StepsEnum); + parameters->FindParam(&indices,&nindices,IndicesEnum); + + /*Get rank:*/ + int my_rank=IssmComm::GetRank(); + + /*Open files and read them complelety, in a distributed way:*/ + range=DetermineLocalSize(nsamples,IssmComm::GetComm()); + GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm()); + + /*Initialize arrays:*/ + xs=xNew(nfields*nsteps); + xtype=xNew(nfields*nsteps); + xsize=xNew(nfields*nsteps); + + /*Start opening files:*/ + for (int i=(lower_row+1);i<=upper_row;i++){ + _printf0_("reading file #: " << i << "\n"); + char file[1000]; + long int length; + char* buffer=NULL; + + /*string:*/ + sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i); + + /*open file: */ + _printf0_(" opening file:\n"); + FILE* fid=fopen(file,"rb"); + + /*figure out size of file, and read the whole thing:*/ + _printf0_(" reading file:\n"); + fseek (fid, 0, SEEK_END); + length = ftell (fid); + fseek (fid, 0, SEEK_SET); + buffer = xNew(length); + fread (buffer, sizeof(char), length, fid); + + /*close file:*/ + fclose (fid); + + /*create a memory stream with this buffer:*/ + _printf0_(" processing file:\n"); + fid=fmemopen(buffer, length, "rb"); + + /*start reading data from the buffer directly:*/ + for (int f=0;f(range); + x[0]=scalar; + xs[counter]=x; + xsize[counter]=range; + } + else if (xtype[counter]==3){ + x=xNew(nindices*range); + for(int k=0;k(buffer); + } + ISSM_MPI_Barrier(IssmComm::GetComm()); + _printf0_("Done reading files, now assembling time series.\n"); + + for (int f=0;f(nsamples); + MPI_Gather(x, range, ISSM_MPI_PDOUBLE,allx, range, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm()); + /*add to results:*/ + if(my_rank==0){ + char fieldname[1000]; + + sprintf(fieldname,"%s%s",fields[f],"Samples"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,allx,nsamples,1,j+1,0)); + } + } + else{ + /*we are broadcasting double arrays:*/ + x=xs[counter]; + allx=xNew(nsamples*nindices); + + MPI_Gather(x, range*nindices, ISSM_MPI_PDOUBLE,allx, range*nindices, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm()); + + /*add to results:*/ + if(my_rank==0){ + char fieldname[1000]; + sprintf(fieldname,"%s%s",fields[f],"Samples"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,allx,nsamples,nindices,j+1,0)); + } + } + } + } + + +} /*}}}*/ +int OutputStatistics(Parameters* parameters,Results* results){ /*{{{*/ + + char outputfilename[1000]; + char* directory=NULL; + char* model=NULL; + char* method=NULL; + int nsamples; + int* steps=NULL; + int nsteps; + + FemModel* femmodel=new FemModel(); + + /*Some parameters that will allow us to use the OutputResultsx module:*/ + parameters->AddObject(new BoolParam(QmuIsdakotaEnum,false)); + parameters->AddObject(new BoolParam(SettingsIoGatherEnum,true)); + + parameters->FindParam(&directory,DirectoryNameEnum); + parameters->FindParam(&model,InputFileNameEnum); + parameters->FindParam(&nsamples,QmuNsampleEnum); + parameters->FindParam(&steps,&nsteps,StepsEnum); + + sprintf(outputfilename,"%s/%s.stats",directory,model); + parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename)); + + /*Call OutputResults module:*/ + femmodel->parameters=parameters; + femmodel->results=results; + + OutputResultsx(femmodel); +} /*}}}*/ Index: /issm/branches/trunk-larour-SLPS2020/src/c/modules/QmuStatisticsx/QmuStatisticsx.h =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/modules/QmuStatisticsx/QmuStatisticsx.h (revision 25596) +++ /issm/branches/trunk-larour-SLPS2020/src/c/modules/QmuStatisticsx/QmuStatisticsx.h (revision 25596) @@ -0,0 +1,17 @@ +/*!\file: QmuStatisticsx.h + */ + +#ifndef _QMU_STATISTCSX_H_ +#define _QMU_STATISTCSX_H_ + +#include "../../classes/classes.h" + +int ComputeMeanVariance(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm); +int ComputeSampleSeries(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm); +int OutputStatistics(Parameters* parameters,Results* results); +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); + +/* local prototypes: */ + +#endif /* _QMU_STATISTCSX_H_ */ Index: /issm/branches/trunk-larour-SLPS2020/src/c/modules/modules.h =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/modules/modules.h (revision 25595) +++ /issm/branches/trunk-larour-SLPS2020/src/c/modules/modules.h (revision 25596) @@ -81,4 +81,5 @@ #include "./PointCloudFindNeighborsx/PointCloudFindNeighborsx.h" #include "./PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h" +#include "./QmuStatisticsx/QmuStatisticsx.h" #include "./Reduceloadx/Reduceloadx.h" #include "./Reducevectorgtofx/Reducevectorgtofx.h" Index: /issm/branches/trunk-larour-SLPS2020/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/shared/Enum/Enum.vim (revision 25595) +++ /issm/branches/trunk-larour-SLPS2020/src/c/shared/Enum/Enum.vim (revision 25596) @@ -314,7 +314,14 @@ syn keyword cConstant QmuResponsePartitionsEnum syn keyword cConstant QmuResponsePartitionsNpartEnum +syn keyword cConstant QmuStatisticsEnum +syn keyword cConstant QmuNumstatisticsEnum +syn keyword cConstant QmuNdirectoriesEnum +syn keyword cConstant QmuNfilesPerDirectoryEnum +syn keyword cConstant QmuStatisticsMethodEnum +syn keyword cConstant QmuMethodsEnum syn keyword cConstant RestartFileNameEnum syn keyword cConstant ResultsEnum syn keyword cConstant RootPathEnum +syn keyword cConstant ModelnameEnum syn keyword cConstant SaveResultsEnum syn keyword cConstant SolidearthPlanetRadiusEnum @@ -1458,4 +1465,5 @@ syn keyword cType PowerVariogram syn keyword cType Profiler +syn keyword cType QmuStatisticsMethod syn keyword cType Quadtree syn keyword cType Radar Index: /issm/branches/trunk-larour-SLPS2020/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/shared/Enum/EnumDefinitions.h (revision 25595) +++ /issm/branches/trunk-larour-SLPS2020/src/c/shared/Enum/EnumDefinitions.h (revision 25596) @@ -308,7 +308,14 @@ QmuResponsePartitionsEnum, QmuResponsePartitionsNpartEnum, + QmuStatisticsEnum, + QmuNumstatisticsEnum, + QmuNdirectoriesEnum, + QmuNfilesPerDirectoryEnum, + QmuStatisticsMethodEnum, + QmuMethodsEnum, RestartFileNameEnum, ResultsEnum, RootPathEnum, + ModelnameEnum, SaveResultsEnum, SolidearthPlanetRadiusEnum, Index: /issm/branches/trunk-larour-SLPS2020/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/shared/Enum/EnumToStringx.cpp (revision 25595) +++ /issm/branches/trunk-larour-SLPS2020/src/c/shared/Enum/EnumToStringx.cpp (revision 25596) @@ -316,7 +316,14 @@ case QmuResponsePartitionsEnum : return "QmuResponsePartitions"; case QmuResponsePartitionsNpartEnum : return "QmuResponsePartitionsNpart"; + case QmuStatisticsEnum : return "QmuStatistics"; + case QmuNumstatisticsEnum : return "QmuNumstatistics"; + case QmuNdirectoriesEnum : return "QmuNdirectories"; + case QmuNfilesPerDirectoryEnum : return "QmuNfilesPerDirectory"; + case QmuStatisticsMethodEnum : return "QmuStatisticsMethod"; + case QmuMethodsEnum : return "QmuMethods"; case RestartFileNameEnum : return "RestartFileName"; case ResultsEnum : return "Results"; case RootPathEnum : return "RootPath"; + case ModelnameEnum : return "Modelname"; case SaveResultsEnum : return "SaveResults"; case SolidearthPlanetRadiusEnum : return "SolidearthPlanetRadius"; Index: /issm/branches/trunk-larour-SLPS2020/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/branches/trunk-larour-SLPS2020/src/c/shared/Enum/StringToEnumx.cpp (revision 25595) +++ /issm/branches/trunk-larour-SLPS2020/src/c/shared/Enum/StringToEnumx.cpp (revision 25596) @@ -322,7 +322,14 @@ else if (strcmp(name,"QmuResponsePartitions")==0) return QmuResponsePartitionsEnum; else if (strcmp(name,"QmuResponsePartitionsNpart")==0) return QmuResponsePartitionsNpartEnum; + else if (strcmp(name,"QmuStatistics")==0) return QmuStatisticsEnum; + else if (strcmp(name,"QmuNumstatistics")==0) return QmuNumstatisticsEnum; + else if (strcmp(name,"QmuNdirectories")==0) return QmuNdirectoriesEnum; + else if (strcmp(name,"QmuNfilesPerDirectory")==0) return QmuNfilesPerDirectoryEnum; + else if (strcmp(name,"QmuStatisticsMethod")==0) return QmuStatisticsMethodEnum; + else if (strcmp(name,"QmuMethods")==0) return QmuMethodsEnum; else if (strcmp(name,"RestartFileName")==0) return RestartFileNameEnum; else if (strcmp(name,"Results")==0) return ResultsEnum; else if (strcmp(name,"RootPath")==0) return RootPathEnum; + else if (strcmp(name,"Modelname")==0) return ModelnameEnum; else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum; else if (strcmp(name,"SolidearthPlanetRadius")==0) return SolidearthPlanetRadiusEnum; @@ -376,5 +383,8 @@ else if (strcmp(name,"SmbDsnowIdx")==0) return SmbDsnowIdxEnum; else if (strcmp(name,"SmbCldFrac")==0) return SmbCldFracEnum; - else if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum; + else stage=4; + } + if(stage==4){ + if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum; else if (strcmp(name,"SmbDelta18oSurface")==0) return SmbDelta18oSurfaceEnum; else if (strcmp(name,"SmbDenIdx")==0) return SmbDenIdxEnum; @@ -383,8 +393,5 @@ else if (strcmp(name,"SmbF")==0) return SmbFEnum; else if (strcmp(name,"SmbInitDensityScaling")==0) return SmbInitDensityScalingEnum; - else stage=4; - } - if(stage==4){ - if (strcmp(name,"SmbIsaccumulation")==0) return SmbIsaccumulationEnum; + else if (strcmp(name,"SmbIsaccumulation")==0) return SmbIsaccumulationEnum; else if (strcmp(name,"SmbIsalbedo")==0) return SmbIsalbedoEnum; else if (strcmp(name,"SmbIsclimatology")==0) return SmbIsclimatologyEnum; @@ -499,5 +506,8 @@ else if (strcmp(name,"BalancethicknessOmega")==0) return BalancethicknessOmegaEnum; else if (strcmp(name,"BalancethicknessThickeningRate")==0) return BalancethicknessThickeningRateEnum; - else if (strcmp(name,"BasalCrevasse")==0) return BasalCrevasseEnum; + else stage=5; + } + if(stage==5){ + if (strcmp(name,"BasalCrevasse")==0) return BasalCrevasseEnum; else if (strcmp(name,"BasalforcingsFloatingiceMeltingRate")==0) return BasalforcingsFloatingiceMeltingRateEnum; else if (strcmp(name,"BasalforcingsGeothermalflux")==0) return BasalforcingsGeothermalfluxEnum; @@ -506,8 +516,5 @@ else if (strcmp(name,"BasalforcingsIsmip6BasinId")==0) return BasalforcingsIsmip6BasinIdEnum; else if (strcmp(name,"BasalforcingsIsmip6Tf")==0) return BasalforcingsIsmip6TfEnum; - else stage=5; - } - if(stage==5){ - if (strcmp(name,"BasalforcingsIsmip6TfShelf")==0) return BasalforcingsIsmip6TfShelfEnum; + else if (strcmp(name,"BasalforcingsIsmip6TfShelf")==0) return BasalforcingsIsmip6TfShelfEnum; else if (strcmp(name,"BasalforcingsIsmip6MeltAnomaly")==0) return BasalforcingsIsmip6MeltAnomalyEnum; else if (strcmp(name,"BasalforcingsOceanSalinity")==0) return BasalforcingsOceanSalinityEnum; @@ -622,5 +629,8 @@ else if (strcmp(name,"UGia")==0) return UGiaEnum; else if (strcmp(name,"UGiaRate")==0) return UGiaRateEnum; - else if (strcmp(name,"Gradient")==0) return GradientEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"Gradient")==0) return GradientEnum; else if (strcmp(name,"GroundinglineHeight")==0) return GroundinglineHeightEnum; else if (strcmp(name,"HydraulicPotential")==0) return HydraulicPotentialEnum; @@ -629,8 +639,5 @@ else if (strcmp(name,"HydrologyBumpHeight")==0) return HydrologyBumpHeightEnum; else if (strcmp(name,"HydrologyBumpSpacing")==0) return HydrologyBumpSpacingEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"HydrologydcBasalMoulinInput")==0) return HydrologydcBasalMoulinInputEnum; + else if (strcmp(name,"HydrologydcBasalMoulinInput")==0) return HydrologydcBasalMoulinInputEnum; else if (strcmp(name,"HydrologydcEplThickness")==0) return HydrologydcEplThicknessEnum; else if (strcmp(name,"HydrologydcEplThicknessOld")==0) return HydrologydcEplThicknessOldEnum; @@ -745,5 +752,8 @@ else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum; else if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum; - else if (strcmp(name,"SmbBMin")==0) return SmbBMinEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"SmbBMin")==0) return SmbBMinEnum; else if (strcmp(name,"SmbBNeg")==0) return SmbBNegEnum; else if (strcmp(name,"SmbBPos")==0) return SmbBPosEnum; @@ -752,8 +762,5 @@ else if (strcmp(name,"SmbDailyairdensity")==0) return SmbDailyairdensityEnum; else if (strcmp(name,"SmbDailyairhumidity")==0) return SmbDailyairhumidityEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"SmbDailydlradiation")==0) return SmbDailydlradiationEnum; + else if (strcmp(name,"SmbDailydlradiation")==0) return SmbDailydlradiationEnum; else if (strcmp(name,"SmbDailydsradiation")==0) return SmbDailydsradiationEnum; else if (strcmp(name,"SmbDailypressure")==0) return SmbDailypressureEnum; @@ -868,5 +875,8 @@ else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum; else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum; - else if (strcmp(name,"Temperature")==0) return TemperatureEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"Temperature")==0) return TemperatureEnum; else if (strcmp(name,"TemperaturePDD")==0) return TemperaturePDDEnum; else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum; @@ -875,8 +885,5 @@ else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum; else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum; + else if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum; else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum; else if (strcmp(name,"Thickness")==0) return ThicknessEnum; @@ -991,5 +998,8 @@ else if (strcmp(name,"Outputdefinition87")==0) return Outputdefinition87Enum; else if (strcmp(name,"Outputdefinition88")==0) return Outputdefinition88Enum; - else if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum; else if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum; else if (strcmp(name,"Outputdefinition90")==0) return Outputdefinition90Enum; @@ -998,8 +1008,5 @@ else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum; else if (strcmp(name,"Outputdefinition94")==0) return Outputdefinition94Enum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"Outputdefinition95")==0) return Outputdefinition95Enum; + else if (strcmp(name,"Outputdefinition95")==0) return Outputdefinition95Enum; else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum; else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum; @@ -1114,5 +1121,8 @@ else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum; else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum; - else if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum; else if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum; else if (strcmp(name,"Fset")==0) return FsetEnum; @@ -1121,8 +1131,5 @@ else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum; else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum; + else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum; else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum; else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum; @@ -1237,5 +1244,8 @@ else if (strcmp(name,"Moulin")==0) return MoulinEnum; else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum; - else if (strcmp(name,"Mpi")==0) return MpiEnum; + else stage=11; + } + if(stage==11){ + if (strcmp(name,"Mpi")==0) return MpiEnum; else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum; else if (strcmp(name,"Mumps")==0) return MumpsEnum; @@ -1244,8 +1254,5 @@ else if (strcmp(name,"Nodal")==0) return NodalEnum; else if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum; - else stage=11; - } - if(stage==11){ - if (strcmp(name,"NodeSId")==0) return NodeSIdEnum; + else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum; else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum; else if (strcmp(name,"None")==0) return NoneEnum; @@ -1360,5 +1367,8 @@ else if (strcmp(name,"TotalGroundedBmbScaled")==0) return TotalGroundedBmbScaledEnum; else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum; - else if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum; + else stage=12; + } + if(stage==12){ + if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum; else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum; else if (strcmp(name,"TransientInput")==0) return TransientInputEnum; @@ -1367,8 +1377,5 @@ else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum; else if (strcmp(name,"Tria")==0) return TriaEnum; - else stage=12; - } - if(stage==12){ - if (strcmp(name,"TriaInput")==0) return TriaInputEnum; + else if (strcmp(name,"TriaInput")==0) return TriaInputEnum; else if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum; else if (strcmp(name,"VectorParam")==0) return VectorParamEnum;