Changeset 25596


Ignore:
Timestamp:
09/25/20 18:31:59 (5 years ago)
Author:
Eric.Larour
Message:

CHG: diverse

Location:
issm/branches/trunk-larour-SLPS2020/src/c
Files:
3 added
15 edited

Legend:

Unmodified
Added
Removed

  • issm/branches/trunk-larour-SLPS2020/src/c/Makefile.am

    r25349 r25596  
    635635        ./modules/NodeConnectivityx/NodeConnectivityx.cpp \
    636636        ./modules/ElementConnectivityx/ElementConnectivityx.cpp \
    637         ./modules/PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.cpp
     637        ./modules/PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.cpp \
     638        ./modules/QmuStatisticsx/QmuStatisticsx.cpp
    638639
    639640if CHACO

  • issm/branches/trunk-larour-SLPS2020/src/c/classes/FemModel.cpp

    r25534 r25596  
    7070        char *restartfilename  = NULL;
    7171        char *rootpath       = NULL;
     72        char *modelname       = NULL;
    7273
    7374        /*First things first, store the communicator, and set it as a global variable: */
    … …  
    8586
    8687        /*From command line arguments, retrieve different filenames needed to create the FemModel: */
    87         ProcessArguments(&solution_type,&binfilename,&outbinfilename,&petscfilename,&lockfilename,&restartfilename,&rootpath,argc,argv);
     88        ProcessArguments(&solution_type,&binfilename,&outbinfilename,&petscfilename,&lockfilename,&restartfilename,&rootpath,&modelname,argc,argv);
    8889
    8990        /*Create femmodel from input files: */
    9091        profiler->Start(MPROCESSOR);
    91         this->InitFromFiles(rootpath,binfilename,outbinfilename,petscfilename,lockfilename,restartfilename, solution_type,trace,NULL);
     92        this->InitFromFiles(rootpath,binfilename,outbinfilename,petscfilename,lockfilename,restartfilename, modelname, solution_type,trace,NULL);
    9293        profiler->Stop(MPROCESSOR);
    9394
    … …  
    148149}
    149150/*}}}*/
    150 FemModel::FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X){ /*{{{*/
     151FemModel::FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, char* modelname, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X){ /*{{{*/
    151152
    152153        bool traceon=true;
    … …  
    159160        /*Create femmodel from input files, with trace activated: */
    160161        profiler->Start(MPROCESSOR);
    161         this->InitFromFiles(rootpath,inputfilename,outputfilename,toolkitsfilename,lockfilename,restartfilename, solution_type,traceon,X);
     162        this->InitFromFiles(rootpath,inputfilename,outputfilename,toolkitsfilename,lockfilename,restartfilename, modelname,solution_type,traceon,X);
    162163        profiler->Stop(MPROCESSOR);
    163164
    … …  
    394395}
    395396/*}}}*/
    396 void FemModel::InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, const int in_solution_type,bool trace,IssmPDouble* X){/*{{{*/
     397void 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){/*{{{*/
    397398
    398399        /*intermediary*/
    … …  
    418419        /*Now save all of these file names into parameters, you never know when you might need them: */
    419420        this->parameters->AddObject(new StringParam(ToolkitsFileNameEnum,toolkitsfilename));
     421        this->parameters->AddObject(new StringParam(ModelnameEnum,modelname));
    420422        this->parameters->AddObject(new StringParam(RootPathEnum,rootpath));
    421423        this->parameters->AddObject(new StringParam(InputFileNameEnum,inputfilename));

  • issm/branches/trunk-larour-SLPS2020/src/c/classes/FemModel.h

    r25534 r25596  
    6969                FemModel(void);
    7070                FemModel(int argc,char** argv,ISSM_MPI_Comm comm_init,bool trace=false);
    71                 FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X);
     71                FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, char* modelname, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X);
    7272                ~FemModel();
    7373
    … …  
    7979                void Echo();
    8080                int  GetElementsWidth(){return 3;};//just tria elements in this first version
    81                 void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, char* restartfilename, const int solution_type,bool trace,IssmPDouble* X=NULL);
     81                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);
    8282                void InitFromFids(char* rootpath, FILE* IOMODEL, FILE* toolkitsoptionsfid, int in_solution_type, bool trace, IssmPDouble* X=NULL);
    8383                void Marshall(char** pmarshalled_data, int* pmarshalled_data_size, int marshall_direction);

  • issm/branches/trunk-larour-SLPS2020/src/c/cores/ProcessArguments.cpp

    r20939 r25596  
    88#include "../shared/shared.h"
    99
    10 void ProcessArguments(int* solution_type,char** pbinfilename,char** poutbinfilename,char** ptoolkitsfilename,char** plockfilename,char** prestartfilename, char** prootpath, int argc,char **argv){
     10void ProcessArguments(int* solution_type,char** pbinfilename,char** poutbinfilename,char** ptoolkitsfilename,char** plockfilename,char** prestartfilename, char** prootpath, char** pmodelname, int argc,char **argv){
    1111
    1212        /*Check input arguments*/
    … …  
    1818        *solution_type = StringToEnumx(argv[1]);
    1919        char* rootpatharg = argv[2];
    20         char* modelname   = argv[3];
     20        char* modelname      = xNew<char>(strlen(argv[3])+1);
     21        xMemCpy<char>(modelname,argv[3],strlen(argv[3])+1);
    2122
    2223        /*Recover myrank and length of string "my_rank" */
    … …  
    4243        *prestartfilename=restartfilename;
    4344        *prootpath=rootpath;
     45        *pmodelname=modelname;
    4446
    4547}

  • issm/branches/trunk-larour-SLPS2020/src/c/cores/cores.h

    r25304 r25596  
    7272
    7373//diverse
    74 void ProcessArguments(int* solution,char** pbinname,char** poutbinname,char** ptoolkitsname,char** plockname,char** prestartname, char** prootpath,int argc,char **argv);
     74void ProcessArguments(int* solution,char** pbinname,char** poutbinname,char** ptoolkitsname,char** plockname,char** prestartname, char** prootpath,char** pmodelname, int argc,char **argv);
    7575void WriteLockFile(char* filename);
    7676void ResetBoundaryConditions(FemModel* femmodel, int analysis_type);

  • issm/branches/trunk-larour-SLPS2020/src/c/main/issm_dakota.cpp

    r23066 r25596  
    44
    55#include "./issm.h"
     6#include <sys/stat.h>
    67
    78/*Dakota includes: */
    … …  
    1415#endif
    1516
    16 int main(int argc,char **argv){
     17/*prototypes:*/
     18int dirstructure(int argc,char** argv);
     19int issm_dakota_statistics(int argc,char** argv);
     20
     21int main(int argc,char **argv){ /*{{{*/
    1722
    1823        #if defined(_HAVE_DAKOTA_) && _DAKOTA_MAJOR_ >= 6
    … …  
    3843        dakota_error_file=xNew<char>((strlen(argv[2])+strlen(argv[3])+strlen(".qmu.err")+2));
    3944        sprintf(dakota_error_file,"%s/%s%s",argv[2],argv[3],".qmu.err");
     45
     46        /*Create directory structure for model outputs:*/
     47        dirstructure(argc,argv);
    4048
    4149        /* Parse input and construct Dakota LibraryEnvironment, performing input data checks*/
    … …  
    8391        env.execute();
    8492
     93        /* Run statistics if requested:*/
     94        issm_dakota_statistics(argc,argv);
     95
     96        /*free allocations:*/
    8597        xDelete<char>(dakota_input_file);
    8698        xDelete<char>(dakota_output_file);
    … …  
    94106        #endif
    95107
    96 }
     108} /*}}}*/
     109int dirstructure(int argc,char** argv){ /*{{{*/
     110
     111        char* input_file;
     112        FILE* fid;
     113        IoModel* iomodel=NULL;
     114        int check;
     115
     116        //qmu statistics
     117        bool statistics    = false;
     118        int  numdirectories = 0;
     119
     120        /*First things first, set the communicator as a global variable: */
     121        IssmComm::SetComm(MPI_COMM_WORLD);
     122
     123        /*Barrier:*/
     124        ISSM_MPI_Barrier(IssmComm::GetComm());
     125        _printf0_("Preparing directory structure for model outputs:" << "\n");
     126
     127        //open model input file for reading
     128        input_file=xNew<char>((strlen(argv[2])+strlen(argv[3])+strlen(".bin")+2));
     129        sprintf(input_file,"%s/%s%s",argv[2],argv[3],".bin");
     130        fid=fopen(input_file,"rb");
     131        if (fid==NULL) Cerr << "dirstructure error message: could not open model " << input_file << " to retrieve qmu statistics parameters" << std::endl;
     132
     133        //initialize IoModel, but light version, we just need it to fetch one constant:
     134        iomodel=new IoModel();
     135        iomodel->fid=fid;
     136        iomodel->FetchConstants();
     137
     138        //early return if statistics not requested:
     139        iomodel->FindConstant(&statistics,"md.qmu.statistics");
     140        if(!statistics){
     141                delete iomodel;
     142                fclose(fid);
     143                return 0;
     144        }
     145
     146        iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories");
     147
     148        /*Ok, we have everything we need to create the directory structure:*/
     149        if(IssmComm::GetRank()==0){
     150                for (int i=0;i<numdirectories;i++){
     151                        char directory[1000];
     152                        sprintf(directory,"./%i",i+1);
     153
     154                        check = mkdir(directory,ACCESSPERMS);
     155                        if (check) _error_("dirstructure error message: could not create directory " << directory << "\n");
     156                }
     157        }
     158
     159        //close model file:
     160        fclose(fid);
     161} /*}}}*/
     162int issm_dakota_statistics(int argc,char** argv){ /*{{{*/
     163
     164        char* input_file;
     165        FILE* fid;
     166        IoModel* iomodel=NULL;
     167        ISSM_MPI_Comm statcomm;
     168        int my_rank;
     169
     170        //qmu statistics
     171        bool statistics    = false;
     172        int  numstatistics = 0;
     173        int  numdirectories = 0;
     174        int  nfilesperdirectory = 0;
     175        char string[1000];
     176        char* name = NULL;
     177        char** fields = NULL;
     178        int    nfields;
     179        int*   steps=NULL;
     180        int    nsteps;
     181        int    nbins;
     182        int*   indices=NULL;
     183        int    nindices;
     184        int    nsamples;
     185        int    dummy;
     186        char*  directory=NULL;
     187        char*  model=NULL;
     188        Results* results=NULL;
     189        Parameters* parameters=NULL;
     190        int color;
     191
     192        /*First things first, set the communicator as a global variable: */
     193        IssmComm::SetComm(MPI_COMM_WORLD);
     194        my_rank=IssmComm::GetRank();
     195
     196        /*Barrier:*/
     197        ISSM_MPI_Barrier(IssmComm::GetComm());
     198        _printf0_("Dakota Statistic Computation" << "\n");
     199
     200        //open model input file for reading
     201        input_file=xNew<char>((strlen(argv[2])+strlen(argv[3])+strlen(".bin")+2));
     202        sprintf(input_file,"%s/%s%s",argv[2],argv[3],".bin");
     203        fid=fopen(input_file,"rb");
     204        if (fid==NULL) Cerr << "issm_dakota_statistics error message: could not open model " << input_file << " to retrieve qmu statistics parameters" << std::endl;
     205
     206        //initialize IoModel, but light version, we'll need it to fetch constants:
     207        iomodel=new IoModel();
     208        iomodel->fid=fid;
     209        iomodel->FetchConstants();
     210
     211        //early return if statistics not requested:
     212        iomodel->FindConstant(&statistics,"md.qmu.statistics");
     213        if(!statistics){
     214                delete iomodel;
     215                fclose(fid);
     216                return 0;
     217        }
     218
     219        //create parameters datasets with al the qmu statistics settings we need:
     220        if(statistics){
     221
     222                /*Initialize parameters and results:*/
     223                results   = new Results();
     224                parameters=new Parameters();
     225       
     226                //root  directory
     227                directory=xNew<char>(strlen(argv[2])+1);
     228                xMemCpy<char>(directory,argv[2],strlen(argv[2])+1);
     229                parameters->AddObject(new StringParam(DirectoryNameEnum,directory));
     230
     231                //model  name
     232                model=xNew<char>(strlen(argv[3])+1);
     233                xMemCpy<char>(model,argv[3],strlen(argv[3])+1);
     234                parameters->AddObject(new StringParam(InputFileNameEnum,model));
     235
     236                //nsamples
     237                iomodel->FindConstant(&nsamples,"md.qmu.method.params.samples");
     238                parameters->AddObject(new IntParam(QmuNsampleEnum,nsamples));
     239
     240                //ndirectories
     241                iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories");
     242                parameters->AddObject(new IntParam(QmuNdirectoriesEnum,numdirectories));
     243
     244                //nfiles per directory
     245                iomodel->FindConstant(&nfilesperdirectory,"md.qmu.statistics.nfiles_per_directory");
     246                parameters->AddObject(new IntParam(QmuNfilesPerDirectoryEnum,nfilesperdirectory));
     247
     248                //At this point, we don't want to go forward any longer, we want to create an MPI
     249                //communicator on which to carry out the computations:
     250                if ((my_rank+1)*nfilesperdirectory>nsamples)color=MPI_UNDEFINED;
     251                else color=0;
     252                ISSM_MPI_Comm_split(ISSM_MPI_COMM_WORLD,color, my_rank, &statcomm);
     253
     254
     255                iomodel->FindConstant(&numstatistics,"md.qmu.statistics.numstatistics");
     256                for (int i=1;i<=numstatistics;i++){
     257
     258                        char* directory=NULL;
     259                        char* model=NULL;
     260                        int   nsamples;
     261                        _printf0_("Dealing with qmu statistical computation #" << i << "\n");
     262               
     263                        sprintf(string,"md.qmu.statistics.method(%i).name",i);
     264                        iomodel->FindConstant(&name,string);
     265
     266                        sprintf(string,"md.qmu.statistics.method(%i).fields",i);
     267                        iomodel->FindConstant(&fields,&nfields,string);
     268                        parameters->AddObject(new StringArrayParam(FieldsEnum,fields,nfields));
     269
     270                        sprintf(string,"md.qmu.statistics.method(%i).steps",i);
     271                        iomodel->FetchData(&steps,&nsteps,&dummy,string);
     272                        parameters->AddObject(new IntVecParam(StepsEnum,steps,nsteps));
     273
     274                        if (strcmp(name,"Histogram")==0){
     275                                /*fetch nbins: */
     276                                sprintf(string,"md.qmu.statistics.method(%i).nbins",i);
     277                                iomodel->FindConstant(&nbins,string);
     278                                parameters->AddObject(new IntParam(NbinsEnum,nbins));
     279                                ComputeHistogram(parameters,results,color,statcomm);
     280                        }
     281                        else if (strcmp(name,"SampleSeries")==0){
     282                                /*fetch indices: */
     283                                sprintf(string,"md.qmu.statistics.method(%i).indices",i);
     284                                iomodel->FetchData(&indices,&nindices,&dummy,string);
     285                                parameters->AddObject(new IntVecParam(IndicesEnum,indices,nindices));
     286               
     287                                ComputeSampleSeries(parameters,results,color,statcomm);
     288                        }
     289                        else if (strcmp(name,"MeanVariance")==0){
     290                                ComputeMeanVariance(parameters,results,color,statcomm);
     291                        }
     292                        else _error_(" error creating qmu statistics methods parameters: unsupported method " << name);
     293                }
     294        }
     295        //close model file:
     296        fclose(fid);
     297
     298        /*output results:*/
     299        ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); _printf0_("Output file.\n");
     300        OutputStatistics(parameters,results);
     301
     302        /*Delete ressources:*/
     303        delete parameters;
     304        delete results;
     305
     306
     307       
     308} /*}}}*/

  • issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp

    r25566 r25596  
    44/*includes and prototypes:*/
    55#include "./issm.h"
    6 int ComputeMeanVariance(Parameters* parameters,Results* results);
    7 int ComputeSampleSeries(Parameters* parameters,Results* results);
    8 int readdata(IssmDouble** pdoublemat, int* pdoublematsize,
    9                      IssmDouble* pdouble, FILE* fid,char* field,int step);
    10 int OutputStatistics(Parameters* parameters,Results* results);
    11 int ComputeHistogram(Parameters* parameters,Results* results);
    126
    137int main(int argc,char **argv){ /*{{{*/
    … …  
    146140
    147141} /*}}}*/
    148 int ComputeMeanVariance(Parameters* parameters,Results* results){ /*{{{*/
    149 
    150         int nsamples;
    151         char* directory=NULL;
    152         char* model=NULL;
    153         char** fields=NULL;
    154         int* steps=NULL;
    155         int nsteps;
    156         int nfields;
    157         int range,lower_row,upper_row;
    158        
    159         /*intermediary:*/
    160         IssmDouble* doublemat=NULL;
    161         int         doublematsize;
    162         IssmDouble scalar;
    163 
    164         /*computation of average and variance itself:*/
    165         IssmDouble*  x = NULL;
    166         IssmDouble*  x2 = NULL;
    167         IssmDouble** xs = NULL;
    168         IssmDouble** xs2 = NULL;
    169         int*         xtype=NULL;
    170         int*         xsize=NULL;
    171 
    172         IssmDouble** meanx=NULL;
    173         IssmDouble** meanx2=NULL;
    174         int*         meantype=NULL;
    175         int*         meansize=NULL;
    176 
    177         /*Retrieve parameters:*/
    178         parameters->FindParam(&nsamples,QmuNsampleEnum);
    179         parameters->FindParam(&directory,DirectoryNameEnum);
    180         parameters->FindParam(&model,InputFileNameEnum);
    181         parameters->FindParam(&fields,&nfields,FieldsEnum);
    182         parameters->FindParam(&steps,&nsteps,StepsEnum);
    183 
    184         /*Get rank:*/
    185         int my_rank=IssmComm::GetRank();
    186 
    187         /*Open files and read them complelety, in a distributed way:*/
    188         range=DetermineLocalSize(nsamples,IssmComm::GetComm());
    189         GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
    190 
    191         /*Initialize arrays:*/
    192         xs=xNew<IssmDouble*>(nfields*nsteps);
    193         xs2=xNew<IssmDouble*>(nfields*nsteps);
    194         xtype=xNew<int>(nfields*nsteps);
    195         xsize=xNew<int>(nfields*nsteps);
    196 
    197         meantype=xNew<int>(nfields);
    198         meansize=xNew<int>(nfields);
    199         meanx=xNew<IssmDouble*>(nfields);
    200         meanx2=xNew<IssmDouble*>(nfields);
    201 
    202         /*Start opening files:*/
    203         for (int i=(lower_row+1);i<=upper_row;i++){
    204                 _printf0_("reading file #: " << i << "\n");
    205                 char file[1000];
    206                 long int  length;
    207                 char* buffer=NULL;
    208 
    209                 /*string:*/
    210                 sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
    211 
    212                 /*open file: */
    213                 _printf0_("    opening file: " << file << "\n");
    214                 FILE* fid=fopen(file,"rb");
    215                 if(fid==NULL) _error_("    could not open file: " << file << "\n");
    216 
    217                 /*figure out size of file, and read the whole thing:*/
    218                 _printf0_("    reading file:\n");
    219                 fseek (fid, 0, SEEK_END);
    220                 length = ftell (fid);
    221                 fseek (fid, 0, SEEK_SET);
    222                 buffer = xNew<char>(length);
    223                 fread (buffer, sizeof(char), length, fid);
    224 
    225                 /*close file:*/
    226                 fclose (fid);
    227 
    228                 /*create a memory stream with this buffer:*/
    229                 _printf0_("    processing file:\n");
    230                 fid=fmemopen(buffer, length, "rb");
    231 
    232                 /*start reading data from the buffer directly:*/
    233                 for (int f=0;f<nfields;f++){
    234                         char* field=fields[f];
    235                         fseek(fid,0,SEEK_SET);
    236                         for (int j=0;j<nsteps;j++){
    237                                 int counter=f*nsteps+j;
    238                                 xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    239                                 if(i==(lower_row+1)){
    240                                         if(xtype[counter]==1){
    241                                                 xs[counter]=xNew<IssmDouble>(1);
    242                                                 xs2[counter]=xNew<IssmDouble>(1);
    243                                                 *xs[counter]=scalar;
    244                                                 *xs2[counter]=pow(scalar,2.0);
    245                                                 xsize[counter]=1;
    246                                         }
    247                                         else if (xtype[counter]==3){
    248                                                 IssmDouble* doublemat2=xNew<IssmDouble>(doublematsize);
    249                                                 for(int k=0;k<doublematsize;k++)doublemat2[k]=pow(doublemat[k],2.0);
    250                                                 xs[counter]=doublemat;
    251                                                 xs2[counter]=doublemat2;
    252                                                 xsize[counter]=doublematsize;
    253                                         }
    254                                         else _error_("cannot carry out statistics on type " << xtype[counter]);
    255                                 }
    256                                 else{
    257                                         if(xtype[counter]==1){
    258                                                 *xs[counter]+=scalar;
    259                                                 *xs2[counter]+=pow(scalar,2.0);
    260                                         }
    261                                         else if (xtype[counter]==3){
    262                                                 IssmDouble* newdoublemat=xs[counter];
    263                                                 IssmDouble* newdoublemat2=xs2[counter];
    264                                                 for(int k=0;k<doublematsize;k++){
    265                                                         newdoublemat[k]+=doublemat[k];
    266                                                         newdoublemat2[k]+=pow(doublemat[k],2.0);
    267                                                 }
    268                                                 xs[counter]=newdoublemat;
    269                                                 xs2[counter]=newdoublemat2;
    270                                         }
    271                                         else _error_("cannot carry out statistics on type " << xtype[counter]);
    272                                 }
    273                         }
    274                 }
    275                
    276                 /*Deal with time mean: */
    277                 for (int f=0;f<nfields;f++){
    278                         char* field=fields[f];
    279                         fseek(fid,0,SEEK_SET);
    280                         meantype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
    281                         if(i==(lower_row+1)){
    282                                 if(meantype[f]==1){
    283                                         meanx[f]=xNewZeroInit<IssmDouble>(1);
    284                                         meanx2[f]=xNewZeroInit<IssmDouble>(1);
    285                                         meansize[f]=1;
    286                                 }
    287                                 else{
    288                                         meanx[f]=xNewZeroInit<IssmDouble>(doublematsize);
    289                                         meanx2[f]=xNewZeroInit<IssmDouble>(doublematsize);
    290                                         meansize[f]=doublematsize;
    291                                 }
    292                         }
    293                         fseek(fid,0,SEEK_SET);
    294                         if(meantype[f]==1){
    295                                 IssmDouble sc=0;
    296                                 IssmDouble sc2=0;
    297                                 for(int j=0;j<nsteps;j++){
    298                                         readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    299                                         sc+=scalar/nsteps;
    300                                 }
    301                                 sc2+=pow(sc,2.0);
    302                                 *meanx[f]+=sc;
    303                                 *meanx2[f]+=sc2;
    304                         }
    305                         else{
    306                                 IssmDouble* sc=meanx[f];
    307                                 IssmDouble* sc2=meanx2[f];
    308                                 IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
    309                                 IssmDouble* timemean2=xNewZeroInit<IssmDouble>(doublematsize);
    310 
    311                                 for(int j=0;j<nsteps;j++){
    312                                         readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    313                                         for (int k=0;k<doublematsize;k++){
    314                                                 timemean[k]+=doublemat[k]/nsteps;
    315                                         }
    316                                 }
    317                                 for (int k=0;k<doublematsize;k++){
    318                                         timemean2[k]=pow(timemean[k],2.0);
    319                                 }
    320                                 for (int k=0;k<doublematsize;k++){
    321                                         sc[k]+=timemean[k];
    322                                         sc2[k]+=timemean2[k];
    323                                 }
    324 
    325                         }
    326 
    327                 }
    328                 fclose(fid);
    329 
    330                 /*delete buffer:*/
    331                 xDelete<char>(buffer);
    332         }
    333         ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD);
    334         _printf0_("Done reading files, now computing mean and variance.\n");
    335 
    336         /*We have agregated x and x^2 across the cluster, now gather across the cluster onto
    337          *cpu0 and then compute statistics:*/
    338         for (int f=0;f<nfields;f++){
    339                 int counter0=f*nsteps+0;
    340                 if (xtype[counter0]==1){ /*deal with scalars {{{*/
    341                         IssmDouble mean,stddev;
    342                         for (int j=0;j<nsteps;j++){
    343                                 int counter=f*nsteps+j;
    344 
    345                                 /*we are broadcasting doubles:*/
    346                                 IssmDouble scalar=*xs[counter];
    347                                 IssmDouble scalar2=*xs2[counter];
    348                                 IssmDouble sumscalar;
    349                                 IssmDouble sumscalar2;
    350 
    351                                 ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    352                                 ISSM_MPI_Reduce(&scalar2,&sumscalar2,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    353                                 /*Build average and standard deviation. For standard deviation, use the
    354                                  *following formula: sigma^2=E(x^2)-mu^2:*/
    355                                 mean=sumscalar/nsamples;
    356                                 stddev=sqrt(sumscalar2/nsamples-pow(mean,2.0));
    357 
    358                                 /*add to results:*/
    359                                 if(my_rank==0){
    360                                         char fieldname[1000];
    361 
    362                                         sprintf(fieldname,"%s%s",fields[f],"Mean");
    363                                         results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,j+1,0));
    364                                         sprintf(fieldname,"%s%s",fields[f],"Stddev");
    365                                         results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,j+1,0));
    366                                 }
    367 
    368                         }
    369                 } /*}}}*/
    370                 else{ /*deal with arrays:{{{*/
    371 
    372                         int size=xsize[counter0];
    373 
    374                         IssmDouble*  mean=xNew<IssmDouble>(size);
    375                         IssmDouble*  stddev=xNew<IssmDouble>(size);
    376 
    377                         for (int j=0;j<nsteps;j++){
    378                                 int counter=f*nsteps+j;
    379 
    380                                 /*we are broadcasting double arrays:*/
    381                                 x=xs[counter];
    382                                 x2=xs2[counter];
    383 
    384                                 IssmDouble*  sumx=xNew<IssmDouble>(size);
    385                                 IssmDouble*  sumx2=xNew<IssmDouble>(size);
    386                                
    387                                 ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    388                                 ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    389 
    390                                 /*Build average and standard deviation. For standard deviation, use the
    391                                  *following formula: sigma^2=E(x^2)-mu^2:*/
    392                                 for (int k=0;k<size;k++){
    393                                         mean[k]=sumx[k]/nsamples;
    394                                         stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
    395                                 }
    396 
    397                                 /*add to results:*/
    398                                 if(my_rank==0){
    399                                         char fieldname[1000];
    400 
    401                                         sprintf(fieldname,"%s%s",fields[f],"Mean");
    402                                         results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,j+1,0));
    403                                         sprintf(fieldname,"%s%s",fields[f],"Stddev");
    404                                         results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,j+1,0));
    405                                 }
    406                         }
    407                 } /*}}}*/
    408         }
    409         /*Do the same but for the time mean:*/
    410         for (int f=0;f<nfields;f++){
    411                 if (meantype[f]==1){ /*deal with scalars {{{*/
    412                         IssmDouble mean,stddev;
    413 
    414                         /*we are broadcasting doubles:*/
    415                         IssmDouble scalar=*meanx[f];
    416                         IssmDouble scalar2=*meanx2[f];
    417                         IssmDouble sumscalar;
    418                         IssmDouble sumscalar2;
    419 
    420                         ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    421                         ISSM_MPI_Reduce(&scalar2,&sumscalar2,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    422                         /*Build average and standard deviation. For standard deviation, use the
    423                          *following formula: sigma^2=E(x^2)-mu^2:*/
    424                         mean=sumscalar/nsamples;
    425                         stddev=sqrt(sumscalar2/nsamples-pow(mean,2.0));
    426 
    427                         /*add to results:*/
    428                         if(my_rank==0){
    429                                 char fieldname[1000];
    430 
    431                                 sprintf(fieldname,"%s%s",fields[f],"TimeMean");
    432                                 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,1,0));
    433                                 sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
    434                                 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,1,0));
    435                         }
    436                 } /*}}}*/
    437                 else{ /*deal with arrays:{{{*/
    438 
    439                         int size=meansize[f];
    440                         IssmDouble*  mean=xNew<IssmDouble>(size);
    441                         IssmDouble*  stddev=xNew<IssmDouble>(size);
    442 
    443                         /*we are broadcasting double arrays:*/
    444                         x=meanx[f];
    445                         x2=meanx2[f];
    446 
    447                         IssmDouble*  sumx=xNew<IssmDouble>(size);
    448                         IssmDouble*  sumx2=xNew<IssmDouble>(size);
    449                                
    450                         ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    451                         ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    452 
    453                         /*Build average and standard deviation. For standard deviation, use the
    454                          *following formula: sigma^2=E(x^2)-mu^2:*/
    455                         for (int k=0;k<size;k++){
    456                                 mean[k]=sumx[k]/nsamples;
    457                                 stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
    458                         }
    459 
    460                         /*add to results:*/
    461                         if(my_rank==0){
    462                                 char fieldname[1000];
    463 
    464                                 sprintf(fieldname,"%s%s",fields[f],"TimeMean");
    465                                 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,1,0));
    466                                 sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
    467                                 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,1,0));
    468                         }
    469                 } /*}}}*/
    470         }
    471 
    472 
    473 } /*}}}*/
    474 int ComputeSampleSeries(Parameters* parameters,Results* results){ /*{{{*/
    475 
    476         int nsamples;
    477         char* directory=NULL;
    478         char* model=NULL;
    479         char** fields=NULL;
    480         int* steps=NULL;
    481         int nsteps;
    482         int nfields;
    483         int range,lower_row,upper_row;
    484         int* indices=NULL;
    485         int  nindices;
    486        
    487         /*intermediary:*/
    488         IssmDouble* doublemat=NULL;
    489         int         doublematsize;
    490         IssmDouble scalar;
    491 
    492         /*computation of average and variance itself:*/
    493         IssmDouble*  x = NULL;
    494         IssmDouble*  allx=NULL;
    495         IssmDouble** xs = NULL;
    496         int*         xtype=NULL;
    497         int*         xsize=NULL;
    498 
    499         /*Retrieve parameters:*/
    500         parameters->FindParam(&nsamples,QmuNsampleEnum);
    501         parameters->FindParam(&directory,DirectoryNameEnum);
    502         parameters->FindParam(&model,InputFileNameEnum);
    503         parameters->FindParam(&fields,&nfields,FieldsEnum);
    504         parameters->FindParam(&steps,&nsteps,StepsEnum);
    505         parameters->FindParam(&indices,&nindices,IndicesEnum);
    506 
    507         /*Get rank:*/
    508         int my_rank=IssmComm::GetRank();
    509 
    510         /*Open files and read them complelety, in a distributed way:*/
    511         range=DetermineLocalSize(nsamples,IssmComm::GetComm());
    512         GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
    513 
    514         /*Initialize arrays:*/
    515         xs=xNew<IssmDouble*>(nfields*nsteps);
    516         xtype=xNew<int>(nfields*nsteps);
    517         xsize=xNew<int>(nfields*nsteps);
    518 
    519         /*Start opening files:*/
    520         for (int i=(lower_row+1);i<=upper_row;i++){
    521                 _printf0_("reading file #: " << i << "\n");
    522                 char file[1000];
    523                 long int  length;
    524                 char* buffer=NULL;
    525 
    526                 /*string:*/
    527                 sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
    528 
    529                 /*open file: */
    530                 _printf0_("    opening file:\n");
    531                 FILE* fid=fopen(file,"rb");
    532 
    533                 /*figure out size of file, and read the whole thing:*/
    534                 _printf0_("    reading file:\n");
    535                 fseek (fid, 0, SEEK_END);
    536                 length = ftell (fid);
    537                 fseek (fid, 0, SEEK_SET);
    538                 buffer = xNew<char>(length);
    539                 fread (buffer, sizeof(char), length, fid);
    540 
    541                 /*close file:*/
    542                 fclose (fid);
    543 
    544                 /*create a memory stream with this buffer:*/
    545                 _printf0_("    processing file:\n");
    546                 fid=fmemopen(buffer, length, "rb");
    547 
    548                 /*start reading data from the buffer directly:*/
    549                 for (int f=0;f<nfields;f++){
    550                         fseek(fid,0,SEEK_SET);
    551                         char* field=fields[f];
    552                         for (int j=0;j<nsteps;j++){
    553                                 int counter=f*nsteps+j;
    554                                 xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    555                                 if(i==(lower_row+1)){
    556                                         if(xtype[counter]==1){
    557                                                 x=xNew<IssmDouble>(range);
    558                                                 x[0]=scalar;
    559                                                 xs[counter]=x;
    560                                                 xsize[counter]=range;
    561                                         }
    562                                         else if (xtype[counter]==3){
    563                                                 x=xNew<IssmDouble>(nindices*range);
    564                                                 for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
    565                                                 xs[counter]=x;
    566                                                 xsize[counter]=range*nindices;
    567                                         }
    568                                         else _error_("cannot carry out statistics on type " << xtype[counter]);
    569                                 }
    570                                 else{
    571                                         if(xtype[counter]==1){
    572                                                 x=xs[counter];
    573                                                 x[i-(lower_row+1)]=scalar;
    574                                                 xs[counter]=x;
    575                                         }
    576                                         else if (xtype[counter]==3){
    577                                                 x=xs[counter];
    578                                                 for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
    579                                                 xs[counter]=x;
    580                                         }
    581                                         else _error_("cannot carry out statistics on type " << xtype[counter]);
    582                                 }
    583                         }
    584                 }
    585                 fclose(fid);
    586 
    587                 /*delete buffer:*/
    588                 xDelete<char>(buffer);
    589         }
    590         ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD);
    591         _printf0_("Done reading files, now assembling time series.\n");
    592 
    593         for (int f=0;f<nfields;f++){
    594                 for (int j=0;j<nsteps;j++){
    595                         int counter=f*nsteps+j;
    596                         if (xtype[counter]==1){
    597                                 /*we are broadcasting range times doubles:*/
    598                                 x=xs[counter];
    599                                 allx=xNew<IssmDouble>(nsamples);
    600                                 MPI_Gather(x, range, ISSM_MPI_PDOUBLE,allx, range, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
    601                                 /*add to results:*/
    602                                 if(my_rank==0){
    603                                         char fieldname[1000];
    604                                        
    605                                         sprintf(fieldname,"%s%s",fields[f],"Samples");
    606                                         results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,1,j+1,0));
    607                                 }
    608                         }
    609                         else{
    610                                 /*we are broadcasting double arrays:*/
    611                                 x=xs[counter];
    612                                 allx=xNew<IssmDouble>(nsamples*nindices);
    613                                
    614                                 MPI_Gather(x, range*nindices, ISSM_MPI_PDOUBLE,allx, range*nindices, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
    615                                
    616                                 /*add to results:*/
    617                                 if(my_rank==0){
    618                                         char fieldname[1000];
    619                                         sprintf(fieldname,"%s%s",fields[f],"Samples");
    620                                         results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,nindices,j+1,0));
    621                                 }
    622                         }
    623                 }
    624         }
    625 
    626 
    627 } /*}}}*/
    628 int ComputeHistogram(Parameters* parameters,Results* results){ /*{{{*/
    629 
    630         int nsamples;
    631         char* directory=NULL;
    632         char* model=NULL;
    633         char** fields=NULL;
    634         int* steps=NULL;
    635         int nsteps;
    636         int nfields;
    637         int nbins;
    638         int range,lower_row,upper_row;
    639        
    640         /*intermediary:*/
    641         IssmDouble* doublemat=NULL;
    642         int         doublematsize;
    643         IssmDouble scalar;
    644 
    645         /*computation of average and variance itself:*/
    646         IssmDouble** maxxs = NULL;
    647         IssmDouble** minxs = NULL;
    648         int*         xtype=NULL;
    649         int*         xsize=NULL;
    650 
    651         IssmDouble** maxmeans=NULL;
    652         IssmDouble** minmeans=NULL;
    653         int*         meanxtype=NULL;
    654         int*         meanxsize=NULL;
    655 
    656         /*Retrieve parameters:*/
    657         parameters->FindParam(&nsamples,QmuNsampleEnum);
    658         parameters->FindParam(&directory,DirectoryNameEnum);
    659         parameters->FindParam(&model,InputFileNameEnum);
    660         parameters->FindParam(&fields,&nfields,FieldsEnum);
    661         parameters->FindParam(&steps,&nsteps,StepsEnum);
    662         parameters->FindParam(&nbins,NbinsEnum);
    663 
    664         /*Get rank:*/
    665         int my_rank=IssmComm::GetRank();
    666 
    667         /*Open files and read them complelety, in a distributed way:*/
    668         range=DetermineLocalSize(nsamples,IssmComm::GetComm());
    669         GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
    670 
    671         /*Initialize arrays:*/
    672 
    673         maxmeans=xNew<IssmDouble*>(nfields);
    674         minmeans=xNew<IssmDouble*>(nfields);
    675         meanxtype=xNew<int>(nfields);
    676         meanxsize=xNew<int>(nfields);
    677 
    678         maxxs=xNew<IssmDouble*>(nfields*nsteps);
    679         minxs=xNew<IssmDouble*>(nfields*nsteps);
    680         xtype=xNew<int>(nfields*nsteps);
    681         xsize=xNew<int>(nfields*nsteps);
    682 
    683         /*Start opening files:*/
    684         for (int i=(lower_row+1);i<=upper_row;i++){
    685                 _printf0_("reading file #: " << i << "\n");
    686                 char file[1000];
    687                 long int  length;
    688                 char* buffer=NULL;
    689 
    690                 /*string:*/
    691                 sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
    692 
    693                 /*open file: */
    694                 _printf0_("    opening file: " << file << "\n");
    695                 FILE* fid=fopen(file,"rb");
    696                 if(fid==NULL)_error_("cound not open file: " << file << "\n");
    697 
    698                 /*figure out size of file, and read the whole thing:*/
    699                 _printf0_("    reading file:\n");
    700                 fseek (fid, 0, SEEK_END);
    701                 length = ftell (fid);
    702                 fseek (fid, 0, SEEK_SET);
    703                 buffer = xNew<char>(length);
    704                 fread (buffer, sizeof(char), length, fid);
    705 
    706                 /*close file:*/
    707                 fclose (fid);
    708 
    709                 /*create a memory stream with this buffer:*/
    710                 _printf0_("    processing file:\n");
    711                 fid=fmemopen(buffer, length, "rb");
    712 
    713                 /*start reading data from the buffer directly:*/
    714                 for (int f=0;f<nfields;f++){
    715                         char* field=fields[f];
    716                         fseek(fid,0,SEEK_SET);
    717                         for (int j=0;j<nsteps;j++){
    718                                 int counter=f*nsteps+j;
    719                                 xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    720                                 if(i==(lower_row+1)){
    721                                         if(xtype[counter]==1){
    722                                                 maxxs[counter]=xNew<IssmDouble>(1);
    723                                                 minxs[counter]=xNew<IssmDouble>(1);
    724                                                 *maxxs[counter]=scalar;
    725                                                 *minxs[counter]=scalar;
    726                                                 xsize[counter]=1;
    727                                         }
    728                                         else if (xtype[counter]==3){
    729                                                 maxxs[counter]=xNew<IssmDouble>(doublematsize);
    730                                                 xMemCpy<IssmDouble>(maxxs[counter],doublemat,doublematsize);
    731                                                 minxs[counter]=xNew<IssmDouble>(doublematsize);
    732                                                 xMemCpy<IssmDouble>(minxs[counter],doublemat,doublematsize);
    733                                                 xsize[counter]=doublematsize;
    734                                                 xDelete<IssmDouble>(doublemat);
    735                                         }
    736                                         else _error_("cannot carry out statistics on type " << xtype[counter]);
    737                                 }
    738                                 else{
    739                                         if(xtype[counter]==1){
    740                                                 *maxxs[counter]=max(*maxxs[counter],scalar);
    741                                                 *minxs[counter]=min(*minxs[counter],scalar);
    742                                         }
    743                                         else if (xtype[counter]==3){
    744                                                 IssmDouble* newmax=maxxs[counter];
    745                                                 IssmDouble* newmin=minxs[counter];
    746                                                 for(int k=0;k<doublematsize;k++){
    747                                                         if(doublemat[k]>newmax[k])newmax[k]=doublemat[k];
    748                                                         if(doublemat[k]<newmin[k])newmin[k]=doublemat[k];
    749                                                 }
    750                                                 xDelete<IssmDouble>(doublemat);
    751                                         }
    752                                         else _error_("cannot carry out statistics on type " << xtype[counter]);
    753                                 }
    754                         }
    755                 }
    756                 _printf0_("    average in time:\n");
    757 
    758                 /*Deal with average in time: */
    759                 for (int f=0;f<nfields;f++){
    760                         fseek(fid,0,SEEK_SET);
    761                         char* field=fields[f];
    762                         meanxtype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
    763                        
    764                         if(meanxtype[f]==1){
    765                                 meanxsize[f]=1;
    766                                 IssmDouble timemean=0;
    767                                 fseek(fid,0,SEEK_SET);
    768                                 for (int j=0;j<nsteps;j++){
    769                                         readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    770                                         timemean+=scalar/nsteps;
    771                                 }
    772                                        
    773                                 /*Figure out max and min of time means: */
    774                                 if(i==(lower_row+1)){
    775                                         maxmeans[f]=xNewZeroInit<IssmDouble>(1);
    776                                         minmeans[f]=xNewZeroInit<IssmDouble>(1);
    777                                         *maxmeans[f]=timemean;
    778                                         *minmeans[f]=timemean;
    779                                 }
    780                                 else{
    781                                         *maxmeans[f]=max(*maxmeans[f],timemean);
    782                                         *minmeans[f]=min(*minmeans[f],timemean);
    783                                 }
    784                         }
    785                         else{
    786                                 meanxsize[f]=doublematsize;
    787                                 fseek(fid,0,SEEK_SET);
    788                                 IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
    789                                 for (int j=0;j<nsteps;j++){
    790                                         readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    791                                         for (int k=0;k<doublematsize;k++){
    792                                                 timemean[k]+=doublemat[k]/nsteps;
    793                                         }
    794                                         xDelete<IssmDouble>(doublemat);
    795                                 }
    796 
    797                                 if(i==(lower_row+1)){
    798                                         maxmeans[f]=xNew<IssmDouble>(doublematsize);
    799                                         xMemCpy<IssmDouble>(maxmeans[f],timemean,doublematsize);
    800                                         minmeans[f]=xNew<IssmDouble>(doublematsize);
    801                                         xMemCpy<IssmDouble>(minmeans[f],timemean,doublematsize);
    802                                 }
    803                                 else{
    804                                         IssmDouble* maxx=maxmeans[f];
    805                                         IssmDouble* minx=minmeans[f];
    806 
    807                                         for(int k=0;k<doublematsize;k++){
    808                                                 maxx[k]=max(maxx[k],timemean[k]);
    809                                                 minx[k]=min(minx[k],timemean[k]);
    810                                         }
    811                                         maxmeans[f]=maxx;
    812                                         minmeans[f]=minx;
    813                                 }
    814                         }
    815                 }
    816                 fclose(fid);
    817 
    818                 /*delete buffer:*/
    819                 xDelete<char>(buffer);
    820         }
    821         ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD);
    822         _printf0_("Done reading files, now computing min and max.\n");
    823 
    824         /*We have agregated minx and max across the cluster, now gather across the cluster onto
    825          *cpu0 and then compute statistics:*/
    826         for (int f=0;f<nfields;f++){
    827                 int counter0=f*nsteps+0;
    828                 if (xtype[counter0]==1){ /*deal with scalars {{{*/
    829                         for (int j=0;j<nsteps;j++){
    830                                 int counter=f*nsteps+j;
    831 
    832                                 /*we are broadcasting doubles:*/
    833                                 IssmDouble maxscalar=*maxxs[counter];
    834                                 IssmDouble minscalar=*minxs[counter];
    835                                 IssmDouble allmaxscalar;
    836                                 IssmDouble allminscalar;
    837                                 IssmDouble sumscalar_alltimes=0;
    838 
    839                                 ISSM_MPI_Allreduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
    840                                 ISSM_MPI_Allreduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
    841 
    842                                 /*Store broadcasted value for later computation of histograms:*/
    843                                 *maxxs[counter]=allmaxscalar;
    844                                 *minxs[counter]=allminscalar;
    845 
    846                         }
    847                 } /*}}}*/
    848                 else{ /*deal with arrays:{{{*/
    849 
    850                         int size=xsize[counter0];
    851                         for (int j=0;j<nsteps;j++){
    852                                 int counter=f*nsteps+j;
    853 
    854                                 /*we are broadcasting double arrays:*/
    855                                 IssmDouble* maxx=maxxs[counter];
    856                                 IssmDouble* minx=minxs[counter];
    857 
    858                                 IssmDouble*  allmax=xNew<IssmDouble>(size);
    859                                 IssmDouble*  allmin=xNew<IssmDouble>(size);
    860                                
    861                                 ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
    862                                 ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
    863 
    864                                 /*Store broadcasted value for later computation of histograms:*/
    865                                 maxxs[counter]=allmax;
    866                                 minxs[counter]=allmin;
    867                         }
    868                 } /*}}}*/
    869         }
    870        
    871         /*Now do the same for the time mean fields:*/
    872         for (int f=0;f<nfields;f++){
    873                 if (meanxtype[f]==1){ /*deal with scalars {{{*/
    874 
    875                         /*we are broadcasting doubles:*/
    876                         IssmDouble maxscalar=*maxmeans[f];
    877                         IssmDouble minscalar=*minmeans[f];
    878                         IssmDouble allmaxscalar;
    879                         IssmDouble allminscalar;
    880 
    881                         ISSM_MPI_Allreduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
    882                         ISSM_MPI_Allreduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
    883 
    884                         /*Store for later use in histogram computation:*/
    885                         *maxmeans[f]=allmaxscalar;
    886                         *minmeans[f]=allminscalar;
    887 
    888                 } /*}}}*/
    889                 else{ /*deal with arrays:{{{*/
    890 
    891                         int size=meanxsize[f];
    892 
    893                         /*we are broadcasting double arrays:*/
    894                         IssmDouble* maxx=maxmeans[f];
    895                         IssmDouble* minx=minmeans[f];
    896 
    897                         IssmDouble*  allmax=xNew<IssmDouble>(size);
    898                         IssmDouble*  allmin=xNew<IssmDouble>(size);
    899 
    900                         ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
    901                         ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
    902 
    903                        
    904                         /*Store for later use in histogram computation:*/
    905                         maxmeans[f]=allmax;
    906                         minmeans[f]=allmin;
    907 
    908                 } /*}}}*/
    909         }
    910 
    911         /*Now that we have the min and max, we can start binning. First allocate
    912          * histograms, then start filling them:*/
    913         IssmDouble** histogram=xNew<IssmDouble*>(nfields*nsteps);
    914         IssmDouble** timehistogram=xNew<IssmDouble*>(nfields);
    915 
    916         _printf0_("Start reading files again, this time binning values in the histogram:\n");
    917         /*Start opening files:*/
    918         for (int i=(lower_row+1);i<=upper_row;i++){
    919                 _printf0_("reading file #: " << i << "\n");
    920                 char file[1000];
    921                 long int  length;
    922                 char* buffer=NULL;
    923 
    924                 /*string:*/
    925                 sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
    926 
    927                 /*open file: */
    928                 _printf0_("    opening file:\n");
    929                 FILE* fid=fopen(file,"rb");
    930                 if(fid==NULL)_error_("cound not open file: " << file << "\n");
    931 
    932                 /*figure out size of file, and read the whole thing:*/
    933                 _printf0_("    reading file:\n");
    934                 fseek (fid, 0, SEEK_END);
    935                 length = ftell (fid);
    936                 fseek (fid, 0, SEEK_SET);
    937                 buffer = xNew<char>(length);
    938                 fread (buffer, sizeof(char), length, fid);
    939 
    940                 /*close file:*/
    941                 fclose (fid);
    942 
    943                 /*create a memory stream with this buffer:*/
    944                 _printf0_("    processing file:\n");
    945                 fid=fmemopen(buffer, length, "rb");
    946 
    947                 /*start reading data from the buffer directly:*/
    948                 for (int f=0;f<nfields;f++){
    949                         char* field=fields[f];
    950                         fseek(fid,0,SEEK_SET);
    951                         for (int j=0;j<nsteps;j++){
    952                                 int counter=f*nsteps+j;
    953                                 xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    954                                 if(i==(lower_row+1)){
    955                                         if(xtype[counter]==1){
    956                                                 IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(nbins);
    957                                                 IssmDouble ma=*maxxs[counter];
    958                                                 IssmDouble mi=*minxs[counter];
    959                                                 int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index--;
    960                                                 localhistogram[index]++;
    961                                                 histogram[counter]=localhistogram;
    962                                         }
    963                                         else if (xtype[counter]==3){
    964                                                 IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(doublematsize*nbins);
    965                                                 IssmDouble* ma=maxxs[counter];
    966                                                 IssmDouble* mi=minxs[counter];
    967                                                 for (int k=0;k<doublematsize;k++){
    968                                                         IssmDouble scalar=doublemat[k];
    969                                                         int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index--;
    970                                                         _assert_(scalar<=ma[k]); _assert_(scalar>=mi[k]); _assert_(index<nbins);
    971                                                         localhistogram[k*nbins+index]++;
    972                                                 }
    973                                                 histogram[counter]=localhistogram;
    974                                                 xDelete<IssmDouble>(doublemat);
    975                                         }
    976                                         else _error_("cannot carry out statistics on type " << xtype[counter]);
    977                                 }
    978                                 else{
    979                                         if(xtype[counter]==1){
    980                                                 IssmDouble* localhistogram=histogram[counter];
    981                                                 IssmDouble ma=*maxxs[counter];
    982                                                 IssmDouble mi=*minxs[counter];
    983                                                 int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
    984                                                 localhistogram[index]++;
    985                                         }
    986                                         else if (xtype[counter]==3){
    987                                                 IssmDouble* localhistogram=histogram[counter];
    988                                                 IssmDouble* ma=maxxs[counter];
    989                                                 IssmDouble* mi=minxs[counter];
    990                                                 for (int k=0;k<doublematsize;k++){
    991                                                         IssmDouble scalar=doublemat[k];
    992                                                         int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
    993                                                         localhistogram[k*nbins+index]++;
    994                                                 }
    995                                                 xDelete<IssmDouble>(doublemat);
    996                                         }
    997                                         else _error_("cannot carry out statistics on type " << xtype[counter]);
    998                                 }
    999                         }
    1000                 }
    1001                 _printf0_("    average in time:\n");
    1002 
    1003                 /*Deal with average in time: */
    1004                 for (int f=0;f<nfields;f++){
    1005                         fseek(fid,0,SEEK_SET);
    1006                         char* field=fields[f];
    1007                         meanxtype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
    1008                        
    1009                         if(meanxtype[f]==1){
    1010                                 IssmDouble timemean=0;
    1011                                 fseek(fid,0,SEEK_SET);
    1012                                 for (int j=0;j<nsteps;j++){
    1013                                         readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    1014                                         timemean+=scalar/nsteps;
    1015                                 }
    1016                                        
    1017                                 /*Figure out max and min of time means: */
    1018                                 if(i==(lower_row+1)){
    1019                                         IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(nbins);
    1020                                         IssmDouble ma=*maxmeans[f];
    1021                                         IssmDouble mi=*minmeans[f];
    1022                                         int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
    1023                                         localhistogram[index]++;
    1024                                         timehistogram[f]=localhistogram;
    1025                                 }
    1026                                 else{
    1027                                         IssmDouble* localhistogram=timehistogram[f];
    1028                                         IssmDouble ma=*maxmeans[f];
    1029                                         IssmDouble mi=*minmeans[f];
    1030                                         int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
    1031                                         localhistogram[index]++;
    1032                                 }
    1033                         }
    1034                         else{
    1035                                 fseek(fid,0,SEEK_SET);
    1036                                 IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
    1037                                 for (int j=0;j<nsteps;j++){
    1038                                         readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    1039                                         for (int k=0;k<doublematsize;k++){
    1040                                                 timemean[k]+=doublemat[k]/nsteps;
    1041                                         }
    1042                                         xDelete<IssmDouble>(doublemat);
    1043                                 }
    1044 
    1045                                 if(i==(lower_row+1)){
    1046                                         IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(doublematsize*nbins);
    1047                                         IssmDouble* ma=maxmeans[f];
    1048                                         IssmDouble* mi=minmeans[f];
    1049                                        
    1050                                         for (int k=0;k<doublematsize;k++){
    1051                                                 IssmDouble scalar=timemean[k];
    1052                                                 int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
    1053                                                 localhistogram[k*nbins+index]++;
    1054                                         }
    1055                                         timehistogram[f]=localhistogram;
    1056                                 }
    1057                                 else{
    1058                                        
    1059                                         IssmDouble* localhistogram=timehistogram[f];
    1060                                         IssmDouble* ma=maxmeans[f];
    1061                                         IssmDouble* mi=minmeans[f];
    1062                                        
    1063                                         for (int k=0;k<doublematsize;k++){
    1064                                                 IssmDouble scalar=timemean[k];
    1065                                                 int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
    1066                                                 localhistogram[k*nbins+index]++;
    1067                                         }
    1068                                 }
    1069                         }
    1070                 }
    1071                 fclose(fid);
    1072 
    1073                 /*delete buffer:*/
    1074                 xDelete<char>(buffer);
    1075         }
    1076         _printf0_("Start aggregating histogram:\n");
    1077 
    1078         /*We have agregated histograms across the cluster, now gather them across  the cluster onto
    1079          *cpu0: */
    1080         for (int f=0;f<nfields;f++){
    1081                 int counter0=f*nsteps+0;
    1082                 if (xtype[counter0]==1){ /*deal with scalars {{{*/
    1083                         for (int j=0;j<nsteps;j++){
    1084                                 int counter=f*nsteps+j;
    1085 
    1086                                 /*we are broadcasting doubles:*/
    1087                                 IssmDouble* histo=histogram[counter]; //size nbins
    1088                                 IssmDouble* allhisto=xNewZeroInit<IssmDouble>(nbins);
    1089 
    1090                                 ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
    1091 
    1092                                 /*add to results:*/
    1093                                 if(my_rank==0){
    1094                                         char fieldname[1000];
    1095 
    1096                                         sprintf(fieldname,"%s%s",fields[f],"Histogram");
    1097                                         results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,1,nbins,j+1,0));
    1098 
    1099                                         sprintf(fieldname,"%s%s",fields[f],"Max");
    1100                                         results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*maxxs[counter],j+1,0));
    1101                                         sprintf(fieldname,"%s%s",fields[f],"Min");
    1102                                         results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*minxs[counter],j+1,0));
    1103                                 }
    1104                         }
    1105                 } /*}}}*/
    1106                 else{ /*deal with arrays:{{{*/
    1107 
    1108                         int size=xsize[counter0];
    1109                         for (int j=0;j<nsteps;j++){
    1110                                 int counter=f*nsteps+j;
    1111 
    1112                                 /*we are broadcasting double arrays:*/
    1113                                 IssmDouble* histo=histogram[counter];
    1114                                 IssmDouble* allhisto=xNew<IssmDouble>(size*nbins);
    1115                                
    1116                                 ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
    1117                                 xDelete<IssmDouble>(histo);
    1118 
    1119                                 /*add to results:*/
    1120                                 if(my_rank==0){
    1121                                         char fieldname[1000];
    1122 
    1123                                         sprintf(fieldname,"%s%s",fields[f],"Histogram");
    1124                                         results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,size,nbins,j+1,0));
    1125 
    1126                                         sprintf(fieldname,"%s%s",fields[f],"Max");
    1127                                         results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,maxxs[counter],size,1,j+1,0));
    1128                                         sprintf(fieldname,"%s%s",fields[f],"Min");
    1129                                         results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,minxs[counter],size,1,j+1,0));
    1130                                 }
    1131                         }
    1132                 } /*}}}*/
    1133         }
    1134         _printf0_("Start aggregating time mean histogram:\n");
    1135        
    1136         /*Now do the same for the time mean fields:*/
    1137         for (int f=0;f<nfields;f++){
    1138                 if (meanxtype[f]==1){ /*deal with scalars {{{*/
    1139 
    1140                         /*we are broadcasting doubles:*/
    1141                         IssmDouble* histo=timehistogram[f];
    1142                         IssmDouble* allhisto=xNewZeroInit<IssmDouble>(nbins);
    1143 
    1144                         ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
    1145 
    1146                         /*add to results at time step 1:*/
    1147                         if(my_rank==0){
    1148                                 char fieldname[1000];
    1149 
    1150                                 sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram");
    1151                                 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,1,nbins,1,0));
    1152 
    1153                                 sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax");
    1154                                 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*maxmeans[f],1,0));
    1155                                 sprintf(fieldname,"%s%s",fields[f],"TimeMeaMin");
    1156                                 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*minmeans[f],1,0));
    1157                         }
    1158                 } /*}}}*/
    1159                 else{ /*deal with arrays:{{{*/
    1160 
    1161                         int size=meanxsize[f];
    1162 
    1163                         /*we are broadcasting double arrays:*/
    1164                         IssmDouble* histo=timehistogram[f];
    1165                         IssmDouble* allhisto=xNewZeroInit<IssmDouble>(size*nbins);
    1166 
    1167                         ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
    1168                         xDelete<IssmDouble>(histo);
    1169                         /*add to results at step 1:*/
    1170                         if(my_rank==0){
    1171                                 char fieldname[1000];
    1172 
    1173                                 sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram");
    1174                                 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,size,nbins,1,0));
    1175                                 sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax");
    1176                                 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,maxmeans[f],size,1,1,0));
    1177                                 sprintf(fieldname,"%s%s",fields[f],"TimeMeanMin");
    1178                                 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,minmeans[f],size,1,1,0));
    1179                         }
    1180                 } /*}}}*/
    1181         }
    1182 } /*}}}*/
    1183 int OutputStatistics(Parameters* parameters,Results* results){ /*{{{*/
    1184        
    1185         char   outputfilename[1000];
    1186         char* directory=NULL;
    1187         char* model=NULL;
    1188         char* method=NULL;
    1189         int   nsamples;
    1190         int* steps=NULL;
    1191         int nsteps;
    1192 
    1193         FemModel* femmodel=new FemModel();
    1194        
    1195         /*Some parameters that will allow us to use the OutputResultsx module:*/
    1196         parameters->AddObject(new BoolParam(QmuIsdakotaEnum,false));
    1197         parameters->AddObject(new BoolParam(SettingsIoGatherEnum,true));
    1198 
    1199         parameters->FindParam(&directory,DirectoryNameEnum);
    1200         parameters->FindParam(&model,InputFileNameEnum);
    1201         parameters->FindParam(&method,AnalysisTypeEnum);
    1202         parameters->FindParam(&nsamples,QmuNsampleEnum);
    1203         parameters->FindParam(&steps,&nsteps,StepsEnum);
    1204 
    1205         sprintf(outputfilename,"%s/%s-%s.stats-samp%i-time%i:%i",directory,model,method,nsamples,steps[0],steps[nsteps-1]);
    1206         parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
    1207 
    1208         /*Call OutputResults module:*/
    1209         femmodel->parameters=parameters;
    1210         femmodel->results=results;
    1211 
    1212         OutputResultsx(femmodel);
    1213 } /*}}}*/
    1214 int readdata(IssmDouble** pdoublemat, int* pdoublematsize, IssmDouble* pdouble, FILE* fid,char* field,int step){ /*{{{*/
    1215 
    1216         int length;
    1217         char fieldname[1000];
    1218         int   fieldname_size;
    1219         IssmDouble   rtime;
    1220         int          rstep;
    1221         int M,N;
    1222 
    1223         //fields that we retrive:
    1224         IssmDouble  dfield;
    1225         char*       sfield    = NULL;
    1226         IssmDouble* dmatfield = NULL;
    1227         int*        imatfield = NULL;
    1228                        
    1229         //type of the returned field:
    1230         int type;
    1231         int found=0;
    1232 
    1233         while(1){
    1234 
    1235                 size_t ret_code = fread(&fieldname_size, sizeof(int), 1, fid);
    1236                 if(ret_code != 1) break; //we are done.
    1237 
    1238                 fread(fieldname, sizeof(char), fieldname_size, fid);
    1239                 //_printf0_("fieldname: " << fieldname << "\n");
    1240 
    1241                 fread(&rtime, sizeof(IssmDouble), 1, fid);
    1242                 fread(&rstep, sizeof(int), 1, fid);
    1243 
    1244                 //check on field:
    1245                 if ((step==rstep) && (strcmp(field,fieldname)==0)){
    1246 
    1247                         //ok, go read the result really:
    1248                         fread(&type,sizeof(int),1,fid);
    1249                         fread(&M,sizeof(int),1,fid);
    1250                         if (type==1){
    1251                                 fread(&dfield,sizeof(IssmDouble),1,fid);
    1252                         }
    1253                         else if (type==2){
    1254                                 fread(&M,sizeof(int),1,fid);
    1255                                 sfield=xNew<char>(M);
    1256                                 fread(sfield,sizeof(char),M,fid);
    1257                         }
    1258                         else if (type==3){
    1259                                 fread(&N,sizeof(int),1,fid);
    1260                                 dmatfield=xNew<IssmDouble>(M*N);
    1261                                 fread(dmatfield,sizeof(IssmDouble),M*N,fid);
    1262                         }
    1263                         else if (type==4){
    1264                                 fread(&N,sizeof(int),1,fid);
    1265                                 imatfield=xNew<int>(M*N);
    1266                                 fread(imatfield,sizeof(int),M*N,fid);
    1267                         }
    1268                         else _error_("cannot read data of type " << type << "\n");
    1269                         found=1;
    1270                         break;
    1271                 }
    1272                 else{
    1273                         //just skim to next results.
    1274                         fread(&type,sizeof(int),1,fid);
    1275                         fread(&M,sizeof(int),1,fid);
    1276                         if (type==1){
    1277                                 fseek(fid,sizeof(IssmDouble),SEEK_CUR);
    1278                         }
    1279                         else if(type==2){
    1280                                 fseek(fid,M*sizeof(char),SEEK_CUR);
    1281                         }
    1282                         else if(type==3){
    1283                                 fread(&N,sizeof(int),1,fid);
    1284                                 fseek(fid,M*N*sizeof(IssmDouble),SEEK_CUR);
    1285                         }
    1286                         else if(type==4){
    1287                                 fread(&N,sizeof(int),1,fid);
    1288                                 fseek(fid,M*N*sizeof(int),SEEK_CUR);
    1289                         }
    1290                         else _error_("cannot read data of type " << type << "\n");
    1291                 }
    1292         }
    1293         if(found==0)_error_("cound not find " << field << " at step " << step  << "\n");
    1294 
    1295         /*assign output pointers:*/
    1296         *pdoublemat=dmatfield;
    1297         *pdoublematsize=M*N;
    1298         *pdouble=dfield;
    1299 
    1300         /*return:*/
    1301         return type;
    1302 
    1303 }
    1304 /*}}}*/

  • issm/branches/trunk-larour-SLPS2020/src/c/modules/ModelProcessorx/Dakota/CreateParametersDakota.cpp

    r25257 r25596  
    4141        int M,N;
    4242
     43        //qmu statistics
     44        bool statistics    = false;
     45        int  numdirectories = 0;
     46        int  nfilesperdirectory = 0;
     47
    4348        /*recover parameters: */
    4449        iomodel->FindConstant(&dakota_analysis,"md.qmu.isdakota");
    … …  
    7277                /*Ok, we have all the response descriptors. Build a parameter with it: */
    7378                parameters->AddObject(new StringArrayParam(QmuResponsedescriptorsEnum,responsedescriptors,numresponsedescriptors));
     79
     80                /*Deal with statistics: */
     81                iomodel->FindConstant(&statistics,"md.qmu.statistics");
     82                parameters->AddObject(new BoolParam(QmuStatisticsEnum,statistics));
     83                if(statistics){
     84                        iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories");
     85                        parameters->AddObject(new IntParam(QmuNdirectoriesEnum,numdirectories));
     86
     87                        iomodel->FindConstant(&nfilesperdirectory,"md.qmu.statistics.nfiles_per_directory");
     88                        parameters->AddObject(new IntParam(QmuNfilesPerDirectoryEnum,nfilesperdirectory));
     89                }
    7490
    7591                /*Load partitioning vectors specific to variables:*/
    … …  
    109125                /*}}}*/
    110126
    111 
    112127                /*Deal with data needed because of qmu variables*/
    113128                DataSet* dataset_variable_descriptors = new DataSet(QmuVariableDescriptorsEnum);
    … …  
    143158                delete dataset_variable_descriptors;
    144159
    145                 /*clean-up*/
     160                /*clean-up {{{*/
    146161                for(i=0;i<numresponsedescriptors;i++){
    147162                        descriptor=responsedescriptors[i];
    … …  
    157172                xDelete<char>(qmuerrname);
    158173                xDelete<char>(qmuoutname);
    159 
    160                
    161                
     174                xDelete<char>(name);
     175                /*}}}*/
    162176        }
    163177
    164         /*Free data*/
    165         xDelete<char>(name);
    166178}

  • issm/branches/trunk-larour-SLPS2020/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp

    r14999 r25596  
    1919
    2020        int i,j,n;
    21         const int maxels=50;
     21        const int maxels=100;
    2222        const int width=maxels+1;
    2323

  • issm/branches/trunk-larour-SLPS2020/src/c/modules/OutputResultsx/OutputResultsx.cpp

    r25335 r25596  
    6363                else{
    6464                        if(my_rank==0){
     65                                /*a little bit complicated. Either statistic computations are requested, which means we
     66                                 * put our outbin files in subidirectories with numbers, or we don't, and we dump our
     67                                 * outbins directly in the current directory:*/
    6568                                int currEvalId ;
     69                                int nfilesperdirectory;
     70                                bool statistics=false;
     71                                char* root=NULL;
     72                                char* modelname=NULL;
     73                               
    6674                                femmodel->parameters->FindParam(&currEvalId,QmuCurrEvalIdEnum);
    67                                 sprintf(outputfilename2,"%s.%i",outputfilename,currEvalId);
     75                                femmodel->parameters->FindParam(&statistics,QmuStatisticsEnum);
     76                                femmodel->parameters->FindParam(&nfilesperdirectory,QmuNfilesPerDirectoryEnum);
     77
     78                                if(statistics){
     79                                        femmodel->parameters->FindParam(&root,RootPathEnum);
     80                                        femmodel->parameters->FindParam(&modelname,ModelnameEnum);
     81                                        sprintf(outputfilename2,"%s/%i/%s.outbin.%i",root,(int)(floor((currEvalId-1)/nfilesperdirectory)+1),modelname,currEvalId);
     82                                }
     83                                else{
     84                                        sprintf(outputfilename2,"%s.%i",outputfilename,currEvalId);
     85                                }
    6886                                fid=pfopen0(outputfilename2,"ab+");
    6987                        }

  • issm/branches/trunk-larour-SLPS2020/src/c/modules/modules.h

    r23992 r25596  
    8181#include "./PointCloudFindNeighborsx/PointCloudFindNeighborsx.h"
    8282#include "./PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h"
     83#include "./QmuStatisticsx/QmuStatisticsx.h"
    8384#include "./Reduceloadx/Reduceloadx.h"
    8485#include "./Reducevectorgtofx/Reducevectorgtofx.h"

  • issm/branches/trunk-larour-SLPS2020/src/c/shared/Enum/Enum.vim

    r25535 r25596  
    314314syn keyword cConstant QmuResponsePartitionsEnum
    315315syn keyword cConstant QmuResponsePartitionsNpartEnum
     316syn keyword cConstant QmuStatisticsEnum
     317syn keyword cConstant QmuNumstatisticsEnum
     318syn keyword cConstant QmuNdirectoriesEnum
     319syn keyword cConstant QmuNfilesPerDirectoryEnum
     320syn keyword cConstant QmuStatisticsMethodEnum
     321syn keyword cConstant QmuMethodsEnum
    316322syn keyword cConstant RestartFileNameEnum
    317323syn keyword cConstant ResultsEnum
    318324syn keyword cConstant RootPathEnum
     325syn keyword cConstant ModelnameEnum
    319326syn keyword cConstant SaveResultsEnum
    320327syn keyword cConstant SolidearthPlanetRadiusEnum
    … …  
    14581465syn keyword cType PowerVariogram
    14591466syn keyword cType Profiler
     1467syn keyword cType QmuStatisticsMethod
    14601468syn keyword cType Quadtree
    14611469syn keyword cType Radar

  • issm/branches/trunk-larour-SLPS2020/src/c/shared/Enum/EnumDefinitions.h

    r25535 r25596  
    308308        QmuResponsePartitionsEnum,
    309309        QmuResponsePartitionsNpartEnum,
     310        QmuStatisticsEnum,
     311        QmuNumstatisticsEnum,
     312        QmuNdirectoriesEnum,
     313        QmuNfilesPerDirectoryEnum,
     314        QmuStatisticsMethodEnum,
     315        QmuMethodsEnum,
    310316        RestartFileNameEnum,
    311317        ResultsEnum,
    312318        RootPathEnum,
     319        ModelnameEnum,
    313320        SaveResultsEnum,
    314321        SolidearthPlanetRadiusEnum,

  • issm/branches/trunk-larour-SLPS2020/src/c/shared/Enum/EnumToStringx.cpp

    r25535 r25596  
    316316                case QmuResponsePartitionsEnum : return "QmuResponsePartitions";
    317317                case QmuResponsePartitionsNpartEnum : return "QmuResponsePartitionsNpart";
     318                case QmuStatisticsEnum : return "QmuStatistics";
     319                case QmuNumstatisticsEnum : return "QmuNumstatistics";
     320                case QmuNdirectoriesEnum : return "QmuNdirectories";
     321                case QmuNfilesPerDirectoryEnum : return "QmuNfilesPerDirectory";
     322                case QmuStatisticsMethodEnum : return "QmuStatisticsMethod";
     323                case QmuMethodsEnum : return "QmuMethods";
    318324                case RestartFileNameEnum : return "RestartFileName";
    319325                case ResultsEnum : return "Results";
    320326                case RootPathEnum : return "RootPath";
     327                case ModelnameEnum : return "Modelname";
    321328                case SaveResultsEnum : return "SaveResults";
    322329                case SolidearthPlanetRadiusEnum : return "SolidearthPlanetRadius";

  • issm/branches/trunk-larour-SLPS2020/src/c/shared/Enum/StringToEnumx.cpp

    r25535 r25596  
    322322              else if (strcmp(name,"QmuResponsePartitions")==0) return QmuResponsePartitionsEnum;
    323323              else if (strcmp(name,"QmuResponsePartitionsNpart")==0) return QmuResponsePartitionsNpartEnum;
     324              else if (strcmp(name,"QmuStatistics")==0) return QmuStatisticsEnum;
     325              else if (strcmp(name,"QmuNumstatistics")==0) return QmuNumstatisticsEnum;
     326              else if (strcmp(name,"QmuNdirectories")==0) return QmuNdirectoriesEnum;
     327              else if (strcmp(name,"QmuNfilesPerDirectory")==0) return QmuNfilesPerDirectoryEnum;
     328              else if (strcmp(name,"QmuStatisticsMethod")==0) return QmuStatisticsMethodEnum;
     329              else if (strcmp(name,"QmuMethods")==0) return QmuMethodsEnum;
    324330              else if (strcmp(name,"RestartFileName")==0) return RestartFileNameEnum;
    325331              else if (strcmp(name,"Results")==0) return ResultsEnum;
    326332              else if (strcmp(name,"RootPath")==0) return RootPathEnum;
     333              else if (strcmp(name,"Modelname")==0) return ModelnameEnum;
    327334              else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
    328335              else if (strcmp(name,"SolidearthPlanetRadius")==0) return SolidearthPlanetRadiusEnum;
    … …  
    376383              else if (strcmp(name,"SmbDsnowIdx")==0) return SmbDsnowIdxEnum;
    377384              else if (strcmp(name,"SmbCldFrac")==0) return SmbCldFracEnum;
    378               else if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum;
     385         else stage=4;
     386   }
     387   if(stage==4){
     388              if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum;
    379389              else if (strcmp(name,"SmbDelta18oSurface")==0) return SmbDelta18oSurfaceEnum;
    380390              else if (strcmp(name,"SmbDenIdx")==0) return SmbDenIdxEnum;
    … …  
    383393              else if (strcmp(name,"SmbF")==0) return SmbFEnum;
    384394              else if (strcmp(name,"SmbInitDensityScaling")==0) return SmbInitDensityScalingEnum;
    385          else stage=4;
    386    }
    387    if(stage==4){
    388               if (strcmp(name,"SmbIsaccumulation")==0) return SmbIsaccumulationEnum;
     395              else if (strcmp(name,"SmbIsaccumulation")==0) return SmbIsaccumulationEnum;
    389396              else if (strcmp(name,"SmbIsalbedo")==0) return SmbIsalbedoEnum;
    390397              else if (strcmp(name,"SmbIsclimatology")==0) return SmbIsclimatologyEnum;
    … …  
    499506              else if (strcmp(name,"BalancethicknessOmega")==0) return BalancethicknessOmegaEnum;
    500507              else if (strcmp(name,"BalancethicknessThickeningRate")==0) return BalancethicknessThickeningRateEnum;
    501               else if (strcmp(name,"BasalCrevasse")==0) return BasalCrevasseEnum;
     508         else stage=5;
     509   }
     510   if(stage==5){
     511              if (strcmp(name,"BasalCrevasse")==0) return BasalCrevasseEnum;
    502512              else if (strcmp(name,"BasalforcingsFloatingiceMeltingRate")==0) return BasalforcingsFloatingiceMeltingRateEnum;
    503513              else if (strcmp(name,"BasalforcingsGeothermalflux")==0) return BasalforcingsGeothermalfluxEnum;
    … …  
    506516              else if (strcmp(name,"BasalforcingsIsmip6BasinId")==0) return BasalforcingsIsmip6BasinIdEnum;
    507517              else if (strcmp(name,"BasalforcingsIsmip6Tf")==0) return BasalforcingsIsmip6TfEnum;
    508          else stage=5;
    509    }
    510    if(stage==5){
    511               if (strcmp(name,"BasalforcingsIsmip6TfShelf")==0) return BasalforcingsIsmip6TfShelfEnum;
     518              else if (strcmp(name,"BasalforcingsIsmip6TfShelf")==0) return BasalforcingsIsmip6TfShelfEnum;
    512519              else if (strcmp(name,"BasalforcingsIsmip6MeltAnomaly")==0) return BasalforcingsIsmip6MeltAnomalyEnum;
    513520              else if (strcmp(name,"BasalforcingsOceanSalinity")==0) return BasalforcingsOceanSalinityEnum;
    … …  
    622629              else if (strcmp(name,"UGia")==0) return UGiaEnum;
    623630              else if (strcmp(name,"UGiaRate")==0) return UGiaRateEnum;
    624               else if (strcmp(name,"Gradient")==0) return GradientEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"Gradient")==0) return GradientEnum;
    625635              else if (strcmp(name,"GroundinglineHeight")==0) return GroundinglineHeightEnum;
    626636              else if (strcmp(name,"HydraulicPotential")==0) return HydraulicPotentialEnum;
    … …  
    629639              else if (strcmp(name,"HydrologyBumpHeight")==0) return HydrologyBumpHeightEnum;
    630640              else if (strcmp(name,"HydrologyBumpSpacing")==0) return HydrologyBumpSpacingEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"HydrologydcBasalMoulinInput")==0) return HydrologydcBasalMoulinInputEnum;
     641              else if (strcmp(name,"HydrologydcBasalMoulinInput")==0) return HydrologydcBasalMoulinInputEnum;
    635642              else if (strcmp(name,"HydrologydcEplThickness")==0) return HydrologydcEplThicknessEnum;
    636643              else if (strcmp(name,"HydrologydcEplThicknessOld")==0) return HydrologydcEplThicknessOldEnum;
    … …  
    745752              else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
    746753              else if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum;
    747               else if (strcmp(name,"SmbBMin")==0) return SmbBMinEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"SmbBMin")==0) return SmbBMinEnum;
    748758              else if (strcmp(name,"SmbBNeg")==0) return SmbBNegEnum;
    749759              else if (strcmp(name,"SmbBPos")==0) return SmbBPosEnum;
    … …  
    752762              else if (strcmp(name,"SmbDailyairdensity")==0) return SmbDailyairdensityEnum;
    753763              else if (strcmp(name,"SmbDailyairhumidity")==0) return SmbDailyairhumidityEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"SmbDailydlradiation")==0) return SmbDailydlradiationEnum;
     764              else if (strcmp(name,"SmbDailydlradiation")==0) return SmbDailydlradiationEnum;
    758765              else if (strcmp(name,"SmbDailydsradiation")==0) return SmbDailydsradiationEnum;
    759766              else if (strcmp(name,"SmbDailypressure")==0) return SmbDailypressureEnum;
    … …  
    868875              else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum;
    869876              else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum;
    870               else if (strcmp(name,"Temperature")==0) return TemperatureEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"Temperature")==0) return TemperatureEnum;
    871881              else if (strcmp(name,"TemperaturePDD")==0) return TemperaturePDDEnum;
    872882              else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
    … …  
    875885              else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum;
    876886              else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum;
     887              else if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum;
    881888              else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum;
    882889              else if (strcmp(name,"Thickness")==0) return ThicknessEnum;
    … …  
    991998              else if (strcmp(name,"Outputdefinition87")==0) return Outputdefinition87Enum;
    992999              else if (strcmp(name,"Outputdefinition88")==0) return Outputdefinition88Enum;
    993               else if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum;
    9941004              else if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum;
    9951005              else if (strcmp(name,"Outputdefinition90")==0) return Outputdefinition90Enum;
    … …  
    9981008              else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum;
    9991009              else if (strcmp(name,"Outputdefinition94")==0) return Outputdefinition94Enum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"Outputdefinition95")==0) return Outputdefinition95Enum;
     1010              else if (strcmp(name,"Outputdefinition95")==0) return Outputdefinition95Enum;
    10041011              else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum;
    10051012              else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum;
    … …  
    11141121              else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
    11151122              else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
    1116               else if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum;
    11171127              else if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum;
    11181128              else if (strcmp(name,"Fset")==0) return FsetEnum;
    … …  
    11211131              else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum;
    11221132              else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum;
     1133              else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum;
    11271134              else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
    11281135              else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
    … …  
    12371244              else if (strcmp(name,"Moulin")==0) return MoulinEnum;
    12381245              else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;
    1239               else if (strcmp(name,"Mpi")==0) return MpiEnum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"Mpi")==0) return MpiEnum;
    12401250              else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum;
    12411251              else if (strcmp(name,"Mumps")==0) return MumpsEnum;
    … …  
    12441254              else if (strcmp(name,"Nodal")==0) return NodalEnum;
    12451255              else if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
     1256              else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
    12501257              else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
    12511258              else if (strcmp(name,"None")==0) return NoneEnum;
    … …  
    13601367              else if (strcmp(name,"TotalGroundedBmbScaled")==0) return TotalGroundedBmbScaledEnum;
    13611368              else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum;
    1362               else if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum;
     1369         else stage=12;
     1370   }
     1371   if(stage==12){
     1372              if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum;
    13631373              else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum;
    13641374              else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
    … …  
    13671377              else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
    13681378              else if (strcmp(name,"Tria")==0) return TriaEnum;
    1369          else stage=12;
    1370    }
    1371    if(stage==12){
    1372               if (strcmp(name,"TriaInput")==0) return TriaInputEnum;
     1379              else if (strcmp(name,"TriaInput")==0) return TriaInputEnum;
    13731380              else if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum;
    13741381              else if (strcmp(name,"VectorParam")==0) return VectorParamEnum;
Note: See TracChangeset for help on using the changeset viewer.