diff --git a/src/java/org/apache/commons/math/estimation/AbstractEstimator.java b/src/java/org/apache/commons/math/estimation/AbstractEstimator.java new file mode 100644 index 000000000..d33e1e5a5 --- /dev/null +++ b/src/java/org/apache/commons/math/estimation/AbstractEstimator.java @@ -0,0 +1,314 @@ +/* + * 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.estimation; + +import java.util.Arrays; + +import org.apache.commons.math.linear.InvalidMatrixException; +import org.apache.commons.math.linear.MatrixUtils; +import org.apache.commons.math.linear.RealMatrix; +import org.apache.commons.math.linear.decomposition.LUDecompositionImpl; + +/** + * Base class for implementing estimators. + *

This base class handles the boilerplates methods associated to thresholds + * settings, jacobian and error estimation.

+ * @version $Revision$ $Date$ + * @since 1.2 + * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has + * been deprecated and replaced by package org.apache.commons.math.optimization.general + * + */ +@Deprecated +public abstract class AbstractEstimator implements Estimator { + + /** Default maximal number of cost evaluations allowed. */ + public static final int DEFAULT_MAX_COST_EVALUATIONS = 100; + + /** + * Build an abstract estimator for least squares problems. + *

The maximal number of cost evaluations allowed is set + * to its default value {@link #DEFAULT_MAX_COST_EVALUATIONS}.

+ */ + protected AbstractEstimator() { + setMaxCostEval(DEFAULT_MAX_COST_EVALUATIONS); + } + + /** + * Set the maximal number of cost evaluations allowed. + * + * @param maxCostEval maximal number of cost evaluations allowed + * @see #estimate + */ + public final void setMaxCostEval(int maxCostEval) { + this.maxCostEval = maxCostEval; + } + + /** + * Get the number of cost evaluations. + * + * @return number of cost evaluations + * */ + public final int getCostEvaluations() { + return costEvaluations; + } + + /** + * Get the number of jacobian evaluations. + * + * @return number of jacobian evaluations + * */ + public final int getJacobianEvaluations() { + return jacobianEvaluations; + } + + /** + * Update the jacobian matrix. + */ + protected void updateJacobian() { + incrementJacobianEvaluationsCounter(); + Arrays.fill(jacobian, 0); + for (int i = 0, index = 0; i < rows; i++) { + WeightedMeasurement wm = measurements[i]; + double factor = -Math.sqrt(wm.getWeight()); + for (int j = 0; j < cols; ++j) { + jacobian[index++] = factor * wm.getPartial(parameters[j]); + } + } + } + + /** + * Increment the jacobian evaluations counter. + */ + protected final void incrementJacobianEvaluationsCounter() { + ++jacobianEvaluations; + } + + /** + * Update the residuals array and cost function value. + * @exception EstimationException if the number of cost evaluations + * exceeds the maximum allowed + */ + protected void updateResidualsAndCost() + throws EstimationException { + + if (++costEvaluations > maxCostEval) { + throw new EstimationException("maximal number of evaluations exceeded ({0})", + maxCostEval); + } + + cost = 0; + for (int i = 0, index = 0; i < rows; i++, index += cols) { + WeightedMeasurement wm = measurements[i]; + double residual = wm.getResidual(); + residuals[i] = Math.sqrt(wm.getWeight()) * residual; + cost += wm.getWeight() * residual * residual; + } + cost = Math.sqrt(cost); + + } + + /** + * Get the Root Mean Square value. + * Get the Root Mean Square value, i.e. the root of the arithmetic + * mean of the square of all weighted residuals. This is related to the + * criterion that is minimized by the estimator as follows: if + * c if the criterion, and n is the number of + * measurements, then the RMS is sqrt (c/n). + * + * @param problem estimation problem + * @return RMS value + */ + public double getRMS(EstimationProblem problem) { + WeightedMeasurement[] wm = problem.getMeasurements(); + double criterion = 0; + for (int i = 0; i < wm.length; ++i) { + double residual = wm[i].getResidual(); + criterion += wm[i].getWeight() * residual * residual; + } + return Math.sqrt(criterion / wm.length); + } + + /** + * Get the Chi-Square value. + * @param problem estimation problem + * @return chi-square value + */ + public double getChiSquare(EstimationProblem problem) { + WeightedMeasurement[] wm = problem.getMeasurements(); + double chiSquare = 0; + for (int i = 0; i < wm.length; ++i) { + double residual = wm[i].getResidual(); + chiSquare += residual * residual / wm[i].getWeight(); + } + return chiSquare; + } + + /** + * Get the covariance matrix of unbound estimated parameters. + * @param problem estimation problem + * @return covariance matrix + * @exception EstimationException if the covariance matrix + * cannot be computed (singular problem) + */ + public double[][] getCovariances(EstimationProblem problem) + throws EstimationException { + + // set up the jacobian + updateJacobian(); + + // compute transpose(J).J, avoiding building big intermediate matrices + final int rows = problem.getMeasurements().length; + final int cols = problem.getUnboundParameters().length; + final int max = cols * rows; + double[][] jTj = new double[cols][cols]; + for (int i = 0; i < cols; ++i) { + for (int j = i; j < cols; ++j) { + double sum = 0; + for (int k = 0; k < max; k += cols) { + sum += jacobian[k + i] * jacobian[k + j]; + } + jTj[i][j] = sum; + jTj[j][i] = sum; + } + } + + try { + // compute the covariances matrix + RealMatrix inverse = + new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse(); + return inverse.getData(); + } catch (InvalidMatrixException ime) { + throw new EstimationException("unable to compute covariances: singular problem"); + } + + } + + /** + * Guess the errors in unbound estimated parameters. + *

Guessing is covariance-based, it only gives rough order of magnitude.

+ * @param problem estimation problem + * @return errors in estimated parameters + * @exception EstimationException if the covariances matrix cannot be computed + * or the number of degrees of freedom is not positive (number of measurements + * lesser or equal to number of parameters) + */ + public double[] guessParametersErrors(EstimationProblem problem) + throws EstimationException { + int m = problem.getMeasurements().length; + int p = problem.getUnboundParameters().length; + if (m <= p) { + throw new EstimationException( + "no degrees of freedom ({0} measurements, {1} parameters)", + m, p); + } + double[] errors = new double[problem.getUnboundParameters().length]; + final double c = Math.sqrt(getChiSquare(problem) / (m - p)); + double[][] covar = getCovariances(problem); + for (int i = 0; i < errors.length; ++i) { + errors[i] = Math.sqrt(covar[i][i]) * c; + } + return errors; + } + + /** + * Initialization of the common parts of the estimation. + *

This method must be called at the start + * of the {@link #estimate(EstimationProblem) estimate} + * method.

+ * @param problem estimation problem to solve + */ + protected void initializeEstimate(EstimationProblem problem) { + + // reset counters + costEvaluations = 0; + jacobianEvaluations = 0; + + // retrieve the equations and the parameters + measurements = problem.getMeasurements(); + parameters = problem.getUnboundParameters(); + + // arrays shared with the other private methods + rows = measurements.length; + cols = parameters.length; + jacobian = new double[rows * cols]; + residuals = new double[rows]; + + cost = Double.POSITIVE_INFINITY; + + } + + /** + * Solve an estimation problem. + * + *

The method should set the parameters of the problem to several + * trial values until it reaches convergence. If this method returns + * normally (i.e. without throwing an exception), then the best + * estimate of the parameters can be retrieved from the problem + * itself, through the {@link EstimationProblem#getAllParameters + * EstimationProblem.getAllParameters} method.

+ * + * @param problem estimation problem to solve + * @exception EstimationException if the problem cannot be solved + * + */ + public abstract void estimate(EstimationProblem problem) + throws EstimationException; + + /** Array of measurements. */ + protected WeightedMeasurement[] measurements; + + /** Array of parameters. */ + protected EstimatedParameter[] parameters; + + /** + * Jacobian matrix. + *

This matrix is in canonical form just after the calls to + * {@link #updateJacobian()}, but may be modified by the solver + * in the derived class (the {@link LevenbergMarquardtEstimator + * Levenberg-Marquardt estimator} does this).

+ */ + protected double[] jacobian; + + /** Number of columns of the jacobian matrix. */ + protected int cols; + + /** Number of rows of the jacobian matrix. */ + protected int rows; + + /** Residuals array. + *

This array is in canonical form just after the calls to + * {@link #updateJacobian()}, but may be modified by the solver + * in the derived class (the {@link LevenbergMarquardtEstimator + * Levenberg-Marquardt estimator} does this).

+ */ + protected double[] residuals; + + /** Cost value (square root of the sum of the residuals). */ + protected double cost; + + /** Maximal allowed number of cost evaluations. */ + private int maxCostEval; + + /** Number of cost evaluations. */ + private int costEvaluations; + + /** Number of jacobian evaluations. */ + private int jacobianEvaluations; + +} \ No newline at end of file diff --git a/src/java/org/apache/commons/math/estimation/EstimatedParameter.java b/src/java/org/apache/commons/math/estimation/EstimatedParameter.java new file mode 100644 index 000000000..b193d6701 --- /dev/null +++ b/src/java/org/apache/commons/math/estimation/EstimatedParameter.java @@ -0,0 +1,126 @@ +/* + * 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.estimation; + +import java.io.Serializable; + +/** This class represents the estimated parameters of an estimation problem. + * + *

The parameters of an estimation problem have a name, a value and + * a bound flag. The value of bound parameters is considered trusted + * and the solvers should not adjust them. On the other hand, the + * solvers should adjust the value of unbounds parameters until they + * satisfy convergence criterions specific to each solver.

+ * + * @version $Revision$ $Date$ + * @since 1.2 + * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has + * been deprecated and replaced by package org.apache.commons.math.optimization.general + * + */ +@Deprecated +public class EstimatedParameter + implements Serializable { + + /** Simple constructor. + * Build an instance from a first estimate of the parameter, + * initially considered unbound. + * @param name name of the parameter + * @param firstEstimate first estimate of the parameter + */ + public EstimatedParameter(String name, double firstEstimate) { + this.name = name; + estimate = firstEstimate; + bound = false; + } + + /** Simple constructor. + * Build an instance from a first estimate of the parameter and a + * bound flag + * @param name name of the parameter + * @param firstEstimate first estimate of the parameter + * @param bound flag, should be true if the parameter is bound + */ + public EstimatedParameter(String name, + double firstEstimate, + boolean bound) { + this.name = name; + estimate = firstEstimate; + this.bound = bound; + } + + /** Copy constructor. + * Build a copy of a parameter + * @param parameter instance to copy + */ + public EstimatedParameter(EstimatedParameter parameter) { + name = parameter.name; + estimate = parameter.estimate; + bound = parameter.bound; + } + + /** Set a new estimated value for the parameter. + * @param estimate new estimate for the parameter + */ + public void setEstimate(double estimate) { + this.estimate = estimate; + } + + /** Get the current estimate of the parameter + * @return current estimate + */ + public double getEstimate() { + return estimate; + } + + /** get the name of the parameter + * @return parameter name + */ + public String getName() { + return name; + } + + /** Set the bound flag of the parameter + * @param bound this flag should be set to true if the parameter is + * bound (i.e. if it should not be adjusted by the solver). + */ + public void setBound(boolean bound) { + this.bound = bound; + } + + /** Check if the parameter is bound + * @return true if the parameter is bound */ + public boolean isBound() { + return bound; + } + + /** Name of the parameter */ + private String name; + + /** Current value of the parameter */ + protected double estimate; + + /** Indicator for bound parameters + * (ie parameters that should not be estimated) + */ + private boolean bound; + + /** Serializable version identifier */ + private static final long serialVersionUID = -555440800213416949L; + +} diff --git a/src/java/org/apache/commons/math/estimation/EstimationException.java b/src/java/org/apache/commons/math/estimation/EstimationException.java new file mode 100644 index 000000000..f3fda25d8 --- /dev/null +++ b/src/java/org/apache/commons/math/estimation/EstimationException.java @@ -0,0 +1,48 @@ +/* + * 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.estimation; + +import org.apache.commons.math.MathException; + +/** + * This class represents exceptions thrown by the estimation solvers. + * + * @version $Revision$ $Date$ + * @since 1.2 + * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has + * been deprecated and replaced by package org.apache.commons.math.optimization.general + * + */ +@Deprecated +public class EstimationException +extends MathException { + + /** Serializable version identifier. */ + private static final long serialVersionUID = -573038581493881337L; + + /** + * Simple constructor. + * Build an exception by translating and formating a message + * @param specifier format specifier (to be translated) + * @param parts to insert in the format (no translation) + */ + public EstimationException(String specifier, Object ... parts) { + super(specifier, parts); + } + +} diff --git a/src/java/org/apache/commons/math/estimation/EstimationProblem.java b/src/java/org/apache/commons/math/estimation/EstimationProblem.java new file mode 100644 index 000000000..7ce7007cd --- /dev/null +++ b/src/java/org/apache/commons/math/estimation/EstimationProblem.java @@ -0,0 +1,67 @@ +/* + * 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.estimation; + +/** + * This interface represents an estimation problem. + * + *

This interface should be implemented by all real estimation + * problems before they can be handled by the estimators through the + * {@link Estimator#estimate Estimator.estimate} method.

+ * + *

An estimation problem, as seen by a solver is a set of + * parameters and a set of measurements. The parameters are adjusted + * during the estimation through the {@link #getUnboundParameters + * getUnboundParameters} and {@link EstimatedParameter#setEstimate + * EstimatedParameter.setEstimate} methods. The measurements both have + * a measured value which is generally fixed at construction and a + * theoretical value which depends on the model and hence varies as + * the parameters are adjusted. The purpose of the solver is to reduce + * the residual between these values, it can retrieve the measurements + * through the {@link #getMeasurements getMeasurements} method.

+ * + * @see Estimator + * @see WeightedMeasurement + * + * @version $Revision$ $Date$ + * @since 1.2 + * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has + * been deprecated and replaced by package org.apache.commons.math.optimization.general + * + */ +@Deprecated +public interface EstimationProblem { + /** + * Get the measurements of an estimation problem. + * @return measurements + */ + public WeightedMeasurement[] getMeasurements(); + + /** + * Get the unbound parameters of the problem. + * @return unbound parameters + */ + public EstimatedParameter[] getUnboundParameters(); + + /** + * Get all the parameters of the problem. + * @return parameters + */ + public EstimatedParameter[] getAllParameters(); + +} diff --git a/src/java/org/apache/commons/math/estimation/Estimator.java b/src/java/org/apache/commons/math/estimation/Estimator.java new file mode 100644 index 000000000..57c851171 --- /dev/null +++ b/src/java/org/apache/commons/math/estimation/Estimator.java @@ -0,0 +1,93 @@ +/* + * 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.estimation; + +/** + * This interface represents solvers for estimation problems. + * + *

The classes which are devoted to solve estimation problems + * should implement this interface. The problems which can be handled + * should implement the {@link EstimationProblem} interface which + * gather all the information needed by the solver.

+ * + *

The interface is composed only of the {@link #estimate estimate} + * method.

+ * + * @see EstimationProblem + * + * @version $Revision$ $Date$ + * @since 1.2 + * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has + * been deprecated and replaced by package org.apache.commons.math.optimization.general + * + */ +@Deprecated +public interface Estimator { + + /** + * Solve an estimation problem. + * + *

The method should set the parameters of the problem to several + * trial values until it reaches convergence. If this method returns + * normally (i.e. without throwing an exception), then the best + * estimate of the parameters can be retrieved from the problem + * itself, through the {@link EstimationProblem#getAllParameters + * EstimationProblem.getAllParameters} method.

+ * + * @param problem estimation problem to solve + * @exception EstimationException if the problem cannot be solved + * + */ + public void estimate(EstimationProblem problem) + throws EstimationException; + + /** + * Get the Root Mean Square value. + * Get the Root Mean Square value, i.e. the root of the arithmetic + * mean of the square of all weighted residuals. This is related to the + * criterion that is minimized by the estimator as follows: if + * c is the criterion, and n is the number of + * measurements, then the RMS is sqrt (c/n). + * @see #guessParametersErrors(EstimationProblem) + * + * @param problem estimation problem + * @return RMS value + */ + public double getRMS(EstimationProblem problem); + + /** + * Get the covariance matrix of estimated parameters. + * @param problem estimation problem + * @return covariance matrix + * @exception EstimationException if the covariance matrix + * cannot be computed (singular problem) + */ + public double[][] getCovariances(EstimationProblem problem) + throws EstimationException; + + /** + * Guess the errors in estimated parameters. + * @see #getRMS(EstimationProblem) + * @param problem estimation problem + * @return errors in estimated parameters + * @exception EstimationException if the error cannot be guessed + */ + public double[] guessParametersErrors(EstimationProblem problem) + throws EstimationException; + +} diff --git a/src/java/org/apache/commons/math/estimation/GaussNewtonEstimator.java b/src/java/org/apache/commons/math/estimation/GaussNewtonEstimator.java new file mode 100644 index 000000000..eb76e0e99 --- /dev/null +++ b/src/java/org/apache/commons/math/estimation/GaussNewtonEstimator.java @@ -0,0 +1,228 @@ +/* + * 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.estimation; + +import java.io.Serializable; + +import org.apache.commons.math.linear.InvalidMatrixException; +import org.apache.commons.math.linear.MatrixUtils; +import org.apache.commons.math.linear.RealMatrix; +import org.apache.commons.math.linear.RealVector; +import org.apache.commons.math.linear.RealVectorImpl; +import org.apache.commons.math.linear.decomposition.LUDecompositionImpl; + +/** + * This class implements a solver for estimation problems. + * + *

This class solves estimation problems using a weighted least + * squares criterion on the measurement residuals. It uses a + * Gauss-Newton algorithm.

+ * + * @version $Revision$ $Date$ + * @since 1.2 + * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has + * been deprecated and replaced by package org.apache.commons.math.optimization.general + * + */ +@Deprecated +public class GaussNewtonEstimator extends AbstractEstimator implements Serializable { + + /** Serializable version identifier */ + private static final long serialVersionUID = 5485001826076289109L; + + /** Default threshold for cost steady state detection. */ + private static final double DEFAULT_STEADY_STATE_THRESHOLD = 1.0e-6; + + /** Default threshold for cost convergence. */ + private static final double DEFAULT_CONVERGENCE = 1.0e-6; + + /** Threshold for cost steady state detection. */ + private double steadyStateThreshold; + + /** Threshold for cost convergence. */ + private double convergence; + + /** Simple constructor with default settings. + *

+ * The estimator is built with default values for all settings. + *

+ * @see #DEFAULT_STEADY_STATE_THRESHOLD + * @see #DEFAULT_CONVERGENCE + * @see AbstractEstimator#DEFAULT_MAX_COST_EVALUATIONS + */ + public GaussNewtonEstimator() { + this.steadyStateThreshold = DEFAULT_STEADY_STATE_THRESHOLD; + this.convergence = DEFAULT_CONVERGENCE; + } + + /** + * Simple constructor. + * + *

This constructor builds an estimator and stores its convergence + * characteristics.

+ * + *

An estimator is considered to have converged whenever either + * the criterion goes below a physical threshold under which + * improvements are considered useless or when the algorithm is + * unable to improve it (even if it is still high). The first + * condition that is met stops the iterations.

+ * + *

The fact an estimator has converged does not mean that the + * model accurately fits the measurements. It only means no better + * solution can be found, it does not mean this one is good. Such an + * analysis is left to the caller.

+ * + *

If neither conditions are fulfilled before a given number of + * iterations, the algorithm is considered to have failed and an + * {@link EstimationException} is thrown.

+ * + * @param maxCostEval maximal number of cost evaluations allowed + * @param convergence criterion threshold below which we do not need + * to improve the criterion anymore + * @param steadyStateThreshold steady state detection threshold, the + * problem has converged has reached a steady state if + * Math.abs(Jn - Jn-1) < + * Jn × convergence, where Jn + * and Jn-1 are the current and preceding criterion + * values (square sum of the weighted residuals of considered measurements). + */ + public GaussNewtonEstimator(final int maxCostEval, final double convergence, + final double steadyStateThreshold) { + setMaxCostEval(maxCostEval); + this.steadyStateThreshold = steadyStateThreshold; + this.convergence = convergence; + } + + /** + * Set the convergence criterion threshold. + * @param convergence criterion threshold below which we do not need + * to improve the criterion anymore + */ + public void setConvergence(final double convergence) { + this.convergence = convergence; + } + + /** + * Set the steady state detection threshold. + *

+ * The problem has converged has reached a steady state if + * Math.abs(Jn - Jn-1) < + * Jn × convergence, where Jn + * and Jn-1 are the current and preceding criterion + * values (square sum of the weighted residuals of considered measurements). + *

+ * @param steadyStateThreshold steady state detection threshold + */ + public void setSteadyStateThreshold(final double steadyStateThreshold) { + this.steadyStateThreshold = steadyStateThreshold; + } + + /** + * Solve an estimation problem using a least squares criterion. + * + *

This method set the unbound parameters of the given problem + * starting from their current values through several iterations. At + * each step, the unbound parameters are changed in order to + * minimize a weighted least square criterion based on the + * measurements of the problem.

+ * + *

The iterations are stopped either when the criterion goes + * below a physical threshold under which improvement are considered + * useless or when the algorithm is unable to improve it (even if it + * is still high). The first condition that is met stops the + * iterations. If the convergence it not reached before the maximum + * number of iterations, an {@link EstimationException} is + * thrown.

+ * + * @param problem estimation problem to solve + * @exception EstimationException if the problem cannot be solved + * + * @see EstimationProblem + * + */ + public void estimate(EstimationProblem problem) + throws EstimationException { + + initializeEstimate(problem); + + // work matrices + double[] grad = new double[parameters.length]; + RealVectorImpl bDecrement = new RealVectorImpl(parameters.length); + double[] bDecrementData = bDecrement.getDataRef(); + RealMatrix wGradGradT = MatrixUtils.createRealMatrix(parameters.length, parameters.length); + + // iterate until convergence is reached + double previous = Double.POSITIVE_INFINITY; + do { + + // build the linear problem + incrementJacobianEvaluationsCounter(); + RealVector b = new RealVectorImpl(parameters.length); + RealMatrix a = MatrixUtils.createRealMatrix(parameters.length, parameters.length); + for (int i = 0; i < measurements.length; ++i) { + if (! measurements [i].isIgnored()) { + + double weight = measurements[i].getWeight(); + double residual = measurements[i].getResidual(); + + // compute the normal equation + for (int j = 0; j < parameters.length; ++j) { + grad[j] = measurements[i].getPartial(parameters[j]); + bDecrementData[j] = weight * residual * grad[j]; + } + + // build the contribution matrix for measurement i + for (int k = 0; k < parameters.length; ++k) { + double gk = grad[k]; + for (int l = 0; l < parameters.length; ++l) { + wGradGradT.setEntry(k, l, weight * gk * grad[l]); + } + } + + // update the matrices + a = a.add(wGradGradT); + b = b.add(bDecrement); + + } + } + + try { + + // solve the linearized least squares problem + RealVector dX = new LUDecompositionImpl(a).getSolver().solve(b); + + // update the estimated parameters + for (int i = 0; i < parameters.length; ++i) { + parameters[i].setEstimate(parameters[i].getEstimate() + dX.getEntry(i)); + } + + } catch(InvalidMatrixException e) { + throw new EstimationException("unable to solve: singular problem"); + } + + + previous = cost; + updateResidualsAndCost(); + + } while ((getCostEvaluations() < 2) || + (Math.abs(previous - cost) > (cost * steadyStateThreshold) && + (Math.abs(cost) > convergence))); + + } + +} diff --git a/src/java/org/apache/commons/math/estimation/LevenbergMarquardtEstimator.java b/src/java/org/apache/commons/math/estimation/LevenbergMarquardtEstimator.java new file mode 100644 index 000000000..cf9139699 --- /dev/null +++ b/src/java/org/apache/commons/math/estimation/LevenbergMarquardtEstimator.java @@ -0,0 +1,874 @@ +/* + * 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.estimation; + +import java.io.Serializable; +import java.util.Arrays; + + +/** + * This class solves a least squares problem. + * + *

This implementation should work even for over-determined systems + * (i.e. systems having more variables than equations). Over-determined systems + * are solved by ignoring the variables which have the smallest impact according + * to their jacobian column norm. Only the rank of the matrix and some loop bounds + * are changed to implement this.

+ * + *

The resolution engine is a simple translation of the MINPACK lmder routine with minor + * changes. The changes include the over-determined resolution and the Q.R. + * decomposition which has been rewritten following the algorithm described in the + * P. Lascaux and R. Theodor book Analyse numérique matricielle + * appliquée à l'art de l'ingénieur, Masson 1986. The + * redistribution policy for MINPACK is available here, for convenience, it + * is reproduced below.

+ * + * + * + * + *
+ * Minpack Copyright Notice (1999) University of Chicago. + * All rights reserved + *
+ * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + *
    + *
  1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer.
    2. + *
  2. Redistributions in binary form must reproduce the above + * copyright notice, this list of conditions and the following + * disclaimer in the documentation and/or other materials provided + * with the distribution.
    3. + *
  3. The end-user documentation included with the redistribution, if any, + * must include the following acknowledgment: + * This product includes software developed by the University of + * Chicago, as Operator of Argonne National Laboratory. + * Alternately, this acknowledgment may appear in the software itself, + * if and wherever such third-party acknowledgments normally appear.
    4. + *
  4. WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS" + * WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE + * UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND + * THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES + * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE + * OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY + * OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR + * USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF + * THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4) + * DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION + * UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL + * BE CORRECTED.
    5. + *
  5. LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT + * HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF + * ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT, + * INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF + * ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF + * PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER + * SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT + * (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE, + * EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE + * POSSIBILITY OF SUCH LOSS OR DAMAGES.
    6. + *
    + + * @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran) + * @author Burton S. Garbow (original fortran) + * @author Kenneth E. Hillstrom (original fortran) + * @author Jorge J. More (original fortran) + + * @version $Revision$ $Date$ + * @since 1.2 + * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has + * been deprecated and replaced by package org.apache.commons.math.optimization.general + * + */ +@Deprecated +public class LevenbergMarquardtEstimator extends AbstractEstimator implements Serializable { + + /** + * Build an estimator for least squares problems. + *

    The default values for the algorithm settings are: + *

    + *

    + */ + public LevenbergMarquardtEstimator() { + + // set up the superclass with a default max cost evaluations setting + setMaxCostEval(1000); + + // default values for the tuning parameters + setInitialStepBoundFactor(100.0); + setCostRelativeTolerance(1.0e-10); + setParRelativeTolerance(1.0e-10); + setOrthoTolerance(1.0e-10); + + } + + /** + * Set the positive input variable used in determining the initial step bound. + * This bound is set to the product of initialStepBoundFactor and the euclidean norm of diag*x if nonzero, + * or else to initialStepBoundFactor itself. In most cases factor should lie + * in the interval (0.1, 100.0). 100.0 is a generally recommended value + * + * @param initialStepBoundFactor initial step bound factor + * @see #estimate + */ + public void setInitialStepBoundFactor(double initialStepBoundFactor) { + this.initialStepBoundFactor = initialStepBoundFactor; + } + + /** + * Set the desired relative error in the sum of squares. + * + * @param costRelativeTolerance desired relative error in the sum of squares + * @see #estimate + */ + public void setCostRelativeTolerance(double costRelativeTolerance) { + this.costRelativeTolerance = costRelativeTolerance; + } + + /** + * Set the desired relative error in the approximate solution parameters. + * + * @param parRelativeTolerance desired relative error + * in the approximate solution parameters + * @see #estimate + */ + public void setParRelativeTolerance(double parRelativeTolerance) { + this.parRelativeTolerance = parRelativeTolerance; + } + + /** + * Set the desired max cosine on the orthogonality. + * + * @param orthoTolerance desired max cosine on the orthogonality + * between the function vector and the columns of the jacobian + * @see #estimate + */ + public void setOrthoTolerance(double orthoTolerance) { + this.orthoTolerance = orthoTolerance; + } + + /** + * Solve an estimation problem using the Levenberg-Marquardt algorithm. + *

    The algorithm used is a modified Levenberg-Marquardt one, based + * on the MINPACK lmder + * routine. The algorithm settings must have been set up before this method + * is called with the {@link #setInitialStepBoundFactor}, + * {@link #setMaxCostEval}, {@link #setCostRelativeTolerance}, + * {@link #setParRelativeTolerance} and {@link #setOrthoTolerance} methods. + * If these methods have not been called, the default values set up by the + * {@link #LevenbergMarquardtEstimator() constructor} will be used.

    + *

    The authors of the original fortran function are:

    + * + *

    Luc Maisonobe did the Java translation.

    + * + * @param problem estimation problem to solve + * @exception EstimationException if convergence cannot be + * reached with the specified algorithm settings or if there are more variables + * than equations + * @see #setInitialStepBoundFactor + * @see #setCostRelativeTolerance + * @see #setParRelativeTolerance + * @see #setOrthoTolerance + */ + public void estimate(EstimationProblem problem) + throws EstimationException { + + initializeEstimate(problem); + + // arrays shared with the other private methods + solvedCols = Math.min(rows, cols); + diagR = new double[cols]; + jacNorm = new double[cols]; + beta = new double[cols]; + permutation = new int[cols]; + lmDir = new double[cols]; + + // local variables + double delta = 0, xNorm = 0; + double[] diag = new double[cols]; + double[] oldX = new double[cols]; + double[] oldRes = new double[rows]; + double[] work1 = new double[cols]; + double[] work2 = new double[cols]; + double[] work3 = new double[cols]; + + // evaluate the function at the starting point and calculate its norm + updateResidualsAndCost(); + + // outer loop + lmPar = 0; + boolean firstIteration = true; + while (true) { + + // compute the Q.R. decomposition of the jacobian matrix + updateJacobian(); + qrDecomposition(); + + // compute Qt.res + qTy(residuals); + + // now we don't need Q anymore, + // so let jacobian contain the R matrix with its diagonal elements + for (int k = 0; k < solvedCols; ++k) { + int pk = permutation[k]; + jacobian[k * cols + pk] = diagR[pk]; + } + + if (firstIteration) { + + // scale the variables according to the norms of the columns + // of the initial jacobian + xNorm = 0; + for (int k = 0; k < cols; ++k) { + double dk = jacNorm[k]; + if (dk == 0) { + dk = 1.0; + } + double xk = dk * parameters[k].getEstimate(); + xNorm += xk * xk; + diag[k] = dk; + } + xNorm = Math.sqrt(xNorm); + + // initialize the step bound delta + delta = (xNorm == 0) ? initialStepBoundFactor : (initialStepBoundFactor * xNorm); + + } + + // check orthogonality between function vector and jacobian columns + double maxCosine = 0; + if (cost != 0) { + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + double s = jacNorm[pj]; + if (s != 0) { + double sum = 0; + for (int i = 0, index = pj; i <= j; ++i, index += cols) { + sum += jacobian[index] * residuals[i]; + } + maxCosine = Math.max(maxCosine, Math.abs(sum) / (s * cost)); + } + } + } + if (maxCosine <= orthoTolerance) { + return; + } + + // rescale if necessary + for (int j = 0; j < cols; ++j) { + diag[j] = Math.max(diag[j], jacNorm[j]); + } + + // inner loop + for (double ratio = 0; ratio < 1.0e-4;) { + + // save the state + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + oldX[pj] = parameters[pj].getEstimate(); + } + double previousCost = cost; + double[] tmpVec = residuals; + residuals = oldRes; + oldRes = tmpVec; + + // determine the Levenberg-Marquardt parameter + determineLMParameter(oldRes, delta, diag, work1, work2, work3); + + // compute the new point and the norm of the evolution direction + double lmNorm = 0; + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + lmDir[pj] = -lmDir[pj]; + parameters[pj].setEstimate(oldX[pj] + lmDir[pj]); + double s = diag[pj] * lmDir[pj]; + lmNorm += s * s; + } + lmNorm = Math.sqrt(lmNorm); + + // on the first iteration, adjust the initial step bound. + if (firstIteration) { + delta = Math.min(delta, lmNorm); + } + + // evaluate the function at x + p and calculate its norm + updateResidualsAndCost(); + + // compute the scaled actual reduction + double actRed = -1.0; + if (0.1 * cost < previousCost) { + double r = cost / previousCost; + actRed = 1.0 - r * r; + } + + // compute the scaled predicted reduction + // and the scaled directional derivative + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + double dirJ = lmDir[pj]; + work1[j] = 0; + for (int i = 0, index = pj; i <= j; ++i, index += cols) { + work1[i] += jacobian[index] * dirJ; + } + } + double coeff1 = 0; + for (int j = 0; j < solvedCols; ++j) { + coeff1 += work1[j] * work1[j]; + } + double pc2 = previousCost * previousCost; + coeff1 = coeff1 / pc2; + double coeff2 = lmPar * lmNorm * lmNorm / pc2; + double preRed = coeff1 + 2 * coeff2; + double dirDer = -(coeff1 + coeff2); + + // ratio of the actual to the predicted reduction + ratio = (preRed == 0) ? 0 : (actRed / preRed); + + // update the step bound + if (ratio <= 0.25) { + double tmp = + (actRed < 0) ? (0.5 * dirDer / (dirDer + 0.5 * actRed)) : 0.5; + if ((0.1 * cost >= previousCost) || (tmp < 0.1)) { + tmp = 0.1; + } + delta = tmp * Math.min(delta, 10.0 * lmNorm); + lmPar /= tmp; + } else if ((lmPar == 0) || (ratio >= 0.75)) { + delta = 2 * lmNorm; + lmPar *= 0.5; + } + + // test for successful iteration. + if (ratio >= 1.0e-4) { + // successful iteration, update the norm + firstIteration = false; + xNorm = 0; + for (int k = 0; k < cols; ++k) { + double xK = diag[k] * parameters[k].getEstimate(); + xNorm += xK * xK; + } + xNorm = Math.sqrt(xNorm); + } else { + // failed iteration, reset the previous values + cost = previousCost; + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + parameters[pj].setEstimate(oldX[pj]); + } + tmpVec = residuals; + residuals = oldRes; + oldRes = tmpVec; + } + + // tests for convergence. + if (((Math.abs(actRed) <= costRelativeTolerance) && + (preRed <= costRelativeTolerance) && + (ratio <= 2.0)) || + (delta <= parRelativeTolerance * xNorm)) { + return; + } + + // tests for termination and stringent tolerances + // (2.2204e-16 is the machine epsilon for IEEE754) + if ((Math.abs(actRed) <= 2.2204e-16) && (preRed <= 2.2204e-16) && (ratio <= 2.0)) { + throw new EstimationException("cost relative tolerance is too small ({0})," + + " no further reduction in the" + + " sum of squares is possible", + costRelativeTolerance); + } else if (delta <= 2.2204e-16 * xNorm) { + throw new EstimationException("parameters relative tolerance is too small" + + " ({0}), no further improvement in" + + " the approximate solution is possible", + parRelativeTolerance); + } else if (maxCosine <= 2.2204e-16) { + throw new EstimationException("orthogonality tolerance is too small ({0})," + + " solution is orthogonal to the jacobian", + orthoTolerance); + } + + } + + } + + } + + /** + * Determine the Levenberg-Marquardt parameter. + *

    This implementation is a translation in Java of the MINPACK + * lmpar + * routine.

    + *

    This method sets the lmPar and lmDir attributes.

    + *

    The authors of the original fortran function are:

    + * + *

    Luc Maisonobe did the Java translation.

    + * + * @param qy array containing qTy + * @param delta upper bound on the euclidean norm of diagR * lmDir + * @param diag diagonal matrix + * @param work1 work array + * @param work2 work array + * @param work3 work array + */ + private void determineLMParameter(double[] qy, double delta, double[] diag, + double[] work1, double[] work2, double[] work3) { + + // compute and store in x the gauss-newton direction, if the + // jacobian is rank-deficient, obtain a least squares solution + for (int j = 0; j < rank; ++j) { + lmDir[permutation[j]] = qy[j]; + } + for (int j = rank; j < cols; ++j) { + lmDir[permutation[j]] = 0; + } + for (int k = rank - 1; k >= 0; --k) { + int pk = permutation[k]; + double ypk = lmDir[pk] / diagR[pk]; + for (int i = 0, index = pk; i < k; ++i, index += cols) { + lmDir[permutation[i]] -= ypk * jacobian[index]; + } + lmDir[pk] = ypk; + } + + // evaluate the function at the origin, and test + // for acceptance of the Gauss-Newton direction + double dxNorm = 0; + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + double s = diag[pj] * lmDir[pj]; + work1[pj] = s; + dxNorm += s * s; + } + dxNorm = Math.sqrt(dxNorm); + double fp = dxNorm - delta; + if (fp <= 0.1 * delta) { + lmPar = 0; + return; + } + + // if the jacobian is not rank deficient, the Newton step provides + // a lower bound, parl, for the zero of the function, + // otherwise set this bound to zero + double sum2, parl = 0; + if (rank == solvedCols) { + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + work1[pj] *= diag[pj] / dxNorm; + } + sum2 = 0; + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + double sum = 0; + for (int i = 0, index = pj; i < j; ++i, index += cols) { + sum += jacobian[index] * work1[permutation[i]]; + } + double s = (work1[pj] - sum) / diagR[pj]; + work1[pj] = s; + sum2 += s * s; + } + parl = fp / (delta * sum2); + } + + // calculate an upper bound, paru, for the zero of the function + sum2 = 0; + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + double sum = 0; + for (int i = 0, index = pj; i <= j; ++i, index += cols) { + sum += jacobian[index] * qy[i]; + } + sum /= diag[pj]; + sum2 += sum * sum; + } + double gNorm = Math.sqrt(sum2); + double paru = gNorm / delta; + if (paru == 0) { + // 2.2251e-308 is the smallest positive real for IEE754 + paru = 2.2251e-308 / Math.min(delta, 0.1); + } + + // if the input par lies outside of the interval (parl,paru), + // set par to the closer endpoint + lmPar = Math.min(paru, Math.max(lmPar, parl)); + if (lmPar == 0) { + lmPar = gNorm / dxNorm; + } + + for (int countdown = 10; countdown >= 0; --countdown) { + + // evaluate the function at the current value of lmPar + if (lmPar == 0) { + lmPar = Math.max(2.2251e-308, 0.001 * paru); + } + double sPar = Math.sqrt(lmPar); + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + work1[pj] = sPar * diag[pj]; + } + determineLMDirection(qy, work1, work2, work3); + + dxNorm = 0; + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + double s = diag[pj] * lmDir[pj]; + work3[pj] = s; + dxNorm += s * s; + } + dxNorm = Math.sqrt(dxNorm); + double previousFP = fp; + fp = dxNorm - delta; + + // if the function is small enough, accept the current value + // of lmPar, also test for the exceptional cases where parl is zero + if ((Math.abs(fp) <= 0.1 * delta) || + ((parl == 0) && (fp <= previousFP) && (previousFP < 0))) { + return; + } + + // compute the Newton correction + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + work1[pj] = work3[pj] * diag[pj] / dxNorm; + } + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + work1[pj] /= work2[j]; + double tmp = work1[pj]; + for (int i = j + 1; i < solvedCols; ++i) { + work1[permutation[i]] -= jacobian[i * cols + pj] * tmp; + } + } + sum2 = 0; + for (int j = 0; j < solvedCols; ++j) { + double s = work1[permutation[j]]; + sum2 += s * s; + } + double correction = fp / (delta * sum2); + + // depending on the sign of the function, update parl or paru. + if (fp > 0) { + parl = Math.max(parl, lmPar); + } else if (fp < 0) { + paru = Math.min(paru, lmPar); + } + + // compute an improved estimate for lmPar + lmPar = Math.max(parl, lmPar + correction); + + } + } + + /** + * Solve a*x = b and d*x = 0 in the least squares sense. + *

    This implementation is a translation in Java of the MINPACK + * qrsolv + * routine.

    + *

    This method sets the lmDir and lmDiag attributes.

    + *

    The authors of the original fortran function are:

    + * + *

    Luc Maisonobe did the Java translation.

    + * + * @param qy array containing qTy + * @param diag diagonal matrix + * @param lmDiag diagonal elements associated with lmDir + * @param work work array + */ + private void determineLMDirection(double[] qy, double[] diag, + double[] lmDiag, double[] work) { + + // copy R and Qty to preserve input and initialize s + // in particular, save the diagonal elements of R in lmDir + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + for (int i = j + 1; i < solvedCols; ++i) { + jacobian[i * cols + pj] = jacobian[j * cols + permutation[i]]; + } + lmDir[j] = diagR[pj]; + work[j] = qy[j]; + } + + // eliminate the diagonal matrix d using a Givens rotation + for (int j = 0; j < solvedCols; ++j) { + + // prepare the row of d to be eliminated, locating the + // diagonal element using p from the Q.R. factorization + int pj = permutation[j]; + double dpj = diag[pj]; + if (dpj != 0) { + Arrays.fill(lmDiag, j + 1, lmDiag.length, 0); + } + lmDiag[j] = dpj; + + // the transformations to eliminate the row of d + // modify only a single element of Qty + // beyond the first n, which is initially zero. + double qtbpj = 0; + for (int k = j; k < solvedCols; ++k) { + int pk = permutation[k]; + + // determine a Givens rotation which eliminates the + // appropriate element in the current row of d + if (lmDiag[k] != 0) { + + double sin, cos; + double rkk = jacobian[k * cols + pk]; + if (Math.abs(rkk) < Math.abs(lmDiag[k])) { + double cotan = rkk / lmDiag[k]; + sin = 1.0 / Math.sqrt(1.0 + cotan * cotan); + cos = sin * cotan; + } else { + double tan = lmDiag[k] / rkk; + cos = 1.0 / Math.sqrt(1.0 + tan * tan); + sin = cos * tan; + } + + // compute the modified diagonal element of R and + // the modified element of (Qty,0) + jacobian[k * cols + pk] = cos * rkk + sin * lmDiag[k]; + double temp = cos * work[k] + sin * qtbpj; + qtbpj = -sin * work[k] + cos * qtbpj; + work[k] = temp; + + // accumulate the tranformation in the row of s + for (int i = k + 1; i < solvedCols; ++i) { + double rik = jacobian[i * cols + pk]; + temp = cos * rik + sin * lmDiag[i]; + lmDiag[i] = -sin * rik + cos * lmDiag[i]; + jacobian[i * cols + pk] = temp; + } + + } + } + + // store the diagonal element of s and restore + // the corresponding diagonal element of R + int index = j * cols + permutation[j]; + lmDiag[j] = jacobian[index]; + jacobian[index] = lmDir[j]; + + } + + // solve the triangular system for z, if the system is + // singular, then obtain a least squares solution + int nSing = solvedCols; + for (int j = 0; j < solvedCols; ++j) { + if ((lmDiag[j] == 0) && (nSing == solvedCols)) { + nSing = j; + } + if (nSing < solvedCols) { + work[j] = 0; + } + } + if (nSing > 0) { + for (int j = nSing - 1; j >= 0; --j) { + int pj = permutation[j]; + double sum = 0; + for (int i = j + 1; i < nSing; ++i) { + sum += jacobian[i * cols + pj] * work[i]; + } + work[j] = (work[j] - sum) / lmDiag[j]; + } + } + + // permute the components of z back to components of lmDir + for (int j = 0; j < lmDir.length; ++j) { + lmDir[permutation[j]] = work[j]; + } + + } + + /** + * Decompose a matrix A as A.P = Q.R using Householder transforms. + *

    As suggested in the P. Lascaux and R. Theodor book + * Analyse numérique matricielle appliquée à + * l'art de l'ingénieur (Masson, 1986), instead of representing + * the Householder transforms with uk unit vectors such that: + * 

    +   * Hk = I - 2uk.ukt
+   *
    + * we use k non-unit vectors such that: + * 
    +   * Hk = I - betakvk.vkt
+   *
    + * where vk = ak - alphak ek. + * The betak coefficients are provided upon exit as recomputing + * them from the vk vectors would be costly.

    + *

    This decomposition handles rank deficient cases since the tranformations + * are performed in non-increasing columns norms order thanks to columns + * pivoting. The diagonal elements of the R matrix are therefore also in + * non-increasing absolute values order.

    + * @exception EstimationException if the decomposition cannot be performed + */ + private void qrDecomposition() throws EstimationException { + + // initializations + for (int k = 0; k < cols; ++k) { + permutation[k] = k; + double norm2 = 0; + for (int index = k; index < jacobian.length; index += cols) { + double akk = jacobian[index]; + norm2 += akk * akk; + } + jacNorm[k] = Math.sqrt(norm2); + } + + // transform the matrix column after column + for (int k = 0; k < cols; ++k) { + + // select the column with the greatest norm on active components + int nextColumn = -1; + double ak2 = Double.NEGATIVE_INFINITY; + for (int i = k; i < cols; ++i) { + double norm2 = 0; + int iDiag = k * cols + permutation[i]; + for (int index = iDiag; index < jacobian.length; index += cols) { + double aki = jacobian[index]; + norm2 += aki * aki; + } + if (Double.isInfinite(norm2) || Double.isNaN(norm2)) { + throw new EstimationException( + "unable to perform Q.R decomposition on the {0}x{1} jacobian matrix", + rows, cols); + } + if (norm2 > ak2) { + nextColumn = i; + ak2 = norm2; + } + } + if (ak2 == 0) { + rank = k; + return; + } + int pk = permutation[nextColumn]; + permutation[nextColumn] = permutation[k]; + permutation[k] = pk; + + // choose alpha such that Hk.u = alpha ek + int kDiag = k * cols + pk; + double akk = jacobian[kDiag]; + double alpha = (akk > 0) ? -Math.sqrt(ak2) : Math.sqrt(ak2); + double betak = 1.0 / (ak2 - akk * alpha); + beta[pk] = betak; + + // transform the current column + diagR[pk] = alpha; + jacobian[kDiag] -= alpha; + + // transform the remaining columns + for (int dk = cols - 1 - k; dk > 0; --dk) { + int dkp = permutation[k + dk] - pk; + double gamma = 0; + for (int index = kDiag; index < jacobian.length; index += cols) { + gamma += jacobian[index] * jacobian[index + dkp]; + } + gamma *= betak; + for (int index = kDiag; index < jacobian.length; index += cols) { + jacobian[index + dkp] -= gamma * jacobian[index]; + } + } + + } + + rank = solvedCols; + + } + + /** + * Compute the product Qt.y for some Q.R. decomposition. + * + * @param y vector to multiply (will be overwritten with the result) + */ + private void qTy(double[] y) { + for (int k = 0; k < cols; ++k) { + int pk = permutation[k]; + int kDiag = k * cols + pk; + double gamma = 0; + for (int i = k, index = kDiag; i < rows; ++i, index += cols) { + gamma += jacobian[index] * y[i]; + } + gamma *= beta[pk]; + for (int i = k, index = kDiag; i < rows; ++i, index += cols) { + y[i] -= gamma * jacobian[index]; + } + } + } + + /** Number of solved variables. */ + private int solvedCols; + + /** Diagonal elements of the R matrix in the Q.R. decomposition. */ + private double[] diagR; + + /** Norms of the columns of the jacobian matrix. */ + private double[] jacNorm; + + /** Coefficients of the Householder transforms vectors. */ + private double[] beta; + + /** Columns permutation array. */ + private int[] permutation; + + /** Rank of the jacobian matrix. */ + private int rank; + + /** Levenberg-Marquardt parameter. */ + private double lmPar; + + /** Parameters evolution direction associated with lmPar. */ + private double[] lmDir; + + /** Positive input variable used in determining the initial step bound. */ + private double initialStepBoundFactor; + + /** Desired relative error in the sum of squares. */ + private double costRelativeTolerance; + + /** Desired relative error in the approximate solution parameters. */ + private double parRelativeTolerance; + + /** Desired max cosine on the orthogonality between the function vector + * and the columns of the jacobian. */ + private double orthoTolerance; + + /** Serializable version identifier */ + private static final long serialVersionUID = -5705952631533171019L; + +} diff --git a/src/java/org/apache/commons/math/estimation/SimpleEstimationProblem.java b/src/java/org/apache/commons/math/estimation/SimpleEstimationProblem.java new file mode 100644 index 000000000..b5f0a2bdd --- /dev/null +++ b/src/java/org/apache/commons/math/estimation/SimpleEstimationProblem.java @@ -0,0 +1,111 @@ +/* + * 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.estimation; + +import java.util.ArrayList; +import java.util.List; + +/** + * Simple implementation of the {@link EstimationProblem + * EstimationProblem} interface for boilerplate data handling. + *

    This class only handles parameters and measurements + * storage and unbound parameters filtering. It does not compute + * anything by itself. It should either be used with measurements + * implementation that are smart enough to know about the + * various parameters in order to compute the partial derivatives + * appropriately. Since the problem-specific logic is mainly related to + * the various measurements models, the simplest way to use this class + * is by extending it and using one internal class extending + * {@link WeightedMeasurement WeightedMeasurement} for each measurement + * type. The instances of the internal classes would have access to the + * various parameters and their current estimate.

    + + * @version $Revision$ $Date$ + * @since 1.2 + * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has + * been deprecated and replaced by package org.apache.commons.math.optimization.general + + */ +@Deprecated +public class SimpleEstimationProblem implements EstimationProblem { + + /** + * Build an empty instance without parameters nor measurements. + */ + public SimpleEstimationProblem() { + parameters = new ArrayList(); + measurements = new ArrayList(); + } + + /** + * Get all the parameters of the problem. + * @return parameters + */ + public EstimatedParameter[] getAllParameters() { + return (EstimatedParameter[]) parameters.toArray(new EstimatedParameter[parameters.size()]); + } + + /** + * Get the unbound parameters of the problem. + * @return unbound parameters + */ + public EstimatedParameter[] getUnboundParameters() { + + // filter the unbound parameters + List unbound = new ArrayList(parameters.size()); + for (EstimatedParameter p : parameters) { + if (! p.isBound()) { + unbound.add(p); + } + } + + // convert to an array + return (EstimatedParameter[]) unbound.toArray(new EstimatedParameter[unbound.size()]); + + } + + /** + * Get the measurements of an estimation problem. + * @return measurements + */ + public WeightedMeasurement[] getMeasurements() { + return (WeightedMeasurement[]) measurements.toArray(new WeightedMeasurement[measurements.size()]); + } + + /** Add a parameter to the problem. + * @param p parameter to add + */ + protected void addParameter(EstimatedParameter p) { + parameters.add(p); + } + + /** + * Add a new measurement to the set. + * @param m measurement to add + */ + protected void addMeasurement(WeightedMeasurement m) { + measurements.add(m); + } + + /** Estimated parameters. */ + private final List parameters; + + /** Measurements. */ + private final List measurements; + +} diff --git a/src/java/org/apache/commons/math/estimation/WeightedMeasurement.java b/src/java/org/apache/commons/math/estimation/WeightedMeasurement.java new file mode 100644 index 000000000..de191806a --- /dev/null +++ b/src/java/org/apache/commons/math/estimation/WeightedMeasurement.java @@ -0,0 +1,172 @@ +/* + * 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.estimation; + +import java.io.Serializable; + +/** + * This class represents measurements in estimation problems. + * + *

    This abstract class implements all the methods needed to handle + * measurements in a general way. It defines neither the {@link + * #getTheoreticalValue getTheoreticalValue} nor the {@link + * #getPartial getPartial} methods, which should be defined by + * sub-classes according to the specific problem.

    + * + *

    The {@link #getTheoreticalValue getTheoreticalValue} and {@link + * #getPartial getPartial} methods must always use the current + * estimate of the parameters set by the solver in the problem. These + * parameters can be retrieved through the {@link + * EstimationProblem#getAllParameters + * EstimationProblem.getAllParameters} method if the measurements are + * independent of the problem, or directly if they are implemented as + * inner classes of the problem.

    + * + *

    The instances for which the ignored flag is set + * through the {@link #setIgnored setIgnored} method are ignored by the + * solvers. This can be used to reject wrong measurements at some + * steps of the estimation.

    + * + * @see EstimationProblem + * + * @version $Revision$ $Date$ + * @since 1.2 + * @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has + * been deprecated and replaced by package org.apache.commons.math.optimization.general + */ + +@Deprecated +public abstract class WeightedMeasurement implements Serializable { + + /** Serializable version identifier. */ + private static final long serialVersionUID = 4360046376796901941L; + + /** + * Simple constructor. + * Build a measurement with the given parameters, and set its ignore + * flag to false. + * @param weight weight of the measurement in the least squares problem + * (two common choices are either to use 1.0 for all measurements, or to + * use a value proportional to the inverse of the variance of the measurement + * type) + * + * @param measuredValue measured value + */ + public WeightedMeasurement(double weight, double measuredValue) { + this.weight = weight; + this.measuredValue = measuredValue; + ignored = false; + } + + /** Simple constructor. + * + * Build a measurement with the given parameters + * + * @param weight weight of the measurement in the least squares problem + * @param measuredValue measured value + * @param ignored true if the measurement should be ignored + */ + public WeightedMeasurement(double weight, double measuredValue, + boolean ignored) { + this.weight = weight; + this.measuredValue = measuredValue; + this.ignored = ignored; + } + + /** + * Get the weight of the measurement in the least squares problem + * + * @return weight + */ + public double getWeight() { + return weight; + } + + /** + * Get the measured value + * + * @return measured value + */ + public double getMeasuredValue() { + return measuredValue; + } + + /** + * Get the residual for this measurement + * The residual is the measured value minus the theoretical value. + * + * @return residual + */ + public double getResidual() { + return measuredValue - getTheoreticalValue(); + } + + /** + * Get the theoretical value expected for this measurement + *

    The theoretical value is the value expected for this measurement + * if the model and its parameter were all perfectly known.

    + *

    The value must be computed using the current estimate of the parameters + * set by the solver in the problem.

    + * + * @return theoretical value + */ + public abstract double getTheoreticalValue(); + + /** + * Get the partial derivative of the {@link #getTheoreticalValue + * theoretical value} according to the parameter. + *

    The value must be computed using the current estimate of the parameters + * set by the solver in the problem.

    + * + * @param parameter parameter against which the partial derivative + * should be computed + * @return partial derivative of the {@link #getTheoreticalValue + * theoretical value} + */ + public abstract double getPartial(EstimatedParameter parameter); + + /** + * Set the ignore flag to the specified value + * Setting the ignore flag to true allow to reject wrong + * measurements, which sometimes can be detected only rather late. + * + * @param ignored value for the ignore flag + */ + public void setIgnored(boolean ignored) { + this.ignored = ignored; + } + + /** + * Check if this measurement should be ignored + * + * @return true if the measurement should be ignored + */ + public boolean isIgnored() { + return ignored; + } + + /** Measurement weight. */ + private final double weight; + + /** Value of the measurements. */ + private final double measuredValue; + + /** Ignore measurement indicator. */ + private boolean ignored; + +} diff --git a/src/java/org/apache/commons/math/estimation/package.html b/src/java/org/apache/commons/math/estimation/package.html new file mode 100644 index 000000000..9181d3690 --- /dev/null +++ b/src/java/org/apache/commons/math/estimation/package.html @@ -0,0 +1,25 @@ + + + + +This package provided classes to solve estimation problems, it is deprecated since 2.0. + +

    This package has been deprecated as of 2.0. It is replaced by the optimization.general package.

    + + + diff --git a/src/test/org/apache/commons/math/estimation/EstimatedParameterTest.java b/src/test/org/apache/commons/math/estimation/EstimatedParameterTest.java new file mode 100644 index 000000000..6331a2156 --- /dev/null +++ b/src/test/org/apache/commons/math/estimation/EstimatedParameterTest.java @@ -0,0 +1,77 @@ +/* + * 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.estimation; + +import org.apache.commons.math.estimation.EstimatedParameter; + +import junit.framework.*; + +@Deprecated +public class EstimatedParameterTest + extends TestCase { + + public EstimatedParameterTest(String name) { + super(name); + } + + public void testConstruction() { + + EstimatedParameter p1 = new EstimatedParameter("p1", 1.0); + assertTrue(p1.getName().equals("p1")); + checkValue(p1.getEstimate(), 1.0); + assertTrue(! p1.isBound()); + + EstimatedParameter p2 = new EstimatedParameter("p2", 2.0, true); + assertTrue(p2.getName().equals("p2")); + checkValue(p2.getEstimate(), 2.0); + assertTrue(p2.isBound()); + + } + + public void testBound() { + + EstimatedParameter p = new EstimatedParameter("p", 0.0); + assertTrue(! p.isBound()); + p.setBound(true); + assertTrue(p.isBound()); + p.setBound(false); + assertTrue(! p.isBound()); + + } + + public void testEstimate() { + + EstimatedParameter p = new EstimatedParameter("p", 0.0); + checkValue(p.getEstimate(), 0.0); + + for (double e = 0.0; e < 10.0; e += 0.5) { + p.setEstimate(e); + checkValue(p.getEstimate(), e); + } + + } + + public static Test suite() { + return new TestSuite(EstimatedParameterTest.class); + } + + private void checkValue(double value, double expected) { + assertTrue(Math.abs(value - expected) < 1.0e-10); + } + +} diff --git a/src/test/org/apache/commons/math/estimation/GaussNewtonEstimatorTest.java b/src/test/org/apache/commons/math/estimation/GaussNewtonEstimatorTest.java new file mode 100644 index 000000000..e150ce5c7 --- /dev/null +++ b/src/test/org/apache/commons/math/estimation/GaussNewtonEstimatorTest.java @@ -0,0 +1,732 @@ +/* + * 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.estimation; + +import java.util.ArrayList; +import java.util.HashSet; + +import junit.framework.Test; +import junit.framework.TestCase; +import junit.framework.TestSuite; + +/** + *

    Some of the unit tests are re-implementations of the MINPACK file17 and file22 test files. + * The redistribution policy for MINPACK is available here, for + * convenience, it is reproduced below.

    + + * + * + * + *
    + * Minpack Copyright Notice (1999) University of Chicago. + * All rights reserved + *
    + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + *
      + *
    1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer.
      2. + *
    2. Redistributions in binary form must reproduce the above + * copyright notice, this list of conditions and the following + * disclaimer in the documentation and/or other materials provided + * with the distribution.
      3. + *
    3. The end-user documentation included with the redistribution, if any, + * must include the following acknowledgment: + * This product includes software developed by the University of + * Chicago, as Operator of Argonne National Laboratory. + * Alternately, this acknowledgment may appear in the software itself, + * if and wherever such third-party acknowledgments normally appear.
      4. + *
    4. WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS" + * WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE + * UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND + * THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES + * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE + * OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY + * OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR + * USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF + * THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4) + * DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION + * UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL + * BE CORRECTED.
      5. + *
    5. LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT + * HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF + * ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT, + * INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF + * ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF + * PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER + * SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT + * (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE, + * EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE + * POSSIBILITY OF SUCH LOSS OR DAMAGES.
      6. + *
      + + * @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran minpack tests) + * @author Burton S. Garbow (original fortran minpack tests) + * @author Kenneth E. Hillstrom (original fortran minpack tests) + * @author Jorge J. More (original fortran minpack tests) + * @author Luc Maisonobe (non-minpack tests and minpack tests Java translation) + */ +@Deprecated +public class GaussNewtonEstimatorTest + extends TestCase { + + public GaussNewtonEstimatorTest(String name) { + super(name); + } + + public void testTrivial() throws EstimationException { + LinearProblem problem = + new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] {2}, + new EstimatedParameter[] { + new EstimatedParameter("p0", 0) + }, 3.0) + }); + GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); + estimator.estimate(problem); + assertEquals(0, estimator.getRMS(problem), 1.0e-10); + assertEquals(1.5, + problem.getUnboundParameters()[0].getEstimate(), + 1.0e-10); + } + + public void testQRColumnsPermutation() throws EstimationException { + + EstimatedParameter[] x = { + new EstimatedParameter("p0", 0), new EstimatedParameter("p1", 0) + }; + LinearProblem problem = new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] { 1.0, -1.0 }, + new EstimatedParameter[] { x[0], x[1] }, + 4.0), + new LinearMeasurement(new double[] { 2.0 }, + new EstimatedParameter[] { x[1] }, + 6.0), + new LinearMeasurement(new double[] { 1.0, -2.0 }, + new EstimatedParameter[] { x[0], x[1] }, + 1.0) + }); + + GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); + estimator.estimate(problem); + assertEquals(0, estimator.getRMS(problem), 1.0e-10); + assertEquals(7.0, x[0].getEstimate(), 1.0e-10); + assertEquals(3.0, x[1].getEstimate(), 1.0e-10); + + } + + public void testNoDependency() throws EstimationException { + EstimatedParameter[] p = new EstimatedParameter[] { + new EstimatedParameter("p0", 0), + new EstimatedParameter("p1", 0), + new EstimatedParameter("p2", 0), + new EstimatedParameter("p3", 0), + new EstimatedParameter("p4", 0), + new EstimatedParameter("p5", 0) + }; + LinearProblem problem = new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[0] }, 0.0), + new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[1] }, 1.1), + new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[2] }, 2.2), + new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[3] }, 3.3), + new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[4] }, 4.4), + new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[5] }, 5.5) + }); + GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); + estimator.estimate(problem); + assertEquals(0, estimator.getRMS(problem), 1.0e-10); + for (int i = 0; i < p.length; ++i) { + assertEquals(0.55 * i, p[i].getEstimate(), 1.0e-10); + } +} + + public void testOneSet() throws EstimationException { + + EstimatedParameter[] p = { + new EstimatedParameter("p0", 0), + new EstimatedParameter("p1", 0), + new EstimatedParameter("p2", 0) + }; + LinearProblem problem = new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] { 1.0 }, + new EstimatedParameter[] { p[0] }, + 1.0), + new LinearMeasurement(new double[] { -1.0, 1.0 }, + new EstimatedParameter[] { p[0], p[1] }, + 1.0), + new LinearMeasurement(new double[] { -1.0, 1.0 }, + new EstimatedParameter[] { p[1], p[2] }, + 1.0) + }); + + GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); + estimator.estimate(problem); + assertEquals(0, estimator.getRMS(problem), 1.0e-10); + assertEquals(1.0, p[0].getEstimate(), 1.0e-10); + assertEquals(2.0, p[1].getEstimate(), 1.0e-10); + assertEquals(3.0, p[2].getEstimate(), 1.0e-10); + + } + + public void testTwoSets() throws EstimationException { + EstimatedParameter[] p = { + new EstimatedParameter("p0", 0), + new EstimatedParameter("p1", 1), + new EstimatedParameter("p2", 2), + new EstimatedParameter("p3", 3), + new EstimatedParameter("p4", 4), + new EstimatedParameter("p5", 5) + }; + + double epsilon = 1.0e-7; + LinearProblem problem = new LinearProblem(new LinearMeasurement[] { + + // 4 elements sub-problem + new LinearMeasurement(new double[] { 2.0, 1.0, 4.0 }, + new EstimatedParameter[] { p[0], p[1], p[3] }, + 2.0), + new LinearMeasurement(new double[] { -4.0, -2.0, 3.0, -7.0 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + -9.0), + new LinearMeasurement(new double[] { 4.0, 1.0, -2.0, 8.0 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + 2.0), + new LinearMeasurement(new double[] { -3.0, -12.0, -1.0 }, + new EstimatedParameter[] { p[1], p[2], p[3] }, + 2.0), + + // 2 elements sub-problem + new LinearMeasurement(new double[] { epsilon, 1.0 }, + new EstimatedParameter[] { p[4], p[5] }, + 1.0 + epsilon * epsilon), + new LinearMeasurement(new double[] { 1.0, 1.0 }, + new EstimatedParameter[] { p[4], p[5] }, + 2.0) + + }); + + GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); + estimator.estimate(problem); + assertEquals(0, estimator.getRMS(problem), 1.0e-10); + assertEquals( 3.0, p[0].getEstimate(), 1.0e-10); + assertEquals( 4.0, p[1].getEstimate(), 1.0e-10); + assertEquals(-1.0, p[2].getEstimate(), 1.0e-10); + assertEquals(-2.0, p[3].getEstimate(), 1.0e-10); + assertEquals( 1.0 + epsilon, p[4].getEstimate(), 1.0e-10); + assertEquals( 1.0 - epsilon, p[5].getEstimate(), 1.0e-10); + + } + + public void testNonInversible() throws EstimationException { + + EstimatedParameter[] p = { + new EstimatedParameter("p0", 0), + new EstimatedParameter("p1", 0), + new EstimatedParameter("p2", 0) + }; + LinearMeasurement[] m = new LinearMeasurement[] { + new LinearMeasurement(new double[] { 1.0, 2.0, -3.0 }, + new EstimatedParameter[] { p[0], p[1], p[2] }, + 1.0), + new LinearMeasurement(new double[] { 2.0, 1.0, 3.0 }, + new EstimatedParameter[] { p[0], p[1], p[2] }, + 1.0), + new LinearMeasurement(new double[] { -3.0, -9.0 }, + new EstimatedParameter[] { p[0], p[2] }, + 1.0) + }; + LinearProblem problem = new LinearProblem(m); + + GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); + try { + estimator.estimate(problem); + fail("an exception should have been caught"); + } catch (EstimationException ee) { + // expected behavior + } catch (Exception e) { + fail("wrong exception type caught"); + } + } + + public void testIllConditioned() throws EstimationException { + EstimatedParameter[] p = { + new EstimatedParameter("p0", 0), + new EstimatedParameter("p1", 1), + new EstimatedParameter("p2", 2), + new EstimatedParameter("p3", 3) + }; + + LinearProblem problem1 = new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] { 10.0, 7.0, 8.0, 7.0 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + 32.0), + new LinearMeasurement(new double[] { 7.0, 5.0, 6.0, 5.0 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + 23.0), + new LinearMeasurement(new double[] { 8.0, 6.0, 10.0, 9.0 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + 33.0), + new LinearMeasurement(new double[] { 7.0, 5.0, 9.0, 10.0 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + 31.0) + }); + GaussNewtonEstimator estimator1 = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); + estimator1.estimate(problem1); + assertEquals(0, estimator1.getRMS(problem1), 1.0e-10); + assertEquals(1.0, p[0].getEstimate(), 1.0e-10); + assertEquals(1.0, p[1].getEstimate(), 1.0e-10); + assertEquals(1.0, p[2].getEstimate(), 1.0e-10); + assertEquals(1.0, p[3].getEstimate(), 1.0e-10); + + LinearProblem problem2 = new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] { 10.0, 7.0, 8.1, 7.2 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + 32.0), + new LinearMeasurement(new double[] { 7.08, 5.04, 6.0, 5.0 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + 23.0), + new LinearMeasurement(new double[] { 8.0, 5.98, 9.89, 9.0 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + 33.0), + new LinearMeasurement(new double[] { 6.99, 4.99, 9.0, 9.98 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + 31.0) + }); + GaussNewtonEstimator estimator2 = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); + estimator2.estimate(problem2); + assertEquals(0, estimator2.getRMS(problem2), 1.0e-10); + assertEquals(-81.0, p[0].getEstimate(), 1.0e-8); + assertEquals(137.0, p[1].getEstimate(), 1.0e-8); + assertEquals(-34.0, p[2].getEstimate(), 1.0e-8); + assertEquals( 22.0, p[3].getEstimate(), 1.0e-8); + + } + + public void testMoreEstimatedParametersSimple() throws EstimationException { + + EstimatedParameter[] p = { + new EstimatedParameter("p0", 7), + new EstimatedParameter("p1", 6), + new EstimatedParameter("p2", 5), + new EstimatedParameter("p3", 4) + }; + LinearProblem problem = new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] { 3.0, 2.0 }, + new EstimatedParameter[] { p[0], p[1] }, + 7.0), + new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 }, + new EstimatedParameter[] { p[1], p[2], p[3] }, + 3.0), + new LinearMeasurement(new double[] { 2.0, 1.0 }, + new EstimatedParameter[] { p[0], p[2] }, + 5.0) + }); + + GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); + try { + estimator.estimate(problem); + fail("an exception should have been caught"); + } catch (EstimationException ee) { + // expected behavior + } catch (Exception e) { + fail("wrong exception type caught"); + } + + } + + public void testMoreEstimatedParametersUnsorted() throws EstimationException { + EstimatedParameter[] p = { + new EstimatedParameter("p0", 2), + new EstimatedParameter("p1", 2), + new EstimatedParameter("p2", 2), + new EstimatedParameter("p3", 2), + new EstimatedParameter("p4", 2), + new EstimatedParameter("p5", 2) + }; + LinearProblem problem = new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] { 1.0, 1.0 }, + new EstimatedParameter[] { p[0], p[1] }, + 3.0), + new LinearMeasurement(new double[] { 1.0, 1.0, 1.0 }, + new EstimatedParameter[] { p[2], p[3], p[4] }, + 12.0), + new LinearMeasurement(new double[] { 1.0, -1.0 }, + new EstimatedParameter[] { p[4], p[5] }, + -1.0), + new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 }, + new EstimatedParameter[] { p[3], p[2], p[5] }, + 7.0), + new LinearMeasurement(new double[] { 1.0, -1.0 }, + new EstimatedParameter[] { p[4], p[3] }, + 1.0) + }); + + GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); + try { + estimator.estimate(problem); + fail("an exception should have been caught"); + } catch (EstimationException ee) { + // expected behavior + } catch (Exception e) { + fail("wrong exception type caught"); + } + + } + + public void testRedundantEquations() throws EstimationException { + EstimatedParameter[] p = { + new EstimatedParameter("p0", 1), + new EstimatedParameter("p1", 1) + }; + LinearProblem problem = new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] { 1.0, 1.0 }, + new EstimatedParameter[] { p[0], p[1] }, + 3.0), + new LinearMeasurement(new double[] { 1.0, -1.0 }, + new EstimatedParameter[] { p[0], p[1] }, + 1.0), + new LinearMeasurement(new double[] { 1.0, 3.0 }, + new EstimatedParameter[] { p[0], p[1] }, + 5.0) + }); + + GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); + estimator.estimate(problem); + assertEquals(0, estimator.getRMS(problem), 1.0e-10); + EstimatedParameter[] all = problem.getAllParameters(); + for (int i = 0; i < all.length; ++i) { + assertEquals(all[i].getName().equals("p0") ? 2.0 : 1.0, + all[i].getEstimate(), 1.0e-10); + } + + } + + public void testInconsistentEquations() throws EstimationException { + EstimatedParameter[] p = { + new EstimatedParameter("p0", 1), + new EstimatedParameter("p1", 1) + }; + LinearProblem problem = new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] { 1.0, 1.0 }, + new EstimatedParameter[] { p[0], p[1] }, + 3.0), + new LinearMeasurement(new double[] { 1.0, -1.0 }, + new EstimatedParameter[] { p[0], p[1] }, + 1.0), + new LinearMeasurement(new double[] { 1.0, 3.0 }, + new EstimatedParameter[] { p[0], p[1] }, + 4.0) + }); + + GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); + estimator.estimate(problem); + assertTrue(estimator.getRMS(problem) > 0.1); + + } + + public void testBoundParameters() throws EstimationException { + EstimatedParameter[] p = { + new EstimatedParameter("unbound0", 2, false), + new EstimatedParameter("unbound1", 2, false), + new EstimatedParameter("bound", 2, true) + }; + LinearProblem problem = new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] { 1.0, 1.0, 1.0 }, + new EstimatedParameter[] { p[0], p[1], p[2] }, + 3.0), + new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 }, + new EstimatedParameter[] { p[0], p[1], p[2] }, + 1.0), + new LinearMeasurement(new double[] { 1.0, 3.0, 2.0 }, + new EstimatedParameter[] { p[0], p[1], p[2] }, + 7.0) + }); + + GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); + estimator.estimate(problem); + assertTrue(estimator.getRMS(problem) < 1.0e-10); + double[][] covariances = estimator.getCovariances(problem); + int i0 = 0, i1 = 1; + if (problem.getUnboundParameters()[0].getName().endsWith("1")) { + i0 = 1; + i1 = 0; + } + assertEquals(11.0 / 24, covariances[i0][i0], 1.0e-10); + assertEquals(-3.0 / 24, covariances[i0][i1], 1.0e-10); + assertEquals(-3.0 / 24, covariances[i1][i0], 1.0e-10); + assertEquals( 3.0 / 24, covariances[i1][i1], 1.0e-10); + + double[] errors = estimator.guessParametersErrors(problem); + assertEquals(0, errors[i0], 1.0e-10); + assertEquals(0, errors[i1], 1.0e-10); + + } + + public void testMaxIterations() { + Circle circle = new Circle(98.680, 47.345); + circle.addPoint( 30.0, 68.0); + circle.addPoint( 50.0, -6.0); + circle.addPoint(110.0, -20.0); + circle.addPoint( 35.0, 15.0); + circle.addPoint( 45.0, 97.0); + try { + GaussNewtonEstimator estimator = new GaussNewtonEstimator(4, 1.0e-14, 1.0e-14); + estimator.estimate(circle); + fail("an exception should have been caught"); + } catch (EstimationException ee) { + // expected behavior + } catch (Exception e) { + fail("wrong exception type caught"); + } + } + + public void testCircleFitting() throws EstimationException { + Circle circle = new Circle(98.680, 47.345); + circle.addPoint( 30.0, 68.0); + circle.addPoint( 50.0, -6.0); + circle.addPoint(110.0, -20.0); + circle.addPoint( 35.0, 15.0); + circle.addPoint( 45.0, 97.0); + GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-10, 1.0e-10); + estimator.estimate(circle); + double rms = estimator.getRMS(circle); + assertEquals(1.768262623567235, Math.sqrt(circle.getM()) * rms, 1.0e-10); + assertEquals(69.96016176931406, circle.getRadius(), 1.0e-10); + assertEquals(96.07590211815305, circle.getX(), 1.0e-10); + assertEquals(48.13516790438953, circle.getY(), 1.0e-10); + } + + public void testCircleFittingBadInit() throws EstimationException { + Circle circle = new Circle(-12, -12); + double[][] points = new double[][] { + {-0.312967, 0.072366}, {-0.339248, 0.132965}, {-0.379780, 0.202724}, + {-0.390426, 0.260487}, {-0.361212, 0.328325}, {-0.346039, 0.392619}, + {-0.280579, 0.444306}, {-0.216035, 0.470009}, {-0.149127, 0.493832}, + {-0.075133, 0.483271}, {-0.007759, 0.452680}, { 0.060071, 0.410235}, + { 0.103037, 0.341076}, { 0.118438, 0.273884}, { 0.131293, 0.192201}, + { 0.115869, 0.129797}, { 0.072223, 0.058396}, { 0.022884, 0.000718}, + {-0.053355, -0.020405}, {-0.123584, -0.032451}, {-0.216248, -0.032862}, + {-0.278592, -0.005008}, {-0.337655, 0.056658}, {-0.385899, 0.112526}, + {-0.405517, 0.186957}, {-0.415374, 0.262071}, {-0.387482, 0.343398}, + {-0.347322, 0.397943}, {-0.287623, 0.458425}, {-0.223502, 0.475513}, + {-0.135352, 0.478186}, {-0.061221, 0.483371}, { 0.003711, 0.422737}, + { 0.065054, 0.375830}, { 0.108108, 0.297099}, { 0.123882, 0.222850}, + { 0.117729, 0.134382}, { 0.085195, 0.056820}, { 0.029800, -0.019138}, + {-0.027520, -0.072374}, {-0.102268, -0.091555}, {-0.200299, -0.106578}, + {-0.292731, -0.091473}, {-0.356288, -0.051108}, {-0.420561, 0.014926}, + {-0.471036, 0.074716}, {-0.488638, 0.182508}, {-0.485990, 0.254068}, + {-0.463943, 0.338438}, {-0.406453, 0.404704}, {-0.334287, 0.466119}, + {-0.254244, 0.503188}, {-0.161548, 0.495769}, {-0.075733, 0.495560}, + { 0.001375, 0.434937}, { 0.082787, 0.385806}, { 0.115490, 0.323807}, + { 0.141089, 0.223450}, { 0.138693, 0.131703}, { 0.126415, 0.049174}, + { 0.066518, -0.010217}, {-0.005184, -0.070647}, {-0.080985, -0.103635}, + {-0.177377, -0.116887}, {-0.260628, -0.100258}, {-0.335756, -0.056251}, + {-0.405195, -0.000895}, {-0.444937, 0.085456}, {-0.484357, 0.175597}, + {-0.472453, 0.248681}, {-0.438580, 0.347463}, {-0.402304, 0.422428}, + {-0.326777, 0.479438}, {-0.247797, 0.505581}, {-0.152676, 0.519380}, + {-0.071754, 0.516264}, { 0.015942, 0.472802}, { 0.076608, 0.419077}, + { 0.127673, 0.330264}, { 0.159951, 0.262150}, { 0.153530, 0.172681}, + { 0.140653, 0.089229}, { 0.078666, 0.024981}, { 0.023807, -0.037022}, + {-0.048837, -0.077056}, {-0.127729, -0.075338}, {-0.221271, -0.067526} + }; + for (int i = 0; i < points.length; ++i) { + circle.addPoint(points[i][0], points[i][1]); + } + GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); + try { + estimator.estimate(circle); + fail("an exception should have been caught"); + } catch (EstimationException ee) { + // expected behavior + } catch (Exception e) { + fail("wrong exception type caught"); + } +} + + private static class LinearProblem extends SimpleEstimationProblem { + + public LinearProblem(LinearMeasurement[] measurements) { + HashSet set = new HashSet(); + for (int i = 0; i < measurements.length; ++i) { + addMeasurement(measurements[i]); + EstimatedParameter[] parameters = measurements[i].getParameters(); + for (int j = 0; j < parameters.length; ++j) { + set.add(parameters[j]); + } + } + for (EstimatedParameter p : set) { + addParameter(p); + } + } + + } + + private static class LinearMeasurement extends WeightedMeasurement { + + public LinearMeasurement(double[] factors, EstimatedParameter[] parameters, + double setPoint) { + super(1.0, setPoint, true); + this.factors = factors; + this.parameters = parameters; + setIgnored(false); + } + + public double getTheoreticalValue() { + double v = 0; + for (int i = 0; i < factors.length; ++i) { + v += factors[i] * parameters[i].getEstimate(); + } + return v; + } + + public double getPartial(EstimatedParameter parameter) { + for (int i = 0; i < parameters.length; ++i) { + if (parameters[i] == parameter) { + return factors[i]; + } + } + return 0; + } + + public EstimatedParameter[] getParameters() { + return parameters; + } + + private double[] factors; + private EstimatedParameter[] parameters; + private static final long serialVersionUID = -3922448707008868580L; + + } + + private static class Circle implements EstimationProblem { + + public Circle(double cx, double cy) { + this.cx = new EstimatedParameter("cx", cx); + this.cy = new EstimatedParameter(new EstimatedParameter("cy", cy)); + points = new ArrayList(); + } + + public void addPoint(double px, double py) { + points.add(new PointModel(px, py)); + } + + public int getM() { + return points.size(); + } + + public WeightedMeasurement[] getMeasurements() { + return (WeightedMeasurement[]) points.toArray(new PointModel[points.size()]); + } + + public EstimatedParameter[] getAllParameters() { + return new EstimatedParameter[] { cx, cy }; + } + + public EstimatedParameter[] getUnboundParameters() { + return new EstimatedParameter[] { cx, cy }; + } + + public double getPartialRadiusX() { + double dRdX = 0; + for (PointModel point : points) { + dRdX += point.getPartialDiX(); + } + return dRdX / points.size(); + } + + public double getPartialRadiusY() { + double dRdY = 0; + for (PointModel point : points) { + dRdY += point.getPartialDiY(); + } + return dRdY / points.size(); + } + + public double getRadius() { + double r = 0; + for (PointModel point : points) { + r += point.getCenterDistance(); + } + return r / points.size(); + } + + public double getX() { + return cx.getEstimate(); + } + + public double getY() { + return cy.getEstimate(); + } + + private class PointModel extends WeightedMeasurement { + + public PointModel(double px, double py) { + super(1.0, 0.0); + this.px = px; + this.py = py; + } + + public double getPartial(EstimatedParameter parameter) { + if (parameter == cx) { + return getPartialDiX() - getPartialRadiusX(); + } else if (parameter == cy) { + return getPartialDiY() - getPartialRadiusY(); + } + return 0; + } + + public double getCenterDistance() { + double dx = px - cx.getEstimate(); + double dy = py - cy.getEstimate(); + return Math.sqrt(dx * dx + dy * dy); + } + + public double getPartialDiX() { + return (cx.getEstimate() - px) / getCenterDistance(); + } + + public double getPartialDiY() { + return (cy.getEstimate() - py) / getCenterDistance(); + } + + public double getTheoreticalValue() { + return getCenterDistance() - getRadius(); + } + + private double px; + private double py; + private static final long serialVersionUID = 1L; + + } + + private EstimatedParameter cx; + private EstimatedParameter cy; + private ArrayList points; + + } + + public static Test suite() { + return new TestSuite(GaussNewtonEstimatorTest.class); + } + +} diff --git a/src/test/org/apache/commons/math/estimation/LevenbergMarquardtEstimatorTest.java b/src/test/org/apache/commons/math/estimation/LevenbergMarquardtEstimatorTest.java new file mode 100644 index 000000000..512455572 --- /dev/null +++ b/src/test/org/apache/commons/math/estimation/LevenbergMarquardtEstimatorTest.java @@ -0,0 +1,844 @@ +/* + * 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.estimation; + +import java.util.ArrayList; +import java.util.HashSet; + +import junit.framework.Test; +import junit.framework.TestCase; +import junit.framework.TestSuite; + +/** + *

      Some of the unit tests are re-implementations of the MINPACK file17 and file22 test files. + * The redistribution policy for MINPACK is available here, for + * convenience, it is reproduced below.

      + + * + * + * + *
      + * Minpack Copyright Notice (1999) University of Chicago. + * All rights reserved + *
      + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + *
        + *
      1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer.
        2. + *
      2. Redistributions in binary form must reproduce the above + * copyright notice, this list of conditions and the following + * disclaimer in the documentation and/or other materials provided + * with the distribution.
        3. + *
      3. The end-user documentation included with the redistribution, if any, + * must include the following acknowledgment: + * This product includes software developed by the University of + * Chicago, as Operator of Argonne National Laboratory. + * Alternately, this acknowledgment may appear in the software itself, + * if and wherever such third-party acknowledgments normally appear.
        4. + *
      4. WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS" + * WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE + * UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND + * THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES + * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE + * OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY + * OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR + * USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF + * THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4) + * DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION + * UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL + * BE CORRECTED.
        5. + *
      5. LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT + * HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF + * ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT, + * INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF + * ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF + * PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER + * SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT + * (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE, + * EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE + * POSSIBILITY OF SUCH LOSS OR DAMAGES.
        6. + *
        + + * @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran minpack tests) + * @author Burton S. Garbow (original fortran minpack tests) + * @author Kenneth E. Hillstrom (original fortran minpack tests) + * @author Jorge J. More (original fortran minpack tests) + * @author Luc Maisonobe (non-minpack tests and minpack tests Java translation) + */ +@Deprecated +public class LevenbergMarquardtEstimatorTest + extends TestCase { + + public LevenbergMarquardtEstimatorTest(String name) { + super(name); + } + + public void testTrivial() throws EstimationException { + LinearProblem problem = + new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] {2}, + new EstimatedParameter[] { + new EstimatedParameter("p0", 0) + }, 3.0) + }); + LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); + estimator.estimate(problem); + assertEquals(0, estimator.getRMS(problem), 1.0e-10); + try { + estimator.guessParametersErrors(problem); + fail("an exception should have been thrown"); + } catch (EstimationException ee) { + // expected behavior + } catch (Exception e) { + fail("wrong exception caught"); + } + assertEquals(1.5, + problem.getUnboundParameters()[0].getEstimate(), + 1.0e-10); + } + + public void testQRColumnsPermutation() throws EstimationException { + + EstimatedParameter[] x = { + new EstimatedParameter("p0", 0), new EstimatedParameter("p1", 0) + }; + LinearProblem problem = new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] { 1.0, -1.0 }, + new EstimatedParameter[] { x[0], x[1] }, + 4.0), + new LinearMeasurement(new double[] { 2.0 }, + new EstimatedParameter[] { x[1] }, + 6.0), + new LinearMeasurement(new double[] { 1.0, -2.0 }, + new EstimatedParameter[] { x[0], x[1] }, + 1.0) + }); + + LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); + estimator.estimate(problem); + assertEquals(0, estimator.getRMS(problem), 1.0e-10); + assertEquals(7.0, x[0].getEstimate(), 1.0e-10); + assertEquals(3.0, x[1].getEstimate(), 1.0e-10); + + } + + public void testNoDependency() throws EstimationException { + EstimatedParameter[] p = new EstimatedParameter[] { + new EstimatedParameter("p0", 0), + new EstimatedParameter("p1", 0), + new EstimatedParameter("p2", 0), + new EstimatedParameter("p3", 0), + new EstimatedParameter("p4", 0), + new EstimatedParameter("p5", 0) + }; + LinearProblem problem = new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[0] }, 0.0), + new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[1] }, 1.1), + new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[2] }, 2.2), + new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[3] }, 3.3), + new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[4] }, 4.4), + new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[5] }, 5.5) + }); + LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); + estimator.estimate(problem); + assertEquals(0, estimator.getRMS(problem), 1.0e-10); + for (int i = 0; i < p.length; ++i) { + assertEquals(0.55 * i, p[i].getEstimate(), 1.0e-10); + } +} + + public void testOneSet() throws EstimationException { + + EstimatedParameter[] p = { + new EstimatedParameter("p0", 0), + new EstimatedParameter("p1", 0), + new EstimatedParameter("p2", 0) + }; + LinearProblem problem = new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] { 1.0 }, + new EstimatedParameter[] { p[0] }, + 1.0), + new LinearMeasurement(new double[] { -1.0, 1.0 }, + new EstimatedParameter[] { p[0], p[1] }, + 1.0), + new LinearMeasurement(new double[] { -1.0, 1.0 }, + new EstimatedParameter[] { p[1], p[2] }, + 1.0) + }); + + LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); + estimator.estimate(problem); + assertEquals(0, estimator.getRMS(problem), 1.0e-10); + assertEquals(1.0, p[0].getEstimate(), 1.0e-10); + assertEquals(2.0, p[1].getEstimate(), 1.0e-10); + assertEquals(3.0, p[2].getEstimate(), 1.0e-10); + + } + + public void testTwoSets() throws EstimationException { + EstimatedParameter[] p = { + new EstimatedParameter("p0", 0), + new EstimatedParameter("p1", 1), + new EstimatedParameter("p2", 2), + new EstimatedParameter("p3", 3), + new EstimatedParameter("p4", 4), + new EstimatedParameter("p5", 5) + }; + + double epsilon = 1.0e-7; + LinearProblem problem = new LinearProblem(new LinearMeasurement[] { + + // 4 elements sub-problem + new LinearMeasurement(new double[] { 2.0, 1.0, 4.0 }, + new EstimatedParameter[] { p[0], p[1], p[3] }, + 2.0), + new LinearMeasurement(new double[] { -4.0, -2.0, 3.0, -7.0 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + -9.0), + new LinearMeasurement(new double[] { 4.0, 1.0, -2.0, 8.0 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + 2.0), + new LinearMeasurement(new double[] { -3.0, -12.0, -1.0 }, + new EstimatedParameter[] { p[1], p[2], p[3] }, + 2.0), + + // 2 elements sub-problem + new LinearMeasurement(new double[] { epsilon, 1.0 }, + new EstimatedParameter[] { p[4], p[5] }, + 1.0 + epsilon * epsilon), + new LinearMeasurement(new double[] { 1.0, 1.0 }, + new EstimatedParameter[] { p[4], p[5] }, + 2.0) + + }); + + LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); + estimator.estimate(problem); + assertEquals(0, estimator.getRMS(problem), 1.0e-10); + assertEquals( 3.0, p[0].getEstimate(), 1.0e-10); + assertEquals( 4.0, p[1].getEstimate(), 1.0e-10); + assertEquals(-1.0, p[2].getEstimate(), 1.0e-10); + assertEquals(-2.0, p[3].getEstimate(), 1.0e-10); + assertEquals( 1.0 + epsilon, p[4].getEstimate(), 1.0e-10); + assertEquals( 1.0 - epsilon, p[5].getEstimate(), 1.0e-10); + + } + + public void testNonInversible() throws EstimationException { + + EstimatedParameter[] p = { + new EstimatedParameter("p0", 0), + new EstimatedParameter("p1", 0), + new EstimatedParameter("p2", 0) + }; + LinearMeasurement[] m = new LinearMeasurement[] { + new LinearMeasurement(new double[] { 1.0, 2.0, -3.0 }, + new EstimatedParameter[] { p[0], p[1], p[2] }, + 1.0), + new LinearMeasurement(new double[] { 2.0, 1.0, 3.0 }, + new EstimatedParameter[] { p[0], p[1], p[2] }, + 1.0), + new LinearMeasurement(new double[] { -3.0, -9.0 }, + new EstimatedParameter[] { p[0], p[2] }, + 1.0) + }; + LinearProblem problem = new LinearProblem(m); + + LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); + double initialCost = estimator.getRMS(problem); + estimator.estimate(problem); + assertTrue(estimator.getRMS(problem) < initialCost); + assertTrue(Math.sqrt(m.length) * estimator.getRMS(problem) > 0.6); + try { + estimator.getCovariances(problem); + fail("an exception should have been thrown"); + } catch (EstimationException ee) { + // expected behavior + } catch (Exception e) { + fail("wrong exception caught"); + } + double dJ0 = 2 * (m[0].getResidual() * m[0].getPartial(p[0]) + + m[1].getResidual() * m[1].getPartial(p[0]) + + m[2].getResidual() * m[2].getPartial(p[0])); + double dJ1 = 2 * (m[0].getResidual() * m[0].getPartial(p[1]) + + m[1].getResidual() * m[1].getPartial(p[1])); + double dJ2 = 2 * (m[0].getResidual() * m[0].getPartial(p[2]) + + m[1].getResidual() * m[1].getPartial(p[2]) + + m[2].getResidual() * m[2].getPartial(p[2])); + assertEquals(0, dJ0, 1.0e-10); + assertEquals(0, dJ1, 1.0e-10); + assertEquals(0, dJ2, 1.0e-10); + + } + + public void testIllConditioned() throws EstimationException { + EstimatedParameter[] p = { + new EstimatedParameter("p0", 0), + new EstimatedParameter("p1", 1), + new EstimatedParameter("p2", 2), + new EstimatedParameter("p3", 3) + }; + + LinearProblem problem1 = new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] { 10.0, 7.0, 8.0, 7.0 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + 32.0), + new LinearMeasurement(new double[] { 7.0, 5.0, 6.0, 5.0 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + 23.0), + new LinearMeasurement(new double[] { 8.0, 6.0, 10.0, 9.0 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + 33.0), + new LinearMeasurement(new double[] { 7.0, 5.0, 9.0, 10.0 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + 31.0) + }); + LevenbergMarquardtEstimator estimator1 = new LevenbergMarquardtEstimator(); + estimator1.estimate(problem1); + assertEquals(0, estimator1.getRMS(problem1), 1.0e-10); + assertEquals(1.0, p[0].getEstimate(), 1.0e-10); + assertEquals(1.0, p[1].getEstimate(), 1.0e-10); + assertEquals(1.0, p[2].getEstimate(), 1.0e-10); + assertEquals(1.0, p[3].getEstimate(), 1.0e-10); + + LinearProblem problem2 = new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] { 10.0, 7.0, 8.1, 7.2 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + 32.0), + new LinearMeasurement(new double[] { 7.08, 5.04, 6.0, 5.0 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + 23.0), + new LinearMeasurement(new double[] { 8.0, 5.98, 9.89, 9.0 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + 33.0), + new LinearMeasurement(new double[] { 6.99, 4.99, 9.0, 9.98 }, + new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, + 31.0) + }); + LevenbergMarquardtEstimator estimator2 = new LevenbergMarquardtEstimator(); + estimator2.estimate(problem2); + assertEquals(0, estimator2.getRMS(problem2), 1.0e-10); + assertEquals(-81.0, p[0].getEstimate(), 1.0e-8); + assertEquals(137.0, p[1].getEstimate(), 1.0e-8); + assertEquals(-34.0, p[2].getEstimate(), 1.0e-8); + assertEquals( 22.0, p[3].getEstimate(), 1.0e-8); + + } + + public void testMoreEstimatedParametersSimple() throws EstimationException { + + EstimatedParameter[] p = { + new EstimatedParameter("p0", 7), + new EstimatedParameter("p1", 6), + new EstimatedParameter("p2", 5), + new EstimatedParameter("p3", 4) + }; + LinearProblem problem = new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] { 3.0, 2.0 }, + new EstimatedParameter[] { p[0], p[1] }, + 7.0), + new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 }, + new EstimatedParameter[] { p[1], p[2], p[3] }, + 3.0), + new LinearMeasurement(new double[] { 2.0, 1.0 }, + new EstimatedParameter[] { p[0], p[2] }, + 5.0) + }); + + LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); + estimator.estimate(problem); + assertEquals(0, estimator.getRMS(problem), 1.0e-10); + + } + + public void testMoreEstimatedParametersUnsorted() throws EstimationException { + EstimatedParameter[] p = { + new EstimatedParameter("p0", 2), + new EstimatedParameter("p1", 2), + new EstimatedParameter("p2", 2), + new EstimatedParameter("p3", 2), + new EstimatedParameter("p4", 2), + new EstimatedParameter("p5", 2) + }; + LinearProblem problem = new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] { 1.0, 1.0 }, + new EstimatedParameter[] { p[0], p[1] }, + 3.0), + new LinearMeasurement(new double[] { 1.0, 1.0, 1.0 }, + new EstimatedParameter[] { p[2], p[3], p[4] }, + 12.0), + new LinearMeasurement(new double[] { 1.0, -1.0 }, + new EstimatedParameter[] { p[4], p[5] }, + -1.0), + new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 }, + new EstimatedParameter[] { p[3], p[2], p[5] }, + 7.0), + new LinearMeasurement(new double[] { 1.0, -1.0 }, + new EstimatedParameter[] { p[4], p[3] }, + 1.0) + }); + + LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); + estimator.estimate(problem); + assertEquals(0, estimator.getRMS(problem), 1.0e-10); + assertEquals(3.0, p[2].getEstimate(), 1.0e-10); + assertEquals(4.0, p[3].getEstimate(), 1.0e-10); + assertEquals(5.0, p[4].getEstimate(), 1.0e-10); + assertEquals(6.0, p[5].getEstimate(), 1.0e-10); + + } + + public void testRedundantEquations() throws EstimationException { + EstimatedParameter[] p = { + new EstimatedParameter("p0", 1), + new EstimatedParameter("p1", 1) + }; + LinearProblem problem = new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] { 1.0, 1.0 }, + new EstimatedParameter[] { p[0], p[1] }, + 3.0), + new LinearMeasurement(new double[] { 1.0, -1.0 }, + new EstimatedParameter[] { p[0], p[1] }, + 1.0), + new LinearMeasurement(new double[] { 1.0, 3.0 }, + new EstimatedParameter[] { p[0], p[1] }, + 5.0) + }); + + LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); + estimator.estimate(problem); + assertEquals(0, estimator.getRMS(problem), 1.0e-10); + assertEquals(2.0, p[0].getEstimate(), 1.0e-10); + assertEquals(1.0, p[1].getEstimate(), 1.0e-10); + + } + + public void testInconsistentEquations() throws EstimationException { + EstimatedParameter[] p = { + new EstimatedParameter("p0", 1), + new EstimatedParameter("p1", 1) + }; + LinearProblem problem = new LinearProblem(new LinearMeasurement[] { + new LinearMeasurement(new double[] { 1.0, 1.0 }, + new EstimatedParameter[] { p[0], p[1] }, + 3.0), + new LinearMeasurement(new double[] { 1.0, -1.0 }, + new EstimatedParameter[] { p[0], p[1] }, + 1.0), + new LinearMeasurement(new double[] { 1.0, 3.0 }, + new EstimatedParameter[] { p[0], p[1] }, + 4.0) + }); + + LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); + estimator.estimate(problem); + assertTrue(estimator.getRMS(problem) > 0.1); + + } + + public void testControlParameters() throws EstimationException { + Circle circle = new Circle(98.680, 47.345); + circle.addPoint( 30.0, 68.0); + circle.addPoint( 50.0, -6.0); + circle.addPoint(110.0, -20.0); + circle.addPoint( 35.0, 15.0); + circle.addPoint( 45.0, 97.0); + checkEstimate(circle, 0.1, 10, 1.0e-14, 1.0e-16, 1.0e-10, false); + checkEstimate(circle, 0.1, 10, 1.0e-15, 1.0e-17, 1.0e-10, true); + checkEstimate(circle, 0.1, 5, 1.0e-15, 1.0e-16, 1.0e-10, true); + circle.addPoint(300, -300); + checkEstimate(circle, 0.1, 20, 1.0e-18, 1.0e-16, 1.0e-10, true); + } + + private void checkEstimate(EstimationProblem problem, + double initialStepBoundFactor, int maxCostEval, + double costRelativeTolerance, double parRelativeTolerance, + double orthoTolerance, boolean shouldFail) { + try { + LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); + estimator.setInitialStepBoundFactor(initialStepBoundFactor); + estimator.setMaxCostEval(maxCostEval); + estimator.setCostRelativeTolerance(costRelativeTolerance); + estimator.setParRelativeTolerance(parRelativeTolerance); + estimator.setOrthoTolerance(orthoTolerance); + estimator.estimate(problem); + assertTrue(! shouldFail); + } catch (EstimationException ee) { + assertTrue(shouldFail); + } catch (Exception e) { + fail("wrong exception type caught"); + } + } + + public void testCircleFitting() throws EstimationException { + Circle circle = new Circle(98.680, 47.345); + circle.addPoint( 30.0, 68.0); + circle.addPoint( 50.0, -6.0); + circle.addPoint(110.0, -20.0); + circle.addPoint( 35.0, 15.0); + circle.addPoint( 45.0, 97.0); + LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); + estimator.estimate(circle); + assertTrue(estimator.getCostEvaluations() < 10); + assertTrue(estimator.getJacobianEvaluations() < 10); + double rms = estimator.getRMS(circle); + assertEquals(1.768262623567235, Math.sqrt(circle.getM()) * rms, 1.0e-10); + assertEquals(69.96016176931406, circle.getRadius(), 1.0e-10); + assertEquals(96.07590211815305, circle.getX(), 1.0e-10); + assertEquals(48.13516790438953, circle.getY(), 1.0e-10); + double[][] cov = estimator.getCovariances(circle); + assertEquals(1.839, cov[0][0], 0.001); + assertEquals(0.731, cov[0][1], 0.001); + assertEquals(cov[0][1], cov[1][0], 1.0e-14); + assertEquals(0.786, cov[1][1], 0.001); + double[] errors = estimator.guessParametersErrors(circle); + assertEquals(1.384, errors[0], 0.001); + assertEquals(0.905, errors[1], 0.001); + + // add perfect measurements and check errors are reduced + double cx = circle.getX(); + double cy = circle.getY(); + double r = circle.getRadius(); + for (double d= 0; d < 2 * Math.PI; d += 0.01) { + circle.addPoint(cx + r * Math.cos(d), cy + r * Math.sin(d)); + } + estimator = new LevenbergMarquardtEstimator(); + estimator.estimate(circle); + cov = estimator.getCovariances(circle); + assertEquals(0.004, cov[0][0], 0.001); + assertEquals(6.40e-7, cov[0][1], 1.0e-9); + assertEquals(cov[0][1], cov[1][0], 1.0e-14); + assertEquals(0.003, cov[1][1], 0.001); + errors = estimator.guessParametersErrors(circle); + assertEquals(0.004, errors[0], 0.001); + assertEquals(0.004, errors[1], 0.001); + + } + + public void testCircleFittingBadInit() throws EstimationException { + Circle circle = new Circle(-12, -12); + double[][] points = new double[][] { + {-0.312967, 0.072366}, {-0.339248, 0.132965}, {-0.379780, 0.202724}, + {-0.390426, 0.260487}, {-0.361212, 0.328325}, {-0.346039, 0.392619}, + {-0.280579, 0.444306}, {-0.216035, 0.470009}, {-0.149127, 0.493832}, + {-0.075133, 0.483271}, {-0.007759, 0.452680}, { 0.060071, 0.410235}, + { 0.103037, 0.341076}, { 0.118438, 0.273884}, { 0.131293, 0.192201}, + { 0.115869, 0.129797}, { 0.072223, 0.058396}, { 0.022884, 0.000718}, + {-0.053355, -0.020405}, {-0.123584, -0.032451}, {-0.216248, -0.032862}, + {-0.278592, -0.005008}, {-0.337655, 0.056658}, {-0.385899, 0.112526}, + {-0.405517, 0.186957}, {-0.415374, 0.262071}, {-0.387482, 0.343398}, + {-0.347322, 0.397943}, {-0.287623, 0.458425}, {-0.223502, 0.475513}, + {-0.135352, 0.478186}, {-0.061221, 0.483371}, { 0.003711, 0.422737}, + { 0.065054, 0.375830}, { 0.108108, 0.297099}, { 0.123882, 0.222850}, + { 0.117729, 0.134382}, { 0.085195, 0.056820}, { 0.029800, -0.019138}, + {-0.027520, -0.072374}, {-0.102268, -0.091555}, {-0.200299, -0.106578}, + {-0.292731, -0.091473}, {-0.356288, -0.051108}, {-0.420561, 0.014926}, + {-0.471036, 0.074716}, {-0.488638, 0.182508}, {-0.485990, 0.254068}, + {-0.463943, 0.338438}, {-0.406453, 0.404704}, {-0.334287, 0.466119}, + {-0.254244, 0.503188}, {-0.161548, 0.495769}, {-0.075733, 0.495560}, + { 0.001375, 0.434937}, { 0.082787, 0.385806}, { 0.115490, 0.323807}, + { 0.141089, 0.223450}, { 0.138693, 0.131703}, { 0.126415, 0.049174}, + { 0.066518, -0.010217}, {-0.005184, -0.070647}, {-0.080985, -0.103635}, + {-0.177377, -0.116887}, {-0.260628, -0.100258}, {-0.335756, -0.056251}, + {-0.405195, -0.000895}, {-0.444937, 0.085456}, {-0.484357, 0.175597}, + {-0.472453, 0.248681}, {-0.438580, 0.347463}, {-0.402304, 0.422428}, + {-0.326777, 0.479438}, {-0.247797, 0.505581}, {-0.152676, 0.519380}, + {-0.071754, 0.516264}, { 0.015942, 0.472802}, { 0.076608, 0.419077}, + { 0.127673, 0.330264}, { 0.159951, 0.262150}, { 0.153530, 0.172681}, + { 0.140653, 0.089229}, { 0.078666, 0.024981}, { 0.023807, -0.037022}, + {-0.048837, -0.077056}, {-0.127729, -0.075338}, {-0.221271, -0.067526} + }; + for (int i = 0; i < points.length; ++i) { + circle.addPoint(points[i][0], points[i][1]); + } + LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); + estimator.estimate(circle); + assertTrue(estimator.getCostEvaluations() < 15); + assertTrue(estimator.getJacobianEvaluations() < 10); + assertEquals( 0.030184491196225207, estimator.getRMS(circle), 1.0e-9); + assertEquals( 0.2922350065939634, circle.getRadius(), 1.0e-9); + assertEquals(-0.15173845023862165, circle.getX(), 1.0e-8); + assertEquals( 0.20750021499570379, circle.getY(), 1.0e-8); + } + + public void testMath199() { + try { + QuadraticProblem problem = new QuadraticProblem(); + problem.addPoint (0, -3.182591015485607, 0.0); + problem.addPoint (1, -2.5581184967730577, 4.4E-323); + problem.addPoint (2, -2.1488478161387325, 1.0); + problem.addPoint (3, -1.9122489313410047, 4.4E-323); + problem.addPoint (4, 1.7785661310051026, 0.0); + new LevenbergMarquardtEstimator().estimate(problem); + fail("an exception should have been thrown"); + } catch (EstimationException ee) { + // expected behavior + } + + } + + private static class LinearProblem implements EstimationProblem { + + public LinearProblem(LinearMeasurement[] measurements) { + this.measurements = measurements; + } + + public WeightedMeasurement[] getMeasurements() { + return measurements; + } + + public EstimatedParameter[] getUnboundParameters() { + return getAllParameters(); + } + + public EstimatedParameter[] getAllParameters() { + HashSet set = new HashSet(); + for (int i = 0; i < measurements.length; ++i) { + EstimatedParameter[] parameters = measurements[i].getParameters(); + for (int j = 0; j < parameters.length; ++j) { + set.add(parameters[j]); + } + } + return (EstimatedParameter[]) set.toArray(new EstimatedParameter[set.size()]); + } + + private LinearMeasurement[] measurements; + + } + + private static class LinearMeasurement extends WeightedMeasurement { + + public LinearMeasurement(double[] factors, EstimatedParameter[] parameters, + double setPoint) { + super(1.0, setPoint); + this.factors = factors; + this.parameters = parameters; + } + + public double getTheoreticalValue() { + double v = 0; + for (int i = 0; i < factors.length; ++i) { + v += factors[i] * parameters[i].getEstimate(); + } + return v; + } + + public double getPartial(EstimatedParameter parameter) { + for (int i = 0; i < parameters.length; ++i) { + if (parameters[i] == parameter) { + return factors[i]; + } + } + return 0; + } + + public EstimatedParameter[] getParameters() { + return parameters; + } + + private double[] factors; + private EstimatedParameter[] parameters; + private static final long serialVersionUID = -3922448707008868580L; + + } + + private static class Circle implements EstimationProblem { + + public Circle(double cx, double cy) { + this.cx = new EstimatedParameter("cx", cx); + this.cy = new EstimatedParameter("cy", cy); + points = new ArrayList(); + } + + public void addPoint(double px, double py) { + points.add(new PointModel(px, py)); + } + + public int getM() { + return points.size(); + } + + public WeightedMeasurement[] getMeasurements() { + return (WeightedMeasurement[]) points.toArray(new PointModel[points.size()]); + } + + public EstimatedParameter[] getAllParameters() { + return new EstimatedParameter[] { cx, cy }; + } + + public EstimatedParameter[] getUnboundParameters() { + return new EstimatedParameter[] { cx, cy }; + } + + public double getPartialRadiusX() { + double dRdX = 0; + for (PointModel point : points) { + dRdX += point.getPartialDiX(); + } + return dRdX / points.size(); + } + + public double getPartialRadiusY() { + double dRdY = 0; + for (PointModel point : points) { + dRdY += point.getPartialDiY(); + } + return dRdY / points.size(); + } + + public double getRadius() { + double r = 0; + for (PointModel point : points) { + r += point.getCenterDistance(); + } + return r / points.size(); + } + + public double getX() { + return cx.getEstimate(); + } + + public double getY() { + return cy.getEstimate(); + } + + private class PointModel extends WeightedMeasurement { + + public PointModel(double px, double py) { + super(1.0, 0.0); + this.px = px; + this.py = py; + } + + public double getPartial(EstimatedParameter parameter) { + if (parameter == cx) { + return getPartialDiX() - getPartialRadiusX(); + } else if (parameter == cy) { + return getPartialDiY() - getPartialRadiusY(); + } + return 0; + } + + public double getCenterDistance() { + double dx = px - cx.getEstimate(); + double dy = py - cy.getEstimate(); + return Math.sqrt(dx * dx + dy * dy); + } + + public double getPartialDiX() { + return (cx.getEstimate() - px) / getCenterDistance(); + } + + public double getPartialDiY() { + return (cy.getEstimate() - py) / getCenterDistance(); + } + + public double getTheoreticalValue() { + return getCenterDistance() - getRadius(); + } + + private double px; + private double py; + private static final long serialVersionUID = 1L; + + } + + private EstimatedParameter cx; + private EstimatedParameter cy; + private ArrayList points; + + } + + private static class QuadraticProblem extends SimpleEstimationProblem { + + private EstimatedParameter a; + private EstimatedParameter b; + private EstimatedParameter c; + + public QuadraticProblem() { + a = new EstimatedParameter("a", 0.0); + b = new EstimatedParameter("b", 0.0); + c = new EstimatedParameter("c", 0.0); + addParameter(a); + addParameter(b); + addParameter(c); + } + + public void addPoint(double x, double y, double w) { + addMeasurement(new LocalMeasurement(x, y, w)); + } + + public double getA() { + return a.getEstimate(); + } + + public double getB() { + return b.getEstimate(); + } + + public double getC() { + return c.getEstimate(); + } + + public double theoreticalValue(double x) { + return ( (a.getEstimate() * x + b.getEstimate() ) * x + c.getEstimate()); + } + + private double partial(double x, EstimatedParameter parameter) { + if (parameter == a) { + return x * x; + } else if (parameter == b) { + return x; + } else { + return 1.0; + } + } + + private class LocalMeasurement extends WeightedMeasurement { + + private static final long serialVersionUID = 1555043155023729130L; + private final double x; + + // constructor + public LocalMeasurement(double x, double y, double w) { + super(w, y); + this.x = x; + } + + public double getTheoreticalValue() { + return theoreticalValue(x); + } + + public double getPartial(EstimatedParameter parameter) { + return partial(x, parameter); + } + + } + } + + public static Test suite() { + return new TestSuite(LevenbergMarquardtEstimatorTest.class); + } + +} diff --git a/src/test/org/apache/commons/math/estimation/MinpackTest.java b/src/test/org/apache/commons/math/estimation/MinpackTest.java new file mode 100644 index 000000000..f17df0a94 --- /dev/null +++ b/src/test/org/apache/commons/math/estimation/MinpackTest.java @@ -0,0 +1,1518 @@ +/* + * 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.estimation; + +import java.util.Arrays; + +import org.apache.commons.math.estimation.EstimatedParameter; +import org.apache.commons.math.estimation.EstimationException; +import org.apache.commons.math.estimation.EstimationProblem; +import org.apache.commons.math.estimation.LevenbergMarquardtEstimator; +import org.apache.commons.math.estimation.WeightedMeasurement; + +import junit.framework.*; + +/** + *

        Some of the unit tests are re-implementations of the MINPACK file17 and file22 test files. + * The redistribution policy for MINPACK is available here, for + * convenience, it is reproduced below.

        + + * + * + * + *
        + * Minpack Copyright Notice (1999) University of Chicago. + * All rights reserved + *
        + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + *
          + *
        1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer.
          2. + *
        2. Redistributions in binary form must reproduce the above + * copyright notice, this list of conditions and the following + * disclaimer in the documentation and/or other materials provided + * with the distribution.
          3. + *
        3. The end-user documentation included with the redistribution, if any, + * must include the following acknowledgment: + * This product includes software developed by the University of + * Chicago, as Operator of Argonne National Laboratory. + * Alternately, this acknowledgment may appear in the software itself, + * if and wherever such third-party acknowledgments normally appear.
          4. + *
        4. WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS" + * WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE + * UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND + * THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR + * IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES + * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE + * OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY + * OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR + * USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF + * THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4) + * DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION + * UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL + * BE CORRECTED.
          5. + *
        5. LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT + * HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF + * ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT, + * INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF + * ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF + * PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER + * SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT + * (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE, + * EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE + * POSSIBILITY OF SUCH LOSS OR DAMAGES.
          6. + *
          + + * @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran minpack tests) + * @author Burton S. Garbow (original fortran minpack tests) + * @author Kenneth E. Hillstrom (original fortran minpack tests) + * @author Jorge J. More (original fortran minpack tests) + * @author Luc Maisonobe (non-minpack tests and minpack tests Java translation) + */ +@Deprecated +public class MinpackTest + extends TestCase { + + public MinpackTest(String name) { + super(name); + } + + public void testMinpackLinearFullRank() + throws EstimationException { + minpackTest(new LinearFullRankFunction(10, 5, 1.0, + 5.0, 2.23606797749979), false); + minpackTest(new LinearFullRankFunction(50, 5, 1.0, + 8.06225774829855, 6.70820393249937), false); + } + + public void testMinpackLinearRank1() + throws EstimationException { + minpackTest(new LinearRank1Function(10, 5, 1.0, + 291.521868819476, 1.4638501094228), false); + minpackTest(new LinearRank1Function(50, 5, 1.0, + 3101.60039334535, 3.48263016573496), false); + } + + public void testMinpackLinearRank1ZeroColsAndRows() + throws EstimationException { + minpackTest(new LinearRank1ZeroColsAndRowsFunction(10, 5, 1.0), false); + minpackTest(new LinearRank1ZeroColsAndRowsFunction(50, 5, 1.0), false); + } + + public void testMinpackRosenbrok() + throws EstimationException { + minpackTest(new RosenbrockFunction(new double[] { -1.2, 1.0 }, + Math.sqrt(24.2)), false); + minpackTest(new RosenbrockFunction(new double[] { -12.0, 10.0 }, + Math.sqrt(1795769.0)), false); + minpackTest(new RosenbrockFunction(new double[] { -120.0, 100.0 }, + 11.0 * Math.sqrt(169000121.0)), false); + } + + public void testMinpackHelicalValley() + throws EstimationException { + minpackTest(new HelicalValleyFunction(new double[] { -1.0, 0.0, 0.0 }, + 50.0), false); + minpackTest(new HelicalValleyFunction(new double[] { -10.0, 0.0, 0.0 }, + 102.95630140987), false); + minpackTest(new HelicalValleyFunction(new double[] { -100.0, 0.0, 0.0}, + 991.261822123701), false); + } + + public void testMinpackPowellSingular() + throws EstimationException { + minpackTest(new PowellSingularFunction(new double[] { 3.0, -1.0, 0.0, 1.0 }, + 14.6628782986152), false); + minpackTest(new PowellSingularFunction(new double[] { 30.0, -10.0, 0.0, 10.0 }, + 1270.9838708654), false); + minpackTest(new PowellSingularFunction(new double[] { 300.0, -100.0, 0.0, 100.0 }, + 126887.903284750), false); + } + + public void testMinpackFreudensteinRoth() + throws EstimationException { + minpackTest(new FreudensteinRothFunction(new double[] { 0.5, -2.0 }, + 20.0124960961895, 6.99887517584575, + new double[] { + 11.4124844654993, + -0.896827913731509 + }), false); + minpackTest(new FreudensteinRothFunction(new double[] { 5.0, -20.0 }, + 12432.833948863, 6.9988751744895, + new double[] { + 11.4130046614746, + -0.896796038685958 + }), false); + minpackTest(new FreudensteinRothFunction(new double[] { 50.0, -200.0 }, + 11426454.595762, 6.99887517242903, + new double[] { + 11.4127817857886, + -0.89680510749204 + }), false); + } + + public void testMinpackBard() + throws EstimationException { + minpackTest(new BardFunction(1.0, 6.45613629515967, 0.0906359603390466, + new double[] { + 0.0824105765758334, + 1.1330366534715, + 2.34369463894115 + }), false); + minpackTest(new BardFunction(10.0, 36.1418531596785, 4.17476870138539, + new double[] { + 0.840666673818329, + -158848033.259565, + -164378671.653535 + }), false); + minpackTest(new BardFunction(100.0, 384.114678637399, 4.17476870135969, + new double[] { + 0.840666673867645, + -158946167.205518, + -164464906.857771 + }), false); + } + + public void testMinpackKowalikOsborne() + throws EstimationException { + minpackTest(new KowalikOsborneFunction(new double[] { 0.25, 0.39, 0.415, 0.39 }, + 0.0728915102882945, + 0.017535837721129, + new double[] { + 0.192807810476249, + 0.191262653354071, + 0.123052801046931, + 0.136053221150517 + }), false); + minpackTest(new KowalikOsborneFunction(new double[] { 2.5, 3.9, 4.15, 3.9 }, + 2.97937007555202, + 0.032052192917937, + new double[] { + 728675.473768287, + -14.0758803129393, + -32977797.7841797, + -20571594.1977912 + }), false); + minpackTest(new KowalikOsborneFunction(new double[] { 25.0, 39.0, 41.5, 39.0 }, + 29.9590617016037, + 0.0175364017658228, + new double[] { + 0.192948328597594, + 0.188053165007911, + 0.122430604321144, + 0.134575665392506 + }), true); + } + + public void testMinpackMeyer() + throws EstimationException { + minpackTest(new MeyerFunction(new double[] { 0.02, 4000.0, 250.0 }, + 41153.4665543031, 9.37794514651874, + new double[] { + 0.00560963647102661, + 6181.34634628659, + 345.223634624144 + }), false); + minpackTest(new MeyerFunction(new double[] { 0.2, 40000.0, 2500.0 }, + 4168216.89130846, 792.917871779501, + new double[] { + 1.42367074157994e-11, + 33695.7133432541, + 901.268527953801 + }), true); + } + + public void testMinpackWatson() + throws EstimationException { + + minpackTest(new WatsonFunction(6, 0.0, + 5.47722557505166, 0.0478295939097601, + new double[] { + -0.0157249615083782, 1.01243488232965, + -0.232991722387673, 1.26043101102818, + -1.51373031394421, 0.99299727291842 + }), false); + minpackTest(new WatsonFunction(6, 10.0, + 6433.12578950026, 0.0478295939096951, + new double[] { + -0.0157251901386677, 1.01243485860105, + -0.232991545843829, 1.26042932089163, + -1.51372776706575, 0.99299573426328 + }), false); + minpackTest(new WatsonFunction(6, 100.0, + 674256.040605213, 0.047829593911544, + new double[] { + -0.0157247019712586, 1.01243490925658, + -0.232991922761641, 1.26043292929555, + -1.51373320452707, 0.99299901922322 + }), false); + + minpackTest(new WatsonFunction(9, 0.0, + 5.47722557505166, 0.00118311459212420, + new double[] { + -0.153070644166722e-4, 0.999789703934597, + 0.0147639634910978, 0.146342330145992, + 1.00082109454817, -2.61773112070507, + 4.10440313943354, -3.14361226236241, + 1.05262640378759 + }), false); + minpackTest(new WatsonFunction(9, 10.0, + 12088.127069307, 0.00118311459212513, + new double[] { + -0.153071334849279e-4, 0.999789703941234, + 0.0147639629786217, 0.146342334818836, + 1.00082107321386, -2.61773107084722, + 4.10440307655564, -3.14361222178686, + 1.05262639322589 + }), false); + minpackTest(new WatsonFunction(9, 100.0, + 1269109.29043834, 0.00118311459212384, + new double[] { + -0.153069523352176e-4, 0.999789703958371, + 0.0147639625185392, 0.146342341096326, + 1.00082104729164, -2.61773101573645, + 4.10440301427286, -3.14361218602503, + 1.05262638516774 + }), false); + + minpackTest(new WatsonFunction(12, 0.0, + 5.47722557505166, 0.217310402535861e-4, + new double[] { + -0.660266001396382e-8, 1.00000164411833, + -0.000563932146980154, 0.347820540050756, + -0.156731500244233, 1.05281515825593, + -3.24727109519451, 7.2884347837505, + -10.271848098614, 9.07411353715783, + -4.54137541918194, 1.01201187975044 + }), false); + minpackTest(new WatsonFunction(12, 10.0, + 19220.7589790951, 0.217310402518509e-4, + new double[] { + -0.663710223017410e-8, 1.00000164411787, + -0.000563932208347327, 0.347820540486998, + -0.156731503955652, 1.05281517654573, + -3.2472711515214, 7.28843489430665, + -10.2718482369638, 9.07411364383733, + -4.54137546533666, 1.01201188830857 + }), false); + minpackTest(new WatsonFunction(12, 100.0, + 2018918.04462367, 0.217310402539845e-4, + new double[] { + -0.663806046485249e-8, 1.00000164411786, + -0.000563932210324959, 0.347820540503588, + -0.156731504091375, 1.05281517718031, + -3.24727115337025, 7.28843489775302, + -10.2718482410813, 9.07411364688464, + -4.54137546660822, 1.0120118885369 + }), false); + + } + + public void testMinpackBox3Dimensional() + throws EstimationException { + minpackTest(new Box3DimensionalFunction(10, new double[] { 0.0, 10.0, 20.0 }, + 32.1115837449572), false); + } + + public void testMinpackJennrichSampson() + throws EstimationException { + minpackTest(new JennrichSampsonFunction(10, new double[] { 0.3, 0.4 }, + 64.5856498144943, 11.1517793413499, + new double[] { + 0.257819926636811, 0.257829976764542 + }), false); + } + + public void testMinpackBrownDennis() + throws EstimationException { + minpackTest(new BrownDennisFunction(20, + new double[] { 25.0, 5.0, -5.0, -1.0 }, + 2815.43839161816, 292.954288244866, + new double[] { + -11.59125141003, 13.2024883984741, + -0.403574643314272, 0.236736269844604 + }), false); + minpackTest(new BrownDennisFunction(20, + new double[] { 250.0, 50.0, -50.0, -10.0 }, + 555073.354173069, 292.954270581415, + new double[] { + -11.5959274272203, 13.2041866926242, + -0.403417362841545, 0.236771143410386 + }), false); + minpackTest(new BrownDennisFunction(20, + new double[] { 2500.0, 500.0, -500.0, -100.0 }, + 61211252.2338581, 292.954306151134, + new double[] { + -11.5902596937374, 13.2020628854665, + -0.403688070279258, 0.236665033746463 + }), false); + } + + public void testMinpackChebyquad() + throws EstimationException { + minpackTest(new ChebyquadFunction(1, 8, 1.0, + 1.88623796907732, 1.88623796907732, + new double[] { 0.5 }), false); + minpackTest(new ChebyquadFunction(1, 8, 10.0, + 5383344372.34005, 1.88424820499951, + new double[] { 0.9817314924684 }), false); + minpackTest(new ChebyquadFunction(1, 8, 100.0, + 0.118088726698392e19, 1.88424820499347, + new double[] { 0.9817314852934 }), false); + minpackTest(new ChebyquadFunction(8, 8, 1.0, + 0.196513862833975, 0.0593032355046727, + new double[] { + 0.0431536648587336, 0.193091637843267, + 0.266328593812698, 0.499999334628884, + 0.500000665371116, 0.733671406187302, + 0.806908362156733, 0.956846335141266 + }), false); + minpackTest(new ChebyquadFunction(9, 9, 1.0, + 0.16994993465202, 0.0, + new double[] { + 0.0442053461357828, 0.199490672309881, + 0.23561910847106, 0.416046907892598, + 0.5, 0.583953092107402, + 0.764380891528940, 0.800509327690119, + 0.955794653864217 + }), false); + minpackTest(new ChebyquadFunction(10, 10, 1.0, + 0.183747831178711, 0.0806471004038253, + new double[] { + 0.0596202671753563, 0.166708783805937, + 0.239171018813509, 0.398885290346268, + 0.398883667870681, 0.601116332129320, + 0.60111470965373, 0.760828981186491, + 0.833291216194063, 0.940379732824644 + }), false); + } + + public void testMinpackBrownAlmostLinear() + throws EstimationException { + minpackTest(new BrownAlmostLinearFunction(10, 0.5, + 16.5302162063499, 0.0, + new double[] { + 0.979430303349862, 0.979430303349862, + 0.979430303349862, 0.979430303349862, + 0.979430303349862, 0.979430303349862, + 0.979430303349862, 0.979430303349862, + 0.979430303349862, 1.20569696650138 + }), false); + minpackTest(new BrownAlmostLinearFunction(10, 5.0, + 9765624.00089211, 0.0, + new double[] { + 0.979430303349865, 0.979430303349865, + 0.979430303349865, 0.979430303349865, + 0.979430303349865, 0.979430303349865, + 0.979430303349865, 0.979430303349865, + 0.979430303349865, 1.20569696650135 + }), false); + minpackTest(new BrownAlmostLinearFunction(10, 50.0, + 0.9765625e17, 0.0, + new double[] { + 1.0, 1.0, 1.0, 1.0, 1.0, + 1.0, 1.0, 1.0, 1.0, 1.0 + }), false); + minpackTest(new BrownAlmostLinearFunction(30, 0.5, + 83.476044467848, 0.0, + new double[] { + 0.997754216442807, 0.997754216442807, + 0.997754216442807, 0.997754216442807, + 0.997754216442807, 0.997754216442807, + 0.997754216442807, 0.997754216442807, + 0.997754216442807, 0.997754216442807, + 0.997754216442807, 0.997754216442807, + 0.997754216442807, 0.997754216442807, + 0.997754216442807, 0.997754216442807, + 0.997754216442807, 0.997754216442807, + 0.997754216442807, 0.997754216442807, + 0.997754216442807, 0.997754216442807, + 0.997754216442807, 0.997754216442807, + 0.997754216442807, 0.997754216442807, + 0.997754216442807, 0.997754216442807, + 0.997754216442807, 1.06737350671578 + }), false); + minpackTest(new BrownAlmostLinearFunction(40, 0.5, + 128.026364472323, 0.0, + new double[] { + 1.00000000000002, 1.00000000000002, + 1.00000000000002, 1.00000000000002, + 1.00000000000002, 1.00000000000002, + 1.00000000000002, 1.00000000000002, + 1.00000000000002, 1.00000000000002, + 1.00000000000002, 1.00000000000002, + 1.00000000000002, 1.00000000000002, + 1.00000000000002, 1.00000000000002, + 1.00000000000002, 1.00000000000002, + 1.00000000000002, 1.00000000000002, + 1.00000000000002, 1.00000000000002, + 1.00000000000002, 1.00000000000002, + 1.00000000000002, 1.00000000000002, + 1.00000000000002, 1.00000000000002, + 1.00000000000002, 1.00000000000002, + 1.00000000000002, 1.00000000000002, + 1.00000000000002, 1.00000000000002, + 0.999999999999121 + }), false); + } + + public void testMinpackOsborne1() + throws EstimationException { + minpackTest(new Osborne1Function(new double[] { 0.5, 1.5, -1.0, 0.01, 0.02, }, + 0.937564021037838, 0.00739249260904843, + new double[] { + 0.375410049244025, 1.93584654543108, + -1.46468676748716, 0.0128675339110439, + 0.0221227011813076 + }), false); + } + + public void testMinpackOsborne2() + throws EstimationException { + + minpackTest(new Osborne2Function(new double[] { + 1.3, 0.65, 0.65, 0.7, 0.6, + 3.0, 5.0, 7.0, 2.0, 4.5, 5.5 + }, + 1.44686540984712, 0.20034404483314, + new double[] { + 1.30997663810096, 0.43155248076, + 0.633661261602859, 0.599428560991695, + 0.754179768272449, 0.904300082378518, + 1.36579949521007, 4.82373199748107, + 2.39868475104871, 4.56887554791452, + 5.67534206273052 + }), false); + } + + private void minpackTest(MinpackFunction function, boolean exceptionExpected) { + LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); + estimator.setMaxCostEval(100 * (function.getN() + 1)); + estimator.setCostRelativeTolerance(Math.sqrt(2.22044604926e-16)); + estimator.setParRelativeTolerance(Math.sqrt(2.22044604926e-16)); + estimator.setOrthoTolerance(2.22044604926e-16); + assertTrue(function.checkTheoreticalStartCost(estimator.getRMS(function))); + try { + estimator.estimate(function); + assertFalse(exceptionExpected); + } catch (EstimationException lsse) { + assertTrue(exceptionExpected); + } + assertTrue(function.checkTheoreticalMinCost(estimator.getRMS(function))); + assertTrue(function.checkTheoreticalMinParams()); + } + + private static abstract class MinpackFunction implements EstimationProblem { + + protected MinpackFunction(int m, + double[] startParams, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + this.m = m; + this.n = startParams.length; + parameters = new EstimatedParameter[n]; + for (int i = 0; i < n; ++i) { + parameters[i] = new EstimatedParameter("p" + i, startParams[i]); + } + this.theoreticalStartCost = theoreticalStartCost; + this.theoreticalMinCost = theoreticalMinCost; + this.theoreticalMinParams = theoreticalMinParams; + this.costAccuracy = 1.0e-8; + this.paramsAccuracy = 1.0e-5; + } + + protected static double[] buildArray(int n, double x) { + double[] array = new double[n]; + Arrays.fill(array, x); + return array; + } + + protected void setCostAccuracy(double costAccuracy) { + this.costAccuracy = costAccuracy; + } + + protected void setParamsAccuracy(double paramsAccuracy) { + this.paramsAccuracy = paramsAccuracy; + } + + public int getN() { + return parameters.length; + } + + public boolean checkTheoreticalStartCost(double rms) { + double threshold = costAccuracy * (1.0 + theoreticalStartCost); + return Math.abs(Math.sqrt(m) * rms - theoreticalStartCost) <= threshold; + } + + public boolean checkTheoreticalMinCost(double rms) { + double threshold = costAccuracy * (1.0 + theoreticalMinCost); + return Math.abs(Math.sqrt(m) * rms - theoreticalMinCost) <= threshold; + } + + public boolean checkTheoreticalMinParams() { + if (theoreticalMinParams != null) { + for (int i = 0; i < theoreticalMinParams.length; ++i) { + double mi = theoreticalMinParams[i]; + double vi = parameters[i].getEstimate(); + if (Math.abs(mi - vi) > (paramsAccuracy * (1.0 + Math.abs(mi)))) { + return false; + } + } + } + return true; + } + + public WeightedMeasurement[] getMeasurements() { + WeightedMeasurement[] measurements = new WeightedMeasurement[m]; + for (int i = 0; i < m; ++i) { + measurements[i] = new MinpackMeasurement(i); + } + return measurements; + } + + public EstimatedParameter[] getUnboundParameters() { + return parameters; + } + + public EstimatedParameter[] getAllParameters() { + return parameters; + } + + protected abstract double[][] getJacobian(); + + protected abstract double[] getResiduals(); + + private class MinpackMeasurement extends WeightedMeasurement { + + public MinpackMeasurement(int index) { + super(1.0, 0.0); + this.index = index; + } + + public double getTheoreticalValue() { + // this is obviously NOT efficient as we recompute the whole vector + // each time we need only one element, but it is only for test + // purposes and is simpler to check. + // This implementation should NOT be taken as an example, it is ugly! + return getResiduals()[index]; + } + + public double getPartial(EstimatedParameter parameter) { + // this is obviously NOT efficient as we recompute the whole jacobian + // each time we need only one element, but it is only for test + // purposes and is simpler to check. + // This implementation should NOT be taken as an example, it is ugly! + for (int j = 0; j < n; ++j) { + if (parameter == parameters[j]) { + return getJacobian()[index][j]; + } + } + return 0; + } + + private int index; + private static final long serialVersionUID = 1L; + + } + + protected int n; + protected int m; + protected EstimatedParameter[] parameters; + protected double theoreticalStartCost; + protected double theoreticalMinCost; + protected double[] theoreticalMinParams; + protected double costAccuracy; + protected double paramsAccuracy; + + } + + private static class LinearFullRankFunction extends MinpackFunction { + + public LinearFullRankFunction(int m, int n, double x0, + double theoreticalStartCost, + double theoreticalMinCost) { + super(m, buildArray(n, x0), theoreticalStartCost, + theoreticalMinCost, buildArray(n, -1.0)); + } + + protected double[][] getJacobian() { + double t = 2.0 / m; + double[][] jacobian = new double[m][]; + for (int i = 0; i < m; ++i) { + jacobian[i] = new double[n]; + for (int j = 0; j < n; ++j) { + jacobian[i][j] = (i == j) ? (1 - t) : -t; + } + } + return jacobian; + } + + protected double[] getResiduals() { + double sum = 0; + for (int i = 0; i < n; ++i) { + sum += parameters[i].getEstimate(); + } + double t = 1 + 2 * sum / m; + double[] f = new double[m]; + for (int i = 0; i < n; ++i) { + f[i] = parameters[i].getEstimate() - t; + } + Arrays.fill(f, n, m, -t); + return f; + } + + } + + private static class LinearRank1Function extends MinpackFunction { + + public LinearRank1Function(int m, int n, double x0, + double theoreticalStartCost, + double theoreticalMinCost) { + super(m, buildArray(n, x0), theoreticalStartCost, theoreticalMinCost, null); + } + + protected double[][] getJacobian() { + double[][] jacobian = new double[m][]; + for (int i = 0; i < m; ++i) { + jacobian[i] = new double[n]; + for (int j = 0; j < n; ++j) { + jacobian[i][j] = (i + 1) * (j + 1); + } + } + return jacobian; + } + + protected double[] getResiduals() { + double[] f = new double[m]; + double sum = 0; + for (int i = 0; i < n; ++i) { + sum += (i + 1) * parameters[i].getEstimate(); + } + for (int i = 0; i < m; ++i) { + f[i] = (i + 1) * sum - 1; + } + return f; + } + + } + + private static class LinearRank1ZeroColsAndRowsFunction extends MinpackFunction { + + public LinearRank1ZeroColsAndRowsFunction(int m, int n, double x0) { + super(m, buildArray(n, x0), + Math.sqrt(m + (n+1)*(n-2)*(m-2)*(m-1) * ((n+1)*(n-2)*(2*m-3) - 12) / 24.0), + Math.sqrt((m * (m + 3) - 6) / (2.0 * (2 * m - 3))), + null); + } + + protected double[][] getJacobian() { + double[][] jacobian = new double[m][]; + for (int i = 0; i < m; ++i) { + jacobian[i] = new double[n]; + jacobian[i][0] = 0; + for (int j = 1; j < (n - 1); ++j) { + if (i == 0) { + jacobian[i][j] = 0; + } else if (i != (m - 1)) { + jacobian[i][j] = i * (j + 1); + } else { + jacobian[i][j] = 0; + } + } + jacobian[i][n - 1] = 0; + } + return jacobian; + } + + protected double[] getResiduals() { + double[] f = new double[m]; + double sum = 0; + for (int i = 1; i < (n - 1); ++i) { + sum += (i + 1) * parameters[i].getEstimate(); + } + for (int i = 0; i < (m - 1); ++i) { + f[i] = i * sum - 1; + } + f[m - 1] = -1; + return f; + } + + } + + private static class RosenbrockFunction extends MinpackFunction { + + public RosenbrockFunction(double[] startParams, double theoreticalStartCost) { + super(2, startParams, theoreticalStartCost, 0.0, buildArray(2, 1.0)); + } + + protected double[][] getJacobian() { + double x1 = parameters[0].getEstimate(); + return new double[][] { { -20 * x1, 10 }, { -1, 0 } }; + } + + protected double[] getResiduals() { + double x1 = parameters[0].getEstimate(); + double x2 = parameters[1].getEstimate(); + return new double[] { 10 * (x2 - x1 * x1), 1 - x1 }; + } + + } + + private static class HelicalValleyFunction extends MinpackFunction { + + public HelicalValleyFunction(double[] startParams, + double theoreticalStartCost) { + super(3, startParams, theoreticalStartCost, 0.0, + new double[] { 1.0, 0.0, 0.0 }); + } + + protected double[][] getJacobian() { + double x1 = parameters[0].getEstimate(); + double x2 = parameters[1].getEstimate(); + double tmpSquare = x1 * x1 + x2 * x2; + double tmp1 = twoPi * tmpSquare; + double tmp2 = Math.sqrt(tmpSquare); + return new double[][] { + { 100 * x2 / tmp1, -100 * x1 / tmp1, 10 }, + { 10 * x1 / tmp2, 10 * x2 / tmp2, 0 }, + { 0, 0, 1 } + }; + } + + protected double[] getResiduals() { + double x1 = parameters[0].getEstimate(); + double x2 = parameters[1].getEstimate(); + double x3 = parameters[2].getEstimate(); + double tmp1; + if (x1 == 0) { + tmp1 = (x2 >= 0) ? 0.25 : -0.25; + } else { + tmp1 = Math.atan(x2 / x1) / twoPi; + if (x1 < 0) { + tmp1 += 0.5; + } + } + double tmp2 = Math.sqrt(x1 * x1 + x2 * x2); + return new double[] { + 10.0 * (x3 - 10 * tmp1), + 10.0 * (tmp2 - 1), + x3 + }; + } + + private static final double twoPi = 2.0 * Math.PI; + + } + + private static class PowellSingularFunction extends MinpackFunction { + + public PowellSingularFunction(double[] startParams, + double theoreticalStartCost) { + super(4, startParams, theoreticalStartCost, 0.0, buildArray(4, 0.0)); + } + + protected double[][] getJacobian() { + double x1 = parameters[0].getEstimate(); + double x2 = parameters[1].getEstimate(); + double x3 = parameters[2].getEstimate(); + double x4 = parameters[3].getEstimate(); + return new double[][] { + { 1, 10, 0, 0 }, + { 0, 0, sqrt5, -sqrt5 }, + { 0, 2 * (x2 - 2 * x3), -4 * (x2 - 2 * x3), 0 }, + { 2 * sqrt10 * (x1 - x4), 0, 0, -2 * sqrt10 * (x1 - x4) } + }; + } + + protected double[] getResiduals() { + double x1 = parameters[0].getEstimate(); + double x2 = parameters[1].getEstimate(); + double x3 = parameters[2].getEstimate(); + double x4 = parameters[3].getEstimate(); + return new double[] { + x1 + 10 * x2, + sqrt5 * (x3 - x4), + (x2 - 2 * x3) * (x2 - 2 * x3), + sqrt10 * (x1 - x4) * (x1 - x4) + }; + } + + private static final double sqrt5 = Math.sqrt( 5.0); + private static final double sqrt10 = Math.sqrt(10.0); + + } + + private static class FreudensteinRothFunction extends MinpackFunction { + + public FreudensteinRothFunction(double[] startParams, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(2, startParams, theoreticalStartCost, + theoreticalMinCost, theoreticalMinParams); + } + + protected double[][] getJacobian() { + double x2 = parameters[1].getEstimate(); + return new double[][] { + { 1, x2 * (10 - 3 * x2) - 2 }, + { 1, x2 * ( 2 + 3 * x2) - 14, } + }; + } + + protected double[] getResiduals() { + double x1 = parameters[0].getEstimate(); + double x2 = parameters[1].getEstimate(); + return new double[] { + -13.0 + x1 + ((5.0 - x2) * x2 - 2.0) * x2, + -29.0 + x1 + ((1.0 + x2) * x2 - 14.0) * x2 + }; + } + + } + + private static class BardFunction extends MinpackFunction { + + public BardFunction(double x0, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(15, buildArray(3, x0), theoreticalStartCost, + theoreticalMinCost, theoreticalMinParams); + } + + protected double[][] getJacobian() { + double x2 = parameters[1].getEstimate(); + double x3 = parameters[2].getEstimate(); + double[][] jacobian = new double[m][]; + for (int i = 0; i < m; ++i) { + double tmp1 = i + 1; + double tmp2 = 15 - i; + double tmp3 = (i <= 7) ? tmp1 : tmp2; + double tmp4 = x2 * tmp2 + x3 * tmp3; + tmp4 *= tmp4; + jacobian[i] = new double[] { -1, tmp1 * tmp2 / tmp4, tmp1 * tmp3 / tmp4 }; + } + return jacobian; + } + + protected double[] getResiduals() { + double x1 = parameters[0].getEstimate(); + double x2 = parameters[1].getEstimate(); + double x3 = parameters[2].getEstimate(); + double[] f = new double[m]; + for (int i = 0; i < m; ++i) { + double tmp1 = i + 1; + double tmp2 = 15 - i; + double tmp3 = (i <= 7) ? tmp1 : tmp2; + f[i] = y[i] - (x1 + tmp1 / (x2 * tmp2 + x3 * tmp3)); + } + return f; + } + + private static final double[] y = { + 0.14, 0.18, 0.22, 0.25, 0.29, + 0.32, 0.35, 0.39, 0.37, 0.58, + 0.73, 0.96, 1.34, 2.10, 4.39 + }; + + } + + private static class KowalikOsborneFunction extends MinpackFunction { + + public KowalikOsborneFunction(double[] startParams, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(11, startParams, theoreticalStartCost, + theoreticalMinCost, theoreticalMinParams); + if (theoreticalStartCost > 20.0) { + setCostAccuracy(2.0e-4); + setParamsAccuracy(5.0e-3); + } + } + + protected double[][] getJacobian() { + double x1 = parameters[0].getEstimate(); + double x2 = parameters[1].getEstimate(); + double x3 = parameters[2].getEstimate(); + double x4 = parameters[3].getEstimate(); + double[][] jacobian = new double[m][]; + for (int i = 0; i < m; ++i) { + double tmp = v[i] * (v[i] + x3) + x4; + double j1 = -v[i] * (v[i] + x2) / tmp; + double j2 = -v[i] * x1 / tmp; + double j3 = j1 * j2; + double j4 = j3 / v[i]; + jacobian[i] = new double[] { j1, j2, j3, j4 }; + } + return jacobian; + } + + protected double[] getResiduals() { + double x1 = parameters[0].getEstimate(); + double x2 = parameters[1].getEstimate(); + double x3 = parameters[2].getEstimate(); + double x4 = parameters[3].getEstimate(); + double[] f = new double[m]; + for (int i = 0; i < m; ++i) { + f[i] = y[i] - x1 * (v[i] * (v[i] + x2)) / (v[i] * (v[i] + x3) + x4); + } + return f; + } + + private static final double[] v = { + 4.0, 2.0, 1.0, 0.5, 0.25, 0.167, 0.125, 0.1, 0.0833, 0.0714, 0.0625 + }; + + private static final double[] y = { + 0.1957, 0.1947, 0.1735, 0.1600, 0.0844, 0.0627, + 0.0456, 0.0342, 0.0323, 0.0235, 0.0246 + }; + + } + + private static class MeyerFunction extends MinpackFunction { + + public MeyerFunction(double[] startParams, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(16, startParams, theoreticalStartCost, + theoreticalMinCost, theoreticalMinParams); + if (theoreticalStartCost > 1.0e6) { + setCostAccuracy(7.0e-3); + setParamsAccuracy(2.0e-2); + } + } + + protected double[][] getJacobian() { + double x1 = parameters[0].getEstimate(); + double x2 = parameters[1].getEstimate(); + double x3 = parameters[2].getEstimate(); + double[][] jacobian = new double[m][]; + for (int i = 0; i < m; ++i) { + double temp = 5.0 * (i + 1) + 45.0 + x3; + double tmp1 = x2 / temp; + double tmp2 = Math.exp(tmp1); + double tmp3 = x1 * tmp2 / temp; + jacobian[i] = new double[] { tmp2, tmp3, -tmp1 * tmp3 }; + } + return jacobian; + } + + protected double[] getResiduals() { + double x1 = parameters[0].getEstimate(); + double x2 = parameters[1].getEstimate(); + double x3 = parameters[2].getEstimate(); + double[] f = new double[m]; + for (int i = 0; i < m; ++i) { + f[i] = x1 * Math.exp(x2 / (5.0 * (i + 1) + 45.0 + x3)) - y[i]; + } + return f; + } + + private static final double[] y = { + 34780.0, 28610.0, 23650.0, 19630.0, + 16370.0, 13720.0, 11540.0, 9744.0, + 8261.0, 7030.0, 6005.0, 5147.0, + 4427.0, 3820.0, 3307.0, 2872.0 + }; + + } + + private static class WatsonFunction extends MinpackFunction { + + public WatsonFunction(int n, double x0, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(31, buildArray(n, x0), theoreticalStartCost, + theoreticalMinCost, theoreticalMinParams); + } + + protected double[][] getJacobian() { + + double[][] jacobian = new double[m][]; + + for (int i = 0; i < (m - 2); ++i) { + double div = (i + 1) / 29.0; + double s2 = 0.0; + double dx = 1.0; + for (int j = 0; j < n; ++j) { + s2 += dx * parameters[j].getEstimate(); + dx *= div; + } + double temp= 2 * div * s2; + dx = 1.0 / div; + jacobian[i] = new double[n]; + for (int j = 0; j < n; ++j) { + jacobian[i][j] = dx * (j - temp); + dx *= div; + } + } + + jacobian[m - 2] = new double[n]; + jacobian[m - 2][0] = 1; + + jacobian[m - 1] = new double[n]; + jacobian[m - 1][0]= -2 * parameters[0].getEstimate(); + jacobian[m - 1][1]= 1; + + return jacobian; + + } + + protected double[] getResiduals() { + double[] f = new double[m]; + for (int i = 0; i < (m - 2); ++i) { + double div = (i + 1) / 29.0; + double s1 = 0; + double dx = 1; + for (int j = 1; j < n; ++j) { + s1 += j * dx * parameters[j].getEstimate(); + dx *= div; + } + double s2 =0; + dx =1; + for (int j = 0; j < n; ++j) { + s2 += dx * parameters[j].getEstimate(); + dx *= div; + } + f[i] = s1 - s2 * s2 - 1; + } + + double x1 = parameters[0].getEstimate(); + double x2 = parameters[1].getEstimate(); + f[m - 2] = x1; + f[m - 1] = x2 - x1 * x1 - 1; + + return f; + + } + + } + + private static class Box3DimensionalFunction extends MinpackFunction { + + public Box3DimensionalFunction(int m, double[] startParams, + double theoreticalStartCost) { + super(m, startParams, theoreticalStartCost, + 0.0, new double[] { 1.0, 10.0, 1.0 }); + } + + protected double[][] getJacobian() { + double x1 = parameters[0].getEstimate(); + double x2 = parameters[1].getEstimate(); + double[][] jacobian = new double[m][]; + for (int i = 0; i < m; ++i) { + double tmp = (i + 1) / 10.0; + jacobian[i] = new double[] { + -tmp * Math.exp(-tmp * x1), + tmp * Math.exp(-tmp * x2), + Math.exp(-i - 1) - Math.exp(-tmp) + }; + } + return jacobian; + } + + protected double[] getResiduals() { + double x1 = parameters[0].getEstimate(); + double x2 = parameters[1].getEstimate(); + double x3 = parameters[2].getEstimate(); + double[] f = new double[m]; + for (int i = 0; i < m; ++i) { + double tmp = (i + 1) / 10.0; + f[i] = Math.exp(-tmp * x1) - Math.exp(-tmp * x2) + + (Math.exp(-i - 1) - Math.exp(-tmp)) * x3; + } + return f; + } + + } + + private static class JennrichSampsonFunction extends MinpackFunction { + + public JennrichSampsonFunction(int m, double[] startParams, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(m, startParams, theoreticalStartCost, + theoreticalMinCost, theoreticalMinParams); + } + + protected double[][] getJacobian() { + double x1 = parameters[0].getEstimate(); + double x2 = parameters[1].getEstimate(); + double[][] jacobian = new double[m][]; + for (int i = 0; i < m; ++i) { + double t = i + 1; + jacobian[i] = new double[] { -t * Math.exp(t * x1), -t * Math.exp(t * x2) }; + } + return jacobian; + } + + protected double[] getResiduals() { + double x1 = parameters[0].getEstimate(); + double x2 = parameters[1].getEstimate(); + double[] f = new double[m]; + for (int i = 0; i < m; ++i) { + double temp = i + 1; + f[i] = 2 + 2 * temp - Math.exp(temp * x1) - Math.exp(temp * x2); + } + return f; + } + + } + + private static class BrownDennisFunction extends MinpackFunction { + + public BrownDennisFunction(int m, double[] startParams, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(m, startParams, theoreticalStartCost, + theoreticalMinCost, theoreticalMinParams); + } + + protected double[][] getJacobian() { + double x1 = parameters[0].getEstimate(); + double x2 = parameters[1].getEstimate(); + double x3 = parameters[2].getEstimate(); + double x4 = parameters[3].getEstimate(); + double[][] jacobian = new double[m][]; + for (int i = 0; i < m; ++i) { + double temp = (i + 1) / 5.0; + double ti = Math.sin(temp); + double tmp1 = x1 + temp * x2 - Math.exp(temp); + double tmp2 = x3 + ti * x4 - Math.cos(temp); + jacobian[i] = new double[] { + 2 * tmp1, 2 * temp * tmp1, 2 * tmp2, 2 * ti * tmp2 + }; + } + return jacobian; + } + + protected double[] getResiduals() { + double x1 = parameters[0].getEstimate(); + double x2 = parameters[1].getEstimate(); + double x3 = parameters[2].getEstimate(); + double x4 = parameters[3].getEstimate(); + double[] f = new double[m]; + for (int i = 0; i < m; ++i) { + double temp = (i + 1) / 5.0; + double tmp1 = x1 + temp * x2 - Math.exp(temp); + double tmp2 = x3 + Math.sin(temp) * x4 - Math.cos(temp); + f[i] = tmp1 * tmp1 + tmp2 * tmp2; + } + return f; + } + + } + + private static class ChebyquadFunction extends MinpackFunction { + + private static double[] buildChebyquadArray(int n, double factor) { + double[] array = new double[n]; + double inv = factor / (n + 1); + for (int i = 0; i < n; ++i) { + array[i] = (i + 1) * inv; + } + return array; + } + + public ChebyquadFunction(int n, int m, double factor, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(m, buildChebyquadArray(n, factor), theoreticalStartCost, + theoreticalMinCost, theoreticalMinParams); + } + + protected double[][] getJacobian() { + + double[][] jacobian = new double[m][]; + for (int i = 0; i < m; ++i) { + jacobian[i] = new double[n]; + } + + double dx = 1.0 / n; + for (int j = 0; j < n; ++j) { + double tmp1 = 1; + double tmp2 = 2 * parameters[j].getEstimate() - 1; + double temp = 2 * tmp2; + double tmp3 = 0; + double tmp4 = 2; + for (int i = 0; i < m; ++i) { + jacobian[i][j] = dx * tmp4; + double ti = 4 * tmp2 + temp * tmp4 - tmp3; + tmp3 = tmp4; + tmp4 = ti; + ti = temp * tmp2 - tmp1; + tmp1 = tmp2; + tmp2 = ti; + } + } + + return jacobian; + + } + + protected double[] getResiduals() { + + double[] f = new double[m]; + + for (int j = 0; j < n; ++j) { + double tmp1 = 1; + double tmp2 = 2 * parameters[j].getEstimate() - 1; + double temp = 2 * tmp2; + for (int i = 0; i < m; ++i) { + f[i] += tmp2; + double ti = temp * tmp2 - tmp1; + tmp1 = tmp2; + tmp2 = ti; + } + } + + double dx = 1.0 / n; + boolean iev = false; + for (int i = 0; i < m; ++i) { + f[i] *= dx; + if (iev) { + f[i] += 1.0 / (i * (i + 2)); + } + iev = ! iev; + } + + return f; + + } + + } + + private static class BrownAlmostLinearFunction extends MinpackFunction { + + public BrownAlmostLinearFunction(int m, double factor, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(m, buildArray(m, factor), theoreticalStartCost, + theoreticalMinCost, theoreticalMinParams); + } + + protected double[][] getJacobian() { + double[][] jacobian = new double[m][]; + for (int i = 0; i < m; ++i) { + jacobian[i] = new double[n]; + } + + double prod = 1; + for (int j = 0; j < n; ++j) { + prod *= parameters[j].getEstimate(); + for (int i = 0; i < n; ++i) { + jacobian[i][j] = 1; + } + jacobian[j][j] = 2; + } + + for (int j = 0; j < n; ++j) { + EstimatedParameter vj = parameters[j]; + double temp = vj.getEstimate(); + if (temp == 0) { + temp = 1; + prod = 1; + for (int k = 0; k < n; ++k) { + if (k != j) { + prod *= parameters[k].getEstimate(); + } + } + } + jacobian[n - 1][j] = prod / temp; + } + + return jacobian; + + } + + protected double[] getResiduals() { + double[] f = new double[m]; + double sum = -(n + 1); + double prod = 1; + for (int j = 0; j < n; ++j) { + sum += parameters[j].getEstimate(); + prod *= parameters[j].getEstimate(); + } + for (int i = 0; i < n; ++i) { + f[i] = parameters[i].getEstimate() + sum; + } + f[n - 1] = prod - 1; + return f; + } + + } + + private static class Osborne1Function extends MinpackFunction { + + public Osborne1Function(double[] startParams, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(33, startParams, theoreticalStartCost, + theoreticalMinCost, theoreticalMinParams); + } + + protected double[][] getJacobian() { + double x2 = parameters[1].getEstimate(); + double x3 = parameters[2].getEstimate(); + double x4 = parameters[3].getEstimate(); + double x5 = parameters[4].getEstimate(); + double[][] jacobian = new double[m][]; + for (int i = 0; i < m; ++i) { + double temp = 10.0 * i; + double tmp1 = Math.exp(-temp * x4); + double tmp2 = Math.exp(-temp * x5); + jacobian[i] = new double[] { + -1, -tmp1, -tmp2, temp * x2 * tmp1, temp * x3 * tmp2 + }; + } + return jacobian; + } + + protected double[] getResiduals() { + double x1 = parameters[0].getEstimate(); + double x2 = parameters[1].getEstimate(); + double x3 = parameters[2].getEstimate(); + double x4 = parameters[3].getEstimate(); + double x5 = parameters[4].getEstimate(); + double[] f = new double[m]; + for (int i = 0; i < m; ++i) { + double temp = 10.0 * i; + double tmp1 = Math.exp(-temp * x4); + double tmp2 = Math.exp(-temp * x5); + f[i] = y[i] - (x1 + x2 * tmp1 + x3 * tmp2); + } + return f; + } + + private static final double[] y = { + 0.844, 0.908, 0.932, 0.936, 0.925, 0.908, 0.881, 0.850, 0.818, 0.784, 0.751, + 0.718, 0.685, 0.658, 0.628, 0.603, 0.580, 0.558, 0.538, 0.522, 0.506, 0.490, + 0.478, 0.467, 0.457, 0.448, 0.438, 0.431, 0.424, 0.420, 0.414, 0.411, 0.406 + }; + + } + + private static class Osborne2Function extends MinpackFunction { + + public Osborne2Function(double[] startParams, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(65, startParams, theoreticalStartCost, + theoreticalMinCost, theoreticalMinParams); + } + + protected double[][] getJacobian() { + double x01 = parameters[0].getEstimate(); + double x02 = parameters[1].getEstimate(); + double x03 = parameters[2].getEstimate(); + double x04 = parameters[3].getEstimate(); + double x05 = parameters[4].getEstimate(); + double x06 = parameters[5].getEstimate(); + double x07 = parameters[6].getEstimate(); + double x08 = parameters[7].getEstimate(); + double x09 = parameters[8].getEstimate(); + double x10 = parameters[9].getEstimate(); + double x11 = parameters[10].getEstimate(); + double[][] jacobian = new double[m][]; + for (int i = 0; i < m; ++i) { + double temp = i / 10.0; + double tmp1 = Math.exp(-x05 * temp); + double tmp2 = Math.exp(-x06 * (temp - x09) * (temp - x09)); + double tmp3 = Math.exp(-x07 * (temp - x10) * (temp - x10)); + double tmp4 = Math.exp(-x08 * (temp - x11) * (temp - x11)); + jacobian[i] = new double[] { + -tmp1, + -tmp2, + -tmp3, + -tmp4, + temp * x01 * tmp1, + x02 * (temp - x09) * (temp - x09) * tmp2, + x03 * (temp - x10) * (temp - x10) * tmp3, + x04 * (temp - x11) * (temp - x11) * tmp4, + -2 * x02 * x06 * (temp - x09) * tmp2, + -2 * x03 * x07 * (temp - x10) * tmp3, + -2 * x04 * x08 * (temp - x11) * tmp4 + }; + } + return jacobian; + } + + protected double[] getResiduals() { + double x01 = parameters[0].getEstimate(); + double x02 = parameters[1].getEstimate(); + double x03 = parameters[2].getEstimate(); + double x04 = parameters[3].getEstimate(); + double x05 = parameters[4].getEstimate(); + double x06 = parameters[5].getEstimate(); + double x07 = parameters[6].getEstimate(); + double x08 = parameters[7].getEstimate(); + double x09 = parameters[8].getEstimate(); + double x10 = parameters[9].getEstimate(); + double x11 = parameters[10].getEstimate(); + double[] f = new double[m]; + for (int i = 0; i < m; ++i) { + double temp = i / 10.0; + double tmp1 = Math.exp(-x05 * temp); + double tmp2 = Math.exp(-x06 * (temp - x09) * (temp - x09)); + double tmp3 = Math.exp(-x07 * (temp - x10) * (temp - x10)); + double tmp4 = Math.exp(-x08 * (temp - x11) * (temp - x11)); + f[i] = y[i] - (x01 * tmp1 + x02 * tmp2 + x03 * tmp3 + x04 * tmp4); + } + return f; + } + + private static final double[] y = { + 1.366, 1.191, 1.112, 1.013, 0.991, + 0.885, 0.831, 0.847, 0.786, 0.725, + 0.746, 0.679, 0.608, 0.655, 0.616, + 0.606, 0.602, 0.626, 0.651, 0.724, + 0.649, 0.649, 0.694, 0.644, 0.624, + 0.661, 0.612, 0.558, 0.533, 0.495, + 0.500, 0.423, 0.395, 0.375, 0.372, + 0.391, 0.396, 0.405, 0.428, 0.429, + 0.523, 0.562, 0.607, 0.653, 0.672, + 0.708, 0.633, 0.668, 0.645, 0.632, + 0.591, 0.559, 0.597, 0.625, 0.739, + 0.710, 0.729, 0.720, 0.636, 0.581, + 0.428, 0.292, 0.162, 0.098, 0.054 + }; + + } + + public static Test suite() { + return new TestSuite(MinpackTest.class); + } + +} diff --git a/src/test/org/apache/commons/math/estimation/WeightedMeasurementTest.java b/src/test/org/apache/commons/math/estimation/WeightedMeasurementTest.java new file mode 100644 index 000000000..e2b9749ce --- /dev/null +++ b/src/test/org/apache/commons/math/estimation/WeightedMeasurementTest.java @@ -0,0 +1,127 @@ +/* + * 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.estimation; + +import org.apache.commons.math.estimation.EstimatedParameter; +import org.apache.commons.math.estimation.WeightedMeasurement; + +import junit.framework.*; + +@Deprecated +public class WeightedMeasurementTest + extends TestCase { + + public WeightedMeasurementTest(String name) { + super(name); + p1 = null; + p2 = null; + } + + public void testConstruction() { + WeightedMeasurement m = new MyMeasurement(3.0, theoretical() + 0.1, this); + checkValue(m.getWeight(), 3.0); + checkValue(m.getMeasuredValue(), theoretical() + 0.1); + } + + public void testIgnored() { + WeightedMeasurement m = new MyMeasurement(3.0, theoretical() + 0.1, this); + assertTrue(!m.isIgnored()); + m.setIgnored(true); + assertTrue(m.isIgnored()); + m.setIgnored(false); + assertTrue(!m.isIgnored()); + } + + public void testTheory() { + WeightedMeasurement m = new MyMeasurement(3.0, theoretical() + 0.1, this); + checkValue(m.getTheoreticalValue(), theoretical()); + checkValue(m.getResidual(), 0.1); + + double oldP1 = p1.getEstimate(); + p1.setEstimate(oldP1 + m.getResidual() / m.getPartial(p1)); + checkValue(m.getResidual(), 0.0); + p1.setEstimate(oldP1); + checkValue(m.getResidual(), 0.1); + + double oldP2 = p2.getEstimate(); + p2.setEstimate(oldP2 + m.getResidual() / m.getPartial(p2)); + checkValue(m.getResidual(), 0.0); + p2.setEstimate(oldP2); + checkValue(m.getResidual(), 0.1); + + } + + public static Test suite() { + return new TestSuite(WeightedMeasurementTest.class); + } + + public void setUp() { + p1 = new EstimatedParameter("p1", 1.0); + p2 = new EstimatedParameter("p2", 2.0); + } + + public void tearDown() { + p1 = null; + p2 = null; + } + + private void checkValue(double value, double expected) { + assertTrue(Math.abs(value - expected) < 1.0e-10); + } + + private double theoretical() { + return 3 * p1.getEstimate() - p2.getEstimate(); + } + + private double partial(EstimatedParameter p) { + if (p == p1) { + return 3.0; + } else if (p == p2) { + return -1.0; + } else { + return 0.0; + } + } + + private static class MyMeasurement + extends WeightedMeasurement { + + public MyMeasurement(double weight, double measuredValue, + WeightedMeasurementTest testInstance) { + super(weight, measuredValue); + this.testInstance = testInstance; + } + + public double getTheoreticalValue() { + return testInstance.theoretical(); + } + + public double getPartial(EstimatedParameter p) { + return testInstance.partial(p); + } + + private transient WeightedMeasurementTest testInstance; + + private static final long serialVersionUID = -246712922500792332L; + + } + + private EstimatedParameter p1; + private EstimatedParameter p2; + +}