/*!\file:  QmuStatisticsx routines
 */ 
/*includes and prototypes:*/
#include <sys/stat.h>
#include "./QmuStatisticsx.h"
#include "../OutputResultsx/OutputResultsx.h"

int DakotaStatistics(int argc,char** argv){ /*{{{*/

	/*Variables:{{{*/
	char* input_file; 
	FILE* fid;
	IoModel* iomodel=NULL;
	ISSM_MPI_Comm statcomm;
	int my_rank;

	//qmu statistics
	bool statistics    = false;
	int  numstatistics = 0;
	int  numdirectories = 0;
	int  nfilesperdirectory = 0;
	char string[1000];
	char* name = NULL;
	char** fields = NULL;
	int    nfields; 
	int*   steps=NULL;
	int    nsteps;
	int    nbins;
	int*   indices=NULL;
	int    nindices;
	int    nsamples;
	int    dummy;
	char*  directory=NULL;
	char*  model=NULL;
	Results* results=NULL;
	Parameters* parameters=NULL;
	int color;
	/*}}}*/
	//First things first, set the communicator as a global variable and be sure we are all here: {{{
	IssmComm::SetComm(MPI_COMM_WORLD);
	my_rank=IssmComm::GetRank();

	/*Barrier:*/
	ISSM_MPI_Barrier(IssmComm::GetComm());
	_printf0_("Dakota Statistic Computation" << "\n");
	/*}}}*/
	//Open model input file for reading {{{
	input_file=xNew<char>((strlen(argv[2])+strlen(argv[3])+strlen(".bin")+2));
	sprintf(input_file,"%s/%s%s",argv[2],argv[3],".bin");
	fid=fopen(input_file,"rb");
	if (fid==NULL) Cerr << "issm_dakota_statistics error message: could not open model " << input_file << " to retrieve qmu statistics parameters" << std::endl;
	//}}}
	//Initialize IoModel, but light version, we'll need it to fetch constants: {{{
	iomodel=new IoModel();
	iomodel->fid=fid;
	iomodel->FetchConstants();
	/*}}}*/
	//Early return if statistics not requested:  {{{
	iomodel->FindConstant(&statistics,"md.qmu.statistics");
	if(!statistics){
		delete iomodel;
		xDelete<char>(input_file);
		fclose(fid); 
		return 0;
	}
	/*}}}*/
	//Create parameters datasets with al the qmu statistics settings we need:  {{{

	/*Initialize parameters and results:*/
	results   = new Results();
	parameters=new Parameters();

	//solution type: 
	parameters->AddObject(new IntParam(SolutionTypeEnum,StatisticsSolutionEnum));

	//root  directory
	directory=xNew<char>(strlen(argv[2])+1);
	xMemCpy<char>(directory,argv[2],strlen(argv[2])+1);
	parameters->AddObject(new StringParam(DirectoryNameEnum,directory));

	//model  name
	model=xNew<char>(strlen(argv[3])+1);
	xMemCpy<char>(model,argv[3],strlen(argv[3])+1);
	parameters->AddObject(new StringParam(InputFileNameEnum,model));

	//nsamples
	iomodel->FindConstant(&nsamples,"md.qmu.method.params.samples");
	parameters->AddObject(new IntParam(QmuNsampleEnum,nsamples));

	//ndirectories
	iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories");
	parameters->AddObject(new IntParam(QmuNdirectoriesEnum,numdirectories));

	//nfiles per directory
	iomodel->FindConstant(&nfilesperdirectory,"md.qmu.statistics.nfiles_per_directory");
	parameters->AddObject(new IntParam(QmuNfilesPerDirectoryEnum,nfilesperdirectory));
	/*}}}*/
	/*Create MPI world: {{{*/
	//At this point, we don't want to go forward any longer, we want to create an MPI 
	//communicator on which to carry out the computations:
	if ((my_rank+1)*nfilesperdirectory>nsamples)color=MPI_UNDEFINED;
	else color=0;
	ISSM_MPI_Comm_split(ISSM_MPI_COMM_WORLD,color, my_rank, &statcomm);
	/*}}}*/

	iomodel->FindConstant(&numstatistics,"md.qmu.statistics.numstatistics");
	for (int i=1;i<=numstatistics;i++){
	//for (int i=9;i<=9;i++){

		char* directory=NULL;
		char* model=NULL;
		int   nsamples;
		_printf0_("Dealing with qmu statistical computation #" << i << "\n");

		sprintf(string,"md.qmu.statistics.method(%i).name",i);
		iomodel->FindConstant(&name,string);

		sprintf(string,"md.qmu.statistics.method(%i).fields",i);
		iomodel->FindConstant(&fields,&nfields,string);
		parameters->AddObject(new StringArrayParam(FieldsEnum,fields,nfields));

		sprintf(string,"md.qmu.statistics.method(%i).steps",i);
		iomodel->FetchData(&steps,&dummy,&nsteps,string);
		parameters->AddObject(new IntVecParam(StepsEnum,steps,nsteps));

		if (strcmp(name,"Histogram")==0){
			/*fetch nbins: */
			sprintf(string,"md.qmu.statistics.method(%i).nbins",i);
			iomodel->FindConstant(&nbins,string);
			parameters->AddObject(new IntParam(NbinsEnum,nbins));
			ComputeHistogram(parameters,results,color,statcomm);
		}
		else if (strcmp(name,"SampleSeries")==0){
			/*fetch indices: */
			sprintf(string,"md.qmu.statistics.method(%i).indices",i);
			iomodel->FetchData(&indices,&nindices,&dummy,string);
			parameters->AddObject(new IntVecParam(IndicesEnum,indices,nindices));

			ComputeSampleSeries(parameters,results,color,statcomm);
		}
		else if (strcmp(name,"MeanVariance")==0){
			ComputeMeanVariance(parameters,results,color,statcomm);
		}
		else _error_(" error creating qmu statistics methods parameters: unsupported method " << name);
	}

	/*Delete resources:*/
	xDelete<char>(input_file);
	delete iomodel;

	//close model file: 
	fclose(fid);

	/*output results:*/
	OutputStatistics(parameters,results,color,statcomm);

	/*all meet here: */
	ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); _printf0_("Output file.\n");

	/*Delete resources:*/
	delete parameters; 
	delete results;

	return 1;
} /*}}}*/
int ComputeHistogram(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){  /*{{{*/

	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;
	int nfilesperdirectory;

	/*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;

	/*only work on the statistical communicator: */
	if (color==MPI_UNDEFINED)return 0;

	/*Retrieve parameters:*/
	parameters->FindParam(&nfilesperdirectory,QmuNfilesPerDirectoryEnum);
	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 from the stat comm communicator:*/
	IssmComm::SetComm(statcomm);
	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<int>(nfields);
	meanxsize=xNew<int>(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");
		/*First read file to figure out size of it in order to create memory buffer mapping into the file.  {{{
		 *This makes it much more efficient to read files without lag.:*/
		char file[1000];
		long int  length;
		char* buffer=NULL;

		/*string:*/
		sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);

		/*open file: */
		_printf0_("    opening file: " << 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 which will be use to read the files:*/
		_printf0_("    processing file:\n");
		fid=fmemopen(buffer, length, "rb");
		/*}}}*/
		/*Figure out for each field, each time step, arrays on each cpu holwing min anx max values:{{{*/
		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){
						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]=xNew<IssmDouble>(doublematsize); 
						xMemCpy<IssmDouble>(maxxs[counter],doublemat,doublematsize);
						minxs[counter]=xNew<IssmDouble>(doublematsize); 
						xMemCpy<IssmDouble>(minxs[counter],doublemat,doublematsize);
						xsize[counter]=doublematsize;
						xDelete<IssmDouble>(doublemat);
					}
					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=maxxs[counter];
						IssmDouble* newmin=minxs[counter];
						for(int k=0;k<doublematsize;k++){
							if(doublemat[k]>newmax[k])newmax[k]=doublemat[k];
							if(doublemat[k]<newmin[k])newmin[k]=doublemat[k];
						}
						xDelete<IssmDouble>(doublemat);
					}
					else _error_("cannot carry out statistics on type " << xtype[counter]); 
				}
			}
		}
		/*}}}*/
		/*Same thing for average in time:{{{*/
		_printf0_("    average in time:\n");

		/*Deal with average in time: */
		for (int f=0;f<nfields;f++){
			fseek(fid,0,SEEK_SET);
			char* field=fields[f];
			meanxtype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);

			if(meanxtype[f]==1){
				meanxsize[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)){
					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{
				meanxsize[f]=doublematsize;
				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;k<doublematsize;k++){
						timemean[k]+=doublemat[k]/nsteps;
					}
					xDelete<IssmDouble>(doublemat);
				}

				if(i==(lower_row+1)){
					maxmeans[f]=xNew<IssmDouble>(doublematsize);
					xMemCpy<IssmDouble>(maxmeans[f],timemean,doublematsize);
					minmeans[f]=xNew<IssmDouble>(doublematsize);
					xMemCpy<IssmDouble>(minmeans[f],timemean,doublematsize);
				}
				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;
				}
			}
		}
		/*}}}*/
		/*Done reading files, close buffer and free memory:{{{*/
		fclose(fid);
		xDelete<char>(buffer);
		/*}}}*/
	}
	ISSM_MPI_Barrier(IssmComm::GetComm());
	_printf0_("Done reading files, now computing min and max.\n"); 

	/*We have collected 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_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());

				/*Store broadcasted value for later computation of histograms:*/
				*maxxs[counter]=allmaxscalar;
				*minxs[counter]=allminscalar;

			}
		} /*}}}*/
		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_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());

				/*Store broadcasted value for later computation of histograms:*/
				maxxs[counter]=allmax;
				minxs[counter]=allmin;
			}
		} /*}}}*/
	}

	/*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_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());

			/*Store for later use in histogram computation:*/
			*maxmeans[f]=allmaxscalar;
			*minmeans[f]=allminscalar;

		} /*}}}*/
		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_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());

			/*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");
		/*read file and make a buffer:{{{*/
		char file[1000];
		long 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");
		/*}}}*/
		/*read data and fill up the histogram using the min and max values from before:{{{*/
		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--;
						if(ma==mi)index=0;
						//_printf_( index << "|" << scalar << "|" << mi << "|" << ma << "|" << nbins << "\n");
						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--;
							if (mi[k]==ma[k])index=0;
							_assert_(scalar<=ma[k]); _assert_(scalar>=mi[k]); _assert_(index<nbins);
							localhistogram[k*nbins+index]++;
						}
						histogram[counter]=localhistogram;
						xDelete<IssmDouble>(doublemat);
					}
					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;
						if(ma==mi)index=0;
						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;
							if (mi[k]==ma[k])index=0;

							localhistogram[k*nbins+index]++;
						}
						xDelete<IssmDouble>(doublemat);
					}
					else _error_("cannot carry out statistics on type " << xtype[counter]); 
				}
			}
		}
		/*}}}*/
		/*Deal with average in time: {{{*/
		_printf0_("    average in time:\n");
		for (int f=0;f<nfields;f++){
			fseek(fid,0,SEEK_SET);
			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;
					if(ma==mi)index=0;
					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;
					if(ma==mi)index=0;
					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;k<doublematsize;k++){
						timemean[k]+=doublemat[k]/nsteps;
					}
					xDelete<IssmDouble>(doublemat);
				}

				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;
						if (mi[k]==ma[k])index=0;
						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;
						if (mi[k]==ma[k])index=0;

						localhistogram[k*nbins+index]++;
					}
				}
			}
		}
		/*}}}*/
		/*close file and delete allocation:{{{*/
		fclose(fid);
		xDelete<char>(buffer);
		/*}}}*/
	}


	/*We have agregated histograms across the cluster, now gather them across  the cluster onto
	 *cpu0: */
	_printf0_("Collect histograms on cpu 0 and save to results:\n");
	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());
				xDelete<IssmDouble>(histo);

				/*add to results while deallocating as much as possible:*/
				char fieldname[1000];
				if(my_rank==0){
					sprintf(fieldname,"%s%s",fields[f],"Histogram"); results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,1,nbins,j+1,0));
				}
				xDelete<IssmDouble>(allhisto); 
				if(my_rank==0){
					sprintf(fieldname,"%s%s",fields[f],"Max"); results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*maxxs[counter],steps[j],0));
					sprintf(fieldname,"%s%s",fields[f],"Min"); results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*minxs[counter],steps[j],0));
				}
				xDelete<IssmDouble>(maxxs[counter]); 
				xDelete<IssmDouble>(minxs[counter]); 
			}
		} /*}}}*/
		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_SUM,IssmComm::GetComm());
				xDelete<IssmDouble>(histo);

				/*add to results:*/
				char fieldname[1000];
				if(my_rank==0){
					sprintf(fieldname,"%s%s",fields[f],"Histogram"); results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,size,nbins,j+1,0));
				}
				xDelete<IssmDouble>(allhisto);
				if(my_rank==0){
					sprintf(fieldname,"%s%s",fields[f],"Max"); results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,maxxs[counter],size,1,steps[j],0));
					sprintf(fieldname,"%s%s",fields[f],"Min"); results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,minxs[counter],size,1,steps[j],0));
				}
				xDelete<IssmDouble>(maxxs[counter]);
				xDelete<IssmDouble>(minxs[counter]);
			}
		} /*}}}*/
	}
	/*Now do the same for the time mean fields:*/
	_printf0_("Collect time mean histograms  on cpu 0 and save to results:\n");
	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());
			xDelete<IssmDouble>(histo);

			/*add to results at time step 1:*/
			char fieldname[1000];
			if(my_rank==0){
				sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram"); results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,1,nbins,steps[0],0));
			}
			xDelete<IssmDouble>(allhisto);
			if(my_rank==0){
				sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax"); results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*maxmeans[f],steps[0],0));
				sprintf(fieldname,"%s%s",fields[f],"TimeMeaMin"); results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*minmeans[f],steps[0],0));
			}
			xDelete<IssmDouble>(maxmeans[f]);
			xDelete<IssmDouble>(minmeans[f]);
		} /*}}}*/
		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_SUM,IssmComm::GetComm());
			xDelete<IssmDouble>(histo);
			
			/*add to results at step 1:*/
			char fieldname[1000];
			if(my_rank==0){
				sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram"); results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,size,nbins,steps[0],0));
			}
			xDelete<IssmDouble>(allhisto);
			if(my_rank==0){
				sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax"); results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,maxmeans[f],size,1,steps[0],0));
				sprintf(fieldname,"%s%s",fields[f],"TimeMeanMin"); results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,minmeans[f],size,1,steps[0],0));
			}
			xDelete<IssmDouble>(maxmeans[f]);
			xDelete<IssmDouble>(minmeans[f]);
		} /*}}}*/
	}
	_printf0_("Done aggregating time mean histogram:\n");
	IssmComm::SetComm(ISSM_MPI_COMM_WORLD);
	
	/*Free allocations:*/
	xDelete<char>(directory);
	xDelete<char>(model);
	for (int i=0;i<nfields;i++)xDelete<char>(fields[i]);
	xDelete<char*>(fields);
	xDelete<int>(steps);
	xDelete<IssmDouble*>(maxxs); 
	xDelete<IssmDouble*>(minxs);
	xDelete<IssmDouble*>(maxmeans); 
	xDelete<IssmDouble*>(minmeans);
	xDelete<int>(xtype);
	xDelete<int>(xsize);

	return 1;
}
/*}}}*/
int ComputeMeanVariance(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){  /*{{{*/

	/*variables: {{{*/
	int nsamples; 
	char* directory=NULL;
	char* model=NULL;
	char** fields=NULL;
	int* steps=NULL;
	int nsteps;
	int nfields;
	int range,lower_row,upper_row;
	int nfilesperdirectory;

	/*intermediary:*/
	IssmDouble* doublemat=NULL;
	int         doublematsize;
	IssmDouble scalar;

	/*computation of average and variance itself:*/
	IssmDouble*  x = NULL;
	IssmDouble*  x2 = NULL;
	IssmDouble** xs = NULL;
	IssmDouble** xs2 = NULL;
	int*         xtype=NULL;
	int*         xsize=NULL;

	IssmDouble** meanx=NULL;
	IssmDouble** meanx2=NULL;
	int*         meantype=NULL;
	int*         meansize=NULL;
	/*}}}*/

	/*only work on the statistical communicator: */
	if (color==MPI_UNDEFINED)return 0;

	/*Retrieve parameters:*/
	parameters->FindParam(&nfilesperdirectory,QmuNfilesPerDirectoryEnum);
	parameters->FindParam(&nsamples,QmuNsampleEnum);
	parameters->FindParam(&directory,DirectoryNameEnum);
	parameters->FindParam(&model,InputFileNameEnum);
	parameters->FindParam(&fields,&nfields,FieldsEnum);
	parameters->FindParam(&steps,&nsteps,StepsEnum);

	/*Get rank from the stat comm communicator:*/
	IssmComm::SetComm(statcomm);
	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:{{{*/
	xs=xNew<IssmDouble*>(nfields*nsteps);
	xs2=xNew<IssmDouble*>(nfields*nsteps);
	xtype=xNew<int>(nfields*nsteps);
	xsize=xNew<int>(nfields*nsteps);

	meantype=xNew<int>(nfields);
	meansize=xNew<int>(nfields);
	meanx=xNew<IssmDouble*>(nfields);
	meanx2=xNew<IssmDouble*>(nfields);
	/*}}}*/

	/*Start opening files:*/
	for (int i=(lower_row+1);i<=upper_row;i++){
		_printf0_("reading file #: " << i << "\n");
		/*open buffer linked to file: {{{*/
		char file[1000];
		long int  length;
		char* buffer=NULL;

		/*string:*/
		sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);

		/*open file: */
		_printf0_("    opening file: " << file << "\n");
		FILE* fid=fopen(file,"rb");
		if(fid==NULL) _error_("    could 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");
		/*}}}*/
		/*read x and x^2 values for each time step:{{{*/
		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){
						xs[counter]=xNew<IssmDouble>(1); 
						xs2[counter]=xNew<IssmDouble>(1); 
						*xs[counter]=scalar;
						*xs2[counter]=pow(scalar,2.0);
						xsize[counter]=1;
					}
					else if (xtype[counter]==3){
						IssmDouble* doublemat2=xNew<IssmDouble>(doublematsize);
						for(int k=0;k<doublematsize;k++)doublemat2[k]=pow(doublemat[k],2.0);
						xs[counter]=doublemat;
						xs2[counter]=doublemat2;
						xsize[counter]=doublematsize;
					}
					else _error_("cannot carry out statistics on type " << xtype[counter]); 
				}
				else{
					if(xtype[counter]==1){
						*xs[counter]+=scalar;
						*xs2[counter]+=pow(scalar,2.0);
					}
					else if (xtype[counter]==3){
						IssmDouble* newdoublemat=xs[counter];
						IssmDouble* newdoublemat2=xs2[counter];
						for(int k=0;k<doublematsize;k++){
							newdoublemat[k]+=doublemat[k];
							newdoublemat2[k]+=pow(doublemat[k],2.0);
						}
						xs[counter]=newdoublemat;
						xs2[counter]=newdoublemat2;
					}
					else _error_("cannot carry out statistics on type " << xtype[counter]); 
				}
			}
		}/*}}}*/
		/*Same but for time mean: {{{*/
		for (int f=0;f<nfields;f++){
			char* field=fields[f];
			fseek(fid,0,SEEK_SET);
			meantype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
			if(i==(lower_row+1)){
				if(meantype[f]==1){
					meanx[f]=xNewZeroInit<IssmDouble>(1);
					meanx2[f]=xNewZeroInit<IssmDouble>(1);
					meansize[f]=1;
				}
				else{
					meanx[f]=xNewZeroInit<IssmDouble>(doublematsize);
					meanx2[f]=xNewZeroInit<IssmDouble>(doublematsize);
					meansize[f]=doublematsize;
					xDelete<IssmDouble>(doublemat);
				}
			}
			fseek(fid,0,SEEK_SET); 
			if(meantype[f]==1){
				IssmDouble sc=0;
				IssmDouble sc2=0;
				for(int j=0;j<nsteps;j++){
					readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
					sc+=scalar/nsteps;
				}
				sc2+=pow(sc,2.0);
				*meanx[f]+=sc;
				*meanx2[f]+=sc2;
			}
			else{
				IssmDouble* sc=meanx[f];
				IssmDouble* sc2=meanx2[f];
				IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
				IssmDouble* timemean2=xNewZeroInit<IssmDouble>(doublematsize);

				for(int j=0;j<nsteps;j++){
					readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
					for (int k=0;k<doublematsize;k++){
						timemean[k]+=doublemat[k]/nsteps;
					}
					xDelete<IssmDouble>(doublemat);
				}
				for (int k=0;k<doublematsize;k++){
					timemean2[k]=pow(timemean[k],2.0);
				}
				for (int k=0;k<doublematsize;k++){
					sc[k]+=timemean[k];
					sc2[k]+=timemean2[k];
				}

			}

		} /*}}}*/
		/*cleanup:{{{*/
		fclose(fid);
		xDelete<char>(buffer);
		/*}}}*/
	}
	ISSM_MPI_Barrier(IssmComm::GetComm());
	_printf0_("Done reading files, now computing mean and variance.\n"); 

	/*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++){
		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];
				IssmDouble scalar2=*xs2[counter];
				IssmDouble sumscalar;
				IssmDouble sumscalar2;

				ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
				ISSM_MPI_Reduce(&scalar2,&sumscalar2,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
				
				xDelete<IssmDouble>(xs[counter]);
				xDelete<IssmDouble>(xs2[counter]);

				/*Build average and standard deviation. For standard deviation, use the 
				 *following formula: sigma^2=E(x^2)-mu^2:*/
				mean=sumscalar/nsamples;
				stddev=sqrt(sumscalar2/nsamples-pow(mean,2.0));

				/*add to results:*/
				if(my_rank==0){
					char fieldname[1000];

					sprintf(fieldname,"%s%s",fields[f],"Mean");
					results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,steps[j],0));
					sprintf(fieldname,"%s%s",fields[f],"Stddev");
					results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,steps[j],0));
				}

			}
		} /*}}}*/
		else{ /*deal with arrays:{{{*/

			int size=xsize[counter0];

			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];

				IssmDouble*  sumx=xNew<IssmDouble>(size);
				IssmDouble*  sumx2=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());
				
				xDelete<IssmDouble>(xs[counter]);
				xDelete<IssmDouble>(xs2[counter]);

				/*Build average and standard deviation. For standard deviation, use the 
				 *following formula: sigma^2=E(x^2)-mu^2:*/
				for (int k=0;k<size;k++){
					mean[k]=sumx[k]/nsamples;
					stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
				}
				xDelete<IssmDouble>(sumx);
				xDelete<IssmDouble>(sumx2);

				/*add to results:*/
				if(my_rank==0){
					char fieldname[1000];

					sprintf(fieldname,"%s%s",fields[f],"Mean");
					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,steps[j],0));
					sprintf(fieldname,"%s%s",fields[f],"Stddev");
					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,steps[j],0));
				}
			}
			/*Deallocate:*/
			xDelete<IssmDouble>(mean);
			xDelete<IssmDouble>(stddev);

		} /*}}}*/
	}
	/*Do the same but for the time mean:*/
	for (int f=0;f<nfields;f++){
		if (meantype[f]==1){ /*deal with scalars {{{*/
			IssmDouble mean,stddev;

			/*we are broadcasting doubles:*/
			IssmDouble scalar=*meanx[f];
			IssmDouble scalar2=*meanx2[f];
			IssmDouble sumscalar;
			IssmDouble sumscalar2;

			xDelete<IssmDouble>(meanx[f]);
			xDelete<IssmDouble>(meanx2[f]);

			ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
			ISSM_MPI_Reduce(&scalar2,&sumscalar2,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
			/*Build average and standard deviation. For standard deviation, use the 
			 *following formula: sigma^2=E(x^2)-mu^2:*/
			mean=sumscalar/nsamples;
			stddev=sqrt(sumscalar2/nsamples-pow(mean,2.0));

			/*add to results:*/
			if(my_rank==0){
				char fieldname[1000];

				sprintf(fieldname,"%s%s",fields[f],"TimeMean");
				results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,steps[0],0));
				sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
				results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,steps[0],0));
			}
		} /*}}}*/
		else{ /*deal with arrays:{{{*/

			int size=meansize[f];
			IssmDouble*  mean=xNew<IssmDouble>(size);
			IssmDouble*  stddev=xNew<IssmDouble>(size);

			/*we are broadcasting double arrays:*/
			x=meanx[f];
			x2=meanx2[f];

			IssmDouble*  sumx=xNew<IssmDouble>(size);
			IssmDouble*  sumx2=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());

			xDelete<IssmDouble>(meanx[f]);
			xDelete<IssmDouble>(meanx2[f]);

			/*Build average and standard deviation. For standard deviation, use the 
			 *following formula: sigma^2=E(x^2)-mu^2:*/
			for (int k=0;k<size;k++){
				mean[k]=sumx[k]/nsamples;
				stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
			}

			/*add to results:*/
			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,steps[0],0));
				sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
				results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,steps[0],0));
			}
			/*Deallocate:*/
			xDelete<IssmDouble>(sumx);
			xDelete<IssmDouble>(sumx2);
			xDelete<IssmDouble>(mean);
			xDelete<IssmDouble>(stddev);
		} /*}}}*/
	}
	_printf0_("Done with MeanVariance:\n");
	IssmComm::SetComm(ISSM_MPI_COMM_WORLD);
	return 1;
} /*}}}*/
int ComputeSampleSeries(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){ /*{{{*/

	int nsamples; 
	char* directory=NULL;
	char* model=NULL;
	char** fields=NULL;
	int* steps=NULL;
	int nsteps;
	int nfields;
	int range,lower_row,upper_row;
	int nfilesperdirectory;
	int* indices=NULL;
	int  nindices;

	/*intermediary:*/
	IssmDouble* doublemat=NULL;
	int         doublematsize;
	IssmDouble scalar;

	/*computation of average and variance itself:*/
	IssmDouble*  x = NULL;
	IssmDouble*  allx=NULL;
	IssmDouble** xs = NULL;
	int*         xtype=NULL;
	int*         xsize=NULL;

	/*only work on the statistical communicator: */
	if (color==MPI_UNDEFINED)return 0;

	/*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(&indices,&nindices,IndicesEnum);

	/*Get rank from the stat comm communicator:*/
	IssmComm::SetComm(statcomm);
	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:*/
	xs=xNew<IssmDouble*>(nfields*nsteps);
	xtype=xNew<int>(nfields*nsteps);
	xsize=xNew<int>(nfields*nsteps);

	/*_printf0_("nindices: " << nindices << "\n");
	for (int i=0;i<nindices;i++){
		_printf0_(indices[i] << " " );
	}
	_printf0_("\n");

	_printf0_("nsteps: " << nsteps << "\n");
	for (int i=0;i<nsteps;i++){
		_printf0_(steps[i] << " " );
	}
	_printf0_("\n");*/



	/*Start opening files:*/
	for (int i=(lower_row+1);i<=upper_row;i++){
		_printf0_("reading file #: " << i << "\n");
		/*Create memory buffer for file, to speed things up: {{{*/
		char file[1000];
		long 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");
		/*}}}*/
		/*Retrieve the values for all fields and time steps:{{{*/
		for (int f=0;f<nfields;f++){
			fseek(fid,0,SEEK_SET);
			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){
						x=xNew<IssmDouble>(range);
						x[0]=scalar;
						xs[counter]=x;
						xsize[counter]=range;
					}
					else if (xtype[counter]==3){
						x=xNew<IssmDouble>(nindices*range);
						for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
						xs[counter]=x;
						xsize[counter]=range*nindices;
					}
					else _error_("cannot carry out statistics on type " << xtype[counter]); 
				}
				else{
					if(xtype[counter]==1){
						x=xs[counter]; 
						x[i-(lower_row+1)]=scalar;
						xs[counter]=x;
					}
					else if (xtype[counter]==3){
						x=xs[counter];
						for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
						xs[counter]=x;
					}
					else _error_("cannot carry out statistics on type " << xtype[counter]); 
				}
			}
		}
		fclose(fid);
		xDelete<char>(buffer);
		/*}}}*/
	}
	ISSM_MPI_Barrier(IssmComm::GetComm());
	_printf0_("Done reading files, now assembling time series.\n"); //{{{

	for (int f=0;f<nfields;f++){
		for (int j=0;j<nsteps;j++){
			int counter=f*nsteps+j;
			if (xtype[counter]==1){
				/*we are broadcasting range times doubles:*/
				x=xs[counter]; 
				allx=xNew<IssmDouble>(nsamples);
				MPI_Gather(x, range, ISSM_MPI_PDOUBLE,allx, range, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
				/*add to results:*/
				if(my_rank==0){
					char fieldname[1000];

					sprintf(fieldname,"%s%s",fields[f],"Samples");
					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,1,j+1,0));
				}
			}
			else{
				/*we are broadcasting double arrays:*/
				x=xs[counter];
				allx=xNew<IssmDouble>(nsamples*nindices);

				MPI_Gather(x, range*nindices, ISSM_MPI_PDOUBLE,allx, range*nindices, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());

				/*add to results:*/
				if(my_rank==0){
					char fieldname[1000];
					sprintf(fieldname,"%s%s",fields[f],"Samples");
					results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,nindices,j+1,0));
				}
			}
		}
	} //}}}
	_printf0_("Done with SampleSeries:\n");
	IssmComm::SetComm(ISSM_MPI_COMM_WORLD);

	return 1;
} /*}}}*/
int OutputStatistics(Parameters* parameters,Results* results,int color,ISSM_MPI_Comm statcomm){ /*{{{*/

	char   outputfilename[1000];
	char* directory=NULL;
	char* model=NULL;
	char* method=NULL;
	int   nsamples;
	int* steps=NULL;
	int nsteps;

	/*only work on the statistical communicator: */
	if (color==MPI_UNDEFINED)return 0;

	FemModel* femmodel=new FemModel();

	/*Some parameters that will allow us to use the OutputResultsx module:*/
	parameters->AddObject(new BoolParam(QmuIsdakotaEnum,false));
	parameters->AddObject(new BoolParam(SettingsIoGatherEnum,true));

	parameters->FindParam(&directory,DirectoryNameEnum);
	parameters->FindParam(&model,InputFileNameEnum);
	parameters->FindParam(&nsamples,QmuNsampleEnum);
	parameters->FindParam(&steps,&nsteps,StepsEnum);

	sprintf(outputfilename,"%s/%s.stats",directory,model);
	parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));

	/*Call OutputResults module:*/
	femmodel->parameters=parameters;
	femmodel->results=results;

	OutputResultsx(femmodel);

	return 1;
} /*}}}*/
int readdata(IssmDouble** pdoublemat, int* pdoublematsize, IssmDouble* pdouble, FILE* fid,char* field,int step){ /*{{{*/

	int length;
	char fieldname[1000];
	int   fieldname_size;
	IssmDouble   rtime;
	int          rstep;
	int M,N;

	//fields that we retrive: 
	IssmDouble  dfield; 
	char*       sfield    = NULL;
	IssmDouble* dmatfield = NULL; 
	int*        imatfield = NULL; 

	//type of the returned field: 
	int type;
	int found=0;

	while(1){

		size_t ret_code = fread(&fieldname_size, sizeof(int), 1, fid); 
		if(ret_code != 1) break; //we are done.

		fread(fieldname, sizeof(char), fieldname_size, fid); 
		//_printf0_("fieldname: " << fieldname << "\n");

		fread(&rtime, sizeof(IssmDouble), 1, fid); 
		fread(&rstep, sizeof(int), 1, fid); 

		//check on field: 
		if ((step==rstep) && (strcmp(field,fieldname)==0)){

			//ok, go read the result really: 
			fread(&type,sizeof(int),1,fid);
			fread(&M,sizeof(int),1,fid);
			if (type==1){
				fread(&dfield,sizeof(IssmDouble),1,fid);
			}
			else if (type==2){
				fread(&M,sizeof(int),1,fid);
				sfield=xNew<char>(M);
				fread(sfield,sizeof(char),M,fid);
			}
			else if (type==3){
				fread(&N,sizeof(int),1,fid);
				dmatfield=xNew<IssmDouble>(M*N);
				fread(dmatfield,sizeof(IssmDouble),M*N,fid);
			}
			else if (type==4){
				fread(&N,sizeof(int),1,fid);
				imatfield=xNew<int>(M*N);
				fread(imatfield,sizeof(int),M*N,fid);
			}
			else _error_("cannot read data of type " << type << "\n");
			found=1;
			break;
		}
		else{
			//just skim to next results.
			fread(&type,sizeof(int),1,fid);
			fread(&M,sizeof(int),1,fid);
			if (type==1){
				fseek(fid,sizeof(IssmDouble),SEEK_CUR);
			}
			else if(type==2){
				fseek(fid,M*sizeof(char),SEEK_CUR);
			}
			else if(type==3){
				fread(&N,sizeof(int),1,fid);
				fseek(fid,M*N*sizeof(IssmDouble),SEEK_CUR);
			}
			else if(type==4){
				fread(&N,sizeof(int),1,fid);
				fseek(fid,M*N*sizeof(int),SEEK_CUR);
			}
			else _error_("cannot read data of type " << type << "\n");
		}
	}
	if(found==0)_error_("cound not find " << field << " at step " << step  << "\n");

	/*assign output pointers:*/
	*pdoublemat=dmatfield;
	*pdoublematsize=M*N;
	*pdouble=dfield;

	/*return:*/
	return type;

}
/*}}}*/
bool DakotaDirStructure(int argc,char** argv){ /*{{{*/

	char* input_file; 
	FILE* fid;
	IoModel* iomodel=NULL;
	int check;

	//qmu statistics
	bool statistics    = false;
	int  numdirectories = 0;

	/*First things first, set the communicator as a global variable: */
	IssmComm::SetComm(MPI_COMM_WORLD);

	/*Barrier:*/
	ISSM_MPI_Barrier(IssmComm::GetComm());
	_printf0_("Preparing directory structure for model outputs:" << "\n");

	//open model input file for reading
	input_file=xNew<char>((strlen(argv[2])+strlen(argv[3])+strlen(".bin")+2));
	sprintf(input_file,"%s/%s%s",argv[2],argv[3],".bin");
	fid=fopen(input_file,"rb");
	if (fid==NULL) Cerr << "dirstructure error message: could not open model " << input_file << " to retrieve qmu statistics parameters" << std::endl;

	//initialize IoModel, but light version, we just need it to fetch one constant: 
	iomodel=new IoModel();
	iomodel->fid=fid;
	iomodel->FetchConstants();

	//early return if statistics not requested: 
	iomodel->FindConstant(&statistics,"md.qmu.statistics");
	if(!statistics){
		delete iomodel;
		xDelete<char>(input_file);
		fclose(fid);
		return false; //important return value!
	}

	iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories");

	/*Ok, we have everything we need to create the directory structure:*/
	if(IssmComm::GetRank()==0){
		for (int i=0;i<numdirectories;i++){
			char directory[1000];
			sprintf(directory,"./%i",i+1);

			check = mkdir(directory,ACCESSPERMS);
			if (check) _error_("dirstructure error message: could not create directory " << directory << "\n");
		}
	}

	/*Delete resources:*/
	delete iomodel;
	xDelete<char>(input_file);

	//close model file: 
	fclose(fid);

	//return value: 
	return true; //statistics computation on!
} /*}}}*/