@@ -334,18 +334,6 @@ static privates* calculateTubes(double *x, double *y, size_t length, double r)
334334 return priv ;
335335}
336336
337- /* Count the number of errors */
338- static int validate (int n , double * low , double * high , double * val )
339- {
340- int iErrors = 0 , i ;
341- for (i = 0 ; i < n ; i ++ ) {
342- if (val [i ] < low [i ] || val [i ] > high [i ]) {
343- iErrors ++ ; //Count as error if current value is bigger than value of high tube and smaller than value of lowtube
344- }
345- }
346- return iErrors ;
347- }
348-
349337static inline double linearInterpolation (double x , double x0 , double x1 , double y0 , double y1 , double xabstol )
350338{
351339 if (almostEqualRelativeAndAbs (x0 ,x ,0 ,xabstol )) { //prevent NaN -> division by zero
@@ -506,7 +494,7 @@ static addTargetEventTimesRes mergeTimelines(addTargetEventTimesRes ref, addTarg
506494 return res ;
507495}
508496
509- static addTargetEventTimesRes removeUneventfulPoints (addTargetEventTimesRes in , int removeNonEvents , double reltol , double xabstol )
497+ static addTargetEventTimesRes removeUneventfulPoints (addTargetEventTimesRes in , double reltol , double xabstol )
510498{
511499 int i ,iter ;
512500 addTargetEventTimesRes res ;
@@ -532,7 +520,7 @@ static addTargetEventTimesRes removeUneventfulPoints(addTargetEventTimesRes in,
532520 double x1 = in .time [i + 1 ];
533521 double y1 = in .values [i + 1 ];
534522 int isLocalMinOrMax = (x > x1 && x > x0 ) || (x < x1 && x < x0 );
535- if (!isLocalMinOrMax && almostEqualRelativeAndAbs (y ,linearInterpolation (x ,x0 ,x1 ,y0 ,y1 ,xabstol ),reltol ,0 ) && !isEvent && removeNonEvents ) {
523+ if (!isLocalMinOrMax && almostEqualRelativeAndAbs (y ,linearInterpolation (x ,x0 ,x1 ,y0 ,y1 ,xabstol ),reltol ,0 ) && !isEvent ) {
536524 /* The point can be reconstructed using linear interpolation and is not a local minimum or maximum */
537525 res .time [res .size ] = x1 ;
538526 res .values [res .size ] = y1 ;
@@ -553,26 +541,15 @@ static addTargetEventTimesRes removeUneventfulPoints(addTargetEventTimesRes in,
553541 res .time [res .size ] = in .time [in .size - 1 ];
554542 res .size ++ ;
555543 }
556-
557544 if (iter ) {
558- /* swap buffers for in and res */
559- double * tmp ;
560- tmp = in .time ;
561- in .time = res .time ;
562- res .time = tmp ;
563- tmp = in .values ;
564- in .values = res .values ;
565- res .values = tmp ;
566- /* set size of in, res will be overwritten anyways*/
567- in .size = res .size ;
568- /* set new tolerance */
569- reltol /= 10.0 ;
545+ /* 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) */
546+ in = res ;
547+ reltol = 0 ;
570548 } else {
571549 /* no more recursion */
572550 break ;
573551 }
574552 } while (1 );
575-
576553 return res ;
577554}
578555
@@ -602,9 +579,47 @@ static void assertMonotonic(addTargetEventTimesRes series)
602579 }
603580}
604581
582+ /* Return NULL if there were no errors */
583+ static double * validate (int n , addTargetEventTimesRes ref , double * low , double * high , double * calibrated_values , double reltol , double abstol , double xabstol )
584+ {
585+ double * error = GC_malloc_atomic (n * sizeof (double ));
586+ int isdifferent = 0 ;
587+ int i ,lastStepError = 1 ;
588+ for (i = 0 ; i < n ; i ++ ) {
589+ int thisStepError = 0 ;
590+ 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 ));
591+ if (isEvent ) {
592+ double refv = ref .values [i ];
593+ double val = calibrated_values [i ];
594+ double tol = fmax (abstol * 10 ,fmax (fabs (refv ),fabs (val ))* reltol * 10 );
595+ /* If there was no error in the last step before the event, give a nice tight curve around both values */
596+ high [i ] = (lastStepError ? refv : fmax (refv ,val )) + tol ;
597+ low [i ] = (lastStepError ? refv : fmin (refv ,val )) - tol ;
598+ error [i ] = NAN ;
599+ } else {
600+ error [i ] = 0 ;
601+ thisStepError = lastStepError ;
602+ if (calibrated_values [i ] < low [i ]) {
603+ error [i ] = low [i ]- calibrated_values [i ];
604+ isdifferent ++ ;
605+ thisStepError = 1 ;
606+ } else if (calibrated_values [i ] > high [i ]) {
607+ error [i ] = calibrated_values [i ]- high [i ];
608+ isdifferent ++ ;
609+ thisStepError = 1 ;
610+ }
611+ }
612+ lastStepError = thisStepError ;
613+ }
614+ if (isdifferent ) {
615+ return error ;
616+ }
617+ return NULL ;
618+ }
619+
605620static 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 )
606621{
607- int i , isdifferent = 0 ;
622+ int i ;
608623 FILE * fout = NULL ;
609624 char * fname = NULL ;
610625 char * html ;
@@ -621,8 +636,8 @@ static unsigned int cmpDataTubes(int isResultCmp, char* varname, DataField *time
621636 actual = actualoriginal ;
622637 /* assertMonotonic(ref); */
623638 /* assertMonotonic(actual); */
624- ref = removeUneventfulPoints (ref , 1 , reltol * reltol , xabstol );
625- actual = removeUneventfulPoints (actual , 1 , reltol * reltol , xabstol );
639+ ref = removeUneventfulPoints (ref , reltol * reltol , xabstol );
640+ actual = removeUneventfulPoints (actual , reltol * reltol , xabstol );
626641 /* assertMonotonic(ref); */
627642 /* assertMonotonic(actual); */
628643 privates * priv = calculateTubes (ref .time ,ref .values ,ref .size ,rangeDelta );
@@ -639,6 +654,7 @@ static unsigned int cmpDataTubes(int isResultCmp, char* varname, DataField *time
639654 double abstol = fabs ((priv -> max - priv -> min )* reltolDiffMaxMin );
640655 addRelativeTolerance (high ,ref .values ,n ,reltol ,abstol ,1 );
641656 addRelativeTolerance (low ,ref .values ,n ,reltol ,abstol ,-1 );
657+ double * error = validate (n ,ref ,low ,high ,calibrated_values ,reltol ,abstol ,xabstol );
642658 if (isHtml ) {
643659#if _XOPEN_SOURCE >= 700 || _POSIX_C_SOURCE >= 200809L
644660 fout = open_memstream (& html , & html_size );
@@ -686,7 +702,7 @@ static unsigned int cmpDataTubes(int isResultCmp, char* varname, DataField *time
686702 reltolDiffMaxMin ,
687703 rangeDelta
688704);
689- } else if (!isResultCmp ) {
705+ } else if (!isResultCmp && ( error || keepEqualResults ) ) {
690706 fname = (char * ) GC_malloc_atomic (25 + strlen (prefix ) + strlen (varname ));
691707 sprintf (fname , "%s.%s.csv" , prefix , varname );
692708 fout = fopen (fname ,"w" );
@@ -705,28 +721,10 @@ static unsigned int cmpDataTubes(int isResultCmp, char* varname, DataField *time
705721 int thisStepError = 0 ;
706722 fprintf (fout , "%s%.15g,%.15g," ,lbracket ,ref .time [i ],ref .values [i ]);
707723 if (i < n ) {
708- 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 ));
709- if (isEvent ) {
710- double refv = ref .values [i ];
711- double val = calibrated_values [i ];
712- double tol = fmax (abstol * 10 ,fmax (fabs (refv ),fabs (val ))* reltol * 10 );
713- /* If there was no error in the last step before the event, give a nice tight curve around both values */
714- double high = (lastStepError ? refv : fmax (refv ,val )) + tol ;
715- double low = (lastStepError ? refv : fmin (refv ,val )) - tol ;
716- fprintf (fout , "%.15g,%.15g,%.15g,%s" ,val ,high ,low ,empty );
724+ if (!error || isnan (error [i ])) {
725+ fprintf (fout , "%.15g,%.15g,%.15g,%s" ,calibrated_values [i ],high [i ],low [i ],empty );
717726 } else {
718- double error = 0 ;
719- thisStepError = lastStepError ;
720- if (calibrated_values [i ] < low [i ]) {
721- error = low [i ]- calibrated_values [i ];
722- isdifferent ++ ;
723- thisStepError = 1 ;
724- } else if (calibrated_values [i ] > high [i ]) {
725- error = calibrated_values [i ]- high [i ];
726- isdifferent ++ ;
727- thisStepError = 1 ;
728- }
729- fprintf (fout , "%.15g,%.15g,%.15g,%.15g" ,calibrated_values [i ],high [i ],low [i ],error );
727+ fprintf (fout , "%.15g,%.15g,%.15g,%.15g" ,calibrated_values [i ],high [i ],low [i ],error [i ]);
730728 }
731729 if (ref .time [i ] == actualoriginal .time [j ] && j < actualoriginal .size ) {
732730 fprintf (fout , ",%.15g%s\n" ,actualoriginal .values [j ++ ],rbracket );
@@ -743,10 +741,8 @@ static unsigned int cmpDataTubes(int isResultCmp, char* varname, DataField *time
743741 lastStepError = thisStepError ;
744742 }
745743 fputs (isHtml ? "],\n" : "\n" , fout );
746- } else {
747- isdifferent = validate (n ,low ,high ,calibrated_values );
748744 }
749- if (isdifferent ) {
745+ if (error ) {
750746 cmpdiffvars [vardiffindx ] = varname ;
751747 vardiffindx ++ ;
752748 if (!isResultCmp ) {
@@ -779,8 +775,22 @@ fprintf(fout, "{title: '%s',\n"
779775 free (html );
780776 }
781777 }
782- if (!isdifferent && 0 == keepEqualResults && 0 == isResultCmp ) {
783- SystemImpl__removeFile (fname );
784- }
778+ /* Tell the GC some variables have been free'd */
779+ if (error ) GC_free (error );
780+ if (fname ) GC_free (fname );
781+ GC_free (low );
782+ GC_free (high );
783+ GC_free (priv -> mh );
784+ GC_free (priv -> i0h );
785+ GC_free (priv -> i1h );
786+ GC_free (priv -> ml );
787+ GC_free (priv -> i0l );
788+ GC_free (priv -> i1l );
789+ GC_free (priv -> xHigh );
790+ GC_free (priv -> xLow );
791+ GC_free (priv -> yHigh );
792+ GC_free (priv -> yLow );
793+ GC_free (priv );
794+ GC_free (calibrated_values );
785795 return vardiffindx ;
786796}
0 commit comments