You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by tn...@apache.org on 2015/02/25 23:02:45 UTC
[3/4] [math] Remove deprecated interpolation and fitter classes.
http://git-wip-us.apache.org/repos/asf/commons-math/blob/0a5cd113/src/main/java/org/apache/commons/math4/fitting/CurveFitter.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math4/fitting/CurveFitter.java b/src/main/java/org/apache/commons/math4/fitting/CurveFitter.java
deleted file mode 100644
index 8cce426..0000000
--- a/src/main/java/org/apache/commons/math4/fitting/CurveFitter.java
+++ /dev/null
@@ -1,233 +0,0 @@
-/*
- * 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.math4.fitting;
-
-import java.util.ArrayList;
-import java.util.List;
-
-import org.apache.commons.math4.analysis.MultivariateMatrixFunction;
-import org.apache.commons.math4.analysis.MultivariateVectorFunction;
-import org.apache.commons.math4.analysis.ParametricUnivariateFunction;
-import org.apache.commons.math4.optim.InitialGuess;
-import org.apache.commons.math4.optim.MaxEval;
-import org.apache.commons.math4.optim.PointVectorValuePair;
-import org.apache.commons.math4.optim.nonlinear.vector.ModelFunction;
-import org.apache.commons.math4.optim.nonlinear.vector.ModelFunctionJacobian;
-import org.apache.commons.math4.optim.nonlinear.vector.MultivariateVectorOptimizer;
-import org.apache.commons.math4.optim.nonlinear.vector.Target;
-import org.apache.commons.math4.optim.nonlinear.vector.Weight;
-
-/**
- * Fitter for parametric univariate real functions y = f(x).
- * <br/>
- * When a univariate real function y = f(x) does depend on some
- * unknown parameters p<sub>0</sub>, p<sub>1</sub> ... p<sub>n-1</sub>,
- * this class can be used to find these parameters. It does this
- * by <em>fitting</em> the curve so it remains very close to a set of
- * observed points (x<sub>0</sub>, y<sub>0</sub>), (x<sub>1</sub>,
- * y<sub>1</sub>) ... (x<sub>k-1</sub>, y<sub>k-1</sub>). This fitting
- * is done by finding the parameters values that minimizes the objective
- * function ∑(y<sub>i</sub>-f(x<sub>i</sub>))<sup>2</sup>. This is
- * really a least squares problem.
- *
- * @param <T> Function to use for the fit.
- *
- * @since 2.0
- * @deprecated As of 3.3. Please use {@link AbstractCurveFitter} and
- * {@link WeightedObservedPoints} instead.
- */
-@Deprecated
-public class CurveFitter<T extends ParametricUnivariateFunction> {
- /** Optimizer to use for the fitting. */
- private final MultivariateVectorOptimizer optimizer;
- /** Observed points. */
- private final List<WeightedObservedPoint> observations;
-
- /**
- * Simple constructor.
- *
- * @param optimizer Optimizer to use for the fitting.
- * @since 3.1
- */
- public CurveFitter(final MultivariateVectorOptimizer optimizer) {
- this.optimizer = optimizer;
- observations = new ArrayList<WeightedObservedPoint>();
- }
-
- /** Add an observed (x,y) point to the sample with unit weight.
- * <p>Calling this method is equivalent to call
- * {@code addObservedPoint(1.0, x, y)}.</p>
- * @param x abscissa of the point
- * @param y observed value of the point at x, after fitting we should
- * have f(x) as close as possible to this value
- * @see #addObservedPoint(double, double, double)
- * @see #addObservedPoint(WeightedObservedPoint)
- * @see #getObservations()
- */
- public void addObservedPoint(double x, double y) {
- addObservedPoint(1.0, x, y);
- }
-
- /** Add an observed weighted (x,y) point to the sample.
- * @param weight weight of the observed point in the fit
- * @param x abscissa of the point
- * @param y observed value of the point at x, after fitting we should
- * have f(x) as close as possible to this value
- * @see #addObservedPoint(double, double)
- * @see #addObservedPoint(WeightedObservedPoint)
- * @see #getObservations()
- */
- public void addObservedPoint(double weight, double x, double y) {
- observations.add(new WeightedObservedPoint(weight, x, y));
- }
-
- /** Add an observed weighted (x,y) point to the sample.
- * @param observed observed point to add
- * @see #addObservedPoint(double, double)
- * @see #addObservedPoint(double, double, double)
- * @see #getObservations()
- */
- public void addObservedPoint(WeightedObservedPoint observed) {
- observations.add(observed);
- }
-
- /** Get the observed points.
- * @return observed points
- * @see #addObservedPoint(double, double)
- * @see #addObservedPoint(double, double, double)
- * @see #addObservedPoint(WeightedObservedPoint)
- */
- public WeightedObservedPoint[] getObservations() {
- return observations.toArray(new WeightedObservedPoint[observations.size()]);
- }
-
- /**
- * Remove all observations.
- */
- public void clearObservations() {
- observations.clear();
- }
-
- /**
- * Fit a curve.
- * This method compute the coefficients of the curve that best
- * fit the sample of observed points previously given through calls
- * to the {@link #addObservedPoint(WeightedObservedPoint)
- * addObservedPoint} method.
- *
- * @param f parametric function to fit.
- * @param initialGuess first guess of the function parameters.
- * @return the fitted parameters.
- * @throws org.apache.commons.math4.exception.DimensionMismatchException
- * if the start point dimension is wrong.
- */
- public double[] fit(T f, final double[] initialGuess) {
- return fit(Integer.MAX_VALUE, f, initialGuess);
- }
-
- /**
- * Fit a curve.
- * This method compute the coefficients of the curve that best
- * fit the sample of observed points previously given through calls
- * to the {@link #addObservedPoint(WeightedObservedPoint)
- * addObservedPoint} method.
- *
- * @param f parametric function to fit.
- * @param initialGuess first guess of the function parameters.
- * @param maxEval Maximum number of function evaluations.
- * @return the fitted parameters.
- * @throws org.apache.commons.math4.exception.TooManyEvaluationsException
- * if the number of allowed evaluations is exceeded.
- * @throws org.apache.commons.math4.exception.DimensionMismatchException
- * if the start point dimension is wrong.
- * @since 3.0
- */
- public double[] fit(int maxEval, T f,
- final double[] initialGuess) {
- // Prepare least squares problem.
- double[] target = new double[observations.size()];
- double[] weights = new double[observations.size()];
- int i = 0;
- for (WeightedObservedPoint point : observations) {
- target[i] = point.getY();
- weights[i] = point.getWeight();
- ++i;
- }
-
- // Input to the optimizer: the model and its Jacobian.
- final TheoreticalValuesFunction model = new TheoreticalValuesFunction(f);
-
- // Perform the fit.
- final PointVectorValuePair optimum
- = optimizer.optimize(new MaxEval(maxEval),
- model.getModelFunction(),
- model.getModelFunctionJacobian(),
- new Target(target),
- new Weight(weights),
- new InitialGuess(initialGuess));
- // Extract the coefficients.
- return optimum.getPointRef();
- }
-
- /** Vectorial function computing function theoretical values. */
- private class TheoreticalValuesFunction {
- /** Function to fit. */
- private final ParametricUnivariateFunction f;
-
- /**
- * @param f function to fit.
- */
- public TheoreticalValuesFunction(final ParametricUnivariateFunction f) {
- this.f = f;
- }
-
- /**
- * @return the model function values.
- */
- public ModelFunction getModelFunction() {
- return new ModelFunction(new MultivariateVectorFunction() {
- /** {@inheritDoc} */
- public double[] value(double[] point) {
- // compute the residuals
- final double[] values = new double[observations.size()];
- int i = 0;
- for (WeightedObservedPoint observed : observations) {
- values[i++] = f.value(observed.getX(), point);
- }
-
- return values;
- }
- });
- }
-
- /**
- * @return the model function Jacobian.
- */
- public ModelFunctionJacobian getModelFunctionJacobian() {
- return new ModelFunctionJacobian(new MultivariateMatrixFunction() {
- public double[][] value(double[] point) {
- final double[][] jacobian = new double[observations.size()][];
- int i = 0;
- for (WeightedObservedPoint observed : observations) {
- jacobian[i++] = f.gradient(observed.getX(), point);
- }
- return jacobian;
- }
- });
- }
- }
-}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/0a5cd113/src/main/java/org/apache/commons/math4/fitting/GaussianFitter.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math4/fitting/GaussianFitter.java b/src/main/java/org/apache/commons/math4/fitting/GaussianFitter.java
deleted file mode 100644
index 285a467..0000000
--- a/src/main/java/org/apache/commons/math4/fitting/GaussianFitter.java
+++ /dev/null
@@ -1,365 +0,0 @@
-/*
- * 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.math4.fitting;
-
-import java.util.Arrays;
-import java.util.Comparator;
-
-import org.apache.commons.math4.analysis.function.Gaussian;
-import org.apache.commons.math4.exception.NotStrictlyPositiveException;
-import org.apache.commons.math4.exception.NullArgumentException;
-import org.apache.commons.math4.exception.NumberIsTooSmallException;
-import org.apache.commons.math4.exception.OutOfRangeException;
-import org.apache.commons.math4.exception.ZeroException;
-import org.apache.commons.math4.exception.util.LocalizedFormats;
-import org.apache.commons.math4.optim.nonlinear.vector.MultivariateVectorOptimizer;
-import org.apache.commons.math4.util.FastMath;
-
-/**
- * Fits points to a {@link
- * org.apache.commons.math4.analysis.function.Gaussian.Parametric Gaussian} function.
- * <p>
- * Usage example:
- * <pre>
- * GaussianFitter fitter = new GaussianFitter(
- * new LevenbergMarquardtOptimizer());
- * fitter.addObservedPoint(4.0254623, 531026.0);
- * fitter.addObservedPoint(4.03128248, 984167.0);
- * fitter.addObservedPoint(4.03839603, 1887233.0);
- * fitter.addObservedPoint(4.04421621, 2687152.0);
- * fitter.addObservedPoint(4.05132976, 3461228.0);
- * fitter.addObservedPoint(4.05326982, 3580526.0);
- * fitter.addObservedPoint(4.05779662, 3439750.0);
- * fitter.addObservedPoint(4.0636168, 2877648.0);
- * fitter.addObservedPoint(4.06943698, 2175960.0);
- * fitter.addObservedPoint(4.07525716, 1447024.0);
- * fitter.addObservedPoint(4.08237071, 717104.0);
- * fitter.addObservedPoint(4.08366408, 620014.0);
- * double[] parameters = fitter.fit();
- * </pre>
- *
- * @since 2.2
- * @deprecated As of 3.3. Please use {@link GaussianCurveFitter} and
- * {@link WeightedObservedPoints} instead.
- */
-@Deprecated
-public class GaussianFitter extends CurveFitter<Gaussian.Parametric> {
- /**
- * Constructs an instance using the specified optimizer.
- *
- * @param optimizer Optimizer to use for the fitting.
- */
- public GaussianFitter(MultivariateVectorOptimizer optimizer) {
- super(optimizer);
- }
-
- /**
- * Fits a Gaussian function to the observed points.
- *
- * @param initialGuess First guess values in the following order:
- * <ul>
- * <li>Norm</li>
- * <li>Mean</li>
- * <li>Sigma</li>
- * </ul>
- * @return the parameters of the Gaussian function that best fits the
- * observed points (in the same order as above).
- * @since 3.0
- */
- public double[] fit(double[] initialGuess) {
- final Gaussian.Parametric f = new Gaussian.Parametric() {
- @Override
- public double value(double x, double ... p) {
- double v = Double.POSITIVE_INFINITY;
- try {
- v = super.value(x, p);
- } catch (NotStrictlyPositiveException e) { // NOPMD
- // Do nothing.
- }
- return v;
- }
-
- @Override
- public double[] gradient(double x, double ... p) {
- double[] v = { Double.POSITIVE_INFINITY,
- Double.POSITIVE_INFINITY,
- Double.POSITIVE_INFINITY };
- try {
- v = super.gradient(x, p);
- } catch (NotStrictlyPositiveException e) { // NOPMD
- // Do nothing.
- }
- return v;
- }
- };
-
- return fit(f, initialGuess);
- }
-
- /**
- * Fits a Gaussian function to the observed points.
- *
- * @return the parameters of the Gaussian function that best fits the
- * observed points (in the same order as above).
- */
- public double[] fit() {
- final double[] guess = (new ParameterGuesser(getObservations())).guess();
- return fit(guess);
- }
-
- /**
- * Guesses the parameters {@code norm}, {@code mean}, and {@code sigma}
- * of a {@link org.apache.commons.math4.analysis.function.Gaussian.Parametric}
- * based on the specified observed points.
- */
- public static class ParameterGuesser {
- /** Normalization factor. */
- private final double norm;
- /** Mean. */
- private final double mean;
- /** Standard deviation. */
- private final double sigma;
-
- /**
- * Constructs instance with the specified observed points.
- *
- * @param observations Observed points from which to guess the
- * parameters of the Gaussian.
- * @throws NullArgumentException if {@code observations} is
- * {@code null}.
- * @throws NumberIsTooSmallException if there are less than 3
- * observations.
- */
- public ParameterGuesser(WeightedObservedPoint[] observations) {
- if (observations == null) {
- throw new NullArgumentException(LocalizedFormats.INPUT_ARRAY);
- }
- if (observations.length < 3) {
- throw new NumberIsTooSmallException(observations.length, 3, true);
- }
-
- final WeightedObservedPoint[] sorted = sortObservations(observations);
- final double[] params = basicGuess(sorted);
-
- norm = params[0];
- mean = params[1];
- sigma = params[2];
- }
-
- /**
- * Gets an estimation of the parameters.
- *
- * @return the guessed parameters, in the following order:
- * <ul>
- * <li>Normalization factor</li>
- * <li>Mean</li>
- * <li>Standard deviation</li>
- * </ul>
- */
- public double[] guess() {
- return new double[] { norm, mean, sigma };
- }
-
- /**
- * Sort the observations.
- *
- * @param unsorted Input observations.
- * @return the input observations, sorted.
- */
- private WeightedObservedPoint[] sortObservations(WeightedObservedPoint[] unsorted) {
- final WeightedObservedPoint[] observations = unsorted.clone();
- final Comparator<WeightedObservedPoint> cmp
- = new Comparator<WeightedObservedPoint>() {
- public int compare(WeightedObservedPoint p1,
- WeightedObservedPoint p2) {
- if (p1 == null && p2 == null) {
- return 0;
- }
- if (p1 == null) {
- return -1;
- }
- if (p2 == null) {
- return 1;
- }
- if (p1.getX() < p2.getX()) {
- return -1;
- }
- if (p1.getX() > p2.getX()) {
- return 1;
- }
- if (p1.getY() < p2.getY()) {
- return -1;
- }
- if (p1.getY() > p2.getY()) {
- return 1;
- }
- if (p1.getWeight() < p2.getWeight()) {
- return -1;
- }
- if (p1.getWeight() > p2.getWeight()) {
- return 1;
- }
- return 0;
- }
- };
-
- Arrays.sort(observations, cmp);
- return observations;
- }
-
- /**
- * Guesses the parameters based on the specified observed points.
- *
- * @param points Observed points, sorted.
- * @return the guessed parameters (normalization factor, mean and
- * sigma).
- */
- private double[] basicGuess(WeightedObservedPoint[] points) {
- final int maxYIdx = findMaxY(points);
- final double n = points[maxYIdx].getY();
- final double m = points[maxYIdx].getX();
-
- double fwhmApprox;
- try {
- final double halfY = n + ((m - n) / 2);
- final double fwhmX1 = interpolateXAtY(points, maxYIdx, -1, halfY);
- final double fwhmX2 = interpolateXAtY(points, maxYIdx, 1, halfY);
- fwhmApprox = fwhmX2 - fwhmX1;
- } catch (OutOfRangeException e) {
- // TODO: Exceptions should not be used for flow control.
- fwhmApprox = points[points.length - 1].getX() - points[0].getX();
- }
- final double s = fwhmApprox / (2 * FastMath.sqrt(2 * FastMath.log(2)));
-
- return new double[] { n, m, s };
- }
-
- /**
- * Finds index of point in specified points with the largest Y.
- *
- * @param points Points to search.
- * @return the index in specified points array.
- */
- private int findMaxY(WeightedObservedPoint[] points) {
- int maxYIdx = 0;
- for (int i = 1; i < points.length; i++) {
- if (points[i].getY() > points[maxYIdx].getY()) {
- maxYIdx = i;
- }
- }
- return maxYIdx;
- }
-
- /**
- * Interpolates using the specified points to determine X at the
- * specified Y.
- *
- * @param points Points to use for interpolation.
- * @param startIdx Index within points from which to start the search for
- * interpolation bounds points.
- * @param idxStep Index step for searching interpolation bounds points.
- * @param y Y value for which X should be determined.
- * @return the value of X for the specified Y.
- * @throws ZeroException if {@code idxStep} is 0.
- * @throws OutOfRangeException if specified {@code y} is not within the
- * range of the specified {@code points}.
- */
- private double interpolateXAtY(WeightedObservedPoint[] points,
- int startIdx,
- int idxStep,
- double y)
- throws OutOfRangeException {
- if (idxStep == 0) {
- throw new ZeroException();
- }
- final WeightedObservedPoint[] twoPoints
- = getInterpolationPointsForY(points, startIdx, idxStep, y);
- final WeightedObservedPoint p1 = twoPoints[0];
- final WeightedObservedPoint p2 = twoPoints[1];
- if (p1.getY() == y) {
- return p1.getX();
- }
- if (p2.getY() == y) {
- return p2.getX();
- }
- return p1.getX() + (((y - p1.getY()) * (p2.getX() - p1.getX())) /
- (p2.getY() - p1.getY()));
- }
-
- /**
- * Gets the two bounding interpolation points from the specified points
- * suitable for determining X at the specified Y.
- *
- * @param points Points to use for interpolation.
- * @param startIdx Index within points from which to start search for
- * interpolation bounds points.
- * @param idxStep Index step for search for interpolation bounds points.
- * @param y Y value for which X should be determined.
- * @return the array containing two points suitable for determining X at
- * the specified Y.
- * @throws ZeroException if {@code idxStep} is 0.
- * @throws OutOfRangeException if specified {@code y} is not within the
- * range of the specified {@code points}.
- */
- private WeightedObservedPoint[] getInterpolationPointsForY(WeightedObservedPoint[] points,
- int startIdx,
- int idxStep,
- double y)
- throws OutOfRangeException {
- if (idxStep == 0) {
- throw new ZeroException();
- }
- for (int i = startIdx;
- idxStep < 0 ? i + idxStep >= 0 : i + idxStep < points.length;
- i += idxStep) {
- final WeightedObservedPoint p1 = points[i];
- final WeightedObservedPoint p2 = points[i + idxStep];
- if (isBetween(y, p1.getY(), p2.getY())) {
- if (idxStep < 0) {
- return new WeightedObservedPoint[] { p2, p1 };
- } else {
- return new WeightedObservedPoint[] { p1, p2 };
- }
- }
- }
-
- // Boundaries are replaced by dummy values because the raised
- // exception is caught and the message never displayed.
- // TODO: Exceptions should not be used for flow control.
- throw new OutOfRangeException(y,
- Double.NEGATIVE_INFINITY,
- Double.POSITIVE_INFINITY);
- }
-
- /**
- * Determines whether a value is between two other values.
- *
- * @param value Value to test whether it is between {@code boundary1}
- * and {@code boundary2}.
- * @param boundary1 One end of the range.
- * @param boundary2 Other end of the range.
- * @return {@code true} if {@code value} is between {@code boundary1} and
- * {@code boundary2} (inclusive), {@code false} otherwise.
- */
- private boolean isBetween(double value,
- double boundary1,
- double boundary2) {
- return (value >= boundary1 && value <= boundary2) ||
- (value >= boundary2 && value <= boundary1);
- }
- }
-}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/0a5cd113/src/main/java/org/apache/commons/math4/fitting/HarmonicFitter.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math4/fitting/HarmonicFitter.java b/src/main/java/org/apache/commons/math4/fitting/HarmonicFitter.java
deleted file mode 100644
index e74d0ac..0000000
--- a/src/main/java/org/apache/commons/math4/fitting/HarmonicFitter.java
+++ /dev/null
@@ -1,384 +0,0 @@
-/*
- * 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.math4.fitting;
-
-import org.apache.commons.math4.analysis.function.HarmonicOscillator;
-import org.apache.commons.math4.exception.MathIllegalStateException;
-import org.apache.commons.math4.exception.NumberIsTooSmallException;
-import org.apache.commons.math4.exception.ZeroException;
-import org.apache.commons.math4.exception.util.LocalizedFormats;
-import org.apache.commons.math4.optim.nonlinear.vector.MultivariateVectorOptimizer;
-import org.apache.commons.math4.util.FastMath;
-
-/**
- * Class that implements a curve fitting specialized for sinusoids.
- *
- * Harmonic fitting is a very simple case of curve fitting. The
- * estimated coefficients are the amplitude a, the pulsation ω and
- * the phase φ: <code>f (t) = a cos (ω t + φ)</code>. They are
- * searched by a least square estimator initialized with a rough guess
- * based on integrals.
- *
- * @since 2.0
- * @deprecated As of 3.3. Please use {@link HarmonicCurveFitter} and
- * {@link WeightedObservedPoints} instead.
- */
-@Deprecated
-public class HarmonicFitter extends CurveFitter<HarmonicOscillator.Parametric> {
- /**
- * Simple constructor.
- * @param optimizer Optimizer to use for the fitting.
- */
- public HarmonicFitter(final MultivariateVectorOptimizer optimizer) {
- super(optimizer);
- }
-
- /**
- * Fit an harmonic function to the observed points.
- *
- * @param initialGuess First guess values in the following order:
- * <ul>
- * <li>Amplitude</li>
- * <li>Angular frequency</li>
- * <li>Phase</li>
- * </ul>
- * @return the parameters of the harmonic function that best fits the
- * observed points (in the same order as above).
- */
- public double[] fit(double[] initialGuess) {
- return fit(new HarmonicOscillator.Parametric(), initialGuess);
- }
-
- /**
- * Fit an harmonic function to the observed points.
- * An initial guess will be automatically computed.
- *
- * @return the parameters of the harmonic function that best fits the
- * observed points (see the other {@link #fit(double[]) fit} method.
- * @throws NumberIsTooSmallException if the sample is too short for the
- * the first guess to be computed.
- * @throws ZeroException if the first guess cannot be computed because
- * the abscissa range is zero.
- */
- public double[] fit() {
- return fit((new ParameterGuesser(getObservations())).guess());
- }
-
- /**
- * This class guesses harmonic coefficients from a sample.
- * <p>The algorithm used to guess the coefficients is as follows:</p>
- *
- * <p>We know f (t) at some sampling points t<sub>i</sub> and want to find a,
- * ω and φ such that f (t) = a cos (ω t + φ).
- * </p>
- *
- * <p>From the analytical expression, we can compute two primitives :
- * <pre>
- * If2 (t) = ∫ f<sup>2</sup> = a<sup>2</sup> × [t + S (t)] / 2
- * If'2 (t) = ∫ f'<sup>2</sup> = a<sup>2</sup> ω<sup>2</sup> × [t - S (t)] / 2
- * where S (t) = sin (2 (ω t + φ)) / (2 ω)
- * </pre>
- * </p>
- *
- * <p>We can remove S between these expressions :
- * <pre>
- * If'2 (t) = a<sup>2</sup> ω<sup>2</sup> t - ω<sup>2</sup> If2 (t)
- * </pre>
- * </p>
- *
- * <p>The preceding expression shows that If'2 (t) is a linear
- * combination of both t and If2 (t): If'2 (t) = A × t + B × If2 (t)
- * </p>
- *
- * <p>From the primitive, we can deduce the same form for definite
- * integrals between t<sub>1</sub> and t<sub>i</sub> for each t<sub>i</sub> :
- * <pre>
- * If2 (t<sub>i</sub>) - If2 (t<sub>1</sub>) = A × (t<sub>i</sub> - t<sub>1</sub>) + B × (If2 (t<sub>i</sub>) - If2 (t<sub>1</sub>))
- * </pre>
- * </p>
- *
- * <p>We can find the coefficients A and B that best fit the sample
- * to this linear expression by computing the definite integrals for
- * each sample points.
- * </p>
- *
- * <p>For a bilinear expression z (x<sub>i</sub>, y<sub>i</sub>) = A × x<sub>i</sub> + B × y<sub>i</sub>, the
- * coefficients A and B that minimize a least square criterion
- * ∑ (z<sub>i</sub> - z (x<sub>i</sub>, y<sub>i</sub>))<sup>2</sup> are given by these expressions:</p>
- * <pre>
- *
- * ∑y<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>z<sub>i</sub> - ∑x<sub>i</sub>y<sub>i</sub> ∑y<sub>i</sub>z<sub>i</sub>
- * A = ------------------------
- * ∑x<sub>i</sub>x<sub>i</sub> ∑y<sub>i</sub>y<sub>i</sub> - ∑x<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>y<sub>i</sub>
- *
- * ∑x<sub>i</sub>x<sub>i</sub> ∑y<sub>i</sub>z<sub>i</sub> - ∑x<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>z<sub>i</sub>
- * B = ------------------------
- * ∑x<sub>i</sub>x<sub>i</sub> ∑y<sub>i</sub>y<sub>i</sub> - ∑x<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>y<sub>i</sub>
- * </pre>
- * </p>
- *
- *
- * <p>In fact, we can assume both a and ω are positive and
- * compute them directly, knowing that A = a<sup>2</sup> ω<sup>2</sup> and that
- * B = - ω<sup>2</sup>. The complete algorithm is therefore:</p>
- * <pre>
- *
- * for each t<sub>i</sub> from t<sub>1</sub> to t<sub>n-1</sub>, compute:
- * f (t<sub>i</sub>)
- * f' (t<sub>i</sub>) = (f (t<sub>i+1</sub>) - f(t<sub>i-1</sub>)) / (t<sub>i+1</sub> - t<sub>i-1</sub>)
- * x<sub>i</sub> = t<sub>i</sub> - t<sub>1</sub>
- * y<sub>i</sub> = ∫ f<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub>
- * z<sub>i</sub> = ∫ f'<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub>
- * update the sums ∑x<sub>i</sub>x<sub>i</sub>, ∑y<sub>i</sub>y<sub>i</sub>, ∑x<sub>i</sub>y<sub>i</sub>, ∑x<sub>i</sub>z<sub>i</sub> and ∑y<sub>i</sub>z<sub>i</sub>
- * end for
- *
- * |--------------------------
- * \ | ∑y<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>z<sub>i</sub> - ∑x<sub>i</sub>y<sub>i</sub> ∑y<sub>i</sub>z<sub>i</sub>
- * a = \ | ------------------------
- * \| ∑x<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>z<sub>i</sub> - ∑x<sub>i</sub>x<sub>i</sub> ∑y<sub>i</sub>z<sub>i</sub>
- *
- *
- * |--------------------------
- * \ | ∑x<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>z<sub>i</sub> - ∑x<sub>i</sub>x<sub>i</sub> ∑y<sub>i</sub>z<sub>i</sub>
- * ω = \ | ------------------------
- * \| ∑x<sub>i</sub>x<sub>i</sub> ∑y<sub>i</sub>y<sub>i</sub> - ∑x<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>y<sub>i</sub>
- *
- * </pre>
- * </p>
- *
- * <p>Once we know ω, we can compute:
- * <pre>
- * fc = ω f (t) cos (ω t) - f' (t) sin (ω t)
- * fs = ω f (t) sin (ω t) + f' (t) cos (ω t)
- * </pre>
- * </p>
- *
- * <p>It appears that <code>fc = a ω cos (φ)</code> and
- * <code>fs = -a ω sin (φ)</code>, so we can use these
- * expressions to compute φ. The best estimate over the sample is
- * given by averaging these expressions.
- * </p>
- *
- * <p>Since integrals and means are involved in the preceding
- * estimations, these operations run in O(n) time, where n is the
- * number of measurements.</p>
- */
- public static class ParameterGuesser {
- /** Amplitude. */
- private final double a;
- /** Angular frequency. */
- private final double omega;
- /** Phase. */
- private final double phi;
-
- /**
- * Simple constructor.
- *
- * @param observations Sampled observations.
- * @throws NumberIsTooSmallException if the sample is too short.
- * @throws ZeroException if the abscissa range is zero.
- * @throws MathIllegalStateException when the guessing procedure cannot
- * produce sensible results.
- */
- public ParameterGuesser(WeightedObservedPoint[] observations) {
- if (observations.length < 4) {
- throw new NumberIsTooSmallException(LocalizedFormats.INSUFFICIENT_OBSERVED_POINTS_IN_SAMPLE,
- observations.length, 4, true);
- }
-
- final WeightedObservedPoint[] sorted = sortObservations(observations);
-
- final double aOmega[] = guessAOmega(sorted);
- a = aOmega[0];
- omega = aOmega[1];
-
- phi = guessPhi(sorted);
- }
-
- /**
- * Gets an estimation of the parameters.
- *
- * @return the guessed parameters, in the following order:
- * <ul>
- * <li>Amplitude</li>
- * <li>Angular frequency</li>
- * <li>Phase</li>
- * </ul>
- */
- public double[] guess() {
- return new double[] { a, omega, phi };
- }
-
- /**
- * Sort the observations with respect to the abscissa.
- *
- * @param unsorted Input observations.
- * @return the input observations, sorted.
- */
- private WeightedObservedPoint[] sortObservations(WeightedObservedPoint[] unsorted) {
- final WeightedObservedPoint[] observations = unsorted.clone();
-
- // Since the samples are almost always already sorted, this
- // method is implemented as an insertion sort that reorders the
- // elements in place. Insertion sort is very efficient in this case.
- WeightedObservedPoint curr = observations[0];
- for (int j = 1; j < observations.length; ++j) {
- WeightedObservedPoint prec = curr;
- curr = observations[j];
- if (curr.getX() < prec.getX()) {
- // the current element should be inserted closer to the beginning
- int i = j - 1;
- WeightedObservedPoint mI = observations[i];
- while ((i >= 0) && (curr.getX() < mI.getX())) {
- observations[i + 1] = mI;
- if (i-- != 0) {
- mI = observations[i];
- }
- }
- observations[i + 1] = curr;
- curr = observations[j];
- }
- }
-
- return observations;
- }
-
- /**
- * Estimate a first guess of the amplitude and angular frequency.
- * This method assumes that the {@link #sortObservations(WeightedObservedPoint[])} method
- * has been called previously.
- *
- * @param observations Observations, sorted w.r.t. abscissa.
- * @throws ZeroException if the abscissa range is zero.
- * @throws MathIllegalStateException when the guessing procedure cannot
- * produce sensible results.
- * @return the guessed amplitude (at index 0) and circular frequency
- * (at index 1).
- */
- private double[] guessAOmega(WeightedObservedPoint[] observations) {
- final double[] aOmega = new double[2];
-
- // initialize the sums for the linear model between the two integrals
- double sx2 = 0;
- double sy2 = 0;
- double sxy = 0;
- double sxz = 0;
- double syz = 0;
-
- double currentX = observations[0].getX();
- double currentY = observations[0].getY();
- double f2Integral = 0;
- double fPrime2Integral = 0;
- final double startX = currentX;
- for (int i = 1; i < observations.length; ++i) {
- // one step forward
- final double previousX = currentX;
- final double previousY = currentY;
- currentX = observations[i].getX();
- currentY = observations[i].getY();
-
- // update the integrals of f<sup>2</sup> and f'<sup>2</sup>
- // considering a linear model for f (and therefore constant f')
- final double dx = currentX - previousX;
- final double dy = currentY - previousY;
- final double f2StepIntegral =
- dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3;
- final double fPrime2StepIntegral = dy * dy / dx;
-
- final double x = currentX - startX;
- f2Integral += f2StepIntegral;
- fPrime2Integral += fPrime2StepIntegral;
-
- sx2 += x * x;
- sy2 += f2Integral * f2Integral;
- sxy += x * f2Integral;
- sxz += x * fPrime2Integral;
- syz += f2Integral * fPrime2Integral;
- }
-
- // compute the amplitude and pulsation coefficients
- double c1 = sy2 * sxz - sxy * syz;
- double c2 = sxy * sxz - sx2 * syz;
- double c3 = sx2 * sy2 - sxy * sxy;
- if ((c1 / c2 < 0) || (c2 / c3 < 0)) {
- final int last = observations.length - 1;
- // Range of the observations, assuming that the
- // observations are sorted.
- final double xRange = observations[last].getX() - observations[0].getX();
- if (xRange == 0) {
- throw new ZeroException();
- }
- aOmega[1] = 2 * Math.PI / xRange;
-
- double yMin = Double.POSITIVE_INFINITY;
- double yMax = Double.NEGATIVE_INFINITY;
- for (int i = 1; i < observations.length; ++i) {
- final double y = observations[i].getY();
- if (y < yMin) {
- yMin = y;
- }
- if (y > yMax) {
- yMax = y;
- }
- }
- aOmega[0] = 0.5 * (yMax - yMin);
- } else {
- if (c2 == 0) {
- // In some ill-conditioned cases (cf. MATH-844), the guesser
- // procedure cannot produce sensible results.
- throw new MathIllegalStateException(LocalizedFormats.ZERO_DENOMINATOR);
- }
-
- aOmega[0] = FastMath.sqrt(c1 / c2);
- aOmega[1] = FastMath.sqrt(c2 / c3);
- }
-
- return aOmega;
- }
-
- /**
- * Estimate a first guess of the phase.
- *
- * @param observations Observations, sorted w.r.t. abscissa.
- * @return the guessed phase.
- */
- private double guessPhi(WeightedObservedPoint[] observations) {
- // initialize the means
- double fcMean = 0;
- double fsMean = 0;
-
- double currentX = observations[0].getX();
- double currentY = observations[0].getY();
- for (int i = 1; i < observations.length; ++i) {
- // one step forward
- final double previousX = currentX;
- final double previousY = currentY;
- currentX = observations[i].getX();
- currentY = observations[i].getY();
- final double currentYPrime = (currentY - previousY) / (currentX - previousX);
-
- double omegaX = omega * currentX;
- double cosine = FastMath.cos(omegaX);
- double sine = FastMath.sin(omegaX);
- fcMean += omega * currentY * cosine - currentYPrime * sine;
- fsMean += omega * currentY * sine + currentYPrime * cosine;
- }
-
- return FastMath.atan2(-fsMean, fcMean);
- }
- }
-}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/0a5cd113/src/main/java/org/apache/commons/math4/fitting/PolynomialFitter.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math4/fitting/PolynomialFitter.java b/src/main/java/org/apache/commons/math4/fitting/PolynomialFitter.java
deleted file mode 100644
index 38ebe91..0000000
--- a/src/main/java/org/apache/commons/math4/fitting/PolynomialFitter.java
+++ /dev/null
@@ -1,72 +0,0 @@
-/*
- * 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.math4.fitting;
-
-import org.apache.commons.math4.analysis.polynomials.PolynomialFunction;
-import org.apache.commons.math4.optim.nonlinear.vector.MultivariateVectorOptimizer;
-
-/**
- * Polynomial fitting is a very simple case of {@link CurveFitter curve fitting}.
- * The estimated coefficients are the polynomial coefficients (see the
- * {@link #fit(double[]) fit} method).
- *
- * @since 2.0
- * @deprecated As of 3.3. Please use {@link PolynomialCurveFitter} and
- * {@link WeightedObservedPoints} instead.
- */
-@Deprecated
-public class PolynomialFitter extends CurveFitter<PolynomialFunction.Parametric> {
- /**
- * Simple constructor.
- *
- * @param optimizer Optimizer to use for the fitting.
- */
- public PolynomialFitter(MultivariateVectorOptimizer optimizer) {
- super(optimizer);
- }
-
- /**
- * Get the coefficients of the polynomial fitting the weighted data points.
- * The degree of the fitting polynomial is {@code guess.length - 1}.
- *
- * @param guess First guess for the coefficients. They must be sorted in
- * increasing order of the polynomial's degree.
- * @param maxEval Maximum number of evaluations of the polynomial.
- * @return the coefficients of the polynomial that best fits the observed points.
- * @throws org.apache.commons.math4.exception.TooManyEvaluationsException if
- * the number of evaluations exceeds {@code maxEval}.
- * @throws org.apache.commons.math4.exception.ConvergenceException
- * if the algorithm failed to converge.
- */
- public double[] fit(int maxEval, double[] guess) {
- return fit(maxEval, new PolynomialFunction.Parametric(), guess);
- }
-
- /**
- * Get the coefficients of the polynomial fitting the weighted data points.
- * The degree of the fitting polynomial is {@code guess.length - 1}.
- *
- * @param guess First guess for the coefficients. They must be sorted in
- * increasing order of the polynomial's degree.
- * @return the coefficients of the polynomial that best fits the observed points.
- * @throws org.apache.commons.math4.exception.ConvergenceException
- * if the algorithm failed to converge.
- */
- public double[] fit(double[] guess) {
- return fit(new PolynomialFunction.Parametric(), guess);
- }
-}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/0a5cd113/src/test/java/org/apache/commons/math4/analysis/interpolation/BicubicSplineInterpolatingFunctionTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math4/analysis/interpolation/BicubicSplineInterpolatingFunctionTest.java b/src/test/java/org/apache/commons/math4/analysis/interpolation/BicubicSplineInterpolatingFunctionTest.java
deleted file mode 100644
index 92bf82b..0000000
--- a/src/test/java/org/apache/commons/math4/analysis/interpolation/BicubicSplineInterpolatingFunctionTest.java
+++ /dev/null
@@ -1,670 +0,0 @@
-/*
- * 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.math4.analysis.interpolation;
-
-import org.apache.commons.math4.analysis.BivariateFunction;
-import org.apache.commons.math4.analysis.interpolation.BicubicSplineFunction;
-import org.apache.commons.math4.analysis.interpolation.BicubicSplineInterpolatingFunction;
-import org.apache.commons.math4.distribution.UniformRealDistribution;
-import org.apache.commons.math4.exception.DimensionMismatchException;
-import org.apache.commons.math4.exception.MathIllegalArgumentException;
-import org.apache.commons.math4.exception.OutOfRangeException;
-import org.apache.commons.math4.random.RandomGenerator;
-import org.apache.commons.math4.random.Well19937c;
-import org.junit.Assert;
-import org.junit.Test;
-import org.junit.Ignore;
-
-/**
- * Test case for the bicubic function.
- *
- */
-public final class BicubicSplineInterpolatingFunctionTest {
- /**
- * Test preconditions.
- */
- @Test
- public void testPreconditions() {
- double[] xval = new double[] {3, 4, 5, 6.5};
- double[] yval = new double[] {-4, -3, -1, 2.5};
- double[][] zval = new double[xval.length][yval.length];
-
- @SuppressWarnings("unused")
- BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval,
- zval, zval, zval);
-
- double[] wxval = new double[] {3, 2, 5, 6.5};
- try {
- bcf = new BicubicSplineInterpolatingFunction(wxval, yval, zval, zval, zval, zval);
- Assert.fail("an exception should have been thrown");
- } catch (MathIllegalArgumentException e) {
- // Expected
- }
- double[] wyval = new double[] {-4, -1, -1, 2.5};
- try {
- bcf = new BicubicSplineInterpolatingFunction(xval, wyval, zval, zval, zval, zval);
- Assert.fail("an exception should have been thrown");
- } catch (MathIllegalArgumentException e) {
- // Expected
- }
- double[][] wzval = new double[xval.length][yval.length - 1];
- try {
- bcf = new BicubicSplineInterpolatingFunction(xval, yval, wzval, zval, zval, zval);
- Assert.fail("an exception should have been thrown");
- } catch (DimensionMismatchException e) {
- // Expected
- }
- try {
- bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, wzval, zval, zval);
- Assert.fail("an exception should have been thrown");
- } catch (DimensionMismatchException e) {
- // Expected
- }
- try {
- bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, zval, wzval, zval);
- Assert.fail("an exception should have been thrown");
- } catch (DimensionMismatchException e) {
- // Expected
- }
- try {
- bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, zval, zval, wzval);
- Assert.fail("an exception should have been thrown");
- } catch (DimensionMismatchException e) {
- // Expected
- }
-
- wzval = new double[xval.length - 1][yval.length];
- try {
- bcf = new BicubicSplineInterpolatingFunction(xval, yval, wzval, zval, zval, zval);
- Assert.fail("an exception should have been thrown");
- } catch (DimensionMismatchException e) {
- // Expected
- }
- try {
- bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, wzval, zval, zval);
- Assert.fail("an exception should have been thrown");
- } catch (DimensionMismatchException e) {
- // Expected
- }
- try {
- bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, zval, wzval, zval);
- Assert.fail("an exception should have been thrown");
- } catch (DimensionMismatchException e) {
- // Expected
- }
- try {
- bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval, zval, zval, wzval);
- Assert.fail("an exception should have been thrown");
- } catch (DimensionMismatchException e) {
- // Expected
- }
- }
-
- /**
- * Test for a plane.
- * <p>
- * z = 2 x - 3 y + 5
- */
- @Ignore@Test
- public void testPlane() {
- double[] xval = new double[] {3, 4, 5, 6.5};
- double[] yval = new double[] {-4, -3, -1, 2, 2.5};
- // Function values
- BivariateFunction f = new BivariateFunction() {
- public double value(double x, double y) {
- return 2 * x - 3 * y + 5;
- }
- };
- double[][] zval = new double[xval.length][yval.length];
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- zval[i][j] = f.value(xval[i], yval[j]);
- }
- }
- // Partial derivatives with respect to x
- double[][] dZdX = new double[xval.length][yval.length];
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- dZdX[i][j] = 2;
- }
- }
- // Partial derivatives with respect to y
- double[][] dZdY = new double[xval.length][yval.length];
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- dZdY[i][j] = -3;
- }
- }
- // Partial cross-derivatives
- double[][] dZdXdY = new double[xval.length][yval.length];
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- dZdXdY[i][j] = 0;
- }
- }
-
- BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval,
- dZdX, dZdY, dZdXdY);
- double x, y;
- double expected, result;
-
- x = 4;
- y = -3;
- expected = f.value(x, y);
- result = bcf.value(x, y);
- Assert.assertEquals("On sample point",
- expected, result, 1e-15);
-
- x = 4.5;
- y = -1.5;
- expected = f.value(x, y);
- result = bcf.value(x, y);
- Assert.assertEquals("Half-way between sample points (middle of the patch)",
- expected, result, 0.3);
-
- x = 3.5;
- y = -3.5;
- expected = f.value(x, y);
- result = bcf.value(x, y);
- Assert.assertEquals("Half-way between sample points (border of the patch)",
- expected, result, 0.3);
- }
-
- /**
- * Test for a paraboloid.
- * <p>
- * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
- */
- @Ignore@Test
- public void testParaboloid() {
- double[] xval = new double[] {3, 4, 5, 6.5};
- double[] yval = new double[] {-4, -3, -1, 2, 2.5};
- // Function values
- BivariateFunction f = new BivariateFunction() {
- public double value(double x, double y) {
- return 2 * x * x - 3 * y * y + 4 * x * y - 5;
- }
- };
- double[][] zval = new double[xval.length][yval.length];
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- zval[i][j] = f.value(xval[i], yval[j]);
- }
- }
- // Partial derivatives with respect to x
- double[][] dZdX = new double[xval.length][yval.length];
- BivariateFunction dfdX = new BivariateFunction() {
- public double value(double x, double y) {
- return 4 * (x + y);
- }
- };
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- dZdX[i][j] = dfdX.value(xval[i], yval[j]);
- }
- }
- // Partial derivatives with respect to y
- double[][] dZdY = new double[xval.length][yval.length];
- BivariateFunction dfdY = new BivariateFunction() {
- public double value(double x, double y) {
- return 4 * x - 6 * y;
- }
- };
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- dZdY[i][j] = dfdY.value(xval[i], yval[j]);
- }
- }
- // Partial cross-derivatives
- double[][] dZdXdY = new double[xval.length][yval.length];
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- dZdXdY[i][j] = 4;
- }
- }
-
- BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval,
- dZdX, dZdY, dZdXdY);
- double x, y;
- double expected, result;
-
- x = 4;
- y = -3;
- expected = f.value(x, y);
- result = bcf.value(x, y);
- Assert.assertEquals("On sample point",
- expected, result, 1e-15);
-
- x = 4.5;
- y = -1.5;
- expected = f.value(x, y);
- result = bcf.value(x, y);
- Assert.assertEquals("Half-way between sample points (middle of the patch)",
- expected, result, 2);
-
- x = 3.5;
- y = -3.5;
- expected = f.value(x, y);
- result = bcf.value(x, y);
- Assert.assertEquals("Half-way between sample points (border of the patch)",
- expected, result, 2);
- }
-
- /**
- * Test for partial derivatives of {@link BicubicSplineFunction}.
- * <p>
- * f(x, y) = Σ<sub>i</sub>Σ<sub>j</sub> (i+1) (j+2) x<sup>i</sup> y<sup>j</sup>
- */
- @Ignore@Test
- public void testSplinePartialDerivatives() {
- final int N = 4;
- final double[] coeff = new double[16];
-
- for (int i = 0; i < N; i++) {
- for (int j = 0; j < N; j++) {
- coeff[i + N * j] = (i + 1) * (j + 2);
- }
- }
-
- final BicubicSplineFunction f = new BicubicSplineFunction(coeff);
- BivariateFunction derivative;
- final double x = 0.435;
- final double y = 0.776;
- final double tol = 1e-13;
-
- derivative = new BivariateFunction() {
- public double value(double x, double y) {
- final double x2 = x * x;
- final double y2 = y * y;
- final double y3 = y2 * y;
- final double yFactor = 2 + 3 * y + 4 * y2 + 5 * y3;
- return yFactor * (2 + 6 * x + 12 * x2);
- }
- };
- Assert.assertEquals("dFdX", derivative.value(x, y),
- f.partialDerivativeX().value(x, y), tol);
-
- derivative = new BivariateFunction() {
- public double value(double x, double y) {
- final double x2 = x * x;
- final double x3 = x2 * x;
- final double y2 = y * y;
- final double xFactor = 1 + 2 * x + 3 * x2 + 4 * x3;
- return xFactor * (3 + 8 * y + 15 * y2);
- }
- };
- Assert.assertEquals("dFdY", derivative.value(x, y),
- f.partialDerivativeY().value(x, y), tol);
-
- derivative = new BivariateFunction() {
- public double value(double x, double y) {
- final double y2 = y * y;
- final double y3 = y2 * y;
- final double yFactor = 2 + 3 * y + 4 * y2 + 5 * y3;
- return yFactor * (6 + 24 * x);
- }
- };
- Assert.assertEquals("d2FdX2", derivative.value(x, y),
- f.partialDerivativeXX().value(x, y), tol);
-
- derivative = new BivariateFunction() {
- public double value(double x, double y) {
- final double x2 = x * x;
- final double x3 = x2 * x;
- final double xFactor = 1 + 2 * x + 3 * x2 + 4 * x3;
- return xFactor * (8 + 30 * y);
- }
- };
- Assert.assertEquals("d2FdY2", derivative.value(x, y),
- f.partialDerivativeYY().value(x, y), tol);
-
- derivative = new BivariateFunction() {
- public double value(double x, double y) {
- final double x2 = x * x;
- final double y2 = y * y;
- final double yFactor = 3 + 8 * y + 15 * y2;
- return yFactor * (2 + 6 * x + 12 * x2);
- }
- };
- Assert.assertEquals("d2FdXdY", derivative.value(x, y),
- f.partialDerivativeXY().value(x, y), tol);
- }
-
- /**
- * Test that the partial derivatives computed from a
- * {@link BicubicSplineInterpolatingFunction} match the input data.
- * <p>
- * f(x, y) = 5
- * - 3 x + 2 y
- * - x y + 2 x<sup>2</sup> - 3 y<sup>2</sup>
- * + 4 x<sup>2</sup> y - x y<sup>2</sup> - 3 x<sup>3</sup> + y<sup>3</sup>
- */
- @Ignore@Test
- public void testMatchingPartialDerivatives() {
- final int sz = 21;
- double[] val = new double[sz];
- // Coordinate values
- final double delta = 1d / (sz - 1);
- for (int i = 0; i < sz; i++) {
- val[i] = i * delta;
- }
- // Function values
- BivariateFunction f = new BivariateFunction() {
- public double value(double x, double y) {
- final double x2 = x * x;
- final double x3 = x2 * x;
- final double y2 = y * y;
- final double y3 = y2 * y;
-
- return 5
- - 3 * x + 2 * y
- - x * y + 2 * x2 - 3 * y2
- + 4 * x2 * y - x * y2 - 3 * x3 + y3;
- }
- };
- double[][] fval = new double[sz][sz];
- for (int i = 0; i < sz; i++) {
- for (int j = 0; j < sz; j++) {
- fval[i][j] = f.value(val[i], val[j]);
- }
- }
- // Partial derivatives with respect to x
- double[][] dFdX = new double[sz][sz];
- BivariateFunction dfdX = new BivariateFunction() {
- public double value(double x, double y) {
- final double x2 = x * x;
- final double y2 = y * y;
- return - 3 - y + 4 * x + 8 * x * y - y2 - 9 * x2;
- }
- };
- for (int i = 0; i < sz; i++) {
- for (int j = 0; j < sz; j++) {
- dFdX[i][j] = dfdX.value(val[i], val[j]);
- }
- }
- // Partial derivatives with respect to y
- double[][] dFdY = new double[sz][sz];
- BivariateFunction dfdY = new BivariateFunction() {
- public double value(double x, double y) {
- final double x2 = x * x;
- final double y2 = y * y;
- return 2 - x - 6 * y + 4 * x2 - 2 * x * y + 3 * y2;
- }
- };
- for (int i = 0; i < sz; i++) {
- for (int j = 0; j < sz; j++) {
- dFdY[i][j] = dfdY.value(val[i], val[j]);
- }
- }
- // Partial cross-derivatives
- double[][] d2FdXdY = new double[sz][sz];
- BivariateFunction d2fdXdY = new BivariateFunction() {
- public double value(double x, double y) {
- return -1 + 8 * x - 2 * y;
- }
- };
- for (int i = 0; i < sz; i++) {
- for (int j = 0; j < sz; j++) {
- d2FdXdY[i][j] = d2fdXdY.value(val[i], val[j]);
- }
- }
-
- BicubicSplineInterpolatingFunction bcf
- = new BicubicSplineInterpolatingFunction(val, val, fval, dFdX, dFdY, d2FdXdY);
-
- double x, y;
- double expected, result;
-
- final double tol = 1e-12;
- for (int i = 0; i < sz; i++) {
- x = val[i];
- for (int j = 0; j < sz; j++) {
- y = val[j];
-
- expected = dfdX.value(x, y);
- result = bcf.partialDerivativeX(x, y);
- Assert.assertEquals(x + " " + y + " dFdX", expected, result, tol);
-
- expected = dfdY.value(x, y);
- result = bcf.partialDerivativeY(x, y);
- Assert.assertEquals(x + " " + y + " dFdY", expected, result, tol);
-
- expected = d2fdXdY.value(x, y);
- result = bcf.partialDerivativeXY(x, y);
- Assert.assertEquals(x + " " + y + " d2FdXdY", expected, result, tol);
- }
- }
- }
-
- /**
- * Interpolating a plane.
- * <p>
- * z = 2 x - 3 y + 5
- */
- @Test
- public void testInterpolation1() {
- final int sz = 21;
- double[] xval = new double[sz];
- double[] yval = new double[sz];
- // Coordinate values
- final double delta = 1d / (sz - 1);
- for (int i = 0; i < sz; i++) {
- xval[i] = -1 + 15 * i * delta;
- yval[i] = -20 + 30 * i * delta;
- }
-
- // Function values
- BivariateFunction f = new BivariateFunction() {
- public double value(double x, double y) {
- return 2 * x - 3 * y + 5;
- }
- };
- double[][] zval = new double[xval.length][yval.length];
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- zval[i][j] = f.value(xval[i], yval[j]);
- }
- }
- // Partial derivatives with respect to x
- double[][] dZdX = new double[xval.length][yval.length];
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- dZdX[i][j] = 2;
- }
- }
- // Partial derivatives with respect to y
- double[][] dZdY = new double[xval.length][yval.length];
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- dZdY[i][j] = -3;
- }
- }
- // Partial cross-derivatives
- double[][] dZdXdY = new double[xval.length][yval.length];
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- dZdXdY[i][j] = 0;
- }
- }
-
- final BivariateFunction bcf
- = new BicubicSplineInterpolatingFunction(xval, yval, zval,
- dZdX, dZdY, dZdXdY);
- double x, y;
-
- final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
- final UniformRealDistribution distX
- = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
- final UniformRealDistribution distY
- = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);
-
- final int numSamples = 50;
- final double tol = 6;
- for (int i = 0; i < numSamples; i++) {
- x = distX.sample();
- for (int j = 0; j < numSamples; j++) {
- y = distY.sample();
-// System.out.println(x + " " + y + " " + f.value(x, y) + " " + bcf.value(x, y));
- Assert.assertEquals(f.value(x, y), bcf.value(x, y), tol);
- }
-// System.out.println();
- }
- }
-
- /**
- * Interpolating a paraboloid.
- * <p>
- * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
- */
- @Test
- public void testInterpolation2() {
- final int sz = 21;
- double[] xval = new double[sz];
- double[] yval = new double[sz];
- // Coordinate values
- final double delta = 1d / (sz - 1);
- for (int i = 0; i < sz; i++) {
- xval[i] = -1 + 15 * i * delta;
- yval[i] = -20 + 30 * i * delta;
- }
-
- // Function values
- BivariateFunction f = new BivariateFunction() {
- public double value(double x, double y) {
- return 2 * x * x - 3 * y * y + 4 * x * y - 5;
- }
- };
- double[][] zval = new double[xval.length][yval.length];
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- zval[i][j] = f.value(xval[i], yval[j]);
- }
- }
- // Partial derivatives with respect to x
- double[][] dZdX = new double[xval.length][yval.length];
- BivariateFunction dfdX = new BivariateFunction() {
- public double value(double x, double y) {
- return 4 * (x + y);
- }
- };
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- dZdX[i][j] = dfdX.value(xval[i], yval[j]);
- }
- }
- // Partial derivatives with respect to y
- double[][] dZdY = new double[xval.length][yval.length];
- BivariateFunction dfdY = new BivariateFunction() {
- public double value(double x, double y) {
- return 4 * x - 6 * y;
- }
- };
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- dZdY[i][j] = dfdY.value(xval[i], yval[j]);
- }
- }
- // Partial cross-derivatives
- double[][] dZdXdY = new double[xval.length][yval.length];
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- dZdXdY[i][j] = 4;
- }
- }
-
- BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval,
- dZdX, dZdY, dZdXdY);
- double x, y;
-
- final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
- final UniformRealDistribution distX
- = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
- final UniformRealDistribution distY
- = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);
-
- final double tol = 224;
- for (int i = 0; i < sz; i++) {
- x = distX.sample();
- for (int j = 0; j < sz; j++) {
- y = distY.sample();
-// System.out.println(x + " " + y + " " + f.value(x, y) + " " + bcf.value(x, y));
- Assert.assertEquals(f.value(x, y), bcf.value(x, y), tol);
- }
-// System.out.println();
- }
- }
-
- @Test
- public void testIsValidPoint() {
- final double xMin = -12;
- final double xMax = 34;
- final double yMin = 5;
- final double yMax = 67;
- final double[] xval = new double[] { xMin, xMax };
- final double[] yval = new double[] { yMin, yMax };
- final double[][] f = new double[][] { { 1, 2 },
- { 3, 4 } };
- final double[][] dFdX = f;
- final double[][] dFdY = f;
- final double[][] dFdXdY = f;
-
- final BicubicSplineInterpolatingFunction bcf
- = new BicubicSplineInterpolatingFunction(xval, yval, f,
- dFdX, dFdY, dFdXdY);
-
- double x, y;
-
- x = xMin;
- y = yMin;
- Assert.assertTrue(bcf.isValidPoint(x, y));
- // Ensure that no exception is thrown.
- bcf.value(x, y);
-
- x = xMax;
- y = yMax;
- Assert.assertTrue(bcf.isValidPoint(x, y));
- // Ensure that no exception is thrown.
- bcf.value(x, y);
-
- final double xRange = xMax - xMin;
- final double yRange = yMax - yMin;
- x = xMin + xRange / 3.4;
- y = yMin + yRange / 1.2;
- Assert.assertTrue(bcf.isValidPoint(x, y));
- // Ensure that no exception is thrown.
- bcf.value(x, y);
-
- final double small = 1e-8;
- x = xMin - small;
- y = yMax;
- Assert.assertFalse(bcf.isValidPoint(x, y));
- // Ensure that an exception would have been thrown.
- try {
- bcf.value(x, y);
- Assert.fail("OutOfRangeException expected");
- } catch (OutOfRangeException expected) {}
-
- x = xMin;
- y = yMax + small;
- Assert.assertFalse(bcf.isValidPoint(x, y));
- // Ensure that an exception would have been thrown.
- try {
- bcf.value(x, y);
- Assert.fail("OutOfRangeException expected");
- } catch (OutOfRangeException expected) {}
- }
-}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/0a5cd113/src/test/java/org/apache/commons/math4/analysis/interpolation/BicubicSplineInterpolatorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math4/analysis/interpolation/BicubicSplineInterpolatorTest.java b/src/test/java/org/apache/commons/math4/analysis/interpolation/BicubicSplineInterpolatorTest.java
deleted file mode 100644
index 91f1f66..0000000
--- a/src/test/java/org/apache/commons/math4/analysis/interpolation/BicubicSplineInterpolatorTest.java
+++ /dev/null
@@ -1,186 +0,0 @@
-/*
- * 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.math4.analysis.interpolation;
-
-import org.apache.commons.math4.analysis.BivariateFunction;
-import org.apache.commons.math4.analysis.interpolation.BicubicSplineInterpolator;
-import org.apache.commons.math4.analysis.interpolation.BivariateGridInterpolator;
-import org.apache.commons.math4.distribution.UniformRealDistribution;
-import org.apache.commons.math4.exception.DimensionMismatchException;
-import org.apache.commons.math4.exception.MathIllegalArgumentException;
-import org.apache.commons.math4.random.RandomGenerator;
-import org.apache.commons.math4.random.Well19937c;
-import org.junit.Assert;
-import org.junit.Test;
-
-/**
- * Test case for the bicubic interpolator.
- *
- */
-public final class BicubicSplineInterpolatorTest {
- /**
- * Test preconditions.
- */
- @Test
- public void testPreconditions() {
- double[] xval = new double[] {3, 4, 5, 6.5};
- double[] yval = new double[] {-4, -3, -1, 2.5};
- double[][] zval = new double[xval.length][yval.length];
-
- BivariateGridInterpolator interpolator = new BicubicSplineInterpolator();
-
- @SuppressWarnings("unused")
- BivariateFunction p = interpolator.interpolate(xval, yval, zval);
-
- double[] wxval = new double[] {3, 2, 5, 6.5};
- try {
- p = interpolator.interpolate(wxval, yval, zval);
- Assert.fail("an exception should have been thrown");
- } catch (MathIllegalArgumentException e) {
- // Expected
- }
-
- double[] wyval = new double[] {-4, -3, -1, -1};
- try {
- p = interpolator.interpolate(xval, wyval, zval);
- Assert.fail("an exception should have been thrown");
- } catch (MathIllegalArgumentException e) {
- // Expected
- }
-
- double[][] wzval = new double[xval.length][yval.length + 1];
- try {
- p = interpolator.interpolate(xval, yval, wzval);
- Assert.fail("an exception should have been thrown");
- } catch (DimensionMismatchException e) {
- // Expected
- }
- wzval = new double[xval.length - 1][yval.length];
- try {
- p = interpolator.interpolate(xval, yval, wzval);
- Assert.fail("an exception should have been thrown");
- } catch (DimensionMismatchException e) {
- // Expected
- }
- }
-
- /**
- * Interpolating a plane.
- * <p>
- * z = 2 x - 3 y + 5
- */
- @Test
- public void testInterpolation1() {
- final int sz = 21;
- double[] xval = new double[sz];
- double[] yval = new double[sz];
- // Coordinate values
- final double delta = 1d / (sz - 1);
- for (int i = 0; i < sz; i++) {
- xval[i] = -1 + 15 * i * delta;
- yval[i] = -20 + 30 * i * delta;
- }
-
- // Function values
- BivariateFunction f = new BivariateFunction() {
- public double value(double x, double y) {
- return 2 * x - 3 * y + 5;
- }
- };
- double[][] zval = new double[xval.length][yval.length];
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- zval[i][j] = f.value(xval[i], yval[j]);
- }
- }
-
- BivariateGridInterpolator interpolator = new BicubicSplineInterpolator();
- BivariateFunction p = interpolator.interpolate(xval, yval, zval);
- double x, y;
-
- final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
- final UniformRealDistribution distX
- = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
- final UniformRealDistribution distY
- = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);
-
- final int numSamples = 50;
- final double tol = 6;
- for (int i = 0; i < numSamples; i++) {
- x = distX.sample();
- for (int j = 0; j < numSamples; j++) {
- y = distY.sample();
-// System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
- Assert.assertEquals(f.value(x, y), p.value(x, y), tol);
- }
-// System.out.println();
- }
- }
-
- /**
- * Interpolating a paraboloid.
- * <p>
- * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
- */
- @Test
- public void testInterpolation2() {
- final int sz = 21;
- double[] xval = new double[sz];
- double[] yval = new double[sz];
- // Coordinate values
- final double delta = 1d / (sz - 1);
- for (int i = 0; i < sz; i++) {
- xval[i] = -1 + 15 * i * delta;
- yval[i] = -20 + 30 * i * delta;
- }
-
- // Function values
- BivariateFunction f = new BivariateFunction() {
- public double value(double x, double y) {
- return 2 * x * x - 3 * y * y + 4 * x * y - 5;
- }
- };
- double[][] zval = new double[xval.length][yval.length];
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- zval[i][j] = f.value(xval[i], yval[j]);
- }
- }
-
- BivariateGridInterpolator interpolator = new BicubicSplineInterpolator();
- BivariateFunction p = interpolator.interpolate(xval, yval, zval);
- double x, y;
-
- final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
- final UniformRealDistribution distX
- = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
- final UniformRealDistribution distY
- = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);
-
- final int numSamples = 50;
- final double tol = 251;
- for (int i = 0; i < numSamples; i++) {
- x = distX.sample();
- for (int j = 0; j < numSamples; j++) {
- y = distY.sample();
-// System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
- Assert.assertEquals(f.value(x, y), p.value(x, y), tol);
- }
-// System.out.println();
- }
- }
-}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/0a5cd113/src/test/java/org/apache/commons/math4/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolatorTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math4/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolatorTest.java b/src/test/java/org/apache/commons/math4/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolatorTest.java
deleted file mode 100644
index 9cdd888..0000000
--- a/src/test/java/org/apache/commons/math4/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolatorTest.java
+++ /dev/null
@@ -1,181 +0,0 @@
-/*
- * 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.math4.analysis.interpolation;
-
-import org.apache.commons.math4.analysis.BivariateFunction;
-import org.apache.commons.math4.analysis.interpolation.BivariateGridInterpolator;
-import org.apache.commons.math4.analysis.interpolation.SmoothingPolynomialBicubicSplineInterpolator;
-import org.apache.commons.math4.exception.DimensionMismatchException;
-import org.apache.commons.math4.exception.MathIllegalArgumentException;
-import org.apache.commons.math4.util.FastMath;
-import org.junit.Assert;
-import org.junit.Test;
-
-/**
- * Test case for the smoothing bicubic interpolator.
- *
- */
-public final class SmoothingPolynomialBicubicSplineInterpolatorTest {
- /**
- * Test preconditions.
- */
- @Test
- public void testPreconditions() {
- double[] xval = new double[] {3, 4, 5, 6.5, 7.5};
- double[] yval = new double[] {-4, -3, -1, 2.5, 3};
- double[][] zval = new double[xval.length][yval.length];
-
- BivariateGridInterpolator interpolator = new SmoothingPolynomialBicubicSplineInterpolator(0);
-
- @SuppressWarnings("unused")
- BivariateFunction p = interpolator.interpolate(xval, yval, zval);
-
- double[] wxval = new double[] {3, 2, 5, 6.5};
- try {
- p = interpolator.interpolate(wxval, yval, zval);
- Assert.fail("an exception should have been thrown");
- } catch (MathIllegalArgumentException e) {
- // Expected
- }
-
- double[] wyval = new double[] {-4, -3, -1, -1};
- try {
- p = interpolator.interpolate(xval, wyval, zval);
- Assert.fail("an exception should have been thrown");
- } catch (MathIllegalArgumentException e) {
- // Expected
- }
-
- double[][] wzval = new double[xval.length][yval.length + 1];
- try {
- p = interpolator.interpolate(xval, yval, wzval);
- Assert.fail("an exception should have been thrown");
- } catch (DimensionMismatchException e) {
- // Expected
- }
- wzval = new double[xval.length - 1][yval.length];
- try {
- p = interpolator.interpolate(xval, yval, wzval);
- Assert.fail("an exception should have been thrown");
- } catch (DimensionMismatchException e) {
- // Expected
- }
- wzval = new double[xval.length][yval.length - 1];
- try {
- p = interpolator.interpolate(xval, yval, wzval);
- Assert.fail("an exception should have been thrown");
- } catch (DimensionMismatchException e) {
- // Expected
- }
- }
-
- /**
- * Test of interpolator for a plane.
- * <p>
- * z = 2 x - 3 y + 5
- */
- @Test
- public void testPlane() {
- BivariateFunction f = new BivariateFunction() {
- public double value(double x, double y) {
- return 2 * x - 3 * y + 5
- + ((int) (FastMath.abs(5 * x + 3 * y)) % 2 == 0 ? 1 : -1);
- }
- };
-
- BivariateGridInterpolator interpolator = new SmoothingPolynomialBicubicSplineInterpolator(1);
-
- double[] xval = new double[] {3, 4, 5, 6.5, 7.5};
- double[] yval = new double[] {-4, -3, -1, 2, 2.5, 3.5};
- double[][] zval = new double[xval.length][yval.length];
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- zval[i][j] = f.value(xval[i], yval[j]);
- }
- }
-
- BivariateFunction p = interpolator.interpolate(xval, yval, zval);
- double x, y;
- double expected, result;
-
- x = 4;
- y = -3;
- expected = f.value(x, y);
- result = p.value(x, y);
- Assert.assertEquals("On sample point", expected, result, 2);
-
- x = 4.5;
- y = -1.5;
- expected = f.value(x, y);
- result = p.value(x, y);
- Assert.assertEquals("half-way between sample points (middle of the patch)", expected, result, 2);
-
- x = 3.5;
- y = -3.5;
- expected = f.value(x, y);
- result = p.value(x, y);
- Assert.assertEquals("half-way between sample points (border of the patch)", expected, result, 2);
- }
-
- /**
- * Test of interpolator for a paraboloid.
- * <p>
- * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
- */
- @Test
- public void testParaboloid() {
- BivariateFunction f = new BivariateFunction() {
- public double value(double x, double y) {
- return 2 * x * x - 3 * y * y + 4 * x * y - 5
- + ((int) (FastMath.abs(5 * x + 3 * y)) % 2 == 0 ? 1 : -1);
- }
- };
-
- BivariateGridInterpolator interpolator = new SmoothingPolynomialBicubicSplineInterpolator(4);
-
- double[] xval = new double[] {3, 4, 5, 6.5, 7.5, 8};
- double[] yval = new double[] {-4, -3, -2, -1, 0.5, 2.5};
- double[][] zval = new double[xval.length][yval.length];
- for (int i = 0; i < xval.length; i++) {
- for (int j = 0; j < yval.length; j++) {
- zval[i][j] = f.value(xval[i], yval[j]);
- }
- }
-
- BivariateFunction p = interpolator.interpolate(xval, yval, zval);
- double x, y;
- double expected, result;
-
- x = 5;
- y = 0.5;
- expected = f.value(x, y);
- result = p.value(x, y);
- Assert.assertEquals("On sample point", expected, result, 2);
-
- x = 4.5;
- y = -1.5;
- expected = f.value(x, y);
- result = p.value(x, y);
- Assert.assertEquals("half-way between sample points (middle of the patch)", expected, result, 2);
-
- x = 3.5;
- y = -3.5;
- expected = f.value(x, y);
- result = p.value(x, y);
- Assert.assertEquals("half-way between sample points (border of the patch)", expected, result, 2);
- }
-}