Changeset 27268


Ignore:
Timestamp:
09/07/22 19:46:13 (3 years ago)
Author:
Eric.Larour
Message:

CHG: better optimization for memory, fixed leaks.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/branches/trunk-larour-SLPS2022/src/c/modules/QmuStatisticsx/QmuStatisticsx.cpp

    r27267 r27268  
    66#include "../OutputResultsx/OutputResultsx.h"
    77
    8 int readdata(IssmDouble** pdoublemat, int* pdoublematsize, IssmDouble* pdouble, FILE* fid,char* field,int step){ /*{{{*/
    9 
    10         int length;
    11         char fieldname[1000];
    12         int   fieldname_size;
    13         IssmDouble   rtime;
    14         int          rstep;
    15         int M,N;
    16 
    17         //fields that we retrive:
    18         IssmDouble  dfield;
    19         char*       sfield    = NULL;
    20         IssmDouble* dmatfield = NULL;
    21         int*        imatfield = NULL;
    22 
    23         //type of the returned field:
    24         int type;
    25         int found=0;
    26 
    27         while(1){
    28 
    29                 size_t ret_code = fread(&fieldname_size, sizeof(int), 1, fid);
    30                 if(ret_code != 1) break; //we are done.
    31 
    32                 fread(fieldname, sizeof(char), fieldname_size, fid);
    33                 //_printf0_("fieldname: " << fieldname << "\n");
    34 
    35                 fread(&rtime, sizeof(IssmDouble), 1, fid);
    36                 fread(&rstep, sizeof(int), 1, fid);
    37 
    38                 //check on field:
    39                 if ((step==rstep) && (strcmp(field,fieldname)==0)){
    40 
    41                         //ok, go read the result really:
    42                         fread(&type,sizeof(int),1,fid);
    43                         fread(&M,sizeof(int),1,fid);
    44                         if (type==1){
    45                                 fread(&dfield,sizeof(IssmDouble),1,fid);
    46                         }
    47                         else if (type==2){
    48                                 fread(&M,sizeof(int),1,fid);
    49                                 sfield=xNew<char>(M);
    50                                 fread(sfield,sizeof(char),M,fid);
    51                         }
    52                         else if (type==3){
    53                                 fread(&N,sizeof(int),1,fid);
    54                                 dmatfield=xNew<IssmDouble>(M*N);
    55                                 fread(dmatfield,sizeof(IssmDouble),M*N,fid);
    56                         }
    57                         else if (type==4){
    58                                 fread(&N,sizeof(int),1,fid);
    59                                 imatfield=xNew<int>(M*N);
    60                                 fread(imatfield,sizeof(int),M*N,fid);
    61                         }
    62                         else _error_("cannot read data of type " << type << "\n");
    63                         found=1;
    64                         break;
    65                 }
    66                 else{
    67                         //just skim to next results.
    68                         fread(&type,sizeof(int),1,fid);
    69                         fread(&M,sizeof(int),1,fid);
    70                         if (type==1){
    71                                 fseek(fid,sizeof(IssmDouble),SEEK_CUR);
    72                         }
    73                         else if(type==2){
    74                                 fseek(fid,M*sizeof(char),SEEK_CUR);
    75                         }
    76                         else if(type==3){
    77                                 fread(&N,sizeof(int),1,fid);
    78                                 fseek(fid,M*N*sizeof(IssmDouble),SEEK_CUR);
    79                         }
    80                         else if(type==4){
    81                                 fread(&N,sizeof(int),1,fid);
    82                                 fseek(fid,M*N*sizeof(int),SEEK_CUR);
    83                         }
    84                         else _error_("cannot read data of type " << type << "\n");
    85                 }
    86         }
    87         if(found==0)_error_("cound not find " << field << " at step " << step  << "\n");
    88 
    89         /*assign output pointers:*/
    90         *pdoublemat=dmatfield;
    91         *pdoublematsize=M*N;
    92         *pdouble=dfield;
    93 
    94         /*return:*/
    95         return type;
    96 
    97 }
    98 /*}}}*/
    99 int ComputeHistogram(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){  /*{{{*/
    100 
    101         int nsamples;
    102         char* directory=NULL;
    103         char* model=NULL;
    104         char** fields=NULL;
    105         int* steps=NULL;
    106         int nsteps;
    107         int nfields;
    108         int nbins;
    109         int range,lower_row,upper_row;
    110         int nfilesperdirectory;
    111 
    112         /*intermediary:*/
    113         IssmDouble* doublemat=NULL;
    114         int         doublematsize;
    115         IssmDouble scalar;
    116 
    117         /*computation of average and variance itself:*/
    118         IssmDouble** maxxs = NULL;
    119         IssmDouble** minxs = NULL;
    120         int*         xtype=NULL;
    121         int*         xsize=NULL;
    122 
    123         IssmDouble** maxmeans=NULL;
    124         IssmDouble** minmeans=NULL;
    125         int*         meanxtype=NULL;
    126         int*         meanxsize=NULL;
    127 
    128         /*only work on the statistical communicator: */
    129         if (color==MPI_UNDEFINED)return 0;
    130 
    131         /*Retrieve parameters:*/
    132         parameters->FindParam(&nfilesperdirectory,QmuNfilesPerDirectoryEnum);
    133         parameters->FindParam(&nsamples,QmuNsampleEnum);
    134         parameters->FindParam(&directory,DirectoryNameEnum);
    135         parameters->FindParam(&model,InputFileNameEnum);
    136         parameters->FindParam(&fields,&nfields,FieldsEnum);
    137         parameters->FindParam(&steps,&nsteps,StepsEnum);
    138         parameters->FindParam(&nbins,NbinsEnum);
    139 
    140         /*Get rank from the stat comm communicator:*/
    141         IssmComm::SetComm(statcomm);
    142         int my_rank=IssmComm::GetRank();
    143 
    144         /*Open files and read them complelety, in a distributed way:*/
    145         range=DetermineLocalSize(nsamples,IssmComm::GetComm());
    146         GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
    147 
    148         /*Initialize arrays:*/
    149         maxmeans=xNew<IssmDouble*>(nfields);
    150         minmeans=xNew<IssmDouble*>(nfields);
    151         meanxtype=xNew<int>(nfields);
    152         meanxsize=xNew<int>(nfields);
    153 
    154         maxxs=xNew<IssmDouble*>(nfields*nsteps);
    155         minxs=xNew<IssmDouble*>(nfields*nsteps);
    156         xtype=xNew<int>(nfields*nsteps);
    157         xsize=xNew<int>(nfields*nsteps);
    158 
    159         /*Start opening files:*/
    160         for(int i=(lower_row+1);i<=upper_row;i++){
    161                 _printf0_("reading file #: " << i << "\n");
    162                 char file[1000];
    163                 long int  length;
    164                 char* buffer=NULL;
    165 
    166                 /*string:*/
    167                 sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
    168 
    169                 /*open file: */
    170                 _printf0_("    opening file: " << file << "\n");
    171                 FILE* fid=fopen(file,"rb");
    172                 if(fid==NULL)_error_("cound not open file: " << file << "\n");
    173 
    174                 /*figure out size of file, and read the whole thing:*/
    175                 _printf0_("    reading file:\n");
    176                 fseek(fid, 0, SEEK_END);
    177                 length = ftell (fid);
    178                 fseek(fid, 0, SEEK_SET);
    179                 buffer = xNew<char>(length);
    180                 fread(buffer, sizeof(char), length, fid);
    181 
    182                 /*close file:*/
    183                 fclose(fid);
    184 
    185                 /*create a memory stream with this buffer:*/
    186                 _printf0_("    processing file:\n");
    187                 fid=fmemopen(buffer, length, "rb");
    188 
    189                 /*start reading data from the buffer directly:*/
    190                 for (int f=0;f<nfields;f++){
    191                         char* field=fields[f];
    192                         fseek(fid,0,SEEK_SET);
    193                         for (int j=0;j<nsteps;j++){
    194                                 int counter=f*nsteps+j;
    195                                 xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    196                                 if(i==(lower_row+1)){
    197                                         if(xtype[counter]==1){
    198                                                 maxxs[counter]=xNew<IssmDouble>(1);
    199                                                 minxs[counter]=xNew<IssmDouble>(1);
    200                                                 *maxxs[counter]=scalar;
    201                                                 *minxs[counter]=scalar;
    202                                                 xsize[counter]=1;
    203                                         }
    204                                         else if (xtype[counter]==3){
    205                                                 maxxs[counter]=xNew<IssmDouble>(doublematsize);
    206                                                 xMemCpy<IssmDouble>(maxxs[counter],doublemat,doublematsize);
    207                                                 minxs[counter]=xNew<IssmDouble>(doublematsize);
    208                                                 xMemCpy<IssmDouble>(minxs[counter],doublemat,doublematsize);
    209                                                 xsize[counter]=doublematsize;
    210                                                 xDelete<IssmDouble>(doublemat);
    211                                         }
    212                                         else _error_("cannot carry out statistics on type " << xtype[counter]);
    213                                 }
    214                                 else{
    215                                         if(xtype[counter]==1){
    216                                                 *maxxs[counter]=max(*maxxs[counter],scalar);
    217                                                 *minxs[counter]=min(*minxs[counter],scalar);
    218                                         }
    219                                         else if (xtype[counter]==3){
    220                                                 IssmDouble* newmax=maxxs[counter];
    221                                                 IssmDouble* newmin=minxs[counter];
    222                                                 for(int k=0;k<doublematsize;k++){
    223                                                         if(doublemat[k]>newmax[k])newmax[k]=doublemat[k];
    224                                                         if(doublemat[k]<newmin[k])newmin[k]=doublemat[k];
    225                                                 }
    226                                                 xDelete<IssmDouble>(doublemat);
    227                                         }
    228                                         else _error_("cannot carry out statistics on type " << xtype[counter]);
    229                                 }
    230                         }
    231                 }
    232                 _printf0_("    average in time:\n");
    233 
    234                 /*Deal with average in time: */
    235                 for (int f=0;f<nfields;f++){
    236                         fseek(fid,0,SEEK_SET);
    237                         char* field=fields[f];
    238                         meanxtype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
    239 
    240                         if(meanxtype[f]==1){
    241                                 meanxsize[f]=1;
    242                                 IssmDouble timemean=0;
    243                                 fseek(fid,0,SEEK_SET);
    244                                 for (int j=0;j<nsteps;j++){
    245                                         readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    246                                         timemean+=scalar/nsteps;
    247                                 }
    248 
    249                                 /*Figure out max and min of time means: */
    250                                 if(i==(lower_row+1)){
    251                                         maxmeans[f]=xNewZeroInit<IssmDouble>(1);
    252                                         minmeans[f]=xNewZeroInit<IssmDouble>(1);
    253                                         *maxmeans[f]=timemean;
    254                                         *minmeans[f]=timemean;
    255                                 }
    256                                 else{
    257                                         *maxmeans[f]=max(*maxmeans[f],timemean);
    258                                         *minmeans[f]=min(*minmeans[f],timemean);
    259                                 }
    260                         }
    261                         else{
    262                                 meanxsize[f]=doublematsize;
    263                                 fseek(fid,0,SEEK_SET);
    264                                 IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
    265                                 for (int j=0;j<nsteps;j++){
    266                                         readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    267                                         for (int k=0;k<doublematsize;k++){
    268                                                 timemean[k]+=doublemat[k]/nsteps;
    269                                         }
    270                                         xDelete<IssmDouble>(doublemat);
    271                                 }
    272 
    273                                 if(i==(lower_row+1)){
    274                                         maxmeans[f]=xNew<IssmDouble>(doublematsize);
    275                                         xMemCpy<IssmDouble>(maxmeans[f],timemean,doublematsize);
    276                                         minmeans[f]=xNew<IssmDouble>(doublematsize);
    277                                         xMemCpy<IssmDouble>(minmeans[f],timemean,doublematsize);
    278                                 }
    279                                 else{
    280                                         IssmDouble* maxx=maxmeans[f];
    281                                         IssmDouble* minx=minmeans[f];
    282 
    283                                         for(int k=0;k<doublematsize;k++){
    284                                                 maxx[k]=max(maxx[k],timemean[k]);
    285                                                 minx[k]=min(minx[k],timemean[k]);
    286                                         }
    287                                         maxmeans[f]=maxx;
    288                                         minmeans[f]=minx;
    289                                 }
    290                         }
    291                 }
    292                 fclose(fid);
    293 
    294                 /*delete buffer:*/
    295                 xDelete<char>(buffer);
    296         }
    297         ISSM_MPI_Barrier(IssmComm::GetComm());
    298         _printf0_("Done reading files, now computing min and max.\n");
    299 
    300         /*We have agregated minx and max across the cluster, now gather across the cluster onto
    301          *cpu0 and then compute statistics:*/
    302         for (int f=0;f<nfields;f++){
    303                 int counter0=f*nsteps+0;
    304                 if (xtype[counter0]==1){ /*deal with scalars {{{*/
    305                         for (int j=0;j<nsteps;j++){
    306                                 int counter=f*nsteps+j;
    307 
    308                                 /*we are broadcasting doubles:*/
    309                                 IssmDouble maxscalar=*maxxs[counter];
    310                                 IssmDouble minscalar=*minxs[counter];
    311                                 IssmDouble allmaxscalar;
    312                                 IssmDouble allminscalar;
    313                                 IssmDouble sumscalar_alltimes=0;
    314 
    315                                 ISSM_MPI_Allreduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
    316                                 ISSM_MPI_Allreduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
    317 
    318                                 /*Store broadcasted value for later computation of histograms:*/
    319                                 *maxxs[counter]=allmaxscalar;
    320                                 *minxs[counter]=allminscalar;
    321 
    322                         }
    323                 } /*}}}*/
    324                 else{ /*deal with arrays:{{{*/
    325 
    326                         int size=xsize[counter0];
    327                         for (int j=0;j<nsteps;j++){
    328                                 int counter=f*nsteps+j;
    329 
    330                                 /*we are broadcasting double arrays:*/
    331                                 IssmDouble* maxx=maxxs[counter];
    332                                 IssmDouble* minx=minxs[counter];
    333 
    334                                 IssmDouble*  allmax=xNew<IssmDouble>(size);
    335                                 IssmDouble*  allmin=xNew<IssmDouble>(size);
    336 
    337                                 ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
    338                                 ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
    339 
    340                                 /*Store broadcasted value for later computation of histograms:*/
    341                                 maxxs[counter]=allmax;
    342                                 minxs[counter]=allmin;
    343                         }
    344                 } /*}}}*/
    345         }
    346 
    347         /*Now do the same for the time mean fields:*/
    348         for (int f=0;f<nfields;f++){
    349                 if (meanxtype[f]==1){ /*deal with scalars {{{*/
    350 
    351                         /*we are broadcasting doubles:*/
    352                         IssmDouble maxscalar=*maxmeans[f];
    353                         IssmDouble minscalar=*minmeans[f];
    354                         IssmDouble allmaxscalar;
    355                         IssmDouble allminscalar;
    356 
    357                         ISSM_MPI_Allreduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
    358                         ISSM_MPI_Allreduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
    359 
    360                         /*Store for later use in histogram computation:*/
    361                         *maxmeans[f]=allmaxscalar;
    362                         *minmeans[f]=allminscalar;
    363 
    364                 } /*}}}*/
    365                 else{ /*deal with arrays:{{{*/
    366 
    367                         int size=meanxsize[f];
    368 
    369                         /*we are broadcasting double arrays:*/
    370                         IssmDouble* maxx=maxmeans[f];
    371                         IssmDouble* minx=minmeans[f];
    372 
    373                         IssmDouble*  allmax=xNew<IssmDouble>(size);
    374                         IssmDouble*  allmin=xNew<IssmDouble>(size);
    375 
    376                         ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
    377                         ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
    378 
    379                         /*Store for later use in histogram computation:*/
    380                         maxmeans[f]=allmax;
    381                         minmeans[f]=allmin;
    382 
    383                 } /*}}}*/
    384         }
    385 
    386         /*Now that we have the min and max, we can start binning. First allocate
    387          * histograms, then start filling them:*/
    388         IssmDouble** histogram=xNew<IssmDouble*>(nfields*nsteps);
    389         IssmDouble** timehistogram=xNew<IssmDouble*>(nfields);
    390 
    391         _printf0_("Start reading files again, this time binning values in the histogram:\n");
    392         /*Start opening files:*/
    393         for (int i=(lower_row+1);i<=upper_row;i++){
    394                 _printf0_("reading file #: " << i << "\n");
    395                 char file[1000];
    396                 long int  length;
    397                 char* buffer=NULL;
    398 
    399                 /*string:*/
    400                 sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
    401 
    402                 /*open file: */
    403                 _printf0_("    opening file:\n");
    404                 FILE* fid=fopen(file,"rb");
    405                 if(fid==NULL)_error_("cound not open file: " << file << "\n");
    406 
    407                 /*figure out size of file, and read the whole thing:*/
    408                 _printf0_("    reading file:\n");
    409                 fseek (fid, 0, SEEK_END);
    410                 length = ftell (fid);
    411                 fseek (fid, 0, SEEK_SET);
    412                 buffer = xNew<char>(length);
    413                 fread (buffer, sizeof(char), length, fid);
    414 
    415                 /*close file:*/
    416                 fclose (fid);
    417 
    418                 /*create a memory stream with this buffer:*/
    419                 _printf0_("    processing file:\n");
    420                 fid=fmemopen(buffer, length, "rb");
    421 
    422                 /*start reading data from the buffer directly:*/
    423                 for (int f=0;f<nfields;f++){
    424                         char* field=fields[f];
    425                         fseek(fid,0,SEEK_SET);
    426                         for (int j=0;j<nsteps;j++){
    427                                 int counter=f*nsteps+j;
    428                                 xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    429                                 if(i==(lower_row+1)){
    430                                         if(xtype[counter]==1){
    431                                                 IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(nbins);
    432                                                 IssmDouble ma=*maxxs[counter];
    433                                                 IssmDouble mi=*minxs[counter];
    434                                                 int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index--;
    435                                                 if(ma==mi)index=0;
    436                                                 //_printf_( index << "|" << scalar << "|" << mi << "|" << ma << "|" << nbins << "\n");
    437                                                 localhistogram[index]++;
    438                                                 histogram[counter]=localhistogram;
    439                                         }
    440                                         else if (xtype[counter]==3){
    441                                                 IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(doublematsize*nbins);
    442                                                 IssmDouble* ma=maxxs[counter];
    443                                                 IssmDouble* mi=minxs[counter];
    444                                                 for (int k=0;k<doublematsize;k++){
    445                                                         IssmDouble scalar=doublemat[k];
    446                                                         int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index--;
    447                                                         if (mi[k]==ma[k])index=0;
    448                                                         _assert_(scalar<=ma[k]); _assert_(scalar>=mi[k]); _assert_(index<nbins);
    449                                                         localhistogram[k*nbins+index]++;
    450                                                 }
    451                                                 histogram[counter]=localhistogram;
    452                                                 xDelete<IssmDouble>(doublemat);
    453                                         }
    454                                         else _error_("cannot carry out statistics on type " << xtype[counter]);
    455                                 }
    456                                 else{
    457                                         if(xtype[counter]==1){
    458                                                 IssmDouble* localhistogram=histogram[counter];
    459                                                 IssmDouble ma=*maxxs[counter];
    460                                                 IssmDouble mi=*minxs[counter];
    461                                                 int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
    462                                                 if(ma==mi)index=0;
    463                                                 localhistogram[index]++;
    464                                         }
    465                                         else if (xtype[counter]==3){
    466                                                 IssmDouble* localhistogram=histogram[counter];
    467                                                 IssmDouble* ma=maxxs[counter];
    468                                                 IssmDouble* mi=minxs[counter];
    469                                                 for (int k=0;k<doublematsize;k++){
    470                                                         IssmDouble scalar=doublemat[k];
    471                                                         int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
    472                                                         if (mi[k]==ma[k])index=0;
    473 
    474                                                         localhistogram[k*nbins+index]++;
    475                                                 }
    476                                                 xDelete<IssmDouble>(doublemat);
    477                                         }
    478                                         else _error_("cannot carry out statistics on type " << xtype[counter]);
    479                                 }
    480                         }
    481                 }
    482                 _printf0_("    average in time:\n");
    483 
    484                 /*Deal with average in time: */
    485                 for (int f=0;f<nfields;f++){
    486                         fseek(fid,0,SEEK_SET);
    487                         char* field=fields[f];
    488                         meanxtype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
    489 
    490                         if(meanxtype[f]==1){
    491                                 IssmDouble timemean=0;
    492                                 fseek(fid,0,SEEK_SET);
    493                                 for (int j=0;j<nsteps;j++){
    494                                         readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    495                                         timemean+=scalar/nsteps;
    496                                 }
    497 
    498                                 /*Figure out max and min of time means: */
    499                                 if(i==(lower_row+1)){
    500                                         IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(nbins);
    501                                         IssmDouble ma=*maxmeans[f];
    502                                         IssmDouble mi=*minmeans[f];
    503                                         int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
    504                                         if(ma==mi)index=0;
    505                                         localhistogram[index]++;
    506                                         timehistogram[f]=localhistogram;
    507                                 }
    508                                 else{
    509                                         IssmDouble* localhistogram=timehistogram[f];
    510                                         IssmDouble ma=*maxmeans[f];
    511                                         IssmDouble mi=*minmeans[f];
    512                                         int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
    513                                         if(ma==mi)index=0;
    514                                         localhistogram[index]++;
    515                                 }
    516                         }
    517                         else{
    518                                 fseek(fid,0,SEEK_SET);
    519                                 IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
    520                                 for (int j=0;j<nsteps;j++){
    521                                         readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    522                                         for (int k=0;k<doublematsize;k++){
    523                                                 timemean[k]+=doublemat[k]/nsteps;
    524                                         }
    525                                         xDelete<IssmDouble>(doublemat);
    526                                 }
    527 
    528                                 if(i==(lower_row+1)){
    529                                         IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(doublematsize*nbins);
    530                                         IssmDouble* ma=maxmeans[f];
    531                                         IssmDouble* mi=minmeans[f];
    532 
    533                                         for (int k=0;k<doublematsize;k++){
    534                                                 IssmDouble scalar=timemean[k];
    535                                                 int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
    536                                                 if (mi[k]==ma[k])index=0;
    537                                                 localhistogram[k*nbins+index]++;
    538                                         }
    539                                         timehistogram[f]=localhistogram;
    540                                 }
    541                                 else{
    542 
    543                                         IssmDouble* localhistogram=timehistogram[f];
    544                                         IssmDouble* ma=maxmeans[f];
    545                                         IssmDouble* mi=minmeans[f];
    546 
    547                                         for (int k=0;k<doublematsize;k++){
    548                                                 IssmDouble scalar=timemean[k];
    549                                                 int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
    550                                                 if (mi[k]==ma[k])index=0;
    551 
    552                                                 localhistogram[k*nbins+index]++;
    553                                         }
    554                                 }
    555                         }
    556                 }
    557                 fclose(fid);
    558 
    559                 /*delete buffer:*/
    560                 xDelete<char>(buffer);
    561         }
    562         _printf0_("Start aggregating histogram:\n");
    563 
    564         /*We have agregated histograms across the cluster, now gather them across  the cluster onto
    565          *cpu0: */
    566         for (int f=0;f<nfields;f++){
    567                 int counter0=f*nsteps+0;
    568                 if (xtype[counter0]==1){ /*deal with scalars {{{*/
    569                         for (int j=0;j<nsteps;j++){
    570                                 int counter=f*nsteps+j;
    571 
    572                                 /*we are broadcasting doubles:*/
    573                                 IssmDouble* histo=histogram[counter]; //size nbins
    574                                 IssmDouble* allhisto=xNewZeroInit<IssmDouble>(nbins);
    575 
    576                                 ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
    577 
    578                                 /*add to results:*/
    579                                 if(my_rank==0){
    580                                         char fieldname[1000];
    581 
    582                                         sprintf(fieldname,"%s%s",fields[f],"Histogram");
    583                                         results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,1,nbins,j+1,0));
    584 
    585                                         sprintf(fieldname,"%s%s",fields[f],"Max");
    586                                         results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*maxxs[counter],j+1,0));
    587                                         sprintf(fieldname,"%s%s",fields[f],"Min");
    588                                         results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*minxs[counter],j+1,0));
    589                                 }
    590                         }
    591                 } /*}}}*/
    592                 else{ /*deal with arrays:{{{*/
    593 
    594                         int size=xsize[counter0];
    595                         for (int j=0;j<nsteps;j++){
    596                                 int counter=f*nsteps+j;
    597 
    598                                 /*we are broadcasting double arrays:*/
    599                                 IssmDouble* histo=histogram[counter];
    600                                 IssmDouble* allhisto=xNew<IssmDouble>(size*nbins);
    601 
    602                                 ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
    603                                 xDelete<IssmDouble>(histo);
    604 
    605                                 /*add to results:*/
    606                                 if(my_rank==0){
    607                                         char fieldname[1000];
    608 
    609                                         sprintf(fieldname,"%s%s",fields[f],"Histogram");
    610                                         results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,size,nbins,j+1,0));
    611 
    612                                         sprintf(fieldname,"%s%s",fields[f],"Max");
    613                                         results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,maxxs[counter],size,1,j+1,0));
    614                                         sprintf(fieldname,"%s%s",fields[f],"Min");
    615                                         results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,minxs[counter],size,1,j+1,0));
    616                                 }
    617                         }
    618                 } /*}}}*/
    619         }
    620         _printf0_("Start aggregating time mean histogram:\n");
    621 
    622         /*Now do the same for the time mean fields:*/
    623         for (int f=0;f<nfields;f++){
    624                 if (meanxtype[f]==1){ /*deal with scalars {{{*/
    625 
    626                         /*we are broadcasting doubles:*/
    627                         IssmDouble* histo=timehistogram[f];
    628                         IssmDouble* allhisto=xNewZeroInit<IssmDouble>(nbins);
    629 
    630                         ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
    631 
    632                         /*add to results at time step 1:*/
    633                         if(my_rank==0){
    634                                 char fieldname[1000];
    635 
    636                                 sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram");
    637                                 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,1,nbins,steps[0],0));
    638 
    639                                 sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax");
    640                                 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*maxmeans[f],steps[0],0));
    641                                 sprintf(fieldname,"%s%s",fields[f],"TimeMeaMin");
    642                                 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*minmeans[f],steps[0],0));
    643                         }
    644                 } /*}}}*/
    645                 else{ /*deal with arrays:{{{*/
    646 
    647                         int size=meanxsize[f];
    648 
    649                         /*we are broadcasting double arrays:*/
    650                         IssmDouble* histo=timehistogram[f];
    651                         IssmDouble* allhisto=xNewZeroInit<IssmDouble>(size*nbins);
    652 
    653                         ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
    654                         xDelete<IssmDouble>(histo);
    655                         /*add to results at step 1:*/
    656                         if(my_rank==0){
    657                                 char fieldname[1000];
    658 
    659                                 sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram");
    660                                 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,size,nbins,steps[0],0));
    661                                 sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax");
    662                                 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,maxmeans[f],size,1,steps[0],0));
    663                                 sprintf(fieldname,"%s%s",fields[f],"TimeMeanMin");
    664                                 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,minmeans[f],size,1,steps[0],0));
    665                         }
    666                 } /*}}}*/
    667         }
    668         _printf0_("Done aggregating time mean histogram:\n");
    669         IssmComm::SetComm(ISSM_MPI_COMM_WORLD);
    670 
    671         return 1;
    672 }
    673 /*}}}*/
    674 int ComputeMeanVariance(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){  /*{{{*/
    675 
    676         int nsamples;
    677         char* directory=NULL;
    678         char* model=NULL;
    679         char** fields=NULL;
    680         int* steps=NULL;
    681         int nsteps;
    682         int nfields;
    683         int range,lower_row,upper_row;
    684         int nfilesperdirectory;
    685 
    686         /*intermediary:*/
    687         IssmDouble* doublemat=NULL;
    688         int         doublematsize;
    689         IssmDouble scalar;
    690 
    691         /*computation of average and variance itself:*/
    692         IssmDouble*  x = NULL;
    693         IssmDouble*  x2 = NULL;
    694         IssmDouble** xs = NULL;
    695         IssmDouble** xs2 = NULL;
    696         int*         xtype=NULL;
    697         int*         xsize=NULL;
    698 
    699         IssmDouble** meanx=NULL;
    700         IssmDouble** meanx2=NULL;
    701         int*         meantype=NULL;
    702         int*         meansize=NULL;
    703 
    704         /*only work on the statistical communicator: */
    705         if (color==MPI_UNDEFINED)return 0;
    706 
    707         /*Retrieve parameters:*/
    708         parameters->FindParam(&nfilesperdirectory,QmuNfilesPerDirectoryEnum);
    709         parameters->FindParam(&nsamples,QmuNsampleEnum);
    710         parameters->FindParam(&directory,DirectoryNameEnum);
    711         parameters->FindParam(&model,InputFileNameEnum);
    712         parameters->FindParam(&fields,&nfields,FieldsEnum);
    713         parameters->FindParam(&steps,&nsteps,StepsEnum);
    714 
    715         /*Get rank from the stat comm communicator:*/
    716         IssmComm::SetComm(statcomm);
    717         int my_rank=IssmComm::GetRank();
    718 
    719         /*Open files and read them complelety, in a distributed way:*/
    720         range=DetermineLocalSize(nsamples,IssmComm::GetComm());
    721         GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
    722 
    723         /*Initialize arrays:*/
    724         xs=xNew<IssmDouble*>(nfields*nsteps);
    725         xs2=xNew<IssmDouble*>(nfields*nsteps);
    726         xtype=xNew<int>(nfields*nsteps);
    727         xsize=xNew<int>(nfields*nsteps);
    728 
    729         meantype=xNew<int>(nfields);
    730         meansize=xNew<int>(nfields);
    731         meanx=xNew<IssmDouble*>(nfields);
    732         meanx2=xNew<IssmDouble*>(nfields);
    733 
    734         /*Start opening files:*/
    735         for (int i=(lower_row+1);i<=upper_row;i++){
    736                 _printf0_("reading file #: " << i << "\n");
    737                 char file[1000];
    738                 long int  length;
    739                 char* buffer=NULL;
    740 
    741                 /*string:*/
    742                 sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
    743 
    744                 /*open file: */
    745                 _printf0_("    opening file: " << file << "\n");
    746                 FILE* fid=fopen(file,"rb");
    747                 if(fid==NULL) _error_("    could not open file: " << file << "\n");
    748 
    749                 /*figure out size of file, and read the whole thing:*/
    750                 _printf0_("    reading file:\n");
    751                 fseek (fid, 0, SEEK_END);
    752                 length = ftell (fid);
    753                 fseek (fid, 0, SEEK_SET);
    754                 buffer = xNew<char>(length);
    755                 fread (buffer, sizeof(char), length, fid);
    756 
    757                 /*close file:*/
    758                 fclose (fid);
    759 
    760                 /*create a memory stream with this buffer:*/
    761                 _printf0_("    processing file:\n");
    762                 fid=fmemopen(buffer, length, "rb");
    763 
    764                 /*start reading data from the buffer directly:*/
    765                 for (int f=0;f<nfields;f++){
    766                         char* field=fields[f];
    767                         fseek(fid,0,SEEK_SET);
    768                         for (int j=0;j<nsteps;j++){
    769                                 int counter=f*nsteps+j;
    770                                 xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    771                                 if(i==(lower_row+1)){
    772                                         if(xtype[counter]==1){
    773                                                 xs[counter]=xNew<IssmDouble>(1);
    774                                                 xs2[counter]=xNew<IssmDouble>(1);
    775                                                 *xs[counter]=scalar;
    776                                                 *xs2[counter]=pow(scalar,2.0);
    777                                                 xsize[counter]=1;
    778                                         }
    779                                         else if (xtype[counter]==3){
    780                                                 IssmDouble* doublemat2=xNew<IssmDouble>(doublematsize);
    781                                                 for(int k=0;k<doublematsize;k++)doublemat2[k]=pow(doublemat[k],2.0);
    782                                                 xs[counter]=doublemat;
    783                                                 xs2[counter]=doublemat2;
    784                                                 xsize[counter]=doublematsize;
    785                                         }
    786                                         else _error_("cannot carry out statistics on type " << xtype[counter]);
    787                                 }
    788                                 else{
    789                                         if(xtype[counter]==1){
    790                                                 *xs[counter]+=scalar;
    791                                                 *xs2[counter]+=pow(scalar,2.0);
    792                                         }
    793                                         else if (xtype[counter]==3){
    794                                                 IssmDouble* newdoublemat=xs[counter];
    795                                                 IssmDouble* newdoublemat2=xs2[counter];
    796                                                 for(int k=0;k<doublematsize;k++){
    797                                                         newdoublemat[k]+=doublemat[k];
    798                                                         newdoublemat2[k]+=pow(doublemat[k],2.0);
    799                                                 }
    800                                                 xs[counter]=newdoublemat;
    801                                                 xs2[counter]=newdoublemat2;
    802                                         }
    803                                         else _error_("cannot carry out statistics on type " << xtype[counter]);
    804                                 }
    805                         }
    806                 }
    807 
    808                 /*Deal with time mean: */
    809                 for (int f=0;f<nfields;f++){
    810                         char* field=fields[f];
    811                         fseek(fid,0,SEEK_SET);
    812                         meantype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
    813                         if(i==(lower_row+1)){
    814                                 if(meantype[f]==1){
    815                                         meanx[f]=xNewZeroInit<IssmDouble>(1);
    816                                         meanx2[f]=xNewZeroInit<IssmDouble>(1);
    817                                         meansize[f]=1;
    818                                 }
    819                                 else{
    820                                         meanx[f]=xNewZeroInit<IssmDouble>(doublematsize);
    821                                         meanx2[f]=xNewZeroInit<IssmDouble>(doublematsize);
    822                                         meansize[f]=doublematsize;
    823                                 }
    824                         }
    825                         fseek(fid,0,SEEK_SET);
    826                         if(meantype[f]==1){
    827                                 IssmDouble sc=0;
    828                                 IssmDouble sc2=0;
    829                                 for(int j=0;j<nsteps;j++){
    830                                         readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    831                                         sc+=scalar/nsteps;
    832                                 }
    833                                 sc2+=pow(sc,2.0);
    834                                 *meanx[f]+=sc;
    835                                 *meanx2[f]+=sc2;
    836                         }
    837                         else{
    838                                 IssmDouble* sc=meanx[f];
    839                                 IssmDouble* sc2=meanx2[f];
    840                                 IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
    841                                 IssmDouble* timemean2=xNewZeroInit<IssmDouble>(doublematsize);
    842 
    843                                 for(int j=0;j<nsteps;j++){
    844                                         readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    845                                         for (int k=0;k<doublematsize;k++){
    846                                                 timemean[k]+=doublemat[k]/nsteps;
    847                                         }
    848                                 }
    849                                 for (int k=0;k<doublematsize;k++){
    850                                         timemean2[k]=pow(timemean[k],2.0);
    851                                 }
    852                                 for (int k=0;k<doublematsize;k++){
    853                                         sc[k]+=timemean[k];
    854                                         sc2[k]+=timemean2[k];
    855                                 }
    856 
    857                         }
    858 
    859                 }
    860                 fclose(fid);
    861 
    862                 /*delete buffer:*/
    863                 xDelete<char>(buffer);
    864         }
    865         ISSM_MPI_Barrier(IssmComm::GetComm());
    866         _printf0_("Done reading files, now computing mean and variance.\n");
    867 
    868         /*We have agregated x and x^2 across the cluster, now gather across the cluster onto
    869          *cpu0 and then compute statistics:*/
    870         for (int f=0;f<nfields;f++){
    871                 int counter0=f*nsteps+0;
    872                 if (xtype[counter0]==1){ /*deal with scalars {{{*/
    873                         IssmDouble mean,stddev;
    874                         for (int j=0;j<nsteps;j++){
    875                                 int counter=f*nsteps+j;
    876 
    877                                 /*we are broadcasting doubles:*/
    878                                 IssmDouble scalar=*xs[counter];
    879                                 IssmDouble scalar2=*xs2[counter];
    880                                 IssmDouble sumscalar;
    881                                 IssmDouble sumscalar2;
    882 
    883                                 ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    884                                 ISSM_MPI_Reduce(&scalar2,&sumscalar2,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    885                                 /*Build average and standard deviation. For standard deviation, use the
    886                                  *following formula: sigma^2=E(x^2)-mu^2:*/
    887                                 mean=sumscalar/nsamples;
    888                                 stddev=sqrt(sumscalar2/nsamples-pow(mean,2.0));
    889 
    890                                 /*add to results:*/
    891                                 if(my_rank==0){
    892                                         char fieldname[1000];
    893 
    894                                         sprintf(fieldname,"%s%s",fields[f],"Mean");
    895                                         results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,j+1,0));
    896                                         sprintf(fieldname,"%s%s",fields[f],"Stddev");
    897                                         results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,j+1,0));
    898                                 }
    899 
    900                         }
    901                 } /*}}}*/
    902                 else{ /*deal with arrays:{{{*/
    903 
    904                         int size=xsize[counter0];
    905 
    906                         IssmDouble*  mean=xNew<IssmDouble>(size);
    907                         IssmDouble*  stddev=xNew<IssmDouble>(size);
    908 
    909                         for (int j=0;j<nsteps;j++){
    910                                 int counter=f*nsteps+j;
    911 
    912                                 /*we are broadcasting double arrays:*/
    913                                 x=xs[counter];
    914                                 x2=xs2[counter];
    915 
    916                                 IssmDouble*  sumx=xNew<IssmDouble>(size);
    917                                 IssmDouble*  sumx2=xNew<IssmDouble>(size);
    918 
    919                                 ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    920                                 ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    921 
    922                                 /*Build average and standard deviation. For standard deviation, use the
    923                                  *following formula: sigma^2=E(x^2)-mu^2:*/
    924                                 for (int k=0;k<size;k++){
    925                                         mean[k]=sumx[k]/nsamples;
    926                                         stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
    927                                 }
    928 
    929                                 /*add to results:*/
    930                                 if(my_rank==0){
    931                                         char fieldname[1000];
    932 
    933                                         sprintf(fieldname,"%s%s",fields[f],"Mean");
    934                                         results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,j+1,0));
    935                                         sprintf(fieldname,"%s%s",fields[f],"Stddev");
    936                                         results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,j+1,0));
    937                                 }
    938                         }
    939                 } /*}}}*/
    940         }
    941         /*Do the same but for the time mean:*/
    942         for (int f=0;f<nfields;f++){
    943                 if (meantype[f]==1){ /*deal with scalars {{{*/
    944                         IssmDouble mean,stddev;
    945 
    946                         /*we are broadcasting doubles:*/
    947                         IssmDouble scalar=*meanx[f];
    948                         IssmDouble scalar2=*meanx2[f];
    949                         IssmDouble sumscalar;
    950                         IssmDouble sumscalar2;
    951 
    952                         ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    953                         ISSM_MPI_Reduce(&scalar2,&sumscalar2,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    954                         /*Build average and standard deviation. For standard deviation, use the
    955                          *following formula: sigma^2=E(x^2)-mu^2:*/
    956                         mean=sumscalar/nsamples;
    957                         stddev=sqrt(sumscalar2/nsamples-pow(mean,2.0));
    958 
    959                         /*add to results:*/
    960                         if(my_rank==0){
    961                                 char fieldname[1000];
    962 
    963                                 sprintf(fieldname,"%s%s",fields[f],"TimeMean");
    964                                 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,steps[0],0));
    965                                 sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
    966                                 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,steps[0],0));
    967                         }
    968                 } /*}}}*/
    969                 else{ /*deal with arrays:{{{*/
    970 
    971                         int size=meansize[f];
    972                         IssmDouble*  mean=xNew<IssmDouble>(size);
    973                         IssmDouble*  stddev=xNew<IssmDouble>(size);
    974 
    975                         /*we are broadcasting double arrays:*/
    976                         x=meanx[f];
    977                         x2=meanx2[f];
    978 
    979                         IssmDouble*  sumx=xNew<IssmDouble>(size);
    980                         IssmDouble*  sumx2=xNew<IssmDouble>(size);
    981 
    982                         ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    983                         ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
    984 
    985                         /*Build average and standard deviation. For standard deviation, use the
    986                          *following formula: sigma^2=E(x^2)-mu^2:*/
    987                         for (int k=0;k<size;k++){
    988                                 mean[k]=sumx[k]/nsamples;
    989                                 stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
    990                         }
    991 
    992                         /*add to results:*/
    993                         if(my_rank==0){
    994                                 char fieldname[1000];
    995 
    996                                 sprintf(fieldname,"%s%s",fields[f],"TimeMean");
    997                                 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,steps[0],0));
    998                                 sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
    999                                 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,steps[0],0));
    1000                         }
    1001                 } /*}}}*/
    1002         }
    1003 
    1004         _printf0_("Done with MeanVariance:\n");
    1005         IssmComm::SetComm(ISSM_MPI_COMM_WORLD);
    1006 
    1007         return 1;
    1008 } /*}}}*/
    1009 int ComputeSampleSeries(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){ /*{{{*/
    1010 
    1011         int nsamples;
    1012         char* directory=NULL;
    1013         char* model=NULL;
    1014         char** fields=NULL;
    1015         int* steps=NULL;
    1016         int nsteps;
    1017         int nfields;
    1018         int range,lower_row,upper_row;
    1019         int nfilesperdirectory;
    1020         int* indices=NULL;
    1021         int  nindices;
    1022 
    1023         /*intermediary:*/
    1024         IssmDouble* doublemat=NULL;
    1025         int         doublematsize;
    1026         IssmDouble scalar;
    1027 
    1028         /*computation of average and variance itself:*/
    1029         IssmDouble*  x = NULL;
    1030         IssmDouble*  allx=NULL;
    1031         IssmDouble** xs = NULL;
    1032         int*         xtype=NULL;
    1033         int*         xsize=NULL;
    1034 
    1035         /*only work on the statistical communicator: */
    1036         if (color==MPI_UNDEFINED)return 0;
    1037 
    1038         /*Retrieve parameters:*/
    1039         parameters->FindParam(&nsamples,QmuNsampleEnum);
    1040         parameters->FindParam(&directory,DirectoryNameEnum);
    1041         parameters->FindParam(&model,InputFileNameEnum);
    1042         parameters->FindParam(&fields,&nfields,FieldsEnum);
    1043         parameters->FindParam(&steps,&nsteps,StepsEnum);
    1044         parameters->FindParam(&indices,&nindices,IndicesEnum);
    1045 
    1046         /*Get rank from the stat comm communicator:*/
    1047         IssmComm::SetComm(statcomm);
    1048         int my_rank=IssmComm::GetRank();
    1049 
    1050         /*Open files and read them complelety, in a distributed way:*/
    1051         range=DetermineLocalSize(nsamples,IssmComm::GetComm());
    1052         GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
    1053 
    1054         /*Initialize arrays:*/
    1055         xs=xNew<IssmDouble*>(nfields*nsteps);
    1056         xtype=xNew<int>(nfields*nsteps);
    1057         xsize=xNew<int>(nfields*nsteps);
    1058 
    1059         /*Start opening files:*/
    1060         for (int i=(lower_row+1);i<=upper_row;i++){
    1061                 _printf0_("reading file #: " << i << "\n");
    1062                 char file[1000];
    1063                 long int  length;
    1064                 char* buffer=NULL;
    1065 
    1066                 /*string:*/
    1067                 sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
    1068 
    1069                 /*open file: */
    1070                 _printf0_("    opening file:\n");
    1071                 FILE* fid=fopen(file,"rb");
    1072 
    1073                 /*figure out size of file, and read the whole thing:*/
    1074                 _printf0_("    reading file:\n");
    1075                 fseek (fid, 0, SEEK_END);
    1076                 length = ftell (fid);
    1077                 fseek (fid, 0, SEEK_SET);
    1078                 buffer = xNew<char>(length);
    1079                 fread (buffer, sizeof(char), length, fid);
    1080 
    1081                 /*close file:*/
    1082                 fclose (fid);
    1083 
    1084                 /*create a memory stream with this buffer:*/
    1085                 _printf0_("    processing file:\n");
    1086                 fid=fmemopen(buffer, length, "rb");
    1087 
    1088                 /*start reading data from the buffer directly:*/
    1089                 for (int f=0;f<nfields;f++){
    1090                         fseek(fid,0,SEEK_SET);
    1091                         char* field=fields[f];
    1092                         for (int j=0;j<nsteps;j++){
    1093                                 int counter=f*nsteps+j;
    1094                                 xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
    1095                                 if(i==(lower_row+1)){
    1096                                         if(xtype[counter]==1){
    1097                                                 x=xNew<IssmDouble>(range);
    1098                                                 x[0]=scalar;
    1099                                                 xs[counter]=x;
    1100                                                 xsize[counter]=range;
    1101                                         }
    1102                                         else if (xtype[counter]==3){
    1103                                                 x=xNew<IssmDouble>(nindices*range);
    1104                                                 for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
    1105                                                 xs[counter]=x;
    1106                                                 xsize[counter]=range*nindices;
    1107                                         }
    1108                                         else _error_("cannot carry out statistics on type " << xtype[counter]);
    1109                                 }
    1110                                 else{
    1111                                         if(xtype[counter]==1){
    1112                                                 x=xs[counter];
    1113                                                 x[i-(lower_row+1)]=scalar;
    1114                                                 xs[counter]=x;
    1115                                         }
    1116                                         else if (xtype[counter]==3){
    1117                                                 x=xs[counter];
    1118                                                 for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
    1119                                                 xs[counter]=x;
    1120                                         }
    1121                                         else _error_("cannot carry out statistics on type " << xtype[counter]);
    1122                                 }
    1123                         }
    1124                 }
    1125                 fclose(fid);
    1126 
    1127                 /*delete buffer:*/
    1128                 xDelete<char>(buffer);
    1129         }
    1130         ISSM_MPI_Barrier(IssmComm::GetComm());
    1131         _printf0_("Done reading files, now assembling time series.\n");
    1132 
    1133         for (int f=0;f<nfields;f++){
    1134                 for (int j=0;j<nsteps;j++){
    1135                         int counter=f*nsteps+j;
    1136                         if (xtype[counter]==1){
    1137                                 /*we are broadcasting range times doubles:*/
    1138                                 x=xs[counter];
    1139                                 allx=xNew<IssmDouble>(nsamples);
    1140                                 MPI_Gather(x, range, ISSM_MPI_PDOUBLE,allx, range, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
    1141                                 /*add to results:*/
    1142                                 if(my_rank==0){
    1143                                         char fieldname[1000];
    1144 
    1145                                         sprintf(fieldname,"%s%s",fields[f],"Samples");
    1146                                         results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,1,j+1,0));
    1147                                 }
    1148                         }
    1149                         else{
    1150                                 /*we are broadcasting double arrays:*/
    1151                                 x=xs[counter];
    1152                                 allx=xNew<IssmDouble>(nsamples*nindices);
    1153 
    1154                                 MPI_Gather(x, range*nindices, ISSM_MPI_PDOUBLE,allx, range*nindices, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
    1155 
    1156                                 /*add to results:*/
    1157                                 if(my_rank==0){
    1158                                         char fieldname[1000];
    1159                                         sprintf(fieldname,"%s%s",fields[f],"Samples");
    1160                                         results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,nindices,j+1,0));
    1161                                 }
    1162                         }
    1163                 }
    1164         }
    1165         _printf0_("Done with SampleSeries:\n");
    1166         IssmComm::SetComm(ISSM_MPI_COMM_WORLD);
    1167 
    1168         return 1;
    1169 } /*}}}*/
    1170 int OutputStatistics(Parameters* parameters,Results* results,int color,ISSM_MPI_Comm statcomm){ /*{{{*/
    1171 
    1172         char   outputfilename[1000];
    1173         char* directory=NULL;
    1174         char* model=NULL;
    1175         char* method=NULL;
    1176         int   nsamples;
    1177         int* steps=NULL;
    1178         int nsteps;
    1179 
    1180         /*only work on the statistical communicator: */
    1181         if (color==MPI_UNDEFINED)return 0;
    1182 
    1183         FemModel* femmodel=new FemModel();
    1184 
    1185         /*Some parameters that will allow us to use the OutputResultsx module:*/
    1186         parameters->AddObject(new BoolParam(QmuIsdakotaEnum,false));
    1187         parameters->AddObject(new BoolParam(SettingsIoGatherEnum,true));
    1188 
    1189         parameters->FindParam(&directory,DirectoryNameEnum);
    1190         parameters->FindParam(&model,InputFileNameEnum);
    1191         parameters->FindParam(&nsamples,QmuNsampleEnum);
    1192         parameters->FindParam(&steps,&nsteps,StepsEnum);
    1193 
    1194         sprintf(outputfilename,"%s/%s.stats",directory,model);
    1195         parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
    1196 
    1197         /*Call OutputResults module:*/
    1198         femmodel->parameters=parameters;
    1199         femmodel->results=results;
    1200 
    1201         OutputResultsx(femmodel);
    1202 
    1203         return 1;
    1204 } /*}}}*/
    1205 bool DakotaDirStructure(int argc,char** argv){ /*{{{*/
    1206 
    1207         char* input_file;
    1208         FILE* fid;
    1209         IoModel* iomodel=NULL;
    1210         int check;
    1211 
    1212         //qmu statistics
    1213         bool statistics    = false;
    1214         int  numdirectories = 0;
    1215 
    1216         /*First things first, set the communicator as a global variable: */
    1217         IssmComm::SetComm(MPI_COMM_WORLD);
    1218 
    1219         /*Barrier:*/
    1220         ISSM_MPI_Barrier(IssmComm::GetComm());
    1221         _printf0_("Preparing directory structure for model outputs:" << "\n");
    1222 
    1223         //open model input file for reading
    1224         input_file=xNew<char>((strlen(argv[2])+strlen(argv[3])+strlen(".bin")+2));
    1225         sprintf(input_file,"%s/%s%s",argv[2],argv[3],".bin");
    1226         fid=fopen(input_file,"rb");
    1227         if (fid==NULL) Cerr << "dirstructure error message: could not open model " << input_file << " to retrieve qmu statistics parameters" << std::endl;
    1228 
    1229         //initialize IoModel, but light version, we just need it to fetch one constant:
    1230         iomodel=new IoModel();
    1231         iomodel->fid=fid;
    1232         iomodel->FetchConstants();
    1233 
    1234         //early return if statistics not requested:
    1235         iomodel->FindConstant(&statistics,"md.qmu.statistics");
    1236         if(!statistics){
    1237                 delete iomodel;
    1238                 xDelete<char>(input_file);
    1239                 fclose(fid);
    1240                 return false; //important return value!
    1241         }
    1242 
    1243         iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories");
    1244 
    1245         /*Ok, we have everything we need to create the directory structure:*/
    1246         if(IssmComm::GetRank()==0){
    1247                 for (int i=0;i<numdirectories;i++){
    1248                         char directory[1000];
    1249                         sprintf(directory,"./%i",i+1);
    1250 
    1251                         check = mkdir(directory,ACCESSPERMS);
    1252                         if (check) _error_("dirstructure error message: could not create directory " << directory << "\n");
    1253                 }
    1254         }
    1255 
    1256         /*Delete resources:*/
    1257         delete iomodel;
    1258         xDelete<char>(input_file);
    1259 
    1260         //close model file:
    1261         fclose(fid);
    1262 
    1263         //return value:
    1264         return true; //statistics computation on!
    1265 } /*}}}*/
    12668int DakotaStatistics(int argc,char** argv){ /*{{{*/
    12679
     10        /*Variables:{{{*/
    126811        char* input_file;
    126912        FILE* fid;
    … …  
    129336        Parameters* parameters=NULL;
    129437        int color;
    1295 
    1296         /*First things first, set the communicator as a global variable: */
     38        /*}}}*/
     39        //First things first, set the communicator as a global variable and be sure we are all here: {{{
    129740        IssmComm::SetComm(MPI_COMM_WORLD);
    129841        my_rank=IssmComm::GetRank();
    … …  
    130144        ISSM_MPI_Barrier(IssmComm::GetComm());
    130245        _printf0_("Dakota Statistic Computation" << "\n");
    1303 
    1304         //open model input file for reading
     46        /*}}}*/
     47        //Open model input file for reading {{{
    130548        input_file=xNew<char>((strlen(argv[2])+strlen(argv[3])+strlen(".bin")+2));
    130649        sprintf(input_file,"%s/%s%s",argv[2],argv[3],".bin");
    130750        fid=fopen(input_file,"rb");
    130851        if (fid==NULL) Cerr << "issm_dakota_statistics error message: could not open model " << input_file << " to retrieve qmu statistics parameters" << std::endl;
    1309 
    1310         //initialize IoModel, but light version, we'll need it to fetch constants:
     52        //}}}
     53        //Initialize IoModel, but light version, we'll need it to fetch constants: {{{
    131154        iomodel=new IoModel();
    131255        iomodel->fid=fid;
    131356        iomodel->FetchConstants();
    1314 
    1315         //early return if statistics not requested:
     57        /*}}}*/
     58        //Early return if statistics not requested:  {{{
    131659        iomodel->FindConstant(&statistics,"md.qmu.statistics");
    131760        if(!statistics){
    … …  
    132063                fclose(fid);
    132164                return 0;
    1322         }else{
    1323                 //create parameters datasets with al the qmu statistics settings we need:
    1324 
    1325                 /*Initialize parameters and results:*/
    1326                 results   = new Results();
    1327                 parameters=new Parameters();
    1328 
    1329                 //solution type:
    1330                 parameters->AddObject(new IntParam(SolutionTypeEnum,StatisticsSolutionEnum));
    1331 
    1332                 //root  directory
    1333                 directory=xNew<char>(strlen(argv[2])+1);
    1334                 xMemCpy<char>(directory,argv[2],strlen(argv[2])+1);
    1335                 parameters->AddObject(new StringParam(DirectoryNameEnum,directory));
    1336 
    1337                 //model  name
    1338                 model=xNew<char>(strlen(argv[3])+1);
    1339                 xMemCpy<char>(model,argv[3],strlen(argv[3])+1);
    1340                 parameters->AddObject(new StringParam(InputFileNameEnum,model));
    1341 
    1342                 //nsamples
    1343                 iomodel->FindConstant(&nsamples,"md.qmu.method.params.samples");
    1344                 parameters->AddObject(new IntParam(QmuNsampleEnum,nsamples));
    1345 
    1346                 //ndirectories
    1347                 iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories");
    1348                 parameters->AddObject(new IntParam(QmuNdirectoriesEnum,numdirectories));
    1349 
    1350                 //nfiles per directory
    1351                 iomodel->FindConstant(&nfilesperdirectory,"md.qmu.statistics.nfiles_per_directory");
    1352                 parameters->AddObject(new IntParam(QmuNfilesPerDirectoryEnum,nfilesperdirectory));
    1353 
    1354                 //At this point, we don't want to go forward any longer, we want to create an MPI
    1355                 //communicator on which to carry out the computations:
    1356                 if ((my_rank+1)*nfilesperdirectory>nsamples)color=MPI_UNDEFINED;
    1357                 else color=0;
    1358                 ISSM_MPI_Comm_split(ISSM_MPI_COMM_WORLD,color, my_rank, &statcomm);
    1359 
    1360                 iomodel->FindConstant(&numstatistics,"md.qmu.statistics.numstatistics");
    1361                 for (int i=1;i<=numstatistics;i++){
    1362 
    1363                         char* directory=NULL;
    1364                         char* model=NULL;
    1365                         int   nsamples;
    1366                         _printf0_("Dealing with qmu statistical computation #" << i << "\n");
    1367 
    1368                         sprintf(string,"md.qmu.statistics.method(%i).name",i);
    1369                         iomodel->FindConstant(&name,string);
    1370 
    1371                         sprintf(string,"md.qmu.statistics.method(%i).fields",i);
    1372                         iomodel->FindConstant(&fields,&nfields,string);
    1373                         parameters->AddObject(new StringArrayParam(FieldsEnum,fields,nfields));
    1374 
    1375                         sprintf(string,"md.qmu.statistics.method(%i).steps",i);
    1376                         iomodel->FetchData(&steps,&dummy,&nsteps,string);
    1377                         parameters->AddObject(new IntVecParam(StepsEnum,steps,nsteps));
    1378 
    1379                         if (strcmp(name,"Histogram")==0){
    1380                                 /*fetch nbins: */
    1381                                 sprintf(string,"md.qmu.statistics.method(%i).nbins",i);
    1382                                 iomodel->FindConstant(&nbins,string);
    1383                                 parameters->AddObject(new IntParam(NbinsEnum,nbins));
    1384                                 ComputeHistogram(parameters,results,color,statcomm);
    1385                         }
    1386                         else if (strcmp(name,"SampleSeries")==0){
    1387                                 /*fetch indices: */
    1388                                 sprintf(string,"md.qmu.statistics.method(%i).indices",i);
    1389                                 iomodel->FetchData(&indices,&dummy,&nindices,string);
    1390                                 parameters->AddObject(new IntVecParam(IndicesEnum,indices,nindices));
    1391 
    1392                                 ComputeSampleSeries(parameters,results,color,statcomm);
    1393                         }
    1394                         else if (strcmp(name,"MeanVariance")==0){
    1395                                 ComputeMeanVariance(parameters,results,color,statcomm);
    1396                         }
    1397                         else _error_(" error creating qmu statistics methods parameters: unsupported method " << name);
     65        }
     66        /*}}}*/
     67        //Create parameters datasets with al the qmu statistics settings we need:  {{{
     68
     69        /*Initialize parameters and results:*/
     70        results   = new Results();
     71        parameters=new Parameters();
     72
     73        //solution type:
     74        parameters->AddObject(new IntParam(SolutionTypeEnum,StatisticsSolutionEnum));
     75
     76        //root  directory
     77        directory=xNew<char>(strlen(argv[2])+1);
     78        xMemCpy<char>(directory,argv[2],strlen(argv[2])+1);
     79        parameters->AddObject(new StringParam(DirectoryNameEnum,directory));
     80
     81        //model  name
     82        model=xNew<char>(strlen(argv[3])+1);
     83        xMemCpy<char>(model,argv[3],strlen(argv[3])+1);
     84        parameters->AddObject(new StringParam(InputFileNameEnum,model));
     85
     86        //nsamples
     87        iomodel->FindConstant(&nsamples,"md.qmu.method.params.samples");
     88        parameters->AddObject(new IntParam(QmuNsampleEnum,nsamples));
     89
     90        //ndirectories
     91        iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories");
     92        parameters->AddObject(new IntParam(QmuNdirectoriesEnum,numdirectories));
     93
     94        //nfiles per directory
     95        iomodel->FindConstant(&nfilesperdirectory,"md.qmu.statistics.nfiles_per_directory");
     96        parameters->AddObject(new IntParam(QmuNfilesPerDirectoryEnum,nfilesperdirectory));
     97        /*}}}*/
     98        /*Create MPI world: {{{*/
     99        //At this point, we don't want to go forward any longer, we want to create an MPI
     100        //communicator on which to carry out the computations:
     101        if ((my_rank+1)*nfilesperdirectory>nsamples)color=MPI_UNDEFINED;
     102        else color=0;
     103        ISSM_MPI_Comm_split(ISSM_MPI_COMM_WORLD,color, my_rank, &statcomm);
     104        /*}}}*/
     105
     106        iomodel->FindConstant(&numstatistics,"md.qmu.statistics.numstatistics");
     107        for (int i=1;i<=numstatistics;i++){
     108
     109                char* directory=NULL;
     110                char* model=NULL;
     111                int   nsamples;
     112                _printf0_("Dealing with qmu statistical computation #" << i << "\n");
     113
     114                sprintf(string,"md.qmu.statistics.method(%i).name",i);
     115                iomodel->FindConstant(&name,string);
     116
     117                sprintf(string,"md.qmu.statistics.method(%i).fields",i);
     118                iomodel->FindConstant(&fields,&nfields,string);
     119                parameters->AddObject(new StringArrayParam(FieldsEnum,fields,nfields));
     120
     121                sprintf(string,"md.qmu.statistics.method(%i).steps",i);
     122                iomodel->FetchData(&steps,&dummy,&nsteps,string);
     123                parameters->AddObject(new IntVecParam(StepsEnum,steps,nsteps));
     124
     125                if (strcmp(name,"Histogram")==0){
     126                        /*fetch nbins: */
     127                        sprintf(string,"md.qmu.statistics.method(%i).nbins",i);
     128                        iomodel->FindConstant(&nbins,string);
     129                        parameters->AddObject(new IntParam(NbinsEnum,nbins));
     130                        ComputeHistogram(parameters,results,color,statcomm);
    1398131                }
    1399 
    1400                 /*Delete resources:*/
    1401                 xDelete<char>(input_file);
    1402                 delete iomodel;
    1403         }
     132                else if (strcmp(name,"SampleSeries")==0){
     133                        /*fetch indices: */
     134                        sprintf(string,"md.qmu.statistics.method(%i).indices",i);
     135                        iomodel->FetchData(&indices,&dummy,&nindices,string);
     136                        parameters->AddObject(new IntVecParam(IndicesEnum,indices,nindices));
     137
     138                        ComputeSampleSeries(parameters,results,color,statcomm);
     139                }
     140                else if (strcmp(name,"MeanVariance")==0){
     141                        ComputeMeanVariance(parameters,results,color,statcomm);
     142                }
     143                else _error_(" error creating qmu statistics methods parameters: unsupported method " << name);
     144        }
     145
     146        /*Delete resources:*/
     147        xDelete<char>(input_file);
     148        delete iomodel;
    1404149
    1405150        //close model file:
    … …  
    1418163        return 1;
    1419164} /*}}}*/
     165int ComputeHistogram(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){  /*{{{*/
     166
     167        int nsamples;
     168        char* directory=NULL;
     169        char* model=NULL;
     170        char** fields=NULL;
     171        int* steps=NULL;
     172        int nsteps;
     173        int nfields;
     174        int nbins;
     175        int range,lower_row,upper_row;
     176        int nfilesperdirectory;
     177
     178        /*intermediary:*/
     179        IssmDouble* doublemat=NULL;
     180        int         doublematsize;
     181        IssmDouble scalar;
     182
     183        /*computation of average and variance itself:*/
     184        IssmDouble** maxxs = NULL;
     185        IssmDouble** minxs = NULL;
     186        int*         xtype=NULL;
     187        int*         xsize=NULL;
     188
     189        IssmDouble** maxmeans=NULL;
     190        IssmDouble** minmeans=NULL;
     191        int*         meanxtype=NULL;
     192        int*         meanxsize=NULL;
     193
     194        /*only work on the statistical communicator: */
     195        if (color==MPI_UNDEFINED)return 0;
     196
     197        /*Retrieve parameters:*/
     198        parameters->FindParam(&nfilesperdirectory,QmuNfilesPerDirectoryEnum);
     199        parameters->FindParam(&nsamples,QmuNsampleEnum);
     200        parameters->FindParam(&directory,DirectoryNameEnum);
     201        parameters->FindParam(&model,InputFileNameEnum);
     202        parameters->FindParam(&fields,&nfields,FieldsEnum);
     203        parameters->FindParam(&steps,&nsteps,StepsEnum);
     204        parameters->FindParam(&nbins,NbinsEnum);
     205
     206        /*Get rank from the stat comm communicator:*/
     207        IssmComm::SetComm(statcomm);
     208        int my_rank=IssmComm::GetRank();
     209
     210        /*Open files and read them complelety, in a distributed way:*/
     211        range=DetermineLocalSize(nsamples,IssmComm::GetComm());
     212        GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
     213
     214        /*Initialize arrays:*/
     215        maxmeans=xNew<IssmDouble*>(nfields);
     216        minmeans=xNew<IssmDouble*>(nfields);
     217        meanxtype=xNew<int>(nfields);
     218        meanxsize=xNew<int>(nfields);
     219
     220        maxxs=xNew<IssmDouble*>(nfields*nsteps);
     221        minxs=xNew<IssmDouble*>(nfields*nsteps);
     222        xtype=xNew<int>(nfields*nsteps);
     223        xsize=xNew<int>(nfields*nsteps);
     224
     225        /*Start opening files:*/
     226        for(int i=(lower_row+1);i<=upper_row;i++){
     227                _printf0_("reading file #: " << i << "\n");
     228                /*First read file to figure out size of it in order to create memory buffer mapping into the file.  {{{
     229                 *This makes it much more efficient to read files without lag.:*/
     230                char file[1000];
     231                long int  length;
     232                char* buffer=NULL;
     233
     234                /*string:*/
     235                sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
     236
     237                /*open file: */
     238                _printf0_("    opening file: " << file << "\n");
     239                FILE* fid=fopen(file,"rb");
     240                if(fid==NULL)_error_("cound not open file: " << file << "\n");
     241
     242                /*figure out size of file, and read the whole thing:*/
     243                _printf0_("    reading file:\n");
     244                fseek(fid, 0, SEEK_END);
     245                length = ftell (fid);
     246                fseek(fid, 0, SEEK_SET);
     247                buffer = xNew<char>(length);
     248                fread(buffer, sizeof(char), length, fid);
     249
     250                /*close file:*/
     251                fclose(fid);
     252
     253                /*create a memory stream with this buffer which will be use to read the files:*/
     254                _printf0_("    processing file:\n");
     255                fid=fmemopen(buffer, length, "rb");
     256                /*}}}*/
     257                /*Figure out for each field, each time step, arrays on each cpu holwing min anx max values:{{{*/
     258                for (int f=0;f<nfields;f++){
     259                        char* field=fields[f];
     260                        fseek(fid,0,SEEK_SET);
     261                        for (int j=0;j<nsteps;j++){
     262                                int counter=f*nsteps+j;
     263                                xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
     264                                if(i==(lower_row+1)){
     265                                        if(xtype[counter]==1){
     266                                                maxxs[counter]=xNew<IssmDouble>(1);
     267                                                minxs[counter]=xNew<IssmDouble>(1);
     268                                                *maxxs[counter]=scalar;
     269                                                *minxs[counter]=scalar;
     270                                                xsize[counter]=1;
     271                                        }
     272                                        else if (xtype[counter]==3){
     273                                                maxxs[counter]=xNew<IssmDouble>(doublematsize);
     274                                                xMemCpy<IssmDouble>(maxxs[counter],doublemat,doublematsize);
     275                                                minxs[counter]=xNew<IssmDouble>(doublematsize);
     276                                                xMemCpy<IssmDouble>(minxs[counter],doublemat,doublematsize);
     277                                                xsize[counter]=doublematsize;
     278                                                xDelete<IssmDouble>(doublemat);
     279                                        }
     280                                        else _error_("cannot carry out statistics on type " << xtype[counter]);
     281                                }
     282                                else{
     283                                        if(xtype[counter]==1){
     284                                                *maxxs[counter]=max(*maxxs[counter],scalar);
     285                                                *minxs[counter]=min(*minxs[counter],scalar);
     286                                        }
     287                                        else if (xtype[counter]==3){
     288                                                IssmDouble* newmax=maxxs[counter];
     289                                                IssmDouble* newmin=minxs[counter];
     290                                                for(int k=0;k<doublematsize;k++){
     291                                                        if(doublemat[k]>newmax[k])newmax[k]=doublemat[k];
     292                                                        if(doublemat[k]<newmin[k])newmin[k]=doublemat[k];
     293                                                }
     294                                                xDelete<IssmDouble>(doublemat);
     295                                        }
     296                                        else _error_("cannot carry out statistics on type " << xtype[counter]);
     297                                }
     298                        }
     299                }
     300                /*}}}*/
     301                /*Same thing for average in time:{{{*/
     302                _printf0_("    average in time:\n");
     303
     304                /*Deal with average in time: */
     305                for (int f=0;f<nfields;f++){
     306                        fseek(fid,0,SEEK_SET);
     307                        char* field=fields[f];
     308                        meanxtype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
     309
     310                        if(meanxtype[f]==1){
     311                                meanxsize[f]=1;
     312                                IssmDouble timemean=0;
     313                                fseek(fid,0,SEEK_SET);
     314                                for (int j=0;j<nsteps;j++){
     315                                        readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
     316                                        timemean+=scalar/nsteps;
     317                                }
     318
     319                                /*Figure out max and min of time means: */
     320                                if(i==(lower_row+1)){
     321                                        maxmeans[f]=xNewZeroInit<IssmDouble>(1);
     322                                        minmeans[f]=xNewZeroInit<IssmDouble>(1);
     323                                        *maxmeans[f]=timemean;
     324                                        *minmeans[f]=timemean;
     325                                }
     326                                else{
     327                                        *maxmeans[f]=max(*maxmeans[f],timemean);
     328                                        *minmeans[f]=min(*minmeans[f],timemean);
     329                                }
     330                        }
     331                        else{
     332                                meanxsize[f]=doublematsize;
     333                                fseek(fid,0,SEEK_SET);
     334                                IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
     335                                for (int j=0;j<nsteps;j++){
     336                                        readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
     337                                        for (int k=0;k<doublematsize;k++){
     338                                                timemean[k]+=doublemat[k]/nsteps;
     339                                        }
     340                                        xDelete<IssmDouble>(doublemat);
     341                                }
     342
     343                                if(i==(lower_row+1)){
     344                                        maxmeans[f]=xNew<IssmDouble>(doublematsize);
     345                                        xMemCpy<IssmDouble>(maxmeans[f],timemean,doublematsize);
     346                                        minmeans[f]=xNew<IssmDouble>(doublematsize);
     347                                        xMemCpy<IssmDouble>(minmeans[f],timemean,doublematsize);
     348                                }
     349                                else{
     350                                        IssmDouble* maxx=maxmeans[f];
     351                                        IssmDouble* minx=minmeans[f];
     352
     353                                        for(int k=0;k<doublematsize;k++){
     354                                                maxx[k]=max(maxx[k],timemean[k]);
     355                                                minx[k]=min(minx[k],timemean[k]);
     356                                        }
     357                                        maxmeans[f]=maxx;
     358                                        minmeans[f]=minx;
     359                                }
     360                        }
     361                }
     362                /*}}}*/
     363                /*Done reading files, close buffer and free memory:{{{*/
     364                fclose(fid);
     365                xDelete<char>(buffer);
     366                /*}}}*/
     367        }
     368        ISSM_MPI_Barrier(IssmComm::GetComm());
     369        _printf0_("Done reading files, now computing min and max.\n");
     370
     371        /*We have collected minx and max across the cluster, now gather across the cluster onto
     372         *cpu0 and then compute statistics:*/
     373        for (int f=0;f<nfields;f++){
     374                int counter0=f*nsteps+0;
     375                if (xtype[counter0]==1){ /*deal with scalars {{{*/
     376                        for (int j=0;j<nsteps;j++){
     377                                int counter=f*nsteps+j;
     378
     379                                /*we are broadcasting doubles:*/
     380                                IssmDouble maxscalar=*maxxs[counter];
     381                                IssmDouble minscalar=*minxs[counter];
     382                                IssmDouble allmaxscalar;
     383                                IssmDouble allminscalar;
     384                                IssmDouble sumscalar_alltimes=0;
     385
     386                                ISSM_MPI_Allreduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
     387                                ISSM_MPI_Allreduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
     388
     389                                /*Store broadcasted value for later computation of histograms:*/
     390                                *maxxs[counter]=allmaxscalar;
     391                                *minxs[counter]=allminscalar;
     392
     393                        }
     394                } /*}}}*/
     395                else{ /*deal with arrays:{{{*/
     396
     397                        int size=xsize[counter0];
     398                        for (int j=0;j<nsteps;j++){
     399                                int counter=f*nsteps+j;
     400
     401                                /*we are broadcasting double arrays:*/
     402                                IssmDouble* maxx=maxxs[counter];
     403                                IssmDouble* minx=minxs[counter];
     404
     405                                IssmDouble*  allmax=xNew<IssmDouble>(size);
     406                                IssmDouble*  allmin=xNew<IssmDouble>(size);
     407
     408                                ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
     409                                ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
     410
     411                                /*Store broadcasted value for later computation of histograms:*/
     412                                maxxs[counter]=allmax;
     413                                minxs[counter]=allmin;
     414                        }
     415                } /*}}}*/
     416        }
     417
     418        /*Now do the same for the time mean fields:*/
     419        for (int f=0;f<nfields;f++){
     420                if (meanxtype[f]==1){ /*deal with scalars {{{*/
     421
     422                        /*we are broadcasting doubles:*/
     423                        IssmDouble maxscalar=*maxmeans[f];
     424                        IssmDouble minscalar=*minmeans[f];
     425                        IssmDouble allmaxscalar;
     426                        IssmDouble allminscalar;
     427
     428                        ISSM_MPI_Allreduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
     429                        ISSM_MPI_Allreduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
     430
     431                        /*Store for later use in histogram computation:*/
     432                        *maxmeans[f]=allmaxscalar;
     433                        *minmeans[f]=allminscalar;
     434
     435                } /*}}}*/
     436                else{ /*deal with arrays:{{{*/
     437
     438                        int size=meanxsize[f];
     439
     440                        /*we are broadcasting double arrays:*/
     441                        IssmDouble* maxx=maxmeans[f];
     442                        IssmDouble* minx=minmeans[f];
     443
     444                        IssmDouble*  allmax=xNew<IssmDouble>(size);
     445                        IssmDouble*  allmin=xNew<IssmDouble>(size);
     446
     447                        ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
     448                        ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
     449
     450                        /*Store for later use in histogram computation:*/
     451                        maxmeans[f]=allmax;
     452                        minmeans[f]=allmin;
     453
     454                } /*}}}*/
     455        }
     456
     457        /*Now that we have the min and max, we can start binning. First allocate
     458         * histograms, then start filling them:*/
     459        IssmDouble** histogram=xNew<IssmDouble*>(nfields*nsteps);
     460        IssmDouble** timehistogram=xNew<IssmDouble*>(nfields);
     461        _printf0_("Start reading files again, this time binning values in the histogram:\n");
     462        /*Start opening files:*/
     463        for (int i=(lower_row+1);i<=upper_row;i++){
     464                _printf0_("reading file #: " << i << "\n");
     465                /*read file and make a buffer:{{{*/
     466                char file[1000];
     467                long int  length;
     468                char* buffer=NULL;
     469
     470                /*string:*/
     471                sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
     472
     473                /*open file: */
     474                _printf0_("    opening file:\n");
     475                FILE* fid=fopen(file,"rb");
     476                if(fid==NULL)_error_("cound not open file: " << file << "\n");
     477
     478                /*figure out size of file, and read the whole thing:*/
     479                _printf0_("    reading file:\n");
     480                fseek (fid, 0, SEEK_END);
     481                length = ftell (fid);
     482                fseek (fid, 0, SEEK_SET);
     483                buffer = xNew<char>(length);
     484                fread (buffer, sizeof(char), length, fid);
     485
     486                /*close file:*/
     487                fclose (fid);
     488
     489                /*create a memory stream with this buffer:*/
     490                _printf0_("    processing file:\n");
     491                fid=fmemopen(buffer, length, "rb");
     492                /*}}}*/
     493                /*read data and fill up the histogram using the min and max values from before:{{{*/
     494                for (int f=0;f<nfields;f++){
     495                        char* field=fields[f];
     496                        fseek(fid,0,SEEK_SET);
     497                        for (int j=0;j<nsteps;j++){
     498                                int counter=f*nsteps+j;
     499                                xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
     500                                if(i==(lower_row+1)){
     501                                        if(xtype[counter]==1){
     502                                                IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(nbins);
     503                                                IssmDouble ma=*maxxs[counter];
     504                                                IssmDouble mi=*minxs[counter];
     505                                                int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index--;
     506                                                if(ma==mi)index=0;
     507                                                //_printf_( index << "|" << scalar << "|" << mi << "|" << ma << "|" << nbins << "\n");
     508                                                localhistogram[index]++;
     509                                                histogram[counter]=localhistogram;
     510                                        }
     511                                        else if (xtype[counter]==3){
     512                                                IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(doublematsize*nbins);
     513                                                IssmDouble* ma=maxxs[counter];
     514                                                IssmDouble* mi=minxs[counter];
     515                                                for (int k=0;k<doublematsize;k++){
     516                                                        IssmDouble scalar=doublemat[k];
     517                                                        int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index--;
     518                                                        if (mi[k]==ma[k])index=0;
     519                                                        _assert_(scalar<=ma[k]); _assert_(scalar>=mi[k]); _assert_(index<nbins);
     520                                                        localhistogram[k*nbins+index]++;
     521                                                }
     522                                                histogram[counter]=localhistogram;
     523                                                xDelete<IssmDouble>(doublemat);
     524                                        }
     525                                        else _error_("cannot carry out statistics on type " << xtype[counter]);
     526                                }
     527                                else{
     528                                        if(xtype[counter]==1){
     529                                                IssmDouble* localhistogram=histogram[counter];
     530                                                IssmDouble ma=*maxxs[counter];
     531                                                IssmDouble mi=*minxs[counter];
     532                                                int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
     533                                                if(ma==mi)index=0;
     534                                                localhistogram[index]++;
     535                                        }
     536                                        else if (xtype[counter]==3){
     537                                                IssmDouble* localhistogram=histogram[counter];
     538                                                IssmDouble* ma=maxxs[counter];
     539                                                IssmDouble* mi=minxs[counter];
     540                                                for (int k=0;k<doublematsize;k++){
     541                                                        IssmDouble scalar=doublemat[k];
     542                                                        int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
     543                                                        if (mi[k]==ma[k])index=0;
     544
     545                                                        localhistogram[k*nbins+index]++;
     546                                                }
     547                                                xDelete<IssmDouble>(doublemat);
     548                                        }
     549                                        else _error_("cannot carry out statistics on type " << xtype[counter]);
     550                                }
     551                        }
     552                }
     553                /*}}}*/
     554                /*Deal with average in time: {{{*/
     555                _printf0_("    average in time:\n");
     556                for (int f=0;f<nfields;f++){
     557                        fseek(fid,0,SEEK_SET);
     558                        char* field=fields[f];
     559                        meanxtype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
     560
     561                        if(meanxtype[f]==1){
     562                                IssmDouble timemean=0;
     563                                fseek(fid,0,SEEK_SET);
     564                                for (int j=0;j<nsteps;j++){
     565                                        readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
     566                                        timemean+=scalar/nsteps;
     567                                }
     568
     569                                /*Figure out max and min of time means: */
     570                                if(i==(lower_row+1)){
     571                                        IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(nbins);
     572                                        IssmDouble ma=*maxmeans[f];
     573                                        IssmDouble mi=*minmeans[f];
     574                                        int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
     575                                        if(ma==mi)index=0;
     576                                        localhistogram[index]++;
     577                                        timehistogram[f]=localhistogram;
     578                                }
     579                                else{
     580                                        IssmDouble* localhistogram=timehistogram[f];
     581                                        IssmDouble ma=*maxmeans[f];
     582                                        IssmDouble mi=*minmeans[f];
     583                                        int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
     584                                        if(ma==mi)index=0;
     585                                        localhistogram[index]++;
     586                                }
     587                        }
     588                        else{
     589                                fseek(fid,0,SEEK_SET);
     590                                IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
     591                                for (int j=0;j<nsteps;j++){
     592                                        readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
     593                                        for (int k=0;k<doublematsize;k++){
     594                                                timemean[k]+=doublemat[k]/nsteps;
     595                                        }
     596                                        xDelete<IssmDouble>(doublemat);
     597                                }
     598
     599                                if(i==(lower_row+1)){
     600                                        IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(doublematsize*nbins);
     601                                        IssmDouble* ma=maxmeans[f];
     602                                        IssmDouble* mi=minmeans[f];
     603
     604                                        for (int k=0;k<doublematsize;k++){
     605                                                IssmDouble scalar=timemean[k];
     606                                                int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
     607                                                if (mi[k]==ma[k])index=0;
     608                                                localhistogram[k*nbins+index]++;
     609                                        }
     610                                        timehistogram[f]=localhistogram;
     611                                }
     612                                else{
     613
     614                                        IssmDouble* localhistogram=timehistogram[f];
     615                                        IssmDouble* ma=maxmeans[f];
     616                                        IssmDouble* mi=minmeans[f];
     617
     618                                        for (int k=0;k<doublematsize;k++){
     619                                                IssmDouble scalar=timemean[k];
     620                                                int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
     621                                                if (mi[k]==ma[k])index=0;
     622
     623                                                localhistogram[k*nbins+index]++;
     624                                        }
     625                                }
     626                        }
     627                }
     628                /*}}}*/
     629                /*close file and delete allocation:{{{*/
     630                fclose(fid);
     631                xDelete<char>(buffer);
     632                /*}}}*/
     633        }
     634
     635
     636        /*We have agregated histograms across the cluster, now gather them across  the cluster onto
     637         *cpu0: */
     638        _printf0_("Collect histograms on cpu 0 and save to results:\n");
     639        for (int f=0;f<nfields;f++){
     640                int counter0=f*nsteps+0;
     641                if (xtype[counter0]==1){ /*deal with scalars {{{*/
     642                        for (int j=0;j<nsteps;j++){
     643                                int counter=f*nsteps+j;
     644
     645                                /*we are broadcasting doubles:*/
     646                                IssmDouble* histo=histogram[counter]; //size nbins
     647                                IssmDouble* allhisto=xNewZeroInit<IssmDouble>(nbins);
     648
     649                                ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
     650                                xDelete<IssmDouble>(histo);
     651
     652                                /*add to results while deallocating as much as possible:*/
     653                                char fieldname[1000];
     654                                if(my_rank==0){
     655                                        sprintf(fieldname,"%s%s",fields[f],"Histogram"); results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,1,nbins,j+1,0));
     656                                }
     657                                xDelete<IssmDouble>(allhisto);
     658                                if(my_rank==0){
     659                                        sprintf(fieldname,"%s%s",fields[f],"Max"); results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*maxxs[counter],j+1,0));
     660                                        sprintf(fieldname,"%s%s",fields[f],"Min"); results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*minxs[counter],j+1,0));
     661                                }
     662                                xDelete<IssmDouble>(maxxs[counter]);
     663                                xDelete<IssmDouble>(minxs[counter]);
     664                        }
     665                } /*}}}*/
     666                else{ /*deal with arrays:{{{*/
     667
     668                        int size=xsize[counter0];
     669                        for (int j=0;j<nsteps;j++){
     670                                int counter=f*nsteps+j;
     671
     672                                /*we are broadcasting double arrays:*/
     673                                IssmDouble* histo=histogram[counter];
     674                                IssmDouble* allhisto=xNew<IssmDouble>(size*nbins);
     675
     676                                ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
     677                                xDelete<IssmDouble>(histo);
     678
     679                                /*add to results:*/
     680                                char fieldname[1000];
     681                                if(my_rank==0){
     682                                        sprintf(fieldname,"%s%s",fields[f],"Histogram"); results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,size,nbins,j+1,0));
     683                                }
     684                                xDelete<IssmDouble>(allhisto);
     685                                if(my_rank==0){
     686                                        sprintf(fieldname,"%s%s",fields[f],"Max"); results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,maxxs[counter],size,1,j+1,0));
     687                                        sprintf(fieldname,"%s%s",fields[f],"Min"); results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,minxs[counter],size,1,j+1,0));
     688                                }
     689                                xDelete<IssmDouble>(maxxs[counter]);
     690                                xDelete<IssmDouble>(minxs[counter]);
     691                        }
     692                } /*}}}*/
     693        }
     694        /*Now do the same for the time mean fields:*/
     695        _printf0_("Collect time mean histograms  on cpu 0 and save to results:\n");
     696        for (int f=0;f<nfields;f++){
     697                if (meanxtype[f]==1){ /*deal with scalars {{{*/
     698
     699                        /*we are broadcasting doubles:*/
     700                        IssmDouble* histo=timehistogram[f];
     701                        IssmDouble* allhisto=xNewZeroInit<IssmDouble>(nbins);
     702
     703                        ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
     704                        xDelete<IssmDouble>(histo);
     705
     706                        /*add to results at time step 1:*/
     707                        char fieldname[1000];
     708                        if(my_rank==0){
     709                                sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram"); results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,1,nbins,steps[0],0));
     710                        }
     711                        xDelete<IssmDouble>(allhisto);
     712                        if(my_rank==0){
     713                                sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax"); results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*maxmeans[f],steps[0],0));
     714                                sprintf(fieldname,"%s%s",fields[f],"TimeMeaMin"); results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*minmeans[f],steps[0],0));
     715                        }
     716                        xDelete<IssmDouble>(maxmeans[f]);
     717                        xDelete<IssmDouble>(minmeans[f]);
     718                } /*}}}*/
     719                else{ /*deal with arrays:{{{*/
     720
     721                        int size=meanxsize[f];
     722
     723                        /*we are broadcasting double arrays:*/
     724                        IssmDouble* histo=timehistogram[f];
     725                        IssmDouble* allhisto=xNewZeroInit<IssmDouble>(size*nbins);
     726
     727                        ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
     728                        xDelete<IssmDouble>(histo);
     729                       
     730                        /*add to results at step 1:*/
     731                        char fieldname[1000];
     732                        if(my_rank==0){
     733                                sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram"); results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,size,nbins,steps[0],0));
     734                        }
     735                        xDelete<IssmDouble>(allhisto);
     736                        if(my_rank==0){
     737                                sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax"); results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,maxmeans[f],size,1,steps[0],0));
     738                                sprintf(fieldname,"%s%s",fields[f],"TimeMeanMin"); results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,minmeans[f],size,1,steps[0],0));
     739                        }
     740                        xDelete<IssmDouble>(maxmeans[f]);
     741                        xDelete<IssmDouble>(minmeans[f]);
     742                } /*}}}*/
     743        }
     744        _printf0_("Done aggregating time mean histogram:\n");
     745        IssmComm::SetComm(ISSM_MPI_COMM_WORLD);
     746       
     747        /*Free allocations:*/
     748        xDelete<char>(directory);
     749        xDelete<char>(model);
     750        for (int i=0;i<nfields;i++)xDelete<char>(fields[i]);
     751        xDelete<char*>(fields);
     752        xDelete<int>(steps);
     753        xDelete<IssmDouble*>(maxxs);
     754        xDelete<IssmDouble*>(minxs);
     755        xDelete<IssmDouble*>(maxmeans);
     756        xDelete<IssmDouble*>(minmeans);
     757        xDelete<int>(xtype);
     758        xDelete<int>(xsize);
     759
     760        return 1;
     761}
     762/*}}}*/
     763int ComputeMeanVariance(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){  /*{{{*/
     764
     765        int nsamples;
     766        char* directory=NULL;
     767        char* model=NULL;
     768        char** fields=NULL;
     769        int* steps=NULL;
     770        int nsteps;
     771        int nfields;
     772        int range,lower_row,upper_row;
     773        int nfilesperdirectory;
     774
     775        /*intermediary:*/
     776        IssmDouble* doublemat=NULL;
     777        int         doublematsize;
     778        IssmDouble scalar;
     779
     780        /*computation of average and variance itself:*/
     781        IssmDouble*  x = NULL;
     782        IssmDouble*  x2 = NULL;
     783        IssmDouble** xs = NULL;
     784        IssmDouble** xs2 = NULL;
     785        int*         xtype=NULL;
     786        int*         xsize=NULL;
     787
     788        IssmDouble** meanx=NULL;
     789        IssmDouble** meanx2=NULL;
     790        int*         meantype=NULL;
     791        int*         meansize=NULL;
     792
     793        /*only work on the statistical communicator: */
     794        if (color==MPI_UNDEFINED)return 0;
     795
     796        /*Retrieve parameters:*/
     797        parameters->FindParam(&nfilesperdirectory,QmuNfilesPerDirectoryEnum);
     798        parameters->FindParam(&nsamples,QmuNsampleEnum);
     799        parameters->FindParam(&directory,DirectoryNameEnum);
     800        parameters->FindParam(&model,InputFileNameEnum);
     801        parameters->FindParam(&fields,&nfields,FieldsEnum);
     802        parameters->FindParam(&steps,&nsteps,StepsEnum);
     803
     804        /*Get rank from the stat comm communicator:*/
     805        IssmComm::SetComm(statcomm);
     806        int my_rank=IssmComm::GetRank();
     807
     808        /*Open files and read them complelety, in a distributed way:*/
     809        range=DetermineLocalSize(nsamples,IssmComm::GetComm());
     810        GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
     811
     812        /*Initialize arrays:*/
     813        xs=xNew<IssmDouble*>(nfields*nsteps);
     814        xs2=xNew<IssmDouble*>(nfields*nsteps);
     815        xtype=xNew<int>(nfields*nsteps);
     816        xsize=xNew<int>(nfields*nsteps);
     817
     818        meantype=xNew<int>(nfields);
     819        meansize=xNew<int>(nfields);
     820        meanx=xNew<IssmDouble*>(nfields);
     821        meanx2=xNew<IssmDouble*>(nfields);
     822
     823        /*Start opening files:*/
     824        for (int i=(lower_row+1);i<=upper_row;i++){
     825                _printf0_("reading file #: " << i << "\n");
     826                char file[1000];
     827                long int  length;
     828                char* buffer=NULL;
     829
     830                /*string:*/
     831                sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
     832
     833                /*open file: */
     834                _printf0_("    opening file: " << file << "\n");
     835                FILE* fid=fopen(file,"rb");
     836                if(fid==NULL) _error_("    could not open file: " << file << "\n");
     837
     838                /*figure out size of file, and read the whole thing:*/
     839                _printf0_("    reading file:\n");
     840                fseek (fid, 0, SEEK_END);
     841                length = ftell (fid);
     842                fseek (fid, 0, SEEK_SET);
     843                buffer = xNew<char>(length);
     844                fread (buffer, sizeof(char), length, fid);
     845
     846                /*close file:*/
     847                fclose (fid);
     848
     849                /*create a memory stream with this buffer:*/
     850                _printf0_("    processing file:\n");
     851                fid=fmemopen(buffer, length, "rb");
     852
     853                /*start reading data from the buffer directly:*/
     854                for (int f=0;f<nfields;f++){
     855                        char* field=fields[f];
     856                        fseek(fid,0,SEEK_SET);
     857                        for (int j=0;j<nsteps;j++){
     858                                int counter=f*nsteps+j;
     859                                xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
     860                                if(i==(lower_row+1)){
     861                                        if(xtype[counter]==1){
     862                                                xs[counter]=xNew<IssmDouble>(1);
     863                                                xs2[counter]=xNew<IssmDouble>(1);
     864                                                *xs[counter]=scalar;
     865                                                *xs2[counter]=pow(scalar,2.0);
     866                                                xsize[counter]=1;
     867                                        }
     868                                        else if (xtype[counter]==3){
     869                                                IssmDouble* doublemat2=xNew<IssmDouble>(doublematsize);
     870                                                for(int k=0;k<doublematsize;k++)doublemat2[k]=pow(doublemat[k],2.0);
     871                                                xs[counter]=doublemat;
     872                                                xs2[counter]=doublemat2;
     873                                                xsize[counter]=doublematsize;
     874                                        }
     875                                        else _error_("cannot carry out statistics on type " << xtype[counter]);
     876                                }
     877                                else{
     878                                        if(xtype[counter]==1){
     879                                                *xs[counter]+=scalar;
     880                                                *xs2[counter]+=pow(scalar,2.0);
     881                                        }
     882                                        else if (xtype[counter]==3){
     883                                                IssmDouble* newdoublemat=xs[counter];
     884                                                IssmDouble* newdoublemat2=xs2[counter];
     885                                                for(int k=0;k<doublematsize;k++){
     886                                                        newdoublemat[k]+=doublemat[k];
     887                                                        newdoublemat2[k]+=pow(doublemat[k],2.0);
     888                                                }
     889                                                xs[counter]=newdoublemat;
     890                                                xs2[counter]=newdoublemat2;
     891                                        }
     892                                        else _error_("cannot carry out statistics on type " << xtype[counter]);
     893                                }
     894                        }
     895                }
     896
     897                /*Deal with time mean: */
     898                for (int f=0;f<nfields;f++){
     899                        char* field=fields[f];
     900                        fseek(fid,0,SEEK_SET);
     901                        meantype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
     902                        if(i==(lower_row+1)){
     903                                if(meantype[f]==1){
     904                                        meanx[f]=xNewZeroInit<IssmDouble>(1);
     905                                        meanx2[f]=xNewZeroInit<IssmDouble>(1);
     906                                        meansize[f]=1;
     907                                }
     908                                else{
     909                                        meanx[f]=xNewZeroInit<IssmDouble>(doublematsize);
     910                                        meanx2[f]=xNewZeroInit<IssmDouble>(doublematsize);
     911                                        meansize[f]=doublematsize;
     912                                }
     913                        }
     914                        fseek(fid,0,SEEK_SET);
     915                        if(meantype[f]==1){
     916                                IssmDouble sc=0;
     917                                IssmDouble sc2=0;
     918                                for(int j=0;j<nsteps;j++){
     919                                        readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
     920                                        sc+=scalar/nsteps;
     921                                }
     922                                sc2+=pow(sc,2.0);
     923                                *meanx[f]+=sc;
     924                                *meanx2[f]+=sc2;
     925                        }
     926                        else{
     927                                IssmDouble* sc=meanx[f];
     928                                IssmDouble* sc2=meanx2[f];
     929                                IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
     930                                IssmDouble* timemean2=xNewZeroInit<IssmDouble>(doublematsize);
     931
     932                                for(int j=0;j<nsteps;j++){
     933                                        readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
     934                                        for (int k=0;k<doublematsize;k++){
     935                                                timemean[k]+=doublemat[k]/nsteps;
     936                                        }
     937                                }
     938                                for (int k=0;k<doublematsize;k++){
     939                                        timemean2[k]=pow(timemean[k],2.0);
     940                                }
     941                                for (int k=0;k<doublematsize;k++){
     942                                        sc[k]+=timemean[k];
     943                                        sc2[k]+=timemean2[k];
     944                                }
     945
     946                        }
     947
     948                }
     949                fclose(fid);
     950
     951                /*delete buffer:*/
     952                xDelete<char>(buffer);
     953        }
     954        ISSM_MPI_Barrier(IssmComm::GetComm());
     955        _printf0_("Done reading files, now computing mean and variance.\n");
     956
     957        /*We have agregated x and x^2 across the cluster, now gather across the cluster onto
     958         *cpu0 and then compute statistics:*/
     959        for (int f=0;f<nfields;f++){
     960                int counter0=f*nsteps+0;
     961                if (xtype[counter0]==1){ /*deal with scalars {{{*/
     962                        IssmDouble mean,stddev;
     963                        for (int j=0;j<nsteps;j++){
     964                                int counter=f*nsteps+j;
     965
     966                                /*we are broadcasting doubles:*/
     967                                IssmDouble scalar=*xs[counter];
     968                                IssmDouble scalar2=*xs2[counter];
     969                                IssmDouble sumscalar;
     970                                IssmDouble sumscalar2;
     971
     972                                ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
     973                                ISSM_MPI_Reduce(&scalar2,&sumscalar2,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
     974                                /*Build average and standard deviation. For standard deviation, use the
     975                                 *following formula: sigma^2=E(x^2)-mu^2:*/
     976                                mean=sumscalar/nsamples;
     977                                stddev=sqrt(sumscalar2/nsamples-pow(mean,2.0));
     978
     979                                /*add to results:*/
     980                                if(my_rank==0){
     981                                        char fieldname[1000];
     982
     983                                        sprintf(fieldname,"%s%s",fields[f],"Mean");
     984                                        results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,j+1,0));
     985                                        sprintf(fieldname,"%s%s",fields[f],"Stddev");
     986                                        results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,j+1,0));
     987                                }
     988
     989                        }
     990                } /*}}}*/
     991                else{ /*deal with arrays:{{{*/
     992
     993                        int size=xsize[counter0];
     994
     995                        IssmDouble*  mean=xNew<IssmDouble>(size);
     996                        IssmDouble*  stddev=xNew<IssmDouble>(size);
     997
     998                        for (int j=0;j<nsteps;j++){
     999                                int counter=f*nsteps+j;
     1000
     1001                                /*we are broadcasting double arrays:*/
     1002                                x=xs[counter];
     1003                                x2=xs2[counter];
     1004
     1005                                IssmDouble*  sumx=xNew<IssmDouble>(size);
     1006                                IssmDouble*  sumx2=xNew<IssmDouble>(size);
     1007
     1008                                ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
     1009                                ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
     1010
     1011                                /*Build average and standard deviation. For standard deviation, use the
     1012                                 *following formula: sigma^2=E(x^2)-mu^2:*/
     1013                                for (int k=0;k<size;k++){
     1014                                        mean[k]=sumx[k]/nsamples;
     1015                                        stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
     1016                                }
     1017
     1018                                /*add to results:*/
     1019                                if(my_rank==0){
     1020                                        char fieldname[1000];
     1021
     1022                                        sprintf(fieldname,"%s%s",fields[f],"Mean");
     1023                                        results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,j+1,0));
     1024                                        sprintf(fieldname,"%s%s",fields[f],"Stddev");
     1025                                        results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,j+1,0));
     1026                                }
     1027                        }
     1028                } /*}}}*/
     1029        }
     1030        /*Do the same but for the time mean:*/
     1031        for (int f=0;f<nfields;f++){
     1032                if (meantype[f]==1){ /*deal with scalars {{{*/
     1033                        IssmDouble mean,stddev;
     1034
     1035                        /*we are broadcasting doubles:*/
     1036                        IssmDouble scalar=*meanx[f];
     1037                        IssmDouble scalar2=*meanx2[f];
     1038                        IssmDouble sumscalar;
     1039                        IssmDouble sumscalar2;
     1040
     1041                        ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
     1042                        ISSM_MPI_Reduce(&scalar2,&sumscalar2,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
     1043                        /*Build average and standard deviation. For standard deviation, use the
     1044                         *following formula: sigma^2=E(x^2)-mu^2:*/
     1045                        mean=sumscalar/nsamples;
     1046                        stddev=sqrt(sumscalar2/nsamples-pow(mean,2.0));
     1047
     1048                        /*add to results:*/
     1049                        if(my_rank==0){
     1050                                char fieldname[1000];
     1051
     1052                                sprintf(fieldname,"%s%s",fields[f],"TimeMean");
     1053                                results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,steps[0],0));
     1054                                sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
     1055                                results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,steps[0],0));
     1056                        }
     1057                } /*}}}*/
     1058                else{ /*deal with arrays:{{{*/
     1059
     1060                        int size=meansize[f];
     1061                        IssmDouble*  mean=xNew<IssmDouble>(size);
     1062                        IssmDouble*  stddev=xNew<IssmDouble>(size);
     1063
     1064                        /*we are broadcasting double arrays:*/
     1065                        x=meanx[f];
     1066                        x2=meanx2[f];
     1067
     1068                        IssmDouble*  sumx=xNew<IssmDouble>(size);
     1069                        IssmDouble*  sumx2=xNew<IssmDouble>(size);
     1070
     1071                        ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
     1072                        ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
     1073
     1074                        /*Build average and standard deviation. For standard deviation, use the
     1075                         *following formula: sigma^2=E(x^2)-mu^2:*/
     1076                        for (int k=0;k<size;k++){
     1077                                mean[k]=sumx[k]/nsamples;
     1078                                stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
     1079                        }
     1080
     1081                        /*add to results:*/
     1082                        if(my_rank==0){
     1083                                char fieldname[1000];
     1084
     1085                                sprintf(fieldname,"%s%s",fields[f],"TimeMean");
     1086                                results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,steps[0],0));
     1087                                sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
     1088                                results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,steps[0],0));
     1089                        }
     1090                } /*}}}*/
     1091        }
     1092
     1093        _printf0_("Done with MeanVariance:\n");
     1094        IssmComm::SetComm(ISSM_MPI_COMM_WORLD);
     1095
     1096        return 1;
     1097} /*}}}*/
     1098int ComputeSampleSeries(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){ /*{{{*/
     1099
     1100        int nsamples;
     1101        char* directory=NULL;
     1102        char* model=NULL;
     1103        char** fields=NULL;
     1104        int* steps=NULL;
     1105        int nsteps;
     1106        int nfields;
     1107        int range,lower_row,upper_row;
     1108        int nfilesperdirectory;
     1109        int* indices=NULL;
     1110        int  nindices;
     1111
     1112        /*intermediary:*/
     1113        IssmDouble* doublemat=NULL;
     1114        int         doublematsize;
     1115        IssmDouble scalar;
     1116
     1117        /*computation of average and variance itself:*/
     1118        IssmDouble*  x = NULL;
     1119        IssmDouble*  allx=NULL;
     1120        IssmDouble** xs = NULL;
     1121        int*         xtype=NULL;
     1122        int*         xsize=NULL;
     1123
     1124        /*only work on the statistical communicator: */
     1125        if (color==MPI_UNDEFINED)return 0;
     1126
     1127        /*Retrieve parameters:*/
     1128        parameters->FindParam(&nsamples,QmuNsampleEnum);
     1129        parameters->FindParam(&directory,DirectoryNameEnum);
     1130        parameters->FindParam(&model,InputFileNameEnum);
     1131        parameters->FindParam(&fields,&nfields,FieldsEnum);
     1132        parameters->FindParam(&steps,&nsteps,StepsEnum);
     1133        parameters->FindParam(&indices,&nindices,IndicesEnum);
     1134
     1135        /*Get rank from the stat comm communicator:*/
     1136        IssmComm::SetComm(statcomm);
     1137        int my_rank=IssmComm::GetRank();
     1138
     1139        /*Open files and read them complelety, in a distributed way:*/
     1140        range=DetermineLocalSize(nsamples,IssmComm::GetComm());
     1141        GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
     1142
     1143        /*Initialize arrays:*/
     1144        xs=xNew<IssmDouble*>(nfields*nsteps);
     1145        xtype=xNew<int>(nfields*nsteps);
     1146        xsize=xNew<int>(nfields*nsteps);
     1147
     1148        /*Start opening files:*/
     1149        for (int i=(lower_row+1);i<=upper_row;i++){
     1150                _printf0_("reading file #: " << i << "\n");
     1151                char file[1000];
     1152                long int  length;
     1153                char* buffer=NULL;
     1154
     1155                /*string:*/
     1156                sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
     1157
     1158                /*open file: */
     1159                _printf0_("    opening file:\n");
     1160                FILE* fid=fopen(file,"rb");
     1161
     1162                /*figure out size of file, and read the whole thing:*/
     1163                _printf0_("    reading file:\n");
     1164                fseek (fid, 0, SEEK_END);
     1165                length = ftell (fid);
     1166                fseek (fid, 0, SEEK_SET);
     1167                buffer = xNew<char>(length);
     1168                fread (buffer, sizeof(char), length, fid);
     1169
     1170                /*close file:*/
     1171                fclose (fid);
     1172
     1173                /*create a memory stream with this buffer:*/
     1174                _printf0_("    processing file:\n");
     1175                fid=fmemopen(buffer, length, "rb");
     1176
     1177                /*start reading data from the buffer directly:*/
     1178                for (int f=0;f<nfields;f++){
     1179                        fseek(fid,0,SEEK_SET);
     1180                        char* field=fields[f];
     1181                        for (int j=0;j<nsteps;j++){
     1182                                int counter=f*nsteps+j;
     1183                                xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
     1184                                if(i==(lower_row+1)){
     1185                                        if(xtype[counter]==1){
     1186                                                x=xNew<IssmDouble>(range);
     1187                                                x[0]=scalar;
     1188                                                xs[counter]=x;
     1189                                                xsize[counter]=range;
     1190                                        }
     1191                                        else if (xtype[counter]==3){
     1192                                                x=xNew<IssmDouble>(nindices*range);
     1193                                                for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
     1194                                                xs[counter]=x;
     1195                                                xsize[counter]=range*nindices;
     1196                                        }
     1197                                        else _error_("cannot carry out statistics on type " << xtype[counter]);
     1198                                }
     1199                                else{
     1200                                        if(xtype[counter]==1){
     1201                                                x=xs[counter];
     1202                                                x[i-(lower_row+1)]=scalar;
     1203                                                xs[counter]=x;
     1204                                        }
     1205                                        else if (xtype[counter]==3){
     1206                                                x=xs[counter];
     1207                                                for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
     1208                                                xs[counter]=x;
     1209                                        }
     1210                                        else _error_("cannot carry out statistics on type " << xtype[counter]);
     1211                                }
     1212                        }
     1213                }
     1214                fclose(fid);
     1215
     1216                /*delete buffer:*/
     1217                xDelete<char>(buffer);
     1218        }
     1219        ISSM_MPI_Barrier(IssmComm::GetComm());
     1220        _printf0_("Done reading files, now assembling time series.\n");
     1221
     1222        for (int f=0;f<nfields;f++){
     1223                for (int j=0;j<nsteps;j++){
     1224                        int counter=f*nsteps+j;
     1225                        if (xtype[counter]==1){
     1226                                /*we are broadcasting range times doubles:*/
     1227                                x=xs[counter];
     1228                                allx=xNew<IssmDouble>(nsamples);
     1229                                MPI_Gather(x, range, ISSM_MPI_PDOUBLE,allx, range, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
     1230                                /*add to results:*/
     1231                                if(my_rank==0){
     1232                                        char fieldname[1000];
     1233
     1234                                        sprintf(fieldname,"%s%s",fields[f],"Samples");
     1235                                        results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,1,j+1,0));
     1236                                }
     1237                        }
     1238                        else{
     1239                                /*we are broadcasting double arrays:*/
     1240                                x=xs[counter];
     1241                                allx=xNew<IssmDouble>(nsamples*nindices);
     1242
     1243                                MPI_Gather(x, range*nindices, ISSM_MPI_PDOUBLE,allx, range*nindices, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
     1244
     1245                                /*add to results:*/
     1246                                if(my_rank==0){
     1247                                        char fieldname[1000];
     1248                                        sprintf(fieldname,"%s%s",fields[f],"Samples");
     1249                                        results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,nindices,j+1,0));
     1250                                }
     1251                        }
     1252                }
     1253        }
     1254        _printf0_("Done with SampleSeries:\n");
     1255        IssmComm::SetComm(ISSM_MPI_COMM_WORLD);
     1256
     1257        return 1;
     1258} /*}}}*/
     1259int OutputStatistics(Parameters* parameters,Results* results,int color,ISSM_MPI_Comm statcomm){ /*{{{*/
     1260
     1261        char   outputfilename[1000];
     1262        char* directory=NULL;
     1263        char* model=NULL;
     1264        char* method=NULL;
     1265        int   nsamples;
     1266        int* steps=NULL;
     1267        int nsteps;
     1268
     1269        /*only work on the statistical communicator: */
     1270        if (color==MPI_UNDEFINED)return 0;
     1271
     1272        FemModel* femmodel=new FemModel();
     1273
     1274        /*Some parameters that will allow us to use the OutputResultsx module:*/
     1275        parameters->AddObject(new BoolParam(QmuIsdakotaEnum,false));
     1276        parameters->AddObject(new BoolParam(SettingsIoGatherEnum,true));
     1277
     1278        parameters->FindParam(&directory,DirectoryNameEnum);
     1279        parameters->FindParam(&model,InputFileNameEnum);
     1280        parameters->FindParam(&nsamples,QmuNsampleEnum);
     1281        parameters->FindParam(&steps,&nsteps,StepsEnum);
     1282
     1283        sprintf(outputfilename,"%s/%s.stats",directory,model);
     1284        parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
     1285
     1286        /*Call OutputResults module:*/
     1287        femmodel->parameters=parameters;
     1288        femmodel->results=results;
     1289
     1290        OutputResultsx(femmodel);
     1291
     1292        return 1;
     1293} /*}}}*/
     1294int readdata(IssmDouble** pdoublemat, int* pdoublematsize, IssmDouble* pdouble, FILE* fid,char* field,int step){ /*{{{*/
     1295
     1296        int length;
     1297        char fieldname[1000];
     1298        int   fieldname_size;
     1299        IssmDouble   rtime;
     1300        int          rstep;
     1301        int M,N;
     1302
     1303        //fields that we retrive:
     1304        IssmDouble  dfield;
     1305        char*       sfield    = NULL;
     1306        IssmDouble* dmatfield = NULL;
     1307        int*        imatfield = NULL;
     1308
     1309        //type of the returned field:
     1310        int type;
     1311        int found=0;
     1312
     1313        while(1){
     1314
     1315                size_t ret_code = fread(&fieldname_size, sizeof(int), 1, fid);
     1316                if(ret_code != 1) break; //we are done.
     1317
     1318                fread(fieldname, sizeof(char), fieldname_size, fid);
     1319                //_printf0_("fieldname: " << fieldname << "\n");
     1320
     1321                fread(&rtime, sizeof(IssmDouble), 1, fid);
     1322                fread(&rstep, sizeof(int), 1, fid);
     1323
     1324                //check on field:
     1325                if ((step==rstep) && (strcmp(field,fieldname)==0)){
     1326
     1327                        //ok, go read the result really:
     1328                        fread(&type,sizeof(int),1,fid);
     1329                        fread(&M,sizeof(int),1,fid);
     1330                        if (type==1){
     1331                                fread(&dfield,sizeof(IssmDouble),1,fid);
     1332                        }
     1333                        else if (type==2){
     1334                                fread(&M,sizeof(int),1,fid);
     1335                                sfield=xNew<char>(M);
     1336                                fread(sfield,sizeof(char),M,fid);
     1337                        }
     1338                        else if (type==3){
     1339                                fread(&N,sizeof(int),1,fid);
     1340                                dmatfield=xNew<IssmDouble>(M*N);
     1341                                fread(dmatfield,sizeof(IssmDouble),M*N,fid);
     1342                        }
     1343                        else if (type==4){
     1344                                fread(&N,sizeof(int),1,fid);
     1345                                imatfield=xNew<int>(M*N);
     1346                                fread(imatfield,sizeof(int),M*N,fid);
     1347                        }
     1348                        else _error_("cannot read data of type " << type << "\n");
     1349                        found=1;
     1350                        break;
     1351                }
     1352                else{
     1353                        //just skim to next results.
     1354                        fread(&type,sizeof(int),1,fid);
     1355                        fread(&M,sizeof(int),1,fid);
     1356                        if (type==1){
     1357                                fseek(fid,sizeof(IssmDouble),SEEK_CUR);
     1358                        }
     1359                        else if(type==2){
     1360                                fseek(fid,M*sizeof(char),SEEK_CUR);
     1361                        }
     1362                        else if(type==3){
     1363                                fread(&N,sizeof(int),1,fid);
     1364                                fseek(fid,M*N*sizeof(IssmDouble),SEEK_CUR);
     1365                        }
     1366                        else if(type==4){
     1367                                fread(&N,sizeof(int),1,fid);
     1368                                fseek(fid,M*N*sizeof(int),SEEK_CUR);
     1369                        }
     1370                        else _error_("cannot read data of type " << type << "\n");
     1371                }
     1372        }
     1373        if(found==0)_error_("cound not find " << field << " at step " << step  << "\n");
     1374
     1375        /*assign output pointers:*/
     1376        *pdoublemat=dmatfield;
     1377        *pdoublematsize=M*N;
     1378        *pdouble=dfield;
     1379
     1380        /*return:*/
     1381        return type;
     1382
     1383}
     1384/*}}}*/
     1385bool DakotaDirStructure(int argc,char** argv){ /*{{{*/
     1386
     1387        char* input_file;
     1388        FILE* fid;
     1389        IoModel* iomodel=NULL;
     1390        int check;
     1391
     1392        //qmu statistics
     1393        bool statistics    = false;
     1394        int  numdirectories = 0;
     1395
     1396        /*First things first, set the communicator as a global variable: */
     1397        IssmComm::SetComm(MPI_COMM_WORLD);
     1398
     1399        /*Barrier:*/
     1400        ISSM_MPI_Barrier(IssmComm::GetComm());
     1401        _printf0_("Preparing directory structure for model outputs:" << "\n");
     1402
     1403        //open model input file for reading
     1404        input_file=xNew<char>((strlen(argv[2])+strlen(argv[3])+strlen(".bin")+2));
     1405        sprintf(input_file,"%s/%s%s",argv[2],argv[3],".bin");
     1406        fid=fopen(input_file,"rb");
     1407        if (fid==NULL) Cerr << "dirstructure error message: could not open model " << input_file << " to retrieve qmu statistics parameters" << std::endl;
     1408
     1409        //initialize IoModel, but light version, we just need it to fetch one constant:
     1410        iomodel=new IoModel();
     1411        iomodel->fid=fid;
     1412        iomodel->FetchConstants();
     1413
     1414        //early return if statistics not requested:
     1415        iomodel->FindConstant(&statistics,"md.qmu.statistics");
     1416        if(!statistics){
     1417                delete iomodel;
     1418                xDelete<char>(input_file);
     1419                fclose(fid);
     1420                return false; //important return value!
     1421        }
     1422
     1423        iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories");
     1424
     1425        /*Ok, we have everything we need to create the directory structure:*/
     1426        if(IssmComm::GetRank()==0){
     1427                for (int i=0;i<numdirectories;i++){
     1428                        char directory[1000];
     1429                        sprintf(directory,"./%i",i+1);
     1430
     1431                        check = mkdir(directory,ACCESSPERMS);
     1432                        if (check) _error_("dirstructure error message: could not create directory " << directory << "\n");
     1433                }
     1434        }
     1435
     1436        /*Delete resources:*/
     1437        delete iomodel;
     1438        xDelete<char>(input_file);
     1439
     1440        //close model file:
     1441        fclose(fid);
     1442
     1443        //return value:
     1444        return true; //statistics computation on!
     1445} /*}}}*/
Note: See TracChangeset for help on using the changeset viewer.