Index: /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp
===================================================================
--- /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp	(revision 25493)
+++ /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp	(revision 25494)
@@ -10,4 +10,5 @@
 	   	     IssmDouble* pdouble, FILE* fid,char* field,int step);
 int OutputStatistics(Parameters* parameters,Results* results);
+int ComputeHistogram(Parameters* parameters,Results* results);
 
 int main(int argc,char **argv){ /*{{{*/
@@ -576,6 +577,6 @@
 	maxmeans=xNew<IssmDouble*>(nfields);
 	minmeans=xNew<IssmDouble*>(nfields);
-	meanxtype=xNew<IssmDouble*>(nfields);
-	meanxsize=xNew<IssmDouble*>(nfields);
+	meanxtype=xNew<int>(nfields);
+	meanxsize=xNew<int>(nfields);
 
 	maxxs=xNew<IssmDouble*>(nfields*nsteps);
@@ -597,4 +598,5 @@
 		_printf0_("    opening file:\n");
 		FILE* fid=fopen(file,"rb");
+		if(fid==NULL)_error_("cound not open file: " << file << "\n");
 
 		/*figure out size of file, and read the whole thing:*/
@@ -616,4 +618,5 @@
 		for (int f=0;f<nfields;f++){
 			char* field=fields[f];
+			fseek(fid,0,SEEK_SET);
 			for (int j=0;j<nsteps;j++){
 				int counter=f*nsteps+j;
@@ -653,42 +656,58 @@
 			}
 		}
-		
+		_printf0_("    average in time:\n");
+
 		/*Deal with average in time: */
+		fseek(fid,0,SEEK_SET);
 		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;
+				fseek(fid,0,SEEK_SET);
 				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);
+				if(i==(lower_row+1)){
+					maxmeans[f]=xNewZeroInit<IssmDouble>(1); 
+					minmeans[f]=xNewZeroInit<IssmDouble>(1); 
+					*maxmeans[f]=timemean;
+					*minmeans[f]=timemean;
+				}
+				else{
+					*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;
-				}
+				fseek(fid,0,SEEK_SET);
+				IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
+				for (int j=0;j<nsteps;j++){
+					readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
+					for (int k=0;i<doublematsize;k++){
+						timemean[k]+=doublemat[k]/nsteps;
+					}
+				}
+
+				if(i==(lower_row+1)){
+					maxmeans[f]=timemean;
+					minmeans[f]=timemean;
+				}
+				else{
+					IssmDouble* maxx=maxmeans[f];
+					IssmDouble* minx=minmeans[f];
+
+					for(int k=0;k<doublematsize;k++){
+						maxx[k]=max(maxx[k],timemean[k]);
+						minx[k]=min(minx[k],timemean[k]);
+					}
+					maxmeans[f]=maxx;
+					minmeans[f]=minx;
+				}
+			}
 		}
 		fclose(fid);
@@ -715,6 +734,6 @@
 				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());
+				ISSM_MPI_Allreduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
+				ISSM_MPI_Allreduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
 
 				/*add to results:*/
@@ -727,4 +746,9 @@
 					results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,allminscalar,j+1,0));
 				}
+
+				/*Store broadcasted value for later computation of histograms:*/
+				*maxxs[counter]=allmaxscalar;
+				*minxs[counter]=allminscalar;
+
 			}
 		} /*}}}*/
@@ -736,12 +760,12 @@
 
 				/*we are broadcasting double arrays:*/
-				IssmDouble* maxx==maxxs[counter];
-				IssmDouble* minx==minxs[counter];
+				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());
+				ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
+				ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
 
 				/*add to results:*/
@@ -754,4 +778,7 @@
 					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allmin,size,1,j+1,0));
 				}
+				/*Store broadcasted value for later computation of histograms:*/
+				maxxs[counter]=allmax;
+				minxs[counter]=allmin;
 			}
 		} /*}}}*/
@@ -768,6 +795,6 @@
 			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());
+			ISSM_MPI_Allreduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
+			ISSM_MPI_Allreduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
 
 			/*add to results at time step 1:*/
@@ -780,4 +807,8 @@
 				results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,allminscalar,1,0));
 			}
+			/*Store for later use in histogram computation:*/
+			*maxmeans[f]=allmaxscalar;
+			*minmeans[f]=allminscalar;
+
 		} /*}}}*/
 		else{ /*deal with arrays:{{{*/
@@ -786,12 +817,12 @@
 
 			/*we are broadcasting double arrays:*/
-			IssmDouble* maxx==maxmeans[f];
-			IssmDouble* minx==minmeans[f];
+			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());
+			ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
+			ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
 
 			/*add to results at step 1:*/
@@ -800,12 +831,260 @@
 
 				sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax");
-				results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allmax,size,1,j,0));
+				results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allmax,size,1,1,0));
 				sprintf(fieldname,"%s%s",fields[f],"TimeMeanMin");
-				results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allmin,size,1,j,0));
-			}
+				results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allmin,size,1,1,0));
+			}
+			
+			/*Store for later use in histogram computation:*/
+			maxmeans[f]=allmax;
+			minmeans[f]=allmin;
+
 		} /*}}}*/
 	}
 
-
+	/*Now that we have the min and max, we can start binning. First allocate 
+	 * histograms, then start filling them:*/
+	IssmDouble** histogram=xNew<IssmDouble*>(nfields*nsteps);
+	IssmDouble** timehistogram=xNew<IssmDouble*>(nfields);
+
+	_printf0_("Start reading files again, this time binning values in the histogram:\n");
+	/*Start opening files:*/
+	for (int i=(lower_row+1);i<=upper_row;i++){
+		_printf0_("reading file #: " << i << "\n");
+		char file[1000];
+		int  length;
+		char* buffer=NULL;
+
+		/*string:*/
+		sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
+
+		/*open file: */
+		_printf0_("    opening file:\n");
+		FILE* fid=fopen(file,"rb");
+		if(fid==NULL)_error_("cound not open file: " << file << "\n");
+
+		/*figure out size of file, and read the whole thing:*/
+		_printf0_("    reading file:\n");
+		fseek (fid, 0, SEEK_END);
+		length = ftell (fid);
+		fseek (fid, 0, SEEK_SET);
+		buffer = xNew<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];
+			fseek(fid,0,SEEK_SET);
+			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){
+						IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(nbins);
+						IssmDouble ma=*maxxs[counter];
+						IssmDouble mi=*minxs[counter];
+						int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
+						localhistogram[index]++;
+						histogram[counter]=localhistogram;
+					}
+					else if (xtype[counter]==3){
+						IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(doublematsize*nbins);
+						IssmDouble* ma=maxxs[counter];
+						IssmDouble* mi=minxs[counter];
+						for (int k=0;k<doublematsize;k++){
+							IssmDouble scalar=doublemat[k];
+							int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
+							localhistogram[k*nbins+index]++;
+						}
+						histogram[counter]=localhistogram;
+					}
+					else _error_("cannot carry out statistics on type " << xtype[counter]); 
+				}
+				else{
+					if(xtype[counter]==1){
+						IssmDouble* localhistogram=histogram[counter];
+						IssmDouble ma=*maxxs[counter];
+						IssmDouble mi=*minxs[counter];
+						int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
+						localhistogram[index]++;
+					}
+					else if (xtype[counter]==3){
+						IssmDouble* localhistogram=histogram[counter];
+						IssmDouble* ma=maxxs[counter];
+						IssmDouble* mi=minxs[counter];
+						for (int k=0;k<doublematsize;k++){
+							IssmDouble scalar=doublemat[k];
+							int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
+							localhistogram[k*nbins+index]++;
+						}
+					}
+					else _error_("cannot carry out statistics on type " << xtype[counter]); 
+				}
+			}
+		}
+		_printf0_("    average in time:\n");
+
+		/*Deal with average in time: */
+		fseek(fid,0,SEEK_SET);
+		for (int f=0;f<nfields;f++){
+			char* field=fields[f];
+			meanxtype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
+			
+			if(meanxtype[f]==1){
+				IssmDouble timemean=0;
+				fseek(fid,0,SEEK_SET);
+				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: */
+				if(i==(lower_row+1)){
+					IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(nbins); 
+					IssmDouble ma=*maxmeans[f];
+					IssmDouble mi=*minmeans[f];
+					int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
+					localhistogram[index]++;
+					timehistogram[f]=localhistogram;
+				}
+				else{
+					IssmDouble* localhistogram=timehistogram[f];
+					IssmDouble ma=*maxmeans[f];
+					IssmDouble mi=*minmeans[f];
+					int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
+					localhistogram[index]++;
+				}
+			}
+			else{
+				fseek(fid,0,SEEK_SET);
+				IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
+				for (int j=0;j<nsteps;j++){
+					readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
+					for (int k=0;i<doublematsize;k++){
+						timemean[k]+=doublemat[k]/nsteps;
+					}
+				}
+
+				if(i==(lower_row+1)){
+					IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(doublematsize*nbins);
+					IssmDouble* ma=maxmeans[f];
+					IssmDouble* mi=minmeans[f];
+					
+					for (int k=0;k<doublematsize;k++){
+						IssmDouble scalar=timemean[k];
+						int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
+						localhistogram[k*nbins+index]++;
+					}
+					timehistogram[f]=localhistogram;
+				}
+				else{
+					
+					IssmDouble* localhistogram=timehistogram[f];
+					IssmDouble* ma=maxmeans[f];
+					IssmDouble* mi=minmeans[f];
+					
+					for (int k=0;k<doublematsize;k++){
+						IssmDouble scalar=timemean[k];
+						int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
+						localhistogram[k*nbins+index]++;
+					}
+				}
+			}
+		}
+		fclose(fid);
+
+		/*delete buffer:*/
+		xDelete<char>(buffer);
+	}
+
+	/*We have agregated histograms across the cluster, now gather them across  the cluster onto
+	 *cpu0: */
+	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* histo=histogram[counter]; //size nbins
+				IssmDouble* allhisto=xNewZeroInit<IssmDouble>(nbins);
+
+				ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
+
+				/*add to results:*/
+				if(my_rank==0){
+					char fieldname[1000];
+
+					sprintf(fieldname,"%s%s",fields[f],"Histogram");
+					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,1,nbins,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* histo=histogram[counter];
+				IssmDouble*  allhisto=xNew<IssmDouble>(size*nbins);
+				
+				ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
+
+				/*add to results:*/
+				if(my_rank==0){
+					char fieldname[1000];
+
+					sprintf(fieldname,"%s%s",fields[f],"Histogram");
+					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,size,nbins,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* histo=timehistogram[f];
+			IssmDouble* allhisto=xNewZeroInit<IssmDouble>(nbins);
+
+			ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
+
+			/*add to results at time step 1:*/
+			if(my_rank==0){
+				char fieldname[1000];
+
+				sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram");
+				results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,1,nbins,1,0));
+			}
+		} /*}}}*/
+		else{ /*deal with arrays:{{{*/
+
+			int size=meanxsize[f];
+
+			/*we are broadcasting double arrays:*/
+			IssmDouble* histo=timehistogram[f];
+			IssmDouble* allhisto=xNewZeroInit<IssmDouble>(size*nbins);
+
+			ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
+			/*add to results at step 1:*/
+			if(my_rank==0){
+				char fieldname[1000];
+
+				sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram");
+				results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,size,nbins,1,0));
+			}
+		} /*}}}*/
+	}
 } /*}}}*/
 int OutputStatistics(Parameters* parameters,Results* results){ /*{{{*/