From 1cc14ced95db30d4df82a66908c107d07e9868f2 Mon Sep 17 00:00:00 2001 From: Candy Chiu Date: Tue, 15 Jul 2014 15:34:28 -0400 Subject: [PATCH] Added a Step function interpolator. --- src/Numerics/Interpolate.cs | 21 +++ .../Interpolation/StepInterpolation.cs | 146 ++++++++++++++++++ src/Numerics/Numerics.csproj | 1 + .../StepInterpolationTest.cs | 89 +++++++++++ src/UnitTests/UnitTests.csproj | 1 + 5 files changed, 258 insertions(+) create mode 100644 src/Numerics/Interpolation/StepInterpolation.cs create mode 100644 src/UnitTests/InterpolationTests/StepInterpolationTest.cs diff --git a/src/Numerics/Interpolate.cs b/src/Numerics/Interpolate.cs index f074b6cb..9049f55f 100644 --- a/src/Numerics/Interpolate.cs +++ b/src/Numerics/Interpolate.cs @@ -223,5 +223,26 @@ namespace MathNet.Numerics { return Interpolation.CubicSpline.InterpolateHermite(points, values, firstDerivatives); } + + /// + /// Create log linear spline interpolation based on arbitrary points. + /// + /// The sample points t. Optimized for arrays. + /// The sample point values x(t). Optimized for arrays. + /// + /// An interpolation scheme optimized for the given sample points and values, + /// which can then be used to compute interpolations and extrapolations + /// on arbitrary points. + /// + /// + /// The value pairs do not have to be sorted, but if they are not sorted ascendingly + /// and the passed x and y arguments are arrays, they will be sorted inplace and thus modified. + /// + /// If the values are passed as an array, they will be modified inplace, even it is already sorted. + /// + public static IInterpolation Step(IEnumerable points, IEnumerable values) + { + return new StepInterpolation(points, values); + } } } diff --git a/src/Numerics/Interpolation/StepInterpolation.cs b/src/Numerics/Interpolation/StepInterpolation.cs new file mode 100644 index 00000000..18892568 --- /dev/null +++ b/src/Numerics/Interpolation/StepInterpolation.cs @@ -0,0 +1,146 @@ +// +// 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-2014 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. +// + +using System; +using System.Collections.Generic; +using System.Linq; +using MathNet.Numerics.Properties; + +namespace MathNet.Numerics.Interpolation +{ + /// + /// A step function where the start of each segment is included, and last is excluded. Segment i is [x_i, x_i+1). + /// The domain of the function is all real numbers, such that y = 0 where x < x_0 or x gt; x_n + /// + public class StepInterpolation : IInterpolation + { + readonly double[] _x; + readonly double[] _y; + readonly Lazy _indefiniteIntegral; + + /// Sample points (N+1) sorted in ascending order + /// Functional value (N) of each segment. + public StepInterpolation(IEnumerable x, IEnumerable y) + { + var xx = (x as double[]) ?? x.ToArray(); + var yy = (y as double[]) ?? y.ToArray(); + + if (xx.Length != yy.Length + 1) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + _x = xx; + _y = yy; + _indefiniteIntegral = new Lazy(ComputeIndefiniteIntegral); + } + + double[] ComputeIndefiniteIntegral() + { + var integral = new double[_x.Length]; + for (int i = 0; i < integral.Length - 1; i++) + { + integral[i + 1] = integral[i] + (_x[i + 1] - _x[i])*_y[i]; + } + return integral; + } + + public bool SupportsDifferentiation + { + get { return true; } + } + + public bool SupportsIntegration + { + get { return true; } + } + + public double Differentiate(double t) + { + int index = Array.BinarySearch(_x, t); + if (index >= 0) + return double.NaN; + return 0d; + } + + public double Differentiate2(double t) + { + return Differentiate(t); + } + + /// + /// Indefinite integral at point t. + /// + /// Point t to integrate at. + public double Integrate(double t) + { + if (t <= _x[0]) + return 0.0; + int last = _x.Length - 1; + if (t >= _x[last]) + return _indefiniteIntegral.Value[last]; + + int k = LeftBracketIndex(t); + var x = (t - _x[k]); + return _indefiniteIntegral.Value[k] + x*_y[k]; + } + + /// + /// Definite integral between points a and b. + /// + /// Left bound of the integration interval [a,b]. + /// Right bound of the integration interval [a,b]. + public double Integrate(double a, double b) + { + return Integrate(b) - Integrate(a); + } + + public double Interpolate(double t) + { + if (t < _x[0] || t >= _x[_x.Length - 1]) + return 0.0; + + int k = LeftBracketIndex(t); + return _y[k]; + } + + /// + /// Find the index of the greatest sample point smaller than t. + /// + int LeftBracketIndex(double t) + { + int index = Array.BinarySearch(_x, t); + if (index >= 0) + return index; + index = ~index; + return index - 1; + } + } +} diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index 106dec4f..b0600c03 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -92,6 +92,7 @@ + diff --git a/src/UnitTests/InterpolationTests/StepInterpolationTest.cs b/src/UnitTests/InterpolationTests/StepInterpolationTest.cs new file mode 100644 index 00000000..d84ff19f --- /dev/null +++ b/src/UnitTests/InterpolationTests/StepInterpolationTest.cs @@ -0,0 +1,89 @@ +// +// 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) 2002-2014 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. +// + +using MathNet.Numerics.Interpolation; +using NUnit.Framework; + +namespace MathNet.Numerics.UnitTests.InterpolationTests +{ + [TestFixture, Category("Interpolation")] + public class StepInterpolationTest + { + readonly double[] _t = { -2.0, -1.0, 0.0, 1.0, 2.0, 3.0 }; + readonly double[] _y = { 1.0, 2.0, -1.0, 0.0, 1.0 }; + + [Test] + public void FirstDerivative() + { + IInterpolation ip = new StepInterpolation(_t, _y); + Assert.That(ip.Differentiate(-3.0), Is.EqualTo(0.0)); + Assert.That(ip.Differentiate(-2.0), Is.EqualTo(double.NaN)); + Assert.That(ip.Differentiate(-1.5), Is.EqualTo(0.0)); + Assert.That(ip.Differentiate(-1.0), Is.EqualTo(double.NaN)); + Assert.That(ip.Differentiate(-0.5), Is.EqualTo(0.0)); + Assert.That(ip.Differentiate(0.0), Is.EqualTo(double.NaN)); + Assert.That(ip.Differentiate(0.5), Is.EqualTo(0.0)); + Assert.That(ip.Differentiate(1.0), Is.EqualTo(double.NaN)); + Assert.That(ip.Differentiate(2.0), Is.EqualTo(double.NaN)); + Assert.That(ip.Differentiate(3.0), Is.EqualTo(double.NaN)); + } + + [Test] + public void DefiniteIntegral() + { + IInterpolation ip = new StepInterpolation(_t, _y); + Assert.That(ip.Integrate(-3.0, -2.0), Is.EqualTo(0.0)); + Assert.That(ip.Integrate(-2.0, -1.0), Is.EqualTo(1.0)); + Assert.That(ip.Integrate(-1.0, 0.0), Is.EqualTo(2.0)); + Assert.That(ip.Integrate(0.0, 1.0), Is.EqualTo(-1.0)); + Assert.That(ip.Integrate(1.0, 2.0), Is.EqualTo(0.0)); + Assert.That(ip.Integrate(2.0, 3.0), Is.EqualTo(1.0)); + Assert.That(ip.Integrate(3.0, 4.0), Is.EqualTo(0.0)); + Assert.That(ip.Integrate(0.0, 4.0), Is.EqualTo(0.0)); + Assert.That(ip.Integrate(-3.0, -1.0), Is.EqualTo(1.0)); + Assert.That(ip.Integrate(-3.0, 4.0), Is.EqualTo(3.0)); + Assert.That(ip.Integrate(0.5, 1.5), Is.EqualTo(-0.5)); + Assert.That(ip.Integrate(-1.5, -0.5), Is.EqualTo(1.5)); + } + + /// + /// Verifies that the interpolation matches the given value at all the provided sample points. + /// + [Test] + public void FitsAtSamplePoints() + { + IInterpolation ip = new StepInterpolation(_t, _y); + for (int i = 0; i < _y.Length; i++) + { + Assert.AreEqual(_y[i], ip.Interpolate(_t[i]), "A Exact Point " + i); + } + } + } +} diff --git a/src/UnitTests/UnitTests.csproj b/src/UnitTests/UnitTests.csproj index 96b31628..fe2b9319 100644 --- a/src/UnitTests/UnitTests.csproj +++ b/src/UnitTests/UnitTests.csproj @@ -147,6 +147,7 @@ +