/*
* 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.direct;
import java.util.Arrays;
import java.util.Comparator;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.MaxEvaluationsExceededException;
import org.apache.commons.math.MaxIterationsExceededException;
import org.apache.commons.math.analysis.MultivariateRealFunction;
import org.apache.commons.math.exception.util.LocalizedFormats;
import org.apache.commons.math.optimization.GoalType;
import org.apache.commons.math.optimization.MultivariateRealOptimizer;
import org.apache.commons.math.optimization.OptimizationException;
import org.apache.commons.math.optimization.RealConvergenceChecker;
import org.apache.commons.math.optimization.RealPointValuePair;
import org.apache.commons.math.optimization.SimpleScalarValueChecker;
/**
* This class implements simplex-based direct search optimization
* algorithms.
*
* <p>Direct search methods only use objective function values, they don't
* need derivatives and don't either try to compute approximation of
* the derivatives. According to a 1996 paper by Margaret H. Wright
* (<a href="http://cm.bell-labs.com/cm/cs/doc/96/4-02.ps.gz">Direct
* Search Methods: Once Scorned, Now Respectable</a>), they are used
* when either the computation of the derivative is impossible (noisy
* functions, unpredictable discontinuities) or difficult (complexity,
* computation cost). In the first cases, rather than an optimum, a
* <em>not too bad</em> point is desired. In the latter cases, an
* optimum is desired but cannot be reasonably found. In all cases
* direct search methods can be useful.</p>
*
* <p>Simplex-based direct search methods are based on comparison of
* the objective function values at the vertices of a simplex (which is a
* set of n+1 points in dimension n) that is updated by the algorithms
* steps.<p>
*
* <p>The initial configuration of the simplex can be set using either
* {@link #setStartConfiguration(double[])} or {@link
* #setStartConfiguration(double[][])}. If neither method has been called
* before optimization is attempted, an explicit call to the first method
* with all steps set to +1 is triggered, thus building a default
* configuration from a unit hypercube. Each call to {@link
* #optimize(MultivariateRealFunction, GoalType, double[]) optimize} will reuse
* the current start configuration and move it such that its first vertex
* is at the provided start point of the optimization. If the {@code optimize}
* method is called to solve a different problem and the number of parameters
* change, the start configuration will be reset to a default one with the
* appropriate dimensions.</p>
*
* <p>If {@link #setConvergenceChecker(RealConvergenceChecker)} is not called,
* a default {@link SimpleScalarValueChecker} is used.</p>
*
* <p>Convergence is checked by providing the <em>worst</em> points of
* previous and current simplex to the convergence checker, not the best ones.</p>
*
* <p>This class is the base class performing the boilerplate simplex
* initialization and handling. The simplex update by itself is
* performed by the derived classes according to the implemented
* algorithms.</p>
*
* implements MultivariateRealOptimizer since 2.0
*
* @see MultivariateRealFunction
* @see NelderMead
* @see MultiDirectional
* @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 févr. 2011) $
* @since 1.2
*/
public abstract class DirectSearchOptimizer implements MultivariateRealOptimizer {
/** Simplex. */
protected RealPointValuePair[] simplex;
/** Objective function. */
private MultivariateRealFunction f;
/** Convergence checker. */
private RealConvergenceChecker checker;
/** Maximal number of iterations allowed. */
private int maxIterations;
/** Number of iterations already performed. */
private int iterations;
/** Maximal number of evaluations allowed. */
private int maxEvaluations;
/** Number of evaluations already performed. */
private int evaluations;
/** Start simplex configuration. */
private double[][] startConfiguration;
/** Simple constructor.
*/
protected DirectSearchOptimizer() {
setConvergenceChecker(new SimpleScalarValueChecker());
setMaxIterations(Integer.MAX_VALUE);
setMaxEvaluations(Integer.MAX_VALUE);
}
/** Set start configuration for simplex.
* <p>The start configuration for simplex is built from a box parallel to
* the canonical axes of the space. The simplex is the subset of vertices
* of a box parallel to the canonical axes. It is built as the path followed
* while traveling from one vertex of the box to the diagonally opposite
* vertex moving only along the box edges. The first vertex of the box will
* be located at the start point of the optimization.</p>
* <p>As an example, in dimension 3 a simplex has 4 vertices. Setting the
* steps to (1, 10, 2) and the start point to (1, 1, 1) would imply the
* start simplex would be: { (1, 1, 1), (2, 1, 1), (2, 11, 1), (2, 11, 3) }.
* The first vertex would be set to the start point at (1, 1, 1) and the
* last vertex would be set to the diagonally opposite vertex at (2, 11, 3).</p>
* @param steps steps along the canonical axes representing box edges,
* they may be negative but not null
* @exception IllegalArgumentException if one step is null
*/
public void setStartConfiguration(final double[] steps)
throws IllegalArgumentException {
// only the relative position of the n final vertices with respect
// to the first one are stored
final int n = steps.length;
startConfiguration = new double[n][n];
for (int i = 0; i < n; ++i) {
final double[] vertexI = startConfiguration[i];
for (int j = 0; j < i + 1; ++j) {
if (steps[j] == 0.0) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.EQUAL_VERTICES_IN_SIMPLEX, j, j + 1);
}
System.arraycopy(steps, 0, vertexI, 0, j + 1);
}
}
}
/** Set start configuration for simplex.
* <p>The real initial simplex will be set up by moving the reference
* simplex such that its first point is located at the start point of the
* optimization.</p>
* @param referenceSimplex reference simplex
* @exception IllegalArgumentException if the reference simplex does not
* contain at least one point, or if there is a dimension mismatch
* in the reference simplex or if one of its vertices is duplicated
*/
public void setStartConfiguration(final double[][] referenceSimplex)
throws IllegalArgumentException {
// only the relative position of the n final vertices with respect
// to the first one are stored
final int n = referenceSimplex.length - 1;
if (n < 0) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.SIMPLEX_NEED_ONE_POINT);
}
startConfiguration = new double[n][n];
final double[] ref0 = referenceSimplex[0];
// vertices loop
for (int i = 0; i < n + 1; ++i) {
final double[] refI = referenceSimplex[i];
// safety checks
if (refI.length != n) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, refI.length, n);
}
for (int j = 0; j < i; ++j) {
final double[] refJ = referenceSimplex[j];
boolean allEquals = true;
for (int k = 0; k < n; ++k) {
if (refI[k] != refJ[k]) {
allEquals = false;
break;
}
}
if (allEquals) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.EQUAL_VERTICES_IN_SIMPLEX, i, j);
}
}
// store vertex i position relative to vertex 0 position
if (i > 0) {
final double[] confI = startConfiguration[i - 1];
for (int k = 0; k < n; ++k) {
confI[k] = refI[k] - ref0[k];
}
}
}
}
/** {@inheritDoc} */
public void setMaxIterations(int maxIterations) {
this.maxIterations = maxIterations;
}
/** {@inheritDoc} */
public int getMaxIterations() {
return maxIterations;
}
/** {@inheritDoc} */
public void setMaxEvaluations(int maxEvaluations) {
this.maxEvaluations = maxEvaluations;
}
/** {@inheritDoc} */
public int getMaxEvaluations() {
return maxEvaluations;
}
/** {@inheritDoc} */
public int getIterations() {
return iterations;
}
/** {@inheritDoc} */
public int getEvaluations() {
return evaluations;
}
/** {@inheritDoc} */
public void setConvergenceChecker(RealConvergenceChecker convergenceChecker) {
this.checker = convergenceChecker;
}
/** {@inheritDoc} */
public RealConvergenceChecker getConvergenceChecker() {
return checker;
}
/** {@inheritDoc} */
public RealPointValuePair optimize(final MultivariateRealFunction function,
final GoalType goalType,
final double[] startPoint)
throws FunctionEvaluationException, OptimizationException, IllegalArgumentException {
if ((startConfiguration == null) ||
(startConfiguration.length != startPoint.length)) {
// no initial configuration has been set up for simplex
// build a default one from a unit hypercube
final double[] unit = new double[startPoint.length];
Arrays.fill(unit, 1.0);
setStartConfiguration(unit);
}
this.f = function;
final Comparator<RealPointValuePair> comparator =
new Comparator<RealPointValuePair>() {
public int compare(final RealPointValuePair o1,
final RealPointValuePair o2) {
final double v1 = o1.getValue();
final double v2 = o2.getValue();
return (goalType == GoalType.MINIMIZE) ?
Double.compare(v1, v2) : Double.compare(v2, v1);
}
};
// initialize search
iterations = 0;
evaluations = 0;
buildSimplex(startPoint);
evaluateSimplex(comparator);
RealPointValuePair[] previous = new RealPointValuePair[simplex.length];
while (true) {
if (iterations > 0) {
boolean converged = true;
for (int i = 0; i < simplex.length; ++i) {
converged &= checker.converged(iterations, previous[i], simplex[i]);
}
if (converged) {
// we have found an optimum
return simplex[0];
}
}
// we still need to search
System.arraycopy(simplex, 0, previous, 0, simplex.length);
iterateSimplex(comparator);
}
}
/** Increment the iterations counter by 1.
* @exception OptimizationException if the maximal number
* of iterations is exceeded
*/
protected void incrementIterationsCounter()
throws OptimizationException {
if (++iterations > maxIterations) {
throw new OptimizationException(new MaxIterationsExceededException(maxIterations));
}
}
/** Compute the next simplex of the algorithm.
* @param comparator comparator to use to sort simplex vertices from best to worst
* @exception FunctionEvaluationException if the function cannot be evaluated at
* some point
* @exception OptimizationException if the algorithm fails to converge
* @exception IllegalArgumentException if the start point dimension is wrong
*/
protected abstract void iterateSimplex(final Comparator<RealPointValuePair> comparator)
throws FunctionEvaluationException, OptimizationException, IllegalArgumentException;
/** Evaluate the objective function on one point.
* <p>A side effect of this method is to count the number of
* function evaluations</p>
* @param x point on which the objective function should be evaluated
* @return objective function value at the given point
* @exception FunctionEvaluationException if no value can be computed for the
* parameters or if the maximal number of evaluations is exceeded
* @exception IllegalArgumentException if the start point dimension is wrong
*/
protected double evaluate(final double[] x)
throws FunctionEvaluationException, IllegalArgumentException {
if (++evaluations > maxEvaluations) {
throw new FunctionEvaluationException(new MaxEvaluationsExceededException(maxEvaluations), x);
}
return f.value(x);
}
/** Build an initial simplex.
* @param startPoint the start point for optimization
* @exception IllegalArgumentException if the start point does not match
* simplex dimension
*/
private void buildSimplex(final double[] startPoint)
throws IllegalArgumentException {
final int n = startPoint.length;
if (n != startConfiguration.length) {
throw MathRuntimeException.createIllegalArgumentException(
LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, n, startConfiguration.length);
}
// set first vertex
simplex = new RealPointValuePair[n + 1];
simplex[0] = new RealPointValuePair(startPoint, Double.NaN);
// set remaining vertices
for (int i = 0; i < n; ++i) {
final double[] confI = startConfiguration[i];
final double[] vertexI = new double[n];
for (int k = 0; k < n; ++k) {
vertexI[k] = startPoint[k] + confI[k];
}
simplex[i + 1] = new RealPointValuePair(vertexI, Double.NaN);
}
}
/** Evaluate all the non-evaluated points of the simplex.
* @param comparator comparator to use to sort simplex vertices from best to worst
* @exception FunctionEvaluationException if no value can be computed for the parameters
* @exception OptimizationException if the maximal number of evaluations is exceeded
*/
protected void evaluateSimplex(final Comparator<RealPointValuePair> comparator)
throws FunctionEvaluationException, OptimizationException {
// evaluate the objective function at all non-evaluated simplex points
for (int i = 0; i < simplex.length; ++i) {
final RealPointValuePair vertex = simplex[i];
final double[] point = vertex.getPointRef();
if (Double.isNaN(vertex.getValue())) {
simplex[i] = new RealPointValuePair(point, evaluate(point), false);
}
}
// sort the simplex from best to worst
Arrays.sort(simplex, comparator);
}
/** Replace the worst point of the simplex by a new point.
* @param pointValuePair point to insert
* @param comparator comparator to use to sort simplex vertices from best to worst
*/
protected void replaceWorstPoint(RealPointValuePair pointValuePair,
final Comparator<RealPointValuePair> comparator) {
int n = simplex.length - 1;
for (int i = 0; i < n; ++i) {
if (comparator.compare(simplex[i], pointValuePair) > 0) {
RealPointValuePair tmp = simplex[i];
simplex[i] = pointValuePair;
pointValuePair = tmp;
}
}
simplex[n] = pointValuePair;
}
}