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

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

CHG: better optimization for memory, fixed leaks.

File size: 44.6 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],j+1,0));
660 sprintf(fieldname,"%s%s",fields[f],"Min"); results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*minxs[counter],j+1,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,j+1,0));
687 sprintf(fieldname,"%s%s",fields[f],"Min"); results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,minxs[counter],size,1,j+1,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 int nsamples;
766 char* directory=NULL;
767 char* model=NULL;
768 char** fields=NULL;
769 int* steps=NULL;
770 int nsteps;
771 int nfields;
772 int range,lower_row,upper_row;
773 int nfilesperdirectory;
774
775 /*intermediary:*/
776 IssmDouble* doublemat=NULL;
777 int doublematsize;
778 IssmDouble scalar;
779
780 /*computation of average and variance itself:*/
781 IssmDouble* x = NULL;
782 IssmDouble* x2 = NULL;
783 IssmDouble** xs = NULL;
784 IssmDouble** xs2 = NULL;
785 int* xtype=NULL;
786 int* xsize=NULL;
787
788 IssmDouble** meanx=NULL;
789 IssmDouble** meanx2=NULL;
790 int* meantype=NULL;
791 int* meansize=NULL;
792
793 /*only work on the statistical communicator: */
794 if (color==MPI_UNDEFINED)return 0;
795
796 /*Retrieve parameters:*/
797 parameters->FindParam(&nfilesperdirectory,QmuNfilesPerDirectoryEnum);
798 parameters->FindParam(&nsamples,QmuNsampleEnum);
799 parameters->FindParam(&directory,DirectoryNameEnum);
800 parameters->FindParam(&model,InputFileNameEnum);
801 parameters->FindParam(&fields,&nfields,FieldsEnum);
802 parameters->FindParam(&steps,&nsteps,StepsEnum);
803
804 /*Get rank from the stat comm communicator:*/
805 IssmComm::SetComm(statcomm);
806 int my_rank=IssmComm::GetRank();
807
808 /*Open files and read them complelety, in a distributed way:*/
809 range=DetermineLocalSize(nsamples,IssmComm::GetComm());
810 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
811
812 /*Initialize arrays:*/
813 xs=xNew<IssmDouble*>(nfields*nsteps);
814 xs2=xNew<IssmDouble*>(nfields*nsteps);
815 xtype=xNew<int>(nfields*nsteps);
816 xsize=xNew<int>(nfields*nsteps);
817
818 meantype=xNew<int>(nfields);
819 meansize=xNew<int>(nfields);
820 meanx=xNew<IssmDouble*>(nfields);
821 meanx2=xNew<IssmDouble*>(nfields);
822
823 /*Start opening files:*/
824 for (int i=(lower_row+1);i<=upper_row;i++){
825 _printf0_("reading file #: " << i << "\n");
826 char file[1000];
827 long int length;
828 char* buffer=NULL;
829
830 /*string:*/
831 sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
832
833 /*open file: */
834 _printf0_(" opening file: " << file << "\n");
835 FILE* fid=fopen(file,"rb");
836 if(fid==NULL) _error_(" could not open file: " << file << "\n");
837
838 /*figure out size of file, and read the whole thing:*/
839 _printf0_(" reading file:\n");
840 fseek (fid, 0, SEEK_END);
841 length = ftell (fid);
842 fseek (fid, 0, SEEK_SET);
843 buffer = xNew<char>(length);
844 fread (buffer, sizeof(char), length, fid);
845
846 /*close file:*/
847 fclose (fid);
848
849 /*create a memory stream with this buffer:*/
850 _printf0_(" processing file:\n");
851 fid=fmemopen(buffer, length, "rb");
852
853 /*start reading data from the buffer directly:*/
854 for (int f=0;f<nfields;f++){
855 char* field=fields[f];
856 fseek(fid,0,SEEK_SET);
857 for (int j=0;j<nsteps;j++){
858 int counter=f*nsteps+j;
859 xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
860 if(i==(lower_row+1)){
861 if(xtype[counter]==1){
862 xs[counter]=xNew<IssmDouble>(1);
863 xs2[counter]=xNew<IssmDouble>(1);
864 *xs[counter]=scalar;
865 *xs2[counter]=pow(scalar,2.0);
866 xsize[counter]=1;
867 }
868 else if (xtype[counter]==3){
869 IssmDouble* doublemat2=xNew<IssmDouble>(doublematsize);
870 for(int k=0;k<doublematsize;k++)doublemat2[k]=pow(doublemat[k],2.0);
871 xs[counter]=doublemat;
872 xs2[counter]=doublemat2;
873 xsize[counter]=doublematsize;
874 }
875 else _error_("cannot carry out statistics on type " << xtype[counter]);
876 }
877 else{
878 if(xtype[counter]==1){
879 *xs[counter]+=scalar;
880 *xs2[counter]+=pow(scalar,2.0);
881 }
882 else if (xtype[counter]==3){
883 IssmDouble* newdoublemat=xs[counter];
884 IssmDouble* newdoublemat2=xs2[counter];
885 for(int k=0;k<doublematsize;k++){
886 newdoublemat[k]+=doublemat[k];
887 newdoublemat2[k]+=pow(doublemat[k],2.0);
888 }
889 xs[counter]=newdoublemat;
890 xs2[counter]=newdoublemat2;
891 }
892 else _error_("cannot carry out statistics on type " << xtype[counter]);
893 }
894 }
895 }
896
897 /*Deal with time mean: */
898 for (int f=0;f<nfields;f++){
899 char* field=fields[f];
900 fseek(fid,0,SEEK_SET);
901 meantype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
902 if(i==(lower_row+1)){
903 if(meantype[f]==1){
904 meanx[f]=xNewZeroInit<IssmDouble>(1);
905 meanx2[f]=xNewZeroInit<IssmDouble>(1);
906 meansize[f]=1;
907 }
908 else{
909 meanx[f]=xNewZeroInit<IssmDouble>(doublematsize);
910 meanx2[f]=xNewZeroInit<IssmDouble>(doublematsize);
911 meansize[f]=doublematsize;
912 }
913 }
914 fseek(fid,0,SEEK_SET);
915 if(meantype[f]==1){
916 IssmDouble sc=0;
917 IssmDouble sc2=0;
918 for(int j=0;j<nsteps;j++){
919 readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
920 sc+=scalar/nsteps;
921 }
922 sc2+=pow(sc,2.0);
923 *meanx[f]+=sc;
924 *meanx2[f]+=sc2;
925 }
926 else{
927 IssmDouble* sc=meanx[f];
928 IssmDouble* sc2=meanx2[f];
929 IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
930 IssmDouble* timemean2=xNewZeroInit<IssmDouble>(doublematsize);
931
932 for(int j=0;j<nsteps;j++){
933 readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
934 for (int k=0;k<doublematsize;k++){
935 timemean[k]+=doublemat[k]/nsteps;
936 }
937 }
938 for (int k=0;k<doublematsize;k++){
939 timemean2[k]=pow(timemean[k],2.0);
940 }
941 for (int k=0;k<doublematsize;k++){
942 sc[k]+=timemean[k];
943 sc2[k]+=timemean2[k];
944 }
945
946 }
947
948 }
949 fclose(fid);
950
951 /*delete buffer:*/
952 xDelete<char>(buffer);
953 }
954 ISSM_MPI_Barrier(IssmComm::GetComm());
955 _printf0_("Done reading files, now computing mean and variance.\n");
956
957 /*We have agregated x and x^2 across the cluster, now gather across the cluster onto
958 *cpu0 and then compute statistics:*/
959 for (int f=0;f<nfields;f++){
960 int counter0=f*nsteps+0;
961 if (xtype[counter0]==1){ /*deal with scalars {{{*/
962 IssmDouble mean,stddev;
963 for (int j=0;j<nsteps;j++){
964 int counter=f*nsteps+j;
965
966 /*we are broadcasting doubles:*/
967 IssmDouble scalar=*xs[counter];
968 IssmDouble scalar2=*xs2[counter];
969 IssmDouble sumscalar;
970 IssmDouble sumscalar2;
971
972 ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
973 ISSM_MPI_Reduce(&scalar2,&sumscalar2,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
974 /*Build average and standard deviation. For standard deviation, use the
975 *following formula: sigma^2=E(x^2)-mu^2:*/
976 mean=sumscalar/nsamples;
977 stddev=sqrt(sumscalar2/nsamples-pow(mean,2.0));
978
979 /*add to results:*/
980 if(my_rank==0){
981 char fieldname[1000];
982
983 sprintf(fieldname,"%s%s",fields[f],"Mean");
984 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,j+1,0));
985 sprintf(fieldname,"%s%s",fields[f],"Stddev");
986 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,j+1,0));
987 }
988
989 }
990 } /*}}}*/
991 else{ /*deal with arrays:{{{*/
992
993 int size=xsize[counter0];
994
995 IssmDouble* mean=xNew<IssmDouble>(size);
996 IssmDouble* stddev=xNew<IssmDouble>(size);
997
998 for (int j=0;j<nsteps;j++){
999 int counter=f*nsteps+j;
1000
1001 /*we are broadcasting double arrays:*/
1002 x=xs[counter];
1003 x2=xs2[counter];
1004
1005 IssmDouble* sumx=xNew<IssmDouble>(size);
1006 IssmDouble* sumx2=xNew<IssmDouble>(size);
1007
1008 ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
1009 ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
1010
1011 /*Build average and standard deviation. For standard deviation, use the
1012 *following formula: sigma^2=E(x^2)-mu^2:*/
1013 for (int k=0;k<size;k++){
1014 mean[k]=sumx[k]/nsamples;
1015 stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
1016 }
1017
1018 /*add to results:*/
1019 if(my_rank==0){
1020 char fieldname[1000];
1021
1022 sprintf(fieldname,"%s%s",fields[f],"Mean");
1023 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,j+1,0));
1024 sprintf(fieldname,"%s%s",fields[f],"Stddev");
1025 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,j+1,0));
1026 }
1027 }
1028 } /*}}}*/
1029 }
1030 /*Do the same but for the time mean:*/
1031 for (int f=0;f<nfields;f++){
1032 if (meantype[f]==1){ /*deal with scalars {{{*/
1033 IssmDouble mean,stddev;
1034
1035 /*we are broadcasting doubles:*/
1036 IssmDouble scalar=*meanx[f];
1037 IssmDouble scalar2=*meanx2[f];
1038 IssmDouble sumscalar;
1039 IssmDouble sumscalar2;
1040
1041 ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
1042 ISSM_MPI_Reduce(&scalar2,&sumscalar2,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
1043 /*Build average and standard deviation. For standard deviation, use the
1044 *following formula: sigma^2=E(x^2)-mu^2:*/
1045 mean=sumscalar/nsamples;
1046 stddev=sqrt(sumscalar2/nsamples-pow(mean,2.0));
1047
1048 /*add to results:*/
1049 if(my_rank==0){
1050 char fieldname[1000];
1051
1052 sprintf(fieldname,"%s%s",fields[f],"TimeMean");
1053 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,steps[0],0));
1054 sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
1055 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,steps[0],0));
1056 }
1057 } /*}}}*/
1058 else{ /*deal with arrays:{{{*/
1059
1060 int size=meansize[f];
1061 IssmDouble* mean=xNew<IssmDouble>(size);
1062 IssmDouble* stddev=xNew<IssmDouble>(size);
1063
1064 /*we are broadcasting double arrays:*/
1065 x=meanx[f];
1066 x2=meanx2[f];
1067
1068 IssmDouble* sumx=xNew<IssmDouble>(size);
1069 IssmDouble* sumx2=xNew<IssmDouble>(size);
1070
1071 ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
1072 ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
1073
1074 /*Build average and standard deviation. For standard deviation, use the
1075 *following formula: sigma^2=E(x^2)-mu^2:*/
1076 for (int k=0;k<size;k++){
1077 mean[k]=sumx[k]/nsamples;
1078 stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
1079 }
1080
1081 /*add to results:*/
1082 if(my_rank==0){
1083 char fieldname[1000];
1084
1085 sprintf(fieldname,"%s%s",fields[f],"TimeMean");
1086 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,steps[0],0));
1087 sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
1088 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,steps[0],0));
1089 }
1090 } /*}}}*/
1091 }
1092
1093 _printf0_("Done with MeanVariance:\n");
1094 IssmComm::SetComm(ISSM_MPI_COMM_WORLD);
1095
1096 return 1;
1097} /*}}}*/
1098int ComputeSampleSeries(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){ /*{{{*/
1099
1100 int nsamples;
1101 char* directory=NULL;
1102 char* model=NULL;
1103 char** fields=NULL;
1104 int* steps=NULL;
1105 int nsteps;
1106 int nfields;
1107 int range,lower_row,upper_row;
1108 int nfilesperdirectory;
1109 int* indices=NULL;
1110 int nindices;
1111
1112 /*intermediary:*/
1113 IssmDouble* doublemat=NULL;
1114 int doublematsize;
1115 IssmDouble scalar;
1116
1117 /*computation of average and variance itself:*/
1118 IssmDouble* x = NULL;
1119 IssmDouble* allx=NULL;
1120 IssmDouble** xs = NULL;
1121 int* xtype=NULL;
1122 int* xsize=NULL;
1123
1124 /*only work on the statistical communicator: */
1125 if (color==MPI_UNDEFINED)return 0;
1126
1127 /*Retrieve parameters:*/
1128 parameters->FindParam(&nsamples,QmuNsampleEnum);
1129 parameters->FindParam(&directory,DirectoryNameEnum);
1130 parameters->FindParam(&model,InputFileNameEnum);
1131 parameters->FindParam(&fields,&nfields,FieldsEnum);
1132 parameters->FindParam(&steps,&nsteps,StepsEnum);
1133 parameters->FindParam(&indices,&nindices,IndicesEnum);
1134
1135 /*Get rank from the stat comm communicator:*/
1136 IssmComm::SetComm(statcomm);
1137 int my_rank=IssmComm::GetRank();
1138
1139 /*Open files and read them complelety, in a distributed way:*/
1140 range=DetermineLocalSize(nsamples,IssmComm::GetComm());
1141 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
1142
1143 /*Initialize arrays:*/
1144 xs=xNew<IssmDouble*>(nfields*nsteps);
1145 xtype=xNew<int>(nfields*nsteps);
1146 xsize=xNew<int>(nfields*nsteps);
1147
1148 /*Start opening files:*/
1149 for (int i=(lower_row+1);i<=upper_row;i++){
1150 _printf0_("reading file #: " << i << "\n");
1151 char file[1000];
1152 long int length;
1153 char* buffer=NULL;
1154
1155 /*string:*/
1156 sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
1157
1158 /*open file: */
1159 _printf0_(" opening file:\n");
1160 FILE* fid=fopen(file,"rb");
1161
1162 /*figure out size of file, and read the whole thing:*/
1163 _printf0_(" reading file:\n");
1164 fseek (fid, 0, SEEK_END);
1165 length = ftell (fid);
1166 fseek (fid, 0, SEEK_SET);
1167 buffer = xNew<char>(length);
1168 fread (buffer, sizeof(char), length, fid);
1169
1170 /*close file:*/
1171 fclose (fid);
1172
1173 /*create a memory stream with this buffer:*/
1174 _printf0_(" processing file:\n");
1175 fid=fmemopen(buffer, length, "rb");
1176
1177 /*start reading data from the buffer directly:*/
1178 for (int f=0;f<nfields;f++){
1179 fseek(fid,0,SEEK_SET);
1180 char* field=fields[f];
1181 for (int j=0;j<nsteps;j++){
1182 int counter=f*nsteps+j;
1183 xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
1184 if(i==(lower_row+1)){
1185 if(xtype[counter]==1){
1186 x=xNew<IssmDouble>(range);
1187 x[0]=scalar;
1188 xs[counter]=x;
1189 xsize[counter]=range;
1190 }
1191 else if (xtype[counter]==3){
1192 x=xNew<IssmDouble>(nindices*range);
1193 for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
1194 xs[counter]=x;
1195 xsize[counter]=range*nindices;
1196 }
1197 else _error_("cannot carry out statistics on type " << xtype[counter]);
1198 }
1199 else{
1200 if(xtype[counter]==1){
1201 x=xs[counter];
1202 x[i-(lower_row+1)]=scalar;
1203 xs[counter]=x;
1204 }
1205 else if (xtype[counter]==3){
1206 x=xs[counter];
1207 for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
1208 xs[counter]=x;
1209 }
1210 else _error_("cannot carry out statistics on type " << xtype[counter]);
1211 }
1212 }
1213 }
1214 fclose(fid);
1215
1216 /*delete buffer:*/
1217 xDelete<char>(buffer);
1218 }
1219 ISSM_MPI_Barrier(IssmComm::GetComm());
1220 _printf0_("Done reading files, now assembling time series.\n");
1221
1222 for (int f=0;f<nfields;f++){
1223 for (int j=0;j<nsteps;j++){
1224 int counter=f*nsteps+j;
1225 if (xtype[counter]==1){
1226 /*we are broadcasting range times doubles:*/
1227 x=xs[counter];
1228 allx=xNew<IssmDouble>(nsamples);
1229 MPI_Gather(x, range, ISSM_MPI_PDOUBLE,allx, range, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
1230 /*add to results:*/
1231 if(my_rank==0){
1232 char fieldname[1000];
1233
1234 sprintf(fieldname,"%s%s",fields[f],"Samples");
1235 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,1,j+1,0));
1236 }
1237 }
1238 else{
1239 /*we are broadcasting double arrays:*/
1240 x=xs[counter];
1241 allx=xNew<IssmDouble>(nsamples*nindices);
1242
1243 MPI_Gather(x, range*nindices, ISSM_MPI_PDOUBLE,allx, range*nindices, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
1244
1245 /*add to results:*/
1246 if(my_rank==0){
1247 char fieldname[1000];
1248 sprintf(fieldname,"%s%s",fields[f],"Samples");
1249 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,nindices,j+1,0));
1250 }
1251 }
1252 }
1253 }
1254 _printf0_("Done with SampleSeries:\n");
1255 IssmComm::SetComm(ISSM_MPI_COMM_WORLD);
1256
1257 return 1;
1258} /*}}}*/
1259int OutputStatistics(Parameters* parameters,Results* results,int color,ISSM_MPI_Comm statcomm){ /*{{{*/
1260
1261 char outputfilename[1000];
1262 char* directory=NULL;
1263 char* model=NULL;
1264 char* method=NULL;
1265 int nsamples;
1266 int* steps=NULL;
1267 int nsteps;
1268
1269 /*only work on the statistical communicator: */
1270 if (color==MPI_UNDEFINED)return 0;
1271
1272 FemModel* femmodel=new FemModel();
1273
1274 /*Some parameters that will allow us to use the OutputResultsx module:*/
1275 parameters->AddObject(new BoolParam(QmuIsdakotaEnum,false));
1276 parameters->AddObject(new BoolParam(SettingsIoGatherEnum,true));
1277
1278 parameters->FindParam(&directory,DirectoryNameEnum);
1279 parameters->FindParam(&model,InputFileNameEnum);
1280 parameters->FindParam(&nsamples,QmuNsampleEnum);
1281 parameters->FindParam(&steps,&nsteps,StepsEnum);
1282
1283 sprintf(outputfilename,"%s/%s.stats",directory,model);
1284 parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
1285
1286 /*Call OutputResults module:*/
1287 femmodel->parameters=parameters;
1288 femmodel->results=results;
1289
1290 OutputResultsx(femmodel);
1291
1292 return 1;
1293} /*}}}*/
1294int readdata(IssmDouble** pdoublemat, int* pdoublematsize, IssmDouble* pdouble, FILE* fid,char* field,int step){ /*{{{*/
1295
1296 int length;
1297 char fieldname[1000];
1298 int fieldname_size;
1299 IssmDouble rtime;
1300 int rstep;
1301 int M,N;
1302
1303 //fields that we retrive:
1304 IssmDouble dfield;
1305 char* sfield = NULL;
1306 IssmDouble* dmatfield = NULL;
1307 int* imatfield = NULL;
1308
1309 //type of the returned field:
1310 int type;
1311 int found=0;
1312
1313 while(1){
1314
1315 size_t ret_code = fread(&fieldname_size, sizeof(int), 1, fid);
1316 if(ret_code != 1) break; //we are done.
1317
1318 fread(fieldname, sizeof(char), fieldname_size, fid);
1319 //_printf0_("fieldname: " << fieldname << "\n");
1320
1321 fread(&rtime, sizeof(IssmDouble), 1, fid);
1322 fread(&rstep, sizeof(int), 1, fid);
1323
1324 //check on field:
1325 if ((step==rstep) && (strcmp(field,fieldname)==0)){
1326
1327 //ok, go read the result really:
1328 fread(&type,sizeof(int),1,fid);
1329 fread(&M,sizeof(int),1,fid);
1330 if (type==1){
1331 fread(&dfield,sizeof(IssmDouble),1,fid);
1332 }
1333 else if (type==2){
1334 fread(&M,sizeof(int),1,fid);
1335 sfield=xNew<char>(M);
1336 fread(sfield,sizeof(char),M,fid);
1337 }
1338 else if (type==3){
1339 fread(&N,sizeof(int),1,fid);
1340 dmatfield=xNew<IssmDouble>(M*N);
1341 fread(dmatfield,sizeof(IssmDouble),M*N,fid);
1342 }
1343 else if (type==4){
1344 fread(&N,sizeof(int),1,fid);
1345 imatfield=xNew<int>(M*N);
1346 fread(imatfield,sizeof(int),M*N,fid);
1347 }
1348 else _error_("cannot read data of type " << type << "\n");
1349 found=1;
1350 break;
1351 }
1352 else{
1353 //just skim to next results.
1354 fread(&type,sizeof(int),1,fid);
1355 fread(&M,sizeof(int),1,fid);
1356 if (type==1){
1357 fseek(fid,sizeof(IssmDouble),SEEK_CUR);
1358 }
1359 else if(type==2){
1360 fseek(fid,M*sizeof(char),SEEK_CUR);
1361 }
1362 else if(type==3){
1363 fread(&N,sizeof(int),1,fid);
1364 fseek(fid,M*N*sizeof(IssmDouble),SEEK_CUR);
1365 }
1366 else if(type==4){
1367 fread(&N,sizeof(int),1,fid);
1368 fseek(fid,M*N*sizeof(int),SEEK_CUR);
1369 }
1370 else _error_("cannot read data of type " << type << "\n");
1371 }
1372 }
1373 if(found==0)_error_("cound not find " << field << " at step " << step << "\n");
1374
1375 /*assign output pointers:*/
1376 *pdoublemat=dmatfield;
1377 *pdoublematsize=M*N;
1378 *pdouble=dfield;
1379
1380 /*return:*/
1381 return type;
1382
1383}
1384/*}}}*/
1385bool DakotaDirStructure(int argc,char** argv){ /*{{{*/
1386
1387 char* input_file;
1388 FILE* fid;
1389 IoModel* iomodel=NULL;
1390 int check;
1391
1392 //qmu statistics
1393 bool statistics = false;
1394 int numdirectories = 0;
1395
1396 /*First things first, set the communicator as a global variable: */
1397 IssmComm::SetComm(MPI_COMM_WORLD);
1398
1399 /*Barrier:*/
1400 ISSM_MPI_Barrier(IssmComm::GetComm());
1401 _printf0_("Preparing directory structure for model outputs:" << "\n");
1402
1403 //open model input file for reading
1404 input_file=xNew<char>((strlen(argv[2])+strlen(argv[3])+strlen(".bin")+2));
1405 sprintf(input_file,"%s/%s%s",argv[2],argv[3],".bin");
1406 fid=fopen(input_file,"rb");
1407 if (fid==NULL) Cerr << "dirstructure error message: could not open model " << input_file << " to retrieve qmu statistics parameters" << std::endl;
1408
1409 //initialize IoModel, but light version, we just need it to fetch one constant:
1410 iomodel=new IoModel();
1411 iomodel->fid=fid;
1412 iomodel->FetchConstants();
1413
1414 //early return if statistics not requested:
1415 iomodel->FindConstant(&statistics,"md.qmu.statistics");
1416 if(!statistics){
1417 delete iomodel;
1418 xDelete<char>(input_file);
1419 fclose(fid);
1420 return false; //important return value!
1421 }
1422
1423 iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories");
1424
1425 /*Ok, we have everything we need to create the directory structure:*/
1426 if(IssmComm::GetRank()==0){
1427 for (int i=0;i<numdirectories;i++){
1428 char directory[1000];
1429 sprintf(directory,"./%i",i+1);
1430
1431 check = mkdir(directory,ACCESSPERMS);
1432 if (check) _error_("dirstructure error message: could not create directory " << directory << "\n");
1433 }
1434 }
1435
1436 /*Delete resources:*/
1437 delete iomodel;
1438 xDelete<char>(input_file);
1439
1440 //close model file:
1441 fclose(fid);
1442
1443 //return value:
1444 return true; //statistics computation on!
1445} /*}}}*/
Note: See TracBrowser for help on using the repository browser.