From 01ed5cf41f81f658f9c306995f9ce86feb7d1698 Mon Sep 17 00:00:00 2001
From: Christoph Ruegg <git@cdrnet.ch>
Date: Fri, 20 Dec 2013 21:46:37 +0100
Subject: [PATCH] Interpolation: add QuadraticSpline (without solver yet)

---
 src/Numerics/Interpolation/QuadraticSpline.cs | 166 ++++++++++++++++++
 src/Numerics/Numerics.csproj                  |   1 +
 2 files changed, 167 insertions(+)
 create mode 100644 src/Numerics/Interpolation/QuadraticSpline.cs

diff --git a/src/Numerics/Interpolation/QuadraticSpline.cs b/src/Numerics/Interpolation/QuadraticSpline.cs
new file mode 100644
index 00000000..948fe518
--- /dev/null
+++ b/src/Numerics/Interpolation/QuadraticSpline.cs
@@ -0,0 +1,166 @@
+// <copyright file="QuadraticSpline.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+// http://mathnetnumerics.codeplex.com
+//
+// Copyright (c) 2009-2013 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using System;
+
+namespace MathNet.Numerics.Interpolation
+{
+    /// <summary>
+    /// Quadratic Spline Interpolation.
+    /// </summary>
+    /// <remarks>Supports both differentiation and integration.</remarks>
+    public class QuadraticSpline : IInterpolation
+    {
+        readonly double[] _x;
+        readonly double[] _c0;
+        readonly double[] _c1;
+        readonly double[] _c2;
+        readonly Lazy<double[]> _indefiniteIntegral;
+
+        /// <param name="x">sample points (N+1), sorted ascending</param>
+        /// <param name="c0">Zero order spline coefficients (N)</param>
+        /// <param name="c1">First order spline coefficients (N)</param>
+        /// <param name="c2">second order spline coefficients (N)</param>
+        public QuadraticSpline(double[] x, double[] c0, double[] c1, double[] c2)
+        {
+            _x = x;
+            _c0 = c0;
+            _c1 = c1;
+            _c2 = c2;
+            _indefiniteIntegral = new Lazy<double[]>(ComputeIndefiniteIntegral);
+        }
+
+        /// <summary>
+        /// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative).
+        /// </summary>
+        public bool SupportsDifferentiation
+        {
+            get { return true; }
+        }
+
+        /// <summary>
+        /// Gets a value indicating whether the algorithm supports integration (interpolated quadrature).
+        /// </summary>
+        public bool SupportsIntegration
+        {
+            get { return true; }
+        }
+
+        /// <summary>
+        /// Interpolate at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated value x(t).</returns>
+        public double Interpolate(double t)
+        {
+            int k = LeftBracketIndex(t);
+            var x = (t - _x[k]);
+            return _c0[k] + x*(_c1[k] + x*_c2[k]);
+        }
+
+        /// <summary>
+        /// Differentiate at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated first derivative at point t.</returns>
+        public double Differentiate(double t)
+        {
+            int k = LeftBracketIndex(t);
+            return _c1[k] + (t - _x[k])*2*_c2[k];
+        }
+
+        /// <summary>
+        /// Differentiate twice at point t.
+        /// </summary>
+        /// <param name="t">Point t to interpolate at.</param>
+        /// <returns>Interpolated second derivative at point t.</returns>
+        public double Differentiate2(double t)
+        {
+            int k = LeftBracketIndex(t);
+            return 2*_c2[k];
+        }
+
+        /// <summary>
+        /// Indefinite integral at point t.
+        /// </summary>
+        /// <param name="t">Point t to integrate at.</param>
+        public double Integrate(double t)
+        {
+            int k = LeftBracketIndex(t);
+            var x = (t - _x[k]);
+            return _indefiniteIntegral.Value[k] + x*(_c0[k] + x*(_c1[k]/2 + x*_c2[k]/3));
+        }
+
+        /// <summary>
+        /// Definite integral between points a and b.
+        /// </summary>
+        /// <param name="a">Left bound of the integration interval [a,b].</param>
+        /// <param name="b">Right bound of the integration interval [a,b].</param>
+        public double Integrate(double a, double b)
+        {
+            return Integrate(b) - Integrate(a);
+        }
+
+        double[] ComputeIndefiniteIntegral()
+        {
+            var integral = new double[_c1.Length];
+            for (int i = 0; i < integral.Length - 1; i++)
+            {
+                double w = _x[i + 1] - _x[i];
+                integral[i + 1] = integral[i] + w*(_c0[i] + w*(_c1[i]/2 + w*_c2[i]/3));
+            }
+            return integral;
+        }
+
+        /// <summary>
+        /// Find the index of the greatest sample point smaller than t.
+        /// </summary>
+        int LeftBracketIndex(double t)
+        {
+            // Binary search in the [ t[0], ..., t[n-2] ] (t[n-1] is not included)
+            int low = 0;
+            int high = _x.Length - 1;
+            while (low != high - 1)
+            {
+                int middle = (low + high)/2;
+                if (_x[middle] > t)
+                {
+                    high = middle;
+                }
+                else
+                {
+                    low = middle;
+                }
+            }
+
+            return low;
+        }
+    }
+}
diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj
index 37cbe59c..ff8ce74b 100644
--- a/src/Numerics/Numerics.csproj
+++ b/src/Numerics/Numerics.csproj
@@ -90,6 +90,7 @@
     <Compile Include="Euclid.cs" />
     <Compile Include="Generate.cs" />
     <Compile Include="GoodnessOfFit.cs" />
+    <Compile Include="Interpolation\QuadraticSpline.cs" />
     <Compile Include="Precision.Comparison.cs" />
     <Compile Include="Precision.Equality.cs" />
     <Compile Include="Distributions\Bernoulli.cs" />