Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

QmuStatisticsx.cpp

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

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

Context Navigation

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

Download in other formats: