From a89fcd09ac3ca39d72d7b4e3e66cb0890672107d Mon Sep 17 00:00:00 2001 From: Eduardo Rojas Date: Wed, 15 Dec 2021 11:42:19 -0300 Subject: [PATCH 1/3] Added function go get points where cubic spline derivative is 0 (GetHorizontalPoints) and also to get the maximum value of the cubic spline function within its domain (GetMaxPoints) --- src/Numerics/Interpolation/CubicSpline.cs | 70 +++++++++++++++++++++++ 1 file changed, 70 insertions(+) diff --git a/src/Numerics/Interpolation/CubicSpline.cs b/src/Numerics/Interpolation/CubicSpline.cs index 1830da86..f3446667 100644 --- a/src/Numerics/Interpolation/CubicSpline.cs +++ b/src/Numerics/Interpolation/CubicSpline.cs @@ -623,5 +623,75 @@ namespace MathNet.Numerics.Interpolation return Math.Min(Math.Max(index, 0), _x.Length - 2); } + + /// + /// Gets all the points (x values) where the derivative is 0 + /// + /// An array of X values where the derivative of the spline is 0 + public double[] GetHorizontalPoints() + { + List points = new List(); + for (int index = 0; index < _x.Length-1; index++) + { + double a = 6 * _c3[index]; //derive ax^3 and multiply by 2 + double b = 2 * _c2[index]; //derive bx^2 + double c = _c1[index];//derive cx + double d = b * b - 2 * a * c; + //first check if a is 0, if so its a linear function, this happens with quadratic condition + if (a.AlmostEqual(0)) + { + double x = _x[index] - c / b; + //check if the result is in the domain + if (_x[index] <= x && x <= _x[index + 1]) points.Add(x); + } + else if (d.AlmostEqual(0))//its a quadratic with a single solution + { + double x = _x[index] - b / a; + if (_x[index] <= x && x <= _x[index + 1]) points.Add(x); + } + else if (d > 0)//only has a solution if d is greater than 0 + { + d = (double)System.Math.Sqrt(d); + //apply quadratic equation + double x1 = _x[index]+ (-b + d) / a; + double x2 = _x[index]+ (-b - d) / a; + //Add any solution points that fall within the domain to the list + if ((_x[index] <= x1) && (x1 <= _x[index + 1])) points.Add(x1); + if ((_x[index] <= x2) && (x2 <= _x[index + 1])) points.Add(x2); + } + } + return points.ToArray(); + } + + /// + /// Returns the minimum value for the spline within its domain. + /// + /// + public double GetMaxPoint() + { + double max = double.MinValue; + double x=0; + //go through the functions points to check if one of them has a higher value + foreach (double p in _x) + { + double y = Interpolate(p); + if (y > max) + { + max = y; + x = p; + } + } + //go through the inflexion, local minimums and local maximums + foreach (double p in GetHorizontalPoints()) + { + double y = Interpolate(p); + if (y > max) + { + max = y; + x = p; + } + } + return x; + } } } From c6d10b2e18a5bba8f1bea3ab46108b7b568f6380 Mon Sep 17 00:00:00 2001 From: Eduardo Rojas <30539773+NetPro69@users.noreply.github.com> Date: Fri, 17 Dec 2021 12:25:40 -0300 Subject: [PATCH 2/3] Update CubicSpline.cs Fixed function documentation error, changed minimum for maximum --- src/Numerics/Interpolation/CubicSpline.cs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Numerics/Interpolation/CubicSpline.cs b/src/Numerics/Interpolation/CubicSpline.cs index f3446667..cb9d0bdf 100644 --- a/src/Numerics/Interpolation/CubicSpline.cs +++ b/src/Numerics/Interpolation/CubicSpline.cs @@ -664,7 +664,7 @@ namespace MathNet.Numerics.Interpolation } /// - /// Returns the minimum value for the spline within its domain. + /// Returns the maximum value for the spline within its domain. /// /// public double GetMaxPoint() From ea875ad6b9c186f3f1e7f805f39b05f554055de5 Mon Sep 17 00:00:00 2001 From: Eduardo Rojas Date: Fri, 24 Dec 2021 16:28:47 -0300 Subject: [PATCH 3/3] Renamed GetHorizontalPoints to GetHorizontalDerivativeTValues and renamed and extended GetMaxPoint to GetMinMaxTValues. Aslo added unit tests for both functions. --- .../InterpolationTests/CubicSplineTest.cs | 37 ++++++++++- src/Numerics/Interpolation/CubicSpline.cs | 66 +++++++++++-------- 2 files changed, 75 insertions(+), 28 deletions(-) diff --git a/src/Numerics.Tests/InterpolationTests/CubicSplineTest.cs b/src/Numerics.Tests/InterpolationTests/CubicSplineTest.cs index 2ba52584..bdeb210d 100644 --- a/src/Numerics.Tests/InterpolationTests/CubicSplineTest.cs +++ b/src/Numerics.Tests/InterpolationTests/CubicSplineTest.cs @@ -39,7 +39,10 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests { readonly double[] _t = { -2.0, -1.0, 0.0, 1.0, 2.0 }; readonly double[] _y = { 1.0, 2.0, -1.0, 0.0, 1.0 }; - + //test data for min max values + readonly double[] _x = { -4, -3, -2, -1, 0, 1, 2, 3, 4 }; + readonly double[] _z = { -7, 2, 5, 0, -3, -1, -4, 0, 6 }; + /// /// Verifies that the interpolation matches the given value at all the provided sample points. /// @@ -202,5 +205,37 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests CollectionAssert.DoesNotContain(yipol, Double.NaN); } + + /// + /// Tests that the cubic spline returns the correct t values where the derivative is 0 + /// + [Test] + public void NaturalSplineGetHorizontalDerivativeTValues() + { + CubicSpline it = CubicSpline.InterpolateBoundaries(_x, _z, SplineBoundaryCondition.SecondDerivative, -4.0, SplineBoundaryCondition.SecondDerivative, 4.0); + var horizontalDerivatives = it.GetHorizontalDerivativeTValues(); + //readonly double[] _x = { -4, -3, -2, -1, 0, 1, 2, 3, 4 }; + //readonly double[] _z = { -7, 2, 5, 0, -3, -1, -4, 0, 6 }; + Assert.AreEqual(4, horizontalDerivatives.Length, "Incorrect number of points with derivative value equal to 0"); + Assert.IsTrue(horizontalDerivatives[0]>=-3 && horizontalDerivatives[0] <= -2,"Spline returns wrong t value: "+horizontalDerivatives[0]+" for first point"); + Assert.IsTrue(horizontalDerivatives[1] >= -1 && horizontalDerivatives[1] <= 0, "Spline returns wrong t value: " + horizontalDerivatives[1] + " for second point"); + Assert.IsTrue(horizontalDerivatives[2] >= 0 && horizontalDerivatives[2] <= 1, "Spline returns wrong t value: " + horizontalDerivatives[2] + " for third point"); + Assert.IsTrue(horizontalDerivatives[3] >= 2 && horizontalDerivatives[3] <= 3, "Spline returns wrong t value: " + horizontalDerivatives[3] + " for fourth point"); + Console.WriteLine("GetHorizontalDerivativeTValues checked out ok for cubic spline."); + } + + /// + /// Tests that the min and max values for the natural spline are correct + /// + [Test] + public void NaturalSplineGetMinMaxTvalues() + { + CubicSpline it = CubicSpline.InterpolateBoundaries(_x, _z, SplineBoundaryCondition.SecondDerivative, -4.0, SplineBoundaryCondition.SecondDerivative, 4.0); + var minMax = it.GetMinMaxTValues(); + Assert.AreEqual(-4, minMax.Item1, "Spline returns wrong t value for global minimum"); + Assert.AreEqual(4, minMax.Item2, "Spline returns wrong t value for global maximum"); + Console.WriteLine("GetMinMaxTValues checked out ok for cubic spline."); + + } } } diff --git a/src/Numerics/Interpolation/CubicSpline.cs b/src/Numerics/Interpolation/CubicSpline.cs index cb9d0bdf..2d00bce9 100644 --- a/src/Numerics/Interpolation/CubicSpline.cs +++ b/src/Numerics/Interpolation/CubicSpline.cs @@ -1,4 +1,4 @@ -// +// // Math.NET Numerics, part of the Math.NET Project // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics @@ -625,13 +625,13 @@ namespace MathNet.Numerics.Interpolation } /// - /// Gets all the points (x values) where the derivative is 0 + /// Gets all the t values where the derivative is 0 /// - /// An array of X values where the derivative of the spline is 0 - public double[] GetHorizontalPoints() + /// An array of t values (in the domain of teh function) where the derivative of the spline is 0 + public double[] GetHorizontalDerivativeTValues() { List points = new List(); - for (int index = 0; index < _x.Length-1; index++) + for (int index = 0; index < _x.Length - 1; index++) { double a = 6 * _c3[index]; //derive ax^3 and multiply by 2 double b = 2 * _c2[index]; //derive bx^2 @@ -640,58 +640,70 @@ namespace MathNet.Numerics.Interpolation //first check if a is 0, if so its a linear function, this happens with quadratic condition if (a.AlmostEqual(0)) { - double x = _x[index] - c / b; - //check if the result is in the domain - if (_x[index] <= x && x <= _x[index + 1]) points.Add(x); + double x = _x[index] - c / b; + //check if the result is in the domain + if (_x[index] <= x && x <= _x[index + 1]) points.Add(x); } else if (d.AlmostEqual(0))//its a quadratic with a single solution { - double x = _x[index] - b / a; - if (_x[index] <= x && x <= _x[index + 1]) points.Add(x); + double x = _x[index] - b / a; + if (_x[index] <= x && x <= _x[index + 1]) points.Add(x); } else if (d > 0)//only has a solution if d is greater than 0 { - d = (double)System.Math.Sqrt(d); - //apply quadratic equation - double x1 = _x[index]+ (-b + d) / a; - double x2 = _x[index]+ (-b - d) / a; - //Add any solution points that fall within the domain to the list - if ((_x[index] <= x1) && (x1 <= _x[index + 1])) points.Add(x1); - if ((_x[index] <= x2) && (x2 <= _x[index + 1])) points.Add(x2); + d = (double)System.Math.Sqrt(d); + //apply quadratic equation + double x1 = _x[index] + (-b + d) / a; + double x2 = _x[index] + (-b - d) / a; + //Add any solution points that fall within the domain to the list + if ((_x[index] <= x1) && (x1 <= _x[index + 1])) points.Add(x1); + if ((_x[index] <= x2) && (x2 <= _x[index + 1])) points.Add(x2); } } return points.ToArray(); } /// - /// Returns the maximum value for the spline within its domain. + /// Returns the t values in the domain of the spline for which it takes the minimum and maximum value. /// - /// - public double GetMaxPoint() + /// A tuple containing the t value for which the spline is minimum in the first component and maximum in the second component + public Tuple GetMinMaxTValues() { double max = double.MinValue; - double x=0; + double min = double.MaxValue; + double minT = 0; + double maxT = 0; //go through the functions points to check if one of them has a higher value foreach (double p in _x) { double y = Interpolate(p); if (y > max) { - max = y; - x = p; + max = y; + maxT = p; + } + if (y < min) + { + min = y; + minT = p; } } //go through the inflexion, local minimums and local maximums - foreach (double p in GetHorizontalPoints()) + foreach (double p in GetHorizontalDerivativeTValues()) { double y = Interpolate(p); if (y > max) { - max = y; - x = p; + max = y; + maxT = p; + } + if (y < min) + { + min = y; + minT = p; } } - return x; + return new Tuple(minT, maxT); } } }