/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.optim.nonlinear.scalar.noderiv;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.MathArrays;
import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.exception.NotStrictlyPositiveException;
import org.apache.commons.math3.exception.MathUnsupportedOperationException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
import org.apache.commons.math3.optim.MaxEval;
import org.apache.commons.math3.optim.PointValuePair;
import org.apache.commons.math3.optim.ConvergenceChecker;
import org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer;
import org.apache.commons.math3.optim.univariate.BracketFinder;
import org.apache.commons.math3.optim.univariate.BrentOptimizer;
import org.apache.commons.math3.optim.univariate.UnivariatePointValuePair;
import org.apache.commons.math3.optim.univariate.SimpleUnivariateValueChecker;
import org.apache.commons.math3.optim.univariate.SearchInterval;
import org.apache.commons.math3.optim.univariate.UnivariateObjectiveFunction;
/**
* Powell's algorithm.
* This code is translated and adapted from the Python version of this
* algorithm (as implemented in module {@code optimize.py} v0.5 of
* <em>SciPy</em>).
* <br/>
* The default stopping criterion is based on the differences of the
* function value between two successive iterations. It is however possible
* to define a custom convergence checker that might terminate the algorithm
* earlier.
* <br/>
* The internal line search optimizer is a {@link BrentOptimizer} with a
* convergence checker set to {@link SimpleUnivariateValueChecker}.
* <br/>
* Constraints are not supported: the call to
* {@link #optimize(OptimizationData[]) optimize} will throw
* {@link MathUnsupportedOperationException} if bounds are passed to it.
* In order to impose simple constraints, the objective function must be
* wrapped in an adapter like
* {@link org.apache.commons.math3.optim.nonlinear.scalar.MultivariateFunctionMappingAdapter
* MultivariateFunctionMappingAdapter} or
* {@link org.apache.commons.math3.optim.nonlinear.scalar.MultivariateFunctionPenaltyAdapter
* MultivariateFunctionPenaltyAdapter}.
*
* @version $Id: PowellOptimizer.java 1462503 2013-03-29 15:48:27Z luc $
* @since 2.2
*/
public class PowellOptimizer
extends MultivariateOptimizer {
/**
* Minimum relative tolerance.
*/
private static final double MIN_RELATIVE_TOLERANCE = 2 * FastMath.ulp(1d);
/**
* Relative threshold.
*/
private final double relativeThreshold;
/**
* Absolute threshold.
*/
private final double absoluteThreshold;
/**
* Line search.
*/
private final LineSearch line;
/**
* This constructor allows to specify a user-defined convergence checker,
* in addition to the parameters that control the default convergence
* checking procedure.
* <br/>
* The internal line search tolerances are set to the square-root of their
* corresponding value in the multivariate optimizer.
*
* @param rel Relative threshold.
* @param abs Absolute threshold.
* @param checker Convergence checker.
* @throws NotStrictlyPositiveException if {@code abs <= 0}.
* @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
*/
public PowellOptimizer(double rel,
double abs,
ConvergenceChecker<PointValuePair> checker) {
this(rel, abs, FastMath.sqrt(rel), FastMath.sqrt(abs), checker);
}
/**
* This constructor allows to specify a user-defined convergence checker,
* in addition to the parameters that control the default convergence
* checking procedure and the line search tolerances.
*
* @param rel Relative threshold for this optimizer.
* @param abs Absolute threshold for this optimizer.
* @param lineRel Relative threshold for the internal line search optimizer.
* @param lineAbs Absolute threshold for the internal line search optimizer.
* @param checker Convergence checker.
* @throws NotStrictlyPositiveException if {@code abs <= 0}.
* @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
*/
public PowellOptimizer(double rel,
double abs,
double lineRel,
double lineAbs,
ConvergenceChecker<PointValuePair> checker) {
super(checker);
if (rel < MIN_RELATIVE_TOLERANCE) {
throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true);
}
if (abs <= 0) {
throw new NotStrictlyPositiveException(abs);
}
relativeThreshold = rel;
absoluteThreshold = abs;
// Create the line search optimizer.
line = new LineSearch(lineRel,
lineAbs);
}
/**
* The parameters control the default convergence checking procedure.
* <br/>
* The internal line search tolerances are set to the square-root of their
* corresponding value in the multivariate optimizer.
*
* @param rel Relative threshold.
* @param abs Absolute threshold.
* @throws NotStrictlyPositiveException if {@code abs <= 0}.
* @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
*/
public PowellOptimizer(double rel,
double abs) {
this(rel, abs, null);
}
/**
* Builds an instance with the default convergence checking procedure.
*
* @param rel Relative threshold.
* @param abs Absolute threshold.
* @param lineRel Relative threshold for the internal line search optimizer.
* @param lineAbs Absolute threshold for the internal line search optimizer.
* @throws NotStrictlyPositiveException if {@code abs <= 0}.
* @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
*/
public PowellOptimizer(double rel,
double abs,
double lineRel,
double lineAbs) {
this(rel, abs, lineRel, lineAbs, null);
}
/** {@inheritDoc} */
@Override
protected PointValuePair doOptimize() {
checkParameters();
final GoalType goal = getGoalType();
final double[] guess = getStartPoint();
final int n = guess.length;
final double[][] direc = new double[n][n];
for (int i = 0; i < n; i++) {
direc[i][i] = 1;
}
final ConvergenceChecker<PointValuePair> checker
= getConvergenceChecker();
double[] x = guess;
double fVal = computeObjectiveValue(x);
double[] x1 = x.clone();
while (true) {
incrementIterationCount();
double fX = fVal;
double fX2 = 0;
double delta = 0;
int bigInd = 0;
double alphaMin = 0;
for (int i = 0; i < n; i++) {
final double[] d = MathArrays.copyOf(direc[i]);
fX2 = fVal;
final UnivariatePointValuePair optimum = line.search(x, d);
fVal = optimum.getValue();
alphaMin = optimum.getPoint();
final double[][] result = newPointAndDirection(x, d, alphaMin);
x = result[0];
if ((fX2 - fVal) > delta) {
delta = fX2 - fVal;
bigInd = i;
}
}
// Default convergence check.
boolean stop = 2 * (fX - fVal) <=
(relativeThreshold * (FastMath.abs(fX) + FastMath.abs(fVal)) +
absoluteThreshold);
final PointValuePair previous = new PointValuePair(x1, fX);
final PointValuePair current = new PointValuePair(x, fVal);
if (!stop && checker != null) { // User-defined stopping criteria.
stop = checker.converged(getIterations(), previous, current);
}
if (stop) {
if (goal == GoalType.MINIMIZE) {
return (fVal < fX) ? current : previous;
} else {
return (fVal > fX) ? current : previous;
}
}
final double[] d = new double[n];
final double[] x2 = new double[n];
for (int i = 0; i < n; i++) {
d[i] = x[i] - x1[i];
x2[i] = 2 * x[i] - x1[i];
}
x1 = x.clone();
fX2 = computeObjectiveValue(x2);
if (fX > fX2) {
double t = 2 * (fX + fX2 - 2 * fVal);
double temp = fX - fVal - delta;
t *= temp * temp;
temp = fX - fX2;
t -= delta * temp * temp;
if (t < 0.0) {
final UnivariatePointValuePair optimum = line.search(x, d);
fVal = optimum.getValue();
alphaMin = optimum.getPoint();
final double[][] result = newPointAndDirection(x, d, alphaMin);
x = result[0];
final int lastInd = n - 1;
direc[bigInd] = direc[lastInd];
direc[lastInd] = result[1];
}
}
}
}
/**
* Compute a new point (in the original space) and a new direction
* vector, resulting from the line search.
*
* @param p Point used in the line search.
* @param d Direction used in the line search.
* @param optimum Optimum found by the line search.
* @return a 2-element array containing the new point (at index 0) and
* the new direction (at index 1).
*/
private double[][] newPointAndDirection(double[] p,
double[] d,
double optimum) {
final int n = p.length;
final double[] nP = new double[n];
final double[] nD = new double[n];
for (int i = 0; i < n; i++) {
nD[i] = d[i] * optimum;
nP[i] = p[i] + nD[i];
}
final double[][] result = new double[2][];
result[0] = nP;
result[1] = nD;
return result;
}
/**
* Class for finding the minimum of the objective function along a given
* direction.
*/
private class LineSearch extends BrentOptimizer {
/**
* Value that will pass the precondition check for {@link BrentOptimizer}
* but will not pass the convergence check, so that the custom checker
* will always decide when to stop the line search.
*/
private static final double REL_TOL_UNUSED = 1e-15;
/**
* Value that will pass the precondition check for {@link BrentOptimizer}
* but will not pass the convergence check, so that the custom checker
* will always decide when to stop the line search.
*/
private static final double ABS_TOL_UNUSED = Double.MIN_VALUE;
/**
* Automatic bracketing.
*/
private final BracketFinder bracket = new BracketFinder();
/**
* The "BrentOptimizer" default stopping criterion uses the tolerances
* to check the domain (point) values, not the function values.
* We thus create a custom checker to use function values.
*
* @param rel Relative threshold.
* @param abs Absolute threshold.
*/
LineSearch(double rel,
double abs) {
super(REL_TOL_UNUSED,
ABS_TOL_UNUSED,
new SimpleUnivariateValueChecker(rel, abs));
}
/**
* Find the minimum of the function {@code f(p + alpha * d)}.
*
* @param p Starting point.
* @param d Search direction.
* @return the optimum.
* @throws org.apache.commons.math3.exception.TooManyEvaluationsException
* if the number of evaluations is exceeded.
*/
public UnivariatePointValuePair search(final double[] p, final double[] d) {
final int n = p.length;
final UnivariateFunction f = new UnivariateFunction() {
public double value(double alpha) {
final double[] x = new double[n];
for (int i = 0; i < n; i++) {
x[i] = p[i] + alpha * d[i];
}
final double obj = PowellOptimizer.this.computeObjectiveValue(x);
return obj;
}
};
final GoalType goal = PowellOptimizer.this.getGoalType();
bracket.search(f, goal, 0, 1);
// Passing "MAX_VALUE" as a dummy value because it is the enclosing
// class that counts the number of evaluations (and will eventually
// generate the exception).
return optimize(new MaxEval(Integer.MAX_VALUE),
new UnivariateObjectiveFunction(f),
goal,
new SearchInterval(bracket.getLo(),
bracket.getHi(),
bracket.getMid()));
}
}
/**
* @throws MathUnsupportedOperationException if bounds were passed to the
* {@link #optimize(OptimizationData[]) optimize} method.
*/
private void checkParameters() {
if (getLowerBound() != null ||
getUpperBound() != null) {
throw new MathUnsupportedOperationException(LocalizedFormats.CONSTRAINT);
}
}
}