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

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

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