src/main/java/org/apache/commons/math/optimization/direct/DirectSearchOptimizer.java

/*
 * 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;
    }

}