Bug Summary

File:OMCompiler/Compiler/runtime/SimulationResults_omc.c
Warning:line 391, column 31
2nd function call argument is an uninitialized value

Annotated Source Code

[?] Use j/k keys for keyboard navigation

SimulationResults_omc.c

1/*
2 * This file is part of OpenModelica.
3 *
4 * Copyright (c) 1998-2010, Linköpings University,
5 * 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 THIS OSMC PUBLIC
11 * LICENSE (OSMC-PL). ANY USE, REPRODUCTION OR DISTRIBUTION OF
12 * THIS PROGRAM CONSTITUTES RECIPIENT'S ACCEPTANCE OF THE OSMC
13 * PUBLIC LICENSE.
14 *
15 * The OpenModelica software and the Open Source Modelica
16 * Consortium (OSMC) Public License (OSMC-PL) are obtained
17 * from Linköpings University, either from the above address,
18 * from the URL: http://www.ida.liu.se/projects/OpenModelica
19 * and in the OpenModelica distribution.
20 *
21 * This program is distributed WITHOUT ANY WARRANTY; without
22 * even the implied warranty of MERCHANTABILITY or FITNESS
23 * FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY SET FORTH
24 * IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS
25 * OF OSMC-PL.
26 *
27 * See the full OSMC Public License conditions for more details.
28 *
29 */
30
31#if defined(_MSC_VER) || defined(__MINGW32__)
32 #define WIN32_LEAN_AND_MEAN
33 #include <windows.h>
34#endif
35
36#include <stdio.h>
37#include <stdlib.h>
38
39
40#include "meta_modelica.h"
41#define ADD_METARECORD_DEFINITIONSstatic static
42#include "OpenModelicaBootstrappingHeader.h"
43
44#include "SimulationResults.c"
45#include "SimulationResultsCmp.c"
46
47void* SimulationResults_readVariables(const char *filename, int readParameters, int omcStyle)
48{
49 return SimulationResultsImpl__readVars(filename, readParameters, omcStyle, &simresglob);
50}
51
52extern void* _ValuesUtil_reverseMatrix(void*);
53void* SimulationResults_readDataset(const char *filename, void *vars, int datasize)
54{
55 void *res = SimulationResultsImpl__readDataset(filename,vars,datasize,0,&simresglob,0);
56 if (res == NULL((void*)0)) MMC_THROW(){longjmp(*((threadData_t*)pthread_getspecific(mmc_thread_data_key
))->mmc_jumper,1);}
;
57 return res;
58}
59
60int SimulationResults_readSimulationResultSize(const char *filename)
61{
62 return SimulationResultsImpl__readSimulationResultSize(filename,&simresglob);
63}
64
65double SimulationResults_val(const char *filename, const char *varname, double timeStamp)
66{
67 return SimulationResultsImpl__val(filename,varname,timeStamp,&simresglob);
68}
69
70void* SimulationResults_cmpSimulationResults(int runningTestsuite, const char *filename,const char *reffilename,const char *logfilename, double refTol, double absTol, void *vars)
71{
72 return SimulationResultsCmp_compareResults(1,runningTestsuite,filename,reffilename,logfilename,refTol,absTol,0,0,vars,0,NULL((void*)0),0,NULL((void*)0));
73}
74
75double SimulationResults_deltaSimulationResults(const char *filename,const char *reffilename, const char *methodname, void *vars)
76{
77 double res = SimulationResultsCmp_deltaResults(filename,reffilename,methodname,vars);
78 return res;
79}
80
81void* SimulationResults_diffSimulationResults(int runningTestsuite, const char *filename,const char *reffilename,const char *logfilename, double refTol, double reltolDiffMaxMin, double rangeDelta, void *vars, int keepEqualResults, int *success)
82{
83 return SimulationResultsCmp_compareResults(0,runningTestsuite,filename,reffilename,logfilename,refTol,0,reltolDiffMaxMin,rangeDelta,vars,keepEqualResults,success,0,NULL((void*)0));
84}
85
86const char* SimulationResults_diffSimulationResultsHtml(int runningTestsuite, const char *var, const char *filename,const char *reffilename, double refTol, double reltolDiffMaxMin, double rangeDelta)
87{
88 char *res = "";
89 SimulationResultsCmp_compareResults(0,runningTestsuite,filename,reffilename,"",0,refTol,reltolDiffMaxMin,rangeDelta,mmc_mk_cons(mmc_mk_scon(var),mmc_mk_nil()((void*)((char*)(&(mmc_nil).header) + 3))),0,NULL((void*)0),1,&res);
1
Calling 'SimulationResultsCmp_compareResults'
90 return res;
91}
92
93void SimulationResults_close()
94{
95 SimulationResultsImpl__close(&simresglob);
96}

./SimulationResultsCmp.c

1/*
2* This file is part of OpenModelica.
3*
4* Copyright (c) 1998-2010, Linköpings University,
5* 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 THIS OSMC PUBLIC
11* LICENSE (OSMC-PL). ANY USE, REPRODUCTION OR DISTRIBUTION OF
12* THIS PROGRAM CONSTITUTES RECIPIENT'S ACCEPTANCE OF THE OSMC
13* PUBLIC LICENSE.
14*
15* The OpenModelica software and the Open Source Modelica
16* Consortium (OSMC) Public License (OSMC-PL) are obtained
17* from Linköpings University, either from the above address,
18* from the URL: http://www.ida.liu.se/projects/OpenModelica
19* and in the OpenModelica distribution.
20*
21* This program is distributed WITHOUT ANY WARRANTY; without
22* even the implied warranty of MERCHANTABILITY or FITNESS
23* FOR A PARTICULAR PURPOSE, EXCEPT AS EXPRESSLY SET FORTH
24* IN THE BY RECIPIENT SELECTED SUBSIDIARY LICENSE CONDITIONS
25* OF OSMC-PL.
26*
27* See the full OSMC Public License conditions for more details.
28*
29*/
30
31#include <stdio.h>
32#include <stdlib.h>
33#include <stdint.h>
34#include <errno(*__errno_location ()).h>
35#include <string.h>
36#include <assert.h>
37
38#include "systemimpl.h"
39
40/* Size of the buffer for warnings and other messages */
41#define WARNINGBUFFSIZE4096 4096
42
43typedef struct {
44 double *data;
45 unsigned int n;
46} DataField;
47
48typedef struct {
49 char *name;
50 double data;
51 double dataref;
52 double time;
53 double timeref;
54 char interpolate;
55} DiffData;
56
57typedef struct {
58 DiffData *data;
59 unsigned int n;
60 unsigned int n_max;
61} DiffDataField;
62
63typedef enum {
64 NORM1,
65 NORM2,
66 MAX_ERR
67} ErrorMethod;
68
69#define DOUBLEEQUAL_TOTAL0.0000000001 0.0000000001
70#define DOUBLEEQUAL_REL0.00001 0.00001
71
72static SimulationResult_Globals simresglob_c = {UNKNOWN_PLOT,0};
73static SimulationResult_Globals simresglob_ref = {UNKNOWN_PLOT,0};
74
75static char ** getVars(void *vars, unsigned int* nvars)
76{
77 char **cmpvars = NULL((void*)0);
78 unsigned int ncmpvars = 0;
79 unsigned int i = 0;
80 void *v;
81
82 /* Count the number of variables in the list. */
83 for(v = vars; MMC_NILHDR(((0) << 10) + (((0) & 255) << 2)) != MMC_GETHDR(v)(*(mmc_uint_t*)((void*)((char*)(v) - 3))); v = MMC_CDR(v)(*(void**)(((void*)((void**)(((void*)((char*)(v) - 3))) + (2)
))))
) {
84 ++ncmpvars;
85 }
86
87 /* Allocate a new string array to contain the variable names. */
88 cmpvars = (char**)omc_alloc_interface.malloc(sizeof(char*)*(ncmpvars));
89
90 /* Copy the variable names from the MetaModelica list to the string array. */
91 for(; MMC_NILHDR(((0) << 10) + (((0) & 255) << 2)) != MMC_GETHDR(vars)(*(mmc_uint_t*)((void*)((char*)(vars) - 3))); vars = MMC_CDR(vars)(*(void**)(((void*)((void**)(((void*)((char*)(vars) - 3))) + (
2)))))
) {
92 cmpvars[i++] = omc_alloc_interface.malloc_strdup(MMC_STRINGDATA(MMC_CAR(vars))(((struct mmc_string*)((void*)((char*)((*(void**)(((void*)((void
**)(((void*)((char*)(vars) - 3))) + (1)))))) - 3)))->data)
);
93 }
94
95 *nvars = ncmpvars;
96 return cmpvars;
97}
98
99static DataField getData(const char *varname,const char *filename, unsigned int size, int suggestRealAll, SimulationResult_Globals* srg, int runningTestsuite)
100{
101 DataField res;
102 void *cmpvar,*dataset,*lst,*datasetBackup;
103 unsigned int i;
104 res.n = 0;
105 res.data = NULL((void*)0);
106
107 /* fprintf(stderr, "getData of Var: %s from file %s\n", varname,filename); */
108 cmpvar = mmc_mk_nil()((void*)((char*)(&(mmc_nil).header) + 3));
109 cmpvar = mmc_mk_cons(mmc_mk_scon(varname),cmpvar);
110 dataset = SimulationResultsImpl__readDataset(filename,cmpvar,size,suggestRealAll,srg,runningTestsuite);
111 if (dataset==NULL((void*)0)) {
112 /* fprintf(stderr, "getData of Var: %s failed!\n",varname); */
113 return res;
114 }
115
116 /* fprintf(stderr, "Data of Var: %s\n", varname); */
117 /* First calculate the length of the matrix */
118 datasetBackup = dataset;
119 while (MMC_NILHDR(((0) << 10) + (((0) & 255) << 2)) != MMC_GETHDR(dataset)(*(mmc_uint_t*)((void*)((char*)(dataset) - 3)))) {
120 lst = MMC_CAR(dataset)(*(void**)(((void*)((void**)(((void*)((char*)(dataset) - 3)))
+ (1)))))
;
121 while (MMC_NILHDR(((0) << 10) + (((0) & 255) << 2)) != MMC_GETHDR(lst)(*(mmc_uint_t*)((void*)((char*)(lst) - 3)))) {
122 res.n++;
123 lst = MMC_CDR(lst)(*(void**)(((void*)((void**)(((void*)((char*)(lst) - 3))) + (
2)))))
;
124 }
125 dataset = MMC_CDR(dataset)(*(void**)(((void*)((void**)(((void*)((char*)(dataset) - 3)))
+ (2)))))
;
126 }
127 if (res.n == 0) return res;
128
129 /* The allocate and read the values */
130 dataset = datasetBackup;
131 i = res.n;
132 res.data = (double*) malloc(sizeof(double)*res.n);
133 while (MMC_NILHDR(((0) << 10) + (((0) & 255) << 2)) != MMC_GETHDR(dataset)(*(mmc_uint_t*)((void*)((char*)(dataset) - 3)))) {
134 lst = MMC_CAR(dataset)(*(void**)(((void*)((void**)(((void*)((char*)(dataset) - 3)))
+ (1)))))
;
135 while (MMC_NILHDR(((0) << 10) + (((0) & 255) << 2)) != MMC_GETHDR(lst)(*(mmc_uint_t*)((void*)((char*)(lst) - 3)))) {
136 res.data[--i] = mmc_prim_get_real(MMC_CAR(lst)(*(void**)(((void*)((void**)(((void*)((char*)(lst) - 3))) + (
1)))))
);
137 lst = MMC_CDR(lst)(*(void**)(((void*)((void**)(((void*)((char*)(lst) - 3))) + (
2)))))
;
138 }
139 dataset = MMC_CDR(dataset)(*(void**)(((void*)((void**)(((void*)((char*)(dataset) - 3)))
+ (2)))))
;
140 }
141 assert(i == 0)((void) sizeof ((i == 0) ? 1 : 0), __extension__ ({ if (i == 0
) ; else __assert_fail ("i == 0", "./SimulationResultsCmp.c",
141, __extension__ __PRETTY_FUNCTION__); }))
;
142
143 /* for (i=0;i<res.n;i++)
144 fprintf(stderr, "%d: %.15g\n", i, res.data[i]); */
145
146 return res;
147}
148
149/* see http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/ */
150static char almostEqualRelativeAndAbs(double a, double b, double reltol, double abstol)
151{
152 /* Check if the numbers are really close -- needed when comparing numbers near zero. */
153 double diff = fabs(a - b);
154 if (diff <= abstol) {
155 return 1;
156 }
157 if (diff <= fmax(fabs(a),fabs(b)) * reltol) {
158 return 1;
159 }
160 return 0;
161}
162
163static char almostEqualWithDefaultTolerance(double a, double b)
164{
165 return almostEqualRelativeAndAbs(a,b,DOUBLEEQUAL_REL0.00001,DOUBLEEQUAL_TOTAL0.0000000001);
166}
167
168static double deltaData(ErrorMethod errMethod, DataField *time, DataField *reftime, DataField *data, DataField *refdata)
169{
170 unsigned int i, iRef, i2;
171 double res,res0, val, valRef, t;
172 res = 0;
173 res0 = 0;
174 i = 0;
175 for (iRef=0;iRef < reftime->n;iRef++){
176 valRef = refdata->data[iRef];
177 t = reftime->data[iRef];
178 //get left and right reference point to interpolate
179 while((time->data[i] < t))
180 {
181 i++;
182 }
183 i2 = i+1;
184 //there is a result value at time t
185 if (fabs(time->data[i]-t) <= fmax((0.0001*time->data[time->n]),1e-12))
186 {
187 val = data->data[i];
188 }
189 //interpolate result value at time t
190 else
191 {
192 //printf("xl %f xr %f tl %f tr %f \n", data->data[i], data->data[i2], time->data[i], time->data[i2]);
193 val = (data->data[i2] - data->data[i])/(time->data[i2] - time->data[i])*(t - time->data[i])+data->data[i];
194 }
195 switch(errMethod){
196 case NORM1:
197 res0 += fabs(valRef-val); break;
198 case NORM2:
199 res0 += pow((valRef-val),2); break;
200 case MAX_ERR:
201 res0 = fmax(fabs(valRef-val),res0); break;
202 default:
203 res0 += fabs(valRef-val); break;
204 }
205 //fprintf(stderr, "at time %f, val: %f and valRef: %f and res %f\n",t,val,valRef, res0);
206 }
207 switch(errMethod){
208 case NORM2:
209 res = sqrt(res0); break;
210 default:
211 res = res0; break;
212 }
213 return res;
214}
215
216
217static unsigned int cmpData(int isResultCmp, char* varname, DataField *time, DataField *reftime, DataField *data, DataField *refdata, double reltol, double abstol, DiffDataField *ddf, char **cmpdiffvars, unsigned int vardiffindx, int keepEqualResults, void **diffLst, const char *prefix)
218{
219 unsigned int i,j,k,j_event;
220 double t,tr,d,dr,err,d_left,d_right,dr_left,dr_right,t_event;
221 char increased = 0;
222 char interpolate = 0;
223 char isdifferent = 0;
224 char refevent = 0;
225 double average=0;
226 FILE *fout = NULL((void*)0);
227 char *fname = NULL((void*)0);
228 if (!isResultCmp) {
229 fname = (char*) malloc(25 + strlen(prefix) + strlen(varname));
230 sprintf(fname, "%s.%s.csv", prefix, varname);
231 fout = fopen(fname,"w");
232 if (fout) {
233 fprintf(fout, "time,reference,actual,err,relerr,threshold\n");
234 }
235 }
236 for (i=0;i<refdata->n;i++){
237 average += fabs(refdata->data[i]);
238 }
239 average = average/((double)refdata->n);
240#ifdef DEBUGOUTPUT
241 fprintf(stderrstderr, "average: %.15g\n",average);
242#endif
243 average = reltol*fabs(average)+abstol;
244 j = 0;
245 tr = reftime->data[j];
246 dr = refdata->data[j];
247#ifdef DEBUGOUTPUT
248 fprintf(stderrstderr, "compare: %s\n",varname);
249#endif
250 for (i=0;i<data->n;i++){
251 t = time->data[i];
252 d = data->data[i];
253 increased = 0;
254#ifdef DEBUGOUTPUT
255 fprintf(stderrstderr, "i: %d t: %.15g d:%.15g\n",i,t,d);
256#endif
257 while(tr < t){
258 if (j +1< reftime->n) {
259 j += 1;
260 tr = reftime->data[j];
261 increased = 1;
262 if (tr == t) {
263 break;
264 }
265#ifdef DEBUGOUTPUT
266 fprintf(stderrstderr, "j: %d tr:%.15g\n",j,tr);
267#endif
268 }
269 else
270 break;
271 }
272 if (increased==1) {
273 if ( (fabs((t-tr)/tr) > reltol) || (fabs(t-tr) > fabs(t-reftime->data[j-1]))) {
274 j = j- 1;
275 tr = reftime->data[j];
276 }
277 }
278#ifdef DEBUGOUTPUT
279 fprintf(stderrstderr, "i: %d t: %.15g d:%.15g j: %d tr:%.15g\n",i,t,d,j,tr);
280#endif
281 /* events, in case of an event compare only the left and right values of the absolute event time range,
282 * this means ta_left = min(t_left,tr_left) and
283 * ta_right = max(t_right,ta_right) */
284 if(i<time->n) {
285#ifdef DEBUGOUTPUT
286 fprintf(stderrstderr, "check event: %.15g - %.15g = %.15g\n",t,time->data[i+1],fabs(t-time->data[i+1]));
287#endif
288 /* an event */
289 if (almostEqualWithDefaultTolerance(t,time->data[i+1])) {
290#ifdef DEBUGOUTPUT
291 fprintf(stderrstderr, "event: %.15g %d %.15g\n",t,i,d);
292#endif
293 /* left value */
294 d_left = d;
295#ifdef DEBUGOUTPUT
296 fprintf(stderrstderr, "left value: %.15g %d %.15g\n",t,i,d_left);
297#endif
298 /* right value */
299 if (i+1<data->n) {
300 while (almostEqualWithDefaultTolerance(t,time->data[i+1])) {
301 i +=1;
302 if (i+1>=data->n) break;
303 }
304 }
305 t = time->data[i];
306 d_right = data->data[i];
307#ifdef DEBUGOUTPUT
308 fprintf(stderrstderr, "right value: %.15g %d %.15g\n",t,i,d_right);
309#endif
310 /* search event in reference forwards */
311 refevent = 0;
312 t_event = t + t*reltol*0.1;
313 /* do not exceed next time step */
314 if (i+1<=data->n) {
315 t_event = (t_event > time->data[i])?time->data[i]:t_event;
316 }else{
317 t_event = (t_event > time->data[i+1])?time->data[i+1]:t_event;
318 }
319 j_event = j;
320 while(tr < t_event) {
321 if (j+1<reftime->n) {
322 if (almostEqualWithDefaultTolerance(tr,reftime->data[j+1])) {
323 dr_left = refdata->data[j];
324#ifdef DEBUGOUTPUT
325 fprintf(stderrstderr, "ref left value: %.15g %d %.15g\n",tr,j,dr_left);
326#endif
327 refevent = 1;
328
329 do {
330 j +=1;
331 if (j+1>=reftime->n) break;
332 } while (almostEqualWithDefaultTolerance(tr,reftime->data[j+1]));
333 }
334 }
335 if (refevent == 0) {
336 j += 1;
337 if (j >= reftime->n)
338 break;
339 tr = reftime->data[j];
340 }
341 else {
342 tr = reftime->data[j];
343 break;
344 }
345 }
346 if (refevent==1) {
347 tr = reftime->data[j];
348 dr_right = refdata->data[j];
349#ifdef DEBUGOUTPUT
350 fprintf(stderrstderr, "ref right value: %.15g %d %.15g\n",tr,j,dr_right);
351#endif
352
353 err = fabs(d_left-dr_left);
354#ifdef DEBUGOUTPUT
355 fprintf(stderrstderr, "delta:%.15g reltol:%.15g\n",err,average);
356#endif
357 if ( err < average){
358 err = fabs(d_right-dr_right);
359#ifdef DEBUGOUTPUT
360 fprintf(stderrstderr, "delta:%.15g reltol:%.15g\n",err,average);
361#endif
362 if ( err < average ) {
363 continue;
364 }
365 }
366 }
367 else {
368 /* search event in reference backwards */
369 j = j_event;
370 tr = reftime->data[j];
371 refevent = 0;
372 t_event = t - t*reltol*0.1;
373 while(tr > t_event) {
374 if (j-1>0) {
375 if (almostEqualWithDefaultTolerance(tr,reftime->data[j-1])) {
376 dr_right = refdata->data[j];
377#ifdef DEBUGOUTPUT
378 fprintf(stderrstderr, "ref right value: %.15g %d %.15g\n",tr,j,dr_right);
379#endif
380 refevent = 1;
381
382 do {
383 j -=1;
384 if (j-1<=0) break;
385 } while (almostEqualWithDefaultTolerance(tr,reftime->data[j-1]));
386 }
387 }
388 if (refevent == 0) {
389 j -= 1;
390 if (j == 0)
391 break;
392 tr = reftime->data[j];
393 }
394 else {
395 tr = reftime->data[j];
396 break;
397 }
398 }
399 if (refevent==1) {
400 tr = reftime->data[j];
401 dr_left = refdata->data[j];
402#ifdef DEBUGOUTPUT
403 fprintf(stderrstderr, "ref left value: %.15g %d %.15g\n",tr,j,dr_left);
404#endif
405 err = fabs(d_left-dr_left);
406#ifdef DEBUGOUTPUT
407 fprintf(stderrstderr, "delta:%.15g reltol:%.15g\n",err,average);
408#endif
409 if ( err < average){
410 err = fabs(d_right-dr_right);
411#ifdef DEBUGOUTPUT
412 fprintf(stderrstderr, "delta:%.15g reltol:%.15g\n",err,average);
413#endif
414 if ( err < average){
415 j = j_event;
416 tr = reftime->data[j];
417 continue;
418 }
419 }
420 }
421 j = j_event;
422 tr = reftime->data[j];
423 }
424 }
425 }
426
427 interpolate = 0;
428#ifdef DEBUGOUTPUT
429 fprintf(stderrstderr, "interpolate? %d %.15g:%.15g %.15g:%.15g\n",i,t,tr,fabs((t-tr)/tr),abstol);
430#endif
431 if (fabs(t-tr) > 0.00001) {
432 interpolate = 1;
433 }
434
435 dr = refdata->data[j];
436 if (interpolate==1){
437#ifdef DEBUGOUTPUT
438 fprintf(stderrstderr, "interpolate %.15g:%.15g %.15g:%.15g %d",t,d,tr,dr,j);
439#endif
440 unsigned int jj = j;
441 /* look for interpolation partner */
442 if (tr > t) {
443 if (j-1 > 0) {
444 jj = j-1;
445 increased = 0;
446 if (reftime->data[jj] == tr){
447 increased = 1;
448 do {
449 jj -= 1;
450 if (jj<=0) break;
451 } while (reftime->data[jj] == tr);
452 }
453 }
454#ifdef DEBUGOUTPUT
455 fprintf(stderrstderr, "-> %d %.15g %.15g\n",jj,reftime->data[jj],refdata->data[jj]);
456#endif
457 if (reftime->data[jj] != tr){
458 dr = refdata->data[jj] + ((dr-refdata->data[jj])/(tr-reftime->data[jj]))*(t-reftime->data[jj]);
459 }
460#ifdef DEBUGOUTPUT
461 fprintf(stderrstderr, "-> dr:%.15g\n",dr);
462#endif
463 }
464 else {
465 if (j+1<reftime->n) {
466 jj = j+1;
467 increased = 0;
468 if (reftime->data[jj] == tr){
469 increased = 1;
470 do {
471 jj += 1;
472 if (jj>=reftime->n) break;
473 } while (reftime->data[jj] == tr);
474 }
475 }
476#ifdef DEBUGOUTPUT
477 fprintf(stderrstderr, "-> %d %.15g %.15g\n",jj,reftime->data[jj],tr);
478#endif
479 if (reftime->data[jj] != tr){
480 dr = dr + ((refdata->data[jj] - dr)/(reftime->data[jj] - tr))*(t-tr);
481 }
482#ifdef DEBUGOUTPUT
483 fprintf(stderrstderr, "-> dr:%.15g\n",dr);
484#endif
485 }
486 }
487#ifdef DEBUGOUTPUT
488 fprintf(stderrstderr, "j: %d tr: %.15g dr:%.15g t:%.15g d:%.15g\n",j,tr,dr,t,d);
489#endif
490 err = fabs(d-dr);
491#ifdef DEBUGOUTPUT
492 fprintf(stderrstderr, "delta:%.15g reltol:%.15g\n",err,average);
493#endif
494 if (fout) {
495 fprintf(fout, "%.15g,%.15g,%.15g,%.15g,%.15g,%.15g\n",tr,dr,d,err,almostEqualWithDefaultTolerance(d,0) ? err/average : fabs(err/d),average);
496 }
497 if ( err > average){
498 if (j+1<reftime->n) {
499 if (reftime->data[j+1] == tr) {
500 dr = refdata->data[j+1];
501 err = fabs(d-dr);
502 }
503 }
504
505 if (err < average){
506 continue;
507 }
508
509 isdifferent = 1;
510 if (isResultCmp) { /* If we produce the full diff, this data has already been output */
511 if (ddf->n >= ddf->n_max) {
512 DiffData *newData;
513 ddf->n_max = ddf->n_max ? ddf->n_max*2 : 1024;
514 newData = (DiffData*) realloc(ddf->data, sizeof(DiffData)*(ddf->n_max));
515 if (!newData) continue; /* realloc failed... pretty bad, but let's continue */
516 ddf->data = newData;
517 }
518 ddf->data[ddf->n].name = varname;
519 ddf->data[ddf->n].data = d;
520 ddf->data[ddf->n].dataref = dr;
521 ddf->data[ddf->n].time = t;
522 ddf->data[ddf->n].timeref = tr;
523 ddf->data[ddf->n].interpolate = interpolate?'1':'0';
524 ddf->n +=1;
525 }
526 }
527 }
528 if (isdifferent) {
529 cmpdiffvars[vardiffindx] = varname;
530 vardiffindx++;
531 if (!isResultCmp) {
532 *diffLst = mmc_mk_cons(mmc_mk_scon(varname),*diffLst);
533 }
534 }
535 if (fout) {
536 fclose(fout);
537 }
538 if (!isdifferent && 0==keepEqualResults && 0==isResultCmp) {
539 SystemImpl__removeFile(fname);
540 }
541 if (fname) {
542 free(fname);
543 }
544 return vardiffindx;
545}
546
547static int writeLogFile(const char *filename,DiffDataField *ddf,const char *f,const char *reff,double reltol,double abstol)
548{
549 FILE* fout;
550 unsigned int i;
551 /* fprintf(stderr, "writeLogFile: %s\n",filename); */
552 fout = fopen(filename, "w");
553 if (!fout)
554 return -1;
555
556 fprintf(fout, "\"Generated by OpenModelica\";;;;;\n");
557 fprintf(fout, "\"Compared Files\";;;\"absolute tolerance\";%.15g;relative tolerance;%.15g\n",abstol,reltol);
558 fprintf(fout, "\"%s\";;;;;;\n",f);
559 fprintf(fout, "\"%s\";;;;;;\n",reff);
560 fprintf(fout, "\"Name\";\"Time\";\"DataPoint\";\"RefTime\";\"RefDataPoint\";\"absolute error\";\"relative error\";interpolate;\n");
561 for (i=0;i<ddf->n;i++){
562 fprintf(fout, "%s;%.15g;%.15g;%.15g;%.15g;%.15g;%.15g;%c;\n",ddf->data[i].name,ddf->data[i].time,ddf->data[i].data,ddf->data[i].timeref,ddf->data[i].dataref,
563 fabs(ddf->data[i].data-ddf->data[i].dataref),fabs((ddf->data[i].data-ddf->data[i].dataref)/ddf->data[i].dataref),ddf->data[i].interpolate);
564 }
565 fclose(fout);
566 /* fprintf(stderr, "writeLogFile: %s finished\n",filename); */
567 return 0;
568}
569
570static const char* getTimeVarName(void *vars) {
571 const char *res = "time";
572 while (MMC_NILHDR(((0) << 10) + (((0) & 255) << 2)) != MMC_GETHDR(vars)(*(mmc_uint_t*)((void*)((char*)(vars) - 3)))) {
573 char *var = MMC_STRINGDATA(MMC_CAR(vars))(((struct mmc_string*)((void*)((char*)((*(void**)(((void*)((void
**)(((void*)((char*)(vars) - 3))) + (1)))))) - 3)))->data)
;
574 if (0==strcmp("time",var)) break;
575 if (0==strcmp("Time",var)) {
576 res="Time";
577 break;
578 }
579 vars = MMC_CDR(vars)(*(void**)(((void*)((void**)(((void*)((char*)(vars) - 3))) + (
2)))))
;
580 }
581 return res;
582}
583
584#include "SimulationResultsCmpTubes.c"
585
586/* Common, huge function, for both result comparison and result diff */
587void* SimulationResultsCmp_compareResults(int isResultCmp, int runningTestsuite, const char *filename, const char *reffilename, const char *resultfilename, double reltol, double abstol, double reltolDiffMaxMin, double rangeDelta, void *vars, int keepEqualResults, int *success, int isHtml, char **htmlOut)
588{
589 char **cmpvars=NULL((void*)0);
590 char **cmpdiffvars=NULL((void*)0);
591 unsigned int vardiffindx=0;
592 unsigned int ncmpvars = 0;
593 unsigned int ngetfailedvars = 0;
594 void *allvars,*allvarsref,*res;
595 unsigned int i,size,size_ref,len,j,k;
596 char *var,*var1,*var2;
597 DataField time,timeref,data,dataref;
598 DiffDataField ddf;
599 const char *msg[2] = {"",""};
600 const char *timeVarName, *timeVarNameRef;
601 int suggestReadAll=0;
602 ddf.data=NULL((void*)0);
603 ddf.n=0;
604 ddf.n_max=0;
605 len = 1;
606 int offset, offsetRef;
607
608 /* open files */
609 /* fprintf(stderr, "Open File %s\n", filename); */
610 if (UNKNOWN_PLOT == SimulationResultsImpl__openFile(filename,&simresglob_c)) {
2
Taking false branch
611 char *str = (char*) omc_alloc_interface.malloc(25+strlen(filename));
612 void *res = NULL((void*)0);
613 *str = 0;
614 strcat(strcat(str,"Error opening file: "), filename);
615 res = mmc_mk_scon(str);
616 GC_free(str);
617 return mmc_mk_cons(res,mmc_mk_nil()((void*)((char*)(&(mmc_nil).header) + 3)));
618 }
619 /* fprintf(stderr, "Open File %s\n", reffilename); */
620 if (UNKNOWN_PLOT == SimulationResultsImpl__openFile(reffilename,&simresglob_ref)) {
3
Taking false branch
621 char *str = (char*) omc_alloc_interface.malloc(35+strlen(reffilename));
622 void *res = NULL((void*)0);
623 *str = 0;
624 strcat(strcat(str,"Error opening reference file: "), reffilename);
625 res = mmc_mk_scon(str);
626 GC_free(str);
627 return mmc_mk_cons(res,mmc_mk_nil()((void*)((char*)(&(mmc_nil).header) + 3)));
628 }
629
630 size = SimulationResultsImpl__readSimulationResultSize(filename,&simresglob_c);
631 /* fprintf(stderr, "Read size of File %s size= %d\n", filename,size); */
632 size_ref = SimulationResultsImpl__readSimulationResultSize(reffilename,&simresglob_ref);
633 /* fprintf(stderr, "Read size of File %s size= %d\n", reffilename,size_ref); */
634
635 /* get vars to compare */
636 cmpvars = getVars(vars,&ncmpvars);
637 /* if no var compare all vars */
638 allvars = SimulationResultsImpl__readVarsFilterAliases(filename,&simresglob_c);
639 allvarsref = SimulationResultsImpl__readVarsFilterAliases(reffilename,&simresglob_ref);
640 if (ncmpvars==0) {
4
Taking false branch
641 suggestReadAll = 1;
642 cmpvars = getVars(allvarsref,&ncmpvars);
643 if (ncmpvars==0) return mmc_mk_cons(mmc_mk_scon("Error Get Vars!"),mmc_mk_nil()((void*)((char*)(&(mmc_nil).header) + 3)));
644 }
645#ifdef DEBUGOUTPUT
646 fprintf(stderrstderr, "Compare Vars:\n");
647 for(i=0;i<ncmpvars;i++)
648 fprintf(stderrstderr, "Var: %s\n", cmpvars[i]);
649#endif
650 /* get time */
651 /* fprintf(stderr, "get time\n"); */
652 timeVarName = getTimeVarName(allvars);
653 timeVarNameRef = getTimeVarName(allvarsref);
654 time = getData(timeVarName,filename,size,suggestReadAll,&simresglob_c,runningTestsuite);
655 if (time.n==0) {
5
Assuming the condition is false
6
Taking false branch
656 return mmc_mk_cons(mmc_mk_scon("Error get time!"),mmc_mk_nil()((void*)((char*)(&(mmc_nil).header) + 3)));
657 }
658 /* fprintf(stderr, "get reftime\n"); */
659 timeref = getData(timeVarNameRef,reffilename,size_ref,suggestReadAll,&simresglob_ref,runningTestsuite);
660 if (timeref.n==0) {
7
Assuming the condition is false
8
Taking false branch
661 return mmc_mk_cons(mmc_mk_scon("Error get ref time!"),mmc_mk_nil()((void*)((char*)(&(mmc_nil).header) + 3)));
662 }
663 cmpdiffvars = (char**)omc_alloc_interface.malloc(sizeof(char*)*(ncmpvars));
664 /* check if time is larger or less reftime */
665 res = mmc_mk_nil()((void*)((char*)(&(mmc_nil).header) + 3));
666 if (fabs(time.data[time.n-1]-timeref.data[timeref.n-1]) > reltol*fabs(timeref.data[timeref.n-1])) {
9
Taking false branch
667 char buf[WARNINGBUFFSIZE4096];
668#ifdef DEBUGOUTPUT
669 fprintf(stderrstderr, "max time value=%.15g ref max time value: %.15g\n",time.data[time.n-1],timeref.data[timeref.n-1]);
670#endif
671 snprintf(buf,WARNINGBUFFSIZE4096,"Resultfile and Reference have different end time points!\n"
672 "Reffile[%d]=%f\n"
673 "File[%d]=%f\n",timeref.n,timeref.data[timeref.n-1],time.n,time.data[time.n-1]);
674 c_add_message(NULL((void*)0),-1, ErrorType_scripting, ErrorLevel_warning, buf, NULL((void*)0), 0);
675 }
676 /* calculate offsets */
677 for(offset=0; offset<time.n-1 && time.data[offset] == time.data[offset+1]; ++offset);
10
Assuming the condition is false
678 for(offsetRef=0; offsetRef<timeref.n-1 && timeref.data[offsetRef] == timeref.data[offsetRef+1]; ++offsetRef);
11
Assuming the condition is false
679 var1=NULL((void*)0);
680 var2=NULL((void*)0);
681 /* compare vars */
682 /* fprintf(stderr, "compare vars\n"); */
683 for (i=0;i<ncmpvars;i++) {
12
Loop condition is true. Entering loop body
684 var = cmpvars[i];
685 len = strlen(var);
686 if (var1) {
13
Taking false branch
687 free(var1);
688 var1 = NULL((void*)0);
689 }
690 var1 = (char*) omc_alloc_interface.malloc(len+10);
691 k = 0;
692 for (j=0;j<len;j++) {
14
Assuming 'j' is >= 'len'
15
Loop condition is false. Execution continues on line 698
693 if (var[j] !='\"' ) {
694 var1[k] = var[j];
695 k +=1;
696 }
697 }
698 var1[k] = 0;
699 /* fprintf(stderr, "compare var: %s\n",var); */
700 /* check if in ref_file */
701 dataref = getData(var1,reffilename,size_ref,suggestReadAll,&simresglob_ref,runningTestsuite);
702 if (dataref.n==0) {
16
Assuming the condition is false
17
Taking false branch
703 if (dataref.data) {
704 free(dataref.data);
705 }
706 if (var1) {
707 GC_free(var1);
708 var1 = NULL((void*)0);
709 }
710 msg[0] = runningTestsuite ? SystemImpl__basename(reffilename) : reffilename;
711 msg[1] = var;
712 c_add_message(NULL((void*)0),-1, ErrorType_scripting, ErrorLevel_warning, gettext("Get data of variable %s from file %s failed!\n")dcgettext (((void*)0), "Get data of variable %s from file %s failed!\n"
, 5)
, msg, 2);
713 ngetfailedvars++;
714 continue;
715 }
716 /* check if in file */
717 data = getData(var1,filename,size,suggestReadAll,&simresglob_c,runningTestsuite);
718 if (data.n==0) {
18
Assuming the condition is false
19
Taking false branch
719 if (data.data) {
720 free(data.data);
721 }
722 if (var1) {
723 GC_free(var1);
724 var1 = NULL((void*)0);
725 }
726 msg[0] = runningTestsuite ? SystemImpl__basename(filename) : filename;
727 msg[1] = var;
728 c_add_message(NULL((void*)0),-1, ErrorType_scripting, ErrorLevel_warning, gettext("Get data of variable %s from file %s failed!\n")dcgettext (((void*)0), "Get data of variable %s from file %s failed!\n"
, 5)
, msg, 2);
729 ngetfailedvars++;
730 continue;
731 }
732 /* adjust initial data points */
733 for(j=offset; j>0; j--)
20
Loop condition is false. Execution continues on line 735
734 data.data[j-1] = data.data[j];
735 for(j=offsetRef; j>0; j--)
21
Loop condition is false. Execution continues on line 738
736 dataref.data[j-1] = dataref.data[j];
737 /* compare */
738 if (isHtml) {
22
Taking true branch
739 vardiffindx = cmpDataTubes(isResultCmp,var,&time,&timeref,&data,&dataref,reltol,rangeDelta,reltolDiffMaxMin,&ddf,cmpdiffvars,vardiffindx,keepEqualResults,&res,resultfilename,1,htmlOut);
23
Calling 'cmpDataTubes'
740 } else if (isResultCmp) {
741 vardiffindx = cmpData(isResultCmp,var,&time,&timeref,&data,&dataref,reltol,abstol,&ddf,cmpdiffvars,vardiffindx,keepEqualResults,&res,resultfilename);
742 } else {
743 vardiffindx = cmpDataTubes(isResultCmp,var,&time,&timeref,&data,&dataref,reltol,rangeDelta,reltolDiffMaxMin,&ddf,cmpdiffvars,vardiffindx,keepEqualResults,&res,resultfilename,0,0);
744 }
745 /* free */
746 if (dataref.data) {
747 free(dataref.data);
748 }
749 if (data.data) {
750 free(data.data);
751 }
752 if (var1) {
753 GC_free(var1);
754 var1 = NULL((void*)0);
755 }
756 }
757
758 if (isResultCmp) {
759 if (writeLogFile(resultfilename,&ddf,filename,reffilename,reltol,abstol)) {
760 c_add_message(NULL((void*)0),-1, ErrorType_scripting, ErrorLevel_warning, gettext("Cannot write to the difference (.csv) file!\n")dcgettext (((void*)0), "Cannot write to the difference (.csv) file!\n"
, 5)
, msg, 0);
761 }
762
763 if ((ddf.n > 0) || (ngetfailedvars > 0) || vardiffindx > 0){
764 /* fprintf(stderr, "diff: %d\n",ddf.n); */
765 /* for (i=0;i<vardiffindx;i++)
766 fprintf(stderr, "diffVar: %s\n",cmpdiffvars[i]); */
767 for (i=0;i<vardiffindx;i++){
768 res = (void*)mmc_mk_cons(mmc_mk_scon(cmpdiffvars[i]),res);
769 }
770 res = mmc_mk_cons(mmc_mk_scon("Files not Equal!"),res);
771 c_add_message(NULL((void*)0),-1, ErrorType_scripting, ErrorLevel_warning, gettext("Files not Equal\n")dcgettext (((void*)0), "Files not Equal\n", 5), msg, 0);
772 } else {
773 res = mmc_mk_cons(mmc_mk_scon("Files Equal!"),res);
774 }
775 } else {
776 if (success) {
777 *success = ((ddf.n == 0) && (vardiffindx == 0));
778 }
779 }
780
781 // if (var1) free(var1);
782 // if (var2) free(var2);
783 if (ddf.data) free(ddf.data);
784 if (cmpvars) GC_free(cmpvars);
785 if (time.data) free(time.data);
786 if (timeref.data) free(timeref.data);
787 if (cmpdiffvars) GC_free(cmpdiffvars);
788 /* close files */
789 SimulationResultsImpl__close(&simresglob_c);
790 SimulationResultsImpl__close(&simresglob_ref);
791
792 return res;
793}
794
795
796/* Common, huge function, for both result comparison and result diff */
797double SimulationResultsCmp_deltaResults(const char *filename, const char *reffilename, const char *methodname, void *vars)
798{
799 double res = 0;
800 unsigned int i,size,size_ref, len, k, j;
801 unsigned int ncmpvars = 0;
802 void *allvars,*allvarsref;
803 char *var,*var1;
804 char **cmpvars=NULL((void*)0);
805 int suggestReadAll=0;
806 DataField time,timeref,data,dataref;
807 const char *timeVarName, *timeVarNameRef;
808 int offset, offsetRef;
809 const char *msg[2] = {"",""};
810
811 /* choose error method */
812 ErrorMethod errMethod;
813
814 if (0 == strcmp(methodname, "1norm")){
815 errMethod = NORM1;
816 }
817 else if (0 == strcmp(methodname, "2norm")){
818 errMethod = NORM2;
819 }
820 else if (0 == strcmp(methodname, "maxerr")){
821 errMethod = MAX_ERR;
822 }
823 else {
824 msg[0] = methodname;
825 c_add_message(NULL((void*)0),-1, ErrorType_scripting, ErrorLevel_warning, gettext("Unknown method string: %s. 1-Norm is chosen.")dcgettext (((void*)0), "Unknown method string: %s. 1-Norm is chosen."
, 5)
, msg, 1);
826 errMethod = NORM1;
827 }
828
829 /* open files */
830 /* fprintf(stderr, "Open File %s\n", filename); */
831
832 SimulationResultsImpl__close(&simresglob_c);
833 SimulationResultsImpl__close(&simresglob_ref);
834
835 if (UNKNOWN_PLOT == SimulationResultsImpl__openFile(filename,&simresglob_c)) {
836 msg[0] = filename;
837 c_add_message(NULL((void*)0),-1, ErrorType_scripting, ErrorLevel_error, gettext("Error opening file: %s.")dcgettext (((void*)0), "Error opening file: %s.", 5), msg, 1);
838 return -1;
839 }
840 /* fprintf(stderr, "Open File %s\n", reffilename); */
841 if (UNKNOWN_PLOT == SimulationResultsImpl__openFile(reffilename,&simresglob_ref)) {
842 msg[0] = filename;
843 c_add_message(NULL((void*)0),-1, ErrorType_scripting, ErrorLevel_error, gettext("Error opening reference file: %s.")dcgettext (((void*)0), "Error opening reference file: %s.", 5
)
, msg, 1);
844 return -1;
845 }
846
847 size = SimulationResultsImpl__readSimulationResultSize(filename,&simresglob_c);
848 /*fprintf(stderr, "Read size of File %s size= %d\n", filename,size);*/
849 size_ref = SimulationResultsImpl__readSimulationResultSize(reffilename,&simresglob_ref);
850 /*fprintf(stderr, "Read size of File %s size= %d\n", reffilename,size_ref);*/
851
852 /* get vars to compare */
853 cmpvars = getVars(vars,&ncmpvars);
854 /* if no var compare all vars */
855 allvars = SimulationResultsImpl__readVarsFilterAliases(filename,&simresglob_c);
856 allvarsref = SimulationResultsImpl__readVarsFilterAliases(reffilename,&simresglob_ref);
857 if (ncmpvars==0) {
858 suggestReadAll = 1;
859 cmpvars = getVars(allvarsref,&ncmpvars);
860 if (ncmpvars==0){
861 c_add_message(NULL((void*)0),-1, ErrorType_scripting, ErrorLevel_error, gettext("Error Getting Vars.")dcgettext (((void*)0), "Error Getting Vars.", 5), msg, 1);
862 return -1;
863 }
864 }
865 #ifdef DEBUGOUTPUT
866 fprintf(stderrstderr, "Compare Vars:\n");
867 for(i=0;i<ncmpvars;i++)
868 fprintf(stderrstderr, "Var: %s\n", cmpvars[i]);
869 #endif
870
871 /* get time */
872 /*fprintf(stderr, "get time\n");*/
873 timeVarName = getTimeVarName(allvars);
874 timeVarNameRef = getTimeVarName(allvarsref);
875 time = getData(timeVarName,filename,size,suggestReadAll,&simresglob_c,0);
876 if (time.n==0) {
877 c_add_message(NULL((void*)0),-1, ErrorType_scripting, ErrorLevel_error, gettext("Error get time!")dcgettext (((void*)0), "Error get time!", 5), msg, 0);
878 return -1;
879 }
880 /*fprintf(stderr, "get reftime\n");*/
881 timeref = getData(timeVarNameRef,reffilename,size_ref,suggestReadAll,&simresglob_ref,0);
882 if (timeref.n==0) {
883 c_add_message(NULL((void*)0),-1, ErrorType_scripting, ErrorLevel_error, gettext("Error get reference time!")dcgettext (((void*)0), "Error get reference time!", 5), msg, 0);
884 return -1;
885 }
886
887 /* check if time is larger or less reftime */
888 if (time.data[0] > timeref.data[0]) {
889 c_add_message(NULL((void*)0),-1, ErrorType_scripting, ErrorLevel_error, gettext("The result file starts before the reference file.")dcgettext (((void*)0), "The result file starts before the reference file."
, 5)
, msg, 0);
890 return -1;
891 }
892 if(time.data[time.n-1] < timeref.data[timeref.n-1]){
893 c_add_message(NULL((void*)0),-1, ErrorType_scripting, ErrorLevel_error, gettext("The result file ends before the reference file.")dcgettext (((void*)0), "The result file ends before the reference file."
, 5)
, msg, 0);
894 return -1;
895 }
896
897 /* calculate offsets */
898 for(offset=0; offset<time.n-1 && time.data[offset] == time.data[offset+1]; ++offset);
899 for(offsetRef=0; offsetRef<timeref.n-1 && timeref.data[offsetRef] == timeref.data[offsetRef+1]; ++offsetRef);
900 var1=NULL((void*)0);
901 /* compare vars */
902 /* fprintf(stderr, "compare vars\n"); */
903 for (i=0;i<ncmpvars;i++) {
904 var = cmpvars[i];
905 len = strlen(var);
906 if (var1) {
907 free(var1);
908 var1 = NULL((void*)0);
909 }
910 var1 = (char*) omc_alloc_interface.malloc(len+10);
911 k = 0;
912 for (j=0;j<len;j++) {
913 if (var[j] !='\"' ) {
914 var1[k] = var[j];
915 k +=1;
916 }
917 }
918 var1[k] = 0;
919 /* fprintf(stderr, "compare var: %s\n",var); */
920 /* check if in ref_file */
921 dataref = getData(var1,reffilename,size_ref,suggestReadAll,&simresglob_ref,0);
922 if (dataref.n==0) {
923 if (dataref.data) {
924 free(dataref.data);
925 }
926 if (var1) {
927 GC_free(var1);
928 var1 = NULL((void*)0);
929 }
930 continue;
931 }
932 /* check if in file */
933 data = getData(var1,filename,size,suggestReadAll,&simresglob_c,0);
934 if (data.n==0) {
935 if (data.data) {
936 free(data.data);
937 }
938 if (var1) {
939 GC_free(var1);
940 var1 = NULL((void*)0);
941 }
942 continue;
943 }
944 /* adjust initial data points */
945 for(j=offset; j>0; j--)
946 data.data[j-1] = data.data[j];
947 for(j=offsetRef; j>0; j--)
948 dataref.data[j-1] = dataref.data[j];
949
950 /* calulate delta */
951 res += deltaData(errMethod,&time,&timeref,&data,&dataref);
952
953 /* free */
954 if (dataref.data) {
955 free(dataref.data);
956 }
957 if (data.data) {
958 free(data.data);
959 }
960 if (var1) {
961 GC_free(var1);
962 var1 = NULL((void*)0);
963 }
964 }
965 if (cmpvars) GC_free(cmpvars);
966 if (time.data) free(time.data);
967 if (timeref.data) free(timeref.data);
968
969 /* close files */
970 SimulationResultsImpl__close(&simresglob_c);
971 SimulationResultsImpl__close(&simresglob_ref);
972 return res;
973}
974

./SimulationResultsCmpTubes.c

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;
28
'x0' declared without an initial value
377 size_t n;
378
379 if (0 == nsource) {
29
Taking false branch
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++) {
30
Loop condition is true. Entering loop body
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
31
Taking true branch
391 interpolatedValues[i] = linearInterpolation(x,x0,x1,y0,y1,xabstol);
32
2nd function call argument is an uninitialized value
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);
24
Assuming 'withTubes' is not equal to 0
25
'?' condition is true
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);
26
'?' condition is true
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);
27
Calling 'calibrateValues'
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);
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}