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

MATH-439 

Refactored the "solvers" package. Implementations refer to number of
evaluation of the objective function (instead of the number of iterations).
New interfaces and base classes.
"NewtonSolver" fits in the design without resorting to a cast.
Created class "MullerSolver2" to contain the code of the method named "solve2"
in class "MullerSolver".
Removed "UnivariateRealSolverFactory" and "UnivariateRealSolverFactoryImpl".
Default solver in "UnivariateRealSolverUtils" is explicitely instantiated.
"AbstractContinuousDistribution": Type of exception thrown changed in
"UnivariateRealSolverUtils".
Factored out duplicate code (in "GaussNewtonOptimizerTest" and
"LevenbergMarquardtOptimizerTest"): class "Circle" is now called
"CircleVectorial". Also factored out the "Circle" class from
"NonLinearConjugateGradientOptimizerTest": class is named "CircleScalar".
Created "SecantSolverTest", moving there all the tests for the class
"SecantSolver" that were located in class "BrentSolverTest".
Created new interface and base class for polynomial functions solvers
("LaguerreSolver") so that the function type is now checked at compile time.
Removed deprecated exceptions (MATH-441).
Javadoc clean-up.
Lowered tolerance values in some unit tests.
Tests upgraded to Junit 4 (MATH-423).



git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1039083 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Gilles Sadowski 2010-11-25 16:22:00 +00:00
parent c32b8042d8
commit b9559bfce9
50 changed files with 2810 additions and 2207 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
src/main/java/org/apache/commons/math/analysis/FunctionUtils.java
View File

@ -26,6 +26,11 @@ import org.apache.commons.math.analysis.function.Identity;
* @since 3.0
*/
public class FunctionUtils {
/**
* Class only contains static methods.
*/
private FunctionUtils() {}
/**
* Compose functions.
*
@ -142,6 +147,7 @@ public class FunctionUtils {
*
* @param f Binary function.
* @param fixed Value to which the first argument of {@code f} is set.
* @return a unary function.
*/
public static UnivariateRealFunction fix1stArgument(final BivariateRealFunction f,
final double fixed) {
@ -157,6 +163,7 @@ public class FunctionUtils {
*
* @param f Binary function.
* @param fixed Value to which the second argument of {@code f} is set.
* @return a unary function.
*/
public static UnivariateRealFunction fix2ndArgument(final BivariateRealFunction f,
final double fixed) {

10
src/main/java/org/apache/commons/math/analysis/polynomials/PolynomialFunctionNewtonForm.java
View File

@ -16,7 +16,6 @@
*/
package org.apache.commons.math.analysis.polynomials;
import org.apache.commons.math.exception.NullArgumentException;
import org.apache.commons.math.exception.NoDataException;
import org.apache.commons.math.exception.DimensionMismatchException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
@ -69,7 +68,8 @@ public class PolynomialFunctionNewtonForm implements UnivariateRealFunction {
*
* @param a Coefficients in Newton form formula.
* @param c Centers.
* @throws NullArgumentException if any argument is {@code null}.
* @throws org.apache.commons.math.exception.NullArgumentException if
* any argument is {@code null}.
* @throws NoDataException if any array has zero length.
* @throws DimensionMismatchException if the size difference between
* {@code a} and {@code c} is not equal to 1.
@ -154,7 +154,8 @@ public class PolynomialFunctionNewtonForm implements UnivariateRealFunction {
* @param c Centers.
* @param z Point at which the function value is to be computed.
* @return the function value.
* @throws NullArgumentException if any argument is {@code null}.
* @throws org.apache.commons.math.exception.NullArgumentException if
* any argument is {@code null}.
* @throws NoDataException if any array has zero length.
* @throws DimensionMismatchException if the size difference between
* {@code a} and {@code c} is not equal to 1.
@ -202,7 +203,8 @@ public class PolynomialFunctionNewtonForm implements UnivariateRealFunction {
*
* @param a the coefficients in Newton form formula
* @param c the centers
* @throws NullArgumentException if any argument is {@code null}.
* @throws org.apache.commons.math.exception.NullArgumentException if
* any argument is {@code null}.
* @throws NoDataException if any array has zero length.
* @throws DimensionMismatchException if the size difference between
* {@code a} and {@code c} is not equal to 1.

81
src/main/java/org/apache/commons/math/analysis/solvers/AbstractDifferentiableUnivariateRealSolver.java Normal file
View File

@ -0,0 +1,81 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.analysis.DifferentiableUnivariateRealFunction;
import org.apache.commons.math.analysis.UnivariateRealFunction;
/**
* Provide a default implementation for several functions useful to generic
* solvers.
*
* @version $Revision$ $Date$
* @since 3.0
*/
public abstract class AbstractDifferentiableUnivariateRealSolver
extends BaseAbstractUnivariateRealSolver<DifferentiableUnivariateRealFunction>
implements DifferentiableUnivariateRealSolver {
/** Derivative of the function to solve. */
private UnivariateRealFunction functionDerivative;
/**
* Construct a solver with given absolute accuracy.
*
* @param absoluteAccuracy Maximum absolute error.
*/
protected AbstractDifferentiableUnivariateRealSolver(final double absoluteAccuracy) {
super(absoluteAccuracy);
}
/**
* Construct a solver with given accuracies.
*
* @param relativeAccuracy Maximum relative error.
* @param absoluteAccuracy Maximum absolute error.
* @param functionValueAccuracy Maximum function value error.
*/
protected AbstractDifferentiableUnivariateRealSolver(final double relativeAccuracy,
final double absoluteAccuracy,
final double functionValueAccuracy) {
super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
}
/**
* Compute the objective function value.
*
* @param point Point at which the objective function must be evaluated.
* @return the objective function value at specified point.
* @throws org.apache.commons.math.exception.TooManyEvaluationsException
* if the maximal number of evaluations is exceeded.
*/
protected double computeDerivativeObjectiveValue(double point) {
incrementEvaluationCount();
return functionDerivative.value(point);
}
/**
* {@inheritDoc}
*/
@Override
protected void setup(DifferentiableUnivariateRealFunction f,
double min, double max,
double startValue) {
super.setup(f, min, max, startValue);
functionDerivative = f.derivative();
}
}

82
src/main/java/org/apache/commons/math/analysis/solvers/AbstractPolynomialSolver.java Normal file
View File

@ -0,0 +1,82 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
/**
* Base class for solvers.
*
* @version $Revision$ $Date$
* @since 3.0
*/
public abstract class AbstractPolynomialSolver
extends BaseAbstractUnivariateRealSolver<PolynomialFunction>
implements PolynomialSolver {
/** Function. */
private PolynomialFunction polynomialFunction;
/**
* Construct a solver with given absolute accuracy.
*
* @param absoluteAccuracy Maximum absolute error.
*/
protected AbstractPolynomialSolver(final double absoluteAccuracy) {
super(absoluteAccuracy);
}
/**
* Construct a solver with given accuracies.
*
* @param relativeAccuracy Maximum relative error.
* @param absoluteAccuracy Maximum absolute error.
*/
protected AbstractPolynomialSolver(final double relativeAccuracy,
final double absoluteAccuracy) {
super(relativeAccuracy, absoluteAccuracy);
}
/**
* Construct a solver with given accuracies.
*
* @param relativeAccuracy Maximum relative error.
* @param absoluteAccuracy Maximum absolute error.
* @param functionValueAccuracy Maximum function value error.
*/
protected AbstractPolynomialSolver(final double relativeAccuracy,
final double absoluteAccuracy,
final double functionValueAccuracy) {
super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
}
/**
* {@inheritDoc}
*/
@Override
protected void setup(PolynomialFunction f,
double min, double max,
double startValue) {
super.setup(f, min, max, startValue);
polynomialFunction = f;
}
/**
* @return the coefficients of the polynomial function.
*/
protected double[] getCoefficients() {
return polynomialFunction.getCoefficients();
}
}

61
src/main/java/org/apache/commons/math/analysis/solvers/AbstractUnivariateRealSolver.java Normal file
View File

@ -0,0 +1,61 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.analysis.UnivariateRealFunction;
/**
* Base class for solvers.
*
* @version $Revision$ $Date$
* @since 3.0
*/
public abstract class AbstractUnivariateRealSolver
extends BaseAbstractUnivariateRealSolver<UnivariateRealFunction>
implements UnivariateRealSolver {
/**
* Construct a solver with given absolute accuracy.
*
* @param absoluteAccuracy Maximum absolute error.
*/
protected AbstractUnivariateRealSolver(final double absoluteAccuracy) {
super(absoluteAccuracy);
}
/**
* Construct a solver with given accuracies.
*
* @param relativeAccuracy Maximum relative error.
* @param absoluteAccuracy Maximum absolute error.
*/
protected AbstractUnivariateRealSolver(final double relativeAccuracy,
final double absoluteAccuracy) {
super(relativeAccuracy, absoluteAccuracy);
}
/**
* Construct a solver with given accuracies.
*
* @param relativeAccuracy Maximum relative error.
* @param absoluteAccuracy Maximum absolute error.
* @param functionValueAccuracy Maximum function value error.
*/
protected AbstractUnivariateRealSolver(final double relativeAccuracy,
final double absoluteAccuracy,
final double functionValueAccuracy) {
super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
}
}

299
src/main/java/org/apache/commons/math/analysis/solvers/BaseAbstractUnivariateRealSolver.java Normal file
View File

@ -0,0 +1,299 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.util.Incrementor;
import org.apache.commons.math.exception.MaxCountExceededException;
import org.apache.commons.math.exception.TooManyEvaluationsException;
import org.apache.commons.math.exception.NullArgumentException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
/**
* Provide a default implementation for several functions useful to generic
* solvers.
*
* @param <FUNC> Type of function to solve.
*
* @version $Revision: 1030464 $ $Date: 2010-11-03 14:46:04 +0100 (Wed, 03 Nov 2010) $
* @since 2.0
*/
public abstract class BaseAbstractUnivariateRealSolver<FUNC extends UnivariateRealFunction>
implements BaseUnivariateRealSolver<FUNC> {
/** Default absolute accuracy */
public static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
/** Default relative accuracy. */
public static final double DEFAULT_RELATIVE_ACCURACY = 1e-14;
/** Default function value accuracy. */
public static final double DEFAULT_FUNCTION_VALUE_ACCURACY = 1e-15;
/** Function value accuracy. */
private final double functionValueAccuracy;
/** Absolute accuracy. */
private final double absoluteAccuracy;
/** Relative accuracy. */
private final double relativeAccuracy;
/** Evaluations counter. */
private final Incrementor evaluations = new Incrementor();
/** Lower end of search interval. */
private double searchMin;
/** Higher end of search interval. */
private double searchMax;
/** Initial guess. */
private double searchStart;
/** Function to solve. */
private FUNC function;
/**
* Construct a solver with given absolute accuracy.
*
* @param absoluteAccuracy Maximum absolute error.
*/
protected BaseAbstractUnivariateRealSolver(final double absoluteAccuracy) {
this(DEFAULT_RELATIVE_ACCURACY,
absoluteAccuracy,
DEFAULT_FUNCTION_VALUE_ACCURACY);
}
/**
* Construct a solver with given accuracies.
*
* @param relativeAccuracy Maximum relative error.
* @param absoluteAccuracy Maximum absolute error.
*/
protected BaseAbstractUnivariateRealSolver(final double relativeAccuracy,
final double absoluteAccuracy) {
this(relativeAccuracy,
absoluteAccuracy,
DEFAULT_FUNCTION_VALUE_ACCURACY);
}
/**
* Construct a solver with given accuracies.
*
* @param relativeAccuracy Maximum relative error.
* @param absoluteAccuracy Maximum absolute error.
* @param functionValueAccuracy Maximum function value error.
*/
protected BaseAbstractUnivariateRealSolver(final double relativeAccuracy,
final double absoluteAccuracy,
final double functionValueAccuracy) {
this.absoluteAccuracy = absoluteAccuracy;
this.relativeAccuracy = relativeAccuracy;
this.functionValueAccuracy = functionValueAccuracy;
}
/** {@inheritDoc} */
public void setMaxEvaluations(int maxEvaluations) {
evaluations.setMaximalCount(maxEvaluations);
}
/** {@inheritDoc} */
public int getMaxEvaluations() {
return evaluations.getMaximalCount();
}
/** {@inheritDoc} */
public int getEvaluations() {
return evaluations.getCount();
}
/**
* @return the lower end of the search interval.
*/
public double getMin() {
return searchMin;
}
/**
* @return the higher end of the search interval.
*/
public double getMax() {
return searchMax;
}
/**
* @return the initial guess.
*/
public double getStartValue() {
return searchStart;
}
/**
* {@inheritDoc}
*/
public double getAbsoluteAccuracy() {
return absoluteAccuracy;
}
/**
* {@inheritDoc}
*/
public double getRelativeAccuracy() {
return relativeAccuracy;
}
/**
* {@inheritDoc}
*/
public double getFunctionValueAccuracy() {
return functionValueAccuracy;
}
/**
* Compute the objective function value.
*
* @param point Point at which the objective function must be evaluated.
* @return the objective function value at specified point.
* @throws TooManyEvaluationsException if the maximal number of evaluations
* is exceeded.
*/
protected double computeObjectiveValue(double point) {
incrementEvaluationCount();
return function.value(point);
}
/**
* Prepare for computation.
* Subclasses must call this method if they override any of the
* {@code solve} methods.
*
* @param f Function to solve.
* @param min Lower bound for the interval.
* @param max Upper bound for the interval.
* @param startValue Start value to use.
*/
protected void setup(FUNC f,
double min, double max,
double startValue) {
// Checks.
if (f == null) {
throw new NullArgumentException();
}
// Reset.
searchMin = min;
searchMax = max;
searchStart = startValue;
function = f;
evaluations.resetCount();
}
/** {@inheritDoc} */
public double solve(FUNC f, double min, double max, double startValue) {
// Initialization.
setup(f, min, max, startValue);
// Perform computation.
return doSolve();
}
/** {@inheritDoc} */
public double solve(FUNC f, double min, double max) {
return solve(f, min, max, min + 0.5 * (max - min));
}
/** {@inheritDoc} */
public double solve(FUNC f, double startValue) {
return solve(f, Double.NaN, Double.NaN, startValue);
}
/**
* Method for implementing actual optimization algorithms in derived
* classes.
*
* @return the root.
* @throws TooManyEvaluationsException if the maximal number of evaluations
* is exceeded.
*/
protected abstract double doSolve();
/**
* Check whether the function takes opposite signs at the endpoints.
*
* @param lower Lower endpoint.
* @param upper Upper endpoint.
* @return {@code true} if the function values have opposite signs at the
* given points.
*/
protected boolean isBracketing(final double lower,
final double upper) {
return UnivariateRealSolverUtils.isBracketing(function, lower, upper);
}
/**
* Check whether the arguments form a (strictly) increasing sequence.
*
* @param start First number.
* @param mid Second number.
* @param end Third number.
* @return {@code true} if the arguments form an increasing sequence.
*/
protected boolean isSequence(final double start,
final double mid,
final double end) {
return UnivariateRealSolverUtils.isSequence(start, mid, end);
}
/**
* Check that the endpoints specify an interval.
*
* @param lower Lower endpoint.
* @param upper Upper endpoint.
* @throws org.apache.commons.math.exception.NumberIsTooLargeException
* if {@code lower >= upper}.
*/
protected void verifyInterval(final double lower,
final double upper) {
UnivariateRealSolverUtils.verifyInterval(lower, upper);
}
/**
* Check that {@code lower < initial < upper}.
*
* @param lower Lower endpoint.
* @param initial Initial value.
* @param upper Upper endpoint.
* @throws org.apache.commons.math.exception.NumberIsTooLargeException
* if {@code lower >= initial} or {@code initial >= upper}.
*/
protected void verifySequence(final double lower,
final double initial,
final double upper) {
UnivariateRealSolverUtils.verifySequence(lower, initial, upper);
}
/**
* Check that the endpoints specify an interval and the function takes
* opposite signs at the endpoints.
*
* @param lower Lower endpoint.
* @param upper Upper endpoint.
* @throws org.apache.commons.math.exception.NoBracketingException if
* the function has the same sign at the endpoints.
*/
protected void verifyBracketing(final double lower,
final double upper) {
UnivariateRealSolverUtils.verifyBracketing(function, lower, upper);
}
/**
* Increment the evaluation count by one.
* Method {@link #computeObjectiveValue(double)} calls this method internally.
* It is provided for subclasses that do not exclusively use
* {@code computeObjectiveValue} to solve the function.
* See e.g. {@link AbstractDifferentiableUnivariateRealSolver}.
*/
protected void incrementEvaluationCount() {
try {
evaluations.incrementCount();
} catch (MaxCountExceededException e) {
throw new TooManyEvaluationsException(e.getMax());
}
}
}

114
src/main/java/org/apache/commons/math/analysis/solvers/BaseUnivariateRealSolver.java Normal file
View File

@ -0,0 +1,114 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.analysis.UnivariateRealFunction;
/**
* Interface for (univariate real) rootfinding algorithms.
* Implementations will search for only one zero in the given interval.
*
* @param <FUNC> Type of function to solve.
*
* @version $Revision$ $Date$
* @since 3.0
*/
public interface BaseUnivariateRealSolver<FUNC extends UnivariateRealFunction> {
/**
* Set the maximal number of function evaluations.
*
* @param maxEvaluations Maximal number of function evaluations.
*/
void setMaxEvaluations(int maxEvaluations);
/**
* Get the maximal number of function evaluations.
*
* @return the maximal number of function evaluations.
*/
int getMaxEvaluations();
/**
* Get the number of evaluations of the objective function.
* The number of evaluations corresponds to the last call to the
* {@code optimize} method. It is 0 if the method has not been
* called yet.
*
* @return the number of evaluations of the objective function.
*/
int getEvaluations();
/**
* @return the absolute accuracy.
*/
double getAbsoluteAccuracy();
/**
* @return the relative accuracy.
*/
double getRelativeAccuracy();
/**
* @return the function value accuracy.
*/
double getFunctionValueAccuracy();
/**
* Solve for a zero root in the given interval.
* A solver may require that the interval brackets a single zero root.
* Solvers that do require bracketing should be able to handle the case
* where one of the endpoints is itself a root.
*
* @param f Function to solve.
* @param min Lower bound for the interval.
* @param max Upper bound for the interval.
* @return a value where the function is zero.
* @throws IllegalArgumentException if {@code min > max} or the endpoints
* do not satisfy the requirements specified by the solver.
* @since 2.0
*/
double solve(FUNC f, double min, double max);
/**
* Solve for a zero in the given interval, start at {@code startValue}.
* A solver may require that the interval brackets a single zero root.
* Solvers that do require bracketing should be able to handle the case
* where one of the endpoints is itself a root.
*
* @param f Function to solve.
* @param min Lower bound for the interval.
* @param max Upper bound for the interval.
* @param startValue Start value to use.
* @return a value where the function is zero.
* @throws IllegalArgumentException if {@code min > max} or the arguments
* do not satisfy the requirements specified by the solver.
* @since 2.0
*/
double solve(FUNC f, double min, double max, double startValue);
/**
* Solve for a zero in the vicinity of {@code startValue}.
* A solver may require that the interval brackets a single zero root.
*
* @param f Function to solve.
* @param startValue Start value to use.
* @return a value where the function is zero.
* @throws IllegalArgumentException if {@code min > max} or the arguments
* do not satisfy the requirements specified by the solver.
* @since 2.0
*/
double solve(FUNC f, double startValue);
}

62
src/main/java/org/apache/commons/math/analysis/solvers/BisectionSolver.java
View File

@ -16,9 +16,6 @@
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.util.FastMath;
/**
@ -29,39 +26,54 @@ import org.apache.commons.math.util.FastMath;
*
* @version $Revision$ $Date$
*/
public class BisectionSolver extends UnivariateRealSolverImpl {
public class BisectionSolver extends AbstractUnivariateRealSolver {
/** Default absolute accuracy. */
public static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
/**
* Construct a solver with default accuracy.
*/
public BisectionSolver() {
this(DEFAULT_ABSOLUTE_ACCURACY);
}
/**
* Construct a solver.
*
* @param absoluteAccuracy Absolute accuracy.
*/
public BisectionSolver() {
super(100, 1E-6);
public BisectionSolver(double absoluteAccuracy) {
super(absoluteAccuracy);
}
/**
* Construct a solver.
*
* @param relativeAccuracy Relative accuracy.
* @param absoluteAccuracy Absolute accuracy.
*/
public BisectionSolver(double relativeAccuracy,
double absoluteAccuracy) {
super(relativeAccuracy, absoluteAccuracy);
}
/** {@inheritDoc} */
public double solve(final UnivariateRealFunction f, double min, double max, double initial)
throws MaxIterationsExceededException, MathUserException {
return solve(f, min, max);
}
/** {@inheritDoc} */
public double solve(final UnivariateRealFunction f, double min, double max)
throws MaxIterationsExceededException, MathUserException {
clearResult();
verifyInterval(min,max);
/**
* {@inheritDoc}
*/
@Override
protected double doSolve() {
double min = getMin();
double max = getMax();
verifyInterval(min, max);
final double absoluteAccuracy = getAbsoluteAccuracy();
double m;
double fm;
double fmin;
int i = 0;
while (i < maximalIterationCount) {
while (true) {
m = UnivariateRealSolverUtils.midpoint(min, max);
fmin = f.value(min);
fm = f.value(m);
fmin = computeObjectiveValue(min);
fm = computeObjectiveValue(m);
if (fm * fmin > 0.0) {
if (fm * fmin > 0) {
// max and m bracket the root.
min = m;
} else {
@ -71,12 +83,8 @@ public class BisectionSolver extends UnivariateRealSolverImpl {
if (FastMath.abs(max - min) <= absoluteAccuracy) {
m = UnivariateRealSolverUtils.midpoint(min, max);
setResult(m, i);
return m;
}
++i;
}
throw new MaxIterationsExceededException(maximalIterationCount);
}
}

365
src/main/java/org/apache/commons/math/analysis/solvers/BrentSolver.java
View File

@ -17,297 +17,214 @@
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.exception.NoBracketingException;
import org.apache.commons.math.util.FastMath;
import org.apache.commons.math.util.MathUtils;
/**
* Implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
* Brent algorithm</a> for finding zeros of real univariate functions.
* <p>
* The function should be continuous but not necessarily smooth.</p>
* This class implements the <a href="http://mathworld.wolfram.com/BrentsMethod.html">
* Brent algorithm</a> for finding zeros of real univariate functions.
* The function should be continuous but not necessarily smooth.
* The {@code solve} method returns a zero {@code x} of the function {@code f}
* in the given interval {@code [a, b]} to within a tolerance
* {@code 6 eps abs(x) + t} where {@code eps} is the relative accuracy and
* {@code t} is the absolute accuracy.
* The given interval must bracket the root.
*
* @version $Revision:670469 $ $Date:2008-06-23 10:01:38 +0200 (lun., 23 juin 2008) $
*/
public class BrentSolver extends UnivariateRealSolverImpl {
/**
* Default absolute accuracy
* @since 2.1
*/
public static final double DEFAULT_ABSOLUTE_ACCURACY = 1E-6;
/** Default maximum number of iterations
* @since 2.1
*/
public static final int DEFAULT_MAXIMUM_ITERATIONS = 100;
public class BrentSolver extends AbstractUnivariateRealSolver {
/** Serializable version identifier */
private static final long serialVersionUID = 7694577816772532779L;
/** Default absolute accuracy. */
public static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
/**
* Construct a solver with default properties.
* Construct a solver with default accuracies.
*/
public BrentSolver() {
super(DEFAULT_MAXIMUM_ITERATIONS, DEFAULT_ABSOLUTE_ACCURACY);
this(DEFAULT_ABSOLUTE_ACCURACY);
}
/**
* Construct a solver with the given absolute accuracy.
* Construct a solver.
*
* @param absoluteAccuracy lower bound for absolute accuracy of solutions returned by the solver
* @since 2.1
* @param absoluteAccuracy Absolute accuracy.
*/
public BrentSolver(double absoluteAccuracy) {
super(DEFAULT_MAXIMUM_ITERATIONS, absoluteAccuracy);
super(absoluteAccuracy);
}
/**
* Construct a solver.
*
* @param relativeAccuracy Relative accuracy.
* @param absoluteAccuracy Absolute accuracy.
*/
public BrentSolver(double relativeAccuracy,
double absoluteAccuracy) {
super(relativeAccuracy, absoluteAccuracy);
}
/**
* Construct a solver.
*
* @param relativeAccuracy Relative accuracy.
* @param absoluteAccuracy Absolute accuracy.
* @param functionValueAccuracy Function value accuracy.
*/
public BrentSolver(double relativeAccuracy,
double absoluteAccuracy,
double functionValueAccuracy) {
super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
}
/**
* Contstruct a solver with the given maximum iterations and absolute accuracy.
*
* @param maximumIterations maximum number of iterations
* @param absoluteAccuracy lower bound for absolute accuracy of solutions returned by the solver
* @since 2.1
* {@inheritDoc}
*/
public BrentSolver(int maximumIterations, double absoluteAccuracy) {
super(maximumIterations, absoluteAccuracy);
}
@Override
protected double doSolve() {
double min = getMin();
double max = getMax();
final double initial = getStartValue();
final double functionValueAccuracy = getFunctionValueAccuracy();
/**
* Find a zero in the given interval with an initial guess.
* <p>Throws <code>IllegalArgumentException</code> if the values of the
* function at the three points have the same sign (note that it is
* allowed to have endpoints with the same sign if the initial point has
* opposite sign function-wise).</p>
*
* @param f function to solve.
* @param min the lower bound for the interval.
* @param max the upper bound for the interval.
* @param initial the start value to use (must be set to min if no
* initial point is known).
* @return the value where the function is zero
* @throws MaxIterationsExceededException the maximum iteration count is exceeded
* @throws MathUserException if an error occurs evaluating the function
* @throws IllegalArgumentException if initial is not between min and max
* (even if it <em>is</em> a root)
*/
public double solve(final UnivariateRealFunction f,
final double min, final double max, final double initial)
throws MaxIterationsExceededException, MathUserException {
verifySequence(min, initial, max);
clearResult();
if ((initial < min) || (initial > max)) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.INVALID_INTERVAL_INITIAL_VALUE_PARAMETERS,
min, initial, max);
}
// return the initial guess if it is good enough
double yInitial = f.value(initial);
// Return the initial guess if it is good enough.
double yInitial = computeObjectiveValue(initial);
if (FastMath.abs(yInitial) <= functionValueAccuracy) {
setResult(initial, 0);
return result;
return initial;
}
// return the first endpoint if it is good enough
double yMin = f.value(min);
// Return the first endpoint if it is good enough.
double yMin = computeObjectiveValue(min);
if (FastMath.abs(yMin) <= functionValueAccuracy) {
setResult(min, 0);
return result;
return min;
}
// reduce interval if min and initial bracket the root
// Reduce interval if min and initial bracket the root.
if (yInitial * yMin < 0) {
return solve(f, min, yMin, initial, yInitial, min, yMin);
return brent(min, initial, yMin, yInitial);
}
// return the second endpoint if it is good enough
double yMax = f.value(max);
// Return the second endpoint if it is good enough.
double yMax = computeObjectiveValue(max);
if (FastMath.abs(yMax) <= functionValueAccuracy) {
setResult(max, 0);
return result;
return max;
}
// reduce interval if initial and max bracket the root
// Reduce interval if initial and max bracket the root.
if (yInitial * yMax < 0) {
return solve(f, initial, yInitial, max, yMax, initial, yInitial);
return brent(initial, max, yInitial, yMax);
}
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.SAME_SIGN_AT_ENDPOINTS, min, max, yMin, yMax);
throw new NoBracketingException(min, max, yMin, yMax);
}
/**
* Find a zero in the given interval.
* <p>
* Requires that the values of the function at the endpoints have opposite
* signs. An <code>IllegalArgumentException</code> is thrown if this is not
* the case.</p>
* Search for a zero inside the provided interval.
* This implemenation is based on the algorithm described at page 58 of
* the book
* <quote>
* <b>Algorithms for Minimization Without Derivatives</b>
* <it>Richard P. Brent</it>
* Dover 0-486-41998-3
* </quote>
*
* @param f the function to solve
* @param min the lower bound for the interval.
* @param max the upper bound for the interval.
* @return the value where the function is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* @throws MathUserException if an error occurs evaluating the function
* @throws IllegalArgumentException if min is not less than max or the
* signs of the values of the function at the endpoints are not opposites
* @param lo Lower bound of the search interval.
* @param hi Higher bound of the search interval.
* @param fLo Function value at the lower bound of the search interval.
* @param fHi Function value at the higher bound of the search interval.
* @return the value where the function is zero.
*/
public double solve(final UnivariateRealFunction f,
final double min, final double max)
throws MaxIterationsExceededException, MathUserException {
private double brent(double lo, double hi,
double fLo, double fHi) {
double a = lo;
double fa = fLo;
double b = hi;
double fb = fHi;
double c = a;
double fc = fa;
double d = b - a;
double e = d;
clearResult();
verifyInterval(min, max);
final double t = getAbsoluteAccuracy();
final double eps = getRelativeAccuracy();
double ret = Double.NaN;
double yMin = f.value(min);
double yMax = f.value(max);
// Verify bracketing
double sign = yMin * yMax;
if (sign > 0) {
// check if either value is close to a zero
if (FastMath.abs(yMin) <= functionValueAccuracy) {
setResult(min, 0);
ret = min;
} else if (FastMath.abs(yMax) <= functionValueAccuracy) {
setResult(max, 0);
ret = max;
} else {
// neither value is close to zero and min and max do not bracket root.
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.SAME_SIGN_AT_ENDPOINTS, min, max, yMin, yMax);
while (true) {
if (FastMath.abs(fc) < FastMath.abs(fb)) {
a = b;
b = c;
c = a;
fa = fb;
fb = fc;
fc = fa;
}
} else if (sign < 0){
// solve using only the first endpoint as initial guess
ret = solve(f, min, yMin, max, yMax, min, yMin);
} else {
// either min or max is a root
if (yMin == 0.0) {
ret = min;
} else {
ret = max;
}
}
return ret;
}
final double tol = 2 * eps * FastMath.abs(b) + t;
final double m = 0.5 * (c - b);
/**
* Find a zero starting search according to the three provided points.
* @param f the function to solve
* @param x0 old approximation for the root
* @param y0 function value at the approximation for the root
* @param x1 last calculated approximation for the root
* @param y1 function value at the last calculated approximation
* for the root
* @param x2 bracket point (must be set to x0 if no bracket point is
* known, this will force starting with linear interpolation)
* @param y2 function value at the bracket point.
* @return the value where the function is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* @throws MathUserException if an error occurs evaluating the function
*/
private double solve(final UnivariateRealFunction f,
double x0, double y0,
double x1, double y1,
double x2, double y2)
throws MaxIterationsExceededException, MathUserException {
double delta = x1 - x0;
double oldDelta = delta;
int i = 0;
while (i < maximalIterationCount) {
if (FastMath.abs(y2) < FastMath.abs(y1)) {
// use the bracket point if is better than last approximation
x0 = x1;
x1 = x2;
x2 = x0;
y0 = y1;
y1 = y2;
y2 = y0;
if (FastMath.abs(m) <= tol ||
MathUtils.equals(fb, 0)) {
return b;
}
if (FastMath.abs(y1) <= functionValueAccuracy) {
// Avoid division by very small values. Assume
// the iteration has converged (the problem may
// still be ill conditioned)
setResult(x1, i);
return result;
}
double dx = x2 - x1;
double tolerance =
FastMath.max(relativeAccuracy * FastMath.abs(x1), absoluteAccuracy);
if (FastMath.abs(dx) <= tolerance) {
setResult(x1, i);
return result;
}
if ((FastMath.abs(oldDelta) < tolerance) ||
(FastMath.abs(y0) <= FastMath.abs(y1))) {
if (FastMath.abs(e) < tol ||
FastMath.abs(fa) <= FastMath.abs(fb)) {
// Force bisection.
delta = 0.5 * dx;
oldDelta = delta;
d = m;
e = d;
} else {
double r3 = y1 / y0;
double s = fb / fa;
double p;
double p1;
// the equality test (x0 == x2) is intentional,
// it is part of the original Brent's method,
// it should NOT be replaced by proximity test
if (x0 == x2) {
double q;
// The equality test (a == c) is intentional,
// it is part of the original Brent's method and
// it should NOT be replaced by proximity test.
if (a == c) {
// Linear interpolation.
p = dx * r3;
p1 = 1.0 - r3;
p = 2 * m * s;
q = 1 - s;
} else {
// Inverse quadratic interpolation.
double r1 = y0 / y2;
double r2 = y1 / y2;
p = r3 * (dx * r1 * (r1 - r2) - (x1 - x0) * (r2 - 1.0));
p1 = (r1 - 1.0) * (r2 - 1.0) * (r3 - 1.0);
q = fa / fc;
final double r = fb / fc;
p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
q = (q - 1) * (r - 1) * (s - 1);
}
if (p > 0.0) {
p1 = -p1;
if (p > 0) {
q = -q;
} else {
p = -p;
}
if (2.0 * p >= 1.5 * dx * p1 - FastMath.abs(tolerance * p1) ||
p >= FastMath.abs(0.5 * oldDelta * p1)) {
s = e;
e = d;
if (p >= 1.5 * m * q - FastMath.abs(tol * q) ||
p >= FastMath.abs(0.5 * s * q)) {
// Inverse quadratic interpolation gives a value
// in the wrong direction, or progress is slow.
// Fall back to bisection.
delta = 0.5 * dx;
oldDelta = delta;
d = m;
e = d;
} else {
oldDelta = delta;
delta = p / p1;
d = p / q;
}
}
// Save old X1, Y1
x0 = x1;
y0 = y1;
// Compute new X1, Y1
if (FastMath.abs(delta) > tolerance) {
x1 = x1 + delta;
} else if (dx > 0.0) {
x1 = x1 + 0.5 * tolerance;
} else if (dx <= 0.0) {
x1 = x1 - 0.5 * tolerance;
a = b;
fa = fb;
if (FastMath.abs(d) > tol) {
b += d;
} else if (m > 0) {
b += tol;
} else {
b -= tol;
}
y1 = f.value(x1);
if ((y1 > 0) == (y2 > 0)) {
x2 = x0;
y2 = y0;
delta = x1 - x0;
oldDelta = delta;
fb = computeObjectiveValue(b);
if ((fb > 0 && fc > 0) ||
(fb <= 0 && fc <= 0)) {
c = a;
fc = fa;
d = b - a;
e = d;
}
i++;
}
throw new MaxIterationsExceededException(maximalIterationCount);
}
}

29
src/main/java/org/apache/commons/math/analysis/solvers/DifferentiableUnivariateRealSolver.java Normal file
View File

@ -0,0 +1,29 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.analysis.DifferentiableUnivariateRealFunction;
/**
* Interface for (univariate real) rootfinding algorithms.
* Implementations will search for only one zero in the given interval.
*
* @version $Revision: 1034896 $ $Date: 2010-11-13 23:27:34 +0100 (Sat, 13 Nov 2010) $
*/
public interface DifferentiableUnivariateRealSolver
extends BaseUnivariateRealSolver<DifferentiableUnivariateRealFunction> {}

525
src/main/java/org/apache/commons/math/analysis/solvers/LaguerreSolver.java
View File

@ -16,304 +16,335 @@
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math.complex.Complex;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.exception.NoBracketingException;
import org.apache.commons.math.exception.NullArgumentException;
import org.apache.commons.math.exception.NoDataException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.util.FastMath;
/**
* Implements the <a href="http://mathworld.wolfram.com/LaguerresMethod.html">
* Laguerre's Method</a> for root finding of real coefficient polynomials.
* For reference, see <b>A First Course in Numerical Analysis</b>,
* ISBN 048641454X, chapter 8.
* <p>
* For reference, see
* <quote>
* <b>A First Course in Numerical Analysis</b>
* ISBN 048641454X, chapter 8.
* </quote>
* Laguerre's method is global in the sense that it can start with any initial
* approximation and be able to solve all roots from that point.</p>
* approximation and be able to solve all roots from that point.
* The algorithm requires a bracketing condition.
*
* @version $Revision$ $Date$
* @since 1.2
*/
public class LaguerreSolver extends UnivariateRealSolverImpl {
public class LaguerreSolver extends AbstractPolynomialSolver {
/** Default absolute accuracy. */
public static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
/** Complex solver. */
protected ComplexSolver complexSolver = new ComplexSolver();
/**
* Construct a solver.
* Construct a solver with default accuracies.
*/
public LaguerreSolver() {
super(100, 1E-6);
this(DEFAULT_ABSOLUTE_ACCURACY);
}
/**
* Construct a solver.
*
* @param absoluteAccuracy Absolute accuracy.
*/
public LaguerreSolver(double absoluteAccuracy) {
super(absoluteAccuracy);
}
/**
* Construct a solver.
*
* @param relativeAccuracy Relative accuracy.
* @param absoluteAccuracy Absolute accuracy.
*/
public LaguerreSolver(double relativeAccuracy,
double absoluteAccuracy) {
super(relativeAccuracy, absoluteAccuracy);
}
/**
* Construct a solver.
*
* @param relativeAccuracy Relative accuracy.
* @param absoluteAccuracy Absolute accuracy.
* @param functionValueAccuracy Function value accuracy.
*/
public LaguerreSolver(double relativeAccuracy,
double absoluteAccuracy,
double functionValueAccuracy) {
super(relativeAccuracy, absoluteAccuracy, functionValueAccuracy);
}
/**
* Find a real root in the given interval with initial value.
* <p>
* Requires bracketing condition.</p>
*
* @param f function to solve (must be polynomial)
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @param initial the start value to use
* @return the point at which the function value is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws MathUserException if an error occurs evaluating the function
* @throws IllegalArgumentException if any parameters are invalid
* {@inheritDoc}
*/
public double solve(final UnivariateRealFunction f,
final double min, final double max, final double initial)
throws ConvergenceException, MathUserException {
@Override
public double doSolve() {
double min = getMin();
double max = getMax();
double initial = getStartValue();
final double functionValueAccuracy = getFunctionValueAccuracy();
// check for zeros before verifying bracketing
if (f.value(min) == 0.0) {
return min;
}
if (f.value(max) == 0.0) {
return max;
}
if (f.value(initial) == 0.0) {
verifySequence(min, initial, max);
// Return the initial guess if it is good enough.
double yInitial = computeObjectiveValue(initial);
if (FastMath.abs(yInitial) <= functionValueAccuracy) {
return initial;
}
verifyBracketing(min, max, f);
verifySequence(min, initial, max);
if (isBracketing(min, initial, f)) {
return solve(f, min, initial);
} else {
return solve(f, initial, max);
// Return the first endpoint if it is good enough.
double yMin = computeObjectiveValue(min);
if (FastMath.abs(yMin) <= functionValueAccuracy) {
return min;
}
// Reduce interval if min and initial bracket the root.
if (yInitial * yMin < 0) {
return laguerre(min, initial, yMin, yInitial);
}
// Return the second endpoint if it is good enough.
double yMax = computeObjectiveValue(max);
if (FastMath.abs(yMax) <= functionValueAccuracy) {
return max;
}
// Reduce interval if initial and max bracket the root.
if (yInitial * yMax < 0) {
return laguerre(initial, max, yInitial, yMax);
}
throw new NoBracketingException(min, max, yMin, yMax);
}
/**
* Find a real root in the given interval.
* <p>
* Despite the bracketing condition, the root returned by solve(Complex[],
* Complex) may not be a real zero inside [min, max]. For example,
* p(x) = x^3 + 1, min = -2, max = 2, initial = 0. We can either try
* another initial value, or, as we did here, call solveAll() to obtain
* all roots and pick up the one that we're looking for.</p>
*
* @param f the function to solve
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the point at which the function value is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws MathUserException if an error occurs evaluating the function
* @throws IllegalArgumentException if any parameters are invalid
* Despite the bracketing condition, the root returned by
* {@link LaguerreSolver.ComplexSolver#solve(Complex[],Complex)} may
* not be a real zero inside {@code [min, max]}.
* For example, <code>p(x) = x<sup>3</sup> + 1,</code>
* with {@code min = -2}, {@code max = 2}, {@code initial = 0}.
* When it occurs, this code calls
* {@link LaguerreSolver.ComplexSolver#solveAll(Complex[],Complex)}
* in order to obtain all roots and picks up one real root.
*
* @param lo Lower bound of the search interval.
* @param hi Higher bound of the search interval.
* @param fLo Function value at the lower bound of the search interval.
* @param fHi Function value at the higher bound of the search interval.
* @return the point at which the function value is zero.
*/
public double solve(final UnivariateRealFunction f,
final double min, final double max)
throws ConvergenceException, MathUserException {
// check function type
if (!(f instanceof PolynomialFunction)) {
throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.FUNCTION_NOT_POLYNOMIAL);
}
// check for zeros before verifying bracketing
if (f.value(min) == 0.0) { return min; }
if (f.value(max) == 0.0) { return max; }
verifyBracketing(min, max, f);
double coefficients[] = ((PolynomialFunction) f).getCoefficients();
public double laguerre(double lo, double hi,
double fLo, double fHi) {
double result = Double.NaN;
double coefficients[] = getCoefficients();
Complex c[] = new Complex[coefficients.length];
for (int i = 0; i < coefficients.length; i++) {
c[i] = new Complex(coefficients[i], 0.0);
c[i] = new Complex(coefficients[i], 0);
}
Complex initial = new Complex(0.5 * (min + max), 0.0);
Complex z = solve(c, initial);
if (isRootOK(min, max, z)) {
setResult(z.getReal(), iterationCount);
return result;
}
// solve all roots and select the one we're seeking
Complex[] root = solveAll(c, initial);
for (int i = 0; i < root.length; i++) {
if (isRootOK(min, max, root[i])) {
setResult(root[i].getReal(), iterationCount);
return result;
Complex initial = new Complex(0.5 * (lo + hi), 0);
Complex z = complexSolver.solve(c, initial);
if (complexSolver.isRoot(lo, hi, z)) {
return z.getReal();
} else {
double r = Double.NaN;
// Solve all roots and select the one we are seeking.
Complex[] root = complexSolver.solveAll(c, initial);
for (int i = 0; i < root.length; i++) {
if (complexSolver.isRoot(lo, hi, root[i])) {
r = root[i].getReal();
break;
}
}
return r;
}
// should never happen
throw new ConvergenceException();
}
/**
* Returns true iff the given complex root is actually a real zero
* in the given interval, within the solver tolerance level.
*
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @param z the complex root
* @return true iff z is the sought-after real zero
* Class for searching all (complex) roots.
*/
protected boolean isRootOK(double min, double max, Complex z) {
double tolerance = FastMath.max(relativeAccuracy * z.abs(), absoluteAccuracy);
return (isSequence(min, z.getReal(), max)) &&
(FastMath.abs(z.getImaginary()) <= tolerance ||
z.abs() <= functionValueAccuracy);
}
/**
* Find all complex roots for the polynomial with the given coefficients,
* starting from the given initial value.
*
* @param coefficients the polynomial coefficients array
* @param initial the start value to use
* @return the point at which the function value is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws MathUserException if an error occurs evaluating the function
* @throws IllegalArgumentException if any parameters are invalid
*/
public Complex[] solveAll(double coefficients[], double initial) throws
ConvergenceException, MathUserException {
Complex c[] = new Complex[coefficients.length];
Complex z = new Complex(initial, 0.0);
for (int i = 0; i < c.length; i++) {
c[i] = new Complex(coefficients[i], 0.0);
}
return solveAll(c, z);
}
/**
* Find all complex roots for the polynomial with the given coefficients,
* starting from the given initial value.
*
* @param coefficients the polynomial coefficients array
* @param initial the start value to use
* @return the point at which the function value is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws MathUserException if an error occurs evaluating the function
* @throws IllegalArgumentException if any parameters are invalid
*/
public Complex[] solveAll(Complex coefficients[], Complex initial) throws
MaxIterationsExceededException, MathUserException {
int n = coefficients.length - 1;
int iterationCount = 0;
if (n < 1) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.NON_POSITIVE_POLYNOMIAL_DEGREE, n);
}
Complex c[] = new Complex[n+1]; // coefficients for deflated polynomial
for (int i = 0; i <= n; i++) {
c[i] = coefficients[i];
private class ComplexSolver {
/**
* Check whether the given complex root is actually a real zero
* in the given interval, within the solver tolerance level.
*
* @param min Lower bound for the interval.
* @param max Upper bound for the interval.
* @param z Complex root.
* @return {@code true} if z is a real zero.
*/
public boolean isRoot(double min, double max, Complex z) {
double tolerance = FastMath.max(getRelativeAccuracy() * z.abs(), getAbsoluteAccuracy());
return (isSequence(min, z.getReal(), max)) &&
(FastMath.abs(z.getImaginary()) <= tolerance ||
z.abs() <= getFunctionValueAccuracy());
}
// solve individual root successively
Complex root[] = new Complex[n];
for (int i = 0; i < n; i++) {
Complex subarray[] = new Complex[n-i+1];
System.arraycopy(c, 0, subarray, 0, subarray.length);
root[i] = solve(subarray, initial);
// polynomial deflation using synthetic division
Complex newc = c[n-i];
Complex oldc = null;
for (int j = n-i-1; j >= 0; j--) {
oldc = c[j];
c[j] = newc;
newc = oldc.add(newc.multiply(root[i]));
/**
* Find all complex roots for the polynomial with the given
* coefficients, starting from the given initial value.
*
* @param coefficients Polynomial coefficients.
* @param initial Start value.
* @return the point at which the function value is zero.
* @throws org.apache.commons.math.exception.TooManyEvaluationsException
* if the maximum number of evaluations is exceeded.
* @throws NullArgumentException if the {@code coefficients} is
* {@code null}.
* @throws NoDataException if the {@code coefficients} array is empty.
*/
public Complex[] solveAll(double coefficients[], double initial) {
if (coefficients == null) {
throw new NullArgumentException();
}
iterationCount += this.iterationCount;
Complex c[] = new Complex[coefficients.length];
Complex z = new Complex(initial, 0);
for (int i = 0; i < c.length; i++) {
c[i] = new Complex(coefficients[i], 0);
}
return solveAll(c, z);
}
resultComputed = true;
this.iterationCount = iterationCount;
return root;
}
/**
* Find a complex root for the polynomial with the given coefficients,
* starting from the given initial value.
*
* @param coefficients the polynomial coefficients array
* @param initial the start value to use
* @return the point at which the function value is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws MathUserException if an error occurs evaluating the function
* @throws IllegalArgumentException if any parameters are invalid
*/
public Complex solve(Complex coefficients[], Complex initial) throws
MaxIterationsExceededException, MathUserException {
int n = coefficients.length - 1;
if (n < 1) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.NON_POSITIVE_POLYNOMIAL_DEGREE, n);
}
Complex N = new Complex(n, 0.0);
Complex N1 = new Complex(n - 1, 0.0);
int i = 1;
Complex pv = null;
Complex dv = null;
Complex d2v = null;
Complex G = null;
Complex G2 = null;
Complex H = null;
Complex delta = null;
Complex denominator = null;
Complex z = initial;
Complex oldz = new Complex(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY);
while (i <= maximalIterationCount) {
// Compute pv (polynomial value), dv (derivative value), and
// d2v (second derivative value) simultaneously.
pv = coefficients[n];
dv = Complex.ZERO;
d2v = Complex.ZERO;
for (int j = n-1; j >= 0; j--) {
d2v = dv.add(z.multiply(d2v));
dv = pv.add(z.multiply(dv));
pv = coefficients[j].add(z.multiply(pv));
/**
* Find all complex roots for the polynomial with the given
* coefficients, starting from the given initial value.
*
* @param coefficients Polynomial coefficients.
* @param initial Start value.
* @return the point at which the function value is zero.
* @throws org.apache.commons.math.exception.TooManyEvaluationsException
* if the maximum number of evaluations is exceeded.
* @throws NullArgumentException if the {@code coefficients} is
* {@code null}.
* @throws NoDataException if the {@code coefficients} array is empty.
*/
public Complex[] solveAll(Complex coefficients[], Complex initial) {
if (coefficients == null) {
throw new NullArgumentException();
}
d2v = d2v.multiply(new Complex(2.0, 0.0));
// check for convergence
double tolerance = FastMath.max(relativeAccuracy * z.abs(),
absoluteAccuracy);
if ((z.subtract(oldz)).abs() <= tolerance) {
resultComputed = true;
iterationCount = i;
return z;
int n = coefficients.length - 1;
if (n == 0) {
throw new NoDataException(LocalizedFormats.POLYNOMIAL);
}
if (pv.abs() <= functionValueAccuracy) {
resultComputed = true;
iterationCount = i;
return z;
// Coefficients for deflated polynomial.
Complex c[] = new Complex[n + 1];
for (int i = 0; i <= n; i++) {
c[i] = coefficients[i];
}
// now pv != 0, calculate the new approximation
G = dv.divide(pv);
G2 = G.multiply(G);
H = G2.subtract(d2v.divide(pv));
delta = N1.multiply((N.multiply(H)).subtract(G2));
// choose a denominator larger in magnitude
Complex deltaSqrt = delta.sqrt();
Complex dplus = G.add(deltaSqrt);
Complex dminus = G.subtract(deltaSqrt);
denominator = dplus.abs() > dminus.abs() ? dplus : dminus;
// Perturb z if denominator is zero, for instance,
// p(x) = x^3 + 1, z = 0.
if (denominator.equals(new Complex(0.0, 0.0))) {
z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));
oldz = new Complex(Double.POSITIVE_INFINITY,
Double.POSITIVE_INFINITY);
} else {
oldz = z;
z = z.subtract(N.divide(denominator));
// Solve individual roots successively.
Complex root[] = new Complex[n];
for (int i = 0; i < n; i++) {
Complex subarray[] = new Complex[n - i + 1];
System.arraycopy(c, 0, subarray, 0, subarray.length);
root[i] = solve(subarray, initial);
// Polynomial deflation using synthetic division.
Complex newc = c[n - i];
Complex oldc = null;
for (int j = n - i - 1; j >= 0; j--) {
oldc = c[j];
c[j] = newc;
newc = oldc.add(newc.multiply(root[i]));
}
}
return root;
}
/**
* Find a complex root for the polynomial with the given coefficients,
* starting from the given initial value.
*
* @param coefficients Polynomial coefficients.
* @param initial Start value.
* @return the point at which the function value is zero.
* @throws org.apache.commons.math.exception.TooManyEvaluationsException
* if the maximum number of evaluations is exceeded.
* @throws NullArgumentException if the {@code coefficients} is
* {@code null}.
* @throws NoDataException if the {@code coefficients} array is empty.
*/
public Complex solve(Complex coefficients[], Complex initial) {
if (coefficients == null) {
throw new NullArgumentException();
}
int n = coefficients.length - 1;
if (n == 0) {
throw new NoDataException(LocalizedFormats.POLYNOMIAL);
}
final double absoluteAccuracy = getAbsoluteAccuracy();
final double relativeAccuracy = getRelativeAccuracy();
final double functionValueAccuracy = getFunctionValueAccuracy();
Complex N = new Complex(n, 0.0);
Complex N1 = new Complex(n - 1, 0.0);
Complex pv = null;
Complex dv = null;
Complex d2v = null;
Complex G = null;
Complex G2 = null;
Complex H = null;
Complex delta = null;
Complex denominator = null;
Complex z = initial;
Complex oldz = new Complex(Double.POSITIVE_INFINITY,
Double.POSITIVE_INFINITY);
while (true) {
// Compute pv (polynomial value), dv (derivative value), and
// d2v (second derivative value) simultaneously.
pv = coefficients[n];
dv = Complex.ZERO;
d2v = Complex.ZERO;
for (int j = n-1; j >= 0; j--) {
d2v = dv.add(z.multiply(d2v));
dv = pv.add(z.multiply(dv));
pv = coefficients[j].add(z.multiply(pv));
}
d2v = d2v.multiply(new Complex(2.0, 0.0));
// check for convergence
double tolerance = FastMath.max(relativeAccuracy * z.abs(),
absoluteAccuracy);
if ((z.subtract(oldz)).abs() <= tolerance) {
return z;
}
if (pv.abs() <= functionValueAccuracy) {
return z;
}
// now pv != 0, calculate the new approximation
G = dv.divide(pv);
G2 = G.multiply(G);
H = G2.subtract(d2v.divide(pv));
delta = N1.multiply((N.multiply(H)).subtract(G2));
// choose a denominator larger in magnitude
Complex deltaSqrt = delta.sqrt();
Complex dplus = G.add(deltaSqrt);
Complex dminus = G.subtract(deltaSqrt);
denominator = dplus.abs() > dminus.abs() ? dplus : dminus;
// Perturb z if denominator is zero, for instance,
// p(x) = x^3 + 1, z = 0.
if (denominator.equals(new Complex(0.0, 0.0))) {
z = z.add(new Complex(absoluteAccuracy, absoluteAccuracy));
oldz = new Complex(Double.POSITIVE_INFINITY,
Double.POSITIVE_INFINITY);
} else {
oldz = z;
z = z.subtract(N.divide(denominator));
}
incrementEvaluationCount();
}
i++;
}
throw new MaxIterationsExceededException(maximalIterationCount);
}
}

270
src/main/java/org/apache/commons/math/analysis/solvers/MullerSolver.java
View File

@ -16,95 +16,116 @@
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.util.FastMath;
import org.apache.commons.math.util.MathUtils;
/**
* Implements the <a href="http://mathworld.wolfram.com/MullersMethod.html">
* This class implements the <a href="http://mathworld.wolfram.com/MullersMethod.html">
* Muller's Method</a> for root finding of real univariate functions. For
* reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477,
* chapter 3.
* <p>
* Muller's method applies to both real and complex functions, but here we
* restrict ourselves to real functions. Methods solve() and solve2() find
* real zeros, using different ways to bypass complex arithmetics.</p>
* restrict ourselves to real functions.
* This class differs from {@link MullerSolver} in the way it avoids complex
* operations.</p>
* Muller's original method would have function evaluation at complex point.
* Since our f(x) is real, we have to find ways to avoid that. Bracketing
* condition is one way to go: by requiring bracketing in every iteration,
* the newly computed approximation is guaranteed to be real.</p>
* <p>
* Normally Muller's method converges quadratically in the vicinity of a
* zero, however it may be very slow in regions far away from zeros. For
* example, f(x) = exp(x) - 1, min = -50, max = 100. In such case we use
* bisection as a safety backup if it performs very poorly.</p>
* <p>
* The formulas here use divided differences directly.</p>
*
* @version $Revision$ $Date$
* @since 1.2
* @see MullerSolver2
*/
public class MullerSolver extends UnivariateRealSolverImpl {
public class MullerSolver extends AbstractUnivariateRealSolver {
/** Serializable version identifier */
private static final long serialVersionUID = 7694577816772532779L;
/** Default absolute accuracy. */
public static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
/**
* Construct a solver.
* Construct a solver with default accuracies.
*/
public MullerSolver() {
super(100, 1E-6);
this(DEFAULT_ABSOLUTE_ACCURACY);
}
/**
* Construct a solver.
*
* @param absoluteAccuracy Absolute accuracy.
*/
public MullerSolver(double absoluteAccuracy) {
super(absoluteAccuracy);
}
/**
* Construct a solver.
*
* @param relativeAccuracy Relative accuracy.
* @param absoluteAccuracy Absolute accuracy.
*/
public MullerSolver(double relativeAccuracy,
double absoluteAccuracy) {
super(relativeAccuracy, absoluteAccuracy);
}
/**
* Find a real root in the given interval with initial value.
* <p>
* Requires bracketing condition.</p>
*
* @param f the function to solve
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @param initial the start value to use
* @return the point at which the function value is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws MathUserException if an error occurs evaluating the function
* @throws IllegalArgumentException if any parameters are invalid
* {@inheritDoc}
*/
public double solve(final UnivariateRealFunction f,
final double min, final double max, final double initial)
throws MaxIterationsExceededException, MathUserException {
@Override
protected double doSolve() {
final double min = getMin();
final double max = getMax();
final double initial = getStartValue();
final double functionValueAccuracy = getFunctionValueAccuracy();
verifySequence(min, initial, max);
// check for zeros before verifying bracketing
if (f.value(min) == 0.0) { return min; }
if (f.value(max) == 0.0) { return max; }
if (f.value(initial) == 0.0) { return initial; }
final double fMin = computeObjectiveValue(min);
if (FastMath.abs(fMin) < functionValueAccuracy) {
return min;
}
final double fMax = computeObjectiveValue(max);
if (FastMath.abs(fMax) < functionValueAccuracy) {
return max;
}
final double fInitial = computeObjectiveValue(initial);
if (FastMath.abs(fInitial) < functionValueAccuracy) {
return initial;
}
verifyBracketing(min, max, f);
verifySequence(min, initial, max);
if (isBracketing(min, initial, f)) {
return solve(f, min, initial);
verifyBracketing(min, max);
if (isBracketing(min, initial)) {
return solve(min, initial, fMin, fInitial);
} else {
return solve(f, initial, max);
return solve(initial, max, fInitial, fMax);
}
}
/**
* Find a real root in the given interval.
* <p>
* Original Muller's method would have function evaluation at complex point.
* Since our f(x) is real, we have to find ways to avoid that. Bracketing
* condition is one way to go: by requiring bracketing in every iteration,
* the newly computed approximation is guaranteed to be real.</p>
* <p>
* Normally Muller's method converges quadratically in the vicinity of a
* zero, however it may be very slow in regions far away from zeros. For
* example, f(x) = exp(x) - 1, min = -50, max = 100. In such case we use
* bisection as a safety backup if it performs very poorly.</p>
* <p>
* The formulas here use divided differences directly.</p>
*
* @param f the function to solve
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the point at which the function value is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws MathUserException if an error occurs evaluating the function
* @throws IllegalArgumentException if any parameters are invalid
* @param min Lower bound for the interval.
* @param max Upper bound for the interval.
* @param fMin function value at the lower bound.
* @param fMax function value at the upper bound.
* @return the point at which the function value is zero.
*/
public double solve(final UnivariateRealFunction f,
final double min, final double max)
throws MaxIterationsExceededException, MathUserException {
private double solve(double min, double max,
double fMin, double fMax) {
final double relativeAccuracy = getRelativeAccuracy();
final double absoluteAccuracy = getAbsoluteAccuracy();
final double functionValueAccuracy = getFunctionValueAccuracy();
// [x0, x2] is the bracketing interval in each iteration
// x1 is the last approximation and an interpolation point in (x0, x2)
@ -112,23 +133,14 @@ public class MullerSolver extends UnivariateRealSolverImpl {
// d01, d12, d012 are divided differences
double x0 = min;
double y0 = f.value(x0);
double y0 = fMin;
double x2 = max;
double y2 = f.value(x2);
double y2 = fMax;
double x1 = 0.5 * (x0 + x2);
double y1 = f.value(x1);
// check for zeros before verifying bracketing
if (y0 == 0.0) {
return min;
}
if (y2 == 0.0) {
return max;
}
verifyBracketing(min, max, f);
double y1 = computeObjectiveValue(x1);
double oldx = Double.POSITIVE_INFINITY;
for (int i = 1; i <= maximalIterationCount; ++i) {
while (true) {
// Muller's method employs quadratic interpolation through
// x0, x1, x2 and x is the zero of the interpolating parabola.
// Due to bracketing condition, this parabola must have two
@ -143,17 +155,13 @@ public class MullerSolver extends UnivariateRealSolverImpl {
// xplus and xminus are two roots of parabola and at least
// one of them should lie in (x0, x2)
final double x = isSequence(x0, xplus, x2) ? xplus : xminus;
final double y = f.value(x);
final double y = computeObjectiveValue(x);
// check for convergence
final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
if (FastMath.abs(x - oldx) <= tolerance) {
setResult(x, i);
return result;
}
if (FastMath.abs(y) <= functionValueAccuracy) {
setResult(x, i);
return result;
if (FastMath.abs(x - oldx) <= tolerance ||
FastMath.abs(y) <= functionValueAccuracy) {
return x;
}
// Bisect if convergence is too slow. Bisection would waste
@ -173,118 +181,16 @@ public class MullerSolver extends UnivariateRealSolverImpl {
oldx = x;
} else {
double xm = 0.5 * (x0 + x2);
double ym = f.value(xm);
double ym = computeObjectiveValue(xm);
if (MathUtils.sign(y0) + MathUtils.sign(ym) == 0.0) {
x2 = xm; y2 = ym;
} else {
x0 = xm; y0 = ym;
}
x1 = 0.5 * (x0 + x2);
y1 = f.value(x1);
y1 = computeObjectiveValue(x1);
oldx = Double.POSITIVE_INFINITY;
}
}
throw new MaxIterationsExceededException(maximalIterationCount);
}
/**
* Find a real root in the given interval.
* <p>
* solve2() differs from solve() in the way it avoids complex operations.
* Except for the initial [min, max], solve2() does not require bracketing
* condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex
* number arises in the computation, we simply use its modulus as real
* approximation.</p>
* <p>
* Because the interval may not be bracketing, bisection alternative is
* not applicable here. However in practice our treatment usually works
* well, especially near real zeros where the imaginary part of complex
* approximation is often negligible.</p>
* <p>
* The formulas here do not use divided differences directly.</p>
*
* @param f the function to solve
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the point at which the function value is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise
* @throws MathUserException if an error occurs evaluating the function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double solve2(final UnivariateRealFunction f,
final double min, final double max)
throws MaxIterationsExceededException, MathUserException {
// x2 is the last root approximation
// x is the new approximation and new x2 for next round
// x0 < x1 < x2 does not hold here
double x0 = min;
double y0 = f.value(x0);
double x1 = max;
double y1 = f.value(x1);
double x2 = 0.5 * (x0 + x1);
double y2 = f.value(x2);
// check for zeros before verifying bracketing
if (y0 == 0.0) { return min; }
if (y1 == 0.0) { return max; }
verifyBracketing(min, max, f);
double oldx = Double.POSITIVE_INFINITY;
for (int i = 1; i <= maximalIterationCount; ++i) {
// quadratic interpolation through x0, x1, x2
final double q = (x2 - x1) / (x1 - x0);
final double a = q * (y2 - (1 + q) * y1 + q * y0);
final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
final double c = (1 + q) * y2;
final double delta = b * b - 4 * a * c;
double x;
final double denominator;
if (delta >= 0.0) {
// choose a denominator larger in magnitude
double dplus = b + FastMath.sqrt(delta);
double dminus = b - FastMath.sqrt(delta);
denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus;
} else {
// take the modulus of (B +/- FastMath.sqrt(delta))
denominator = FastMath.sqrt(b * b - delta);
}
if (denominator != 0) {
x = x2 - 2.0 * c * (x2 - x1) / denominator;
// perturb x if it exactly coincides with x1 or x2
// the equality tests here are intentional
while (x == x1 || x == x2) {
x += absoluteAccuracy;
}
} else {
// extremely rare case, get a random number to skip it
x = min + FastMath.random() * (max - min);
oldx = Double.POSITIVE_INFINITY;
}
final double y = f.value(x);
// check for convergence
final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
if (FastMath.abs(x - oldx) <= tolerance) {
setResult(x, i);
return result;
}
if (FastMath.abs(y) <= functionValueAccuracy) {
setResult(x, i);
return result;
}
// prepare the next iteration
x0 = x1;
y0 = y1;
x1 = x2;
y1 = y2;
x2 = x;
y2 = y;
oldx = x;
}
throw new MaxIterationsExceededException(maximalIterationCount);
}
}

165
src/main/java/org/apache/commons/math/analysis/solvers/MullerSolver2.java Normal file
View File

@ -0,0 +1,165 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.exception.NoBracketingException;
import org.apache.commons.math.util.FastMath;
/**
* This class implements the <a href="http://mathworld.wolfram.com/MullersMethod.html">
* Muller's Method</a> for root finding of real univariate functions. For
* reference, see <b>Elementary Numerical Analysis</b>, ISBN 0070124477,
* chapter 3.
* <p>
* Muller's method applies to both real and complex functions, but here we
* restrict ourselves to real functions.<
* This class differs from {@link MullerSolver} in the way it avoids complex
* operations.</p>
* Except for the initial [min, max], it does not require bracketing
* condition, e.g. f(x0), f(x1), f(x2) can have the same sign. If complex
* number arises in the computation, we simply use its modulus as real
* approximation.</p>
* <p>
* Because the interval may not be bracketing, bisection alternative is
* not applicable here. However in practice our treatment usually works
* well, especially near real zeroes where the imaginary part of complex
* approximation is often negligible.</p>
* <p>
* The formulas here do not use divided differences directly.</p>
*
* @version $Revision: 1034896 $ $Date: 2010-11-13 23:27:34 +0100 (Sat, 13 Nov 2010) $
* @since 1.2
* @see MullerSolver
*/
public class MullerSolver2 extends AbstractUnivariateRealSolver {
/** Serializable version identifier */
private static final long serialVersionUID = 7694577816772532779L;
/** Default absolute accuracy. */
public static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
/**
* Construct a solver with default accuracies.
*/
public MullerSolver2() {
this(DEFAULT_ABSOLUTE_ACCURACY);
}
/**
* Construct a solver.
*
* @param absoluteAccuracy Absolute accuracy.
*/
public MullerSolver2(double absoluteAccuracy) {
super(absoluteAccuracy);
}
/**
* Construct a solver.
*
* @param relativeAccuracy Relative accuracy.
* @param absoluteAccuracy Absolute accuracy.
*/
public MullerSolver2(double relativeAccuracy,
double absoluteAccuracy) {
super(relativeAccuracy, absoluteAccuracy);
}
/**
* {@inheritDoc}
*/
@Override
protected double doSolve() {
final double min = getMin();
final double max = getMax();
verifyInterval(min, max);
final double relativeAccuracy = getRelativeAccuracy();
final double absoluteAccuracy = getAbsoluteAccuracy();
final double functionValueAccuracy = getFunctionValueAccuracy();
// x2 is the last root approximation
// x is the new approximation and new x2 for next round
// x0 < x1 < x2 does not hold here
double x0 = min;
double y0 = computeObjectiveValue(x0);
if (FastMath.abs(y0) < functionValueAccuracy) {
return x0;
}
double x1 = max;
double y1 = computeObjectiveValue(x1);
if (FastMath.abs(y1) < functionValueAccuracy) {
return x1;
}
if(y0 * y1 > 0) {
throw new NoBracketingException(x0, x1, y0, y1);
}
double x2 = 0.5 * (x0 + x1);
double y2 = computeObjectiveValue(x2);
double oldx = Double.POSITIVE_INFINITY;
while (true) {
// quadratic interpolation through x0, x1, x2
final double q = (x2 - x1) / (x1 - x0);
final double a = q * (y2 - (1 + q) * y1 + q * y0);
final double b = (2 * q + 1) * y2 - (1 + q) * (1 + q) * y1 + q * q * y0;
final double c = (1 + q) * y2;
final double delta = b * b - 4 * a * c;
double x;
final double denominator;
if (delta >= 0.0) {
// choose a denominator larger in magnitude
double dplus = b + FastMath.sqrt(delta);
double dminus = b - FastMath.sqrt(delta);
denominator = FastMath.abs(dplus) > FastMath.abs(dminus) ? dplus : dminus;
} else {
// take the modulus of (B +/- FastMath.sqrt(delta))
denominator = FastMath.sqrt(b * b - delta);
}
if (denominator != 0) {
x = x2 - 2.0 * c * (x2 - x1) / denominator;
// perturb x if it exactly coincides with x1 or x2
// the equality tests here are intentional
while (x == x1 || x == x2) {
x += absoluteAccuracy;
}
} else {
// extremely rare case, get a random number to skip it
x = min + FastMath.random() * (max - min);
oldx = Double.POSITIVE_INFINITY;
}
final double y = computeObjectiveValue(x);
// check for convergence
final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
if (FastMath.abs(x - oldx) <= tolerance ||
FastMath.abs(y) <= functionValueAccuracy) {
return x;
}
// prepare the next iteration
x0 = x1;
y0 = y1;
x1 = x2;
y1 = y2;
x2 = x;
y2 = y;
oldx = x;
}
}
}

93
src/main/java/org/apache/commons/math/analysis/solvers/NewtonSolver.java
View File

@ -17,12 +17,7 @@
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.analysis.DifferentiableUnivariateRealFunction;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.util.FastMath;
/**
@ -33,76 +28,58 @@ import org.apache.commons.math.util.FastMath;
*
* @version $Revision$ $Date$
*/
public class NewtonSolver extends UnivariateRealSolverImpl {
public class NewtonSolver extends AbstractDifferentiableUnivariateRealSolver {
/** Default absolute accuracy. */
public static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
/**
* Construct a solver.
*/
public NewtonSolver() {
super(100, 1E-6);
this(DEFAULT_ABSOLUTE_ACCURACY);
}
/**
* Construct a solver.
*
* @param absoluteAccuracy Absolute accuracy.
*/
public NewtonSolver(double absoluteAccuracy) {
super(absoluteAccuracy);
}
/**
* Find a zero near the midpoint of <code>min</code> and <code>max</code>.
* Find a zero near the midpoint of {@code min} and {@code max}.
*
* @param f the function to solve
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the value where the function is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* @throws MathUserException if an error occurs evaluating the function or derivative
* @throws IllegalArgumentException if min is not less than max
* @param f Function to solve.
* @param min Lower bound for the interval?
* @param max Upper bound for the interval.
* @return the value where the function is zero.
* @throws org.apache.commons.math.exception.TooManyEvaluationsException
* if the maximum evaluation count is exceeded.
* @throws IllegalArgumentException if {@code min >= max}.
*/
public double solve(final UnivariateRealFunction f,
final double min, final double max)
throws MaxIterationsExceededException, MathUserException {
return solve(f, min, max, UnivariateRealSolverUtils.midpoint(min, max));
public double solve(final DifferentiableUnivariateRealFunction f,
final double min, final double max) {
return super.solve(f, UnivariateRealSolverUtils.midpoint(min, max));
}
/**
* Find a zero near the value <code>startValue</code>.
*
* @param f the function to solve
* @param min the lower bound for the interval (ignored).
* @param max the upper bound for the interval (ignored).
* @param startValue the start value to use.
* @return the value where the function is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* @throws MathUserException if an error occurs evaluating the function or derivative
* @throws IllegalArgumentException if startValue is not between min and max or
* if function is not a {@link DifferentiableUnivariateRealFunction} instance
* {@inheritDoc}
*/
public double solve(final UnivariateRealFunction f,
final double min, final double max, final double startValue)
throws MaxIterationsExceededException, MathUserException {
@Override
protected double doSolve() {
final double startValue = getStartValue();
final double absoluteAccuracy = getAbsoluteAccuracy();
try {
final UnivariateRealFunction derivative =
((DifferentiableUnivariateRealFunction) f).derivative();
clearResult();
verifySequence(min, startValue, max);
double x0 = startValue;
double x1;
int i = 0;
while (i < maximalIterationCount) {
x1 = x0 - (f.value(x0) / derivative.value(x0));
if (FastMath.abs(x1 - x0) <= absoluteAccuracy) {
setResult(x1, i);
return x1;
}
x0 = x1;
++i;
double x0 = startValue;
double x1;
while (true) {
x1 = x0 - (computeObjectiveValue(x0) / computeDerivativeObjectiveValue(x0));
if (FastMath.abs(x1 - x0) <= absoluteAccuracy) {
return x1;
}
throw new MaxIterationsExceededException(maximalIterationCount);
} catch (ClassCastException cce) {
throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.FUNCTION_NOT_DIFFERENTIABLE);
x0 = x1;
}
}
}

29
src/main/java/org/apache/commons/math/analysis/solvers/PolynomialSolver.java Normal file
View File

@ -0,0 +1,29 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
/**
* Interface for (polynomial) root-finding algorithms.
* Implementations will search for only one zero in the given interval.
*
* @version $Revision$ $Date$
* @since 3.0
*/
public interface PolynomialSolver
extends BaseUnivariateRealSolver<PolynomialFunction> {}

114
src/main/java/org/apache/commons/math/analysis/solvers/RiddersSolver.java
View File

@ -16,10 +16,6 @@
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.util.FastMath;
import org.apache.commons.math.util.MathUtils;
@ -35,106 +31,84 @@ import org.apache.commons.math.util.MathUtils;
* @version $Revision$ $Date$
* @since 1.2
*/
public class RiddersSolver extends UnivariateRealSolverImpl {
public class RiddersSolver extends AbstractUnivariateRealSolver {
/** Default absolute accuracy. */
public static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
/**
* Construct a solver.
* Construct a solver with default accuracy.
*/
public RiddersSolver() {
super(100, 1E-6);
this(DEFAULT_ABSOLUTE_ACCURACY);
}
/**
* Construct a solver.
*
* @param absoluteAccuracy Absolute accuracy.
*/
public RiddersSolver(double absoluteAccuracy) {
super(absoluteAccuracy);
}
/**
* Construct a solver.
*
* @param relativeAccuracy Relative accuracy.
* @param absoluteAccuracy Absolute accuracy.
*/
public RiddersSolver(double relativeAccuracy,
double absoluteAccuracy) {
super(relativeAccuracy, absoluteAccuracy);
}
/**
* Find a root in the given interval with initial value.
* <p>
* Requires bracketing condition.</p>
*
* @param f the function to solve
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @param initial the start value to use
* @return the point at which the function value is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* @throws MathUserException if an error occurs evaluating the function
* @throws IllegalArgumentException if any parameters are invalid
* {@inheritDoc}
*/
public double solve(final UnivariateRealFunction f,
final double min, final double max, final double initial)
throws MaxIterationsExceededException, MathUserException {
// check for zeros before verifying bracketing
if (f.value(min) == 0.0) { return min; }
if (f.value(max) == 0.0) { return max; }
if (f.value(initial) == 0.0) { return initial; }
verifyBracketing(min, max, f);
verifySequence(min, initial, max);
if (isBracketing(min, initial, f)) {
return solve(f, min, initial);
} else {
return solve(f, initial, max);
}
}
/**
* Find a root in the given interval.
* <p>
* Requires bracketing condition.</p>
*
* @param f the function to solve
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @return the point at which the function value is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* @throws MathUserException if an error occurs evaluating the function
* @throws IllegalArgumentException if any parameters are invalid
*/
public double solve(final UnivariateRealFunction f,
final double min, final double max)
throws MaxIterationsExceededException, MathUserException {
@Override
protected double doSolve() {
double min = getMin();
double max = getMax();
// [x1, x2] is the bracketing interval in each iteration
// x3 is the midpoint of [x1, x2]
// x is the new root approximation and an endpoint of the new interval
double x1 = min;
double y1 = f.value(x1);
double y1 = computeObjectiveValue(x1);
double x2 = max;
double y2 = f.value(x2);
double y2 = computeObjectiveValue(x2);
// check for zeros before verifying bracketing
if (y1 == 0.0) {
if (y1 == 0) {
return min;
}
if (y2 == 0.0) {
if (y2 == 0) {
return max;
}
verifyBracketing(min, max, f);
verifyBracketing(min, max);
final double absoluteAccuracy = getAbsoluteAccuracy();
final double functionValueAccuracy = getFunctionValueAccuracy();
final double relativeAccuracy = getRelativeAccuracy();
int i = 1;
double oldx = Double.POSITIVE_INFINITY;
while (i <= maximalIterationCount) {
while (true) {
// calculate the new root approximation
final double x3 = 0.5 * (x1 + x2);
final double y3 = f.value(x3);
final double y3 = computeObjectiveValue(x3);
if (FastMath.abs(y3) <= functionValueAccuracy) {
setResult(x3, i);
return result;
return x3;
}
final double delta = 1 - (y1 * y2) / (y3 * y3); // delta > 1 due to bracketing
final double correction = (MathUtils.sign(y2) * MathUtils.sign(y3)) *
(x3 - x1) / FastMath.sqrt(delta);
final double x = x3 - correction; // correction != 0
final double y = f.value(x);
final double y = computeObjectiveValue(x);
// check for convergence
final double tolerance = FastMath.max(relativeAccuracy * FastMath.abs(x), absoluteAccuracy);
if (FastMath.abs(x - oldx) <= tolerance) {
setResult(x, i);
return result;
return x;
}
if (FastMath.abs(y) <= functionValueAccuracy) {
setResult(x, i);
return result;
return x;
}
// prepare the new interval for next iteration
@ -161,8 +135,6 @@ public class RiddersSolver extends UnivariateRealSolverImpl {
}
}
oldx = x;
i++;
}
throw new MaxIterationsExceededException(maximalIterationCount);
}
}

103
src/main/java/org/apache/commons/math/analysis/solvers/SecantSolver.java
View File

@ -16,12 +16,7 @@
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.exception.NoBracketingException;
import org.apache.commons.math.util.FastMath;
@ -41,52 +36,46 @@ import org.apache.commons.math.util.FastMath;
*
* @version $Revision$ $Date$
*/
public class SecantSolver extends UnivariateRealSolverImpl {
public class SecantSolver extends AbstractUnivariateRealSolver {
/** Default absolute accuracy. */
public static final double DEFAULT_ABSOLUTE_ACCURACY = 1e-6;
/**
* Construct a solver.
* Construct a solver with default accuracy.
*/
public SecantSolver() {
super(100, 1E-6);
this(DEFAULT_ABSOLUTE_ACCURACY);
}
/**
* Find a zero in the given interval.
* Construct a solver.
*
* @param f the function to solve
* @param min the lower bound for the interval
* @param max the upper bound for the interval
* @param initial the start value to use (ignored)
* @return the value where the function is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* @throws MathUserException if an error occurs evaluating the function
* @throws IllegalArgumentException if min is not less than max or the
* signs of the values of the function at the endpoints are not opposites
* @param absoluteAccuracy Absolute accuracy.
*/
public double solve(final UnivariateRealFunction f,
final double min, final double max, final double initial)
throws MaxIterationsExceededException, MathUserException {
return solve(f, min, max);
public SecantSolver(double absoluteAccuracy) {
super(absoluteAccuracy);
}
/**
* Construct a solver.
*
* @param relativeAccuracy Relative accuracy.
* @param absoluteAccuracy Absolute accuracy.
*/
public SecantSolver(double relativeAccuracy,
double absoluteAccuracy) {
super(relativeAccuracy, absoluteAccuracy);
}
/**
* Find a zero in the given interval.
* @param f the function to solve
* @param min the lower bound for the interval.
* @param max the upper bound for the interval.
* @return the value where the function is zero
* @throws MaxIterationsExceededException if the maximum iteration count is exceeded
* @throws MathUserException if an error occurs evaluating the function
* @throws IllegalArgumentException if min is not less than max or the
* signs of the values of the function at the endpoints are not opposites
* {@inheritDoc}
*/
public double solve(final UnivariateRealFunction f,
final double min, final double max)
throws MaxIterationsExceededException, MathUserException {
clearResult();
@Override
protected double doSolve() {
double min = getMin();
double max = getMax();
verifyInterval(min, max);
final double functionValueAccuracy = getFunctionValueAccuracy();
// Index 0 is the old approximation for the root.
// Index 1 is the last calculated approximation for the root.
// Index 2 is a bracket for the root with respect to x0.
@ -94,20 +83,31 @@ public class SecantSolver extends UnivariateRealSolverImpl {
// iteration.
double x0 = min;
double x1 = max;
double y0 = f.value(x0);
double y1 = f.value(x1);
double y0 = computeObjectiveValue(x0);
// return the first endpoint if it is good enough
if (FastMath.abs(y0) <= functionValueAccuracy) {
return x0;
}
// return the second endpoint if it is good enough
double y1 = computeObjectiveValue(x1);
if (FastMath.abs(y1) <= functionValueAccuracy) {
return x1;
}
// Verify bracketing
if (y0 * y1 >= 0) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.SAME_SIGN_AT_ENDPOINTS, min, max, y0, y1);
throw new NoBracketingException(min, max, y0, y1);
}
final double absoluteAccuracy = getAbsoluteAccuracy();
final double relativeAccuracy = getRelativeAccuracy();
double x2 = x0;
double y2 = y0;
double oldDelta = x2 - x1;
int i = 0;
while (i < maximalIterationCount) {
while (true) {
if (FastMath.abs(y2) < FastMath.abs(y1)) {
x0 = x1;
x1 = x2;
@ -117,13 +117,11 @@ public class SecantSolver extends UnivariateRealSolverImpl {
y2 = y0;
}
if (FastMath.abs(y1) <= functionValueAccuracy) {
setResult(x1, i);
return result;
return x1;
}
if (FastMath.abs(oldDelta) <
FastMath.max(relativeAccuracy * FastMath.abs(x1), absoluteAccuracy)) {
setResult(x1, i);
return result;
if (FastMath.abs(oldDelta) < FastMath.max(relativeAccuracy * FastMath.abs(x1),
absoluteAccuracy)) {
return x1;
}
double delta;
if (FastMath.abs(y1) > FastMath.abs(y0)) {
@ -140,16 +138,13 @@ public class SecantSolver extends UnivariateRealSolverImpl {
x0 = x1;
y0 = y1;
x1 = x1 + delta;
y1 = f.value(x1);
y1 = computeObjectiveValue(x1);
if ((y1 > 0) == (y2 > 0)) {
// New bracket is (x0,x1).
x2 = x0;
y2 = y0;
}
oldDelta = x2 - x1;
i++;
}
throw new MaxIterationsExceededException(maximalIterationCount);
}
}

98
src/main/java/org/apache/commons/math/analysis/solvers/UnivariateRealSolver.java
View File

@ -16,104 +16,14 @@
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.ConvergingAlgorithm;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
/**
* Interface for (univariate real) rootfinding algorithms.
* <p>
* Implementations will search for only one zero in the given interval.</p>
* Interface for (univariate real) root-finding algorithms.
* Implementations will search for only one zero in the given interval.
*
* @version $Revision$ $Date$
*/
public interface UnivariateRealSolver extends ConvergingAlgorithm {
/**
* Set the function value accuracy.
* <p>
* This is used to determine when an evaluated function value or some other
* value which is used as divisor is zero.</p>
* <p>
* This is a safety guard and it shouldn't be necessary to change this in
* general.</p>
*
* @param accuracy the accuracy.
* @throws IllegalArgumentException if the accuracy can't be achieved by
* the solver or is otherwise deemed unreasonable.
*/
void setFunctionValueAccuracy(double accuracy);
/**
* Get the actual function value accuracy.
* @return the accuracy
*/
double getFunctionValueAccuracy();
/**
* Reset the actual function accuracy to the default.
* The default value is provided by the solver implementation.
*/
void resetFunctionValueAccuracy();
/**
* Solve for a zero root in the given interval.
* <p>A solver may require that the interval brackets a single zero root.
* Solvers that do require bracketing should be able to handle the case
* where one of the endpoints is itself a root.</p>
*
* @param f the function to solve.
* @param min the lower bound for the interval.
* @param max the upper bound for the interval.
* @return a value where the function is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise.
* @throws MathUserException if an error occurs evaluating the function
* @throws IllegalArgumentException if min > max or the endpoints do not
* satisfy the requirements specified by the solver
* @since 2.0
*/
double solve(UnivariateRealFunction f, double min, double max)
throws ConvergenceException, MathUserException;
/**
* Solve for a zero in the given interval, start at startValue.
* <p>A solver may require that the interval brackets a single zero root.
* Solvers that do require bracketing should be able to handle the case
* where one of the endpoints is itself a root.</p>
*
* @param f the function to solve.
* @param min the lower bound for the interval.
* @param max the upper bound for the interval.
* @param startValue the start value to use
* @return a value where the function is zero
* @throws ConvergenceException if the maximum iteration count is exceeded
* or the solver detects convergence problems otherwise.
* @throws MathUserException if an error occurs evaluating the function
* @throws IllegalArgumentException if min > max or the arguments do not
* satisfy the requirements specified by the solver
* @since 2.0
*/
double solve(UnivariateRealFunction f, double min, double max, double startValue)
throws ConvergenceException, MathUserException, IllegalArgumentException;
/**
* Get the result of the last run of the solver.
*
* @return the last result.
* @throws IllegalStateException if there is no result available, either
* because no result was yet computed or the last attempt failed.
*/
double getResult();
/**
* Get the result of the last run of the solver.
*
* @return the value of the function at the last result.
* @throws IllegalStateException if there is no result available, either
* because no result was yet computed or the last attempt failed.
*/
double getFunctionValue();
}
public interface UnivariateRealSolver
extends BaseUnivariateRealSolver<UnivariateRealFunction> {}

90
src/main/java/org/apache/commons/math/analysis/solvers/UnivariateRealSolverFactory.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.analysis.solvers;
/**
* Abstract factory class used to create {@link UnivariateRealSolver} instances.
* <p>
* Solvers implementing the following algorithms are supported:
* <ul>
* <li>Bisection</li>
* <li>Brent's method</li>
* <li>Secant method</li>
* </ul>
* Concrete factories extending this class also specify a default solver, instances of which
* are returned by <code>newDefaultSolver()</code>.</p>
* <p>
* Common usage:<pre>
* SolverFactory factory = UnivariateRealSolverFactory.newInstance();</p>
*
* // create a Brent solver to use
* BrentSolver solver = factory.newBrentSolver();
* </pre>
*
* @version $Revision$ $Date$
*/
public abstract class UnivariateRealSolverFactory {
/**
* Default constructor.
*/
protected UnivariateRealSolverFactory() {
}
/**
* Create a new factory.
* @return a new factory.
*/
public static UnivariateRealSolverFactory newInstance() {
return new UnivariateRealSolverFactoryImpl();
}
/**
* Create a new {@link UnivariateRealSolver}. The
* actual solver returned is determined by the underlying factory.
* @return the new solver.
*/
public abstract UnivariateRealSolver newDefaultSolver();
/**
* Create a new {@link UnivariateRealSolver}. The
* solver is an implementation of the bisection method.
* @return the new solver.
*/
public abstract UnivariateRealSolver newBisectionSolver();
/**
* Create a new {@link UnivariateRealSolver}. The
* solver is an implementation of the Brent method.
* @return the new solver.
*/
public abstract UnivariateRealSolver newBrentSolver();
/**
* Create a new {@link UnivariateRealSolver}. The
* solver is an implementation of Newton's Method.
* @return the new solver.
*/
public abstract UnivariateRealSolver newNewtonSolver();
/**
* Create a new {@link UnivariateRealSolver}. The
* solver is an implementation of the secant method.
* @return the new solver.
*/
public abstract UnivariateRealSolver newSecantSolver();
}

64
src/main/java/org/apache/commons/math/analysis/solvers/UnivariateRealSolverFactoryImpl.java
View File

@ -1,64 +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.analysis.solvers;
/**
* A concrete {@link UnivariateRealSolverFactory}. This is the default solver factory
* used by commons-math.
* <p>
* The default solver returned by this factory is a {@link BrentSolver}.</p>
*
* @version $Revision$ $Date$
*/
public class UnivariateRealSolverFactoryImpl extends UnivariateRealSolverFactory {
/**
* Default constructor.
*/
public UnivariateRealSolverFactoryImpl() {
}
/** {@inheritDoc} */
@Override
public UnivariateRealSolver newDefaultSolver() {
return newBrentSolver();
}
/** {@inheritDoc} */
@Override
public UnivariateRealSolver newBisectionSolver() {
return new BisectionSolver();
}
/** {@inheritDoc} */
@Override
public UnivariateRealSolver newBrentSolver() {
return new BrentSolver();
}
/** {@inheritDoc} */
@Override
public UnivariateRealSolver newNewtonSolver() {
return new NewtonSolver();
}
/** {@inheritDoc} */
@Override
public UnivariateRealSolver newSecantSolver() {
return new SecantSolver();
}
}

7
src/main/java/org/apache/commons/math/analysis/solvers/UnivariateRealSolverImpl.java
View File

@ -22,16 +22,17 @@ import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.exception.NullArgumentException;
/**
* Provide a default implementation for several functions useful to generic
* solvers.
*
* @version $Revision$ $Date$
* @deprecated in 2.2 (to be removed in 3.0). Please use
* {@link AbstractUnivariateRealSolver} instead.
*/
public abstract class UnivariateRealSolverImpl
extends ConvergingAlgorithmImpl implements UnivariateRealSolver {
@Deprecated
public abstract class UnivariateRealSolverImpl extends ConvergingAlgorithmImpl {
/** Maximum error of function. */
protected double functionValueAccuracy;

244
src/main/java/org/apache/commons/math/analysis/solvers/UnivariateRealSolverUtils.java
View File

@ -16,12 +16,12 @@
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.exception.NullArgumentException;
import org.apache.commons.math.exception.NoBracketingException;
import org.apache.commons.math.exception.NumberIsTooLargeException;
import org.apache.commons.math.exception.NotStrictlyPositiveException;
import org.apache.commons.math.util.FastMath;
/**
@ -30,56 +30,53 @@ import org.apache.commons.math.util.FastMath;
* @version $Revision$ $Date$
*/
public class UnivariateRealSolverUtils {
/**
* Default constructor.
* Class contains only static methods.
*/
private UnivariateRealSolverUtils() {
super();
}
private UnivariateRealSolverUtils() {}
/**
* Convenience method to find a zero of a univariate real function. A default
* solver is used.
*
* @param f the function.
* @param x0 the lower bound for the interval.
* @param x1 the upper bound for the interval.
* @param function Function.
* @param x0 Lower bound for the interval.
* @param x1 Upper bound for the interval.
* @return a value where the function is zero.
* @throws ConvergenceException if the iteration count was exceeded
* @throws MathUserException if an error occurs evaluating the function
* @throws IllegalArgumentException if f is null or the endpoints do not
* specify a valid interval
* specify a valid interval.
*/
public static double solve(UnivariateRealFunction f, double x0, double x1)
throws ConvergenceException, MathUserException {
setup(f);
return LazyHolder.FACTORY.newDefaultSolver().solve(f, x0, x1);
public static double solve(UnivariateRealFunction function, double x0, double x1) {
if (function == null) {
throw new NullArgumentException(LocalizedFormats.FUNCTION);
}
final UnivariateRealSolver solver = new BrentSolver();
solver.setMaxEvaluations(Integer.MAX_VALUE);
return solver.solve(function, x0, x1);
}
/**
* Convenience method to find a zero of a univariate real function. A default
* solver is used.
*
* @param f the function
* @param x0 the lower bound for the interval
* @param x1 the upper bound for the interval
* @param absoluteAccuracy the accuracy to be used by the solver
* @return a value where the function is zero
* @throws ConvergenceException if the iteration count is exceeded
* @throws MathUserException if an error occurs evaluating the function
* @throws IllegalArgumentException if f is null, the endpoints do not
* specify a valid interval, or the absoluteAccuracy is not valid for the
* default solver
* @param function Function.
* @param x0 Lower bound for the interval.
* @param x1 Upper bound for the interval.
* @param absoluteAccuracy Accuracy to be used by the solver.
* @return a value where the function is zero.
* @throws IllegalArgumentException if {@code function} is {@code null},
* the endpoints do not specify a valid interval, or the absolute accuracy
* is not valid for the default solver.
*/
public static double solve(UnivariateRealFunction f, double x0, double x1,
double absoluteAccuracy) throws ConvergenceException,
MathUserException {
setup(f);
UnivariateRealSolver solver = LazyHolder.FACTORY.newDefaultSolver();
solver.setAbsoluteAccuracy(absoluteAccuracy);
return solver.solve(f, x0, x1);
public static double solve(UnivariateRealFunction function,
double x0, double x1,
double absoluteAccuracy) {
if (function == null) {
throw new NullArgumentException(LocalizedFormats.FUNCTION);
}
final UnivariateRealSolver solver = new BrentSolver(absoluteAccuracy);
solver.setMaxEvaluations(Integer.MAX_VALUE);
return solver.solve(function, x0, x1);
}
/**
@ -110,23 +107,21 @@ public class UnivariateRealSolverUtils {
* {@link #bracket(UnivariateRealFunction, double, double, double, int)},
* explicitly specifying the maximum number of iterations.</p>
*
* @param function the function
* @param initial initial midpoint of interval being expanded to
* bracket a root
* @param lowerBound lower bound (a is never lower than this value)
* @param upperBound upper bound (b never is greater than this
* value)
* @return a two element array holding {a, b}
* @throws ConvergenceException if a root can not be bracketted
* @throws MathUserException if an error occurs evaluating the function
* @param function Function.
* @param initial Initial midpoint of interval being expanded to
* bracket a root.
* @param lowerBound Lower bound (a is never lower than this value)
* @param upperBound Upper bound (b never is greater than this
* value).
* @return a two-element array holding a and b.
* @throws NoBracketingException if a root cannot be bracketted.
* @throws IllegalArgumentException if function is null, maximumIterations
* is not positive, or initial is not between lowerBound and upperBound
* is not positive, or initial is not between lowerBound and upperBound.
*/
public static double[] bracket(UnivariateRealFunction function,
double initial, double lowerBound, double upperBound)
throws ConvergenceException, MathUserException {
return bracket( function, initial, lowerBound, upperBound,
Integer.MAX_VALUE ) ;
double initial,
double lowerBound, double upperBound) {
return bracket(function, initial, lowerBound, upperBound, Integer.MAX_VALUE);
}
/**
@ -148,42 +143,36 @@ public class UnivariateRealSolverUtils {
* <li> <code> maximumIterations</code> iterations elapse
* -- ConvergenceException </li></ul></p>
*
* @param function the function
* @param initial initial midpoint of interval being expanded to
* bracket a root
* @param lowerBound lower bound (a is never lower than this value)
* @param upperBound upper bound (b never is greater than this
* value)
* @param maximumIterations maximum number of iterations to perform
* @return a two element array holding {a, b}.
* @throws ConvergenceException if the algorithm fails to find a and b
* satisfying the desired conditions
* @throws MathUserException if an error occurs evaluating the function
* @param function Function.
* @param initial Initial midpoint of interval being expanded to
* bracket a root.
* @param lowerBound Lower bound (a is never lower than this value).
* @param upperBound Upper bound (b never is greater than this
* value).
* @param maximumIterations Maximum number of iterations to perform
* @return a two element array holding a and b.
* @throws NoBracketingException if the algorithm fails to find a and b
* satisfying the desired conditions.
* @throws IllegalArgumentException if function is null, maximumIterations
* is not positive, or initial is not between lowerBound and upperBound
* is not positive, or initial is not between lowerBound and upperBound.
*/
public static double[] bracket(UnivariateRealFunction function,
double initial, double lowerBound, double upperBound,
int maximumIterations) throws ConvergenceException,
MathUserException {
double initial,
double lowerBound, double upperBound,
int maximumIterations) {
if (function == null) {
throw new NullArgumentException(LocalizedFormats.FUNCTION);
}
if (maximumIterations <= 0) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.INVALID_MAX_ITERATIONS, maximumIterations);
}
if (initial < lowerBound || initial > upperBound || lowerBound >= upperBound) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.INVALID_BRACKETING_PARAMETERS,
lowerBound, initial, upperBound);
throw new NotStrictlyPositiveException(LocalizedFormats.INVALID_MAX_ITERATIONS, maximumIterations);
}
verifySequence(lowerBound, initial, upperBound);
double a = initial;
double b = initial;
double fa;
double fb;
int numIterations = 0 ;
int numIterations = 0;
do {
a = FastMath.max(a - 1.0, lowerBound);
@ -191,18 +180,18 @@ public class UnivariateRealSolverUtils {
fa = function.value(a);
fb = function.value(b);
numIterations++ ;
++numIterations;
} while ((fa * fb > 0.0) && (numIterations < maximumIterations) &&
((a > lowerBound) || (b < upperBound)));
if (fa * fb > 0.0 ) {
throw new ConvergenceException(
LocalizedFormats.FAILED_BRACKETING,
numIterations, maximumIterations, initial,
lowerBound, upperBound, a, b, fa, fb);
if (fa * fb > 0.0) {
throw new NoBracketingException(LocalizedFormats.FAILED_BRACKETING,
a, b, fa, fb,
numIterations, maximumIterations, initial,
lowerBound, upperBound);
}
return new double[]{a, b};
return new double[] {a, b};
}
/**
@ -213,28 +202,95 @@ public class UnivariateRealSolverUtils {
* @return the midpoint.
*/
public static double midpoint(double a, double b) {
return (a + b) * .5;
return (a + b) * 0.5;
}
/**
* Checks to see if f is null, throwing IllegalArgumentException if so.
* @param f input function
* @throws IllegalArgumentException if f is null
* Check whether the function takes opposite signs at the endpoints.
*
* @param function Function.
* @param lower Lower endpoint.
* @param upper Upper endpoint.
* @return {@code true} if the function values have opposite signs at the
* given points.
*/
private static void setup(UnivariateRealFunction f) {
if (f == null) {
public static boolean isBracketing(UnivariateRealFunction function,
final double lower,
final double upper) {
if (function == null) {
throw new NullArgumentException(LocalizedFormats.FUNCTION);
}
final double fLo = function.value(lower);
final double fHi = function.value(upper);
return (fLo > 0 && fHi < 0) || (fLo < 0 && fHi > 0);
}
/**
* Check whether the arguments form a (strictly) increasing sequence.
*
* @param start First number.
* @param mid Second number.
* @param end Third number.
* @return {@code true} if the arguments form an increasing sequence.
*/
public static boolean isSequence(final double start,
final double mid,
final double end) {
return (start < mid) && (mid < end);
}
/**
* Check that the endpoints specify an interval.
*
* @param lower Lower endpoint.
* @param upper Upper endpoint.
* @throws NumberIsTooLargeException if {@code lower >= upper}.
*/
public static void verifyInterval(final double lower,
final double upper) {
if (lower >= upper) {
throw new NumberIsTooLargeException(LocalizedFormats.ENDPOINTS_NOT_AN_INTERVAL,
lower, upper, false);
}
}
// CHECKSTYLE: stop HideUtilityClassConstructor
/** Holder for the factory.
* <p>We use here the Initialization On Demand Holder Idiom.</p>
/**
* Check that {@code lower < initial < upper}.
*
* @param lower Lower endpoint.
* @param initial Initial value.
* @param upper Upper endpoint.
* @throws NumberIsTooLargeException if {@code lower >= initial} or
* {@code initial >= upper}.
*/
private static class LazyHolder {
/** Cached solver factory */
private static final UnivariateRealSolverFactory FACTORY = UnivariateRealSolverFactory.newInstance();
public static void verifySequence(final double lower,
final double initial,
final double upper) {
verifyInterval(lower, initial);
verifyInterval(initial, upper);
}
// CHECKSTYLE: resume HideUtilityClassConstructor
/**
* Check that the endpoints specify an interval and the function takes
* opposite signs at the endpoints.
*
* @param function Function.
* @param lower Lower endpoint.
* @param upper Upper endpoint.
* @throws NoBracketingException if function has the same sign at the
* endpoints.
*/
public static void verifyBracketing(UnivariateRealFunction function,
final double lower,
final double upper) {
if (function == null) {
throw new NullArgumentException(LocalizedFormats.FUNCTION);
}
verifyInterval(lower, upper);
if (!isBracketing(function, lower, upper)) {
throw new NoBracketingException(lower, upper,
function.value(lower),
function.value(upper));
}
}
}

4
src/main/java/org/apache/commons/math/distribution/AbstractContinuousDistribution.java
View File

@ -18,7 +18,6 @@ package org.apache.commons.math.distribution;
import java.io.Serializable;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.MathException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.analysis.solvers.BrentSolver;
@ -27,6 +26,7 @@ import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.exception.NotStrictlyPositiveException;
import org.apache.commons.math.exception.OutOfRangeException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.exception.NumberIsTooLargeException;
import org.apache.commons.math.random.RandomDataImpl;
import org.apache.commons.math.util.FastMath;
@ -106,7 +106,7 @@ public abstract class AbstractContinuousDistribution
bracket = UnivariateRealSolverUtils.bracket(
rootFindingFunction, getInitialDomain(p),
lowerBound, upperBound);
} catch (ConvergenceException ex) {
} catch (NumberIsTooLargeException ex) {
/*
* Check domain endpoints to see if one gives value that is within
* the default solver's defaultAbsoluteAccuracy of 0 (will be the

124
src/main/java/org/apache/commons/math/exception/NoBracketingException.java Normal file
View File

@ -0,0 +1,124 @@
/*
* 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.Localizable;
import org.apache.commons.math.exception.util.LocalizedFormats;
/**
* Exception to be thrown when function values have the same sign at both
* ends of an interval.
*
* @since 3.0
* @version $Revision$ $Date$
*/
public class NoBracketingException extends MathIllegalArgumentException {
/** Serializable version Id. */
private static final long serialVersionUID = -3629324471511904459L;
/** Lower end of the interval. */
private final double lo;
/** Higher end of the interval. */
private final double hi;
/** Value at lower end of the interval. */
private final double fLo;
/** Value at higher end of the interval. */
private final double fHi;
/**
* Construct the exception.
*
* @param lo Lower end of the interval.
* @param hi Higher end of the interval.
* @param fLo Value at lower end of the interval.
* @param fHi Value at higher end of the interval.
*/
public NoBracketingException(double lo, double hi,
double fLo, double fHi) {
this(null, lo, hi, fLo, fHi);
}
/**
* Construct the exception with a specific context.
*
* @param specific Contextual information on what caused the exception.
* @param lo Lower end of the interval.
* @param hi Higher end of the interval.
* @param fLo Value at lower end of the interval.
* @param fHi Value at higher end of the interval.
*/
public NoBracketingException(Localizable specific,
double lo, double hi,
double fLo, double fHi) {
super(specific, LocalizedFormats.SAME_SIGN_AT_ENDPOINTS, lo, hi, fLo, fHi);
this.lo = lo;
this.hi = hi;
this.fLo = fLo;
this.fHi = fHi;
}
/**
* Construct the exception with a specific context.
*
* @param specific Contextual information on what caused the exception.
* @param lo Lower end of the interval.
* @param hi Higher end of the interval.
* @param fLo Value at lower end of the interval.
* @param fHi Value at higher end of the interval.
* @param args Additional arguments.
*/
public NoBracketingException(Localizable specific,
double lo, double hi,
double fLo, double fHi,
Object ... args) {
super(specific, LocalizedFormats.SAME_SIGN_AT_ENDPOINTS, lo, hi, fLo, fHi, args);
this.lo = lo;
this.hi = hi;
this.fLo = fLo;
this.fHi = fHi;
}
/**
* Get the lower end of the interval.
*
* @return the lower end.
*/
public double getLo() {
return lo;
}
/**
* Get the higher end of the interval.
*
* @return the higher end.
*/
public double getHi() {
return hi;
}
/**
* Get the value at the lower end of the interval.
*
* @return the value at the lower end.
*/
public double getFLo() {
return fLo;
}
/**
* Get the value at the higher end of the interval.
*
* @return the value at the higher end.
*/
public double getFHi() {
return fHi;
}
}

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

@ -98,7 +98,7 @@ public enum LocalizedFormats implements Localizable {
EVALUATION("evaluation"), /* keep */
EXPANSION_FACTOR_SMALLER_THAN_ONE("expansion factor smaller than one ({0})"),
FACTORIAL_NEGATIVE_PARAMETER("must have n >= 0 for n!, got n = {0}"),
FAILED_BRACKETING("number of iterations={0}, maximum iterations={1}, initial={2}, lower bound={3}, upper bound={4}, final a value={5}, final b value={6}, f(a)={7}, f(b)={8}"),
FAILED_BRACKETING("number of iterations={4}, maximum iterations={5}, initial={6}, lower bound={7}, upper bound={8}, final a value={0}, final b value={1}, f(a)={2}, f(b)={3}"),
FAILED_FRACTION_CONVERSION("Unable to convert {0} to fraction after {1} iterations"),
FIRST_COLUMNS_NOT_INITIALIZED_YET("first {0} columns are not initialized yet"),
FIRST_ELEMENT_NOT_ZERO("first element is not 0: {0}"),
@ -263,6 +263,7 @@ public enum LocalizedFormats implements Localizable {
PERCENTILE_IMPLEMENTATION_CANNOT_ACCESS_METHOD("cannot access {0} method in percentile implementation {1}"),
PERCENTILE_IMPLEMENTATION_UNSUPPORTED_METHOD("percentile implementation {0} does not support {1}"),
PERMUTATION_EXCEEDS_N("permutation size ({0}) exceeds permuation domain ({1})"), /* keep */
POLYNOMIAL("polynomial"), /* keep */
POLYNOMIAL_INTERPOLANTS_MISMATCH_SEGMENTS("number of polynomial interpolants must match the number of segments ({0} != {1} - 1)"),
POPULATION_LIMIT_NOT_POSITIVE("population limit has to be positive"),
POSITION_SIZE_MISMATCH_INPUT_ARRAY("position {0} and size {1} don't fit to the size of the input array {2}"),

3
src/main/java/org/apache/commons/math/ode/events/EventState.java
View File

@ -252,7 +252,8 @@ public class EventState {
}
}
};
final BrentSolver solver = new BrentSolver(maxIterationCount, convergence);
final BrentSolver solver = new BrentSolver(convergence);
solver.setMaxEvaluations(maxIterationCount);
final double root = (ta <= tb) ? solver.solve(f, ta, tb) : solver.solve(f, tb, ta);
if ((FastMath.abs(root - ta) <= convergence) &&
(FastMath.abs(root - previousEventTime) <= convergence)) {

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

@ -109,8 +109,7 @@ public class LeastSquaresConverter implements MultivariateRealFunction {
* the {@link #value(double[])} method is called)
*/
public LeastSquaresConverter(final MultivariateVectorialFunction function,
final double[] observations, final double[] weights)
throws IllegalArgumentException {
final double[] observations, final double[] weights) {
if (observations.length != weights.length) {
throw new DimensionMismatchException(observations.length, weights.length);
}

4
src/main/java/org/apache/commons/math/optimization/fitting/GaussianFitter.java
View File

@ -18,7 +18,6 @@
package org.apache.commons.math.optimization.fitting;
import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer;
import org.apache.commons.math.optimization.OptimizationException;
import org.apache.commons.math.optimization.fitting.CurveFitter;
import org.apache.commons.math.optimization.fitting.WeightedObservedPoint;
@ -86,7 +85,8 @@ public class GaussianFitter {
/**
* Fits Gaussian function to the observed points.
* It will call {@link CurveFitter#fit()}.
* It will call the base class
* {@link CurveFitter#fit(ParametricRealFunction,double[]) fit} method.
*
* @return the Gaussian function that best fits the observed points.
* @see CurveFitter

3
src/main/java/org/apache/commons/math/optimization/fitting/HarmonicFitter.java
View File

@ -77,6 +77,7 @@ public class HarmonicFitter {
* @return harmonic Function that best fits the observed points.
* @throws NumberIsTooSmallException if the sample is too short or if
* the first guess cannot be computed.
* @throws OptimizationException
*/
public HarmonicFunction fit() throws OptimizationException {
// shall we compute the first guess of the parameters ourselves ?
@ -93,7 +94,7 @@ public class HarmonicFitter {
guesser.getGuessedAmplitude(),
guesser.getGuessedPulsation(),
guesser.getGuessedPhase()
};
};
}
double[] fitted = fitter.fit(new ParametricHarmonicFunction(), parameters);

6
src/main/java/org/apache/commons/math/optimization/fitting/ParametricGaussianFunction.java
View File

@ -97,11 +97,9 @@ public class ParametricGaussianFunction implements ParametricRealFunction, Seria
* respect to {@code c}, and the partial derivative of {@code f(a, b, c,
* d)} with respect to {@code d}.
*
* @param x {@code x} value to be used as constant in {@code f(a, b, c, d)}.
* @param parameters values of {@code a}, {@code b}, {@code c}, and
* {@code d} for computation of gradient vector of {@code f(a, b, c, d)}.
* @return the gradient vector of {@code f(a, b, c, d)}.
* @param x Value to be used as constant in {@code f(x, a, b, c, d)}.
* @param parameters Values of {@code a}, {@code b}, {@code c}, and {@code d}.
* @return the gradient vector of {@code f(a, b, c, d)}.
* @throws NullArgumentException if {@code parameters} is {@code null}.
* @throws DimensionMismatchException if the size of {@code parameters} is
* not 4.

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

@ -17,7 +17,6 @@
package org.apache.commons.math.optimization.general;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.exception.NumberIsTooSmallException;
import org.apache.commons.math.exception.DimensionMismatchException;
import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction;
@ -98,9 +97,10 @@ public abstract class AbstractLeastSquaresOptimizer
*
* @throws DimensionMismatchException if the Jacobian dimension does not
* match problem dimension.
* @throws MathUserException if users jacobian function throws one
* @throws org.apache.commons.math.exception.MathUserException if the jacobian
* function throws one.
*/
protected void updateJacobian() throws MathUserException {
protected void updateJacobian() {
++jacobianEvaluations;
weightedResidualJacobian = jF.value(point);
if (weightedResidualJacobian.length != rows) {
@ -126,7 +126,7 @@ public abstract class AbstractLeastSquaresOptimizer
* @throws org.apache.commons.math.exception.TooManyEvaluationsException
* if the maximal number of evaluations is exceeded.
*/
protected void updateResidualsAndCost() throws MathUserException {
protected void updateResidualsAndCost() {
objective = computeObjectiveValue(point);
if (objective.length != rows) {
throw new DimensionMismatchException(objective.length, rows);
@ -176,9 +176,10 @@ public abstract class AbstractLeastSquaresOptimizer
* @return the covariance matrix.
* @throws org.apache.commons.math.exception.SingularMatrixException
* if the covariance matrix cannot be computed (singular problem).
* @throws MathUserException if jacobian function throws one
* @throws org.apache.commons.math.exception.MathUserException if the jacobian
* function throws one.
*/
public double[][] getCovariances() throws MathUserException {
public double[][] getCovariances() {
// set up the jacobian
updateJacobian();
@ -211,9 +212,10 @@ public abstract class AbstractLeastSquaresOptimizer
* @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.
* @throws MathUserException if jacobian function throws one
* @throws org.apache.commons.math.exception.MathUserException if the jacobian
* function throws one.
*/
public double[] guessParametersErrors() throws MathUserException {
public double[] guessParametersErrors() {
if (rows <= cols) {
throw new NumberIsTooSmallException(LocalizedFormats.NO_DEGREES_OF_FREEDOM,
rows, cols, false);
@ -231,7 +233,7 @@ public abstract class AbstractLeastSquaresOptimizer
@Override
public VectorialPointValuePair optimize(final DifferentiableMultivariateVectorialFunction f,
final double[] target, final double[] weights,
final double[] startPoint) throws MathUserException {
final double[] startPoint) {
// Reset counter.
jacobianEvaluations = 0;

27
src/main/java/org/apache/commons/math/optimization/general/NonLinearConjugateGradientOptimizer.java
View File

@ -18,7 +18,6 @@
package org.apache.commons.math.optimization.general;
import org.apache.commons.math.exception.MathIllegalStateException;
import org.apache.commons.math.exception.ConvergenceException;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.analysis.solvers.BrentSolver;
@ -40,7 +39,6 @@ import org.apache.commons.math.util.FastMath;
* @since 2.0
*
*/
public class NonLinearConjugateGradientOptimizer
extends AbstractScalarDifferentiableOptimizer {
/** Update formula for the beta parameter. */
@ -86,7 +84,8 @@ public class NonLinearConjugateGradientOptimizer
* default {@link BrentSolver Brent solver}.
*/
public void setLineSearchSolver(final UnivariateRealSolver lineSearchSolver) {
this.solver = lineSearchSolver;
solver = lineSearchSolver;
solver.setMaxEvaluations(getMaxEvaluations());
}
/**
@ -116,6 +115,7 @@ public class NonLinearConjugateGradientOptimizer
}
if (solver == null) {
solver = new BrentSolver();
solver.setMaxEvaluations(getMaxEvaluations());
}
point = getStartPoint();
final GoalType goal = getGoalType();
@ -158,15 +158,15 @@ public class NonLinearConjugateGradientOptimizer
// Find the optimal step in the search direction.
final UnivariateRealFunction lsf = new LineSearchFunction(searchDirection);
try {
final double step = solver.solve(lsf, 0, findUpperBound(lsf, 0, initialStep));
final double uB = findUpperBound(lsf, 0, initialStep);
// XXX Last parameters is set to a value clode to zero in order to
// work around the divergence problem in the "testCircleFitting"
// unit test (see MATH-439).
final double step = solver.solve(lsf, 0, uB, 1e-15);
// Validate new point.
for (int i = 0; i < point.length; ++i) {
point[i] += step * searchDirection[i];
}
} catch (org.apache.commons.math.ConvergenceException e) {
throw new ConvergenceException(); // XXX ugly workaround.
// Validate new point.
for (int i = 0; i < point.length; ++i) {
point[i] += step * searchDirection[i];
}
r = computeObjectiveGradient(point);
@ -242,7 +242,6 @@ public class NonLinearConjugateGradientOptimizer
public double[] precondition(double[] variables, double[] r) {
return r.clone();
}
}
/** Internal class for line search.
@ -267,7 +266,6 @@ public class NonLinearConjugateGradientOptimizer
/** {@inheritDoc} */
public double value(double x) throws MathUserException {
// current point in the search direction
final double[] shiftedPoint = point.clone();
for (int i = 0; i < shiftedPoint.length; ++i) {
@ -284,9 +282,6 @@ public class NonLinearConjugateGradientOptimizer
}
return dotProduct;
}
}
}

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

@ -70,7 +70,7 @@ EVALUATION_FAILED = erreur d''\u00e9valuation pour l''argument {0}
EVALUATION = \u00e9valuation
EXPANSION_FACTOR_SMALLER_THAN_ONE = facteur d''extension inf\u00e9rieur \u00e0 un ({0})
FACTORIAL_NEGATIVE_PARAMETER = n doit \u00eatre positif pour le calcul de n!, or n = {0}
FAILED_BRACKETING = nombre d''it\u00e9rations = {0}, it\u00e9rations maximum = {1}, valeur initiale = {2}, borne inf\u00e9rieure = {3}, borne sup\u00e9rieure = {4}, valeur a finale = {5}, valeur b finale = {6}, f(a) = {7}, f(b) = {8}
FAILED_BRACKETING = nombre d''it\u00e9rations = {4}, it\u00e9rations maximum = {5}, valeur initiale = {6}, borne inf\u00e9rieure = {7}, borne sup\u00e9rieure = {8}, valeur a finale = {0}, valeur b finale = {1}, f(a) = {2}, f(b) = {3}
FAILED_FRACTION_CONVERSION = Impossible de convertir {0} en fraction apr\u00e8s {1} it\u00e9rations
FIRST_COLUMNS_NOT_INITIALIZED_YET = les {0} premi\u00e8res colonnes ne sont pas encore initialis\u00e9es
FIRST_ELEMENT_NOT_ZERO = le premier \u00e9l\u00e9ment n''est pas nul : {0}
@ -235,6 +235,7 @@ OVERFLOW_IN_SUBTRACTION = d\u00e9passement de capacit\u00e9 pour la soustraction
PERCENTILE_IMPLEMENTATION_CANNOT_ACCESS_METHOD = acc\u00e8s impossible \u00e0 la m\u00e9thode {0}
PERCENTILE_IMPLEMENTATION_UNSUPPORTED_METHOD = l''implantation de pourcentage {0} ne dispose pas de la m\u00e9thode {1}
PERMUTATION_EXCEEDS_N = la taille de la permutation ({0}) d\u00e9passe le domaine de la permutation ({1})
POLYNOMIAL = polyn\u00f4me
POLYNOMIAL_INTERPOLANTS_MISMATCH_SEGMENTS = le nombre d''interpolants polyn\u00f4miaux doit correspondre au nombre de segments ({0} != {1} - 1)
POPULATION_LIMIT_NOT_POSITIVE = la limite de population doit \u00eatre positive
POSITION_SIZE_MISMATCH_INPUT_ARRAY = la position {0} et la taille {1} sont incompatibles avec la taille du tableau d''entr\u00e9e {2}

144
src/test/java/org/apache/commons/math/analysis/solvers/BisectionSolverTest.java
View File

@ -16,166 +16,90 @@
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.MathException;
import org.apache.commons.math.analysis.QuinticFunction;
import org.apache.commons.math.analysis.SinFunction;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.util.FastMath;
import junit.framework.TestCase;
import org.junit.Assert;
import org.junit.Test;
/**
* @version $Revision$ $Date$
*/
public final class BisectionSolverTest extends TestCase {
public void testSinZero() throws MathException {
public final class BisectionSolverTest {
@Test
public void testSinZero() {
UnivariateRealFunction f = new SinFunction();
double result;
UnivariateRealSolver solver = new BisectionSolver();
BisectionSolver solver = new BisectionSolver();
solver.setMaxEvaluations(50);
result = solver.solve(f, 3, 4);
assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
result = solver.solve(f, 1, 4);
assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
}
public void testQuinticZero() throws MathException {
@Test
public void testQuinticZero() {
UnivariateRealFunction f = new QuinticFunction();
double result;
UnivariateRealSolver solver = new BisectionSolver();
BisectionSolver solver = new BisectionSolver();
solver.setMaxEvaluations(50);
result = solver.solve(f, -0.2, 0.2);
assertEquals(result, 0, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 0, solver.getAbsoluteAccuracy());
result = solver.solve(f, -0.1, 0.3);
assertEquals(result, 0, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 0, solver.getAbsoluteAccuracy());
result = solver.solve(f, -0.3, 0.45);
assertEquals(result, 0, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 0, solver.getAbsoluteAccuracy());
result = solver.solve(f, 0.3, 0.7);
assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
result = solver.solve(f, 0.2, 0.6);
assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
result = solver.solve(f, 0.05, 0.95);
assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
result = solver.solve(f, 0.85, 1.25);
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
result = solver.solve(f, 0.8, 1.2);
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
result = solver.solve(f, 0.85, 1.75);
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
result = solver.solve(f, 0.55, 1.45);
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
result = solver.solve(f, 0.85, 5);
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
assertEquals(result, solver.getResult(), 0);
assertTrue(solver.getIterationCount() > 0);
Assert.assertTrue(solver.getEvaluations() > 0);
}
@Test
public void testMath369() throws Exception {
UnivariateRealFunction f = new SinFunction();
UnivariateRealSolver solver = new BisectionSolver();
assertEquals(FastMath.PI, solver.solve(f, 3.0, 3.2, 3.1), solver.getAbsoluteAccuracy());
BisectionSolver solver = new BisectionSolver();
solver.setMaxEvaluations(40);
Assert.assertEquals(FastMath.PI, solver.solve(f, 3.0, 3.2, 3.1), solver.getAbsoluteAccuracy());
}
/**
*
*/
public void testSetFunctionValueAccuracy(){
double expected = 1.0e-2;
UnivariateRealSolver solver = new BisectionSolver();
solver.setFunctionValueAccuracy(expected);
assertEquals(expected, solver.getFunctionValueAccuracy(), 1.0e-2);
}
/**
*
*/
public void testResetFunctionValueAccuracy(){
double newValue = 1.0e-2;
UnivariateRealSolver solver = new BisectionSolver();
double oldValue = solver.getFunctionValueAccuracy();
solver.setFunctionValueAccuracy(newValue);
solver.resetFunctionValueAccuracy();
assertEquals(oldValue, solver.getFunctionValueAccuracy(), 1.0e-2);
}
/**
*
*/
public void testSetAbsoluteAccuracy(){
double expected = 1.0e-2;
UnivariateRealSolver solver = new BisectionSolver();
solver.setAbsoluteAccuracy(expected);
assertEquals(expected, solver.getAbsoluteAccuracy(), 1.0e-2);
}
/**
*
*/
public void testResetAbsoluteAccuracy(){
double newValue = 1.0e-2;
UnivariateRealSolver solver = new BisectionSolver();
double oldValue = solver.getAbsoluteAccuracy();
solver.setAbsoluteAccuracy(newValue);
solver.resetAbsoluteAccuracy();
assertEquals(oldValue, solver.getAbsoluteAccuracy(), 1.0e-2);
}
/**
*
*/
public void testSetMaximalIterationCount(){
@Test
public void testSetMaximalEvaluationCount(){
int expected = 100;
UnivariateRealSolver solver = new BisectionSolver();
solver.setMaximalIterationCount(expected);
assertEquals(expected, solver.getMaximalIterationCount());
BisectionSolver solver = new BisectionSolver();
solver.setMaxEvaluations(expected);
Assert.assertEquals(expected, solver.getMaxEvaluations());
}
/**
*
*/
public void testResetMaximalIterationCount(){
int newValue = 10000;
UnivariateRealSolver solver = new BisectionSolver();
int oldValue = solver.getMaximalIterationCount();
solver.setMaximalIterationCount(newValue);
solver.resetMaximalIterationCount();
assertEquals(oldValue, solver.getMaximalIterationCount());
}
/**
*
*/
public void testSetRelativeAccuracy(){
double expected = 1.0e-2;
UnivariateRealSolver solver = new BisectionSolver();
solver.setRelativeAccuracy(expected);
assertEquals(expected, solver.getRelativeAccuracy(), 1.0e-2);
}
/**
*
*/
public void testResetRelativeAccuracy(){
double newValue = 1.0e-2;
UnivariateRealSolver solver = new BisectionSolver();
double oldValue = solver.getRelativeAccuracy();
solver.setRelativeAccuracy(newValue);
solver.resetRelativeAccuracy();
assertEquals(oldValue, solver.getRelativeAccuracy(), 1.0e-2);
}
}

279
src/test/java/org/apache/commons/math/analysis/solvers/BrentSolverTest.java
View File

@ -16,17 +16,16 @@
*/
package org.apache.commons.math.analysis.solvers;
import junit.framework.TestCase;
import org.apache.commons.math.MathException;
import org.apache.commons.math.analysis.MonitoredFunction;
import org.apache.commons.math.analysis.QuinticFunction;
import org.apache.commons.math.analysis.SinFunction;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
/**
* Testcase for UnivariateRealSolver.
* Testcase for {@link BrentSolver Brent} solver.
* Because Brent-Dekker is guaranteed to converge in less than the default
* maximum iteration count due to bisection fallback, it is quite hard to
* debug. I include measured iteration counts plus one in order to detect
@ -36,50 +35,32 @@ import org.apache.commons.math.util.FastMath;
*
* @version $Revision:670469 $ $Date:2008-06-23 10:01:38 +0200 (lun., 23 juin 2008) $
*/
public final class BrentSolverTest extends TestCase {
public BrentSolverTest(String name) {
super(name);
}
public void testSinZero() throws MathException {
// The sinus function is behaved well around the root at #pi. The second
public final class BrentSolverTest {
@Test
public void testSinZero() {
// The sinus function is behaved well around the root at pi. The second
// order derivative is zero, which means linar approximating methods will
// still converge quadratically.
UnivariateRealFunction f = new SinFunction();
double result;
UnivariateRealSolver solver = new BrentSolver();
solver.setMaxEvaluations(10);
// Somewhat benign interval. The function is monotone.
result = solver.solve(f, 3, 4);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
// 4 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 5);
// System.out.println(
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 7);
// Larger and somewhat less benign interval. The function is grows first.
result = solver.solve(f, 1, 4);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
// 5 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 6);
solver = new SecantSolver();
result = solver.solve(f, 3, 4);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
// 4 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 5);
result = solver.solve(f, 1, 4);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
// 5 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 6);
assertEquals(result, solver.getResult(), 0);
// System.out.println(
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 8);
}
public void testQuinticZero() throws MathException {
@Test
public void testQuinticZero() {
// The quintic function has zeros at 0, +-0.5 and +-1.
// Around the root of 0 the function is well behaved, with a second derivative
// of zero a 0.
@ -91,238 +72,144 @@ public final class BrentSolverTest extends TestCase {
double result;
// Brent-Dekker solver.
UnivariateRealSolver solver = new BrentSolver();
solver.setMaxEvaluations(20);
// Symmetric bracket around 0. Test whether solvers can handle hitting
// the root in the first iteration.
result = solver.solve(f, -0.2, 0.2);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 0, solver.getAbsoluteAccuracy());
assertTrue(solver.getIterationCount() <= 2);
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 0, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 3);
// 1 iterations on i586 JDK 1.4.1.
// Asymmetric bracket around 0, just for fun. Contains extremum.
result = solver.solve(f, -0.1, 0.3);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 0, solver.getAbsoluteAccuracy());
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 0, solver.getAbsoluteAccuracy());
// 5 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 6);
Assert.assertTrue(solver.getEvaluations() <= 7);
// Large bracket around 0. Contains two extrema.
result = solver.solve(f, -0.3, 0.45);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 0, solver.getAbsoluteAccuracy());
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 0, solver.getAbsoluteAccuracy());
// 6 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 7);
Assert.assertTrue(solver.getEvaluations() <= 8);
// Benign bracket around 0.5, function is monotonous.
result = solver.solve(f, 0.3, 0.7);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
// 6 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 7);
Assert.assertTrue(solver.getEvaluations() <= 9);
// Less benign bracket around 0.5, contains one extremum.
result = solver.solve(f, 0.2, 0.6);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
// 6 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 7);
// System.out.println(
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 10);
// Large, less benign bracket around 0.5, contains both extrema.
result = solver.solve(f, 0.05, 0.95);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
// 8 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 9);
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 11);
// Relatively benign bracket around 1, function is monotonous. Fast growth for x>1
// is still a problem.
result = solver.solve(f, 0.85, 1.25);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
// 8 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 9);
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 11);
// Less benign bracket around 1 with extremum.
result = solver.solve(f, 0.8, 1.2);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
// 8 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 9);
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 11);
// Large bracket around 1. Monotonous.
result = solver.solve(f, 0.85, 1.75);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
// 10 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 11);
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 13);
// Large bracket around 1. Interval contains extremum.
result = solver.solve(f, 0.55, 1.45);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
// 7 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 8);
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 10);
// Very large bracket around 1 for testing fast growth behaviour.
result = solver.solve(f, 0.85, 5);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
// 12 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 13);
// Secant solver.
solver = new SecantSolver();
result = solver.solve(f, -0.2, 0.2);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 0, solver.getAbsoluteAccuracy());
// 1 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 2);
result = solver.solve(f, -0.1, 0.3);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 0, solver.getAbsoluteAccuracy());
// 5 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 6);
result = solver.solve(f, -0.3, 0.45);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 0, solver.getAbsoluteAccuracy());
// 6 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 7);
result = solver.solve(f, 0.3, 0.7);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
// 7 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 8);
result = solver.solve(f, 0.2, 0.6);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
// 6 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 7);
result = solver.solve(f, 0.05, 0.95);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
// 8 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 9);
result = solver.solve(f, 0.85, 1.25);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
// 10 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 11);
result = solver.solve(f, 0.8, 1.2);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
// 8 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 9);
result = solver.solve(f, 0.85, 1.75);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
// 14 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 15);
// The followig is especially slow because the solver first has to reduce
// the bracket to exclude the extremum. After that, convergence is rapide.
result = solver.solve(f, 0.55, 1.45);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
// 7 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 8);
result = solver.solve(f, 0.85, 5);
//System.out.println(
// "Root: " + result + " Iterations: " + solver.getIterationCount());
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
// 14 iterations on i586 JDK 1.4.1.
assertTrue(solver.getIterationCount() <= 15);
// Static solve method
result = UnivariateRealSolverUtils.solve(f, -0.2, 0.2);
assertEquals(result, 0, solver.getAbsoluteAccuracy());
result = UnivariateRealSolverUtils.solve(f, -0.1, 0.3);
assertEquals(result, 0, 1E-8);
result = UnivariateRealSolverUtils.solve(f, -0.3, 0.45);
assertEquals(result, 0, 1E-6);
result = UnivariateRealSolverUtils.solve(f, 0.3, 0.7);
assertEquals(result, 0.5, 1E-6);
result = UnivariateRealSolverUtils.solve(f, 0.2, 0.6);
assertEquals(result, 0.5, 1E-6);
result = UnivariateRealSolverUtils.solve(f, 0.05, 0.95);
assertEquals(result, 0.5, 1E-6);
result = UnivariateRealSolverUtils.solve(f, 0.85, 1.25);
assertEquals(result, 1.0, 1E-6);
result = UnivariateRealSolverUtils.solve(f, 0.8, 1.2);
assertEquals(result, 1.0, 1E-6);
result = UnivariateRealSolverUtils.solve(f, 0.85, 1.75);
assertEquals(result, 1.0, 1E-6);
result = UnivariateRealSolverUtils.solve(f, 0.55, 1.45);
assertEquals(result, 1.0, 1E-6);
result = UnivariateRealSolverUtils.solve(f, 0.85, 5);
assertEquals(result, 1.0, 1E-6);
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 15);
}
public void testRootEndpoints() throws Exception {
@Test
public void testRootEndpoints() {
UnivariateRealFunction f = new SinFunction();
UnivariateRealSolver solver = new BrentSolver();
BrentSolver solver = new BrentSolver();
solver.setMaxEvaluations(10);
// endpoint is root
double result = solver.solve(f, FastMath.PI, 4);
assertEquals(FastMath.PI, result, solver.getAbsoluteAccuracy());
Assert.assertEquals(FastMath.PI, result, solver.getAbsoluteAccuracy());
result = solver.solve(f, 3, FastMath.PI);
assertEquals(FastMath.PI, result, solver.getAbsoluteAccuracy());
Assert.assertEquals(FastMath.PI, result, solver.getAbsoluteAccuracy());
result = solver.solve(f, FastMath.PI, 4, 3.5);
assertEquals(FastMath.PI, result, solver.getAbsoluteAccuracy());
Assert.assertEquals(FastMath.PI, result, solver.getAbsoluteAccuracy());
result = solver.solve(f, 3, FastMath.PI, 3.07);
assertEquals(FastMath.PI, result, solver.getAbsoluteAccuracy());
Assert.assertEquals(FastMath.PI, result, solver.getAbsoluteAccuracy());
}
public void testBadEndpoints() throws Exception {
@Test
public void testBadEndpoints() {
UnivariateRealFunction f = new SinFunction();
UnivariateRealSolver solver = new BrentSolver();
BrentSolver solver = new BrentSolver();
solver.setMaxEvaluations(10);
try { // bad interval
solver.solve(f, 1, -1);
fail("Expecting IllegalArgumentException - bad interval");
Assert.fail("Expecting IllegalArgumentException - bad interval");
} catch (IllegalArgumentException ex) {
// expected
}
try { // no bracket
solver.solve(f, 1, 1.5);
fail("Expecting IllegalArgumentException - non-bracketing");
Assert.fail("Expecting IllegalArgumentException - non-bracketing");
} catch (IllegalArgumentException ex) {
// expected
}
try { // no bracket
solver.solve(f, 1, 1.5, 1.2);
fail("Expecting IllegalArgumentException - non-bracketing");
Assert.fail("Expecting IllegalArgumentException - non-bracketing");
} catch (IllegalArgumentException ex) {
// expected
}
}
public void testInitialGuess() throws MathException {
@Test
public void testInitialGuess() {
MonitoredFunction f = new MonitoredFunction(new QuinticFunction());
UnivariateRealSolver solver = new BrentSolver();
BrentSolver solver = new BrentSolver();
solver.setMaxEvaluations(20);
double result;
// no guess
result = solver.solve(f, 0.6, 7.0);
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
int referenceCallsCount = f.getCallsCount();
assertTrue(referenceCallsCount >= 13);
Assert.assertTrue(referenceCallsCount >= 13);
// invalid guess (it *is* a root, but outside of the range)
try {
result = solver.solve(f, 0.6, 7.0, 0.0);
fail("an IllegalArgumentException was expected");
Assert.fail("an IllegalArgumentException was expected");
} catch (IllegalArgumentException iae) {
// expected behaviour
}
@ -330,22 +217,20 @@ public final class BrentSolverTest extends TestCase {
// bad guess
f.setCallsCount(0);
result = solver.solve(f, 0.6, 7.0, 0.61);
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
assertTrue(f.getCallsCount() > referenceCallsCount);
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertTrue(f.getCallsCount() > referenceCallsCount);
// good guess
f.setCallsCount(0);
result = solver.solve(f, 0.6, 7.0, 0.999999);
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
assertTrue(f.getCallsCount() < referenceCallsCount);
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertTrue(f.getCallsCount() < referenceCallsCount);
// perfect guess
f.setCallsCount(0);
result = solver.solve(f, 0.6, 7.0, 1.0);
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
assertEquals(0, solver.getIterationCount());
assertEquals(1, f.getCallsCount());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertEquals(1, solver.getEvaluations());
Assert.assertEquals(1, f.getCallsCount());
}
}

112
src/test/java/org/apache/commons/math/analysis/solvers/LaguerreSolverTest.java
View File

@ -22,8 +22,9 @@ import org.apache.commons.math.analysis.SinFunction;
import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math.complex.Complex;
import org.apache.commons.math.util.FastMath;
import junit.framework.TestCase;
import org.junit.Assert;
import org.junit.Test;
import org.junit.Ignore;
/**
* Testcase for Laguerre solver.
@ -35,144 +36,145 @@ import junit.framework.TestCase;
*
* @version $Revision$ $Date$
*/
public final class LaguerreSolverTest extends TestCase {
public final class LaguerreSolverTest {
/**
* Test of solver for the linear function.
*/
public void testLinearFunction() throws MathException {
@Test
public void testLinearFunction() {
double min, max, expected, result, tolerance;
// p(x) = 4x - 1
double coefficients[] = { -1.0, 4.0 };
PolynomialFunction f = new PolynomialFunction(coefficients);
UnivariateRealSolver solver = new LaguerreSolver();
LaguerreSolver solver = new LaguerreSolver();
solver.setMaxEvaluations(10);
min = 0.0; max = 1.0; expected = 0.25;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the quadratic function.
*/
public void testQuadraticFunction() throws MathException {
@Test
public void testQuadraticFunction() {
double min, max, expected, result, tolerance;
// p(x) = 2x^2 + 5x - 3 = (x+3)(2x-1)
double coefficients[] = { -3.0, 5.0, 2.0 };
PolynomialFunction f = new PolynomialFunction(coefficients);
UnivariateRealSolver solver = new LaguerreSolver();
LaguerreSolver solver = new LaguerreSolver();
solver.setMaxEvaluations(10);
min = 0.0; max = 2.0; expected = 0.5;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
min = -4.0; max = -1.0; expected = -3.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the quintic function.
*/
public void testQuinticFunction() throws MathException {
@Test
public void testQuinticFunction() {
double min, max, expected, result, tolerance;
// p(x) = x^5 - x^4 - 12x^3 + x^2 - x - 12 = (x+1)(x+3)(x-4)(x^2-x+1)
double coefficients[] = { -12.0, -1.0, 1.0, -12.0, -1.0, 1.0 };
PolynomialFunction f = new PolynomialFunction(coefficients);
UnivariateRealSolver solver = new LaguerreSolver();
LaguerreSolver solver = new LaguerreSolver();
solver.setMaxEvaluations(10);
min = -2.0; max = 2.0; expected = -1.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
min = -5.0; max = -2.5; expected = -3.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
min = 3.0; max = 6.0; expected = 4.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the quintic function using solveAll().
* XXX commented out because "solveAll" is not part of the API.
*/
public void testQuinticFunction2() throws MathException {
double initial = 0.0, tolerance;
Complex expected, result[];
// public void testQuinticFunction2() {
// double initial = 0.0, tolerance;
// Complex expected, result[];
// p(x) = x^5 + 4x^3 + x^2 + 4 = (x+1)(x^2-x+1)(x^2+4)
double coefficients[] = { 4.0, 0.0, 1.0, 4.0, 0.0, 1.0 };
LaguerreSolver solver = new LaguerreSolver();
result = solver.solveAll(coefficients, initial);
// // p(x) = x^5 + 4x^3 + x^2 + 4 = (x+1)(x^2-x+1)(x^2+4)
// double coefficients[] = { 4.0, 0.0, 1.0, 4.0, 0.0, 1.0 };
// LaguerreSolver solver = new LaguerreSolver();
// result = solver.solveAll(coefficients, initial);
expected = new Complex(0.0, -2.0);
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected.abs() * solver.getRelativeAccuracy()));
TestUtils.assertContains(result, expected, tolerance);
// expected = new Complex(0.0, -2.0);
// tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
// FastMath.abs(expected.abs() * solver.getRelativeAccuracy()));
// TestUtils.assertContains(result, expected, tolerance);
expected = new Complex(0.0, 2.0);
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected.abs() * solver.getRelativeAccuracy()));
TestUtils.assertContains(result, expected, tolerance);
// expected = new Complex(0.0, 2.0);
// tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
// FastMath.abs(expected.abs() * solver.getRelativeAccuracy()));
// TestUtils.assertContains(result, expected, tolerance);
expected = new Complex(0.5, 0.5 * FastMath.sqrt(3.0));
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected.abs() * solver.getRelativeAccuracy()));
TestUtils.assertContains(result, expected, tolerance);
// expected = new Complex(0.5, 0.5 * FastMath.sqrt(3.0));
// tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
// FastMath.abs(expected.abs() * solver.getRelativeAccuracy()));
// TestUtils.assertContains(result, expected, tolerance);
expected = new Complex(-1.0, 0.0);
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected.abs() * solver.getRelativeAccuracy()));
TestUtils.assertContains(result, expected, tolerance);
// expected = new Complex(-1.0, 0.0);
// tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
// FastMath.abs(expected.abs() * solver.getRelativeAccuracy()));
// TestUtils.assertContains(result, expected, tolerance);
expected = new Complex(0.5, -0.5 * FastMath.sqrt(3.0));
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected.abs() * solver.getRelativeAccuracy()));
TestUtils.assertContains(result, expected, tolerance);
}
// expected = new Complex(0.5, -0.5 * FastMath.sqrt(3.0));
// tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
// FastMath.abs(expected.abs() * solver.getRelativeAccuracy()));
// TestUtils.assertContains(result, expected, tolerance);
// }
/**
* Test of parameters for the solver.
*/
public void testParameters() throws Exception {
@Test
public void testParameters() {
double coefficients[] = { -3.0, 5.0, 2.0 };
PolynomialFunction f = new PolynomialFunction(coefficients);
UnivariateRealSolver solver = new LaguerreSolver();
LaguerreSolver solver = new LaguerreSolver();
solver.setMaxEvaluations(10);
try {
// bad interval
solver.solve(f, 1, -1);
fail("Expecting IllegalArgumentException - bad interval");
Assert.fail("Expecting IllegalArgumentException - bad interval");
} catch (IllegalArgumentException ex) {
// expected
}
try {
// no bracketing
solver.solve(f, 2, 3);
fail("Expecting IllegalArgumentException - no bracketing");
} catch (IllegalArgumentException ex) {
// expected
}
try {
// bad function
solver.solve(new SinFunction(), -1, 1);
fail("Expecting IllegalArgumentException - bad function");
Assert.fail("Expecting IllegalArgumentException - no bracketing");
} catch (IllegalArgumentException ex) {
// expected
}

150
src/test/java/org/apache/commons/math/analysis/solvers/MullerSolver2Test.java Normal file
View File

@ -0,0 +1,150 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.analysis.Expm1Function;
import org.apache.commons.math.analysis.QuinticFunction;
import org.apache.commons.math.analysis.SinFunction;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.exception.NumberIsTooLargeException;
import org.apache.commons.math.exception.NoBracketingException;
import org.apache.commons.math.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
/**
* Testcase for {@link MullerSolver2 Muller} solver.
* <p>
* Muller's method converges almost quadratically near roots, but it can
* be very slow in regions far away from zeros. Test runs show that for
* reasonably good initial values, for a default absolute accuracy of 1E-6,
* it generally takes 5 to 10 iterations for the solver to converge.
* <p>
* Tests for the exponential function illustrate the situations where
* Muller solver performs poorly.
*
* @version $Revision: 1034896 $ $Date: 2010-11-13 23:27:34 +0100 (Sat, 13 Nov 2010) $
*/
public final class MullerSolver2Test {
/**
* Test of solver for the sine function.
*/
@Test
public void testSinFunction() {
UnivariateRealFunction f = new SinFunction();
UnivariateRealSolver solver = new MullerSolver2();
solver.setMaxEvaluations(10);
double min, max, expected, result, tolerance;
min = 3.0; max = 4.0; expected = FastMath.PI;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
Assert.assertEquals(expected, result, tolerance);
min = -1.0; max = 1.5; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
Assert.assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the quintic function.
*/
@Test
public void testQuinticFunction() {
UnivariateRealFunction f = new QuinticFunction();
UnivariateRealSolver solver = new MullerSolver2();
solver.setMaxEvaluations(10);
double min, max, expected, result, tolerance;
min = -0.4; max = 0.2; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
Assert.assertEquals(expected, result, tolerance);
min = 0.75; max = 1.5; expected = 1.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
Assert.assertEquals(expected, result, tolerance);
min = -0.9; max = -0.2; expected = -0.5;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
Assert.assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the exponential function.
* <p>
* It takes 25 to 50 iterations for the last two tests to converge.
*/
@Test
public void testExpm1Function() {
UnivariateRealFunction f = new Expm1Function();
UnivariateRealSolver solver = new MullerSolver2();
solver.setMaxEvaluations(55);
double min, max, expected, result, tolerance;
min = -1.0; max = 2.0; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
Assert.assertEquals(expected, result, tolerance);
min = -20.0; max = 10.0; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
Assert.assertEquals(expected, result, tolerance);
min = -50.0; max = 100.0; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
Assert.assertEquals(expected, result, tolerance);
}
/**
* Test of parameters for the solver.
*/
@Test
public void testParameters() throws Exception {
UnivariateRealFunction f = new SinFunction();
UnivariateRealSolver solver = new MullerSolver2();
solver.setMaxEvaluations(10);
try {
// bad interval
solver.solve(f, 1, -1);
Assert.fail("Expecting IllegalArgumentException - bad interval");
} catch (NumberIsTooLargeException ex) {
// expected
}
try {
// no bracketing
solver.solve(f, 2, 3);
Assert.fail("Expecting IllegalArgumentException - no bracketing");
} catch (NoBracketingException ex) {
// expected
}
}
}

130
src/test/java/org/apache/commons/math/analysis/solvers/MullerSolverTest.java
View File

@ -16,17 +16,18 @@
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.MathException;
import org.apache.commons.math.analysis.Expm1Function;
import org.apache.commons.math.analysis.QuinticFunction;
import org.apache.commons.math.analysis.SinFunction;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.exception.NumberIsTooLargeException;
import org.apache.commons.math.exception.NoBracketingException;
import org.apache.commons.math.util.FastMath;
import junit.framework.TestCase;
import org.junit.Assert;
import org.junit.Test;
/**
* Testcase for Muller solver.
* Testcase for {@link MullerSolver Muller} solver.
* <p>
* Muller's method converges almost quadratically near roots, but it can
* be very slow in regions far away from zeros. Test runs show that for
@ -38,102 +39,57 @@ import junit.framework.TestCase;
*
* @version $Revision$ $Date$
*/
public final class MullerSolverTest extends TestCase {
public final class MullerSolverTest {
/**
* Test of solver for the sine function.
*/
public void testSinFunction() throws MathException {
@Test
public void testSinFunction() {
UnivariateRealFunction f = new SinFunction();
UnivariateRealSolver solver = new MullerSolver();
solver.setMaxEvaluations(10);
double min, max, expected, result, tolerance;
min = 3.0; max = 4.0; expected = FastMath.PI;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
min = -1.0; max = 1.5; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the sine function using solve2().
*/
public void testSinFunction2() throws MathException {
UnivariateRealFunction f = new SinFunction();
MullerSolver solver = new MullerSolver();
double min, max, expected, result, tolerance;
min = 3.0; max = 4.0; expected = FastMath.PI;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve2(f, min, max);
assertEquals(expected, result, tolerance);
min = -1.0; max = 1.5; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve2(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the quintic function.
*/
public void testQuinticFunction() throws MathException {
@Test
public void testQuinticFunction() {
UnivariateRealFunction f = new QuinticFunction();
UnivariateRealSolver solver = new MullerSolver();
solver.setMaxEvaluations(15);
double min, max, expected, result, tolerance;
min = -0.4; max = 0.2; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
min = 0.75; max = 1.5; expected = 1.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
min = -0.9; max = -0.2; expected = -0.5;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the quintic function using solve2().
*/
public void testQuinticFunction2() throws MathException {
UnivariateRealFunction f = new QuinticFunction();
MullerSolver solver = new MullerSolver();
double min, max, expected, result, tolerance;
min = -0.4; max = 0.2; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve2(f, min, max);
assertEquals(expected, result, tolerance);
min = 0.75; max = 1.5; expected = 1.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve2(f, min, max);
assertEquals(expected, result, tolerance);
min = -0.9; max = -0.2; expected = -0.5;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve2(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
}
/**
@ -143,78 +99,54 @@ public final class MullerSolverTest extends TestCase {
* In fact, if not for the bisection alternative, the solver would
* exceed the default maximal iteration of 100.
*/
public void testExpm1Function() throws MathException {
@Test
public void testExpm1Function() {
UnivariateRealFunction f = new Expm1Function();
UnivariateRealSolver solver = new MullerSolver();
solver.setMaxEvaluations(25);
double min, max, expected, result, tolerance;
min = -1.0; max = 2.0; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
min = -20.0; max = 10.0; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
min = -50.0; max = 100.0; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the exponential function using solve2().
* <p>
* It takes 25 to 50 iterations for the last two tests to converge.
*/
public void testExpm1Function2() throws MathException {
UnivariateRealFunction f = new Expm1Function();
MullerSolver solver = new MullerSolver();
double min, max, expected, result, tolerance;
min = -1.0; max = 2.0; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve2(f, min, max);
assertEquals(expected, result, tolerance);
min = -20.0; max = 10.0; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve2(f, min, max);
assertEquals(expected, result, tolerance);
min = -50.0; max = 100.0; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve2(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
}
/**
* Test of parameters for the solver.
*/
@Test
public void testParameters() throws Exception {
UnivariateRealFunction f = new SinFunction();
UnivariateRealSolver solver = new MullerSolver();
solver.setMaxEvaluations(10);
try {
// bad interval
solver.solve(f, 1, -1);
fail("Expecting IllegalArgumentException - bad interval");
} catch (IllegalArgumentException ex) {
double root = solver.solve(f, 1, -1);
System.out.println("root=" + root);
Assert.fail("Expecting IllegalArgumentException - bad interval");
} catch (NumberIsTooLargeException ex) {
// expected
}
try {
// no bracketing
solver.solve(f, 2, 3);
fail("Expecting IllegalArgumentException - no bracketing");
} catch (IllegalArgumentException ex) {
Assert.fail("Expecting IllegalArgumentException - no bracketing");
} catch (NoBracketingException ex) {
// expected
}
}

50
src/test/java/org/apache/commons/math/analysis/solvers/NewtonSolverTest.java
View File

@ -16,78 +16,78 @@
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.MathException;
import org.apache.commons.math.analysis.DifferentiableUnivariateRealFunction;
import org.apache.commons.math.analysis.QuinticFunction;
import org.apache.commons.math.analysis.SinFunction;
import org.apache.commons.math.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
import junit.framework.TestCase;
/**
* @version $Revision$ $Date$
*/
public final class NewtonSolverTest extends TestCase {
public final class NewtonSolverTest {
/**
*
*/
public void testSinZero() throws MathException {
@Test
public void testSinZero() {
DifferentiableUnivariateRealFunction f = new SinFunction();
double result;
UnivariateRealSolver solver = new NewtonSolver();
NewtonSolver solver = new NewtonSolver();
solver.setMaxEvaluations(10);
result = solver.solve(f, 3, 4);
assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
result = solver.solve(f, 1, 4);
assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
assertEquals(result, solver.getResult(), 0);
assertTrue(solver.getIterationCount() > 0);
Assert.assertTrue(solver.getEvaluations() > 0);
}
/**
*
*/
public void testQuinticZero() throws MathException {
@Test
public void testQuinticZero() {
DifferentiableUnivariateRealFunction f = new QuinticFunction();
double result;
UnivariateRealSolver solver = new NewtonSolver();
NewtonSolver solver = new NewtonSolver();
solver.setMaxEvaluations(30);
result = solver.solve(f, -0.2, 0.2);
assertEquals(result, 0, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 0, solver.getAbsoluteAccuracy());
result = solver.solve(f, -0.1, 0.3);
assertEquals(result, 0, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 0, solver.getAbsoluteAccuracy());
result = solver.solve(f, -0.3, 0.45);
assertEquals(result, 0, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 0, solver.getAbsoluteAccuracy());
result = solver.solve(f, 0.3, 0.7);
assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
result = solver.solve(f, 0.2, 0.6);
assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
result = solver.solve(f, 0.05, 0.95);
assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
result = solver.solve(f, 0.85, 1.25);
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
result = solver.solve(f, 0.8, 1.2);
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
result = solver.solve(f, 0.85, 1.75);
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
result = solver.solve(f, 0.55, 1.45);
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
result = solver.solve(f, 0.85, 5);
assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
}
}

52
src/test/java/org/apache/commons/math/analysis/solvers/RiddersSolverTest.java
View File

@ -16,17 +16,18 @@
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.MathException;
import org.apache.commons.math.analysis.Expm1Function;
import org.apache.commons.math.analysis.QuinticFunction;
import org.apache.commons.math.analysis.SinFunction;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.util.FastMath;
import junit.framework.TestCase;
import org.apache.commons.math.exception.NoBracketingException;
import org.apache.commons.math.exception.NumberIsTooLargeException;
import org.junit.Assert;
import org.junit.Test;
/**
* Testcase for Ridders solver.
* Testcase for {@link RiddersSolver Ridders} solver.
* <p>
* Ridders' method converges superlinearly, more specific, its rate of
* convergence is sqrt(2). Test runs show that for a default absolute
@ -36,102 +37,109 @@ import junit.framework.TestCase;
*
* @version $Revision$ $Date$
*/
public final class RiddersSolverTest extends TestCase {
public final class RiddersSolverTest {
/**
* Test of solver for the sine function.
*/
public void testSinFunction() throws MathException {
@Test
public void testSinFunction() {
UnivariateRealFunction f = new SinFunction();
UnivariateRealSolver solver = new RiddersSolver();
solver.setMaxEvaluations(10);
double min, max, expected, result, tolerance;
min = 3.0; max = 4.0; expected = FastMath.PI;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
min = -1.0; max = 1.5; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the quintic function.
*/
public void testQuinticFunction() throws MathException {
@Test
public void testQuinticFunction() {
UnivariateRealFunction f = new QuinticFunction();
UnivariateRealSolver solver = new RiddersSolver();
solver.setMaxEvaluations(15);
double min, max, expected, result, tolerance;
min = -0.4; max = 0.2; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
min = 0.75; max = 1.5; expected = 1.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
min = -0.9; max = -0.2; expected = -0.5;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
}
/**
* Test of solver for the exponential function.
*/
public void testExpm1Function() throws MathException {
@Test
public void testExpm1Function() {
UnivariateRealFunction f = new Expm1Function();
UnivariateRealSolver solver = new RiddersSolver();
solver.setMaxEvaluations(20);
double min, max, expected, result, tolerance;
min = -1.0; max = 2.0; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
min = -20.0; max = 10.0; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
min = -50.0; max = 100.0; expected = 0.0;
tolerance = FastMath.max(solver.getAbsoluteAccuracy(),
FastMath.abs(expected * solver.getRelativeAccuracy()));
result = solver.solve(f, min, max);
assertEquals(expected, result, tolerance);
Assert.assertEquals(expected, result, tolerance);
}
/**
* Test of parameters for the solver.
*/
public void testParameters() throws Exception {
@Test
public void testParameters() {
UnivariateRealFunction f = new SinFunction();
UnivariateRealSolver solver = new RiddersSolver();
solver.setMaxEvaluations(10);
try {
// bad interval
solver.solve(f, 1, -1);
fail("Expecting IllegalArgumentException - bad interval");
} catch (IllegalArgumentException ex) {
Assert.fail("Expecting IllegalArgumentException - bad interval");
} catch (NumberIsTooLargeException ex) {
// expected
}
try {
// no bracketing
solver.solve(f, 2, 3);
fail("Expecting IllegalArgumentException - no bracketing");
} catch (IllegalArgumentException ex) {
Assert.fail("Expecting IllegalArgumentException - no bracketing");
} catch (NoBracketingException ex) {
// expected
}
}

173
src/test/java/org/apache/commons/math/analysis/solvers/SecantSolverTest.java Normal file
View File

@ -0,0 +1,173 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.analysis.solvers;
import org.apache.commons.math.analysis.MonitoredFunction;
import org.apache.commons.math.analysis.QuinticFunction;
import org.apache.commons.math.analysis.SinFunction;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
/**
* Testcase for {@link SecantSolver}.
*
* @version $Revision:670469 $ $Date:2008-06-23 10:01:38 +0200 (lun., 23 juin 2008) $
*/
public final class SecantSolverTest {
@Test
public void testSinZero() {
// The sinus function is behaved well around the root at pi. The second
// order derivative is zero, which means linar approximating methods will
// still converge quadratically.
UnivariateRealFunction f = new SinFunction();
double result;
UnivariateRealSolver solver = new SecantSolver();
solver.setMaxEvaluations(10);
result = solver.solve(f, 3, 4);
//System.out.println(
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 6);
result = solver.solve(f, 1, 4);
//System.out.println(
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, FastMath.PI, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 7);
}
@Test
public void testQuinticZero() {
// The quintic function has zeros at 0, +-0.5 and +-1.
// Around the root of 0 the function is well behaved, with a second derivative
// of zero a 0.
// The other roots are less well to find, in particular the root at 1, because
// the function grows fast for x>1.
// The function has extrema (first derivative is zero) at 0.27195613 and 0.82221643,
// intervals containing these values are harder for the solvers.
UnivariateRealFunction f = new QuinticFunction();
double result;
// Brent-Dekker solver.
UnivariateRealSolver solver = new SecantSolver();
solver.setMaxEvaluations(20);
result = solver.solve(f, -0.2, 0.2);
//System.out.println(
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 0, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 3);
result = solver.solve(f, -0.1, 0.3);
//System.out.println(
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 0, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 7);
result = solver.solve(f, -0.3, 0.45);
//System.out.println(
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 0, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 8);
result = solver.solve(f, 0.3, 0.7);
//System.out.println(
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 9);
result = solver.solve(f, 0.2, 0.6);
//System.out.println(
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 8);
result = solver.solve(f, 0.05, 0.95);
//System.out.println(
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 0.5, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 10);
result = solver.solve(f, 0.85, 1.25);
//System.out.println(
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 12);
result = solver.solve(f, 0.8, 1.2);
//System.out.println(
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 10);
result = solver.solve(f, 0.85, 1.75);
//System.out.println(
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 16);
// The followig is especially slow because the solver first has to reduce
// the bracket to exclude the extremum. After that, convergence is rapide.
result = solver.solve(f, 0.55, 1.45);
//System.out.println(
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 9);
result = solver.solve(f, 0.85, 5);
//System.out.println(
// "Root: " + result + " Evaluations: " + solver.getEvaluations());
Assert.assertEquals(result, 1.0, solver.getAbsoluteAccuracy());
Assert.assertTrue(solver.getEvaluations() <= 16);
}
@Test
public void testRootEndpoints() {
UnivariateRealFunction f = new SinFunction();
SecantSolver solver = new SecantSolver();
solver.setMaxEvaluations(10);
// endpoint is root
double result = solver.solve(f, FastMath.PI, 4);
Assert.assertEquals(FastMath.PI, result, solver.getAbsoluteAccuracy());
result = solver.solve(f, 3, FastMath.PI);
Assert.assertEquals(FastMath.PI, result, solver.getAbsoluteAccuracy());
result = solver.solve(f, FastMath.PI, 4, 3.5);
Assert.assertEquals(FastMath.PI, result, solver.getAbsoluteAccuracy());
result = solver.solve(f, 3, FastMath.PI, 3.07);
Assert.assertEquals(FastMath.PI, result, solver.getAbsoluteAccuracy());
}
@Test
public void testBadEndpoints() {
UnivariateRealFunction f = new SinFunction();
SecantSolver solver = new SecantSolver();
solver.setMaxEvaluations(10);
try { // bad interval
solver.solve(f, 1, -1);
Assert.fail("Expecting IllegalArgumentException - bad interval");
} catch (IllegalArgumentException ex) {
// expected
}
try { // no bracket
solver.solve(f, 1, 1.5);
Assert.fail("Expecting IllegalArgumentException - non-bracketing");
} catch (IllegalArgumentException ex) {
// expected
}
try { // no bracket
solver.solve(f, 1, 1.5, 1.2);
Assert.fail("Expecting IllegalArgumentException - non-bracketing");
} catch (IllegalArgumentException ex) {
// expected
}
}
}

75
src/test/java/org/apache/commons/math/analysis/solvers/UnivariateRealSolverFactoryImplTest.java
View File

@ -1,75 +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.analysis.solvers;
import junit.framework.TestCase;
/**
* @version $Revision$ $Date$
*/
public class UnivariateRealSolverFactoryImplTest extends TestCase {
/** solver factory */
private UnivariateRealSolverFactory factory;
/**
* @throws java.lang.Exception
* @see junit.framework.TestCase#tearDown()
*/
@Override
protected void setUp() throws Exception {
super.setUp();
factory = new UnivariateRealSolverFactoryImpl();
}
/**
* @throws java.lang.Exception
* @see junit.framework.TestCase#tearDown()
*/
@Override
protected void tearDown() throws Exception {
factory = null;
super.tearDown();
}
public void testNewBisectionSolverValid() {
UnivariateRealSolver solver = factory.newBisectionSolver();
assertNotNull(solver);
assertTrue(solver instanceof BisectionSolver);
}
public void testNewNewtonSolverValid() {
UnivariateRealSolver solver = factory.newNewtonSolver();
assertNotNull(solver);
assertTrue(solver instanceof NewtonSolver);
}
public void testNewBrentSolverValid() {
UnivariateRealSolver solver = factory.newBrentSolver();
assertNotNull(solver);
assertTrue(solver instanceof BrentSolver);
}
public void testNewSecantSolverValid() {
UnivariateRealSolver solver = factory.newSecantSolver();
assertNotNull(solver);
assertTrue(solver instanceof SecantSolver);
}
}

107
src/test/java/org/apache/commons/math/analysis/solvers/UnivariateRealSolverUtilsTest.java
View File

@ -17,125 +17,168 @@
package org.apache.commons.math.analysis.solvers;
import junit.framework.TestCase;
import org.apache.commons.math.MathException;
import org.apache.commons.math.analysis.SinFunction;
import org.apache.commons.math.analysis.QuinticFunction;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
/**
* @version $Revision$ $Date$
*/
public class UnivariateRealSolverUtilsTest extends TestCase {
public class UnivariateRealSolverUtilsTest {
protected UnivariateRealFunction sin = new SinFunction();
public void testSolveNull() throws MathException {
@Test
public void testSolveNull() {
try {
UnivariateRealSolverUtils.solve(null, 0.0, 4.0);
fail();
Assert.fail();
} catch(IllegalArgumentException ex){
// success
}
}
public void testSolveBadEndpoints() throws MathException {
@Test
public void testSolveBadEndpoints() {
try { // bad endpoints
UnivariateRealSolverUtils.solve(sin, -0.1, 4.0, 4.0);
fail("Expecting IllegalArgumentException");
double root = UnivariateRealSolverUtils.solve(sin, 4.0, -0.1, 1e-6);
System.out.println("root=" + root);
Assert.fail("Expecting IllegalArgumentException");
} catch (IllegalArgumentException ex) {
// expected
}
}
public void testSolveBadAccuracy() throws MathException {
@Test
public void testSolveBadAccuracy() {
try { // bad accuracy
UnivariateRealSolverUtils.solve(sin, 0.0, 4.0, 0.0);
// fail("Expecting IllegalArgumentException"); // TODO needs rework since convergence behaviour was changed
// Assert.fail("Expecting IllegalArgumentException"); // TODO needs rework since convergence behaviour was changed
} catch (IllegalArgumentException ex) {
// expected
}
}
public void testSolveSin() throws MathException {
@Test
public void testSolveSin() {
double x = UnivariateRealSolverUtils.solve(sin, 1.0, 4.0);
assertEquals(FastMath.PI, x, 1.0e-4);
Assert.assertEquals(FastMath.PI, x, 1.0e-4);
}
public void testSolveAccuracyNull() throws MathException {
@Test
public void testSolveAccuracyNull() {
try {
double accuracy = 1.0e-6;
UnivariateRealSolverUtils.solve(null, 0.0, 4.0, accuracy);
fail();
Assert.fail();
} catch(IllegalArgumentException ex){
// success
}
}
public void testSolveAccuracySin() throws MathException {
@Test
public void testSolveAccuracySin() {
double accuracy = 1.0e-6;
double x = UnivariateRealSolverUtils.solve(sin, 1.0,
4.0, accuracy);
assertEquals(FastMath.PI, x, accuracy);
Assert.assertEquals(FastMath.PI, x, accuracy);
}
public void testSolveNoRoot() throws MathException {
@Test
public void testSolveNoRoot() {
try {
UnivariateRealSolverUtils.solve(sin, 1.0, 1.5);
fail("Expecting IllegalArgumentException ");
Assert.fail("Expecting IllegalArgumentException ");
} catch (IllegalArgumentException ex) {
// expected
}
}
public void testBracketSin() throws MathException {
@Test
public void testBracketSin() {
double[] result = UnivariateRealSolverUtils.bracket(sin,
0.0, -2.0, 2.0);
assertTrue(sin.value(result[0]) < 0);
assertTrue(sin.value(result[1]) > 0);
Assert.assertTrue(sin.value(result[0]) < 0);
Assert.assertTrue(sin.value(result[1]) > 0);
}
public void testBracketEndpointRoot() throws MathException {
@Test
public void testBracketEndpointRoot() {
double[] result = UnivariateRealSolverUtils.bracket(sin, 1.5, 0, 2.0);
assertEquals(0.0, sin.value(result[0]), 1.0e-15);
assertTrue(sin.value(result[1]) > 0);
Assert.assertEquals(0.0, sin.value(result[0]), 1.0e-15);
Assert.assertTrue(sin.value(result[1]) > 0);
}
public void testNullFunction() throws MathException {
@Test
public void testNullFunction() {
try { // null function
UnivariateRealSolverUtils.bracket(null, 1.5, 0, 2.0);
fail("Expecting IllegalArgumentException");
Assert.fail("Expecting IllegalArgumentException");
} catch (IllegalArgumentException ex) {
// expected
}
}
public void testBadInitial() throws MathException {
@Test
public void testBadInitial() {
try { // initial not between endpoints
UnivariateRealSolverUtils.bracket(sin, 2.5, 0, 2.0);
fail("Expecting IllegalArgumentException");
Assert.fail("Expecting IllegalArgumentException");
} catch (IllegalArgumentException ex) {
// expected
}
}
public void testBadEndpoints() throws MathException {
@Test
public void testBadEndpoints() {
try { // endpoints not valid
UnivariateRealSolverUtils.bracket(sin, 1.5, 2.0, 1.0);
fail("Expecting IllegalArgumentException");
Assert.fail("Expecting IllegalArgumentException");
} catch (IllegalArgumentException ex) {
// expected
}
}
public void testBadMaximumIterations() throws MathException {
@Test
public void testBadMaximumIterations() {
try { // bad maximum iterations
UnivariateRealSolverUtils.bracket(sin, 1.5, 0, 2.0, 0);
fail("Expecting IllegalArgumentException");
Assert.fail("Expecting IllegalArgumentException");
} catch (IllegalArgumentException ex) {
// expected
}
}
@Test
public void testMisc() {
UnivariateRealFunction f = new QuinticFunction();
double result;
// Static solve method
result = UnivariateRealSolverUtils.solve(f, -0.2, 0.2);
Assert.assertEquals(result, 0, 1E-8);
result = UnivariateRealSolverUtils.solve(f, -0.1, 0.3);
Assert.assertEquals(result, 0, 1E-8);
result = UnivariateRealSolverUtils.solve(f, -0.3, 0.45);
Assert.assertEquals(result, 0, 1E-6);
result = UnivariateRealSolverUtils.solve(f, 0.3, 0.7);
Assert.assertEquals(result, 0.5, 1E-6);
result = UnivariateRealSolverUtils.solve(f, 0.2, 0.6);
Assert.assertEquals(result, 0.5, 1E-6);
result = UnivariateRealSolverUtils.solve(f, 0.05, 0.95);
Assert.assertEquals(result, 0.5, 1E-6);
result = UnivariateRealSolverUtils.solve(f, 0.85, 1.25);
Assert.assertEquals(result, 1.0, 1E-6);
result = UnivariateRealSolverUtils.solve(f, 0.8, 1.2);
Assert.assertEquals(result, 1.0, 1E-6);
result = UnivariateRealSolverUtils.solve(f, 0.85, 1.75);
Assert.assertEquals(result, 1.0, 1E-6);
result = UnivariateRealSolverUtils.solve(f, 0.55, 1.45);
Assert.assertEquals(result, 1.0, 1E-6);
result = UnivariateRealSolverUtils.solve(f, 0.85, 5);
Assert.assertEquals(result, 1.0, 1E-6);
}
}

2
src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince853IntegratorTest.java
View File

@ -235,7 +235,7 @@ public class DormandPrince853IntegratorTest
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
assertTrue(handler.getMaximalValueError() < 5.0e-8);
assertEquals(0, handler.getMaximalValueError(), 1.1e-7);
assertEquals(0, handler.getMaximalTimeError(), 1.0e-12);
assertEquals(12.0, handler.getLastTime(), 1.0e-8 * maxStep);
integ.clearEventHandlers();

8
src/test/java/org/apache/commons/math/ode/nonstiff/HighamHall54IntegratorTest.java
View File

@ -20,8 +20,9 @@ package org.apache.commons.math.ode.nonstiff;
import junit.framework.TestCase;
import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.exception.TooManyEvaluationsException;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math.ode.FirstOrderIntegrator;
import org.apache.commons.math.ode.IntegratorException;
@ -268,9 +269,8 @@ public class HighamHall54IntegratorTest
pb.getInitialTime(), pb.getInitialState(),
pb.getFinalTime(), new double[pb.getDimension()]);
fail("an exception should have been thrown");
} catch (IntegratorException ie) {
assertTrue(ie.getCause() != null);
assertTrue(ie.getCause() instanceof ConvergenceException);
} catch (TooManyEvaluationsException tmee) {
// Expected.
}
}

97
src/test/java/org/apache/commons/math/optimization/general/CircleScalar.java Normal file
View File

@ -0,0 +1,97 @@
/*
* 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.optimization.general;
import java.awt.geom.Point2D;
import java.util.ArrayList;
import org.apache.commons.math.analysis.DifferentiableMultivariateRealFunction;
import org.apache.commons.math.analysis.MultivariateRealFunction;
import org.apache.commons.math.analysis.MultivariateVectorialFunction;
/**
* Class used in the tests.
*/
class CircleScalar implements DifferentiableMultivariateRealFunction {
private ArrayList<Point2D.Double> points;
public CircleScalar() {
points = new ArrayList<Point2D.Double>();
}
public void addPoint(double px, double py) {
points.add(new Point2D.Double(px, py));
}
public double getRadius(Point2D.Double center) {
double r = 0;
for (Point2D.Double point : points) {
r += point.distance(center);
}
return r / points.size();
}
private double[] gradient(double[] point) {
// optimal radius
Point2D.Double center = new Point2D.Double(point[0], point[1]);
double radius = getRadius(center);
// gradient of the sum of squared residuals
double dJdX = 0;
double dJdY = 0;
for (Point2D.Double pk : points) {
double dk = pk.distance(center);
dJdX += (center.x - pk.x) * (dk - radius) / dk;
dJdY += (center.y - pk.y) * (dk - radius) / dk;
}
dJdX *= 2;
dJdY *= 2;
return new double[] { dJdX, dJdY };
}
public double value(double[] variables) {
Point2D.Double center = new Point2D.Double(variables[0], variables[1]);
double radius = getRadius(center);
double sum = 0;
for (Point2D.Double point : points) {
double di = point.distance(center) - radius;
sum += di * di;
}
return sum;
}
public MultivariateVectorialFunction gradient() {
return new MultivariateVectorialFunction() {
private static final long serialVersionUID = 3174909643301201710L;
public double[] value(double[] point) {
return gradient(point);
}
};
}
public MultivariateRealFunction partialDerivative(final int k) {
return new MultivariateRealFunction() {
private static final long serialVersionUID = 3073956364104833888L;
public double value(double[] point) {
return gradient(point)[k];
}
};
}
}

98
src/test/java/org/apache/commons/math/optimization/general/CircleVectorial.java Normal file
View File

@ -0,0 +1,98 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math.optimization.general;
import java.awt.geom.Point2D;
import java.util.ArrayList;
import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction;
import org.apache.commons.math.analysis.MultivariateMatrixFunction;
/**
* Class used in the tests.
*/
class CircleVectorial implements DifferentiableMultivariateVectorialFunction {
private ArrayList<Point2D.Double> points;
public CircleVectorial() {
points = new ArrayList<Point2D.Double>();
}
public void addPoint(double px, double py) {
points.add(new Point2D.Double(px, py));
}
public int getN() {
return points.size();
}
public double getRadius(Point2D.Double center) {
double r = 0;
for (Point2D.Double point : points) {
r += point.distance(center);
}
return r / points.size();
}
private double[][] jacobian(double[] point) {
int n = points.size();
Point2D.Double center = new Point2D.Double(point[0], point[1]);
// gradient of the optimal radius
double dRdX = 0;
double dRdY = 0;
for (Point2D.Double pk : points) {
double dk = pk.distance(center);
dRdX += (center.x - pk.x) / dk;
dRdY += (center.y - pk.y) / dk;
}
dRdX /= n;
dRdY /= n;
// jacobian of the radius residuals
double[][] jacobian = new double[n][2];
for (int i = 0; i < n; ++i) {
Point2D.Double pi = points.get(i);
double di = pi.distance(center);
jacobian[i][0] = (center.x - pi.x) / di - dRdX;
jacobian[i][1] = (center.y - pi.y) / di - dRdY;
}
return jacobian;
}
public double[] value(double[] variables) {
Point2D.Double center = new Point2D.Double(variables[0], variables[1]);
double radius = getRadius(center);
double[] residuals = new double[points.size()];
for (int i = 0; i < residuals.length; ++i) {
residuals[i] = points.get(i).distance(center) - radius;
}
return residuals;
}
public MultivariateMatrixFunction jacobian() {
return new MultivariateMatrixFunction() {
private static final long serialVersionUID = -4340046230875165095L;
public double[][] value(double[] point) {
return jacobian(point);
}
};
}
}

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

@ -368,7 +368,7 @@ extends TestCase {
}
public void testMaxEvaluations() throws Exception {
Circle circle = new Circle();
CircleVectorial circle = new CircleVectorial();
circle.addPoint( 30.0, 68.0);
circle.addPoint( 50.0, -6.0);
circle.addPoint(110.0, -20.0);
@ -388,7 +388,7 @@ extends TestCase {
}
public void testCircleFitting() throws MathUserException {
Circle circle = new Circle();
CircleVectorial circle = new CircleVectorial();
circle.addPoint( 30.0, 68.0);
circle.addPoint( 50.0, -6.0);
circle.addPoint(110.0, -20.0);
@ -409,7 +409,7 @@ extends TestCase {
}
public void testCircleFittingBadInit() throws MathUserException {
Circle circle = new Circle();
CircleVectorial circle = new CircleVectorial();
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},
@ -488,86 +488,5 @@ extends TestCase {
}
};
}
}
private static class Circle implements DifferentiableMultivariateVectorialFunction, Serializable {
private static final long serialVersionUID = -7165774454925027042L;
private ArrayList<Point2D.Double> points;
public Circle() {
points = new ArrayList<Point2D.Double>();
}
public void addPoint(double px, double py) {
points.add(new Point2D.Double(px, py));
}
public int getN() {
return points.size();
}
public double getRadius(Point2D.Double center) {
double r = 0;
for (Point2D.Double point : points) {
r += point.distance(center);
}
return r / points.size();
}
private double[][] jacobian(double[] variables) {
int n = points.size();
Point2D.Double center = new Point2D.Double(variables[0], variables[1]);
// gradient of the optimal radius
double dRdX = 0;
double dRdY = 0;
for (Point2D.Double pk : points) {
double dk = pk.distance(center);
dRdX += (center.x - pk.x) / dk;
dRdY += (center.y - pk.y) / dk;
}
dRdX /= n;
dRdY /= n;
// jacobian of the radius residuals
double[][] jacobian = new double[n][2];
for (int i = 0; i < n; ++i) {
Point2D.Double pi = points.get(i);
double di = pi.distance(center);
jacobian[i][0] = (center.x - pi.x) / di - dRdX;
jacobian[i][1] = (center.y - pi.y) / di - dRdY;
}
return jacobian;
}
public double[] value(double[] variables) {
Point2D.Double center = new Point2D.Double(variables[0], variables[1]);
double radius = getRadius(center);
double[] residuals = new double[points.size()];
for (int i = 0; i < residuals.length; ++i) {
residuals[i] = points.get(i).distance(center) - radius;
}
return residuals;
}
public MultivariateMatrixFunction jacobian() {
return new MultivariateMatrixFunction() {
private static final long serialVersionUID = -4340046230875165095L;
public double[][] value(double[] point) {
return jacobian(point);
}
};
}
}
}

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

@ -25,20 +25,19 @@ import java.util.List;
import junit.framework.TestCase;
import org.apache.commons.math.analysis.DifferentiableMultivariateVectorialFunction;
import org.apache.commons.math.analysis.MultivariateMatrixFunction;
import org.apache.commons.math.exception.SingularMatrixException;
import org.apache.commons.math.exception.ConvergenceException;
import org.apache.commons.math.exception.DimensionMismatchException;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.exception.NumberIsTooSmallException;
import org.apache.commons.math.exception.SingularMatrixException;
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;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.optimization.SimpleVectorialValueChecker;
import org.apache.commons.math.optimization.VectorialPointValuePair;
import org.apache.commons.math.util.FastMath;
import org.apache.commons.math.util.MathUtils;
import org.apache.commons.math.util.FastMath;
/**
* <p>Some of the unit tests are re-implementations of the MINPACK <a
@ -109,7 +108,7 @@ public class LevenbergMarquardtOptimizerTest
super(name);
}
public void testTrivial() throws MathUserException {
public void testTrivial() {
LinearProblem problem =
new LinearProblem(new double[][] { { 2 } }, new double[] { 3 });
LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
@ -126,7 +125,7 @@ public class LevenbergMarquardtOptimizerTest
assertEquals(3.0, optimum.getValue()[0], 1.0e-10);
}
public void testQRColumnsPermutation() throws MathUserException {
public void testQRColumnsPermutation() {
LinearProblem problem =
new LinearProblem(new double[][] { { 1.0, -1.0 }, { 0.0, 2.0 }, { 1.0, -2.0 } },
@ -143,7 +142,7 @@ public class LevenbergMarquardtOptimizerTest
assertEquals(1.0, optimum.getValue()[2], 1.0e-10);
}
public void testNoDependency() throws MathUserException {
public void testNoDependency() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 2, 0, 0, 0, 0, 0 },
{ 0, 2, 0, 0, 0, 0 },
@ -162,7 +161,7 @@ public class LevenbergMarquardtOptimizerTest
}
}
public void testOneSet() throws MathUserException {
public void testOneSet() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 1, 0, 0 },
@ -178,7 +177,7 @@ public class LevenbergMarquardtOptimizerTest
assertEquals(3.0, optimum.getPoint()[2], 1.0e-10);
}
public void testTwoSets() throws MathUserException {
public void testTwoSets() {
double epsilon = 1.0e-7;
LinearProblem problem = new LinearProblem(new double[][] {
{ 2, 1, 0, 4, 0, 0 },
@ -202,7 +201,7 @@ public class LevenbergMarquardtOptimizerTest
assertEquals( 1.0 - epsilon, optimum.getPoint()[5], 1.0e-10);
}
public void testNonInvertible() throws MathUserException {
public void testNonInvertible() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 1, 2, -3 },
@ -221,7 +220,7 @@ public class LevenbergMarquardtOptimizerTest
}
}
public void testIllConditioned() throws MathUserException {
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 },
@ -254,7 +253,7 @@ public class LevenbergMarquardtOptimizerTest
assertEquals( 22.0, optimum2.getPoint()[3], 1.0e-8);
}
public void testMoreEstimatedParametersSimple() throws MathUserException {
public void testMoreEstimatedParametersSimple() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 3.0, 2.0, 0.0, 0.0 },
@ -268,7 +267,7 @@ public class LevenbergMarquardtOptimizerTest
assertEquals(0, optimizer.getRMS(), 1.0e-10);
}
public void testMoreEstimatedParametersUnsorted() throws MathUserException {
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 },
@ -288,7 +287,7 @@ public class LevenbergMarquardtOptimizerTest
assertEquals(6.0, optimum.getPointRef()[5], 1.0e-10);
}
public void testRedundantEquations() throws MathUserException {
public void testRedundantEquations() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 1.0, 1.0 },
{ 1.0, -1.0 },
@ -304,7 +303,7 @@ public class LevenbergMarquardtOptimizerTest
assertEquals(1.0, optimum.getPointRef()[1], 1.0e-10);
}
public void testInconsistentEquations() throws MathUserException {
public void testInconsistentEquations() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 1.0, 1.0 },
{ 1.0, -1.0 },
@ -316,7 +315,7 @@ public class LevenbergMarquardtOptimizerTest
assertTrue(optimizer.getRMS() > 0.1);
}
public void testInconsistentSizes() throws MathUserException {
public void testInconsistentSizes() {
LinearProblem problem =
new LinearProblem(new double[][] { { 1, 0 }, { 0, 1 } }, new double[] { -1, 1 });
LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
@ -346,8 +345,8 @@ public class LevenbergMarquardtOptimizerTest
}
}
public void testControlParameters() throws MathUserException {
Circle circle = new Circle();
public void testControlParameters() {
CircleVectorial circle = new CircleVectorial();
circle.addPoint( 30.0, 68.0);
circle.addPoint( 50.0, -6.0);
circle.addPoint(110.0, -20.0);
@ -363,7 +362,7 @@ public class LevenbergMarquardtOptimizerTest
private void checkEstimate(DifferentiableMultivariateVectorialFunction problem,
double initialStepBoundFactor, int maxCostEval,
double costRelativeTolerance, double parRelativeTolerance,
double orthoTolerance, boolean shouldFail) throws MathUserException {
double orthoTolerance, boolean shouldFail) {
try {
LevenbergMarquardtOptimizer optimizer
= new LevenbergMarquardtOptimizer(initialStepBoundFactor,
@ -382,8 +381,8 @@ public class LevenbergMarquardtOptimizerTest
}
}
public void testCircleFitting() throws MathUserException {
Circle circle = new Circle();
public void testCircleFitting() {
CircleVectorial circle = new CircleVectorial();
circle.addPoint( 30.0, 68.0);
circle.addPoint( 50.0, -6.0);
circle.addPoint(110.0, -20.0);
@ -430,8 +429,8 @@ public class LevenbergMarquardtOptimizerTest
assertEquals(0.004, errors[1], 0.001);
}
public void testCircleFittingBadInit() throws MathUserException {
Circle circle = new Circle();
public void testCircleFittingBadInit() {
CircleVectorial circle = new CircleVectorial();
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},
@ -483,7 +482,7 @@ public class LevenbergMarquardtOptimizerTest
assertEquals( 0.2075001, center.y, 1.0e-6);
}
public void testMath199() throws MathUserException {
public void testMath199() {
try {
QuadraticProblem problem = new QuadraticProblem();
problem.addPoint (0, -3.182591015485607);
@ -527,83 +526,6 @@ public class LevenbergMarquardtOptimizerTest
}
}
private static class Circle implements DifferentiableMultivariateVectorialFunction, Serializable {
private static final long serialVersionUID = -4711170319243817874L;
private ArrayList<Point2D.Double> points;
public Circle() {
points = new ArrayList<Point2D.Double>();
}
public void addPoint(double px, double py) {
points.add(new Point2D.Double(px, py));
}
public int getN() {
return points.size();
}
public double getRadius(Point2D.Double center) {
double r = 0;
for (Point2D.Double point : points) {
r += point.distance(center);
}
return r / points.size();
}
private double[][] jacobian(double[] point) {
int n = points.size();
Point2D.Double center = new Point2D.Double(point[0], point[1]);
// gradient of the optimal radius
double dRdX = 0;
double dRdY = 0;
for (Point2D.Double pk : points) {
double dk = pk.distance(center);
dRdX += (center.x - pk.x) / dk;
dRdY += (center.y - pk.y) / dk;
}
dRdX /= n;
dRdY /= n;
// jacobian of the radius residuals
double[][] jacobian = new double[n][2];
for (int i = 0; i < n; ++i) {
Point2D.Double pi = points.get(i);
double di = pi.distance(center);
jacobian[i][0] = (center.x - pi.x) / di - dRdX;
jacobian[i][1] = (center.y - pi.y) / di - dRdY;
}
return jacobian;
}
public double[] value(double[] variables) {
Point2D.Double center = new Point2D.Double(variables[0], variables[1]);
double radius = getRadius(center);
double[] residuals = new double[points.size()];
for (int i = 0; i < residuals.length; ++i) {
residuals[i] = points.get(i).distance(center) - radius;
}
return residuals;
}
public MultivariateMatrixFunction jacobian() {
return new MultivariateMatrixFunction() {
private static final long serialVersionUID = -4340046230875165095L;
public double[][] value(double[] point) {
return jacobian(point);
}
};
}
}
private static class QuadraticProblem implements DifferentiableMultivariateVectorialFunction, Serializable {
private static final long serialVersionUID = 7072187082052755854L;

213
src/test/java/org/apache/commons/math/optimization/general/NonLinearConjugateGradientOptimizerTest.java
View File

@ -21,18 +21,18 @@ import java.awt.geom.Point2D;
import java.io.Serializable;
import java.util.ArrayList;
import junit.framework.TestCase;
import org.apache.commons.math.analysis.DifferentiableMultivariateRealFunction;
import org.apache.commons.math.analysis.MultivariateRealFunction;
import org.apache.commons.math.analysis.MultivariateVectorialFunction;
import org.apache.commons.math.analysis.solvers.UnivariateRealSolver;
import org.apache.commons.math.analysis.solvers.BrentSolver;
import org.apache.commons.math.exception.MathUserException;
import org.apache.commons.math.linear.BlockRealMatrix;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.optimization.GoalType;
import org.apache.commons.math.optimization.RealPointValuePair;
import org.apache.commons.math.optimization.SimpleScalarValueChecker;
import org.junit.Assert;
import org.junit.Test;
/**
* <p>Some of the unit tests are re-implementations of the MINPACK <a
@ -96,14 +96,9 @@ import org.apache.commons.math.optimization.SimpleScalarValueChecker;
* @author Jorge J. More (original fortran minpack tests)
* @author Luc Maisonobe (non-minpack tests and minpack tests Java translation)
*/
public class NonLinearConjugateGradientOptimizerTest
extends TestCase {
public NonLinearConjugateGradientOptimizerTest(String name) {
super(name);
}
public void testTrivial() throws MathUserException {
public class NonLinearConjugateGradientOptimizerTest {
@Test
public void testTrivial() {
LinearProblem problem =
new LinearProblem(new double[][] { { 2 } }, new double[] { 3 });
NonLinearConjugateGradientOptimizer optimizer =
@ -112,12 +107,12 @@ extends TestCase {
optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-6, 1.0e-6));
RealPointValuePair optimum =
optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 0 });
assertEquals(1.5, optimum.getPoint()[0], 1.0e-10);
assertEquals(0.0, optimum.getValue(), 1.0e-10);
Assert.assertEquals(1.5, optimum.getPoint()[0], 1.0e-10);
Assert.assertEquals(0.0, optimum.getValue(), 1.0e-10);
}
public void testColumnsPermutation() throws MathUserException {
@Test
public void testColumnsPermutation() {
LinearProblem problem =
new LinearProblem(new double[][] { { 1.0, -1.0 }, { 0.0, 2.0 }, { 1.0, -2.0 } },
new double[] { 4.0, 6.0, 1.0 });
@ -128,13 +123,14 @@ extends TestCase {
optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-6, 1.0e-6));
RealPointValuePair optimum =
optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 0, 0 });
assertEquals(7.0, optimum.getPoint()[0], 1.0e-10);
assertEquals(3.0, optimum.getPoint()[1], 1.0e-10);
assertEquals(0.0, optimum.getValue(), 1.0e-10);
Assert.assertEquals(7.0, optimum.getPoint()[0], 1.0e-10);
Assert.assertEquals(3.0, optimum.getPoint()[1], 1.0e-10);
Assert.assertEquals(0.0, optimum.getValue(), 1.0e-10);
}
public void testNoDependency() throws MathUserException {
@Test
public void testNoDependency() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 2, 0, 0, 0, 0, 0 },
{ 0, 2, 0, 0, 0, 0 },
@ -150,12 +146,12 @@ extends TestCase {
RealPointValuePair optimum =
optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 0, 0, 0, 0, 0, 0 });
for (int i = 0; i < problem.target.length; ++i) {
assertEquals(0.55 * i, optimum.getPoint()[i], 1.0e-10);
Assert.assertEquals(0.55 * i, optimum.getPoint()[i], 1.0e-10);
}
}
public void testOneSet() throws MathUserException {
@Test
public void testOneSet() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 1, 0, 0 },
{ -1, 1, 0 },
@ -167,13 +163,14 @@ extends TestCase {
optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-6, 1.0e-6));
RealPointValuePair optimum =
optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 0, 0, 0 });
assertEquals(1.0, optimum.getPoint()[0], 1.0e-10);
assertEquals(2.0, optimum.getPoint()[1], 1.0e-10);
assertEquals(3.0, optimum.getPoint()[2], 1.0e-10);
Assert.assertEquals(1.0, optimum.getPoint()[0], 1.0e-10);
Assert.assertEquals(2.0, optimum.getPoint()[1], 1.0e-10);
Assert.assertEquals(3.0, optimum.getPoint()[2], 1.0e-10);
}
public void testTwoSets() throws MathUserException {
@Test
public void testTwoSets() {
final double epsilon = 1.0e-7;
LinearProblem problem = new LinearProblem(new double[][] {
{ 2, 1, 0, 4, 0, 0 },
@ -203,17 +200,17 @@ extends TestCase {
RealPointValuePair optimum =
optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 0, 0, 0, 0, 0, 0 });
assertEquals( 3.0, optimum.getPoint()[0], 1.0e-10);
assertEquals( 4.0, optimum.getPoint()[1], 1.0e-10);
assertEquals(-1.0, optimum.getPoint()[2], 1.0e-10);
assertEquals(-2.0, optimum.getPoint()[3], 1.0e-10);
assertEquals( 1.0 + epsilon, optimum.getPoint()[4], 1.0e-10);
assertEquals( 1.0 - epsilon, optimum.getPoint()[5], 1.0e-10);
Assert.assertEquals( 3.0, optimum.getPoint()[0], 1.0e-10);
Assert.assertEquals( 4.0, optimum.getPoint()[1], 1.0e-10);
Assert.assertEquals(-1.0, optimum.getPoint()[2], 1.0e-10);
Assert.assertEquals(-2.0, optimum.getPoint()[3], 1.0e-10);
Assert.assertEquals( 1.0 + epsilon, optimum.getPoint()[4], 1.0e-10);
Assert.assertEquals( 1.0 - epsilon, optimum.getPoint()[5], 1.0e-10);
}
public void testNonInversible() throws MathUserException {
@Test
public void testNonInversible() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 1, 2, -3 },
{ 2, 1, 3 },
@ -225,10 +222,11 @@ extends TestCase {
optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-6, 1.0e-6));
RealPointValuePair optimum =
optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 0, 0, 0 });
assertTrue(optimum.getValue() > 0.5);
Assert.assertTrue(optimum.getValue() > 0.5);
}
public void testIllConditioned() throws MathUserException {
@Test
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 },
@ -239,16 +237,14 @@ extends TestCase {
new NonLinearConjugateGradientOptimizer(ConjugateGradientFormula.POLAK_RIBIERE);
optimizer.setMaxEvaluations(100);
optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-13, 1.0e-13));
BrentSolver solver = new BrentSolver();
solver.setAbsoluteAccuracy(1.0e-15);
solver.setRelativeAccuracy(1.0e-15);
BrentSolver solver = new BrentSolver(1e-15, 1e-15);
optimizer.setLineSearchSolver(solver);
RealPointValuePair optimum1 =
optimizer.optimize(problem1, GoalType.MINIMIZE, new double[] { 0, 1, 2, 3 });
assertEquals(1.0, optimum1.getPoint()[0], 1.0e-5);
assertEquals(1.0, optimum1.getPoint()[1], 1.0e-5);
assertEquals(1.0, optimum1.getPoint()[2], 1.0e-5);
assertEquals(1.0, optimum1.getPoint()[3], 1.0e-5);
Assert.assertEquals(1.0, optimum1.getPoint()[0], 1.0e-4);
Assert.assertEquals(1.0, optimum1.getPoint()[1], 1.0e-4);
Assert.assertEquals(1.0, optimum1.getPoint()[2], 1.0e-4);
Assert.assertEquals(1.0, optimum1.getPoint()[3], 1.0e-4);
LinearProblem problem2 = new LinearProblem(new double[][] {
{ 10.00, 7.00, 8.10, 7.20 },
@ -258,16 +254,15 @@ extends TestCase {
}, new double[] { 32, 23, 33, 31 });
RealPointValuePair optimum2 =
optimizer.optimize(problem2, GoalType.MINIMIZE, new double[] { 0, 1, 2, 3 });
assertEquals(-81.0, optimum2.getPoint()[0], 1.0e-1);
assertEquals(137.0, optimum2.getPoint()[1], 1.0e-1);
assertEquals(-34.0, optimum2.getPoint()[2], 1.0e-1);
assertEquals( 22.0, optimum2.getPoint()[3], 1.0e-1);
Assert.assertEquals(-81.0, optimum2.getPoint()[0], 1.0e-1);
Assert.assertEquals(137.0, optimum2.getPoint()[1], 1.0e-1);
Assert.assertEquals(-34.0, optimum2.getPoint()[2], 1.0e-1);
Assert.assertEquals( 22.0, optimum2.getPoint()[3], 1.0e-1);
}
public void testMoreEstimatedParametersSimple()
throws MathUserException {
@Test
public void testMoreEstimatedParametersSimple() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 3.0, 2.0, 0.0, 0.0 },
{ 0.0, 1.0, -1.0, 1.0 },
@ -280,12 +275,12 @@ extends TestCase {
optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-6, 1.0e-6));
RealPointValuePair optimum =
optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 7, 6, 5, 4 });
assertEquals(0, optimum.getValue(), 1.0e-10);
Assert.assertEquals(0, optimum.getValue(), 1.0e-10);
}
public void testMoreEstimatedParametersUnsorted()
throws MathUserException {
@Test
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 },
@ -299,10 +294,11 @@ extends TestCase {
optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-6, 1.0e-6));
RealPointValuePair optimum =
optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 2, 2, 2, 2, 2, 2 });
assertEquals(0, optimum.getValue(), 1.0e-10);
Assert.assertEquals(0, optimum.getValue(), 1.0e-10);
}
public void testRedundantEquations() throws MathUserException {
@Test
public void testRedundantEquations() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 1.0, 1.0 },
{ 1.0, -1.0 },
@ -315,12 +311,13 @@ extends TestCase {
optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-6, 1.0e-6));
RealPointValuePair optimum =
optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 1, 1 });
assertEquals(2.0, optimum.getPoint()[0], 1.0e-8);
assertEquals(1.0, optimum.getPoint()[1], 1.0e-8);
Assert.assertEquals(2.0, optimum.getPoint()[0], 1.0e-8);
Assert.assertEquals(1.0, optimum.getPoint()[1], 1.0e-8);
}
public void testInconsistentEquations() throws MathUserException {
@Test
public void testInconsistentEquations() {
LinearProblem problem = new LinearProblem(new double[][] {
{ 1.0, 1.0 },
{ 1.0, -1.0 },
@ -333,12 +330,13 @@ extends TestCase {
optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-6, 1.0e-6));
RealPointValuePair optimum =
optimizer.optimize(problem, GoalType.MINIMIZE, new double[] { 1, 1 });
assertTrue(optimum.getValue() > 0.1);
Assert.assertTrue(optimum.getValue() > 0.1);
}
public void testCircleFitting() throws MathUserException {
Circle circle = new Circle();
@Test
public void testCircleFitting() {
CircleScalar circle = new CircleScalar();
circle.addPoint( 30.0, 68.0);
circle.addPoint( 50.0, -6.0);
circle.addPoint(110.0, -20.0);
@ -348,16 +346,14 @@ extends TestCase {
new NonLinearConjugateGradientOptimizer(ConjugateGradientFormula.POLAK_RIBIERE);
optimizer.setMaxEvaluations(100);
optimizer.setConvergenceChecker(new SimpleScalarValueChecker(1.0e-30, 1.0e-30));
BrentSolver solver = new BrentSolver();
solver.setAbsoluteAccuracy(1.0e-13);
solver.setRelativeAccuracy(1.0e-15);
UnivariateRealSolver solver = new BrentSolver(1e-15, 1e-13);
optimizer.setLineSearchSolver(solver);
RealPointValuePair optimum =
optimizer.optimize(circle, GoalType.MINIMIZE, new double[] { 98.680, 47.345 });
Point2D.Double center = new Point2D.Double(optimum.getPointRef()[0], optimum.getPointRef()[1]);
assertEquals(69.960161753, circle.getRadius(center), 1.0e-8);
assertEquals(96.075902096, center.x, 1.0e-8);
assertEquals(48.135167894, center.y, 1.0e-8);
Assert.assertEquals(69.960161753, circle.getRadius(center), 1.0e-8);
Assert.assertEquals(96.075902096, center.x, 1.0e-8);
Assert.assertEquals(48.135167894, center.y, 1.0e-8);
}
private static class LinearProblem implements DifferentiableMultivariateRealFunction, Serializable {
@ -382,7 +378,7 @@ extends TestCase {
return p;
}
public double value(double[] variables) throws MathUserException {
public double value(double[] variables) {
double[] y = factors.operate(variables);
double sum = 0;
for (int i = 0; i < y.length; ++i) {
@ -409,86 +405,5 @@ extends TestCase {
}
};
}
}
private static class Circle implements DifferentiableMultivariateRealFunction, Serializable {
private static final long serialVersionUID = -4711170319243817874L;
private ArrayList<Point2D.Double> points;
public Circle() {
points = new ArrayList<Point2D.Double>();
}
public void addPoint(double px, double py) {
points.add(new Point2D.Double(px, py));
}
public double getRadius(Point2D.Double center) {
double r = 0;
for (Point2D.Double point : points) {
r += point.distance(center);
}
return r / points.size();
}
private double[] gradient(double[] point) {
// optimal radius
Point2D.Double center = new Point2D.Double(point[0], point[1]);
double radius = getRadius(center);
// gradient of the sum of squared residuals
double dJdX = 0;
double dJdY = 0;
for (Point2D.Double pk : points) {
double dk = pk.distance(center);
dJdX += (center.x - pk.x) * (dk - radius) / dk;
dJdY += (center.y - pk.y) * (dk - radius) / dk;
}
dJdX *= 2;
dJdY *= 2;
return new double[] { dJdX, dJdY };
}
public double value(double[] variables)
throws IllegalArgumentException, MathUserException {
Point2D.Double center = new Point2D.Double(variables[0], variables[1]);
double radius = getRadius(center);
double sum = 0;
for (Point2D.Double point : points) {
double di = point.distance(center) - radius;
sum += di * di;
}
return sum;
}
public MultivariateVectorialFunction gradient() {
return new MultivariateVectorialFunction() {
private static final long serialVersionUID = 3174909643301201710L;
public double[] value(double[] point) {
return gradient(point);
}
};
}
public MultivariateRealFunction partialDerivative(final int k) {
return new MultivariateRealFunction() {
private static final long serialVersionUID = 3073956364104833888L;
public double value(double[] point) {
return gradient(point)[k];
}
};
}
}
}