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

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

CHG: diverse
File size: 36.2 KB
Line 
1/*!\file: QmuStatisticsx routines
2 */
3/*includes and prototypes:*/
4#include "./QmuStatisticsx.h"
5#include "../OutputResultsx/OutputResultsx.h"
6
7int readdata(IssmDouble** pdoublemat, int* pdoublematsize, IssmDouble* pdouble, FILE* fid,char* field,int step){ /*{{{*/
8
9 int length;
10 char fieldname[1000];
11 int fieldname_size;
12 IssmDouble rtime;
13 int rstep;
14 int M,N;
15
16 //fields that we retrive:
17 IssmDouble dfield;
18 char* sfield = NULL;
19 IssmDouble* dmatfield = NULL;
20 int* imatfield = NULL;
21
22 //type of the returned field:
23 int type;
24 int found=0;
25
26 while(1){
27
28 size_t ret_code = fread(&fieldname_size, sizeof(int), 1, fid);
29 if(ret_code != 1) break; //we are done.
30
31 fread(fieldname, sizeof(char), fieldname_size, fid);
32 //_printf0_("fieldname: " << fieldname << "\n");
33
34 fread(&rtime, sizeof(IssmDouble), 1, fid);
35 fread(&rstep, sizeof(int), 1, fid);
36
37 //check on field:
38 if ((step==rstep) && (strcmp(field,fieldname)==0)){
39
40 //ok, go read the result really:
41 fread(&type,sizeof(int),1,fid);
42 fread(&M,sizeof(int),1,fid);
43 if (type==1){
44 fread(&dfield,sizeof(IssmDouble),1,fid);
45 }
46 else if (type==2){
47 fread(&M,sizeof(int),1,fid);
48 sfield=xNew<char>(M);
49 fread(sfield,sizeof(char),M,fid);
50 }
51 else if (type==3){
52 fread(&N,sizeof(int),1,fid);
53 dmatfield=xNew<IssmDouble>(M*N);
54 fread(dmatfield,sizeof(IssmDouble),M*N,fid);
55 }
56 else if (type==4){
57 fread(&N,sizeof(int),1,fid);
58 imatfield=xNew<int>(M*N);
59 fread(imatfield,sizeof(int),M*N,fid);
60 }
61 else _error_("cannot read data of type " << type << "\n");
62 found=1;
63 break;
64 }
65 else{
66 //just skim to next results.
67 fread(&type,sizeof(int),1,fid);
68 fread(&M,sizeof(int),1,fid);
69 if (type==1){
70 fseek(fid,sizeof(IssmDouble),SEEK_CUR);
71 }
72 else if(type==2){
73 fseek(fid,M*sizeof(char),SEEK_CUR);
74 }
75 else if(type==3){
76 fread(&N,sizeof(int),1,fid);
77 fseek(fid,M*N*sizeof(IssmDouble),SEEK_CUR);
78 }
79 else if(type==4){
80 fread(&N,sizeof(int),1,fid);
81 fseek(fid,M*N*sizeof(int),SEEK_CUR);
82 }
83 else _error_("cannot read data of type " << type << "\n");
84 }
85 }
86 if(found==0)_error_("cound not find " << field << " at step " << step << "\n");
87
88 /*assign output pointers:*/
89 *pdoublemat=dmatfield;
90 *pdoublematsize=M*N;
91 *pdouble=dfield;
92
93 /*return:*/
94 return type;
95
96}
97/*}}}*/
98int ComputeHistogram(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){ /*{{{*/
99
100 int nsamples;
101 char* directory=NULL;
102 char* model=NULL;
103 char** fields=NULL;
104 int* steps=NULL;
105 int nsteps;
106 int nfields;
107 int nbins;
108 int range,lower_row,upper_row;
109 int nfilesperdirectory;
110
111 /*intermediary:*/
112 IssmDouble* doublemat=NULL;
113 int doublematsize;
114 IssmDouble scalar;
115
116 /*computation of average and variance itself:*/
117 IssmDouble** maxxs = NULL;
118 IssmDouble** minxs = NULL;
119 int* xtype=NULL;
120 int* xsize=NULL;
121
122 IssmDouble** maxmeans=NULL;
123 IssmDouble** minmeans=NULL;
124 int* meanxtype=NULL;
125 int* meanxsize=NULL;
126
127 /*only work on the statistical communicator: */
128 if (color==MPI_UNDEFINED)return 0;
129
130 /*Retrieve parameters:*/
131 parameters->FindParam(&nfilesperdirectory,QmuNfilesPerDirectoryEnum);
132 parameters->FindParam(&nsamples,QmuNsampleEnum);
133 parameters->FindParam(&directory,DirectoryNameEnum);
134 parameters->FindParam(&model,InputFileNameEnum);
135 parameters->FindParam(&fields,&nfields,FieldsEnum);
136 parameters->FindParam(&steps,&nsteps,StepsEnum);
137 parameters->FindParam(&nbins,NbinsEnum);
138
139 /*Get rank from the stat comm communicator:*/
140 IssmComm::SetComm(statcomm);
141 int my_rank=IssmComm::GetRank();
142
143 /*Open files and read them complelety, in a distributed way:*/
144 range=DetermineLocalSize(nsamples,IssmComm::GetComm());
145 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
146
147 /*Initialize arrays:*/
148 maxmeans=xNew<IssmDouble*>(nfields);
149 minmeans=xNew<IssmDouble*>(nfields);
150 meanxtype=xNew<int>(nfields);
151 meanxsize=xNew<int>(nfields);
152
153 maxxs=xNew<IssmDouble*>(nfields*nsteps);
154 minxs=xNew<IssmDouble*>(nfields*nsteps);
155 xtype=xNew<int>(nfields*nsteps);
156 xsize=xNew<int>(nfields*nsteps);
157
158 /*Start opening files:*/
159 for(int i=(lower_row+1);i<=upper_row;i++){
160 _printf0_("reading file #: " << i << "\n");
161 char file[1000];
162 long int length;
163 char* buffer=NULL;
164
165 /*string:*/
166 sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
167
168 /*open file: */
169 _printf0_(" opening file: " << file << "\n");
170 FILE* fid=fopen(file,"rb");
171 if(fid==NULL)_error_("cound not open file: " << file << "\n");
172
173 /*figure out size of file, and read the whole thing:*/
174 _printf0_(" reading file:\n");
175 fseek (fid, 0, SEEK_END);
176 length = ftell (fid);
177 fseek (fid, 0, SEEK_SET);
178 buffer = xNew<char>(length);
179 fread (buffer, sizeof(char), length, fid);
180
181 /*close file:*/
182 fclose (fid);
183
184 /*create a memory stream with this buffer:*/
185 _printf0_(" processing file:\n");
186 fid=fmemopen(buffer, length, "rb");
187
188 /*start reading data from the buffer directly:*/
189 for (int f=0;f<nfields;f++){
190 char* field=fields[f];
191 fseek(fid,0,SEEK_SET);
192 for (int j=0;j<nsteps;j++){
193 int counter=f*nsteps+j;
194 xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
195 if(i==(lower_row+1)){
196 if(xtype[counter]==1){
197 maxxs[counter]=xNew<IssmDouble>(1);
198 minxs[counter]=xNew<IssmDouble>(1);
199 *maxxs[counter]=scalar;
200 *minxs[counter]=scalar;
201 xsize[counter]=1;
202 }
203 else if (xtype[counter]==3){
204 maxxs[counter]=xNew<IssmDouble>(doublematsize);
205 xMemCpy<IssmDouble>(maxxs[counter],doublemat,doublematsize);
206 minxs[counter]=xNew<IssmDouble>(doublematsize);
207 xMemCpy<IssmDouble>(minxs[counter],doublemat,doublematsize);
208 xsize[counter]=doublematsize;
209 xDelete<IssmDouble>(doublemat);
210 }
211 else _error_("cannot carry out statistics on type " << xtype[counter]);
212 }
213 else{
214 if(xtype[counter]==1){
215 *maxxs[counter]=max(*maxxs[counter],scalar);
216 *minxs[counter]=min(*minxs[counter],scalar);
217 }
218 else if (xtype[counter]==3){
219 IssmDouble* newmax=maxxs[counter];
220 IssmDouble* newmin=minxs[counter];
221 for(int k=0;k<doublematsize;k++){
222 if(doublemat[k]>newmax[k])newmax[k]=doublemat[k];
223 if(doublemat[k]<newmin[k])newmin[k]=doublemat[k];
224 }
225 xDelete<IssmDouble>(doublemat);
226 }
227 else _error_("cannot carry out statistics on type " << xtype[counter]);
228 }
229 }
230 }
231 _printf0_(" average in time:\n");
232
233 /*Deal with average in time: */
234 for (int f=0;f<nfields;f++){
235 fseek(fid,0,SEEK_SET);
236 char* field=fields[f];
237 meanxtype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
238
239 if(meanxtype[f]==1){
240 meanxsize[f]=1;
241 IssmDouble timemean=0;
242 fseek(fid,0,SEEK_SET);
243 for (int j=0;j<nsteps;j++){
244 readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
245 timemean+=scalar/nsteps;
246 }
247
248 /*Figure out max and min of time means: */
249 if(i==(lower_row+1)){
250 maxmeans[f]=xNewZeroInit<IssmDouble>(1);
251 minmeans[f]=xNewZeroInit<IssmDouble>(1);
252 *maxmeans[f]=timemean;
253 *minmeans[f]=timemean;
254 }
255 else{
256 *maxmeans[f]=max(*maxmeans[f],timemean);
257 *minmeans[f]=min(*minmeans[f],timemean);
258 }
259 }
260 else{
261 meanxsize[f]=doublematsize;
262 fseek(fid,0,SEEK_SET);
263 IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
264 for (int j=0;j<nsteps;j++){
265 readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
266 for (int k=0;k<doublematsize;k++){
267 timemean[k]+=doublemat[k]/nsteps;
268 }
269 xDelete<IssmDouble>(doublemat);
270 }
271
272 if(i==(lower_row+1)){
273 maxmeans[f]=xNew<IssmDouble>(doublematsize);
274 xMemCpy<IssmDouble>(maxmeans[f],timemean,doublematsize);
275 minmeans[f]=xNew<IssmDouble>(doublematsize);
276 xMemCpy<IssmDouble>(minmeans[f],timemean,doublematsize);
277 }
278 else{
279 IssmDouble* maxx=maxmeans[f];
280 IssmDouble* minx=minmeans[f];
281
282 for(int k=0;k<doublematsize;k++){
283 maxx[k]=max(maxx[k],timemean[k]);
284 minx[k]=min(minx[k],timemean[k]);
285 }
286 maxmeans[f]=maxx;
287 minmeans[f]=minx;
288 }
289 }
290 }
291 fclose(fid);
292
293 /*delete buffer:*/
294 xDelete<char>(buffer);
295 }
296 ISSM_MPI_Barrier(IssmComm::GetComm());
297 _printf0_("Done reading files, now computing min and max.\n");
298
299 /*We have agregated minx and max across the cluster, now gather across the cluster onto
300 *cpu0 and then compute statistics:*/
301 for (int f=0;f<nfields;f++){
302 int counter0=f*nsteps+0;
303 if (xtype[counter0]==1){ /*deal with scalars {{{*/
304 for (int j=0;j<nsteps;j++){
305 int counter=f*nsteps+j;
306
307 /*we are broadcasting doubles:*/
308 IssmDouble maxscalar=*maxxs[counter];
309 IssmDouble minscalar=*minxs[counter];
310 IssmDouble allmaxscalar;
311 IssmDouble allminscalar;
312 IssmDouble sumscalar_alltimes=0;
313
314 ISSM_MPI_Allreduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
315 ISSM_MPI_Allreduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
316
317 /*Store broadcasted value for later computation of histograms:*/
318 *maxxs[counter]=allmaxscalar;
319 *minxs[counter]=allminscalar;
320
321 }
322 } /*}}}*/
323 else{ /*deal with arrays:{{{*/
324
325 int size=xsize[counter0];
326 for (int j=0;j<nsteps;j++){
327 int counter=f*nsteps+j;
328
329 /*we are broadcasting double arrays:*/
330 IssmDouble* maxx=maxxs[counter];
331 IssmDouble* minx=minxs[counter];
332
333 IssmDouble* allmax=xNew<IssmDouble>(size);
334 IssmDouble* allmin=xNew<IssmDouble>(size);
335
336 ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
337 ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
338
339 /*Store broadcasted value for later computation of histograms:*/
340 maxxs[counter]=allmax;
341 minxs[counter]=allmin;
342 }
343 } /*}}}*/
344 }
345
346 /*Now do the same for the time mean fields:*/
347 for (int f=0;f<nfields;f++){
348 if (meanxtype[f]==1){ /*deal with scalars {{{*/
349
350 /*we are broadcasting doubles:*/
351 IssmDouble maxscalar=*maxmeans[f];
352 IssmDouble minscalar=*minmeans[f];
353 IssmDouble allmaxscalar;
354 IssmDouble allminscalar;
355
356 ISSM_MPI_Allreduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
357 ISSM_MPI_Allreduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
358
359 /*Store for later use in histogram computation:*/
360 *maxmeans[f]=allmaxscalar;
361 *minmeans[f]=allminscalar;
362
363 } /*}}}*/
364 else{ /*deal with arrays:{{{*/
365
366 int size=meanxsize[f];
367
368 /*we are broadcasting double arrays:*/
369 IssmDouble* maxx=maxmeans[f];
370 IssmDouble* minx=minmeans[f];
371
372 IssmDouble* allmax=xNew<IssmDouble>(size);
373 IssmDouble* allmin=xNew<IssmDouble>(size);
374
375 ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
376 ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm());
377
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 localhistogram[index]++;
436 histogram[counter]=localhistogram;
437 }
438 else if (xtype[counter]==3){
439 IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(doublematsize*nbins);
440 IssmDouble* ma=maxxs[counter];
441 IssmDouble* mi=minxs[counter];
442 for (int k=0;k<doublematsize;k++){
443 IssmDouble scalar=doublemat[k];
444 int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index--;
445 _assert_(scalar<=ma[k]); _assert_(scalar>=mi[k]); _assert_(index<nbins);
446 localhistogram[k*nbins+index]++;
447 }
448 histogram[counter]=localhistogram;
449 xDelete<IssmDouble>(doublemat);
450 }
451 else _error_("cannot carry out statistics on type " << xtype[counter]);
452 }
453 else{
454 if(xtype[counter]==1){
455 IssmDouble* localhistogram=histogram[counter];
456 IssmDouble ma=*maxxs[counter];
457 IssmDouble mi=*minxs[counter];
458 int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
459 localhistogram[index]++;
460 }
461 else if (xtype[counter]==3){
462 IssmDouble* localhistogram=histogram[counter];
463 IssmDouble* ma=maxxs[counter];
464 IssmDouble* mi=minxs[counter];
465 for (int k=0;k<doublematsize;k++){
466 IssmDouble scalar=doublemat[k];
467 int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
468 localhistogram[k*nbins+index]++;
469 }
470 xDelete<IssmDouble>(doublemat);
471 }
472 else _error_("cannot carry out statistics on type " << xtype[counter]);
473 }
474 }
475 }
476 _printf0_(" average in time:\n");
477
478 /*Deal with average in time: */
479 for (int f=0;f<nfields;f++){
480 fseek(fid,0,SEEK_SET);
481 char* field=fields[f];
482 meanxtype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
483
484 if(meanxtype[f]==1){
485 IssmDouble timemean=0;
486 fseek(fid,0,SEEK_SET);
487 for (int j=0;j<nsteps;j++){
488 readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
489 timemean+=scalar/nsteps;
490 }
491
492 /*Figure out max and min of time means: */
493 if(i==(lower_row+1)){
494 IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(nbins);
495 IssmDouble ma=*maxmeans[f];
496 IssmDouble mi=*minmeans[f];
497 int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
498 localhistogram[index]++;
499 timehistogram[f]=localhistogram;
500 }
501 else{
502 IssmDouble* localhistogram=timehistogram[f];
503 IssmDouble ma=*maxmeans[f];
504 IssmDouble mi=*minmeans[f];
505 int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1;
506 localhistogram[index]++;
507 }
508 }
509 else{
510 fseek(fid,0,SEEK_SET);
511 IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
512 for (int j=0;j<nsteps;j++){
513 readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
514 for (int k=0;k<doublematsize;k++){
515 timemean[k]+=doublemat[k]/nsteps;
516 }
517 xDelete<IssmDouble>(doublemat);
518 }
519
520 if(i==(lower_row+1)){
521 IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(doublematsize*nbins);
522 IssmDouble* ma=maxmeans[f];
523 IssmDouble* mi=minmeans[f];
524
525 for (int k=0;k<doublematsize;k++){
526 IssmDouble scalar=timemean[k];
527 int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
528 localhistogram[k*nbins+index]++;
529 }
530 timehistogram[f]=localhistogram;
531 }
532 else{
533
534 IssmDouble* localhistogram=timehistogram[f];
535 IssmDouble* ma=maxmeans[f];
536 IssmDouble* mi=minmeans[f];
537
538 for (int k=0;k<doublematsize;k++){
539 IssmDouble scalar=timemean[k];
540 int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1;
541 localhistogram[k*nbins+index]++;
542 }
543 }
544 }
545 }
546 fclose(fid);
547
548 /*delete buffer:*/
549 xDelete<char>(buffer);
550 }
551 _printf0_("Start aggregating histogram:\n");
552
553 /*We have agregated histograms across the cluster, now gather them across the cluster onto
554 *cpu0: */
555 for (int f=0;f<nfields;f++){
556 int counter0=f*nsteps+0;
557 if (xtype[counter0]==1){ /*deal with scalars {{{*/
558 for (int j=0;j<nsteps;j++){
559 int counter=f*nsteps+j;
560
561 /*we are broadcasting doubles:*/
562 IssmDouble* histo=histogram[counter]; //size nbins
563 IssmDouble* allhisto=xNewZeroInit<IssmDouble>(nbins);
564
565 ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
566
567 /*add to results:*/
568 if(my_rank==0){
569 char fieldname[1000];
570
571 sprintf(fieldname,"%s%s",fields[f],"Histogram");
572 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,1,nbins,j+1,0));
573
574 sprintf(fieldname,"%s%s",fields[f],"Max");
575 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*maxxs[counter],j+1,0));
576 sprintf(fieldname,"%s%s",fields[f],"Min");
577 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*minxs[counter],j+1,0));
578 }
579 }
580 } /*}}}*/
581 else{ /*deal with arrays:{{{*/
582
583 int size=xsize[counter0];
584 for (int j=0;j<nsteps;j++){
585 int counter=f*nsteps+j;
586
587 /*we are broadcasting double arrays:*/
588 IssmDouble* histo=histogram[counter];
589 IssmDouble* allhisto=xNew<IssmDouble>(size*nbins);
590
591 ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
592 xDelete<IssmDouble>(histo);
593
594 /*add to results:*/
595 if(my_rank==0){
596 char fieldname[1000];
597
598 sprintf(fieldname,"%s%s",fields[f],"Histogram");
599 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,size,nbins,j+1,0));
600
601 sprintf(fieldname,"%s%s",fields[f],"Max");
602 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,maxxs[counter],size,1,j+1,0));
603 sprintf(fieldname,"%s%s",fields[f],"Min");
604 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,minxs[counter],size,1,j+1,0));
605 }
606 }
607 } /*}}}*/
608 }
609 _printf0_("Start aggregating time mean histogram:\n");
610
611 /*Now do the same for the time mean fields:*/
612 for (int f=0;f<nfields;f++){
613 if (meanxtype[f]==1){ /*deal with scalars {{{*/
614
615 /*we are broadcasting doubles:*/
616 IssmDouble* histo=timehistogram[f];
617 IssmDouble* allhisto=xNewZeroInit<IssmDouble>(nbins);
618
619 ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
620
621 /*add to results at time step 1:*/
622 if(my_rank==0){
623 char fieldname[1000];
624
625 sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram");
626 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,1,nbins,1,0));
627
628 sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax");
629 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*maxmeans[f],1,0));
630 sprintf(fieldname,"%s%s",fields[f],"TimeMeaMin");
631 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,*minmeans[f],1,0));
632 }
633 } /*}}}*/
634 else{ /*deal with arrays:{{{*/
635
636 int size=meanxsize[f];
637
638 /*we are broadcasting double arrays:*/
639 IssmDouble* histo=timehistogram[f];
640 IssmDouble* allhisto=xNewZeroInit<IssmDouble>(size*nbins);
641
642 ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
643 xDelete<IssmDouble>(histo);
644 /*add to results at step 1:*/
645 if(my_rank==0){
646 char fieldname[1000];
647
648 sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram");
649 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,size,nbins,1,0));
650 sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax");
651 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,maxmeans[f],size,1,1,0));
652 sprintf(fieldname,"%s%s",fields[f],"TimeMeanMin");
653 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,minmeans[f],size,1,1,0));
654 }
655 } /*}}}*/
656 }
657 _printf0_("Done aggregating time mean histogram:\n");
658 IssmComm::SetComm(ISSM_MPI_COMM_WORLD);
659}
660/*}}}*/
661int ComputeMeanVariance(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){ /*{{{*/
662
663 int nsamples;
664 char* directory=NULL;
665 char* model=NULL;
666 char** fields=NULL;
667 int* steps=NULL;
668 int nsteps;
669 int nfields;
670 int range,lower_row,upper_row;
671 int nfilesperdirectory;
672
673 /*intermediary:*/
674 IssmDouble* doublemat=NULL;
675 int doublematsize;
676 IssmDouble scalar;
677
678 /*computation of average and variance itself:*/
679 IssmDouble* x = NULL;
680 IssmDouble* x2 = NULL;
681 IssmDouble** xs = NULL;
682 IssmDouble** xs2 = NULL;
683 int* xtype=NULL;
684 int* xsize=NULL;
685
686 IssmDouble** meanx=NULL;
687 IssmDouble** meanx2=NULL;
688 int* meantype=NULL;
689 int* meansize=NULL;
690
691 /*only work on the statistical communicator: */
692 if (color==MPI_UNDEFINED)return 0;
693
694 /*Retrieve parameters:*/
695 parameters->FindParam(&nfilesperdirectory,QmuNfilesPerDirectoryEnum);
696 parameters->FindParam(&nsamples,QmuNsampleEnum);
697 parameters->FindParam(&directory,DirectoryNameEnum);
698 parameters->FindParam(&model,InputFileNameEnum);
699 parameters->FindParam(&fields,&nfields,FieldsEnum);
700 parameters->FindParam(&steps,&nsteps,StepsEnum);
701
702 /*Get rank from the stat comm communicator:*/
703 IssmComm::SetComm(statcomm);
704 int my_rank=IssmComm::GetRank();
705
706 /*Open files and read them complelety, in a distributed way:*/
707 range=DetermineLocalSize(nsamples,IssmComm::GetComm());
708 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
709
710 /*Initialize arrays:*/
711 xs=xNew<IssmDouble*>(nfields*nsteps);
712 xs2=xNew<IssmDouble*>(nfields*nsteps);
713 xtype=xNew<int>(nfields*nsteps);
714 xsize=xNew<int>(nfields*nsteps);
715
716 meantype=xNew<int>(nfields);
717 meansize=xNew<int>(nfields);
718 meanx=xNew<IssmDouble*>(nfields);
719 meanx2=xNew<IssmDouble*>(nfields);
720
721 /*Start opening files:*/
722 for (int i=(lower_row+1);i<=upper_row;i++){
723 _printf0_("reading file #: " << i << "\n");
724 char file[1000];
725 long int length;
726 char* buffer=NULL;
727
728 /*string:*/
729 sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
730
731 /*open file: */
732 _printf0_(" opening file: " << file << "\n");
733 FILE* fid=fopen(file,"rb");
734 if(fid==NULL) _error_(" could not open file: " << file << "\n");
735
736 /*figure out size of file, and read the whole thing:*/
737 _printf0_(" reading file:\n");
738 fseek (fid, 0, SEEK_END);
739 length = ftell (fid);
740 fseek (fid, 0, SEEK_SET);
741 buffer = xNew<char>(length);
742 fread (buffer, sizeof(char), length, fid);
743
744 /*close file:*/
745 fclose (fid);
746
747 /*create a memory stream with this buffer:*/
748 _printf0_(" processing file:\n");
749 fid=fmemopen(buffer, length, "rb");
750
751 /*start reading data from the buffer directly:*/
752 for (int f=0;f<nfields;f++){
753 char* field=fields[f];
754 fseek(fid,0,SEEK_SET);
755 for (int j=0;j<nsteps;j++){
756 int counter=f*nsteps+j;
757 xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
758 if(i==(lower_row+1)){
759 if(xtype[counter]==1){
760 xs[counter]=xNew<IssmDouble>(1);
761 xs2[counter]=xNew<IssmDouble>(1);
762 *xs[counter]=scalar;
763 *xs2[counter]=pow(scalar,2.0);
764 xsize[counter]=1;
765 }
766 else if (xtype[counter]==3){
767 IssmDouble* doublemat2=xNew<IssmDouble>(doublematsize);
768 for(int k=0;k<doublematsize;k++)doublemat2[k]=pow(doublemat[k],2.0);
769 xs[counter]=doublemat;
770 xs2[counter]=doublemat2;
771 xsize[counter]=doublematsize;
772 }
773 else _error_("cannot carry out statistics on type " << xtype[counter]);
774 }
775 else{
776 if(xtype[counter]==1){
777 *xs[counter]+=scalar;
778 *xs2[counter]+=pow(scalar,2.0);
779 }
780 else if (xtype[counter]==3){
781 IssmDouble* newdoublemat=xs[counter];
782 IssmDouble* newdoublemat2=xs2[counter];
783 for(int k=0;k<doublematsize;k++){
784 newdoublemat[k]+=doublemat[k];
785 newdoublemat2[k]+=pow(doublemat[k],2.0);
786 }
787 xs[counter]=newdoublemat;
788 xs2[counter]=newdoublemat2;
789 }
790 else _error_("cannot carry out statistics on type " << xtype[counter]);
791 }
792 }
793 }
794
795 /*Deal with time mean: */
796 for (int f=0;f<nfields;f++){
797 char* field=fields[f];
798 fseek(fid,0,SEEK_SET);
799 meantype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
800 if(i==(lower_row+1)){
801 if(meantype[f]==1){
802 meanx[f]=xNewZeroInit<IssmDouble>(1);
803 meanx2[f]=xNewZeroInit<IssmDouble>(1);
804 meansize[f]=1;
805 }
806 else{
807 meanx[f]=xNewZeroInit<IssmDouble>(doublematsize);
808 meanx2[f]=xNewZeroInit<IssmDouble>(doublematsize);
809 meansize[f]=doublematsize;
810 }
811 }
812 fseek(fid,0,SEEK_SET);
813 if(meantype[f]==1){
814 IssmDouble sc=0;
815 IssmDouble sc2=0;
816 for(int j=0;j<nsteps;j++){
817 readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
818 sc+=scalar/nsteps;
819 }
820 sc2+=pow(sc,2.0);
821 *meanx[f]+=sc;
822 *meanx2[f]+=sc2;
823 }
824 else{
825 IssmDouble* sc=meanx[f];
826 IssmDouble* sc2=meanx2[f];
827 IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
828 IssmDouble* timemean2=xNewZeroInit<IssmDouble>(doublematsize);
829
830 for(int j=0;j<nsteps;j++){
831 readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
832 for (int k=0;k<doublematsize;k++){
833 timemean[k]+=doublemat[k]/nsteps;
834 }
835 }
836 for (int k=0;k<doublematsize;k++){
837 timemean2[k]=pow(timemean[k],2.0);
838 }
839 for (int k=0;k<doublematsize;k++){
840 sc[k]+=timemean[k];
841 sc2[k]+=timemean2[k];
842 }
843
844 }
845
846 }
847 fclose(fid);
848
849 /*delete buffer:*/
850 xDelete<char>(buffer);
851 }
852 ISSM_MPI_Barrier(IssmComm::GetComm());
853 _printf0_("Done reading files, now computing mean and variance.\n");
854
855 /*We have agregated x and x^2 across the cluster, now gather across the cluster onto
856 *cpu0 and then compute statistics:*/
857 for (int f=0;f<nfields;f++){
858 int counter0=f*nsteps+0;
859 if (xtype[counter0]==1){ /*deal with scalars {{{*/
860 IssmDouble mean,stddev;
861 for (int j=0;j<nsteps;j++){
862 int counter=f*nsteps+j;
863
864 /*we are broadcasting doubles:*/
865 IssmDouble scalar=*xs[counter];
866 IssmDouble scalar2=*xs2[counter];
867 IssmDouble sumscalar;
868 IssmDouble sumscalar2;
869
870 ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
871 ISSM_MPI_Reduce(&scalar2,&sumscalar2,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
872 /*Build average and standard deviation. For standard deviation, use the
873 *following formula: sigma^2=E(x^2)-mu^2:*/
874 mean=sumscalar/nsamples;
875 stddev=sqrt(sumscalar2/nsamples-pow(mean,2.0));
876
877 /*add to results:*/
878 if(my_rank==0){
879 char fieldname[1000];
880
881 sprintf(fieldname,"%s%s",fields[f],"Mean");
882 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,j+1,0));
883 sprintf(fieldname,"%s%s",fields[f],"Stddev");
884 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,j+1,0));
885 }
886
887 }
888 } /*}}}*/
889 else{ /*deal with arrays:{{{*/
890
891 int size=xsize[counter0];
892
893 IssmDouble* mean=xNew<IssmDouble>(size);
894 IssmDouble* stddev=xNew<IssmDouble>(size);
895
896 for (int j=0;j<nsteps;j++){
897 int counter=f*nsteps+j;
898
899 /*we are broadcasting double arrays:*/
900 x=xs[counter];
901 x2=xs2[counter];
902
903 IssmDouble* sumx=xNew<IssmDouble>(size);
904 IssmDouble* sumx2=xNew<IssmDouble>(size);
905
906 ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
907 ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
908
909 /*Build average and standard deviation. For standard deviation, use the
910 *following formula: sigma^2=E(x^2)-mu^2:*/
911 for (int k=0;k<size;k++){
912 mean[k]=sumx[k]/nsamples;
913 stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
914 }
915
916 /*add to results:*/
917 if(my_rank==0){
918 char fieldname[1000];
919
920 sprintf(fieldname,"%s%s",fields[f],"Mean");
921 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,j+1,0));
922 sprintf(fieldname,"%s%s",fields[f],"Stddev");
923 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,j+1,0));
924 }
925 }
926 } /*}}}*/
927 }
928 /*Do the same but for the time mean:*/
929 for (int f=0;f<nfields;f++){
930 if (meantype[f]==1){ /*deal with scalars {{{*/
931 IssmDouble mean,stddev;
932
933 /*we are broadcasting doubles:*/
934 IssmDouble scalar=*meanx[f];
935 IssmDouble scalar2=*meanx2[f];
936 IssmDouble sumscalar;
937 IssmDouble sumscalar2;
938
939 ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
940 ISSM_MPI_Reduce(&scalar2,&sumscalar2,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
941 /*Build average and standard deviation. For standard deviation, use the
942 *following formula: sigma^2=E(x^2)-mu^2:*/
943 mean=sumscalar/nsamples;
944 stddev=sqrt(sumscalar2/nsamples-pow(mean,2.0));
945
946 /*add to results:*/
947 if(my_rank==0){
948 char fieldname[1000];
949
950 sprintf(fieldname,"%s%s",fields[f],"TimeMean");
951 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,1,0));
952 sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
953 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,1,0));
954 }
955 } /*}}}*/
956 else{ /*deal with arrays:{{{*/
957
958 int size=meansize[f];
959 IssmDouble* mean=xNew<IssmDouble>(size);
960 IssmDouble* stddev=xNew<IssmDouble>(size);
961
962 /*we are broadcasting double arrays:*/
963 x=meanx[f];
964 x2=meanx2[f];
965
966 IssmDouble* sumx=xNew<IssmDouble>(size);
967 IssmDouble* sumx2=xNew<IssmDouble>(size);
968
969 ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
970 ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
971
972 /*Build average and standard deviation. For standard deviation, use the
973 *following formula: sigma^2=E(x^2)-mu^2:*/
974 for (int k=0;k<size;k++){
975 mean[k]=sumx[k]/nsamples;
976 stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
977 }
978
979 /*add to results:*/
980 if(my_rank==0){
981 char fieldname[1000];
982
983 sprintf(fieldname,"%s%s",fields[f],"TimeMean");
984 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,1,0));
985 sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
986 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,1,0));
987 }
988 } /*}}}*/
989 }
990
991 _printf0_("Done with MeanVariance :\n");
992 IssmComm::SetComm(ISSM_MPI_COMM_WORLD);
993
994} /*}}}*/
995int ComputeSampleSeries(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){ /*{{{*/
996
997 int nsamples;
998 char* directory=NULL;
999 char* model=NULL;
1000 char** fields=NULL;
1001 int* steps=NULL;
1002 int nsteps;
1003 int nfields;
1004 int range,lower_row,upper_row;
1005 int nfilesperdirectory;
1006 int* indices=NULL;
1007 int nindices;
1008
1009 /*intermediary:*/
1010 IssmDouble* doublemat=NULL;
1011 int doublematsize;
1012 IssmDouble scalar;
1013
1014 /*computation of average and variance itself:*/
1015 IssmDouble* x = NULL;
1016 IssmDouble* allx=NULL;
1017 IssmDouble** xs = NULL;
1018 int* xtype=NULL;
1019 int* xsize=NULL;
1020
1021 /*only work on the statistical communicator: */
1022 if (color==MPI_UNDEFINED)return 0;
1023
1024 /*Retrieve parameters:*/
1025 parameters->FindParam(&nsamples,QmuNsampleEnum);
1026 parameters->FindParam(&directory,DirectoryNameEnum);
1027 parameters->FindParam(&model,InputFileNameEnum);
1028 parameters->FindParam(&fields,&nfields,FieldsEnum);
1029 parameters->FindParam(&steps,&nsteps,StepsEnum);
1030 parameters->FindParam(&indices,&nindices,IndicesEnum);
1031
1032 /*Get rank from the stat comm communicator:*/
1033 IssmComm::SetComm(statcomm);
1034 int my_rank=IssmComm::GetRank();
1035
1036 /*Open files and read them complelety, in a distributed way:*/
1037 range=DetermineLocalSize(nsamples,IssmComm::GetComm());
1038 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
1039
1040 /*Initialize arrays:*/
1041 xs=xNew<IssmDouble*>(nfields*nsteps);
1042 xtype=xNew<int>(nfields*nsteps);
1043 xsize=xNew<int>(nfields*nsteps);
1044
1045 /*Start opening files:*/
1046 for (int i=(lower_row+1);i<=upper_row;i++){
1047 _printf0_("reading file #: " << i << "\n");
1048 char file[1000];
1049 long int length;
1050 char* buffer=NULL;
1051
1052 /*string:*/
1053 sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
1054
1055 /*open file: */
1056 _printf0_(" opening file:\n");
1057 FILE* fid=fopen(file,"rb");
1058
1059 /*figure out size of file, and read the whole thing:*/
1060 _printf0_(" reading file:\n");
1061 fseek (fid, 0, SEEK_END);
1062 length = ftell (fid);
1063 fseek (fid, 0, SEEK_SET);
1064 buffer = xNew<char>(length);
1065 fread (buffer, sizeof(char), length, fid);
1066
1067 /*close file:*/
1068 fclose (fid);
1069
1070 /*create a memory stream with this buffer:*/
1071 _printf0_(" processing file:\n");
1072 fid=fmemopen(buffer, length, "rb");
1073
1074 /*start reading data from the buffer directly:*/
1075 for (int f=0;f<nfields;f++){
1076 fseek(fid,0,SEEK_SET);
1077 char* field=fields[f];
1078 for (int j=0;j<nsteps;j++){
1079 int counter=f*nsteps+j;
1080 xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
1081 if(i==(lower_row+1)){
1082 if(xtype[counter]==1){
1083 x=xNew<IssmDouble>(range);
1084 x[0]=scalar;
1085 xs[counter]=x;
1086 xsize[counter]=range;
1087 }
1088 else if (xtype[counter]==3){
1089 x=xNew<IssmDouble>(nindices*range);
1090 for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
1091 xs[counter]=x;
1092 xsize[counter]=range*nindices;
1093 }
1094 else _error_("cannot carry out statistics on type " << xtype[counter]);
1095 }
1096 else{
1097 if(xtype[counter]==1){
1098 x=xs[counter];
1099 x[i-(lower_row+1)]=scalar;
1100 xs[counter]=x;
1101 }
1102 else if (xtype[counter]==3){
1103 x=xs[counter];
1104 for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
1105 xs[counter]=x;
1106 }
1107 else _error_("cannot carry out statistics on type " << xtype[counter]);
1108 }
1109 }
1110 }
1111 fclose(fid);
1112
1113 /*delete buffer:*/
1114 xDelete<char>(buffer);
1115 }
1116 ISSM_MPI_Barrier(IssmComm::GetComm());
1117 _printf0_("Done reading files, now assembling time series.\n");
1118
1119 for (int f=0;f<nfields;f++){
1120 for (int j=0;j<nsteps;j++){
1121 int counter=f*nsteps+j;
1122 if (xtype[counter]==1){
1123 /*we are broadcasting range times doubles:*/
1124 x=xs[counter];
1125 allx=xNew<IssmDouble>(nsamples);
1126 MPI_Gather(x, range, ISSM_MPI_PDOUBLE,allx, range, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
1127 /*add to results:*/
1128 if(my_rank==0){
1129 char fieldname[1000];
1130
1131 sprintf(fieldname,"%s%s",fields[f],"Samples");
1132 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,1,j+1,0));
1133 }
1134 }
1135 else{
1136 /*we are broadcasting double arrays:*/
1137 x=xs[counter];
1138 allx=xNew<IssmDouble>(nsamples*nindices);
1139
1140 MPI_Gather(x, range*nindices, ISSM_MPI_PDOUBLE,allx, range*nindices, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
1141
1142 /*add to results:*/
1143 if(my_rank==0){
1144 char fieldname[1000];
1145 sprintf(fieldname,"%s%s",fields[f],"Samples");
1146 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,nindices,j+1,0));
1147 }
1148 }
1149 }
1150 }
1151 _printf0_("Done with SampleSeries :\n");
1152 IssmComm::SetComm(ISSM_MPI_COMM_WORLD);
1153
1154} /*}}}*/
1155int OutputStatistics(Parameters* parameters,Results* results,int color,ISSM_MPI_Comm statcomm){ /*{{{*/
1156
1157 char outputfilename[1000];
1158 char* directory=NULL;
1159 char* model=NULL;
1160 char* method=NULL;
1161 int nsamples;
1162 int* steps=NULL;
1163 int nsteps;
1164
1165 /*only work on the statistical communicator: */
1166 if (color==MPI_UNDEFINED)return 0;
1167
1168 FemModel* femmodel=new FemModel();
1169
1170 /*Some parameters that will allow us to use the OutputResultsx module:*/
1171 parameters->AddObject(new BoolParam(QmuIsdakotaEnum,false));
1172 parameters->AddObject(new BoolParam(SettingsIoGatherEnum,true));
1173
1174 parameters->FindParam(&directory,DirectoryNameEnum);
1175 parameters->FindParam(&model,InputFileNameEnum);
1176 parameters->FindParam(&nsamples,QmuNsampleEnum);
1177 parameters->FindParam(&steps,&nsteps,StepsEnum);
1178
1179 sprintf(outputfilename,"%s/%s.stats",directory,model);
1180 parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
1181
1182 /*Call OutputResults module:*/
1183 femmodel->parameters=parameters;
1184 femmodel->results=results;
1185
1186 OutputResultsx(femmodel);
1187} /*}}}*/
Note: See TracBrowser for help on using the repository browser.