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

Last change on this file since 27273 was 27273, checked in by Eric.Larour, 3 years ago

CHG: fixing issues with time indexing of results. fixed lots of memory leaks.

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