Improve tube calculation

Modelica_ResultCompare/CsvFile.cs

@@ -265,7 +265,69 @@ public Report CompareFiles(Log log, CsvFile csvBase, ref Options options)
             else
                 return CompareFiles(log, csvBase, null, ref options);
         }
+
+        static Curve Trim_curve(Curve curve,
+                       double start_time,
+                       double end_time)
+        {
+            int start_idx = 0;
+            int end_idx = curve.Count;
+
+            /* Find first valid sample */
+            while (start_idx < curve.Count &&
+                   curve.X[start_idx] < start_time)
+            {
+                start_idx++;
+            }
+
+            /* Find last valid sample */
+            while (end_idx > start_idx &&
+                   curve.X[end_idx - 1] > end_time)
+            {
+                end_idx--;
+            }
+
+            int new_size = end_idx - start_idx;
+            double[] TargetValues = new double[new_size];
+            double[] TargetTime = new double[new_size];
+            /* Shift data to beginning */
+            for (int i = 0; i < new_size; ++i)
+            {
+                TargetTime[i] = curve.X[start_idx + i];
+                TargetValues[i] = curve.Y[start_idx + i];
+            }
+            Curve NewCurve = new Curve(curve.Name, TargetTime, TargetValues);
+            return NewCurve; 
+        }
+
+        public double[] GetCommonInterval(Curve reference, Curve compareCurve)
+        {
+            double[] common_interval = new double[2] { 0.0, 0.0 };
 
+            if (compareCurve.Count == 0 || reference.Count == 0)
+            {
+                return common_interval;
+            }
+            // Determine common interval
+            double common_start =
+            (reference.X[0] > compareCurve.X[0])
+            ? reference.X[0]
+            : compareCurve.X[0];
+
+            double common_stop =
+                (reference.X[reference.X.Length - 1] <
+                 compareCurve.X[compareCurve.X.Length - 1])
+                    ? reference.X[reference.X.Length - 1]
+                    : compareCurve.X[compareCurve.X.Length - 1];
+            if (common_start > common_stop)
+            {
+                return common_interval;
+            }
+            common_interval[0] = common_start;
+            common_interval[1] = common_stop;
+            return common_interval;
+        }
+
         public Report CompareFiles(Log log, CsvFile csvBase, string sReportPath, ref Options options)
         {
             int iInvalids = 0;
@@ -276,8 +338,11 @@ public Report CompareFiles(Log log, CsvFile csvBase, string sReportPath, ref Opt
             rep.BaseFile = csvBase.ToString();
             rep.CompareFile = _fileName;
 
+            double[] common_interval =  new double [2];
             Curve reference = new Curve();
             Curve compareCurve = new Curve();
+            Curve trimmedReference = new Curve();
+            Curve trimmedCompareCurve = new Curve();
             TubeReport report = new TubeReport();
             TubeSize size = null;
             Tube tube = new Tube(size);
@@ -313,13 +378,26 @@ public Report CompareFiles(Log log, CsvFile csvBase, string sReportPath, ref Opt
                     else
                         log.WriteLine(LogLevel.Debug, "The resolution of the base x-axis is good.");
 
+                    common_interval = GetCommonInterval(compareCurve, reference);
+
+                    if (common_interval != null && common_interval[0] != common_interval[1])
+                    {
+                        trimmedReference = Trim_curve(reference, common_interval[0], common_interval[1]);
+                        trimmedCompareCurve = Trim_curve(compareCurve, common_interval[0], common_interval[1]);
+                    }
+                    else
+                    {
+                        trimmedReference = reference;
+                        trimmedCompareCurve = compareCurve;
+                    }
+
                     // The actual nominal attribute should be used, but is unfortunately unavailable in the CSV files.
                     // A default nominal value of 0.001 was chosen as a compromise between having many false negatives
                     // and passing wrong result files.
                     // See discussion in https://github.com/modelica/ModelicaStandardLibrary/issues/4421
                     const double defaultNominalValue = 0.001;
-                    const bool useLegacyBaseAndRatio = true;
-                    size = new TubeSize(reference, defaultNominalValue, useLegacyBaseAndRatio);
+                    const bool useLegacyBaseAndRatio = false;
+                    size = new TubeSize(trimmedReference, defaultNominalValue, useLegacyBaseAndRatio);
                     size.Calculate(_dRangeDelta, Axes.X, Relativity.Relative);
                     tube = new Tube(size);
                     var calcResult = tube.Calculate(reference);
@@ -331,7 +409,7 @@ public Report CompareFiles(Log log, CsvFile csvBase, string sReportPath, ref Opt
                         continue;
                     }
                     report = calcResult.Item1;
-                    bool validationSuccess = Tube.Validate(compareCurve, report);
+                    bool validationSuccess = Tube.Validate(trimmedCompareCurve, report);
                     if (!validationSuccess)
                     {
                         log.Error("Error in the validation of the tube. Skipping {0}", res.Key);
@@ -350,7 +428,7 @@ public Report CompareFiles(Log log, CsvFile csvBase, string sReportPath, ref Opt
                     }
                 }
                 if (null != report) //No charts for missing reports
-                    PrepareCharts(reference, compareCurve, report.Errors, rep, report, res, options.UseBitmapPlots);
+                    PrepareCharts(reference, compareCurve, trimmedCompareCurve, report.Errors, rep, report, res, options.UseBitmapPlots);
             }
             rep.Tolerance = _dRangeDelta;
 
@@ -391,10 +469,10 @@ public Report CompareFiles(Log log, CsvFile csvBase, string sReportPath, ref Opt
 
         private void PrepareCharts(Report rep, Curve compare)//Draw result only
         {
-            PrepareCharts(compare, null, null, rep, null, new KeyValuePair<string, List<double>>(compare.Name, null), false);
+            PrepareCharts(compare, null, null, null, rep, null, new KeyValuePair<string, List<double>>(compare.Name, null), false);
         }
 
-        private void PrepareCharts(Curve reference, Curve compare, Curve error, Report rep, TubeReport tubeReport, KeyValuePair<string, List<double>> res, bool bDrawBitmapPlots)
+        private void PrepareCharts(Curve reference, Curve compare, Curve trimmedCompare, Curve error, Report rep, TubeReport tubeReport, KeyValuePair<string, List<double>> res, bool bDrawBitmapPlots)
         {
 
             Chart ch = new Chart()
@@ -456,15 +534,15 @@ private void PrepareCharts(Curve reference, Curve compare, Curve error, Report r
             if (null != error && null != error.X && error.X.Length > 0)
             {
                 //Get complete error curve as "error" only holds error points
-                Curve curveErrors = new Curve("ERRORS", new double[compare.X.Length], new double[compare.X.Length]);
+                Curve curveErrors = new Curve("ERRORS", new double[trimmedCompare.X.Length], new double[trimmedCompare.X.Length]);
                 int j = 0;
-                for (int i = 0; i < compare.X.Length - 1; i++)
+                for (int i = 0; i <= trimmedCompare.X.Length - 1; i++)
                 {
-                    curveErrors.X[i] = compare.X[i];
-                    if (error.X.Contains(compare.X[i]))
+                    curveErrors.X[i] = trimmedCompare.X[i];
+                    if (error.X.Contains(trimmedCompare.X[i]))
                     {
                         curveErrors.Y[i] = (this._bShowRelativeErrors) ? error.Y[j] : 1;
-                        if (compare.X[i + 1] > compare.X[i])
+                        if (i == trimmedCompare.X.Length - 1 || trimmedCompare.X[i + 1] > trimmedCompare.X[i])
                         {
                             j++;
                         }
@@ -487,24 +565,24 @@ private void PrepareCharts(Curve reference, Curve compare, Curve error, Report r
                 //Calculate delta error
                 List<double> lDeltas = new List<double>();
                 j = 0;
-                for (int i = 1; i < compare.X.Length - 1; i++)
+                for (int i = 1; i < trimmedCompare.X.Length - 1; i++)
                 {
                     if (j < error.X.Length)
                     {
-                        while (compare.X[i] < error.X[j])
+                        while (trimmedCompare.X[i] < error.X[j])
                         {
                             i++;
                             continue;
                         }
 
-                        if (i < compare.X.Length - 1)
-                            lDeltas.Add((Math.Abs(error.Y[j]) * ((Math.Abs(compare.X[i] - compare.X[i - 1])) + (Math.Abs(compare.X[i + 1] - compare.X[i])))) / 2);
+                        if (i < trimmedCompare.X.Length - 1)
+                            lDeltas.Add((Math.Abs(error.Y[j]) * ((Math.Abs(trimmedCompare.X[i] - trimmedCompare.X[i - 1])) + (Math.Abs(trimmedCompare.X[i + 1] - trimmedCompare.X[i])))) / 2);
                         else // handle errors in the last point (there is no i+1)
-                            lDeltas.Add((Math.Abs(error.Y[j]) * ((Math.Abs(compare.X[i] - compare.X[i - 1])))) / 2);
+                            lDeltas.Add((Math.Abs(error.Y[j]) * ((Math.Abs(trimmedCompare.X[i] - trimmedCompare.X[i - 1])))) / 2);
                         j++;
                     }
                 }
-                ch.DeltaError = lDeltas.Sum() / (1e-3 + compare.Y.Max(x => Math.Abs(x)));
+                ch.DeltaError = lDeltas.Sum() / (1e-3 + trimmedCompare.Y.Max(x => Math.Abs(x)));
             }
             if (null != tubeReport && tubeReport.Lower.X.ToList<double>().Count > 2)//Remember Start and Stop values for graph scaling
             {

Modelica_ResultCompare/CurveCompare/Algorithms/Rectangle.cs

@@ -81,12 +81,25 @@ private static Curve CalculateLower(Curve reference, TubeSize size)
 
             // ignore identical point at the beginning
             b = 0;
+
+            // calculate slope at the beginning
             while (b + 1 < reference.Count && (reference.X[b] - reference.X[b + 1] == 0) && (reference.Y[b] - reference.Y[b + 1] == 0))
                 b++;
 
+            if (reference.X.Length == b + 1)
+                b--;
+
+            if (reference.X[b + 1] != reference.X[b])
+                m0 = (reference.Y[b + 1] - reference.Y[b]) / (reference.X[b + 1] - reference.X[b]);
+            else
+                m0 = 0;
+
             // add point down left
             LX.Add(reference.X[b] - size.X);
-            LY.Add(reference.Y[b] - size.Y);
+            if (m0 > 0)
+                LY.Add(reference.Y[b] - size.Y - 2 * m0 * size.X);
+            else
+                LY.Add(reference.Y[b] - size.Y);
 
             if (b + 1 < reference.Count)
             {
@@ -195,9 +208,17 @@ private static Curve CalculateLower(Curve reference, TubeSize size)
                 }
             }
 
+            // calculate slope at the end
+            if (reference.X[reference.Count - 1] != reference.X[reference.Count - 2])
+                m0 = (reference.Y[reference.Count - 1] - reference.Y[reference.Count - 2]) / (reference.X[reference.Count - 1] - reference.X[reference.Count - 2]);
+            else
+                m0 = 0;
             // add point down right
             LX.Add(reference.X[reference.Count - 1] + size.X);
-            LY.Add(reference.Y[reference.Count - 1] - size.Y);
+            if (m0 < 0)
+                LY.Add(reference.Y[reference.Count - 1] - size.Y + 2 * m0 * size.X);
+            else
+                LY.Add(reference.Y[reference.Count - 1] - size.Y);
 
             // -------------------------------------------------------------------------------------------------------------
             // -------------- 2. Remove points and add intersection points in case of backward order -----------------------
@@ -239,9 +260,21 @@ private static Curve CalculateUpper(Curve reference, TubeSize size)
             while (b + 1 < reference.Count && (reference.X[b] - reference.X[b + 1] == 0) && (reference.Y[b] - reference.Y[b + 1] == 0))
                 b++;
 
+            if (reference.X.Length == b + 1)
+                b--;
+
+            // calculate slope at the beginning
+            if (reference.X[b + 1] != reference.X[b])
+                m0 = (reference.Y[b + 1] - reference.Y[b]) / (reference.X[b + 1] - reference.X[b]);
+            else
+                m0 = 0;
+
             // add point top left
             UX.Add(reference.X[b] - size.X);
-            UY.Add(reference.Y[b] + size.Y);
+            if (m0 < 0)
+                UY.Add(reference.Y[b] + size.Y - 2 * m0 * size.X);
+            else
+                UY.Add(reference.Y[b] + size.Y);
 
             if (b + 1 < reference.Count)
             {
@@ -350,10 +383,18 @@ private static Curve CalculateUpper(Curve reference, TubeSize size)
                 }
             }
 
-            // add point top right
-            UX.Add(reference.X[reference.Count - 1] + size.X);
-            UY.Add(reference.Y[reference.Count - 1] + size.Y);
 
+            // calculate slope at the end
+            if (reference.X[reference.Count - 1] != reference.X[reference.Count - 2])
+                m0 = (reference.Y[reference.Count - 1] - reference.Y[reference.Count - 2]) / (reference.X[reference.Count - 1] - reference.X[reference.Count - 2]);
+            else
+                m0 = 0;
+            // add point down right
+            UX.Add(reference.X[reference.Count - 1] + size.X);
+            if (m0 > 0)
+                UY.Add(reference.Y[reference.Count - 1] + size.Y + 2 * m0 * size.X);
+            else
+                UY.Add(reference.Y[reference.Count - 1] + size.Y);
             // ---------------------------------------------------------------------------------------------------------
             // -------------- 2. Remove points and add intersection points in case of backward order -------------------
             // ---------------------------------------------------------------------------------------------------------
@@ -433,7 +474,7 @@ private static int RemoveLoop(List<double> X, List<double> Y, bool lower)
                         k++;
                         //while ((X[i] < X[k] || (X[i] == X[k] && Y[i] < Y[k])) && i < j)
                         while ((X[i] < X[k] || (lower && X[i] == X[k] && Y[i] < Y[k] && !(k + 1 < X.Count && X[k] == X[k + 1] && Y[k + 1] < Y[k])) || (!lower && X[i] == X[k] && Y[i] > Y[k] && !(k + 1 < X.Count && X[k] == X[k + 1] && Y[k + 1] > Y[k]))) && i < j)
-                                i++;
+                            i++;
                         // it holds X[i - 1] < X[k] <= X[i], particularly X[i] != X[i - 1]
                         // for i < j and X[i - 1] < X[k] it holds X[i - 1] < X[k] <= X[i], particularly X[i] != X[i - 1]
                         // linear interpolation of (x, y) = (X[k], y) on segment (i - 1, i)

Modelica_ResultCompare/CurveCompare/TubeSize.cs

@@ -14,7 +14,7 @@ namespace CurveCompare
     /// </summary>
     public class TubeSize
     {
-        private double x, y, baseX, baseY, ratio;
+        private double x, y, baseX, baseY, ratio, mean, width;
         private Curve reference;
         bool successful;
 
@@ -81,7 +81,8 @@ public TubeSize(Curve reference, double nominalValue, bool formerBaseAndRatio)
             if (formerBaseAndRatio)
                 SetFormerBaseAndRatio(nominalValue);
             else
-                SetStandardBaseAndRatio(nominalValue);
+                // SetStandardBaseAndRatio(nominalValue);
+                SetNewBaseAndRatio(nominalValue);
             successful = false;
         }
         /// <summary>
@@ -104,6 +105,27 @@ private void SetStandardBaseAndRatio(double nominalValue)
             else
                 ratio = 0;
         }
+
+        private void SetNewBaseAndRatio(double nominalValue)
+        {
+            // set baseX
+            baseX = reference.X.Max() - reference.X.Min(); //reference.X.Max() - reference.X.Min() + Math.Abs(reference.X.Min());
+            if (baseX == 0) // nonsense case, no data
+                baseX = Math.Abs(reference.X.Max());
+            if (baseX == 0) // nonsense case, no data
+                baseX = 1;
+            // set baseY
+            width = Math.Abs((reference.Y.Max() - reference.Y.Min()) / 2); // half the range/width of the curve in y- direction
+            mean = Math.Abs((reference.Y.Max() + reference.Y.Min()) / 4); // half the median of the curve in y- direction
+
+
+            baseY = Math.Max(width + mean, nominalValue);
+            // set ratio
+            if (baseX != 0)
+                ratio = baseY / baseX;
+            else
+                ratio = 0;
+        }
         /// <summary>
         /// Calculates former standard values for BaseX , BaseY and Ratio.
         /// </summary>
@@ -127,7 +149,7 @@ private void SetFormerBaseAndRatio(double nominalValue)
         /// Calculation fails, if (Ratio = 0 or BaseX = 0 or BaseY = 0).<para/>
         /// If calculation fails, [set Ratio and BaseX and BaseY != 0] or [call Calculate(double x, double y, Relativity relativity) with parameter Relativity.Absolute]</para></remarks>
         /// <exception cref="ArgumentOutOfRangeException">Relative value is out of expected range [0,1].</exception>
-        public void Calculate(double value, Axes axes, Relativity relativity)
+        public void Calculate(double value, double valueT, Axes axes, Relativity relativity)
         {
             successful = false;
 
@@ -138,16 +160,19 @@ public void Calculate(double value, Axes axes, Relativity relativity)
                     if ((value < 0) || (value > 1))
                         throw new ArgumentOutOfRangeException("Relative value is out of expected range [0,1].");
 
+                    if ((valueT < 0) || (valueT > 1))
+                        throw new ArgumentOutOfRangeException("Relative value is out of expected range [0,1].");
+
                     if (axes == Axes.Y && baseY > 0)
                     {
                         y = value * baseY;
-                        x = y / ratio;
+                        x = valueT * baseX;
                         successful = true;
                     }
                     else if (axes == Axes.X && baseX > 0)
                     {
-                        x = value * baseX;
-                        y = ratio * x;
+                        x = valueT * baseX;
+                        y = value * baseY;
                         successful = true;
                     }
                 }

Modelica_ResultCompare/Options.cs

@@ -49,7 +49,7 @@ public class Options
         [Option('i', "inline", DefaultValue = false, HelpText = "If set, javascript and style sheet files are inserted as inline text in every html output file")]
         public bool InlineScripts { get; set; }
 
-        [Option('t', "tolerance", Required = false, DefaultValue = "0.002", HelpText = "Set the width of the tube at discontinuity in x-direction [Default is 0.002].")]
+        [Option('t', "tolerance", Required = false, DefaultValue = "0.002", HelpText = "Set the width of the tube at discontinuity [Default is 0.002].")]
         public string Tolerance { get; set; }
 
         [Option('v', "verbosity", DefaultValue = 4, Required = false, HelpText = "Sets the verbosity of the output (1 most to 4[Default] less verbose).")]