1 /* simi.c
2
3 Copyright (c) 1991-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 Contains input routines for simulation.
21 This file checks whether the SUNDIALS CVODES libraries are installed.
22 If so, the symbols HAVE_LIBSUNDIALS_CVODES and HAVE_LIBSUNDIALS_NVECSERIAL
23 should be defined in config.h (or in the makefile).
24 Otherwise, the corresponding features are disabled (and the program
25 switches back to the Lsodes integrator and issues a warming message).
26
27 */
28
29 #include <assert.h>
30 #include <ctype.h>
31 #include <float.h>
32 #include <stdio.h>
33 #include <stdlib.h>
34 #include <string.h>
35 #include <limits.h>
36
37 #include "lexerr.h"
38 #include "list.h"
39 #include "simi.h"
40 #include "siminit.h"
41 #include "simmonte.h"
42 #include "config.h"
43
44 PSTRLEX vrgszlexArgs[ARGS_MAX];
45
46 /* Keyword Map Structure */
47
48 typedef struct tagKM {
49 PSTR szKeyword;
50 int iKWCode; /* Enumeration code of Keyword KM_* */
51 WORD fContext; /* Bit flags of valid context for KW */
52 } KM, *PKM; /* Keyword Map */
53
54 KM vrgkmKeywordMap[] = { /* Global Keyword - code map */
55
56 /* Simulation syntax */
57
58 {"Level", KM_LEVEL, CN_GLOBAL},
59
60 {"Experiment", KM_EXPERIMENT, CN_GLOBAL}, /* Obsolete!! */
61 {"Simulation", KM_EXPERIMENT, CN_GLOBAL},
62
63 {"OutputFile", KM_OUTPUTFILE, CN_GLOBAL}, /* Output file name spec */
64
65 {"MCMC", KM_MCMC, CN_GLOBAL}, /* Metropolis spec */
66
67 {"OptimalDesign", KM_OPTDESIGN, CN_GLOBAL}, /* Optimal design spec */
68 {"MonteCarlo", KM_MONTECARLO, CN_GLOBAL}, /* Monte Carlo spec */
69
70 {"Distrib", KM_MCVARY, CN_GLOBAL}, /* MC Variable mod */
71 {"Likelihood", KM_MCVARY, CN_GLOBAL}, /* special MC Variable */
72 {"Density", KM_MCVARY, CN_GLOBAL}, /* special MC Variable */
73 {"MCVary", KM_MCVARY, CN_GLOBAL}, /* Obsolete!! */
74 {"InvTemperature",KM_TEMPERATURE, CN_GLOBAL}, /* Obsolete!! */
75 {"Perks", KM_TEMPERATURE, CN_GLOBAL}, /* tempered MCMC option */
76
77 {"SetPoints", KM_SETPOINTS, CN_GLOBAL}, /* Forced point runs*/
78
79 {"Integrate", KM_INTEGRATE, CN_GLOBAL | CN_EXPERIMENT},
80 {"Simulate", KM_SIMULATE, CN_GLOBAL | CN_EXPERIMENT}, /* obsolete */
81 {"StartTime", KM_STARTTIME, CN_GLOBAL | CN_EXPERIMENT},
82
83 {"Print", KM_PRINT, CN_EXPERIMENT},
84 {"Prediction", KM_PRINT, CN_EXPERIMENT},
85 {"PrintStep", KM_PRINTSTEP, CN_EXPERIMENT},
86 {"Data", KM_DATA, CN_EXPERIMENT},
87
88 /* If a type is not seen before the first section is found, that section
89 becomes the type. Subsequent statements are ignored. */
90 {"SimType", KM_SIMTYPE, CN_GLOBAL}, /* Optional SimType */
91
92 {"End", KM_END, CN_GLOBAL}, /* Optional End statement */
93 {"END", KM_END, CN_GLOBAL}, /* Optional End statement */
94
95 /* Function arguments */
96
97 {"DefaultSim", KM_DEFAULTSIM, CN_FUNCARG}, /* For SimType() only */
98
99 {"No", KM_NO, CN_FUNCARG}, /* Use YesNoFromLex() */
100 {"Yes", KM_YES, CN_FUNCARG},
101
102 {"Beta", KM_BETA, CN_FUNCARG}, /* Use McvFromLex() */
103 {"Binomial", KM_BINOMIAL, CN_FUNCARG},
104 {"BinomialBeta", KM_BINOMIALBETA, CN_FUNCARG},
105 {"Cauchy", KM_CAUCHY, CN_FUNCARG},
106 {"Chi2", KM_CHI2, CN_FUNCARG},
107 {"Exponential", KM_EXPONENTIAL, CN_FUNCARG},
108 {"Gamma", KM_GGAMMA, CN_FUNCARG},
109 {"GenLogNormal", KM_GENLOGNORMAL, CN_FUNCARG},
110 {"HalfCauchy", KM_HALFCAUCHY, CN_FUNCARG},
111 {"HalfNormal", KM_HALFNORMAL, CN_FUNCARG},
112 {"InvGamma", KM_INVGGAMMA, CN_FUNCARG},
113 {"LogNormal", KM_LOGNORMAL, CN_FUNCARG},
114 {"LogNormal_v", KM_LOGNORMALV, CN_FUNCARG},
115 {"LogUniform", KM_LOGUNIFORM, CN_FUNCARG},
116 {"NegativeBinomial",KM_NEGATIVEBINOM, CN_FUNCARG},
117 {"Normal", KM_NORMAL, CN_FUNCARG},
118 {"Normal_cv", KM_NORMALCV, CN_FUNCARG},
119 {"Normal_v", KM_NORMALV, CN_FUNCARG},
120 {"Piecewise", KM_PIECEWISE, CN_FUNCARG},
121 {"Poisson", KM_POISSON, CN_FUNCARG},
122 {"StudentT", KM_STUDENTT, CN_FUNCARG},
123 {"TruncInvGamma", KM_TRUNCINVGGAMMA, CN_FUNCARG},
124 {"TruncLogNormal", KM_TRUNCLOGNORMAL, CN_FUNCARG},
125 {"TruncLogNormal_v",KM_TRUNCLOGNORMALV, CN_FUNCARG},
126 {"TruncNormal", KM_TRUNCNORMAL, CN_FUNCARG},
127 {"TruncNormal_cv", KM_TRUNCNORMALCV, CN_FUNCARG},
128 {"TruncNormal_v", KM_TRUNCNORMALV, CN_FUNCARG},
129 {"Uniform", KM_UNIFORM, CN_FUNCARG},
130 {"UserSpecifiedLL", KM_USERLL, CN_FUNCARG},
131
132 {"Prediction", KM_PREDICTION, CN_FUNCARG},
133 {"Data", KM_DATA, CN_FUNCARG},
134
135 {"Euler", KM_EULER, CN_FUNCARG},
136 {"Lsodes", KM_LSODES, CN_FUNCARG},
137 {"LSODES", KM_LSODES, CN_FUNCARG},
138 {"Cvodes", KM_CVODES, CN_FUNCARG},
139 {"CVODES", KM_CVODES, CN_FUNCARG},
140
141 {"Forward", KM_FORWARD, CN_FUNCARG},
142 {"Backward", KM_BACKWARD, CN_FUNCARG},
143
144 {"Replace", KM_REPLACE, CN_FUNCARG},
145 {"Add", KM_ADD, CN_FUNCARG},
146 {"Multiply", KM_MULTIPLY, CN_FUNCARG},
147
148 /* Variables names valid in all CN_ */
149
150 {"", 0, CN_ALL} /* End flag */
151
152 }; /* vrgkmKeywordMap[] */
153
154
155 /* ----------------------------------------------------------------------------
156 GetKeywordCode
157
158 Returns the code of the szKeyword given. If the string is not
159 a valid keyword or abbreviation, returns 0.
160
161 If pfContext is non-NULL, contexts in which the code is valid is
162 set too.
163 */
GetKeywordCode(PSTR szKeyword,PINT pfContext)164 int GetKeywordCode (PSTR szKeyword, PINT pfContext)
165 {
166 PKM pkm = &vrgkmKeywordMap[0];
167
168 while (*pkm->szKeyword && MyStrcmp(szKeyword, pkm->szKeyword))
169 pkm++;
170
171 if (pfContext)
172 *pfContext = pkm->fContext; /* Set iContext flag */
173
174 return(pkm->iKWCode); /* Return Keyword Code or 0 */
175
176 } /* GetKeywordCode */
177
178
179 /* ----------------------------------------------------------------------------
180 GetKeywordCode_in_context
181
182 Returns the code of the szKeyword given, in a given context.
183 If the string is not a valid keyword or abbreviation in context, returns 0.
184
185 */
GetKeywordCode_in_context(PSTR szKeyword,WORD fContext)186 int GetKeywordCode_in_context (PSTR szKeyword, WORD fContext)
187 {
188 PKM pkm = &vrgkmKeywordMap[0];
189
190 while (*pkm->szKeyword && !((pkm->fContext == fContext) &&
191 !MyStrcmp(szKeyword, pkm->szKeyword)))
192 pkm++;
193
194 return(pkm->iKWCode); /* Return Keyword Code or 0 */
195
196 } /* GetKeywordCode_in_context */
197
198
199 /* ----------------------------------------------------------------------------
200 GetKeyword
201
202 Returns the first string of the KM_ keyword map code given. If the
203 code is not a valid keyword code, returns NULL.
204 */
GetKeyword(int iKWCode)205 PSTR GetKeyword (int iKWCode)
206 {
207 PKM pkm = &vrgkmKeywordMap[0];
208
209 while (*pkm->szKeyword && iKWCode != pkm->iKWCode)
210 pkm++;
211
212 return(pkm->szKeyword); /* Return Keyword Code or 0 */
213
214 } /* GetKeyword */
215
216
217 /* ----------------------------------------------------------------------------
218 YesNoFromLex
219
220 Converts an string input argument into a Boolean,
221 Yes being TRUE and No being FALSE. Also, a numeric argument
222 is converted to Yes if it is non-zero.
223 */
YesNoFromLex(PSTR szLex)224 BOOL YesNoFromLex (PSTR szLex)
225 {
226 int ikwcode = GetKeywordCode(szLex, NULL);
227 BOOL bReturn;
228
229 bReturn = (!isalpha(szLex[0]) ? atoi(szLex)
230 : ikwcode == KM_YES ? TRUE
231 : ikwcode == KM_NO ? FALSE
232 : FALSE);
233
234 return bReturn;
235
236 } /* YesNoFromLex */
237
238
239 /* ----------------------------------------------------------------------------
240 ImFromLex
241
242 Converts an string input argument into the correct IM_
243 integration method.
244 */
ImFromLex(PSTR szLex)245 long ImFromLex (PSTR szLex)
246 {
247 int ikwcode = GetKeywordCode(szLex, NULL);
248 long lReturn;
249
250 lReturn = (!isalpha(szLex[0]) ? atoi(szLex)
251 : ikwcode == KM_LSODES ? IAL_LSODES
252 : ikwcode == KM_CVODES ? IAL_CVODES
253 : ikwcode == KM_EULER ? IAL_EULER
254 : 0);
255
256 return(lReturn);
257
258 } /* ImFromLex */
259
260
261 /* ----------------------------------------------------------------------------
262 McvFromLex
263
264 Converts a string input argument into the correct MCV_
265 Monte Carlo variation distribution type.
266 */
McvFromLex(PSTR szLex)267 int McvFromLex (PSTR szLex)
268 {
269 int ikwcode = GetKeywordCode(szLex, NULL);
270 int iReturn;
271
272 iReturn = (ikwcode == KM_UNIFORM ? MCV_UNIFORM
273 : ikwcode == KM_LOGUNIFORM ? MCV_LOGUNIFORM
274 : ikwcode == KM_BETA ? MCV_BETA
275 : ikwcode == KM_NORMAL ? MCV_NORMAL
276 : ikwcode == KM_HALFNORMAL ? MCV_HALFNORMAL
277 : ikwcode == KM_LOGNORMAL ? MCV_LOGNORMAL
278 : ikwcode == KM_TRUNCNORMAL ? MCV_TRUNCNORMAL
279 : ikwcode == KM_TRUNCLOGNORMAL ? MCV_TRUNCLOGNORMAL
280 : ikwcode == KM_CHI2 ? MCV_CHI2
281 : ikwcode == KM_BINOMIAL ? MCV_BINOMIAL
282 : ikwcode == KM_PIECEWISE ? MCV_PIECEWISE
283 : ikwcode == KM_EXPONENTIAL ? MCV_EXPONENTIAL
284 : ikwcode == KM_GGAMMA ? MCV_GGAMMA
285 : ikwcode == KM_POISSON ? MCV_POISSON
286 : ikwcode == KM_INVGGAMMA ? MCV_INVGGAMMA
287 : ikwcode == KM_TRUNCINVGGAMMA ? MCV_TRUNCINVGGAMMA
288 : ikwcode == KM_NORMALV ? MCV_NORMALV
289 : ikwcode == KM_NORMALCV ? MCV_NORMALCV
290 : ikwcode == KM_LOGNORMALV ? MCV_LOGNORMALV
291 : ikwcode == KM_TRUNCNORMALV ? MCV_TRUNCNORMALV
292 : ikwcode == KM_TRUNCNORMALCV ? MCV_TRUNCNORMALCV
293 : ikwcode == KM_TRUNCLOGNORMALV? MCV_TRUNCLOGNORMALV
294 : ikwcode == KM_BINOMIALBETA ? MCV_BINOMIALBETA
295 : ikwcode == KM_GENLOGNORMAL ? MCV_GENLOGNORMAL
296 : ikwcode == KM_STUDENTT ? MCV_STUDENTT
297 : ikwcode == KM_CAUCHY ? MCV_CAUCHY
298 : ikwcode == KM_HALFCAUCHY ? MCV_HALFCAUCHY
299 : ikwcode == KM_USERLL ? MCV_USERLL
300 : ikwcode == KM_NEGATIVEBINOM ? MCV_NEGATIVEBINOM
301 : (-1));
302
303 return iReturn;
304
305 } /* McvFromLex */
306
307
308 /* ----------------------------------------------------------------------------
309 GetTerminator
310
311 Tries to read a statement terminator. Reports Errors.
312 */
GetTerminator(PINPUTBUF pibIn,PSTR szLex)313 int GetTerminator (PINPUTBUF pibIn, PSTR szLex)
314 {
315 int iErr;
316
317 if ((iErr = !GetPunct(pibIn, szLex, CH_STMTTERM))) {
318 szLex[1] = CH_STMTTERM;
319 ReportError(pibIn, RE_EXPECTED, szLex, NULL);
320 }
321
322 return(iErr);
323
324 } /* GetTerminator */
325
326
327 /* ----------------------------------------------------------------------------
328 GetSimTYpe
329
330 Get the type of simulation to perform.
331 */
GetSimType(PINPUTBUF pibIn)332 BOOL GetSimType (PINPUTBUF pibIn)
333 {
334 #define NAT_ARGS 1 /* the type */
335
336 static int vrgiAtArgTypes[NAT_ARGS] = {LX_IDENTIFIER};
337
338 PANALYSIS panal = (PANALYSIS) pibIn->pInfo;
339
340 int iAT = AT_DEFAULTSIM;
341 int iKwCode = 0;
342 BOOL bErr=!GetFuncArgs(pibIn, NAT_ARGS, vrgiAtArgTypes, vrgszlexArgs[0]);
343
344 if (!bErr) {
345 iKwCode = GetKeywordCode(vrgszlexArgs[0], NULL);
346 switch (iKwCode) {
347
348 case KM_MONTECARLO:
349 iAT = AT_MONTECARLO;
350 break;
351
352 case KM_SETPOINTS:
353 iAT = AT_SETPOINTS;
354 break;
355
356 case KM_MCMC:
357 iAT = AT_MCMC;
358 break;
359
360 case KM_OPTDESIGN:
361 iAT = AT_OPTDESIGN;
362 break;
363
364 case KM_DEFAULTSIM:
365 iAT = AT_DEFAULTSIM;
366 break;
367
368 default:
369 ReportError(pibIn, RE_SPECERR | RE_FATAL, "Unknown SimType ",
370 vrgszlexArgs[0]);
371 break;
372 } /* switch */
373 } /* if */
374 else
375 printf("Syntax: %s (Normal | MonteCarlo | SetPoints | MCMC)\n"
376 " -- if not specified, the first spec section will be used.\n\n",
377 GetKeyword(KM_SIMTYPE));
378
379 if (!bErr)
380 panal->iType = iAT;
381
382 return(bErr);
383
384 } /* GetSimType */
385
386
387 /* ----------------------------------------------------------------------------
388 GetPerks
389
390 In case of tempered MCMC, read in a user-defined inverse
391 temperature schedule.
392 Return TRUE if the structure is defined, FALSE on error. The
393 command syntax includes the number of inverse temperatures then the list of
394 them (numbers in [0,inf] interval).
395 */
396
GetPerks(PINPUTBUF pibIn,PSTR szLex,PGIBBSDATA pgd)397 BOOL GetPerks (PINPUTBUF pibIn, PSTR szLex, PGIBBSDATA pgd)
398 {
399 int iType;
400 int i;
401 BOOL bOK = TRUE;
402 BOOL bErr = FALSE; /* Return value flags error condition */
403
404 if ((bErr = EGetPunct(pibIn, szLex, CH_LPAREN)))
405 goto Exit_GetPerks;
406
407 if ((bErr = ENextLex(pibIn, szLex, LX_INTEGER)))
408 goto Exit_GetPerks;
409
410 pgd->nPerks = atoi(szLex);
411
412 if ((bErr = (pgd->nPerks <= 0))) {
413 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, "positive-integer", szLex);
414 goto Exit_GetPerks;
415 }
416
417 pgd->endT = pgd->nPerks - 1;
418
419 /* allocate inverse temperature array */
420 if (!(pgd->rgdPerks = InitdVector(pgd->nPerks)))
421 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetPerks", NULL);
422
423 /* allocate working arrays */
424 if (!(pgd->rgdlnPi = InitdVector(pgd->nPerks)) ||
425 !(pgd->rglCount = InitlVector(pgd->nPerks)))
426 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetPerks", NULL);
427
428 /* try to get temperatures list: */
429
430 for (i = 0; i < pgd->nPerks && bOK; i++) {
431
432 /* get comma */
433 if (!(bOK = GetOptPunct(pibIn, szLex, ','))) {
434 szLex[0] = *pibIn->pbufCur++;
435 szLex[1] = '\0';
436 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, ",", szLex);
437 } /* if */
438
439 /* get floating point number */
440 NextLex(pibIn, szLex, &iType);
441 if (!(bOK &= (iType & LX_NUMBER) > 0))
442 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, "number", szLex);
443 /* assign */
444 pgd->rgdPerks[i] = atof(szLex);
445 /* initialize */
446 pgd->rgdlnPi[i] = 0; /* perks */
447 pgd->rglCount[i] = 0; /* temperature counts */
448 if (pgd->rgdPerks[i] < 0)
449 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL,
450 "positive inverse temperature", szLex);
451 if ((i > 0) && (pgd->rgdPerks[i] <= pgd->rgdPerks[i-1]))
452 ReportError(pibIn, RE_SPECERR | RE_FATAL,
453 "Inverse temperatures out of order", NULL);
454 } /* for */
455
456 bErr = EGetPunct(pibIn, szLex, CH_RPAREN);
457
458 Exit_GetPerks:
459
460 if (bErr)
461 printf("Syntax: Inverse temperatures (nPerks, "
462 "<n increasing inverse temperature values >= 0>)\n\n");
463
464 return(!bErr);
465
466 } /* GetPerks */
467
468
469 /* ----------------------------------------------------------------------------
470 GetLsodesOptions
471 */
GetLsodesOptions(PINPUTBUF pibIn,PSTR szLex,PINTSPEC pis)472 BOOL GetLsodesOptions (PINPUTBUF pibIn, PSTR szLex, PINTSPEC pis)
473 {
474 BOOL bErr = FALSE;
475
476 if (!(GetPunct(pibIn, szLex, ',')))
477 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, ",", szLex);
478
479 if ((bErr = ENextLex(pibIn, szLex, LX_NUMBER)))
480 goto Exit_GetLsodes;
481 pis->dRtol = atof(szLex);
482
483 if (!(GetPunct(pibIn, szLex, ',')))
484 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, ",", szLex);
485
486 if ((bErr = ENextLex(pibIn, szLex, LX_NUMBER)))
487 goto Exit_GetLsodes;
488 pis->dAtol = atof(szLex);
489
490 if (!(GetPunct(pibIn, szLex, ',')))
491 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, ",", szLex);
492
493 if ((bErr = ENextLex(pibIn, szLex, LX_INTEGER)))
494 goto Exit_GetLsodes;
495 pis->iMf = atoi(szLex);
496
497 /* translate input iMf to original lsodes mf codes */
498 switch (pis->iMf) {
499
500 case 0:
501 pis->iMf = 10;
502 break;
503
504 case 1:
505 pis->iMf = 222;
506 break;
507
508 case 2:
509 pis->iMf = 121;
510 break;
511
512 default:
513 printf("Error: method flag must be 0, 1 or 2 for Lsodes - ");
514 printf("Exiting\n\n");
515 exit(0);
516 break;
517
518 } /* switch */
519
520 pis->iDSFlag = 1;
521
522 Exit_GetLsodes:
523
524 if (bErr) {
525 printf("Lsodes options are: relative tolerance, absolute tolerance, "
526 "method.\n\n");
527 exit(0);
528 }
529
530 return(bErr);
531
532 } /* GetLsodesOptions */
533
534
535 /* ----------------------------------------------------------------------------
536 GetCvodesOptions
537 */
GetCvodesOptions(PINPUTBUF pibIn,PSTR szLex,PINTSPEC pis)538 BOOL GetCvodesOptions (PINPUTBUF pibIn, PSTR szLex, PINTSPEC pis)
539 {
540 BOOL bErr = FALSE;
541
542 if (!(GetPunct(pibIn, szLex, ',')))
543 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, ",", szLex);
544
545 if ((bErr = ENextLex(pibIn, szLex, LX_NUMBER)))
546 goto Exit_GetCvodes;
547 pis->dRtol = atof(szLex);
548
549 if (!(GetPunct(pibIn, szLex, ',')))
550 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, ",", szLex);
551
552 if ((bErr = ENextLex(pibIn, szLex, LX_NUMBER)))
553 goto Exit_GetCvodes;
554 pis->dAtol = atof(szLex);
555
556 if (!(GetPunct(pibIn, szLex, ',')))
557 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, ",", szLex);
558
559 if ((bErr = ENextLex(pibIn, szLex, LX_INTEGER)))
560 goto Exit_GetCvodes;
561 pis->maxsteps = atoi(szLex);
562
563 if (!(GetPunct(pibIn, szLex, ',')))
564 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, ",", szLex);
565
566 if ((bErr = ENextLex(pibIn, szLex, LX_INTEGER)))
567 goto Exit_GetCvodes;
568 pis->maxnef = atoi(szLex);
569
570 if (!(GetPunct(pibIn, szLex, ',')))
571 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, ",", szLex);
572
573 if ((bErr = ENextLex(pibIn, szLex, LX_INTEGER)))
574 goto Exit_GetCvodes;
575 pis->maxcor = atoi(szLex);
576
577 if (!(GetPunct(pibIn, szLex, ',')))
578 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, ",", szLex);
579
580 if ((bErr = ENextLex(pibIn, szLex, LX_INTEGER)))
581 goto Exit_GetCvodes;
582 pis->maxncf = atoi(szLex);
583
584 if (!(GetPunct(pibIn, szLex, ',')))
585 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, ",", szLex);
586
587 if ((bErr = ENextLex(pibIn, szLex, LX_NUMBER)))
588 goto Exit_GetCvodes;
589 pis->nlscoef = atof(szLex);
590
591 Exit_GetCvodes:
592
593 if (bErr) {
594 printf("Cvodes options are: relative tolerance, absolute tolerance, "
595 "maxsteps, maxnef, maxcor, maxncf, nlscoef.\n\n");
596 exit(0);
597 }
598
599 return(bErr);
600
601 } /* GetCvodesOptions */
602
603
604 /* ----------------------------------------------------------------------------
605 GetEulerOptions
606 */
GetEulerOptions(PINPUTBUF pibIn,PSTR szLex,PINTSPEC pis)607 BOOL GetEulerOptions (PINPUTBUF pibIn, PSTR szLex, PINTSPEC pis)
608 {
609 BOOL bErr = FALSE;
610
611 if (!(GetPunct(pibIn, szLex, ',')))
612 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, ",", szLex);
613
614 if ((bErr = ENextLex(pibIn, szLex, LX_NUMBER)))
615 goto Exit_GetEuler;
616 pis->dTStep = atof(szLex);
617
618 if (pis->dTStep <= 0) {
619 printf("Error: Time step specified is null or negative - Exiting\n\n");
620 exit(0);
621 }
622
623 Exit_GetEuler:
624
625 if (bErr) {
626 printf("Euler has one option: time-step.\n\n");
627 exit(0);
628 }
629
630 return(bErr);
631
632 } /* GetEulerOptions */
633
634
635 /* ----------------------------------------------------------------------------
636 GetIntegrate
637 */
GetIntegrate(PINPUTBUF pibIn,PSTR szLex,PINTSPEC pis)638 BOOL GetIntegrate (PINPUTBUF pibIn, PSTR szLex, PINTSPEC pis)
639 {
640 BOOL bErr = FALSE;
641
642 if ((bErr = EGetPunct(pibIn, szLex, CH_LPAREN)))
643 goto Exit_GetIntegrate;
644
645 if ((bErr = ENextLex(pibIn, szLex, LX_IDENTIFIER)))
646 goto Exit_GetIntegrate;
647
648 pis->iAlgo = ImFromLex(szLex);
649
650 switch (pis->iAlgo) {
651 case IAL_LSODES:
652 GetLsodesOptions(pibIn, szLex, pis);
653 break;
654
655 case IAL_CVODES:
656 #ifndef HAVE_LIBSUNDIALS_CVODES
657 #ifndef HAVE_LIBSUNDIALS_NVECSERIAL
658 /* switch to lsodes if cvodes is not installed */
659 printf("Warning: Cvodes libraries are not available -\n"
660 " Switching to Lsodes with default options\n\n");
661 pis->iAlgo = IAL_DEFAULT;
662 goto Exit_GetIntegrate;
663 #endif
664 #endif
665 GetCvodesOptions(pibIn, szLex, pis);
666 break;
667
668 case IAL_EULER:
669 GetEulerOptions(pibIn, szLex, pis);
670 break;
671
672 default:
673 printf("Error: Unknown integration method: %s - Exiting\n\n",
674 vrgszlexArgs[0]);
675 bErr = TRUE;
676 goto Exit_GetIntegrate;
677 break;
678 }
679
680 bErr = EGetPunct(pibIn, szLex, CH_RPAREN);
681
682 Exit_GetIntegrate:
683
684 if (bErr) {
685 printf("Syntax: %s([Lsodes | Cvodes | Euler], [OPTIONS]);\n\n",
686 GetKeyword(KM_INTEGRATE));
687 exit(0);
688 }
689
690 return(!bErr);
691
692 } /* GetIntegrate */
693
694
695 /* ----------------------------------------------------------------------------
696 OneDToArray
697
698 Copies one double from the list to the newly formed array.
699 Increments the info pointer which is the pointer into the array.
700 */
OneDToArray(PVOID pData,PVOID pInfo)701 int OneDToArray (PVOID pData, PVOID pInfo)
702 {
703 PDOUBLE *ppdArrayVal = (PDOUBLE *) pInfo;
704
705 *(*ppdArrayVal)++ = *(PDOUBLE) pData;
706
707 return 0;
708
709 } /* OneDToArray */
710
711
712 /* ----------------------------------------------------------------------------
713 DListToArray
714
715 Converts a list a doubles to an array of doubles. *pcDouble is
716 the count of doubles in the array, and *ppDouble is the array
717 pointer.
718 */
DListToArray(PLIST plist,PLONG pcDouble,PDOUBLE * ppDouble)719 void DListToArray (PLIST plist, PLONG pcDouble, PDOUBLE *ppDouble)
720 {
721 PDOUBLE pdTmp; /* Temp pointer to get incremented */
722
723 *pcDouble = ListLength(plist);
724
725 if (!(pdTmp = *ppDouble = InitdVector(*pcDouble)))
726 ReportError(NULL, RE_OUTOFMEM | RE_FATAL, "DListToArray", NULL);
727
728 ForAllList(plist, &OneDToArray, (PVOID) &pdTmp);
729
730 } /* DListToArray */
731
732
733 /* ----------------------------------------------------------------------------
734 GetListOfTimes
735
736 Reads an arbitrary length list of times and closing parenthesis.
737 Defines the count and array of times in the PRINTREC structure.
738 */
GetListOfTimes(PINPUTBUF pibIn,int nRecs,PPRINTREC * ppr,PSTR szLex)739 BOOL GetListOfTimes (PINPUTBUF pibIn, int nRecs, PPRINTREC *ppr, PSTR szLex)
740 {
741 PLIST plistTimes = InitList();
742 PDOUBLE pdTmp;
743 int iNLI, i, j;
744 BOOL bErr;
745
746 do {
747 if ( !(pdTmp = InitdVector(1)))
748 ReportError(NULL, RE_OUTOFMEM | RE_FATAL, "GetListOfTimes", NULL);
749
750 *pdTmp = atof(szLex);
751 QueueListItem(plistTimes, (PVOID) pdTmp);
752 } while ((iNLI = NextListItem(pibIn, szLex, LX_NUMBER, 1, CH_RPAREN)) > 0);
753
754 if (!iNLI) /* List terminator */
755 bErr = EGetPunct(pibIn, szLex, CH_RPAREN) || !ListLength(plistTimes);
756 else {
757 bErr = TRUE;
758 ReportError(pibIn, RE_LEXEXPECTED, "number", szLex);
759 }
760
761 if (!bErr)
762 for (i = 0; i < nRecs; ++i)
763 DListToArray(plistTimes, &ppr[i]->cTimes, &ppr[i]->pdTimes);
764
765 FreeList(&plistTimes, NULL, TRUE); /* Free list and cells */
766
767 for (i = 1; i < ppr[0]->cTimes && !bErr; i++) /* Verify Times */
768 if ((bErr = (*(ppr[0]->pdTimes+i) <= *(ppr[0]->pdTimes+i-1)))) {
769 for (j = 0; j < nRecs; ++j)
770 free(ppr[j]->pdTimes);
771 ReportError(pibIn, RE_SPECERR | RE_FATAL, "Times out of order", NULL);
772 } /* if */
773
774 return(bErr);
775
776 } /* GetListOfTimes */
777
778
779 /* ----------------------------------------------------------------------------
780 GetListOfData
781
782 Reads an arbitrary length list of data and closing parenthesis.
783 Defines the count and array of data in the DATAREC structure.
784 */
GetListOfData(PINPUTBUF pibIn,PDATAREC pda,PSTR szLex)785 BOOL GetListOfData (PINPUTBUF pibIn, PDATAREC pda, PSTR szLex)
786 {
787 PLIST plistData = InitList();
788 PDOUBLE pdTmp;
789 int iNLI;
790 BOOL bErr;
791
792 while ((iNLI = NextListItem(pibIn, szLex, LX_NUMBER, 1, CH_RPAREN))
793 > 0) {
794 if ( !(pdTmp = InitdVector(1)))
795 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetListOfData", NULL);
796
797 *pdTmp = atof(szLex);
798 QueueListItem(plistData, (PVOID) pdTmp);
799
800 } /* while */
801
802 if (!iNLI) /* List terminator */
803 bErr = EGetPunct(pibIn, szLex, CH_RPAREN)
804 || !ListLength(plistData);
805 else {
806 bErr = TRUE;
807 ReportError(pibIn, RE_LEXEXPECTED, "number", szLex);
808 } /* else */
809
810 if (!bErr) DListToArray(plistData, &pda->cData, &pda->pdData);
811
812 FreeList(&plistData, NULL, TRUE); /* Free list and cells */
813 return(bErr);
814
815 } /* GetListOfData */
816
817
818 /* ----------------------------------------------------------------------------
819 GetPrint
820
821 Gets the arguments to a Print() statement. Put them in
822 a list plistPrintRecs of PRINTREC structures
823 */
824 BOOL bGavePrintUsage = FALSE; /* prevent multiple diagnostics */
825
GetPrint(PINPUTBUF pibIn,PSTR szLex,POUTSPEC pos)826 BOOL GetPrint (PINPUTBUF pibIn, PSTR szLex, POUTSPEC pos)
827 {
828 PPRINTREC pprintrec[MAX_PRINT_VARS];
829 BOOL bErr = FALSE;
830 HVAR hvar;
831 int nVars = 0, n, iLex;
832 long i, iLB, iUB;
833 PSTRLEX szTmp;
834
835 /* get the opening parenthesis and list of variables to print */
836 if (!(bErr = EGetPunct(pibIn, szLex, CH_LPAREN))) {
837 for (;;) {
838 NextLex(pibIn, szLex, &iLex);
839
840 if (iLex != LX_IDENTIFIER)
841 break; /* hitting the list of times, stop reading variable names */
842
843 /* check if it is a scalar or an array and act accordingly */
844 iLB = iUB = -1;
845 if (GetPunct(pibIn, szTmp, '[')) /* array found, read bounds */
846 GetArrayBounds(pibIn, &iLB, &iUB);
847
848 if (iUB == -1) { /* scalar */
849
850 if (nVars == MAX_PRINT_VARS)
851 ReportError(pibIn, RE_TOOMANYPVARS | RE_FATAL, "GetPrint", NULL);
852
853 if ((bErr = !(hvar = GetVarHandle(szLex))))
854 ReportError(pibIn, RE_UNDEFINED | RE_FATAL, szLex, NULL);
855 else {
856 if ( !(pprintrec[nVars] = (PPRINTREC) malloc(sizeof(PRINTREC))) ||
857 !(pprintrec[nVars]->szOutputName =
858 (PSTR) malloc(MyStrlen(szLex)+1)))
859 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetPrint", NULL);
860 MyStrcpy(pprintrec[nVars]->szOutputName, szLex);
861 pprintrec[nVars]->hvar = hvar;
862 assert(pprintrec[nVars]);
863 ++nVars;
864 }
865 GetOptPunct(pibIn, szLex, ',');
866 }
867 else { /* array */
868
869 for (i = iLB; i < iUB; i++) {
870 sprintf(szTmp, "%s_%ld", szLex, i); /* create names */
871
872 if (nVars == MAX_PRINT_VARS)
873 ReportError(pibIn, RE_TOOMANYPVARS | RE_FATAL, "GetPrint", NULL);
874
875 if ((bErr = !(hvar = GetVarHandle(szTmp))))
876 ReportError(pibIn, RE_UNDEFINED | RE_FATAL, szTmp, NULL);
877 else {
878 if ( !(pprintrec[nVars] = (PPRINTREC) malloc(sizeof(PRINTREC))) ||
879 !(pprintrec[nVars]->szOutputName =
880 (PSTR) malloc(strlen(szTmp)+1))) {
881 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetPrint", NULL);
882 }
883 strcpy(pprintrec[nVars]->szOutputName, szTmp);
884 pprintrec[nVars]->hvar = hvar;
885 assert(pprintrec[nVars]);
886 ++nVars;
887 }
888
889 if (i == (iUB - 1))
890 GetOptPunct(pibIn, szTmp, ',');
891 } /* end for i */
892 } /* end else */
893 } /* end for ever */
894
895 /* check the number of variables read in */
896 if (nVars < 1)
897 ReportError(pibIn, RE_LEXEXPECTED, "identifier", szLex);
898
899 /* get list of output times */
900 bErr = GetListOfTimes(pibIn, nVars, pprintrec, szLex);
901
902 if (bErr) {
903 for (n = 0; n < nVars; ++n) {
904 free(pprintrec[n]->szOutputName);
905 free(pprintrec[n]);
906 }
907 }
908 else
909 for (n = 0; n < nVars; ++n)
910 QueueListItem(pos->plistPrintRecs, (PVOID) pprintrec[n]);
911 } /* if */
912
913 if (!bErr) bErr = GetTerminator(pibIn, szLex);
914 else {
915 if (!bGavePrintUsage) {
916 printf("Syntax: %s (<Identifiers>, Time1, Time2, ...)\n\n",
917 GetKeyword(KM_PRINT));
918 bGavePrintUsage = TRUE;
919 }
920 }
921
922 return(bErr);
923
924 } /* GetPrint */
925
926
927 /* ----------------------------------------------------------------------------
928 GetPrintStep
929
930 Gets the arguments to a PrintStep() statement. They are: a list of
931 identifiers, a start time, an end time, a time step.
932 If the time period is not congruent with the time step the last step
933 will be shortened.
934 */
935 BOOL bGavePrintStepUsage = FALSE; /* prevent multiple diagnostics */
936
GetPrintStep(PINPUTBUF pibIn,PSTR szLex,POUTSPEC pos)937 BOOL GetPrintStep (PINPUTBUF pibIn, PSTR szLex, POUTSPEC pos)
938 {
939 PPRINTREC pprintrec[MAX_PRINT_VARS];
940 BOOL bErr = FALSE, bOK = TRUE;
941 HVAR hvar = 0;
942 int nVars = 0, n, iLex;
943 long i, iLB, iUB;
944 double dStart = 0, dEnd = 0, dStep = 0, dTmp;
945 PSTRLEX szTmp;
946
947 /* get the opening parenthesis and list of variables to print */
948 bErr = EGetPunct(pibIn, szLex, CH_LPAREN);
949 if (bErr)
950 goto Exit_GetPrintStep;
951
952 for (;;) {
953 NextLex(pibIn, szLex, &iLex);
954
955 if (iLex != LX_IDENTIFIER)
956 break; /* hitting the list of times, stop reading variable names */
957
958 /* check if it is a scalar or an array and act accordingly */
959 iLB = iUB = -1;
960 if (GetPunct(pibIn, szTmp, '[')) /* array found, read bounds */
961 GetArrayBounds(pibIn, &iLB, &iUB);
962
963 if (iUB == -1) { /* scalar */
964
965 if (nVars == MAX_PRINT_VARS)
966 ReportError(pibIn, RE_TOOMANYPVARS | RE_FATAL, "GetPrint", NULL);
967
968 if ((bErr = !(hvar = GetVarHandle(szLex))))
969 ReportError(pibIn, RE_UNDEFINED | RE_FATAL, szLex, NULL);
970 else {
971 if ( !(pprintrec[nVars] = (PPRINTREC) malloc(sizeof(PRINTREC))) ||
972 !(pprintrec[nVars]->szOutputName =
973 (PSTR) malloc(MyStrlen(szLex)+1)))
974 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetPrint", NULL);
975 MyStrcpy(pprintrec[nVars]->szOutputName, szLex);
976 pprintrec[nVars]->hvar = hvar;
977 assert(pprintrec[nVars]);
978 ++nVars;
979 }
980 GetOptPunct(pibIn, szLex, ',');
981 }
982 else { /* array */
983
984 for (i = iLB; i < iUB; i++) {
985 sprintf(szTmp, "%s_%ld", szLex, i); /* create names */
986
987 if (nVars == MAX_PRINT_VARS)
988 ReportError(pibIn, RE_TOOMANYPVARS | RE_FATAL, "GetPrint", NULL);
989
990 if ((bErr = !(hvar = GetVarHandle(szTmp))))
991 ReportError(pibIn, RE_UNDEFINED | RE_FATAL, szTmp, NULL);
992 else {
993 if ( !(pprintrec[nVars] = (PPRINTREC) malloc(sizeof(PRINTREC))) ||
994 !(pprintrec[nVars]->szOutputName =
995 (PSTR) malloc(strlen(szTmp)+1))) {
996 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetPrint", NULL);
997 }
998 strcpy(pprintrec[nVars]->szOutputName, szTmp);
999 pprintrec[nVars]->hvar = hvar;
1000 assert(pprintrec[nVars]);
1001 ++nVars;
1002 }
1003
1004 if (i == (iUB - 1))
1005 GetOptPunct(pibIn, szTmp, ',');
1006 } /* end for i*/
1007 } /* end else */
1008 } /* end for ever */
1009
1010 /* check the number of variables read in */
1011 if (nVars < 1) {
1012 ReportError(pibIn, RE_LEXEXPECTED, "identifier", szLex);
1013 bErr = TRUE;
1014 goto Exit_GetPrintStep;
1015 }
1016
1017 /* the starting output time has already been read */
1018 dStart = atof(szLex);
1019 if (!(bOK = GetPunct(pibIn, szLex, ',')))
1020 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, ",", szLex);
1021
1022 /* get ending output time */
1023 if ((bErr = ENextLex(pibIn, szLex, LX_NUMBER)))
1024 goto Exit_GetPrintStep;
1025 dEnd = atof(szLex);
1026 if (!(bOK = GetPunct(pibIn, szLex, ',')))
1027 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, ",", szLex);
1028
1029 /* get output time step */
1030 if ((bErr = ENextLex(pibIn, szLex, LX_NUMBER)))
1031 goto Exit_GetPrintStep;
1032 dStep = atof(szLex);
1033
1034 /* get closing parenthesis */
1035 if (!(bOK = GetPunct(pibIn, szLex, CH_RPAREN)))
1036 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, ")", szLex);
1037
1038 /* check times for consistency */
1039 if ((bErr = (dEnd <= dStart))) {
1040 ReportError(pibIn, RE_SPECERR, "End_time must be > Start_time", NULL);
1041 goto Exit_GetPrintStep;
1042 }
1043 else if ((bErr = (dStep > (dEnd - dStart)))) {
1044 ReportError(pibIn, RE_SPECERR, "Time_step too large", NULL);
1045 goto Exit_GetPrintStep;
1046 }
1047
1048 /* fix an unlikely overflow */
1049 dTmp = 1 + ceil((dEnd - dStart) / dStep);
1050 for (n = 0; n < nVars; ++n) {
1051 if (dTmp < LONG_MAX)
1052 pprintrec[n]->cTimes = (long) dTmp;
1053 else
1054 pprintrec[n]->cTimes = LONG_MAX;
1055 }
1056
1057 for (n = 0; n < nVars; ++n) {
1058 if ( !(pprintrec[n]->pdTimes = InitdVector(pprintrec[n]->cTimes)))
1059 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetPrintStep", NULL);
1060
1061 for (i = 0; i < pprintrec[n]->cTimes - 1; i++)
1062 pprintrec[n]->pdTimes[i] = dStart + (i * dStep);
1063
1064 pprintrec[n]->pdTimes[pprintrec[n]->cTimes - 1] = dEnd;
1065
1066 QueueListItem(pos->plistPrintRecs, (PVOID) pprintrec[n]);
1067 }
1068
1069 Exit_GetPrintStep:
1070
1071 if (bErr)
1072 if (!bGavePrintStepUsage) {
1073 printf("Syntax: %s (<Identifiers>, Start_time, End_time, Time_step)\n\n",
1074 GetKeyword(KM_PRINTSTEP));
1075 bGavePrintStepUsage = TRUE;
1076 }
1077
1078 return(bErr);
1079
1080 } /* GetPrintStep */
1081
1082
1083 /* ----------------------------------------------------------------------------
1084 GetData
1085
1086 Gets the arguments to a Data() statement
1087 */
1088 BOOL bGaveDataUsage = FALSE; /* prevent multiple diagnostics */
1089
GetData(PINPUTBUF pibIn,PSTR szLex,POUTSPEC pos)1090 BOOL GetData (PINPUTBUF pibIn, PSTR szLex, POUTSPEC pos)
1091 {
1092 PDATAREC pdatarec;
1093 BOOL bErr = FALSE;
1094 HVAR hvar;
1095
1096 if (!(bErr = EGetPunct(pibIn, szLex, CH_LPAREN))) {
1097 if (!(bErr = ENextLex(pibIn, szLex, LX_IDENTIFIER))) {
1098
1099 if ((bErr = !(hvar = GetVarHandle(szLex))))
1100 ReportError(pibIn, RE_UNDEFINED, szLex, NULL);
1101
1102 else {
1103 if ( !(pdatarec = (PDATAREC) malloc(sizeof(DATAREC))))
1104 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetData", NULL);
1105
1106 if ( !(pdatarec->szDataName = (PSTR) malloc(MyStrlen(szLex)+1)))
1107 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetData", NULL);
1108
1109 MyStrcpy(pdatarec->szDataName, szLex);
1110 assert(pdatarec);
1111
1112 pdatarec->hvar = hvar;
1113
1114 bErr = GetListOfData(pibIn, pdatarec, szLex);
1115
1116 if (bErr) {
1117 free(pdatarec->szDataName);
1118 free(pdatarec);
1119 } /* if */
1120 else
1121 QueueListItem(pos->plistDataRecs, (PVOID) pdatarec);
1122 } /* else */
1123 } /* if */
1124 } /* if */
1125
1126 if (!bErr) bErr = GetTerminator(pibIn, szLex);
1127 else {
1128 if (!bGaveDataUsage) {
1129 printf ("Syntax: %s (identifier, Time1, Time2, ...)\n\n",
1130 GetKeyword(KM_DATA));
1131 bGaveDataUsage = TRUE;
1132 } /* if */
1133 } /* else */
1134
1135 return(bErr);
1136
1137 } /* GetData */
1138
1139
1140 /* ----------------------------------------------------------------------------
1141 GetStringArg
1142
1143 tries to read a string argument from pibIn and assign it to
1144 *pszArg. If pszArg is NULL, the argument is read, but no
1145 assigment is made. If pszArg is not NULL, space is allocated for
1146 argument read. szLex is a workspace. If bDelim is TRUE, a
1147 delimiter is skipped in the input buffer.
1148
1149 The return value is TRUE for error. Errors are reported.
1150 */
GetStringArg(PINPUTBUF pibIn,PSTR * pszArg,PSTR szLex,BOOL bDelim)1151 BOOL GetStringArg (PINPUTBUF pibIn, PSTR *pszArg, PSTR szLex, BOOL bDelim)
1152 {
1153 BOOL bErr;
1154
1155 assert(szLex); /* Workspace must be given */
1156
1157 if (bDelim)
1158 GetOptPunct(pibIn, szLex, ',');
1159
1160 bErr = ENextLex(pibIn, szLex, LX_STRING);
1161
1162 if (!bErr) {
1163 if (szLex[0]) {
1164 /* Allocate and copy the string */
1165 if ( !(*pszArg = (PSTR) malloc(MyStrlen(szLex) + 1)))
1166 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetStringArg", NULL);
1167
1168 MyStrcpy(*pszArg, szLex);
1169 }
1170 else
1171 *pszArg = NULL; /* No string given */
1172 } /* if */
1173
1174 return(bErr);
1175
1176 } /* GetStringArg */
1177
1178
1179 /* ----------------------------------------------------------------------------
1180 GetOutputFile
1181
1182 Use a name different from the default for the regular output
1183 */
GetOutputFile(PINPUTBUF pibIn,PSTR szLex,PANALYSIS panal)1184 BOOL GetOutputFile (PINPUTBUF pibIn, PSTR szLex, PANALYSIS panal)
1185 {
1186 BOOL bErr = FALSE;
1187
1188 bErr = EGetPunct(pibIn, szLex, CH_LPAREN)
1189 || GetStringArg(pibIn, &panal->szOutfilename, szLex, FALSE);
1190
1191 if (!bErr) {
1192 panal->bAllocatedFileName = TRUE;
1193 bErr = EGetPunct(pibIn, szLex, CH_RPAREN);
1194 }
1195
1196 if (!bErr)
1197 bErr = GetTerminator(pibIn, szLex);
1198 else
1199 printf("Syntax: %s (szOutputFilename)\n\n", GetKeyword(KM_OUTPUTFILE));
1200
1201 return(bErr);
1202
1203 } /* GetOutputFile */
1204
1205
1206 /* ----------------------------------------------------------------------------
1207 GetSimulate: obsolete, give an ignore warning
1208 */
1209 BOOL bGaveSimulateUsage = FALSE; /* prevent multiple diagnostics */
1210
GetSimulate()1211 BOOL GetSimulate ()
1212 {
1213
1214 if (!bGaveSimulateUsage) {
1215 printf ("Warning: %s statements are obsolete and ignored.\n\n",
1216 GetKeyword(KM_SIMULATE));
1217 bGaveSimulateUsage = TRUE;
1218 }
1219
1220 return(1);
1221
1222 } /* GetSimulate */
1223
1224
1225 /* ----------------------------------------------------------------------------
1226 GetStartTime
1227 Reads the starting time for the integration (one argument).
1228
1229 * If the argument is a number, the start time takes that value.
1230
1231 * If the argument is an identifier which is a valid model
1232 parameter, the start time is defined as being dependent on that parameter.
1233 The symbolic parameter will be replaced by its value after model
1234 initialization in DoOneExperiment.
1235 */
1236 BOOL bGaveSrtTUsage = FALSE; /* prevent multiple diagnostics */
1237
GetStartTime(PINPUTBUF pibIn,PEXPERIMENT pexp)1238 BOOL GetStartTime (PINPUTBUF pibIn, PEXPERIMENT pexp)
1239 {
1240 static int vrgiSimArgTypes[1] = {LX_NUMBER | LX_IDENTIFIER};
1241
1242 BOOL bErr=!GetFuncArgs(pibIn, 1, vrgiSimArgTypes, vrgszlexArgs[0]);
1243
1244 if (!bErr) {
1245 if (!DefDepParm(vrgszlexArgs[0], &pexp->dT0, &pexp->hT0))
1246 ReportError(pibIn, RE_EXPECTED, "StartTime spec", NULL);
1247 }
1248 else {
1249 if (!bGaveSrtTUsage) {
1250 printf("Syntax: %s (InitialTime)\n\n", GetKeyword(KM_STARTTIME));
1251 bGaveSrtTUsage = TRUE;
1252 }
1253 }
1254
1255 return(bErr);
1256
1257 } /* GetStartTime */
1258
1259
1260 /* ----------------------------------------------------------------------------
1261 GetMCMCSpec
1262
1263 get the MCMC specification.
1264 */
GetMCMCSpec(PINPUTBUF pibIn,PEXPERIMENT pexp)1265 BOOL GetMCMCSpec (PINPUTBUF pibIn, PEXPERIMENT pexp)
1266 {
1267 #define NMCMC_ARGS 8 /* # Function arguments for MCMC spec */
1268
1269 static int vrgiGibbsArgTypes[NMCMC_ARGS] = {LX_STRING, LX_STRING, LX_STRING,
1270 LX_INTEGER, LX_INTEGER, LX_INTEGER,
1271 LX_INTEGER, LX_NUMBER};
1272
1273 PANALYSIS panal = (PANALYSIS) pibIn->pInfo;
1274
1275 BOOL bErr= !GetFuncArgs(pibIn, NMCMC_ARGS,
1276 vrgiGibbsArgTypes, vrgszlexArgs[0]);
1277
1278 static char vszGibbsOutDefault[] = "MCMC.default.out";
1279
1280 if (!bErr) {
1281 if (*vrgszlexArgs[0]) { /* Get output Filename */
1282 if ( !(panal->gd.szGout = (PSTR) malloc(MyStrlen(vrgszlexArgs[0]) + 1)))
1283 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetMCMCSpec", NULL);
1284
1285 MyStrcpy(panal->gd.szGout, vrgszlexArgs[0]);
1286 panal->bAllocatedFileName = TRUE;
1287 }
1288 else panal->gd.szGout = vszGibbsOutDefault;
1289
1290 if (*vrgszlexArgs[1]) { /* Get restart file */
1291 if ( !(panal->gd.szGrestart =
1292 (PSTR) malloc(MyStrlen(vrgszlexArgs[1]) + 1)))
1293 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetMCMCSpec", NULL);
1294
1295 MyStrcpy(panal->gd.szGrestart, vrgszlexArgs[1]);
1296 }
1297
1298 if (panal->gd.szGrestart != NULL &&
1299 !strcmp(panal->gd.szGout, panal->gd.szGrestart))
1300 ReportError(pibIn, RE_OUTISRESTART | RE_FATAL, "GetMCMCSpec", NULL);
1301
1302 if (*vrgszlexArgs[2]) { /* get external data file name */
1303 if ( !(panal->gd.szGdata = (PSTR) malloc(MyStrlen(vrgszlexArgs[2]) + 1)))
1304 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetMCMCSpec", NULL);
1305
1306 MyStrcpy(panal->gd.szGdata, vrgszlexArgs[2]);
1307 }
1308
1309 panal->gd.nMaxIter = atol(vrgszlexArgs[3]);
1310 panal->gd.nSimTypeFlag = atol(vrgszlexArgs[4]); /* should be 0 to 5 */
1311 panal->gd.nPrintFreq = atol(vrgszlexArgs[5]);
1312 panal->gd.nPrintIter = atol(vrgszlexArgs[6]);
1313
1314 panal->dSeed = atof(vrgszlexArgs[7]);
1315
1316 if (((panal->gd.nSimTypeFlag==1) && (panal->gd.szGrestart == NULL)) ||
1317 ((panal->gd.nSimTypeFlag==2) && (panal->gd.szGrestart == NULL))) {
1318 printf ("Error: if simTypeFlag is one or two a restart file must be "
1319 "given - Exiting\n\n");
1320 exit(0);
1321 }
1322 } /* if */
1323 else {
1324 printf ("Syntax:\n%s (szOut, szRestart, szData, "
1325 "nMaxIters, simTypeFlag, nPrintFreq,\n"
1326 " nIterToPrint, dSeed)\nExiting.\n\n",
1327 GetKeyword(KM_MCMC));
1328 exit(0);
1329 }
1330
1331 if (!bErr)
1332 panal->iType = AT_MCMC;
1333
1334 return(!bErr);
1335
1336 } /* GetMCMCSpec */
1337
1338
1339 /* ----------------------------------------------------------------------------
1340 GetOptDSpec
1341
1342 get the optimal design specification. It is based on the GetSetPointsSpec
1343 routine.
1344 The modification list is kept in MCVAR variation records, although this is
1345 not really a Monte Carlo analysis. This structure should eventually be
1346 changed to reflect a more general variation specification.
1347 */
GetOptDSpec(PINPUTBUF pibIn,PANALYSIS panal,PSTR szLex)1348 BOOL GetOptDSpec (PINPUTBUF pibIn, PANALYSIS panal, PSTR szLex)
1349 {
1350 PMCVAR pMCVar;
1351 HVAR hvar;
1352 int iErr = 0;
1353 int iNLI;
1354 int ikwcode;
1355
1356 /* Try to get open paren and filenames */
1357 if ((iErr = EGetPunct(pibIn, szLex, CH_LPAREN) ||
1358 GetStringArg(pibIn, &panal->gd.szGout, szLex, FALSE) ||
1359 GetStringArg(pibIn, &panal->gd.szGrestart, szLex, TRUE))) {
1360 goto Exit_GetOptDSpec;
1361 }
1362 else {
1363 panal->bAllocatedFileName = TRUE;
1364 }
1365
1366 /* There has to be a restart file */
1367 if (!panal->gd.szGrestart)
1368 ReportError(pibIn, RE_SPECERR | RE_FATAL, "Missing restart file", NULL);
1369
1370 /* Try to get number of parameter samples to read in */
1371 GetOptPunct(pibIn, szLex, ',');
1372 if ((iErr = ENextLex(pibIn, szLex, LX_INTEGER)))
1373 goto Exit_GetOptDSpec;
1374 panal->mc.nRuns = atol(szLex);
1375
1376 /* Try to get the random seed */
1377 GetOptPunct(pibIn, szLex, ',');
1378 if ((iErr = ENextLex(pibIn, szLex, LX_NUMBER)))
1379 goto Exit_GetOptDSpec;
1380 panal->dSeed = atof(szLex);
1381
1382 /* Try to get the style (Forward or Backward) */
1383 GetOptPunct(pibIn, szLex, ',');
1384 if ((iErr = ENextLex(pibIn, szLex, LX_IDENTIFIER)))
1385 goto Exit_GetOptDSpec;
1386 ikwcode = GetKeywordCode(szLex, NULL);
1387 if (ikwcode == KM_FORWARD)
1388 panal->mc.style = forward;
1389 else if (ikwcode == KM_BACKWARD)
1390 panal->mc.style = backward;
1391 else {
1392 iErr = TRUE;
1393 goto Exit_GetOptDSpec;
1394 }
1395
1396 /* Try to get identifier list */
1397 while ((iNLI=NextListItem(pibIn, szLex, LX_IDENTIFIER, 1, CH_RPAREN)) > 0) {
1398 hvar = GetVarHandle(szLex);
1399 if ((iErr = (!hvar || IsInput(hvar))))
1400 break; /* Is this reported ? */
1401
1402 if ( !(pMCVar = (PMCVAR) malloc(sizeof(MCVAR))))
1403 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetOptDSpec", NULL);
1404
1405 pMCVar->hvar = hvar;
1406 pMCVar->iType = MCV_SETPOINTS;
1407 pMCVar->dParm[2] = pMCVar->dParm[3] = 0.0;
1408
1409 QueueListItem(panal->mc.plistMCVars, pMCVar);
1410
1411 } /* while */
1412
1413 panal->mc.nSetParms = ListLength(panal->mc.plistMCVars);
1414
1415 if (panal->mc.nSetParms == 0) {
1416 iErr = TRUE;
1417 printf(
1418 "\nError: you must specify a list of parameters to read.\n\n");
1419 goto Exit_GetOptDSpec;
1420 }
1421
1422 if (!iNLI) /* List terminator */
1423 iErr = EGetPunct(pibIn, szLex, CH_RPAREN);
1424 else {
1425 iErr = TRUE;
1426 ReportError(pibIn, RE_LEXEXPECTED, "identifier", szLex);
1427 } /* else */
1428
1429 Exit_GetOptDSpec: ;
1430
1431 if (iErr) {
1432 printf ("Syntax:\n"
1433 "%s (\"Output_File\", \"Param_Sample_File\", nSamples, "
1434 "random_seed, <Forward or Backward>, "
1435 "<param-id-list...>)\n\n", GetKeyword(KM_OPTDESIGN));
1436 printf("Exiting...\n");
1437 exit(0);
1438 }
1439 else
1440 panal->iType = AT_OPTDESIGN; /* Flag SetPoints anal if not chosen */
1441
1442 return(iErr);
1443
1444 } /* GetOptDSpec */
1445
1446
1447 /* ----------------------------------------------------------------------------
1448 GetDistribSpec
1449
1450 reads in a Distrib statement (previously MCVary statement).
1451 */
1452 BOOL bGaveMCVaryUsage = FALSE; /* prevent multiple diagnostics */
1453
GetDistribSpec(PINPUTBUF pibIn,PSTR szLex,PANALYSIS panal)1454 int GetDistribSpec (PINPUTBUF pibIn, PSTR szLex, PANALYSIS panal)
1455 {
1456 PLIST plist;
1457 PMCVAR pMCVar = NULL;
1458 HVAR hvar;
1459 int n, iErr = 0;
1460 PSTRLEX szDummy;
1461
1462 /* if we are not doing Monte Carlo, SetPoints, Optimal Design or
1463 MCMC sampling, then ignore this specification */
1464 if (panal->iType && !((panal->iType == AT_MONTECARLO) ||
1465 (panal->iType == AT_SETPOINTS) ||
1466 (panal->iType == AT_OPTDESIGN) ||
1467 (panal->iType == AT_MCMC))) {
1468 EatStatement(pibIn); /* Ignore this Distrib() */
1469 goto Exit_MCVarySpec;
1470 }
1471
1472 /* Get the Distrib() spec. Check syntax at each element. */
1473
1474 /* Get the parameter to be varied */
1475 if ((iErr = (EGetPunct(pibIn, szLex, CH_LPAREN) ||
1476 ENextLex(pibIn, szLex, LX_IDENTIFIER))))
1477 goto Done_GetMCVary;
1478
1479 if (GetKeywordCode(szLex, NULL) == KM_DATA) { /* Data keyword used */
1480 /* try to get opening parenthesis */
1481 if (EGetPunct(pibIn, szLex, CH_LPAREN))
1482 exit(0); /* error */ /* this should be more graceful */
1483
1484 /* get the name of variable predicted */
1485 ENextLex(pibIn, szLex, LX_IDENTIFIER);
1486
1487 /* Data variable must exist and must be an input, state or output */
1488 if (!(hvar = GetVarHandle(szLex)) || IsParm(hvar))
1489 ReportError (pibIn, RE_LEXEXPECTED | RE_FATAL,
1490 "input, output or state variable", szLex);
1491
1492 /* try to read off closing parenthesis */
1493 if (EGetPunct(pibIn, szDummy, CH_RPAREN))
1494 exit(0); /* error */ /* this should be more graceful */
1495 }
1496 else
1497 if ((iErr = (!(hvar = GetVarHandle(szLex)) || /* Invalid variable name? */
1498 IsInput(hvar))))
1499 {
1500 ReportError(pibIn, RE_LEXEXPECTED, "state, output or parameter", szLex);
1501 goto Done_GetMCVary;
1502 } /* if */
1503
1504 /* Find the list in which pMCVar will be queued */
1505 if (panal->iCurrentDepth == 0) /* not an MCMC simulation */
1506 plist = panal->mc.plistMCVars;
1507 else { /* an MCMC simulation */
1508 if (!IsParm(hvar)) /* Distrib is for a likelihood definition */
1509 plist = panal->pCurrentLevel[panal->iCurrentDepth-1]->plistLikes;
1510 else
1511 plist = panal->pCurrentLevel[panal->iCurrentDepth-1]->plistMCVars;
1512 }
1513
1514 if (!(pMCVar = (PMCVAR) malloc(sizeof(MCVAR))))
1515 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetDistribSpec", NULL);
1516
1517 if(!(pMCVar->pszName = (PSTR) malloc(strlen(szLex)+1)))
1518 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetDistribSpec", NULL);
1519
1520 /* Initialize pMCVar */
1521 strcpy(pMCVar->pszName, szLex);
1522 pMCVar->hvar = hvar;
1523 pMCVar->pdVal = &(pMCVar->dVal);
1524 pMCVar->iDepth = panal->iCurrentDepth - 1;
1525 pMCVar->plistDependents = InitList();
1526 pMCVar->bExptIsDep = pMCVar->bIsFixed = FALSE;
1527 pMCVar->lJumps = pMCVar->lCount = 0;
1528 pMCVar->dKernelSD = INIT_KERNELSD;
1529 pMCVar->bGibbs = FALSE;
1530 for (n = 0; n < 4; n++) {
1531 pMCVar->hParm[n] = 0;
1532 pMCVar->pMCVParent[n] = NULL;
1533 pMCVar->pdParm[n] = &(pMCVar->dParm[n]);
1534 pMCVar->iParmType[n] = 0;
1535 }
1536
1537 /* Get the distribution type */
1538 GetOptPunct(pibIn, szLex, ',');
1539 iErr |= ENextLex(pibIn, szLex, LX_IDENTIFIER);
1540 pMCVar->iType = McvFromLex(szLex);
1541 if (iErr |= pMCVar->iType < 0) {
1542 ReportError(pibIn, RE_LEXEXPECTED, "distribution-type", szLex);
1543 goto Done_GetMCVary;
1544 }
1545
1546 /* Get parameters of the distribution. These vary by distribution type.
1547 No value checking is made because assignements can be symbolic
1548 */
1549 switch (pMCVar->iType) {
1550 /* ----------------------------------------------------------------------*/
1551 case MCV_UNIFORM:
1552 case MCV_LOGUNIFORM: /* 2 parameters: min and max */
1553
1554 if ((iErr = GetDistribParam(pibIn, szLex, plist, 0, pMCVar)))
1555 goto Done_GetMCVary;
1556
1557 if ((iErr = GetDistribParam(pibIn, szLex, plist, 1, pMCVar)))
1558 goto Done_GetMCVary;
1559
1560 pMCVar->dParm[2] = -DBL_MAX * 0.5;
1561 pMCVar->dParm[3] = DBL_MAX * 0.5;
1562
1563 break;
1564
1565 /* ----------------------------------------------------------------------*/
1566 case MCV_NORMAL: /* 2 parameters, mean and SD */
1567 case MCV_LOGNORMAL: /* ~ */
1568 case MCV_NORMALCV: /* 2 parameters, mean and coefficient of variation */
1569 case MCV_NORMALV: /* 2 parameters, mean and VARIANCE */
1570 case MCV_LOGNORMALV: /* ~ */
1571
1572 if ((iErr = GetDistribParam(pibIn, szLex, plist, 0, pMCVar)))
1573 goto Done_GetMCVary;
1574
1575 if ((iErr = GetDistribParam(pibIn, szLex, plist, 1, pMCVar)))
1576 goto Done_GetMCVary;
1577
1578 /* set the range */
1579 if ((pMCVar->iType == MCV_NORMAL) || (pMCVar->iType == MCV_NORMALCV) ||
1580 (pMCVar->iType == MCV_NORMALV)) {
1581 pMCVar->dParm[2] = -DBL_MAX * 0.5;
1582 pMCVar->dParm[3] = DBL_MAX * 0.5;
1583 }
1584 else {
1585 pMCVar->dParm[2] = 0.0;
1586 pMCVar->dParm[3] = DBL_MAX;
1587 }
1588
1589 break;
1590
1591 /* ----------------------------------------------------------------------*/
1592 case MCV_HALFNORMAL: /* 1 parameter: SD; mean is set to zero */
1593
1594 if ((iErr = GetDistribParam(pibIn, szLex, plist, 1, pMCVar)))
1595 goto Done_GetMCVary;
1596
1597 /* Set the mean */
1598 pMCVar->dParm[0] = 0.0;
1599
1600 /* Set the range */
1601 pMCVar->dParm[2] = 0.0;
1602 pMCVar->dParm[3] = DBL_MAX;
1603
1604 break;
1605
1606 /* ----------------------------------------------------------------------*/
1607 case MCV_BETA: /* 2 or 4 parameters */
1608 case MCV_TRUNCNORMAL: /* 4 parameters, the last 2 are min and max */
1609 case MCV_TRUNCLOGNORMAL: /* ~ */
1610 case MCV_TRUNCNORMALCV: /* Coefficient of variation instead of SD */
1611 case MCV_TRUNCNORMALV: /* VARIANCE instead of SD */
1612 case MCV_TRUNCLOGNORMALV: /* ~ */
1613
1614 if ((iErr = GetDistribParam(pibIn, szLex, plist, 0, pMCVar)))
1615 goto Done_GetMCVary;
1616
1617 if ((iErr = GetDistribParam(pibIn, szLex, plist, 1, pMCVar)))
1618 goto Done_GetMCVary;
1619
1620 /* Set min-max range defaults */
1621 pMCVar->dParm[2] = 0.0; /* Standard range for beta */
1622 pMCVar->dParm[3] = 1.0;
1623 if ((pMCVar->iType == MCV_TRUNCNORMAL) ||
1624 (pMCVar->iType == MCV_TRUNCNORMALCV) ||
1625 (pMCVar->iType == MCV_TRUNCNORMALV)) {
1626 pMCVar->dParm[2] = -DBL_MAX * 0.5;
1627 pMCVar->dParm[3] = DBL_MAX * 0.5;
1628 }
1629 else if ((pMCVar->iType == MCV_TRUNCLOGNORMAL) ||
1630 (pMCVar->iType == MCV_TRUNCLOGNORMALV))
1631 pMCVar->dParm[3] = DBL_MAX;
1632
1633 /* Look if a min-max range is included. For truncated types
1634 it is required. */
1635 SkipWhitespace(pibIn);
1636 if ((pMCVar->iType == MCV_BETA) && NextChar(pibIn) == CH_RPAREN)
1637 break; /* The spec is finished */
1638
1639 /* Get the min and max */
1640 if ((iErr = GetDistribParam(pibIn, szLex, plist, 2, pMCVar)))
1641 goto Done_GetMCVary;
1642
1643 if ((iErr = GetDistribParam(pibIn, szLex, plist, 3, pMCVar)))
1644 goto Done_GetMCVary;
1645
1646 break;
1647
1648 /* ----------------------------------------------------------------------*/
1649 case MCV_CHI2: /* only one parameter: degrees of freedom */
1650
1651 if ((iErr = GetDistribParam(pibIn, szLex, plist, 0, pMCVar)))
1652 goto Done_GetMCVary;
1653
1654 /* set the range */
1655 pMCVar->dParm[2] = 0.0;
1656 pMCVar->dParm[3] = DBL_MAX;
1657
1658 break;
1659
1660 /* ----------------------------------------------------------------------*/
1661 case MCV_BINOMIAL: /* 2 parameters, p and N */
1662
1663 if ((iErr = GetDistribParam(pibIn, szLex, plist, 0, pMCVar)))
1664 goto Done_GetMCVary;
1665
1666 if((iErr = GetDistribParam(pibIn, szLex, plist, 1, pMCVar)))
1667 goto Done_GetMCVary;
1668
1669 /* set the range */
1670 pMCVar->dParm[2] = 0.0; /* minimum */
1671 if (pMCVar->iParmType[1] != MCVP_FIXD)
1672 /* N symbolic: could be anything: assign DBL_MAX as maximum */
1673 pMCVar->dParm[3] = DBL_MAX; /* FB 18/07/97 */
1674 else
1675 /* N numeric: assign it as maximum */
1676 pMCVar->dParm[3] = pMCVar->dParm[1];
1677
1678 break;
1679
1680 /* ----------------------------------------------------------------------*/
1681 case MCV_NEGATIVEBINOM: /* 2 parameters, r and p */
1682
1683 if ((iErr = GetDistribParam(pibIn, szLex, plist, 0, pMCVar)))
1684 goto Done_GetMCVary;
1685
1686 if((iErr = GetDistribParam(pibIn, szLex, plist, 1, pMCVar)))
1687 goto Done_GetMCVary;
1688
1689 /* set the range */
1690 pMCVar->dParm[2] = 0.0; /* minimum */
1691 pMCVar->dParm[3] = DBL_MAX; /* maximum */
1692
1693 break;
1694
1695 /* ----------------------------------------------------------------------*/
1696 case MCV_PIECEWISE: /* 4 parameters, note the particular order */
1697
1698 if ((iErr = GetDistribParam(pibIn, szLex, plist, 2, pMCVar)))
1699 goto Done_GetMCVary;
1700
1701 if ((iErr = GetDistribParam(pibIn, szLex, plist, 0, pMCVar)))
1702 goto Done_GetMCVary;
1703
1704 if ((iErr = GetDistribParam(pibIn, szLex, plist, 1, pMCVar)))
1705 goto Done_GetMCVary;
1706
1707 if ((iErr = GetDistribParam(pibIn, szLex, plist, 3, pMCVar)))
1708 goto Done_GetMCVary;
1709
1710 break;
1711
1712 /* ----------------------------------------------------------------------*/
1713 case MCV_EXPONENTIAL: /* only one parameter: inverse scale */
1714
1715 if ((iErr = GetDistribParam(pibIn, szLex, plist, 0, pMCVar)))
1716 goto Done_GetMCVary;
1717
1718 /* set the range */
1719 pMCVar->dParm[2] = 0.0;
1720 pMCVar->dParm[3] = DBL_MAX;
1721
1722 break;
1723
1724 /* ----------------------------------------------------------------------*/
1725 case MCV_GGAMMA:
1726 case MCV_INVGGAMMA: /* 2 parameter: shape and inverse scale */
1727
1728 if ((iErr = GetDistribParam(pibIn, szLex, plist, 0, pMCVar)))
1729 goto Done_GetMCVary;
1730
1731 if ((iErr = GetDistribParam(pibIn, szLex, plist, 1, pMCVar)))
1732 goto Done_GetMCVary;
1733
1734 /* set the range */
1735 pMCVar->dParm[2] = 0.0;
1736 pMCVar->dParm[3] = DBL_MAX;
1737
1738 break;
1739
1740 /* ----------------------------------------------------------------------*/
1741 case MCV_TRUNCINVGGAMMA: /* 4 parameter: shape, inverse scale and bounds */
1742
1743 #ifndef HAVE_LIBGSL
1744 printf("Warning: The truncated inverse gamma density cannot be\n");
1745 printf(" used in MCMC simulations if the GNU Scientific\n");
1746 printf(" Library is not installed.\n");
1747 #endif
1748
1749 if ((iErr = GetDistribParam(pibIn, szLex, plist, 0, pMCVar)))
1750 goto Done_GetMCVary;
1751
1752 if ((iErr = GetDistribParam(pibIn, szLex, plist, 1, pMCVar)))
1753 goto Done_GetMCVary;
1754
1755 /* Get the min and max */
1756 if ((iErr = GetDistribParam(pibIn, szLex, plist, 2, pMCVar)))
1757 goto Done_GetMCVary;
1758
1759 if ((iErr = GetDistribParam(pibIn, szLex, plist, 3, pMCVar)))
1760 goto Done_GetMCVary;
1761
1762 break;
1763
1764 /* ----------------------------------------------------------------------*/
1765 case MCV_POISSON: /* 1 parameter: rate */
1766
1767 if ((iErr = GetDistribParam(pibIn, szLex, plist, 0, pMCVar)))
1768 goto Done_GetMCVary;
1769
1770 /* set the range */
1771 pMCVar->dParm[2] = 0.0;
1772 pMCVar->dParm[3] = DBL_MAX;
1773
1774 break;
1775
1776 /* ----------------------------------------------------------------------*/
1777 case MCV_BINOMIALBETA: /* 3 parameter: mean, alpha, beta */
1778 case MCV_GENLOGNORMAL: /* 3 parameter: mean, sdnorm, sdlognorm */
1779 case MCV_STUDENTT: /* 3 parameter: dof, location, scale */
1780
1781 if ((iErr = GetDistribParam(pibIn, szLex, plist, 0, pMCVar)))
1782 goto Done_GetMCVary;
1783
1784 if ((iErr = GetDistribParam(pibIn, szLex, plist, 1, pMCVar)))
1785 goto Done_GetMCVary;
1786
1787 if ((iErr = GetDistribParam(pibIn, szLex, plist, 2, pMCVar)))
1788 goto Done_GetMCVary;
1789
1790 /* set the last parameter, unused */
1791 pMCVar->dParm[3] = DBL_MAX;
1792
1793 break;
1794
1795 /* ----------------------------------------------------------------------*/
1796 case MCV_CAUCHY: /* one parameter: scale */
1797
1798 if ((iErr = GetDistribParam(pibIn, szLex, plist, 0, pMCVar)))
1799 goto Done_GetMCVary;
1800
1801 pMCVar->dParm[1] = DBL_MAX;
1802 pMCVar->dParm[2] = -DBL_MAX * 0.5;
1803 pMCVar->dParm[3] = DBL_MAX * 0.5;
1804
1805 break;
1806
1807 /* ----------------------------------------------------------------------*/
1808 case MCV_HALFCAUCHY: /* one parameter: scale */
1809
1810 if ((iErr = GetDistribParam(pibIn, szLex, plist, 0, pMCVar)))
1811 goto Done_GetMCVary;
1812
1813 pMCVar->dParm[1] = DBL_MAX * 0.5;
1814 pMCVar->dParm[2] = 0;
1815 pMCVar->dParm[3] = DBL_MAX;
1816
1817 break;
1818
1819 /* ----------------------------------------------------------------------*/
1820 case MCV_USERLL: /* one parameter: the model computed likelihood */
1821
1822 if (!panal->pCurrentLevel[panal->iCurrentDepth-1] ||
1823 (plist != panal->pCurrentLevel[panal->iCurrentDepth-1]->plistLikes)) {
1824 printf("UserSpecifiefLL can only be used in Likelihood().");
1825 iErr = 1;
1826 goto Done_GetMCVary;
1827 }
1828
1829 if ((iErr = GetDistribParam(pibIn, szLex, plist, 0, pMCVar)))
1830 goto Done_GetMCVary;
1831
1832 pMCVar->dParm[1] = DBL_MAX * 0.5;
1833 pMCVar->dParm[2] = 0;
1834 pMCVar->dParm[3] = DBL_MAX;
1835
1836 break;
1837
1838 /* ----------------------------------------------------------------------*/
1839 default:
1840 ReportRunTimeError(panal, RE_UNKNOWNDIST | RE_FATAL, "GetDistribSpec");
1841 break;
1842
1843 } /* switch */
1844
1845 EGetPunct(pibIn, szLex, CH_RPAREN);
1846
1847 /* Check for a range error if the bounds are numeric. If there is a problem,
1848 correct it, but issue a warning in case this is wrong. */
1849 if ((pMCVar->iParmType[2] == MCVP_FIXD) &&
1850 (pMCVar->iParmType[3] == MCVP_FIXD) &&
1851 (pMCVar->dParm[3] < pMCVar->dParm[2])) {
1852 double dTmp = pMCVar->dParm[3]; /* Swap ranges */
1853 pMCVar->dParm[3] = pMCVar->dParm[2];
1854 pMCVar->dParm[2] = dTmp;
1855 ReportError(pibIn, RE_MAXMIN_RANGE | RE_WARNING, NULL, NULL);
1856 }
1857
1858 /* If there's no error at this point, queue the variation in
1859 the Monte Carlo record(s). */
1860 if (!iErr) {
1861 QueueListItem(plist, pMCVar);
1862 } /* if */
1863
1864 Done_GetMCVary: ;
1865
1866 if (iErr) {
1867 if (pMCVar) free(pMCVar);
1868
1869 if (!bGaveMCVaryUsage) {
1870 printf("\nSyntax: Check the syntax of %s.\n", GetKeyword (KM_MCVARY));
1871 bGaveMCVaryUsage = TRUE;
1872 }
1873
1874 ReportError(pibIn, RE_SYNTAXERR | RE_FATAL, NULL, NULL);
1875 }
1876
1877 Exit_MCVarySpec: ;
1878
1879 return(iErr);
1880
1881 } /* GetDistribSpec */
1882
1883
1884 /* ----------------------------------------------------------------------------
1885 CheckDistribParam
1886
1887 If the nth distribution parameter for hvar2 is a parameter and the same as
1888 the variable hvar1 for which the Distrib statement is specified,
1889 check false. Used to check that you don't have parameter self-dependency
1890 at level 0.
1891 */
CheckDistribParam(PLIST plist,HVAR hvar1,HVAR hvar2)1892 BOOL CheckDistribParam (PLIST plist, HVAR hvar1, HVAR hvar2) {
1893 int n;
1894 PLISTELEM p = plist->pleHead;
1895 PMCVAR pMCVar;
1896
1897 if (plist == NULL) return TRUE;
1898
1899 for (n = 0; n < plist->iSize; ++n) {
1900 pMCVar = (PMCVAR) p->pData;
1901 if (hvar2 == pMCVar->hvar) {
1902 if ((pMCVar->iParmType[0] == MCVP_PARM) && (hvar1 == pMCVar->hParm[0]))
1903 return FALSE;
1904 if ((pMCVar->iParmType[1] == MCVP_PARM) && (hvar1 == pMCVar->hParm[1]))
1905 return FALSE;
1906 if ((pMCVar->iParmType[2] == MCVP_PARM) && (hvar1 == pMCVar->hParm[2]))
1907 return FALSE;
1908 if ((pMCVar->iParmType[3] == MCVP_PARM) && (hvar1 == pMCVar->hParm[3]))
1909 return FALSE;
1910 }
1911 p = p->pleNext;
1912 }
1913
1914 return TRUE;
1915
1916 } /* CheckDistribParam */
1917
1918
1919 /* ----------------------------------------------------------------------------
1920 GetDistribParam
1921
1922 Determine if argument `n' of the Distrib statement is a variable name or
1923 a number; set the parameter accordingly.
1924
1925 If the argument is a variable, set a pointer to its MC structure
1926 */
GetDistribParam(PINPUTBUF pibIn,PSTR szLex,PLIST plist,int n,PMCVAR pMCVar)1927 int GetDistribParam (PINPUTBUF pibIn, PSTR szLex, PLIST plist, int n,
1928 PMCVAR pMCVar) {
1929 PANALYSIS panal = (PANALYSIS) pibIn->pInfo;
1930 int iLex, iCode;
1931 HVAR hvar;
1932
1933 GetOptPunct(pibIn, szLex, ',');
1934 if (n != 3)
1935 NextLex(pibIn, szLex, &iLex);
1936 else {
1937 SkipWhitespace(pibIn);
1938 iLex = LX_NULL;
1939 if (NextChar(pibIn) != CH_RPAREN)
1940 NextLex(pibIn, szLex, &iLex);
1941 }
1942
1943 if (iLex == LX_IDENTIFIER) { /* symbol used for parameter */
1944
1945 iCode = GetKeywordCode_in_context(szLex, CN_FUNCARG);
1946
1947 if ((iCode == KM_PREDICTION) || (iCode == KM_DATA)) {
1948 /* Prediction or Data keywords used */
1949
1950 /* Only inputs, states and outputs can have predicted or data
1951 parameters */
1952 if (IsParm(pMCVar->hvar))
1953 ReportError(pibIn, RE_BADCONTEXT | RE_FATAL, szLex, NULL);
1954
1955 /* Try to get opening parenthesis */
1956 if (EGetPunct(pibIn, szLex, CH_LPAREN))
1957 return 1; /* error */
1958
1959 /* Get the name of variable specified */
1960 NextLex(pibIn, szLex, &iLex);
1961
1962 /* Specified variable must exist and be input, state or output */
1963 if (!(hvar = GetVarHandle(szLex)) || IsParm(hvar))
1964 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL,
1965 "input, output or state variable", szLex);
1966
1967 /* Try to get closing parenthesis */
1968 if (EGetPunct(pibIn, szLex, CH_RPAREN))
1969 return 1; /* error */
1970
1971 } /* end of Prediction/Data case */
1972 else {
1973 /* No keyword used, a regular symbol: that symbol should
1974 be declared and should not be an input or an output or a state
1975 (because those should use keywords); so it should be
1976 a parameter - FB 18/10/98 */
1977 if (!(hvar = GetVarHandle(szLex)) || !IsParm(hvar))
1978 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, "parameter", szLex);
1979 }
1980
1981 /* No self-dependency at level 0 allowed, except in Optimal design */
1982 if (!(panal->iType == AT_OPTDESIGN) &&
1983 ((panal->iCurrentDepth == 0 && hvar == pMCVar->hvar) ||
1984 !CheckDistribParam(plist, pMCVar->hvar, hvar)))
1985 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, "valid parameter", szLex);
1986
1987 /* Declare it a symbolic parameter */
1988 if (iCode == KM_PREDICTION)
1989 pMCVar->iParmType[n] = MCVP_PRED;
1990 else if (iCode == KM_DATA)
1991 pMCVar->iParmType[n] = MCVP_DATA;
1992 else
1993 pMCVar->iParmType[n] = MCVP_PARM;
1994
1995 /* Attach handle */
1996 pMCVar->hParm[n] = hvar;
1997
1998 }
1999 else /* i.e. not an identifier */
2000 if (iLex == LX_FLOAT || iLex == LX_INTEGER) {
2001 pMCVar->iParmType[n] = MCVP_FIXD;
2002 pMCVar->dParm[n] = atof(szLex);
2003 }
2004 else {
2005 /* allow max to be absent - set to default */
2006 if (n == 3) {
2007 pMCVar->iParmType[n] = MCVP_FIXD;
2008 pMCVar->dParm[n] = DBL_MAX;
2009 }
2010 else
2011 return 1; /* error */
2012 }
2013
2014 return 0; /* OK */
2015
2016 } /* GetDistribParam */
2017
2018
2019 /* ----------------------------------------------------------------------------
2020 GetSetPointsSpec
2021
2022 Reads the SetPoints() arguments. The modification list is kept
2023 in MCVAR variation records, although this is not really a Monte
2024 Carlo analysis. This structure should eventually be changed to
2025 reflect a more general variation specification.
2026 */
GetSetPointsSpec(PINPUTBUF pibIn,PANALYSIS panal,PSTR szLex)2027 int GetSetPointsSpec (PINPUTBUF pibIn, PANALYSIS panal, PSTR szLex)
2028 {
2029 PMCVAR pMCVar;
2030 PSTRLEX szTmp;
2031 HVAR hvar;
2032 int iErr = 0;
2033 int iNLI;
2034 long j, iLB, iUB;
2035
2036 /* MonteCarlo sampling can be mixed with SetPoints sampling if Distrib
2037 specs appear after the SetPoint spec */
2038 if (ListLength(panal->mc.plistMCVars) > 0) {
2039 printf ("Error: Distrib() statements can only appear after the SetPoints()"
2040 "specification, not before - Exiting\n\n");
2041 exit(0);
2042 }
2043
2044 /* Try to get open paren and filenames */
2045 if ((iErr = EGetPunct(pibIn, szLex, CH_LPAREN) ||
2046 GetStringArg(pibIn, &panal->mc.szMCOutfilename, szLex, FALSE) ||
2047 GetStringArg(pibIn, &panal->mc.szSetPointsFilename, szLex,
2048 TRUE))) {
2049 goto Exit_GetSetPointsSpec;
2050 }
2051 else {
2052 panal->bAllocatedFileName = TRUE;
2053 }
2054
2055 /* There has to be a restart file */
2056 if (!panal->mc.szSetPointsFilename)
2057 ReportError(pibIn, RE_SPECERR | RE_FATAL, "Missing setpoints file", NULL);
2058
2059 /* Output file and setpoint file have to be different */
2060 if (!MyStrcmp(panal->mc.szMCOutfilename, panal->mc.szSetPointsFilename))
2061 ReportError(pibIn, RE_SPECERR | RE_FATAL, "Same name for 2 files", NULL);
2062
2063 /* Try to get number of runs */
2064 GetOptPunct(pibIn, szLex, ',');
2065 if ((iErr = ENextLex(pibIn, szLex, LX_INTEGER)))
2066 goto Exit_GetSetPointsSpec;
2067
2068 panal->mc.nRuns = atol(szLex);
2069
2070 /* Try to get identifier list */
2071 while ((iNLI=NextListItem(pibIn, szLex, LX_IDENTIFIER, 1, CH_RPAREN)) > 0) {
2072 /* check if szLex is a scalar or an array and act accordingly */
2073 iLB = iUB = -1;
2074 if (GetPunct(pibIn, szTmp, '[')) /* array found, read bounds */
2075 GetArrayBounds(pibIn, &iLB, &iUB);
2076
2077 if (iUB == -1) { /* scalar, store and continue */
2078 hvar = GetVarHandle(szLex);
2079 if ((iErr = (!hvar || IsInput(hvar))))
2080 break; /* will generate an error message */
2081
2082 if ( !(pMCVar = (PMCVAR) malloc(sizeof(MCVAR))))
2083 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetSetPointsSpec", NULL);
2084
2085 pMCVar->hvar = hvar;
2086 pMCVar->iType = MCV_SETPOINTS;
2087 pMCVar->dParm[2] = pMCVar->dParm[3] = 0.0;
2088
2089 QueueListItem(panal->mc.plistMCVars, pMCVar);
2090 }
2091 else { /* array */
2092
2093 for (j = iLB; j < iUB; j++) {
2094 sprintf(szTmp, "%s_%ld", szLex, j); /* create names */
2095
2096 hvar = GetVarHandle(szTmp);
2097 if ((iErr = (!hvar || IsInput(hvar))))
2098 break; /* will generate an error message */
2099
2100 if ( !(pMCVar = (PMCVAR) malloc(sizeof(MCVAR))))
2101 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetSetPointsSpec", NULL);
2102
2103 pMCVar->hvar = hvar;
2104 pMCVar->iType = MCV_SETPOINTS;
2105 pMCVar->dParm[2] = pMCVar->dParm[3] = 0.0;
2106
2107 QueueListItem(panal->mc.plistMCVars, pMCVar);
2108 }
2109 } /* end else */
2110
2111 } /* while */
2112
2113 panal->mc.nSetParms = ListLength(panal->mc.plistMCVars);
2114
2115 /* FB 19 nov 96 */
2116 if (panal->mc.nSetParms == 0) {
2117 iErr = TRUE;
2118 printf(
2119 "\nError: you must specify a list of parameters to read.\n\n");
2120 goto Exit_GetSetPointsSpec;
2121 }
2122
2123 if (!iNLI) /* List terminator */
2124 iErr = ((szTmp[0] != CH_RPAREN) && (EGetPunct(pibIn, szLex, CH_RPAREN))) ||
2125 InitSetPoints(&panal->mc);
2126 else {
2127 iErr = TRUE;
2128 ReportError(pibIn, RE_LEXEXPECTED, "identifier", szLex);
2129 } /* else */
2130
2131 Exit_GetSetPointsSpec: ;
2132
2133 if (iErr) {
2134 printf ("Syntax:\n"
2135 "%s (\"OutputFile\", \"SetPtsFile\", nRuns, "
2136 "<param-id-list...>)\n\n", GetKeyword(KM_SETPOINTS));
2137 printf("Exiting...\n");
2138 exit(0);
2139 }
2140 else
2141 panal->iType = AT_SETPOINTS; /* Flag SetPoints anal */
2142
2143 return(iErr);
2144
2145 } /* GetSetPointsSpec */
2146
2147
2148 /* ----------------------------------------------------------------------------
2149 GetMonteCarloSpec
2150 */
GetMonteCarloSpec(PINPUTBUF pibIn,PANALYSIS panal,PSTR szLex)2151 int GetMonteCarloSpec (PINPUTBUF pibIn, PANALYSIS panal, PSTR szLex)
2152 {
2153 #define NMC_ARGS 3 /* 3 MonteCarlo Args */
2154
2155 static int vrgiMCArgTypes[NMC_ARGS] = {LX_STRING, LX_INTEGER, LX_NUMBER};
2156
2157 int iErr = 0;
2158
2159 iErr = !GetFuncArgs(pibIn, NMC_ARGS, vrgiMCArgTypes, vrgszlexArgs[0]);
2160
2161 if (!iErr) {
2162 if (*vrgszlexArgs[0]) {
2163 if ( !(panal->mc.szMCOutfilename =
2164 (PSTR) malloc(MyStrlen(vrgszlexArgs[0]) + 1)))
2165 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetMonteCarloSpec", NULL);
2166
2167 MyStrcpy(panal->mc.szMCOutfilename, vrgszlexArgs[0]);
2168 panal->bAllocatedFileName = TRUE;
2169 }
2170
2171 panal->mc.nRuns = atol(vrgszlexArgs[1]);
2172 panal->dSeed = atof(vrgszlexArgs[2]);
2173 } /* if */
2174 else
2175 printf ("Syntax: %s (szOutfilename, nRuns, dSeed)\n\n",
2176 GetKeyword(KM_MONTECARLO));
2177
2178 if (!iErr)
2179 panal->iType = AT_MONTECARLO; /* Flag as MC */
2180
2181 return(iErr);
2182
2183 } /* GetMonteCarloSpec */
2184
2185
2186 /* ----------------------------------------------------------------------------
2187 GetParmMod
2188
2189 Reads a variable assignment and store the assigned value (or link if
2190 symbolic assignment)
2191 */
GetParmMod(PINPUTBUF pibIn,PSTRLEX szLex)2192 BOOL GetParmMod (PINPUTBUF pibIn, PSTRLEX szLex)
2193 {
2194 HVAR hvar = GetVarHandle(szLex);
2195 PANALYSIS panal = (PANALYSIS) pibIn->pInfo;
2196 PEXPERIMENT pexp = panal->pexpCurrent;
2197
2198 PSTRLEX szPunct;
2199 int iErr;
2200 PVARMOD pvarmod; /* Pointer to the variable modification */
2201
2202 if ((iErr = !hvar))
2203 ReportError(pibIn, RE_LEXEXPECTED, "model-variable", szLex);
2204
2205 else {
2206 /* Allocate space and initialize modification */
2207
2208 if ( !(pvarmod = (PVARMOD) malloc(sizeof(VARMODIFICATION))))
2209 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetParmMod", NULL);
2210
2211 pvarmod->hvar = hvar; /* The variable handle */
2212
2213 if (!GetOptPunct(pibIn, szPunct, '=')) { /* Try to get '=' */
2214 iErr = szPunct[1] = '=';
2215 ReportError(pibIn, RE_EXPECTED, szPunct, NULL);
2216 }
2217
2218 else if (IsInput(hvar)) { /* Process INPUT */
2219 if ( !(pvarmod->uvar.pifn = (PIFN) malloc(sizeof(IFN))))
2220 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetParmMod", NULL);
2221
2222 iErr = !pvarmod->uvar.pifn
2223 || !GetInputFn(pibIn, NULL, pvarmod->uvar.pifn);
2224 if (iErr) {
2225 free(pvarmod->uvar.pifn); /* Cleanup if error */
2226 pvarmod->uvar.pifn = NULL;
2227 }
2228 } /* else if */
2229 else { /* PARM, STATE, etc */
2230 if (!(iErr = ENextLex(pibIn, szLex, LX_NUMBER)))
2231 pvarmod->uvar.dVal = atof(szLex);
2232 }
2233
2234 if (!iErr) { /* No errors, add mod to list */
2235 if(panal->iCurrentDepth == 0 || panal->wContext == CN_EXPERIMENT)
2236 QueueListItem(pexp->plistParmMods, pvarmod);
2237 else
2238 QueueListItem(panal->pCurrentLevel[panal->iCurrentDepth-1]->plistVars,
2239 pvarmod);
2240 iErr = GetTerminator(pibIn, szLex);
2241 } /* if */
2242 else /* Invalid mod, cleanup */
2243 free(pvarmod);
2244
2245 } /* else valid id */
2246
2247 return((BOOL) iErr);
2248
2249 } /* GetParmMod */
2250
2251
2252 /* ----------------------------------------------------------------------------
2253 GetParmMod2
2254
2255 Store the assigned value (or link if symbolic assignment) passed by szeqn.
2256 */
GetParmMod2(PINPUTBUF pibIn,PSTRLEX szLex,PSTREQN szEqn)2257 BOOL GetParmMod2 (PINPUTBUF pibIn, PSTRLEX szLex, PSTREQN szEqn)
2258 {
2259 int iErr;
2260 PVARMOD pvarmod; /* Pointer to the variable modification */
2261
2262 HVAR hvar = GetVarHandle(szLex);
2263 PANALYSIS panal = (PANALYSIS) pibIn->pInfo;
2264 PEXPERIMENT pexp = panal->pexpCurrent;
2265
2266 if ((iErr = !hvar))
2267 ReportError(pibIn, RE_LEXEXPECTED, "model-variable", szLex);
2268
2269 else { /* valid ID */
2270 /* Allocate space and initialize modification */
2271
2272 if ( !(pvarmod = (PVARMOD) malloc(sizeof(VARMODIFICATION))))
2273 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetParmMod", NULL);
2274
2275 pvarmod->hvar = hvar; /* The variable handle */
2276
2277 if (IsInput(hvar)) { /* Process INPUT */
2278 if ( !(pvarmod->uvar.pifn = (PIFN) malloc(sizeof(IFN))))
2279 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "GetParmMod", NULL);
2280
2281 iErr = !pvarmod->uvar.pifn
2282 || !GetInputFn(pibIn, NULL, pvarmod->uvar.pifn);
2283 if (iErr) {
2284 free(pvarmod->uvar.pifn); /* Cleanup if error */
2285 pvarmod->uvar.pifn = NULL;
2286 }
2287 }
2288 else { /* PARM, STATE, etc */
2289 pvarmod->uvar.dVal = atof(szEqn);
2290 }
2291
2292 if (!iErr) { /* No errors, add mod to list */
2293 if(panal->iCurrentDepth == 0 || panal->wContext == CN_EXPERIMENT)
2294 QueueListItem(pexp->plistParmMods, pvarmod);
2295 else
2296 QueueListItem(panal->pCurrentLevel[panal->iCurrentDepth-1]->plistVars,
2297 pvarmod);
2298 }
2299 else /* Invalid mod, cleanup */
2300 free(pvarmod);
2301
2302 } /* else valid id */
2303
2304 return((BOOL) iErr);
2305
2306 } /* GetParmMod2 */
2307
2308
2309 /* ----------------------------------------------------------------------------
2310 NewExperiment
2311
2312 creates a new experiment in the analysis and copies global defaults.
2313 */
NewExperiment(PINPUTBUF pibIn)2314 void NewExperiment (PINPUTBUF pibIn)
2315 {
2316 PANALYSIS panal = (PANALYSIS)pibIn->pInfo;
2317 PLEVEL plevel;
2318 int n;
2319
2320 /* Allocate new experiment and assign list and current pointers */
2321
2322 if (panal->iCurrentDepth < 0) { /* something is real wrong - FB 03/08/97 */
2323 ReportError (pibIn, RE_LEXEXPECTED | RE_FATAL, "Level statement",
2324 "Simulation");
2325 }
2326
2327 if (panal->iCurrentDepth == 0) {
2328 panal->expGlobal.iExp++; /* Increment number of experiment */
2329 panal->pexpCurrent = panal->rgpExps[panal->expGlobal.iExp - 1] =
2330 (PEXPERIMENT) malloc(sizeof(EXPERIMENT));
2331 if (!panal->pexpCurrent)
2332 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "NewExperiment()", NULL);
2333
2334 if (panal->rank == 0)
2335 printf("Reading experiment %d.\n", panal->expGlobal.iExp);
2336 }
2337 else {
2338 plevel = panal->pLevels[panal->iInstances - 1];
2339 for (n = 0; n < panal->iCurrentDepth-1; ++n) {
2340 plevel = plevel->pLevels[plevel->iInstances - 1];
2341 }
2342 if (plevel->iInstances == MAX_INSTANCES - 1)
2343 ReportError(pibIn, RE_TOOMANYINST | RE_FATAL, "NewExperiment", NULL);
2344 n = panal->pCurrentLevel[panal->iCurrentDepth - 1]->iInstances++;
2345 if (!(plevel = plevel->pLevels[n] = (PLEVEL)malloc(sizeof(LEVEL))))
2346 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "NewExperiment", NULL);
2347 plevel->iInstances = 0;
2348 plevel->iSequence = n + 1;
2349 plevel->iDepth = panal->iCurrentDepth;
2350 panal->pCurrentLevel[panal->iCurrentDepth++] = plevel;
2351 if (panal->iDepth < panal->iCurrentDepth)
2352 panal->iDepth = panal->iCurrentDepth;
2353
2354 plevel->nMCVars = plevel->nFixedVars = plevel->nLikes = 0;
2355 plevel->plistVars = InitList();
2356 plevel->plistMCVars = InitList();
2357 plevel->plistLikes = InitList();
2358
2359 if (!(plevel->pexpt = (PEXPERIMENT) malloc(sizeof(EXPERIMENT))))
2360 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "NewExperiment", NULL);
2361
2362 panal->pexpCurrent = plevel->pexpt;
2363 panal->pexpCurrent->iExp = panal->expGlobal.iExp = ++panal->iExpts;
2364 panal->wContext = CN_EXPERIMENT;
2365
2366 if (panal->rank == 0)
2367 printf ("Simulation %d - depth %d, instance %d\n", panal->iExpts,
2368 panal->iCurrentDepth,
2369 panal->pCurrentLevel[panal->iCurrentDepth-2]->iInstances);
2370 }
2371
2372 memcpy(panal->pexpCurrent, &panal->expGlobal, sizeof(EXPERIMENT));
2373 panal->wContext = CN_EXPERIMENT;
2374 panal->pexpCurrent->plistParmMods = InitList(); /* Local mods */
2375 panal->pexpCurrent->os.plistPrintRecs = InitList();
2376 panal->pexpCurrent->os.plistDataRecs = InitList();
2377
2378 } /* NewExperiment */
2379
2380
2381 /* ----------------------------------------------------------------------------
2382 EndExperiment
2383
2384 cleans up at the end of defining a new experiment section.
2385 */
EndExperiment(PINPUTBUF pibIn,PANALYSIS panal)2386 BOOL EndExperiment (PINPUTBUF pibIn, PANALYSIS panal)
2387 {
2388 BOOL bReturn;
2389
2390 bReturn = !ErrorsReported(pibIn);
2391
2392 if (!bReturn) {
2393 /* Experiment had errors. Cleanup this space and continue */
2394 ReportError(pibIn, RE_ERRORSINEXP | RE_FATAL,
2395 (PSTR) &panal->pexpCurrent->iExp, NULL);
2396 ClearErrors(pibIn);
2397 panal->rgpExps[--panal->expGlobal.iExp] = NULL;
2398 free(panal->pexpCurrent);
2399 } /* if */
2400
2401 else {
2402 /* Create space for outputs and data */
2403 PrepareOutSpec(panal->pexpCurrent);
2404 }
2405
2406 /* Reset current exp to global context. */
2407
2408 panal->pexpCurrent = &panal->expGlobal;
2409 panal->wContext = CN_GLOBAL;
2410
2411 if (panal->iType == AT_MCMC && panal->iCurrentDepth-- == 0)
2412 return FALSE;
2413
2414 return(bReturn);
2415
2416 } /* EndExperiment */
2417
2418
2419 /* ----------------------------------------------------------------------------
2420 SetLevel
2421
2422 Encountered `Level' keyword, increments number of levels, allocates
2423 structure, initializes
2424 */
SetLevel(PINPUTBUF pibIn)2425 int SetLevel(PINPUTBUF pibIn)
2426 {
2427 PSTRLEX szPunct;
2428 PANALYSIS panal = (PANALYSIS)pibIn->pInfo;
2429 PLEVEL plevel;
2430 int n;
2431
2432 if (panal->iType != AT_MCMC)
2433 ReportError(pibIn, RE_TYPENOTMCMC | RE_FATAL, "SetLevel", NULL);
2434
2435 if (panal->iCurrentDepth == MAX_LEVELS)
2436 ReportError(pibIn, RE_TOOMANYLEVELS | RE_FATAL, "SetLevel", NULL);
2437
2438 if (panal->wContext == CN_EXPERIMENT)
2439 ReportError(pibIn, RE_LEVINEXPT | RE_FATAL, "SetLevel", NULL);
2440
2441 if (EGetPunct(pibIn, szPunct, CH_LBRACE))
2442 return 1;
2443
2444 if (panal->iCurrentDepth == 0) {
2445
2446 plevel = panal->pLevels[panal->iInstances++]
2447 = (PLEVEL) malloc(sizeof(LEVEL));
2448 if (plevel == NULL)
2449 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "SetLevel", NULL);
2450
2451 /* we do not want more than 1 top level, exit */
2452 if (panal->iInstances > 1) {
2453 printf("Error: only one top level is allowed - Exiting.\n");
2454 exit(0);
2455 }
2456
2457 plevel->iSequence = panal->iInstances;
2458 if (panal->rank == 0)
2459 printf("New level - depth 1, instance %d\n", panal->iInstances);
2460
2461 }
2462 else {
2463
2464 plevel = panal->pLevels[panal->iInstances-1];
2465
2466 for (n = 0; n < panal->iCurrentDepth-1; ++n)
2467 plevel = plevel->pLevels[plevel->iInstances - 1];
2468
2469 if (plevel->iInstances == MAX_INSTANCES - 1)
2470 ReportError(pibIn, RE_TOOMANYINST | RE_FATAL, "SetLevel", NULL);
2471
2472 n = panal->pCurrentLevel[panal->iCurrentDepth-1]->iInstances++;
2473
2474 plevel = plevel->pLevels[n]
2475 = (PLEVEL)malloc(sizeof(LEVEL));
2476
2477 if (plevel == NULL)
2478 ReportError(pibIn, RE_OUTOFMEM | RE_FATAL, "SetLevel", NULL);
2479
2480 plevel->iSequence = n + 1;
2481 if (panal->rank == 0)
2482 printf("New level - depth %d, instance %d\n", panal->iCurrentDepth + 1,
2483 panal->pCurrentLevel[panal->iCurrentDepth-1]->iInstances);
2484 }
2485
2486 plevel->iInstances = 0;
2487 plevel->iDepth = panal->iCurrentDepth;
2488 panal->pCurrentLevel[panal->iCurrentDepth++] = plevel;
2489
2490 if (panal->iDepth < panal->iCurrentDepth)
2491 panal->iDepth = panal->iCurrentDepth;
2492
2493 plevel->nMCVars = plevel->nFixedVars = plevel->nLikes = 0;
2494 plevel->plistVars = InitList();
2495 plevel->plistMCVars = InitList();
2496 plevel->plistLikes = InitList();
2497 plevel->pexpt = NULL;
2498
2499 return 0;
2500
2501 } /* SetLevel */
2502
2503
2504 /* ----------------------------------------------------------------------------
2505 EndLevel
2506 */
EndLevel(PANALYSIS panal)2507 BOOL EndLevel (PANALYSIS panal)
2508 {
2509 if(panal->iCurrentDepth-- == 0)
2510 return FALSE;
2511 return TRUE;
2512
2513 } /* EndLevel */
2514
2515
2516 /* ----------------------------------------------------------------------------
2517 FreeLevels
2518 */
FreeLevels(PANALYSIS panal)2519 void FreeLevels (PANALYSIS panal)
2520 {
2521 int n;
2522
2523 for (n = 0; n < panal->iInstances; n++)
2524 if (panal->pLevels[n] != NULL)
2525 FreeOneLevel(panal->pLevels[n]);
2526
2527 if (panal->bAllocatedFileName)
2528 free(panal->gd.szGout);
2529
2530 FreeList(&panal->mc.plistMCVars, NULL, TRUE);
2531 FreeList(&panal->expGlobal.plistParmMods, NULL, TRUE);
2532
2533 free(panal->expGlobal.is.iwork);
2534 free(panal->expGlobal.is.rwork);
2535
2536 free(panal->modelinfo.pStateHvar);
2537 free(panal);
2538
2539 } /* FreeLevels */
2540
2541
2542 /* ----------------------------------------------------------------------------
2543 FreeMCVar
2544 */
FreeMCVar(PVOID pData,PVOID pUserInfo)2545 int FreeMCVar (PVOID pData, PVOID pUserInfo)
2546 {
2547 PMCVAR pMCVar = (PMCVAR) pData;
2548 FreeList(&pMCVar->plistDependents, NULL, TRUE);
2549 free(pMCVar->pszName);
2550 /* dummy */
2551 return 0;
2552
2553 } /* FreeMCVar */
2554
2555
2556 /* ----------------------------------------------------------------------------
2557 FreeDataRec
2558
2559 Partial freeing (pdData is left in memory)
2560 */
FreeDataRec(PVOID pData,PVOID pUserInfo)2561 int FreeDataRec (PVOID pData, PVOID pUserInfo)
2562 {
2563 PDATAREC pDataRecord = (PDATAREC) pData;
2564
2565 /* do not free(pDataRecord->pdData); */
2566 free(pDataRecord->szDataName);
2567 free(pDataRecord);
2568 return 0;
2569
2570 } /* FreeDataRec */
2571
2572
2573 /* ----------------------------------------------------------------------------
2574 FreePrintRec
2575 */
FreePrintRec(PVOID pData,PVOID pUserInfo)2576 int FreePrintRec (PVOID pData, PVOID pUserInfo)
2577 {
2578 PPRINTREC pPrintRecord = (PPRINTREC) pData;
2579
2580 free(pPrintRecord->pdTimes);
2581 free(pPrintRecord->szOutputName);
2582 free(pPrintRecord);
2583 return 0;
2584
2585 } /* FreePrintRec */
2586
2587
2588 /* ----------------------------------------------------------------------------
2589 FreeOneLevel (recursive)
2590 */
FreeOneLevel(PLEVEL plevel)2591 void FreeOneLevel (PLEVEL plevel)
2592 {
2593 int n;
2594
2595 for (n = 0; n < plevel->iInstances; n++)
2596 if (plevel->pLevels[n] != NULL)
2597 FreeOneLevel(plevel->pLevels[n]);
2598
2599 FreeList(&plevel->plistVars, NULL, TRUE);
2600
2601 ForAllList(plevel->plistMCVars, &FreeMCVar, NULL);
2602 FreeList(&plevel->plistMCVars, NULL, TRUE);
2603
2604 ForAllList(plevel->plistLikes, &FreeMCVar, NULL);
2605 FreeList(&plevel->plistLikes, NULL, TRUE);
2606
2607 if (plevel->pexpt != NULL) {
2608 FreeList(&plevel->pexpt->plistParmMods, &FreeVarMod, FALSE);
2609
2610 POUTSPEC pos = &plevel->pexpt->os;
2611 free(pos->pszOutputNames);
2612 free(pos->phvar_out);
2613 free(pos->pcOutputTimes);
2614 free(pos->piCurrentOut);
2615 free(pos->prgdOutputTimes);
2616 for (n = 0; n < pos->nOutputs; n++)
2617 free(pos->prgdOutputVals[n]);
2618 free(pos->prgdOutputVals);
2619 free(pos->rgdDistinctTimes);
2620 ForAllList(pos->plistPrintRecs, &FreePrintRec, NULL);
2621 FreeList(&pos->plistPrintRecs, NULL, FALSE);
2622 free(pos->plistPrintRecs);
2623
2624 free(pos->pcData);
2625 free(pos->phvar_dat);
2626 free(pos->pszDataNames);
2627 for (n = 0; n < pos->nData; n++)
2628 free(pos->prgdDataVals[n]);
2629 free(pos->prgdDataVals);
2630 //ForAllList(pos->plistDataRecs, &FreeDataRec, NULL);
2631 //FreeList(&pos->plistDataRecs, NULL, FALSE);
2632 //free(pos->plistDataRecs);
2633
2634 free(plevel->pexpt);
2635 }
2636 if (plevel->nFixedVars > 0)
2637 free(plevel->rgpFixedVars);
2638 /* for (n = 0; n < plevel->nMCVars; n++) */
2639 /* FreeList(&plevel->rgpMCVars[n]->plistDependents, NULL, TRUE); */
2640 if (plevel->nMCVars > 0)
2641 free(plevel->rgpMCVars);
2642 //if (plevel->nLikes > 0) {
2643 //printf("plevel->nLikes %ld\n", plevel->nLikes);
2644 //for (n = 0; n < plevel->nLikes; n++) {
2645 //printf("plevel->rgpLikes[n]-> %lu\n", &plevel->rgpLikes[n]);
2646 /* clean array of MCVars */
2647 //if (&plevel->rgpLikes[n] != NULL)
2648 //FreeMCVar ((PVOID) (plevel->rgpLikes[n]), NULL);
2649 //}
2650 free(plevel->rgpLikes);
2651 //}
2652
2653 free(plevel);
2654
2655 } /* FreeOneLevel */
2656
2657
2658 /* ----------------------------------------------------------------------------
2659 ProcessWord
2660
2661 processes the word szLex.
2662
2663 This is the main loop of the interpreter. It is a big switch that
2664 recognizes keywords that are specifications and takes the
2665 appropriate action.
2666
2667 If the word szLex is not a keyword, ProcessWord() attempts to
2668 define a parameter specification.
2669 */
ProcessWord(PINPUTBUF pibIn,PSTR szLex,PSTR szEqn)2670 void ProcessWord (PINPUTBUF pibIn, PSTR szLex, PSTR szEqn)
2671 {
2672 int iErr = 0;
2673 int iKWCode, fContext;
2674 long i, iLB, iUB;
2675 PSTREQN szEqnU;
2676 PSTRLEX szTmp;
2677 PANALYSIS panal;
2678
2679 if (!pibIn || !szLex || !szLex[0] || !szEqn)
2680 return;
2681
2682 panal = (PANALYSIS) pibIn->pInfo;
2683
2684 iKWCode = GetKeywordCode(szLex, &fContext);
2685
2686 assert(panal->wContext != CN_END);
2687
2688 if ((iErr =
2689 (iKWCode /* Is a keyword */
2690 && !(fContext & panal->wContext)))) /* In invalid context */
2691 ReportError(pibIn, RE_BADCONTEXT, szLex, NULL);
2692
2693 else {
2694 switch (iKWCode) {
2695
2696 default: /* If a keyword is not found, try to get an assignment */
2697
2698 /* check if it is an array */
2699 iLB = iUB = -1;
2700 if (GetPunct(pibIn, szTmp, '[')) /* array found, read bounds */
2701 GetArrayBounds(pibIn, &iLB, &iUB);
2702 else
2703 if (!strcmp(szTmp, "=")) /* put it back... */
2704 pibIn->pbufCur--;
2705
2706 if (iUB == -1) { /* scalar */
2707 iErr = GetParmMod(pibIn, szLex);
2708 }
2709 else { /* array */
2710
2711 if (GetPunct(pibIn, szTmp, '=')) { /* read assignment */
2712 GetStatement(pibIn, szEqn);
2713 for (i = iLB; i < iUB; i++) {
2714 sprintf(szTmp, "%s_%ld", szLex, i); /* create names */
2715 UnrollEquation(pibIn, i, szEqn, szEqnU);
2716 iErr = GetParmMod2(pibIn, szTmp, szEqnU);
2717 if (iErr)
2718 break;
2719 }
2720 }
2721 else {
2722 ReportError(pibIn, RE_LEXEXPECTED | RE_FATAL, "= or [", NULL);
2723 }
2724 }
2725 break;
2726
2727 /* Otherwise process the following keywords */
2728
2729 case KM_PRINT:
2730 iErr = GetPrint(pibIn, szLex, &panal->pexpCurrent->os);
2731 break;
2732
2733 case KM_PRINTSTEP:
2734 iErr = GetPrintStep(pibIn, szLex, &panal->pexpCurrent->os);
2735 break;
2736
2737 case KM_EXPERIMENT:
2738 if (!(iErr = EGetPunct(pibIn, szTmp, CH_LBRACE)))
2739 NewExperiment(pibIn);
2740 break;
2741
2742 case KM_LEVEL:
2743 iErr = SetLevel(pibIn);
2744 break;
2745
2746 case KM_MCVARY:
2747 iErr = GetDistribSpec(pibIn, szLex, panal);
2748 break;
2749
2750 case KM_OUTPUTFILE:
2751 if (panal->szOutfilename)
2752 EatStatement(pibIn);
2753 else
2754 iErr = GetOutputFile(pibIn, szLex, panal);
2755 break;
2756
2757 case KM_DATA:
2758 iErr = GetData(pibIn, szLex, &panal->pexpCurrent->os);
2759 break;
2760
2761 case KM_INTEGRATE:
2762 iErr = GetIntegrate(pibIn, szLex, &panal->pexpCurrent->is);
2763 break;
2764
2765 case KM_MCMC:
2766 iErr = GetMCMCSpec(pibIn, panal->pexpCurrent);
2767 break;
2768
2769 case KM_OPTDESIGN:
2770 iErr = GetOptDSpec(pibIn, panal, szLex);
2771 break;
2772
2773 case KM_MONTECARLO:
2774 iErr = GetMonteCarloSpec(pibIn, panal, szLex);
2775 break;
2776
2777 case KM_SETPOINTS:
2778 iErr = GetSetPointsSpec(pibIn, panal, szLex);
2779 break;
2780
2781 case KM_SIMULATE:
2782 iErr = GetSimulate();
2783 break;
2784
2785 case KM_STARTTIME:
2786 iErr = GetStartTime(pibIn, panal->pexpCurrent);
2787 break;
2788
2789 case KM_SIMTYPE:
2790 iErr = GetSimType(pibIn);
2791 break;
2792
2793 case KM_TEMPERATURE:
2794 iErr = GetPerks(pibIn, szLex, &panal->gd);
2795 break;
2796
2797 case KM_END:
2798 panal->wContext = CN_END;
2799 break;
2800
2801 } /* switch */
2802 } /* else */
2803
2804 if (iErr)
2805 EatStatement(pibIn);
2806
2807 } /* ProcessWord */
2808
2809
2810 /* ----------------------------------------------------------------------------
2811 ReadAnalysis
2812
2813 Core routine for input file parsing. Called from sim.c in particular.
2814 */
ReadAnalysis(PINPUTBUF pibIn)2815 BOOL ReadAnalysis (PINPUTBUF pibIn)
2816 {
2817 PSTRLEX szLex; /* Lex elem of MAX_LEX length */
2818 PSTREQN szEqn; /* Equation buffer of MAX_EQN length */
2819 int iLexType;
2820
2821 BOOL bReturn = TRUE;
2822 PANALYSIS panal;
2823
2824 if (!pibIn) return(FALSE);
2825
2826 panal = (PANALYSIS) pibIn->pInfo;
2827 panal->iDepth = panal->iCurrentDepth = panal->iInstances = 0;
2828 panal->mc.plistMCVars = InitList();
2829
2830 do {
2831 /* State machine for parsing syntax */
2832 NextLex(pibIn, szLex, &iLexType);
2833
2834 switch (iLexType) {
2835
2836 case LX_NULL:
2837 if (panal->wContext != CN_GLOBAL)
2838 ReportError(pibIn, RE_WARNING, NULL, "Unexpected end of file");
2839
2840 if (panal->wContext == CN_EXPERIMENT)
2841 bReturn &= EndExperiment(pibIn, panal);
2842 panal->wContext = CN_END;
2843 break;
2844
2845 case LX_IDENTIFIER:
2846 ProcessWord(pibIn, szLex, szEqn);
2847 break;
2848
2849 case LX_PUNCT:
2850 if (szLex[0] == CH_STMTTERM)
2851 break;
2852 else if (szLex[0] == CH_RBRACE) {
2853 if (panal->wContext & CN_EXPERIMENT) {
2854 bReturn &= EndExperiment(pibIn, panal);
2855 break;
2856 }
2857 else {
2858 bReturn &= EndLevel(panal);
2859 break;
2860 }
2861 }
2862 else
2863 if (szLex[0] == CH_COMMENT) {
2864 SkipComment(pibIn);
2865 break;
2866 }
2867
2868 /* else -- fall through! */
2869
2870 default:
2871 ReportError(pibIn, RE_UNEXPECTED, szLex, "* Ignoring");
2872 break;
2873
2874 case LX_INTEGER:
2875 case LX_FLOAT:
2876 ReportError(pibIn, RE_UNEXPNUMBER, szLex, "* Ignoring");
2877 break;
2878
2879 } /* switch */
2880
2881 } while (panal->wContext != CN_END
2882 && (*pibIn->pbufCur || FillBuffer(pibIn) != EOF));
2883
2884 if(panal->iCurrentDepth != 0)
2885 ReportError(pibIn, RE_OPENLEVEL | RE_FATAL, "ReadAnalysis", NULL);
2886 return(bReturn);
2887
2888 } /* ReadAnalysis */
2889
2890 /* End */
2891