Christoph Ruegg
13 years ago
5 changed files with 193 additions and 214 deletions
@ -0,0 +1,182 @@ |
|||
// <copyright file="LinearSpline.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; |
|||
using System.Collections.Generic; |
|||
using System.Linq; |
|||
using MathNet.Numerics.Properties; |
|||
|
|||
namespace MathNet.Numerics.Interpolation |
|||
{ |
|||
/// <summary>
|
|||
/// Linear Spline Interpolation.
|
|||
/// </summary>
|
|||
/// <remarks>Supports both differentiation and integration.</remarks>
|
|||
public class LinearSpline : IInterpolation |
|||
{ |
|||
readonly double[] _x; |
|||
readonly double[] _c0; |
|||
readonly double[] _c1; |
|||
|
|||
/// <param name="x">Sample points (N+1), sorted ascending</param>
|
|||
/// <param name="c0">Sample values (N or N+1) at the corresponding points; intercept, zero order coefficients</param>
|
|||
/// <param name="c1">Slopes (N) at the sample points (first order coefficients): N</param>
|
|||
public LinearSpline(double[] x, double[] c0, double[] c1) |
|||
{ |
|||
_x = x; |
|||
_c0 = c0; |
|||
_c1 = c1; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Create a linear spline interpolation from a set of (x,y) value pairs.
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// 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.
|
|||
/// </remarks>
|
|||
public static LinearSpline Interpolate(IEnumerable<double> x, IEnumerable<double> y) |
|||
{ |
|||
var xx = (x as double[]) ?? x.ToArray(); |
|||
var yy = (y as double[]) ?? y.ToArray(); |
|||
|
|||
if (xx.Length != yy.Length) |
|||
{ |
|||
throw new ArgumentException(Resources.ArgumentVectorsSameLength); |
|||
} |
|||
|
|||
Sorting.Sort(xx, yy); |
|||
|
|||
var c1 = new double[xx.Length - 1]; |
|||
for (int i = 0; i < c1.Length; i++) |
|||
{ |
|||
c1[i] = (yy[i + 1] - yy[i])/(xx[i + 1] - xx[i]); |
|||
} |
|||
|
|||
return new LinearSpline(xx, yy, c1); |
|||
} |
|||
|
|||
/// <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); |
|||
return _c0[k] + (t - _x[k])*_c1[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]; |
|||
} |
|||
|
|||
/// <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) |
|||
{ |
|||
return 0d; |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Integrate up to point t.
|
|||
/// </summary>
|
|||
/// <param name="t">Right bound of the integration interval [a,t].</param>
|
|||
/// <returns>Interpolated definite integral over the interval [a,t].</returns>
|
|||
public double Integrate(double t) |
|||
{ |
|||
int k = LeftBracketIndex(t); |
|||
var x = (t - _x[k]); |
|||
return x*(_c0[k] + x*0.5*_c1[k]); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Interpolate, differentiate and 2nd differentiate at point t.
|
|||
/// </summary>
|
|||
/// <param name="t">Point t to interpolate at.</param>
|
|||
/// <returns>Interpolated first derivative at point t.</returns>
|
|||
public Tuple<double, double, double> DifferentiateAll(double t) |
|||
{ |
|||
int k = LeftBracketIndex(t); |
|||
var x = (t - _x[k]); |
|||
return new Tuple<double, double, double>(_c0[k] + x*_c1[k], _c1[k], 0d); |
|||
} |
|||
|
|||
/// <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; |
|||
} |
|||
} |
|||
} |
|||
@ -1,188 +0,0 @@ |
|||
// <copyright file="LinearSplineInterpolation.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-2010 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; |
|||
using System.Collections.Generic; |
|||
using MathNet.Numerics.Properties; |
|||
|
|||
namespace MathNet.Numerics.Interpolation |
|||
{ |
|||
/// <summary>
|
|||
/// Linear Spline Interpolation Algorithm.
|
|||
/// </summary>
|
|||
/// <remarks>
|
|||
/// This algorithm supports both differentiation and integration.
|
|||
/// </remarks>
|
|||
public class LinearSplineInterpolation : IInterpolation |
|||
{ |
|||
/// <summary>
|
|||
/// Internal Spline Interpolation
|
|||
/// </summary>
|
|||
readonly SplineInterpolation _spline; |
|||
|
|||
/// <summary>
|
|||
/// Initializes a new instance of the LinearSplineInterpolation class.
|
|||
/// </summary>
|
|||
public LinearSplineInterpolation() |
|||
{ |
|||
_spline = new SplineInterpolation(); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Initializes a new instance of the LinearSplineInterpolation class.
|
|||
/// </summary>
|
|||
/// <param name="samplePoints">Sample Points t, sorted ascending.</param>
|
|||
/// <param name="sampleValues">Sample Values x(t)</param>
|
|||
public LinearSplineInterpolation(IList<double> samplePoints, IList<double> sampleValues) |
|||
{ |
|||
_spline = new SplineInterpolation(); |
|||
Initialize(samplePoints, sampleValues); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets a value indicating whether the algorithm supports differentiation (interpolated derivative).
|
|||
/// </summary>
|
|||
/// <seealso cref="Differentiate(double)"/>
|
|||
/// <seealso cref="DifferentiateAll(double)"/>
|
|||
bool IInterpolation.SupportsDifferentiation |
|||
{ |
|||
get { return true; } |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Gets a value indicating whether the algorithm supports integration (interpolated quadrature).
|
|||
/// </summary>
|
|||
/// <seealso cref="Integrate"/>
|
|||
bool IInterpolation.SupportsIntegration |
|||
{ |
|||
get { return true; } |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Initialize the interpolation method with the given spline coefficients (sorted by the sample points t).
|
|||
/// </summary>
|
|||
/// <param name="samplePoints">Sample Points t, sorted ascending.</param>
|
|||
/// <param name="sampleValues">Sample Values x(t)</param>
|
|||
public void Initialize(IList<double> samplePoints, IList<double> sampleValues) |
|||
{ |
|||
double[] coefficients = EvaluateSplineCoefficients(samplePoints, sampleValues); |
|||
_spline.Initialize(samplePoints, coefficients); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Evaluate the spline coefficients as used
|
|||
/// internally by this interpolation algorithm.
|
|||
/// </summary>
|
|||
/// <param name="samplePoints">Sample Points t, sorted ascending.</param>
|
|||
/// <param name="sampleValues">Sample Values x(t)</param>
|
|||
/// <returns>Spline Coefficient Vector</returns>
|
|||
public static double[] EvaluateSplineCoefficients(IList<double> samplePoints, IList<double> 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(Resources.ArgumentVectorsSameLength); |
|||
} |
|||
|
|||
for (var i = 1; i < samplePoints.Count; ++i) |
|||
if (samplePoints[i] <= samplePoints[i - 1]) |
|||
throw new ArgumentException(Resources.Interpolation_Initialize_SamplePointsNotStrictlyAscendingOrder, "samplePoints"); |
|||
|
|||
var coefficients = new double[4*(samplePoints.Count - 1)]; |
|||
for (int i = 0, j = 0; i < samplePoints.Count - 1; i++, j += 4) |
|||
{ |
|||
coefficients[j] = sampleValues[i]; |
|||
coefficients[j + 1] = (sampleValues[i + 1] - sampleValues[i])/(samplePoints[i + 1] - samplePoints[i]); |
|||
coefficients[j + 2] = 0; |
|||
coefficients[j + 3] = 0; |
|||
} |
|||
return coefficients; |
|||
} |
|||
|
|||
/// <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) |
|||
{ |
|||
return _spline.Interpolate(t); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Differentiate at point t.
|
|||
/// </summary>
|
|||
/// <param name="t">Point t to interpolate at.</param>
|
|||
/// <returns>Interpolated first derivative at point t.</returns>
|
|||
/// <seealso cref="IInterpolation.SupportsDifferentiation"/>
|
|||
/// <seealso cref="DifferentiateAll(double)"/>
|
|||
public double Differentiate(double t) |
|||
{ |
|||
return _spline.Differentiate(t); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Interpolate, differentiate and 2nd differentiate at point t.
|
|||
/// </summary>
|
|||
/// <param name="t">Point t to interpolate at.</param>
|
|||
/// <returns>Interpolated first derivative at point t.</returns>
|
|||
/// <seealso cref="IInterpolation.SupportsDifferentiation"/>
|
|||
/// <seealso cref="Differentiate(double)"/>
|
|||
public Tuple<double, double, double> DifferentiateAll(double t) |
|||
{ |
|||
return _spline.DifferentiateAll(t); |
|||
} |
|||
|
|||
/// <summary>
|
|||
/// Integrate up to point t.
|
|||
/// </summary>
|
|||
/// <param name="t">Right bound of the integration interval [a,t].</param>
|
|||
/// <returns>Interpolated definite integral over the interval [a,t].</returns>
|
|||
/// <seealso cref="IInterpolation.SupportsIntegration"/>
|
|||
public double Integrate(double t) |
|||
{ |
|||
return _spline.Integrate(t); |
|||
} |
|||
} |
|||
} |
|||
Loading…
Reference in new issue