Use InterpolateNaturalSorted as data is sorted
MODIFY
WalletWasabi.Fluent/MathNet/CubicSpline.cs
WalletWasabi.Fluent/MathNet/CubicSpline.cs
// </copyright>
using System;
-using System.Collections.Generic;
-using System.Linq;
namespace MathNet.Numerics.Interpolation
{
@@ -104,216 +102,6 @@ namespace MathNet.Numerics.Interpolation
return new CubicSpline(x, c0, c1, c2, c3);
}
- /// <summary>
- /// Create a Hermite cubic spline interpolation from an unsorted set of (x,y) value pairs and their slope (first derivative).
- /// WARNING: Works in-place and can thus causes the data array to be reordered.
- /// </summary>
- public static CubicSpline InterpolateHermiteInplace(double[] x, double[] y, double[] firstDerivatives)
- {
- if (x.Length != y.Length || x.Length != firstDerivatives.Length)
- {
- throw new ArgumentException("All vectors must have the same dimensionality.");
- }
-
- if (x.Length < 2)
- {
- throw new ArgumentException("The given array is too small. It must be at least 2 long.", nameof(x));
- }
-
- Sorting.Sort(x, y, firstDerivatives);
- return InterpolateHermiteSorted(x, y, firstDerivatives);
- }
-
- /// <summary>
- /// Create a Hermite cubic spline interpolation from an unsorted set of (x,y) value pairs and their slope (first derivative).
- /// </summary>
- public static CubicSpline InterpolateHermite(IEnumerable<double> x, IEnumerable<double> y, IEnumerable<double> firstDerivatives)
- {
- // note: we must make a copy, even if the input was arrays already
- return InterpolateHermiteInplace(x.ToArray(), y.ToArray(), firstDerivatives.ToArray());
- }
-
- /// <summary>
- /// Create an Akima cubic spline interpolation from a set of (x,y) value pairs, sorted ascendingly by x.
- /// Akima splines are robust to outliers.
- /// </summary>
- public static CubicSpline InterpolateAkimaSorted(double[] x, double[] y)
- {
- if (x.Length != y.Length)
- {
- throw new ArgumentException("All vectors must have the same dimensionality.");
- }
-
- if (x.Length < 5)
- {
- throw new ArgumentException("The given array is too small. It must be at least 5 long.", nameof(x));
- }
-
- /* Prepare divided differences (diff) and weights (w) */
-
- var diff = new double[x.Length - 1];
- var weights = new double[x.Length - 1];
-
- for (int i = 0; i < diff.Length; i++)
- {
- diff[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
- }
-
- for (int i = 1; i < weights.Length; i++)
- {
- weights[i] = Math.Abs(diff[i] - diff[i - 1]);
- }
-
- /* Prepare Hermite interpolation scheme */
-
- var dd = new double[x.Length];
-
- for (int i = 2; i < dd.Length - 2; i++)
- {
- dd[i] = weights[i - 1].AlmostEqual(0.0) && weights[i + 1].AlmostEqual(0.0)
- ? (((x[i + 1] - x[i]) * diff[i - 1]) + ((x[i] - x[i - 1]) * diff[i])) / (x[i + 1] - x[i - 1])
- : ((weights[i + 1] * diff[i - 1]) + (weights[i - 1] * diff[i])) / (weights[i + 1] + weights[i - 1]);
- }
-
- dd[0] = DifferentiateThreePoint(x, y, 0, 0, 1, 2);
- dd[1] = DifferentiateThreePoint(x, y, 1, 0, 1, 2);
- dd[x.Length - 2] = DifferentiateThreePoint(x, y, x.Length - 2, x.Length - 3, x.Length - 2, x.Length - 1);
- dd[x.Length - 1] = DifferentiateThreePoint(x, y, x.Length - 1, x.Length - 3, x.Length - 2, x.Length - 1);
-
- /* Build Akima spline using Hermite interpolation scheme */
-
- return InterpolateHermiteSorted(x, y, dd);
- }
-
- /// <summary>
- /// Create an Akima cubic spline interpolation from an unsorted set of (x,y) value pairs.
- /// Akima splines are robust to outliers.
- /// WARNING: Works in-place and can thus causes the data array to be reordered.
- /// </summary>
- public static CubicSpline InterpolateAkimaInplace(double[] x, double[] y)
- {
- if (x.Length != y.Length)
- {
- throw new ArgumentException("All vectors must have the same dimensionality.");
- }
-
- Sorting.Sort(x, y);
- return InterpolateAkimaSorted(x, y);
- }
-
- /// <summary>
- /// Create an Akima cubic spline interpolation from an unsorted set of (x,y) value pairs.
- /// Akima splines are robust to outliers.
- /// </summary>
- public static CubicSpline InterpolateAkima(IEnumerable<double> x, IEnumerable<double> y)
- {
- // note: we must make a copy, even if the input was arrays already
- return InterpolateAkimaInplace(x.ToArray(), y.ToArray());
- }
-
- /// <summary>
- /// Create a piecewise cubic Hermite interpolating polynomial from an unsorted set of (x,y) value pairs.
- /// Monotone-preserving interpolation with continuous first derivative.
- /// </summary>
- public static CubicSpline InterpolatePchipSorted(double[] x, double[] y)
- {
- // Implementation based on "Numerical Computing with Matlab" (Moler, 2004).
-
- if (x.Length != y.Length)
- {
- throw new ArgumentException("All vectors must have the same dimensionality.");
- }
-
- if (x.Length < 3)
- {
- throw new ArgumentException("The given array is too small. It must be at least 3 long.", nameof(x));
- }
-
- var m = new double[x.Length - 1];
-
- for (int i = 0; i < m.Length; i++)
- {
- m[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
- }
-
- var dd = new double[x.Length];
- var hPrev = x[1] - x[0];
-
- // This check is quite costly as it usually involves a Math.Pow().
- var mPrevIs0 = m[0].AlmostEqual(0.0);
-
- for (var i = 1; i < x.Length - 1; ++i)
- {
- var h = x[i + 1] - x[i];
- var mIs0 = m[i].AlmostEqual(0.0);
-
- if (mIs0 || mPrevIs0 || Math.Sign(m[i]) != Math.Sign(m[i - 1]))
- {
- dd[i] = 0;
- }
- else
- {
- // Weighted harmonic mean of each slope.
- var w1 = 2 * h + hPrev;
- var w2 = h + 2 * hPrev;
- dd[i] = (w1 + w2) / (w1 / m[i - 1] + w2 / m[i]);
- }
-
- hPrev = h;
- mPrevIs0 = mIs0;
- }
-
- // Special case end-points.
- dd[0] = PchipEndPoints(x[1] - x[0], x[2] - x[1], m[0], m[1]);
- dd[dd.Length - 1] = PchipEndPoints(x[x.Length - 1] - x[x.Length - 2], x[x.Length - 2] - x[x.Length - 3], m[m.Length - 1], m[m.Length - 2]);
-
- return InterpolateHermiteSorted(x, y, dd);
- }
-
- private static double PchipEndPoints(double h0, double h1, double m0, double m1)
- {
- // One-sided, shape-preserving, three-point estimate for the derivative.
- var d = ((2 * h0 + h1) * m0 - h0 * m1) / (h0 + h1);
-
- if (Math.Sign(d) != Math.Sign(m0))
- {
- return 0.0;
- }
-
- if (Math.Sign(m0) != Math.Sign(m1) && (Math.Abs(d) > 3 * Math.Abs(m0)))
- {
- return 3 * m0;
- }
-
- return d;
- }
-
- /// <summary>
- /// Create a piecewise cubic Hermite interpolating polynomial from an unsorted set of (x,y) value pairs.
- /// Monotone-preserving interpolation with continuous first derivative.
- /// WARNING: Works in-place and can thus causes the data array to be reordered.
- /// </summary>
- public static CubicSpline InterpolatePchipInplace(double[] x, double[] y)
- {
- if (x.Length != y.Length)
- {
- throw new ArgumentException("All vectors must have the same dimensionality.");
- }
-
- Sorting.Sort(x, y);
- return InterpolatePchipSorted(x, y);
- }
-
- /// <summary>
- /// Create a piecewise cubic Hermite interpolating polynomial from an unsorted set of (x,y) value pairs.
- /// Monotone-preserving interpolation with continuous first derivative.
- /// </summary>
- public static CubicSpline InterpolatePchip(IEnumerable<double> x, IEnumerable<double> y)
- {
- // note: we must make a copy, even if the input was arrays already
- return InterpolatePchipInplace(x.ToArray(), y.ToArray());
- }
-
/// <summary>
/// Create a cubic spline interpolation from a set of (x,y) value pairs, sorted ascendingly by x,
/// and custom boundary/termination conditions.
@@ -426,30 +214,6 @@ namespace MathNet.Numerics.Interpolation
return InterpolateHermiteSorted(x, y, dd);
}
- /// <summary>
- /// Create a cubic spline interpolation from an unsorted set of (x,y) value pairs and custom boundary/termination conditions.
- /// WARNING: Works in-place and can thus causes the data array to be reordered.
- /// </summary>
- public static CubicSpline InterpolateBoundariesInplace(double[] x, double[] y, SplineBoundaryCondition leftBoundaryCondition, double leftBoundary, SplineBoundaryCondition rightBoundaryCondition, double rightBoundary)
- {
- if (x.Length != y.Length)
- {
- throw new ArgumentException("All vectors must have the same dimensionality.");
- }
-
- Sorting.Sort(x, y);
- return InterpolateBoundariesSorted(x, y, leftBoundaryCondition, leftBoundary, rightBoundaryCondition, rightBoundary);
- }
-
- /// <summary>
- /// Create a cubic spline interpolation from an unsorted set of (x,y) value pairs and custom boundary/termination conditions.
- /// </summary>
- public static CubicSpline InterpolateBoundaries(IEnumerable<double> x, IEnumerable<double> y, SplineBoundaryCondition leftBoundaryCondition, double leftBoundary, SplineBoundaryCondition rightBoundaryCondition, double rightBoundary)
- {
- // note: we must make a copy, even if the input was arrays already
- return InterpolateBoundariesInplace(x.ToArray(), y.ToArray(), leftBoundaryCondition, leftBoundary, rightBoundaryCondition, rightBoundary);
- }
-
/// <summary>
/// Create a natural cubic spline interpolation from a set of (x,y) value pairs
/// and zero second derivatives at the two boundaries, sorted ascendingly by x.
@@ -459,50 +223,6 @@ namespace MathNet.Numerics.Interpolation
return InterpolateBoundariesSorted(x, y, SplineBoundaryCondition.SecondDerivative, 0.0, SplineBoundaryCondition.SecondDerivative, 0.0);
}
- /// <summary>
- /// Create a natural cubic spline interpolation from an unsorted set of (x,y) value pairs
- /// and zero second derivatives at the two boundaries.
- /// WARNING: Works in-place and can thus causes the data array to be reordered.
- /// </summary>
- public static CubicSpline InterpolateNaturalInplace(double[] x, double[] y)
- {
- return InterpolateBoundariesInplace(x, y, SplineBoundaryCondition.SecondDerivative, 0.0, SplineBoundaryCondition.SecondDerivative, 0.0);
- }
-
- /// <summary>
- /// Create a natural cubic spline interpolation from an unsorted set of (x,y) value pairs
- /// and zero second derivatives at the two boundaries.
- /// </summary>
- public static CubicSpline InterpolateNatural(IEnumerable<double> x, IEnumerable<double> y)
- {
- return InterpolateBoundaries(x, y, SplineBoundaryCondition.SecondDerivative, 0.0, SplineBoundaryCondition.SecondDerivative, 0.0);
- }
-
- /// <summary>
- /// Three-Point Differentiation Helper.
- /// </summary>
- /// <param name="xx">Sample Points t.</param>
- /// <param name="yy">Sample Values x(t).</param>
- /// <param name="indexT">Index of the point of the differentiation.</param>
- /// <param name="index0">Index of the first sample.</param>
- /// <param name="index1">Index of the second sample.</param>
- /// <param name="index2">Index of the third sample.</param>
- /// <returns>The derivative approximation.</returns>
- private static double DifferentiateThreePoint(double[] xx, double[] yy, int indexT, int index0, int index1, int index2)
- {
- double x0 = yy[index0];
- double x1 = yy[index1];
- double x2 = yy[index2];
-
- double t = xx[indexT] - xx[index0];
- double t1 = xx[index1] - xx[index0];
- double t2 = xx[index2] - xx[index0];
-
- double a = (x2 - x0 - (t2 / t1 * (x1 - x0))) / (t2 * (t2 - t1));
- double b = (x1 - x0 - a * t1 * t1) / t1;
- return (2 * a * t) + b;
- }
-
/// <summary>
/// Tridiagonal Solve Helper.
/// </summary>
@@ -542,51 +262,6 @@ namespace MathNet.Numerics.Interpolation
return _c0[k] + x * (_c1[k] + x * (_c2[k] + x * _c3[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 = LeftSegmentIndex(t);
- var x = t - _x[k];
- return _c1[k] + x * (2 * _c2[k] + x * 3 * _c3[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 = LeftSegmentIndex(t);
- var x = t - _x[k];
- return 2 * _c2[k] + x * 6 * _c3[k];
- }
-
- /// <summary>
- /// Indefinite integral at point t.
- /// </summary>
- /// <param name="t">Point t to integrate at.</param>
- public double Integrate(double t)
- {
- int k = LeftSegmentIndex(t);
- var x = t - _x[k];
- return _indefiniteIntegral.Value[k] + x * (_c0[k] + x * (_c1[k] / 2 + x * (_c2[k] / 3 + x * _c3[k] / 4)));
- }
-
- /// <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);
- }
-
private double[] ComputeIndefiniteIntegral()
{
var integral = new double[_c1.Length];

MODIFY
WalletWasabi.Fluent/ViewModels/TestLineChartViewModel.cs
WalletWasabi.Fluent/ViewModels/TestLineChartViewModel.cs
@@ -65,7 +65,7 @@ namespace WalletWasabi.Fluent.ViewModels
y[i] = Math.Sin(x[i]);
Console.WriteLine($"{x[i]} {y[i]}");
}
- var spline = CubicSpline.InterpolateNatural(x, y);
+ var spline = CubicSpline.InterpolateNaturalSorted(x, y);
double testX = 1.5 * Math.PI / 12;
double testY = spline.Interpolate(testX);
Console.WriteLine($"Interpolated:");

Use InterpolateNaturalSorted as data is sorted