1 /* sim.c
2
3 Copyright (c) 1993-2017 Free Software Foundation, Inc.
4
5 This file is part of GNU MCSim.
6
7 GNU MCSim is free software; you can redistribute it and/or
8 modify it under the terms of the GNU General Public License
9 as published by the Free Software Foundation; either version 3
10 of the License, or (at your option) any later version.
11
12 GNU MCSim is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with GNU MCSim; if not, see <http://www.gnu.org/licenses/>
19
20 Entry point and main simulation routines for 'sim' program.
21
22 This file can use the SUNDIALS CVODES libraries if they are installed.
23 If so, the symbols HAVE_LIBSUNDIALS_CVODES and HAVE_LIBSUNDIALS_NVECSERIAL
24 should be defined in config.h (or in the makefile).
25 Otherwise, the corresponding features are disabled.
26
27 */
28
29 #include <stdio.h>
30 #include <stdlib.h>
31 #include <string.h>
32 #include <ctype.h>
33 #include <math.h>
34 #include <assert.h>
35
36 #include "delays.h"
37 #include "yourcode.h"
38 #include "getopt.h"
39 #include "mh.h"
40 #include "optdsign.h"
41 #include "lexerr.h"
42 #include "lsodes.h"
43 #include "sim.h"
44 #include "simi.h"
45 #include "siminit.h"
46 #include "simo.h"
47 #include "simmonte.h"
48 #include "strutil.h"
49 #include "config.h"
50
51 /* CVODES specific includes and routines */
52
53 #ifdef HAVE_LIBSUNDIALS_CVODES
54
55 #include <nvector/nvector_serial.h> /* serial N_Vector types, fcts., macros */
56
57 #include <cvodes/cvodes.h> /* prototypes CVODE fcts. and consts. */
58 #include <cvodes/cvodes_band.h> /* prototype for CVBand */
59 #include <cvodes/cvodes_dense.h> /* prototype for CVDense */
60 #include <sundials/sundials_dense.h> /* definitions DlsMat DENSE_ELEM */
61 #include <sundials/sundials_types.h> /* definition of type realtype */
62 #include <sundials/sundials_math.h> /* definition of ABS and EXP */
63
64
65 typedef struct {
66 int nVars; /* number of state and output variables */
67 } UserData;
68
69 /* Check function return value...
70 opt == 0 means SUNDIALS function allocates memory so check if
71 returned NULL pointer
72 opt == 1 means SUNDIALS function returns a flag so check if
73 flag >= 0
74 opt == 2 means function allocates memory so check if returned
75 NULL pointer */
check_flag(void * flagvalue,char * funcname,int opt)76 static int check_flag(void *flagvalue, char *funcname, int opt)
77 {
78 int *errflag;
79
80 /* Check if SUNDIALS function returned NULL pointer - no memory allocated */
81
82 if (opt == 0 && flagvalue == NULL) {
83 fprintf(stderr, "[SUNDIALS ERROR] %s() failed - returned NULL pointer\n\n",
84 funcname);
85 return(1); }
86
87 /* Check if flag < 0 */
88
89 else if (opt == 1) {
90 errflag = (int *) flagvalue;
91 if (*errflag < 0) {
92 fprintf(stderr, "[SUNDIALS ERROR] %s() failed with flag = %d\n\n",
93 funcname, *errflag);
94 return(1); }}
95
96 /* Check if function returned NULL pointer - no memory allocated */
97
98 else if (opt == 2 && flagvalue == NULL) {
99 fprintf(stderr, "[MEMORY ERROR] %s() failed - returned NULL pointer\n\n",
100 funcname);
101 return(1); }
102
103 return(0);
104 }
105
106
107 /* ----------------------------------------------------------------------------
108 f_for_cvodes
109
110 stupid routine interfacing cvodes derivative function call and CalcDeriv.
111 */
f_for_cvodes(realtype t,N_Vector u,N_Vector udot,void * user_data)112 static int f_for_cvodes(realtype t, N_Vector u, N_Vector udot, void *user_data)
113 {
114 static realtype *rgMVars = NULL;
115 static int nStates, nVars;
116 int i;
117
118 /* problem with the output variables, derivatives have the proper length
119 and do not need to be translated */
120 if (rgMVars == NULL) { /* initialize */
121
122 nStates = NV_LENGTH_S(udot);
123
124 /* Extract number of outputs from user_data */
125 nVars = ((UserData *) user_data)->nVars;
126
127 /* state and output vector */
128 rgMVars = (realtype *) malloc(nVars * sizeof(realtype));
129
130 if (/*!dudata ||*/ !rgMVars)
131 ReportError (NULL, RE_OUTOFMEM | RE_FATAL, "f_for_cvodes", NULL);
132 }
133
134 /* copy u to rgMVars start */
135 for (i = 0; i < nStates; i++) {
136 rgMVars[i] = NV_Ith_S(u, i);
137 }
138
139 CalcDeriv ((PDOUBLE) rgMVars, (PDOUBLE) NV_DATA_S(udot), (PDOUBLE) &t);
140
141 return(0);
142
143 } /* f_for_cvodes */
144
145 #endif
146
147 /* MPI specific includes */
148
149 #ifdef USEMPI
150 #include "mpi.h"
151 #endif
152
153
154 /* ----------------------------------------------------------------------------
155 CorrectInputToTransition
156
157 resets the integrator and inputs when an input transition occurs.
158
159 returns the simulation time pexp->dTime and input values to
160 the input discontinuity, or transition point *pdTtrans.
161
162 The inputs are updated to reflect their state just after the
163 transition. The integrator is initialized for a new segment.
164
165 This does NOT affect state and output definitions.
166 */
CorrectInputToTransition(PEXPERIMENT pexp,PDOUBLE pdTtrans)167 void CorrectInputToTransition (PEXPERIMENT pexp, PDOUBLE pdTtrans)
168 {
169 pexp->dTime = *pdTtrans;
170 UpdateInputs (&pexp->dTime, pdTtrans);
171
172 } /* CorrectInputToTransition */
173
174
175 /* ----------------------------------------------------------------------------
176 Euler
177
178 Simple Euler integrator.
179 */
Euler(long neq,double * y,double * t,double tout,double dTStep)180 int Euler (long neq, double *y, double *t, double tout, double dTStep)
181 {
182 static PDOUBLE rgdDeriv;
183 double dTmp_step;
184 long i;
185
186 if (!(rgdDeriv))
187 if ( !(rgdDeriv = InitdVector (neq)))
188 ReportError (NULL, RE_OUTOFMEM | RE_FATAL, "Euler", NULL);
189
190 /* Iterate through time out to prescrebed output */
191 while (*t < tout) {
192
193 /* Compute derivatives at current time */
194 CalcDeriv (y, rgdDeriv, t);
195
196 /* Update time */
197 *t = *t + dTStep;
198
199 /* But do not exceed prescribed tout */
200 if (*t > tout) {
201 dTmp_step = tout - (*t - dTStep);
202 *t = tout;
203 }
204 else
205 dTmp_step = dTStep;
206
207 /* Update the state variables */
208 for (i = 0; i < neq; i++)
209 y[i] = y[i] + dTmp_step * rgdDeriv[i];
210
211 if (bDelays)
212 StoreDelayed(*t);
213
214 DoStep_by_Step();
215 }
216
217 /* Calculate the derivatives at t = tout for cleanliness */
218 CalcDeriv (y, rgdDeriv, t);
219
220 return (0);
221
222 } /* Euler */
223
224
225 /* ----------------------------------------------------------------------------
226 FreeVarMod
227
228 Frees a VARMOD structure
229 */
FreeVarMod(PVOID pData)230 void FreeVarMod (PVOID pData)
231 {
232 PVARMOD pvarmod = (PVARMOD) pData;
233
234 if (IsInput (pvarmod->hvar))
235 if (pvarmod->uvar.pifn)
236 free (pvarmod->uvar.pifn);
237
238 free (pvarmod);
239
240 /* we might have to free those too... */
241 /* PDOUBLE rgT0s; /\* Array of start times *\/ */
242 /* PDOUBLE rgMags; /\* Array of magnitudes *\/ */
243 /* HANDLE *rghT0s; /\* Handles to start times *\/ */
244 /* HANDLE *rghMags; /\* Handles to magnitudes *\/ */
245 /* PINT rgOper; /\* Array of operation types *\/ */
246
247 } /* FreeVarMod */
248
249
250 /* ----------------------------------------------------------------------------
251 ModifyOneParm
252
253 Callback function for ModifyParms.
254 */
255
ModifyOneParm(PVOID pData,PVOID pNullInfo)256 int ModifyOneParm (PVOID pData, PVOID pNullInfo)
257 {
258 PVARMOD pvarmod = (PVARMOD) pData;
259
260 if (IsInput(pvarmod->hvar))
261 SetInput (pvarmod->hvar, pvarmod->uvar.pifn);
262 else
263 SetVar (pvarmod->hvar, pvarmod->uvar.dVal);
264
265 return 0;
266
267 } /* ModifyOneParm */
268
269
270 /* ----------------------------------------------------------------------------
271 ModifyParms
272
273 Modifies the parameters in the plistParmMods LIST of the experiment
274 spec by call ForAllList to increment through the list.
275 */
ModifyParms(PLIST plistParmMods)276 void ModifyParms (PLIST plistParmMods)
277 {
278
279 assert (plistParmMods);
280 ForAllList (plistParmMods, &ModifyOneParm, NULL);
281
282 } /* ModifyParms */
283
284
285 /* ----------------------------------------------------------------------------
286 DoOneExperiment
287
288 Runs one experiment - return 1 on success and 0 in case of errors
289 */
DoOneExperiment(PEXPERIMENT pexp)290 int DoOneExperiment (PEXPERIMENT pexp)
291 {
292
293 double dTout; /* next output time */
294 double dTtrans; /* next exposure transition time */
295 double dTup; /* the smaller one of dTout or dTtrans*/
296 int iOut; /* index to next output time */
297 PMODELINFO pmod; /* pointer to the current model info */
298 PINTSPEC pis; /* pointer to the integrator specs */
299
300 #ifdef HAVE_LIBSUNDIALS_CVODES
301 /* CVODES specific variables */
302 static N_Vector u = NULL;
303 static UserData user_data;
304 static void *cvode_mem = NULL;
305 int flag, i;
306 #endif
307
308 if (!pexp) return 0;
309
310 pmod = pexp->pmodelinfo;
311 pis = &(pexp->is);
312
313 if (!InitOutputs (pexp, &iOut, &dTout)) return 0;
314
315 /* Resolve dependency for initial time, eventually */
316 if (pexp->hT0)
317 pexp->dT0 = GetVarValue (pexp->hT0);
318
319 /* Resolve dependent inputs, which calls ScaleModel */
320 UpdateInputs (&pexp->dT0, &dTtrans);
321
322 if (bDelays)
323 InitDelays(pexp->hT0);
324
325 if (pexp->dT0 > dTtrans) {
326 printf ("\nError: starting time is greater than first discontinuity,"
327 " check your inputs - Exiting.\n\n");
328 exit (0);
329 }
330
331 if (pexp->dT0 > dTout) {
332 printf ("\nError: starting time is greater than first output time,"
333 " check your outputs - Exiting.\n\n");
334 exit (0);
335 }
336
337 pexp->dTime = pexp->dT0;
338
339 /* integrator initializations */
340 if (pis->iAlgo == IAL_LSODES) { /* Lsodes algorithm */
341 /* set lsodes return flag to 1 for first call */
342 pis->iDSFlag = 1;
343 }
344 else {
345 if (pis->iAlgo == IAL_CVODES) { /* Sundials CVODE algorithm */
346
347 #ifdef HAVE_LIBSUNDIALS_CVODES
348
349 if (1 || u == NULL) { /* always done for now */
350
351 /* create a serial vector for state variables */
352 u = N_VNew_Serial(pmod->nStates); /* Allocate u vector */
353 if (check_flag((void*)u, "N_VNew_Serial", 0)) return(1);
354
355 /* call CVodeCreate to create the solver memory and specify the
356 Backward Differentiation Formula and the use of a Newton iteration */
357 cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
358 if (check_flag((void *)cvode_mem, "CVodeCreate", 0)) return(1);
359
360 /* set initial state values */
361 for (i = 0; i < pmod->nStates; i++)
362 NV_Ith_S(u, i) = pmod->pdModelVars[i];
363
364 /* initialize the integrator memory and specify the
365 user's right hand side function in u'=f(t,u), the inital time T0,
366 and the initial dependent variable vector u. */
367 flag = CVodeInit(cvode_mem, f_for_cvodes, pexp->hT0, u);
368 if (check_flag(&flag, "CVodeInit", 1)) return(1);
369
370 /* set the scalar relative tolerance and scalar absolute tolerance */
371 flag = CVodeSStolerances(cvode_mem,
372 RCONST(pis->dRtol), RCONST(pis->dAtol));
373 if (check_flag(&flag, "CVodeSStolerances", 1)) return(1);
374
375 /* set the pointer to user-defined data used to store the total
376 number of state and output variables */
377 user_data.nVars = pmod->nModelVars;
378 flag = CVodeSetUserData(cvode_mem, &user_data);
379 if (check_flag(&flag, "CVodeSetUserData", 1)) return(1);
380
381 /* set the maximum number of internal steps before t_out */
382 flag = CVodeSetMaxNumSteps(cvode_mem, pis->maxsteps);
383 if (check_flag(&flag, "CVodeSetMaxNumSteps", 1)) return(1);
384
385 /* set the maximum number of error test failures */
386 flag = CVodeSetMaxErrTestFails(cvode_mem, pis->maxnef);
387 if (check_flag(&flag, "CVodeSetMaxErrTestFails", 1)) return(1);
388
389 /* set the maximum number of nonlinear iterations */
390 flag = CVodeSetMaxNonlinIters(cvode_mem, pis->maxcor);
391 if (check_flag(&flag, "CVodeSetMaxNonlinIters", 1)) return(1);
392
393 /* set the maximum number of convergence failures */
394 flag = CVodeSetMaxConvFails(cvode_mem, pis->maxncf);
395 if (check_flag(&flag, "CVodeSetMaxConvFails", 1)) return(1);
396
397 /* set the oefficient in the nonlinear convergence test */
398 flag = CVodeSetNonlinConvCoef(cvode_mem, RCONST(pis->nlscoef));
399 if (check_flag(&flag, "CVodeSetNonlinConvCoef", 1)) return(1);
400
401 /* call CVDense to specify the CVDENSE dense linear solver */
402 flag = CVDense(cvode_mem, pmod->nStates);
403 if (check_flag(&flag, "CVDense", 1)) return(1);
404
405 /* set the user-supplied Jacobian routine Jac, not used now
406 flag = CVDlsSetBandJacFn(cvode_mem, Jac);
407 if(check_flag(&flag, "CVDlsSetBandJacFn", 1)) return(1); */
408
409 }
410 else { /* disabled for now */
411 /* reset initial state values */
412 for (i = 0; i < pmod->nStates; i++)
413 NV_Ith_S(u, i) = pmod->pdModelVars[i];
414
415 flag = CVodeReInit(cvode_mem, pexp->dT0, u);
416 }
417 #endif
418 } /* end if IAL_CVODES */
419 }
420
421 /* Iterate to final time */
422 while (pexp->dTime < pexp->dTfinal) {
423
424 /* If dynamics equations are defined */
425 if (pmod->nStates > 0) {
426
427 /* the upper limit of integration dTup should be either dTout
428 or dTtrans, whichever is smaller */
429 /* F. Bois, 08 April 2007: before that, if the difference between dTout
430 and dTtrans is too small make dTtrans = dTout to avoid problems with
431 the integration */
432 if (fabs(dTout - dTtrans) < DBL_EPSILON * 2.0 *
433 mymax(fabs(dTout), fabs(dTtrans)))
434 dTtrans = dTout;
435
436 dTup = (dTout < dTtrans) ? dTout : dTtrans;
437
438 /* F. Bois, 07 October 2008: same rounding problem fix: */
439 if (fabs(dTup - pexp->dTime) < DBL_EPSILON * 2.0 *
440 mymax(fabs(dTup), fabs(pexp->dTime)))
441 pexp->dTime = dTup;
442
443 if (pis->iAlgo == IAL_LSODES) { /* Lsodes algorithm */
444
445 pis->rwork[0] = dTup; /* do not overshoot dTup - FB 01/07/97 */
446
447 lsodes_(&pmod->nStates, pmod->pdModelVars, &(pexp)->dTime,
448 &dTup, &pis->itol, &pis->dRtol, &pis->dAtol,
449 &pis->itask, &pis->iDSFlag, &pis->iopt, pis->rwork,
450 &pis->lrw, pis->iwork, &pis->liw, &pis->iMf);
451
452 /* Handle error returns : FB 25/11/96 : */
453 if (pis->iDSFlag < 0) {
454 /* We cannot guarantee the accuracy of the results, exit the routine
455 with an error flag */
456 return (0);
457 }
458 }
459 else {
460 if (pis->iAlgo == IAL_CVODES) { /* Sundials CVODE algorithm */
461
462 #ifdef HAVE_LIBSUNDIALS_CVODES
463 if (dTup > (pexp)->dTime) {
464 /* do not overshoot */
465 flag = CVodeSetStopTime(cvode_mem, (realtype) dTup);
466 flag = CVode(cvode_mem, (realtype) dTup, u, &(pexp)->dTime,
467 CV_NORMAL);
468 if (check_flag(&flag, "CVode", 1)) {
469 /* we cannot guarantee the accuracy of the results, exit routine
470 with an error flag */
471 return (0);
472 }
473 /* copy back state values */
474 for (i = 0; i < pmod->nStates; i++) {
475 pmod->pdModelVars[i] = NV_Ith_S(u, i);
476 }
477 }
478 #endif
479 }
480 else if (pis->iAlgo == IAL_EULER) { /* Euler algorithm */
481 Euler(pmod->nStates, pmod->pdModelVars, &(pexp)->dTime, dTup,
482 pis->dTStep);
483 }
484 }
485 }
486 else {
487 /* We still need to advance the time */
488 pexp->dTime = (dTout < dTtrans) ? dTout : dTtrans;
489 }
490
491 if (dTtrans <= dTout) {
492 /* dTime == dTtrans <= dTout: we are at a discontinuity.
493 This point belongs to the NEW period UNLESS we are at
494 the final time */
495 if (dTtrans < dTout) {
496 if (dTtrans < pexp->dTfinal) {
497 CorrectInputToTransition (pexp, &dTtrans);
498 pis->iDSFlag = 1;
499 }
500 }
501 else {
502 /* dTtrans == dTout */
503 if (dTtrans < pexp->dTfinal) {
504 CorrectInputToTransition (pexp, &dTtrans);
505 pis->iDSFlag = 1;
506 }
507 SaveOutputs (pexp, &dTout);
508 NextOutputTime (pexp, &dTout, &iOut);
509 }
510 }
511 else {
512 /* dTime == dTout < dTtrans: */
513 SaveOutputs (pexp, &dTout);
514 NextOutputTime (pexp, &dTout, &iOut);
515 }
516
517 } /* while dTime < final time */
518
519 if (pis->iAlgo == IAL_CVODES) { /* cleanup Sundials CVODE algorithm */
520 #ifdef HAVE_LIBSUNDIALS_CVODES
521 /* Free vector u */
522 N_VDestroy_Serial(u);
523 CVodeFree(&cvode_mem); /* Free the integrator memory */
524 #endif
525 }
526
527 /* success */
528 return 1;
529
530 } /* DoOneExperiment */
531
532
533 /* ----------------------------------------------------------------------------
534 DoOneNormalExp
535
536 Does one AT_DEFAULTSIM simulation.
537
538 Return 1 on success and 0 in case of failure
539 */
DoOneNormalExp(PANALYSIS panal,PEXPERIMENT pexp)540 int DoOneNormalExp (PANALYSIS panal, PEXPERIMENT pexp)
541 {
542 printf("experiment %d\n", pexp->iExp); /* Show what experiment it is */
543
544 InitModel();
545 ModifyParms(panal->expGlobal.plistParmMods); /* Global modifications */
546 ModifyParms(pexp->plistParmMods); /* Mods for this experiment */
547 if (!DoOneExperiment(pexp)) {
548 /* Error */
549 return 0;
550 }
551
552 return (1);
553
554 } /* DoOneNormalExp */
555
556
557 /* ----------------------------------------------------------------------------
558 DoOneMCExp
559
560 Does one AT_MONTECARLO simulation.
561
562 Can maybe merge this with DoOneNormalExp() in the future.
563
564 The major issue is the order of setting parameters. For each
565 experiment in a Monte Carlo run of an analysis, the order must be
566 as follows:
567
568 Each Run
569 calc mc mods
570
571 Each Simulation (old "Experiment")
572 1) Init the model
573 2) Global parm mods
574 3) Monte Carlo mods
575 4) Local mods override everything
576
577 The problem becomes that for the simulation to be started over
578 again, the inputs have to be told to initialize and parm mods for
579 the current experiment must be made (body weight, etc). This
580 currently won't happen unless the model is init'd. Maybe create a
581 ResetInputs() with starting time which will do the funky stuff
582 done by the global variable right now.
583
584 Return 1 on success and 0 in case of failure
585 */
DoOneMCExp(PANALYSIS panal,PEXPERIMENT pexp)586 int DoOneMCExp (PANALYSIS panal, PEXPERIMENT pexp)
587 {
588 register MONTECARLO *pmc = &panal->mc;
589
590 InitModel ();
591 ModifyParms (panal->expGlobal.plistParmMods); /* Global modifications */
592 SetParms (pmc->nParms, pmc->rghvar, pmc->rgdParms); /* MC mods */
593 ModifyParms (pexp->plistParmMods); /* Mods for this experiment */
594 if (!DoOneExperiment (pexp)) {
595 /* Error */
596 return 0;
597 }
598
599 return (1);
600
601 } /* DoOneMCExp */
602
603
604 /* ----------------------------------------------------------------------------
605 DoNormal
606
607 Does a normal analysis
608 */
DoNormal(PANALYSIS panal)609 void DoNormal (PANALYSIS panal)
610 {
611 int nExps = panal->expGlobal.iExp;
612 int i;
613
614 printf ("\nDoing analysis - %d normal experiment%c\n", nExps,
615 (nExps > 1 ? 's' : ' '));
616
617 for (i = 0; i < nExps; i++) {
618 if (DoOneNormalExp (panal, panal->rgpExps[i])) {
619 /* if successfull write out the results */
620 WriteNormalOutput (panal, panal->rgpExps[i]);
621 }
622 else
623 printf("[MCSIM ERROR] Integration failed - No output generated\n\n");
624 }
625
626 } /* DoNormal */
627
628
629 /* ----------------------------------------------------------------------------
630 DoMonteCarlo
631
632 Does Monte Carlo simulations.
633 */
DoMonteCarlo(PANALYSIS panal)634 void DoMonteCarlo (PANALYSIS panal)
635 {
636 int nExps = panal->expGlobal.iExp;
637 long nRuns = panal->mc.nRuns;
638 MCPREDOUT mcpredout;
639 BOOL bOK;
640 long i, j;
641
642 if (panal->rank == 0) {
643 printf("Doing %ld Monte Carlo simulation%c, %d experiment%c%s\n",
644 nRuns, (nRuns != 1 ? 's' : ' '),
645 nExps, (nExps > 1 ? 's' : ' '), (nRuns != 1 ? " each" : "."));
646 if (panal->size > 1)
647 printf("Split between %d processors\n", panal->size);
648 }
649 else
650 printf("\n");
651
652 SetParents(&panal->mc, 0); /* do all requested simulations */
653
654 OpenMCFiles(panal);
655 mcpredout.pred = NULL;
656
657 #ifdef USEMPI
658 MPI_Barrier(MPI_COMM_WORLD);
659 #endif
660
661 /* split the work between processors */
662 for (i = panal->rank; i < nRuns; i += panal->size) {
663
664 /* output to screen, if required as a command-line option */
665 if (i == 0)
666 printf("\n");
667 if (panal->bOutputIter && ((i+1) % panal->nOutputFreq == 0)) {
668 if (panal->size > 1)
669 printf("Processor %d, Iteration %ld\n", panal->rank, i + 1);
670 else
671 printf("Iteration %ld\n", i + 1);
672 }
673
674 panal->mc.lRun = i;
675
676 CalcMCParms (&panal->mc, NULL, 0); /* start at 0, do them all */
677
678 for (j = 0; j < nExps; j++) { /* do all experiments */
679 bOK = DoOneMCExp (panal, panal->rgpExps[j]);
680 if (!bOK)
681 break;
682 }
683
684 if (bOK) { /* if successful write results */
685 TransformPred(panal, &mcpredout); /* transform output run */
686 WriteMCOutput(panal, &mcpredout);
687 }
688 else
689 printf("Warning: Integration failed on iteration %ld, experiment %ld:\n"
690 " No output generated\n", panal->mc.lRun+1, j+1);
691
692 } /* for i */
693
694 CloseMCFiles(panal);
695 if (mcpredout.pred)
696 free(mcpredout.pred);
697
698 } /* DoMonteCarlo */
699
700
701 /* ----------------------------------------------------------------------------
702 DoSetPoints
703
704 Does Set Points analysis.
705
706 If the number of runs (nRuns) is specified to be zero, set points
707 are read from the set points file until end of file is reached.
708 Otherwise, the number of runs explicity stated are read. Not
709 having enough points in the file in this latter case yields an
710 error.
711
712 If nRuns == 0, the test at the end of the while{} loop is not
713 used, and the decision to continue is made by the return value of
714 GetMCMods().
715 */
DoSetPoints(PANALYSIS panal)716 void DoSetPoints (PANALYSIS panal)
717 {
718 int nExps = panal->expGlobal.iExp;
719 long nRuns = panal->mc.nRuns;
720 MCPREDOUT mcpredout;
721 BOOL bOK = FALSE, bNotDone; /* Not finished with analysis */
722 int i;
723
724 mcpredout.pred = NULL;
725
726 OpenMCFiles(panal);
727
728 if (panal->rank == 0) {
729 printf ("Doing analysis - %ld SetPoints run%c... %d experiment%c%s\n",
730 nRuns, (nRuns != 1 ? 's' : ' '), nExps, (nExps > 1 ? 's' : ' '),
731 (nRuns != 1 ? " each" : " "));
732 if (panal->size > 1)
733 printf("Split between %d processors\n", panal->size);
734 }
735 else
736 printf("\n");
737
738 if ((!nRuns) && panal->rank == 0)
739 printf("0 runs specified for SetPoint(). Reading entire file.\n\n");
740
741 SetParents(&panal->mc, panal->mc.nSetParms);
742
743 #ifdef USEMPI
744 MPI_Barrier(MPI_COMM_WORLD);
745 #endif
746
747 panal->mc.lRun = 0; /* first run */
748 bNotDone = TRUE;
749
750 while (bNotDone) {
751
752 bNotDone = GetSPMods (panal, NULL);
753
754 /* split the work between processors */
755 if ((bNotDone) && (panal->mc.lRun % panal->size == panal->rank)) {
756
757 /* output to screen, if required as a command-line option */
758 if (panal->bOutputIter &&
759 ((panal->mc.lRun+1) % panal->nOutputFreq == 0)) {
760 if (panal->size > 1)
761 printf("Processor %d, Iteration %ld\n",
762 panal->rank, panal->mc.lRun + 1);
763 else
764 printf("Iteration %ld\n", panal->mc.lRun + 1);
765 }
766
767 /* do analysis if not finished */
768 for (i = 0; i < nExps; i++) { /* do all experiments */
769 bOK = DoOneMCExp (panal, panal->rgpExps[i]);
770 if (!bOK) break;
771 }
772
773 if (bOK) {
774 /* if successful write results */
775 TransformPred (panal, &mcpredout); /* transform output run */
776 WriteMCOutput (panal, &mcpredout);
777 }
778 else
779 printf ("Warning: Integration failed on iteration %ld, experiment %d:\n"
780 " No output generated\n", panal->mc.lRun+1, i+1);
781 } /* if bNotDone */
782
783 panal->mc.lRun++; /* next run */
784 if (nRuns) /* if a number of runs spec'd... */
785 bNotDone = (panal->mc.lRun < nRuns);
786
787 } /* while */
788
789 CloseMCFiles (panal);
790
791 if (mcpredout.pred) free(mcpredout.pred);
792
793 } /* DoSetPoints */
794
795
796 /* ----------------------------------------------------------------------------
797 DoAnalysis
798
799 Does the analysis in the given specification.
800 */
DoAnalysis(PANALYSIS panal)801 void DoAnalysis (PANALYSIS panal)
802 {
803
804 if (panal->size == 1)
805 InitRandom (panal->rank, panal->dSeed, TRUE);
806 else
807 InitRandom (panal->rank, panal->dSeed + panal->rank, TRUE);
808
809 switch (panal->iType) {
810
811 default:
812 case AT_DEFAULTSIM:
813 if (panal->rank == 0) /* not parallel */
814 DoNormal (panal);
815 break;
816
817 case AT_SETPOINTS:
818 DoSetPoints (panal);
819 break;
820
821 case AT_MONTECARLO:
822 DoMonteCarlo (panal);
823 break;
824
825 case AT_MCMC:
826 DoMarkov (panal);
827 break;
828
829 case AT_OPTDESIGN:
830 if (panal->rank == 0) /* not parallel */
831 DoOptimalDesign (panal);
832 break;
833
834 } /* switch */
835
836 if (panal->pfileOut) {
837 fclose (panal->pfileOut);
838 printf ("Wrote output file \"%s\"\n", panal->szOutfilename);
839 }
840
841 } /* DoAnalysis */
842
843
844 /* ----------------------------------------------------------------------------
845 FreeMemory
846
847 To use in the case of simulations without Levels.
848 */
FreeMemory(PANALYSIS panal)849 void FreeMemory (PANALYSIS panal)
850 {
851 int i, j;
852
853 free(panal->modelinfo.pStateHvar);
854
855 FreeList (&panal->mc.plistMCVars, NULL, TRUE);
856 if (panal->mc.rgdParms) {
857 free (panal->mc.rgdParms);
858 free (panal->mc.rghvar);
859 }
860
861 PINTSPEC pis = &panal->rgpExps[0]->is;
862 free (pis->iwork);
863 free (pis->rwork);
864
865 for (i = 0; i < panal->expGlobal.iExp; i++) {
866 if (panal->rgpExps[i] != NULL) {
867 FreeList (&panal->rgpExps[i]->plistParmMods, &FreeVarMod, TRUE);
868
869 POUTSPEC pos = &panal->rgpExps[i]->os;
870 free (pos->pszOutputNames);
871 free (pos->phvar_out);
872 free (pos->pcOutputTimes);
873 free (pos->piCurrentOut);
874 free (pos->prgdOutputTimes);
875 for (j = 0; j < pos->nOutputs; j++)
876 free(pos->prgdOutputVals[j]);
877 free (pos->prgdOutputVals);
878 free (pos->rgdDistinctTimes);
879 ForAllList (pos->plistPrintRecs, &FreePrintRec, NULL);
880 FreeList (&pos->plistPrintRecs, NULL, FALSE);
881 free (pos->plistPrintRecs);
882 ForAllList (pos->plistDataRecs, &FreeDataRec, NULL);
883 FreeList (&pos->plistDataRecs, NULL, FALSE);
884 free (pos->plistDataRecs);
885 free (panal->rgpExps[i]);
886 }
887 }
888 if (panal->bAllocatedFileName) {
889 if (panal->szOutfilename) free (panal->szOutfilename);
890 if (panal->mc.szMCOutfilename) free (panal->mc.szMCOutfilename);
891 if (panal->gd.szGout) free (panal->gd.szGout);
892 }
893
894 if (panal->mc.szSetPointsFilename) free (panal->mc.szSetPointsFilename);
895 if (panal->gd.szGrestart) free (panal->gd.szGrestart);
896 if (panal->gd.szGdata) free (panal->gd.szGdata);
897
898 FreeList (&panal->expGlobal.plistParmMods, NULL, TRUE);
899 free (panal);
900
901 } /* FreeMemory */
902
903
904 /* ----------------------------------------------------------------------------
905 MCVarListToArray
906
907 converts a list of MCVAR to an array. This must be a callback for
908 ForAllList() since we are making the change here that will let us
909 not to be forced to use list traversal in the future.
910 */
911
912 MCVAR **vrgpMCVar; /* Avoid hairy pointers in here */
913 int viMCVar; /* Index to the array */
914
MCVarListToArray(PVOID pv_pMCVar,PVOID pv_Null)915 int MCVarListToArray (PVOID pv_pMCVar, PVOID pv_Null)
916 {
917
918 vrgpMCVar[viMCVar] = (MCVAR *) pv_pMCVar; /* Copy the pointer and.. */
919 viMCVar++; /* Advance to next element of array */
920 return 1;
921
922 } /* MCVarListToArray */
923
924
925 /* ----------------------------------------------------------------------------
926 PrepAnalysis
927
928 makes the ANALYSIS structure easier to work with in the simulation
929 code. Specifically, changes lists to arrays.
930 */
PrepAnalysis(PANALYSIS panal)931 void PrepAnalysis (PANALYSIS panal)
932 {
933 register MONTECARLO *pmc = &panal->mc;
934 register int l;
935
936 pmc->nParms = ListLength (pmc->plistMCVars);
937 /* avoid zero pmc->nParms which can confuse some implementations of
938 malloc. If pmc->nParms is zero no use is going to be made of these
939 arrays anyway */
940 if (pmc->nParms == 0) return;
941
942 pmc->rgdParms = InitdVector (pmc->nParms);
943 pmc->rgpMCVar = (MCVAR **) malloc((pmc->nParms)*sizeof(MCVAR *));
944 if (!(pmc->rgdParms && pmc->rgpMCVar))
945 ReportError (NULL, RE_OUTOFMEM | RE_FATAL, "PrepAnalysis", NULL);
946
947 /* Address of the first pointer */
948 vrgpMCVar = &pmc->rgpMCVar[0];
949
950 /* Initialize global array index */
951 viMCVar = 0;
952 ForAllList (pmc->plistMCVars, MCVarListToArray, (PVOID) NULL);
953 FreeList (&pmc->plistMCVars, NULL, FALSE);
954
955 /* Make a handle vector */
956 pmc->rghvar = (HVAR *) malloc((pmc->nParms)*sizeof(HVAR));
957 if (pmc->rghvar) {
958 for (l = 0; l < pmc->nParms; l++)
959 pmc->rghvar[l] = pmc->rgpMCVar[l]->hvar;
960 }
961 else
962 ReportError (NULL, RE_OUTOFMEM | RE_FATAL, "PrepAnalysis", NULL);
963
964 } /* PrepAnalysis */
965
966
967 /* ----------------------------------------------------------------------------
968 PromptFilenames
969
970 prompts for both input and output file names. The space allocated
971 for inputting the files is reallocated to their actual size.
972 */
PromptFilenames(PSTR * pszFileIn,PSTR * pszFileOut)973 void PromptFilenames (PSTR *pszFileIn, PSTR *pszFileOut)
974 {
975 *pszFileIn = (PSTR) calloc (1, MAX_FILENAMESIZE);
976 *pszFileOut = (PSTR) calloc (1, MAX_FILENAMESIZE);
977
978 printf ("Input filename? ");
979
980 if (!fgets (*pszFileIn, MAX_FILENAMESIZE, stdin)) {
981 ReportError (NULL, RE_READERROR | RE_FATAL, "stdin", NULL);
982 }
983 else
984 *pszFileIn = strtok (*pszFileIn, " \t\n");
985
986 if (!(*pszFileIn)) /* Nothing entered, quit */
987 return;
988
989 if ((*pszFileIn)[0]) { /* Input file specified */
990 printf ("Output filename? ");
991
992 if (!fgets (*pszFileOut, MAX_FILENAMESIZE, stdin)) {
993 ReportError (NULL, RE_READERROR | RE_FATAL, "stdin", NULL);
994 }
995 else
996 *pszFileOut = strtok (*pszFileOut, " \t\n");
997 }
998
999 if (!(*pszFileOut) || !(*pszFileOut)[0]) { /* If no output specified */
1000 free (*pszFileOut); /* .. use default later */
1001 *pszFileOut = NULL;
1002 }
1003 else {
1004 *pszFileIn = (PSTR) realloc (*pszFileIn, MyStrlen(*pszFileIn) + 1);
1005 *pszFileOut = (PSTR) realloc (*pszFileOut, MyStrlen(*pszFileOut) + 1);
1006 }
1007
1008 } /* PromptFilenames */
1009
1010
1011 /* ----------------------------------------------------------------------------
1012 GetCmdLineArgs
1013
1014 retrieves options and filenames from the command line arguments passed to
1015 the program.
1016
1017 The command line syntax is:
1018
1019 mcsim [-options] [input-file [output-file]]
1020
1021 If the output filename is not given a default is used.
1022 If neither the input, nor output filenames are given, the
1023 program prompts for them both.
1024
1025 The options can appear anywhere in the line and in any order.
1026
1027 The options are parsed with _getopt(). After _getopt() is called,
1028 the args in rgszArg have been permuted so that non-option args are
1029 first, which in this case means the filenames.
1030
1031 Uses the following globals (see getopt.c comments):
1032 char *optarg; -- Contains the string argument to each option in turn
1033
1034 The following options are defined:
1035 -c (without argument) : print MCMC convergence check to stdout,
1036 inhibits -i option if set (redundant)
1037 -h or -H (without argument) : give minimal help
1038 -i (with integer argument) : print one in every x MCMC iteration
1039 counter to stdout; x is given by the
1040 argument
1041 -D (with argument print-hierarchy): print the MCMC Level structure
1042 */
1043 static char vszOptions[] = "c::h::H::i:D:";
1044
GetCmdLineArgs(int cArg,char * const * rgszArg,PSTR * pszFileIn,PSTR * pszFileOut,PANALYSIS panal)1045 void GetCmdLineArgs (int cArg, char *const *rgszArg, PSTR *pszFileIn,
1046 PSTR *pszFileOut, PANALYSIS panal)
1047 {
1048 int c;
1049
1050 *pszFileIn = *pszFileOut = (PSTR) NULL;
1051
1052 while (1) {
1053
1054 c = _getopt (cArg, rgszArg, vszOptions);
1055 if (c == EOF) /* finished with option args */
1056 break;
1057
1058 switch (c) {
1059 case 'c':
1060 #ifdef USEMPI
1061 panal->bPrintConvergence = TRUE;
1062 if (optarg)
1063 printf(">> Option -c argument %s will be ignored\n\n", optarg);
1064 #else
1065 printf(">> MPI parallelization not active: option -c is ignored\n\n");
1066 #endif
1067 break;
1068
1069 case 'D':
1070 /* remove optional leading '=' character */
1071 if (optarg[0] == '=')
1072 optarg++;
1073 if (!strcmp(optarg, "print-hierarchy")) {
1074 printf (">> Debug option %s\n\n", optarg);
1075 panal->bDependents = TRUE;
1076 }
1077 else {
1078 printf(">> A known debugging code must follow -D\nExiting.\n\n");
1079 exit(-1);
1080 }
1081 break;
1082
1083 case 'H':
1084 case 'h':
1085 printf("Usage: %s [options] <input-file> [<output-file>]\n",
1086 rgszArg[0]);
1087 printf("Options:\n");
1088 printf(" -c "
1089 "Display MCMC convergence (if MPI is used)\n");
1090 printf(" -D=print-hierarchy "
1091 "Print out the hierarchy for debugging\n");
1092 printf(" -h "
1093 "Display this information\n");
1094 printf(" -H "
1095 "Display this information\n");
1096 printf(" -i=<arg> "
1097 "Print out every <arg> iteration\n");
1098 printf("\nFor further help on GNU MCSim please see:\n"
1099 "http://www.gnu.org/software/mcsim.\n\n");
1100 exit(-1);
1101 break;
1102
1103 case 'i':
1104 /* remove optional leading '=' character */
1105 if (optarg[0] == '=')
1106 optarg++;
1107 /* convert argument to base 10 */
1108 panal->nOutputFreq = strtol(optarg, NULL, 10);
1109 if (panal->nOutputFreq > 0) {
1110 if (panal->rank == 0)
1111 printf (">> Print iteration frequency %d\n\n", panal->nOutputFreq);
1112 panal->bOutputIter = TRUE;
1113 }
1114 else {
1115 printf(">> Error: An integer print step must "
1116 "follow -i\nExiting.\n\n");
1117 exit(-1);
1118 }
1119 break;
1120
1121 default:
1122 printf("Exiting.\n\n");
1123 exit(-1);
1124
1125 } /* switch */
1126
1127 } /* while */
1128
1129 switch (cArg - optind) { /* Remaining args are filenames */
1130 case 2: /* Output and input file specificed */
1131 *pszFileOut = rgszArg[optind + 1];
1132 /* Fall through! */
1133
1134 case 1: /* Input file specificed */
1135 *pszFileIn = rgszArg[optind];
1136 break;
1137
1138 case 0: /* No file names specified */
1139 PromptFilenames (pszFileIn, pszFileOut);
1140 break;
1141
1142 default:
1143 exit (-1);
1144 break;
1145
1146 } /* switch */
1147
1148 while (*pszFileIn && (*pszFileIn)[0] && /* files specified */
1149 !MyStrcmp(*pszFileIn, *pszFileOut)) { /* but not different */
1150 printf ("\n** Input and output filename must be different.\n");
1151 PromptFilenames (pszFileIn, pszFileOut);
1152 }
1153
1154 if (!(*pszFileIn && (*pszFileIn)[0])) { /* no input name given is an error */
1155 printf ("Error: an input file name must be specified - Exiting.\n\n");
1156 exit (-1);
1157 }
1158
1159 } /* GetCmdLineArgs */
1160
1161
1162 /* ----------------------------------------------------------------------------
1163 */
AnnounceProgram(int iRank)1164 void AnnounceProgram (int iRank)
1165 {
1166 if (iRank == 0) { /* serial */
1167 printf ("\n________________________________________\n");
1168 printf ("\nMCSim " VSZ_VERSION "\n\n");
1169 printf (VSZ_COPYRIGHT "\n\n");
1170
1171 printf ("MCSim comes with ABSOLUTELY NO WARRANTY;\n"
1172 "This is free software, and you are welcome to redistribute it\n"
1173 "under certain conditions; "
1174 "see the GNU General Public License.\n\n");
1175
1176 #ifdef HAVE_LIBGSL
1177 printf ("* Using GNU Scientific Library (GSL)\n\n");
1178 #endif
1179
1180 printf ("* Using `%s' model in file \"%s\" created by %s\n\n",
1181 szModelDescFilename, szModelSourceFilename, szModelGenAndVersion);
1182
1183 } /* end if iRank == 0 */
1184
1185 } /* AnnounceProgram */
1186
1187
1188 /* ----------------------------------------------------------------------------
1189 main
1190
1191 Entry point for simulation and analysis program.
1192 */
main(int nArg,char ** rgszArg)1193 int main (int nArg, char **rgszArg)
1194 {
1195
1196 #ifdef USEMPI
1197 MPI_Init(&nArg,&rgszArg);
1198 #endif
1199
1200 int rank = 0;
1201 int size = 1;
1202
1203 #ifdef USEMPI
1204 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
1205 MPI_Comm_size(MPI_COMM_WORLD, &size);
1206 #endif
1207
1208 PSTR szFileIn, szFileOut;
1209 INPUTBUF ibIn;
1210 PANALYSIS panal = (PANALYSIS) malloc (sizeof(ANALYSIS));
1211
1212 panal->rank = rank;
1213 panal->size = size;
1214
1215 AnnounceProgram(panal->rank);
1216
1217 if (!panal)
1218 ReportError (NULL, RE_OUTOFMEM | RE_FATAL,
1219 "ANALYSIS specification too large", NULL);
1220
1221 InitAnalysis (panal);
1222 GetCmdLineArgs (nArg, rgszArg, &szFileIn, &szFileOut, panal);
1223
1224 /* Define the output file as the global experiment default */
1225 panal->szOutfilename = szFileOut;
1226 szFileOut == NULL ? (panal->bCommandLineSpec = FALSE) :
1227 (panal->bCommandLineSpec = TRUE);
1228
1229 if (!InitBuffer (&ibIn, szFileIn))
1230 ReportError (&ibIn, RE_INIT | RE_FATAL, "ReadInput", NULL);
1231
1232 ibIn.pInfo = (PVOID) panal; /* Attach analysis specification to input */
1233
1234 if (ReadAnalysis (&ibIn)) {
1235 PrepAnalysis (panal);
1236 DoAnalysis (panal);
1237 }
1238
1239 if (panal->rank == 0)
1240 printf("Done.\n\n");
1241
1242 if (panal->iType == AT_MCMC)
1243 FreeLevels (panal);
1244 else {
1245 FreeMemory (panal);
1246 free (ibIn.pbufOrg);
1247 }
1248
1249 #ifdef USEMPI
1250 MPI_Finalize();
1251 #endif
1252
1253 return 0;
1254
1255 } /* main */
1256