Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

QmuStatisticsx.cpp

Last change on this file was 27275, checked in by Eric.Larour, 3 years ago
CHG: indices were not being picked up!
File size: 45.8 KB

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

Note: See TracBrowser for help on using the repository browser.

Context Navigation

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

Download in other formats: