From acd569595ee919d1afed3884cf7584cc609ad9ef Mon Sep 17 00:00:00 2001 From: Gilles Sadowski Date: Tue, 30 Jul 2013 15:04:22 +0000 Subject: [PATCH] MATH-1008 Created package "o.a.c.m.fitting.leastsquares" to contain a modified version of the contents of (to-be-deprecated) "o.a.c.m.optim.nonlinear.vector", the main purpose being to provide a new "fluent" API (cf. "withXxx" methods). Along the way, class "LevenbergMarquardtOptimizer" has been further cleaned up (e.g. removing protected fields and deprecated methods and using local variables instead of instance fields). An additional constructor was necessary in "BaseOptimizer". git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1508481 13f79535-47bb-0310-9956-ffa450edef68 --- .../AbstractLeastSquaresOptimizer.java | 362 ++++ .../leastsquares/GaussNewtonOptimizer.java | 325 ++++ .../LevenbergMarquardtOptimizer.java | 1172 +++++++++++++ .../leastsquares/WithConvergenceChecker.java | 39 + .../leastsquares/WithMaxEvaluations.java | 36 + .../leastsquares/WithMaxIterations.java | 36 + .../leastsquares/WithModelAndJacobian.java | 41 + .../fitting/leastsquares/WithStartPoint.java | 36 + .../fitting/leastsquares/WithTarget.java | 36 + .../fitting/leastsquares/WithWeight.java | 38 + .../fitting/leastsquares/package-info.java | 37 + .../commons/math3/optim/BaseOptimizer.java | 34 +- ...ractLeastSquaresOptimizerAbstractTest.java | 643 +++++++ .../AbstractLeastSquaresOptimizerTest.java | 109 ++ ...ctLeastSquaresOptimizerTestValidation.java | 308 ++++ .../fitting/leastsquares/CircleProblem.java | 175 ++ .../fitting/leastsquares/CircleVectorial.java | 94 ++ .../GaussNewtonOptimizerTest.java | 109 ++ .../LevenbergMarquardtOptimizerTest.java | 356 ++++ .../fitting/leastsquares/MinpackTest.java | 1499 +++++++++++++++++ .../RandomCirclePointGenerator.java | 91 + .../RandomStraightLinePointGenerator.java | 98 ++ .../StatisticalReferenceDataset.java | 370 ++++ .../StatisticalReferenceDatasetFactory.java | 201 +++ .../leastsquares/StraightLineProblem.java | 165 ++ .../math3/fitting/leastsquares/Hahn1.dat | 296 ++++ .../math3/fitting/leastsquares/Kirby2.dat | 211 +++ .../math3/fitting/leastsquares/Lanczos1.dat | 84 + .../math3/fitting/leastsquares/MGH17.dat | 93 + 29 files changed, 7092 insertions(+), 2 deletions(-) create mode 100644 src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizer.java create mode 100644 src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java create mode 100644 src/main/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizer.java create mode 100644 src/main/java/org/apache/commons/math3/fitting/leastsquares/WithConvergenceChecker.java create mode 100644 src/main/java/org/apache/commons/math3/fitting/leastsquares/WithMaxEvaluations.java create mode 100644 src/main/java/org/apache/commons/math3/fitting/leastsquares/WithMaxIterations.java create mode 100644 src/main/java/org/apache/commons/math3/fitting/leastsquares/WithModelAndJacobian.java create mode 100644 src/main/java/org/apache/commons/math3/fitting/leastsquares/WithStartPoint.java create mode 100644 src/main/java/org/apache/commons/math3/fitting/leastsquares/WithTarget.java create mode 100644 src/main/java/org/apache/commons/math3/fitting/leastsquares/WithWeight.java create mode 100644 src/main/java/org/apache/commons/math3/fitting/leastsquares/package-info.java create mode 100644 src/test/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizerAbstractTest.java create mode 100644 src/test/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizerTest.java create mode 100644 src/test/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizerTestValidation.java create mode 100644 src/test/java/org/apache/commons/math3/fitting/leastsquares/CircleProblem.java create mode 100644 src/test/java/org/apache/commons/math3/fitting/leastsquares/CircleVectorial.java create mode 100644 src/test/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizerTest.java create mode 100644 src/test/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizerTest.java create mode 100644 src/test/java/org/apache/commons/math3/fitting/leastsquares/MinpackTest.java create mode 100644 src/test/java/org/apache/commons/math3/fitting/leastsquares/RandomCirclePointGenerator.java create mode 100644 src/test/java/org/apache/commons/math3/fitting/leastsquares/RandomStraightLinePointGenerator.java create mode 100644 src/test/java/org/apache/commons/math3/fitting/leastsquares/StatisticalReferenceDataset.java create mode 100644 src/test/java/org/apache/commons/math3/fitting/leastsquares/StatisticalReferenceDatasetFactory.java create mode 100644 src/test/java/org/apache/commons/math3/fitting/leastsquares/StraightLineProblem.java create mode 100644 src/test/resources/org/apache/commons/math3/fitting/leastsquares/Hahn1.dat create mode 100644 src/test/resources/org/apache/commons/math3/fitting/leastsquares/Kirby2.dat create mode 100644 src/test/resources/org/apache/commons/math3/fitting/leastsquares/Lanczos1.dat create mode 100644 src/test/resources/org/apache/commons/math3/fitting/leastsquares/MGH17.dat diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizer.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizer.java new file mode 100644 index 000000000..2c771db82 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizer.java @@ -0,0 +1,362 @@ +/* + * 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.math3.fitting.leastsquares; + +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.MathUnsupportedOperationException; +import org.apache.commons.math3.analysis.MultivariateVectorFunction; +import org.apache.commons.math3.analysis.MultivariateMatrixFunction; +import org.apache.commons.math3.linear.ArrayRealVector; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.DiagonalMatrix; +import org.apache.commons.math3.linear.DecompositionSolver; +import org.apache.commons.math3.linear.MatrixUtils; +import org.apache.commons.math3.linear.QRDecomposition; +import org.apache.commons.math3.linear.EigenDecomposition; +import org.apache.commons.math3.optim.ConvergenceChecker; +import org.apache.commons.math3.optim.BaseOptimizer; +import org.apache.commons.math3.optim.PointVectorValuePair; +import org.apache.commons.math3.optim.OptimizationData; +import org.apache.commons.math3.util.FastMath; + +/** + * Base class for implementing least-squares optimizers. + * It provides methods for error estimation. + * + * @version $Id$ + * @since 3.3 + */ +public abstract class AbstractLeastSquaresOptimizer + extends BaseOptimizer { + /** Target values for the model function at optimum. */ + private final double[] target; + /** Weight matrix. */ + private final RealMatrix weight; + /** Model function. */ + private final MultivariateVectorFunction model; + /** Jacobian of the model function. */ + private final MultivariateMatrixFunction jacobian; + /** Square-root of the weight matrix. */ + private final RealMatrix weightSqrt; + /** Initial guess. */ + private final double[] start; + + /** + * @param target Observations. + * @param weight Weight of the observations. + * For performance, no defensive copy is performed. + * @param weightSqrt Square-root of the {@code weight} matrix. + * If {@code null}, it will be computed; otherwise it is the caller's + * responsibility that {@code weight} and {@code weightSqrt} are + * consistent. + * No defensive copy is performed. + * @param model ModelFunction. + * @param jacobian Jacobian of the model function. + * @param checker Convergence checker. + * @param start Initial guess. + * @param maxEval Maximum number of evaluations of the model + * function. + * @param maxIter Maximum number of iterations. + */ + protected AbstractLeastSquaresOptimizer(double[] target, + RealMatrix weight, + RealMatrix weightSqrt, + MultivariateVectorFunction model, + MultivariateMatrixFunction jacobian, + ConvergenceChecker checker, + double[] start, + int maxEval, + int maxIter) { + super(checker, maxEval, maxIter); + + this.target = target; + this.weight = weight; + this.model = model; + this.jacobian = jacobian; + this.start = start; + + this.weightSqrt = weightSqrt == null ? + (weight == null ? + null : squareRoot(weight)) : weightSqrt; + } + + /** + * Gets the target values. + * + * @return the target values. + */ + public double[] getTarget() { + return target == null ? null : target.clone(); + } + + /** + * Gets the initial guess. + * + * @return the initial guess values. + */ + public double[] getStart() { + return start == null ? null : start.clone(); + } + + /** + * Gets the square-root of the weight matrix. + * + * @return the square-root of the weight matrix. + */ + public RealMatrix getWeightSquareRoot() { + return weightSqrt == null ? null : weightSqrt.copy(); + } + + /** + * Gets the model function. + * + * @return the model function. + */ + public MultivariateVectorFunction getModel() { + return model; + } + + /** + * Gets the model function's Jacobian. + * + * @return the Jacobian. + */ + public MultivariateMatrixFunction getJacobian() { + return jacobian; + } + + /** + * Get the covariance matrix of the optimized parameters. + *
+ * Note that this operation involves the inversion of the + * JTJ matrix, where {@code J} is the + * Jacobian matrix. + * The {@code threshold} parameter is a way for the caller to specify + * that the result of this computation should be considered meaningless, + * and thus trigger an exception. + * + * @param params Model parameters. + * @param threshold Singularity threshold. + * @return the covariance matrix. + * @throws org.apache.commons.math3.linear.SingularMatrixException + * if the covariance matrix cannot be computed (singular problem). + */ + public double[][] computeCovariances(double[] params, + double threshold) { + // Set up the Jacobian. + final RealMatrix j = computeWeightedJacobian(params); + + // Compute transpose(J)J. + final RealMatrix jTj = j.transpose().multiply(j); + + // Compute the covariances matrix. + final DecompositionSolver solver + = new QRDecomposition(jTj, threshold).getSolver(); + return solver.getInverse().getData(); + } + + /** + * Computes an estimate of the standard deviation of the parameters. The + * returned values are the square root of the diagonal coefficients of the + * covariance matrix, {@code sd(a[i]) ~= sqrt(C[i][i])}, where {@code a[i]} + * is the optimized value of the {@code i}-th parameter, and {@code C} is + * the covariance matrix. + * + * @param params Model parameters. + * @param covarianceSingularityThreshold Singularity threshold (see + * {@link #computeCovariances(double[],double) computeCovariances}). + * @return an estimate of the standard deviation of the optimized parameters + * @throws org.apache.commons.math3.linear.SingularMatrixException + * if the covariance matrix cannot be computed. + */ + public double[] computeSigma(double[] params, + double covarianceSingularityThreshold) { + final int nC = params.length; + final double[] sig = new double[nC]; + final double[][] cov = computeCovariances(params, covarianceSingularityThreshold); + for (int i = 0; i < nC; ++i) { + sig[i] = FastMath.sqrt(cov[i][i]); + } + return sig; + } + + /** + * Gets the weight matrix of the observations. + * + * @return the weight matrix. + */ + public RealMatrix getWeight() { + return weight.copy(); + } + + /** + * Computes the normalized cost. + * It is the square-root of the sum of squared of the residuals, divided + * by the number of measurements. + * + * @param params Model function parameters. + * @return the cost. + */ + public double computeRMS(double[] params) { + final double cost = computeCost(computeResiduals(getModel().value(params))); + return FastMath.sqrt(cost * cost / target.length); + } + + /** + * Calling this method will raise an exception. + * + * @param optData Obsolete. + * @return nothing. + * @throws MathUnsupportedOperationException if called. + * @deprecated Do not use this method. + */ + @Deprecated + @Override + public PointVectorValuePair optimize(OptimizationData... optData) + throws MathUnsupportedOperationException { + throw new MathUnsupportedOperationException(); + } + + /** + * Gets a reference to the corresponding field. + * Altering it could jeopardize the consistency of this class. + * + * @return the reference. + */ + protected double[] getTargetInternal() { + return target; + } + + /** + * Gets a reference to the corresponding field. + * Altering it could jeopardize the consistency of this class. + * + * @return the reference. + */ + protected RealMatrix getWeightInternal() { + return weight; + } + + /** + * Gets a reference to the corresponding field. + * Altering it could jeopardize the consistency of this class. + * + * @return the reference. + */ + protected RealMatrix getWeightSquareRootInternal() { + return weightSqrt; + } + + /** + * Computes the objective function value. + * This method must be called by subclasses to enforce the + * evaluation counter limit. + * + * @param params Point at which the objective function must be evaluated. + * @return the objective function value at the specified point. + * @throws org.apache.commons.math3.exception.TooManyEvaluationsException + * if the maximal number of evaluations (of the model vector function) is + * exceeded. + */ + protected double[] computeObjectiveValue(double[] params) { + super.incrementEvaluationCount(); + return model.value(params); + } + + /** + * Computes the weighted Jacobian matrix. + * + * @param params Model parameters at which to compute the Jacobian. + * @return the weighted Jacobian: W1/2 J. + * @throws DimensionMismatchException if the Jacobian dimension does not + * match problem dimension. + */ + protected RealMatrix computeWeightedJacobian(double[] params) { + return weightSqrt.multiply(MatrixUtils.createRealMatrix(computeJacobian(params))); + } + + /** + * Computes the Jacobian matrix. + * + * @param params Point at which the Jacobian must be evaluated. + * @return the Jacobian at the specified point. + */ + protected double[][] computeJacobian(final double[] params) { + return jacobian.value(params); + } + + /** + * Computes the cost. + * + * @param residuals Residuals. + * @return the cost. + * @see #computeResiduals(double[]) + */ + protected double computeCost(double[] residuals) { + final ArrayRealVector r = new ArrayRealVector(residuals); + return FastMath.sqrt(r.dotProduct(weight.operate(r))); + } + + /** + * Computes the residuals. + * The residual is the difference between the observed (target) + * values and the model (objective function) value. + * There is one residual for each element of the vector-valued + * function. + * + * @param objectiveValue Value of the the objective function. This is + * the value returned from a call to + * {@link #computeObjectiveValue(double[]) computeObjectiveValue} + * (whose array argument contains the model parameters). + * @return the residuals. + * @throws DimensionMismatchException if {@code params} has a wrong + * length. + */ + protected double[] computeResiduals(double[] objectiveValue) { + if (objectiveValue.length != target.length) { + throw new DimensionMismatchException(target.length, + objectiveValue.length); + } + + final double[] residuals = new double[target.length]; + for (int i = 0; i < target.length; i++) { + residuals[i] = target[i] - objectiveValue[i]; + } + + return residuals; + } + + /** + * Computes the square-root of the weight matrix. + * + * @param m Symmetric, positive-definite (weight) matrix. + * @return the square-root of the weight matrix. + */ + private RealMatrix squareRoot(RealMatrix m) { + if (m instanceof DiagonalMatrix) { + final int dim = m.getRowDimension(); + final RealMatrix sqrtM = new DiagonalMatrix(dim); + for (int i = 0; i < dim; i++) { + sqrtM.setEntry(i, i, FastMath.sqrt(m.getEntry(i, i))); + } + return sqrtM; + } else { + final EigenDecomposition dec = new EigenDecomposition(m); + return dec.getSquareRoot(); + } + } +} diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java new file mode 100644 index 000000000..8c6988461 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizer.java @@ -0,0 +1,325 @@ +/* + * 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.math3.fitting.leastsquares; + +import org.apache.commons.math3.analysis.MultivariateVectorFunction; +import org.apache.commons.math3.analysis.MultivariateMatrixFunction; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.ConvergenceException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.exception.MathInternalError; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.linear.ArrayRealVector; +import org.apache.commons.math3.linear.BlockRealMatrix; +import org.apache.commons.math3.linear.DecompositionSolver; +import org.apache.commons.math3.linear.LUDecomposition; +import org.apache.commons.math3.linear.QRDecomposition; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.SingularMatrixException; +import org.apache.commons.math3.optim.ConvergenceChecker; +import org.apache.commons.math3.optim.PointVectorValuePair; + +/** + * Gauss-Newton least-squares solver. + * + *

+ * This class solve a least-square problem by solving the normal equations + * of the linearized problem at each iteration. Either LU decomposition or + * QR decomposition can be used to solve the normal equations. LU decomposition + * is faster but QR decomposition is more robust for difficult problems. + *

+ * + * @version $Id$ + * @since 3.3 + * + */ +public class GaussNewtonOptimizer extends AbstractLeastSquaresOptimizer + implements WithTarget, + WithWeight, + WithModelAndJacobian, + WithConvergenceChecker, + WithStartPoint, + WithMaxIterations, + WithMaxEvaluations { + /** Indicator for using LU decomposition. */ + private final boolean useLU; + + /** + * Constructor called by the various {@code withXxx} methods. + * + * @param target Observations. + * @param weight Weight of the observations. + * For performance, no defensive copy is performed. + * @param weightSqrt Square-root of the {@code weight} matrix. + * If {@code null}, it will be computed; otherwise it is the caller's + * responsibility that {@code weight} and {@code weightSqrt} are + * consistent. + * No defensive copy is performed. + * @param model ModelFunction. + * @param jacobian Jacobian of the model function. + * @param checker Convergence checker. + * @param start Initial guess. + * @param maxEval Maximum number of evaluations of the model + * function. + * @param maxIter Maximum number of iterations. + * @param useLU Whether to use LU decomposition. + */ + private GaussNewtonOptimizer(double[] target, + RealMatrix weight, + RealMatrix weightSqrt, + MultivariateVectorFunction model, + MultivariateMatrixFunction jacobian, + ConvergenceChecker checker, + double[] start, + int maxEval, + int maxIter, + boolean useLU) { + super(target, weight, weightSqrt, model, jacobian, checker, start, maxEval, maxIter); + + this.useLU = useLU; + } + + /** + * Creates a bare-bones instance. + * Several calls to {@code withXxx} methods are necessary to obtain + * an object with all necessary fields set to sensible values. + *
+ * The default for the algorithm is to solve the normal equations + * using LU decomposition. + * + * @return an instance of this class. + */ + public static GaussNewtonOptimizer create() { + return new GaussNewtonOptimizer(null, null, null, null, null, null, null, + 0, 0, true); + } + + /** {@inheritDoc} */ + public GaussNewtonOptimizer withTarget(double[] target) { + return new GaussNewtonOptimizer(target, + getWeightInternal(), + getWeightSquareRootInternal(), + getModel(), + getJacobian(), + getConvergenceChecker(), + getStart(), + getMaxEvaluations(), + getMaxIterations(), + useLU); + } + + /** {@inheritDoc} */ + public GaussNewtonOptimizer withWeight(RealMatrix weight) { + return new GaussNewtonOptimizer(getTargetInternal(), + weight, + null, + getModel(), + getJacobian(), + getConvergenceChecker(), + getStart(), + getMaxEvaluations(), + getMaxIterations(), + useLU); + } + + /** {@inheritDoc} */ + public GaussNewtonOptimizer withModelAndJacobian(MultivariateVectorFunction model, + MultivariateMatrixFunction jacobian) { + return new GaussNewtonOptimizer(getTargetInternal(), + getWeightInternal(), + getWeightSquareRootInternal(), + model, + jacobian, + getConvergenceChecker(), + getStart(), + getMaxEvaluations(), + getMaxIterations(), + useLU); + } + + /** {@inheritDoc} */ + public GaussNewtonOptimizer withConvergenceChecker(ConvergenceChecker checker) { + return new GaussNewtonOptimizer(getTarget(), + getWeightInternal(), + getWeightSquareRootInternal(), + getModel(), + getJacobian(), + checker, + getStart(), + getMaxEvaluations(), + getMaxIterations(), + useLU); + } + + /** {@inheritDoc} */ + public GaussNewtonOptimizer withStartPoint(double[] start) { + return new GaussNewtonOptimizer(getTarget(), + getWeightInternal(), + getWeightSquareRootInternal(), + getModel(), + getJacobian(), + getConvergenceChecker(), + start, + getMaxEvaluations(), + getMaxIterations(), + useLU); + } + + /** {@inheritDoc} */ + public GaussNewtonOptimizer withMaxIterations(int maxIter) { + return new GaussNewtonOptimizer(getTarget(), + getWeightInternal(), + getWeightSquareRootInternal(), + getModel(), + getJacobian(), + getConvergenceChecker(), + getStart(), + getMaxEvaluations(), + maxIter, + useLU); + } + + /** {@inheritDoc} */ + public GaussNewtonOptimizer withMaxEvaluations(int maxEval) { + return new GaussNewtonOptimizer(getTarget(), + getWeightInternal(), + getWeightSquareRootInternal(), + getModel(), + getJacobian(), + getConvergenceChecker(), + getStart(), + maxEval, + getMaxIterations(), + useLU); + } + + /** + * Creates a new instance. + * + * @param withLU Whether to use LU decomposition. + * @return a new instance with all fields identical to this instance except + * for the givens arguments. + */ + public GaussNewtonOptimizer withLU(boolean withLU) { + return new GaussNewtonOptimizer(getTarget(), + getWeightInternal(), + getWeightSquareRootInternal(), + getModel(), + getJacobian(), + getConvergenceChecker(), + getStart(), + getMaxEvaluations(), + getMaxIterations(), + withLU); + } + + /** {@inheritDoc} */ + @Override + public PointVectorValuePair doOptimize() { + final ConvergenceChecker checker + = getConvergenceChecker(); + + // Computation will be useless without a checker (see "for-loop"). + if (checker == null) { + throw new NullArgumentException(); + } + + final double[] targetValues = getTarget(); + final int nR = targetValues.length; // Number of observed data. + + final RealMatrix weightMatrix = getWeight(); + if (weightMatrix.getRowDimension() != nR) { + throw new DimensionMismatchException(weightMatrix.getRowDimension(), nR); + } + if (weightMatrix.getColumnDimension() != nR) { + throw new DimensionMismatchException(weightMatrix.getColumnDimension(), nR); + } + + // Diagonal of the weight matrix. + final double[] residualsWeights = new double[nR]; + for (int i = 0; i < nR; i++) { + residualsWeights[i] = weightMatrix.getEntry(i, i); + } + + final double[] currentPoint = getStart(); + final int nC = currentPoint.length; + + // iterate until convergence is reached + PointVectorValuePair current = null; + for (boolean converged = false; !converged;) { + incrementIterationCount(); + + // evaluate the objective function and its jacobian + PointVectorValuePair previous = current; + // Value of the objective function at "currentPoint". + final double[] currentObjective = computeObjectiveValue(currentPoint); + final double[] currentResiduals = computeResiduals(currentObjective); + final RealMatrix weightedJacobian = computeWeightedJacobian(currentPoint); + current = new PointVectorValuePair(currentPoint, currentObjective); + + // build the linear problem + final double[] b = new double[nC]; + final double[][] a = new double[nC][nC]; + for (int i = 0; i < nR; ++i) { + + final double[] grad = weightedJacobian.getRow(i); + final double weight = residualsWeights[i]; + final double residual = currentResiduals[i]; + + // compute the normal equation + final double wr = weight * residual; + for (int j = 0; j < nC; ++j) { + b[j] += wr * grad[j]; + } + + // build the contribution matrix for measurement i + for (int k = 0; k < nC; ++k) { + double[] ak = a[k]; + double wgk = weight * grad[k]; + for (int l = 0; l < nC; ++l) { + ak[l] += wgk * grad[l]; + } + } + } + + // Check convergence. + if (previous != null) { + converged = checker.converged(getIterations(), previous, current); + if (converged) { + return current; + } + } + + try { + // solve the linearized least squares problem + RealMatrix mA = new BlockRealMatrix(a); + DecompositionSolver solver = useLU ? + new LUDecomposition(mA).getSolver() : + new QRDecomposition(mA).getSolver(); + final double[] dX = solver.solve(new ArrayRealVector(b, false)).toArray(); + // update the estimated parameters + for (int i = 0; i < nC; ++i) { + currentPoint[i] += dX[i]; + } + } catch (SingularMatrixException e) { + throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM); + } + } + // Must never happen. + throw new MathInternalError(); + } +} diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizer.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizer.java new file mode 100644 index 000000000..4614fbc53 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizer.java @@ -0,0 +1,1172 @@ +/* + * 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.math3.fitting.leastsquares; + +import java.util.Arrays; +import org.apache.commons.math3.analysis.MultivariateVectorFunction; +import org.apache.commons.math3.analysis.MultivariateMatrixFunction; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.exception.ConvergenceException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.optim.ConvergenceChecker; +import org.apache.commons.math3.optim.PointVectorValuePair; +import org.apache.commons.math3.util.Precision; +import org.apache.commons.math3.util.FastMath; + + +/** + * This class solves a least-squares problem using the Levenberg-Marquardt + * algorithm. + * + *

This implementation should work even for over-determined systems + * (i.e. systems having more point than equations). Over-determined systems + * are solved by ignoring the point 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, the use of + * inherited convergence checker 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 authors of the original fortran version are: + *

    + *
  • Argonne National Laboratory. MINPACK project. March 1980
    • + *
  • Burton S. Garbow
    • + *
  • Kenneth E. Hillstrom
    • + *
  • Jorge J. More
    • + *
+ * 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. + *
    + * + * @version $Id$ + * @since 2.0 + */ +public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer + implements WithTarget, + WithWeight, + WithModelAndJacobian, + WithConvergenceChecker, + WithStartPoint, + WithMaxIterations, + WithMaxEvaluations { + /** Twice the "epsilon machine". */ + private static final double TWO_EPS = 2 * Precision.EPSILON; + /** Positive input variable used in determining the initial step bound. */ + private final double initialStepBoundFactor; + /** Desired relative error in the sum of squares. */ + private final double costRelativeTolerance; + /** Desired relative error in the approximate solution parameters. */ + private final double parRelativeTolerance; + /** Desired max cosine on the orthogonality between the function vector + * and the columns of the jacobian. */ + private final double orthoTolerance; + /** Threshold for QR ranking. */ + private final double qrRankingThreshold; + /** Levenberg-Marquardt parameter. */ + private double lmPar; + /** Parameters evolution direction associated with lmPar. */ + private double[] lmDir; + + /** + * Constructor called by the various {@code withXxx} methods. + * + * @param target Observations. + * @param weight Weight of the observations. + * For performance, no defensive copy is performed. + * @param weightSqrt Square-root of the {@code weight} matrix. + * If {@code null}, it will be computed; otherwise it is the caller's + * responsibility that {@code weight} and {@code weightSqrt} are + * consistent. + * No defensive copy is performed. + * @param model ModelFunction. + * @param jacobian Jacobian of the model function. + * @param checker Convergence checker. + * @param start Initial guess. + * @param maxEval Maximum number of evaluations of the model + * function. + * @param maxIter Maximum number of iterations. + * @param initialStepBoundFactor Positive input variable used in + * determining the initial step bound. This bound is set to the + * product of initialStepBoundFactor and the euclidean norm of + * {@code diag * x} if non-zero, or else to {@code initialStepBoundFactor} + * itself. In most cases factor should lie in the interval + * {@code (0.1, 100.0)}. {@code 100} is a generally recommended value. + * @param costRelativeTolerance Desired relative error in the sum of + * squares. + * @param parRelativeTolerance Desired relative error in the approximate + * solution parameters. + * @param orthoTolerance Desired max cosine on the orthogonality between + * the function vector and the columns of the Jacobian. + * @param threshold Desired threshold for QR ranking. If the squared norm + * of a column vector is smaller or equal to this threshold during QR + * decomposition, it is considered to be a zero vector and hence the rank + * of the matrix is reduced. + */ + private LevenbergMarquardtOptimizer(double[] target, + RealMatrix weight, + RealMatrix weightSqrt, + MultivariateVectorFunction model, + MultivariateMatrixFunction jacobian, + ConvergenceChecker checker, + double[] start, + int maxEval, + int maxIter, + double initialStepBoundFactor, + double costRelativeTolerance, + double parRelativeTolerance, + double orthoTolerance, + double threshold) { + super(target, weight, weightSqrt, model, jacobian, checker, start, maxEval, maxIter); + + this.initialStepBoundFactor = initialStepBoundFactor; + this.costRelativeTolerance = costRelativeTolerance; + this.parRelativeTolerance = parRelativeTolerance; + this.orthoTolerance = orthoTolerance; + this.qrRankingThreshold = threshold; + } + + /** + * Creates a bare-bones instance. + * Several calls to {@code withXxx} methods are necessary to obtain + * an object with all necessary fields set to sensible values. + *
    + * The default values for the algorithm settings are: + *
      + *
    • Initial step bound factor: 100
      • + *
    • Cost relative tolerance: 1e-10
      • + *
    • Parameters relative tolerance: 1e-10
      • + *
    • Orthogonality tolerance: 1e-10
      • + *
    • QR ranking threshold: {@link Precision#SAFE_MIN}
      • + *
    + * + * @return an instance of this class. + */ + public static LevenbergMarquardtOptimizer create() { + return new LevenbergMarquardtOptimizer(null, null, null, null, null, null, null, + 0, 0, 100, 1e-10, 1e-10, 1e-10, + Precision.SAFE_MIN); + } + + /** {@inheritDoc} */ + public LevenbergMarquardtOptimizer withTarget(double[] target) { + return new LevenbergMarquardtOptimizer(target, + getWeightInternal(), + getWeightSquareRootInternal(), + getModel(), + getJacobian(), + getConvergenceChecker(), + getStart(), + getMaxEvaluations(), + getMaxIterations(), + initialStepBoundFactor, + costRelativeTolerance, + parRelativeTolerance, + orthoTolerance, + qrRankingThreshold); + } + + /** {@inheritDoc} */ + public LevenbergMarquardtOptimizer withWeight(RealMatrix weight) { + return new LevenbergMarquardtOptimizer(getTargetInternal(), + weight, + null, + getModel(), + getJacobian(), + getConvergenceChecker(), + getStart(), + getMaxEvaluations(), + getMaxIterations(), + initialStepBoundFactor, + costRelativeTolerance, + parRelativeTolerance, + orthoTolerance, + qrRankingThreshold); + } + + /** {@inheritDoc} */ + public LevenbergMarquardtOptimizer withModelAndJacobian(MultivariateVectorFunction model, + MultivariateMatrixFunction jacobian) { + return new LevenbergMarquardtOptimizer(getTargetInternal(), + getWeightInternal(), + getWeightSquareRootInternal(), + model, + jacobian, + getConvergenceChecker(), + getStart(), + getMaxEvaluations(), + getMaxIterations(), + initialStepBoundFactor, + costRelativeTolerance, + parRelativeTolerance, + orthoTolerance, + qrRankingThreshold); + } + + /** {@inheritDoc} */ + public LevenbergMarquardtOptimizer withConvergenceChecker(ConvergenceChecker checker) { + return new LevenbergMarquardtOptimizer(getTarget(), + getWeightInternal(), + getWeightSquareRootInternal(), + getModel(), + getJacobian(), + checker, + getStart(), + getMaxEvaluations(), + getMaxIterations(), + initialStepBoundFactor, + costRelativeTolerance, + parRelativeTolerance, + orthoTolerance, + qrRankingThreshold); + } + + /** {@inheritDoc} */ + public LevenbergMarquardtOptimizer withStartPoint(double[] start) { + return new LevenbergMarquardtOptimizer(getTarget(), + getWeightInternal(), + getWeightSquareRootInternal(), + getModel(), + getJacobian(), + getConvergenceChecker(), + start, + getMaxEvaluations(), + getMaxIterations(), + initialStepBoundFactor, + costRelativeTolerance, + parRelativeTolerance, + orthoTolerance, + qrRankingThreshold); + } + + /** {@inheritDoc} */ + public LevenbergMarquardtOptimizer withMaxIterations(int maxIter) { + return new LevenbergMarquardtOptimizer(getTarget(), + getWeightInternal(), + getWeightSquareRootInternal(), + getModel(), + getJacobian(), + getConvergenceChecker(), + getStart(), + getMaxEvaluations(), + maxIter, + initialStepBoundFactor, + costRelativeTolerance, + parRelativeTolerance, + orthoTolerance, + qrRankingThreshold); + } + + /** {@inheritDoc} */ + public LevenbergMarquardtOptimizer withMaxEvaluations(int maxEval) { + return new LevenbergMarquardtOptimizer(getTarget(), + getWeightInternal(), + getWeightSquareRootInternal(), + getModel(), + getJacobian(), + getConvergenceChecker(), + getStart(), + maxEval, + getMaxIterations(), + initialStepBoundFactor, + costRelativeTolerance, + parRelativeTolerance, + orthoTolerance, + qrRankingThreshold); + } + + /** + * Creates a new instance. + * + * @param initStepBoundFactor Positive input variable used in + * determining the initial step bound. This bound is set to the + * product of initialStepBoundFactor and the euclidean norm of + * {@code diag * x} if non-zero, or else to {@code initialStepBoundFactor} + * itself. In most cases factor should lie in the interval + * {@code (0.1, 100.0)}. {@code 100} is a generally recommended value. + * @param costRelTol Desired relative error in the sum of squares. + * @param parRelTol Desired relative error in the approximate solution + * parameters. + * @param orthoTol Desired max cosine on the orthogonality between + * the function vector and the columns of the Jacobian. + * @param threshold Desired threshold for QR ranking. If the squared norm + * of a column vector is smaller or equal to this threshold during QR + * decomposition, it is considered to be a zero vector and hence the rank + * of the matrix is reduced. + * @return a new instance with all fields identical to this instance except + * for the givens arguments. + */ + public LevenbergMarquardtOptimizer withTuningParameters(double initStepBoundFactor, + double costRelTol, + double parRelTol, + double orthoTol, + double threshold) { + return new LevenbergMarquardtOptimizer(getTarget(), + getWeightInternal(), + getWeightSquareRootInternal(), + getModel(), + getJacobian(), + getConvergenceChecker(), + getStart(), + getMaxEvaluations(), + getMaxIterations(), + initStepBoundFactor, + costRelTol, + parRelTol, + orthoTol, + threshold); + } + + /** + * Gets the value of a tuning parameter. + * @see #withTuningParameters(double,double,double,double,double) + * + * @return the parameter's value. + */ + public double getInitialStepBoundFactor() { + return initialStepBoundFactor; + } + + /** + * Gets the value of a tuning parameter. + * @see #withTuningParameters(double,double,double,double,double) + * + * @return the parameter's value. + */ + public double getCostRelativeTolerance() { + return costRelativeTolerance; + } + + /** + * Gets the value of a tuning parameter. + * @see #withTuningParameters(double,double,double,double,double) + * + * @return the parameter's value. + */ + public double getParRelativeTolerance() { + return parRelativeTolerance; + } + + /** + * Gets the value of a tuning parameter. + * @see #withTuningParameters(double,double,double,double,double) + * + * @return the parameter's value. + */ + public double getOrthoTolerance() { + return orthoTolerance; + } + + /** + * Gets the value of a tuning parameter. + * @see #withTuningParameters(double,double,double,double,double) + * + * @return the parameter's value. + */ + public double getRankingThreshold() { + return qrRankingThreshold; + } + + /** {@inheritDoc} */ + @Override + protected PointVectorValuePair doOptimize() { + final int nR = getTarget().length; // Number of observed data. + final double[] currentPoint = getStart(); + final int nC = currentPoint.length; // Number of parameters. + + // arrays shared with the other private methods + final int solvedCols = FastMath.min(nR, nC); + lmDir = new double[nC]; + lmPar = 0; + + // local point + double delta = 0; + double xNorm = 0; + double[] diag = new double[nC]; + double[] oldX = new double[nC]; + double[] oldRes = new double[nR]; + double[] oldObj = new double[nR]; + double[] qtf = new double[nR]; + double[] work1 = new double[nC]; + double[] work2 = new double[nC]; + double[] work3 = new double[nC]; + + final RealMatrix weightMatrixSqrt = getWeightSquareRoot(); + + // Evaluate the function at the starting point and calculate its norm. + double[] currentObjective = computeObjectiveValue(currentPoint); + double[] currentResiduals = computeResiduals(currentObjective); + PointVectorValuePair current = new PointVectorValuePair(currentPoint, currentObjective); + double currentCost = computeCost(currentResiduals); + + // Outer loop. + boolean firstIteration = true; + final ConvergenceChecker checker = getConvergenceChecker(); + while (true) { + incrementIterationCount(); + + final PointVectorValuePair previous = current; + + // QR decomposition of the jacobian matrix + final InternalData internalData = qrDecomposition(computeWeightedJacobian(currentPoint), + solvedCols); + final double[][] weightedJacobian = internalData.weightedJacobian; + final int[] permutation = internalData.permutation; + final double[] diagR = internalData.diagR; + final double[] jacNorm = internalData.jacNorm; + + double[] weightedResidual = weightMatrixSqrt.operate(currentResiduals); + for (int i = 0; i < nR; i++) { + qtf[i] = weightedResidual[i]; + } + + // compute Qt.res + qTy(qtf, internalData); + + // 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]; + weightedJacobian[k][pk] = diagR[pk]; + } + + if (firstIteration) { + // scale the point according to the norms of the columns + // of the initial jacobian + xNorm = 0; + for (int k = 0; k < nC; ++k) { + double dk = jacNorm[k]; + if (dk == 0) { + dk = 1.0; + } + double xk = dk * currentPoint[k]; + xNorm += xk * xk; + diag[k] = dk; + } + xNorm = FastMath.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 (currentCost != 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; i <= j; ++i) { + sum += weightedJacobian[i][pj] * qtf[i]; + } + maxCosine = FastMath.max(maxCosine, FastMath.abs(sum) / (s * currentCost)); + } + } + } + if (maxCosine <= orthoTolerance) { + // Convergence has been reached. + return current; + } + + // rescale if necessary + for (int j = 0; j < nC; ++j) { + diag[j] = FastMath.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] = currentPoint[pj]; + } + final double previousCost = currentCost; + double[] tmpVec = weightedResidual; + weightedResidual = oldRes; + oldRes = tmpVec; + tmpVec = currentObjective; + currentObjective = oldObj; + oldObj = tmpVec; + + // determine the Levenberg-Marquardt parameter + determineLMParameter(qtf, delta, diag, + internalData, solvedCols, + 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]; + currentPoint[pj] = oldX[pj] + lmDir[pj]; + double s = diag[pj] * lmDir[pj]; + lmNorm += s * s; + } + lmNorm = FastMath.sqrt(lmNorm); + // on the first iteration, adjust the initial step bound. + if (firstIteration) { + delta = FastMath.min(delta, lmNorm); + } + + // Evaluate the function at x + p and calculate its norm. + currentObjective = computeObjectiveValue(currentPoint); + currentResiduals = computeResiduals(currentObjective); + current = new PointVectorValuePair(currentPoint, currentObjective); + currentCost = computeCost(currentResiduals); + + // compute the scaled actual reduction + double actRed = -1.0; + if (0.1 * currentCost < previousCost) { + double r = currentCost / 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; i <= j; ++i) { + work1[i] += weightedJacobian[i][pj] * 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 * currentCost >= previousCost) || (tmp < 0.1)) { + tmp = 0.1; + } + delta = tmp * FastMath.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 < nC; ++k) { + double xK = diag[k] * currentPoint[k]; + xNorm += xK * xK; + } + xNorm = FastMath.sqrt(xNorm); + + // tests for convergence. + if (checker != null && checker.converged(getIterations(), previous, current)) { + return current; + } + } else { + // failed iteration, reset the previous values + currentCost = previousCost; + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + currentPoint[pj] = oldX[pj]; + } + tmpVec = weightedResidual; + weightedResidual = oldRes; + oldRes = tmpVec; + tmpVec = currentObjective; + currentObjective = oldObj; + oldObj = tmpVec; + // Reset "current" to previous values. + current = new PointVectorValuePair(currentPoint, currentObjective); + } + + // Default convergence criteria. + if ((FastMath.abs(actRed) <= costRelativeTolerance && + preRed <= costRelativeTolerance && + ratio <= 2.0) || + delta <= parRelativeTolerance * xNorm) { + return current; + } + + // tests for termination and stringent tolerances + if (FastMath.abs(actRed) <= TWO_EPS && + preRed <= TWO_EPS && + ratio <= 2.0) { + throw new ConvergenceException(LocalizedFormats.TOO_SMALL_COST_RELATIVE_TOLERANCE, + costRelativeTolerance); + } else if (delta <= TWO_EPS * xNorm) { + throw new ConvergenceException(LocalizedFormats.TOO_SMALL_PARAMETERS_RELATIVE_TOLERANCE, + parRelativeTolerance); + } else if (maxCosine <= TWO_EPS) { + throw new ConvergenceException(LocalizedFormats.TOO_SMALL_ORTHOGONALITY_TOLERANCE, + orthoTolerance); + } + } + } + } + + /** + * Holds internal data. + * This structure was created so that all optimizer fields can be "final". + * Code should be further refactored in order to not pass around arguments + * that will modified in-place (cf. "work" arrays). + */ + private static class InternalData { + /** Weighted Jacobian. */ + final double[][] weightedJacobian; + /** Columns permutation array. */ + final int[] permutation; + /** Rank of the Jacobian matrix. */ + final int rank; + /** Diagonal elements of the R matrix in the QR decomposition. */ + final double[] diagR; + /** Norms of the columns of the jacobian matrix. */ + final double[] jacNorm; + /** Coefficients of the Householder transforms vectors. */ + final double[] beta; + + /** + * @param weightedJacobian Weighted Jacobian. + * @param permutation Columns permutation array. + * @param rank Rank of the Jacobian matrix. + * @param diagR Diagonal elements of the R matrix in the QR decomposition. + * @param jacNorm Norms of the columns of the jacobian matrix. + * @param beta Coefficients of the Householder transforms vectors. + */ + InternalData(double[][] weightedJacobian, + int[] permutation, + int rank, + double[] diagR, + double[] jacNorm, + double[] beta) { + this.weightedJacobian = weightedJacobian; + this.permutation = permutation; + this.rank = rank; + this.diagR = diagR; + this.jacNorm = jacNorm; + this.beta = beta; + } + } + + /** + * Determines 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:

    + *
      + *
    • Argonne National Laboratory. MINPACK project. March 1980
      • + *
    • Burton S. Garbow
      • + *
    • Kenneth E. Hillstrom
      • + *
    • Jorge J. More
      • + *
    + *

    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 internalData Data (modified in-place in this method). + * @param solvedCols Number of solved point. + * @param work1 work array + * @param work2 work array + * @param work3 work array + */ + private void determineLMParameter(double[] qy, double delta, double[] diag, + InternalData internalData, int solvedCols, + double[] work1, double[] work2, double[] work3) { + final double[][] weightedJacobian = internalData.weightedJacobian; + final int[] permutation = internalData.permutation; + final int rank = internalData.rank; + final double[] diagR = internalData.diagR; + + final int nC = weightedJacobian[0].length; + + // 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 < nC; ++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; i < k; ++i) { + lmDir[permutation[i]] -= ypk * weightedJacobian[i][pk]; + } + 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 = FastMath.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; + double 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; i < j; ++i) { + sum += weightedJacobian[i][pj] * 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; i <= j; ++i) { + sum += weightedJacobian[i][pj] * qy[i]; + } + sum /= diag[pj]; + sum2 += sum * sum; + } + double gNorm = FastMath.sqrt(sum2); + double paru = gNorm / delta; + if (paru == 0) { + paru = Precision.SAFE_MIN / FastMath.min(delta, 0.1); + } + + // if the input par lies outside of the interval (parl,paru), + // set par to the closer endpoint + lmPar = FastMath.min(paru, FastMath.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 = FastMath.max(Precision.SAFE_MIN, 0.001 * paru); + } + double sPar = FastMath.sqrt(lmPar); + for (int j = 0; j < solvedCols; ++j) { + int pj = permutation[j]; + work1[pj] = sPar * diag[pj]; + } + determineLMDirection(qy, work1, work2, internalData, solvedCols, 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 = FastMath.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 (FastMath.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]] -= weightedJacobian[i][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 = FastMath.max(parl, lmPar); + } else if (fp < 0) { + paru = FastMath.min(paru, lmPar); + } + + // compute an improved estimate for lmPar + lmPar = FastMath.max(parl, lmPar + correction); + } + + return; + } + + /** + * 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:

    + *
      + *
    • Argonne National Laboratory. MINPACK project. March 1980
      • + *
    • Burton S. Garbow
      • + *
    • Kenneth E. Hillstrom
      • + *
    • Jorge J. More
      • + *
    + *

    Luc Maisonobe did the Java translation.

    + * + * @param qy array containing qTy + * @param diag diagonal matrix + * @param lmDiag diagonal elements associated with lmDir + * @param internalData Data (modified in-place in this method). + * @param solvedCols Number of sloved point. + * @param work work array + */ + private void determineLMDirection(double[] qy, double[] diag, + double[] lmDiag, + InternalData internalData, + int solvedCols, + double[] work) { + final int[] permutation = internalData.permutation; + final double[][] weightedJacobian = internalData.weightedJacobian; + final double[] diagR = internalData.diagR; + + // 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) { + weightedJacobian[i][pj] = weightedJacobian[j][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) { + + final double sin; + final double cos; + double rkk = weightedJacobian[k][pk]; + if (FastMath.abs(rkk) < FastMath.abs(lmDiag[k])) { + final double cotan = rkk / lmDiag[k]; + sin = 1.0 / FastMath.sqrt(1.0 + cotan * cotan); + cos = sin * cotan; + } else { + final double tan = lmDiag[k] / rkk; + cos = 1.0 / FastMath.sqrt(1.0 + tan * tan); + sin = cos * tan; + } + + // compute the modified diagonal element of R and + // the modified element of (Qty,0) + weightedJacobian[k][pk] = cos * rkk + sin * lmDiag[k]; + final 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 = weightedJacobian[i][pk]; + final double temp2 = cos * rik + sin * lmDiag[i]; + lmDiag[i] = -sin * rik + cos * lmDiag[i]; + weightedJacobian[i][pk] = temp2; + } + } + } + + // store the diagonal element of s and restore + // the corresponding diagonal element of R + lmDiag[j] = weightedJacobian[j][permutation[j]]; + weightedJacobian[j][permutation[j]] = 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 += weightedJacobian[i][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.

    + * + * @param jacobian Weighted Jacobian matrix at the current point. + * @param solvedCols Number of solved point. + * @return data used in other methods of this class. + * @throws ConvergenceException if the decomposition cannot be performed. + */ + private InternalData qrDecomposition(RealMatrix jacobian, + int solvedCols) throws ConvergenceException { + // Code in this class assumes that the weighted Jacobian is -(W^(1/2) J), + // hence the multiplication by -1. + final double[][] weightedJacobian = jacobian.scalarMultiply(-1).getData(); + + final int nR = weightedJacobian.length; + final int nC = weightedJacobian[0].length; + + final int[] permutation = new int[nC]; + final double[] diagR = new double[nC]; + final double[] jacNorm = new double[nC]; + final double[] beta = new double[nC]; + + // initializations + for (int k = 0; k < nC; ++k) { + permutation[k] = k; + double norm2 = 0; + for (int i = 0; i < nR; ++i) { + double akk = weightedJacobian[i][k]; + norm2 += akk * akk; + } + jacNorm[k] = FastMath.sqrt(norm2); + } + + // transform the matrix column after column + for (int k = 0; k < nC; ++k) { + + // select the column with the greatest norm on active components + int nextColumn = -1; + double ak2 = Double.NEGATIVE_INFINITY; + for (int i = k; i < nC; ++i) { + double norm2 = 0; + for (int j = k; j < nR; ++j) { + double aki = weightedJacobian[j][permutation[i]]; + norm2 += aki * aki; + } + if (Double.isInfinite(norm2) || Double.isNaN(norm2)) { + throw new ConvergenceException(LocalizedFormats.UNABLE_TO_PERFORM_QR_DECOMPOSITION_ON_JACOBIAN, + nR, nC); + } + if (norm2 > ak2) { + nextColumn = i; + ak2 = norm2; + } + } + if (ak2 <= qrRankingThreshold) { + return new InternalData(weightedJacobian, permutation, k, diagR, jacNorm, beta); + } + int pk = permutation[nextColumn]; + permutation[nextColumn] = permutation[k]; + permutation[k] = pk; + + // choose alpha such that Hk.u = alpha ek + double akk = weightedJacobian[k][pk]; + double alpha = (akk > 0) ? -FastMath.sqrt(ak2) : FastMath.sqrt(ak2); + double betak = 1.0 / (ak2 - akk * alpha); + beta[pk] = betak; + + // transform the current column + diagR[pk] = alpha; + weightedJacobian[k][pk] -= alpha; + + // transform the remaining columns + for (int dk = nC - 1 - k; dk > 0; --dk) { + double gamma = 0; + for (int j = k; j < nR; ++j) { + gamma += weightedJacobian[j][pk] * weightedJacobian[j][permutation[k + dk]]; + } + gamma *= betak; + for (int j = k; j < nR; ++j) { + weightedJacobian[j][permutation[k + dk]] -= gamma * weightedJacobian[j][pk]; + } + } + } + + return new InternalData(weightedJacobian, permutation, solvedCols, diagR, jacNorm, beta); + } + + /** + * Compute the product Qt.y for some Q.R. decomposition. + * + * @param y vector to multiply (will be overwritten with the result) + * @param internalData Data. + */ + private void qTy(double[] y, + InternalData internalData) { + final double[][] weightedJacobian = internalData.weightedJacobian; + final int[] permutation = internalData.permutation; + final double[] beta = internalData.beta; + + final int nR = weightedJacobian.length; + final int nC = weightedJacobian[0].length; + + for (int k = 0; k < nC; ++k) { + int pk = permutation[k]; + double gamma = 0; + for (int i = k; i < nR; ++i) { + gamma += weightedJacobian[i][pk] * y[i]; + } + gamma *= beta[pk]; + for (int i = k; i < nR; ++i) { + y[i] -= gamma * weightedJacobian[i][pk]; + } + } + } +} diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithConvergenceChecker.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithConvergenceChecker.java new file mode 100644 index 000000000..1b523c5b7 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithConvergenceChecker.java @@ -0,0 +1,39 @@ +/* + * 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.math3.fitting.leastsquares; + +import org.apache.commons.math3.optim.ConvergenceChecker; +import org.apache.commons.math3.optim.PointVectorValuePair; + +/** + * Interface for "fluent-API" that advertizes a capability of the optimizer. + * + * @param Concrete optimizer implementation. + * + * @version $Id$ + * @since 3.3 + */ +public interface WithConvergenceChecker { + /** + * Creates a new instance with the specified parameter. + * + * @param checker Convergence checker. + * @return a new optimizer instance with all fields identical to this + * instance except for the given argument. + */ + T withConvergenceChecker(ConvergenceChecker checker); +} diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithMaxEvaluations.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithMaxEvaluations.java new file mode 100644 index 000000000..eec962665 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithMaxEvaluations.java @@ -0,0 +1,36 @@ +/* + * 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.math3.fitting.leastsquares; + +/** + * Interface for "fluent-API" that advertizes a capability of the optimizer. + * + * @param Concrete optimizer implementation. + * + * @version $Id$ + * @since 3.3 + */ +public interface WithMaxEvaluations { + /** + * Creates a new instance with the specified parameter. + * + * @param maxEval Maximum number of evaluations of the model function. + * @return a new optimizer instance with all fields identical to this + * instance except for the given argument. + */ + T withMaxEvaluations(int maxEval); +} diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithMaxIterations.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithMaxIterations.java new file mode 100644 index 000000000..157819822 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithMaxIterations.java @@ -0,0 +1,36 @@ +/* + * 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.math3.fitting.leastsquares; + +/** + * Interface for "fluent-API" that advertizes a capability of the optimizer. + * + * @param Concrete optimizer implementation. + * + * @version $Id$ + * @since 3.3 + */ +public interface WithMaxIterations { + /** + * Creates a new instance with the specified parameter. + * + * @param maxIter Maximum number of iterations. + * @return a new optimizer instance with all fields identical to this + * instance except for the given argument. + */ + T withMaxIterations(int maxIter); +} diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithModelAndJacobian.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithModelAndJacobian.java new file mode 100644 index 000000000..87508f8a0 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithModelAndJacobian.java @@ -0,0 +1,41 @@ +/* + * 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.math3.fitting.leastsquares; + +import org.apache.commons.math3.analysis.MultivariateVectorFunction; +import org.apache.commons.math3.analysis.MultivariateMatrixFunction; + +/** + * Interface for "fluent-API" that advertizes a capability of the optimizer. + * + * @param Concrete optimizer implementation. + * + * @version $Id$ + * @since 3.3 + */ +public interface WithModelAndJacobian { + /** + * Creates a new instance with the specified parameters. + * + * @param model ModelFunction. + * @param jacobian Jacobian of the model function. + * @return a new optimizer instance with all fields identical to this + * instance except for the given arguments. + */ + T withModelAndJacobian(MultivariateVectorFunction model, + MultivariateMatrixFunction jacobian); +} diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithStartPoint.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithStartPoint.java new file mode 100644 index 000000000..cde5e30c6 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithStartPoint.java @@ -0,0 +1,36 @@ +/* + * 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.math3.fitting.leastsquares; + +/** + * Interface for "fluent-API" that advertizes a capability of the optimizer. + * + * @param Concrete optimizer implementation. + * + * @version $Id$ + * @since 3.3 + */ +public interface WithStartPoint { + /** + * Creates a new instance with the specified parameter. + * + * @param start Initial guess for the parameters of the model function. + * @return a new optimizer instance with all fields identical to this + * instance except for the given argument. + */ + T withStartPoint(double[] start); +} diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithTarget.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithTarget.java new file mode 100644 index 000000000..cb0b33f28 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithTarget.java @@ -0,0 +1,36 @@ +/* + * 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.math3.fitting.leastsquares; + +/** + * Interface for "fluent-API" that advertizes a capability of the optimizer. + * + * @param Concrete optimizer implementation. + * + * @version $Id$ + * @since 3.3 + */ +public interface WithTarget { + /** + * Creates a new instance with the specified parameter. + * + * @param target Objective points of the model function. + * @return a new optimizer instance with all fields identical to this + * instance except for the given argument. + */ + T withTarget(double[] target); +} diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithWeight.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithWeight.java new file mode 100644 index 000000000..7d3cd212d --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/WithWeight.java @@ -0,0 +1,38 @@ +/* + * 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.math3.fitting.leastsquares; + +import org.apache.commons.math3.linear.RealMatrix; + +/** + * Interface for "fluent-API" that advertizes a capability of the optimizer. + * + * @param Concrete optimizer implementation. + * + * @version $Id$ + * @since 3.3 + */ +public interface WithWeight { + /** + * Creates a new instance with the specified parameter. + * + * @param weight Weight matrix of the observations. + * @return a new optimizer instance with all fields identical to this + * instance except for the given argument. + */ + T withWeight(RealMatrix weight); +} diff --git a/src/main/java/org/apache/commons/math3/fitting/leastsquares/package-info.java b/src/main/java/org/apache/commons/math3/fitting/leastsquares/package-info.java new file mode 100644 index 000000000..73c695663 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/leastsquares/package-info.java @@ -0,0 +1,37 @@ +/* + * 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. + */ + +/** + * This package provides algorithms that minimize the residuals + * between observations and model values. + * The {@link org.apache.commons.math3.fitting.leastsquares.AbstractLeastSquaresOptimizer + * non-linear least-squares optimizers} minimize the distance (called cost or + * χ2) between model and observations. + * + *
    + * Algorithms in this category need access to a model function + * (represented by a {@link org.apache.commons.math3.analysis.MultivariateVectorFunction + * MultivariateVectorFunction}). + * Such a model predicts a set of values which the algorithm tries to match with a set + * of given set of {@link WithTarget observed values}. + *
    + * The algorithms implemented in this package also require that the user specifies the + * Jacobian matrix of the model (represented by a + * {@link org.apache.commons.math3.analysis.MultivariateMatrixFunction + * MultivariateMatrixFunction}). + */ +package org.apache.commons.math3.fitting.leastsquares; diff --git a/src/main/java/org/apache/commons/math3/optim/BaseOptimizer.java b/src/main/java/org/apache/commons/math3/optim/BaseOptimizer.java index 9aff1df26..3a98feeb1 100644 --- a/src/main/java/org/apache/commons/math3/optim/BaseOptimizer.java +++ b/src/main/java/org/apache/commons/math3/optim/BaseOptimizer.java @@ -45,10 +45,21 @@ public abstract class BaseOptimizer { * @param checker Convergence checker. */ protected BaseOptimizer(ConvergenceChecker checker) { + this(checker, 0, Integer.MAX_VALUE); + } + + /** + * @param checker Convergence checker. + * @param maxEval Maximum number of objective function evaluations. + * @param maxIter Maximum number of algorithm iterations. + */ + protected BaseOptimizer(ConvergenceChecker checker, + int maxEval, + int maxIter) { this.checker = checker; - evaluations = new Incrementor(0, new MaxEvalCallback()); - iterations = new Incrementor(Integer.MAX_VALUE, new MaxIterCallback()); + evaluations = new Incrementor(maxEval, new MaxEvalCallback()); + iterations = new Incrementor(maxIter, new MaxIterCallback()); } /** @@ -143,6 +154,25 @@ public abstract class BaseOptimizer { return doOptimize(); } + /** + * Performs the optimization. + * + * @return a point/value pair that satifies the convergence criteria. + * @throws TooManyEvaluationsException if the maximal number of + * evaluations is exceeded. + * @throws TooManyIterationsException if the maximal number of + * iterations is exceeded. + */ + public PAIR optimize() + throws TooManyEvaluationsException, + TooManyIterationsException { + // Reset counters. + evaluations.resetCount(); + iterations.resetCount(); + // Perform optimization. + return doOptimize(); + } + /** * Performs the bulk of the optimization algorithm. * diff --git a/src/test/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizerAbstractTest.java b/src/test/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizerAbstractTest.java new file mode 100644 index 000000000..912b67ba6 --- /dev/null +++ b/src/test/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizerAbstractTest.java @@ -0,0 +1,643 @@ +/* + * 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.math3.fitting.leastsquares; + +import java.io.IOException; +import java.io.Serializable; +import java.util.Arrays; +import org.apache.commons.math3.analysis.MultivariateVectorFunction; +import org.apache.commons.math3.analysis.MultivariateMatrixFunction; +import org.apache.commons.math3.exception.ConvergenceException; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.geometry.euclidean.twod.Vector2D; +import org.apache.commons.math3.linear.BlockRealMatrix; +import org.apache.commons.math3.linear.RealMatrix; +import org.apache.commons.math3.linear.DiagonalMatrix; +import org.apache.commons.math3.optim.PointVectorValuePair; +import org.apache.commons.math3.util.FastMath; +import org.junit.Assert; +import org.junit.Test; + +/** + * Some of the unit tests are re-implementations of the MINPACK + * file17 and + * file22 test files. + * The redistribution policy for MINPACK is available + * here. + * + * Concrete implementation of an optimizer. + * + * @version $Id$ + */ +public abstract class AbstractLeastSquaresOptimizerAbstractTest & + WithWeight & + WithModelAndJacobian & + WithConvergenceChecker & + WithStartPoint & + WithMaxIterations & + WithMaxEvaluations> { + /** + * @return a concrete optimizer. + */ + public abstract T createOptimizer(); + + /** + * @return the default number of allowed iterations (which will be + * used when not specified otherwise). + */ + public abstract int getMaxIterations(); + + @Test + public void testGetIterations() { + T optim = createOptimizer() + .withMaxEvaluations(100) + .withMaxIterations(getMaxIterations()) + .withTarget(new double[] { 1 }) + .withWeight(new DiagonalMatrix(new double[] { 1 })) + .withStartPoint(new double[] { 3 }) + .withModelAndJacobian(new MultivariateVectorFunction() { + public double[] value(double[] point) { + return new double[] { + FastMath.pow(point[0], 4) + }; + }}, + new MultivariateMatrixFunction() { + public double[][] value(double[] point) { + return new double[][] { + { 0.25 * FastMath.pow(point[0], 3) } + }; + } + }); + + optim.optimize(); + Assert.assertTrue(optim.getIterations() > 0); + } + + @Test + public void testTrivial() { + LinearProblem problem + = new LinearProblem(new double[][] { { 2 } }, + new double[] { 3 }); + T optimizer = createOptimizer() + .withMaxEvaluations(100) + .withMaxIterations(getMaxIterations()) + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(problem.getTarget()) + .withWeight(new DiagonalMatrix(new double[] { 1 })) + .withStartPoint(new double[] { 0 }); + + PointVectorValuePair optimum = optimizer.optimize(); + + Assert.assertEquals(0, optimizer.computeRMS(optimum.getPoint()), 1e-10); + Assert.assertEquals(1.5, optimum.getPoint()[0], 1e-10); + Assert.assertEquals(3.0, optimum.getValue()[0], 1e-10); + } + + @Test + public void testQRColumnsPermutation() { + LinearProblem problem + = new LinearProblem(new double[][] { { 1, -1 }, { 0, 2 }, { 1, -2 } }, + new double[] { 4, 6, 1 }); + + T optimizer = createOptimizer() + .withMaxEvaluations(100) + .withMaxIterations(getMaxIterations()) + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(problem.getTarget()) + .withWeight(new DiagonalMatrix(new double[] { 1, 1, 1 })) + .withStartPoint(new double[] { 0, 0 }); + + PointVectorValuePair optimum = optimizer.optimize(); + + Assert.assertEquals(0, optimizer.computeRMS(optimum.getPoint()), 1e-10); + Assert.assertEquals(7, optimum.getPoint()[0], 1e-10); + Assert.assertEquals(3, optimum.getPoint()[1], 1e-10); + Assert.assertEquals(4, optimum.getValue()[0], 1e-10); + Assert.assertEquals(6, optimum.getValue()[1], 1e-10); + Assert.assertEquals(1, optimum.getValue()[2], 1e-10); + } + + @Test + public void testNoDependency() { + LinearProblem problem = new LinearProblem(new double[][] { + { 2, 0, 0, 0, 0, 0 }, + { 0, 2, 0, 0, 0, 0 }, + { 0, 0, 2, 0, 0, 0 }, + { 0, 0, 0, 2, 0, 0 }, + { 0, 0, 0, 0, 2, 0 }, + { 0, 0, 0, 0, 0, 2 } + }, new double[] { 0, 1.1, 2.2, 3.3, 4.4, 5.5 }); + T optimizer = createOptimizer() + .withMaxEvaluations(100) + .withMaxIterations(getMaxIterations()) + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(problem.getTarget()) + .withWeight(new DiagonalMatrix(new double[] { 1, 1, 1, 1, 1, 1 })) + .withStartPoint(new double[] { 0, 0, 0, 0, 0, 0 }); + + double[] optimum = optimizer.optimize().getPoint(); + Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10); + for (int i = 0; i < problem.target.length; ++i) { + Assert.assertEquals(0.55 * i, optimum[i], 1e-10); + } + } + + @Test + public void testOneSet() { + LinearProblem problem = new LinearProblem(new double[][] { + { 1, 0, 0 }, + { -1, 1, 0 }, + { 0, -1, 1 } + }, new double[] { 1, 1, 1}); + + T optimizer = createOptimizer() + .withMaxEvaluations(100) + .withMaxIterations(getMaxIterations()) + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(problem.getTarget()) + .withWeight(new DiagonalMatrix(new double[] { 1, 1, 1 })) + .withStartPoint(new double[] { 0, 0, 0 }); + + double[] optimum = optimizer.optimize().getPoint(); + Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10); + Assert.assertEquals(1, optimum[0], 1e-10); + Assert.assertEquals(2, optimum[1], 1e-10); + Assert.assertEquals(3, optimum[2], 1e-10); + } + + @Test + public void testTwoSets() { + double epsilon = 1e-7; + LinearProblem problem = new LinearProblem(new double[][] { + { 2, 1, 0, 4, 0, 0 }, + { -4, -2, 3, -7, 0, 0 }, + { 4, 1, -2, 8, 0, 0 }, + { 0, -3, -12, -1, 0, 0 }, + { 0, 0, 0, 0, epsilon, 1 }, + { 0, 0, 0, 0, 1, 1 } + }, new double[] { 2, -9, 2, 2, 1 + epsilon * epsilon, 2}); + + T optimizer = createOptimizer() + .withMaxEvaluations(100) + .withMaxIterations(getMaxIterations()) + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(problem.getTarget()) + .withWeight(new DiagonalMatrix(new double[] { 1, 1, 1, 1, 1, 1 })) + .withStartPoint(new double[] { 0, 0, 0, 0, 0, 0 }); + + double[] optimum = optimizer.optimize().getPoint(); + + Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10); + Assert.assertEquals(3, optimum[0], 1e-10); + Assert.assertEquals(4, optimum[1], 1e-10); + Assert.assertEquals(-1, optimum[2], 1e-10); + Assert.assertEquals(-2, optimum[3], 1e-10); + Assert.assertEquals(1 + epsilon, optimum[4], 1e-10); + Assert.assertEquals(1 - epsilon, optimum[5], 1e-10); + } + + @Test(expected=ConvergenceException.class) + public void testNonInvertible() throws Exception { + LinearProblem problem = new LinearProblem(new double[][] { + { 1, 2, -3 }, + { 2, 1, 3 }, + { -3, 0, -9 } + }, new double[] { 1, 1, 1 }); + + T optimizer = createOptimizer() + .withMaxEvaluations(100) + .withMaxIterations(getMaxIterations()) + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(problem.getTarget()) + .withWeight(new DiagonalMatrix(new double[] { 1, 1, 1 })) + .withStartPoint(new double[] { 0, 0, 0 }); + + optimizer.optimize(); + } + + @Test + public void testIllConditioned() { + LinearProblem problem1 = new LinearProblem(new double[][] { + { 10, 7, 8, 7 }, + { 7, 5, 6, 5 }, + { 8, 6, 10, 9 }, + { 7, 5, 9, 10 } + }, new double[] { 32, 23, 33, 31 }); + T optimizer = createOptimizer() + .withMaxEvaluations(100) + .withMaxIterations(getMaxIterations()) + .withModelAndJacobian(problem1.getModelFunction(), + problem1.getModelFunctionJacobian()) + .withTarget(problem1.getTarget()) + .withWeight(new DiagonalMatrix(new double[] { 1, 1, 1, 1 })) + .withStartPoint(new double[] { 0, 1, 2, 3 }); + + double[] optimum = optimizer.optimize().getPoint(); + + Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10); + Assert.assertEquals(1, optimum[0], 1e-10); + Assert.assertEquals(1, optimum[1], 1e-10); + Assert.assertEquals(1, optimum[2], 1e-10); + Assert.assertEquals(1, optimum[3], 1e-10); + + LinearProblem problem2 = new LinearProblem(new double[][] { + { 10.00, 7.00, 8.10, 7.20 }, + { 7.08, 5.04, 6.00, 5.00 }, + { 8.00, 5.98, 9.89, 9.00 }, + { 6.99, 4.99, 9.00, 9.98 } + }, new double[] { 32, 23, 33, 31 }); + + optimizer = optimizer + .withModelAndJacobian(problem2.getModelFunction(), + problem2.getModelFunctionJacobian()) + .withTarget(problem2.getTarget()); + + optimum = optimizer.optimize().getPoint(); + + Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10); + Assert.assertEquals(-81, optimum[0], 1e-8); + Assert.assertEquals(137, optimum[1], 1e-8); + Assert.assertEquals(-34, optimum[2], 1e-8); + Assert.assertEquals( 22, optimum[3], 1e-8); + } + + @Test + public void testMoreEstimatedParametersSimple() { + LinearProblem problem = new LinearProblem(new double[][] { + { 3, 2, 0, 0 }, + { 0, 1, -1, 1 }, + { 2, 0, 1, 0 } + }, new double[] { 7, 3, 5 }); + + T optimizer = createOptimizer() + .withMaxEvaluations(100) + .withMaxIterations(getMaxIterations()) + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(problem.getTarget()) + .withWeight(new DiagonalMatrix(new double[] { 1, 1, 1 })) + .withStartPoint(new double[] { 7, 6, 5, 4 }); + + double[] optimum = optimizer.optimize().getPoint(); + Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10); + } + + @Test + public void testMoreEstimatedParametersUnsorted() { + LinearProblem problem = new LinearProblem(new double[][] { + { 1, 1, 0, 0, 0, 0 }, + { 0, 0, 1, 1, 1, 0 }, + { 0, 0, 0, 0, 1, -1 }, + { 0, 0, -1, 1, 0, 1 }, + { 0, 0, 0, -1, 1, 0 } + }, new double[] { 3, 12, -1, 7, 1 }); + + T optimizer = createOptimizer() + .withMaxEvaluations(100) + .withMaxIterations(getMaxIterations()) + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(problem.getTarget()) + .withWeight(new DiagonalMatrix(new double[] { 1, 1, 1, 1, 1 })) + .withStartPoint(new double[] { 2, 2, 2, 2, 2, 2 }); + + double[] optimum = optimizer.optimize().getPoint(); + + Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10); + Assert.assertEquals(3, optimum[2], 1e-10); + Assert.assertEquals(4, optimum[3], 1e-10); + Assert.assertEquals(5, optimum[4], 1e-10); + Assert.assertEquals(6, optimum[5], 1e-10); + } + + @Test + public void testRedundantEquations() { + LinearProblem problem = new LinearProblem(new double[][] { + { 1, 1 }, + { 1, -1 }, + { 1, 3 } + }, new double[] { 3, 1, 5 }); + + T optimizer = createOptimizer() + .withMaxEvaluations(100) + .withMaxIterations(getMaxIterations()) + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(problem.getTarget()) + .withWeight(new DiagonalMatrix(new double[] { 1, 1, 1 })) + .withStartPoint(new double[] { 1, 1 }); + + double[] optimum = optimizer.optimize().getPoint(); + + Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10); + Assert.assertEquals(2, optimum[0], 1e-10); + Assert.assertEquals(1, optimum[1], 1e-10); + } + + @Test + public void testInconsistentEquations() { + LinearProblem problem = new LinearProblem(new double[][] { + { 1, 1 }, + { 1, -1 }, + { 1, 3 } + }, new double[] { 3, 1, 4 }); + + T optimizer = createOptimizer() + .withMaxEvaluations(100) + .withMaxIterations(getMaxIterations()) + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(problem.getTarget()) + .withWeight(new DiagonalMatrix(new double[] { 1, 1, 1 })) + .withStartPoint(new double[] { 1, 1 }); + + double[] optimum = optimizer.optimize().getPoint(); + + Assert.assertTrue(optimizer.computeRMS(optimum) > 0.1); + } + + @Test(expected=DimensionMismatchException.class) + public void testInconsistentSizes1() { + LinearProblem problem + = new LinearProblem(new double[][] { { 1, 0 }, + { 0, 1 } }, + new double[] { -1, 1 }); + T optimizer = createOptimizer() + .withMaxEvaluations(100) + .withMaxIterations(getMaxIterations()) + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(problem.getTarget()) + .withWeight(new DiagonalMatrix(new double[] { 1, 1 })) + .withStartPoint(new double[] { 0, 0 }); + + double[] optimum = optimizer.optimize().getPoint(); + + Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10); + Assert.assertEquals(-1, optimum[0], 1e-10); + Assert.assertEquals(1, optimum[1], 1e-10); + + optimizer.withWeight(new DiagonalMatrix(new double[] { 1 })).optimize(); + } + + @Test(expected=DimensionMismatchException.class) + public void testInconsistentSizes2() { + LinearProblem problem + = new LinearProblem(new double[][] { { 1, 0 }, { 0, 1 } }, + new double[] { -1, 1 }); + + T optimizer = createOptimizer() + .withMaxEvaluations(100) + .withMaxIterations(getMaxIterations()) + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(problem.getTarget()) + .withWeight(new DiagonalMatrix(new double[] { 1, 1 })) + .withStartPoint(new double[] { 0, 0 }); + + double[] optimum = optimizer.optimize().getPoint(); + + Assert.assertEquals(0, optimizer.computeRMS(optimum), 1e-10); + Assert.assertEquals(-1, optimum[0], 1e-10); + Assert.assertEquals(1, optimum[1], 1e-10); + + optimizer + .withTarget(new double[] { 1 }) + .withWeight(new DiagonalMatrix(new double[] { 1 })) + .optimize(); + } + + @Test + public void testCircleFitting() { + CircleVectorial circle = new CircleVectorial(); + circle.addPoint( 30, 68); + circle.addPoint( 50, -6); + circle.addPoint(110, -20); + circle.addPoint( 35, 15); + circle.addPoint( 45, 97); + + T optimizer = createOptimizer() + .withMaxEvaluations(100) + .withMaxIterations(getMaxIterations()) + .withModelAndJacobian(circle.getModelFunction(), + circle.getModelFunctionJacobian()) + .withTarget(new double[] { 0, 0, 0, 0, 0 }) + .withWeight(new DiagonalMatrix(new double[] { 1, 1, 1, 1, 1 })) + .withStartPoint(new double[] { 98.680, 47.345 }); + + double[] optimum = optimizer.optimize().getPoint(); + Assert.assertTrue(optimizer.getEvaluations() < 10); + + double rms = optimizer.computeRMS(optimum); + Assert.assertEquals(1.768262623567235, FastMath.sqrt(circle.getN()) * rms, 1e-10); + + Vector2D center = new Vector2D(optimum[0], optimum[1]); + Assert.assertEquals(69.96016176931406, circle.getRadius(center), 1e-6); + Assert.assertEquals(96.07590211815305, center.getX(), 1e-6); + Assert.assertEquals(48.13516790438953, center.getY(), 1e-6); + + double[][] cov = optimizer.computeCovariances(optimum, 1e-14); + Assert.assertEquals(1.839, cov[0][0], 0.001); + Assert.assertEquals(0.731, cov[0][1], 0.001); + Assert.assertEquals(cov[0][1], cov[1][0], 1e-14); + Assert.assertEquals(0.786, cov[1][1], 0.001); + + // add perfect measurements and check errors are reduced + double r = circle.getRadius(center); + for (double d= 0; d < 2 * FastMath.PI; d += 0.01) { + circle.addPoint(center.getX() + r * FastMath.cos(d), center.getY() + r * FastMath.sin(d)); + } + + double[] target = new double[circle.getN()]; + Arrays.fill(target, 0); + double[] weights = new double[circle.getN()]; + Arrays.fill(weights, 2); + optimizer = optimizer.withTarget(target).withWeight(new DiagonalMatrix(weights)); + optimum = optimizer.optimize().getPoint(); + + cov = optimizer.computeCovariances(optimum, 1e-14); + Assert.assertEquals(0.0016, cov[0][0], 0.001); + Assert.assertEquals(3.2e-7, cov[0][1], 1e-9); + Assert.assertEquals(cov[0][1], cov[1][0], 1e-14); + Assert.assertEquals(0.0016, cov[1][1], 0.001); + } + + @Test + public void testCircleFittingBadInit() { + CircleVectorial circle = new CircleVectorial(); + double[][] points = circlePoints; + double[] target = new double[points.length]; + Arrays.fill(target, 0); + double[] weights = new double[points.length]; + Arrays.fill(weights, 2); + for (int i = 0; i < points.length; ++i) { + circle.addPoint(points[i][0], points[i][1]); + } + T optimizer = createOptimizer() + .withMaxEvaluations(100) + .withMaxIterations(getMaxIterations()) + .withModelAndJacobian(circle.getModelFunction(), + circle.getModelFunctionJacobian()) + .withTarget(target) + .withWeight(new DiagonalMatrix(weights)) + .withStartPoint(new double[] { -12, -12 }); + + double[] optimum = optimizer.optimize().getPoint(); + + Vector2D center = new Vector2D(optimum[0], optimum[1]); + Assert.assertTrue(optimizer.getEvaluations() < 25); + Assert.assertEquals( 0.043, optimizer.computeRMS(optimum), 1e-3); + Assert.assertEquals( 0.292235, circle.getRadius(center), 1e-6); + Assert.assertEquals(-0.151738, center.getX(), 1e-6); + Assert.assertEquals( 0.2075001, center.getY(), 1e-6); + } + + @Test + public void testCircleFittingGoodInit() { + CircleVectorial circle = new CircleVectorial(); + double[][] points = circlePoints; + double[] target = new double[points.length]; + Arrays.fill(target, 0); + double[] weights = new double[points.length]; + Arrays.fill(weights, 2); + for (int i = 0; i < points.length; ++i) { + circle.addPoint(points[i][0], points[i][1]); + } + T optimizer = createOptimizer() + .withMaxEvaluations(100) + .withMaxIterations(getMaxIterations()) + .withModelAndJacobian(circle.getModelFunction(), + circle.getModelFunctionJacobian()) + .withTarget(target) + .withWeight(new DiagonalMatrix(weights)) + .withStartPoint(new double[] { 0, 0 }); + + double[] optimum = optimizer.optimize().getPoint(); + + Assert.assertEquals(-0.1517383071957963, optimum[0], 1e-6); + Assert.assertEquals(0.2074999736353867, optimum[1], 1e-6); + Assert.assertEquals(0.04268731682389561, optimizer.computeRMS(optimum), 1e-8); + } + + private final double[][] circlePoints = 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} + }; + + public void doTestStRD(final StatisticalReferenceDataset dataset, + final double errParams, + final double errParamsSd) { + final double[] w = new double[dataset.getNumObservations()]; + Arrays.fill(w, 1); + + final double[][] data = dataset.getData(); + final double[] initial = dataset.getStartingPoint(0); + final StatisticalReferenceDataset.LeastSquaresProblem problem = dataset.getLeastSquaresProblem(); + + final T optimizer = createOptimizer() + .withMaxEvaluations(100) + .withMaxIterations(getMaxIterations()) + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(data[1]) + .withWeight(new DiagonalMatrix(w)) + .withStartPoint(initial); + + final double[] actual = optimizer.optimize().getPoint(); + for (int i = 0; i < actual.length; i++) { + double expected = dataset.getParameter(i); + double delta = FastMath.abs(errParams * expected); + Assert.assertEquals(dataset.getName() + ", param #" + i, + expected, actual[i], delta); + } + } + + @Test + public void testKirby2() throws IOException { + doTestStRD(StatisticalReferenceDatasetFactory.createKirby2(), 1E-7, 1E-7); + } + + @Test + public void testHahn1() throws IOException { + doTestStRD(StatisticalReferenceDatasetFactory.createHahn1(), 1E-7, 1E-4); + } + + static class LinearProblem { + private final RealMatrix factors; + private final double[] target; + + public LinearProblem(double[][] factors, double[] target) { + this.factors = new BlockRealMatrix(factors); + this.target = target; + } + + public double[] getTarget() { + return target; + } + + public MultivariateVectorFunction getModelFunction() { + return new MultivariateVectorFunction() { + public double[] value(double[] params) { + return factors.operate(params); + } + }; + } + + public MultivariateMatrixFunction getModelFunctionJacobian() { + return new MultivariateMatrixFunction() { + public double[][] value(double[] params) { + return factors.getData(); + } + }; + } + } +} diff --git a/src/test/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizerTest.java b/src/test/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizerTest.java new file mode 100644 index 000000000..e971677ff --- /dev/null +++ b/src/test/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizerTest.java @@ -0,0 +1,109 @@ +/* + * 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.math3.fitting.leastsquares; + +import java.io.IOException; +import java.util.Arrays; +import org.apache.commons.math3.optim.PointVectorValuePair; +import org.apache.commons.math3.linear.DiagonalMatrix; +import org.apache.commons.math3.util.FastMath; +import org.junit.Test; +import org.junit.Assert; +/** + * The only features tested here are utility methods defined + * in {@link AbstractLeastSquaresOptimizer} that compute the + * chi-square and parameters standard-deviations. + */ +public class AbstractLeastSquaresOptimizerTest { + @Test + public void testComputeCost() throws IOException { + final StatisticalReferenceDataset dataset + = StatisticalReferenceDatasetFactory.createKirby2(); + final double[] a = dataset.getParameters(); + final double[] y = dataset.getData()[1]; + final double[] w = new double[y.length]; + Arrays.fill(w, 1d); + + StatisticalReferenceDataset.LeastSquaresProblem problem + = dataset.getLeastSquaresProblem(); + + final LevenbergMarquardtOptimizer optim = LevenbergMarquardtOptimizer.create() + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(y) + .withWeight(new DiagonalMatrix(w)) + .withStartPoint(a); + + final double expected = dataset.getResidualSumOfSquares(); + final double cost = optim.computeCost(optim.computeResiduals(optim.getModel().value(optim.getStart()))); + final double actual = cost * cost; + Assert.assertEquals(dataset.getName(), expected, actual, 1e-11 * expected); + } + + @Test + public void testComputeRMS() throws IOException { + final StatisticalReferenceDataset dataset + = StatisticalReferenceDatasetFactory.createKirby2(); + final double[] a = dataset.getParameters(); + final double[] y = dataset.getData()[1]; + final double[] w = new double[y.length]; + Arrays.fill(w, 1d); + + StatisticalReferenceDataset.LeastSquaresProblem problem + = dataset.getLeastSquaresProblem(); + + final LevenbergMarquardtOptimizer optim = LevenbergMarquardtOptimizer.create() + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(y) + .withWeight(new DiagonalMatrix(w)) + .withStartPoint(a); + + final double expected = FastMath.sqrt(dataset.getResidualSumOfSquares() / + dataset.getNumObservations()); + final double actual = optim.computeRMS(optim.getStart()); + Assert.assertEquals(dataset.getName(), expected, actual, 1e-11 * expected); + } + + @Test + public void testComputeSigma() throws IOException { + final StatisticalReferenceDataset dataset + = StatisticalReferenceDatasetFactory.createKirby2(); + final double[] a = dataset.getParameters(); + final double[] y = dataset.getData()[1]; + final double[] w = new double[y.length]; + Arrays.fill(w, 1d); + + StatisticalReferenceDataset.LeastSquaresProblem problem + = dataset.getLeastSquaresProblem(); + + final LevenbergMarquardtOptimizer optim = LevenbergMarquardtOptimizer.create() + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(y) + .withWeight(new DiagonalMatrix(w)) + .withStartPoint(a); + + final double[] expected = dataset.getParametersStandardDeviations(); + + final double cost = optim.computeCost(optim.computeResiduals(optim.getModel().value(optim.getStart()))); + final double[] sig = optim.computeSigma(optim.getStart(), 1e-14); + final int dof = y.length - a.length; + for (int i = 0; i < sig.length; i++) { + final double actual = FastMath.sqrt(cost * cost / dof) * sig[i]; + Assert.assertEquals(dataset.getName() + ", parameter #" + i, + expected[i], actual, 1e-6 * expected[i]); + } + } +} diff --git a/src/test/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizerTestValidation.java b/src/test/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizerTestValidation.java new file mode 100644 index 000000000..5c298e26d --- /dev/null +++ b/src/test/java/org/apache/commons/math3/fitting/leastsquares/AbstractLeastSquaresOptimizerTestValidation.java @@ -0,0 +1,308 @@ +/* + * 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.math3.fitting.leastsquares; + +import java.util.Arrays; +import java.util.List; +import java.util.ArrayList; +import java.awt.geom.Point2D; +import org.apache.commons.math3.optim.PointVectorValuePair; +import org.apache.commons.math3.stat.descriptive.SummaryStatistics; +import org.apache.commons.math3.stat.descriptive.StatisticalSummary; +import org.apache.commons.math3.linear.DiagonalMatrix; +import org.apache.commons.math3.util.FastMath; +import org.junit.Test; +import org.junit.Assert; + +/** + * This class demonstrates the main functionality of the + * {@link AbstractLeastSquaresOptimizer}, common to the + * optimizer implementations in package + * {@link org.apache.commons.math3.fitting.leastsquares}. + *
    + * Not enabled by default, as the class name does not end with "Test". + *
    + * Invoke by running + * 
    
+ *  mvn test -Dtest=AbstractLeastSquaresOptimizerTestValidation
+ *
    + * or by running + * 
    
+ *  mvn test -Dtest=AbstractLeastSquaresOptimizerTestValidation -DargLine="-DmcRuns=1234 -server"
+ *
    + */ +public class AbstractLeastSquaresOptimizerTestValidation { + /** Number of runs. */ + private static final int MONTE_CARLO_RUNS = Integer.parseInt(System.getProperty("mcRuns", + "100")); + + /** + * Using a Monte-Carlo procedure, this test checks the error estimations + * as provided by the square-root of the diagonal elements of the + * covariance matrix. + *
    + * The test generates sets of observations, each sampled from + * a Gaussian distribution. + *
    + * The optimization problem solved is defined in class + * {@link StraightLineProblem}. + *
    + * The output (on stdout) will be a table summarizing the distribution + * of parameters generated by the Monte-Carlo process and by the direct + * estimation provided by the diagonal elements of the covariance matrix. + */ + @Test + public void testParametersErrorMonteCarloObservations() { + // Error on the observations. + final double yError = 15; + + // True values of the parameters. + final double slope = 123.456; + final double offset = -98.765; + + // Samples generator. + final RandomStraightLinePointGenerator lineGenerator + = new RandomStraightLinePointGenerator(slope, offset, + yError, + -1e3, 1e4, + 138577L); + + // Number of observations. + final int numObs = 100; // XXX Should be a command-line option. + // number of parameters. + final int numParams = 2; + + // Parameters found for each of Monte-Carlo run. + final SummaryStatistics[] paramsFoundByDirectSolution = new SummaryStatistics[numParams]; + // Sigma estimations (square-root of the diagonal elements of the + // covariance matrix), for each Monte-Carlo run. + final SummaryStatistics[] sigmaEstimate = new SummaryStatistics[numParams]; + + // Initialize statistics accumulators. + for (int i = 0; i < numParams; i++) { + paramsFoundByDirectSolution[i] = new SummaryStatistics(); + sigmaEstimate[i] = new SummaryStatistics(); + } + + final double[] init = { slope, offset }; + + // Monte-Carlo (generates many sets of observations). + final int mcRepeat = MONTE_CARLO_RUNS; + int mcCount = 0; + while (mcCount < mcRepeat) { + // Observations. + final Point2D.Double[] obs = lineGenerator.generate(numObs); + + final StraightLineProblem problem = new StraightLineProblem(yError); + for (int i = 0; i < numObs; i++) { + final Point2D.Double p = obs[i]; + problem.addPoint(p.x, p.y); + } + + // Direct solution (using simple regression). + final double[] regress = problem.solve(); + + // Estimation of the standard deviation (diagonal elements of the + // covariance matrix). + // Dummy optimizer (to compute the covariance matrix). + final AbstractLeastSquaresOptimizer optim = LevenbergMarquardtOptimizer.create() + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(problem.target()) + .withWeight(new DiagonalMatrix(problem.weight())); + + final double[] sigma = optim.computeSigma(init, 1e-14); + + // Accumulate statistics. + for (int i = 0; i < numParams; i++) { + paramsFoundByDirectSolution[i].addValue(regress[i]); + sigmaEstimate[i].addValue(sigma[i]); + } + + // Next Monte-Carlo. + ++mcCount; + } + + // Print statistics. + final String line = "--------------------------------------------------------------"; + System.out.println(" True value Mean Std deviation"); + for (int i = 0; i < numParams; i++) { + System.out.println(line); + System.out.println("Parameter #" + i); + + StatisticalSummary s = paramsFoundByDirectSolution[i].getSummary(); + System.out.printf(" %+.6e %+.6e %+.6e\n", + init[i], + s.getMean(), + s.getStandardDeviation()); + + s = sigmaEstimate[i].getSummary(); + System.out.printf("sigma: %+.6e (%+.6e)\n", + s.getMean(), + s.getStandardDeviation()); + } + System.out.println(line); + + // Check the error estimation. + for (int i = 0; i < numParams; i++) { + Assert.assertEquals(paramsFoundByDirectSolution[i].getSummary().getStandardDeviation(), + sigmaEstimate[i].getSummary().getMean(), + 8e-2); + } + } + + /** + * In this test, the set of observations is fixed. + * Using a Monte-Carlo procedure, it generates sets of parameters, + * and determine the parameter change that will result in the + * normalized chi-square becoming larger by one than the value from + * the best fit solution. + *
    + * The optimization problem solved is defined in class + * {@link StraightLineProblem}. + *
    + * The output (on stdout) will be a list of lines containing: + *
      + *
    • slope of the straight line,
      • + *
    • intercept of the straight line,
      • + *
    • chi-square of the solution defined by the above two values.
      • + *
    + * The output is separated into two blocks (with a blank line between + * them); the first block will contain all parameter sets for which + * {@code chi2 < chi2_b + 1} + * and the second block, all sets for which + * {@code chi2 >= chi2_b + 1} + * where {@code chi2_b} is the lowest chi-square (corresponding to the + * best solution). + */ + @Test + public void testParametersErrorMonteCarloParameters() { + // Error on the observations. + final double yError = 15; + + // True values of the parameters. + final double slope = 123.456; + final double offset = -98.765; + + // Samples generator. + final RandomStraightLinePointGenerator lineGenerator + = new RandomStraightLinePointGenerator(slope, offset, + yError, + -1e3, 1e4, + 13839013L); + + // Number of observations. + final int numObs = 10; + // number of parameters. + final int numParams = 2; + + // Create a single set of observations. + final Point2D.Double[] obs = lineGenerator.generate(numObs); + + final StraightLineProblem problem = new StraightLineProblem(yError); + for (int i = 0; i < numObs; i++) { + final Point2D.Double p = obs[i]; + problem.addPoint(p.x, p.y); + } + + // Direct solution (using simple regression). + final double[] regress = problem.solve(); + + // Dummy optimizer (to compute the chi-square). + final AbstractLeastSquaresOptimizer optim = LevenbergMarquardtOptimizer.create() + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(problem.target()) + .withWeight(new DiagonalMatrix(problem.weight())); + + final double[] init = { slope, offset }; + // Get chi-square of the best parameters set for the given set of + // observations. + final double bestChi2N = getChi2N(optim, problem, regress); + final double[] sigma = optim.computeSigma(regress, 1e-14); + + // Monte-Carlo (generates a grid of parameters). + final int mcRepeat = MONTE_CARLO_RUNS; + final int gridSize = (int) FastMath.sqrt(mcRepeat); + + // Parameters found for each of Monte-Carlo run. + // Index 0 = slope + // Index 1 = offset + // Index 2 = normalized chi2 + final List paramsAndChi2 = new ArrayList(gridSize * gridSize); + + final double slopeRange = 10 * sigma[0]; + final double offsetRange = 10 * sigma[1]; + final double minSlope = slope - 0.5 * slopeRange; + final double minOffset = offset - 0.5 * offsetRange; + final double deltaSlope = slopeRange/ gridSize; + final double deltaOffset = offsetRange / gridSize; + for (int i = 0; i < gridSize; i++) { + final double s = minSlope + i * deltaSlope; + for (int j = 0; j < gridSize; j++) { + final double o = minOffset + j * deltaOffset; + final double chi2N = getChi2N(optim, problem, new double[] {s, o}); + + paramsAndChi2.add(new double[] {s, o, chi2N}); + } + } + + // Output (for use with "gnuplot"). + + // Some info. + + // For plotting separately sets of parameters that have a large chi2. + final double chi2NPlusOne = bestChi2N + 1; + int numLarger = 0; + + final String lineFmt = "%+.10e %+.10e %.8e\n"; + + // Point with smallest chi-square. + System.out.printf(lineFmt, regress[0], regress[1], bestChi2N); + System.out.println(); // Empty line. + + // Points within the confidence interval. + for (double[] d : paramsAndChi2) { + if (d[2] <= chi2NPlusOne) { + System.out.printf(lineFmt, d[0], d[1], d[2]); + } + } + System.out.println(); // Empty line. + + // Points outside the confidence interval. + for (double[] d : paramsAndChi2) { + if (d[2] > chi2NPlusOne) { + ++numLarger; + System.out.printf(lineFmt, d[0], d[1], d[2]); + } + } + System.out.println(); // Empty line. + + System.out.println("# sigma=" + Arrays.toString(sigma)); + System.out.println("# " + numLarger + " sets filtered out"); + } + + /** + * @return the normalized chi-square. + */ + private double getChi2N(AbstractLeastSquaresOptimizer optim, + StraightLineProblem problem, + double[] params) { + final double[] t = problem.target(); + final double[] w = problem.weight(); + final double cost = optim.computeCost(optim.computeResiduals(optim.getModel().value(params))); + return cost * cost / (t.length - params.length); + } +} + diff --git a/src/test/java/org/apache/commons/math3/fitting/leastsquares/CircleProblem.java b/src/test/java/org/apache/commons/math3/fitting/leastsquares/CircleProblem.java new file mode 100644 index 000000000..ed637d43b --- /dev/null +++ b/src/test/java/org/apache/commons/math3/fitting/leastsquares/CircleProblem.java @@ -0,0 +1,175 @@ +/* + * 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.math3.fitting.leastsquares; + +import java.util.ArrayList; +import org.apache.commons.math3.analysis.MultivariateVectorFunction; +import org.apache.commons.math3.analysis.MultivariateMatrixFunction; +import org.apache.commons.math3.util.MathUtils; +import org.apache.commons.math3.util.FastMath; + +/** + * Class that models a circle. + * The parameters of problem are: + *
      + *
    • the x-coordinate of the circle center,
      • + *
    • the y-coordinate of the circle center,
      • + *
    • the radius of the circle.
      • + *
    + * The model functions are: + *
      + *
    • for each triplet (cx, cy, r), the (x, y) coordinates of a point on the + * corresponding circle.
      • + *
    + */ +class CircleProblem { + /** Cloud of points assumed to be fitted by a circle. */ + private final ArrayList points; + /** Error on the x-coordinate of the points. */ + private final double xSigma; + /** Error on the y-coordinate of the points. */ + private final double ySigma; + /** Number of points on the circumference (when searching which + model point is closest to a given "observation". */ + private final int resolution; + + /** + * @param xError Assumed error for the x-coordinate of the circle points. + * @param yError Assumed error for the y-coordinate of the circle points. + * @param searchResolution Number of points to try when searching the one + * that is closest to a given "observed" point. + */ + public CircleProblem(double xError, + double yError, + int searchResolution) { + points = new ArrayList(); + xSigma = xError; + ySigma = yError; + resolution = searchResolution; + } + + /** + * @param xError Assumed error for the x-coordinate of the circle points. + * @param yError Assumed error for the y-coordinate of the circle points. + */ + public CircleProblem(double xError, + double yError) { + this(xError, yError, 500); + } + + public void addPoint(double px, double py) { + points.add(new double[] { px, py }); + } + + public double[] target() { + final double[] t = new double[points.size() * 2]; + for (int i = 0; i < points.size(); i++) { + final double[] p = points.get(i); + final int index = i * 2; + t[index] = p[0]; + t[index + 1] = p[1]; + } + + return t; + } + + public double[] weight() { + final double wX = 1 / (xSigma * xSigma); + final double wY = 1 / (ySigma * ySigma); + final double[] w = new double[points.size() * 2]; + for (int i = 0; i < points.size(); i++) { + final int index = i * 2; + w[index] = wX; + w[index + 1] = wY; + } + + return w; + } + + public MultivariateVectorFunction getModelFunction() { + return new MultivariateVectorFunction() { + public double[] value(double[] params) { + final double cx = params[0]; + final double cy = params[1]; + final double r = params[2]; + + final double[] model = new double[points.size() * 2]; + + final double deltaTheta = MathUtils.TWO_PI / resolution; + for (int i = 0; i < points.size(); i++) { + final double[] p = points.get(i); + final double px = p[0]; + final double py = p[1]; + + double bestX = 0; + double bestY = 0; + double dMin = Double.POSITIVE_INFINITY; + + // Find the angle for which the circle passes closest to the + // current point (using a resolution of 100 points along the + // circumference). + for (double theta = 0; theta <= MathUtils.TWO_PI; theta += deltaTheta) { + final double currentX = cx + r * FastMath.cos(theta); + final double currentY = cy + r * FastMath.sin(theta); + final double dX = currentX - px; + final double dY = currentY - py; + final double d = dX * dX + dY * dY; + if (d < dMin) { + dMin = d; + bestX = currentX; + bestY = currentY; + } + } + + final int index = i * 2; + model[index] = bestX; + model[index + 1] = bestY; + } + + return model; + } + }; + } + + public MultivariateMatrixFunction getModelFunctionJacobian() { + return new MultivariateMatrixFunction() { + public double[][] value(double[] point) { + return jacobian(point); + } + }; + } + + private double[][] jacobian(double[] params) { + final double[][] jacobian = new double[points.size() * 2][3]; + + for (int i = 0; i < points.size(); i++) { + final int index = i * 2; + // Partial derivative wrt x-coordinate of center. + jacobian[index][0] = 1; + jacobian[index + 1][0] = 0; + // Partial derivative wrt y-coordinate of center. + jacobian[index][1] = 0; + jacobian[index + 1][1] = 1; + // Partial derivative wrt radius. + final double[] p = points.get(i); + jacobian[index][2] = (p[0] - params[0]) / params[2]; + jacobian[index + 1][2] = (p[1] - params[1]) / params[2]; + } + + return jacobian; + } +} diff --git a/src/test/java/org/apache/commons/math3/fitting/leastsquares/CircleVectorial.java b/src/test/java/org/apache/commons/math3/fitting/leastsquares/CircleVectorial.java new file mode 100644 index 000000000..968bea5e3 --- /dev/null +++ b/src/test/java/org/apache/commons/math3/fitting/leastsquares/CircleVectorial.java @@ -0,0 +1,94 @@ +/* + * 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.math3.fitting.leastsquares; + +import java.util.ArrayList; +import org.apache.commons.math3.analysis.MultivariateVectorFunction; +import org.apache.commons.math3.analysis.MultivariateMatrixFunction; +import org.apache.commons.math3.geometry.euclidean.twod.Vector2D; + +/** + * Class used in the tests. + */ +class CircleVectorial { + private ArrayList points; + + public CircleVectorial() { + points = new ArrayList(); + } + + public void addPoint(double px, double py) { + points.add(new Vector2D(px, py)); + } + + public int getN() { + return points.size(); + } + + public double getRadius(Vector2D center) { + double r = 0; + for (Vector2D point : points) { + r += point.distance(center); + } + return r / points.size(); + } + + public MultivariateVectorFunction getModelFunction() { + return new MultivariateVectorFunction() { + public double[] value(double[] params) { + Vector2D center = new Vector2D(params[0], params[1]); + double radius = getRadius(center); + double[] residuals = new double[points.size()]; + for (int i = 0; i < residuals.length; i++) { + residuals[i] = points.get(i).distance(center) - radius; + } + + return residuals; + } + }; + } + + public MultivariateMatrixFunction getModelFunctionJacobian() { + return new MultivariateMatrixFunction() { + public double[][] value(double[] params) { + final int n = points.size(); + final Vector2D center = new Vector2D(params[0], params[1]); + + double dRdX = 0; + double dRdY = 0; + for (Vector2D pk : points) { + double dk = pk.distance(center); + dRdX += (center.getX() - pk.getX()) / dk; + dRdY += (center.getY() - pk.getY()) / dk; + } + dRdX /= n; + dRdY /= n; + + // Jacobian of the radius residuals. + double[][] jacobian = new double[n][2]; + for (int i = 0; i < n; i++) { + final Vector2D pi = points.get(i); + final double di = pi.distance(center); + jacobian[i][0] = (center.getX() - pi.getX()) / di - dRdX; + jacobian[i][1] = (center.getY() - pi.getY()) / di - dRdY; + } + + return jacobian; + } + }; + } +} diff --git a/src/test/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizerTest.java b/src/test/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizerTest.java new file mode 100644 index 000000000..8b707e023 --- /dev/null +++ b/src/test/java/org/apache/commons/math3/fitting/leastsquares/GaussNewtonOptimizerTest.java @@ -0,0 +1,109 @@ +/* + * 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.math3.fitting.leastsquares; + +import java.io.IOException; +import org.apache.commons.math3.exception.ConvergenceException; +import org.apache.commons.math3.exception.TooManyEvaluationsException; +import org.apache.commons.math3.exception.MathUnsupportedOperationException; +import org.apache.commons.math3.optim.SimpleVectorValueChecker; +import org.apache.commons.math3.linear.DiagonalMatrix; +import org.junit.Test; + +/** + *

    Some of the unit tests are re-implementations of the MINPACK file17 and file22 test files. + * The redistribution policy for MINPACK is available here/ + * + * @version $Id$ + */ +public class GaussNewtonOptimizerTest + extends AbstractLeastSquaresOptimizerAbstractTest { + @Override + public GaussNewtonOptimizer createOptimizer() { + return GaussNewtonOptimizer.create() + .withConvergenceChecker(new SimpleVectorValueChecker(1e-6, 1e-6)); + } + + @Override + public int getMaxIterations() { + return 1000; + } + + @Override + @Test(expected=ConvergenceException.class) + public void testMoreEstimatedParametersSimple() { + /* + * Exception is expected with this optimizer + */ + super.testMoreEstimatedParametersSimple(); + } + + @Override + @Test(expected=ConvergenceException.class) + public void testMoreEstimatedParametersUnsorted() { + /* + * Exception is expected with this optimizer + */ + super.testMoreEstimatedParametersUnsorted(); + } + + @Test(expected=TooManyEvaluationsException.class) + public void testMaxEvaluations() throws Exception { + CircleVectorial circle = new CircleVectorial(); + 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); + + GaussNewtonOptimizer optimizer = createOptimizer() + .withConvergenceChecker(new SimpleVectorValueChecker(1e-30, 1e-30)) + .withMaxIterations(Integer.MAX_VALUE) + .withMaxEvaluations(100) + .withModelAndJacobian(circle.getModelFunction(), + circle.getModelFunctionJacobian()) + .withTarget(new double[] { 0, 0, 0, 0, 0 }) + .withWeight(new DiagonalMatrix(new double[] { 1, 1, 1, 1, 1 })) + .withStartPoint(new double[] { 98.680, 47.345 }); + + optimizer.optimize(); + } + + @Override + @Test(expected=ConvergenceException.class) + public void testCircleFittingBadInit() { + /* + * This test does not converge with this optimizer. + */ + super.testCircleFittingBadInit(); + } + + @Override + @Test(expected=ConvergenceException.class) + public void testHahn1() + throws IOException { + /* + * TODO This test leads to a singular problem with the Gauss-Newton + * optimizer. This should be inquired. + */ + super.testHahn1(); + } +} diff --git a/src/test/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizerTest.java b/src/test/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizerTest.java new file mode 100644 index 000000000..2df665bee --- /dev/null +++ b/src/test/java/org/apache/commons/math3/fitting/leastsquares/LevenbergMarquardtOptimizerTest.java @@ -0,0 +1,356 @@ +/* + * 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.math3.fitting.leastsquares; + +import java.io.Serializable; +import java.util.ArrayList; +import java.util.List; +import org.apache.commons.math3.optim.PointVectorValuePair; +import org.apache.commons.math3.analysis.MultivariateVectorFunction; +import org.apache.commons.math3.analysis.MultivariateMatrixFunction; +import org.apache.commons.math3.exception.ConvergenceException; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.TooManyEvaluationsException; +import org.apache.commons.math3.exception.MathUnsupportedOperationException; +import org.apache.commons.math3.geometry.euclidean.twod.Vector2D; +import org.apache.commons.math3.linear.SingularMatrixException; +import org.apache.commons.math3.linear.DiagonalMatrix; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.util.Precision; +import org.junit.Assert; +import org.junit.Test; +import org.junit.Ignore; + +/** + *

    Some of the unit tests are re-implementations of the MINPACK file17 and file22 test files. + * The redistribution policy for MINPACK is available here. + * + * @version $Id$ + */ +public class LevenbergMarquardtOptimizerTest + extends AbstractLeastSquaresOptimizerAbstractTest { + @Override + public LevenbergMarquardtOptimizer createOptimizer() { + return LevenbergMarquardtOptimizer.create(); + } + + @Override + public int getMaxIterations() { + return 25; + } + + @Override + @Test(expected=SingularMatrixException.class) + public void testNonInvertible() { + /* + * Overrides the method from parent class, since the default singularity + * threshold (1e-14) does not trigger the expected exception. + */ + LinearProblem problem = new LinearProblem(new double[][] { + { 1, 2, -3 }, + { 2, 1, 3 }, + { -3, 0, -9 } + }, new double[] { 1, 1, 1 }); + + final LevenbergMarquardtOptimizer optimizer = createOptimizer() + .withMaxEvaluations(100) + .withMaxIterations(20) + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(problem.getTarget()) + .withWeight(new DiagonalMatrix(new double[] { 1, 1, 1 })) + .withStartPoint(new double[] { 0, 0, 0 }); + + final double[] optimum = optimizer.optimize().getPoint(); + Assert.assertTrue(FastMath.sqrt(optimizer.getTarget().length) * optimizer.computeRMS(optimum) > 0.6); + + optimizer.computeCovariances(optimum, 1.5e-14); + } + + @Test + public void testControlParameters() { + CircleVectorial circle = new CircleVectorial(); + 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.getModelFunction(), + circle.getModelFunctionJacobian(), + 0.1, 10, 1.0e-14, 1.0e-16, 1.0e-10, false); + checkEstimate(circle.getModelFunction(), + circle.getModelFunctionJacobian(), + 0.1, 10, 1.0e-15, 1.0e-17, 1.0e-10, true); + checkEstimate(circle.getModelFunction(), + circle.getModelFunctionJacobian(), + 0.1, 5, 1.0e-15, 1.0e-16, 1.0e-10, true); + circle.addPoint(300, -300); + checkEstimate(circle.getModelFunction(), + circle.getModelFunctionJacobian(), + 0.1, 20, 1.0e-18, 1.0e-16, 1.0e-10, true); + } + + private void checkEstimate(MultivariateVectorFunction problem, + MultivariateMatrixFunction problemJacobian, + double initialStepBoundFactor, int maxCostEval, + double costRelativeTolerance, double parRelativeTolerance, + double orthoTolerance, boolean shouldFail) { + try { + final LevenbergMarquardtOptimizer optimizer = LevenbergMarquardtOptimizer.create() + .withTuningParameters(initialStepBoundFactor, + costRelativeTolerance, + parRelativeTolerance, + orthoTolerance, + Precision.SAFE_MIN) + .withMaxEvaluations(maxCostEval) + .withMaxIterations(100) + .withModelAndJacobian(problem, problemJacobian) + .withTarget(new double[] { 0, 0, 0, 0, 0 }) + .withWeight(new DiagonalMatrix(new double[] { 1, 1, 1, 1, 1 })) + .withStartPoint(new double[] { 98.680, 47.345 }); + + optimizer.optimize(); + + Assert.assertTrue(!shouldFail); + } catch (DimensionMismatchException ee) { + Assert.assertTrue(shouldFail); + } catch (TooManyEvaluationsException ee) { + Assert.assertTrue(shouldFail); + } + } + + /** + * Non-linear test case: fitting of decay curve (from Chapter 8 of + * Bevington's textbook, "Data reduction and analysis for the physical sciences"). + * XXX The expected ("reference") values may not be accurate and the tolerance too + * relaxed for this test to be currently really useful (the issue is under + * investigation). + */ + @Test + public void testBevington() { + final double[][] dataPoints = { + // column 1 = times + { 15, 30, 45, 60, 75, 90, 105, 120, 135, 150, + 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, + 315, 330, 345, 360, 375, 390, 405, 420, 435, 450, + 465, 480, 495, 510, 525, 540, 555, 570, 585, 600, + 615, 630, 645, 660, 675, 690, 705, 720, 735, 750, + 765, 780, 795, 810, 825, 840, 855, 870, 885, }, + // column 2 = measured counts + { 775, 479, 380, 302, 185, 157, 137, 119, 110, 89, + 74, 61, 66, 68, 48, 54, 51, 46, 55, 29, + 28, 37, 49, 26, 35, 29, 31, 24, 25, 35, + 24, 30, 26, 28, 21, 18, 20, 27, 17, 17, + 14, 17, 24, 11, 22, 17, 12, 10, 13, 16, + 9, 9, 14, 21, 17, 13, 12, 18, 10, }, + }; + + final BevingtonProblem problem = new BevingtonProblem(); + + final int len = dataPoints[0].length; + final double[] weights = new double[len]; + for (int i = 0; i < len; i++) { + problem.addPoint(dataPoints[0][i], + dataPoints[1][i]); + + weights[i] = 1 / dataPoints[1][i]; + } + + final LevenbergMarquardtOptimizer optimizer = LevenbergMarquardtOptimizer.create() + .withMaxEvaluations(100) + .withMaxIterations(20) + .withModelAndJacobian(problem.getModelFunction(), + problem.getModelFunctionJacobian()) + .withTarget(dataPoints[1]) + .withWeight(new DiagonalMatrix(weights)) + .withStartPoint(new double[] { 10, 900, 80, 27, 225 }); + + final PointVectorValuePair optimum = optimizer.optimize(); + final double[] solution = optimum.getPoint(); + final double[] expectedSolution = { 10.4, 958.3, 131.4, 33.9, 205.0 }; + + final double[][] covarMatrix = optimizer.computeCovariances(solution, 1e-14); + final double[][] expectedCovarMatrix = { + { 3.38, -3.69, 27.98, -2.34, -49.24 }, + { -3.69, 2492.26, 81.89, -69.21, -8.9 }, + { 27.98, 81.89, 468.99, -44.22, -615.44 }, + { -2.34, -69.21, -44.22, 6.39, 53.80 }, + { -49.24, -8.9, -615.44, 53.8, 929.45 } + }; + + final int numParams = expectedSolution.length; + + // Check that the computed solution is within the reference error range. + for (int i = 0; i < numParams; i++) { + final double error = FastMath.sqrt(expectedCovarMatrix[i][i]); + Assert.assertEquals("Parameter " + i, expectedSolution[i], solution[i], error); + } + + // Check that each entry of the computed covariance matrix is within 10% + // of the reference matrix entry. + for (int i = 0; i < numParams; i++) { + for (int j = 0; j < numParams; j++) { + Assert.assertEquals("Covariance matrix [" + i + "][" + j + "]", + expectedCovarMatrix[i][j], + covarMatrix[i][j], + FastMath.abs(0.1 * expectedCovarMatrix[i][j])); + } + } + } + + @Test + public void testCircleFitting2() { + final double xCenter = 123.456; + final double yCenter = 654.321; + final double xSigma = 10; + final double ySigma = 15; + final double radius = 111.111; + // The test is extremely sensitive to the seed. + final long seed = 59421061L; + final RandomCirclePointGenerator factory + = new RandomCirclePointGenerator(xCenter, yCenter, radius, + xSigma, ySigma, + seed); + final CircleProblem circle = new CircleProblem(xSigma, ySigma); + + final int numPoints = 10; + for (Vector2D p : factory.generate(numPoints)) { + circle.addPoint(p.getX(), p.getY()); + } + + // First guess for the center's coordinates and radius. + final double[] init = { 90, 659, 115 }; + + final LevenbergMarquardtOptimizer optimizer = LevenbergMarquardtOptimizer.create() + .withMaxEvaluations(100) + .withMaxIterations(50) + .withModelAndJacobian(circle.getModelFunction(), + circle.getModelFunctionJacobian()) + .withTarget(circle.target()) + .withWeight(new DiagonalMatrix(circle.weight())) + .withStartPoint(init); + + final PointVectorValuePair optimum = optimizer.optimize(); + final double[] paramFound = optimum.getPoint(); + + // Retrieve errors estimation. + final double[] asymptoticStandardErrorFound = optimizer.computeSigma(paramFound, 1e-14); + + // Check that the parameters are found within the assumed error bars. + Assert.assertEquals(xCenter, paramFound[0], asymptoticStandardErrorFound[0]); + Assert.assertEquals(yCenter, paramFound[1], asymptoticStandardErrorFound[1]); + Assert.assertEquals(radius, paramFound[2], asymptoticStandardErrorFound[2]); + } + + private static class QuadraticProblem { + private List x; + private List y; + + public QuadraticProblem() { + x = new ArrayList(); + y = new ArrayList(); + } + + public void addPoint(double x, double y) { + this.x.add(x); + this.y.add(y); + } + + public MultivariateVectorFunction getModelFunction() { + return new MultivariateVectorFunction() { + public double[] value(double[] variables) { + double[] values = new double[x.size()]; + for (int i = 0; i < values.length; ++i) { + values[i] = (variables[0] * x.get(i) + variables[1]) * x.get(i) + variables[2]; + } + return values; + } + }; + } + + public MultivariateMatrixFunction getModelFunctionJacobian() { + return new MultivariateMatrixFunction() { + public double[][] value(double[] params) { + double[][] jacobian = new double[x.size()][3]; + for (int i = 0; i < jacobian.length; ++i) { + jacobian[i][0] = x.get(i) * x.get(i); + jacobian[i][1] = x.get(i); + jacobian[i][2] = 1.0; + } + return jacobian; + } + }; + } + } + + private static class BevingtonProblem { + private List time; + private List count; + + public BevingtonProblem() { + time = new ArrayList(); + count = new ArrayList(); + } + + public void addPoint(double t, double c) { + time.add(t); + count.add(c); + } + + public MultivariateVectorFunction getModelFunction() { + return new MultivariateVectorFunction() { + public double[] value(double[] params) { + double[] values = new double[time.size()]; + for (int i = 0; i < values.length; ++i) { + final double t = time.get(i); + values[i] = params[0] + + params[1] * Math.exp(-t / params[3]) + + params[2] * Math.exp(-t / params[4]); + } + return values; + } + }; + } + + public MultivariateMatrixFunction getModelFunctionJacobian() { + return new MultivariateMatrixFunction() { + public double[][] value(double[] params) { + double[][] jacobian = new double[time.size()][5]; + + for (int i = 0; i < jacobian.length; ++i) { + final double t = time.get(i); + jacobian[i][0] = 1; + + final double p3 = params[3]; + final double p4 = params[4]; + final double tOp3 = t / p3; + final double tOp4 = t / p4; + jacobian[i][1] = Math.exp(-tOp3); + jacobian[i][2] = Math.exp(-tOp4); + jacobian[i][3] = params[1] * Math.exp(-tOp3) * tOp3 / p3; + jacobian[i][4] = params[2] * Math.exp(-tOp4) * tOp4 / p4; + } + return jacobian; + } + }; + } + } +} diff --git a/src/test/java/org/apache/commons/math3/fitting/leastsquares/MinpackTest.java b/src/test/java/org/apache/commons/math3/fitting/leastsquares/MinpackTest.java new file mode 100644 index 000000000..adc77ebce --- /dev/null +++ b/src/test/java/org/apache/commons/math3/fitting/leastsquares/MinpackTest.java @@ -0,0 +1,1499 @@ +/* + * 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.math3.fitting.leastsquares; + +import java.io.Serializable; +import java.util.Arrays; +import org.apache.commons.math3.exception.TooManyEvaluationsException; +import org.apache.commons.math3.analysis.MultivariateVectorFunction; +import org.apache.commons.math3.analysis.MultivariateMatrixFunction; +import org.apache.commons.math3.optim.PointVectorValuePair; +import org.apache.commons.math3.linear.DiagonalMatrix; +import org.apache.commons.math3.util.FastMath; +import org.junit.Assert; +import org.junit.Test; + +/** + *

    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) + */ +public class MinpackTest { + + @Test + public void testMinpackLinearFullRank() { + minpackTest(new LinearFullRankFunction(10, 5, 1.0, + 5.0, 2.23606797749979), false); + minpackTest(new LinearFullRankFunction(50, 5, 1.0, + 8.06225774829855, 6.70820393249937), false); + } + + @Test + public void testMinpackLinearRank1() { + minpackTest(new LinearRank1Function(10, 5, 1.0, + 291.521868819476, 1.4638501094228), false); + minpackTest(new LinearRank1Function(50, 5, 1.0, + 3101.60039334535, 3.48263016573496), false); + } + + @Test + public void testMinpackLinearRank1ZeroColsAndRows() { + minpackTest(new LinearRank1ZeroColsAndRowsFunction(10, 5, 1.0), false); + minpackTest(new LinearRank1ZeroColsAndRowsFunction(50, 5, 1.0), false); + } + + @Test + public void testMinpackRosenbrok() { + minpackTest(new RosenbrockFunction(new double[] { -1.2, 1.0 }, + FastMath.sqrt(24.2)), false); + minpackTest(new RosenbrockFunction(new double[] { -12.0, 10.0 }, + FastMath.sqrt(1795769.0)), false); + minpackTest(new RosenbrockFunction(new double[] { -120.0, 100.0 }, + 11.0 * FastMath.sqrt(169000121.0)), false); + } + + @Test + public void testMinpackHelicalValley() { + 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); + } + + @Test + public void testMinpackPowellSingular() { + 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); + } + + @Test + public void testMinpackFreudensteinRoth() { + 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.41300466147456, + -0.896796038685959 + }), false); + minpackTest(new FreudensteinRothFunction(new double[] { 50.0, -200.0 }, + 11426454.595762, 6.99887517242903, + new double[] { + 11.412781785788564, + -0.8968051074920405 + }), false); + } + + @Test + public void testMinpackBard() { + 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); + } + + @Test + public void testMinpackKowalikOsborne() { + 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 + }), false); + } + + @Test + public void testMinpackMeyer() { + 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); + } + + @Test + public void testMinpackWatson() { + 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); + } + + @Test + public void testMinpackBox3Dimensional() { + minpackTest(new Box3DimensionalFunction(10, new double[] { 0.0, 10.0, 20.0 }, + 32.1115837449572), false); + } + + @Test + public void testMinpackJennrichSampson() { + minpackTest(new JennrichSampsonFunction(10, new double[] { 0.3, 0.4 }, + 64.5856498144943, 11.1517793413499, + new double[] { +// 0.2578330049, 0.257829976764542 + 0.2578199266368004, 0.25782997676455244 + }), false); + } + + @Test + public void testMinpackBrownDennis() { + 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); + } + + @Test + public void testMinpackChebyquad() { + 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); + } + + @Test + public void testMinpackBrownAlmostLinear() { + 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); + } + + @Test + public void testMinpackOsborne1() { + 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); + } + + @Test + public void testMinpackOsborne2() { + 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) { + final double tol = 2.22044604926e-16; + LevenbergMarquardtOptimizer optimizer = LevenbergMarquardtOptimizer.create(); + optimizer = optimizer + .withTuningParameters(optimizer.getInitialStepBoundFactor(), + FastMath.sqrt(tol), + FastMath.sqrt(tol), + tol, + optimizer.getRankingThreshold()) + .withMaxEvaluations(400 * (function.getN() + 1)) + .withMaxIterations(2000) + .withModelAndJacobian(function.getModelFunction(), + function.getModelFunctionJacobian()) + .withTarget(function.getTarget()) + .withWeight(new DiagonalMatrix(function.getWeight())) + .withStartPoint(function.getStartPoint()); + + try { + final double[] optimum = optimizer.optimize().getPoint(); + Assert.assertFalse(exceptionExpected); + function.checkTheoreticalMinCost(optimizer.computeRMS(optimum)); + function.checkTheoreticalMinParams(optimum); + } catch (TooManyEvaluationsException e) { + Assert.assertTrue(exceptionExpected); + } + } + + private static abstract class MinpackFunction { + protected int n; + protected int m; + protected double[] startParams; + protected double theoreticalMinCost; + protected double[] theoreticalMinParams; + protected double costAccuracy; + protected double paramsAccuracy; + + protected MinpackFunction(int m, double[] startParams, + double theoreticalMinCost, + double[] theoreticalMinParams) { + this.m = m; + this.n = startParams.length; + this.startParams = startParams.clone(); + 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; + } + + public double[] getTarget() { + return buildArray(m, 0.0); + } + + public double[] getWeight() { + return buildArray(m, 1.0); + } + + public double[] getStartPoint() { + return startParams.clone(); + } + + protected void setCostAccuracy(double costAccuracy) { + this.costAccuracy = costAccuracy; + } + + protected void setParamsAccuracy(double paramsAccuracy) { + this.paramsAccuracy = paramsAccuracy; + } + + public int getN() { + return startParams.length; + } + + public void checkTheoreticalMinCost(double rms) { + double threshold = costAccuracy * (1.0 + theoreticalMinCost); + Assert.assertEquals(theoreticalMinCost, FastMath.sqrt(m) * rms, threshold); + } + + public void checkTheoreticalMinParams(double[] params) { + if (theoreticalMinParams != null) { + for (int i = 0; i < theoreticalMinParams.length; ++i) { + double mi = theoreticalMinParams[i]; + double vi = params[i]; + Assert.assertEquals(mi, vi, paramsAccuracy * (1.0 + FastMath.abs(mi))); + } + } + } + + public MultivariateVectorFunction getModelFunction() { + return new MultivariateVectorFunction() { + public double[] value(double[] point) { + return computeValue(point); + } + }; + } + + public MultivariateMatrixFunction getModelFunctionJacobian() { + return new MultivariateMatrixFunction() { + public double[][] value(double[] point) { + return computeJacobian(point); + } + }; + } + + public abstract double[][] computeJacobian(double[] variables); + public abstract double[] computeValue(double[] variables); + } + + private static class LinearFullRankFunction extends MinpackFunction { + private static final long serialVersionUID = -9030323226268039536L; + + public LinearFullRankFunction(int m, int n, double x0, + double theoreticalStartCost, + double theoreticalMinCost) { + super(m, buildArray(n, x0), theoreticalMinCost, + buildArray(n, -1.0)); + } + + @Override + public double[][] computeJacobian(double[] variables) { + 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; + } + + @Override + public double[] computeValue(double[] variables) { + double sum = 0; + for (int i = 0; i < n; ++i) { + sum += variables[i]; + } + double t = 1 + 2 * sum / m; + double[] f = new double[m]; + for (int i = 0; i < n; ++i) { + f[i] = variables[i] - t; + } + Arrays.fill(f, n, m, -t); + return f; + } + } + + private static class LinearRank1Function extends MinpackFunction { + private static final long serialVersionUID = 8494863245104608300L; + + public LinearRank1Function(int m, int n, double x0, + double theoreticalStartCost, + double theoreticalMinCost) { + super(m, buildArray(n, x0), theoreticalMinCost, null); + } + + @Override + public double[][] computeJacobian(double[] variables) { + 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; + } + + @Override + public double[] computeValue(double[] variables) { + double[] f = new double[m]; + double sum = 0; + for (int i = 0; i < n; ++i) { + sum += (i + 1) * variables[i]; + } + for (int i = 0; i < m; ++i) { + f[i] = (i + 1) * sum - 1; + } + return f; + } + } + + private static class LinearRank1ZeroColsAndRowsFunction extends MinpackFunction { + private static final long serialVersionUID = -3316653043091995018L; + + public LinearRank1ZeroColsAndRowsFunction(int m, int n, double x0) { + super(m, buildArray(n, x0), + FastMath.sqrt((m * (m + 3) - 6) / (2.0 * (2 * m - 3))), + null); + } + + @Override + public double[][] computeJacobian(double[] variables) { + 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; + } + + @Override + public double[] computeValue(double[] variables) { + double[] f = new double[m]; + double sum = 0; + for (int i = 1; i < (n - 1); ++i) { + sum += (i + 1) * variables[i]; + } + 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 { + private static final long serialVersionUID = 2893438180956569134L; + public RosenbrockFunction(double[] startParams, double theoreticalStartCost) { + super(2, startParams, 0.0, buildArray(2, 1.0)); + } + + @Override + public double[][] computeJacobian(double[] variables) { + double x1 = variables[0]; + return new double[][] { { -20 * x1, 10 }, { -1, 0 } }; + } + + @Override + public double[] computeValue(double[] variables) { + double x1 = variables[0]; + double x2 = variables[1]; + return new double[] { 10 * (x2 - x1 * x1), 1 - x1 }; + } + } + + private static class HelicalValleyFunction extends MinpackFunction { + private static final long serialVersionUID = 220613787843200102L; + public HelicalValleyFunction(double[] startParams, + double theoreticalStartCost) { + super(3, startParams, 0.0, new double[] { 1.0, 0.0, 0.0 }); + } + + @Override + public double[][] computeJacobian(double[] variables) { + double x1 = variables[0]; + double x2 = variables[1]; + double tmpSquare = x1 * x1 + x2 * x2; + double tmp1 = twoPi * tmpSquare; + double tmp2 = FastMath.sqrt(tmpSquare); + return new double[][] { + { 100 * x2 / tmp1, -100 * x1 / tmp1, 10 }, + { 10 * x1 / tmp2, 10 * x2 / tmp2, 0 }, + { 0, 0, 1 } + }; + } + + @Override + public double[] computeValue(double[] variables) { + double x1 = variables[0]; + double x2 = variables[1]; + double x3 = variables[2]; + double tmp1; + if (x1 == 0) { + tmp1 = (x2 >= 0) ? 0.25 : -0.25; + } else { + tmp1 = FastMath.atan(x2 / x1) / twoPi; + if (x1 < 0) { + tmp1 += 0.5; + } + } + double tmp2 = FastMath.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 * FastMath.PI; + } + + private static class PowellSingularFunction extends MinpackFunction { + private static final long serialVersionUID = 7298364171208142405L; + + public PowellSingularFunction(double[] startParams, + double theoreticalStartCost) { + super(4, startParams, 0.0, buildArray(4, 0.0)); + } + + @Override + public double[][] computeJacobian(double[] variables) { + double x1 = variables[0]; + double x2 = variables[1]; + double x3 = variables[2]; + double x4 = variables[3]; + 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) } + }; + } + + @Override + public double[] computeValue(double[] variables) { + double x1 = variables[0]; + double x2 = variables[1]; + double x3 = variables[2]; + double x4 = variables[3]; + 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 = FastMath.sqrt( 5.0); + private static final double sqrt10 = FastMath.sqrt(10.0); + } + + private static class FreudensteinRothFunction extends MinpackFunction { + private static final long serialVersionUID = 2892404999344244214L; + + public FreudensteinRothFunction(double[] startParams, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(2, startParams, theoreticalMinCost, + theoreticalMinParams); + } + + @Override + public double[][] computeJacobian(double[] variables) { + double x2 = variables[1]; + return new double[][] { + { 1, x2 * (10 - 3 * x2) - 2 }, + { 1, x2 * ( 2 + 3 * x2) - 14, } + }; + } + + @Override + public double[] computeValue(double[] variables) { + double x1 = variables[0]; + double x2 = variables[1]; + 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 { + private static final long serialVersionUID = 5990442612572087668L; + + public BardFunction(double x0, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(15, buildArray(3, x0), theoreticalMinCost, + theoreticalMinParams); + } + + @Override + public double[][] computeJacobian(double[] variables) { + double x2 = variables[1]; + double x3 = variables[2]; + 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; + } + + @Override + public double[] computeValue(double[] variables) { + double x1 = variables[0]; + double x2 = variables[1]; + double x3 = variables[2]; + 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 { + private static final long serialVersionUID = -4867445739880495801L; + + public KowalikOsborneFunction(double[] startParams, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(11, startParams, theoreticalMinCost, + theoreticalMinParams); + if (theoreticalStartCost > 20.0) { + setCostAccuracy(2.0e-4); + setParamsAccuracy(5.0e-3); + } + } + + @Override + public double[][] computeJacobian(double[] variables) { + double x1 = variables[0]; + double x2 = variables[1]; + double x3 = variables[2]; + double x4 = variables[3]; + 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; + } + + @Override + public double[] computeValue(double[] variables) { + double x1 = variables[0]; + double x2 = variables[1]; + double x3 = variables[2]; + double x4 = variables[3]; + 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 { + private static final long serialVersionUID = -838060619150131027L; + + public MeyerFunction(double[] startParams, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(16, startParams, theoreticalMinCost, + theoreticalMinParams); + if (theoreticalStartCost > 1.0e6) { + setCostAccuracy(7.0e-3); + setParamsAccuracy(2.0e-2); + } + } + + @Override + public double[][] computeJacobian(double[] variables) { + double x1 = variables[0]; + double x2 = variables[1]; + double x3 = variables[2]; + 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 = FastMath.exp(tmp1); + double tmp3 = x1 * tmp2 / temp; + jacobian[i] = new double[] { tmp2, tmp3, -tmp1 * tmp3 }; + } + return jacobian; + } + + @Override + public double[] computeValue(double[] variables) { + double x1 = variables[0]; + double x2 = variables[1]; + double x3 = variables[2]; + double[] f = new double[m]; + for (int i = 0; i < m; ++i) { + f[i] = x1 * FastMath.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 { + private static final long serialVersionUID = -9034759294980218927L; + + public WatsonFunction(int n, double x0, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(31, buildArray(n, x0), theoreticalMinCost, + theoreticalMinParams); + } + + @Override + public double[][] computeJacobian(double[] variables) { + 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 * variables[j]; + 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 * variables[0]; + jacobian[m - 1][1]= 1; + + return jacobian; + } + + @Override + public double[] computeValue(double[] variables) { + 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 * variables[j]; + dx *= div; + } + double s2 = 0; + dx = 1; + for (int j = 0; j < n; ++j) { + s2 += dx * variables[j]; + dx *= div; + } + f[i] = s1 - s2 * s2 - 1; + } + + double x1 = variables[0]; + double x2 = variables[1]; + f[m - 2] = x1; + f[m - 1] = x2 - x1 * x1 - 1; + + return f; + } + } + + private static class Box3DimensionalFunction extends MinpackFunction { + private static final long serialVersionUID = 5511403858142574493L; + + public Box3DimensionalFunction(int m, double[] startParams, + double theoreticalStartCost) { + super(m, startParams, 0.0, + new double[] { 1.0, 10.0, 1.0 }); + } + + @Override + public double[][] computeJacobian(double[] variables) { + double x1 = variables[0]; + double x2 = variables[1]; + double[][] jacobian = new double[m][]; + for (int i = 0; i < m; ++i) { + double tmp = (i + 1) / 10.0; + jacobian[i] = new double[] { + -tmp * FastMath.exp(-tmp * x1), + tmp * FastMath.exp(-tmp * x2), + FastMath.exp(-i - 1) - FastMath.exp(-tmp) + }; + } + return jacobian; + } + + @Override + public double[] computeValue(double[] variables) { + double x1 = variables[0]; + double x2 = variables[1]; + double x3 = variables[2]; + double[] f = new double[m]; + for (int i = 0; i < m; ++i) { + double tmp = (i + 1) / 10.0; + f[i] = FastMath.exp(-tmp * x1) - FastMath.exp(-tmp * x2) + + (FastMath.exp(-i - 1) - FastMath.exp(-tmp)) * x3; + } + return f; + } + } + + private static class JennrichSampsonFunction extends MinpackFunction { + private static final long serialVersionUID = -2489165190443352947L; + + public JennrichSampsonFunction(int m, double[] startParams, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(m, startParams, theoreticalMinCost, + theoreticalMinParams); + } + + @Override + public double[][] computeJacobian(double[] variables) { + double x1 = variables[0]; + double x2 = variables[1]; + double[][] jacobian = new double[m][]; + for (int i = 0; i < m; ++i) { + double t = i + 1; + jacobian[i] = new double[] { -t * FastMath.exp(t * x1), -t * FastMath.exp(t * x2) }; + } + return jacobian; + } + + @Override + public double[] computeValue(double[] variables) { + double x1 = variables[0]; + double x2 = variables[1]; + double[] f = new double[m]; + for (int i = 0; i < m; ++i) { + double temp = i + 1; + f[i] = 2 + 2 * temp - FastMath.exp(temp * x1) - FastMath.exp(temp * x2); + } + return f; + } + } + + private static class BrownDennisFunction extends MinpackFunction { + private static final long serialVersionUID = 8340018645694243910L; + + public BrownDennisFunction(int m, double[] startParams, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(m, startParams, theoreticalMinCost, + theoreticalMinParams); + setCostAccuracy(2.5e-8); + } + + @Override + public double[][] computeJacobian(double[] variables) { + double x1 = variables[0]; + double x2 = variables[1]; + double x3 = variables[2]; + double x4 = variables[3]; + double[][] jacobian = new double[m][]; + for (int i = 0; i < m; ++i) { + double temp = (i + 1) / 5.0; + double ti = FastMath.sin(temp); + double tmp1 = x1 + temp * x2 - FastMath.exp(temp); + double tmp2 = x3 + ti * x4 - FastMath.cos(temp); + jacobian[i] = new double[] { + 2 * tmp1, 2 * temp * tmp1, 2 * tmp2, 2 * ti * tmp2 + }; + } + return jacobian; + } + + @Override + public double[] computeValue(double[] variables) { + double x1 = variables[0]; + double x2 = variables[1]; + double x3 = variables[2]; + double x4 = variables[3]; + double[] f = new double[m]; + for (int i = 0; i < m; ++i) { + double temp = (i + 1) / 5.0; + double tmp1 = x1 + temp * x2 - FastMath.exp(temp); + double tmp2 = x3 + FastMath.sin(temp) * x4 - FastMath.cos(temp); + f[i] = tmp1 * tmp1 + tmp2 * tmp2; + } + return f; + } + } + + private static class ChebyquadFunction extends MinpackFunction { + private static final long serialVersionUID = -2394877275028008594L; + + 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), theoreticalMinCost, + theoreticalMinParams); + } + + @Override + public double[][] computeJacobian(double[] variables) { + 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 * variables[j] - 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; + } + + @Override + public double[] computeValue(double[] variables) { + double[] f = new double[m]; + + for (int j = 0; j < n; ++j) { + double tmp1 = 1; + double tmp2 = 2 * variables[j] - 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 { + private static final long serialVersionUID = 8239594490466964725L; + + public BrownAlmostLinearFunction(int m, double factor, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(m, buildArray(m, factor), theoreticalMinCost, + theoreticalMinParams); + } + + @Override + public double[][] computeJacobian(double[] variables) { + 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 *= variables[j]; + for (int i = 0; i < n; ++i) { + jacobian[i][j] = 1; + } + jacobian[j][j] = 2; + } + + for (int j = 0; j < n; ++j) { + double temp = variables[j]; + if (temp == 0) { + temp = 1; + prod = 1; + for (int k = 0; k < n; ++k) { + if (k != j) { + prod *= variables[k]; + } + } + } + jacobian[n - 1][j] = prod / temp; + } + + return jacobian; + } + + @Override + public double[] computeValue(double[] variables) { + double[] f = new double[m]; + double sum = -(n + 1); + double prod = 1; + for (int j = 0; j < n; ++j) { + sum += variables[j]; + prod *= variables[j]; + } + for (int i = 0; i < n; ++i) { + f[i] = variables[i] + sum; + } + f[n - 1] = prod - 1; + return f; + } + } + + private static class Osborne1Function extends MinpackFunction { + private static final long serialVersionUID = 4006743521149849494L; + + public Osborne1Function(double[] startParams, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(33, startParams, theoreticalMinCost, + theoreticalMinParams); + } + + @Override + public double[][] computeJacobian(double[] variables) { + double x2 = variables[1]; + double x3 = variables[2]; + double x4 = variables[3]; + double x5 = variables[4]; + double[][] jacobian = new double[m][]; + for (int i = 0; i < m; ++i) { + double temp = 10.0 * i; + double tmp1 = FastMath.exp(-temp * x4); + double tmp2 = FastMath.exp(-temp * x5); + jacobian[i] = new double[] { + -1, -tmp1, -tmp2, temp * x2 * tmp1, temp * x3 * tmp2 + }; + } + return jacobian; + } + + @Override + public double[] computeValue(double[] variables) { + double x1 = variables[0]; + double x2 = variables[1]; + double x3 = variables[2]; + double x4 = variables[3]; + double x5 = variables[4]; + double[] f = new double[m]; + for (int i = 0; i < m; ++i) { + double temp = 10.0 * i; + double tmp1 = FastMath.exp(-temp * x4); + double tmp2 = FastMath.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 { + private static final long serialVersionUID = -8418268780389858746L; + + public Osborne2Function(double[] startParams, + double theoreticalStartCost, + double theoreticalMinCost, + double[] theoreticalMinParams) { + super(65, startParams, theoreticalMinCost, + theoreticalMinParams); + } + + @Override + public double[][] computeJacobian(double[] variables) { + double x01 = variables[0]; + double x02 = variables[1]; + double x03 = variables[2]; + double x04 = variables[3]; + double x05 = variables[4]; + double x06 = variables[5]; + double x07 = variables[6]; + double x08 = variables[7]; + double x09 = variables[8]; + double x10 = variables[9]; + double x11 = variables[10]; + double[][] jacobian = new double[m][]; + for (int i = 0; i < m; ++i) { + double temp = i / 10.0; + double tmp1 = FastMath.exp(-x05 * temp); + double tmp2 = FastMath.exp(-x06 * (temp - x09) * (temp - x09)); + double tmp3 = FastMath.exp(-x07 * (temp - x10) * (temp - x10)); + double tmp4 = FastMath.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; + } + + @Override + public double[] computeValue(double[] variables) { + double x01 = variables[0]; + double x02 = variables[1]; + double x03 = variables[2]; + double x04 = variables[3]; + double x05 = variables[4]; + double x06 = variables[5]; + double x07 = variables[6]; + double x08 = variables[7]; + double x09 = variables[8]; + double x10 = variables[9]; + double x11 = variables[10]; + double[] f = new double[m]; + for (int i = 0; i < m; ++i) { + double temp = i / 10.0; + double tmp1 = FastMath.exp(-x05 * temp); + double tmp2 = FastMath.exp(-x06 * (temp - x09) * (temp - x09)); + double tmp3 = FastMath.exp(-x07 * (temp - x10) * (temp - x10)); + double tmp4 = FastMath.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 + }; + } +} diff --git a/src/test/java/org/apache/commons/math3/fitting/leastsquares/RandomCirclePointGenerator.java b/src/test/java/org/apache/commons/math3/fitting/leastsquares/RandomCirclePointGenerator.java new file mode 100644 index 000000000..3f9a20d11 --- /dev/null +++ b/src/test/java/org/apache/commons/math3/fitting/leastsquares/RandomCirclePointGenerator.java @@ -0,0 +1,91 @@ +/* + * 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.math3.fitting.leastsquares; + +import org.apache.commons.math3.random.RandomGenerator; +import org.apache.commons.math3.random.Well44497b; +import org.apache.commons.math3.util.MathUtils; +import org.apache.commons.math3.util.FastMath; +import org.apache.commons.math3.distribution.RealDistribution; +import org.apache.commons.math3.distribution.UniformRealDistribution; +import org.apache.commons.math3.distribution.NormalDistribution; +import org.apache.commons.math3.geometry.euclidean.twod.Vector2D; + +/** + * Factory for generating a cloud of points that approximate a circle. + */ +public class RandomCirclePointGenerator { + /** RNG for the x-coordinate of the center. */ + private final RealDistribution cX; + /** RNG for the y-coordinate of the center. */ + private final RealDistribution cY; + /** RNG for the parametric position of the point. */ + private final RealDistribution tP; + /** Radius of the circle. */ + private final double radius; + + /** + * @param x Abscissa of the circle center. + * @param y Ordinate of the circle center. + * @param radius Radius of the circle. + * @param xSigma Error on the x-coordinate of the circumference points. + * @param ySigma Error on the y-coordinate of the circumference points. + * @param seed RNG seed. + */ + public RandomCirclePointGenerator(double x, + double y, + double radius, + double xSigma, + double ySigma, + long seed) { + final RandomGenerator rng = new Well44497b(seed); + this.radius = radius; + cX = new NormalDistribution(rng, x, xSigma, + NormalDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY); + cY = new NormalDistribution(rng, y, ySigma, + NormalDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY); + tP = new UniformRealDistribution(rng, 0, MathUtils.TWO_PI, + UniformRealDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY); + } + + /** + * Point generator. + * + * @param n Number of points to create. + * @return the cloud of {@code n} points. + */ + public Vector2D[] generate(int n) { + final Vector2D[] cloud = new Vector2D[n]; + for (int i = 0; i < n; i++) { + cloud[i] = create(); + } + return cloud; + } + + /** + * Create one point. + * + * @return a point. + */ + private Vector2D create() { + final double t = tP.sample(); + final double pX = cX.sample() + radius * FastMath.cos(t); + final double pY = cY.sample() + radius * FastMath.sin(t); + + return new Vector2D(pX, pY); + } +} diff --git a/src/test/java/org/apache/commons/math3/fitting/leastsquares/RandomStraightLinePointGenerator.java b/src/test/java/org/apache/commons/math3/fitting/leastsquares/RandomStraightLinePointGenerator.java new file mode 100644 index 000000000..1c60cd21a --- /dev/null +++ b/src/test/java/org/apache/commons/math3/fitting/leastsquares/RandomStraightLinePointGenerator.java @@ -0,0 +1,98 @@ +/* + * 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.math3.fitting.leastsquares; + +import java.awt.geom.Point2D; +import org.apache.commons.math3.random.RandomGenerator; +import org.apache.commons.math3.random.Well44497b; +import org.apache.commons.math3.distribution.RealDistribution; +import org.apache.commons.math3.distribution.UniformRealDistribution; +import org.apache.commons.math3.distribution.NormalDistribution; + +/** + * Factory for generating a cloud of points that approximate a straight line. + */ +public class RandomStraightLinePointGenerator { + /** Slope. */ + private final double slope; + /** Intercept. */ + private final double intercept; + /** RNG for the x-coordinate. */ + private final RealDistribution x; + /** RNG for the error on the y-coordinate. */ + private final RealDistribution error; + + /** + * The generator will create a cloud of points whose x-coordinates + * will be randomly sampled between {@code xLo} and {@code xHi}, and + * the corresponding y-coordinates will be computed as + * 
      
+     *  y = a x + b + N(0, error)
+     *
      + * where {@code N(mean, sigma)} is a Gaussian distribution with the + * given mean and standard deviation. + * + * @param a Slope. + * @param b Intercept. + * @param sigma Standard deviation on the y-coordinate of the point. + * @param lo Lowest value of the x-coordinate. + * @param hi Highest value of the x-coordinate. + * @param seed RNG seed. + */ + public RandomStraightLinePointGenerator(double a, + double b, + double sigma, + double lo, + double hi, + long seed) { + final RandomGenerator rng = new Well44497b(seed); + slope = a; + intercept = b; + error = new NormalDistribution(rng, 0, sigma, + NormalDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY); + x = new UniformRealDistribution(rng, lo, hi, + UniformRealDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY); + } + + /** + * Point generator. + * + * @param n Number of points to create. + * @return the cloud of {@code n} points. + */ + public Point2D.Double[] generate(int n) { + final Point2D.Double[] cloud = new Point2D.Double[n]; + for (int i = 0; i < n; i++) { + cloud[i] = create(); + } + return cloud; + } + + /** + * Create one point. + * + * @return a point. + */ + private Point2D.Double create() { + final double abscissa = x.sample(); + final double yModel = slope * abscissa + intercept; + final double ordinate = yModel + error.sample(); + + return new Point2D.Double(abscissa, ordinate); + } +} diff --git a/src/test/java/org/apache/commons/math3/fitting/leastsquares/StatisticalReferenceDataset.java b/src/test/java/org/apache/commons/math3/fitting/leastsquares/StatisticalReferenceDataset.java new file mode 100644 index 000000000..5f3dbba07 --- /dev/null +++ b/src/test/java/org/apache/commons/math3/fitting/leastsquares/StatisticalReferenceDataset.java @@ -0,0 +1,370 @@ +/* + * 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.math3.fitting.leastsquares; + +import java.io.BufferedReader; +import java.io.IOException; +import java.util.ArrayList; +import org.apache.commons.math3.analysis.MultivariateVectorFunction; +import org.apache.commons.math3.analysis.MultivariateMatrixFunction; +import org.apache.commons.math3.util.MathArrays; + +/** + * This class gives access to the statistical reference datasets provided by the + * NIST (available + * here). + * Instances of this class can be created by invocation of the + * {@link StatisticalReferenceDatasetFactory}. + */ +public abstract class StatisticalReferenceDataset { + /** The name of this dataset. */ + private final String name; + /** The total number of observations (data points). */ + private final int numObservations; + /** The total number of parameters. */ + private final int numParameters; + /** The total number of starting points for the optimizations. */ + private final int numStartingPoints; + /** The values of the predictor. */ + private final double[] x; + /** The values of the response. */ + private final double[] y; + /** + * The starting values. {@code startingValues[j][i]} is the value of the + * {@code i}-th parameter in the {@code j}-th set of starting values. + */ + private final double[][] startingValues; + /** The certified values of the parameters. */ + private final double[] a; + /** The certified values of the standard deviation of the parameters. */ + private final double[] sigA; + /** The certified value of the residual sum of squares. */ + private double residualSumOfSquares; + /** The least-squares problem. */ + private final LeastSquaresProblem problem; + + /** + * Creates a new instance of this class from the specified data file. The + * file must follow the StRD format. + * + * @param in the data file + * @throws IOException if an I/O error occurs + */ + public StatisticalReferenceDataset(final BufferedReader in) + throws IOException { + + final ArrayList lines = new ArrayList(); + for (String line = in.readLine(); line != null; line = in.readLine()) { + lines.add(line); + } + int[] index = findLineNumbers("Data", lines); + if (index == null) { + throw new AssertionError("could not find line indices for data"); + } + this.numObservations = index[1] - index[0] + 1; + this.x = new double[this.numObservations]; + this.y = new double[this.numObservations]; + for (int i = 0; i < this.numObservations; i++) { + final String line = lines.get(index[0] + i - 1); + final String[] tokens = line.trim().split(" ++"); + // Data columns are in reverse order!!! + this.y[i] = Double.parseDouble(tokens[0]); + this.x[i] = Double.parseDouble(tokens[1]); + } + + index = findLineNumbers("Starting Values", lines); + if (index == null) { + throw new AssertionError( + "could not find line indices for starting values"); + } + this.numParameters = index[1] - index[0] + 1; + + double[][] start = null; + this.a = new double[numParameters]; + this.sigA = new double[numParameters]; + for (int i = 0; i < numParameters; i++) { + final String line = lines.get(index[0] + i - 1); + final String[] tokens = line.trim().split(" ++"); + if (start == null) { + start = new double[tokens.length - 4][numParameters]; + } + for (int j = 2; j < tokens.length - 2; j++) { + start[j - 2][i] = Double.parseDouble(tokens[j]); + } + this.a[i] = Double.parseDouble(tokens[tokens.length - 2]); + this.sigA[i] = Double.parseDouble(tokens[tokens.length - 1]); + } + if (start == null) { + throw new IOException("could not find starting values"); + } + this.numStartingPoints = start.length; + this.startingValues = start; + + double dummyDouble = Double.NaN; + String dummyString = null; + for (String line : lines) { + if (line.contains("Dataset Name:")) { + dummyString = line + .substring(line.indexOf("Dataset Name:") + 13, + line.indexOf("(")).trim(); + } + if (line.contains("Residual Sum of Squares")) { + final String[] tokens = line.split(" ++"); + dummyDouble = Double.parseDouble(tokens[4].trim()); + } + } + if (Double.isNaN(dummyDouble)) { + throw new IOException( + "could not find certified value of residual sum of squares"); + } + this.residualSumOfSquares = dummyDouble; + + if (dummyString == null) { + throw new IOException("could not find dataset name"); + } + this.name = dummyString; + + this.problem = new LeastSquaresProblem(); + } + + class LeastSquaresProblem { + public MultivariateVectorFunction getModelFunction() { + return new MultivariateVectorFunction() { + public double[] value(final double[] a) { + final int n = getNumObservations(); + final double[] yhat = new double[n]; + for (int i = 0; i < n; i++) { + yhat[i] = getModelValue(getX(i), a); + } + return yhat; + } + }; + } + + public MultivariateMatrixFunction getModelFunctionJacobian() { + return new MultivariateMatrixFunction() { + public double[][] value(final double[] a) + throws IllegalArgumentException { + final int n = getNumObservations(); + final double[][] j = new double[n][]; + for (int i = 0; i < n; i++) { + j[i] = getModelDerivatives(getX(i), a); + } + return j; + } + }; + } + } + + /** + * Returns the name of this dataset. + * + * @return the name of the dataset + */ + public String getName() { + return name; + } + + /** + * Returns the total number of observations (data points). + * + * @return the number of observations + */ + public int getNumObservations() { + return numObservations; + } + + /** + * Returns a copy of the data arrays. The data is laid out as follows
    6. + * {@code data[0][i] = x[i]},
    7. {@code data[1][i] = y[i]},
      8. + * + * @return the array of data points. + */ + public double[][] getData() { + return new double[][] { + MathArrays.copyOf(x), MathArrays.copyOf(y) + }; + } + + /** + * Returns the x-value of the {@code i}-th data point. + * + * @param i the index of the data point + * @return the x-value + */ + public double getX(final int i) { + return x[i]; + } + + /** + * Returns the y-value of the {@code i}-th data point. + * + * @param i the index of the data point + * @return the y-value + */ + public double getY(final int i) { + return y[i]; + } + + /** + * Returns the total number of parameters. + * + * @return the number of parameters + */ + public int getNumParameters() { + return numParameters; + } + + /** + * Returns the certified values of the paramters. + * + * @return the values of the parameters + */ + public double[] getParameters() { + return MathArrays.copyOf(a); + } + + /** + * Returns the certified value of the {@code i}-th parameter. + * + * @param i the index of the parameter + * @return the value of the parameter + */ + public double getParameter(final int i) { + return a[i]; + } + + /** + * Reurns the certified values of the standard deviations of the parameters. + * + * @return the standard deviations of the parameters + */ + public double[] getParametersStandardDeviations() { + return MathArrays.copyOf(sigA); + } + + /** + * Returns the certified value of the standard deviation of the {@code i}-th + * parameter. + * + * @param i the index of the parameter + * @return the standard deviation of the parameter + */ + public double getParameterStandardDeviation(final int i) { + return sigA[i]; + } + + /** + * Returns the certified value of the residual sum of squares. + * + * @return the residual sum of squares + */ + public double getResidualSumOfSquares() { + return residualSumOfSquares; + } + + /** + * Returns the total number of starting points (initial guesses for the + * optimization process). + * + * @return the number of starting points + */ + public int getNumStartingPoints() { + return numStartingPoints; + } + + /** + * Returns the {@code i}-th set of initial values of the parameters. + * + * @param i the index of the starting point + * @return the starting point + */ + public double[] getStartingPoint(final int i) { + return MathArrays.copyOf(startingValues[i]); + } + + /** + * Returns the least-squares problem corresponding to fitting the model to + * the specified data. + * + * @return the least-squares problem + */ + public LeastSquaresProblem getLeastSquaresProblem() { + return problem; + } + + /** + * Returns the value of the model for the specified values of the predictor + * variable and the parameters. + * + * @param x the predictor variable + * @param a the parameters + * @return the value of the model + */ + public abstract double getModelValue(final double x, final double[] a); + + /** + * Returns the values of the partial derivatives of the model with respect + * to the parameters. + * + * @param x the predictor variable + * @param a the parameters + * @return the partial derivatives + */ + public abstract double[] getModelDerivatives(final double x, + final double[] a); + + /** + *

      + * Parses the specified text lines, and extracts the indices of the first + * and last lines of the data defined by the specified {@code key}. This key + * must be one of + *

      + *
        + *
      • {@code "Starting Values"},
        • + *
      • {@code "Certified Values"},
        • + *
      • {@code "Data"}.
        • + *
      + *

      + * In the NIST data files, the line indices are separated by the keywords + * {@code "lines"} and {@code "to"}. + *

      + * + * @param lines the line of text to be parsed + * @return an array of two {@code int}s. First value is the index of the + * first line, second value is the index of the last line. + * {@code null} if the line could not be parsed. + */ + private static int[] findLineNumbers(final String key, + final Iterable lines) { + for (String text : lines) { + boolean flag = text.contains(key) && text.contains("lines") && + text.contains("to") && text.contains(")"); + if (flag) { + final int[] numbers = new int[2]; + final String from = text.substring(text.indexOf("lines") + 5, + text.indexOf("to")); + numbers[0] = Integer.parseInt(from.trim()); + final String to = text.substring(text.indexOf("to") + 2, + text.indexOf(")")); + numbers[1] = Integer.parseInt(to.trim()); + return numbers; + } + } + return null; + } +} diff --git a/src/test/java/org/apache/commons/math3/fitting/leastsquares/StatisticalReferenceDatasetFactory.java b/src/test/java/org/apache/commons/math3/fitting/leastsquares/StatisticalReferenceDatasetFactory.java new file mode 100644 index 000000000..cca7c989c --- /dev/null +++ b/src/test/java/org/apache/commons/math3/fitting/leastsquares/StatisticalReferenceDatasetFactory.java @@ -0,0 +1,201 @@ +/* + * 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.math3.fitting.leastsquares; + +import java.io.BufferedReader; +import java.io.IOException; +import java.io.InputStream; +import java.io.InputStreamReader; +import org.apache.commons.math3.util.FastMath; + +/** + * A factory to create instances of {@link StatisticalReferenceDataset} from + * available resources. + */ +public class StatisticalReferenceDatasetFactory { + + private StatisticalReferenceDatasetFactory() { + // Do nothing + } + + /** + * Creates a new buffered reader from the specified resource name. + * + * @param name the name of the resource + * @return a buffered reader + * @throws IOException if an I/O error occured + */ + public static BufferedReader createBufferedReaderFromResource(final String name) + throws IOException { + final InputStream resourceAsStream; + resourceAsStream = StatisticalReferenceDatasetFactory.class + .getResourceAsStream(name); + if (resourceAsStream == null) { + throw new IOException("could not find resource " + name); + } + return new BufferedReader(new InputStreamReader(resourceAsStream)); + } + + public static StatisticalReferenceDataset createKirby2() + throws IOException { + final BufferedReader in = createBufferedReaderFromResource("Kirby2.dat"); + StatisticalReferenceDataset dataset = null; + try { + dataset = new StatisticalReferenceDataset(in) { + + @Override + public double getModelValue(final double x, final double[] a) { + final double p = a[0] + x * (a[1] + x * a[2]); + final double q = 1.0 + x * (a[3] + x * a[4]); + return p / q; + } + + @Override + public double[] getModelDerivatives(final double x, + final double[] a) { + final double[] dy = new double[5]; + final double p = a[0] + x * (a[1] + x * a[2]); + final double q = 1.0 + x * (a[3] + x * a[4]); + dy[0] = 1.0 / q; + dy[1] = x / q; + dy[2] = x * dy[1]; + dy[3] = -x * p / (q * q); + dy[4] = x * dy[3]; + return dy; + } + }; + } finally { + in.close(); + } + return dataset; + } + + public static StatisticalReferenceDataset createHahn1() + throws IOException { + final BufferedReader in = createBufferedReaderFromResource("Hahn1.dat"); + StatisticalReferenceDataset dataset = null; + try { + dataset = new StatisticalReferenceDataset(in) { + + @Override + public double getModelValue(final double x, final double[] a) { + final double p = a[0] + x * (a[1] + x * (a[2] + x * a[3])); + final double q = 1.0 + x * (a[4] + x * (a[5] + x * a[6])); + return p / q; + } + + @Override + public double[] getModelDerivatives(final double x, + final double[] a) { + final double[] dy = new double[7]; + final double p = a[0] + x * (a[1] + x * (a[2] + x * a[3])); + final double q = 1.0 + x * (a[4] + x * (a[5] + x * a[6])); + dy[0] = 1.0 / q; + dy[1] = x * dy[0]; + dy[2] = x * dy[1]; + dy[3] = x * dy[2]; + dy[4] = -x * p / (q * q); + dy[5] = x * dy[4]; + dy[6] = x * dy[5]; + return dy; + } + }; + } finally { + in.close(); + } + return dataset; + } + + public static StatisticalReferenceDataset createMGH17() + throws IOException { + final BufferedReader in = createBufferedReaderFromResource("MGH17.dat"); + StatisticalReferenceDataset dataset = null; + try { + dataset = new StatisticalReferenceDataset(in) { + + @Override + public double getModelValue(final double x, final double[] a) { + return a[0] + a[1] * FastMath.exp(-a[3] * x) + a[2] * + FastMath.exp(-a[4] * x); + } + + @Override + public double[] getModelDerivatives(final double x, + final double[] a) { + final double[] dy = new double[5]; + dy[0] = 1.0; + dy[1] = FastMath.exp(-x * a[3]); + dy[2] = FastMath.exp(-x * a[4]); + dy[3] = -x * a[1] * dy[1]; + dy[4] = -x * a[2] * dy[2]; + return dy; + } + }; + } finally { + in.close(); + } + return dataset; + } + + public static StatisticalReferenceDataset createLanczos1() + throws IOException { + final BufferedReader in = + createBufferedReaderFromResource("Lanczos1.dat"); + StatisticalReferenceDataset dataset = null; + try { + dataset = new StatisticalReferenceDataset(in) { + + @Override + public double getModelValue(final double x, final double[] a) { + System.out.println(a[0]+", "+a[1]+", "+a[2]+", "+a[3]+", "+a[4]+", "+a[5]); + return a[0] * FastMath.exp(-a[3] * x) + + a[1] * FastMath.exp(-a[4] * x) + + a[2] * FastMath.exp(-a[5] * x); + } + + @Override + public double[] getModelDerivatives(final double x, + final double[] a) { + final double[] dy = new double[6]; + dy[0] = FastMath.exp(-x * a[3]); + dy[1] = FastMath.exp(-x * a[4]); + dy[2] = FastMath.exp(-x * a[5]); + dy[3] = -x * a[0] * dy[0]; + dy[4] = -x * a[1] * dy[1]; + dy[5] = -x * a[2] * dy[2]; + return dy; + } + }; + } finally { + in.close(); + } + return dataset; + } + + /** + * Returns an array with all available reference datasets. + * + * @return the array of datasets + * @throws IOException if an I/O error occurs + */ + public StatisticalReferenceDataset[] createAll() + throws IOException { + return new StatisticalReferenceDataset[] { + createKirby2(), createMGH17() + }; + } +} diff --git a/src/test/java/org/apache/commons/math3/fitting/leastsquares/StraightLineProblem.java b/src/test/java/org/apache/commons/math3/fitting/leastsquares/StraightLineProblem.java new file mode 100644 index 000000000..f17c4616c --- /dev/null +++ b/src/test/java/org/apache/commons/math3/fitting/leastsquares/StraightLineProblem.java @@ -0,0 +1,165 @@ +/* + * 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.math3.fitting.leastsquares; + +import java.util.ArrayList; +import org.apache.commons.math3.analysis.MultivariateVectorFunction; +import org.apache.commons.math3.analysis.MultivariateMatrixFunction; +import org.apache.commons.math3.analysis.UnivariateFunction; +import org.apache.commons.math3.stat.regression.SimpleRegression; + +/** + * Class that models a straight line defined as {@code y = a x + b}. + * The parameters of problem are: + *
        + *
      • {@code a}
        • + *
      • {@code b}
        • + *
      + * The model functions are: + *
        + *
      • for each pair (a, b), the y-coordinate of the line.
        • + *
      + */ +class StraightLineProblem { + /** Cloud of points assumed to be fitted by a straight line. */ + private final ArrayList points; + /** Error (on the y-coordinate of the points). */ + private final double sigma; + + /** + * @param error Assumed error for the y-coordinate. + */ + public StraightLineProblem(double error) { + points = new ArrayList(); + sigma = error; + } + + public void addPoint(double px, double py) { + points.add(new double[] { px, py }); + } + + /** + * @return the list of x-coordinates. + */ + public double[] x() { + final double[] v = new double[points.size()]; + for (int i = 0; i < points.size(); i++) { + final double[] p = points.get(i); + v[i] = p[0]; // x-coordinate. + } + + return v; + } + + /** + * @return the list of y-coordinates. + */ + public double[] y() { + final double[] v = new double[points.size()]; + for (int i = 0; i < points.size(); i++) { + final double[] p = points.get(i); + v[i] = p[1]; // y-coordinate. + } + + return v; + } + + public double[] target() { + return y(); + } + + public double[] weight() { + final double weight = 1 / (sigma * sigma); + final double[] w = new double[points.size()]; + for (int i = 0; i < points.size(); i++) { + w[i] = weight; + } + + return w; + } + + public MultivariateVectorFunction getModelFunction() { + return new MultivariateVectorFunction() { + public double[] value(double[] params) { + final Model line = new Model(params[0], params[1]); + + final double[] model = new double[points.size()]; + for (int i = 0; i < points.size(); i++) { + final double[] p = points.get(i); + model[i] = line.value(p[0]); + } + + return model; + } + }; + } + + public MultivariateMatrixFunction getModelFunctionJacobian() { + return new MultivariateMatrixFunction() { + public double[][] value(double[] point) { + return jacobian(point); + } + }; + } + + /** + * Directly solve the linear problem, using the {@link SimpleRegression} + * class. + */ + public double[] solve() { + final SimpleRegression regress = new SimpleRegression(true); + for (double[] d : points) { + regress.addData(d[0], d[1]); + } + + final double[] result = { regress.getSlope(), regress.getIntercept() }; + return result; + } + + private double[][] jacobian(double[] params) { + final double[][] jacobian = new double[points.size()][2]; + + for (int i = 0; i < points.size(); i++) { + final double[] p = points.get(i); + // Partial derivative wrt "a". + jacobian[i][0] = p[0]; + // Partial derivative wrt "b". + jacobian[i][1] = 1; + } + + return jacobian; + } + + /** + * Linear function. + */ + public static class Model implements UnivariateFunction { + final double a; + final double b; + + public Model(double a, + double b) { + this.a = a; + this.b = b; + } + + public double value(double x) { + return a * x + b; + } + } +} diff --git a/src/test/resources/org/apache/commons/math3/fitting/leastsquares/Hahn1.dat b/src/test/resources/org/apache/commons/math3/fitting/leastsquares/Hahn1.dat new file mode 100644 index 000000000..0e493a4e3 --- /dev/null +++ b/src/test/resources/org/apache/commons/math3/fitting/leastsquares/Hahn1.dat @@ -0,0 +1,296 @@ +NIST/ITL StRD +Dataset Name: Hahn1 (Hahn1.dat) + +File Format: ASCII + Starting Values (lines 41 to 47) + Certified Values (lines 41 to 52) + Data (lines 61 to 296) + +Procedure: Nonlinear Least Squares Regression + +Description: These data are the result of a NIST study involving + the thermal expansion of copper. The response + variable is the coefficient of thermal expansion, and + the predictor variable is temperature in degrees + kelvin. + + +Reference: Hahn, T., NIST (197?). + Copper Thermal Expansion Study. + + + + + +Data: 1 Response (y = coefficient of thermal expansion) + 1 Predictor (x = temperature, degrees kelvin) + 236 Observations + Average Level of Difficulty + Observed Data + +Model: Rational Class (cubic/cubic) + 7 Parameters (b1 to b7) + + y = (b1+b2*x+b3*x**2+b4*x**3) / + (1+b5*x+b6*x**2+b7*x**3) + e + + + Starting values Certified Values + + Start 1 Start 2 Parameter Standard Deviation + b1 = 10 1 1.0776351733E+00 1.7070154742E-01 + b2 = -1 -0.1 -1.2269296921E-01 1.2000289189E-02 + b3 = 0.05 0.005 4.0863750610E-03 2.2508314937E-04 + b4 = -0.00001 -0.000001 -1.4262662514E-06 2.7578037666E-07 + b5 = -0.05 -0.005 -5.7609940901E-03 2.4712888219E-04 + b6 = 0.001 0.0001 2.4053735503E-04 1.0449373768E-05 + b7 = -0.000001 -0.0000001 -1.2314450199E-07 1.3027335327E-08 + +Residual Sum of Squares: 1.5324382854E+00 +Residual Standard Deviation: 8.1803852243E-02 +Degrees of Freedom: 229 +Number of Observations: 236 + + + + + + + +Data: y x + .591E0 24.41E0 + 1.547E0 34.82E0 + 2.902E0 44.09E0 + 2.894E0 45.07E0 + 4.703E0 54.98E0 + 6.307E0 65.51E0 + 7.03E0 70.53E0 + 7.898E0 75.70E0 + 9.470E0 89.57E0 + 9.484E0 91.14E0 + 10.072E0 96.40E0 + 10.163E0 97.19E0 + 11.615E0 114.26E0 + 12.005E0 120.25E0 + 12.478E0 127.08E0 + 12.982E0 133.55E0 + 12.970E0 133.61E0 + 13.926E0 158.67E0 + 14.452E0 172.74E0 + 14.404E0 171.31E0 + 15.190E0 202.14E0 + 15.550E0 220.55E0 + 15.528E0 221.05E0 + 15.499E0 221.39E0 + 16.131E0 250.99E0 + 16.438E0 268.99E0 + 16.387E0 271.80E0 + 16.549E0 271.97E0 + 16.872E0 321.31E0 + 16.830E0 321.69E0 + 16.926E0 330.14E0 + 16.907E0 333.03E0 + 16.966E0 333.47E0 + 17.060E0 340.77E0 + 17.122E0 345.65E0 + 17.311E0 373.11E0 + 17.355E0 373.79E0 + 17.668E0 411.82E0 + 17.767E0 419.51E0 + 17.803E0 421.59E0 + 17.765E0 422.02E0 + 17.768E0 422.47E0 + 17.736E0 422.61E0 + 17.858E0 441.75E0 + 17.877E0 447.41E0 + 17.912E0 448.7E0 + 18.046E0 472.89E0 + 18.085E0 476.69E0 + 18.291E0 522.47E0 + 18.357E0 522.62E0 + 18.426E0 524.43E0 + 18.584E0 546.75E0 + 18.610E0 549.53E0 + 18.870E0 575.29E0 + 18.795E0 576.00E0 + 19.111E0 625.55E0 + .367E0 20.15E0 + .796E0 28.78E0 + 0.892E0 29.57E0 + 1.903E0 37.41E0 + 2.150E0 39.12E0 + 3.697E0 50.24E0 + 5.870E0 61.38E0 + 6.421E0 66.25E0 + 7.422E0 73.42E0 + 9.944E0 95.52E0 + 11.023E0 107.32E0 + 11.87E0 122.04E0 + 12.786E0 134.03E0 + 14.067E0 163.19E0 + 13.974E0 163.48E0 + 14.462E0 175.70E0 + 14.464E0 179.86E0 + 15.381E0 211.27E0 + 15.483E0 217.78E0 + 15.59E0 219.14E0 + 16.075E0 262.52E0 + 16.347E0 268.01E0 + 16.181E0 268.62E0 + 16.915E0 336.25E0 + 17.003E0 337.23E0 + 16.978E0 339.33E0 + 17.756E0 427.38E0 + 17.808E0 428.58E0 + 17.868E0 432.68E0 + 18.481E0 528.99E0 + 18.486E0 531.08E0 + 19.090E0 628.34E0 + 16.062E0 253.24E0 + 16.337E0 273.13E0 + 16.345E0 273.66E0 + 16.388E0 282.10E0 + 17.159E0 346.62E0 + 17.116E0 347.19E0 + 17.164E0 348.78E0 + 17.123E0 351.18E0 + 17.979E0 450.10E0 + 17.974E0 450.35E0 + 18.007E0 451.92E0 + 17.993E0 455.56E0 + 18.523E0 552.22E0 + 18.669E0 553.56E0 + 18.617E0 555.74E0 + 19.371E0 652.59E0 + 19.330E0 656.20E0 + 0.080E0 14.13E0 + 0.248E0 20.41E0 + 1.089E0 31.30E0 + 1.418E0 33.84E0 + 2.278E0 39.70E0 + 3.624E0 48.83E0 + 4.574E0 54.50E0 + 5.556E0 60.41E0 + 7.267E0 72.77E0 + 7.695E0 75.25E0 + 9.136E0 86.84E0 + 9.959E0 94.88E0 + 9.957E0 96.40E0 + 11.600E0 117.37E0 + 13.138E0 139.08E0 + 13.564E0 147.73E0 + 13.871E0 158.63E0 + 13.994E0 161.84E0 + 14.947E0 192.11E0 + 15.473E0 206.76E0 + 15.379E0 209.07E0 + 15.455E0 213.32E0 + 15.908E0 226.44E0 + 16.114E0 237.12E0 + 17.071E0 330.90E0 + 17.135E0 358.72E0 + 17.282E0 370.77E0 + 17.368E0 372.72E0 + 17.483E0 396.24E0 + 17.764E0 416.59E0 + 18.185E0 484.02E0 + 18.271E0 495.47E0 + 18.236E0 514.78E0 + 18.237E0 515.65E0 + 18.523E0 519.47E0 + 18.627E0 544.47E0 + 18.665E0 560.11E0 + 19.086E0 620.77E0 + 0.214E0 18.97E0 + 0.943E0 28.93E0 + 1.429E0 33.91E0 + 2.241E0 40.03E0 + 2.951E0 44.66E0 + 3.782E0 49.87E0 + 4.757E0 55.16E0 + 5.602E0 60.90E0 + 7.169E0 72.08E0 + 8.920E0 85.15E0 + 10.055E0 97.06E0 + 12.035E0 119.63E0 + 12.861E0 133.27E0 + 13.436E0 143.84E0 + 14.167E0 161.91E0 + 14.755E0 180.67E0 + 15.168E0 198.44E0 + 15.651E0 226.86E0 + 15.746E0 229.65E0 + 16.216E0 258.27E0 + 16.445E0 273.77E0 + 16.965E0 339.15E0 + 17.121E0 350.13E0 + 17.206E0 362.75E0 + 17.250E0 371.03E0 + 17.339E0 393.32E0 + 17.793E0 448.53E0 + 18.123E0 473.78E0 + 18.49E0 511.12E0 + 18.566E0 524.70E0 + 18.645E0 548.75E0 + 18.706E0 551.64E0 + 18.924E0 574.02E0 + 19.1E0 623.86E0 + 0.375E0 21.46E0 + 0.471E0 24.33E0 + 1.504E0 33.43E0 + 2.204E0 39.22E0 + 2.813E0 44.18E0 + 4.765E0 55.02E0 + 9.835E0 94.33E0 + 10.040E0 96.44E0 + 11.946E0 118.82E0 + 12.596E0 128.48E0 + 13.303E0 141.94E0 + 13.922E0 156.92E0 + 14.440E0 171.65E0 + 14.951E0 190.00E0 + 15.627E0 223.26E0 + 15.639E0 223.88E0 + 15.814E0 231.50E0 + 16.315E0 265.05E0 + 16.334E0 269.44E0 + 16.430E0 271.78E0 + 16.423E0 273.46E0 + 17.024E0 334.61E0 + 17.009E0 339.79E0 + 17.165E0 349.52E0 + 17.134E0 358.18E0 + 17.349E0 377.98E0 + 17.576E0 394.77E0 + 17.848E0 429.66E0 + 18.090E0 468.22E0 + 18.276E0 487.27E0 + 18.404E0 519.54E0 + 18.519E0 523.03E0 + 19.133E0 612.99E0 + 19.074E0 638.59E0 + 19.239E0 641.36E0 + 19.280E0 622.05E0 + 19.101E0 631.50E0 + 19.398E0 663.97E0 + 19.252E0 646.9E0 + 19.89E0 748.29E0 + 20.007E0 749.21E0 + 19.929E0 750.14E0 + 19.268E0 647.04E0 + 19.324E0 646.89E0 + 20.049E0 746.9E0 + 20.107E0 748.43E0 + 20.062E0 747.35E0 + 20.065E0 749.27E0 + 19.286E0 647.61E0 + 19.972E0 747.78E0 + 20.088E0 750.51E0 + 20.743E0 851.37E0 + 20.83E0 845.97E0 + 20.935E0 847.54E0 + 21.035E0 849.93E0 + 20.93E0 851.61E0 + 21.074E0 849.75E0 + 21.085E0 850.98E0 + 20.935E0 848.23E0 diff --git a/src/test/resources/org/apache/commons/math3/fitting/leastsquares/Kirby2.dat b/src/test/resources/org/apache/commons/math3/fitting/leastsquares/Kirby2.dat new file mode 100644 index 000000000..75cd80fb6 --- /dev/null +++ b/src/test/resources/org/apache/commons/math3/fitting/leastsquares/Kirby2.dat @@ -0,0 +1,211 @@ +NIST/ITL StRD +Dataset Name: Kirby2 (Kirby2.dat) + +File Format: ASCII + Starting Values (lines 41 to 45) + Certified Values (lines 41 to 50) + Data (lines 61 to 211) + +Procedure: Nonlinear Least Squares Regression + +Description: These data are the result of a NIST study involving + scanning electron microscope line with standards. + + +Reference: Kirby, R., NIST (197?). + Scanning electron microscope line width standards. + + + + + + + + +Data: 1 Response (y) + 1 Predictor (x) + 151 Observations + Average Level of Difficulty + Observed Data + +Model: Rational Class (quadratic/quadratic) + 5 Parameters (b1 to b5) + + y = (b1 + b2*x + b3*x**2) / + (1 + b4*x + b5*x**2) + e + + + Starting values Certified Values + + Start 1 Start 2 Parameter Standard Deviation + b1 = 2 1.5 1.6745063063E+00 8.7989634338E-02 + b2 = -0.1 -0.15 -1.3927397867E-01 4.1182041386E-03 + b3 = 0.003 0.0025 2.5961181191E-03 4.1856520458E-05 + b4 = -0.001 -0.0015 -1.7241811870E-03 5.8931897355E-05 + b5 = 0.00001 0.00002 2.1664802578E-05 2.0129761919E-07 + +Residual Sum of Squares: 3.9050739624E+00 +Residual Standard Deviation: 1.6354535131E-01 +Degrees of Freedom: 146 +Number of Observations: 151 + + + + + + + + + +Data: y x + 0.0082E0 9.65E0 + 0.0112E0 10.74E0 + 0.0149E0 11.81E0 + 0.0198E0 12.88E0 + 0.0248E0 14.06E0 + 0.0324E0 15.28E0 + 0.0420E0 16.63E0 + 0.0549E0 18.19E0 + 0.0719E0 19.88E0 + 0.0963E0 21.84E0 + 0.1291E0 24.00E0 + 0.1710E0 26.25E0 + 0.2314E0 28.86E0 + 0.3227E0 31.85E0 + 0.4809E0 35.79E0 + 0.7084E0 40.18E0 + 1.0220E0 44.74E0 + 1.4580E0 49.53E0 + 1.9520E0 53.94E0 + 2.5410E0 58.29E0 + 3.2230E0 62.63E0 + 3.9990E0 67.03E0 + 4.8520E0 71.25E0 + 5.7320E0 75.22E0 + 6.7270E0 79.33E0 + 7.8350E0 83.56E0 + 9.0250E0 87.75E0 + 10.2670E0 91.93E0 + 11.5780E0 96.10E0 + 12.9440E0 100.28E0 + 14.3770E0 104.46E0 + 15.8560E0 108.66E0 + 17.3310E0 112.71E0 + 18.8850E0 116.88E0 + 20.5750E0 121.33E0 + 22.3200E0 125.79E0 + 22.3030E0 125.79E0 + 23.4600E0 128.74E0 + 24.0600E0 130.27E0 + 25.2720E0 133.33E0 + 25.8530E0 134.79E0 + 27.1100E0 137.93E0 + 27.6580E0 139.33E0 + 28.9240E0 142.46E0 + 29.5110E0 143.90E0 + 30.7100E0 146.91E0 + 31.3500E0 148.51E0 + 32.5200E0 151.41E0 + 33.2300E0 153.17E0 + 34.3300E0 155.97E0 + 35.0600E0 157.76E0 + 36.1700E0 160.56E0 + 36.8400E0 162.30E0 + 38.0100E0 165.21E0 + 38.6700E0 166.90E0 + 39.8700E0 169.92E0 + 40.0300E0 170.32E0 + 40.5000E0 171.54E0 + 41.3700E0 173.79E0 + 41.6700E0 174.57E0 + 42.3100E0 176.25E0 + 42.7300E0 177.34E0 + 43.4600E0 179.19E0 + 44.1400E0 181.02E0 + 44.5500E0 182.08E0 + 45.2200E0 183.88E0 + 45.9200E0 185.75E0 + 46.3000E0 186.80E0 + 47.0000E0 188.63E0 + 47.6800E0 190.45E0 + 48.0600E0 191.48E0 + 48.7400E0 193.35E0 + 49.4100E0 195.22E0 + 49.7600E0 196.23E0 + 50.4300E0 198.05E0 + 51.1100E0 199.97E0 + 51.5000E0 201.06E0 + 52.1200E0 202.83E0 + 52.7600E0 204.69E0 + 53.1800E0 205.86E0 + 53.7800E0 207.58E0 + 54.4600E0 209.50E0 + 54.8300E0 210.65E0 + 55.4000E0 212.33E0 + 56.4300E0 215.43E0 + 57.0300E0 217.16E0 + 58.0000E0 220.21E0 + 58.6100E0 221.98E0 + 59.5800E0 225.06E0 + 60.1100E0 226.79E0 + 61.1000E0 229.92E0 + 61.6500E0 231.69E0 + 62.5900E0 234.77E0 + 63.1200E0 236.60E0 + 64.0300E0 239.63E0 + 64.6200E0 241.50E0 + 65.4900E0 244.48E0 + 66.0300E0 246.40E0 + 66.8900E0 249.35E0 + 67.4200E0 251.32E0 + 68.2300E0 254.22E0 + 68.7700E0 256.24E0 + 69.5900E0 259.11E0 + 70.1100E0 261.18E0 + 70.8600E0 264.02E0 + 71.4300E0 266.13E0 + 72.1600E0 268.94E0 + 72.7000E0 271.09E0 + 73.4000E0 273.87E0 + 73.9300E0 276.08E0 + 74.6000E0 278.83E0 + 75.1600E0 281.08E0 + 75.8200E0 283.81E0 + 76.3400E0 286.11E0 + 76.9800E0 288.81E0 + 77.4800E0 291.08E0 + 78.0800E0 293.75E0 + 78.6000E0 295.99E0 + 79.1700E0 298.64E0 + 79.6200E0 300.84E0 + 79.8800E0 302.02E0 + 80.1900E0 303.48E0 + 80.6600E0 305.65E0 + 81.2200E0 308.27E0 + 81.6600E0 310.41E0 + 82.1600E0 313.01E0 + 82.5900E0 315.12E0 + 83.1400E0 317.71E0 + 83.5000E0 319.79E0 + 84.0000E0 322.36E0 + 84.4000E0 324.42E0 + 84.8900E0 326.98E0 + 85.2600E0 329.01E0 + 85.7400E0 331.56E0 + 86.0700E0 333.56E0 + 86.5400E0 336.10E0 + 86.8900E0 338.08E0 + 87.3200E0 340.60E0 + 87.6500E0 342.57E0 + 88.1000E0 345.08E0 + 88.4300E0 347.02E0 + 88.8300E0 349.52E0 + 89.1200E0 351.44E0 + 89.5400E0 353.93E0 + 89.8500E0 355.83E0 + 90.2500E0 358.32E0 + 90.5500E0 360.20E0 + 90.9300E0 362.67E0 + 91.2000E0 364.53E0 + 91.5500E0 367.00E0 + 92.2000E0 371.30E0 diff --git a/src/test/resources/org/apache/commons/math3/fitting/leastsquares/Lanczos1.dat b/src/test/resources/org/apache/commons/math3/fitting/leastsquares/Lanczos1.dat new file mode 100644 index 000000000..d23d5e483 --- /dev/null +++ b/src/test/resources/org/apache/commons/math3/fitting/leastsquares/Lanczos1.dat @@ -0,0 +1,84 @@ +NIST/ITL StRD +Dataset Name: Lanczos1 (Lanczos1.dat) + +File Format: ASCII + Starting Values (lines 41 to 46) + Certified Values (lines 41 to 51) + Data (lines 61 to 84) + +Procedure: Nonlinear Least Squares Regression + +Description: These data are taken from an example discussed in + Lanczos (1956). The data were generated to 14-digits + of accuracy using + f(x) = 0.0951*exp(-x) + 0.8607*exp(-3*x) + + 1.5576*exp(-5*x). + + +Reference: Lanczos, C. (1956). + Applied Analysis. + Englewood Cliffs, NJ: Prentice Hall, pp. 272-280. + + + + +Data: 1 Response (y) + 1 Predictor (x) + 24 Observations + Average Level of Difficulty + Generated Data + +Model: Exponential Class + 6 Parameters (b1 to b6) + + y = b1*exp(-b2*x) + b3*exp(-b4*x) + b5*exp(-b6*x) + e + + + + Starting values Certified Values + + Start 1 Start 2 Parameter Standard Deviation + b1 = 1.2 0.5 9.5100000027E-02 5.3347304234E-11 + b2 = 0.3 0.7 1.0000000001E+00 2.7473038179E-10 + b3 = 5.6 3.6 8.6070000013E-01 1.3576062225E-10 + b4 = 5.5 4.2 3.0000000002E+00 3.3308253069E-10 + b5 = 6.5 4 1.5575999998E+00 1.8815731448E-10 + b6 = 7.6 6.3 5.0000000001E+00 1.1057500538E-10 + +Residual Sum of Squares: 1.4307867721E-25 +Residual Standard Deviation: 8.9156129349E-14 +Degrees of Freedom: 18 +Number of Observations: 24 + + + + + + + + +Data: y x + 2.513400000000E+00 0.000000000000E+00 + 2.044333373291E+00 5.000000000000E-02 + 1.668404436564E+00 1.000000000000E-01 + 1.366418021208E+00 1.500000000000E-01 + 1.123232487372E+00 2.000000000000E-01 + 9.268897180037E-01 2.500000000000E-01 + 7.679338563728E-01 3.000000000000E-01 + 6.388775523106E-01 3.500000000000E-01 + 5.337835317402E-01 4.000000000000E-01 + 4.479363617347E-01 4.500000000000E-01 + 3.775847884350E-01 5.000000000000E-01 + 3.197393199326E-01 5.500000000000E-01 + 2.720130773746E-01 6.000000000000E-01 + 2.324965529032E-01 6.500000000000E-01 + 1.996589546065E-01 7.000000000000E-01 + 1.722704126914E-01 7.500000000000E-01 + 1.493405660168E-01 8.000000000000E-01 + 1.300700206922E-01 8.500000000000E-01 + 1.138119324644E-01 9.000000000000E-01 + 1.000415587559E-01 9.500000000000E-01 + 8.833209084540E-02 1.000000000000E+00 + 7.833544019350E-02 1.050000000000E+00 + 6.976693743449E-02 1.100000000000E+00 + 6.239312536719E-02 1.150000000000E+00 diff --git a/src/test/resources/org/apache/commons/math3/fitting/leastsquares/MGH17.dat b/src/test/resources/org/apache/commons/math3/fitting/leastsquares/MGH17.dat new file mode 100644 index 000000000..584f73c2f --- /dev/null +++ b/src/test/resources/org/apache/commons/math3/fitting/leastsquares/MGH17.dat @@ -0,0 +1,93 @@ +NIST/ITL StRD +Dataset Name: MGH17 (MGH17.dat) + +File Format: ASCII + Starting Values (lines 41 to 45) + Certified Values (lines 41 to 50) + Data (lines 61 to 93) + +Procedure: Nonlinear Least Squares Regression + +Description: This problem was found to be difficult for some very + good algorithms. + + See More, J. J., Garbow, B. S., and Hillstrom, K. E. + (1981). Testing unconstrained optimization software. + ACM Transactions on Mathematical Software. 7(1): + pp. 17-41. + +Reference: Osborne, M. R. (1972). + Some aspects of nonlinear least squares + calculations. In Numerical Methods for Nonlinear + Optimization, Lootsma (Ed). + New York, NY: Academic Press, pp. 171-189. + +Data: 1 Response (y) + 1 Predictor (x) + 33 Observations + Average Level of Difficulty + Generated Data + +Model: Exponential Class + 5 Parameters (b1 to b5) + + y = b1 + b2*exp[-x*b4] + b3*exp[-x*b5] + e + + + + Starting values Certified Values + + Start 1 Start 2 Parameter Standard Deviation + b1 = 50 0.5 3.7541005211E-01 2.0723153551E-03 + b2 = 150 1.5 1.9358469127E+00 2.2031669222E-01 + b3 = -100 -1 -1.4646871366E+00 2.2175707739E-01 + b4 = 1 0.01 1.2867534640E-02 4.4861358114E-04 + b5 = 2 0.02 2.2122699662E-02 8.9471996575E-04 + +Residual Sum of Squares: 5.4648946975E-05 +Residual Standard Deviation: 1.3970497866E-03 +Degrees of Freedom: 28 +Number of Observations: 33 + + + + + + + + + +Data: y x + 8.440000E-01 0.000000E+00 + 9.080000E-01 1.000000E+01 + 9.320000E-01 2.000000E+01 + 9.360000E-01 3.000000E+01 + 9.250000E-01 4.000000E+01 + 9.080000E-01 5.000000E+01 + 8.810000E-01 6.000000E+01 + 8.500000E-01 7.000000E+01 + 8.180000E-01 8.000000E+01 + 7.840000E-01 9.000000E+01 + 7.510000E-01 1.000000E+02 + 7.180000E-01 1.100000E+02 + 6.850000E-01 1.200000E+02 + 6.580000E-01 1.300000E+02 + 6.280000E-01 1.400000E+02 + 6.030000E-01 1.500000E+02 + 5.800000E-01 1.600000E+02 + 5.580000E-01 1.700000E+02 + 5.380000E-01 1.800000E+02 + 5.220000E-01 1.900000E+02 + 5.060000E-01 2.000000E+02 + 4.900000E-01 2.100000E+02 + 4.780000E-01 2.200000E+02 + 4.670000E-01 2.300000E+02 + 4.570000E-01 2.400000E+02 + 4.480000E-01 2.500000E+02 + 4.380000E-01 2.600000E+02 + 4.310000E-01 2.700000E+02 + 4.240000E-01 2.800000E+02 + 4.200000E-01 2.900000E+02 + 4.140000E-01 3.000000E+02 + 4.110000E-01 3.100000E+02 + 4.060000E-01 3.200000E+02