1 /* yourcode.c
2
3 Originally written by Frederic Bois
4
5 Copyright (c) 1993-2017 Free Software Foundation, Inc.
6
7 This file is part of GNU MCSim.
8
9 GNU MCSim is free software; you can redistribute it and/or
10 modify it under the terms of the GNU General Public License
11 as published by the Free Software Foundation; either version 3
12 of the License, or (at your option) any later version.
13
14 GNU MCSim is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18
19 You should have received a copy of the GNU General Public License
20 along with GNU MCSim; if not, see <http://www.gnu.org/licenses/>
21
22 Contains the routines most susceptible to be modified by the user and
23 some utilities.
24 */
25
26 #include <stdio.h>
27 #include <stdlib.h>
28 #include <string.h>
29 #include <math.h>
30
31 #include "lexerr.h"
32 #include "simmonte.h"
33 #include "yourcode.h"
34
35
36 /* ---------------------------------------------------------------------------
37 DoStep_by_Step
38
39 routine called after each successful step of the integrator.
40 It can be used for interupts management, for step by step
41 printing or operations such as finding a maximum etc.
42 */
DoStep_by_Step()43 void DoStep_by_Step (/* your needed parameters here, e.g.:
44 double t, long *neq, double *y */)
45 {
46 /* example of code:
47 static FILE *fout;
48 int i;
49
50 if (!fout) fout = fopen("step.out", "w");
51
52 fprintf(fout, "%g\t", t);
53 for (i = 1; i <= *neq; i++) fprintf(fout, "%g\t", y[i]);
54 fprintf(fout, "\n");
55
56 */
57
58 } /* DoStep_by_Step */
59
60
61 /* ---------------------------------------------------------------------------
62 TransformPred
63
64 At least flattens the model predictions in a simple array after a
65 Monte Carlo or a SetPoint simulation.
66
67 Changing this routine should be avoided and using output variables
68 defined through the model specification file should be preferred.
69
70 If you change it make sure that you allocate the data array of the
71 pMCPredOut structure and specify its length. See the routine
72 OutspecToLinearArray for exemple.
73
74 At most it allows the user to manipulate the data output
75 for creating summaries (e.g. sums of variables) which
76 better relate to the experimental data simulated. Those summaries
77 are placed in the pMCPredOut structure and will be used and printed.
78 */
79
TransformPred(PANALYSIS panal,PMCPREDOUT pMCPredOut)80 void TransformPred (PANALYSIS panal, PMCPREDOUT pMCPredOut)
81 {
82
83 OutspecToLinearArray (panal, pMCPredOut); /* no tranformation */
84
85 } /* TransformPred */
86
87
88 /* ---------------------------------------------------------------------------
89 OutspecToLinearArray
90
91 Flattens the panal nested output arrays on a single array.
92 Allocate the data array of the pMCPredOut structure and sets the
93 dimension pMCPredOut->nbrdy to the length of the data array.
94 */
95
OutspecToLinearArray(PANALYSIS panal,PMCPREDOUT pMCPredOut)96 void OutspecToLinearArray (PANALYSIS panal, PMCPREDOUT pMCPredOut)
97 {
98 POUTSPEC pos;
99 long i, j, k;
100
101 pMCPredOut->nbrdy = 0;
102
103 /* get the size needed for the data array of pMCPredOut
104 there should be one cell for each experiment, variable,
105 and output time
106 */
107 for (i = 0; i < panal->expGlobal.iExp; i++)
108 for (j = 0, pos = &panal->rgpExps[i]->os; j < pos->nOutputs; j++)
109 for (k = 0; k < pos->pcOutputTimes[j]; k++)
110 pMCPredOut->nbrdy++;
111
112 /* allocate data */
113
114 if (!(pMCPredOut->pred))
115 if ( !(pMCPredOut->pred = InitdVector (pMCPredOut->nbrdy)))
116 ReportError (NULL, RE_OUTOFMEM | RE_FATAL, "OutspecToLinearArray", NULL);
117
118 pMCPredOut->nbrdy = 0;
119 /* fill in pred array */
120 for (i = 0; i < panal->expGlobal.iExp; i++)
121 for (j = 0, pos = &panal->rgpExps[i]->os; j < pos->nOutputs; j++)
122 for (k = 0; k < pos->pcOutputTimes[j]; k++)
123 pMCPredOut->pred[pMCPredOut->nbrdy++] = pos->prgdOutputVals[j][k];
124
125 } /* OutspecToLinearArray */
126
127
128 /* ---------------------------------------------------------------------------
129 Definite_Integral
130
131 Performs the integration of an analytic function using Rhomberg algorithm.
132
133 Adapted from the algorithm described in the book Numerical Recipes by
134 Press et al.
135
136 */
137
Definite_Integral(double (* Function)(double),double dFrom,double dTo)138 double Definite_Integral (double (*Function)(double), double dFrom, double dTo)
139 {
140 #define PRECISION 1.0e-6
141 #define MAX_STEPS 25
142 #define MIN_STEP 4
143
144 double dTotal_Area, dDelta_Area;
145 double pdArea[MAX_STEPS+2], pdTemp[MAX_STEPS+2];
146 int i;
147
148 if (dFrom >= dTo) {
149 if (dFrom == dTo)
150 return (0);
151 else {
152 printf ("\nError: inverted integration bounds in Definite_Integral - "
153 "Exiting\n\n");
154 exit (0);
155 }
156 }
157
158 pdTemp[0] = 1;
159
160 for (i = 0; i < MAX_STEPS; i++) {
161
162 pdArea[i] = Trapezes (Function, dFrom, dTo, i+1);
163
164 if (i >= MIN_STEP) {
165 Interpolate_Poly (&pdTemp[i-MIN_STEP], &pdArea[i-MIN_STEP], MIN_STEP + 1,
166 0.0, &dTotal_Area, &dDelta_Area);
167
168 if (dTotal_Area == 0) {
169 if (fabs(dDelta_Area) < PRECISION)
170 return dTotal_Area;
171 }
172 else {
173 if (fabs(dDelta_Area) < PRECISION * fabs(dTotal_Area))
174 return dTotal_Area;
175 }
176 }
177
178 pdArea[i+1] = pdArea[i];
179 pdTemp[i+1] = 0.25 * pdTemp[i];
180
181 }
182
183 printf ("\nError: Too many steps in routine Definite_Integral "
184 "- Exiting\n\n");
185 exit (0);
186
187 #undef PRECISION
188 #undef MAX_STEPS
189 #undef MIN_STEP
190
191 } /* Definite_Integral */
192
193
194 /* ---------------------------------------------------------------------------
195 Interpolate_Poly
196
197 Adapted from the algorithm described in the book Numerical Recipes by
198 Press et al.
199
200 */
201
Interpolate_Poly(double rgdX[],double rgdY[],int n,double x,double * pdY,double * pdDY)202 void Interpolate_Poly (double rgdX[], double rgdY[], int n, double x,
203 double *pdY, double *pdDY)
204 {
205 int i, j, nIndex = 1;
206 double dDenom, dDiff, dTemp1, dTemp2;
207 static PDOUBLE pdTerm1 = NULL, pdTerm2 = NULL;
208
209 if (!pdTerm1)
210 if ( !(pdTerm1 = InitdVector (n+1)) || !(pdTerm2 = InitdVector (n+1)))
211 ReportError (NULL, RE_OUTOFMEM | RE_FATAL, "Interpolate_Poly", NULL);
212
213 dDiff = fabs (x - rgdX[0]);
214 pdTerm1[0] = rgdY[0];
215 pdTerm2[0] = rgdY[0];
216
217 for (i = 1; i < n; i++) {
218 /* find the smallest difference between x and components of x[] */
219 if ((dTemp1 = fabs (x - rgdX[i])) < dDiff) {
220 nIndex = i;
221 dDiff = dTemp1;
222 }
223 pdTerm1[i] = rgdY[i];
224 pdTerm2[i] = rgdY[i];
225 }
226
227 *pdY = rgdY[nIndex--];
228
229 for (j = 1; j < n; j++) {
230 for (i = 0; i < n - j; i++) {
231 dTemp1 = rgdX[i] - x;
232 dTemp2 = rgdX[i+j] - x;
233 if ((dDenom = dTemp1 - dTemp2) == 0) {
234 printf ("\nError: null denominator in Interpolate_Poly - Exiting\n\n");
235 exit (0);
236 }
237 dDenom = (pdTerm1[i+1] - pdTerm2[i]) / dDenom;
238 pdTerm2[i] = dTemp2 * dDenom;
239 pdTerm1[i] = dTemp1 * dDenom;
240 }
241
242 *pdDY = (2 * (nIndex + 1) < (n - j) ?
243 pdTerm1[nIndex+1] : pdTerm2[nIndex--]);
244
245 *pdY = *pdY + *pdDY;
246 }
247
248 /* free(pdTerm1); free(pdTerm2); */
249
250 } /* Interpolate_Poly */
251
252
253 /* ---------------------------------------------------------------------------
254 Trapezes
255
256 Integrates an analytic function using the trapeze method. This is called
257 recursively.
258 The function to integrate is evaluated 2^(nStep - 2) times.
259
260 Adapted from the algorithm described in the book Numerical Recipes by
261 Press et al.
262
263 */
264
Trapezes(double (* Function)(double x),double dFrom,double dTo,int nSteps)265 double Trapezes (double (*Function)(double x), double dFrom, double dTo,
266 int nSteps)
267 {
268
269 double x, dSum, dDelta;
270 static double dStoredArea;
271 int i, j;
272
273 if (nSteps == 1) { /* a single trapeze */
274 return (dStoredArea = 0.5 * (dTo - dFrom) *
275 ((*Function)(dFrom) + (*Function)(dTo)));
276 }
277 else {
278 i = 1;
279 for (j = 0; j < nSteps - 2; j++)
280 i = i << 1;
281
282 dDelta = (dTo - dFrom) / (double) i;
283
284 x = dFrom + 0.5 * dDelta;
285 dSum = 0;
286 while (x < dTo) {
287 dSum = dSum + (*Function)(x);
288 x = x + dDelta;
289 }
290
291 return (dStoredArea = 0.5 * (dStoredArea + dDelta * dSum));
292
293 }
294
295 } /* Trapezes */
296
297
298 /* End */
299