Restore original CubicSpline code
MODIFY
WalletWasabi.Fluent/MathNet/CubicSpline.cs
WalletWasabi.Fluent/MathNet/CubicSpline.cs
@@ -7,23 +7,24 @@ namespace WalletWasabi.Fluent.MathNet
/// <summary>
/// Cubic Spline Interpolation.
/// </summary>
+ /// <remarks>Supports both differentiation and integration.</remarks>
public class CubicSpline
{
- private readonly double[] _x;
- private readonly double[] _c0;
- private readonly double[] _c1;
- private readonly double[] _c2;
- private readonly double[] _c3;
+ readonly double[] _x;
+ readonly double[] _c0;
+ readonly double[] _c1;
+ readonly double[] _c2;
+ readonly double[] _c3;
+ readonly Lazy<double[]> _indefiniteIntegral;
- /// <param name="x">Sample points (N+1), sorted ascending</param>
+ /// <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>
- /// <param name="c3">Third order spline coefficients (N)</param>
- private CubicSpline(double[] x, double[] c0, double[] c1, double[] c2, double[] c3)
+ /// <param name="c2">second order spline coefficients (N)</param>
+ /// <param name="c3">third order spline coefficients (N)</param>
+ public CubicSpline(double[] x, double[] c0, double[] c1, double[] c2, double[] c3)
{
- if (x.Length != c0.Length + 1 || x.Length != c1.Length + 1 || x.Length != c2.Length + 1 ||
- x.Length != c3.Length + 1)
+ if (x.Length != c0.Length + 1 || x.Length != c1.Length + 1 || x.Length != c2.Length + 1 || x.Length != c3.Length + 1)
{
throw new ArgumentException("All vectors must have the same dimensionality.");
}
@@ -38,38 +39,13 @@ namespace WalletWasabi.Fluent.MathNet
_c1 = c1;
_c2 = c2;
_c3 = c3;
- }
-
- /// <summary>
- /// Left and right boundary conditions.
- /// </summary>
- private enum SplineBoundaryCondition
- {
- /// <summary>
- /// Natural Boundary (Zero second derivative).
- /// </summary>
- Natural = 0,
-
- /// <summary>
- /// Parabolically Terminated boundary.
- /// </summary>
- ParabolicallyTerminated,
-
- /// <summary>
- /// Fixed first derivative at the boundary.
- /// </summary>
- FirstDerivative,
-
- /// <summary>
- /// Fixed second derivative at the boundary.
- /// </summary>
- SecondDerivative
+ _indefiniteIntegral = new Lazy<double[]>(ComputeIndefiniteIntegral);
}
/// <summary>
/// Create a Hermite cubic spline interpolation from a set of (x,y) value pairs and their slope (first derivative), sorted ascendingly by x.
/// </summary>
- private static CubicSpline InterpolateHermiteSorted(double[] x, double[] y, double[] firstDerivatives)
+ public static CubicSpline InterpolateHermiteSorted(double[] x, double[] y, double[] firstDerivatives)
{
if (x.Length != y.Length || x.Length != firstDerivatives.Length)
{
@@ -98,11 +74,143 @@ namespace WalletWasabi.Fluent.MathNet
return new CubicSpline(x, c0, c1, c2, c3);
}
+ /// <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 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);
+ }
+
+ 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 cubic spline interpolation from a set of (x,y) value pairs, sorted ascendingly by x,
/// and custom boundary/termination conditions.
/// </summary>
- private static CubicSpline InterpolateBoundariesSorted(double[] x, double[] y, SplineBoundaryCondition leftBoundaryCondition, double leftBoundary, SplineBoundaryCondition rightBoundaryCondition, double rightBoundary)
+ public static CubicSpline InterpolateBoundariesSorted(double[] x, double[] y,
+ SplineBoundaryCondition leftBoundaryCondition, double leftBoundary,
+ SplineBoundaryCondition rightBoundaryCondition, double rightBoundary)
{
if (x.Length != y.Length)
{
@@ -175,8 +283,7 @@ namespace WalletWasabi.Fluent.MathNet
a1[i] = x[i + 1] - x[i];
a2[i] = 2*(x[i + 1] - x[i - 1]);
a3[i] = x[i] - x[i - 1];
- b[i] = (3 * (y[i] - y[i - 1]) / (x[i] - x[i - 1]) * (x[i + 1] - x[i])) +
- (3 * (y[i + 1] - y[i]) / (x[i + 1] - x[i]) * (x[i] - x[i - 1]));
+ b[i] = (3*(y[i] - y[i - 1])/(x[i] - x[i - 1])*(x[i + 1] - x[i])) + (3*(y[i + 1] - y[i])/(x[i + 1] - x[i])*(x[i] - x[i - 1]));
}
// Right Boundary
@@ -198,8 +305,7 @@ namespace WalletWasabi.Fluent.MathNet
a1[n - 1] = 1;
a2[n - 1] = 2;
a3[n - 1] = 0;
- b[n - 1] = (3 * (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2])) +
- (0.5 * rightBoundary * (x[n - 1] - x[n - 2]));
+ b[n - 1] = (3*(y[n - 1] - y[n - 2])/(x[n - 1] - x[n - 2])) + (0.5*rightBoundary*(x[n - 1] - x[n - 2]));
break;
default:
throw new NotSupportedException("Invalid Right Boundary Condition.");
@@ -219,6 +325,31 @@ namespace WalletWasabi.Fluent.MathNet
return InterpolateBoundariesSorted(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>
+ 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>
@@ -227,7 +358,7 @@ namespace WalletWasabi.Fluent.MathNet
/// <param name="c">The c-vector[n].</param>
/// <param name="d">The d-vector[n], will be modified by this function.</param>
/// <returns>The x-vector[n]</returns>
- private static double[] SolveTridiagonal(double[] a, double[] b, double[] c, double[] d)
+ static double[] SolveTridiagonal(double[] a, double[] b, double[] c, double[] d)
{
for (int k = 1; k < a.Length; k++)
{
@@ -258,11 +389,68 @@ namespace WalletWasabi.Fluent.MathNet
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);
+ }
+
+ 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 + w*_c3[i]/4)));
+ }
+
+ return integral;
+ }
+
/// <summary>
/// Find the index of the greatest sample point smaller than t,
/// or the left index of the closest segment for extrapolation.
/// </summary>
- private int LeftSegmentIndex(double t)
+ int LeftSegmentIndex(double t)
{
int index = Array.BinarySearch(_x, t);
if (index < 0)

ADD
WalletWasabi.Fluent/MathNet/Precision.Equality.cs
+// Copyright (c) 2009-2013 Math.NET Taken from http://github.com/mathnet/mathnet-numerics and modified for Wasabi Wallet
+
+using System;
+
+namespace WalletWasabi.Fluent.MathNet
+{
+ public static partial class Precision
+ {
+ /// <summary>
+ /// Compares two doubles and determines if they are equal
+ /// within the specified maximum absolute error.
+ /// </summary>
+ /// <param name="a">The norm of the first value (can be negative).</param>
+ /// <param name="b">The norm of the second value (can be negative).</param>
+ /// <param name="diff">The norm of the difference of the two values (can be negative).</param>
+ /// <param name="maximumAbsoluteError">The absolute accuracy required for being almost equal.</param>
+ /// <returns>True if both doubles are almost equal up to the specified maximum absolute error, false otherwise.</returns>
+ public static bool AlmostEqualNorm(this double a, double b, double diff, double maximumAbsoluteError)
+ {
+ // If A or B are infinity (positive or negative) then
+ // only return true if they are exactly equal to each other -
+ // that is, if they are both infinities of the same sign.
+ if (double.IsInfinity(a) || double.IsInfinity(b))
+ {
+ return a == b;
+ }
+
+ // If A or B are a NAN, return false. NANs are equal to nothing,
+ // not even themselves.
+ if (double.IsNaN(a) || double.IsNaN(b))
+ {
+ return false;
+ }
+
+ return Math.Abs(diff) < maximumAbsoluteError;
+ }
+
+ /// <summary>
+ /// Checks whether two real numbers are almost equal.
+ /// </summary>
+ /// <param name="a">The first number</param>
+ /// <param name="b">The second number</param>
+ /// <returns>true if the two values differ by no more than 10 * 2^(-52); false otherwise.</returns>
+ public static bool AlmostEqual(this double a, double b)
+ {
+ return AlmostEqualNorm(a, b, a - b, DefaultDoubleAccuracy);
+ }
+ }
+}
\ No newline at end of file

ADD
WalletWasabi.Fluent/MathNet/Precision.cs
+// Copyright (c) 2009-2015 Math.NET Taken from http://github.com/mathnet/mathnet-numerics and modified for Wasabi Wallet
+
+using System;
+
+namespace WalletWasabi.Fluent.MathNet
+{
+ /// <summary>
+ /// Utilities for working with floating point numbers.
+ /// </summary>
+ /// <remarks>
+ /// <para>
+ /// Useful links:
+ /// <list type="bullet">
+ /// <item>
+ /// http://docs.sun.com/source/806-3568/ncg_goldberg.html#689 - What every computer scientist should know about floating-point arithmetic
+ /// </item>
+ /// <item>
+ /// http://en.wikipedia.org/wiki/Machine_epsilon - Gives the definition of machine epsilon
+ /// </item>
+ /// </list>
+ /// </para>
+ /// </remarks>
+ public static partial class Precision
+ {
+ /// <summary>
+ /// The number of binary digits used to represent the binary number for a double precision floating
+ /// point value. i.e. there are this many digits used to represent the
+ /// actual number, where in a number as: 0.134556 * 10^5 the digits are 0.134556 and the exponent is 5.
+ /// </summary>
+ private const int DoubleWidth = 53;
+
+ /// <summary>
+ /// Standard epsilon, the maximum relative precision of IEEE 754 double-precision floating numbers (64 bit).
+ /// According to the definition of Prof. Demmel and used in LAPACK and Scilab.
+ /// </summary>
+ public static readonly double DoublePrecision = Math.Pow(2, -DoubleWidth);
+
+ /// <summary>
+ /// Value representing 10 * 2^(-53) = 1.11022302462516E-15
+ /// </summary>
+ private static readonly double DefaultDoubleAccuracy = DoublePrecision * 10;
+ }
+}
\ No newline at end of file

Restore original CubicSpline code