Bug Summary

File:OMCompiler/Compiler/runtime/./SimulationResultsCmpTubes.c
Warning:line 753, column 3
Value stored to 'maxn' is never read

Annotated Source Code

[?] Use j/k keys for keyboard navigation

1/*
2 * This file is part of OpenModelica.
3 *
4 * Copyright (c) 1998-CurrentYear, Open Source Modelica Consortium (OSMC),
5 * c/o Linköpings universitet, Department of Computer and Information Science,
6 * SE-58183 Linköping, Sweden.
7 *
8 * All rights reserved.
9 *
10 * THIS PROGRAM IS PROVIDED UNDER THE TERMS OF GPL VERSION 3 LICENSE OR
11 * THIS OSMC PUBLIC LICENSE (OSMC-PL) VERSION 1.2.
12 * ANY USE, REPRODUCTION OR DISTRIBUTION OF THIS PROGRAM CONSTITUTES
13 * RECIPIENT'S ACCEPTANCE OF THE OSMC PUBLIC LICENSE OR THE GPL VERSION 3,
14 * ACCORDING TO RECIPIENTS CHOICE.
15 *
16 * The OpenModelica software and the Open Source Modelica
17 * Consortium (OSMC) Public License (OSMC-PL) are obtained
18 * from OSMC, either from the above address,
19 * from the URLs: http://www.ida.liu.se/projects/OpenModelica or
20 * http://www.openmodelica.org, and in the OpenModelica distribution.
21 * GNU version 3 is obtained from: http://www.gnu.org/copyleft/gpl.html.
22 *
23 * This program is distributed WITHOUT ANY WARRANTY; without
24 * even the implied warranty of MERCHANTABILITY or FITNESS
25 * FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY SET FORTH
26 * IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS OF OSMC-PL.
27 *
28 * See the full OSMC Public License conditions for more details.
29 *
30 */
31
32/*
33 * The original source:
34 * https://svn.modelica.org/projects/Modelica/branches/tools/csv-compare/Modelica_ResultCompare/Tubes.cs
35 * Comes with the following terms [BSD 3-clause]:
36 *
37 * Copyright (c) 2013, ITI GmbH Dresden
38 * All rights reserved.
39 *
40 * Redistribution and use in source and binary forms, with or without
41 * modification, are permitted provided that the following conditions are met:
42 *
43 * Redistributions of source code must retain the above copyright notice,
44 * this list of conditions and the following disclaimer.
45 * Redistributions in binary form must reproduce the above copyright notice,
46 * this list of conditions and the following disclaimer in the documentation
47 * and/or other materials provided with the distribution.
48 *
49 * Neither the name of the ITI GmbH nor the names of its contributors may be
50 * used to endorse or promote products derived from this software without
51 * specific prior written permission.
52 *
53 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
54 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
55 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
56 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
57 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
58 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
59 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
60 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
61 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
62 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
63 * POSSIBILITY OF SUCH DAMAGE.
64 */
65
66#include <stdlib.h>
67#include <math.h>
68#include <gc.h>
69#include <stdio.h>
70
71typedef struct {
72 double *mh,*ml,*xHigh,*xLow,*yHigh,*yLow;
73 int *i0h,*i1h,*i0l,*i1l;
74 double tStart,tStop,x1,y1,x2,y2,currentSlope,slopeDif,delta,S,xRelEps,xMinStep,min,max;
75 size_t countLow,countHigh,length;
76} privates;
77
78static inline__inline__ int intmax(int a, int b) {
79 return a>b ? a : b;
80}
81
82static void generateHighTube(privates *priv, double *x, double *y)
83{
84 int index = priv->countHigh - 1;
85 double m1 = priv->mh[index];
86 double m2 = priv->mh[index - 1];
87 priv->slopeDif = fabs(m1 - m2); // (3.2.6.2)
88
89 if ((priv->slopeDif == 0) || ((priv->slopeDif < 2e-15 * fmax(fabs(m1), fabs(m2))) && (priv->i0h[priv->countHigh - 1] - priv->i1h[priv->countHigh - 2] < 100))) {
90 /* new accumulated value of the saved interval is the terminal
91 * value of the current interval
92 * after that dismiss the current interval (3.2.6.3.2.2.) */
93 double x3,y3,x4,y4;
94
95 priv->i0h[index-1] = priv->i0h[index]; /* Remove the second to last element */
96 priv->countHigh--; /* Remove the last element for priv->i1h and priv->mh, priv->xHigh and priv->yHigh */
97
98
99 /* calculation of the new slope (3.2.6.3.3) */
100 x3 = x[priv->i0h[index - 1]]; /* = _dX1 */
101 y3 = y[priv->i0h[index - 1]]; /* = _dY1 */
102 x4 = x[priv->i1h[index - 1]]; /* < X3 */
103 y4 = y[priv->i1h[index - 1]];
104
105 /* write slope to the list of slopes */
106 priv->mh[index - 1] = (y3 - y4) / (x3 - x4);
107
108 } else { /* If difference is too big: ( 3.2.6.4) */
109 priv->xHigh[index] = priv->x2 - (priv->delta * (m1 + m2) / (sqrt((m2 * m2) + (priv->S * priv->S)) + sqrt((m1 * m1) + (priv->S * priv->S))));
110 if (m1 * m2 < 0) {
111 priv->yHigh[index] = priv->y2 + (priv->delta * (m1 * sqrt((m2 * m2) + (priv->S * priv->S)) - m2 * sqrt((m1 * m1) + (priv->S * priv->S)))) / (m1 - m2);
112 } else {
113 priv->yHigh[index] = priv->y2 + (priv->S * priv->S * priv->delta * (m1 + m2) / (m1 * sqrt((m2 * m2) + (priv->S * priv->S)) + m2 * sqrt((m1 * m1) + (priv->S * priv->S))));
114 }
115
116 if ((priv->xHigh[index] == priv->xHigh[index - 1]) && (priv->yHigh[index] != priv->yHigh[index - 1])) {
117 priv->xHigh[index] = priv->xHigh[index - 1] + priv->xMinStep;
118 priv->yHigh[index] = priv->y2 + m1 * (priv->xHigh[index] - priv->x2) + priv->delta * sqrt((m1 * m1) + (priv->S * priv->S));
119 priv->mh[index - 1] = (priv->yHigh[index] - priv->yHigh[index - 1]) / priv->xMinStep;
120 }
121
122 /* If the juncture of the current interval is before the last saved one (3.2.6.7) */
123 while (/* Avoid underflow. Should this generate an error instead? */ index>1 && priv->xHigh[index] <= priv->xHigh[index - 1]) {
124 /* consolidating the current and the previous interval (3.2.6.7.3) */
125 priv->i0h[index-1] = priv->i0h[index]; /* Remove the second to last element */
126 priv->i1h[index-1] = priv->i1h[index];
127 priv->mh[index-1] = priv->mh[index];
128 index--;
129 priv->countHigh--; /* remove also from xHigh,yHigh */
130
131 /* if the saved interval is the 1st interval (3.2.6.7.3.5.1) */
132 if (index == 0) {
133 double x3 = x[0];
134 priv->xHigh[index] = x3 - priv->delta;
135 priv->yHigh[index] = priv->y2 + m1 * (priv->xHigh[index] - priv->x2) + priv->delta * sqrt((m1 * m1) + (priv->S * priv->S));
136 } else { /* if it is not the first: (3.2.6.7.3.5.2.) */
137 double x3 = priv->xHigh[index - 1];
138 double y3 = priv->yHigh[index - 1];
139 m2 = priv->mh[index - 1];
140
141 priv->xHigh[index] = (m2 * x3 - m1 * priv->x2 + priv->y2 - y3 + priv->delta * sqrt((m1 * m1) + (priv->S * priv->S))) / (m2 - m1);
142 priv->yHigh[index] = (m2 * m1 * (x3 - priv->x2) + m2 * (priv->y2 + priv->delta * sqrt((m1 * m1) + (priv->S * priv->S))) - m1 * y3) / (m2 - m1);
143 }
144 }
145 }
146}
147
148static void generateLowTube(privates *priv, double *x, double *y)
149{
150 int index = priv->countLow - 1; /* = _li0l.Count - 1 = _li1l.Count - 1 = xLow.Count - 1 = yLow.Count - 1 > 0 */
151 double m1 = priv->ml[index];
152 double m2 = priv->ml[index - 1];
153 priv->slopeDif = fabs(m1 - m2);
154
155 if ((priv->slopeDif == 0) || ((priv->slopeDif < 2e-15 * fmax(fabs(m1), fabs(m2))) && (priv->i0l[priv->countLow - 1] - priv->i1l[priv->countLow - 2] < 100))) {
156 double x3,y3,x4,y4;
157 priv->i0l[index-1] = priv->i0l[index];
158 priv->countLow--;
159
160 x3 = x[priv->i0l[index - 1]]; /* = _dX1 */
161 y3 = y[priv->i0l[index - 1]]; /* = _dY1 */
162 x4 = x[priv->i1l[index - 1]]; /* < X3 */
163 y4 = y[priv->i1l[index - 1]];
164
165 priv->ml[index - 1] = (y3 - y4) / (x3 - x4);
166 } else {
167 priv->xLow[index] = priv->x2 + (priv->delta * (m1 + m2) / (sqrt((m2 * m2) + (priv->S * priv->S)) + sqrt((m1 * m1) + (priv->S * priv->S))));
168 if (m1 * m2 < 0) {
169 priv->yLow[index] = priv->y2 - (priv->delta * (m1 * sqrt((m2 * m2) + (priv->S * priv->S)) - m2 * sqrt((m1 * m1) + (priv->S * priv->S)))) / (m1 - m2);
170 } else {
171 priv->yLow[index] = priv->y2 - (priv->S * priv->S * priv->delta * (m1 + m2) / (m1 * sqrt((m2 * m2) + (priv->S * priv->S)) + m2 * sqrt((m1 * m1) + (priv->S * priv->S))));
172 }
173
174 if ((priv->xLow[index] == priv->xLow[index - 1]) && (priv->yLow[index] != priv->yLow[index - 1])) {
175 priv->xLow[index] = priv->xLow[index - 1] + priv->xMinStep;
176 priv->yLow[index] = priv->y2 + m1 * (priv->xLow[index] - priv->x2) - priv->delta * sqrt((m1 * m1) + (priv->S * priv->S));
177 priv->ml[index - 1] = (priv->yLow[index] - priv->yLow[index - 1]) / priv->xMinStep;
178 }
179
180 while (index>1 /* Avoid underflow. Should this generate an error instead? */ && priv->xLow[index] <= priv->xLow[index - 1]) {
181 priv->i0l[index-1] = priv->i0l[index];
182 priv->i1l[index-1] = priv->i1l[index];
183 priv->ml[index-1] = priv->ml[index];
184 index--;
185 priv->countLow--;
186
187 if (index == 0) {
188 double x3 = x[0];
189 priv->xLow[index] = x3 - priv->delta;
190 priv->yLow[index] = priv->y2 + m1 * (priv->xLow[index] - priv->x2) - priv->delta * sqrt((m1 * m1) + (priv->S * priv->S));
191 } else {
192 double x3 = priv->xLow[index - 1];
193 double y3 = priv->yLow[index - 1];
194 m2 = priv->ml[index - 1];
195 priv->xLow[index] = (m2 * x3 - m1 * priv->x2 + priv->y2 - y3 - priv->delta * sqrt((m1 * m1) + (priv->S * priv->S))) / (m2 - m1);
196 priv->yLow[index] = (m2 * m1 * (x3 - priv->x2) + m2 * (priv->y2 - priv->delta * sqrt((m1 * m1) + (priv->S * priv->S))) - m1 * y3) / (m2 - m1);
197 }
198 }
199 }
200}
201
202static privates* skipCalculateTubes(double *x, double *y, size_t length)
203{
204 privates *priv = (privates*) omc_alloc_interface.malloc(sizeof(privates));
205
206 /* set tStart and tStop */
207 priv->length = length;
208 priv->tStart = x[0];
209 priv->tStop = x[length - 1];
210 priv->xRelEps = 1e-15;
211 priv->xMinStep = ((priv->tStop - priv->tStart) + fabs(priv->tStart)) * priv->xRelEps;
212 priv->countLow = length;
213 priv->countHigh = length;
214 priv->yHigh = (double*)omc_alloc_interface.malloc_atomic(sizeof(double)*length);
215 priv->yLow = (double*)omc_alloc_interface.malloc_atomic(sizeof(double)*length);
216 memcpy(priv->yHigh, y, length * sizeof(double));
217 memcpy(priv->yLow, y, length * sizeof(double));
218 return priv;
219}
220
221/* This method generates tubes around a given curve */
222static privates* calculateTubes(double *x, double *y, size_t length, double r)
223{
224 privates *priv = (privates*) omc_alloc_interface.malloc(sizeof(privates));
225 int i;
226 /* set tStart and tStop */
227 priv->length = length;
228 priv->tStart = x[0];
229 priv->tStop = x[length - 1];
230 priv->xRelEps = 1e-15;
231 priv->xMinStep = ((priv->tStop - priv->tStart) + fabs(priv->tStart)) * priv->xRelEps;
232 priv->countLow = 0;
233 priv->countHigh = 0;
234
235 /* Initialize lists (upper tube) */
236 priv->mh = (double*)omc_alloc_interface.malloc_atomic(sizeof(double)*length);
237 priv->i0h = (int*)omc_alloc_interface.malloc_atomic(sizeof(int)*length);
238 priv->i1h = (int*)omc_alloc_interface.malloc_atomic(sizeof(int)*length);
239 /* Initialize lists (lower tube) */
240 priv->ml = (double*)omc_alloc_interface.malloc_atomic(sizeof(double)*length);
241 priv->i0l = (int*)omc_alloc_interface.malloc_atomic(sizeof(int)*length);
242 priv->i1l = (int*)omc_alloc_interface.malloc_atomic(sizeof(int)*length);
243
244 priv->xHigh = (double*)omc_alloc_interface.malloc_atomic(sizeof(double)*length);
245 priv->xLow = (double*)omc_alloc_interface.malloc_atomic(sizeof(double)*length);
246 priv->yHigh = (double*)omc_alloc_interface.malloc_atomic(sizeof(double)*length);
247 priv->yLow = (double*)omc_alloc_interface.malloc_atomic(sizeof(double)*length);
248
249 /* calculate the tubes delta */
250 priv->delta = r * (priv->tStop - priv->tStart);
251
252 /* calculate S */
253 priv->max = y[0];
254 priv->min = y[0];
255
256 for (i = 1; i < length; i++) {
257 priv->max = fmax(y[i],priv->max);
258 priv->min = fmin(y[i],priv->min);
259 }
260 priv->S = fabs(4 * (priv->max - priv->min) / (fabs(priv->tStop - priv->tStart)));
261
262 if (priv->S < 0.0004 / fabs(priv->tStop - priv->tStart)) {
263 priv->S = 0.0004 / fabs(priv->tStop - priv->tStart);
264 }
265
266 /* Begin calculation for the tubes */
267 for (i = 1; i < length; i++) {
268 /* get current value */
269 priv->x1 = x[i];
270 priv->y1 = y[i];
271 /* get previous value */
272 priv->x2 = x[i - 1];
273 priv->y2 = y[i - 1];
274 /* catch jumps */
275 if ((priv->x1 <= priv->x2) && (priv->y1 == priv->y2) && (priv->countHigh == 0)) {
276 continue;
277 }
278 if ((priv->x1 <= priv->x2) && (priv->y1 == priv->y2)) {
279 priv->x1 = fmax(priv->x1, x[priv->i1l[priv->countLow - 1]] + priv->xMinStep);
280 priv->x1 = fmax(priv->x1, x[priv->i1h[priv->countHigh - 1]] + priv->xMinStep);
281 x[i] = priv->x1;
282 priv->currentSlope = priv->mh[priv->countHigh - 1];
283 } else {
284 if (priv->x1 <= priv->x2) {
285 priv->x1 = priv->x2 + priv->xMinStep;
286 x[i] = priv->x1;
287 }
288 priv->currentSlope = (priv->y1 - priv->y2) / (priv->x1 - priv->x2); /* calculate current slope ( 3.2.6.1) */
289 }
290
291 /* fill lists with new values: values upper tube */
292 priv->i0h[priv->countHigh] = i;
293 priv->i1h[priv->countHigh] = i-1;
294 priv->mh[priv->countHigh] = priv->currentSlope;
295
296 /* fill lists with new values: values lower tube */
297 priv->i0l[priv->countLow] = i;
298 priv->i1l[priv->countLow] = i-1;
299 if ((priv->x1 <= priv->x2) && (priv->y1 == priv->y2)) {
300 priv->currentSlope = priv->ml[priv->countLow - 1];
301 }
302 priv->ml[priv->countLow] = priv->currentSlope;
303
304 if (priv->countHigh == 0) { /* 1st interval (3.2.5) */
305 /* initial values upper tube */
306 priv->xHigh[priv->countHigh] = priv->x2 - priv->delta;
307 priv->yHigh[priv->countHigh] = priv->y2 - priv->currentSlope * priv->delta + priv->delta * sqrt((priv->currentSlope * priv->currentSlope) + (priv->S * priv->S));
308
309 /* initial values lower tube */
310 priv->xLow[priv->countLow] = priv->x2 - priv->delta;
311 priv->yLow[priv->countLow] = priv->y2 - priv->currentSlope * priv->delta - priv->delta * sqrt((priv->currentSlope * priv->currentSlope) + (priv->S * priv->S));
312
313 priv->countHigh++;
314 priv->countLow++;
315 } else { // if not 1st interval (3.2.6)
316 /* fill lists with new values, set X and Y to arbitrary value (3.2.6.1) */
317 priv->xHigh[priv->countHigh] = 1;
318 priv->yHigh[priv->countHigh] = 1;
319 priv->xLow[priv->countLow] = 1;
320 priv->yLow[priv->countLow] = 1;
321
322 priv->countHigh++;
323 priv->countLow++;
324
325 /* begin procedure for upper tube */
326 generateHighTube(priv, x, y);
327 /* begin procedure for lower tube */
328 generateLowTube(priv, x, y);
329 }
330
331 }
332
333 // calculate terminal value
334 // upper tube
335 priv->x1 = priv->xHigh[priv->countHigh - 1];
336 priv->y1 = priv->yHigh[priv->countHigh - 1];
337 priv->x2 = priv->tStop;
338
339 priv->currentSlope = priv->mh[priv->countHigh - 1];
340
341 priv->xHigh[priv->countHigh] = priv->x2 + priv->delta;
342 priv->yHigh[priv->countHigh] = priv->y1 + priv->currentSlope * (priv->x2 + priv->delta - priv->x1);
343 priv->countHigh++;
344
345 // lower tube
346 priv->x1 = priv->xLow[priv->countLow - 1];
347 priv->y1 = priv->yLow[priv->countLow - 1];
348 priv->x2 = priv->tStop;
349
350 priv->currentSlope = priv->ml[priv->countLow - 1];
351
352 priv->xLow[priv->countLow] = priv->x2 + priv->delta;
353 priv->yLow[priv->countLow] = priv->y1 + priv->currentSlope * (priv->x2 + priv->delta - priv->x1);
354 priv->countLow++;
355 return priv;
356}
357
358static inline__inline__ double linearInterpolation(double x, double x0, double x1, double y0, double y1, double xabstol)
359{
360 if (almostEqualRelativeAndAbs(x0,x,0,xabstol)) { //prevent NaN -> division by zero
361 return y0;
362 } else if (almostEqualRelativeAndAbs(x1,x,0,xabstol)) { //prevent NaN -> division by zero
363 return y1;
364 } else if (almostEqualRelativeAndAbs(x1,x0,0,xabstol)) { //prevent NaN -> division by zero
365 return y0;
366 } else {
367 return y0 + (((y1 - y0) / (x1 - x0)) * (x - x0)); // linear interpolation of the source value at the target moment in time
368 }
369}
370
371/* Calibrate the target time+value pair onto the source timeline */
372static double* calibrateValues(double* sourceTimeLine, double* targetTimeLine, double* targetValues, size_t *nsource, size_t ntarget, double xabstol)
373{
374 double* interpolatedValues;
375 int j, i;
376 double x0, x1, y0, y1;
377 size_t n;
378
379 if (0 == nsource) {
380 return NULL((void*)0);
381 }
382
383 n = *nsource;
384 interpolatedValues = (double*) omc_alloc_interface.malloc_atomic(sizeof(double)*n);
385
386 j = 1;
387 for (i = 0; i < n; i++) {
388 double x = sourceTimeLine[i];
389
390 if (targetTimeLine[j] > sourceTimeLine[n - 1] && targetTimeLine[j-1] > sourceTimeLine[n - 1]) { // Avoid extrapolation by cutting the sequence
391 interpolatedValues[i] = linearInterpolation(x,x0,x1,y0,y1,xabstol);
392 *nsource = i+1;
393 break;
394 }
395
396 x1 = targetTimeLine[j];
397 y1 = targetValues[j];
398
399 while ((x1 <= x) && ((j + 1) < ntarget)) { // step source timline to the current moment
400 j++;
401 x1 = targetTimeLine[j];
402 y1 = targetValues[j];
403 if (almostEqualRelativeAndAbs(x1,x,0,xabstol)) {
404 break;
405 }
406 }
407 x0 = targetTimeLine[j - 1];
408 y0 = targetValues[j - 1];
409 if (i && almostEqualRelativeAndAbs(sourceTimeLine[i-1],x0,0,xabstol) && almostEqualRelativeAndAbs(x0,x1,0,xabstol)) {
410 /* Previous value was the left limit of the event; use the right limit! */
411 interpolatedValues[i] = y1;
412 } else {
413 interpolatedValues[i] = linearInterpolation(x,x0,x1,y0,y1,xabstol);
414 }
415 }
416
417 return interpolatedValues;
418}
419
420typedef struct {
421 double *time;
422 double *values;
423 size_t size;
424} addTargetEventTimesRes;
425
426static size_t findNextEvent(size_t i, double *time, size_t n, double xabstol)
427{
428 for (; i < n; i++) {
429 if (almostEqualRelativeAndAbs(time[i-1],time[i],0,xabstol) || (i<n-1 && almostEqualRelativeAndAbs(time[i],time[i+1],0,xabstol))) return i;
430 }
431 return 0; /* Not found */
432}
433
434static addTargetEventTimesRes addTargetEventTimes(double* sourceTimeLine, double* targetTimeLine, double *sourceValues, size_t nsource, size_t ntarget, double xabstol)
435{
436 addTargetEventTimesRes res;
437 size_t i=0,j,count=0;
438 int iter=0;
439 while ((i=findNextEvent(i+1,targetTimeLine,ntarget,xabstol))) {
440 if (targetTimeLine[i] >= sourceTimeLine[nsource-1]) {
441 break;
442 }
443 count++; /* The number of found time events in the target file */
444 }
445 if (count == 0) {
446 res.size = nsource;
447 res.time = sourceTimeLine;
448 res.values = sourceValues;
449 return res;
450 }
451 res.size = nsource+count;
452 res.values = omc_alloc_interface.malloc_atomic(sizeof(double)*res.size);
453 res.time = omc_alloc_interface.malloc_atomic(sizeof(double)*res.size);
454 i=0;
455 count=0;
456 j=findNextEvent(1,targetTimeLine,ntarget,xabstol);
457 while (j) {
458 if (targetTimeLine[j] >= sourceTimeLine[nsource-1]) {
459 break;
460 }
461 while (sourceTimeLine[i] < targetTimeLine[j] && i<nsource) {
462 res.values[count] = sourceValues[i];
463 res.time[count++] = sourceTimeLine[i++];
464 }
465 if (sourceTimeLine[i] == targetTimeLine[j]) {
466 res.size--; /* Filter events at identical times in both files */
467 } else {
468 double x0 = sourceTimeLine[intmax(0,i-1)];
469 double y0 = sourceValues[intmax(0,i-1)];
470 double x1 = sourceTimeLine[i];
471 double y1 = sourceValues[i];
472 double x = targetTimeLine[j];
473 double y = linearInterpolation(x,x0,x1,y0,y1,xabstol);
474 res.values[count] = y;
475 res.time[count++] = x;
476 }
477 assert(count < res.size)((void) sizeof ((count < res.size) ? 1 : 0), __extension__
({ if (count < res.size) ; else __assert_fail ("count < res.size"
, "./SimulationResultsCmpTubes.c", 477, __extension__ __PRETTY_FUNCTION__
); }))
;
478 j=findNextEvent(j+1,targetTimeLine,ntarget,xabstol);
479 iter++;
480 }
481 while (i < nsource) {
482 res.values[count] = sourceValues[i];
483 res.time[count++] = sourceTimeLine[i++];
484 }
485 assert(res.size == count)((void) sizeof ((res.size == count) ? 1 : 0), __extension__ (
{ if (res.size == count) ; else __assert_fail ("res.size == count"
, "./SimulationResultsCmpTubes.c", 485, __extension__ __PRETTY_FUNCTION__
); }))
;
486 return res;
487}
488
489static addTargetEventTimesRes mergeTimelines(addTargetEventTimesRes ref, addTargetEventTimesRes actual, double xabstol)
490{
491 int i=0,j=0,count=0;
492 addTargetEventTimesRes res;
493 res.size = ref.size + actual.size;
494 res.values = omc_alloc_interface.malloc_atomic(sizeof(double)*res.size);
495 res.time = omc_alloc_interface.malloc_atomic(sizeof(double)*res.size);
496 res.values[count] = ref.values[0];
497 res.time[count++] = ref.time[0];
498 for (i=1; i<ref.size; i++) {
499 double x0 = ref.time[i-1];
500 double y0 = ref.values[i-1];
501 double x1 = ref.time[i];
502 double y1 = ref.values[i];
503 while (j<actual.size && actual.time[j] <= x1) {
504 double x = actual.time[j];
505 int isEventInRef = almostEqualRelativeAndAbs(x0,x1,0,xabstol);
506 int isEventInActual = j<actual.size-1 && almostEqualRelativeAndAbs(x,actual.time[j+1],0,xabstol);
507 int isEvent = isEventInRef || isEventInActual;
508 if ((almostEqualRelativeAndAbs(x,x0,0,xabstol) || almostEqualRelativeAndAbs(x,x1,0,xabstol)) && !isEvent) {
509 /* Not an event but we got the same time stamp; skip! */
510 j++;
511 continue;
512 }
513 res.time[count] = x;
514 res.values[count++] = linearInterpolation(x,x0,x1,y0,y1,xabstol);
515 j++;
516 }
517 res.time[count] = x1;
518 res.values[count++] = y1;
519 assert(count <= res.size)((void) sizeof ((count <= res.size) ? 1 : 0), __extension__
({ if (count <= res.size) ; else __assert_fail ("count <= res.size"
, "./SimulationResultsCmpTubes.c", 519, __extension__ __PRETTY_FUNCTION__
); }))
;
520 }
521 res.size = count;
522 return res;
523}
524
525static addTargetEventTimesRes removeUneventfulPoints(addTargetEventTimesRes in, double reltol, double xabstol)
526{
527 int i;
528 addTargetEventTimesRes res;
529 res.values = (double*) omc_alloc_interface.malloc_atomic(in.size * sizeof(double));
530 res.time = (double*) omc_alloc_interface.malloc_atomic(in.size * sizeof(double));
531
532 do {
533 int iter = 0;
534 /* Don't remove first point */
535 res.values[0] = in.values[0];
536 res.time[0] = in.time[0];
537 res.size = 1;
538 for (i=1; i<in.size-1; i++) {
539 /* double x0 = res.time[res.size-1]; */
540 double y0 = res.values[res.size-1];
541 double x = in.time[i];
542 double y = in.values[i];
543 double x1 = in.time[i+1];
544 double y1 = in.values[i+1];
545 if (y0 == y1 && y == y0) {
546 res.time[res.size] = x1;
547 res.values[res.size] = y1;
548 res.size++;
549 i++;
550 if (i < in.size-2) {
551 iter=1;
552 }
553 continue;
554 }
555 res.time[res.size] = x;
556 res.values[res.size] = y;
557 res.size++;
558 }
559 if (in.size > 1) {
560 /* Don't remove last point */
561 res.values[res.size] = in.values[in.size-1];
562 res.time[res.size] = in.time[in.size-1];
563 res.size++;
564 }
565 if (iter) {
566 /* This is ok, because we only remove the current element we iterate over; we could do this in-place but we avoid copying the initial array (which we want to remain as is) */
567 in = res;
568 reltol = 0;
569 } else {
570 /* no more recursion */
571 break;
572 }
573 } while (1);
574 return res;
575}
576
577/* Adds a relative tolerance compared to the reference signal. Overwrites the target values vector. */
578static void addRelativeTolerance(double *targetValues, double *sourceValues, size_t length, double reltol, double abstol, int direction)
579{
580 int i;
581 if (direction > 0) {
582 for (i=0; i<length; i++) {
583 targetValues[i] = fmax(sourceValues[i] + fmax(fabs(sourceValues[i]*reltol),abstol), targetValues[i]);
584 }
585 } else {
586 for (i=0; i<length; i++) {
587 targetValues[i] = fmin(sourceValues[i] - fmax(fabs(sourceValues[i]*reltol),abstol), targetValues[i]);
588 }
589 }
590}
591
592static void assertMonotonic(addTargetEventTimesRes series)
593{
594 int i;
595 for (i=1; i<series.size; i++) {
596 if (series.time[i] < series.time[i-1]) {
597 printf("assertion failed, size %ld: %.15g < %.15g\n", series.size, series.time[i], series.time[i-1]);
598 abort();
599 }
600 }
601}
602
603/* Return NULL if there were no errors */
604static double* validate(int n, addTargetEventTimesRes ref, double *low, double *high, double *calibrated_values, double reltol, double abstol, double xabstol)
605{
606 double *error = omc_alloc_interface.malloc_atomic(n * sizeof(double));
607 int isdifferent = 0;
608 int i,lastStepError = 1;
609 for (i=0; i<n; i++) {
610 int thisStepError = 0;
611 int isEvent = (i && almostEqualRelativeAndAbs(ref.time[i],ref.time[i-1],0,xabstol)) || (i+1<n && almostEqualRelativeAndAbs(ref.time[i],ref.time[i+1],0,xabstol));
612 if (isEvent) {
613 double refv = ref.values[i];
614 double val = calibrated_values[i];
615 double tol = fmax(abstol*10,fmax(fabs(refv),fabs(val))*reltol*10);
616 /* If there was no error in the last step before the event, give a nice tight curve around both values */
617 high[i] = (lastStepError ? refv : fmax(refv,val)) + tol;
618 low[i] = (lastStepError ? refv : fmin(refv,val)) - tol;
619 error[i] = NAN(__builtin_nanf (""));
620 } else {
621 error[i] = 0;
622 thisStepError=lastStepError;
623 if (calibrated_values[i] < low[i]) {
624 error[i] = low[i]-calibrated_values[i];
625 isdifferent++;
626 thisStepError=1;
627 } else if (calibrated_values[i] > high[i]) {
628 error[i] = calibrated_values[i]-high[i];
629 isdifferent++;
630 thisStepError=1;
631 }
632 }
633 lastStepError = thisStepError;
634 }
635 if (isdifferent) {
636 return error;
637 }
638 return NULL((void*)0);
639}
640
641static unsigned int cmpDataTubes(int isResultCmp, char* varname, DataField *time, DataField *reftime, DataField *data, DataField *refdata, double reltol, double rangeDelta, double reltolDiffMaxMin, DiffDataField *ddf, char **cmpdiffvars, unsigned int vardiffindx, int keepEqualResults, void **diffLst, const char *prefix, int isHtml, char **htmlOut)
642{
643 int withTubes = 0 == rangeDelta;
644 FILE *fout = NULL((void*)0);
645 char *fname = NULL((void*)0);
646 char *html;
647 /* The tolerance for detecting events is proportional to the number of output points in the file */
648 double xabstol = (reftime->data[reftime->n-1]-reftime->data[0])*(withTubes ? rangeDelta : 1e-3) / fmax(time->n,reftime->n);
649 /* Calculate the tubes without additional events added */
650 addTargetEventTimesRes ref,actual,actualoriginal;
651 privates *priv=NULL((void*)0);
652 size_t n,maxn,html_size=0;
653 double *calibrated_values=NULL((void*)0), *high=NULL((void*)0), *low=NULL((void*)0), *error=NULL((void*)0),maxPlusTol,minMinusTol,abstol;
654
655 ref.values = refdata->data;
656 ref.time = reftime->data;
657 ref.size = reftime->n;
658 actualoriginal.values = data->data;
659 actualoriginal.time = time->data;
660 actualoriginal.size = time->n;
661 actual = actualoriginal;
662 /* assertMonotonic(ref); */
663 /* assertMonotonic(actual); */
664 /* ref = removeUneventfulPoints(ref, reltol*reltol, xabstol); */
665 /* actual = removeUneventfulPoints(actual, reltol*reltol, xabstol); */
666 /* assertMonotonic(ref); */
667 /* assertMonotonic(actual); */
668 priv = withTubes ? skipCalculateTubes(ref.time,ref.values,ref.size) : calculateTubes(ref.time,ref.values,ref.size,rangeDelta);
669 /* ref = mergeTimelines(ref,actual,xabstol); */
670 /* assertMonotonic(ref); */
671 n = ref.size;
672 calibrated_values = calibrateValues(ref.time,actual.time,actual.values,&n,actual.size,xabstol);
673 maxPlusTol = priv->max + fabs(priv->max) * reltol;
674 minMinusTol = priv->min - fabs(priv->min) * reltol;
675 high = calibrateValues(ref.time,priv->xHigh,priv->yHigh,&n,priv->countHigh,xabstol);
676 low = calibrateValues(ref.time,priv->xLow,priv->yLow,&n,priv->countLow,xabstol);
677 /* If all values in the reference are ~0 (and the same)... Allow reltolDiffMaxMin^2 as tolerance
678 * Maybe we should just treat it differently though
679 * Like not creating a tubes and simply check that the other file also has only identical points close to this
680 */
681 abstol = (priv->max-priv->min == 0 && priv->max < reltolDiffMaxMin*reltolDiffMaxMin) ? reltolDiffMaxMin*reltolDiffMaxMin : fabs((priv->max-priv->min)*reltolDiffMaxMin);
682 addRelativeTolerance(high,ref.values,n,reltol,abstol,1);
683 addRelativeTolerance(low ,ref.values,n,reltol,abstol,-1);
684 error = validate(n,ref,low,high,calibrated_values,reltol,abstol,xabstol);
685 if ( isHtml ) {
686
687#if _XOPEN_SOURCE >= 700 || _POSIX_C_SOURCE200809L >= 200809L
688 html_size=0;
689 fout = open_memstream(&html, &html_size);
690#else
691 fname = (char*) omc_alloc_interface.malloc_atomic(25 + strlen(varname));
692 sprintf(fname, "tmp.%s.html.tmp", varname);
693 fout = fopen(fname, "wb+");
694 if (!fout)
695 {
696 perror("Error opening temp file"); fflush(stderrstderr);
697 }
698#endif
699
700 fprintf(fout, "<html>\n"
701"<head>\n"
702"<script type=\"text/javascript\" src=\"dygraph-combined.js\"></script>\n"
703" <style type=\"text/css\">\n"
704" #graphdiv {\n"
705" position: absolute;\n"
706" left: 10px;\n"
707" right: 10px;\n"
708" top: 40px;\n"
709" bottom: 10px;\n"
710" }\n"
711" </style>\n"
712"</head>\n"
713"<body>\n"
714"<div id=\"graphdiv\"></div>\n"
715"<p>"
716"<input type=checkbox id=\"0\" checked onClick=\"change(this)\">\n"
717"<label for=\"0\">reference</label>\n"
718"<input type=checkbox id=\"1\" checked onClick=\"change(this)\">\n"
719"<label for=\"1\">actual</label>\n"
720"<input type=checkbox id=\"2\" checked onClick=\"change(this)\">\n"
721"<label for=\"2\">high</label>\n"
722"<input type=checkbox id=\"3\" checked onClick=\"change(this)\">\n"
723"<label for=\"3\">low</label>\n"
724"<input type=checkbox id=\"4\" checked onClick=\"change(this)\">\n"
725"<label for=\"4\">error</label>\n"
726"<input type=checkbox id=\"5\" onClick=\"change(this)\">\n"
727"<label for=\"5\">actual (original)</label>\n"
728"Reference time: %.15g to %.15g, actual time: %.15g to %.15g. Parameters used for the comparison: Relative tolerance %.2g. Absolute tolerance %.2g (%.2g relative). Range delta %.2g."
729"</p>\n"
730"<script type=\"text/javascript\">\n"
731"g = new Dygraph(document.getElementById(\"graphdiv\"),\n"
732"[\n",
733 reftime->data[0],
734 reftime->data[reftime->n-1],
735 time->data[0],
736 time->data[time->n-1],
737 reltol,
738 abstol,
739 reltolDiffMaxMin,
740 rangeDelta
741);
742 } else if (!isResultCmp && (error || keepEqualResults)) {
743 fname = (char*) omc_alloc_interface.malloc_atomic(25 + strlen(prefix) + strlen(varname));
744 sprintf(fname, "%s.%s.csv", prefix, varname);
745 fout = fopen(fname,"w");
746 if (!fout)
747 {
748 perror("Error opening file");
749 fprintf(stderrstderr, "File: %s\n", fname);
750 fflush(stderrstderr);
751 }
752 }
753 maxn = intmax(intmax(intmax(ref.size,actual.size),priv->countHigh),priv->countLow);
Value stored to 'maxn' is never read
754 if (fout) {
755 int i;
756 const char *empty = isHtml ? "null" : "";
757 const char *lbracket = isHtml ? "[" : "";
758 const char *rbracket = isHtml ? "]," : "";
759 int j=0;
760 int lastStepError = 1;
761 if (!isHtml) {
762 fputs("time,reference,actual,high,low,error,actual (raw)\n", fout);
763 }
764 for (i=0; i<ref.size; i++) {
765 int thisStepError = 0;
766 fprintf(fout, "%s%.15g,%.15g,",lbracket,ref.time[i],ref.values[i]);
767 if (i < n) {
768 if (!error || isnan(error[i])(sizeof ((error[i])) == sizeof (float) ? __isnanf (error[i]) :
sizeof ((error[i])) == sizeof (double) ? __isnan (error[i]) :
__isnanl (error[i]))
) {
769 fprintf(fout, "%.15g,%.15g,%.15g,%s",calibrated_values[i],high[i],low[i],empty);
770 } else {
771 fprintf(fout, "%.15g,%.15g,%.15g,%.15g",calibrated_values[i],high[i],low[i],error[i]);
772 }
773 if (j < actualoriginal.size && ref.time[i] == actualoriginal.time[j]) {
774 fprintf(fout, ",%.15g%s\n",actualoriginal.values[j++],rbracket);
775 } else {
776 fprintf(fout, ",%s%s\n",empty,rbracket);
777 }
778 } else {
779 fputs(isHtml ? "null,null,null,null,null],\n" : ",,,,\n", fout);
780 }
781 while (j < actualoriginal.size && ref.time[i] > actualoriginal.time[j]) {
782 fprintf(fout, "%s%.15g,%s,%s,%s,%s,%s,%.15g%s\n",lbracket,actualoriginal.time[j],empty,empty,empty,empty,empty,actualoriginal.values[j],rbracket);
783 j++;
784 }
785 lastStepError = thisStepError;
786 }
787 fputs(isHtml ? "],\n" : "\n", fout);
788 }
789 if (error) {
790 cmpdiffvars[vardiffindx] = varname;
791 vardiffindx++;
792 if (!isResultCmp) {
793 *diffLst = mmc_mk_cons(mmc_mk_scon(varname),*diffLst);
794 }
795 }
796 if (fout) {
797 if (isHtml) {
798fprintf(fout, "{title: '%s',\n"
799"legend: 'always',\n"
800"xlabel: ['time'],\n"
801"connectSeparatedPoints: true,\n"
802"labels: ['time','reference','actual','high','low','error','actual (original)'],\n"
803"y2label: ['error'],\n"
804"series : { 'error': { axis: 'y2' } },\n"
805"colors: ['blue','red','teal','lightblue','orange','black'],\n"
806"visibility: [true,true,true,true,true,false]\n"
807"});\n"
808"function change(el) {\n"
809" g.setVisibility(parseInt(el.id), el.checked);\n"
810"}\n"
811"</script>\n"
812"</body>\n"
813"</html>\n", varname);
814 }
815
816 if (isHtml) {
817#if !(_XOPEN_SOURCE >= 700 || _POSIX_C_SOURCE200809L >= 200809L)
818 size_t r;
819 if (fseek(fout, 0, SEEK_END2))
820 {
821 perror("Error on fseek end!");
822 }
823 html_size = ftell(fout);
824 if (fseek(fout, 0, SEEK_SET0))
825 {
826 perror("Error on fseek set!");
827 }
828 html = (char*)malloc((html_size + 1) * sizeof(char));
829 r = fread(html, sizeof(char), html_size, fout);
830 if (r != html_size)
831 {
832 perror("Error on fread!");
833 }
834 html[html_size] = '\0';
835 fclose(fout);
836 unlink(fname);
837#else
838 fclose(fout);
839#endif
840 *htmlOut = omc_alloc_interface.malloc_strdup(html);
841 free(html);
842 }
843 else
844 {
845 fclose(fout);
846 }
847 }
848 /* Tell the GC some variables have been free'd */
849 if (error) GC_free(error);
850 if (fname) GC_free(fname);
851 GC_free(low);
852 GC_free(high);
853 if (withTubes) {
854 GC_free(priv->mh);
855 GC_free(priv->i0h);
856 GC_free(priv->i1h);
857 GC_free(priv->ml);
858 GC_free(priv->i0l);
859 GC_free(priv->i1l);
860 GC_free(priv->xHigh);
861 GC_free(priv->xLow);
862 }
863 GC_free(priv->yHigh);
864 GC_free(priv->yLow);
865 GC_free(priv);
866 GC_free(calibrated_values);
867 return vardiffindx;
868}