1 /* simo.c
2
3 Copyright (c) 1991-2017 Free Software Foundation, Inc.
4
5 This file is part of GNU MCSim.
6
7 GNU MCSim is free software; you can redistribute it and/or
8 modify it under the terms of the GNU General Public License
9 as published by the Free Software Foundation; either version 3
10 of the License, or (at your option) any later version.
11
12 GNU MCSim is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with GNU MCSim; if not, see <http://www.gnu.org/licenses/>
19
20 Output routines for the simulation
21 */
22
23 #include <math.h>
24 #include <stdio.h>
25 #include <stdlib.h>
26 #include <string.h>
27
28 #include "lexerr.h"
29 #include "simo.h"
30 #include "modelu.h"
31
32 static char vszDefOutFilename[] = "sim.out";
33 static char vszDefMCOutFilename[] = "simmc.out";
34
35
36 /* ----------------------------------------------------------------------------
37 SaveOutputs
38
39 Also saves states
40 */
41
SaveOutputs(PEXPERIMENT pexp,PDOUBLE pdTout)42 void SaveOutputs (PEXPERIMENT pexp, PDOUBLE pdTout)
43 {
44 #define SO_EPSILON (1e-100) /* Smaller values are zeroed */
45
46 static PDOUBLE rgdInterpStates, rgdInterpDeriv;
47 int i, j, index;
48 PMODELINFO pmod = pexp->pmodelinfo;
49 POUTSPEC pos = &pexp->os;
50 extern IFN vrgInputs[]; /* Input Function records */
51
52 if (!(rgdInterpStates) || !(rgdInterpDeriv))
53 if ( !(rgdInterpStates = InitdVector (GetNModelVars ())) ||
54 !(rgdInterpDeriv = InitdVector (GetNModelVars ())))
55 ReportError (NULL, RE_OUTOFMEM | RE_FATAL, "SaveOutputs", NULL);
56
57 memcpy (rgdInterpStates, pmod->pdModelVars,
58 pmod->nModelVars*sizeof(double));
59
60 /* Update inputs and outputs defined only in Dynamics */
61 CalcDeriv(rgdInterpStates, rgdInterpDeriv, pdTout);
62
63 /* Update output scaling */
64 CalcOutputs (rgdInterpStates, rgdInterpDeriv, pdTout);
65
66 for (i = 0; i < pos->nOutputs; i++) {
67
68 /* Save interpolated value if there are still times to output
69 for this variable, and if this time is scheduled */
70
71 if (pos->piCurrentOut[i] < pos->pcOutputTimes[i]
72 && *pdTout == pos->prgdOutputTimes[i][pos->piCurrentOut[i]]) {
73 double dTmp = 0;
74
75 if (IsModelVar(pos->phvar_out[i])) /* Use interp'd model value */
76 dTmp = rgdInterpStates[ ModelIndex(pos->phvar_out[i])];
77
78 else { /* Use current parm/input value */
79 index = HINDEX(pos->phvar_out[i]);
80 if (IsInput(pos->phvar_out[i]) &&
81 (vrgInputs[index].iType == IFN_SPIKES)) {
82
83 j = vrgInputs[index].iDoseCur;
84
85 if ((vrgInputs[index].rgT0s[j] == pexp->dTime) &&
86 (j < vrgInputs[index].nDoses))
87 dTmp = vrgInputs[index].rgMags[j];
88 else
89 dTmp = 0;
90
91 }
92 else
93 dTmp = GetVarValue (pos->phvar_out[i]);
94 }
95
96 if (fabs(dTmp) < SO_EPSILON) /* Avoid silly little numbers */
97 dTmp = 0.0;
98
99 pos->prgdOutputVals[i][pos->piCurrentOut[i]++] = dTmp;
100
101 } /* if */
102 } /* for */
103
104 } /* SaveOutputs */
105
106
107 /* -----------------------------------------------------------------------------
108 NextOutputTime
109
110 Returns in pdTout,the next time, pdTout, at which an variable is
111 to be output.
112 */
113
NextOutputTime(PEXPERIMENT pexp,PDOUBLE pdTout,PINT piOut)114 void NextOutputTime (PEXPERIMENT pexp, PDOUBLE pdTout, PINT piOut)
115 {
116 if (pexp->dTime < pexp->dTfinal) {
117 if (++*piOut < pexp->os.cDistinctTimes) {
118 *pdTout = pexp->os.rgdDistinctTimes[*piOut];
119 }
120 else {
121 *pdTout = pexp->dTfinal;
122 }
123 }
124
125 } /* NextOutputTime */
126
127
128 /* -----------------------------------------------------------------------------
129 WriteOneMod
130
131 writes one parameter modification from the list. Inputs are *not*
132 written.
133 */
134
WriteOneMod(PVOID pData,PVOID pInfo)135 int WriteOneMod (PVOID pData, PVOID pInfo)
136 {
137 PMCVAR pmcvar = (PMCVAR) pData;
138 PFILE pfile = (PFILE) pInfo;
139
140 if (!IsInput (pmcvar->hvar))
141 fprintf(pfile, "%g\t", pmcvar->dVal);
142
143 return 0;
144
145 } /* WriteOneMod */
146
147
148 /* -----------------------------------------------------------------------------
149 WriteMCHeader
150
151 Write a tabulated text header with the run number, the list of parameters
152 and outputs.
153 */
154
WriteMCHeader(PFILE pfileOut,PANALYSIS panal)155 void WriteMCHeader (PFILE pfileOut, PANALYSIS panal)
156 {
157 long i, j, k;
158 PMONTECARLO pmc = &panal->mc;
159 OUTSPEC *pos;
160
161 fprintf (pfileOut, "Iter");
162
163 for (i = 0; i < pmc->nParms; i++)
164 fprintf (pfileOut, "\t%s", GetVarName(pmc->rgpMCVar[i]->hvar));
165
166 /* print the outputs as they come with experiment and time code */
167 for (i = 0; i < panal->expGlobal.iExp; i++) {
168 pos = &panal->rgpExps[i]->os;
169 for (j = 0; j < pos->nOutputs; j++) {
170 for (k = 0; k < pos->pcOutputTimes[j]; k++)
171 fprintf (pfileOut, "\t%s_%ld.%ld", pos->pszOutputNames[j], i+1, k+1);
172 } /* for j */
173 } /* for i */
174
175 fprintf (pfileOut, "\n");
176
177 fflush (pfileOut);
178
179 } /* WriteMCHeader */
180
181
182 /* -----------------------------------------------------------------------------
183 OpenMCFiles
184
185 For each processor (if more than one), open the Monte Carlo output
186 file to be written to by WriteMCOutput()
187
188 Report fatal error (will lead to exit) if the file cannot be opened.
189 */
OpenMCFiles(PANALYSIS panal)190 void OpenMCFiles (PANALYSIS panal)
191 {
192 PMONTECARLO pmc = &panal->mc;
193
194 /* Use command line spec if given */
195 if (panal->bCommandLineSpec) {
196 free (pmc->szMCOutfilename);
197 panal->bAllocatedFileName = FALSE;
198 pmc->szMCOutfilename = panal->szOutfilename;
199 }
200 else
201 if (!(pmc->szMCOutfilename)) /* Default if none given */
202 pmc->szMCOutfilename = vszDefMCOutFilename;
203
204 /* prefix the filename with the rank of the process if more than one
205 process is used */
206 if (panal->size > 1) {
207 char* with_rank = malloc(sizeof(char)*(6+strlen(pmc->szMCOutfilename)));
208 sprintf(with_rank, "%04d_%s", panal->rank, pmc->szMCOutfilename);
209 pmc->szMCOutfilename = with_rank;
210 }
211
212 if (!pmc->pfileMCOut
213 && !(pmc->pfileMCOut = fopen (pmc->szMCOutfilename, "w"))) {
214 ReportError (NULL, RE_FATAL | RE_CANNOTOPEN, pmc->szMCOutfilename,
215 "OpenMCFiles()");
216 }
217
218 WriteMCHeader (pmc->pfileMCOut, panal);
219
220 } /* OpenMCFiles */
221
222
223 /* -----------------------------------------------------------------------------
224 CloseMCFiles
225
226 Closes output files associated with Monte Carlo and set points runs
227 */
228
CloseMCFiles(PANALYSIS panal)229 void CloseMCFiles (PANALYSIS panal)
230 {
231 fclose (panal->mc.pfileMCOut);
232 printf ("\nWrote results to \"%s\"\n", panal->mc.szMCOutfilename);
233
234 } /* CloseMCFiles */
235
236
237 /* -----------------------------------------------------------------------------
238 WriteMCOutput
239
240 Output the parameters for this run and the results of the
241 simulation (passed through TransformPred).
242 */
243
WriteMCOutput(PANALYSIS panal,PMCPREDOUT pmcpredout)244 void WriteMCOutput (PANALYSIS panal, PMCPREDOUT pmcpredout)
245 {
246 PFILE pfileMC;
247 PMONTECARLO pmc = &panal->mc;
248
249 pfileMC = pmc->pfileMCOut;
250
251 fprintf (pfileMC, "%ld\t", panal->mc.lRun);
252
253 /* Include parameter values for that run */
254 WriteArray (pfileMC, panal->mc.nParms, panal->mc.rgdParms);
255 fprintf (pfileMC, "\t");
256
257 /* write the flattened and eventually transformed predictions */
258 WriteArray (pfileMC, pmcpredout->nbrdy, pmcpredout->pred);
259 fprintf (pfileMC, "\n");
260
261 fflush (pfileMC);
262
263 } /* WriteMCOutput */
264
265
266 /* -----------------------------------------------------------------------------
267 WriteNormalOutput
268
269 Write the results in the output file. This procedure is
270 called only from time to time in order to save storage space
271 */
WriteNormalOutput(PANALYSIS panal,PEXPERIMENT pexp)272 void WriteNormalOutput (PANALYSIS panal, PEXPERIMENT pexp)
273 {
274 long i, j;
275 PFILE pfile;
276 POUTSPEC pos;
277
278 if (!panal) return;
279
280 pos = &pexp->os;
281
282 if (!panal->szOutfilename)
283 panal->szOutfilename = vszDefOutFilename;
284
285 /* prefix the filename with the rank of the process if more than one
286 process is used */
287 if (panal->size > 1) {
288 char* with_rank = malloc(sizeof(char)*(5+strlen(panal->szOutfilename)));
289 sprintf(with_rank,"%4d%s",panal->rank,panal->szOutfilename);
290 panal->szOutfilename = with_rank;
291 }
292
293 if (!(panal->pfileOut))
294 if (!(panal->pfileOut = fopen (panal->szOutfilename, "w")))
295 ReportError (NULL, RE_CANNOTOPEN | RE_FATAL, panal->szOutfilename, NULL);
296
297 pfile = panal->pfileOut;
298 fprintf (pfile, "Results of Simulation %d\n\n", pexp->iExp);
299
300 /* Vertical output: Formatted Time1 Out_Var1 Out_Var2 ... */
301
302 fprintf (pfile, "Time");
303
304 for (i = 0; i < pos->nOutputs; i++)
305 fprintf (pfile, "\t%s", pos->pszOutputNames[i]);
306 fprintf (pfile, "\n");
307
308 for (j = 0; j < pos->nOutputs; j++)
309 pos->piCurrentOut[j] = 0;
310
311 for (i = 0; i < pos->cDistinctTimes; i++) {
312 fprintf (pfile, "%g", pos->rgdDistinctTimes[i]);
313 for (j = 0; j < pos->nOutputs; j++) {
314
315 if (pos->piCurrentOut[j] < pos->pcOutputTimes[j]
316 && pos->rgdDistinctTimes[i]
317 == pos->prgdOutputTimes[j][pos->piCurrentOut[j]])
318 fprintf (pfile, "\t%g",
319 pos->prgdOutputVals[j][pos->piCurrentOut[j]++]);
320
321 else
322 fprintf (pfile, "\t");
323
324 } /* for */
325
326 fprintf (pfile, "\n");
327
328 } /* for */
329
330 fprintf (pfile, "\n\n");
331
332 } /* WriteNormalOutput */
333