Index: src/test/org/apache/commons/math/analysis/interpolation/LoessInterpolatorTest.java
===================================================================
--- src/test/org/apache/commons/math/analysis/interpolation/LoessInterpolatorTest.java (revision 0)
+++ src/test/org/apache/commons/math/analysis/interpolation/LoessInterpolatorTest.java (revision 0)
@@ -0,0 +1,268 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.analysis.interpolation;
+
+import junit.framework.TestCase;
+import org.apache.commons.math.MathException;
+
+/**
+ * Test of the LoessInterpolator class.
+ */
+public class LoessInterpolatorTest extends TestCase {
+ public void testOnOnePoint() throws MathException {
+ double[] xval = {0.5};
+ double[] yval = {0.7};
+ double[] res = new LoessInterpolator().smooth(xval, yval);
+ assertEquals(1, res.length);
+ assertEquals(0.7, res[0]);
+ }
+
+ public void testOnTwoPoints() throws MathException {
+ double[] xval = {0.5, 0.6};
+ double[] yval = {0.7, 0.8};
+ double[] res = new LoessInterpolator().smooth(xval, yval);
+ assertEquals(2, res.length);
+ assertEquals(0.7, res[0]);
+ assertEquals(0.8, res[1]);
+ }
+
+ public void testOnStraightLine() throws MathException {
+ double[] xval = {1,2,3,4,5};
+ double[] yval = {2,4,6,8,10};
+ LoessInterpolator li = new LoessInterpolator();
+ li.setBandwidth(0.6);
+ double[] res = li.smooth(xval, yval);
+ assertEquals(5, res.length);
+ for(int i = 0; i < 5; ++i) {
+ assertEquals(yval[i], res[i], 1e-8);
+ }
+ }
+
+ public void testOnDistortedSine() throws MathException {
+ int numPoints = 100;
+ double[] xval = new double[numPoints];
+ double[] yval = new double[numPoints];
+ double xnoise = 0.1;
+ double ynoise = 0.2;
+
+ generateSineData(xval, yval, xnoise, ynoise);
+
+ LoessInterpolator li = new LoessInterpolator();
+ li.setBandwidth(0.3);
+ li.setRobustnessIters(4);
+
+ double[] res = li.smooth(xval, yval);
+
+ // Check that the resulting curve differs from
+ // the "real" sine less than the jittered one
+
+ double noisyResidualSum = 0;
+ double fitResidualSum = 0;
+
+ System.out.println();
+ for(int i = 0; i < numPoints; ++i) {
+ double expected = Math.sin(xval[i]);
+ double noisy = yval[i];
+ double fit = res[i];
+
+ noisyResidualSum += Math.pow(noisy - expected, 2);
+ fitResidualSum += Math.pow(fit - expected, 2);
+ }
+
+ assertTrue(fitResidualSum < noisyResidualSum);
+ }
+
+ public void testIncreasingBandwidthIncreasesSmoothness() throws MathException {
+ int numPoints = 100;
+ double[] xval = new double[numPoints];
+ double[] yval = new double[numPoints];
+ double xnoise = 0.1;
+ double ynoise = 0.1;
+
+ generateSineData(xval, yval, xnoise, ynoise);
+
+ LoessInterpolator li = new LoessInterpolator();
+ li.setRobustnessIters(4);
+
+ // Check that variance decreases as bandwidth increases
+
+ double[] bandwidths = {0.1, 0.5, 1.0};
+ double[] variances = new double[bandwidths.length];
+ for (int i = 0; i < bandwidths.length; i++) {
+ double bw = bandwidths[i];
+
+ li.setBandwidth(bw);
+
+ double[] res = li.smooth(xval, yval);
+
+ for (int j = 1; j < res.length; ++j) {
+ variances[i] += Math.pow(res[j] - res[j-1], 2);
+ }
+ }
+
+ for(int i = 1; i < variances.length; ++i) {
+ assertTrue(variances[i] < variances[i-1]);
+ }
+ }
+
+ public void testIncreasingRobustnessItersIncreasesSmoothnessWithOutliers() throws MathException {
+ int numPoints = 100;
+ double[] xval = new double[numPoints];
+ double[] yval = new double[numPoints];
+ double xnoise = 0.1;
+ double ynoise = 0.1;
+
+ generateSineData(xval, yval, xnoise, ynoise);
+
+ // Introduce a couple of outliers
+ yval[numPoints/3] *= 100;
+ yval[2 * numPoints/3] *= -100;
+
+ LoessInterpolator li = new LoessInterpolator();
+ li.setBandwidth(0.3);
+
+ // Check that variance decreases as the number of robustness
+ // iterations increases
+
+ double[] variances = new double[4];
+ for (int i = 0; i < 4; i++) {
+ li.setRobustnessIters(i);
+
+ double[] res = li.smooth(xval, yval);
+
+ for (int j = 1; j < res.length; ++j) {
+ variances[i] += Math.abs(res[j] - res[j-1]);
+ }
+ }
+
+ for(int i = 1; i < variances.length; ++i) {
+ assertTrue(variances[i] < variances[i-1]);
+ }
+ }
+
+ public void testUnequalSizeArguments() {
+ try {
+ new LoessInterpolator().smooth(new double[] {1,2,3}, new double[] {1,2,3,4});
+ fail();
+ } catch(MathException e) {
+ // Expected
+ }
+ }
+
+ public void testEmptyData() {
+ try {
+ new LoessInterpolator().smooth(new double[] {}, new double[] {});
+ fail();
+ } catch(MathException e) {
+ // Expected
+ }
+ }
+
+ public void testNonStrictlyIncreasing() {
+ try {
+ new LoessInterpolator().smooth(new double[] {4,3,1,2}, new double[] {3,4,5,6});
+ fail();
+ } catch(MathException e) {
+ // Expected
+ }
+ try {
+ new LoessInterpolator().smooth(new double[] {1,2,2,3}, new double[] {3,4,5,6});
+ fail();
+ } catch(MathException e) {
+ // Expected
+ }
+ }
+
+ public void testNotAllFiniteReal() {
+ try {
+ new LoessInterpolator().smooth(new double[] {1,2,Double.NaN}, new double[] {3,4,5});
+ fail();
+ } catch(MathException e) {
+ // Expected
+ }
+ try {
+ new LoessInterpolator().smooth(new double[] {1,2,Double.POSITIVE_INFINITY}, new double[] {3,4,5});
+ fail();
+ } catch(MathException e) {
+ // Expected
+ }
+ try {
+ new LoessInterpolator().smooth(new double[] {1,2,Double.NEGATIVE_INFINITY}, new double[] {3,4,5});
+ fail();
+ } catch(MathException e) {
+ // Expected
+ }
+ try {
+ new LoessInterpolator().smooth(new double[] {3,4,5}, new double[] {1,2,Double.NaN});
+ fail();
+ } catch(MathException e) {
+ // Expected
+ }
+ try {
+ new LoessInterpolator().smooth(new double[] {3,4,5}, new double[] {1,2,Double.POSITIVE_INFINITY});
+ fail();
+ } catch(MathException e) {
+ // Expected
+ }
+ try {
+ new LoessInterpolator().smooth(new double[] {3,4,5}, new double[] {1,2,Double.NEGATIVE_INFINITY});
+ fail();
+ } catch(MathException e) {
+ // Expected
+ }
+ }
+
+ public void testInsufficientBandwidth() {
+ try {
+ LoessInterpolator li = new LoessInterpolator();
+ li.setBandwidth(0.1);
+ li.smooth(new double[] {1,2,3,4,5,6,7,8,9,10,11,12}, new double[] {1,2,3,4,5,6,7,8,9,10,11,12});
+ fail();
+ } catch(MathException e) {
+ // Expected
+ }
+ }
+
+ public void testCompletelyIncorrectBandwidth() {
+ try {
+ LoessInterpolator li = new LoessInterpolator();
+ li.setBandwidth(-0.2);
+ li.smooth(new double[] {1,2,3}, new double[] {1,2,3});
+ fail();
+ } catch(MathException e) {
+ // Expected
+ }
+ try {
+ LoessInterpolator li = new LoessInterpolator();
+ li.setBandwidth(1.1);
+ li.smooth(new double[] {1,2,3}, new double[] {1,2,3});
+ fail();
+ } catch(MathException e) {
+ // Expected
+ }
+ }
+
+ private void generateSineData(double[] xval, double[] yval, double xnoise, double ynoise) {
+ double dx = 2 * Math.PI / xval.length;
+ double x = 0;
+ for(int i = 0; i < xval.length; ++i) {
+ xval[i] = x;
+ yval[i] = Math.sin(x) + (2 * Math.random() - 1) * ynoise;
+ x += dx * (1 + (2 * Math.random() - 1) * xnoise);
+ }
+ }
+}
Index: src/java/org/apache/commons/math/analysis/interpolation/LoessInterpolator.java
===================================================================
--- src/java/org/apache/commons/math/analysis/interpolation/LoessInterpolator.java (revision 0)
+++ src/java/org/apache/commons/math/analysis/interpolation/LoessInterpolator.java (revision 0)
@@ -0,0 +1,391 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.analysis.interpolation;
+
+import org.apache.commons.math.analysis.UnivariateRealFunction;
+import org.apache.commons.math.MathException;
+
+import java.io.Serializable;
+import java.util.Arrays;
+
+/**
+ * Implements the
+ * Local Regression Algorithm (also Loess, Lowess) for interpolation of
+ * real univariate functions.
+ *
+ * For reference, see
+ *
+ * William S. Cleveland - Robust Locally Weighted Regression and Smoothing
+ * Scatterplots
+ *
+ * This class implements both the loess method and serves as an interpolation
+ * adapter to it, allowing to build a spline on the obtained loess fit.
+ *
+ * @version $Revision: ??? $ $Date: 2009-06-18 20:05:00 +0300 (Thu, 18 Jun 2009) $
+ */
+public class LoessInterpolator
+ implements UnivariateRealInterpolator, Serializable {
+
+ /** serializable version identifier. */
+ static final long serialVersionUID = 5204927143605193821L;
+
+ /**
+ * Default value of the bandwidth parameter.
+ */
+ public static final double DEFAULT_BANDWIDTH = 0.3;
+ /**
+ * Default value of the number of robustness iterations.
+ */
+ public static final int DEFAULT_ROBUSTNESS_ITERS = 2;
+
+ /**
+ * The bandwidth parameter: when computing the loess fit at
+ * a particular point, this fraction of source points closest
+ * to the current point is taken into account for computing
+ * a least-squares regression.
+ *
+ * A sensible value is usually 0.25 to 0.5.
+ */
+ private double bandwidth = DEFAULT_BANDWIDTH;
+
+ /**
+ * The number of robustness iterations parameter: this many
+ * robustness iterations are done.
+ *
+ * A sensible value is usually 0 (just the initial fit without any
+ * robustness iterations) to 4.
+ */
+ private int robustnessIters = DEFAULT_ROBUSTNESS_ITERS;
+
+ /**
+ * Constructs a new
+ * {@link org.apache.commons.math.analysis.interpolation.LoessInterpolator}
+ * with a bandwidth of {@link #DEFAULT_BANDWIDTH} and
+ * {@link #DEFAULT_ROBUSTNESS_ITERS} robustness iterations.
+ */
+ public LoessInterpolator() {
+ }
+
+ /**
+ * Compute an interpolating function by performing a loess fit
+ * on the data at the original abscissae and then building a cubic spline
+ * with a
+ * {@link org.apache.commons.math.analysis.interpolation.SplineInterpolator}
+ * on the resulting fit.
+ *
+ * @param xval the arguments for the interpolation points
+ * @param yval the values for the interpolation points
+ * @return A cubic spline built upon a loess fit to the data
+ * at the original abscissae
+ * @throws MathException if some of the following conditions are false:
+ * Arguments and values are of the same size that is greater than zero
+ * The arguments are in a strictly increasing order
+ * All arguments and values are finite real numbers
+ */
+ public final UnivariateRealFunction interpolate(
+ final double[] xval, final double[] yval) throws MathException {
+ return new SplineInterpolator().interpolate(xval, smooth(xval, yval));
+ }
+
+ /**
+ * Compute a loess fit on the data at the original abscissae.
+ *
+ * @param xval the arguments for the interpolation points
+ * @param yval the values for the interpolation points
+ * @return values of the loess fit at corresponding original abscissae
+ * @throws MathException if some of the following conditions are false:
+ * Arguments and values are of the same size that is greater than zero
+ * The arguments are in a strictly increasing order
+ * All arguments and values are finite real numbers
+ */
+ public final double[] smooth(final double[] xval, final double[] yval)
+ throws MathException {
+ if (xval.length != yval.length) {
+ throw new MathException(
+ "Loess expects the abscissa and ordinate arrays " +
+ "to be of the same size, " +
+ "but got {0} abscisssae and {1} ordinatae",
+ xval.length, yval.length);
+ }
+
+ int n = xval.length;
+
+ if (n == 0) {
+ throw new MathException(
+ "Loess expects at least 1 point");
+ }
+
+ checkAllFiniteReal(xval, "abscissae");
+ checkAllFiniteReal(yval, "ordinatae");
+
+ checkStrictlyIncreasing(xval);
+
+ if (n == 1) {
+ return new double[]{yval[0]};
+ }
+
+ if (n == 2) {
+ return new double[]{yval[0], yval[1]};
+ }
+
+ int bandwidthInPoints = (int) (bandwidth * n);
+
+ if (bandwidthInPoints < 2) {
+ throw new MathException(
+ "The bandwidth must be large enough to " +
+ "accomodate at least 2 points. There are {0} " +
+ " data points, and bandwidth must be at least {1} " +
+ " but it is only {2}",
+ n, 2.0 / n, bandwidth);
+ }
+
+ double[] res = new double[n];
+
+ double[] residuals = new double[n];
+ double[] sortedResiduals = new double[n];
+
+ double[] robustnessWeights = new double[n];
+
+ // Do an initial fit and 'robustnessIters' robustness iterations.
+ // This is equivalent to doing 'robustnessIters+1' robustness iterations
+ // starting with all robustness weights set to 1.
+ Arrays.fill(robustnessWeights, 1);
+
+ for (int iter = 0; iter <= robustnessIters; ++iter) {
+ int[] bandwidthInterval = {0, bandwidthInPoints-1};
+ // At each x, compute a local weighted linear regression
+ for (int i = 0; i < n; ++i) {
+ double x = xval[i];
+
+ // Find out the interval of source points on which
+ // a regression is to be made.
+ if(i > 0) {
+ updateBandwidthInterval(xval, i, bandwidthInterval);
+ }
+
+ int ileft = bandwidthInterval[0];
+ int iright = bandwidthInterval[1];
+
+ // Compute the point of the bandwidth interval that is
+ // farthest from x
+ int edge;
+ if (xval[i] - xval[ileft] > xval[iright] - xval[i]) {
+ edge = ileft;
+ } else {
+ edge = iright;
+ }
+
+ // Compute a least-squares linear fit weighted by
+ // the product of robustness weights and the tricube
+ // weight function.
+ // See http://en.wikipedia.org/wiki/Linear_regression
+ // (section "Univariate linear case")
+ // and http://en.wikipedia.org/wiki/Weighted_least_squares
+ // (section "Weighted least squares")
+ double sumWeights = 0;
+ double sumX = 0, sumXSquared = 0, sumY = 0, sumXY = 0;
+ double denom = Math.abs(1.0 / (xval[edge] - x));
+ for (int k = ileft; k <= iright; ++k) {
+ double xk = xval[k], yk = yval[k];
+ double dist;
+ if (k < i) {
+ dist = (x - xk);
+ } else {
+ dist = (xk - x);
+ }
+ double w = tricube(dist * denom) * robustnessWeights[k];
+ double xkw = xk * w;
+ sumWeights += w;
+ sumX += xkw;
+ sumXSquared += xk * xkw;
+ sumY += yk * w;
+ sumXY += yk * xkw;
+ }
+
+ double meanX = sumX / sumWeights;
+ double meanY = sumY / sumWeights;
+ double meanXY = sumXY / sumWeights;
+ double meanXSquared = sumXSquared / sumWeights;
+
+ double beta;
+ if(meanXSquared == meanX * meanX) {
+ beta = 0;
+ } else {
+ beta = (meanXY - meanX * meanY) / (meanXSquared - meanX * meanX);
+ }
+
+ double alpha = meanY - beta * meanX;
+
+ res[i] = beta * x + alpha;
+ residuals[i] = Math.abs(yval[i] - res[i]);
+ }
+
+ // No need to recompute the robustness weights at the last
+ // iteration, they won't be needed anymore
+ if(iter == robustnessIters) {
+ break;
+ }
+
+ // Recompute the robustness weights.
+
+ // Find the median residual.
+ // An arraycopy and a sort are completely tractable here,
+ // because the preceding loop is a lot more expensive
+ System.arraycopy(residuals, 0, sortedResiduals, 0, n);
+ Arrays.sort(sortedResiduals);
+ double medianResidual = sortedResiduals[n / 2];
+
+ if(medianResidual == 0) {
+ break;
+ }
+
+ for (int i = 0; i < n; ++i) {
+ double arg = residuals[i] / (6 * medianResidual);
+ robustnessWeights[i] = (arg >= 1) ?
+ 0 :
+ Math.pow(1 - arg * arg, 2);
+ }
+ }
+
+ return res;
+ }
+
+ /**
+ * Given an index interval into xval that embraces a certain number of
+ * points closest to xval[i-1], update the interval so that it embraces
+ * the same number of points closest to xval[i]
+ *
+ * @param xval arguments array
+ * @param i the index around which the new interval should be computed
+ * @param bandwidthInterval a two-element array {left, right} such that:
+ * (left==0 or xval[i] - xval[left-1] > xval[right] - xval[i])
+ * and also
+ * (right==xval.length-1 or xval[right+1] - xval[i] > xval[i] - xval[left]).
+ * The array will be updated.
+ */
+ private static void updateBandwidthInterval(double[] xval, int i, int[] bandwidthInterval) {
+ int left = bandwidthInterval[0];
+ int right = bandwidthInterval[1];
+
+ // 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
+ if(right < xval.length - 1 &&
+ xval[right+1] - xval[i] < xval[i] - xval[left])
+ {
+ bandwidthInterval[0]++;
+ bandwidthInterval[1]++;
+ }
+ }
+
+ /**
+ * Compute the tricube weight function
+ * http://en.wikipedia.org/wiki/Local_regression#Weight_function
+ *
+ * @param x the argument
+ * @return (1-|x|^3)^3
+ */
+ private static double tricube(double x) {
+ double tmp = 1 - x * x * x;
+ return tmp * tmp * tmp;
+ }
+
+ /**
+ * Check that all elements of an array are finite real numbers.
+ *
+ * @param values the values array
+ * @param pluralLabel plural form of a label for the values that will
+ * be used in text of the exception
+ * @throws MathException if one of the values is not
+ * a finite real number
+ */
+ private static void checkAllFiniteReal(double[] values, String pluralLabel) throws MathException {
+ for (int i = 0; i < values.length; i++) {
+ double x = values[i];
+ if (Double.isInfinite(x) || Double.isNaN(x)) {
+ throw new MathException(
+ "All {0} must be finite real numbers, " +
+ "but {1}-th is {2}",
+ pluralLabel, i, x);
+ }
+ }
+ }
+
+ /**
+ * Check that elements of the abscissae array are in a strictly
+ * increasing order.
+ *
+ * @param xval the abscissae array
+ * @throws MathException if the abscissae array
+ * is not in a strictly increasing order
+ */
+ private static void checkStrictlyIncreasing(double[] xval) throws MathException {
+ for (int i = 0; i < xval.length; ++i) {
+ if (i >= 1 && xval[i - 1] >= xval[i]) {
+ throw new MathException(
+ "The abscissae array must be sorted in a strictly " +
+ "increasing order, but the {0}-th element is {1} " +
+ "whereas {2}-th is {3}",
+ i - 1, xval[i - 1], i, xval[i]);
+ }
+ }
+ }
+
+ /**
+ * Set the bandwidth parameter.
+ *
+ * When computing the loess fit at
+ * a particular point, this fraction of source points closest
+ * to the current point is taken into account for computing
+ * a least-squares regression.
+ *
+ * A sensible value is usually 0.25 to 0.5, the default value is
+ * {@link #DEFAULT_BANDWIDTH}.
+ *
+ * @param bandwidth the new value of the bandwidth parameter
+ * @throws MathException if bandwidth does not lie in the interval [0,1]
+ */
+ public void setBandwidth(double bandwidth) throws MathException {
+ if (bandwidth < 0 || bandwidth > 1) {
+ throw new MathException(
+ "Bandwidth must be in the interval [0,1], but got {0}",
+ bandwidth);
+ }
+ this.bandwidth = bandwidth;
+ }
+
+ /**
+ * Set the number of robustness iterations.
+ *
+ * This many robustness iterations are done.
+ *
+ * A sensible value is usually 0 (just the initial fit without any
+ * robustness iterations) to 4.
+ *
+ * @param robustnessIters the new value of the number of robustness
+ * iterations
+ * @throws MathException if robustnessIters is negative
+ */
+ public void setRobustnessIters(int robustnessIters) throws MathException {
+ if (robustnessIters < 0) {
+ throw new MathException(
+ "The number of robustness iterations must " +
+ "be non-negative, but got {0}",
+ robustnessIters);
+ }
+ this.robustnessIters = robustnessIters;
+ }
+}
Index: src/java/org/apache/commons/math/ode/AbstractIntegrator.java
===================================================================
--- src/java/org/apache/commons/math/ode/AbstractIntegrator.java (revision 786076)
+++ src/java/org/apache/commons/math/ode/AbstractIntegrator.java (working copy)
@@ -78,6 +78,8 @@
resetEvaluations();
}
+ /** Build an instance with a name of null.
+ */
protected AbstractIntegrator() {
this(null);
}