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

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

CHG: diverse
File size: 35.5 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
672 /*intermediary:*/
673 IssmDouble* doublemat=NULL;
674 int doublematsize;
675 IssmDouble scalar;
676
677 /*computation of average and variance itself:*/
678 IssmDouble* x = NULL;
679 IssmDouble* x2 = NULL;
680 IssmDouble** xs = NULL;
681 IssmDouble** xs2 = NULL;
682 int* xtype=NULL;
683 int* xsize=NULL;
684
685 IssmDouble** meanx=NULL;
686 IssmDouble** meanx2=NULL;
687 int* meantype=NULL;
688 int* meansize=NULL;
689
690 /*Retrieve parameters:*/
691 parameters->FindParam(&nsamples,QmuNsampleEnum);
692 parameters->FindParam(&directory,DirectoryNameEnum);
693 parameters->FindParam(&model,InputFileNameEnum);
694 parameters->FindParam(&fields,&nfields,FieldsEnum);
695 parameters->FindParam(&steps,&nsteps,StepsEnum);
696
697 /*Get rank:*/
698 int my_rank=IssmComm::GetRank();
699
700 /*Open files and read them complelety, in a distributed way:*/
701 range=DetermineLocalSize(nsamples,IssmComm::GetComm());
702 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
703
704 /*Initialize arrays:*/
705 xs=xNew<IssmDouble*>(nfields*nsteps);
706 xs2=xNew<IssmDouble*>(nfields*nsteps);
707 xtype=xNew<int>(nfields*nsteps);
708 xsize=xNew<int>(nfields*nsteps);
709
710 meantype=xNew<int>(nfields);
711 meansize=xNew<int>(nfields);
712 meanx=xNew<IssmDouble*>(nfields);
713 meanx2=xNew<IssmDouble*>(nfields);
714
715 /*Start opening files:*/
716 for (int i=(lower_row+1);i<=upper_row;i++){
717 _printf0_("reading file #: " << i << "\n");
718 char file[1000];
719 long int length;
720 char* buffer=NULL;
721
722 /*string:*/
723 sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
724
725 /*open file: */
726 _printf0_(" opening file: " << file << "\n");
727 FILE* fid=fopen(file,"rb");
728 if(fid==NULL) _error_(" could not open file: " << file << "\n");
729
730 /*figure out size of file, and read the whole thing:*/
731 _printf0_(" reading file:\n");
732 fseek (fid, 0, SEEK_END);
733 length = ftell (fid);
734 fseek (fid, 0, SEEK_SET);
735 buffer = xNew<char>(length);
736 fread (buffer, sizeof(char), length, fid);
737
738 /*close file:*/
739 fclose (fid);
740
741 /*create a memory stream with this buffer:*/
742 _printf0_(" processing file:\n");
743 fid=fmemopen(buffer, length, "rb");
744
745 /*start reading data from the buffer directly:*/
746 for (int f=0;f<nfields;f++){
747 char* field=fields[f];
748 fseek(fid,0,SEEK_SET);
749 for (int j=0;j<nsteps;j++){
750 int counter=f*nsteps+j;
751 xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
752 if(i==(lower_row+1)){
753 if(xtype[counter]==1){
754 xs[counter]=xNew<IssmDouble>(1);
755 xs2[counter]=xNew<IssmDouble>(1);
756 *xs[counter]=scalar;
757 *xs2[counter]=pow(scalar,2.0);
758 xsize[counter]=1;
759 }
760 else if (xtype[counter]==3){
761 IssmDouble* doublemat2=xNew<IssmDouble>(doublematsize);
762 for(int k=0;k<doublematsize;k++)doublemat2[k]=pow(doublemat[k],2.0);
763 xs[counter]=doublemat;
764 xs2[counter]=doublemat2;
765 xsize[counter]=doublematsize;
766 }
767 else _error_("cannot carry out statistics on type " << xtype[counter]);
768 }
769 else{
770 if(xtype[counter]==1){
771 *xs[counter]+=scalar;
772 *xs2[counter]+=pow(scalar,2.0);
773 }
774 else if (xtype[counter]==3){
775 IssmDouble* newdoublemat=xs[counter];
776 IssmDouble* newdoublemat2=xs2[counter];
777 for(int k=0;k<doublematsize;k++){
778 newdoublemat[k]+=doublemat[k];
779 newdoublemat2[k]+=pow(doublemat[k],2.0);
780 }
781 xs[counter]=newdoublemat;
782 xs2[counter]=newdoublemat2;
783 }
784 else _error_("cannot carry out statistics on type " << xtype[counter]);
785 }
786 }
787 }
788
789 /*Deal with time mean: */
790 for (int f=0;f<nfields;f++){
791 char* field=fields[f];
792 fseek(fid,0,SEEK_SET);
793 meantype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
794 if(i==(lower_row+1)){
795 if(meantype[f]==1){
796 meanx[f]=xNewZeroInit<IssmDouble>(1);
797 meanx2[f]=xNewZeroInit<IssmDouble>(1);
798 meansize[f]=1;
799 }
800 else{
801 meanx[f]=xNewZeroInit<IssmDouble>(doublematsize);
802 meanx2[f]=xNewZeroInit<IssmDouble>(doublematsize);
803 meansize[f]=doublematsize;
804 }
805 }
806 fseek(fid,0,SEEK_SET);
807 if(meantype[f]==1){
808 IssmDouble sc=0;
809 IssmDouble sc2=0;
810 for(int j=0;j<nsteps;j++){
811 readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
812 sc+=scalar/nsteps;
813 }
814 sc2+=pow(sc,2.0);
815 *meanx[f]+=sc;
816 *meanx2[f]+=sc2;
817 }
818 else{
819 IssmDouble* sc=meanx[f];
820 IssmDouble* sc2=meanx2[f];
821 IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
822 IssmDouble* timemean2=xNewZeroInit<IssmDouble>(doublematsize);
823
824 for(int j=0;j<nsteps;j++){
825 readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
826 for (int k=0;k<doublematsize;k++){
827 timemean[k]+=doublemat[k]/nsteps;
828 }
829 }
830 for (int k=0;k<doublematsize;k++){
831 timemean2[k]=pow(timemean[k],2.0);
832 }
833 for (int k=0;k<doublematsize;k++){
834 sc[k]+=timemean[k];
835 sc2[k]+=timemean2[k];
836 }
837
838 }
839
840 }
841 fclose(fid);
842
843 /*delete buffer:*/
844 xDelete<char>(buffer);
845 }
846 ISSM_MPI_Barrier(IssmComm::GetComm());
847 _printf0_("Done reading files, now computing mean and variance.\n");
848
849 /*We have agregated x and x^2 across the cluster, now gather across the cluster onto
850 *cpu0 and then compute statistics:*/
851 for (int f=0;f<nfields;f++){
852 int counter0=f*nsteps+0;
853 if (xtype[counter0]==1){ /*deal with scalars {{{*/
854 IssmDouble mean,stddev;
855 for (int j=0;j<nsteps;j++){
856 int counter=f*nsteps+j;
857
858 /*we are broadcasting doubles:*/
859 IssmDouble scalar=*xs[counter];
860 IssmDouble scalar2=*xs2[counter];
861 IssmDouble sumscalar;
862 IssmDouble sumscalar2;
863
864 ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
865 ISSM_MPI_Reduce(&scalar2,&sumscalar2,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
866 /*Build average and standard deviation. For standard deviation, use the
867 *following formula: sigma^2=E(x^2)-mu^2:*/
868 mean=sumscalar/nsamples;
869 stddev=sqrt(sumscalar2/nsamples-pow(mean,2.0));
870
871 /*add to results:*/
872 if(my_rank==0){
873 char fieldname[1000];
874
875 sprintf(fieldname,"%s%s",fields[f],"Mean");
876 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,j+1,0));
877 sprintf(fieldname,"%s%s",fields[f],"Stddev");
878 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,j+1,0));
879 }
880
881 }
882 } /*}}}*/
883 else{ /*deal with arrays:{{{*/
884
885 int size=xsize[counter0];
886
887 IssmDouble* mean=xNew<IssmDouble>(size);
888 IssmDouble* stddev=xNew<IssmDouble>(size);
889
890 for (int j=0;j<nsteps;j++){
891 int counter=f*nsteps+j;
892
893 /*we are broadcasting double arrays:*/
894 x=xs[counter];
895 x2=xs2[counter];
896
897 IssmDouble* sumx=xNew<IssmDouble>(size);
898 IssmDouble* sumx2=xNew<IssmDouble>(size);
899
900 ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
901 ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
902
903 /*Build average and standard deviation. For standard deviation, use the
904 *following formula: sigma^2=E(x^2)-mu^2:*/
905 for (int k=0;k<size;k++){
906 mean[k]=sumx[k]/nsamples;
907 stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
908 }
909
910 /*add to results:*/
911 if(my_rank==0){
912 char fieldname[1000];
913
914 sprintf(fieldname,"%s%s",fields[f],"Mean");
915 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,j+1,0));
916 sprintf(fieldname,"%s%s",fields[f],"Stddev");
917 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,j+1,0));
918 }
919 }
920 } /*}}}*/
921 }
922 /*Do the same but for the time mean:*/
923 for (int f=0;f<nfields;f++){
924 if (meantype[f]==1){ /*deal with scalars {{{*/
925 IssmDouble mean,stddev;
926
927 /*we are broadcasting doubles:*/
928 IssmDouble scalar=*meanx[f];
929 IssmDouble scalar2=*meanx2[f];
930 IssmDouble sumscalar;
931 IssmDouble sumscalar2;
932
933 ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
934 ISSM_MPI_Reduce(&scalar2,&sumscalar2,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
935 /*Build average and standard deviation. For standard deviation, use the
936 *following formula: sigma^2=E(x^2)-mu^2:*/
937 mean=sumscalar/nsamples;
938 stddev=sqrt(sumscalar2/nsamples-pow(mean,2.0));
939
940 /*add to results:*/
941 if(my_rank==0){
942 char fieldname[1000];
943
944 sprintf(fieldname,"%s%s",fields[f],"TimeMean");
945 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,1,0));
946 sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
947 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,1,0));
948 }
949 } /*}}}*/
950 else{ /*deal with arrays:{{{*/
951
952 int size=meansize[f];
953 IssmDouble* mean=xNew<IssmDouble>(size);
954 IssmDouble* stddev=xNew<IssmDouble>(size);
955
956 /*we are broadcasting double arrays:*/
957 x=meanx[f];
958 x2=meanx2[f];
959
960 IssmDouble* sumx=xNew<IssmDouble>(size);
961 IssmDouble* sumx2=xNew<IssmDouble>(size);
962
963 ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
964 ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
965
966 /*Build average and standard deviation. For standard deviation, use the
967 *following formula: sigma^2=E(x^2)-mu^2:*/
968 for (int k=0;k<size;k++){
969 mean[k]=sumx[k]/nsamples;
970 stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
971 }
972
973 /*add to results:*/
974 if(my_rank==0){
975 char fieldname[1000];
976
977 sprintf(fieldname,"%s%s",fields[f],"TimeMean");
978 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,1,0));
979 sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
980 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,1,0));
981 }
982 } /*}}}*/
983 }
984
985
986} /*}}}*/
987int ComputeSampleSeries(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){ /*{{{*/
988
989 int nsamples;
990 char* directory=NULL;
991 char* model=NULL;
992 char** fields=NULL;
993 int* steps=NULL;
994 int nsteps;
995 int nfields;
996 int range,lower_row,upper_row;
997 int* indices=NULL;
998 int nindices;
999
1000 /*intermediary:*/
1001 IssmDouble* doublemat=NULL;
1002 int doublematsize;
1003 IssmDouble scalar;
1004
1005 /*computation of average and variance itself:*/
1006 IssmDouble* x = NULL;
1007 IssmDouble* allx=NULL;
1008 IssmDouble** xs = NULL;
1009 int* xtype=NULL;
1010 int* xsize=NULL;
1011
1012 /*Retrieve parameters:*/
1013 parameters->FindParam(&nsamples,QmuNsampleEnum);
1014 parameters->FindParam(&directory,DirectoryNameEnum);
1015 parameters->FindParam(&model,InputFileNameEnum);
1016 parameters->FindParam(&fields,&nfields,FieldsEnum);
1017 parameters->FindParam(&steps,&nsteps,StepsEnum);
1018 parameters->FindParam(&indices,&nindices,IndicesEnum);
1019
1020 /*Get rank:*/
1021 int my_rank=IssmComm::GetRank();
1022
1023 /*Open files and read them complelety, in a distributed way:*/
1024 range=DetermineLocalSize(nsamples,IssmComm::GetComm());
1025 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
1026
1027 /*Initialize arrays:*/
1028 xs=xNew<IssmDouble*>(nfields*nsteps);
1029 xtype=xNew<int>(nfields*nsteps);
1030 xsize=xNew<int>(nfields*nsteps);
1031
1032 /*Start opening files:*/
1033 for (int i=(lower_row+1);i<=upper_row;i++){
1034 _printf0_("reading file #: " << i << "\n");
1035 char file[1000];
1036 long int length;
1037 char* buffer=NULL;
1038
1039 /*string:*/
1040 sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
1041
1042 /*open file: */
1043 _printf0_(" opening file:\n");
1044 FILE* fid=fopen(file,"rb");
1045
1046 /*figure out size of file, and read the whole thing:*/
1047 _printf0_(" reading file:\n");
1048 fseek (fid, 0, SEEK_END);
1049 length = ftell (fid);
1050 fseek (fid, 0, SEEK_SET);
1051 buffer = xNew<char>(length);
1052 fread (buffer, sizeof(char), length, fid);
1053
1054 /*close file:*/
1055 fclose (fid);
1056
1057 /*create a memory stream with this buffer:*/
1058 _printf0_(" processing file:\n");
1059 fid=fmemopen(buffer, length, "rb");
1060
1061 /*start reading data from the buffer directly:*/
1062 for (int f=0;f<nfields;f++){
1063 fseek(fid,0,SEEK_SET);
1064 char* field=fields[f];
1065 for (int j=0;j<nsteps;j++){
1066 int counter=f*nsteps+j;
1067 xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
1068 if(i==(lower_row+1)){
1069 if(xtype[counter]==1){
1070 x=xNew<IssmDouble>(range);
1071 x[0]=scalar;
1072 xs[counter]=x;
1073 xsize[counter]=range;
1074 }
1075 else if (xtype[counter]==3){
1076 x=xNew<IssmDouble>(nindices*range);
1077 for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
1078 xs[counter]=x;
1079 xsize[counter]=range*nindices;
1080 }
1081 else _error_("cannot carry out statistics on type " << xtype[counter]);
1082 }
1083 else{
1084 if(xtype[counter]==1){
1085 x=xs[counter];
1086 x[i-(lower_row+1)]=scalar;
1087 xs[counter]=x;
1088 }
1089 else if (xtype[counter]==3){
1090 x=xs[counter];
1091 for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
1092 xs[counter]=x;
1093 }
1094 else _error_("cannot carry out statistics on type " << xtype[counter]);
1095 }
1096 }
1097 }
1098 fclose(fid);
1099
1100 /*delete buffer:*/
1101 xDelete<char>(buffer);
1102 }
1103 ISSM_MPI_Barrier(IssmComm::GetComm());
1104 _printf0_("Done reading files, now assembling time series.\n");
1105
1106 for (int f=0;f<nfields;f++){
1107 for (int j=0;j<nsteps;j++){
1108 int counter=f*nsteps+j;
1109 if (xtype[counter]==1){
1110 /*we are broadcasting range times doubles:*/
1111 x=xs[counter];
1112 allx=xNew<IssmDouble>(nsamples);
1113 MPI_Gather(x, range, ISSM_MPI_PDOUBLE,allx, range, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
1114 /*add to results:*/
1115 if(my_rank==0){
1116 char fieldname[1000];
1117
1118 sprintf(fieldname,"%s%s",fields[f],"Samples");
1119 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,1,j+1,0));
1120 }
1121 }
1122 else{
1123 /*we are broadcasting double arrays:*/
1124 x=xs[counter];
1125 allx=xNew<IssmDouble>(nsamples*nindices);
1126
1127 MPI_Gather(x, range*nindices, ISSM_MPI_PDOUBLE,allx, range*nindices, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
1128
1129 /*add to results:*/
1130 if(my_rank==0){
1131 char fieldname[1000];
1132 sprintf(fieldname,"%s%s",fields[f],"Samples");
1133 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,nindices,j+1,0));
1134 }
1135 }
1136 }
1137 }
1138
1139
1140} /*}}}*/
1141int OutputStatistics(Parameters* parameters,Results* results){ /*{{{*/
1142
1143 char outputfilename[1000];
1144 char* directory=NULL;
1145 char* model=NULL;
1146 char* method=NULL;
1147 int nsamples;
1148 int* steps=NULL;
1149 int nsteps;
1150
1151 FemModel* femmodel=new FemModel();
1152
1153 /*Some parameters that will allow us to use the OutputResultsx module:*/
1154 parameters->AddObject(new BoolParam(QmuIsdakotaEnum,false));
1155 parameters->AddObject(new BoolParam(SettingsIoGatherEnum,true));
1156
1157 parameters->FindParam(&directory,DirectoryNameEnum);
1158 parameters->FindParam(&model,InputFileNameEnum);
1159 parameters->FindParam(&nsamples,QmuNsampleEnum);
1160 parameters->FindParam(&steps,&nsteps,StepsEnum);
1161
1162 sprintf(outputfilename,"%s/%s.stats",directory,model);
1163 parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
1164
1165 /*Call OutputResults module:*/
1166 femmodel->parameters=parameters;
1167 femmodel->results=results;
1168
1169 OutputResultsx(femmodel);
1170} /*}}}*/
Note: See TracBrowser for help on using the repository browser.