From 8968416790a42eebd5583718800ccc4fe55c8cbc Mon Sep 17 00:00:00 2001 From: Gilles Sadowski Date: Mon, 19 Jul 2021 18:40:52 +0200 Subject: [PATCH] MATH-1379: Fix "LoessInterpolator" in case of unevenly spaced samples. Thanks to Richard Wilkinson. --- .../interpolation/LoessInterpolator.java | 10 ++- .../interpolation/LoessInterpolatorTest.java | 76 +++++++++++++++++++ src/changes/changes.xml | 3 + 3 files changed, 86 insertions(+), 3 deletions(-) diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/interpolation/LoessInterpolator.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/interpolation/LoessInterpolator.java index c9481d2c3..63a94f77b 100644 --- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/interpolation/LoessInterpolator.java +++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/interpolation/LoessInterpolator.java @@ -409,7 +409,8 @@ public class LoessInterpolator * {@code (right==xval.length-1 or xval[right+1] - xval[i] > xval[i] - xval[left])}. * The array will be updated. */ - private static void updateBandwidthInterval(final double[] xval, final double[] weights, + private static void updateBandwidthInterval(final double[] xval, + final double[] weights, final int i, final int[] bandwidthInterval) { final int left = bandwidthInterval[0]; @@ -418,10 +419,13 @@ public class LoessInterpolator // The right edge should be adjusted if the next point to the right // is closer to xval[i] than the leftmost point of the current interval int nextRight = nextNonzero(weights, right); - if (nextRight < xval.length && xval[nextRight] - xval[i] < xval[i] - xval[left]) { - int nextLeft = nextNonzero(weights, bandwidthInterval[0]); + int nextLeft = left; + while (nextRight < xval.length && + xval[nextRight] - xval[i] < xval[i] - xval[nextLeft]) { + nextLeft = nextNonzero(weights, bandwidthInterval[0]); bandwidthInterval[0] = nextLeft; bandwidthInterval[1] = nextRight; + nextRight = nextNonzero(weights, nextRight); } } diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/analysis/interpolation/LoessInterpolatorTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/analysis/interpolation/LoessInterpolatorTest.java index 6f47e58b2..5e2f9eb67 100644 --- a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/analysis/interpolation/LoessInterpolatorTest.java +++ b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/analysis/interpolation/LoessInterpolatorTest.java @@ -247,6 +247,82 @@ public class LoessInterpolatorTest { } } + // MATH-1379 + @Test + public void testFitWithUnevenXSpacing() { + final double[] xval = { + 0.1, 0.12, 0.23, 0.4, 0.57, + 0.7, 0.87, 1.3, 1.9, 2.2, + 2.3, 2.65, 3.0, 3.1, 3.5, + 4.6, 4.7, 5.8, 5.95, 6.1 + }; + final double[] yval = { + 0.47, 0.48, 0.55, 0.56, -0.08, + -0.04, -0.07, -0.07, -0.56, -0.46, + -0.56, -0.52, -3.03, -3.08, -3.09, + -3.04, 3.54, 3.46, 3.36, 3.35 + }; + + // Compare with output from R. + // predict(loess(y ~ x, data.frame(x=xval, y=yval), span=0.35, degree=1, family="symmetric", control=loess.control(iterations=1, surface="direct"))) + final double[] yref = { + 0.556184894, 0.541907126, 0.455059334, 0.303681477, 0.142126445, + 0.002615653, -0.031178445, -0.187124310, -0.405235207, -0.535023851, + -0.706801740, -1.466740294, -2.349248503, -2.596576469, -3.354222419, + 0.086206868, 0.320251370, 3.064778450, 3.426179479, 3.783500164 + }; + + final double delta = 1e-8; + // Note R counts all iterations whereas LoessInterpolator robustness iterations are + // in addition to the initial fit. + final LoessInterpolator li = new LoessInterpolator(0.35, 0, 1e-12); + final double[] res = li.smooth(xval, yval); + Assert.assertEquals(xval.length, res.length); + for (int i = 0; i < res.length; ++i) { + Assert.assertEquals(yref[i], res[i], delta); + } + } + + // MATH-1379 + @Test + public void testFitWithVaryingWeightsAndUnevenXSpacing() { + final double[] xval = { + 0.1, 0.12, 0.23, 0.4, 0.57, + 0.7, 0.87, 1.3, 1.9, 2.2, + 2.3, 2.65, 3.0, 3.1, 3.5, + 4.6, 4.7, 5.8, 5.95, 6.1 + }; + final double[] yval = { + 0.47, 0.48, 0.55, 0.56, -0.08, + -0.04, -0.07, -0.07, -0.56, -0.46, + -0.56, -0.52, -3.03, -3.08, -3.09, + -3.04, 3.54, 3.46, 3.36, 3.35}; + final double[] weights = { + 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, + 1, 0.8, 0.5, 0.5, 0.8, + 0.8, 0.5, 0.5, 0.9, 1 + }; + + // Compare with output from R. + // predict(loess(y ~ x, data.frame(x=xval, y=yval), weights, span=0.35, degree=1, family="symmetric", control=loess.control(iterations=1, surface="direct"))) + final double[] yref = { + 0.556184894, 0.541907126, 0.455059334, 0.303681477, 0.142126445, + 0.002615653, -0.031178445, -0.187124310, -0.406403569, -0.531957113, + -0.669978426, -1.411039850, -2.225022609, -2.463874517, -3.298767758, + -0.506116921, -0.231242991, 2.873755984, 3.295897106, 3.715495321}; + + final double delta = 1e-8; + // Note R counts all iterations whereas LoessInterpolator robustness iterations are + // in addition to the initial fit. + final LoessInterpolator li = new LoessInterpolator(0.35, 0, 1e-12); + final double[] res = li.smooth(xval, yval, weights); + Assert.assertEquals(xval.length, res.length); + for (int i = 0; i < res.length; ++i) { + Assert.assertEquals(yref[i], res[i], delta); + } + } + private void generateSineData(double[] xval, double[] yval, double xnoise, double ynoise) { double dx = 2 * AccurateMath.PI / xval.length; double x = 0; diff --git a/src/changes/changes.xml b/src/changes/changes.xml index 55b15bd13..2db01a576 100644 --- a/src/changes/changes.xml +++ b/src/changes/changes.xml @@ -95,6 +95,9 @@ Caveat: nightmare was one of the main reasons for creating more focused components.] "> + + Fix "LoessInterpolator" (in package "o.a.c.m.legacy.analysis.interpolation"). + Class "BigReal": Fix equality check.