Merge ea875ad6b9 into e989ab7302

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 }; 
+                                                                  
         /// <summary>
         /// Verifies that the interpolation matches the given value at all the provided sample points.
         /// </summary>
@@ -202,5 +205,37 @@ namespace MathNet.Numerics.UnitTests.InterpolationTests
 
             CollectionAssert.DoesNotContain(yipol, Double.NaN);
         }
+
+        /// <summary>
+        /// Tests that the cubic spline returns the correct t values where the derivative is 0
+        /// </summary>
+        [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.");
+         }
+
+        /// <summary>
+        /// Tests that the min and max values for the natural spline are correct
+        /// </summary>
+        [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.");
+
+        }
     }
 }

src/Numerics/Interpolation/CubicSpline.cs

@@ -1,4 +1,4 @@
-// <copyright file="CubicSpline.cs" company="Math.NET">
+// <copyright file="CubicSpline.cs" company="Math.NET">
 // Math.NET Numerics, part of the Math.NET Project
 // http://numerics.mathdotnet.com
 // http://github.com/mathnet/mathnet-numerics
@@ -623,5 +623,87 @@ namespace MathNet.Numerics.Interpolation
 
             return Math.Min(Math.Max(index, 0), _x.Length - 2);
         }
+
+        /// <summary>
+        /// Gets all the t values where the derivative is 0
+        /// </summary>
+        /// <returns>An array of t values (in the domain of teh function) where the derivative of the spline is 0</returns>
+        public double[] GetHorizontalDerivativeTValues()
+        {
+            List<double> points = new List<double>();
+            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();
+        }
+
+        /// <summary>
+        /// Returns the t values in the domain of the spline for which it takes the minimum and maximum value.
+        /// </summary>
+        /// <returns>A tuple containing the t value for which the spline is minimum in the first component and maximum in the second component </returns>
+        public Tuple<double, double> GetMinMaxTValues()
+        {
+            double max = double.MinValue;
+            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;
+                    maxT = p;
+                }
+                if (y < min)
+                {
+                    min = y;
+                    minT = p;
+                }
+            }
+            //go through the inflexion, local minimums and local maximums
+            foreach (double p in GetHorizontalDerivativeTValues())
+            {
+                double y = Interpolate(p);
+                if (y > max)
+                {
+                    max = y;
+                    maxT = p;
+                }
+                if (y < min)
+                {
+                    min = y;
+                    minT = p;
+                }
+            }
+            return new Tuple<double, double>(minT, maxT);
+        }
     }
 }