From 9a57ff4894e1bb0fba04481ff193100abaa9ecb1 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Thu, 2 Jul 2009 00:56:45 +0800 Subject: [PATCH] interpolation: rational pole-free, linear spline Signed-off-by: Christoph Ruegg --- .../Algorithms/BarycentricInterpolation.cs | 226 ++++++++++++++++ .../Algorithms/LinearSplineInterpolation.cs | 177 +++++++++++++ .../RationalPoleFreeInterpolation.cs | 248 +++++++++++++++++ .../Algorithms/SplineInterpolation.cs | 249 ++++++++++++++++++ src/Managed/Interpolation/Interpolation.cs | 13 +- src/Managed/Managed.csproj | 4 + src/Native/Native.csproj | 12 + 7 files changed, 925 insertions(+), 4 deletions(-) create mode 100644 src/Managed/Interpolation/Algorithms/BarycentricInterpolation.cs create mode 100644 src/Managed/Interpolation/Algorithms/LinearSplineInterpolation.cs create mode 100644 src/Managed/Interpolation/Algorithms/RationalPoleFreeInterpolation.cs create mode 100644 src/Managed/Interpolation/Algorithms/SplineInterpolation.cs diff --git a/src/Managed/Interpolation/Algorithms/BarycentricInterpolation.cs b/src/Managed/Interpolation/Algorithms/BarycentricInterpolation.cs new file mode 100644 index 00000000..43c9a770 --- /dev/null +++ b/src/Managed/Interpolation/Algorithms/BarycentricInterpolation.cs @@ -0,0 +1,226 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 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. +// + +namespace MathNet.Numerics.Interpolation.Algorithms +{ + using System; + using System.Collections.Generic; + + /// + /// Barycentric Interpolation Algorithm. + /// + /// + /// This algorithm neither supports differentiation nor integration. + /// + public class BarycentricInterpolation : IInterpolation + { + /// + /// Sample Points t. + /// + private IList points; + + /// + /// Sample Values x(t). + /// + private IList values; + + /// + /// Barycentric Weights w(t). + /// + private IList weights; + + /// + /// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative). + /// + /// + /// + bool IInterpolation.SupportsDifferentiation + { + get { return false; } + } + + /// + /// Gets a value indicating whether the algorithm supports integration (interpolated quadrature). + /// + /// + bool IInterpolation.SupportsIntegration + { + get { return false; } + } + + /// + /// Initialize the interpolation method with the given sample set. + /// + /// Sample Points t + /// Sample Values x(t) + /// Barycentric weights w(t) + public + void + Initialize( + IList samplePoints, + IList sampleValues, + IList barycentricWeights) + { + if (null == samplePoints) + { + throw new ArgumentNullException("samplePoints"); + } + + if (null == sampleValues) + { + throw new ArgumentNullException("sampleValues"); + } + + if (null == barycentricWeights) + { + throw new ArgumentNullException("barycentricWeights"); + } + + if (samplePoints.Count < 1) + { + throw new ArgumentOutOfRangeException("samplePoints"); + } + + if (samplePoints.Count != sampleValues.Count) + { + throw new ArgumentException(Properties.Resources.ArgumentVectorsSameLengths); + } + + if (samplePoints.Count != barycentricWeights.Count) + { + throw new ArgumentException(Properties.Resources.ArgumentVectorsSameLengths); + } + + this.points = samplePoints; + this.values = sampleValues; + this.weights = barycentricWeights; + } + + /// + /// Interpolate at point t. + /// + /// Point t to interpolate at. + /// Interpolated value x(t). + public + double + Interpolate(double t) + { + // trivial case: only one sample? + if (this.points.Count == 1) + { + return this.values[0]; + } + + // evaluate closest point and offset from that point + int closestPoint = 0; + double offset = t - this.points[0]; + for (int i = 1; i < this.points.Count; i++) + { + if (Math.Abs(t - this.points[i]) < Math.Abs(offset)) + { + offset = t - this.points[i]; + closestPoint = i; + } + } + + // trivial case: on a known sample point? + // TODO: Number.AlmostZero(offset) instead of == + if (offset == 0.0) + { + return this.values[closestPoint]; + } + + if (Math.Abs(offset) > 1e-150) + { + // no need to guard against overflow, so use fast formula + closestPoint = -1; + offset = 1.0; + } + + double s1 = 0.0; + double s2 = 0.0; + for (int i = 0; i < this.points.Count; i++) + { + if (i != closestPoint) + { + double v = offset * this.weights[i] / (t - this.points[i]); + s1 = s1 + (v * this.values[i]); + s2 = s2 + v; + } + else + { + double v = this.weights[i]; + s1 = s1 + (v * this.values[i]); + s2 = s2 + v; + } + } + + return s1 / s2; + } + + /// + /// Differentiate at point t. + /// + /// Point t to interpolate at. + /// Interpolated first derivative at point t. + /// + /// + double IInterpolation.Differentiate(double t) + { + throw new NotSupportedException(); + } + + /// + /// Differentiate at point t. + /// + /// Point t to interpolate at. + /// Interpolated value x(t) + /// Interpolated second derivative at point t. + /// Interpolated first derivative at point t. + /// + /// + double IInterpolation.Differentiate( + double t, + out double interpolatedValue, + out double secondDerivative) + { + throw new NotSupportedException(); + } + + /// + /// Integrate up to point t. + /// + /// Right bound of the integration interval [a,t]. + /// Interpolated definite integral over the interval [a,t]. + /// + double IInterpolation.Integrate(double t) + { + throw new NotSupportedException(); + } + } +} diff --git a/src/Managed/Interpolation/Algorithms/LinearSplineInterpolation.cs b/src/Managed/Interpolation/Algorithms/LinearSplineInterpolation.cs new file mode 100644 index 00000000..ac1b7f9d --- /dev/null +++ b/src/Managed/Interpolation/Algorithms/LinearSplineInterpolation.cs @@ -0,0 +1,177 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 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. +// + +namespace MathNet.Numerics.Interpolation.Algorithms +{ + using System; + using System.Collections.Generic; + + /// + /// Linear Spline Interpolation Algorithm. + /// + /// + /// This algorithm supports both differentiation and integration. + /// + public class LinearSplineInterpolation : IInterpolation + { + /// + /// Internal Spline Interpolation + /// + private readonly SplineInterpolation spline; + + /// + /// Initializes a new instance of the LinearSplineInterpolation class. + /// + public + LinearSplineInterpolation() + { + this.spline = new SplineInterpolation(); + } + + /// + /// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative). + /// + /// + /// + bool IInterpolation.SupportsDifferentiation + { + get { return true; } + } + + /// + /// Gets a value indicating whether the algorithm supports integration (interpolated quadrature). + /// + /// + bool IInterpolation.SupportsIntegration + { + get { return true; } + } + + /// + /// Initialize the interpolation method with the given spline coefficients. + /// + /// Sample Points t + /// Sample Values x(samplePoints) + public + void + Initialize( + IList samplePoints, + IList sampleValues) + { + if (null == samplePoints) + { + throw new ArgumentNullException("samplePoints"); + } + + if (null == sampleValues) + { + throw new ArgumentNullException("sampleValues"); + } + + if (samplePoints.Count < 2) + { + throw new ArgumentOutOfRangeException("samplePoints"); + } + + if (samplePoints.Count != sampleValues.Count) + { + throw new ArgumentException(Properties.Resources.ArgumentVectorsSameLengths); + } + + double[] coefficients = new double[4 * (samplePoints.Count - 1)]; + double[] sortedPoints = new double[samplePoints.Count]; + samplePoints.CopyTo(sortedPoints, 0); + double[] sortedValues = new double[sampleValues.Count]; + sampleValues.CopyTo(sortedValues, 0); + + // TODO: Sorting.Sort(sortedPoints, sortedValues); + + for (int i = 0, j = 0; i < sortedPoints.Length - 1; i++, j += 4) + { + coefficients[j] = sortedValues[i]; + coefficients[j + 1] = (sortedValues[i + 1] - sortedValues[i]) / (sortedPoints[i + 1] - sortedPoints[i]); + coefficients[j + 2] = 0; + coefficients[j + 3] = 0; + } + + this.spline.Initialize(sortedPoints, coefficients); + } + + /// + /// Interpolate at point t. + /// + /// Point t to interpolate at. + /// Interpolated value x(t). + public + double + Interpolate(double t) + { + return this.spline.Interpolate(t); + } + + /// + /// Differentiate at point t. + /// + /// Point t to interpolate at. + /// Interpolated first derivative at point t. + /// + /// + public double Differentiate(double t) + { + return this.spline.Differentiate(t); + } + + /// + /// Differentiate at point t. + /// + /// Point t to interpolate at. + /// Interpolated value x(t) + /// Interpolated second derivative at point t. + /// Interpolated first derivative at point t. + /// + /// + public double Differentiate( + double t, + out double interpolatedValue, + out double secondDerivative) + { + return this.spline.Differentiate(t, out interpolatedValue, out secondDerivative); + } + + /// + /// Integrate up to point t. + /// + /// Right bound of the integration interval [a,t]. + /// Interpolated definite integral over the interval [a,t]. + /// + public double Integrate(double t) + { + return this.spline.Integrate(t); + } + } +} diff --git a/src/Managed/Interpolation/Algorithms/RationalPoleFreeInterpolation.cs b/src/Managed/Interpolation/Algorithms/RationalPoleFreeInterpolation.cs new file mode 100644 index 00000000..6892dfd5 --- /dev/null +++ b/src/Managed/Interpolation/Algorithms/RationalPoleFreeInterpolation.cs @@ -0,0 +1,248 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 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. +// + +namespace MathNet.Numerics.Interpolation.Algorithms +{ + using System; + using System.Collections.Generic; + + /// + /// Barycentric Rational Interpolation without poles, using Floater and Hormann's Algorithm. + /// + /// + /// This algorithm neither supports differentiation nor integration. + /// + public class RationalPoleFreeInterpolation : IInterpolation + { + /// + /// Internal Barycentric Interpolation + /// + private readonly BarycentricInterpolation barycentric; + + /// + /// Initializes a new instance of the RationalPoleFreeInterpolation class. + /// + public + RationalPoleFreeInterpolation() + { + this.barycentric = new BarycentricInterpolation(); + } + + /// + /// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative). + /// + /// + /// + bool IInterpolation.SupportsDifferentiation + { + get { return false; } + } + + /// + /// Gets a value indicating whether the algorithm supports integration (interpolated quadrature). + /// + /// + bool IInterpolation.SupportsIntegration + { + get { return false; } + } + + /// + /// Initialize the interpolation method with the given sample set. + /// + /// + /// The interpolation scheme order will be set to 3. + /// + /// Sample Points t + /// Sample Values x(t) + public + void + Initialize( + IList samplePoints, + IList sampleValues) + { + this.Initialize(samplePoints, sampleValues, Math.Min(3, samplePoints.Count - 1)); + } + + /// + /// Initialize the interpolation method with the given sample set. + /// + /// Sample Points t + /// Sample Values x(t) + /// + /// Order of the interpolation scheme, 0 <= order <= N. + /// In most cases a value between 3 and 8 gives good results. + /// + public + void + Initialize( + IList samplePoints, + IList sampleValues, + int order) + { + if (null == samplePoints) + { + throw new ArgumentNullException("samplePoints"); + } + + if (null == sampleValues) + { + throw new ArgumentNullException("sampleValues"); + } + + if (samplePoints.Count < 1) + { + throw new ArgumentOutOfRangeException("samplePoints"); + } + + if (samplePoints.Count != sampleValues.Count) + { + throw new ArgumentException(Properties.Resources.ArgumentVectorsSameLengths); + } + + if (0 > order || samplePoints.Count <= order) + { + throw new ArgumentOutOfRangeException("order"); + } + + double[] sortedWeights = new double[sampleValues.Count]; + double[] sortedPoints = new double[samplePoints.Count]; + samplePoints.CopyTo(sortedPoints, 0); + + // order: odd -> negative, even -> positive + double sign = ((order & 0x1) == 0x1) ? -1.0 : 1.0; + + // init permutation vector + int[] perm = new int[sortedWeights.Length]; + for (int i = 0; i < perm.Length; i++) + { + perm[i] = i; + } + + // sort and update permutation vector + for (int i = 0; i < perm.Length - 1; i++) + { + for (int j = i + 1; j < perm.Length; j++) + { + if (sortedPoints[j] < sortedPoints[i]) + { + double s = sortedPoints[i]; + sortedPoints[i] = sortedPoints[j]; + sortedPoints[j] = s; + int k = perm[i]; + perm[i] = perm[j]; + perm[j] = k; + } + } + } + + // compute barycentric weights + for (int k = 0; k < sortedWeights.Length; k++) + { + double s = 0; + for (int i = Math.Max(k - order, 0); i <= Math.Min(k, sortedWeights.Length - 1 - order); i++) + { + double v = 1; + for (int j = i; j <= i + order; j++) + { + if (j != k) + { + v = v / Math.Abs(sortedPoints[k] - sortedPoints[j]); + } + } + + s = s + v; + } + + sortedWeights[k] = sign * s; + sign = -sign; + } + + // reorder back to original order, based on the permutation vector. + double[] weights = new double[sortedWeights.Length]; + for (int i = 0; i < weights.Length; i++) + { + weights[perm[i]] = sortedWeights[i]; + } + + this.barycentric.Initialize(samplePoints, sampleValues, weights); + } + + /// + /// Interpolate at point t. + /// + /// Point t to interpolate at. + /// Interpolated value x(t). + public + double + Interpolate(double t) + { + return this.barycentric.Interpolate(t); + } + + /// + /// Differentiate at point t. + /// + /// Point t to interpolate at. + /// Interpolated first derivative at point t. + /// + /// + double IInterpolation.Differentiate(double t) + { + throw new NotSupportedException(); + } + + /// + /// Differentiate at point t. + /// + /// Point t to interpolate at. + /// Interpolated value x(t) + /// Interpolated second derivative at point t. + /// Interpolated first derivative at point t. + /// + /// + double IInterpolation.Differentiate( + double t, + out double interpolatedValue, + out double secondDerivative) + { + throw new NotSupportedException(); + } + + /// + /// Integrate up to point t. + /// + /// Right bound of the integration interval [a,t]. + /// Interpolated definite integral over the interval [a,t]. + /// + double IInterpolation.Integrate(double t) + { + throw new NotSupportedException(); + } + } +} diff --git a/src/Managed/Interpolation/Algorithms/SplineInterpolation.cs b/src/Managed/Interpolation/Algorithms/SplineInterpolation.cs new file mode 100644 index 00000000..dfbb87af --- /dev/null +++ b/src/Managed/Interpolation/Algorithms/SplineInterpolation.cs @@ -0,0 +1,249 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009 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. +// + +namespace MathNet.Numerics.Interpolation.Algorithms +{ + using System; + using System.Collections.Generic; + + /// + /// Third-Degree Spline Interpolation Algorithm. + /// + /// + /// This algorithm supports both differentiation and integration. + /// + public class SplineInterpolation : IInterpolation + { + /// + /// Sample Points t. + /// + private IList points; + + /// + /// Spline Coefficients c(t). + /// + private IList coefficients; + + /// + /// Number of samples. + /// + private int sampleCount; + + /// + /// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative). + /// + /// + /// + bool IInterpolation.SupportsDifferentiation + { + get { return true; } + } + + /// + /// Gets a value indicating whether the algorithm supports integration (interpolated quadrature). + /// + /// + bool IInterpolation.SupportsIntegration + { + get { return true; } + } + + /// + /// Initialize the interpolation method with the given spline coefficients. + /// + /// Sample Points t (length: N) + /// Spline Coefficients (length: 4*(N-1)) + public + void + Initialize( + IList samplePoints, + IList splineCoefficients) + { + if (null == samplePoints) + { + throw new ArgumentNullException("samplePoints"); + } + + if (null == splineCoefficients) + { + throw new ArgumentNullException("splineCoefficients"); + } + + if (samplePoints.Count < 1) + { + throw new ArgumentOutOfRangeException("samplePoints"); + } + + if (splineCoefficients.Count != 4 * (samplePoints.Count - 1)) + { + throw new ArgumentOutOfRangeException("splineCoefficients"); + } + + this.points = samplePoints; + this.coefficients = splineCoefficients; + this.sampleCount = samplePoints.Count; + } + + /// + /// Interpolate at point t. + /// + /// Point t to interpolate at. + /// Interpolated value x(t). + public + double + Interpolate(double t) + { + // Binary search in the [ t[0], ..., t[n-2] ] (t[n-1] is not included) + int low = 0; + int high = this.sampleCount - 1; + while (low != high - 1) + { + int middle = (low + high) / 2; + if (this.points[middle] > t) + { + high = middle; + } + else + { + low = middle; + } + } + + // Interpolation + t = t - this.points[low]; + int k = low << 2; + return this.coefficients[k] + (t * (this.coefficients[k + 1] + (t * (this.coefficients[k + 2] + (t * this.coefficients[k + 3]))))); + } + + /// + /// Differentiate at point t. + /// + /// Point t to interpolate at. + /// Interpolated first derivative at point t. + /// + /// + public double Differentiate(double t) + { + // Binary search in the [ t[0], ..., t[n-2] ] (t[n-1] is not included) + int low = 0; + int high = this.sampleCount - 1; + while (low != high - 1) + { + int middle = (low + high) / 2; + if (this.points[middle] > t) + { + high = middle; + } + else + { + low = middle; + } + } + + // Differentiation + t = t - this.points[low]; + int k = low << 2; + return this.coefficients[k + 1] + (2 * t * this.coefficients[k + 2]) + (3 * t * t * this.coefficients[k + 3]); + } + + /// + /// Differentiate at point t. + /// + /// Point t to interpolate at. + /// Interpolated value x(t) + /// Interpolated second derivative at point t. + /// Interpolated first derivative at point t. + /// + /// + public double Differentiate( + double t, + out double interpolatedValue, + out double secondDerivative) + { + // Binary search in the [ t[0], ..., t[n-2] ] (t[n-1] is not included) + int low = 0; + int high = this.sampleCount - 1; + while (low != high - 1) + { + int middle = (low + high) / 2; + if (this.points[middle] > t) + { + high = middle; + } + else + { + low = middle; + } + } + + // Differentiation + t = t - this.points[low]; + int k = low << 2; + interpolatedValue = this.coefficients[k] + (t * (this.coefficients[k + 1] + (t * (this.coefficients[k + 2] + (t * this.coefficients[k + 3]))))); + secondDerivative = (2 * this.coefficients[k + 2]) + (6 * t * this.coefficients[k + 3]); + return this.coefficients[k + 1] + (2 * t * this.coefficients[k + 2]) + (3 * t * t * this.coefficients[k + 3]); + } + + /// + /// Integrate up to point t. + /// + /// Right bound of the integration interval [a,t]. + /// Interpolated definite integral over the interval [a,t]. + /// + public double Integrate(double t) + { + // Binary search in the [ t[0], ..., t[n-2] ] (t[n-1] is not included) + int low = 0; + int high = this.sampleCount - 1; + while (low != high - 1) + { + int middle = (low + high) / 2; + if (this.points[middle] > t) + { + high = middle; + } + else + { + low = middle; + } + } + + // Integration + double result = 0; + for (int i = 0, j = 0; i < low; i++, j += 4) + { + double w = this.points[i + 1] - this.points[i]; + result += w * (this.coefficients[j] + ((w * (this.coefficients[j + 1] * 0.5)) + (w * ((this.coefficients[j + 2] / 3) + (w * this.coefficients[j + 3] * 0.25))))); + } + + t = t - this.points[low]; + int k = low << 2; + return result + (t * (this.coefficients[k] + ((t * (this.coefficients[k + 1] * 0.5)) + (t * (this.coefficients[k + 2] / 3)) + (t * this.coefficients[k + 3] * 0.25)))); + } + } +} diff --git a/src/Managed/Interpolation/Interpolation.cs b/src/Managed/Interpolation/Interpolation.cs index e57bcee1..424822c6 100644 --- a/src/Managed/Interpolation/Interpolation.cs +++ b/src/Managed/Interpolation/Interpolation.cs @@ -30,6 +30,7 @@ namespace MathNet.Numerics.Interpolation { using System; using System.Collections.Generic; + using Algorithms; /// /// Interpolation Factory. @@ -37,7 +38,7 @@ namespace MathNet.Numerics.Interpolation public static class Interpolation { /// - /// Create a rational pole-free interpolation based on arbitrary points. This is the default interpolation scheme. + /// Creates an interpolation based on arbitrary points. /// /// The sample points t. Supports both lists and arrays. /// The sample point values x(t). Supports both lists and arrays. @@ -50,7 +51,7 @@ namespace MathNet.Numerics.Interpolation IList points, IList values) { - throw new NotImplementedException(); + return CreateRationalPoleFree(points, values); } /// @@ -67,7 +68,9 @@ namespace MathNet.Numerics.Interpolation IList points, IList values) { - throw new NotImplementedException(); + LinearSplineInterpolation method = new LinearSplineInterpolation(); + method.Initialize(points, values); + return method; } /// @@ -84,7 +87,9 @@ namespace MathNet.Numerics.Interpolation IList points, IList values) { - throw new NotImplementedException(); + RationalPoleFreeInterpolation method = new RationalPoleFreeInterpolation(); + method.Initialize(points, values); + return method; } } } diff --git a/src/Managed/Managed.csproj b/src/Managed/Managed.csproj index 8e7e7de3..d076ff04 100644 --- a/src/Managed/Managed.csproj +++ b/src/Managed/Managed.csproj @@ -54,6 +54,10 @@ + + + + diff --git a/src/Native/Native.csproj b/src/Native/Native.csproj index 31e1919c..5b250a92 100644 --- a/src/Native/Native.csproj +++ b/src/Native/Native.csproj @@ -71,6 +71,18 @@ Distributions\IDistribution.cs + + Interpolation\Algorithms\BarycentricInterpolation.cs + + + Interpolation\Algorithms\LinearSplineInterpolation.cs + + + Interpolation\Algorithms\RationalPoleFreeInterpolation.cs + + + Interpolation\Algorithms\SplineInterpolation.cs + Interpolation\IInterpolation.cs