Changeset 25480

Timestamp:

08/26/20 18:47:35 (5 years ago)

Author:

Eric.Larour

Message:

CHG: starting implementation of a histogram capability in issm_post.

File:

• issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp (modified) (10 diffs)

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

@@ -27 +27 @@
         int    nindices;
         int*   indices=NULL;
+        int    nbins;
         
         /*datasets*/
@@ -92 +93 @@
 
                 ComputeMeanVariance(parameters,results);
+        }
+        else if(strcmp(method,"Histogram")==0){
+                
+                /*Retrieve the size of the histogram (number of bins):*/
+                nbins=atoi(argv[offset]); offset++;
+                parameters->AddObject(new IntParam(NbinsEnum,nbins));
+                
+                ComputeHistogram(parameters,results);
+
         }
         else if(strcmp(method,"SampleSeries")==0){
@@ -244 +254 @@
         _printf0_("Done reading files, now computing mean and variance.\n"); 
 
-        /*We have agregated x and x^2 across the cluster, now consolidate: */
+        /*We have agregated x and x^2 across the cluster, now gather across the cluster onto
+         *cpu0 and then compute statistics:*/
         for (int f=0;f<nfields;f++){
-                for (int j=0;j<nsteps;j++){
-                        int counter=f*nsteps+j;
-                        if (xtype[counter]==1){
+
+                IssmDouble sumscalar_alltimes=0;
+                IssmDouble sumscalar2_alltimes=0;
+                IssmDouble* sumarray_alltimes=NULL;
+                IssmDouble* sumarray2_alltimes=NULL;
+                        
+                int counter0=f*nsteps+0;
+                if (xtype[counter0]==1){ /*deal with scalars {{{*/
+                        IssmDouble mean,stddev;
+                        for (int j=0;j<nsteps;j++){
+                                int counter=f*nsteps+j;
+
                                 /*we are broadcasting doubles:*/
                                 IssmDouble scalar=*xs[counter];
@@ -254 +274 @@
                                 IssmDouble sumscalar;
                                 IssmDouble sumscalar2;
-                                IssmDouble mean,stddev;
 
                                 ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
@@ -266 +285 @@
                                 if(my_rank==0){
                                         char fieldname[1000];
-                                        
+
                                         sprintf(fieldname,"%s%s",fields[f],"Mean");
                                         results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,j+1,0));
@@ -272 +291 @@
                                         results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,j+1,0));
                                 }
-                        }
-                        else{
+
+                                /*We are also looking for the mean and stddev  over all time steps, 
+                                 *so accumulate sums again:*/
+                                sumscalar_alltimes+=sumscalar;
+                                sumscalar2_alltimes+=sumscalar2;
+                        }
+                        /*Build averae and standard deviation for all time steps:*/
+                        mean=sumscalar_alltimes/nsamples/nsteps;
+                        stddev=sqrt(sumscalar2_alltimes/nsamples/nsteps-pow(mean,2.0));
+                        /*add to results at step 1:*/
+                        if(my_rank==0){
+                                char fieldname[1000];
+
+                                sprintf(fieldname,"%s%s",fields[f],"TimeMean");
+                                results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,1,0));
+                                sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
+                                results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,1,0));
+                        }
+                } /*}}}*/
+                else{ /*deal with arrays:{{{*/
+
+                        int size=xsize[counter0];
+                        sumarray_alltimes=xNewZeroInit<IssmDouble>(size);
+                        sumarray2_alltimes=xNewZeroInit<IssmDouble>(size);
+
+                        IssmDouble*  mean=xNew<IssmDouble>(size);
+                        IssmDouble*  stddev=xNew<IssmDouble>(size);
+
+                        for (int j=0;j<nsteps;j++){
+                                int counter=f*nsteps+j;
+
                                 /*we are broadcasting double arrays:*/
                                 x=xs[counter];
                                 x2=xs2[counter];
-                                int size=xsize[counter];
-                                
+
                                 IssmDouble*  sumx=xNew<IssmDouble>(size);
                                 IssmDouble*  sumx2=xNew<IssmDouble>(size);
-                                IssmDouble*  mean=xNew<IssmDouble>(size);
-                                IssmDouble*  stddev=xNew<IssmDouble>(size);
-
+                                
                                 ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
                                 ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
@@ -303 +348 @@
                                         results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,j+1,0));
                                 }
-                        }
-                }
-        }
-
+
+                                /*We want mean and stddev across time too, so accumulate through time: */
+                                for (int k=0;k<size;k++){
+                                        sumarray_alltimes[k]+=sumx[k];
+                                        sumarray2_alltimes[k]+=sumx2[k];
+                                }
+                        }
+
+                        /*build time mean and stddev:*/
+                        for (int k=0;k<size;k++){
+                                mean[k]=sumarray_alltimes[k]/nsamples/nsteps;
+                                stddev[k]=sqrt(sumarray2_alltimes[k]/nsamples/nsteps-pow(mean[k],2.0));
+                        }
+
+                        /*add to results at step 1:*/
+                        if(my_rank==0){
+                                char fieldname[1000];
+
+                                sprintf(fieldname,"%s%s",fields[f],"TimeMean");
+                                results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,1,0));
+                                sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
+                                results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,1,0));
+                        }
+                } /*}}}*/
+
+        }
 
 } /*}}}*/
@@ -462 +529 @@
 
 } /*}}}*/
+int ComputeHistogram(Parameters* parameters,Results* results){ /*{{{*/
+
+        int nsamples; 
+        char* directory=NULL;
+        char* model=NULL;
+        char** fields=NULL;
+        int* steps=NULL;
+        int nsteps;
+        int nfields;
+        int nbins;
+        int range,lower_row,upper_row;
+        
+        /*intermediary:*/
+        IssmDouble* doublemat=NULL;
+        int         doublematsize;
+        IssmDouble scalar;
+
+        /*computation of average and variance itself:*/
+        IssmDouble** maxxs = NULL;
+        IssmDouble** minxs = NULL;
+        int*         xtype=NULL;
+        int*         xsize=NULL;
+
+        IssmDouble** maxmeans=NULL;
+        IssmDouble** minmeans=NULL;
+        int*         meanxtype=NULL;
+        int*         meanxsize=NULL;
+
+        /*Retrieve parameters:*/
+        parameters->FindParam(&nsamples,QmuNsampleEnum);
+        parameters->FindParam(&directory,DirectoryNameEnum);
+        parameters->FindParam(&model,InputFileNameEnum);
+        parameters->FindParam(&fields,&nfields,FieldsEnum);
+        parameters->FindParam(&steps,&nsteps,StepsEnum);
+        parameters->FindParam(&nbins,NbinsEnum);
+
+        /*Get rank:*/
+        int my_rank=IssmComm::GetRank();
+
+        /*Open files and read them complelety, in a distributed way:*/
+        range=DetermineLocalSize(nsamples,IssmComm::GetComm());
+        GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
+
+        /*Initialize arrays:*/
+
+        maxmeans=xNew<IssmDouble*>(nfields);
+        minmeans=xNew<IssmDouble*>(nfields);
+        meanxtype=xNew<IssmDouble*>(nfields);
+        meanxsize=xNew<IssmDouble*>(nfields);
+
+        maxxs=xNew<IssmDouble*>(nfields*nsteps);
+        minxs=xNew<IssmDouble*>(nfields*nsteps);
+        xtype=xNew<int>(nfields*nsteps);
+        xsize=xNew<int>(nfields*nsteps);
+
+        /*Start opening files:*/
+        for (int i=(lower_row+1);i<=upper_row;i++){
+                _printf0_("reading file #: " << i << "\n");
+                char file[1000];
+                int  length;
+                char* buffer=NULL;
+
+                /*string:*/
+                sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
+
+                /*open file: */
+                _printf0_("    opening file:\n");
+                FILE* fid=fopen(file,"rb");
+
+                /*figure out size of file, and read the whole thing:*/
+                _printf0_("    reading file:\n");
+                fseek (fid, 0, SEEK_END);
+                length = ftell (fid);
+                fseek (fid, 0, SEEK_SET);
+                buffer = xNew<char>(length);
+                fread (buffer, sizeof(char), length, fid);
+
+                /*close file:*/
+                fclose (fid);
+
+                /*create a memory stream with this buffer:*/
+                _printf0_("    processing file:\n");
+                fid=fmemopen(buffer, length, "rb");
+
+                /*start reading data from the buffer directly:*/
+                for (int f=0;f<nfields;f++){
+                        char* field=fields[f];
+                        for (int j=0;j<nsteps;j++){
+                                int counter=f*nsteps+j;
+                                xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
+                                if(i==(lower_row+1)){
+                                        if(xtype[counter]==1){
+                                                maxxs[counter]=xNew<IssmDouble>(1); 
+                                                minxs[counter]=xNew<IssmDouble>(1); 
+                                                *maxxs[counter]=scalar;
+                                                *minxs[counter]=scalar;
+                                                xsize[counter]=1;
+                                        }
+                                        else if (xtype[counter]==3){
+                                                maxxs[counter]=doublemat;
+                                                minxs[counter]=doublemat;
+                                                xsize[counter]=doublematsize;
+                                        }
+                                        else _error_("cannot carry out statistics on type " << xtype[counter]); 
+                                }
+                                else{
+                                        if(xtype[counter]==1){
+                                                *maxxs[counter]=max(*maxxs[counter],scalar);
+                                                *minxs[counter]=min(*minxs[counter],scalar);
+                                        }
+                                        else if (xtype[counter]==3){
+                                                IssmDouble* newmax=minxs[counter];
+                                                IssmDouble* newmin=maxxs[counter];
+                                                for(int k=0;k<doublematsize;k++){
+                                                        newmax[k]=max(newmax[k],doublemat[k]);
+                                                        newmin[k]=min(newmin[k],doublemat[k]);
+                                                }
+                                                maxxs[counter]=newmax;
+                                                minxs[counter]=newmin;
+                                        }
+                                        else _error_("cannot carry out statistics on type " << xtype[counter]); 
+                                }
+                        }
+                }
+                
+                /*Deal with average in time: */
+                for (int f=0;f<nfields;f++){
+                        char* field=fields[f];
+                        meanxtype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
+                        if(i==(lower_row+1)){
+                                if(meanxtype[f]==1){
+                                        maxmeans[f]=xNewZeroInit<IssmDouble>(1);
+                                        minmeans[f]=xNewZeroInit<IssmDouble>(1);
+                                        meanxsize[f]=1;
+                                }
+                                else{
+                                        maxmeans[f]=xNewZeroInit<IssmDouble>(doublematsize);
+                                        minmeans[f]=xNewZeroInit<IssmDouble>(doublematsize);
+                                        meanxsize[f]=doublematsize;
+                                }
+
+                        }
+                        if(meanxtype[f]==1){
+                                IssmDouble timemean=0;
+                                for (int j=0;j<nsteps;j++){
+                                        readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
+                                        timemean+=scalar/nsteps;
+                                }
+                                /*Figure out max and min of time means: */
+                                *maxmeans[f]=max(*maxmeans[f],timemean);
+                                *minmeans[f]=min(*minmeans[f],timemean);
+                        }
+                        else{
+                                IssmDouble* maxx=maxmeans[f];
+                                IssmDouble* minx=minmeans[f];
+                                for(int k=0;k<doublematsize;k++){
+                                        maxx[k]=max(maxx[k],doublemat[k]);
+                                        minx[k]=min(minx[k],doublemat[k]);
+                                }
+                                maxmeans[f]=maxx;
+                                minmeans[f]=minx;
+                                }
+                }
+                fclose(fid);
+
+                /*delete buffer:*/
+                xDelete<char>(buffer);
+        }
+        ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); 
+        _printf0_("Done reading files, now computing min and max.\n"); 
+
+        /*We have agregated minx and max across the cluster, now gather across the cluster onto
+         *cpu0 and then compute statistics:*/
+        for (int f=0;f<nfields;f++){
+                int counter0=f*nsteps+0;
+                if (xtype[counter0]==1){ /*deal with scalars {{{*/
+                        for (int j=0;j<nsteps;j++){
+                                int counter=f*nsteps+j;
+
+                                /*we are broadcasting doubles:*/
+                                IssmDouble maxscalar=*maxxs[counter];
+                                IssmDouble minscalar=*minxs[counter];
+                                IssmDouble allmaxscalar;
+                                IssmDouble allminscalar;
+                                IssmDouble sumscalar_alltimes=0;
+
+                                ISSM_MPI_Reduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm());
+                                ISSM_MPI_Reduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,0,IssmComm::GetComm());
+
+                                /*add to results:*/
+                                if(my_rank==0){
+                                        char fieldname[1000];
+
+                                        sprintf(fieldname,"%s%s",fields[f],"Max");
+                                        results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,allmaxscalar,j+1,0));
+                                        sprintf(fieldname,"%s%s",fields[f],"Min");
+                                        results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,allminscalar,j+1,0));
+                                }
+                        }
+                } /*}}}*/
+                else{ /*deal with arrays:{{{*/
+
+                        int size=xsize[counter0];
+                        for (int j=0;j<nsteps;j++){
+                                int counter=f*nsteps+j;
+
+                                /*we are broadcasting double arrays:*/
+                                IssmDouble* maxx==maxxs[counter];
+                                IssmDouble* minx==minxs[counter];
+
+                                IssmDouble*  allmax=xNew<IssmDouble>(size);
+                                IssmDouble*  allmin=xNew<IssmDouble>(size);
+                                
+                                ISSM_MPI_Reduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm());
+                                ISSM_MPI_Reduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,0,IssmComm::GetComm());
+
+                                /*add to results:*/
+                                if(my_rank==0){
+                                        char fieldname[1000];
+
+                                        sprintf(fieldname,"%s%s",fields[f],"Max");
+                                        results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allmax,size,1,j+1,0));
+                                        sprintf(fieldname,"%s%s",fields[f],"Min");
+                                        results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allmin,size,1,j+1,0));
+                                }
+                        }
+                } /*}}}*/
+        }
+        
+        /*Now do the same for the time mean fields:*/
+        for (int f=0;f<nfields;f++){
+                if (meanxtype[f]==1){ /*deal with scalars {{{*/
+
+                        /*we are broadcasting doubles:*/
+                        IssmDouble maxscalar=*maxmeans[f];
+                        IssmDouble minscalar=*minmeans[f];
+                        IssmDouble allmaxscalar;
+                        IssmDouble allminscalar;
+
+                        ISSM_MPI_Reduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm());
+                        ISSM_MPI_Reduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,0,IssmComm::GetComm());
+
+                        /*add to results at time step 1:*/
+                        if(my_rank==0){
+                                char fieldname[1000];
+
+                                sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax");
+                                results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,allmaxscalar,1,0));
+                                sprintf(fieldname,"%s%s",fields[f],"TimeMeaMin");
+                                results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,allminscalar,1,0));
+                        }
+                } /*}}}*/
+                else{ /*deal with arrays:{{{*/
+
+                        int size=meanxsize[f];
+
+                        /*we are broadcasting double arrays:*/
+                        IssmDouble* maxx==maxmeans[f];
+                        IssmDouble* minx==minmeans[f];
+
+                        IssmDouble*  allmax=xNew<IssmDouble>(size);
+                        IssmDouble*  allmin=xNew<IssmDouble>(size);
+
+                        ISSM_MPI_Reduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm());
+                        ISSM_MPI_Reduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,0,IssmComm::GetComm());
+
+                        /*add to results at step 1:*/
+                        if(my_rank==0){
+                                char fieldname[1000];
+
+                                sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax");
+                                results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allmax,size,1,j,0));
+                                sprintf(fieldname,"%s%s",fields[f],"TimeMeanMin");
+                                results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allmin,size,1,j,0));
+                        }
+                } /*}}}*/
+        }
+
+
+} /*}}}*/
 int OutputStatistics(Parameters* parameters,Results* results){ /*{{{*/
         
@@ -469 +816 @@
         char* method=NULL;
         int   nsamples;
+        int* steps=NULL;
+        int nsteps;
+
         FemModel* femmodel=new FemModel();
-
         
         /*Some parameters that will allow us to use the OutputResultsx module:*/
@@ -480 +829 @@
         parameters->FindParam(&method,AnalysisTypeEnum);
         parameters->FindParam(&nsamples,QmuNsampleEnum);
-        sprintf(outputfilename,"%s/%s-%s.stats-%i",directory,model,method,nsamples);
+        parameters->FindParam(&steps,&nsteps,StepsEnum);
+
+        sprintf(outputfilename,"%s/%s-%s.stats-samp%i-time%i:%i",directory,model,method,nsamples,steps[0],steps[nsteps-1]);
         parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));