From 5dca1ce266822d1396f6bd10232d2f353c224515 Mon Sep 17 00:00:00 2001 From: Christoph Ruegg Date: Tue, 7 Jul 2009 01:35:24 +0800 Subject: [PATCH] interpolation: spline refactoring (code review) Signed-off-by: Christoph Ruegg --- .../Algorithms/LinearSplineInterpolation.cs | 2 +- .../Algorithms/SplineInterpolation.cs | 128 +++++++++--------- 2 files changed, 62 insertions(+), 68 deletions(-) diff --git a/src/Managed/Interpolation/Algorithms/LinearSplineInterpolation.cs b/src/Managed/Interpolation/Algorithms/LinearSplineInterpolation.cs index e14f7475..187444b9 100644 --- a/src/Managed/Interpolation/Algorithms/LinearSplineInterpolation.cs +++ b/src/Managed/Interpolation/Algorithms/LinearSplineInterpolation.cs @@ -119,7 +119,7 @@ namespace MathNet.Numerics.Interpolation.Algorithms double[] sortedValues = new double[sampleValues.Count]; sampleValues.CopyTo(sortedValues, 0); - // TODO: Sorting.Sort(sortedPoints, sortedValues); + /* TODO: Sorting.Sort(sortedPoints, sortedValues); */ for (int i = 0, j = 0; i < sortedPoints.Length - 1; i++, j += 4) { diff --git a/src/Managed/Interpolation/Algorithms/SplineInterpolation.cs b/src/Managed/Interpolation/Algorithms/SplineInterpolation.cs index f6801878..586b1fee 100644 --- a/src/Managed/Interpolation/Algorithms/SplineInterpolation.cs +++ b/src/Managed/Interpolation/Algorithms/SplineInterpolation.cs @@ -133,26 +133,16 @@ namespace MathNet.Numerics.Interpolation.Algorithms /// 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; - } - } + int closestLeftIndex = IndexOfClosestPointLeftOf(t); // 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]))))); + double offset = t - this.points[closestLeftIndex]; + int k = closestLeftIndex << 2; + + return this.coefficients[k] + + (offset * (this.coefficients[k + 1] + + (offset * (this.coefficients[k + 2] + + (offset * this.coefficients[k + 3]))))); } /// @@ -164,26 +154,15 @@ namespace MathNet.Numerics.Interpolation.Algorithms /// 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; - } - } + int closestLeftIndex = IndexOfClosestPointLeftOf(t); // 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]); + double offset = t - this.points[closestLeftIndex]; + int k = closestLeftIndex << 2; + + return this.coefficients[k + 1] + + (2 * offset * this.coefficients[k + 2]) + + (3 * offset * offset * this.coefficients[k + 3]); } /// @@ -200,28 +179,23 @@ namespace MathNet.Numerics.Interpolation.Algorithms 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; - } - } + int closestLeftIndex = IndexOfClosestPointLeftOf(t); // 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]); + double offset = t - this.points[closestLeftIndex]; + int k = closestLeftIndex << 2; + + interpolatedValue = this.coefficients[k] + + (offset * (this.coefficients[k + 1] + + (offset * (this.coefficients[k + 2] + + (offset * this.coefficients[k + 3]))))); + + secondDerivative = (2 * this.coefficients[k + 2]) + + (6 * offset * this.coefficients[k + 3]); + + return this.coefficients[k + 1] + + (2 * offset * this.coefficients[k + 2]) + + (3 * offset * offset * this.coefficients[k + 3]); } /// @@ -231,6 +205,36 @@ namespace MathNet.Numerics.Interpolation.Algorithms /// Interpolated definite integral over the interval [a,t]. /// public double Integrate(double t) + { + int closestLeftIndex = IndexOfClosestPointLeftOf(t); + + // Integration + double result = 0; + for (int i = 0, j = 0; i < closestLeftIndex; 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))))); + } + + double offset = t - this.points[closestLeftIndex]; + int k = closestLeftIndex << 2; + + return result + + (offset * (this.coefficients[k] + + ((offset * (this.coefficients[k + 1] * 0.5)) + + (offset * (this.coefficients[k + 2] / 3)) + + (offset * this.coefficients[k + 3] * 0.25)))); + } + + /// + /// Find the index of the greatest sample point smaller than t. + /// + /// The value to look for. + /// The sample point index. + private int IndexOfClosestPointLeftOf(double t) { // Binary search in the [ t[0], ..., t[n-2] ] (t[n-1] is not included) int low = 0; @@ -248,17 +252,7 @@ namespace MathNet.Numerics.Interpolation.Algorithms } } - // 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)))); + return low; } } }