MATH-1529: Modified Akima interpolator.

src/main/java/org/apache/commons/math4/analysis/interpolation/AkimaSplineInterpolator.java

@@ -36,7 +36,13 @@ import org.apache.commons.numbers.core.Precision;
  * <p>
  * This implementation is based on the Akima implementation in the CubicSpline
  * class in the Math.NET Numerics library. The method referenced is
- * CubicSpline.InterpolateAkimaSorted
+ * "CubicSpline.InterpolateAkimaSorted".
+ * </p>
+ * <p>
+ * When the {@link #AkimaSplineInterpolator(boolean) argument to the constructor}
+ * is {@code true}, the weights are modified in order to
+ * <a href="https://nl.mathworks.com/help/matlab/ref/makima.html">smooth out
+ * spurious oscillations</a>.
  * </p>
  * <p>
  * The {@link #interpolate(double[], double[]) interpolate} method returns a
@@ -49,6 +55,24 @@ public class AkimaSplineInterpolator
     implements UnivariateInterpolator {
     /** The minimum number of points that are needed to compute the function. */
     private static final int MINIMUM_NUMBER_POINTS = 5;
+    /** Whether to use the "modified Akima" interpolation. */
+    private final boolean useModified;
+
+    /**
+     * Uses the original Akima algorithm.
+     */
+    public AkimaSplineInterpolator() {
+        this(false);
+    }
+
+    /**
+     * @param useModified Whether to use the
+     * <a href="https://nl.mathworks.com/help/matlab/ref/makima.html">modified Akima</a>
+     * interpolation.
+     */
+    public AkimaSplineInterpolator(boolean useModified) {
+        this.useModified = useModified;
+    }
 
     /**
      * Computes an interpolating function for the data set.
@@ -95,8 +119,16 @@ public class AkimaSplineInterpolator
             differences[i] = (yvals[i + 1] - yvals[i]) / (xvals[i + 1] - xvals[i]);
         }
 
-        for (int i = 1; i < weights.length; i++) {
-            weights[i] = FastMath.abs(differences[i] - differences[i - 1]);
+        if (useModified) {
+            for (int i = 1; i < weights.length; i++) {
+                final double a = differences[i];
+                final double b = differences[i - 1];
+                weights[i] = FastMath.abs(a - b) + 0.5 * FastMath.abs(a + b);
+            }
+        } else {
+            for (int i = 1; i < weights.length; i++) {
+                weights[i] = FastMath.abs(differences[i] - differences[i - 1]);
+            }
         }
 
         // Prepare Hermite interpolation scheme.

src/test/java/org/apache/commons/math4/analysis/interpolation/AkimaSplineInterpolatorTest.java

@@ -187,6 +187,41 @@ public class AkimaSplineInterpolatorTest
                            maxTolerance );
     }
 
+    @Test
+    public void testOriginalVsModified() {
+        final UnivariateFunction f = new UnivariateFunction() {
+            @Override
+            public double value(double x) {
+                return x < -1 ? -1 :
+                    x < 1 ? x : 1;
+            }
+        };
+
+        final double[] xS = new double[] {-1, 0, 1, 2, 3 };
+        final double[] yS = new double[xS.length];
+
+        for (int i = 0; i < xS.length; i++) {
+            yS[i] = f.value(xS[i]);
+        }
+
+        final UnivariateFunction iOriginal = new AkimaSplineInterpolator(false).interpolate(xS, yS);
+        final UnivariateFunction iModified = new AkimaSplineInterpolator(true).interpolate(xS, yS);
+
+        final int n = 100;
+        final double delta = 1d / n;
+        for (int i = 1; i < n - 1; i++) {
+            final double x = 2 - i * delta;
+
+            final double value = f.value(x);
+            final double diffOriginal = Math.abs(iOriginal.value(x) - value);
+            final double diffModified = Math.abs(iModified.value(x) - value);
+
+            // In interval (1, 2), the modified algorithm eliminates interpolation artefacts.
+            Assert.assertTrue(diffOriginal > 0);
+            Assert.assertEquals(0d, diffModified, 0d);
+        }
+    }
+
     private void testInterpolation( double minimumX, double maximumX, int numberOfElements, int numberOfSamples,
                                     UnivariateFunction f, double tolerance, double maxTolerance )
     {