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