1 /* modo.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 Output the model.c file.
21
22 Modified on 07/07/97 by FB - added yourcode.h to WriteIncludes().
23
24 */
25
26 #include <stdio.h>
27 #include <stdlib.h>
28 #include <string.h>
29 #include <ctype.h>
30 #include <assert.h>
31 #include <time.h>
32
33 #include "mod.h"
34 #include "lexfn.h"
35 #include "lexerr.h"
36 #include "modi.h"
37 #include "modd.h"
38 #include "modo.h"
39
40
41 /* Global Variables */
42
43 static char *vszModelFilename = NULL;
44 static char *vszModGenName = NULL;
45
46 static char vszModelArrayName[] = "vrgModelVars";
47 static char vszInputArrayName[] = "vrgInputs";
48 extern char vszHasInitializer[]; /* decl'd in modd.c */
49
50 PVMMAPSTRCT vpvmGloVarList;
51
52 char *vszIFNTypes[] = { /* Must match defines in lexfn.h */
53 "IFN_NULL /* ?? */",
54 "IFN_CONSTANT",
55 "IFN_PERDOSE",
56 "IFN_PERRATE",
57 "IFN_PEREXP",
58 "IFN_NDOSES"
59 }; /* vszIFNTypes[] = */
60
61 int vnStates, vnOutputs, vnInputs, vnParms, vnModelVars;
62
63 BOOL bForR = FALSE;
64 BOOL bForInits = FALSE;
65 BOOL bDelay = TRUE; /* R specific code */
66
67
68 /* ----------------------------------------------------------------------------
69 ForAllVar
70
71 Takes a pfile, a pvm list head and a callback function which is called
72 for all variables in the list if hType == ALL_VARS or only for
73 only those that match hType otherwise.
74 */
ForAllVar(PFILE pfile,PVMMAPSTRCT pvm,PFI_CALLBACK pfiFunc,HANDLE hType,PVOID pinfo)75 int ForAllVar (PFILE pfile, PVMMAPSTRCT pvm, PFI_CALLBACK pfiFunc,
76 HANDLE hType, PVOID pinfo)
77 {
78 int iTotal = 0;
79
80 while (pvm) {
81 if (hType == ALL_VARS /* Do for all ...*/
82 || TYPE(pvm) == hType) { /* ..or do only for vars of hType */
83 if (pfiFunc)
84 iTotal += (*pfiFunc) (pfile, pvm, pinfo);
85 else
86 iTotal++; /* No func! just count */
87 }
88 pvm = pvm->pvmNextVar;
89 } /* while */
90
91 return(iTotal);
92
93 } /* ForAllVar */
94
95
96 /* ----------------------------------------------------------------------------
97 CountOneDecl
98
99 Counts declarations for variables. Note that each variable may only
100 be declared and assigned *once*. This is a problem and should be fixed.
101
102 Callback for ForAllVar().
103 */
CountOneDecl(PFILE pfile,PVMMAPSTRCT pvm,PVOID pInfo)104 int CountOneDecl (PFILE pfile, PVMMAPSTRCT pvm, PVOID pInfo)
105 {
106 if (pvm->szEqn != vszHasInitializer) { /* These appear later */
107 return(1);
108 }
109 return(0);
110
111 } /* CountOneDecl */
112
113
114 /* ----------------------------------------------------------------------------
115 */
WritebDelays(PFILE pfile,BOOL bDelays)116 void WritebDelays (PFILE pfile, BOOL bDelays)
117 {
118
119 fprintf(pfile, "\nBOOL bDelays = %d;\n", bDelays);
120
121 } /* WritebDelays */
122
123
124 /* ----------------------------------------------------------------------------
125 WriteOneName
126
127 Writes a name and nominal value for the given variable of pvm.
128
129 Callback for ForAllVar().
130 */
WriteOneName(PFILE pfile,PVMMAPSTRCT pvm,PVOID pInfo)131 int WriteOneName (PFILE pfile, PVMMAPSTRCT pvm, PVOID pInfo)
132 {
133 if (pvm->szEqn != vszHasInitializer) { /* These appear later */
134 if (bForR) {
135 if (TYPE(pvm) == ID_OUTPUT)
136 fprintf(pfile, " \"%s", pvm->szName);
137 else
138 fprintf(pfile, " %s", pvm->szName);
139
140 if (TYPE(pvm) != ID_INPUT) {
141 if (TYPE(pvm) == ID_OUTPUT)
142 fprintf(pfile, "\",\n");
143 else
144 fprintf(pfile, " = %s,\n", (pvm->szEqn ? pvm->szEqn : "0.0"));
145 }
146 else
147 fprintf(pfile, " (forcing function)\n");
148 }
149 else { /* not for R */
150 fprintf(pfile, " %s", pvm->szName);
151 if (TYPE(pvm) != ID_INPUT)
152 fprintf(pfile, " %s %s;\n", (pvm->szEqn ? "=" : "->"),
153 (pvm->szEqn ? pvm->szEqn : "0.0"));
154 else
155 fprintf(pfile, " (is a function)\n");
156 }
157 return(1);
158 } /* if */
159
160 return(0);
161
162 } /* WriteOneName */
163
164
165 /* ----------------------------------------------------------------------------
166 WriteHeader
167
168 That it does do dashingly. Includes a nice comment at the
169 beginning (hopefully!) of the file replete with filename,
170 timestamp, and a list of model variables.
171 */
WriteHeader(PFILE pfile,PSTR szName,PVMMAPSTRCT pvmGlo)172 void WriteHeader (PFILE pfile, PSTR szName, PVMMAPSTRCT pvmGlo)
173 {
174 time_t ttTime;
175
176 time(&ttTime);
177
178 if (fprintf(pfile, "/* %s\n", szName) < 0)
179 ReportError(NULL, RE_CANNOTOPEN | RE_FATAL,
180 szName, "...in WriteHeader ()");
181
182 fprintf(pfile, " ___________________________________________________\n\n");
183 fprintf(pfile, " Model File: %s\n\n", vszModelFilename);
184 fprintf(pfile, " Date: %s\n", ctime(&ttTime));
185 fprintf(pfile, " Created by: \"%s %s\"\n", vszModGenName, VSZ_VERSION);
186 fprintf(pfile, " -- a model preprocessor by Don Maszle\n");
187 fprintf(pfile, " ___________________________________________________\n\n");
188
189 fprintf(pfile, " " VSZ_COPYRIGHT "\n");
190
191 fprintf(pfile, "\n Model calculations for compartmental model:\n\n");
192
193 if (vnStates == 1) fprintf(pfile, " 1 State:\n");
194 else fprintf(pfile, " %d States:\n", vnStates);
195 ForAllVar(pfile, pvmGlo, &WriteOneName, ID_STATE, NULL);
196
197 if (vnOutputs == 1) fprintf(pfile, "\n 1 Output:\n");
198 else fprintf(pfile, "\n %d Outputs:\n", vnOutputs);
199 ForAllVar(pfile, pvmGlo, &WriteOneName, ID_OUTPUT, NULL);
200
201 if (vnInputs == 1) fprintf(pfile, "\n 1 Input:\n");
202 else fprintf(pfile, "\n %d Inputs:\n", vnInputs);
203 ForAllVar(pfile, pvmGlo, &WriteOneName, ID_INPUT, NULL);
204
205 if (vnParms == 1) fprintf(pfile, "\n 1 Parameter:\n");
206 else fprintf(pfile, "\n %d Parameters:\n", vnParms);
207 ForAllVar(pfile, pvmGlo, &WriteOneName, ID_PARM, NULL);
208
209 fprintf(pfile, "*/\n\n");
210
211 } /* WriteHeader */
212
213
214 /* ----------------------------------------------------------------------------
215 */
WriteIncludes(PFILE pfile)216 void WriteIncludes (PFILE pfile)
217 {
218 fprintf(pfile, "\n#include <stdlib.h>\n");
219 fprintf(pfile, "#include <stdio.h>\n");
220 fprintf(pfile, "#include <math.h>\n");
221 fprintf(pfile, "#include <string.h>\n");
222 fprintf(pfile, "#include <float.h>\n");
223
224 fprintf(pfile, "#include \"modelu.h\"\n");
225 fprintf(pfile, "#include \"random.h\"\n");
226 fprintf(pfile, "#include \"yourcode.h\"\n");
227
228 } /* WriteIncludes */
229
230
231 /* ----------------------------------------------------------------------------
232 WriteOneDecl
233
234 Write one global or local declaration. Callback for ForAllVar().
235 */
WriteOneDecl(PFILE pfile,PVMMAPSTRCT pvm,PVOID pInfo)236 int WriteOneDecl (PFILE pfile, PVMMAPSTRCT pvm, PVOID pInfo)
237 {
238 assert(TYPE(pvm) != ID_INPUT);
239 assert(TYPE(pvm) != ID_OUTPUT);
240 assert(TYPE(pvm) != ID_STATE);
241
242 if (TYPE(pvm) > ID_PARM) {
243 fprintf(pfile, " /* local */ ");
244 }
245
246 fprintf(pfile, "double %s;\n", pvm->szName);
247
248 return(1);
249
250 } /* WriteOneDecl */
251
252
253 /* ----------------------------------------------------------------------------
254 */
WriteOneIndexDefine(PFILE pfile,PVMMAPSTRCT pvm,PVOID pInfo)255 int WriteOneIndexDefine (PFILE pfile, PVMMAPSTRCT pvm, PVOID pInfo)
256 {
257 if (pvm->szEqn != vszHasInitializer) { /* These appear later */
258 fprintf(pfile, "#define ");
259 WriteIndexName(pfile, pvm);
260 if (INDEX(pvm))
261 fprintf(pfile, " 0x%05lx\n", INDEX(pvm));
262 else
263 fprintf(pfile, " 0x00000\n");
264
265 return 1;
266 } /* if */
267
268 return 0;
269
270 } /* WriteOneIndexDefine */
271
272
273 /* ----------------------------------------------------------------------------
274 */
WriteDecls(PFILE pfile,PVMMAPSTRCT pvmGlo)275 void WriteDecls (PFILE pfile, PVMMAPSTRCT pvmGlo)
276 {
277 fprintf(pfile, "\n\n/*----- Indices to Global Variables */\n");
278
279 fprintf(pfile, "\n/* Model variables: States and other outputs */\n");
280 ForAllVar(pfile, pvmGlo, &WriteOneIndexDefine, ID_STATE, NULL);
281 ForAllVar(pfile, pvmGlo, &WriteOneIndexDefine, ID_OUTPUT, NULL);
282
283 fprintf(pfile, "\n/* Inputs */\n");
284 ForAllVar(pfile, pvmGlo, &WriteOneIndexDefine, ID_INPUT, NULL);
285 fprintf(pfile, "\n/* Parameters */\n");
286 ForAllVar(pfile, pvmGlo, &WriteOneIndexDefine, ID_PARM, NULL);
287
288 fprintf(pfile, "\n\n/*----- Global Variables */\n");
289
290 fprintf(pfile, "\n/* For export. Keep track of who we are. */\n");
291 fprintf(pfile, "char szModelDescFilename[] = \"%s\";\n", vszModelFilename);
292 fprintf(pfile, "char szModelSourceFilename[] = __FILE__;\n");
293 fprintf(pfile, "char szModelGenAndVersion[] = \"%s %s\";\n",
294 vszModGenName, VSZ_VERSION);
295
296 fprintf(pfile, "\n/* Externs */\n");
297 fprintf(pfile, "extern BOOL vbModelReinitd;\n");
298
299 fprintf(pfile, "\n/* Model Dimensions */\n");
300 fprintf(pfile, "int vnStates = %d;\n", vnStates);
301 fprintf(pfile, "int vnOutputs = %d;\n", vnOutputs);
302 fprintf(pfile, "int vnModelVars = %d;\n", vnModelVars);
303 fprintf(pfile, "int vnInputs = %d;\n", vnInputs);
304 fprintf(pfile, "int vnParms = %d;\n", vnParms);
305
306 fprintf(pfile, "\n/* States and Outputs*/\n");
307 fprintf(pfile, "double %s[%d];\n", vszModelArrayName, vnModelVars);
308
309 fprintf(pfile, "\n/* Inputs */\n");
310 /* if vnInputs is zero put a dummy 1 for array size */
311 fprintf(pfile, "IFN %s[%d];\n", vszInputArrayName,
312 (vnInputs > 0 ? vnInputs : 1));
313
314 fprintf(pfile, "\n/* Parameters */\n");
315 ForAllVar(pfile, pvmGlo, &WriteOneDecl, ID_PARM, NULL);
316
317 } /* WriteDecls */
318
319
320 /* ----------------------------------------------------------------------------
321 GetName
322
323 returns a string to the name of the pvm variable given. The
324 string is static and must be used immediately or copied. It will
325 be changed on the next call of this function.
326
327 szModelVarName and szDerivName are names to be used for
328 state variables and derivatives arrays, resp.
329
330 The name is determined by hType if hType is non-NULL. If hType is
331 NULL, then the type is taken to be the hType field of pvm.
332 */
GetName(PVMMAPSTRCT pvm,PSTR szModelVarName,PSTR szDerivName,HANDLE hType)333 PSTR GetName (PVMMAPSTRCT pvm, PSTR szModelVarName, PSTR szDerivName,
334 HANDLE hType)
335 {
336 static PSTRLEX vszVarName;
337 HANDLE hTypeToUse = (hType ? hType : TYPE(pvm));
338
339 switch (hTypeToUse) {
340
341 case ID_INPUT:
342 if (bForR)
343 sprintf(vszVarName, "%s", pvm->szName);
344 else
345 sprintf(vszVarName, "vrgInputs[ID_%s]", pvm->szName);
346 break;
347
348 case ID_STATE:
349 if (bForR) {
350 if (bForInits)
351 sprintf(vszVarName, "%s", pvm->szName);
352 else
353 sprintf(vszVarName, "y[ID_%s]", pvm->szName);
354 }
355 else {
356 if (szModelVarName)
357 sprintf(vszVarName, "%s[ID_%s]", szModelVarName, pvm->szName);
358 else
359 sprintf(vszVarName, "vrgModelVars[ID_%s]", pvm->szName);
360 }
361 break;
362
363 case ID_OUTPUT:
364 if (bForR)
365 sprintf(vszVarName, "yout[ID_%s]", pvm->szName);
366 else {
367 if (szModelVarName)
368 sprintf(vszVarName, "%s[ID_%s]", szModelVarName, pvm->szName);
369 else
370 sprintf(vszVarName, "vrgModelVars[ID_%s]", pvm->szName);
371 }
372 break;
373
374 case ID_DERIV:
375 assert(szDerivName);
376 if (bForR)
377 sprintf(vszVarName, "ydot[ID_%s]", pvm->szName);
378 else
379 sprintf(vszVarName, "%s[ID_%s]", szDerivName, pvm->szName);
380 break;
381
382 default: /* Parms and local variables */
383 sprintf(vszVarName, "%s", pvm->szName);
384 break;
385 } /* switch */
386
387 return(vszVarName);
388
389 } /* GetName */
390
391
392 /* ----------------------------------------------------------------------------
393 WriteOneVMEntry
394
395 prints a single entry for the variable map symbol table.
396
397 Callback for ForAllList.
398 */
WriteOneVMEntry(PFILE pfile,PVMMAPSTRCT pvm,PVOID pInfo)399 int WriteOneVMEntry (PFILE pfile, PVMMAPSTRCT pvm, PVOID pInfo)
400 {
401 int iType = TYPE(pvm);
402
403 if (!pvm) {
404 fprintf(pfile, " {\"\", NULL, 0} /* End flag */\n");
405 return 0;
406 }
407
408 assert(iType && \
409 iType != ID_LOCALDYN && iType != ID_LOCALSCALE && iType != ID_LOCALJACOB);
410
411 if (pvm->szEqn != vszHasInitializer) { /* These appear later */
412 fprintf(pfile, " {\"%s\", (PVOID) &%s", pvm->szName,
413 GetName(pvm, vszModelArrayName, NULL, ID_NULL));
414
415 fprintf(pfile, ", ID_%s | ID_%s},\n",
416 (iType == ID_PARM ? "PARM"
417 : (iType == ID_INPUT ? "INPUT"
418 : (iType == ID_OUTPUT ? "OUTPUT"
419 : "STATE"))), pvm->szName);
420
421 return(1);
422
423 } /* if */
424
425 return 0;
426
427 } /* WriteOneVMEntry */
428
429
430 /* ----------------------------------------------------------------------------
431 */
WriteVarMap(PFILE pfile,PVMMAPSTRCT pvmGlo)432 void WriteVarMap (PFILE pfile, PVMMAPSTRCT pvmGlo)
433 {
434 fprintf(pfile, "\n\n/*----- Global Variable Map */\n\n");
435 fprintf(pfile, "VMMAPSTRCT vrgvmGlo[] = {\n");
436 ForAllVar(pfile, pvmGlo, &WriteOneVMEntry, ID_STATE, NULL);
437 ForAllVar(pfile, pvmGlo, &WriteOneVMEntry, ID_OUTPUT, NULL);
438 ForAllVar(pfile, pvmGlo, &WriteOneVMEntry, ID_INPUT, NULL);
439 ForAllVar(pfile, pvmGlo, &WriteOneVMEntry, ID_PARM, NULL);
440 WriteOneVMEntry(pfile, NULL, NULL); /* Include end flag */
441 fprintf(pfile, "}; /* vrgpvmGlo[] */\n");
442
443 } /* WriteVarMap */
444
445
446 /* ----------------------------------------------------------------------------
447 */
WriteOneInit(PFILE pfile,PVMMAPSTRCT pvm,PVOID pInfo)448 int WriteOneInit (PFILE pfile, PVMMAPSTRCT pvm, PVOID pInfo)
449 {
450 IFN ifnNull = {IFN_CONSTANT}; /* Init other fields to zero */
451 PSTR szVarName = GetName(pvm, NULL, NULL, ID_NULL);
452
453 if ((pvm->szEqn != vszHasInitializer && /* These appear later */
454 TYPE(pvm) <= ID_PARM) || /* No Locals! */
455 (TYPE(pvm) == ID_INLINE)) {
456
457 if (TYPE(pvm) == ID_INPUT) {
458 PIFN pifn = (PIFN) pvm->szEqn;
459
460 if (!pifn) /* No eqn, init to zero */
461 pifn = &ifnNull;
462
463 fprintf(pfile, " %s.iType = %s;\n",
464 szVarName, vszIFNTypes[pifn->iType]);
465 fprintf(pfile, " %s.dTStartPeriod = 0;\n", szVarName);
466 fprintf(pfile, " %s.bOn = FALSE;\n", szVarName);
467 fprintf(pfile, " %s.dMag = %f;\n", szVarName, pifn->dMag);
468 fprintf(pfile, " %s.dT0 = %f;\n", szVarName, pifn->dT0);
469 fprintf(pfile, " %s.dTexp = %f;\n", szVarName, pifn->dTexp);
470 fprintf(pfile, " %s.dDecay = %f;\n", szVarName, pifn->dDecay);
471 fprintf(pfile, " %s.dTper = %f;\n", szVarName, pifn->dTper);
472
473 fprintf(pfile, " %s.hMag = %#lx;\n", szVarName, pifn->hMag);
474 fprintf(pfile, " %s.hT0 = %#lx;\n", szVarName, pifn->hT0);
475 fprintf(pfile, " %s.hTexp = %#lx;\n", szVarName, pifn->hTexp);
476 fprintf(pfile, " %s.hDecay = %#lx;\n", szVarName, pifn->hDecay);
477 fprintf(pfile, " %s.hTper = %#lx;\n", szVarName, pifn->hTper);
478
479 fprintf(pfile, " %s.dVal = 0.0;\n", szVarName);
480 fprintf(pfile, " %s.nDoses = 0;\n", szVarName);
481 } /* if */
482 else {
483 if (TYPE(pvm) == ID_INLINE) /* write out just the equation */
484 fprintf(pfile, "\n%s\n", pvm->szEqn);
485 else
486 fprintf(pfile, " %s = %s;\n", szVarName,
487 (pvm->szEqn ? pvm->szEqn : "0.0"));
488 } /* else */
489
490 return(1);
491
492 } /* if */
493
494 return(0);
495
496 } /* WriteOneInit */
497
498
499 /* ----------------------------------------------------------------------------
500 WriteInitModel
501
502 Writes the routine to initialize the model variables.
503 */
WriteInitModel(PFILE pfile,PVMMAPSTRCT pvmGlo)504 void WriteInitModel (PFILE pfile, PVMMAPSTRCT pvmGlo)
505 {
506 fprintf(pfile, "\n\n/*----- InitModel\n");
507 fprintf(pfile, " Should be called to initialize model variables at\n");
508 fprintf(pfile, " the beginning of experiment before reading\n");
509 fprintf(pfile, " variants from the simulation spec file.\n*/\n\n");
510 fprintf(pfile, "void InitModel(void)\n{\n");
511 fprintf(pfile, " /* Initialize things in the order that they appear in\n"
512 " model definition file so that dependencies are\n"
513 " handled correctly. */\n\n");
514 ForAllVar(pfile, pvmGlo, &WriteOneInit, ALL_VARS, NULL);
515 fprintf(pfile, "\n vbModelReinitd = TRUE;\n\n");
516 fprintf(pfile, "} /* InitModel */\n\n\n");
517
518 } /* WriteInitModel */
519
520
521 /* ----------------------------------------------------------------------------
522 TranslateID
523
524 Writes an equation id if it declared and in a valid context.
525 Errors are reported if necessary.
526 */
527
TranslateID(PINPUTBUF pibDum,PFILE pfile,PSTR szLex,int iEqType)528 void TranslateID (PINPUTBUF pibDum, PFILE pfile, PSTR szLex, int iEqType)
529 {
530 int iKWCode, fContext;
531 long iLowerB, iUpperB;
532
533 iKWCode = GetKeywordCode(szLex, &fContext);
534
535 switch (iKWCode) {
536
537 case KM_DXDT: {
538 int iArg = LX_IDENTIFIER;
539 PVMMAPSTRCT pvm = NULL;
540
541 if (GetFuncArgs(pibDum, 1, &iArg, szLex, &iLowerB, &iUpperB)
542 && (pvm = GetVarPTR(vpvmGloVarList, szLex))
543 && TYPE(pvm) == ID_STATE)
544 fprintf(pfile, "%s", GetName(pvm, NULL, "rgDerivs", ID_DERIV));
545 else
546 ReportError(pibDum, RE_BADSTATE | RE_FATAL, (pvm? szLex : NULL), NULL);
547 } /* KM_DXDT block */
548 break;
549
550 case KM_NULL: {
551 PVMMAPSTRCT pvm = GetVarPTR(vpvmGloVarList, szLex);
552
553 /* Handle undeclared ids */
554 if (!pvm) {
555 if ((iEqType == KM_DYNAMICS ||
556 iEqType == KM_SCALE ||
557 iEqType == KM_CALCOUTPUTS) &&
558 !(strcmp(szLex, VSZ_TIME) &&
559 strcmp(szLex, VSZ_TIME_SBML)))
560 /* If this is the time variable, convert to the correct formal arg */
561 fprintf(pfile, "(*pdTime)");
562 else
563 /* otherwise output id exactly as is */
564 fprintf(pfile, "%s", szLex);
565 } /* if */
566
567 else {
568 if (iEqType == KM_SCALE)
569 /* fixed by FB - 1 mar 98. This prints a "vrgModelVars" */
570 fprintf(pfile, "%s", GetName(pvm, NULL, NULL, ID_NULL));
571 else
572 fprintf(pfile, "%s", GetName(pvm, "rgModelVars", NULL, ID_NULL));
573
574 if ((TYPE(pvm) == ID_INPUT) && (!bForR))
575 fprintf(pfile, ".dVal"); /* Use current value */
576 } /* else */
577
578 } /* KM_NULL: */
579 break;
580
581 default: /* No other keywords allowed in equations */
582
583 /* Allow for C keywords here, including math library functions */
584 ReportError(pibDum, RE_BADCONTEXT | RE_FATAL, szLex, NULL);
585 break;
586
587 } /* switch */
588
589 } /* TranslateID */
590
591
592 /* ----------------------------------------------------------------------------
593 TranslateEquation
594
595 Writes an equation to the output file substituting model variable
596 names, inputs names, derivative names, etc.
597
598 Tries to do some hack formatting.
599 */
600
TranslateEquation(PFILE pfile,PSTR szEqn,long iEqType)601 void TranslateEquation (PFILE pfile, PSTR szEqn, long iEqType)
602 {
603 INPUTBUF ibDum;
604 PINPUTBUF pibDum = &ibDum;
605 PSTRLEX szLex;
606 PVMMAPSTRCT pvm = NULL;
607 int iType;
608 BOOL bDelayCall = FALSE;
609
610 MakeStringBuffer(NULL, pibDum, szEqn);
611
612 NextLex(pibDum, szLex, &iType);
613 if (!iType) {
614 fprintf(pfile, "0.0; /* NULL EQN!?? */");
615 return;
616 } /* if */
617
618 do {
619 if (iType == LX_IDENTIFIER) { /* Process Identifier */
620 if (bDelayCall) {
621 /* do not translate the 1st param of CalcDelay but check it */
622 pvm = GetVarPTR(vpvmGloVarList, szLex);
623 if ( ( bForR && ((pvm && (TYPE(pvm) == ID_STATE)))) ||
624 (!bForR && ((pvm && (TYPE(pvm) == ID_STATE)) ||
625 (TYPE(pvm) == ID_OUTPUT)))) {
626 fprintf(pfile, "ID_%s", szLex);
627 fprintf(pfile, ", (*pdTime)"); /* add the automatic time variable */
628 bDelayCall = FALSE; /* turn delay context off */
629 }
630 else
631 ReportError(pibDum, RE_LEXEXPECTED | RE_FATAL,
632 (bForR ? "state" : "state or output") , NULL);
633 }
634 else
635 TranslateID(pibDum, pfile, szLex, iEqType);
636 }
637 else {
638 if ((iType == LX_EQNPUNCT ||
639 iType == LX_PUNCT) &&
640 (*(szLex) == CH_COMMENT)) {
641 while (*pibDum->pbufCur && *pibDum->pbufCur != CH_EOLN)
642 pibDum->pbufCur++;
643
644 fprintf(pfile, "\n");
645 }
646 else /* Spew everything else */
647 fprintf(pfile, "%s", szLex);
648 }
649
650 if (!bDelayCall) { /* check delay context */
651 bDelayCall = (!strcmp("CalcDelay", szLex));
652 bDelay = bDelay || bDelayCall;
653 }
654
655 fprintf(pfile, " ");
656 NextLex(pibDum, szLex, &iType);
657
658 } while (iType);
659
660 if (bForR && bForInits)
661 fprintf(pfile, "\n");
662 else
663 fprintf(pfile, ";\n");
664
665 } /* TranslateEquation */
666
667
668 /* ----------------------------------------------------------------------------
669 WriteOneEquation
670
671 Writes one equation of a function definition. Uses the hType
672 flag for spacing set by DefineEquation().
673
674 Callback function for ForAllVar().
675 */
WriteOneEquation(PFILE pfile,PVMMAPSTRCT pvm,PVOID pInfo)676 int WriteOneEquation (PFILE pfile, PVMMAPSTRCT pvm, PVOID pInfo)
677 {
678 long iType = (long) pInfo;
679
680 if (pvm->hType & ID_SPACEFLAG) /* Flag for space between equations */
681 fprintf(pfile, "\n");
682
683 switch (iType) {
684 default:
685 case KM_SCALE: /* Scale Global Names */
686
687 /* Inputs not allowed anymore in Scale section - FB 7/12/96 */
688 if (TYPE(pvm) == ID_INPUT) {
689 printf("Error: input '%s' used in Scale context.\n",
690 pvm->szName);
691 exit(0);
692 }
693
694 if (TYPE(pvm) != ID_INLINE) { /* do not write "Inline" */
695 if (bForR && bForInits && TYPE(pvm) == ID_STATE)
696 fprintf(pfile, " Y[\"%s\"] <- ",
697 GetName (pvm, NULL, NULL, ID_NULL));
698 else
699 fprintf(pfile, " %s = ", GetName (pvm, NULL, NULL, ID_NULL));
700 }
701 break;
702
703 case KM_CALCOUTPUTS:
704 case KM_DYNAMICS:
705 case KM_JACOB:
706 case KM_EVENTS:
707 case KM_ROOTS:
708 if (TYPE(pvm) != ID_INLINE) /* do not write "Inline" */
709 fprintf(pfile, " %s = ", GetName (pvm, "rgModelVars", "rgDerivs",
710 ID_NULL));
711 break;
712
713 } /* switch */
714
715 if (TYPE(pvm) == ID_INLINE) /* write out the equation */
716 fprintf(pfile, "\n%s\n", pvm->szEqn);
717 else
718 TranslateEquation(pfile, pvm->szEqn, iType);
719
720 return 1;
721
722 } /* WriteOneEquation */
723
724
725 /* ----------------------------------------------------------------------------
726 WriteCalcDeriv
727
728 Writes the CalcDeriv() function. Writes dynamics equations in
729 the order they appeared in the model definition file.
730 */
731
WriteCalcDeriv(PFILE pfile,PVMMAPSTRCT pvmGlo,PVMMAPSTRCT pvmDyn)732 void WriteCalcDeriv (PFILE pfile, PVMMAPSTRCT pvmGlo, PVMMAPSTRCT pvmDyn)
733 {
734 if (!pvmDyn)
735 printf("No Dynamics{} equations.\n\n");
736
737 fprintf(pfile, "/*----- Dynamics section */\n\n");
738 fprintf(pfile, "void CalcDeriv (double rgModelVars[], ");
739 fprintf(pfile, "double rgDerivs[], PDOUBLE pdTime)\n{\n");
740
741 ForAllVar(pfile, pvmGlo, &WriteOneDecl, ID_LOCALDYN, NULL);
742
743 fprintf(pfile, "\n CalcInputs (pdTime); /* Get new input vals */\n\n");
744
745 ForAllVar(pfile, pvmDyn, &WriteOneEquation, ALL_VARS, (PVOID) KM_DYNAMICS);
746
747 fprintf(pfile, "\n} /* CalcDeriv */\n\n\n");
748
749 } /* WriteCalcDeriv */
750
751
752
753 /* ----------------------------------------------------------------------------
754 WriteScale
755 */
WriteScale(PFILE pfile,PVMMAPSTRCT pvmGlo,PVMMAPSTRCT pvmScale)756 void WriteScale (PFILE pfile, PVMMAPSTRCT pvmGlo, PVMMAPSTRCT pvmScale)
757 {
758 if (!pvmScale)
759 printf("No Scale{} equations. Null function defined.\n\n");
760
761 fprintf(pfile, "/*----- Model scaling */\n\n");
762 fprintf(pfile, "void ScaleModel (PDOUBLE pdTime)\n");
763 fprintf(pfile, "{\n");
764 ForAllVar(pfile, pvmGlo, &WriteOneDecl, ID_LOCALSCALE, NULL);
765 ForAllVar(pfile, pvmScale, &WriteOneEquation, ALL_VARS, (PVOID) KM_SCALE);
766 fprintf(pfile, "\n} /* ScaleModel */\n\n\n");
767
768 } /* WriteScale */
769
770
771 /* ----------------------------------------------------------------------------
772 WriteCalcJacob
773 */
WriteCalcJacob(PFILE pfile,PVMMAPSTRCT pvmGlo,PVMMAPSTRCT pvmJacob)774 void WriteCalcJacob (PFILE pfile, PVMMAPSTRCT pvmGlo, PVMMAPSTRCT pvmJacob)
775 {
776 fprintf(pfile, "/*----- Jacobian calculations */\n\n");
777 fprintf(pfile, "void CalcJacob (PDOUBLE pdTime, double rgModelVars[],\n");
778 fprintf(pfile, " long column, double rgdJac[])\n");
779 fprintf(pfile, "{\n");
780 ForAllVar(pfile, pvmGlo, &WriteOneDecl, ID_LOCALJACOB, NULL);
781 ForAllVar(pfile, pvmJacob, &WriteOneEquation, ALL_VARS, (PVOID) KM_JACOB);
782 fprintf(pfile, "\n} /* CalcJacob */\n\n\n");
783
784 } /* WriteCalcJacob */
785
786
787 /* ----------------------------------------------------------------------------
788 WriteCalcOutputs
789 */
WriteCalcOutputs(PFILE pfile,PVMMAPSTRCT pvmGlo,PVMMAPSTRCT pvmCalcOut)790 void WriteCalcOutputs (PFILE pfile, PVMMAPSTRCT pvmGlo, PVMMAPSTRCT pvmCalcOut)
791 {
792 if (!pvmCalcOut)
793 printf("No CalcOutputs{} equations. Null function defined.\n\n");
794
795 fprintf(pfile, "/*----- Outputs calculations */\n\n");
796 fprintf(pfile, "void CalcOutputs (double rgModelVars[], ");
797 fprintf(pfile, "double rgDerivs[], PDOUBLE pdTime)\n{\n");
798 ForAllVar(pfile, pvmGlo, &WriteOneDecl, ID_LOCALCALCOUT, NULL);
799 ForAllVar(pfile, pvmCalcOut, &WriteOneEquation, ALL_VARS,
800 (PVOID) KM_CALCOUTPUTS);
801 fprintf(pfile, "\n} /* CalcOutputs */\n\n\n");
802
803 } /* WriteCalcOutputs */
804
805
806 /* ----------------------------------------------------------------------------
807 IndexOneVar
808
809 ORs the index passed through the info pointer, a PINT,
810 and increments the value of the pointer.
811
812 The handle is corrected for use in the simulation. For this, the
813 TYPE part of the handle is shifted left 4 bits.
814
815 Callback function for ForAllVar().
816 */
IndexOneVar(PFILE pfile,PVMMAPSTRCT pvm,PVOID pInfo)817 int IndexOneVar (PFILE pfile, PVMMAPSTRCT pvm, PVOID pInfo)
818 {
819 if (pvm->szEqn != vszHasInitializer) { /* Handled later */
820 pvm->hType |= (*((PINT) pInfo))++;
821 return 1;
822 } /* if */
823 return 0;
824
825 } /* IndexOneVar */
826
827
828 /* ----------------------------------------------------------------------------
829 IndexVariables
830
831 Creates indices and counts all the model variables.
832
833 The model variables (states and outputs) are indexed
834 together, since they will be stored in one array.
835
836 The indices to the parmeters are into the global var map.
837
838 The indices to inputs are into the array of inputs.
839 */
IndexVariables(PVMMAPSTRCT pvmGlo)840 void IndexVariables (PVMMAPSTRCT pvmGlo)
841 {
842 int iIndex, iMax = MAX_VARS;
843
844 /* Get counts */
845 vnStates = ForAllVar(NULL, pvmGlo, &CountOneDecl, ID_STATE, NULL);
846 vnOutputs = ForAllVar(NULL, pvmGlo, &CountOneDecl, ID_OUTPUT, NULL);
847 vnInputs = ForAllVar(NULL, pvmGlo, &CountOneDecl, ID_INPUT, NULL);
848 vnParms = ForAllVar(NULL, pvmGlo, &CountOneDecl, ID_PARM, NULL);
849 vnModelVars = vnStates + vnOutputs;
850
851 /* Report all errors */
852 if (vnStates > MAX_VARS)
853 ReportError(NULL, RE_TOOMANYVARS, "state", (PSTR) &iMax);
854 if (vnOutputs > MAX_VARS)
855 ReportError(NULL, RE_TOOMANYVARS, "input", (PSTR) &iMax);
856 if (vnInputs > MAX_VARS)
857 ReportError(NULL, RE_TOOMANYVARS, "output", (PSTR) &iMax);
858 if (vnParms > (iMax = MAX_VARS - vnModelVars))
859 ReportError(NULL, RE_TOOMANYVARS, "parameter", (PSTR) &iMax);
860
861 if (vnStates > MAX_VARS
862 || vnInputs > MAX_VARS
863 || vnOutputs > MAX_VARS
864 || vnParms > iMax)
865 ReportError(NULL, RE_FATAL, NULL, NULL); /* Abort generation */
866
867 /* Set indices */
868 iIndex = 0;
869 ForAllVar(NULL, pvmGlo, &IndexOneVar, ID_STATE, (PVOID) &iIndex);
870 ForAllVar(NULL, pvmGlo, &IndexOneVar, ID_OUTPUT, (PVOID) &iIndex);
871
872 iIndex = 0;
873 ForAllVar(NULL, pvmGlo, &IndexOneVar, ID_INPUT, (PVOID) &iIndex);
874
875 iIndex = vnStates + vnOutputs + vnInputs;
876 ForAllVar(NULL, pvmGlo, &IndexOneVar, ID_PARM, (PVOID) &iIndex);
877
878 } /* IndexVariables */
879
880
881 /* ----------------------------------------------------------------------------
882 AdjustOneVar
883
884 Increments the dependent parameter handles of an input by the
885 iOffset given through the info pointer.
886
887 Callback function for ForAllVar().
888 */
AdjustOneVar(PFILE pfile,PVMMAPSTRCT pvm,PVOID pInfo)889 int AdjustOneVar (PFILE pfile, PVMMAPSTRCT pvm, PVOID pInfo)
890 {
891 PIFN pifn = (PIFN) pvm->szEqn;
892 WORD wOffset = *(PWORD) pInfo;
893
894 if (!pifn)
895 return 1; /* No eqn! No dependent parm! */
896
897 if (pifn->hMag)
898 pifn->hMag += wOffset;
899 if (pifn->hTper)
900 pifn->hTper += wOffset;
901 if (pifn->hT0)
902 pifn->hT0 += wOffset;
903 if (pifn->hTexp)
904 pifn->hTexp += wOffset;
905 if (pifn->hDecay)
906 pifn->hDecay += wOffset;
907
908 return 1;
909
910 } /* AdjustOneVar */
911
912
913 /* ----------------------------------------------------------------------------
914 AdjustVarHandles
915
916 Adjusts the variable handles on input definitions. They must be
917 incremented by the beginning offset of the parameter section of
918 the global variable map.
919 */
AdjustVarHandles(PVMMAPSTRCT pvmGlo)920 void AdjustVarHandles (PVMMAPSTRCT pvmGlo)
921 {
922 WORD wOffset = (WORD) vnInputs + vnStates + vnOutputs;
923
924 ForAllVar(NULL, pvmGlo, &AdjustOneVar, ID_INPUT, (PVOID) &wOffset);
925
926 } /* AdjustVarHandles */
927
928
929 /* ----------------------------------------------------------------------------
930 ReversePointers
931
932 Flips the pointer on the var-list so that the head is the tail
933 and the tail is the head and things aren't what they seem to be.
934
935 The dynamic equation list must be reversed so that the equations
936 appear in the right order since they were created as a stack.
937 */
ReversePointers(PVMMAPSTRCT * ppvm)938 void ReversePointers (PVMMAPSTRCT *ppvm)
939 {
940 PVMMAPSTRCT pvmPrev, pvmNext;
941
942 if (!ppvm || !(*ppvm)
943 || !(*ppvm)->pvmNextVar) /* List of one is already reversed! */
944 return;
945
946 pvmPrev = NULL;
947 while ((pvmNext = (*ppvm)->pvmNextVar)) {
948 (*ppvm)->pvmNextVar = pvmPrev;
949 pvmPrev = (*ppvm);
950 *ppvm = pvmNext;
951 } /* while */
952
953 (*ppvm)->pvmNextVar = pvmPrev; /* Link new head to reversed list */
954
955 } /* ReversePointers */
956
957
958 /* ----------------------------------------------------------------------------
959 AssertExistsEqn
960
961 Check that equation exists.
962
963 Uses info pointer of ForAllVar's callback as a second eqn list.
964
965 (1) If the second list is NULL, checks for initialization
966 of pvm.
967
968 (2) If the second list in non-NULL, checks it for definition.
969
970 Note that case (1) and (2) apply to inputs and dynamics eqns respectively.
971 Further assertion will have to be handled otherwise. Parms are
972 assured to be initialized when declared, and outputs can be initialized
973 to 0 automatically.
974
975 The errors are not reported as fatal so that all errors can
976 be discovered. They will cause exit subsequently.
977 */
AssertExistsEqn(PFILE pfile,PVMMAPSTRCT pvm,PVOID pInfo)978 int AssertExistsEqn (PFILE pfile, PVMMAPSTRCT pvm, PVOID pInfo)
979 {
980 int iReturn = 0;
981 PVMMAPSTRCT pvmDyn = (PVMMAPSTRCT) pInfo;
982
983 if (pvm->szEqn != vszHasInitializer) { /* Don't count these! */
984 if (pvmDyn) {
985 if (!(iReturn = (GetVarPTR(pvmDyn, pvm->szName) != NULL)))
986 ReportError(NULL, RE_NODYNEQN, pvm->szName, NULL);
987 }
988 else
989 if (!(iReturn = (pvm->szEqn != NULL)))
990 ReportError(NULL, RE_NOINPDEF, pvm->szName, NULL);
991 } /* if */
992
993 /* If HasInitializer, no error to report */
994 return(iReturn ? 1 : 0);
995
996 } /* AssertExistsEqn */
997
998
999 /* ----------------------------------------------------------------------------
1000 VerifyEqns
1001
1002 Check that each state variable has an associated differential equation
1003 defined.
1004
1005 Calls AssertExistsEqn on each and return a fatal error if one or more
1006 equation is missing.
1007 */
VerifyEqns(PVMMAPSTRCT pvmGlo,PVMMAPSTRCT pvmDyn)1008 void VerifyEqns (PVMMAPSTRCT pvmGlo, PVMMAPSTRCT pvmDyn)
1009 {
1010 BOOL bStatesOK;
1011
1012 bStatesOK = (vnStates == ForAllVar(NULL, pvmGlo, &AssertExistsEqn,
1013 ID_STATE, (PVOID) pvmDyn));
1014
1015 if (!bStatesOK)
1016 ReportError(NULL, RE_FATAL, NULL, "State equations missing.\n");
1017
1018 } /* VerifyEqns */
1019
1020
1021 /* ----------------------------------------------------------------------------
1022 AssertExistsOutputEqn
1023
1024 Check that an equation exists in either the Dynamics section or the
1025 CalcOutputs section for each output variable.
1026
1027 The errors are not reported as fatal so that all errors can
1028 be discovered. They will cause exit subsequently.
1029 */
AssertExistsOutputEqn(PFILE pfile,PVMMAPSTRCT pvm,PVOID pInfo)1030 int AssertExistsOutputEqn (PFILE pfile, PVMMAPSTRCT pvm, PVOID pInfo)
1031 {
1032 int iReturn = 0;
1033 PINPUTINFO pinfo = (PINPUTINFO) pInfo;
1034 PVMMAPSTRCT pvmDyn = (PVMMAPSTRCT) pinfo->pvmDynEqns;
1035 PVMMAPSTRCT pvmOut = (PVMMAPSTRCT) pinfo->pvmCalcOutEqns;
1036
1037 if (pvm->szEqn != vszHasInitializer) { /* Don't count these! */
1038
1039 if ((GetVarPTR(pvmDyn, pvm->szName) == NULL) &&
1040 (GetVarPTR(pvmOut, pvm->szName) == NULL)) {
1041 ReportError(NULL, RE_NOOUTPUTEQN, pvm->szName, NULL);
1042 iReturn = 0;
1043 }
1044 else
1045 iReturn = 1;
1046 } /* if */
1047
1048 /* If HasInitializer, no error to report */
1049 return(iReturn ? 1 : 0);
1050
1051 } /* AssertExistsOutputEqn */
1052
1053
1054 /* ----------------------------------------------------------------------------
1055 VerifyOutputEqns
1056
1057 Check that each output variable has an associated equation defined.
1058
1059 Calls AssertExistsOutputEqn on each and return a fatal error if one or more
1060 equation is missing.
1061 */
VerifyOutputEqns(PINPUTINFO pInfo)1062 void VerifyOutputEqns (PINPUTINFO pInfo)
1063 {
1064 BOOL bOutputsOK;
1065
1066 bOutputsOK = (vnOutputs == ForAllVar(NULL, pInfo->pvmGloVars,
1067 &AssertExistsOutputEqn, ID_OUTPUT,
1068 (PVOID) pInfo));
1069
1070 if (!bOutputsOK)
1071 ReportError(NULL, RE_FATAL, NULL, "Output equations missing.\n");
1072
1073 } /* VerifyOutputEqns */
1074
1075
1076 /* ----------------------------------------------------------------------------
1077 WriteModel
1078
1079 Writes the model calculation file szOutFilename for the parameters
1080 and dynamic equations given.
1081 */
WriteModel(PINPUTINFO pinfo,PSTR szFileOut)1082 void WriteModel (PINPUTINFO pinfo, PSTR szFileOut)
1083 {
1084 PFILE pfile;
1085
1086 if (!pinfo->pvmGloVars || (!pinfo->pvmDynEqns && !pinfo->pvmCalcOutEqns)) {
1087 printf("Error: No Dynamics, no outputs or no global variables defined\n");
1088 return;
1089 }
1090
1091 ReversePointers(&pinfo->pvmGloVars);
1092 ReversePointers(&pinfo->pvmDynEqns);
1093 ReversePointers(&pinfo->pvmScaleEqns);
1094 ReversePointers(&pinfo->pvmCalcOutEqns);
1095 ReversePointers(&pinfo->pvmJacobEqns);
1096 vpvmGloVarList = pinfo->pvmGloVars;
1097
1098 IndexVariables(pinfo->pvmGloVars);
1099 AdjustVarHandles(pinfo->pvmGloVars);
1100 VerifyEqns(pinfo->pvmGloVars, pinfo->pvmDynEqns);
1101
1102 VerifyOutputEqns(pinfo);
1103
1104 pfile = fopen(szFileOut, "w");
1105 if (pfile) {
1106
1107 /* Keep track of the model description file and generator name */
1108 vszModelFilename = pinfo->szInputFilename;
1109 vszModGenName = pinfo->szModGenName;
1110
1111 WriteHeader(pfile, szFileOut, pinfo->pvmGloVars);
1112
1113 WriteIncludes(pfile);
1114 WriteDecls (pfile, pinfo->pvmGloVars);
1115 WritebDelays (pfile, pinfo->bDelays);
1116 WriteVarMap (pfile, pinfo->pvmGloVars);
1117
1118 WriteInitModel (pfile, pinfo->pvmGloVars);
1119 WriteCalcDeriv (pfile, pinfo->pvmGloVars, pinfo->pvmDynEqns);
1120 WriteScale (pfile, pinfo->pvmGloVars, pinfo->pvmScaleEqns);
1121 WriteCalcJacob (pfile, pinfo->pvmGloVars, pinfo->pvmJacobEqns);
1122 WriteCalcOutputs(pfile, pinfo->pvmGloVars, pinfo->pvmCalcOutEqns);
1123
1124 fclose(pfile);
1125
1126 printf("\n* Created model file '%s'.\n\n", szFileOut);
1127
1128 } /* if */
1129 else
1130 ReportError(NULL, RE_CANNOTOPEN | RE_FATAL,
1131 szFileOut, "...in WriteModel ()");
1132
1133 } /* WriteModel */
1134
1135
1136 /* ----------------------------------------------------------------------------
1137 */
WriteOne_R_SODefine(PFILE pfile,PVMMAPSTRCT pvm,PVOID pInfo)1138 int WriteOne_R_SODefine (PFILE pfile, PVMMAPSTRCT pvm, PVOID pInfo)
1139 {
1140 static long iStates = 0; /* indexing states */
1141 static long iOutputs = 0; /* indexing outputs */
1142
1143 if (pvm->szEqn != vszHasInitializer) {
1144 fprintf(pfile, "#define ");
1145 WriteIndexName(pfile, pvm);
1146 if (TYPE(pvm) == ID_STATE) {
1147 fprintf(pfile, " 0x%05lx\n", iStates);
1148 iStates = iStates + 1;
1149 }
1150 else {
1151 fprintf(pfile, " 0x%05lx\n", iOutputs);
1152 iOutputs = iOutputs + 1;
1153 }
1154
1155 return 1;
1156 }
1157
1158 return 0;
1159
1160 } /* WriteOne_R_SODefine */
1161
1162
1163 /* ----------------------------------------------------------------------------
1164 Write_R_Scale
1165 */
Write_R_Scale(PFILE pfile,PVMMAPSTRCT pvmGlo,PVMMAPSTRCT pvmScale)1166 void Write_R_Scale (PFILE pfile, PVMMAPSTRCT pvmGlo, PVMMAPSTRCT pvmScale)
1167 {
1168 fprintf(pfile,
1169 "void getParms (double *inParms, double *out, int *nout) {\n");
1170 fprintf(pfile, "/*----- Model scaling */\n\n");
1171
1172 ForAllVar(pfile, pvmGlo, &WriteOneDecl, ID_LOCALSCALE, NULL);
1173
1174 fprintf(pfile, " int i;\n\n");
1175 fprintf(pfile, " for (i = 0; i < *nout; i++) {\n");
1176 fprintf(pfile, " parms[i] = inParms[i];\n }\n\n");
1177
1178 ForAllVar(pfile, pvmScale, &WriteOneEquation, ID_PARM, (PVOID) KM_SCALE);
1179
1180 fprintf(pfile, "\n for (i = 0; i < *nout; i++) {\n");
1181 fprintf(pfile, " out[i] = parms[i];\n }\n");
1182 fprintf(pfile, " }\n");
1183
1184 } /* Write_R_Scale */
1185
1186
1187 /* ----------------------------------------------------------------------------
1188 Write_R_State_Scale
1189 */
Write_R_State_Scale(PFILE pfile,PVMMAPSTRCT pvmScale)1190 void Write_R_State_Scale (PFILE pfile, PVMMAPSTRCT pvmScale)
1191 {
1192
1193 ForAllVar(pfile, pvmScale, &WriteOneEquation, ID_STATE, (PVOID) KM_SCALE);
1194
1195 ForAllVar(pfile, pvmScale, &WriteOneEquation, ID_INLINE, (PVOID) KM_SCALE);
1196
1197 } /* Write_R_State_Scale */
1198
1199
1200 /* ----------------------------------------------------------------------------
1201 Write_R_CalcDeriv
1202
1203 Writes the CalcDeriv() function for compatibility with R deSolve package.
1204 Writes dynamics equations in the order they appeared in the model definition
1205 file. Appends the CalcOutput equations, if any, at the end.
1206 */
Write_R_CalcDeriv(PFILE pfile,PVMMAPSTRCT pvmGlo,PVMMAPSTRCT pvmDyn,PVMMAPSTRCT pvmCalcOut)1207 void Write_R_CalcDeriv (PFILE pfile, PVMMAPSTRCT pvmGlo, PVMMAPSTRCT pvmDyn,
1208 PVMMAPSTRCT pvmCalcOut)
1209 {
1210 if (!pvmDyn)
1211 printf("No Dynamics{} equations.\n\n");
1212
1213 fprintf(pfile, "/*----- Dynamics section */\n\n");
1214 fprintf(pfile, "void derivs (int *neq, double *pdTime, double *y, ");
1215 fprintf(pfile, "double *ydot, double *yout, int *ip)\n{\n");
1216
1217 ForAllVar(pfile, pvmGlo, &WriteOneDecl, ID_LOCALDYN, NULL);
1218 ForAllVar(pfile, pvmGlo, &WriteOneDecl, ID_LOCALCALCOUT, NULL);
1219
1220 ForAllVar(pfile, pvmDyn, &WriteOneEquation, ALL_VARS, (PVOID) KM_DYNAMICS);
1221
1222 ForAllVar(pfile, pvmCalcOut, &WriteOneEquation, ALL_VARS,
1223 (PVOID) KM_CALCOUTPUTS);
1224
1225 fprintf(pfile, "\n} /* derivs */\n\n\n");
1226
1227 } /* Write_R_CalcDeriv */
1228
1229
1230 /* ----------------------------------------------------------------------------
1231 Write_R_InitModel
1232
1233 Writes the routine to initialize the model variables for R deSolve.
1234 */
Write_R_InitModel(PFILE pfile,PVMMAPSTRCT pvmGlo)1235 void Write_R_InitModel (PFILE pfile, PVMMAPSTRCT pvmGlo)
1236 {
1237 fprintf(pfile, "/*----- Initializers */\n");
1238 fprintf(pfile, "void initmod (void (* odeparms)(int *, double *))\n{\n");
1239 fprintf(pfile, " int N=%d;\n", vnParms);
1240 fprintf(pfile, " odeparms(&N, parms);\n");
1241 fprintf(pfile, "}\n\n");
1242
1243 fprintf(pfile, "void initforc (void (* odeforcs)(int *, double *))\n{\n");
1244 fprintf(pfile, " int N=%d;\n", vnInputs);
1245 fprintf(pfile, " odeforcs(&N, forc);\n");
1246 fprintf(pfile, "}\n\n\n");
1247
1248 if (bDelay) {
1249 fprintf(pfile, "/* Calling R code will ensure that input y has same\n");
1250 fprintf(pfile, " dimension as yini */\n");
1251 fprintf(pfile, "void initState (double *y)\n");
1252 fprintf(pfile, "{\n");
1253 fprintf(pfile, " int i;\n\n");
1254 fprintf(pfile, " for (i = 0; i < sizeof(yini) / sizeof(yini[0]); i++)\n");
1255 fprintf(pfile, " {\n");
1256 fprintf(pfile, " yini[i] = y[i];\n");
1257 fprintf(pfile, " }\n}\n\n");
1258 }
1259
1260 } /* Write_R_InitModel */
1261
1262
1263 /* ----------------------------------------------------------------------------
1264 Write_R_CalcJacob: empty in fact
1265 */
Write_R_CalcJacob(PFILE pfile,PVMMAPSTRCT pvmGlo,PVMMAPSTRCT pvmJacob)1266 void Write_R_CalcJacob (PFILE pfile, PVMMAPSTRCT pvmGlo, PVMMAPSTRCT pvmJacob)
1267 {
1268 fprintf(pfile, "/*----- Jacobian calculations: */\n");
1269 fprintf(pfile, "void jac (int *neq, double *t, double *y, int *ml, ");
1270 fprintf(pfile, "int *mu, ");
1271 fprintf(pfile, "double *pd, int *nrowpd, double *yout, int *ip)\n");
1272 fprintf(pfile, "{\n");
1273 ForAllVar(pfile, pvmGlo, &WriteOneDecl, ID_LOCALJACOB, NULL);
1274 ForAllVar(pfile, pvmJacob, &WriteOneEquation, ALL_VARS, (PVOID) KM_JACOB);
1275 fprintf(pfile, "\n} /* jac */\n\n\n");
1276
1277 } /* Write_R_CalcJacob */
1278
1279
1280 /* ----------------------------------------------------------------------------
1281 Write_R_Events: may contain Inlines
1282 */
Write_R_Events(PFILE pfile,PVMMAPSTRCT pvmGlo,PVMMAPSTRCT pvmEvents)1283 void Write_R_Events (PFILE pfile, PVMMAPSTRCT pvmGlo, PVMMAPSTRCT pvmEvents)
1284 {
1285 fprintf(pfile, "/*----- Events calculations: */\n");
1286 fprintf(pfile, "void event (int *n, double *t, double *y)\n");
1287 fprintf(pfile, "{\n");
1288 ForAllVar(pfile, pvmGlo, &WriteOneDecl, ID_LOCALEVENT, NULL);
1289 ForAllVar(pfile, pvmEvents, &WriteOneEquation, ALL_VARS, (PVOID) KM_EVENTS);
1290 fprintf(pfile, "\n} /* event */\n\n");
1291
1292 } /* Write_R_Events */
1293
1294
1295 /* ----------------------------------------------------------------------------
1296 Write_R_Roots: may contain Inlines
1297 */
Write_R_Roots(PFILE pfile,PVMMAPSTRCT pvmGlo,PVMMAPSTRCT pvmRoots)1298 void Write_R_Roots (PFILE pfile, PVMMAPSTRCT pvmGlo, PVMMAPSTRCT pvmRoots)
1299 {
1300 fprintf(pfile, "/*----- Roots calculations: */\n");
1301 fprintf(pfile, "void root (int *neq, double *t, double *y, ");
1302 fprintf(pfile, "int *ng, double *gout, double *out, int *ip)\n");
1303 fprintf(pfile, "{\n");
1304 ForAllVar(pfile, pvmGlo, &WriteOneDecl, ID_LOCALROOT, NULL);
1305 ForAllVar(pfile, pvmRoots, &WriteOneEquation, ALL_VARS, (PVOID) KM_ROOTS);
1306 fprintf(pfile, "\n} /* root */\n\n");
1307
1308 } /* Write_R_Roots */
1309
1310
1311 /* ----------------------------------------------------------------------------
1312 WriteOne_R_PIDefine
1313
1314 Write one global or local declaration for parameters and inputs (forcing
1315 functions) passed through R.
1316 Increments and prints two parallel counters ("iParms" and "iForcs").
1317 Callback for ForAllVar().
1318 */
WriteOne_R_PIDefine(PFILE pfile,PVMMAPSTRCT pvm,PVOID pInfo)1319 int WriteOne_R_PIDefine (PFILE pfile, PVMMAPSTRCT pvm, PVOID pInfo)
1320 {
1321 static long iParms = 0; /* indexing parameters */
1322 static long iForcs = 0; /* indexing input functions (forcing function) */
1323
1324 assert(TYPE(pvm) != ID_OUTPUT);
1325 assert(TYPE(pvm) != ID_STATE);
1326
1327 if (TYPE(pvm) == ID_PARM) {
1328 fprintf(pfile, "#define %s parms[%ld]\n", pvm->szName, iParms);
1329 iParms = iParms + 1;
1330 }
1331 else {
1332 fprintf(pfile, "#define %s forc[%ld]\n", pvm->szName, iForcs);
1333 iForcs = iForcs + 1;
1334 }
1335
1336 return(1);
1337
1338 } /* WriteOne_R_PIDefine */
1339
1340
1341 /* ----------------------------------------------------------------------------
1342 ForAllVarwSep
1343
1344 Takes a pfile, a pvm list head and a callback function which is called
1345 for all variables in the list if hType == ALL_VARS or only for
1346 only those that match hType otherwise. This version is for writing
1347 lists with separating punctuation. It sets End to -1 for the first
1348 item in the list, 0 for items in the middle of the list, and 1 for the
1349 last item.
1350 */
ForAllVarwSep(PFILE pfile,PVMMAPSTRCT pvm,PFI_CALLBACK pfiFunc,HANDLE hType,PVOID pinfo)1351 int ForAllVarwSep (PFILE pfile, PVMMAPSTRCT pvm, PFI_CALLBACK pfiFunc,
1352 HANDLE hType, PVOID pinfo)
1353 {
1354 int iTotal = 0;
1355 long End = -1;
1356 int iCount = 0;
1357
1358 while (pvm) {
1359 if (hType == ALL_VARS /* Do for all ... */
1360 || TYPE(pvm) == hType) { /* ... or do only for vars of hType */
1361
1362 if (pvm->szEqn != vszHasInitializer) {
1363 if (pfiFunc) {
1364 if (iCount > 0)
1365 End = 0;
1366 iTotal += (*pfiFunc) (pfile, pvm, (PVOID) End);
1367 iCount++;
1368 }
1369 else
1370 iTotal++; /* No func! just count */
1371 }
1372 }
1373 pvm = pvm->pvmNextVar;
1374 } /* while */
1375
1376 End = 1;
1377
1378 (*pfiFunc) (pfile, pvm, (PVOID) End);
1379
1380 return(iTotal);
1381
1382 } /* ForAllVarwSep */
1383
1384
1385 /* ---------------------------------------------------------------------------
1386 Is_numeric
1387
1388 Returns 1 if the argument is interpretable as a double precision number,
1389 0 otherwise.
1390 */
Is_numeric(PSTR str)1391 int Is_numeric(PSTR str)
1392 {
1393 double val;
1394 char *ptr;
1395
1396 if (str) {
1397 /* try to convert str to a number. ptr will hold a pointer
1398 to the part of the string that can't be converted.
1399 If str starts with a variable name, then str and ptr will
1400 be the same, but if str is something like 2 * Vblood,
1401 ptr will point to the '*'. In both these cases, strlen will
1402 return a non-zero value (I think). */
1403 val = strtod((char *)str, &ptr);
1404 if (strlen(ptr) > 0) return(0); else return(1);
1405 }
1406 else {
1407 return(2);
1408 }
1409
1410 } /* Is_numeric */
1411
1412
1413 /* ----------------------------------------------------------------------------
1414 WriteOne_R_ParmDecl
1415
1416 Write an R (comma-separated list) of parameter or state names. pinfo is
1417 used as a beginning - middle -end list flag.
1418 */
WriteOne_R_PSDecl(PFILE pfile,PVMMAPSTRCT pvm,PVOID pInfo)1419 int WriteOne_R_PSDecl (PFILE pfile, PVMMAPSTRCT pvm, PVOID pInfo)
1420 {
1421 PSTR szVarName;
1422 PSTR szZero = "0.0";
1423 long End = (long) pInfo;
1424 PSTR RHS;
1425 int iOut;
1426
1427 if (End < 1) {
1428
1429 szVarName = GetName(pvm, NULL, NULL, ID_NULL);
1430
1431 iOut = Is_numeric(pvm->szEqn);
1432 switch (iOut) {
1433 case 0:
1434 RHS = szZero;
1435 break;
1436 case 1:
1437 RHS = pvm->szEqn;
1438 break;
1439 case 2:
1440 RHS = szZero;
1441 }
1442 } /* end if */
1443
1444 switch (End) {
1445 case -1:
1446 fprintf(pfile, " %s = %s", szVarName, RHS);
1447 break;
1448 case 0:
1449 fprintf(pfile, ",\n %s = %s", szVarName, RHS);
1450 break;
1451 case 1:
1452 fprintf(pfile, "\n");
1453 return(0);
1454 }
1455
1456 return(1);
1457
1458 } /* WriteOne_R_PSDecl */
1459
1460
1461 /* ----------------------------------------------------------------------------
1462 WriteOne_R_ParmInit
1463 */
WriteOne_R_ParmInit(PFILE pfile,PVMMAPSTRCT pvm,PVOID pInfo)1464 int WriteOne_R_ParmInit (PFILE pfile, PVMMAPSTRCT pvm, PVOID pInfo)
1465 {
1466 PSTR szVarName;
1467 int iOut;
1468
1469 if (((long) pInfo) < 1) {
1470 szVarName = GetName(pvm, NULL, NULL, ID_NULL);
1471 iOut = Is_numeric(pvm->szEqn);
1472 if (iOut == 0) {
1473 fprintf(pfile, " %s = %s;\n",szVarName, pvm->szEqn);
1474 }
1475 return(1);
1476 }
1477 else
1478 return(0);
1479
1480 } /* WriteOne_R_ParmInit */
1481
1482
1483 /* ----------------------------------------------------------------------------
1484 WriteOne_R_StateInit
1485 */
WriteOne_R_StateInit(PFILE pfile,PVMMAPSTRCT pvm,PVOID pInfo)1486 int WriteOne_R_StateInit (PFILE pfile, PVMMAPSTRCT pvm, PVOID pInfo)
1487 {
1488 PSTR szVarName;
1489 long End = (long) pInfo;
1490
1491 if (End < 1)
1492 szVarName = GetName(pvm, NULL, NULL, ID_NULL);
1493
1494 switch (End) {
1495 case -1:
1496 fprintf(pfile, " %s = %s", szVarName,
1497 (pvm->szEqn ? pvm->szEqn : "0.0"));
1498 break;
1499
1500 case 0:
1501 fprintf(pfile, ",\n %s = %s", szVarName,
1502 (pvm->szEqn ? pvm->szEqn : "0.0"));
1503 break;
1504
1505 case 1:
1506 fprintf(pfile, "\n");
1507 return(0);
1508 }
1509
1510 return(1);
1511
1512 } /* WriteOne_R_StateInit */
1513
1514
1515 /* ----------------------------------------------------------------------------
1516 WriteOneOutputName
1517
1518 For R only
1519 */
WriteOneOutputName(PFILE pfile,PVMMAPSTRCT pvm,PVOID pInfo)1520 int WriteOneOutputName (PFILE pfile, PVMMAPSTRCT pvm, PVOID pInfo)
1521 {
1522 long END = (long) pInfo;
1523
1524 switch (END) {
1525 case -1:
1526 fprintf(pfile, " \"%s\"", pvm->szName);
1527 break;
1528
1529 case 0:
1530 fprintf(pfile, ",\n \"%s\"", pvm->szName);
1531 break;
1532
1533 case 1:
1534 fprintf(pfile, "\n");
1535 }
1536
1537 return(1);
1538
1539 } /* WriteOneOutputName */
1540
1541
1542 /* ----------------------------------------------------------------------------
1543 This is where we are going...
1544
1545 int WriteOneRForcing (PFILE pfile, PVMMAPSTRCT pvm, PVOID pInfo)
1546 {
1547 PIFN pifn = (PIFN) pvm->szEqn;
1548 }
1549
1550 void Write_R_Forcings(PFILE pfile, PVMMAPSTRCT pvmGlo)
1551 {
1552 ForAllVar(pfile, pvmGlo, &WriteOneRForcing, ID_INIT, NULL);
1553 }
1554 */
1555
1556
1557 /* ----------------------------------------------------------------------------
1558 Write_R_InitPOS
1559
1560 Write a utility R file defining functions to initialize parameters, outputs
1561 and states.
1562 */
Write_R_InitPOS(PFILE pfile,PVMMAPSTRCT pvmGlo,PVMMAPSTRCT pvmScale)1563 void Write_R_InitPOS (PFILE pfile, PVMMAPSTRCT pvmGlo, PVMMAPSTRCT pvmScale)
1564 {
1565 /* write R function initParms */
1566 fprintf(pfile, "initParms <- function(newParms = NULL) {\n");
1567 fprintf(pfile, " parms <- c(\n");
1568
1569 /* We write out all parameters here. If the rhs is a numerical constant,
1570 write it out. If it is an equation, then write out 0.0 */
1571 ForAllVarwSep(pfile, pvmGlo, &WriteOne_R_PSDecl, ID_PARM, NULL);
1572 fprintf(pfile, " )\n\n");
1573 fprintf(pfile, " if (!is.null(newParms)) {\n");
1574 fprintf(pfile, " if (!all(names(newParms) %%in%% c(names(parms)))) {\n");
1575 fprintf(pfile, " stop(\"illegal parameter name\")\n");
1576 fprintf(pfile, " }\n");
1577 fprintf(pfile, " parms[names(newParms)] <- newParms\n");
1578 fprintf(pfile, " }\n\n");
1579 fprintf(pfile, " parms <- within(as.list(parms), {\n");
1580
1581 /* Here, just write out the LHS and RHS for variables with eqns on the RHS */
1582 ForAllVarwSep(pfile, pvmGlo, &WriteOne_R_ParmInit, ID_PARM, NULL);
1583 fprintf(pfile, " })\n");
1584 fprintf(pfile, " out <- .C(\"getParms\", as.double(parms),\n");
1585 fprintf(pfile, " out=double(length(parms)),\n");
1586 fprintf(pfile, " as.integer(length(parms)))$out\n");
1587 fprintf(pfile, " names(out) <- names(parms)\n");
1588 fprintf(pfile, " out\n");
1589 fprintf(pfile, "}\n\n");
1590
1591 /* write R function Outputs */
1592 fprintf(pfile, "Outputs <- c(\n");
1593 ForAllVarwSep(pfile, pvmGlo, &WriteOneOutputName, ID_OUTPUT, NULL);
1594 fprintf(pfile, ")\n\n");
1595
1596 /* write R function initStates */
1597 bForInits = TRUE;
1598 fprintf(pfile, "initStates <- function(parms, newStates = NULL)"
1599 " {\n Y <- c(\n");
1600 ForAllVarwSep(pfile, pvmGlo, &WriteOne_R_PSDecl, ID_STATE, NULL);
1601 fprintf(pfile, " )\n\n");
1602
1603 /* do the next only if needed otherwise Y is reset to null */
1604 if (ForAllVar(pfile, pvmScale, NULL, ID_STATE, NULL) ||
1605 ForAllVar(pfile, pvmScale, NULL, ID_INLINE, NULL)) {
1606 fprintf(pfile, " Y <- within(c(as.list(parms),as.list(Y)), {");
1607 Write_R_State_Scale(pfile, pvmScale);
1608 fprintf(pfile, "\n })$Y\n\n");
1609 }
1610 fprintf(pfile, " if (!is.null(newStates)) {\n");
1611 fprintf(pfile, " if (!all(names(newStates) %%in%% c(names(Y)))) {\n");
1612 fprintf(pfile,
1613 " stop(\"illegal state variable name in newStates\")\n");
1614 fprintf(pfile, " }\n");
1615 fprintf(pfile, " Y[names(newStates)] <- newStates\n }\n\n");
1616
1617 if (bDelay) {
1618 fprintf(pfile, ".C(\"initState\", as.double(Y));\n");
1619 }
1620 fprintf(pfile, "Y\n}\n");
1621
1622 /* Write_R_Forcings: unfinished;
1623 translate mcsim's input information into a default forcing
1624 formulation for deSolve solvers. */
1625
1626 bForInits = FALSE;
1627
1628 } /* Write_R_InitPOS */
1629
1630
1631 /* ----------------------------------------------------------------------------
1632 Write_R_Decls
1633 */
Write_R_Decls(PFILE pfile,PVMMAPSTRCT pvmGlo)1634 void Write_R_Decls (PFILE pfile, PVMMAPSTRCT pvmGlo)
1635 {
1636 fprintf(pfile, "\n/* Model variables: States */\n");
1637 ForAllVar(pfile, pvmGlo, &WriteOne_R_SODefine, ID_STATE, NULL);
1638
1639 fprintf(pfile, "\n/* Model variables: Outputs */\n");
1640 ForAllVar(pfile, pvmGlo, &WriteOne_R_SODefine, ID_OUTPUT, NULL);
1641
1642 fprintf(pfile, "\n/* Parameters */\n");
1643 fprintf(pfile, "static double parms[%d];\n\n", vnParms);
1644 ForAllVar(pfile, pvmGlo, &WriteOne_R_PIDefine, ID_PARM, NULL);
1645
1646 fprintf(pfile, "\n/* Forcing (Input) functions */\n");
1647 fprintf(pfile, "static double forc[%d];\n\n", vnInputs);
1648 ForAllVar(pfile, pvmGlo, &WriteOne_R_PIDefine, ID_INPUT, NULL);
1649 fprintf(pfile, "\n");
1650
1651 if (bDelay) {
1652 fprintf(pfile, "/* Function definitions for delay differential "
1653 "equations */\n\n");
1654 fprintf(pfile, "int Nout=1;\n");
1655 fprintf(pfile, "int nr[1]={0};\n");
1656 fprintf(pfile, "double ytau[1] = {0.0};\n\n");
1657 fprintf(pfile, "static double yini[%d] = {",vnStates);
1658 int i;
1659 for (i = 1; i <= vnStates; i++) {
1660 if (i == vnStates) {
1661 fprintf(pfile, "0.0");
1662 } else{
1663 fprintf(pfile, "0.0, ");
1664 }
1665 }
1666 fprintf(pfile, "}; /*Array of initial state variables*/\n\n");
1667 fprintf(pfile, "void lagvalue(double T, int *nr, int N, double *ytau) "
1668 "{\n");
1669 fprintf(pfile, " static void(*fun)(double, int*, int, double*) = NULL;"
1670 "\n");
1671 fprintf(pfile, " if (fun == NULL)\n");
1672 fprintf(pfile, " fun = (void(*)(double, int*, int, double*))"
1673 "R_GetCCallable(\"deSolve\", \"lagvalue\");\n");
1674 fprintf(pfile, " return fun(T, nr, N, ytau);\n}\n\n");
1675 fprintf(pfile, "double CalcDelay(int hvar, double dTime, double delay) {"
1676 "\n");
1677 fprintf(pfile, " double T = dTime-delay;\n");
1678 fprintf(pfile, " if (dTime > delay){\n");
1679 fprintf(pfile, " nr[0] = hvar;\n");
1680 fprintf(pfile, " lagvalue( T, nr, Nout, ytau );\n}\n");
1681 fprintf(pfile, " else{\n");
1682 fprintf(pfile, " ytau[0] = yini[hvar];\n}\n");
1683 fprintf(pfile, " return(ytau[0]);\n}\n\n");
1684
1685 } /* end if */
1686
1687 } /* Write_R_Decls */
1688
1689
1690 /* ----------------------------------------------------------------------------
1691 */
Write_R_Includes(PFILE pfile)1692 void Write_R_Includes (PFILE pfile)
1693 {
1694 fprintf(pfile, "#include <R.h>\n");
1695
1696 if (bDelay) {
1697 fprintf(pfile, "#include <Rinternals.h>\n");
1698 fprintf(pfile, "#include <Rdefines.h>\n");
1699 fprintf(pfile, "#include <R_ext/Rdynload.h>\n");
1700 }
1701
1702 } /* Write_R_Includes */
1703
1704
1705 /* ----------------------------------------------------------------------------
1706 Write_R_Model
1707
1708 Writes a deSolve (R package) compatible C file "szOutFilename"
1709 corresponding to equations given.
1710 */
Write_R_Model(PINPUTINFO pinfo,PSTR szFileOut)1711 void Write_R_Model (PINPUTINFO pinfo, PSTR szFileOut)
1712 {
1713 static PSTRLEX vszModified_Title;
1714 PFILE pfile;
1715 PSTR Rfile;
1716 PSTR Rappend = "_inits.R";
1717 size_t nRout, nbase;
1718 char * lastdot;
1719
1720 /* set global flag ! */
1721 bForR = TRUE;
1722
1723 if (!pinfo->pvmGloVars || (!pinfo->pvmDynEqns && !pinfo->pvmCalcOutEqns)) {
1724 printf ("Error: No Dynamics, outputs or global variables defined\n");
1725 return;
1726 }
1727
1728 ReversePointers(&pinfo->pvmGloVars);
1729 ReversePointers(&pinfo->pvmDynEqns);
1730 ReversePointers(&pinfo->pvmScaleEqns);
1731 ReversePointers(&pinfo->pvmCalcOutEqns);
1732 ReversePointers(&pinfo->pvmJacobEqns);
1733 ReversePointers(&pinfo->pvmEventEqns);
1734 ReversePointers(&pinfo->pvmRootEqns);
1735 vpvmGloVarList = pinfo->pvmGloVars;
1736
1737 IndexVariables(pinfo->pvmGloVars);
1738 AdjustVarHandles(pinfo->pvmGloVars);
1739 VerifyEqns(pinfo->pvmGloVars, pinfo->pvmDynEqns);
1740
1741 VerifyOutputEqns(pinfo);
1742
1743 pfile = fopen(szFileOut, "w");
1744 if (pfile) {
1745
1746 /* Keep track of the model description file and generator name */
1747 vszModelFilename = pinfo->szInputFilename;
1748 vszModGenName = pinfo->szModGenName;
1749
1750 sprintf(vszModified_Title, "%s %s", szFileOut, "for R deSolve package");
1751 WriteHeader(pfile, vszModified_Title, pinfo->pvmGloVars);
1752
1753 Write_R_Includes(pfile);
1754 Write_R_Decls(pfile, pinfo->pvmGloVars);
1755
1756 Write_R_InitModel(pfile, pinfo->pvmGloVars);
1757 Write_R_Scale (pfile, pinfo->pvmGloVars, pinfo->pvmScaleEqns);
1758 Write_R_CalcDeriv(pfile, pinfo->pvmGloVars, pinfo->pvmDynEqns,
1759 pinfo->pvmCalcOutEqns); /* fold in CaclOutput */
1760 Write_R_CalcJacob(pfile, pinfo->pvmGloVars, pinfo->pvmJacobEqns);
1761 Write_R_Events (pfile, pinfo->pvmGloVars, pinfo->pvmEventEqns);
1762 Write_R_Roots (pfile, pinfo->pvmGloVars, pinfo->pvmRootEqns);
1763
1764 fclose(pfile);
1765
1766 printf("\n* Created C model file '%s'.\n\n", szFileOut);
1767
1768 } /* if */
1769 else
1770 ReportError(NULL, RE_CANNOTOPEN | RE_FATAL,
1771 szFileOut, "in Write_R_Model ()");
1772
1773 /* Write a function to initialize everything */
1774 /* Construct the name of the R file. If szFileOut is
1775 fu.c, R file is fu_inits.R
1776 */
1777 lastdot = strrchr(szFileOut,'.');
1778 if (lastdot != NULL)
1779 *lastdot = '\0';
1780 nbase = strlen(szFileOut);
1781
1782 /* Length of buffer for new file name: includes terminating null */
1783 nRout = nbase + strlen(Rappend) + 1;
1784 Rfile = (PSTR) malloc(nRout);
1785 Rfile = strncpy(Rfile, szFileOut, nbase);
1786 Rfile[nbase] = '\0';
1787 Rfile = strcat(Rfile, Rappend);
1788 pfile = fopen(Rfile, "w");
1789 if (pfile) {
1790 Write_R_InitPOS(pfile, pinfo->pvmGloVars, pinfo->pvmScaleEqns);
1791 fclose(pfile);
1792 printf("\n* Created R parameter initialization file '%s'.\n\n",Rfile);
1793 }
1794 else
1795 ReportError(NULL, RE_CANNOTOPEN | RE_FATAL,
1796 Rfile, "in Write_R_Model ()");
1797
1798 free(Rfile);
1799
1800 } /* Write_R_Model */
1801
1802
1803