diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/ml/distance/CanberraDistance.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/ml/distance/CanberraDistance.java index 8f9295bef0..d1f82a97c0 100644 --- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/ml/distance/CanberraDistance.java +++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/ml/distance/CanberraDistance.java @@ -31,9 +31,16 @@ public double compute(double[] a, double[] b) { MathArrays.checkEqualLength(a, b); double sum = 0; for (int i = 0; i < a.length; i++) { + final double absA = JdkMath.abs(a[i]); + final double absB = JdkMath.abs(b[i]); final double num = JdkMath.abs(a[i] - b[i]); - final double denom = JdkMath.abs(a[i]) + JdkMath.abs(b[i]); - sum += num == 0.0 && denom == 0.0 ? 0.0 : num / denom; + final double denom = absA + absB; + if (Double.isInfinite(denom) && Double.isFinite(a[i]) && Double.isFinite(b[i])) { + // Rescale large finite coordinates to avoid intermediate overflow. + sum += JdkMath.abs(a[i] * 0.5 - b[i] * 0.5) / (absA * 0.5 + absB * 0.5); + } else { + sum += denom == 0.0 ? 0.0 : num / denom; + } } return sum; } diff --git a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/ml/distance/CanberraDistanceTest.java b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/ml/distance/CanberraDistanceTest.java index 85384ba116..8b38df937b 100644 --- a/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/ml/distance/CanberraDistanceTest.java +++ b/commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/ml/distance/CanberraDistanceTest.java @@ -16,6 +16,9 @@ */ package org.apache.commons.math4.legacy.ml.distance; +import java.math.BigDecimal; +import java.math.MathContext; + import org.junit.Assert; import org.junit.Test; @@ -45,4 +48,81 @@ public void test() { Assert.assertEquals(expected, distance.compute(a, b), 0d); Assert.assertEquals(expected, distance.compute(b, a), 0d); } + + @Test + public void testExtremeOppositeSigns() { + // Both |MAX - (-MAX)| and |MAX| + |-MAX| overflow to infinity: the true ratio is 1. + final double max = Double.MAX_VALUE; + Assert.assertEquals(1, distance.compute(new double[] {max}, new double[] {-max}), 0d); + Assert.assertEquals(1, distance.compute(new double[] {-max}, new double[] {max}), 0d); + Assert.assertEquals(1, distance.compute(new double[] {max}, new double[] {-max / 2}), 0d); + Assert.assertEquals(1, distance.compute(new double[] {-max / 2}, new double[] {max}), 0d); + } + + @Test + public void testExtremeSameSign() { + // |MAX - MAX/2| / (|MAX| + |MAX/2|) == 1/3; only the denominator overflows. + final double max = Double.MAX_VALUE; + // 1/3 is not representable, so allow one ulp instead of an exact match. + final double expected = 1.0 / 3.0; + final double tol = Math.ulp(expected); + Assert.assertEquals(expected, distance.compute(new double[] {max}, new double[] {max / 2}), tol); + Assert.assertEquals(expected, distance.compute(new double[] {max / 2}, new double[] {max}), tol); + Assert.assertEquals(expected, distance.compute(new double[] {-max}, new double[] {-max / 2}), tol); + Assert.assertEquals(expected, distance.compute(new double[] {-max / 2}, new double[] {-max}), tol); + // Identical extreme values have distance 0. + Assert.assertEquals(0, distance.compute(new double[] {max}, new double[] {max}), 0d); + Assert.assertEquals(0, distance.compute(new double[] {-max}, new double[] {-max}), 0d); + } + + @Test + public void testExtremeVector() { + final double max = Double.MAX_VALUE; + final double[] a = { max, max / 2, max }; + final double[] b = { -max, max, max }; + // 1 + 1/3 + 0 + Assert.assertEquals(4.0 / 3.0, distance.compute(a, b), 1e-15); + Assert.assertEquals(4.0 / 3.0, distance.compute(b, a), 1e-15); + } + + @Test + public void testExtremeFiniteValuesAgainstBigDecimalReference() { + final double max = Double.MAX_VALUE; + final double[] a = {max, max, 0.75 * max, max, -max, max}; + final double[] b = {max / 2, Math.nextDown(max), Math.nextDown(0.75 * max), -max, -max / 2, max / 16}; + for (int i = 0; i < a.length; i++) { + final double expected = referenceCanberra(a[i], b[i]); + assertWithinOneUlp(expected, distance.compute(new double[] {a[i]}, new double[] {b[i]})); + assertWithinOneUlp(expected, distance.compute(new double[] {b[i]}, new double[] {a[i]})); + } + } + + /** + * Computes a high-precision reference from the exact decimal representations of the input doubles. + */ + private static double referenceCanberra(double a, double b) { + final BigDecimal bdA = new BigDecimal(a); + final BigDecimal bdB = new BigDecimal(b); + final BigDecimal num = bdA.subtract(bdB).abs(); + final BigDecimal den = bdA.abs().add(bdB.abs()); + return num.divide(den, MathContext.DECIMAL128).doubleValue(); + } + + private static void assertWithinOneUlp(double expected, double actual) { + Assert.assertEquals("expected " + expected + " but was " + actual, + expected, actual, Math.ulp(expected)); + } + + @Test + public void testSpecialValuesPreserved() { + final double min = Double.MIN_VALUE; + // Subnormal ratios are unaffected by the extreme-value handling. + Assert.assertEquals(1.0 / 3.0, distance.compute(new double[] {min}, new double[] {2 * min}), 0d); + Assert.assertEquals(1, distance.compute(new double[] {min}, new double[] {-min}), 0d); + // Signed zero: the zero/zero convention contributes 0. + Assert.assertEquals(0, distance.compute(new double[] {0.0}, new double[] {-0.0}), 0d); + // Infinite and NaN inputs keep their historical NaN result. + Assert.assertTrue(Double.isNaN(distance.compute(new double[] {Double.POSITIVE_INFINITY}, new double[] {1}))); + Assert.assertTrue(Double.isNaN(distance.compute(new double[] {Double.NaN}, new double[] {0}))); + } } diff --git a/src/changes/changes.xml b/src/changes/changes.xml index 077104d34d..b111a20cc7 100644 --- a/src/changes/changes.xml +++ b/src/changes/changes.xml @@ -101,6 +101,9 @@ Caveat: Instances not cannot be constructed that would silently overwrite entries due to key duplication from interger overflow. + + "CanberraDistance": Avoid overflow for large finite coordinates. + "SynchronizedDescriptiveStatistics": Synchronize all mutators of the internal state. Updates previously unsynchronized methods