Index: /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp
===================================================================
--- /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp	(revision 25479)
+++ /issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp	(revision 25480)
@@ -27,4 +27,5 @@
 	int    nindices;
 	int*   indices=NULL;
+	int    nbins;
 	
 	/*datasets*/
@@ -92,4 +93,13 @@
 
 		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,9 +254,19 @@
 	_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,5 +274,4 @@
 				IssmDouble sumscalar;
 				IssmDouble sumscalar2;
-				IssmDouble mean,stddev;
 
 				ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
@@ -266,5 +285,5 @@
 				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,16 +291,42 @@
 					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,8 +348,30 @@
 					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,4 +529,284 @@
 
 } /*}}}*/
+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,6 +816,8 @@
 	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,5 +829,7 @@
 	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));