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

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

CHG: time mean.
File size: 43.1 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 /*Store for later use in histogram computation:*/
380 maxmeans[f]=allmax;
381 minmeans[f]=allmin;
382
383 } /*}}}*/
384 }
385
386 /*Now that we have the min and max, we can start binning. First allocate
387 * histograms, then start filling them:*/
388 IssmDouble** histogram=xNew<IssmDouble*>(nfields*nsteps);
389 IssmDouble** timehistogram=xNew<IssmDouble*>(nfields);
390
391 _printf0_("Start reading files again, this time binning values in the histogram:\n");
392 /*Start opening files:*/
393 for (int i=(lower_row+1);i<=upper_row;i++){
394 _printf0_("reading file #: " << i << "\n");
395 char file[1000];
396 long int length;
397 char* buffer=NULL;
398
399 /*string:*/
400 sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
401
402 /*open file: */
403 _printf0_(" opening file:\n");
404 FILE* fid=fopen(file,"rb");
405 if(fid==NULL)_error_("cound not open file: " << file << "\n");
406
407 /*figure out size of file, and read the whole thing:*/
408 _printf0_(" reading file:\n");
409 fseek (fid, 0, SEEK_END);
410 length = ftell (fid);
411 fseek (fid, 0, SEEK_SET);
412 buffer = xNew<char>(length);
413 fread (buffer, sizeof(char), length, fid);
414
415 /*close file:*/
416 fclose (fid);
417
418 /*create a memory stream with this buffer:*/
419 _printf0_(" processing file:\n");
420 fid=fmemopen(buffer, length, "rb");
421
422 /*start reading data from the buffer directly:*/
423 for (int f=0;f<nfields;f++){
424 char* field=fields[f];
425 fseek(fid,0,SEEK_SET);
426 for (int j=0;j<nsteps;j++){
427 int counter=f*nsteps+j;
428 xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
429 if(i==(lower_row+1)){
430 if(xtype[counter]==1){
431 IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(nbins);
432 IssmDouble ma=*maxxs[counter];
433 IssmDouble mi=*minxs[counter];
434 int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index--;
435 if(ma==mi)index=0;
436 //_printf_( index << "|" << scalar << "|" << mi << "|" << ma << "|" << nbins << "\n");
437 localhistogram[index]++;
438 histogram[counter]=localhistogram;
439 }
440 else if (xtype[counter]==3){
441 IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(doublematsize*nbins);
442 IssmDouble* ma=maxxs[counter];
443 IssmDouble* mi=minxs[counter];
444 for (int k=0;k<doublematsize;k++){
445 IssmDouble scalar=doublemat[k];
446 int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index--;
447 if (mi[k]==ma[k])index=0;
448 _assert_(scalar<=ma[k]); _assert_(scalar>=mi[k]); _assert_(index<nbins);
449 localhistogram[k*nbins+index]++;
450 }
451 histogram[counter]=localhistogram;
452 xDelete<IssmDouble>(doublemat);
453 }
454 else _error_("cannot carry out statistics on type " << xtype[counter]);
455 }
456 else{
457 if(xtype[counter]==1){
458 IssmDouble* localhistogram=histogram[counter];
459 IssmDouble ma=*maxxs[counter];
460 IssmDouble mi=*minxs[counter];
461 int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
462 if(ma==mi)index=0;
463 localhistogram[index]++;
464 }
465 else if (xtype[counter]==3){
466 IssmDouble* localhistogram=histogram[counter];
467 IssmDouble* ma=maxxs[counter];
468 IssmDouble* mi=minxs[counter];
469 for (int k=0;k<doublematsize;k++){
470 IssmDouble scalar=doublemat[k];
471 int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
472 if (mi[k]==ma[k])index=0;
473
474 localhistogram[k*nbins+index]++;
475 }
476 xDelete<IssmDouble>(doublemat);
477 }
478 else _error_("cannot carry out statistics on type " << xtype[counter]);
479 }
480 }
481 }
482 _printf0_(" average in time:\n");
483
484 /*Deal with average in time: */
485 for (int f=0;f<nfields;f++){
486 fseek(fid,0,SEEK_SET);
487 char* field=fields[f];
488 meanxtype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
489
490 if(meanxtype[f]==1){
491 IssmDouble timemean=0;
492 fseek(fid,0,SEEK_SET);
493 for (int j=0;j<nsteps;j++){
494 readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
495 timemean+=scalar/nsteps;
496 }
497
498 /*Figure out max and min of time means: */
499 if(i==(lower_row+1)){
500 IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(nbins);
501 IssmDouble ma=*maxmeans[f];
502 IssmDouble mi=*minmeans[f];
503 int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
504 if(ma==mi)index=0;
505 localhistogram[index]++;
506 timehistogram[f]=localhistogram;
507 }
508 else{
509 IssmDouble* localhistogram=timehistogram[f];
510 IssmDouble ma=*maxmeans[f];
511 IssmDouble mi=*minmeans[f];
512 int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
513 if(ma==mi)index=0;
514 localhistogram[index]++;
515 }
516 }
517 else{
518 fseek(fid,0,SEEK_SET);
519 IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
520 for (int j=0;j<nsteps;j++){
521 readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
522 for (int k=0;k<doublematsize;k++){
523 timemean[k]+=doublemat[k]/nsteps;
524 }
525 xDelete<IssmDouble>(doublemat);
526 }
527
528 if(i==(lower_row+1)){
529 IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(doublematsize*nbins);
530 IssmDouble* ma=maxmeans[f];
531 IssmDouble* mi=minmeans[f];
532
533 for (int k=0;k<doublematsize;k++){
534 IssmDouble scalar=timemean[k];
535 int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
536 if (mi[k]==ma[k])index=0;
537 localhistogram[k*nbins+index]++;
538 }
539 timehistogram[f]=localhistogram;
540 }
541 else{
542
543 IssmDouble* localhistogram=timehistogram[f];
544 IssmDouble* ma=maxmeans[f];
545 IssmDouble* mi=minmeans[f];
546
547 for (int k=0;k<doublematsize;k++){
548 IssmDouble scalar=timemean[k];
549 int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
550 if (mi[k]==ma[k])index=0;
551
552 localhistogram[k*nbins+index]++;
553 }
554 }
555 }
556 }
557 fclose(fid);
558
559 /*delete buffer:*/
560 xDelete<char>(buffer);
561 }
562 _printf0_("Start aggregating histogram:\n");
563
564 /*We have agregated histograms across the cluster, now gather them across the cluster onto
565 *cpu0: */
566 for (int f=0;f<nfields;f++){
567 int counter0=f*nsteps+0;
568 if (xtype[counter0]==1){ /*deal with scalars {{{*/
569 for (int j=0;j<nsteps;j++){
570 int counter=f*nsteps+j;
571
572 /*we are broadcasting doubles:*/
573 IssmDouble* histo=histogram[counter]; //size nbins
574 IssmDouble* allhisto=xNewZeroInit<IssmDouble>(nbins);
575
576 ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
577
578 /*add to results:*/
579 if(my_rank==0){
580 char fieldname[1000];
581
582 sprintf(fieldname,"%s%s",fields[f],"Histogram");
583 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,1,nbins,j+1,0));
584
585 sprintf(fieldname,"%s%s",fields[f],"Max");
586 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*maxxs[counter],j+1,0));
587 sprintf(fieldname,"%s%s",fields[f],"Min");
588 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*minxs[counter],j+1,0));
589 }
590 }
591 } /*}}}*/
592 else{ /*deal with arrays:{{{*/
593
594 int size=xsize[counter0];
595 for (int j=0;j<nsteps;j++){
596 int counter=f*nsteps+j;
597
598 /*we are broadcasting double arrays:*/
599 IssmDouble* histo=histogram[counter];
600 IssmDouble* allhisto=xNew<IssmDouble>(size*nbins);
601
602 ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
603 xDelete<IssmDouble>(histo);
604
605 /*add to results:*/
606 if(my_rank==0){
607 char fieldname[1000];
608
609 sprintf(fieldname,"%s%s",fields[f],"Histogram");
610 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,size,nbins,j+1,0));
611
612 sprintf(fieldname,"%s%s",fields[f],"Max");
613 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,maxxs[counter],size,1,j+1,0));
614 sprintf(fieldname,"%s%s",fields[f],"Min");
615 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,minxs[counter],size,1,j+1,0));
616 }
617 }
618 } /*}}}*/
619 }
620 _printf0_("Start aggregating time mean histogram:\n");
621
622 /*Now do the same for the time mean fields:*/
623 for (int f=0;f<nfields;f++){
624 if (meanxtype[f]==1){ /*deal with scalars {{{*/
625
626 /*we are broadcasting doubles:*/
627 IssmDouble* histo=timehistogram[f];
628 IssmDouble* allhisto=xNewZeroInit<IssmDouble>(nbins);
629
630 ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
631
632 /*add to results at time step 1:*/
633 if(my_rank==0){
634 char fieldname[1000];
635
636 sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram");
637 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,1,nbins,steps[0],0));
638
639 sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax");
640 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*maxmeans[f],steps[0],0));
641 sprintf(fieldname,"%s%s",fields[f],"TimeMeaMin");
642 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*minmeans[f],steps[0],0));
643 }
644 } /*}}}*/
645 else{ /*deal with arrays:{{{*/
646
647 int size=meanxsize[f];
648
649 /*we are broadcasting double arrays:*/
650 IssmDouble* histo=timehistogram[f];
651 IssmDouble* allhisto=xNewZeroInit<IssmDouble>(size*nbins);
652
653 ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
654 xDelete<IssmDouble>(histo);
655 /*add to results at step 1:*/
656 if(my_rank==0){
657 char fieldname[1000];
658
659 sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram");
660 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,size,nbins,steps[0],0));
661 sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax");
662 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,maxmeans[f],size,1,steps[0],0));
663 sprintf(fieldname,"%s%s",fields[f],"TimeMeanMin");
664 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,minmeans[f],size,1,steps[0],0));
665 }
666 } /*}}}*/
667 }
668 _printf0_("Done aggregating time mean histogram:\n");
669 IssmComm::SetComm(ISSM_MPI_COMM_WORLD);
670
671 return 1;
672}
673/*}}}*/
674int ComputeMeanVariance(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){ /*{{{*/
675
676 int nsamples;
677 char* directory=NULL;
678 char* model=NULL;
679 char** fields=NULL;
680 int* steps=NULL;
681 int nsteps;
682 int nfields;
683 int range,lower_row,upper_row;
684 int nfilesperdirectory;
685
686 /*intermediary:*/
687 IssmDouble* doublemat=NULL;
688 int doublematsize;
689 IssmDouble scalar;
690
691 /*computation of average and variance itself:*/
692 IssmDouble* x = NULL;
693 IssmDouble* x2 = NULL;
694 IssmDouble** xs = NULL;
695 IssmDouble** xs2 = NULL;
696 int* xtype=NULL;
697 int* xsize=NULL;
698
699 IssmDouble** meanx=NULL;
700 IssmDouble** meanx2=NULL;
701 int* meantype=NULL;
702 int* meansize=NULL;
703
704 /*only work on the statistical communicator: */
705 if (color==MPI_UNDEFINED)return 0;
706
707 /*Retrieve parameters:*/
708 parameters->FindParam(&nfilesperdirectory,QmuNfilesPerDirectoryEnum);
709 parameters->FindParam(&nsamples,QmuNsampleEnum);
710 parameters->FindParam(&directory,DirectoryNameEnum);
711 parameters->FindParam(&model,InputFileNameEnum);
712 parameters->FindParam(&fields,&nfields,FieldsEnum);
713 parameters->FindParam(&steps,&nsteps,StepsEnum);
714
715 /*Get rank from the stat comm communicator:*/
716 IssmComm::SetComm(statcomm);
717 int my_rank=IssmComm::GetRank();
718
719 /*Open files and read them complelety, in a distributed way:*/
720 range=DetermineLocalSize(nsamples,IssmComm::GetComm());
721 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
722
723 /*Initialize arrays:*/
724 xs=xNew<IssmDouble*>(nfields*nsteps);
725 xs2=xNew<IssmDouble*>(nfields*nsteps);
726 xtype=xNew<int>(nfields*nsteps);
727 xsize=xNew<int>(nfields*nsteps);
728
729 meantype=xNew<int>(nfields);
730 meansize=xNew<int>(nfields);
731 meanx=xNew<IssmDouble*>(nfields);
732 meanx2=xNew<IssmDouble*>(nfields);
733
734 /*Start opening files:*/
735 for (int i=(lower_row+1);i<=upper_row;i++){
736 _printf0_("reading file #: " << i << "\n");
737 char file[1000];
738 long int length;
739 char* buffer=NULL;
740
741 /*string:*/
742 sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
743
744 /*open file: */
745 _printf0_(" opening file: " << file << "\n");
746 FILE* fid=fopen(file,"rb");
747 if(fid==NULL) _error_(" could not open file: " << file << "\n");
748
749 /*figure out size of file, and read the whole thing:*/
750 _printf0_(" reading file:\n");
751 fseek (fid, 0, SEEK_END);
752 length = ftell (fid);
753 fseek (fid, 0, SEEK_SET);
754 buffer = xNew<char>(length);
755 fread (buffer, sizeof(char), length, fid);
756
757 /*close file:*/
758 fclose (fid);
759
760 /*create a memory stream with this buffer:*/
761 _printf0_(" processing file:\n");
762 fid=fmemopen(buffer, length, "rb");
763
764 /*start reading data from the buffer directly:*/
765 for (int f=0;f<nfields;f++){
766 char* field=fields[f];
767 fseek(fid,0,SEEK_SET);
768 for (int j=0;j<nsteps;j++){
769 int counter=f*nsteps+j;
770 xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
771 if(i==(lower_row+1)){
772 if(xtype[counter]==1){
773 xs[counter]=xNew<IssmDouble>(1);
774 xs2[counter]=xNew<IssmDouble>(1);
775 *xs[counter]=scalar;
776 *xs2[counter]=pow(scalar,2.0);
777 xsize[counter]=1;
778 }
779 else if (xtype[counter]==3){
780 IssmDouble* doublemat2=xNew<IssmDouble>(doublematsize);
781 for(int k=0;k<doublematsize;k++)doublemat2[k]=pow(doublemat[k],2.0);
782 xs[counter]=doublemat;
783 xs2[counter]=doublemat2;
784 xsize[counter]=doublematsize;
785 }
786 else _error_("cannot carry out statistics on type " << xtype[counter]);
787 }
788 else{
789 if(xtype[counter]==1){
790 *xs[counter]+=scalar;
791 *xs2[counter]+=pow(scalar,2.0);
792 }
793 else if (xtype[counter]==3){
794 IssmDouble* newdoublemat=xs[counter];
795 IssmDouble* newdoublemat2=xs2[counter];
796 for(int k=0;k<doublematsize;k++){
797 newdoublemat[k]+=doublemat[k];
798 newdoublemat2[k]+=pow(doublemat[k],2.0);
799 }
800 xs[counter]=newdoublemat;
801 xs2[counter]=newdoublemat2;
802 }
803 else _error_("cannot carry out statistics on type " << xtype[counter]);
804 }
805 }
806 }
807
808 /*Deal with time mean: */
809 for (int f=0;f<nfields;f++){
810 char* field=fields[f];
811 fseek(fid,0,SEEK_SET);
812 meantype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
813 if(i==(lower_row+1)){
814 if(meantype[f]==1){
815 meanx[f]=xNewZeroInit<IssmDouble>(1);
816 meanx2[f]=xNewZeroInit<IssmDouble>(1);
817 meansize[f]=1;
818 }
819 else{
820 meanx[f]=xNewZeroInit<IssmDouble>(doublematsize);
821 meanx2[f]=xNewZeroInit<IssmDouble>(doublematsize);
822 meansize[f]=doublematsize;
823 }
824 }
825 fseek(fid,0,SEEK_SET);
826 if(meantype[f]==1){
827 IssmDouble sc=0;
828 IssmDouble sc2=0;
829 for(int j=0;j<nsteps;j++){
830 readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
831 sc+=scalar/nsteps;
832 }
833 sc2+=pow(sc,2.0);
834 *meanx[f]+=sc;
835 *meanx2[f]+=sc2;
836 }
837 else{
838 IssmDouble* sc=meanx[f];
839 IssmDouble* sc2=meanx2[f];
840 IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
841 IssmDouble* timemean2=xNewZeroInit<IssmDouble>(doublematsize);
842
843 for(int j=0;j<nsteps;j++){
844 readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
845 for (int k=0;k<doublematsize;k++){
846 timemean[k]+=doublemat[k]/nsteps;
847 }
848 }
849 for (int k=0;k<doublematsize;k++){
850 timemean2[k]=pow(timemean[k],2.0);
851 }
852 for (int k=0;k<doublematsize;k++){
853 sc[k]+=timemean[k];
854 sc2[k]+=timemean2[k];
855 }
856
857 }
858
859 }
860 fclose(fid);
861
862 /*delete buffer:*/
863 xDelete<char>(buffer);
864 }
865 ISSM_MPI_Barrier(IssmComm::GetComm());
866 _printf0_("Done reading files, now computing mean and variance.\n");
867
868 /*We have agregated x and x^2 across the cluster, now gather across the cluster onto
869 *cpu0 and then compute statistics:*/
870 for (int f=0;f<nfields;f++){
871 int counter0=f*nsteps+0;
872 if (xtype[counter0]==1){ /*deal with scalars {{{*/
873 IssmDouble mean,stddev;
874 for (int j=0;j<nsteps;j++){
875 int counter=f*nsteps+j;
876
877 /*we are broadcasting doubles:*/
878 IssmDouble scalar=*xs[counter];
879 IssmDouble scalar2=*xs2[counter];
880 IssmDouble sumscalar;
881 IssmDouble sumscalar2;
882
883 ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
884 ISSM_MPI_Reduce(&scalar2,&sumscalar2,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
885 /*Build average and standard deviation. For standard deviation, use the
886 *following formula: sigma^2=E(x^2)-mu^2:*/
887 mean=sumscalar/nsamples;
888 stddev=sqrt(sumscalar2/nsamples-pow(mean,2.0));
889
890 /*add to results:*/
891 if(my_rank==0){
892 char fieldname[1000];
893
894 sprintf(fieldname,"%s%s",fields[f],"Mean");
895 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,j+1,0));
896 sprintf(fieldname,"%s%s",fields[f],"Stddev");
897 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,j+1,0));
898 }
899
900 }
901 } /*}}}*/
902 else{ /*deal with arrays:{{{*/
903
904 int size=xsize[counter0];
905
906 IssmDouble* mean=xNew<IssmDouble>(size);
907 IssmDouble* stddev=xNew<IssmDouble>(size);
908
909 for (int j=0;j<nsteps;j++){
910 int counter=f*nsteps+j;
911
912 /*we are broadcasting double arrays:*/
913 x=xs[counter];
914 x2=xs2[counter];
915
916 IssmDouble* sumx=xNew<IssmDouble>(size);
917 IssmDouble* sumx2=xNew<IssmDouble>(size);
918
919 ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
920 ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
921
922 /*Build average and standard deviation. For standard deviation, use the
923 *following formula: sigma^2=E(x^2)-mu^2:*/
924 for (int k=0;k<size;k++){
925 mean[k]=sumx[k]/nsamples;
926 stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
927 }
928
929 /*add to results:*/
930 if(my_rank==0){
931 char fieldname[1000];
932
933 sprintf(fieldname,"%s%s",fields[f],"Mean");
934 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,j+1,0));
935 sprintf(fieldname,"%s%s",fields[f],"Stddev");
936 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,j+1,0));
937 }
938 }
939 } /*}}}*/
940 }
941 /*Do the same but for the time mean:*/
942 for (int f=0;f<nfields;f++){
943 if (meantype[f]==1){ /*deal with scalars {{{*/
944 IssmDouble mean,stddev;
945
946 /*we are broadcasting doubles:*/
947 IssmDouble scalar=*meanx[f];
948 IssmDouble scalar2=*meanx2[f];
949 IssmDouble sumscalar;
950 IssmDouble sumscalar2;
951
952 ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
953 ISSM_MPI_Reduce(&scalar2,&sumscalar2,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
954 /*Build average and standard deviation. For standard deviation, use the
955 *following formula: sigma^2=E(x^2)-mu^2:*/
956 mean=sumscalar/nsamples;
957 stddev=sqrt(sumscalar2/nsamples-pow(mean,2.0));
958
959 /*add to results:*/
960 if(my_rank==0){
961 char fieldname[1000];
962
963 sprintf(fieldname,"%s%s",fields[f],"TimeMean");
964 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,steps[0],0));
965 sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
966 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,steps[0],0));
967 }
968 } /*}}}*/
969 else{ /*deal with arrays:{{{*/
970
971 int size=meansize[f];
972 IssmDouble* mean=xNew<IssmDouble>(size);
973 IssmDouble* stddev=xNew<IssmDouble>(size);
974
975 /*we are broadcasting double arrays:*/
976 x=meanx[f];
977 x2=meanx2[f];
978
979 IssmDouble* sumx=xNew<IssmDouble>(size);
980 IssmDouble* sumx2=xNew<IssmDouble>(size);
981
982 ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
983 ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
984
985 /*Build average and standard deviation. For standard deviation, use the
986 *following formula: sigma^2=E(x^2)-mu^2:*/
987 for (int k=0;k<size;k++){
988 mean[k]=sumx[k]/nsamples;
989 stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
990 }
991
992 /*add to results:*/
993 if(my_rank==0){
994 char fieldname[1000];
995
996 sprintf(fieldname,"%s%s",fields[f],"TimeMean");
997 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,steps[0],0));
998 sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
999 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,steps[0],0));
1000 }
1001 } /*}}}*/
1002 }
1003
1004 _printf0_("Done with MeanVariance:\n");
1005 IssmComm::SetComm(ISSM_MPI_COMM_WORLD);
1006
1007 return 1;
1008} /*}}}*/
1009int ComputeSampleSeries(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){ /*{{{*/
1010
1011 int nsamples;
1012 char* directory=NULL;
1013 char* model=NULL;
1014 char** fields=NULL;
1015 int* steps=NULL;
1016 int nsteps;
1017 int nfields;
1018 int range,lower_row,upper_row;
1019 int nfilesperdirectory;
1020 int* indices=NULL;
1021 int nindices;
1022
1023 /*intermediary:*/
1024 IssmDouble* doublemat=NULL;
1025 int doublematsize;
1026 IssmDouble scalar;
1027
1028 /*computation of average and variance itself:*/
1029 IssmDouble* x = NULL;
1030 IssmDouble* allx=NULL;
1031 IssmDouble** xs = NULL;
1032 int* xtype=NULL;
1033 int* xsize=NULL;
1034
1035 /*only work on the statistical communicator: */
1036 if (color==MPI_UNDEFINED)return 0;
1037
1038 /*Retrieve parameters:*/
1039 parameters->FindParam(&nsamples,QmuNsampleEnum);
1040 parameters->FindParam(&directory,DirectoryNameEnum);
1041 parameters->FindParam(&model,InputFileNameEnum);
1042 parameters->FindParam(&fields,&nfields,FieldsEnum);
1043 parameters->FindParam(&steps,&nsteps,StepsEnum);
1044 parameters->FindParam(&indices,&nindices,IndicesEnum);
1045
1046 /*Get rank from the stat comm communicator:*/
1047 IssmComm::SetComm(statcomm);
1048 int my_rank=IssmComm::GetRank();
1049
1050 /*Open files and read them complelety, in a distributed way:*/
1051 range=DetermineLocalSize(nsamples,IssmComm::GetComm());
1052 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
1053
1054 /*Initialize arrays:*/
1055 xs=xNew<IssmDouble*>(nfields*nsteps);
1056 xtype=xNew<int>(nfields*nsteps);
1057 xsize=xNew<int>(nfields*nsteps);
1058
1059 /*Start opening files:*/
1060 for (int i=(lower_row+1);i<=upper_row;i++){
1061 _printf0_("reading file #: " << i << "\n");
1062 char file[1000];
1063 long int length;
1064 char* buffer=NULL;
1065
1066 /*string:*/
1067 sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
1068
1069 /*open file: */
1070 _printf0_(" opening file:\n");
1071 FILE* fid=fopen(file,"rb");
1072
1073 /*figure out size of file, and read the whole thing:*/
1074 _printf0_(" reading file:\n");
1075 fseek (fid, 0, SEEK_END);
1076 length = ftell (fid);
1077 fseek (fid, 0, SEEK_SET);
1078 buffer = xNew<char>(length);
1079 fread (buffer, sizeof(char), length, fid);
1080
1081 /*close file:*/
1082 fclose (fid);
1083
1084 /*create a memory stream with this buffer:*/
1085 _printf0_(" processing file:\n");
1086 fid=fmemopen(buffer, length, "rb");
1087
1088 /*start reading data from the buffer directly:*/
1089 for (int f=0;f<nfields;f++){
1090 fseek(fid,0,SEEK_SET);
1091 char* field=fields[f];
1092 for (int j=0;j<nsteps;j++){
1093 int counter=f*nsteps+j;
1094 xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
1095 if(i==(lower_row+1)){
1096 if(xtype[counter]==1){
1097 x=xNew<IssmDouble>(range);
1098 x[0]=scalar;
1099 xs[counter]=x;
1100 xsize[counter]=range;
1101 }
1102 else if (xtype[counter]==3){
1103 x=xNew<IssmDouble>(nindices*range);
1104 for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
1105 xs[counter]=x;
1106 xsize[counter]=range*nindices;
1107 }
1108 else _error_("cannot carry out statistics on type " << xtype[counter]);
1109 }
1110 else{
1111 if(xtype[counter]==1){
1112 x=xs[counter];
1113 x[i-(lower_row+1)]=scalar;
1114 xs[counter]=x;
1115 }
1116 else if (xtype[counter]==3){
1117 x=xs[counter];
1118 for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
1119 xs[counter]=x;
1120 }
1121 else _error_("cannot carry out statistics on type " << xtype[counter]);
1122 }
1123 }
1124 }
1125 fclose(fid);
1126
1127 /*delete buffer:*/
1128 xDelete<char>(buffer);
1129 }
1130 ISSM_MPI_Barrier(IssmComm::GetComm());
1131 _printf0_("Done reading files, now assembling time series.\n");
1132
1133 for (int f=0;f<nfields;f++){
1134 for (int j=0;j<nsteps;j++){
1135 int counter=f*nsteps+j;
1136 if (xtype[counter]==1){
1137 /*we are broadcasting range times doubles:*/
1138 x=xs[counter];
1139 allx=xNew<IssmDouble>(nsamples);
1140 MPI_Gather(x, range, ISSM_MPI_PDOUBLE,allx, range, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
1141 /*add to results:*/
1142 if(my_rank==0){
1143 char fieldname[1000];
1144
1145 sprintf(fieldname,"%s%s",fields[f],"Samples");
1146 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,1,j+1,0));
1147 }
1148 }
1149 else{
1150 /*we are broadcasting double arrays:*/
1151 x=xs[counter];
1152 allx=xNew<IssmDouble>(nsamples*nindices);
1153
1154 MPI_Gather(x, range*nindices, ISSM_MPI_PDOUBLE,allx, range*nindices, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
1155
1156 /*add to results:*/
1157 if(my_rank==0){
1158 char fieldname[1000];
1159 sprintf(fieldname,"%s%s",fields[f],"Samples");
1160 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,nindices,j+1,0));
1161 }
1162 }
1163 }
1164 }
1165 _printf0_("Done with SampleSeries:\n");
1166 IssmComm::SetComm(ISSM_MPI_COMM_WORLD);
1167
1168 return 1;
1169} /*}}}*/
1170int OutputStatistics(Parameters* parameters,Results* results,int color,ISSM_MPI_Comm statcomm){ /*{{{*/
1171
1172 char outputfilename[1000];
1173 char* directory=NULL;
1174 char* model=NULL;
1175 char* method=NULL;
1176 int nsamples;
1177 int* steps=NULL;
1178 int nsteps;
1179
1180 /*only work on the statistical communicator: */
1181 if (color==MPI_UNDEFINED)return 0;
1182
1183 FemModel* femmodel=new FemModel();
1184
1185 /*Some parameters that will allow us to use the OutputResultsx module:*/
1186 parameters->AddObject(new BoolParam(QmuIsdakotaEnum,false));
1187 parameters->AddObject(new BoolParam(SettingsIoGatherEnum,true));
1188
1189 parameters->FindParam(&directory,DirectoryNameEnum);
1190 parameters->FindParam(&model,InputFileNameEnum);
1191 parameters->FindParam(&nsamples,QmuNsampleEnum);
1192 parameters->FindParam(&steps,&nsteps,StepsEnum);
1193
1194 sprintf(outputfilename,"%s/%s.stats",directory,model);
1195 parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
1196
1197 /*Call OutputResults module:*/
1198 femmodel->parameters=parameters;
1199 femmodel->results=results;
1200
1201 OutputResultsx(femmodel);
1202
1203 return 1;
1204} /*}}}*/
1205bool DakotaDirStructure(int argc,char** argv){ /*{{{*/
1206
1207 char* input_file;
1208 FILE* fid;
1209 IoModel* iomodel=NULL;
1210 int check;
1211
1212 //qmu statistics
1213 bool statistics = false;
1214 int numdirectories = 0;
1215
1216 /*First things first, set the communicator as a global variable: */
1217 IssmComm::SetComm(MPI_COMM_WORLD);
1218
1219 /*Barrier:*/
1220 ISSM_MPI_Barrier(IssmComm::GetComm());
1221 _printf0_("Preparing directory structure for model outputs:" << "\n");
1222
1223 //open model input file for reading
1224 input_file=xNew<char>((strlen(argv[2])+strlen(argv[3])+strlen(".bin")+2));
1225 sprintf(input_file,"%s/%s%s",argv[2],argv[3],".bin");
1226 fid=fopen(input_file,"rb");
1227 if (fid==NULL) Cerr << "dirstructure error message: could not open model " << input_file << " to retrieve qmu statistics parameters" << std::endl;
1228
1229 //initialize IoModel, but light version, we just need it to fetch one constant:
1230 iomodel=new IoModel();
1231 iomodel->fid=fid;
1232 iomodel->FetchConstants();
1233
1234 //early return if statistics not requested:
1235 iomodel->FindConstant(&statistics,"md.qmu.statistics");
1236 if(!statistics){
1237 delete iomodel;
1238 xDelete<char>(input_file);
1239 fclose(fid);
1240 return false; //important return value!
1241 }
1242
1243 iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories");
1244
1245 /*Ok, we have everything we need to create the directory structure:*/
1246 if(IssmComm::GetRank()==0){
1247 for (int i=0;i<numdirectories;i++){
1248 char directory[1000];
1249 sprintf(directory,"./%i",i+1);
1250
1251 check = mkdir(directory,ACCESSPERMS);
1252 if (check) _error_("dirstructure error message: could not create directory " << directory << "\n");
1253 }
1254 }
1255
1256 /*Delete resources:*/
1257 delete iomodel;
1258 xDelete<char>(input_file);
1259
1260 //close model file:
1261 fclose(fid);
1262
1263 //return value:
1264 return true; //statistics computation on!
1265} /*}}}*/
1266int DakotaStatistics(int argc,char** argv){ /*{{{*/
1267
1268 char* input_file;
1269 FILE* fid;
1270 IoModel* iomodel=NULL;
1271 ISSM_MPI_Comm statcomm;
1272 int my_rank;
1273
1274 //qmu statistics
1275 bool statistics = false;
1276 int numstatistics = 0;
1277 int numdirectories = 0;
1278 int nfilesperdirectory = 0;
1279 char string[1000];
1280 char* name = NULL;
1281 char** fields = NULL;
1282 int nfields;
1283 int* steps=NULL;
1284 int nsteps;
1285 int nbins;
1286 int* indices=NULL;
1287 int nindices;
1288 int nsamples;
1289 int dummy;
1290 char* directory=NULL;
1291 char* model=NULL;
1292 Results* results=NULL;
1293 Parameters* parameters=NULL;
1294 int color;
1295
1296 /*First things first, set the communicator as a global variable: */
1297 IssmComm::SetComm(MPI_COMM_WORLD);
1298 my_rank=IssmComm::GetRank();
1299
1300 /*Barrier:*/
1301 ISSM_MPI_Barrier(IssmComm::GetComm());
1302 _printf0_("Dakota Statistic Computation" << "\n");
1303
1304 //open model input file for reading
1305 input_file=xNew<char>((strlen(argv[2])+strlen(argv[3])+strlen(".bin")+2));
1306 sprintf(input_file,"%s/%s%s",argv[2],argv[3],".bin");
1307 fid=fopen(input_file,"rb");
1308 if (fid==NULL) Cerr << "issm_dakota_statistics error message: could not open model " << input_file << " to retrieve qmu statistics parameters" << std::endl;
1309
1310 //initialize IoModel, but light version, we'll need it to fetch constants:
1311 iomodel=new IoModel();
1312 iomodel->fid=fid;
1313 iomodel->FetchConstants();
1314
1315 //early return if statistics not requested:
1316 iomodel->FindConstant(&statistics,"md.qmu.statistics");
1317 if(!statistics){
1318 delete iomodel;
1319 xDelete<char>(input_file);
1320 fclose(fid);
1321 return 0;
1322 }else{
1323 //create parameters datasets with al the qmu statistics settings we need:
1324
1325 /*Initialize parameters and results:*/
1326 results = new Results();
1327 parameters=new Parameters();
1328
1329 //solution type:
1330 parameters->AddObject(new IntParam(SolutionTypeEnum,StatisticsSolutionEnum));
1331
1332 //root directory
1333 directory=xNew<char>(strlen(argv[2])+1);
1334 xMemCpy<char>(directory,argv[2],strlen(argv[2])+1);
1335 parameters->AddObject(new StringParam(DirectoryNameEnum,directory));
1336
1337 //model name
1338 model=xNew<char>(strlen(argv[3])+1);
1339 xMemCpy<char>(model,argv[3],strlen(argv[3])+1);
1340 parameters->AddObject(new StringParam(InputFileNameEnum,model));
1341
1342 //nsamples
1343 iomodel->FindConstant(&nsamples,"md.qmu.method.params.samples");
1344 parameters->AddObject(new IntParam(QmuNsampleEnum,nsamples));
1345
1346 //ndirectories
1347 iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories");
1348 parameters->AddObject(new IntParam(QmuNdirectoriesEnum,numdirectories));
1349
1350 //nfiles per directory
1351 iomodel->FindConstant(&nfilesperdirectory,"md.qmu.statistics.nfiles_per_directory");
1352 parameters->AddObject(new IntParam(QmuNfilesPerDirectoryEnum,nfilesperdirectory));
1353
1354 //At this point, we don't want to go forward any longer, we want to create an MPI
1355 //communicator on which to carry out the computations:
1356 if ((my_rank+1)*nfilesperdirectory>nsamples)color=MPI_UNDEFINED;
1357 else color=0;
1358 ISSM_MPI_Comm_split(ISSM_MPI_COMM_WORLD,color, my_rank, &statcomm);
1359
1360 iomodel->FindConstant(&numstatistics,"md.qmu.statistics.numstatistics");
1361 for (int i=1;i<=numstatistics;i++){
1362
1363 char* directory=NULL;
1364 char* model=NULL;
1365 int nsamples;
1366 _printf0_("Dealing with qmu statistical computation #" << i << "\n");
1367
1368 sprintf(string,"md.qmu.statistics.method(%i).name",i);
1369 iomodel->FindConstant(&name,string);
1370
1371 sprintf(string,"md.qmu.statistics.method(%i).fields",i);
1372 iomodel->FindConstant(&fields,&nfields,string);
1373 parameters->AddObject(new StringArrayParam(FieldsEnum,fields,nfields));
1374
1375 sprintf(string,"md.qmu.statistics.method(%i).steps",i);
1376 iomodel->FetchData(&steps,&dummy,&nsteps,string);
1377 parameters->AddObject(new IntVecParam(StepsEnum,steps,nsteps));
1378
1379 if (strcmp(name,"Histogram")==0){
1380 /*fetch nbins: */
1381 sprintf(string,"md.qmu.statistics.method(%i).nbins",i);
1382 iomodel->FindConstant(&nbins,string);
1383 parameters->AddObject(new IntParam(NbinsEnum,nbins));
1384 ComputeHistogram(parameters,results,color,statcomm);
1385 }
1386 else if (strcmp(name,"SampleSeries")==0){
1387 /*fetch indices: */
1388 sprintf(string,"md.qmu.statistics.method(%i).indices",i);
1389 iomodel->FetchData(&indices,&dummy,&nindices,string);
1390 parameters->AddObject(new IntVecParam(IndicesEnum,indices,nindices));
1391
1392 ComputeSampleSeries(parameters,results,color,statcomm);
1393 }
1394 else if (strcmp(name,"MeanVariance")==0){
1395 ComputeMeanVariance(parameters,results,color,statcomm);
1396 }
1397 else _error_(" error creating qmu statistics methods parameters: unsupported method " << name);
1398 }
1399
1400 /*Delete resources:*/
1401 xDelete<char>(input_file);
1402 delete iomodel;
1403 }
1404
1405 //close model file:
1406 fclose(fid);
1407
1408 /*output results:*/
1409 OutputStatistics(parameters,results,color,statcomm);
1410
1411 /*all meet here: */
1412 ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); _printf0_("Output file.\n");
1413
1414 /*Delete resources:*/
1415 delete parameters;
1416 delete results;
1417
1418 return 1;
1419} /*}}}*/
Note: See TracBrowser for help on using the repository browser.