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

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

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