Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

MATH-425 

Removed usage of exceptions defined in the package "linear". Added
replacements in package "exception". Modified tests to catch the
new exception type. Deleted the deprecated package "estimation" (test
failed because of the changes).
Deprecated "MaxIterationsExceededException".
Changes in package "optimization.general" were necessary to comply with
the new exceptions and to be consistent in reporting the actual cause
of failure.


git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1033486 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Gilles Sadowski 2010-11-10 14:23:36 +00:00
parent 06ec87572a
commit 07f1663320
56 changed files with 495 additions and 5779 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/main/java/org/apache/commons/math/MaxIterationsExceededException.java
View File

@ -28,6 +28,9 @@ import org.apache.commons.math.exception.util.LocalizedFormats;
*
* @version $Revision$ $Date$
* @since 1.2
* @deprecated in 2.2 (to be removed in 3.0). Please use
* {@link org.apache.commons.math.exception.MaxCountExceededException}
* instead.
*/
public class MaxIterationsExceededException extends ConvergenceException {

318
src/main/java/org/apache/commons/math/estimation/AbstractEstimator.java
View File

@ -1,318 +0,0 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.estimation;
import java.util.Arrays;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.linear.InvalidMatrixException;
import org.apache.commons.math.linear.LUDecompositionImpl;
import org.apache.commons.math.linear.MatrixUtils;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.util.FastMath;
/**
* Base class for implementing estimators.
* <p>This base class handles the boilerplates methods associated to thresholds
* settings, jacobian and error estimation.</p>
* @version $Revision$ $Date$
* @since 1.2
* @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
* been deprecated and replaced by package org.apache.commons.math.optimization.general
*
*/
@Deprecated
public abstract class AbstractEstimator implements Estimator {
/** Default maximal number of cost evaluations allowed. */
public static final int DEFAULT_MAX_COST_EVALUATIONS = 100;
/** Array of measurements. */
protected WeightedMeasurement[] measurements;
/** Array of parameters. */
protected EstimatedParameter[] parameters;
/**
* Jacobian matrix.
* <p>This matrix is in canonical form just after the calls to
* {@link #updateJacobian()}, but may be modified by the solver
* in the derived class (the {@link LevenbergMarquardtEstimator
* Levenberg-Marquardt estimator} does this).</p>
*/
protected double[] jacobian;
/** Number of columns of the jacobian matrix. */
protected int cols;
/** Number of rows of the jacobian matrix. */
protected int rows;
/** Residuals array.
* <p>This array is in canonical form just after the calls to
* {@link #updateJacobian()}, but may be modified by the solver
* in the derived class (the {@link LevenbergMarquardtEstimator
* Levenberg-Marquardt estimator} does this).</p>
*/
protected double[] residuals;
/** Cost value (square root of the sum of the residuals). */
protected double cost;
/** Maximal allowed number of cost evaluations. */
private int maxCostEval;
/** Number of cost evaluations. */
private int costEvaluations;
/** Number of jacobian evaluations. */
private int jacobianEvaluations;
/**
* Build an abstract estimator for least squares problems.
* <p>The maximal number of cost evaluations allowed is set
* to its default value {@link #DEFAULT_MAX_COST_EVALUATIONS}.</p>
*/
protected AbstractEstimator() {
setMaxCostEval(DEFAULT_MAX_COST_EVALUATIONS);
}
/**
* Set the maximal number of cost evaluations allowed.
*
* @param maxCostEval maximal number of cost evaluations allowed
* @see #estimate
*/
public final void setMaxCostEval(int maxCostEval) {
this.maxCostEval = maxCostEval;
}
/**
* Get the number of cost evaluations.
*
* @return number of cost evaluations
* */
public final int getCostEvaluations() {
return costEvaluations;
}
/**
* Get the number of jacobian evaluations.
*
* @return number of jacobian evaluations
* */
public final int getJacobianEvaluations() {
return jacobianEvaluations;
}
/**
* Update the jacobian matrix.
*/
protected void updateJacobian() {
incrementJacobianEvaluationsCounter();
Arrays.fill(jacobian, 0);
int index = 0;
for (int i = 0; i < rows; i++) {
WeightedMeasurement wm = measurements[i];
double factor = -FastMath.sqrt(wm.getWeight());
for (int j = 0; j < cols; ++j) {
jacobian[index++] = factor * wm.getPartial(parameters[j]);
}
}
}
/**
* Increment the jacobian evaluations counter.
*/
protected final void incrementJacobianEvaluationsCounter() {
++jacobianEvaluations;
}
/**
* Update the residuals array and cost function value.
* @exception EstimationException if the number of cost evaluations
* exceeds the maximum allowed
*/
protected void updateResidualsAndCost()
throws EstimationException {
if (++costEvaluations > maxCostEval) {
throw new EstimationException(LocalizedFormats.MAX_EVALUATIONS_EXCEEDED,
maxCostEval);
}
cost = 0;
int index = 0;
for (int i = 0; i < rows; i++, index += cols) {
WeightedMeasurement wm = measurements[i];
double residual = wm.getResidual();
residuals[i] = FastMath.sqrt(wm.getWeight()) * residual;
cost += wm.getWeight() * residual * residual;
}
cost = FastMath.sqrt(cost);
}
/**
* Get the Root Mean Square value.
* Get the Root Mean Square value, i.e. the root of the arithmetic
* mean of the square of all weighted residuals. This is related to the
* criterion that is minimized by the estimator as follows: if
* <em>c</em> if the criterion, and <em>n</em> is the number of
* measurements, then the RMS is <em>sqrt (c/n)</em>.
*
* @param problem estimation problem
* @return RMS value
*/
public double getRMS(EstimationProblem problem) {
WeightedMeasurement[] wm = problem.getMeasurements();
double criterion = 0;
for (int i = 0; i < wm.length; ++i) {
double residual = wm[i].getResidual();
criterion += wm[i].getWeight() * residual * residual;
}
return FastMath.sqrt(criterion / wm.length);
}
/**
* Get the Chi-Square value.
* @param problem estimation problem
* @return chi-square value
*/
public double getChiSquare(EstimationProblem problem) {
WeightedMeasurement[] wm = problem.getMeasurements();
double chiSquare = 0;
for (int i = 0; i < wm.length; ++i) {
double residual = wm[i].getResidual();
chiSquare += residual * residual / wm[i].getWeight();
}
return chiSquare;
}
/**
* Get the covariance matrix of unbound estimated parameters.
* @param problem estimation problem
* @return covariance matrix
* @exception EstimationException if the covariance matrix
* cannot be computed (singular problem)
*/
public double[][] getCovariances(EstimationProblem problem)
throws EstimationException {
// set up the jacobian
updateJacobian();
// compute transpose(J).J, avoiding building big intermediate matrices
final int n = problem.getMeasurements().length;
final int m = problem.getUnboundParameters().length;
final int max = m * n;
double[][] jTj = new double[m][m];
for (int i = 0; i < m; ++i) {
for (int j = i; j < m; ++j) {
double sum = 0;
for (int k = 0; k < max; k += m) {
sum += jacobian[k + i] * jacobian[k + j];
}
jTj[i][j] = sum;
jTj[j][i] = sum;
}
}
try {
// compute the covariances matrix
RealMatrix inverse =
new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
return inverse.getData();
} catch (InvalidMatrixException ime) {
throw new EstimationException(LocalizedFormats.UNABLE_TO_COMPUTE_COVARIANCE_SINGULAR_PROBLEM);
}
}
/**
* Guess the errors in unbound estimated parameters.
* <p>Guessing is covariance-based, it only gives rough order of magnitude.</p>
* @param problem estimation problem
* @return errors in estimated parameters
* @exception EstimationException if the covariances matrix cannot be computed
* or the number of degrees of freedom is not positive (number of measurements
* lesser or equal to number of parameters)
*/
public double[] guessParametersErrors(EstimationProblem problem)
throws EstimationException {
int m = problem.getMeasurements().length;
int p = problem.getUnboundParameters().length;
if (m <= p) {
throw new EstimationException(
LocalizedFormats.NO_DEGREES_OF_FREEDOM,
m, p);
}
double[] errors = new double[problem.getUnboundParameters().length];
final double c = FastMath.sqrt(getChiSquare(problem) / (m - p));
double[][] covar = getCovariances(problem);
for (int i = 0; i < errors.length; ++i) {
errors[i] = FastMath.sqrt(covar[i][i]) * c;
}
return errors;
}
/**
* Initialization of the common parts of the estimation.
* <p>This method <em>must</em> be called at the start
* of the {@link #estimate(EstimationProblem) estimate}
* method.</p>
* @param problem estimation problem to solve
*/
protected void initializeEstimate(EstimationProblem problem) {
// reset counters
costEvaluations = 0;
jacobianEvaluations = 0;
// retrieve the equations and the parameters
measurements = problem.getMeasurements();
parameters = problem.getUnboundParameters();
// arrays shared with the other private methods
rows = measurements.length;
cols = parameters.length;
jacobian = new double[rows * cols];
residuals = new double[rows];
cost = Double.POSITIVE_INFINITY;
}
/**
* Solve an estimation problem.
*
* <p>The method should set the parameters of the problem to several
* trial values until it reaches convergence. If this method returns
* normally (i.e. without throwing an exception), then the best
* estimate of the parameters can be retrieved from the problem
* itself, through the {@link EstimationProblem#getAllParameters
* EstimationProblem.getAllParameters} method.</p>
*
* @param problem estimation problem to solve
* @exception EstimationException if the problem cannot be solved
*
*/
public abstract void estimate(EstimationProblem problem)
throws EstimationException;
}

126
src/main/java/org/apache/commons/math/estimation/EstimatedParameter.java
View File

@ -1,126 +0,0 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.estimation;
import java.io.Serializable;
/** This class represents the estimated parameters of an estimation problem.
*
* <p>The parameters of an estimation problem have a name, a value and
* a bound flag. The value of bound parameters is considered trusted
* and the solvers should not adjust them. On the other hand, the
* solvers should adjust the value of unbounds parameters until they
* satisfy convergence criterions specific to each solver.</p>
*
* @version $Revision$ $Date$
* @since 1.2
* @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
* been deprecated and replaced by package org.apache.commons.math.optimization.general
*
*/
@Deprecated
public class EstimatedParameter
implements Serializable {
/** Serializable version identifier */
private static final long serialVersionUID = -555440800213416949L;
/** Current value of the parameter */
protected double estimate;
/** Name of the parameter */
private final String name;
/** Indicator for bound parameters
* (ie parameters that should not be estimated)
*/
private boolean bound;
/** Simple constructor.
* Build an instance from a first estimate of the parameter,
* initially considered unbound.
* @param name name of the parameter
* @param firstEstimate first estimate of the parameter
*/
public EstimatedParameter(String name, double firstEstimate) {
this.name = name;
estimate = firstEstimate;
bound = false;
}
/** Simple constructor.
* Build an instance from a first estimate of the parameter and a
* bound flag
* @param name name of the parameter
* @param firstEstimate first estimate of the parameter
* @param bound flag, should be true if the parameter is bound
*/
public EstimatedParameter(String name,
double firstEstimate,
boolean bound) {
this.name = name;
estimate = firstEstimate;
this.bound = bound;
}
/** Copy constructor.
* Build a copy of a parameter
* @param parameter instance to copy
*/
public EstimatedParameter(EstimatedParameter parameter) {
name = parameter.name;
estimate = parameter.estimate;
bound = parameter.bound;
}
/** Set a new estimated value for the parameter.
* @param estimate new estimate for the parameter
*/
public void setEstimate(double estimate) {
this.estimate = estimate;
}
/** Get the current estimate of the parameter
* @return current estimate
*/
public double getEstimate() {
return estimate;
}
/** get the name of the parameter
* @return parameter name
*/
public String getName() {
return name;
}
/** Set the bound flag of the parameter
* @param bound this flag should be set to true if the parameter is
* bound (i.e. if it should not be adjusted by the solver).
*/
public void setBound(boolean bound) {
this.bound = bound;
}
/** Check if the parameter is bound
* @return true if the parameter is bound */
public boolean isBound() {
return bound;
}
}

61
src/main/java/org/apache/commons/math/estimation/EstimationException.java
View File

@ -1,61 +0,0 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.estimation;
import org.apache.commons.math.MathException;
import org.apache.commons.math.exception.util.DummyLocalizable;
import org.apache.commons.math.exception.util.Localizable;
/**
* This class represents exceptions thrown by the estimation solvers.
*
* @version $Revision$ $Date$
* @since 1.2
* @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
* been deprecated and replaced by package org.apache.commons.math.optimization.general
*
*/
@Deprecated
public class EstimationException
extends MathException {
/** Serializable version identifier. */
private static final long serialVersionUID = -573038581493881337L;
/**
* Simple constructor.
* Build an exception by translating and formating a message
* @param specifier format specifier (to be translated)
* @param parts to insert in the format (no translation)
*/
public EstimationException(String specifier, Object ... parts) {
this(new DummyLocalizable(specifier), parts);
}
/**
* Simple constructor.
* Build an exception by translating and formating a message
* @param specifier format specifier (to be translated)
* @param parts to insert in the format (no translation)
* @since 2.2
*/
public EstimationException(Localizable specifier, Object ... parts) {
super(specifier, parts);
}
}

68
src/main/java/org/apache/commons/math/estimation/EstimationProblem.java
View File

@ -1,68 +0,0 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.estimation;
/**
* This interface represents an estimation problem.
*
* <p>This interface should be implemented by all real estimation
* problems before they can be handled by the estimators through the
* {@link Estimator#estimate Estimator.estimate} method.</p>
*
* <p>An estimation problem, as seen by a solver is a set of
* parameters and a set of measurements. The parameters are adjusted
* during the estimation through the {@link #getUnboundParameters
* getUnboundParameters} and {@link EstimatedParameter#setEstimate
* EstimatedParameter.setEstimate} methods. The measurements both have
* a measured value which is generally fixed at construction and a
* theoretical value which depends on the model and hence varies as
* the parameters are adjusted. The purpose of the solver is to reduce
* the residual between these values, it can retrieve the measurements
* through the {@link #getMeasurements getMeasurements} method.</p>
*
* @see Estimator
* @see WeightedMeasurement
*
* @version $Revision$ $Date$
* @since 1.2
* @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
* been deprecated and replaced by package org.apache.commons.math.optimization.general
*
*/
@Deprecated
public interface EstimationProblem {
/**
* Get the measurements of an estimation problem.
* @return measurements
*/
WeightedMeasurement[] getMeasurements();
/**
* Get the unbound parameters of the problem.
* @return unbound parameters
*/
EstimatedParameter[] getUnboundParameters();
/**
* Get all the parameters of the problem.
* @return parameters
*/
EstimatedParameter[] getAllParameters();
}

90
src/main/java/org/apache/commons/math/estimation/Estimator.java
View File

@ -1,90 +0,0 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.estimation;
/**
* This interface represents solvers for estimation problems.
*
* <p>The classes which are devoted to solve estimation problems
* should implement this interface. The problems which can be handled
* should implement the {@link EstimationProblem} interface which
* gather all the information needed by the solver.</p>
*
* <p>The interface is composed only of the {@link #estimate estimate}
* method.</p>
*
* @see EstimationProblem
*
* @version $Revision$ $Date$
* @since 1.2
* @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
* been deprecated and replaced by package org.apache.commons.math.optimization.general
*
*/
@Deprecated
public interface Estimator {
/**
* Solve an estimation problem.
*
* <p>The method should set the parameters of the problem to several
* trial values until it reaches convergence. If this method returns
* normally (i.e. without throwing an exception), then the best
* estimate of the parameters can be retrieved from the problem
* itself, through the {@link EstimationProblem#getAllParameters
* EstimationProblem.getAllParameters} method.</p>
*
* @param problem estimation problem to solve
* @exception EstimationException if the problem cannot be solved
*
*/
void estimate(EstimationProblem problem) throws EstimationException;
/**
* Get the Root Mean Square value.
* Get the Root Mean Square value, i.e. the root of the arithmetic
* mean of the square of all weighted residuals. This is related to the
* criterion that is minimized by the estimator as follows: if
* <em>c</em> is the criterion, and <em>n</em> is the number of
* measurements, then the RMS is <em>sqrt (c/n)</em>.
* @see #guessParametersErrors(EstimationProblem)
*
* @param problem estimation problem
* @return RMS value
*/
double getRMS(EstimationProblem problem);
/**
* Get the covariance matrix of estimated parameters.
* @param problem estimation problem
* @return covariance matrix
* @exception EstimationException if the covariance matrix
* cannot be computed (singular problem)
*/
double[][] getCovariances(EstimationProblem problem) throws EstimationException;
/**
* Guess the errors in estimated parameters.
* @see #getRMS(EstimationProblem)
* @param problem estimation problem
* @return errors in estimated parameters
* @exception EstimationException if the error cannot be guessed
*/
double[] guessParametersErrors(EstimationProblem problem) throws EstimationException;
}

231
src/main/java/org/apache/commons/math/estimation/GaussNewtonEstimator.java
View File

@ -1,231 +0,0 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.estimation;
import java.io.Serializable;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.linear.InvalidMatrixException;
import org.apache.commons.math.linear.LUDecompositionImpl;
import org.apache.commons.math.linear.MatrixUtils;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.linear.RealVector;
import org.apache.commons.math.linear.ArrayRealVector;
import org.apache.commons.math.util.FastMath;
/**
* This class implements a solver for estimation problems.
*
* <p>This class solves estimation problems using a weighted least
* squares criterion on the measurement residuals. It uses a
* Gauss-Newton algorithm.</p>
*
* @version $Revision$ $Date$
* @since 1.2
* @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
* been deprecated and replaced by package org.apache.commons.math.optimization.general
*
*/
@Deprecated
public class GaussNewtonEstimator extends AbstractEstimator implements Serializable {
/** Serializable version identifier */
private static final long serialVersionUID = 5485001826076289109L;
/** Default threshold for cost steady state detection. */
private static final double DEFAULT_STEADY_STATE_THRESHOLD = 1.0e-6;
/** Default threshold for cost convergence. */
private static final double DEFAULT_CONVERGENCE = 1.0e-6;
/** Threshold for cost steady state detection. */
private double steadyStateThreshold;
/** Threshold for cost convergence. */
private double convergence;
/** Simple constructor with default settings.
* <p>
* The estimator is built with default values for all settings.
* </p>
* @see #DEFAULT_STEADY_STATE_THRESHOLD
* @see #DEFAULT_CONVERGENCE
* @see AbstractEstimator#DEFAULT_MAX_COST_EVALUATIONS
*/
public GaussNewtonEstimator() {
this.steadyStateThreshold = DEFAULT_STEADY_STATE_THRESHOLD;
this.convergence = DEFAULT_CONVERGENCE;
}
/**
* Simple constructor.
*
* <p>This constructor builds an estimator and stores its convergence
* characteristics.</p>
*
* <p>An estimator is considered to have converged whenever either
* the criterion goes below a physical threshold under which
* improvements are considered useless or when the algorithm is
* unable to improve it (even if it is still high). The first
* condition that is met stops the iterations.</p>
*
* <p>The fact an estimator has converged does not mean that the
* model accurately fits the measurements. It only means no better
* solution can be found, it does not mean this one is good. Such an
* analysis is left to the caller.</p>
*
* <p>If neither conditions are fulfilled before a given number of
* iterations, the algorithm is considered to have failed and an
* {@link EstimationException} is thrown.</p>
*
* @param maxCostEval maximal number of cost evaluations allowed
* @param convergence criterion threshold below which we do not need
* to improve the criterion anymore
* @param steadyStateThreshold steady state detection threshold, the
* problem has converged has reached a steady state if
* <code>FastMath.abs(J<sub>n</sub> - J<sub>n-1</sub>) &lt;
* J<sub>n</sub> &times convergence</code>, where <code>J<sub>n</sub></code>
* and <code>J<sub>n-1</sub></code> are the current and preceding criterion
* values (square sum of the weighted residuals of considered measurements).
*/
public GaussNewtonEstimator(final int maxCostEval, final double convergence,
final double steadyStateThreshold) {
setMaxCostEval(maxCostEval);
this.steadyStateThreshold = steadyStateThreshold;
this.convergence = convergence;
}
/**
* Set the convergence criterion threshold.
* @param convergence criterion threshold below which we do not need
* to improve the criterion anymore
*/
public void setConvergence(final double convergence) {
this.convergence = convergence;
}
/**
* Set the steady state detection threshold.
* <p>
* The problem has converged has reached a steady state if
* <code>FastMath.abs(J<sub>n</sub> - J<sub>n-1</sub>) &lt;
* J<sub>n</sub> &times convergence</code>, where <code>J<sub>n</sub></code>
* and <code>J<sub>n-1</sub></code> are the current and preceding criterion
* values (square sum of the weighted residuals of considered measurements).
* </p>
* @param steadyStateThreshold steady state detection threshold
*/
public void setSteadyStateThreshold(final double steadyStateThreshold) {
this.steadyStateThreshold = steadyStateThreshold;
}
/**
* Solve an estimation problem using a least squares criterion.
*
* <p>This method set the unbound parameters of the given problem
* starting from their current values through several iterations. At
* each step, the unbound parameters are changed in order to
* minimize a weighted least square criterion based on the
* measurements of the problem.</p>
*
* <p>The iterations are stopped either when the criterion goes
* below a physical threshold under which improvement are considered
* useless or when the algorithm is unable to improve it (even if it
* is still high). The first condition that is met stops the
* iterations. If the convergence it not reached before the maximum
* number of iterations, an {@link EstimationException} is
* thrown.</p>
*
* @param problem estimation problem to solve
* @exception EstimationException if the problem cannot be solved
*
* @see EstimationProblem
*
*/
@Override
public void estimate(EstimationProblem problem)
throws EstimationException {
initializeEstimate(problem);
// work matrices
double[] grad = new double[parameters.length];
ArrayRealVector bDecrement = new ArrayRealVector(parameters.length);
double[] bDecrementData = bDecrement.getDataRef();
RealMatrix wGradGradT = MatrixUtils.createRealMatrix(parameters.length, parameters.length);
// iterate until convergence is reached
double previous = Double.POSITIVE_INFINITY;
do {
// build the linear problem
incrementJacobianEvaluationsCounter();
RealVector b = new ArrayRealVector(parameters.length);
RealMatrix a = MatrixUtils.createRealMatrix(parameters.length, parameters.length);
for (int i = 0; i < measurements.length; ++i) {
if (! measurements [i].isIgnored()) {
double weight = measurements[i].getWeight();
double residual = measurements[i].getResidual();
// compute the normal equation
for (int j = 0; j < parameters.length; ++j) {
grad[j] = measurements[i].getPartial(parameters[j]);
bDecrementData[j] = weight * residual * grad[j];
}
// build the contribution matrix for measurement i
for (int k = 0; k < parameters.length; ++k) {
double gk = grad[k];
for (int l = 0; l < parameters.length; ++l) {
wGradGradT.setEntry(k, l, weight * gk * grad[l]);
}
}
// update the matrices
a = a.add(wGradGradT);
b = b.add(bDecrement);
}
}
try {
// solve the linearized least squares problem
RealVector dX = new LUDecompositionImpl(a).getSolver().solve(b);
// update the estimated parameters
for (int i = 0; i < parameters.length; ++i) {
parameters[i].setEstimate(parameters[i].getEstimate() + dX.getEntry(i));
}
} catch(InvalidMatrixException e) {
throw new EstimationException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM);
}
previous = cost;
updateResidualsAndCost();
} while ((getCostEvaluations() < 2) ||
(FastMath.abs(previous - cost) > (cost * steadyStateThreshold) &&
(FastMath.abs(cost) > convergence)));
}
}

897
src/main/java/org/apache/commons/math/estimation/LevenbergMarquardtEstimator.java
View File

@ -1,897 +0,0 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.estimation;
import java.io.Serializable;
import java.util.Arrays;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.util.FastMath;
/**
* This class solves a least squares problem.
*
* <p>This implementation <em>should</em> work even for over-determined systems
* (i.e. systems having more variables than equations). Over-determined systems
* are solved by ignoring the variables which have the smallest impact according
* to their jacobian column norm. Only the rank of the matrix and some loop bounds
* are changed to implement this.</p>
*
* <p>The resolution engine is a simple translation of the MINPACK <a
* href="http://www.netlib.org/minpack/lmder.f">lmder</a> routine with minor
* changes. The changes include the over-determined resolution and the Q.R.
* decomposition which has been rewritten following the algorithm described in the
* P. Lascaux and R. Theodor book <i>Analyse num&eacute;rique matricielle
* appliqu&eacute;e &agrave; l'art de l'ing&eacute;nieur</i>, Masson 1986.</p>
* <p>The authors of the original fortran version are:
* <ul>
* <li>Argonne National Laboratory. MINPACK project. March 1980</li>
* <li>Burton S. Garbow</li>
* <li>Kenneth E. Hillstrom</li>
* <li>Jorge J. More</li>
* </ul>
* The redistribution policy for MINPACK is available <a
* href="http://www.netlib.org/minpack/disclaimer">here</a>, for convenience, it
* is reproduced below.</p>
*
* <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
* <tr><td>
* Minpack Copyright Notice (1999) University of Chicago.
* All rights reserved
* </td></tr>
* <tr><td>
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* <ol>
* <li>Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.</li>
* <li>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.</li>
* <li>The end-user documentation included with the redistribution, if any,
* must include the following acknowledgment:
* <code>This product includes software developed by the University of
* Chicago, as Operator of Argonne National Laboratory.</code>
* Alternately, this acknowledgment may appear in the software itself,
* if and wherever such third-party acknowledgments normally appear.</li>
* <li><strong>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.</strong></li>
* <li><strong>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.</strong></li>
* <ol></td></tr>
* </table>
* @version $Revision$ $Date$
* @since 1.2
* @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
* been deprecated and replaced by package org.apache.commons.math.optimization.general
*
*/
@Deprecated
public class LevenbergMarquardtEstimator extends AbstractEstimator implements Serializable {
/** Serializable version identifier */
private static final long serialVersionUID = -5705952631533171019L;
/** Number of solved variables. */
private int solvedCols;
/** Diagonal elements of the R matrix in the Q.R. decomposition. */
private double[] diagR;
/** Norms of the columns of the jacobian matrix. */
private double[] jacNorm;
/** Coefficients of the Householder transforms vectors. */
private double[] beta;
/** Columns permutation array. */
private int[] permutation;
/** Rank of the jacobian matrix. */
private int rank;
/** Levenberg-Marquardt parameter. */
private double lmPar;
/** Parameters evolution direction associated with lmPar. */
private double[] lmDir;
/** Positive input variable used in determining the initial step bound. */
private double initialStepBoundFactor;
/** Desired relative error in the sum of squares. */
private double costRelativeTolerance;
/** Desired relative error in the approximate solution parameters. */
private double parRelativeTolerance;
/** Desired max cosine on the orthogonality between the function vector
* and the columns of the jacobian. */
private double orthoTolerance;
/**
* Build an estimator for least squares problems.
* <p>The default values for the algorithm settings are:
* <ul>
* <li>{@link #setInitialStepBoundFactor initial step bound factor}: 100.0</li>
* <li>{@link #setMaxCostEval maximal cost evaluations}: 1000</li>
* <li>{@link #setCostRelativeTolerance cost relative tolerance}: 1.0e-10</li>
* <li>{@link #setParRelativeTolerance parameters relative tolerance}: 1.0e-10</li>
* <li>{@link #setOrthoTolerance orthogonality tolerance}: 1.0e-10</li>
* </ul>
* </p>
*/
public LevenbergMarquardtEstimator() {
// set up the superclass with a default max cost evaluations setting
setMaxCostEval(1000);
// default values for the tuning parameters
setInitialStepBoundFactor(100.0);
setCostRelativeTolerance(1.0e-10);
setParRelativeTolerance(1.0e-10);
setOrthoTolerance(1.0e-10);
}
/**
* Set the positive input variable used in determining the initial step bound.
* This bound is set to the product of initialStepBoundFactor and the euclidean norm of diag*x if nonzero,
* or else to initialStepBoundFactor itself. In most cases factor should lie
* in the interval (0.1, 100.0). 100.0 is a generally recommended value
*
* @param initialStepBoundFactor initial step bound factor
* @see #estimate
*/
public void setInitialStepBoundFactor(double initialStepBoundFactor) {
this.initialStepBoundFactor = initialStepBoundFactor;
}
/**
* Set the desired relative error in the sum of squares.
*
* @param costRelativeTolerance desired relative error in the sum of squares
* @see #estimate
*/
public void setCostRelativeTolerance(double costRelativeTolerance) {
this.costRelativeTolerance = costRelativeTolerance;
}
/**
* Set the desired relative error in the approximate solution parameters.
*
* @param parRelativeTolerance desired relative error
* in the approximate solution parameters
* @see #estimate
*/
public void setParRelativeTolerance(double parRelativeTolerance) {
this.parRelativeTolerance = parRelativeTolerance;
}
/**
* Set the desired max cosine on the orthogonality.
*
* @param orthoTolerance desired max cosine on the orthogonality
* between the function vector and the columns of the jacobian
* @see #estimate
*/
public void setOrthoTolerance(double orthoTolerance) {
this.orthoTolerance = orthoTolerance;
}
/**
* Solve an estimation problem using the Levenberg-Marquardt algorithm.
* <p>The algorithm used is a modified Levenberg-Marquardt one, based
* on the MINPACK <a href="http://www.netlib.org/minpack/lmder.f">lmder</a>
* routine. The algorithm settings must have been set up before this method
* is called with the {@link #setInitialStepBoundFactor},
* {@link #setMaxCostEval}, {@link #setCostRelativeTolerance},
* {@link #setParRelativeTolerance} and {@link #setOrthoTolerance} methods.
* If these methods have not been called, the default values set up by the
* {@link #LevenbergMarquardtEstimator() constructor} will be used.</p>
* <p>The authors of the original fortran function are:</p>
* <ul>
* <li>Argonne National Laboratory. MINPACK project. March 1980</li>
* <li>Burton S. Garbow</li>
* <li>Kenneth E. Hillstrom</li>
* <li>Jorge J. More</li>
* </ul>
* <p>Luc Maisonobe did the Java translation.</p>
*
* @param problem estimation problem to solve
* @exception EstimationException if convergence cannot be
* reached with the specified algorithm settings or if there are more variables
* than equations
* @see #setInitialStepBoundFactor
* @see #setCostRelativeTolerance
* @see #setParRelativeTolerance
* @see #setOrthoTolerance
*/
@Override
public void estimate(EstimationProblem problem)
throws EstimationException {
initializeEstimate(problem);
// arrays shared with the other private methods
solvedCols = FastMath.min(rows, cols);
diagR = new double[cols];
jacNorm = new double[cols];
beta = new double[cols];
permutation = new int[cols];
lmDir = new double[cols];
// local variables
double delta = 0;
double xNorm = 0;
double[] diag = new double[cols];
double[] oldX = new double[cols];
double[] oldRes = new double[rows];
double[] work1 = new double[cols];
double[] work2 = new double[cols];
double[] work3 = new double[cols];
// evaluate the function at the starting point and calculate its norm
updateResidualsAndCost();
// outer loop
lmPar = 0;
boolean firstIteration = true;
while (true) {
// compute the Q.R. decomposition of the jacobian matrix
updateJacobian();
qrDecomposition();
// compute Qt.res
qTy(residuals);
// now we don't need Q anymore,
// so let jacobian contain the R matrix with its diagonal elements
for (int k = 0; k < solvedCols; ++k) {
int pk = permutation[k];
jacobian[k * cols + pk] = diagR[pk];
}
if (firstIteration) {
// scale the variables according to the norms of the columns
// of the initial jacobian
xNorm = 0;
for (int k = 0; k < cols; ++k) {
double dk = jacNorm[k];
if (dk == 0) {
dk = 1.0;
}
double xk = dk * parameters[k].getEstimate();
xNorm += xk * xk;
diag[k] = dk;
}
xNorm = 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 (cost != 0) {
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
double s = jacNorm[pj];
if (s != 0) {
double sum = 0;
int index = pj;
for (int i = 0; i <= j; ++i) {
sum += jacobian[index] * residuals[i];
index += cols;
}
maxCosine = FastMath.max(maxCosine, FastMath.abs(sum) / (s * cost));
}
}
}
if (maxCosine <= orthoTolerance) {
return;
}
// rescale if necessary
for (int j = 0; j < cols; ++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] = parameters[pj].getEstimate();
}
double previousCost = cost;
double[] tmpVec = residuals;
residuals = oldRes;
oldRes = tmpVec;
// determine the Levenberg-Marquardt parameter
determineLMParameter(oldRes, delta, diag, work1, work2, work3);
// compute the new point and the norm of the evolution direction
double lmNorm = 0;
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
lmDir[pj] = -lmDir[pj];
parameters[pj].setEstimate(oldX[pj] + lmDir[pj]);
double s = diag[pj] * lmDir[pj];
lmNorm += s * s;
}
lmNorm = 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
updateResidualsAndCost();
// compute the scaled actual reduction
double actRed = -1.0;
if (0.1 * cost < previousCost) {
double r = cost / previousCost;
actRed = 1.0 - r * r;
}
// compute the scaled predicted reduction
// and the scaled directional derivative
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
double dirJ = lmDir[pj];
work1[j] = 0;
int index = pj;
for (int i = 0; i <= j; ++i) {
work1[i] += jacobian[index] * dirJ;
index += cols;
}
}
double coeff1 = 0;
for (int j = 0; j < solvedCols; ++j) {
coeff1 += work1[j] * work1[j];
}
double pc2 = previousCost * previousCost;
coeff1 = coeff1 / pc2;
double coeff2 = lmPar * lmNorm * lmNorm / pc2;
double preRed = coeff1 + 2 * coeff2;
double dirDer = -(coeff1 + coeff2);
// ratio of the actual to the predicted reduction
ratio = (preRed == 0) ? 0 : (actRed / preRed);
// update the step bound
if (ratio <= 0.25) {
double tmp =
(actRed < 0) ? (0.5 * dirDer / (dirDer + 0.5 * actRed)) : 0.5;
if ((0.1 * cost >= previousCost) || (tmp < 0.1)) {
tmp = 0.1;
}
delta = tmp * 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 < cols; ++k) {
double xK = diag[k] * parameters[k].getEstimate();
xNorm += xK * xK;
}
xNorm = FastMath.sqrt(xNorm);
} else {
// failed iteration, reset the previous values
cost = previousCost;
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
parameters[pj].setEstimate(oldX[pj]);
}
tmpVec = residuals;
residuals = oldRes;
oldRes = tmpVec;
}
// tests for convergence.
if (((FastMath.abs(actRed) <= costRelativeTolerance) &&
(preRed <= costRelativeTolerance) &&
(ratio <= 2.0)) ||
(delta <= parRelativeTolerance * xNorm)) {
return;
}
// tests for termination and stringent tolerances
// (2.2204e-16 is the machine epsilon for IEEE754)
if ((FastMath.abs(actRed) <= 2.2204e-16) && (preRed <= 2.2204e-16) && (ratio <= 2.0)) {
throw new EstimationException("cost relative tolerance is too small ({0})," +
" no further reduction in the" +
" sum of squares is possible",
costRelativeTolerance);
} else if (delta <= 2.2204e-16 * xNorm) {
throw new EstimationException("parameters relative tolerance is too small" +
" ({0}), no further improvement in" +
" the approximate solution is possible",
parRelativeTolerance);
} else if (maxCosine <= 2.2204e-16) {
throw new EstimationException("orthogonality tolerance is too small ({0})," +
" solution is orthogonal to the jacobian",
orthoTolerance);
}
}
}
}
/**
* Determine the Levenberg-Marquardt parameter.
* <p>This implementation is a translation in Java of the MINPACK
* <a href="http://www.netlib.org/minpack/lmpar.f">lmpar</a>
* routine.</p>
* <p>This method sets the lmPar and lmDir attributes.</p>
* <p>The authors of the original fortran function are:</p>
* <ul>
* <li>Argonne National Laboratory. MINPACK project. March 1980</li>
* <li>Burton S. Garbow</li>
* <li>Kenneth E. Hillstrom</li>
* <li>Jorge J. More</li>
* </ul>
* <p>Luc Maisonobe did the Java translation.</p>
*
* @param qy array containing qTy
* @param delta upper bound on the euclidean norm of diagR * lmDir
* @param diag diagonal matrix
* @param work1 work array
* @param work2 work array
* @param work3 work array
*/
private void determineLMParameter(double[] qy, double delta, double[] diag,
double[] work1, double[] work2, double[] work3) {
// compute and store in x the gauss-newton direction, if the
// jacobian is rank-deficient, obtain a least squares solution
for (int j = 0; j < rank; ++j) {
lmDir[permutation[j]] = qy[j];
}
for (int j = rank; j < cols; ++j) {
lmDir[permutation[j]] = 0;
}
for (int k = rank - 1; k >= 0; --k) {
int pk = permutation[k];
double ypk = lmDir[pk] / diagR[pk];
int index = pk;
for (int i = 0; i < k; ++i) {
lmDir[permutation[i]] -= ypk * jacobian[index];
index += cols;
}
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;
int index = pj;
for (int i = 0; i < j; ++i) {
sum += jacobian[index] * work1[permutation[i]];
index += cols;
}
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;
int index = pj;
for (int i = 0; i <= j; ++i) {
sum += jacobian[index] * qy[i];
index += cols;
}
sum /= diag[pj];
sum2 += sum * sum;
}
double gNorm = FastMath.sqrt(sum2);
double paru = gNorm / delta;
if (paru == 0) {
// 2.2251e-308 is the smallest positive real for IEE754
paru = 2.2251e-308 / 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(2.2251e-308, 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, 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]] -= jacobian[i * cols + pj] * tmp;
}
}
sum2 = 0;
for (int j = 0; j < solvedCols; ++j) {
double s = work1[permutation[j]];
sum2 += s * s;
}
double correction = fp / (delta * sum2);
// depending on the sign of the function, update parl or paru.
if (fp > 0) {
parl = 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);
}
}
/**
* Solve a*x = b and d*x = 0 in the least squares sense.
* <p>This implementation is a translation in Java of the MINPACK
* <a href="http://www.netlib.org/minpack/qrsolv.f">qrsolv</a>
* routine.</p>
* <p>This method sets the lmDir and lmDiag attributes.</p>
* <p>The authors of the original fortran function are:</p>
* <ul>
* <li>Argonne National Laboratory. MINPACK project. March 1980</li>
* <li>Burton S. Garbow</li>
* <li>Kenneth E. Hillstrom</li>
* <li>Jorge J. More</li>
* </ul>
* <p>Luc Maisonobe did the Java translation.</p>
*
* @param qy array containing qTy
* @param diag diagonal matrix
* @param lmDiag diagonal elements associated with lmDir
* @param work work array
*/
private void determineLMDirection(double[] qy, double[] diag,
double[] lmDiag, double[] work) {
// copy R and Qty to preserve input and initialize s
// in particular, save the diagonal elements of R in lmDir
for (int j = 0; j < solvedCols; ++j) {
int pj = permutation[j];
for (int i = j + 1; i < solvedCols; ++i) {
jacobian[i * cols + pj] = jacobian[j * cols + permutation[i]];
}
lmDir[j] = diagR[pj];
work[j] = qy[j];
}
// eliminate the diagonal matrix d using a Givens rotation
for (int j = 0; j < solvedCols; ++j) {
// prepare the row of d to be eliminated, locating the
// diagonal element using p from the Q.R. factorization
int pj = permutation[j];
double dpj = diag[pj];
if (dpj != 0) {
Arrays.fill(lmDiag, j + 1, lmDiag.length, 0);
}
lmDiag[j] = dpj;
// the transformations to eliminate the row of d
// modify only a single element of Qty
// beyond the first n, which is initially zero.
double qtbpj = 0;
for (int k = j; k < solvedCols; ++k) {
int pk = permutation[k];
// determine a Givens rotation which eliminates the
// appropriate element in the current row of d
if (lmDiag[k] != 0) {
final double sin;
final double cos;
double rkk = jacobian[k * cols + 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)
jacobian[k * cols + 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 = jacobian[i * cols + pk];
final double temp2 = cos * rik + sin * lmDiag[i];
lmDiag[i] = -sin * rik + cos * lmDiag[i];
jacobian[i * cols + pk] = temp2;
}
}
}
// store the diagonal element of s and restore
// the corresponding diagonal element of R
int index = j * cols + permutation[j];
lmDiag[j] = jacobian[index];
jacobian[index] = lmDir[j];
}
// solve the triangular system for z, if the system is
// singular, then obtain a least squares solution
int nSing = solvedCols;
for (int j = 0; j < solvedCols; ++j) {
if ((lmDiag[j] == 0) && (nSing == solvedCols)) {
nSing = j;
}
if (nSing < solvedCols) {
work[j] = 0;
}
}
if (nSing > 0) {
for (int j = nSing - 1; j >= 0; --j) {
int pj = permutation[j];
double sum = 0;
for (int i = j + 1; i < nSing; ++i) {
sum += jacobian[i * cols + pj] * work[i];
}
work[j] = (work[j] - sum) / lmDiag[j];
}
}
// permute the components of z back to components of lmDir
for (int j = 0; j < lmDir.length; ++j) {
lmDir[permutation[j]] = work[j];
}
}
/**
* Decompose a matrix A as A.P = Q.R using Householder transforms.
* <p>As suggested in the P. Lascaux and R. Theodor book
* <i>Analyse num&eacute;rique matricielle appliqu&eacute;e &agrave;
* l'art de l'ing&eacute;nieur</i> (Masson, 1986), instead of representing
* the Householder transforms with u<sub>k</sub> unit vectors such that:
* <pre>
* H<sub>k</sub> = I - 2u<sub>k</sub>.u<sub>k</sub><sup>t</sup>
* </pre>
* we use <sub>k</sub> non-unit vectors such that:
* <pre>
* H<sub>k</sub> = I - beta<sub>k</sub>v<sub>k</sub>.v<sub>k</sub><sup>t</sup>
* </pre>
* where v<sub>k</sub> = a<sub>k</sub> - alpha<sub>k</sub> e<sub>k</sub>.
* The beta<sub>k</sub> coefficients are provided upon exit as recomputing
* them from the v<sub>k</sub> vectors would be costly.</p>
* <p>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.</p>
* @exception EstimationException if the decomposition cannot be performed
*/
private void qrDecomposition() throws EstimationException {
// initializations
for (int k = 0; k < cols; ++k) {
permutation[k] = k;
double norm2 = 0;
for (int index = k; index < jacobian.length; index += cols) {
double akk = jacobian[index];
norm2 += akk * akk;
}
jacNorm[k] = FastMath.sqrt(norm2);
}
// transform the matrix column after column
for (int k = 0; k < cols; ++k) {
// select the column with the greatest norm on active components
int nextColumn = -1;
double ak2 = Double.NEGATIVE_INFINITY;
for (int i = k; i < cols; ++i) {
double norm2 = 0;
int iDiag = k * cols + permutation[i];
for (int index = iDiag; index < jacobian.length; index += cols) {
double aki = jacobian[index];
norm2 += aki * aki;
}
if (Double.isInfinite(norm2) || Double.isNaN(norm2)) {
throw new EstimationException(
LocalizedFormats.UNABLE_TO_PERFORM_QR_DECOMPOSITION_ON_JACOBIAN,
rows, cols);
}
if (norm2 > ak2) {
nextColumn = i;
ak2 = norm2;
}
}
if (ak2 == 0) {
rank = k;
return;
}
int pk = permutation[nextColumn];
permutation[nextColumn] = permutation[k];
permutation[k] = pk;
// choose alpha such that Hk.u = alpha ek
int kDiag = k * cols + pk;
double akk = jacobian[kDiag];
double alpha = (akk > 0) ? -FastMath.sqrt(ak2) : FastMath.sqrt(ak2);
double betak = 1.0 / (ak2 - akk * alpha);
beta[pk] = betak;
// transform the current column
diagR[pk] = alpha;
jacobian[kDiag] -= alpha;
// transform the remaining columns
for (int dk = cols - 1 - k; dk > 0; --dk) {
int dkp = permutation[k + dk] - pk;
double gamma = 0;
for (int index = kDiag; index < jacobian.length; index += cols) {
gamma += jacobian[index] * jacobian[index + dkp];
}
gamma *= betak;
for (int index = kDiag; index < jacobian.length; index += cols) {
jacobian[index + dkp] -= gamma * jacobian[index];
}
}
}
rank = solvedCols;
}
/**
* Compute the product Qt.y for some Q.R. decomposition.
*
* @param y vector to multiply (will be overwritten with the result)
*/
private void qTy(double[] y) {
for (int k = 0; k < cols; ++k) {
int pk = permutation[k];
int kDiag = k * cols + pk;
double gamma = 0;
int index = kDiag;
for (int i = k; i < rows; ++i) {
gamma += jacobian[index] * y[i];
index += cols;
}
gamma *= beta[pk];
index = kDiag;
for (int i = k; i < rows; ++i) {
y[i] -= gamma * jacobian[index];
index += cols;
}
}
}
}

111
src/main/java/org/apache/commons/math/estimation/SimpleEstimationProblem.java
View File

@ -1,111 +0,0 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.estimation;
import java.util.ArrayList;
import java.util.List;
/**
* Simple implementation of the {@link EstimationProblem
* EstimationProblem} interface for boilerplate data handling.
* <p>This class <em>only</em> handles parameters and measurements
* storage and unbound parameters filtering. It does not compute
* anything by itself. It should either be used with measurements
* implementation that are smart enough to know about the
* various parameters in order to compute the partial derivatives
* appropriately. Since the problem-specific logic is mainly related to
* the various measurements models, the simplest way to use this class
* is by extending it and using one internal class extending
* {@link WeightedMeasurement WeightedMeasurement} for each measurement
* type. The instances of the internal classes would have access to the
* various parameters and their current estimate.</p>
* @version $Revision$ $Date$
* @since 1.2
* @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
* been deprecated and replaced by package org.apache.commons.math.optimization.general
*/
@Deprecated
public class SimpleEstimationProblem implements EstimationProblem {
/** Estimated parameters. */
private final List<EstimatedParameter> parameters;
/** Measurements. */
private final List<WeightedMeasurement> measurements;
/**
* Build an empty instance without parameters nor measurements.
*/
public SimpleEstimationProblem() {
parameters = new ArrayList<EstimatedParameter>();
measurements = new ArrayList<WeightedMeasurement>();
}
/**
* Get all the parameters of the problem.
* @return parameters
*/
public EstimatedParameter[] getAllParameters() {
return parameters.toArray(new EstimatedParameter[parameters.size()]);
}
/**
* Get the unbound parameters of the problem.
* @return unbound parameters
*/
public EstimatedParameter[] getUnboundParameters() {
// filter the unbound parameters
List<EstimatedParameter> unbound = new ArrayList<EstimatedParameter>(parameters.size());
for (EstimatedParameter p : parameters) {
if (! p.isBound()) {
unbound.add(p);
}
}
// convert to an array
return unbound.toArray(new EstimatedParameter[unbound.size()]);
}
/**
* Get the measurements of an estimation problem.
* @return measurements
*/
public WeightedMeasurement[] getMeasurements() {
return measurements.toArray(new WeightedMeasurement[measurements.size()]);
}
/** Add a parameter to the problem.
* @param p parameter to add
*/
protected void addParameter(EstimatedParameter p) {
parameters.add(p);
}
/**
* Add a new measurement to the set.
* @param m measurement to add
*/
protected void addMeasurement(WeightedMeasurement m) {
measurements.add(m);
}
}

172
src/main/java/org/apache/commons/math/estimation/WeightedMeasurement.java
View File

@ -1,172 +0,0 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.estimation;
import java.io.Serializable;
/**
* This class represents measurements in estimation problems.
*
* <p>This abstract class implements all the methods needed to handle
* measurements in a general way. It defines neither the {@link
* #getTheoreticalValue getTheoreticalValue} nor the {@link
* #getPartial getPartial} methods, which should be defined by
* sub-classes according to the specific problem.</p>
*
* <p>The {@link #getTheoreticalValue getTheoreticalValue} and {@link
* #getPartial getPartial} methods must always use the current
* estimate of the parameters set by the solver in the problem. These
* parameters can be retrieved through the {@link
* EstimationProblem#getAllParameters
* EstimationProblem.getAllParameters} method if the measurements are
* independent of the problem, or directly if they are implemented as
* inner classes of the problem.</p>
*
* <p>The instances for which the <code>ignored</code> flag is set
* through the {@link #setIgnored setIgnored} method are ignored by the
* solvers. This can be used to reject wrong measurements at some
* steps of the estimation.</p>
*
* @see EstimationProblem
*
* @version $Revision$ $Date$
* @since 1.2
* @deprecated as of 2.0, everything in package org.apache.commons.math.estimation has
* been deprecated and replaced by package org.apache.commons.math.optimization.general
*/
@Deprecated
public abstract class WeightedMeasurement implements Serializable {
/** Serializable version identifier. */
private static final long serialVersionUID = 4360046376796901941L;
/** Measurement weight. */
private final double weight;
/** Value of the measurements. */
private final double measuredValue;
/** Ignore measurement indicator. */
private boolean ignored;
/**
* Simple constructor.
* Build a measurement with the given parameters, and set its ignore
* flag to false.
* @param weight weight of the measurement in the least squares problem
* (two common choices are either to use 1.0 for all measurements, or to
* use a value proportional to the inverse of the variance of the measurement
* type)
*
* @param measuredValue measured value
*/
public WeightedMeasurement(double weight, double measuredValue) {
this.weight = weight;
this.measuredValue = measuredValue;
ignored = false;
}
/** Simple constructor.
*
* Build a measurement with the given parameters
*
* @param weight weight of the measurement in the least squares problem
* @param measuredValue measured value
* @param ignored true if the measurement should be ignored
*/
public WeightedMeasurement(double weight, double measuredValue,
boolean ignored) {
this.weight = weight;
this.measuredValue = measuredValue;
this.ignored = ignored;
}
/**
* Get the weight of the measurement in the least squares problem
*
* @return weight
*/
public double getWeight() {
return weight;
}
/**
* Get the measured value
*
* @return measured value
*/
public double getMeasuredValue() {
return measuredValue;
}
/**
* Get the residual for this measurement
* The residual is the measured value minus the theoretical value.
*
* @return residual
*/
public double getResidual() {
return measuredValue - getTheoreticalValue();
}
/**
* Get the theoretical value expected for this measurement
* <p>The theoretical value is the value expected for this measurement
* if the model and its parameter were all perfectly known.</p>
* <p>The value must be computed using the current estimate of the parameters
* set by the solver in the problem.</p>
*
* @return theoretical value
*/
public abstract double getTheoreticalValue();
/**
* Get the partial derivative of the {@link #getTheoreticalValue
* theoretical value} according to the parameter.
* <p>The value must be computed using the current estimate of the parameters
* set by the solver in the problem.</p>
*
* @param parameter parameter against which the partial derivative
* should be computed
* @return partial derivative of the {@link #getTheoreticalValue
* theoretical value}
*/
public abstract double getPartial(EstimatedParameter parameter);
/**
* Set the ignore flag to the specified value
* Setting the ignore flag to true allow to reject wrong
* measurements, which sometimes can be detected only rather late.
*
* @param ignored value for the ignore flag
*/
public void setIgnored(boolean ignored) {
this.ignored = ignored;
}
/**
* Check if this measurement should be ignored
*
* @return true if the measurement should be ignored
*/
public boolean isIgnored() {
return ignored;
}
}

25
src/main/java/org/apache/commons/math/estimation/package.html
View File

@ -1,25 +0,0 @@
<html>
<!--
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.
-->
<!-- $Revision$ -->
<body>
This package provided classes to solve estimation problems, it is deprecated since 2.0.
<p>This package has been deprecated as of 2.0. It is replaced by the optimization.general package.</p>
</body>
</html>

24
src/main/java/org/apache/commons/math/exception/DimensionMismatchException.java
View File

@ -17,6 +17,7 @@
package org.apache.commons.math.exception;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.exception.util.Localizable;
/**
* Exception to be thrown when two dimensions differ.
@ -25,13 +26,27 @@ import org.apache.commons.math.exception.util.LocalizedFormats;
* @version $Revision$ $Date$
*/
public class DimensionMismatchException extends MathIllegalNumberException {
/** Serializable version Id. */
private static final long serialVersionUID = -8415396756375798143L;
/** Correct dimension. */
private final int dimension;
/**
* Construct an exception from the mismatched dimensions.
*
* @param specific Specific context information pattern.
* @param wrong Wrong dimension.
* @param expected Expected dimension.
*/
public DimensionMismatchException(Localizable specific,
int wrong,
int expected) {
super(specific,
LocalizedFormats.DIMENSIONS_MISMATCH,
wrong, expected);
dimension = expected;
}
/**
* Construct an exception from the mismatched dimensions.
*
@ -40,10 +55,7 @@ public class DimensionMismatchException extends MathIllegalNumberException {
*/
public DimensionMismatchException(int wrong,
int expected) {
super(LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
LocalizedFormats.DIMENSIONS_MISMATCH,
wrong, expected);
dimension = expected;
this(LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, wrong, expected);
}
/**

4
src/main/java/org/apache/commons/math/exception/MathIllegalNumberException.java
View File

@ -40,7 +40,7 @@ public class MathIllegalNumberException extends MathIllegalArgumentException {
*
* @param specific Localizable pattern.
* @param general Localizable pattern.
* @param wrong wrong number
* @param wrong Wrong number.
* @param arguments Arguments.
*/
protected MathIllegalNumberException(Localizable specific,
@ -55,7 +55,7 @@ public class MathIllegalNumberException extends MathIllegalArgumentException {
* Construct an exception.
*
* @param general Localizable pattern.
* @param wrong wrong number
* @param wrong Wrong number.
* @param arguments Arguments.
*/
protected MathIllegalNumberException(Localizable general,

6
src/main/java/org/apache/commons/math/exception/MaxCountExceededException.java
View File

@ -28,7 +28,6 @@ import org.apache.commons.math.exception.util.LocalizedFormats;
public class MaxCountExceededException extends MathIllegalStateException {
/** Serializable version Id. */
private static final long serialVersionUID = 4330003017885151975L;
/**
* Maximum number of evaluations.
*/
@ -50,10 +49,7 @@ public class MaxCountExceededException extends MathIllegalStateException {
*/
public MaxCountExceededException(Localizable specific,
Number max) {
super(specific,
LocalizedFormats.MAX_COUNT_EXCEEDED,
max);
super(specific, LocalizedFormats.MAX_COUNT_EXCEEDED, max);
this.max = max;
}

64
src/main/java/org/apache/commons/math/exception/NonPositiveDefiniteMatrixException.java Normal file
View File

@ -0,0 +1,64 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.exception;
import org.apache.commons.math.exception.util.LocalizedFormats;
/**
* Exception to be thrown when a symmetric matrix is expected.
*
* @since 3.0
* @version $Revision$ $Date$
*/
public class NonPositiveDefiniteMatrixException extends MathIllegalArgumentException {
/** Index (diagonal element). */
private final int index;
/** Threshold. */
private final double threshold;
/**
* Construct an exception.
*
* @param index Row (and column) index.
* @param threshold Absolute positivity threshold.
*/
public NonPositiveDefiniteMatrixException(int index,
double threshold) {
super(LocalizedFormats.NON_POSITIVE_DEFINITE_MATRIX, index, threshold);
this.index = index;
this.threshold = threshold;
}
/**
* @return the row index.
*/
public int getRow() {
return index;
}
/**
* @return the column index.
*/
public int getColumn() {
return index;
}
/**
* @return the absolute positivity threshold.
*/
public double getThreshold() {
return threshold;
}
}

38
src/main/java/org/apache/commons/math/exception/NonSquareMatrixException.java Normal file
View File

@ -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.math.exception;
import org.apache.commons.math.exception.util.LocalizedFormats;
/**
* Exception to be thrown when a square matrix is expected.
*
* @since 3.0
* @version $Revision$ $Date$
*/
public class NonSquareMatrixException extends DimensionMismatchException {
/**
* Construct an exception from the mismatched dimensions.
*
* @param wrong Row dimension.
* @param expected Column dimension.
*/
public NonSquareMatrixException(int wrong,
int expected) {
super(LocalizedFormats.NON_SQUARE_MATRIX, wrong, expected);
}
}

69
src/main/java/org/apache/commons/math/exception/NonSymmetricMatrixException.java Normal file
View File

@ -0,0 +1,69 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.exception;
import org.apache.commons.math.exception.util.LocalizedFormats;
/**
* Exception to be thrown when a symmetric matrix is expected.
*
* @since 3.0
* @version $Revision$ $Date$
*/
public class NonSymmetricMatrixException extends MathIllegalArgumentException {
/** Row. */
private final int row;
/** Column. */
private final int column;
/** Threshold. */
private final double threshold;
/**
* Construct an exception.
*
* @param row Row index.
* @param column Column index.
* @param threshold Relative symmetry threshold.
*/
public NonSymmetricMatrixException(int row,
int column,
double threshold) {
super(LocalizedFormats.NON_SYMMETRIC_MATRIX, row, column, threshold);
this.row = row;
this.column = column;
this.threshold = threshold;
}
/**
* @return the row index of the entry.
*/
public int getRow() {
return row;
}
/**
* @return the column index of the entry.
*/
public int getColumn() {
return column;
}
/**
* @return the relative symmetry threshold.
*/
public double getThreshold() {
return threshold;
}
}

34
src/main/java/org/apache/commons/math/exception/SingularMatrixException.java Normal file
View File

@ -0,0 +1,34 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.exception;
import org.apache.commons.math.exception.util.LocalizedFormats;
/**
* Exception to be thrown when a non-singular matrix is expected.
*
* @since 3.0
* @version $Revision$ $Date$
*/
public class SingularMatrixException extends MathIllegalArgumentException {
/**
* Construct an exception.
*/
public SingularMatrixException() {
super(LocalizedFormats.SINGULAR_MATRIX);
}
}

4
src/main/java/org/apache/commons/math/exception/util/LocalizedFormats.java
View File

@ -181,6 +181,7 @@ public enum LocalizedFormats implements Localizable {
BETA("beta"), /* keep */
NOT_POSITIVE_COLUMNDIMENSION("invalid column dimension: {0} (must be positive)"),
NOT_POSITIVE_DEFINITE_MATRIX("not positive definite matrix"),
NON_POSITIVE_DEFINITE_MATRIX("not positive definite matrix: diagonal element at ({0},{0}) is larger than {2}"), /* keep */
DEGREES_OF_FREEDOM("degrees of freedom ({0})"), /* keep */
NOT_POSITIVE_DEGREES_OF_FREEDOM("degrees of freedom must be positive ({0})"),
NOT_POSITIVE_ELEMENT_AT_INDEX("element {0} is not positive: {1}"),
@ -217,6 +218,7 @@ public enum LocalizedFormats implements Localizable {
NOT_STRICTLY_INCREASING_SEQUENCE("points {3} and {2} are not strictly increasing ({1} >= {0})"), /* keep */
NOT_SUBTRACTION_COMPATIBLE_MATRICES("{0}x{1} and {2}x{3} matrices are not subtraction compatible"),
NOT_SYMMETRIC_MATRIX("not symmetric matrix"),
NON_SYMMETRIC_MATRIX("non symmetric matrix: the difference between entries at ({0},{1}) and ({1},{0}) is larger than {2}"), /* keep */
NO_BIN_SELECTED("no bin selected"),
NO_CONVERGENCE_WITH_ANY_START_POINT("none of the {0} start points lead to convergence"), /* keep */
NO_DATA("no data"), /* keep */
@ -276,7 +278,7 @@ public enum LocalizedFormats implements Localizable {
SAMPLE_SIZE_LARGER_THAN_POPULATION_SIZE("sample size ({0}) must be less than or equal to population size ({1})"),
SIMPLEX_NEED_ONE_POINT("simplex must contain at least one point"),
SIMPLE_MESSAGE("{0}"),
SINGULAR_MATRIX("matrix is singular"),
SINGULAR_MATRIX("matrix is singular"), /* keep */
SUBARRAY_ENDS_AFTER_ARRAY_END("subarray ends after array end"),
TOO_LARGE_CUTOFF_SINGULAR_VALUE("cutoff singular value is {0}, should be at most {1}"),
TOO_MANY_ELEMENTS_TO_DISCARD_FROM_ARRAY("cannot discard {0} elements from a {1} elements array"),

4
src/main/java/org/apache/commons/math/linear/AbstractFieldMatrix.java
View File

@ -23,6 +23,7 @@ import java.util.Arrays;
import org.apache.commons.math.Field;
import org.apache.commons.math.FieldElement;
import org.apache.commons.math.exception.MatrixDimensionMismatchException;
import org.apache.commons.math.exception.NonSquareMatrixException;
import org.apache.commons.math.exception.DimensionMismatchException;
import org.apache.commons.math.exception.NoDataException;
import org.apache.commons.math.exception.OutOfRangeException;
@ -620,8 +621,7 @@ public abstract class AbstractFieldMatrix<T extends FieldElement<T>> implements
public abstract int getColumnDimension();
/** {@inheritDoc} */
public T getTrace()
throws NonSquareMatrixException {
public T getTrace() {
final int nRows = getRowDimension();
final int nCols = getColumnDimension();
if (nRows != nCols) {

4
src/main/java/org/apache/commons/math/linear/AbstractRealMatrix.java
View File

@ -20,6 +20,7 @@ package org.apache.commons.math.linear;
import org.apache.commons.math.exception.NoDataException;
import org.apache.commons.math.exception.NotStrictlyPositiveException;
import org.apache.commons.math.exception.DimensionMismatchException;
import org.apache.commons.math.exception.NonSquareMatrixException;
import org.apache.commons.math.exception.MatrixDimensionMismatchException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.util.MathUtils;
@ -576,8 +577,7 @@ public abstract class AbstractRealMatrix implements RealMatrix {
public abstract int getColumnDimension();
/** {@inheritDoc} */
public double getTrace()
throws NonSquareMatrixException {
public double getTrace() {
final int nRows = getRowDimension();
final int nCols = getColumnDimension();
if (nRows != nCols) {

93
src/main/java/org/apache/commons/math/linear/CholeskyDecompositionImpl.java
View File

@ -17,8 +17,10 @@
package org.apache.commons.math.linear;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.exception.DimensionMismatchException;
import org.apache.commons.math.exception.NonSquareMatrixException;
import org.apache.commons.math.exception.NonSymmetricMatrixException;
import org.apache.commons.math.exception.NonPositiveDefiniteMatrixException;
import org.apache.commons.math.util.FastMath;
@ -34,21 +36,20 @@ import org.apache.commons.math.util.FastMath;
* @since 2.0
*/
public class CholeskyDecompositionImpl implements CholeskyDecomposition {
/** Default threshold above which off-diagonal elements are considered too different
* and matrix not symmetric. */
/**
* Default threshold above which off-diagonal elements are considered too different
* and matrix not symmetric.
*/
public static final double DEFAULT_RELATIVE_SYMMETRY_THRESHOLD = 1.0e-15;
/** Default threshold below which diagonal elements are considered null
* and matrix not positive definite. */
/**
* Default threshold below which diagonal elements are considered null
* and matrix not positive definite.
*/
public static final double DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD = 1.0e-10;
/** Row-oriented storage for L<sup>T</sup> matrix data. */
private double[][] lTData;
/** Cached value of L. */
private RealMatrix cachedL;
/** Cached value of LT. */
private RealMatrix cachedLT;
@ -62,17 +63,15 @@ public class CholeskyDecompositionImpl implements CholeskyDecomposition {
* #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD}
* </p>
* @param matrix the matrix to decompose
* @exception NonSquareMatrixException if matrix is not square
* @exception NotSymmetricMatrixException if matrix is not symmetric
* @exception NotPositiveDefiniteMatrixException if the matrix is not
* strictly positive definite
* @throws NonSquareMatrixException if the matrix is not square.
* @throws NonSymmetricMatrixException if the matrix is not symmetric.
* @throws NonPositiveDefiniteMatrixException if the matrix is not
* strictly positive definite.
* @see #CholeskyDecompositionImpl(RealMatrix, double, double)
* @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
* @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
*/
public CholeskyDecompositionImpl(final RealMatrix matrix)
throws NonSquareMatrixException,
NotSymmetricMatrixException, NotPositiveDefiniteMatrixException {
public CholeskyDecompositionImpl(final RealMatrix matrix) {
this(matrix, DEFAULT_RELATIVE_SYMMETRY_THRESHOLD,
DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD);
}
@ -84,20 +83,17 @@ public class CholeskyDecompositionImpl implements CholeskyDecomposition {
* elements are considered too different and matrix not symmetric
* @param absolutePositivityThreshold threshold below which diagonal
* elements are considered null and matrix not positive definite
* @exception NonSquareMatrixException if matrix is not square
* @exception NotSymmetricMatrixException if matrix is not symmetric
* @exception NotPositiveDefiniteMatrixException if the matrix is not
* strictly positive definite
* @throws NonSquareMatrixException if the matrix is not square.
* @throws NonSymmetricMatrixException if the matrix is not symmetric.
* @throws NonPositiveDefiniteMatrixException if the matrix is not
* strictly positive definite.
* @see #CholeskyDecompositionImpl(RealMatrix)
* @see #DEFAULT_RELATIVE_SYMMETRY_THRESHOLD
* @see #DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD
*/
public CholeskyDecompositionImpl(final RealMatrix matrix,
final double relativeSymmetryThreshold,
final double absolutePositivityThreshold)
throws NonSquareMatrixException,
NotSymmetricMatrixException, NotPositiveDefiniteMatrixException {
final double absolutePositivityThreshold) {
if (!matrix.isSquare()) {
throw new NonSquareMatrixException(matrix.getRowDimension(),
matrix.getColumnDimension());
@ -110,7 +106,6 @@ public class CholeskyDecompositionImpl implements CholeskyDecomposition {
// check the matrix before transformation
for (int i = 0; i < order; ++i) {
final double[] lI = lTData[i];
// check off-diagonal elements (and reset them to 0)
@ -121,7 +116,7 @@ public class CholeskyDecompositionImpl implements CholeskyDecomposition {
final double maxDelta =
relativeSymmetryThreshold * FastMath.max(FastMath.abs(lIJ), FastMath.abs(lJI));
if (FastMath.abs(lIJ - lJI) > maxDelta) {
throw new NotSymmetricMatrixException();
throw new NonSymmetricMatrixException(i, j, relativeSymmetryThreshold);
}
lJ[i] = 0;
}
@ -134,7 +129,7 @@ public class CholeskyDecompositionImpl implements CholeskyDecomposition {
// check diagonal element
if (ltI[i] < absolutePositivityThreshold) {
throw new NotPositiveDefiniteMatrixException();
throw new NonPositiveDefiniteMatrixException(i, absolutePositivityThreshold);
}
ltI[i] = FastMath.sqrt(ltI[i]);
@ -147,9 +142,7 @@ public class CholeskyDecompositionImpl implements CholeskyDecomposition {
ltQ[p] -= ltI[q] * ltI[p];
}
}
}
}
/** {@inheritDoc} */
@ -169,7 +162,6 @@ public class CholeskyDecompositionImpl implements CholeskyDecomposition {
// return the cached matrix
return cachedLT;
}
/** {@inheritDoc} */
@ -189,7 +181,6 @@ public class CholeskyDecompositionImpl implements CholeskyDecomposition {
/** Specialized solver. */
private static class Solver implements DecompositionSolver {
/** Row-oriented storage for L<sup>T</sup> matrix data. */
private final double[][] lTData;
@ -208,14 +199,10 @@ public class CholeskyDecompositionImpl implements CholeskyDecomposition {
}
/** {@inheritDoc} */
public double[] solve(double[] b)
throws IllegalArgumentException, InvalidMatrixException {
public double[] solve(double[] b) {
final int m = lTData.length;
if (b.length != m) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.VECTOR_LENGTH_MISMATCH,
b.length, m);
throw new DimensionMismatchException(b.length, m);
}
final double[] x = b.clone();
@ -240,21 +227,17 @@ public class CholeskyDecompositionImpl implements CholeskyDecomposition {
}
return x;
}
/** {@inheritDoc} */
public RealVector solve(RealVector b)
throws IllegalArgumentException, InvalidMatrixException {
public RealVector solve(RealVector b) {
try {
return solve((ArrayRealVector) b);
} catch (ClassCastException cce) {
final int m = lTData.length;
if (b.getDimension() != m) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.VECTOR_LENGTH_MISMATCH,
b.getDimension(), m);
throw new DimensionMismatchException(b.getDimension(), m);
}
final double[] x = b.getData();
@ -279,7 +262,6 @@ public class CholeskyDecompositionImpl implements CholeskyDecomposition {
}
return new ArrayRealVector(x, false);
}
}
@ -287,23 +269,18 @@ public class CholeskyDecompositionImpl implements CholeskyDecomposition {
* <p>The A matrix is implicit here. It is </p>
* @param b right-hand side of the equation A &times; X = B
* @return a vector X such that A &times; X = B
* @exception IllegalArgumentException if matrices dimensions don't match
* @exception InvalidMatrixException if decomposed matrix is singular
* @throws IllegalArgumentException if matrices dimensions don't match
* @throws InvalidMatrixException if decomposed matrix is singular
*/
public ArrayRealVector solve(ArrayRealVector b)
throws IllegalArgumentException, InvalidMatrixException {
public ArrayRealVector solve(ArrayRealVector b) {
return new ArrayRealVector(solve(b.getDataRef()), false);
}
/** {@inheritDoc} */
public RealMatrix solve(RealMatrix b)
throws IllegalArgumentException, InvalidMatrixException {
public RealMatrix solve(RealMatrix b) {
final int m = lTData.length;
if (b.getRowDimension() != m) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
b.getRowDimension(), b.getColumnDimension(), m, "n");
throw new DimensionMismatchException(b.getRowDimension(), m);
}
final int nColB = b.getColumnDimension();
@ -347,10 +324,8 @@ public class CholeskyDecompositionImpl implements CholeskyDecomposition {
}
/** {@inheritDoc} */
public RealMatrix getInverse() throws InvalidMatrixException {
public RealMatrix getInverse() {
return solve(MatrixUtils.createRealIdentityMatrix(lTData.length));
}
}
}

150
src/main/java/org/apache/commons/math/linear/EigenDecompositionImpl.java
View File

@ -17,8 +17,10 @@
package org.apache.commons.math.linear;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.exception.MaxCountExceededException;
import org.apache.commons.math.exception.SingularMatrixException;
import org.apache.commons.math.exception.NonSymmetricMatrixException;
import org.apache.commons.math.exception.DimensionMismatchException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.util.MathUtils;
import org.apache.commons.math.util.FastMath;
@ -88,39 +90,33 @@ public class EigenDecompositionImpl implements EigenDecomposition {
/**
* Calculates the eigen decomposition of the given symmetric matrix.
* @param matrix The <strong>symmetric</strong> matrix to decompose.
* @param splitTolerance dummy parameter, present for backward compatibility only.
* @exception InvalidMatrixException (wrapping a
* {@link org.apache.commons.math.ConvergenceException} if algorithm
* fails to converge
*
* @param matrix Matrix to decompose. It <em>must</em> be symmetric.
* @param splitTolerance Dummy parameter (present for backward
* compatibility only).
* @throws NonSymmetricMatrixException if the matrix is not symmetric.
* @throws MaxCountExceededException if the algorithm fails to converge.
*/
public EigenDecompositionImpl(final RealMatrix matrix,final double splitTolerance)
throws InvalidMatrixException {
if (isSymmetric(matrix)) {
public EigenDecompositionImpl(final RealMatrix matrix,
final double splitTolerance) {
if (isSymmetric(matrix, true)) {
transformToTridiagonal(matrix);
findEigenVectors(transformer.getQ().getData());
} else {
// as of 2.0, non-symmetric matrices (i.e. complex eigenvalues) are
// NOT supported
// see issue https://issues.apache.org/jira/browse/MATH-235
throw new InvalidMatrixException(
LocalizedFormats.ASSYMETRIC_EIGEN_NOT_SUPPORTED);
}
}
/**
* Calculates the eigen decomposition of the symmetric tridiagonal
* matrix. The Householder matrix is assumed to be the identity matrix.
*
* @param main Main diagonal of the symmetric triadiagonal form
* @param secondary Secondary of the tridiagonal form
* @param splitTolerance dummy parameter, present for backward compatibility only.
* @exception InvalidMatrixException (wrapping a
* {@link org.apache.commons.math.ConvergenceException} if algorithm
* fails to converge
* @param splitTolerance Dummy parameter (present for backward
* compatibility only).
* @throws MaxCountExceededException if the algorithm fails to converge.
*/
public EigenDecompositionImpl(final double[] main,final double[] secondary,
final double splitTolerance)
throws InvalidMatrixException {
final double splitTolerance) {
this.main = main.clone();
this.secondary = secondary.clone();
transformer = null;
@ -134,11 +130,16 @@ public class EigenDecompositionImpl implements EigenDecomposition {
/**
* Check if a matrix is symmetric.
* @param matrix
* matrix to check
* @return true if matrix is symmetric
*
* @param matrix Matrix to check.
* @param raiseException If {@code true}, the method will throw an
* exception if {@code matrix} is not symmetric.
* @return {@code true} if {@code matrix} is symmetric.
* @throws NonSymmetricMatrixException if the matrix is not symmetric and
* {@code raiseException} is {@code true}.
*/
private boolean isSymmetric(final RealMatrix matrix) {
private boolean isSymmetric(final RealMatrix matrix,
boolean raiseException) {
final int rows = matrix.getRowDimension();
final int columns = matrix.getColumnDimension();
final double eps = 10 * rows * columns * MathUtils.EPSILON;
@ -148,6 +149,9 @@ public class EigenDecompositionImpl implements EigenDecomposition {
final double mji = matrix.getEntry(j, i);
if (FastMath.abs(mij - mji) >
(FastMath.max(FastMath.abs(mij), FastMath.abs(mji)) * eps)) {
if (raiseException) {
throw new NonSymmetricMatrixException(i, j, eps);
}
return false;
}
}
@ -156,7 +160,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
}
/** {@inheritDoc} */
public RealMatrix getV() throws InvalidMatrixException {
public RealMatrix getV() {
if (cachedV == null) {
final int m = eigenvectors.length;
@ -171,7 +175,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
}
/** {@inheritDoc} */
public RealMatrix getD() throws InvalidMatrixException {
public RealMatrix getD() {
if (cachedD == null) {
// cache the matrix for subsequent calls
cachedD = MatrixUtils.createRealDiagonalMatrix(realEigenvalues);
@ -180,7 +184,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
}
/** {@inheritDoc} */
public RealMatrix getVT() throws InvalidMatrixException {
public RealMatrix getVT() {
if (cachedVt == null) {
final int m = eigenvectors.length;
@ -196,30 +200,27 @@ public class EigenDecompositionImpl implements EigenDecomposition {
}
/** {@inheritDoc} */
public double[] getRealEigenvalues() throws InvalidMatrixException {
public double[] getRealEigenvalues() {
return realEigenvalues.clone();
}
/** {@inheritDoc} */
public double getRealEigenvalue(final int i) throws InvalidMatrixException,
ArrayIndexOutOfBoundsException {
public double getRealEigenvalue(final int i) {
return realEigenvalues[i];
}
/** {@inheritDoc} */
public double[] getImagEigenvalues() throws InvalidMatrixException {
public double[] getImagEigenvalues() {
return imagEigenvalues.clone();
}
/** {@inheritDoc} */
public double getImagEigenvalue(final int i) throws InvalidMatrixException,
ArrayIndexOutOfBoundsException {
public double getImagEigenvalue(final int i) {
return imagEigenvalues[i];
}
/** {@inheritDoc} */
public RealVector getEigenvector(final int i)
throws InvalidMatrixException, ArrayIndexOutOfBoundsException {
public RealVector getEigenvector(final int i) {
return eigenvectors[i].copy();
}
@ -283,8 +284,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
* @exception InvalidMatrixException
* if decomposed matrix is singular
*/
public double[] solve(final double[] b)
throws IllegalArgumentException, InvalidMatrixException {
public double[] solve(final double[] b) {
if (!isNonSingular()) {
throw new SingularMatrixException();
@ -292,9 +292,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
final int m = realEigenvalues.length;
if (b.length != m) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.VECTOR_LENGTH_MISMATCH,
b.length, m);
throw new DimensionMismatchException(b.length, m);
}
final double[] bp = new double[m];
@ -317,26 +315,19 @@ public class EigenDecompositionImpl implements EigenDecomposition {
* This method only find exact linear solutions, i.e. solutions for
* which ||A &times; X - B|| is exactly 0.
* </p>
* @param b
* right-hand side of the equation A &times; X = B
* @return a vector X that minimizes the two norm of A &times; X - B
* @exception IllegalArgumentException
* if matrices dimensions don't match
* @exception InvalidMatrixException
* if decomposed matrix is singular
* @param b Right-hand side of the equation A &times; X = B
* @return a Vector X that minimizes the two norm of A &times; X - B
* @throws DimensionMismatchException if the matrices dimensions do not match.
* @throws SingularMatrixException if the decomposed matrix is singular.
*/
public RealVector solve(final RealVector b)
throws IllegalArgumentException, InvalidMatrixException {
public RealVector solve(final RealVector b) {
if (!isNonSingular()) {
throw new SingularMatrixException();
}
final int m = realEigenvalues.length;
if (b.getDimension() != m) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.VECTOR_LENGTH_MISMATCH, b
.getDimension(), m);
throw new DimensionMismatchException(b.getDimension(), m);
}
final double[] bp = new double[m];
@ -350,7 +341,6 @@ public class EigenDecompositionImpl implements EigenDecomposition {
}
return new ArrayRealVector(bp, false);
}
/**
@ -359,16 +349,12 @@ public class EigenDecompositionImpl implements EigenDecomposition {
* This method only find exact linear solutions, i.e. solutions for
* which ||A &times; X - B|| is exactly 0.
* </p>
* @param b
* right-hand side of the equation A &times; X = B
* @return a matrix X that minimizes the two norm of A &times; X - B
* @exception IllegalArgumentException
* if matrices dimensions don't match
* @exception InvalidMatrixException
* if decomposed matrix is singular
* @param b Right-hand side of the equation A &times; X = B
* @return a Matrix X that minimizes the two norm of A &times; X - B
* @throws DimensionMismatchException if the matrices dimensions do not match.
* @throws SingularMatrixException if the decomposed matrix is singular.
*/
public RealMatrix solve(final RealMatrix b)
throws IllegalArgumentException, InvalidMatrixException {
public RealMatrix solve(final RealMatrix b) {
if (!isNonSingular()) {
throw new SingularMatrixException();
@ -376,11 +362,7 @@ public class EigenDecompositionImpl implements EigenDecomposition {
final int m = realEigenvalues.length;
if (b.getRowDimension() != m) {
throw MathRuntimeException
.createIllegalArgumentException(
LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
b.getRowDimension(), b.getColumnDimension(), m,
"n");
throw new DimensionMismatchException(b.getRowDimension(), m);
}
final int nColB = b.getColumnDimension();
@ -419,12 +401,11 @@ public class EigenDecompositionImpl implements EigenDecomposition {
/**
* Get the inverse of the decomposed matrix.
* @return inverse matrix
* @throws InvalidMatrixException
* if decomposed matrix is singular
*
* @return the inverse matrix.
* @throws SingularMatrixException if the decomposed matrix is singular.
*/
public RealMatrix getInverse() throws InvalidMatrixException {
public RealMatrix getInverse() {
if (!isNonSingular()) {
throw new SingularMatrixException();
}
@ -444,32 +425,28 @@ public class EigenDecompositionImpl implements EigenDecomposition {
}
}
return MatrixUtils.createRealMatrix(invData);
}
}
/**
* Transform matrix to tridiagonal.
* @param matrix
* matrix to transform
*
* @param matrix Matrix to transform.
*/
private void transformToTridiagonal(final RealMatrix matrix) {
// transform the matrix to tridiagonal
transformer = new TriDiagonalTransformer(matrix);
main = transformer.getMainDiagonalRef();
secondary = transformer.getSecondaryDiagonalRef();
}
/**
* Find eigenvalues and eigenvectors (Dubrulle et al., 1971)
*
* @param householderMatrix Householder matrix of the transformation
* to tri-diagonal form.
* to tri-diagonal form.
*/
private void findEigenVectors(double[][] householderMatrix) {
double[][]z = householderMatrix.clone();
final int n = main.length;
realEigenvalues = new double[n];
@ -515,9 +492,10 @@ public class EigenDecompositionImpl implements EigenDecomposition {
}
}
if (m != j) {
if (its == maxIter)
throw new InvalidMatrixException(
new MaxIterationsExceededException(maxIter));
if (its == maxIter) {
throw new MaxCountExceededException(LocalizedFormats.CONVERGENCE_FAILED,
maxIter);
}
its++;
double q = (realEigenvalues[j + 1] - realEigenvalues[j]) / (2 * e[j]);
double t = FastMath.sqrt(1 + q * q);

46
src/main/java/org/apache/commons/math/linear/FieldLUDecompositionImpl.java
View File

@ -21,8 +21,9 @@ import java.lang.reflect.Array;
import org.apache.commons.math.Field;
import org.apache.commons.math.FieldElement;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.exception.DimensionMismatchException;
import org.apache.commons.math.exception.NonSquareMatrixException;
import org.apache.commons.math.exception.SingularMatrixException;
/**
* Calculates the LUP-decomposition of a square matrix.
@ -67,13 +68,12 @@ public class FieldLUDecompositionImpl<T extends FieldElement<T>> implements Fiel
/**
* Calculates the LU-decomposition of the given matrix.
* @param matrix The matrix to decompose.
* @exception NonSquareMatrixException if matrix is not square
* @throws NonSquareMatrixException if matrix is not square
*/
public FieldLUDecompositionImpl(FieldMatrix<T> matrix)
throws NonSquareMatrixException {
public FieldLUDecompositionImpl(FieldMatrix<T> matrix) {
if (!matrix.isSquare()) {
throw new NonSquareMatrixException(matrix.getRowDimension(), matrix.getColumnDimension());
throw new NonSquareMatrixException(matrix.getRowDimension(),
matrix.getColumnDimension());
}
final int m = matrix.getColumnDimension();
@ -258,14 +258,10 @@ public class FieldLUDecompositionImpl<T extends FieldElement<T>> implements Fiel
}
/** {@inheritDoc} */
public T[] solve(T[] b)
throws IllegalArgumentException, InvalidMatrixException {
public T[] solve(T[] b) {
final int m = pivot.length;
if (b.length != m) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.VECTOR_LENGTH_MISMATCH,
b.length, m);
throw new DimensionMismatchException(b.length, m);
}
if (singular) {
throw new SingularMatrixException();
@ -301,17 +297,14 @@ public class FieldLUDecompositionImpl<T extends FieldElement<T>> implements Fiel
}
/** {@inheritDoc} */
public FieldVector<T> solve(FieldVector<T> b)
throws IllegalArgumentException, InvalidMatrixException {
public FieldVector<T> solve(FieldVector<T> b) {
try {
return solve((ArrayFieldVector<T>) b);
} catch (ClassCastException cce) {
final int m = pivot.length;
if (b.getDimension() != m) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.VECTOR_LENGTH_MISMATCH,
b.getDimension(), m);
throw new DimensionMismatchException(b.getDimension(), m);
}
if (singular) {
throw new SingularMatrixException();
@ -351,23 +344,18 @@ public class FieldLUDecompositionImpl<T extends FieldElement<T>> implements Fiel
* <p>The A matrix is implicit here. It is </p>
* @param b right-hand side of the equation A &times; X = B
* @return a vector X such that A &times; X = B
* @exception IllegalArgumentException if matrices dimensions don't match
* @exception InvalidMatrixException if decomposed matrix is singular
* @throws IllegalArgumentException if matrices dimensions don't match
* @throws InvalidMatrixException if decomposed matrix is singular
*/
public ArrayFieldVector<T> solve(ArrayFieldVector<T> b)
throws IllegalArgumentException, InvalidMatrixException {
public ArrayFieldVector<T> solve(ArrayFieldVector<T> b) {
return new ArrayFieldVector<T>(solve(b.getDataRef()), false);
}
/** {@inheritDoc} */
public FieldMatrix<T> solve(FieldMatrix<T> b)
throws IllegalArgumentException, InvalidMatrixException {
public FieldMatrix<T> solve(FieldMatrix<T> b) {
final int m = pivot.length;
if (b.getRowDimension() != m) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
b.getRowDimension(), b.getColumnDimension(), m, "n");
throw new DimensionMismatchException(b.getRowDimension(), m);
}
if (singular) {
throw new SingularMatrixException();
@ -428,7 +416,5 @@ public class FieldLUDecompositionImpl<T extends FieldElement<T>> implements Fiel
}
return solve(identity);
}
}
}

3
src/main/java/org/apache/commons/math/linear/FieldMatrix.java
View File

@ -424,7 +424,8 @@ public interface FieldMatrix<T extends FieldElement<T>> extends AnyMatrix {
* trace</a> of the matrix (the sum of the elements on the main diagonal).
*
* @return trace
* @throws NonSquareMatrixException if the matrix is not square
* @throws org.apache.commons.math.exception.NonSquareMatrixException
* if the matrix is not square.
*/
T getTrace();

58
src/main/java/org/apache/commons/math/linear/LUDecompositionImpl.java
View File

@ -17,8 +17,9 @@
package org.apache.commons.math.linear;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.exception.DimensionMismatchException;
import org.apache.commons.math.exception.NonSquareMatrixException;
import org.apache.commons.math.exception.SingularMatrixException;
import org.apache.commons.math.util.FastMath;
/**
@ -34,35 +35,27 @@ import org.apache.commons.math.util.FastMath;
* @since 2.0
*/
public class LUDecompositionImpl implements LUDecomposition {
/** Default bound to determine effective singularity in LU decomposition */
private static final double DEFAULT_TOO_SMALL = 10E-12;
/** Entries of LU decomposition. */
private double lu[][];
/** Pivot permutation associated with LU decomposition */
private int[] pivot;
/** Parity of the permutation associated with the LU decomposition */
private boolean even;
/** Singularity indicator. */
private boolean singular;
/** Cached value of L. */
private RealMatrix cachedL;
/** Cached value of U. */
private RealMatrix cachedU;
/** Cached value of P. */
private RealMatrix cachedP;
/**
* Calculates the LU-decomposition of the given matrix.
* @param matrix The matrix to decompose.
* @exception InvalidMatrixException if matrix is not square
* @throws InvalidMatrixException if matrix is not square
*/
public LUDecompositionImpl(RealMatrix matrix)
throws InvalidMatrixException {
@ -74,13 +67,12 @@ public class LUDecompositionImpl implements LUDecomposition {
* @param matrix The matrix to decompose.
* @param singularityThreshold threshold (based on partial row norm)
* under which a matrix is considered singular
* @exception NonSquareMatrixException if matrix is not square
* @throws NonSquareMatrixException if matrix is not square
*/
public LUDecompositionImpl(RealMatrix matrix, double singularityThreshold)
throws NonSquareMatrixException {
public LUDecompositionImpl(RealMatrix matrix, double singularityThreshold) {
if (!matrix.isSquare()) {
throw new NonSquareMatrixException(matrix.getRowDimension(), matrix.getColumnDimension());
throw new NonSquareMatrixException(matrix.getRowDimension(),
matrix.getColumnDimension());
}
final int m = matrix.getColumnDimension();
@ -158,7 +150,6 @@ public class LUDecompositionImpl implements LUDecomposition {
lu[row][col] /= luDiag;
}
}
}
/** {@inheritDoc} */
@ -258,13 +249,10 @@ public class LUDecompositionImpl implements LUDecomposition {
}
/** {@inheritDoc} */
public double[] solve(double[] b)
throws IllegalArgumentException, InvalidMatrixException {
public double[] solve(double[] b) {
final int m = pivot.length;
if (b.length != m) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.VECTOR_LENGTH_MISMATCH, b.length, m);
throw new DimensionMismatchException(b.length, m);
}
if (singular) {
throw new SingularMatrixException();
@ -295,20 +283,17 @@ public class LUDecompositionImpl implements LUDecomposition {
}
return bp;
}
/** {@inheritDoc} */
public RealVector solve(RealVector b)
throws IllegalArgumentException, InvalidMatrixException {
public RealVector solve(RealVector b) {
try {
return solve((ArrayRealVector) b);
} catch (ClassCastException cce) {
final int m = pivot.length;
if (b.getDimension() != m) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.VECTOR_LENGTH_MISMATCH, b.getDimension(), m);
throw new DimensionMismatchException(b.getDimension(), m);
}
if (singular) {
throw new SingularMatrixException();
@ -339,7 +324,6 @@ public class LUDecompositionImpl implements LUDecomposition {
}
return new ArrayRealVector(bp, false);
}
}
@ -347,23 +331,20 @@ public class LUDecompositionImpl implements LUDecomposition {
* <p>The A matrix is implicit here. It is </p>
* @param b right-hand side of the equation A &times; X = B
* @return a vector X such that A &times; X = B
* @exception IllegalArgumentException if matrices dimensions don't match
* @exception InvalidMatrixException if decomposed matrix is singular
* @throws DimensionMismatchException if the matrices dimensions
* do not match.
* @throws InvalidMatrixException if decomposed matrix is singular
*/
public ArrayRealVector solve(ArrayRealVector b)
throws IllegalArgumentException, InvalidMatrixException {
public ArrayRealVector solve(ArrayRealVector b) {
return new ArrayRealVector(solve(b.getDataRef()), false);
}
/** {@inheritDoc} */
public RealMatrix solve(RealMatrix b)
throws IllegalArgumentException, InvalidMatrixException {
public RealMatrix solve(RealMatrix b) {
final int m = pivot.length;
if (b.getRowDimension() != m) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
b.getRowDimension(), b.getColumnDimension(), m, "n");
throw new DimensionMismatchException(b.getRowDimension(), m);
}
if (singular) {
throw new SingularMatrixException();
@ -410,14 +391,11 @@ public class LUDecompositionImpl implements LUDecomposition {
}
return new Array2DRowRealMatrix(bp, false);
}
/** {@inheritDoc} */
public RealMatrix getInverse() throws InvalidMatrixException {
return solve(MatrixUtils.createRealIdentityMatrix(pivot.length));
}
}
}

2
src/main/java/org/apache/commons/math/linear/NonSquareMatrixException.java
View File

@ -24,6 +24,8 @@ import org.apache.commons.math.exception.util.LocalizedFormats;
* Thrown when an operation defined only for square matrices is applied to non-square ones.
* @version $Revision$ $Date$
* @since 2.0
* @deprecated in 2.2 (to be removed in 3.0). Please use
* {@link org.apache.commons.math.exception.DimensionMismatchException} instead.
*/
public class NonSquareMatrixException extends InvalidMatrixException {

2
src/main/java/org/apache/commons/math/linear/NotPositiveDefiniteMatrixException.java
View File

@ -26,6 +26,8 @@ import org.apache.commons.math.exception.util.LocalizedFormats;
*
* @since 1.2
* @version $Revision$ $Date$
* @deprecated in 2.2 (to be removed in 3.0). Please use
* {@link org.apache.commons.math.exception.NonPositiveDefiniteMatrixException} instead.
*/
public class NotPositiveDefiniteMatrixException extends MathException {

2
src/main/java/org/apache/commons/math/linear/NotSymmetricMatrixException.java
View File

@ -26,6 +26,8 @@ import org.apache.commons.math.exception.util.LocalizedFormats;
*
* @since 2.0
* @version $Revision$ $Date$
* @deprecated in 2.2 (to be removed in 3.0). Please use
* {@link org.apache.commons.math.exception.NonSymmetricMatrixException} instead.
*/
public class NotSymmetricMatrixException extends MathException {

46
src/main/java/org/apache/commons/math/linear/QRDecompositionImpl.java
View File

@ -19,8 +19,8 @@ package org.apache.commons.math.linear;
import java.util.Arrays;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.exception.DimensionMismatchException;
import org.apache.commons.math.exception.SingularMatrixException;
import org.apache.commons.math.util.FastMath;
@ -162,12 +162,10 @@ public class QRDecompositionImpl implements QRDecomposition {
cachedR.setEntry(row, col, qrt[col][row]);
}
}
}
// return the cached matrix
return cachedR;
}
/** {@inheritDoc} */
@ -180,7 +178,6 @@ public class QRDecompositionImpl implements QRDecomposition {
/** {@inheritDoc} */
public RealMatrix getQT() {
if (cachedQT == null) {
// QT is supposed to be m x m
@ -214,17 +211,14 @@ public class QRDecompositionImpl implements QRDecomposition {
}
}
}
}
// return the cached matrix
return cachedQT;
}
/** {@inheritDoc} */
public RealMatrix getH() {
if (cachedH == null) {
final int n = qrt.length;
@ -235,12 +229,10 @@ public class QRDecompositionImpl implements QRDecomposition {
cachedH.setEntry(i, j, qrt[j][i] / -rDiag[j]);
}
}
}
// return the cached matrix
return cachedH;
}
/** {@inheritDoc} */
@ -281,19 +273,14 @@ public class QRDecompositionImpl implements QRDecomposition {
}
}
return true;
}
/** {@inheritDoc} */
public double[] solve(double[] b)
throws IllegalArgumentException, InvalidMatrixException {
public double[] solve(double[] b) {
final int n = qrt.length;
final int m = qrt[0].length;
if (b.length != m) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.VECTOR_LENGTH_MISMATCH,
b.length, m);
throw new DimensionMismatchException(b.length, m);
}
if (!isNonSingular()) {
throw new SingularMatrixException();
@ -315,7 +302,6 @@ public class QRDecompositionImpl implements QRDecomposition {
for (int row = minor; row < m; row++) {
y[row] += dotProduct * qrtMinor[row];
}
}
// solve triangular system R.x = y
@ -330,12 +316,10 @@ public class QRDecompositionImpl implements QRDecomposition {
}
return x;
}
/** {@inheritDoc} */
public RealVector solve(RealVector b)
throws IllegalArgumentException, InvalidMatrixException {
public RealVector solve(RealVector b) {
try {
return solve((ArrayRealVector) b);
} catch (ClassCastException cce) {
@ -350,21 +334,16 @@ public class QRDecompositionImpl implements QRDecomposition {
* @throws IllegalArgumentException if matrices dimensions don't match
* @throws InvalidMatrixException if decomposed matrix is singular
*/
public ArrayRealVector solve(ArrayRealVector b)
throws IllegalArgumentException, InvalidMatrixException {
public ArrayRealVector solve(ArrayRealVector b) {
return new ArrayRealVector(solve(b.getDataRef()), false);
}
/** {@inheritDoc} */
public RealMatrix solve(RealMatrix b)
throws IllegalArgumentException, InvalidMatrixException {
public RealMatrix solve(RealMatrix b) {
final int n = qrt.length;
final int m = qrt[0].length;
if (b.getRowDimension() != m) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.DIMENSIONS_MISMATCH_2x2,
b.getRowDimension(), b.getColumnDimension(), m, "n");
throw new DimensionMismatchException(b.getRowDimension(), m);
}
if (!isNonSingular()) {
throw new SingularMatrixException();
@ -409,7 +388,6 @@ public class QRDecompositionImpl implements QRDecomposition {
yRow[k] += alpha[k] * d;
}
}
}
// solve triangular system R.x = y
@ -433,21 +411,15 @@ public class QRDecompositionImpl implements QRDecomposition {
yI[k] -= yJ[k] * rIJ;
}
}
}
}
return new BlockRealMatrix(n, columns, xBlocks, false);
}
/** {@inheritDoc} */
public RealMatrix getInverse()
throws InvalidMatrixException {
public RealMatrix getInverse() {
return solve(MatrixUtils.createRealIdentityMatrix(rDiag.length));
}
}
}

5
src/main/java/org/apache/commons/math/linear/RealMatrix.java
View File

@ -466,8 +466,9 @@ public interface RealMatrix extends AnyMatrix {
* Returns the <a href="http://mathworld.wolfram.com/MatrixTrace.html">
* trace</a> of the matrix (the sum of the elements on the main diagonal).
*
* @return trace
* @throws NonSquareMatrixException if the matrix is not square
* @return the trace.
* @throws org.apache.commons.math.exception.NonSquareMatrixException
* if the matrix is not square.
*/
double getTrace();

2
src/main/java/org/apache/commons/math/linear/SingularMatrixException.java
View File

@ -24,6 +24,8 @@ import org.apache.commons.math.exception.util.LocalizedFormats;
* Thrown when a matrix is singular.
* @version $Revision$ $Date$
* @since 2.0
* @deprecated in 2.2 (to be removed in 3.0). Please use
* {@link org.apache.commons.math.exception.SingularMatrixException} instead.
*/
public class SingularMatrixException extends InvalidMatrixException {

27
src/main/java/org/apache/commons/math/linear/TriDiagonalTransformer.java
View File

@ -19,6 +19,7 @@ package org.apache.commons.math.linear;
import java.util.Arrays;
import org.apache.commons.math.exception.NonSquareMatrixException;
import org.apache.commons.math.util.FastMath;
@ -38,22 +39,16 @@ import org.apache.commons.math.util.FastMath;
* @since 2.0
*/
class TriDiagonalTransformer {
/** Householder vectors. */
private final double householderVectors[][];
/** Main diagonal. */
private final double[] main;
/** Secondary diagonal. */
private final double[] secondary;
/** Cached value of Q. */
private RealMatrix cachedQ;
/** Cached value of Qt. */
private RealMatrix cachedQt;
/** Cached value of T. */
private RealMatrix cachedT;
@ -61,11 +56,11 @@ class TriDiagonalTransformer {
* Build the transformation to tridiagonal shape of a symmetrical matrix.
* <p>The specified matrix is assumed to be symmetrical without any check.
* Only the upper triangular part of the matrix is used.</p>
* @param matrix the symmetrical matrix to transform.
* @exception InvalidMatrixException if matrix is not square
*
* @param matrix Symmetrical matrix to transform.
* @exception NonSquareMatrixException if the matrix is not square.
*/
public TriDiagonalTransformer(RealMatrix matrix)
throws InvalidMatrixException {
public TriDiagonalTransformer(RealMatrix matrix) {
if (!matrix.isSquare()) {
throw new NonSquareMatrixException(matrix.getRowDimension(), matrix.getColumnDimension());
}
@ -80,7 +75,6 @@ class TriDiagonalTransformer {
// transform matrix
transform();
}
/**
@ -101,9 +95,7 @@ class TriDiagonalTransformer {
* @return the Q matrix
*/
public RealMatrix getQT() {
if (cachedQt == null) {
final int m = householderVectors.length;
cachedQt = MatrixUtils.createRealMatrix(m, m);
@ -132,12 +124,10 @@ class TriDiagonalTransformer {
}
}
cachedQt.setEntry(0, 0, 1);
}
// return the cached matrix
return cachedQt;
}
/**
@ -145,9 +135,7 @@ class TriDiagonalTransformer {
* @return the T matrix
*/
public RealMatrix getT() {
if (cachedT == null) {
final int m = main.length;
cachedT = MatrixUtils.createRealMatrix(m, m);
for (int i = 0; i < m; ++i) {
@ -164,7 +152,6 @@ class TriDiagonalTransformer {
// return the cached matrix
return cachedT;
}
/**
@ -202,7 +189,6 @@ class TriDiagonalTransformer {
* <p>Transformation is done using Householder transforms.</p>
*/
private void transform() {
final int m = householderVectors.length;
final double[] z = new double[m];
for (int k = 0; k < m - 1; k++) {
@ -260,11 +246,8 @@ class TriDiagonalTransformer {
hI[j] -= hK[i] * z[j] + z[i] * hK[j];
}
}
}
}
main[m - 1] = householderVectors[m - 1][m - 1];
}
}

85
src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java
View File

@ -17,12 +17,11 @@
package org.apache.commons.math.optimization.general;
import org.apache.commons.math.exception.FunctionEvaluationException;
import org.apache.commons.math.exception.ConvergenceException;
import org.apache.commons.math.exception.NumberIsTooSmallException;
import org.apache.commons.math.exception.DimensionMismatchException;
import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction;
import org.apache.commons.math.analysis.MultivariateMatrixFunction;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.linear.InvalidMatrixException;
import org.apache.commons.math.linear.LUDecompositionImpl;
import org.apache.commons.math.linear.MatrixUtils;
import org.apache.commons.math.linear.RealMatrix;
@ -95,15 +94,17 @@ public abstract class AbstractLeastSquaresOptimizer
/**
* Update the jacobian matrix.
* @exception FunctionEvaluationException if the function jacobian
* cannot be evaluated or its dimension doesn't match problem dimension
*
* @throws org.apache.commons.math.exception.FunctionEvaluationException
* if the function Jacobian cannot be evaluated.
* @throws DimensionMismatchException if the Jacobian dimension does not
* match problem dimension.
*/
protected void updateJacobian() throws FunctionEvaluationException {
protected void updateJacobian() {
++jacobianEvaluations;
weightedResidualJacobian = jF.value(point);
if (weightedResidualJacobian.length != rows) {
throw new FunctionEvaluationException(point, LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
weightedResidualJacobian.length, rows);
throw new DimensionMismatchException(weightedResidualJacobian.length, rows);
}
final double[] residualsWeights = getWeightRef();
@ -120,15 +121,17 @@ public abstract class AbstractLeastSquaresOptimizer
/**
* Update the residuals array and cost function value.
* @exception FunctionEvaluationException if the function cannot be evaluated
* or its dimension doesn't match problem dimension or maximal number of
* of evaluations is exceeded
* @throws org.apache.commons.math.exception.FunctionEvaluationException
* if the function cannot be evaluated.
* @throws DimensionMismatchException if the dimension does not match the
* problem dimension.
* @throws org.apache.commons.math.exception.TooManyEvaluationsException
* if the maximal number of evaluations is exceeded.
*/
protected void updateResidualsAndCost() throws FunctionEvaluationException {
protected void updateResidualsAndCost() {
objective = computeObjectiveValue(point);
if (objective.length != rows) {
throw new FunctionEvaluationException(point, LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE,
objective.length, rows);
throw new DimensionMismatchException(objective.length, rows);
}
final double[] targetValues = getTargetRef();
@ -170,16 +173,15 @@ public abstract class AbstractLeastSquaresOptimizer
}
/**
* Get the covariance matrix of optimized parameters.
* @return covariance matrix
* @exception FunctionEvaluationException if the function jacobian cannot
* be evaluated
* @exception ConvergenceException if the covariance matrix
* cannot be computed (singular problem)
* Get the covariance matrix of the optimized parameters.
*
* @return the covariance matrix.
* @throws org.apache.commons.math.exception.FunctionEvaluationException
* if the function Jacobian cannot be evaluated.
* @throws org.apache.commons.math.exception.SingularMatrixException
* if the covariance matrix cannot be computed (singular problem).
*/
public double[][] getCovariances()
throws FunctionEvaluationException {
public double[][] getCovariances() {
// set up the jacobian
updateJacobian();
@ -196,31 +198,29 @@ public abstract class AbstractLeastSquaresOptimizer
}
}
try {
// compute the covariances matrix
RealMatrix inverse =
new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
return inverse.getData();
} catch (InvalidMatrixException ime) {
throw new ConvergenceException(LocalizedFormats.UNABLE_TO_COMPUTE_COVARIANCE_SINGULAR_PROBLEM);
}
// compute the covariances matrix
RealMatrix inverse =
new LUDecompositionImpl(MatrixUtils.createRealMatrix(jTj)).getSolver().getInverse();
return inverse.getData();
}
/**
* Guess the errors in optimized parameters.
* <p>Guessing is covariance-based, it only gives rough order of magnitude.</p>
* Guessing is covariance-based: It only gives a rough order of magnitude.
*
* @return errors in optimized parameters
* @exception FunctionEvaluationException if the function jacobian cannot b evaluated
* @exception ConvergenceException if the covariances matrix cannot be computed
* or the number of degrees of freedom is not positive (number of measurements
* lesser or equal to number of parameters)
* @throws org.apache.commons.math.exception.FunctionEvaluationException
* if the function Jacobian cannot be evaluated.
* @throws org.apache.commons.math.exception.SingularMatrixException if
* the covariances matrix cannot be computed.
* @throws NumberIsTooSmallException if the number of degrees of freedom is not
* positive, i.e. the number of measurements is less or equal to the number of
* parameters.
*/
public double[] guessParametersErrors()
throws FunctionEvaluationException {
public double[] guessParametersErrors() {
if (rows <= cols) {
throw new ConvergenceException(LocalizedFormats.NO_DEGREES_OF_FREEDOM,
rows, cols);
throw new NumberIsTooSmallException(LocalizedFormats.NO_DEGREES_OF_FREEDOM,
rows, cols, false);
}
double[] errors = new double[cols];
final double c = FastMath.sqrt(getChiSquare() / (rows - cols));
@ -235,8 +235,7 @@ public abstract class AbstractLeastSquaresOptimizer
@Override
public VectorialPointValuePair optimize(final DifferentiableMultivariateVectorialFunction f,
final double[] target, final double[] weights,
final double[] startPoint)
throws FunctionEvaluationException {
final double[] startPoint) {
// Reset counter.
jacobianEvaluations = 0;

6
src/main/java/org/apache/commons/math/optimization/general/GaussNewtonOptimizer.java
View File

@ -17,15 +17,15 @@
package org.apache.commons.math.optimization.general;
import org.apache.commons.math.exception.SingularMatrixException;
import org.apache.commons.math.exception.FunctionEvaluationException;
import org.apache.commons.math.exception.ConvergenceException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.linear.BlockRealMatrix;
import org.apache.commons.math.linear.DecompositionSolver;
import org.apache.commons.math.linear.InvalidMatrixException;
import org.apache.commons.math.linear.LUDecompositionImpl;
import org.apache.commons.math.linear.QRDecompositionImpl;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.exception.ConvergenceException;
import org.apache.commons.math.optimization.VectorialPointValuePair;
import org.apache.commons.math.optimization.ConvergenceChecker;
@ -120,7 +120,7 @@ public class GaussNewtonOptimizer extends AbstractLeastSquaresOptimizer {
for (int i = 0; i < cols; ++i) {
point[i] += dX[i];
}
} catch (InvalidMatrixException e) {
} catch (SingularMatrixException e) {
throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM);
}

2
src/main/resources/META-INF/localization/LocalizedFormats_fr.properties
View File

@ -153,6 +153,7 @@ ALPHA = alpha
BETA = beta
NOT_POSITIVE_COLUMNDIMENSION = nombre de colonnes invalide : {0} (doit \u00eatre positif)
NOT_POSITIVE_DEFINITE_MATRIX = matrice non d\u00e9finie positive
NON_POSITIVE_DEFINITE_MATRIX = matrice non d\u00e9finie positive: l''\u00e9l\u00e9ment diagonal ({0},{0}) est sup\u00e9rieur \u0061 {1}
NOT_POSITIVE_DEGREES_OF_FREEDOM = les degr\u00e9s de libert\u00e9 doivent \u00eatre positifs ({0})
DEGREES_OF_FREEDOM = degr\u00e9s de libert\u00e9 ({0})
NOT_POSITIVE_ELEMENT_AT_INDEX = l''\u00e9l\u00e9ment {0} n''est pas positif : {1}
@ -189,6 +190,7 @@ NOT_STRICTLY_INCREASING_NUMBER_OF_POINTS = les points {0} et {1} ne sont pas str
NOT_STRICTLY_INCREASING_SEQUENCE = les points {3} et {2} ne sont pas strictement croissants ({1} >= {0})
NOT_SUBTRACTION_COMPATIBLE_MATRICES = les dimensions {0}x{1} et {2}x{3} sont incompatibles pour la soustraction matricielle
NOT_SYMMETRIC_MATRIX = matrice non symm\u00e9trique
NON_SYMMETRIC_MATRIX = matrice non symm\u00e9trique: la diff\u00e9rence entre les \u00e9l\u00e9ments ({0},{1}) et ({1},{0}) est sup\u00e9rieure \u0061 {2}
NO_BIN_SELECTED = aucun compartiment s\u00e9lectionn\u00e9
NO_CONVERGENCE_WITH_ANY_START_POINT = aucun des {0} points de d\u00e9part n''aboutit \u00e0 une convergence
NO_DATA = aucune donn\u00e9e

74
src/test/java/org/apache/commons/math/estimation/EstimatedParameterTest.java
View File

@ -1,74 +0,0 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.estimation;
import org.apache.commons.math.estimation.EstimatedParameter;
import org.apache.commons.math.util.FastMath;
import junit.framework.*;
@Deprecated
public class EstimatedParameterTest
extends TestCase {
public EstimatedParameterTest(String name) {
super(name);
}
public void testConstruction() {
EstimatedParameter p1 = new EstimatedParameter("p1", 1.0);
assertTrue(p1.getName().equals("p1"));
checkValue(p1.getEstimate(), 1.0);
assertTrue(! p1.isBound());
EstimatedParameter p2 = new EstimatedParameter("p2", 2.0, true);
assertTrue(p2.getName().equals("p2"));
checkValue(p2.getEstimate(), 2.0);
assertTrue(p2.isBound());
}
public void testBound() {
EstimatedParameter p = new EstimatedParameter("p", 0.0);
assertTrue(! p.isBound());
p.setBound(true);
assertTrue(p.isBound());
p.setBound(false);
assertTrue(! p.isBound());
}
public void testEstimate() {
EstimatedParameter p = new EstimatedParameter("p", 0.0);
checkValue(p.getEstimate(), 0.0);
for (double e = 0.0; e < 10.0; e += 0.5) {
p.setEstimate(e);
checkValue(p.getEstimate(), e);
}
}
private void checkValue(double value, double expected) {
assertTrue(FastMath.abs(value - expected) < 1.0e-10);
}
}

724
src/test/java/org/apache/commons/math/estimation/GaussNewtonEstimatorTest.java
View File

@ -1,724 +0,0 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.estimation;
import java.util.ArrayList;
import java.util.HashSet;
import org.apache.commons.math.util.FastMath;
import junit.framework.TestCase;
/**
* <p>Some of the unit tests are re-implementations of the MINPACK <a
* href="http://www.netlib.org/minpack/ex/file17">file17</a> and <a
* href="http://www.netlib.org/minpack/ex/file22">file22</a> test files.
* The redistribution policy for MINPACK is available <a
* href="http://www.netlib.org/minpack/disclaimer">here</a>, for
* convenience, it is reproduced below.</p>
* <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
* <tr><td>
* Minpack Copyright Notice (1999) University of Chicago.
* All rights reserved
* </td></tr>
* <tr><td>
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* <ol>
* <li>Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.</li>
* <li>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.</li>
* <li>The end-user documentation included with the redistribution, if any,
* must include the following acknowledgment:
* <code>This product includes software developed by the University of
* Chicago, as Operator of Argonne National Laboratory.</code>
* Alternately, this acknowledgment may appear in the software itself,
* if and wherever such third-party acknowledgments normally appear.</li>
* <li><strong>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.</strong></li>
* <li><strong>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.</strong></li>
* <ol></td></tr>
* </table>
* @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran minpack tests)
* @author Burton S. Garbow (original fortran minpack tests)
* @author Kenneth E. Hillstrom (original fortran minpack tests)
* @author Jorge J. More (original fortran minpack tests)
* @author Luc Maisonobe (non-minpack tests and minpack tests Java translation)
*/
@Deprecated
public class GaussNewtonEstimatorTest
extends TestCase {
public GaussNewtonEstimatorTest(String name) {
super(name);
}
public void testTrivial() throws EstimationException {
LinearProblem problem =
new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] {2},
new EstimatedParameter[] {
new EstimatedParameter("p0", 0)
}, 3.0)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
assertEquals(1.5,
problem.getUnboundParameters()[0].getEstimate(),
1.0e-10);
}
public void testQRColumnsPermutation() throws EstimationException {
EstimatedParameter[] x = {
new EstimatedParameter("p0", 0), new EstimatedParameter("p1", 0)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { x[0], x[1] },
4.0),
new LinearMeasurement(new double[] { 2.0 },
new EstimatedParameter[] { x[1] },
6.0),
new LinearMeasurement(new double[] { 1.0, -2.0 },
new EstimatedParameter[] { x[0], x[1] },
1.0)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
assertEquals(7.0, x[0].getEstimate(), 1.0e-10);
assertEquals(3.0, x[1].getEstimate(), 1.0e-10);
}
public void testNoDependency() throws EstimationException {
EstimatedParameter[] p = new EstimatedParameter[] {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 0),
new EstimatedParameter("p2", 0),
new EstimatedParameter("p3", 0),
new EstimatedParameter("p4", 0),
new EstimatedParameter("p5", 0)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[0] }, 0.0),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[1] }, 1.1),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[2] }, 2.2),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[3] }, 3.3),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[4] }, 4.4),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[5] }, 5.5)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
for (int i = 0; i < p.length; ++i) {
assertEquals(0.55 * i, p[i].getEstimate(), 1.0e-10);
}
}
public void testOneSet() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 0),
new EstimatedParameter("p2", 0)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0 },
new EstimatedParameter[] { p[0] },
1.0),
new LinearMeasurement(new double[] { -1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1] },
1.0),
new LinearMeasurement(new double[] { -1.0, 1.0 },
new EstimatedParameter[] { p[1], p[2] },
1.0)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
assertEquals(1.0, p[0].getEstimate(), 1.0e-10);
assertEquals(2.0, p[1].getEstimate(), 1.0e-10);
assertEquals(3.0, p[2].getEstimate(), 1.0e-10);
}
public void testTwoSets() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 1),
new EstimatedParameter("p2", 2),
new EstimatedParameter("p3", 3),
new EstimatedParameter("p4", 4),
new EstimatedParameter("p5", 5)
};
double epsilon = 1.0e-7;
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
// 4 elements sub-problem
new LinearMeasurement(new double[] { 2.0, 1.0, 4.0 },
new EstimatedParameter[] { p[0], p[1], p[3] },
2.0),
new LinearMeasurement(new double[] { -4.0, -2.0, 3.0, -7.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
-9.0),
new LinearMeasurement(new double[] { 4.0, 1.0, -2.0, 8.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
2.0),
new LinearMeasurement(new double[] { -3.0, -12.0, -1.0 },
new EstimatedParameter[] { p[1], p[2], p[3] },
2.0),
// 2 elements sub-problem
new LinearMeasurement(new double[] { epsilon, 1.0 },
new EstimatedParameter[] { p[4], p[5] },
1.0 + epsilon * epsilon),
new LinearMeasurement(new double[] { 1.0, 1.0 },
new EstimatedParameter[] { p[4], p[5] },
2.0)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
assertEquals( 3.0, p[0].getEstimate(), 1.0e-10);
assertEquals( 4.0, p[1].getEstimate(), 1.0e-10);
assertEquals(-1.0, p[2].getEstimate(), 1.0e-10);
assertEquals(-2.0, p[3].getEstimate(), 1.0e-10);
assertEquals( 1.0 + epsilon, p[4].getEstimate(), 1.0e-10);
assertEquals( 1.0 - epsilon, p[5].getEstimate(), 1.0e-10);
}
public void testNonInversible() {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 0),
new EstimatedParameter("p2", 0)
};
LinearMeasurement[] m = new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, 2.0, -3.0 },
new EstimatedParameter[] { p[0], p[1], p[2] },
1.0),
new LinearMeasurement(new double[] { 2.0, 1.0, 3.0 },
new EstimatedParameter[] { p[0], p[1], p[2] },
1.0),
new LinearMeasurement(new double[] { -3.0, -9.0 },
new EstimatedParameter[] { p[0], p[2] },
1.0)
};
LinearProblem problem = new LinearProblem(m);
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
try {
estimator.estimate(problem);
fail("an exception should have been caught");
} catch (EstimationException ee) {
// expected behavior
}
}
public void testIllConditioned() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 1),
new EstimatedParameter("p2", 2),
new EstimatedParameter("p3", 3)
};
LinearProblem problem1 = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 10.0, 7.0, 8.0, 7.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
32.0),
new LinearMeasurement(new double[] { 7.0, 5.0, 6.0, 5.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
23.0),
new LinearMeasurement(new double[] { 8.0, 6.0, 10.0, 9.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
33.0),
new LinearMeasurement(new double[] { 7.0, 5.0, 9.0, 10.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
31.0)
});
GaussNewtonEstimator estimator1 = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator1.estimate(problem1);
assertEquals(0, estimator1.getRMS(problem1), 1.0e-10);
assertEquals(1.0, p[0].getEstimate(), 1.0e-10);
assertEquals(1.0, p[1].getEstimate(), 1.0e-10);
assertEquals(1.0, p[2].getEstimate(), 1.0e-10);
assertEquals(1.0, p[3].getEstimate(), 1.0e-10);
LinearProblem problem2 = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 10.0, 7.0, 8.1, 7.2 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
32.0),
new LinearMeasurement(new double[] { 7.08, 5.04, 6.0, 5.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
23.0),
new LinearMeasurement(new double[] { 8.0, 5.98, 9.89, 9.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
33.0),
new LinearMeasurement(new double[] { 6.99, 4.99, 9.0, 9.98 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
31.0)
});
GaussNewtonEstimator estimator2 = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator2.estimate(problem2);
assertEquals(0, estimator2.getRMS(problem2), 1.0e-10);
assertEquals(-81.0, p[0].getEstimate(), 1.0e-8);
assertEquals(137.0, p[1].getEstimate(), 1.0e-8);
assertEquals(-34.0, p[2].getEstimate(), 1.0e-8);
assertEquals( 22.0, p[3].getEstimate(), 1.0e-8);
}
public void testMoreEstimatedParametersSimple() {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 7),
new EstimatedParameter("p1", 6),
new EstimatedParameter("p2", 5),
new EstimatedParameter("p3", 4)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 3.0, 2.0 },
new EstimatedParameter[] { p[0], p[1] },
7.0),
new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 },
new EstimatedParameter[] { p[1], p[2], p[3] },
3.0),
new LinearMeasurement(new double[] { 2.0, 1.0 },
new EstimatedParameter[] { p[0], p[2] },
5.0)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
try {
estimator.estimate(problem);
fail("an exception should have been caught");
} catch (EstimationException ee) {
// expected behavior
}
}
public void testMoreEstimatedParametersUnsorted() {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 2),
new EstimatedParameter("p1", 2),
new EstimatedParameter("p2", 2),
new EstimatedParameter("p3", 2),
new EstimatedParameter("p4", 2),
new EstimatedParameter("p5", 2)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1] },
3.0),
new LinearMeasurement(new double[] { 1.0, 1.0, 1.0 },
new EstimatedParameter[] { p[2], p[3], p[4] },
12.0),
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { p[4], p[5] },
-1.0),
new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 },
new EstimatedParameter[] { p[3], p[2], p[5] },
7.0),
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { p[4], p[3] },
1.0)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
try {
estimator.estimate(problem);
fail("an exception should have been caught");
} catch (EstimationException ee) {
// expected behavior
}
}
public void testRedundantEquations() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 1),
new EstimatedParameter("p1", 1)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1] },
3.0),
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { p[0], p[1] },
1.0),
new LinearMeasurement(new double[] { 1.0, 3.0 },
new EstimatedParameter[] { p[0], p[1] },
5.0)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
EstimatedParameter[] all = problem.getAllParameters();
for (int i = 0; i < all.length; ++i) {
assertEquals(all[i].getName().equals("p0") ? 2.0 : 1.0,
all[i].getEstimate(), 1.0e-10);
}
}
public void testInconsistentEquations() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 1),
new EstimatedParameter("p1", 1)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1] },
3.0),
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { p[0], p[1] },
1.0),
new LinearMeasurement(new double[] { 1.0, 3.0 },
new EstimatedParameter[] { p[0], p[1] },
4.0)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator.estimate(problem);
assertTrue(estimator.getRMS(problem) > 0.1);
}
public void testBoundParameters() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("unbound0", 2, false),
new EstimatedParameter("unbound1", 2, false),
new EstimatedParameter("bound", 2, true)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, 1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1], p[2] },
3.0),
new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1], p[2] },
1.0),
new LinearMeasurement(new double[] { 1.0, 3.0, 2.0 },
new EstimatedParameter[] { p[0], p[1], p[2] },
7.0)
});
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
estimator.estimate(problem);
assertTrue(estimator.getRMS(problem) < 1.0e-10);
double[][] covariances = estimator.getCovariances(problem);
int i0 = 0, i1 = 1;
if (problem.getUnboundParameters()[0].getName().endsWith("1")) {
i0 = 1;
i1 = 0;
}
assertEquals(11.0 / 24, covariances[i0][i0], 1.0e-10);
assertEquals(-3.0 / 24, covariances[i0][i1], 1.0e-10);
assertEquals(-3.0 / 24, covariances[i1][i0], 1.0e-10);
assertEquals( 3.0 / 24, covariances[i1][i1], 1.0e-10);
double[] errors = estimator.guessParametersErrors(problem);
assertEquals(0, errors[i0], 1.0e-10);
assertEquals(0, errors[i1], 1.0e-10);
}
public void testMaxIterations() {
Circle circle = new Circle(98.680, 47.345);
circle.addPoint( 30.0, 68.0);
circle.addPoint( 50.0, -6.0);
circle.addPoint(110.0, -20.0);
circle.addPoint( 35.0, 15.0);
circle.addPoint( 45.0, 97.0);
try {
GaussNewtonEstimator estimator = new GaussNewtonEstimator(4, 1.0e-14, 1.0e-14);
estimator.estimate(circle);
fail("an exception should have been caught");
} catch (EstimationException ee) {
// expected behavior
}
}
public void testCircleFitting() throws EstimationException {
Circle circle = new Circle(98.680, 47.345);
circle.addPoint( 30.0, 68.0);
circle.addPoint( 50.0, -6.0);
circle.addPoint(110.0, -20.0);
circle.addPoint( 35.0, 15.0);
circle.addPoint( 45.0, 97.0);
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-10, 1.0e-10);
estimator.estimate(circle);
double rms = estimator.getRMS(circle);
assertEquals(1.768262623567235, FastMath.sqrt(circle.getM()) * rms, 1.0e-10);
assertEquals(69.96016176931406, circle.getRadius(), 1.0e-10);
assertEquals(96.07590211815305, circle.getX(), 1.0e-10);
assertEquals(48.13516790438953, circle.getY(), 1.0e-10);
}
public void testCircleFittingBadInit() {
Circle circle = new Circle(-12, -12);
double[][] points = new double[][] {
{-0.312967, 0.072366}, {-0.339248, 0.132965}, {-0.379780, 0.202724},
{-0.390426, 0.260487}, {-0.361212, 0.328325}, {-0.346039, 0.392619},
{-0.280579, 0.444306}, {-0.216035, 0.470009}, {-0.149127, 0.493832},
{-0.075133, 0.483271}, {-0.007759, 0.452680}, { 0.060071, 0.410235},
{ 0.103037, 0.341076}, { 0.118438, 0.273884}, { 0.131293, 0.192201},
{ 0.115869, 0.129797}, { 0.072223, 0.058396}, { 0.022884, 0.000718},
{-0.053355, -0.020405}, {-0.123584, -0.032451}, {-0.216248, -0.032862},
{-0.278592, -0.005008}, {-0.337655, 0.056658}, {-0.385899, 0.112526},
{-0.405517, 0.186957}, {-0.415374, 0.262071}, {-0.387482, 0.343398},
{-0.347322, 0.397943}, {-0.287623, 0.458425}, {-0.223502, 0.475513},
{-0.135352, 0.478186}, {-0.061221, 0.483371}, { 0.003711, 0.422737},
{ 0.065054, 0.375830}, { 0.108108, 0.297099}, { 0.123882, 0.222850},
{ 0.117729, 0.134382}, { 0.085195, 0.056820}, { 0.029800, -0.019138},
{-0.027520, -0.072374}, {-0.102268, -0.091555}, {-0.200299, -0.106578},
{-0.292731, -0.091473}, {-0.356288, -0.051108}, {-0.420561, 0.014926},
{-0.471036, 0.074716}, {-0.488638, 0.182508}, {-0.485990, 0.254068},
{-0.463943, 0.338438}, {-0.406453, 0.404704}, {-0.334287, 0.466119},
{-0.254244, 0.503188}, {-0.161548, 0.495769}, {-0.075733, 0.495560},
{ 0.001375, 0.434937}, { 0.082787, 0.385806}, { 0.115490, 0.323807},
{ 0.141089, 0.223450}, { 0.138693, 0.131703}, { 0.126415, 0.049174},
{ 0.066518, -0.010217}, {-0.005184, -0.070647}, {-0.080985, -0.103635},
{-0.177377, -0.116887}, {-0.260628, -0.100258}, {-0.335756, -0.056251},
{-0.405195, -0.000895}, {-0.444937, 0.085456}, {-0.484357, 0.175597},
{-0.472453, 0.248681}, {-0.438580, 0.347463}, {-0.402304, 0.422428},
{-0.326777, 0.479438}, {-0.247797, 0.505581}, {-0.152676, 0.519380},
{-0.071754, 0.516264}, { 0.015942, 0.472802}, { 0.076608, 0.419077},
{ 0.127673, 0.330264}, { 0.159951, 0.262150}, { 0.153530, 0.172681},
{ 0.140653, 0.089229}, { 0.078666, 0.024981}, { 0.023807, -0.037022},
{-0.048837, -0.077056}, {-0.127729, -0.075338}, {-0.221271, -0.067526}
};
for (int i = 0; i < points.length; ++i) {
circle.addPoint(points[i][0], points[i][1]);
}
GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6);
try {
estimator.estimate(circle);
fail("an exception should have been caught");
} catch (EstimationException ee) {
// expected behavior
}
}
private static class LinearProblem extends SimpleEstimationProblem {
public LinearProblem(LinearMeasurement[] measurements) {
HashSet<EstimatedParameter> set = new HashSet<EstimatedParameter>();
for (int i = 0; i < measurements.length; ++i) {
addMeasurement(measurements[i]);
EstimatedParameter[] parameters = measurements[i].getParameters();
for (int j = 0; j < parameters.length; ++j) {
set.add(parameters[j]);
}
}
for (EstimatedParameter p : set) {
addParameter(p);
}
}
}
private static class LinearMeasurement extends WeightedMeasurement {
public LinearMeasurement(double[] factors, EstimatedParameter[] parameters,
double setPoint) {
super(1.0, setPoint, true);
this.factors = factors;
this.parameters = parameters;
setIgnored(false);
}
@Override
public double getTheoreticalValue() {
double v = 0;
for (int i = 0; i < factors.length; ++i) {
v += factors[i] * parameters[i].getEstimate();
}
return v;
}
@Override
public double getPartial(EstimatedParameter parameter) {
for (int i = 0; i < parameters.length; ++i) {
if (parameters[i] == parameter) {
return factors[i];
}
}
return 0;
}
public EstimatedParameter[] getParameters() {
return parameters;
}
private double[] factors;
private EstimatedParameter[] parameters;
private static final long serialVersionUID = -3922448707008868580L;
}
private static class Circle implements EstimationProblem {
public Circle(double cx, double cy) {
this.cx = new EstimatedParameter("cx", cx);
this.cy = new EstimatedParameter(new EstimatedParameter("cy", cy));
points = new ArrayList<PointModel>();
}
public void addPoint(double px, double py) {
points.add(new PointModel(this, px, py));
}
public int getM() {
return points.size();
}
public WeightedMeasurement[] getMeasurements() {
return points.toArray(new PointModel[points.size()]);
}
public EstimatedParameter[] getAllParameters() {
return new EstimatedParameter[] { cx, cy };
}
public EstimatedParameter[] getUnboundParameters() {
return new EstimatedParameter[] { cx, cy };
}
public double getPartialRadiusX() {
double dRdX = 0;
for (PointModel point : points) {
dRdX += point.getPartialDiX();
}
return dRdX / points.size();
}
public double getPartialRadiusY() {
double dRdY = 0;
for (PointModel point : points) {
dRdY += point.getPartialDiY();
}
return dRdY / points.size();
}
public double getRadius() {
double r = 0;
for (PointModel point : points) {
r += point.getCenterDistance();
}
return r / points.size();
}
public double getX() {
return cx.getEstimate();
}
public double getY() {
return cy.getEstimate();
}
private static class PointModel extends WeightedMeasurement {
public PointModel(Circle circle, double px, double py) {
super(1.0, 0.0);
this.px = px;
this.py = py;
this.circle = circle;
}
@Override
public double getPartial(EstimatedParameter parameter) {
if (parameter == circle.cx) {
return getPartialDiX() - circle.getPartialRadiusX();
} else if (parameter == circle.cy) {
return getPartialDiY() - circle.getPartialRadiusY();
}
return 0;
}
public double getCenterDistance() {
double dx = px - circle.cx.getEstimate();
double dy = py - circle.cy.getEstimate();
return FastMath.sqrt(dx * dx + dy * dy);
}
public double getPartialDiX() {
return (circle.cx.getEstimate() - px) / getCenterDistance();
}
public double getPartialDiY() {
return (circle.cy.getEstimate() - py) / getCenterDistance();
}
@Override
public double getTheoreticalValue() {
return getCenterDistance() - circle.getRadius();
}
private double px;
private double py;
private transient final Circle circle;
private static final long serialVersionUID = 1L;
}
private EstimatedParameter cx;
private EstimatedParameter cy;
private ArrayList<PointModel> points;
}
}

832
src/test/java/org/apache/commons/math/estimation/LevenbergMarquardtEstimatorTest.java
View File

@ -1,832 +0,0 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.estimation;
import java.util.ArrayList;
import java.util.HashSet;
import org.apache.commons.math.util.FastMath;
import junit.framework.TestCase;
/**
* <p>Some of the unit tests are re-implementations of the MINPACK <a
* href="http://www.netlib.org/minpack/ex/file17">file17</a> and <a
* href="http://www.netlib.org/minpack/ex/file22">file22</a> test files.
* The redistribution policy for MINPACK is available <a
* href="http://www.netlib.org/minpack/disclaimer">here</a>, for
* convenience, it is reproduced below.</p>
* <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0">
* <tr><td>
* Minpack Copyright Notice (1999) University of Chicago.
* All rights reserved
* </td></tr>
* <tr><td>
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* <ol>
* <li>Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.</li>
* <li>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.</li>
* <li>The end-user documentation included with the redistribution, if any,
* must include the following acknowledgment:
* <code>This product includes software developed by the University of
* Chicago, as Operator of Argonne National Laboratory.</code>
* Alternately, this acknowledgment may appear in the software itself,
* if and wherever such third-party acknowledgments normally appear.</li>
* <li><strong>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.</strong></li>
* <li><strong>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.</strong></li>
* <ol></td></tr>
* </table>
* @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran minpack tests)
* @author Burton S. Garbow (original fortran minpack tests)
* @author Kenneth E. Hillstrom (original fortran minpack tests)
* @author Jorge J. More (original fortran minpack tests)
* @author Luc Maisonobe (non-minpack tests and minpack tests Java translation)
*/
@Deprecated
public class LevenbergMarquardtEstimatorTest
extends TestCase {
public LevenbergMarquardtEstimatorTest(String name) {
super(name);
}
public void testTrivial() throws EstimationException {
LinearProblem problem =
new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] {2},
new EstimatedParameter[] {
new EstimatedParameter("p0", 0)
}, 3.0)
});
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
try {
estimator.guessParametersErrors(problem);
fail("an exception should have been thrown");
} catch (EstimationException ee) {
// expected behavior
}
assertEquals(1.5,
problem.getUnboundParameters()[0].getEstimate(),
1.0e-10);
}
public void testQRColumnsPermutation() throws EstimationException {
EstimatedParameter[] x = {
new EstimatedParameter("p0", 0), new EstimatedParameter("p1", 0)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { x[0], x[1] },
4.0),
new LinearMeasurement(new double[] { 2.0 },
new EstimatedParameter[] { x[1] },
6.0),
new LinearMeasurement(new double[] { 1.0, -2.0 },
new EstimatedParameter[] { x[0], x[1] },
1.0)
});
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
assertEquals(7.0, x[0].getEstimate(), 1.0e-10);
assertEquals(3.0, x[1].getEstimate(), 1.0e-10);
}
public void testNoDependency() throws EstimationException {
EstimatedParameter[] p = new EstimatedParameter[] {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 0),
new EstimatedParameter("p2", 0),
new EstimatedParameter("p3", 0),
new EstimatedParameter("p4", 0),
new EstimatedParameter("p5", 0)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[0] }, 0.0),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[1] }, 1.1),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[2] }, 2.2),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[3] }, 3.3),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[4] }, 4.4),
new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[5] }, 5.5)
});
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
for (int i = 0; i < p.length; ++i) {
assertEquals(0.55 * i, p[i].getEstimate(), 1.0e-10);
}
}
public void testOneSet() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 0),
new EstimatedParameter("p2", 0)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0 },
new EstimatedParameter[] { p[0] },
1.0),
new LinearMeasurement(new double[] { -1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1] },
1.0),
new LinearMeasurement(new double[] { -1.0, 1.0 },
new EstimatedParameter[] { p[1], p[2] },
1.0)
});
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
assertEquals(1.0, p[0].getEstimate(), 1.0e-10);
assertEquals(2.0, p[1].getEstimate(), 1.0e-10);
assertEquals(3.0, p[2].getEstimate(), 1.0e-10);
}
public void testTwoSets() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 1),
new EstimatedParameter("p2", 2),
new EstimatedParameter("p3", 3),
new EstimatedParameter("p4", 4),
new EstimatedParameter("p5", 5)
};
double epsilon = 1.0e-7;
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
// 4 elements sub-problem
new LinearMeasurement(new double[] { 2.0, 1.0, 4.0 },
new EstimatedParameter[] { p[0], p[1], p[3] },
2.0),
new LinearMeasurement(new double[] { -4.0, -2.0, 3.0, -7.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
-9.0),
new LinearMeasurement(new double[] { 4.0, 1.0, -2.0, 8.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
2.0),
new LinearMeasurement(new double[] { -3.0, -12.0, -1.0 },
new EstimatedParameter[] { p[1], p[2], p[3] },
2.0),
// 2 elements sub-problem
new LinearMeasurement(new double[] { epsilon, 1.0 },
new EstimatedParameter[] { p[4], p[5] },
1.0 + epsilon * epsilon),
new LinearMeasurement(new double[] { 1.0, 1.0 },
new EstimatedParameter[] { p[4], p[5] },
2.0)
});
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
assertEquals( 3.0, p[0].getEstimate(), 1.0e-10);
assertEquals( 4.0, p[1].getEstimate(), 1.0e-10);
assertEquals(-1.0, p[2].getEstimate(), 1.0e-10);
assertEquals(-2.0, p[3].getEstimate(), 1.0e-10);
assertEquals( 1.0 + epsilon, p[4].getEstimate(), 1.0e-10);
assertEquals( 1.0 - epsilon, p[5].getEstimate(), 1.0e-10);
}
public void testNonInversible() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 0),
new EstimatedParameter("p2", 0)
};
LinearMeasurement[] m = new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, 2.0, -3.0 },
new EstimatedParameter[] { p[0], p[1], p[2] },
1.0),
new LinearMeasurement(new double[] { 2.0, 1.0, 3.0 },
new EstimatedParameter[] { p[0], p[1], p[2] },
1.0),
new LinearMeasurement(new double[] { -3.0, -9.0 },
new EstimatedParameter[] { p[0], p[2] },
1.0)
};
LinearProblem problem = new LinearProblem(m);
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
double initialCost = estimator.getRMS(problem);
estimator.estimate(problem);
assertTrue(estimator.getRMS(problem) < initialCost);
assertTrue(FastMath.sqrt(m.length) * estimator.getRMS(problem) > 0.6);
try {
estimator.getCovariances(problem);
fail("an exception should have been thrown");
} catch (EstimationException ee) {
// expected behavior
}
double dJ0 = 2 * (m[0].getResidual() * m[0].getPartial(p[0])
+ m[1].getResidual() * m[1].getPartial(p[0])
+ m[2].getResidual() * m[2].getPartial(p[0]));
double dJ1 = 2 * (m[0].getResidual() * m[0].getPartial(p[1])
+ m[1].getResidual() * m[1].getPartial(p[1]));
double dJ2 = 2 * (m[0].getResidual() * m[0].getPartial(p[2])
+ m[1].getResidual() * m[1].getPartial(p[2])
+ m[2].getResidual() * m[2].getPartial(p[2]));
assertEquals(0, dJ0, 1.0e-10);
assertEquals(0, dJ1, 1.0e-10);
assertEquals(0, dJ2, 1.0e-10);
}
public void testIllConditioned() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 0),
new EstimatedParameter("p1", 1),
new EstimatedParameter("p2", 2),
new EstimatedParameter("p3", 3)
};
LinearProblem problem1 = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 10.0, 7.0, 8.0, 7.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
32.0),
new LinearMeasurement(new double[] { 7.0, 5.0, 6.0, 5.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
23.0),
new LinearMeasurement(new double[] { 8.0, 6.0, 10.0, 9.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
33.0),
new LinearMeasurement(new double[] { 7.0, 5.0, 9.0, 10.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
31.0)
});
LevenbergMarquardtEstimator estimator1 = new LevenbergMarquardtEstimator();
estimator1.estimate(problem1);
assertEquals(0, estimator1.getRMS(problem1), 1.0e-10);
assertEquals(1.0, p[0].getEstimate(), 1.0e-10);
assertEquals(1.0, p[1].getEstimate(), 1.0e-10);
assertEquals(1.0, p[2].getEstimate(), 1.0e-10);
assertEquals(1.0, p[3].getEstimate(), 1.0e-10);
LinearProblem problem2 = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 10.0, 7.0, 8.1, 7.2 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
32.0),
new LinearMeasurement(new double[] { 7.08, 5.04, 6.0, 5.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
23.0),
new LinearMeasurement(new double[] { 8.0, 5.98, 9.89, 9.0 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
33.0),
new LinearMeasurement(new double[] { 6.99, 4.99, 9.0, 9.98 },
new EstimatedParameter[] { p[0], p[1], p[2], p[3] },
31.0)
});
LevenbergMarquardtEstimator estimator2 = new LevenbergMarquardtEstimator();
estimator2.estimate(problem2);
assertEquals(0, estimator2.getRMS(problem2), 1.0e-10);
assertEquals(-81.0, p[0].getEstimate(), 1.0e-8);
assertEquals(137.0, p[1].getEstimate(), 1.0e-8);
assertEquals(-34.0, p[2].getEstimate(), 1.0e-8);
assertEquals( 22.0, p[3].getEstimate(), 1.0e-8);
}
public void testMoreEstimatedParametersSimple() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 7),
new EstimatedParameter("p1", 6),
new EstimatedParameter("p2", 5),
new EstimatedParameter("p3", 4)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 3.0, 2.0 },
new EstimatedParameter[] { p[0], p[1] },
7.0),
new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 },
new EstimatedParameter[] { p[1], p[2], p[3] },
3.0),
new LinearMeasurement(new double[] { 2.0, 1.0 },
new EstimatedParameter[] { p[0], p[2] },
5.0)
});
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
}
public void testMoreEstimatedParametersUnsorted() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 2),
new EstimatedParameter("p1", 2),
new EstimatedParameter("p2", 2),
new EstimatedParameter("p3", 2),
new EstimatedParameter("p4", 2),
new EstimatedParameter("p5", 2)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1] },
3.0),
new LinearMeasurement(new double[] { 1.0, 1.0, 1.0 },
new EstimatedParameter[] { p[2], p[3], p[4] },
12.0),
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { p[4], p[5] },
-1.0),
new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 },
new EstimatedParameter[] { p[3], p[2], p[5] },
7.0),
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { p[4], p[3] },
1.0)
});
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
assertEquals(3.0, p[2].getEstimate(), 1.0e-10);
assertEquals(4.0, p[3].getEstimate(), 1.0e-10);
assertEquals(5.0, p[4].getEstimate(), 1.0e-10);
assertEquals(6.0, p[5].getEstimate(), 1.0e-10);
}
public void testRedundantEquations() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 1),
new EstimatedParameter("p1", 1)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1] },
3.0),
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { p[0], p[1] },
1.0),
new LinearMeasurement(new double[] { 1.0, 3.0 },
new EstimatedParameter[] { p[0], p[1] },
5.0)
});
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(problem);
assertEquals(0, estimator.getRMS(problem), 1.0e-10);
assertEquals(2.0, p[0].getEstimate(), 1.0e-10);
assertEquals(1.0, p[1].getEstimate(), 1.0e-10);
}
public void testInconsistentEquations() throws EstimationException {
EstimatedParameter[] p = {
new EstimatedParameter("p0", 1),
new EstimatedParameter("p1", 1)
};
LinearProblem problem = new LinearProblem(new LinearMeasurement[] {
new LinearMeasurement(new double[] { 1.0, 1.0 },
new EstimatedParameter[] { p[0], p[1] },
3.0),
new LinearMeasurement(new double[] { 1.0, -1.0 },
new EstimatedParameter[] { p[0], p[1] },
1.0),
new LinearMeasurement(new double[] { 1.0, 3.0 },
new EstimatedParameter[] { p[0], p[1] },
4.0)
});
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(problem);
assertTrue(estimator.getRMS(problem) > 0.1);
}
public void testControlParameters() {
Circle circle = new Circle(98.680, 47.345);
circle.addPoint( 30.0, 68.0);
circle.addPoint( 50.0, -6.0);
circle.addPoint(110.0, -20.0);
circle.addPoint( 35.0, 15.0);
circle.addPoint( 45.0, 97.0);
checkEstimate(circle, 0.1, 10, 1.0e-14, 1.0e-16, 1.0e-10, false);
checkEstimate(circle, 0.1, 10, 1.0e-15, 1.0e-17, 1.0e-10, true);
checkEstimate(circle, 0.1, 5, 1.0e-15, 1.0e-16, 1.0e-10, true);
circle.addPoint(300, -300);
checkEstimate(circle, 0.1, 20, 1.0e-18, 1.0e-16, 1.0e-10, true);
}
private void checkEstimate(EstimationProblem problem,
double initialStepBoundFactor, int maxCostEval,
double costRelativeTolerance, double parRelativeTolerance,
double orthoTolerance, boolean shouldFail) {
try {
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.setInitialStepBoundFactor(initialStepBoundFactor);
estimator.setMaxCostEval(maxCostEval);
estimator.setCostRelativeTolerance(costRelativeTolerance);
estimator.setParRelativeTolerance(parRelativeTolerance);
estimator.setOrthoTolerance(orthoTolerance);
estimator.estimate(problem);
assertTrue(! shouldFail);
} catch (EstimationException ee) {
assertTrue(shouldFail);
}
}
public void testCircleFitting() throws EstimationException {
Circle circle = new Circle(98.680, 47.345);
circle.addPoint( 30.0, 68.0);
circle.addPoint( 50.0, -6.0);
circle.addPoint(110.0, -20.0);
circle.addPoint( 35.0, 15.0);
circle.addPoint( 45.0, 97.0);
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(circle);
assertTrue(estimator.getCostEvaluations() < 10);
assertTrue(estimator.getJacobianEvaluations() < 10);
double rms = estimator.getRMS(circle);
assertEquals(1.768262623567235, FastMath.sqrt(circle.getM()) * rms, 1.0e-10);
assertEquals(69.96016176931406, circle.getRadius(), 1.0e-10);
assertEquals(96.07590211815305, circle.getX(), 1.0e-10);
assertEquals(48.13516790438953, circle.getY(), 1.0e-10);
double[][] cov = estimator.getCovariances(circle);
assertEquals(1.839, cov[0][0], 0.001);
assertEquals(0.731, cov[0][1], 0.001);
assertEquals(cov[0][1], cov[1][0], 1.0e-14);
assertEquals(0.786, cov[1][1], 0.001);
double[] errors = estimator.guessParametersErrors(circle);
assertEquals(1.384, errors[0], 0.001);
assertEquals(0.905, errors[1], 0.001);
// add perfect measurements and check errors are reduced
double cx = circle.getX();
double cy = circle.getY();
double r = circle.getRadius();
for (double d= 0; d < 2 * FastMath.PI; d += 0.01) {
circle.addPoint(cx + r * FastMath.cos(d), cy + r * FastMath.sin(d));
}
estimator = new LevenbergMarquardtEstimator();
estimator.estimate(circle);
cov = estimator.getCovariances(circle);
assertEquals(0.004, cov[0][0], 0.001);
assertEquals(6.40e-7, cov[0][1], 1.0e-9);
assertEquals(cov[0][1], cov[1][0], 1.0e-14);
assertEquals(0.003, cov[1][1], 0.001);
errors = estimator.guessParametersErrors(circle);
assertEquals(0.004, errors[0], 0.001);
assertEquals(0.004, errors[1], 0.001);
}
public void testCircleFittingBadInit() throws EstimationException {
Circle circle = new Circle(-12, -12);
double[][] points = new double[][] {
{-0.312967, 0.072366}, {-0.339248, 0.132965}, {-0.379780, 0.202724},
{-0.390426, 0.260487}, {-0.361212, 0.328325}, {-0.346039, 0.392619},
{-0.280579, 0.444306}, {-0.216035, 0.470009}, {-0.149127, 0.493832},
{-0.075133, 0.483271}, {-0.007759, 0.452680}, { 0.060071, 0.410235},
{ 0.103037, 0.341076}, { 0.118438, 0.273884}, { 0.131293, 0.192201},
{ 0.115869, 0.129797}, { 0.072223, 0.058396}, { 0.022884, 0.000718},
{-0.053355, -0.020405}, {-0.123584, -0.032451}, {-0.216248, -0.032862},
{-0.278592, -0.005008}, {-0.337655, 0.056658}, {-0.385899, 0.112526},
{-0.405517, 0.186957}, {-0.415374, 0.262071}, {-0.387482, 0.343398},
{-0.347322, 0.397943}, {-0.287623, 0.458425}, {-0.223502, 0.475513},
{-0.135352, 0.478186}, {-0.061221, 0.483371}, { 0.003711, 0.422737},
{ 0.065054, 0.375830}, { 0.108108, 0.297099}, { 0.123882, 0.222850},
{ 0.117729, 0.134382}, { 0.085195, 0.056820}, { 0.029800, -0.019138},
{-0.027520, -0.072374}, {-0.102268, -0.091555}, {-0.200299, -0.106578},
{-0.292731, -0.091473}, {-0.356288, -0.051108}, {-0.420561, 0.014926},
{-0.471036, 0.074716}, {-0.488638, 0.182508}, {-0.485990, 0.254068},
{-0.463943, 0.338438}, {-0.406453, 0.404704}, {-0.334287, 0.466119},
{-0.254244, 0.503188}, {-0.161548, 0.495769}, {-0.075733, 0.495560},
{ 0.001375, 0.434937}, { 0.082787, 0.385806}, { 0.115490, 0.323807},
{ 0.141089, 0.223450}, { 0.138693, 0.131703}, { 0.126415, 0.049174},
{ 0.066518, -0.010217}, {-0.005184, -0.070647}, {-0.080985, -0.103635},
{-0.177377, -0.116887}, {-0.260628, -0.100258}, {-0.335756, -0.056251},
{-0.405195, -0.000895}, {-0.444937, 0.085456}, {-0.484357, 0.175597},
{-0.472453, 0.248681}, {-0.438580, 0.347463}, {-0.402304, 0.422428},
{-0.326777, 0.479438}, {-0.247797, 0.505581}, {-0.152676, 0.519380},
{-0.071754, 0.516264}, { 0.015942, 0.472802}, { 0.076608, 0.419077},
{ 0.127673, 0.330264}, { 0.159951, 0.262150}, { 0.153530, 0.172681},
{ 0.140653, 0.089229}, { 0.078666, 0.024981}, { 0.023807, -0.037022},
{-0.048837, -0.077056}, {-0.127729, -0.075338}, {-0.221271, -0.067526}
};
for (int i = 0; i < points.length; ++i) {
circle.addPoint(points[i][0], points[i][1]);
}
LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator();
estimator.estimate(circle);
assertTrue(estimator.getCostEvaluations() < 15);
assertTrue(estimator.getJacobianEvaluations() < 10);
assertEquals( 0.030184491196225207, estimator.getRMS(circle), 1.0e-9);
assertEquals( 0.2922350065939634, circle.getRadius(), 1.0e-9);
assertEquals(-0.15173845023862165, circle.getX(), 1.0e-8);
assertEquals( 0.20750021499570379, circle.getY(), 1.0e-8);
}
public void testMath199() {
try {
QuadraticProblem problem = new QuadraticProblem();
problem.addPoint (0, -3.182591015485607, 0.0);
problem.addPoint (1, -2.5581184967730577, 4.4E-323);
problem.addPoint (2, -2.1488478161387325, 1.0);
problem.addPoint (3, -1.9122489313410047, 4.4E-323);
problem.addPoint (4, 1.7785661310051026, 0.0);
new LevenbergMarquardtEstimator().estimate(problem);
fail("an exception should have been thrown");
} catch (EstimationException ee) {
// expected behavior
}
}
private static class LinearProblem implements EstimationProblem {
public LinearProblem(LinearMeasurement[] measurements) {
this.measurements = measurements;
}
public WeightedMeasurement[] getMeasurements() {
return measurements;
}
public EstimatedParameter[] getUnboundParameters() {
return getAllParameters();
}
public EstimatedParameter[] getAllParameters() {
HashSet<EstimatedParameter> set = new HashSet<EstimatedParameter>();
for (int i = 0; i < measurements.length; ++i) {
EstimatedParameter[] parameters = measurements[i].getParameters();
for (int j = 0; j < parameters.length; ++j) {
set.add(parameters[j]);
}
}
return set.toArray(new EstimatedParameter[set.size()]);
}
private LinearMeasurement[] measurements;
}
private static class LinearMeasurement extends WeightedMeasurement {
public LinearMeasurement(double[] factors, EstimatedParameter[] parameters,
double setPoint) {
super(1.0, setPoint);
this.factors = factors;
this.parameters = parameters;
}
@Override
public double getTheoreticalValue() {
double v = 0;
for (int i = 0; i < factors.length; ++i) {
v += factors[i] * parameters[i].getEstimate();
}
return v;
}
@Override
public double getPartial(EstimatedParameter parameter) {
for (int i = 0; i < parameters.length; ++i) {
if (parameters[i] == parameter) {
return factors[i];
}
}
return 0;
}
public EstimatedParameter[] getParameters() {
return parameters;
}
private double[] factors;
private EstimatedParameter[] parameters;
private static final long serialVersionUID = -3922448707008868580L;
}
private static class Circle implements EstimationProblem {
public Circle(double cx, double cy) {
this.cx = new EstimatedParameter("cx", cx);
this.cy = new EstimatedParameter("cy", cy);
points = new ArrayList<PointModel>();
}
public void addPoint(double px, double py) {
points.add(new PointModel(this, px, py));
}
public int getM() {
return points.size();
}
public WeightedMeasurement[] getMeasurements() {
return points.toArray(new PointModel[points.size()]);
}
public EstimatedParameter[] getAllParameters() {
return new EstimatedParameter[] { cx, cy };
}
public EstimatedParameter[] getUnboundParameters() {
return new EstimatedParameter[] { cx, cy };
}
public double getPartialRadiusX() {
double dRdX = 0;
for (PointModel point : points) {
dRdX += point.getPartialDiX();
}
return dRdX / points.size();
}
public double getPartialRadiusY() {
double dRdY = 0;
for (PointModel point : points) {
dRdY += point.getPartialDiY();
}
return dRdY / points.size();
}
public double getRadius() {
double r = 0;
for (PointModel point : points) {
r += point.getCenterDistance();
}
return r / points.size();
}
public double getX() {
return cx.getEstimate();
}
public double getY() {
return cy.getEstimate();
}
private static class PointModel extends WeightedMeasurement {
public PointModel(Circle circle, double px, double py) {
super(1.0, 0.0);
this.px = px;
this.py = py;
this.circle = circle;
}
@Override
public double getPartial(EstimatedParameter parameter) {
if (parameter == circle.cx) {
return getPartialDiX() - circle.getPartialRadiusX();
} else if (parameter == circle.cy) {
return getPartialDiY() - circle.getPartialRadiusY();
}
return 0;
}
public double getCenterDistance() {
double dx = px - circle.cx.getEstimate();
double dy = py - circle.cy.getEstimate();
return FastMath.sqrt(dx * dx + dy * dy);
}
public double getPartialDiX() {
return (circle.cx.getEstimate() - px) / getCenterDistance();
}
public double getPartialDiY() {
return (circle.cy.getEstimate() - py) / getCenterDistance();
}
@Override
public double getTheoreticalValue() {
return getCenterDistance() - circle.getRadius();
}
private double px;
private double py;
private transient final Circle circle;
private static final long serialVersionUID = 1L;
}
private EstimatedParameter cx;
private EstimatedParameter cy;
private ArrayList<PointModel> points;
}
private static class QuadraticProblem extends SimpleEstimationProblem {
private EstimatedParameter a;
private EstimatedParameter b;
private EstimatedParameter c;
public QuadraticProblem() {
a = new EstimatedParameter("a", 0.0);
b = new EstimatedParameter("b", 0.0);
c = new EstimatedParameter("c", 0.0);
addParameter(a);
addParameter(b);
addParameter(c);
}
public void addPoint(double x, double y, double w) {
addMeasurement(new LocalMeasurement(this, x, y, w));
}
public double theoreticalValue(double x) {
return ( (a.getEstimate() * x + b.getEstimate() ) * x + c.getEstimate());
}
private double partial(double x, EstimatedParameter parameter) {
if (parameter == a) {
return x * x;
} else if (parameter == b) {
return x;
} else {
return 1.0;
}
}
private static class LocalMeasurement extends WeightedMeasurement {
private static final long serialVersionUID = 1555043155023729130L;
private final double x;
private transient final QuadraticProblem pb;
// constructor
public LocalMeasurement(QuadraticProblem pb, double x, double y, double w) {
super(w, y);
this.x = x;
this.pb = pb;
}
@Override
public double getTheoreticalValue() {
return pb.theoreticalValue(x);
}
@Override
public double getPartial(EstimatedParameter parameter) {
return pb.partial(x, parameter);
}
}
}
}

1538
src/test/java/org/apache/commons/math/estimation/MinpackTest.java
View File

File diff suppressed because it is too large Load Diff

128
src/test/java/org/apache/commons/math/estimation/WeightedMeasurementTest.java
View File

@ -1,128 +0,0 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.estimation;
import org.apache.commons.math.estimation.EstimatedParameter;
import org.apache.commons.math.estimation.WeightedMeasurement;
import org.apache.commons.math.util.FastMath;
import junit.framework.*;
@Deprecated
public class WeightedMeasurementTest
extends TestCase {
public WeightedMeasurementTest(String name) {
super(name);
p1 = null;
p2 = null;
}
public void testConstruction() {
WeightedMeasurement m = new MyMeasurement(3.0, theoretical() + 0.1, this);
checkValue(m.getWeight(), 3.0);
checkValue(m.getMeasuredValue(), theoretical() + 0.1);
}
public void testIgnored() {
WeightedMeasurement m = new MyMeasurement(3.0, theoretical() + 0.1, this);
assertTrue(!m.isIgnored());
m.setIgnored(true);
assertTrue(m.isIgnored());
m.setIgnored(false);
assertTrue(!m.isIgnored());
}
public void testTheory() {
WeightedMeasurement m = new MyMeasurement(3.0, theoretical() + 0.1, this);
checkValue(m.getTheoreticalValue(), theoretical());
checkValue(m.getResidual(), 0.1);
double oldP1 = p1.getEstimate();
p1.setEstimate(oldP1 + m.getResidual() / m.getPartial(p1));
checkValue(m.getResidual(), 0.0);
p1.setEstimate(oldP1);
checkValue(m.getResidual(), 0.1);
double oldP2 = p2.getEstimate();
p2.setEstimate(oldP2 + m.getResidual() / m.getPartial(p2));
checkValue(m.getResidual(), 0.0);
p2.setEstimate(oldP2);
checkValue(m.getResidual(), 0.1);
}
@Override
public void setUp() {
p1 = new EstimatedParameter("p1", 1.0);
p2 = new EstimatedParameter("p2", 2.0);
}
@Override
public void tearDown() {
p1 = null;
p2 = null;
}
private void checkValue(double value, double expected) {
assertTrue(FastMath.abs(value - expected) < 1.0e-10);
}
private double theoretical() {
return 3 * p1.getEstimate() - p2.getEstimate();
}
private double partial(EstimatedParameter p) {
if (p == p1) {
return 3.0;
} else if (p == p2) {
return -1.0;
} else {
return 0.0;
}
}
private static class MyMeasurement
extends WeightedMeasurement {
public MyMeasurement(double weight, double measuredValue,
WeightedMeasurementTest testInstance) {
super(weight, measuredValue);
this.testInstance = testInstance;
}
@Override
public double getTheoreticalValue() {
return testInstance.theoretical();
}
@Override
public double getPartial(EstimatedParameter p) {
return testInstance.partial(p);
}
private transient WeightedMeasurementTest testInstance;
private static final long serialVersionUID = -246712922500792332L;
}
private EstimatedParameter p1;
private EstimatedParameter p2;
}

1
src/test/java/org/apache/commons/math/linear/Array2DRowRealMatrixTest.java
View File

@ -24,6 +24,7 @@ import org.apache.commons.math.exception.MatrixDimensionMismatchException;
import org.apache.commons.math.exception.OutOfRangeException;
import org.apache.commons.math.exception.ZeroException;
import org.apache.commons.math.exception.NumberIsTooSmallException;
import org.apache.commons.math.exception.NonSquareMatrixException;
/**
* Test cases for the {@link Array2DRowRealMatrix} class.

1
src/test/java/org/apache/commons/math/linear/BlockFieldMatrixTest.java
View File

@ -29,6 +29,7 @@ import org.apache.commons.math.exception.NoDataException;
import org.apache.commons.math.exception.OutOfRangeException;
import org.apache.commons.math.exception.NumberIsTooSmallException;
import org.apache.commons.math.exception.NotStrictlyPositiveException;
import org.apache.commons.math.exception.NonSquareMatrixException;
/**
* Test cases for the {@link BlockFieldMatrix} class.

1
src/test/java/org/apache/commons/math/linear/BlockRealMatrixTest.java
View File

@ -27,6 +27,7 @@ import org.apache.commons.math.exception.MatrixDimensionMismatchException;
import org.apache.commons.math.exception.OutOfRangeException;
import org.apache.commons.math.exception.ZeroException;
import org.apache.commons.math.exception.NumberIsTooSmallException;
import org.apache.commons.math.exception.NonSquareMatrixException;
/**
* Test cases for the {@link BlockRealMatrix} class.

16
src/test/java/org/apache/commons/math/linear/CholeskyDecompositionImplTest.java
View File

@ -21,13 +21,9 @@ import static org.junit.Assert.assertEquals;
import static org.junit.Assert.assertTrue;
import org.apache.commons.math.MathException;
import org.apache.commons.math.linear.CholeskyDecomposition;
import org.apache.commons.math.linear.CholeskyDecompositionImpl;
import org.apache.commons.math.linear.MatrixUtils;
import org.apache.commons.math.linear.NonSquareMatrixException;
import org.apache.commons.math.linear.NotPositiveDefiniteMatrixException;
import org.apache.commons.math.linear.NotSymmetricMatrixException;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.exception.NonSquareMatrixException;
import org.apache.commons.math.exception.NonPositiveDefiniteMatrixException;
import org.apache.commons.math.exception.NonSymmetricMatrixException;
import org.junit.Test;
public class CholeskyDecompositionImplTest {
@ -58,7 +54,7 @@ public class CholeskyDecompositionImplTest {
}
/** test non-symmetric matrix */
@Test(expected = NotSymmetricMatrixException.class)
@Test(expected = NonSymmetricMatrixException.class)
public void testNotSymmetricMatrixException() throws MathException {
double[][] changed = testData.clone();
changed[0][changed[0].length - 1] += 1.0e-5;
@ -66,7 +62,7 @@ public class CholeskyDecompositionImplTest {
}
/** test non positive definite matrix */
@Test(expected = NotPositiveDefiniteMatrixException.class)
@Test(expected = NonPositiveDefiniteMatrixException.class)
public void testNotPositiveDefinite() throws MathException {
new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(new double[][] {
{ 14, 11, 13, 15, 24 },
@ -77,7 +73,7 @@ public class CholeskyDecompositionImplTest {
}));
}
@Test(expected = NotPositiveDefiniteMatrixException.class)
@Test(expected = NonPositiveDefiniteMatrixException.class)
public void testMath274() throws MathException {
new CholeskyDecompositionImpl(MatrixUtils.createRealMatrix(new double[][] {
{ 0.40434286, -0.09376327, 0.30328980, 0.04909388 },

3
src/test/java/org/apache/commons/math/linear/EigenSolverTest.java
View File

@ -19,6 +19,7 @@ package org.apache.commons.math.linear;
import java.util.Random;
import org.apache.commons.math.exception.SingularMatrixException;
import junit.framework.TestCase;
import org.apache.commons.math.util.MathUtils;
@ -42,7 +43,7 @@ public class EigenSolverTest extends TestCase {
try {
es.getInverse();
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
} catch (SingularMatrixException ime) {
// expected behavior
}
}

3
src/test/java/org/apache/commons/math/linear/FieldLUDecompositionImplTest.java
View File

@ -22,6 +22,7 @@ import junit.framework.TestCase;
import org.apache.commons.math.TestUtils;
import org.apache.commons.math.fraction.Fraction;
import org.apache.commons.math.fraction.FractionField;
import org.apache.commons.math.exception.NonSquareMatrixException;
public class FieldLUDecompositionImplTest extends TestCase {
private Fraction[][] testData = {
@ -77,7 +78,7 @@ public class FieldLUDecompositionImplTest extends TestCase {
{ Fraction.ZERO, Fraction.ZERO },
{ Fraction.ZERO, Fraction.ZERO }
}));
} catch (InvalidMatrixException ime) {
} catch (NonSquareMatrixException ime) {
// expected behavior
}
}

1
src/test/java/org/apache/commons/math/linear/FieldMatrixImplTest.java
View File

@ -27,6 +27,7 @@ import org.apache.commons.math.exception.OutOfRangeException;
import org.apache.commons.math.exception.NumberIsTooSmallException;
import org.apache.commons.math.exception.NotStrictlyPositiveException;
import org.apache.commons.math.exception.NullArgumentException;
import org.apache.commons.math.exception.NonSquareMatrixException;
/**
* Test cases for the {@link Array2DRowFieldMatrix} class.

5
src/test/java/org/apache/commons/math/linear/LUDecompositionImplTest.java
View File

@ -17,6 +17,7 @@
package org.apache.commons.math.linear;
import org.apache.commons.math.exception.NonSquareMatrixException;
import junit.framework.TestCase;
public class LUDecompositionImplTest extends TestCase {
@ -73,8 +74,8 @@ public class LUDecompositionImplTest extends TestCase {
public void testNonSquare() {
try {
new LUDecompositionImpl(MatrixUtils.createRealMatrix(new double[3][2]));
fail("Expecting InvalidMatrixException");
} catch (InvalidMatrixException ime) {
fail("Expecting NonSquareMatrixException");
} catch (NonSquareMatrixException ime) {
// expected behavior
}
}

9
src/test/java/org/apache/commons/math/linear/LUSolverTest.java
View File

@ -17,6 +17,7 @@
package org.apache.commons.math.linear;
import org.apache.commons.math.exception.SingularMatrixException;
import junit.framework.TestCase;
public class LUSolverTest extends TestCase {
@ -102,25 +103,25 @@ public class LUSolverTest extends TestCase {
try {
solver.solve(b);
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
} catch (SingularMatrixException ime) {
// expected behavior
}
try {
solver.solve(b.getColumn(0));
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
} catch (SingularMatrixException ime) {
// expected behavior
}
try {
solver.solve(b.getColumnVector(0));
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
} catch (SingularMatrixException ime) {
// expected behavior
}
try {
solver.solve(new ArrayRealVectorTest.RealVectorTestImpl(b.getColumn(0)));
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
} catch (SingularMatrixException ime) {
// expected behavior
}
}

7
src/test/java/org/apache/commons/math/linear/QRSolverTest.java
View File

@ -19,6 +19,7 @@ package org.apache.commons.math.linear;
import java.util.Random;
import org.apache.commons.math.exception.SingularMatrixException;
import junit.framework.TestCase;
public class QRSolverTest extends TestCase {
@ -101,19 +102,19 @@ public class QRSolverTest extends TestCase {
try {
solver.solve(b);
fail("an exception should have been thrown");
} catch (InvalidMatrixException iae) {
} catch (SingularMatrixException iae) {
// expected behavior
}
try {
solver.solve(b.getColumn(0));
fail("an exception should have been thrown");
} catch (InvalidMatrixException iae) {
} catch (SingularMatrixException iae) {
// expected behavior
}
try {
solver.solve(b.getColumnVector(0));
fail("an exception should have been thrown");
} catch (InvalidMatrixException iae) {
} catch (SingularMatrixException iae) {
// expected behavior
}
}

1
src/test/java/org/apache/commons/math/linear/SparseFieldMatrixTest.java
View File

@ -26,6 +26,7 @@ import org.apache.commons.math.exception.NoDataException;
import org.apache.commons.math.exception.OutOfRangeException;
import org.apache.commons.math.exception.NumberIsTooSmallException;
import org.apache.commons.math.exception.NullArgumentException;
import org.apache.commons.math.exception.NonSquareMatrixException;
/**
* Test cases for the {@link SparseFieldMatrix} class.

1
src/test/java/org/apache/commons/math/linear/SparseRealMatrixTest.java
View File

@ -22,6 +22,7 @@ import org.apache.commons.math.TestUtils;
import org.apache.commons.math.exception.OutOfRangeException;
import org.apache.commons.math.exception.ZeroException;
import org.apache.commons.math.exception.NumberIsTooSmallException;
import org.apache.commons.math.exception.NonSquareMatrixException;
/**
* Test cases for the {@link OpenMapRealMatrix} class.

3
src/test/java/org/apache/commons/math/linear/TriDiagonalTransformerTest.java
View File

@ -20,6 +20,7 @@ package org.apache.commons.math.linear;
import java.util.Arrays;
import org.apache.commons.math.util.FastMath;
import org.apache.commons.math.exception.NonSquareMatrixException;
import junit.framework.TestCase;
@ -47,7 +48,7 @@ public class TriDiagonalTransformerTest extends TestCase {
try {
new TriDiagonalTransformer(MatrixUtils.createRealMatrix(new double[3][2]));
fail("an exception should have been thrown");
} catch (InvalidMatrixException ime) {
} catch (NonSquareMatrixException ime) {
// expected behavior
}
}

4
src/test/java/org/apache/commons/math/optimization/general/GaussNewtonOptimizerTest.java
View File

@ -24,6 +24,7 @@ import java.util.Arrays;
import junit.framework.TestCase;
import org.apache.commons.math.exception.SingularMatrixException;
import org.apache.commons.math.exception.FunctionEvaluationException;
import org.apache.commons.math.exception.ConvergenceException;
import org.apache.commons.math.exception.TooManyEvaluationsException;
@ -278,7 +279,6 @@ extends TestCase {
} catch (ConvergenceException ee) {
// expected behavior
}
}
public void testMoreEstimatedParametersUnsorted() throws Exception {
@ -362,7 +362,7 @@ extends TestCase {
new double[] { 1 },
new double[] { 0, 0 });
fail("an exception should have been thrown");
} catch (FunctionEvaluationException oe) {
} catch (DimensionMismatchException oe) {
// expected behavior
}

40
src/test/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizerTest.java
View File

@ -25,10 +25,12 @@ import java.util.List;
import junit.framework.TestCase;
import org.apache.commons.math.exception.SingularMatrixException;
import org.apache.commons.math.exception.FunctionEvaluationException;
import org.apache.commons.math.exception.ConvergenceException;
import org.apache.commons.math.exception.DimensionMismatchException;
import org.apache.commons.math.exception.TooManyEvaluationsException;
import org.apache.commons.math.exception.NumberIsTooSmallException;
import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction;
import org.apache.commons.math.analysis.MultivariateMatrixFunction;
import org.apache.commons.math.linear.BlockRealMatrix;
@ -107,7 +109,7 @@ public class LevenbergMarquardtOptimizerTest
super(name);
}
public void testTrivial() throws FunctionEvaluationException {
public void testTrivial() {
LinearProblem problem =
new LinearProblem(new double[][] { { 2 } }, new double[] { 3 });
LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
@ -117,14 +119,14 @@ public class LevenbergMarquardtOptimizerTest
try {
optimizer.guessParametersErrors();
fail("an exception should have been thrown");
} catch (ConvergenceException ee) {
} catch (NumberIsTooSmallException ee) {
// expected behavior
}
assertEquals(1.5, optimum.getPoint()[0], 1.0e-10);
assertEquals(3.0, optimum.getValue()[0], 1.0e-10);
}
public void testQRColumnsPermutation() throws FunctionEvaluationException {
public void testQRColumnsPermutation() {
LinearProblem problem =
new LinearProblem(new double[][] { { 1.0, -1.0 }, { 0.0, 2.0 }, { 1.0, -2.0 } },
@ -141,7 +143,7 @@ public class LevenbergMarquardtOptimizerTest
assertEquals(1.0, optimum.getValue()[2], 1.0e-10);
}
public void testNoDependency() throws FunctionEvaluationException {
public void testNoDependency() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 2, 0, 0, 0, 0, 0 },
{ 0, 2, 0, 0, 0, 0 },
@ -160,7 +162,7 @@ public class LevenbergMarquardtOptimizerTest
}
}
public void testOneSet() throws FunctionEvaluationException {
public void testOneSet() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 1, 0, 0 },
@ -176,7 +178,7 @@ public class LevenbergMarquardtOptimizerTest
assertEquals(3.0, optimum.getPoint()[2], 1.0e-10);
}
public void testTwoSets() throws FunctionEvaluationException {
public void testTwoSets() {
double epsilon = 1.0e-7;
LinearProblem problem = new LinearProblem(new double[][] {
{ 2, 1, 0, 4, 0, 0 },
@ -200,7 +202,7 @@ public class LevenbergMarquardtOptimizerTest
assertEquals( 1.0 - epsilon, optimum.getPoint()[5], 1.0e-10);
}
public void testNonInversible() throws FunctionEvaluationException {
public void testNonInvertible() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 1, 2, -3 },
@ -214,12 +216,12 @@ public class LevenbergMarquardtOptimizerTest
try {
optimizer.getCovariances();
fail("an exception should have been thrown");
} catch (ConvergenceException ee) {
} catch (SingularMatrixException ee) {
// expected behavior
}
}
public void testIllConditioned() throws FunctionEvaluationException {
public void testIllConditioned() {
LinearProblem problem1 = new LinearProblem(new double[][] {
{ 10.0, 7.0, 8.0, 7.0 },
{ 7.0, 5.0, 6.0, 5.0 },
@ -252,7 +254,7 @@ public class LevenbergMarquardtOptimizerTest
assertEquals( 22.0, optimum2.getPoint()[3], 1.0e-8);
}
public void testMoreEstimatedParametersSimple() throws FunctionEvaluationException {
public void testMoreEstimatedParametersSimple() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 3.0, 2.0, 0.0, 0.0 },
@ -266,7 +268,7 @@ public class LevenbergMarquardtOptimizerTest
assertEquals(0, optimizer.getRMS(), 1.0e-10);
}
public void testMoreEstimatedParametersUnsorted() throws FunctionEvaluationException {
public void testMoreEstimatedParametersUnsorted() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 1.0, 1.0, 0.0, 0.0, 0.0, 0.0 },
{ 0.0, 0.0, 1.0, 1.0, 1.0, 0.0 },
@ -286,7 +288,7 @@ public class LevenbergMarquardtOptimizerTest
assertEquals(6.0, optimum.getPointRef()[5], 1.0e-10);
}
public void testRedundantEquations() throws FunctionEvaluationException {
public void testRedundantEquations() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 1.0, 1.0 },
{ 1.0, -1.0 },
@ -302,7 +304,7 @@ public class LevenbergMarquardtOptimizerTest
assertEquals(1.0, optimum.getPointRef()[1], 1.0e-10);
}
public void testInconsistentEquations() throws FunctionEvaluationException {
public void testInconsistentEquations() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 1.0, 1.0 },
{ 1.0, -1.0 },
@ -314,7 +316,7 @@ public class LevenbergMarquardtOptimizerTest
assertTrue(optimizer.getRMS() > 0.1);
}
public void testInconsistentSizes() throws FunctionEvaluationException {
public void testInconsistentSizes() {
LinearProblem problem =
new LinearProblem(new double[][] { { 1, 0 }, { 0, 1 } }, new double[] { -1, 1 });
LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
@ -339,7 +341,7 @@ public class LevenbergMarquardtOptimizerTest
new double[] { 1 },
new double[] { 0, 0 });
fail("an exception should have been thrown");
} catch (FunctionEvaluationException oe) {
} catch (DimensionMismatchException oe) {
// expected behavior
}
}
@ -373,14 +375,14 @@ public class LevenbergMarquardtOptimizerTest
optimizer.optimize(problem, new double[] { 0, 0, 0, 0, 0 }, new double[] { 1, 1, 1, 1, 1 },
new double[] { 98.680, 47.345 });
assertTrue(!shouldFail);
} catch (FunctionEvaluationException ee) {
} catch (DimensionMismatchException ee) {
assertTrue(shouldFail);
} catch (TooManyEvaluationsException ee) {
assertTrue(shouldFail);
}
}
public void testCircleFitting() throws FunctionEvaluationException {
public void testCircleFitting() {
Circle circle = new Circle();
circle.addPoint( 30.0, 68.0);
circle.addPoint( 50.0, -6.0);
@ -428,7 +430,7 @@ public class LevenbergMarquardtOptimizerTest
assertEquals(0.004, errors[1], 0.001);
}
public void testCircleFittingBadInit() throws FunctionEvaluationException {
public void testCircleFittingBadInit() {
Circle circle = new Circle();
double[][] points = new double[][] {
{-0.312967, 0.072366}, {-0.339248, 0.132965}, {-0.379780, 0.202724},
@ -481,7 +483,7 @@ public class LevenbergMarquardtOptimizerTest
assertEquals( 0.2075001, center.y, 1.0e-6);
}
public void testMath199() throws FunctionEvaluationException {
public void testMath199() {
try {
QuadraticProblem problem = new QuadraticProblem();
problem.addPoint (0, -3.182591015485607);