Renamed to StationaryPoints and Extrema, optimized Extrema and added large data set unit test

Author: NetPro69

src/Numerics.Tests/InterpolationTests/CubicSplineTest.cs

@@ -213,7 +213,7 @@ public void InterpolateAkimaSorted_MustBeThreadSafe_GitHub219([Values(8, 32, 256
         public void NaturalSplineGetHorizontalDerivativeTValues()
         {
             CubicSpline it = CubicSpline.InterpolateBoundaries(_x, _z, SplineBoundaryCondition.SecondDerivative, -4.0, SplineBoundaryCondition.SecondDerivative, 4.0);
-            var horizontalDerivatives = it.GetHorizontalDerivativeTValues();
+            var horizontalDerivatives = it.StationaryPoints();
             //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");
@@ -231,11 +231,70 @@ public void NaturalSplineGetHorizontalDerivativeTValues()
         public void NaturalSplineGetMinMaxTvalues()
         {
             CubicSpline it = CubicSpline.InterpolateBoundaries(_x, _z, SplineBoundaryCondition.SecondDerivative, -4.0, SplineBoundaryCondition.SecondDerivative, 4.0);
-            var minMax = it.GetMinMaxTValues();
+            var minMax = it.Extrema();
             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.");
+        }
 
+        /// <summary>
+        /// This tests that the min and max values for a natural spline interpolated on an oscilating decaying function are correct and
+        /// spints out the time it takes to do the calculation for 100 thousand points
+        /// </summary>
+        [Test]
+        public void CheckNaturalSplineMinMaxValuesPerformance()
+        {
+            //first generate the test data and spline
+            double amplitude = 100;
+            double ofset = 0;
+            double period = 10;
+            double decay = 0.1;
+            double minX = -50;
+            double maxX = 50;
+            int points = 100000;
+            var data = GenerateSplineData(amplitude, period, decay, minX, maxX, ofset, points);
+            CubicSpline it = CubicSpline.InterpolateBoundaries(data.Item1, data.Item2, SplineBoundaryCondition.Natural, 0, SplineBoundaryCondition.Natural, 0);
+            //start the time and calculate the min max values
+            var t = DateTime.Now;
+            var minMax = it.Extrema();
+            Assert.IsTrue(minMax.Item2.AlmostEqual(ofset, 0.3), "Expexted max value near ofset.");
+            Assert.IsTrue(minMax.Item1.AlmostEqual(ofset+period/2, 0.3) || minMax.Item1.AlmostEqual(ofset - period / 2, 0.3), "Expexted min value near ofset +- period/2.");
+            //spit out the time it took to calculate
+            Console.WriteLine("Extrema took: " + (DateTime.Now - t).TotalMilliseconds.ToString("000.00")+" ms for "+points.ToString()+" points.");
+            //determine if the values are correct
+            var sp = it.StationaryPoints();
+            foreach (var x in sp)
+            {
+                //check that the stationary point falls roughly at a half period
+                Assert.IsTrue(Math.Abs((Math.Abs((x - ofset) * 2 / period) - Math.Round(Math.Abs(x - ofset) * 2 / period, 0))).AlmostEqual(0,0.3), "Stationary point found outside of period/2 for x="+x.ToString());
+            }
         }
+
+        /// <summary>
+        /// Generates a set of points representing an oscilating decaying function 
+        /// </summary>
+        /// <param name="amplitude">The max amplitude</param>
+        /// <param name="period">The period of oscilation</param>
+        /// <param name="decay">The decaying exponent, the larger the value the faster the functioon decays</param>
+        /// <param name="minX">The min value for X</param>
+        /// <param name="maxX">The max value for X</param>
+        /// <param name="ofset">The x - ofset of the max value of the function</param>
+        /// <param name="points">The number of points to generate, must be greater than 1</param>
+        /// <returns></returns>
+        private Tuple<double[],double[]> GenerateSplineData(double amplitude, double period, double decay, double minX, double maxX,  double ofset, int points)
+        {
+            double delta = (maxX-minX)/(points-1);
+            double[] _x = new double[points];
+            double[] _y = new double[points];   
+            for (int i = 0; i < points; i++)
+            {
+                double x = i * delta + minX;
+                double y = amplitude * Math.Cos((x - ofset) * 2 * Math.PI / period) * Math.Exp(-Math.Abs((x - ofset) * decay));
+                _x[i] = x;
+                _y[i] = y;
+            }
+            return Tuple.Create(_x, _y);
+        }
+
     }
 }

src/Numerics/Interpolation/CubicSpline.cs

@@ -626,9 +626,10 @@ int LeftSegmentIndex(double t)
 
         /// <summary>
         /// Gets all the t values where the derivative is 0
+        /// see: https://mathworld.wolfram.com/StationaryPoint.html
         /// </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()
+        /// <returns>An array of t values (in the domain of the function) where the derivative of the spline is 0</returns>
+        public double[] StationaryPoints()
         {
             List<double> points = new List<double>();
             for (int index = 0; index < _x.Length - 1; index++)
@@ -667,29 +668,31 @@ public double[] GetHorizontalDerivativeTValues()
         /// 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()
+        public Tuple<double, double> Extrema()
         {
-            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())
+            //Check the edges of the domain
+            //set the initial values to the leftmost domain point
+            double t = _x[0];
+            double max = Interpolate(t);
+            double min = max;
+            double minT = t;
+            double maxT = t;
+            //check the rightmost domain point
+            t = _x[_x.Length-1];
+            var ty = Interpolate(t);
+            if (ty > max)
+            {
+                max = ty;
+                maxT = t;
+            }
+            if (ty < min)
+            {
+                min = ty;
+                minT = t;
+            }
+            //check the the inflexion, local minimums and local maximums
+            double[] pointsToCheck = StationaryPoints();
+            foreach (double p in pointsToCheck)
             {
                 double y = Interpolate(p);
                 if (y > max)