commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/SimplexOptimizer.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.math4.legacy.optim.nonlinear.scalar.noderiv;

import java.util.Arrays;
import java.util.List;
import java.util.ArrayList;
import java.util.Comparator;
import java.util.Collections;
import java.util.Objects;
import java.util.function.UnaryOperator;
import java.util.function.IntSupplier;
import java.util.concurrent.CopyOnWriteArrayList;

import org.apache.commons.math4.legacy.core.MathArrays;
import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
import org.apache.commons.math4.legacy.exception.MathUnsupportedOperationException;
import org.apache.commons.math4.legacy.exception.MathInternalError;
import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
import org.apache.commons.math4.legacy.optim.ConvergenceChecker;
import org.apache.commons.math4.legacy.optim.OptimizationData;
import org.apache.commons.math4.legacy.optim.PointValuePair;
import org.apache.commons.math4.legacy.optim.SimpleValueChecker;
import org.apache.commons.math4.legacy.optim.InitialGuess;
import org.apache.commons.math4.legacy.optim.MaxEval;
import org.apache.commons.math4.legacy.optim.nonlinear.scalar.GoalType;
import org.apache.commons.math4.legacy.optim.nonlinear.scalar.MultivariateOptimizer;
import org.apache.commons.math4.legacy.optim.nonlinear.scalar.SimulatedAnnealing;
import org.apache.commons.math4.legacy.optim.nonlinear.scalar.PopulationSize;
import org.apache.commons.math4.legacy.optim.nonlinear.scalar.ObjectiveFunction;

/**
 * This class implements simplex-based direct search optimization.
 *
 * <p>
 * Direct search methods only use objective function values, they do
 * not 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>
 * 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>
 * In addition to those documented in
 * {@link MultivariateOptimizer#optimize(OptimizationData[]) MultivariateOptimizer},
 * an instance of this class will register the following data:
 * <ul>
 *  <li>{@link Simplex}</li>
 *  <li>{@link Simplex.TransformFactory}</li>
 *  <li>{@link SimulatedAnnealing}</li>
 *  <li>{@link PopulationSize}</li>
 * </ul>
 *
 * <p>
 * Each call to {@code optimize} will re-use the start configuration of
 * the current simplex 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 simplex must be re-initialized
 * to one with the appropriate dimensions.
 *
 * <p>
 * Convergence is considered achieved when <em>all</em> the simplex points
 * have converged.
 * <p>
 * Whenever {@link SimulatedAnnealing simulated annealing (SA)} is activated,
 * and the SA phase has completed, convergence has probably not been reached
 * yet; whenever it's the case, an additional (non-SA) search will be performed
 * (using the current best simplex point as a start point).
 * <p>
 * Additional "best list" searches can be requested through setting the
 * {@link PopulationSize} argument of the {@link #optimize(OptimizationData[])
 * optimize} method.
 *
 * <p>
 * This implementation does not directly support constrained optimization
 * with simple bounds.
 * The call to {@link #optimize(OptimizationData[]) optimize} will throw
 * {@link MathUnsupportedOperationException} if bounds are passed to it.
 *
 * @see NelderMeadTransform
 * @see MultiDirectionalTransform
 * @see HedarFukushimaTransform
 */
public class SimplexOptimizer extends MultivariateOptimizer {
    /** Default simplex side length ratio. */
    private static final double SIMPLEX_SIDE_RATIO = 1e-1;
    /** Simplex update function factory. */
    private Simplex.TransformFactory updateRule;
    /** Initial simplex. */
    private Simplex initialSimplex;
    /** Simulated annealing setup (optional). */
    private SimulatedAnnealing simulatedAnnealing = null;
    /** User-defined number of additional optimizations (optional). */
    private int populationSize = 0;
    /** Actual number of additional optimizations. */
    private int additionalSearch = 0;
    /** Callbacks. */
    private final List<Observer> callbacks = new CopyOnWriteArrayList<>();

    /**
     * @param checker Convergence checker.
     */
    public SimplexOptimizer(ConvergenceChecker<PointValuePair> checker) {
        super(checker);
    }

    /**
     * @param rel Relative threshold.
     * @param abs Absolute threshold.
     */
    public SimplexOptimizer(double rel,
                            double abs) {
        this(new SimpleValueChecker(rel, abs));
    }

    /**
     * Callback interface for updating caller's code with the current
     * state of the optimization.
     */
    @FunctionalInterface
    public interface Observer {
        /**
         * Method called after each modification of the {@code simplex}.
         *
         * @param simplex Current simplex.
         * @param isInit {@code true} at the start of a new search (either
         * "main" or "best list"), after the initial simplex's vertices
         * have been evaluated.
         * @param numEval Number of evaluations of the objective function.
         */
        void update(Simplex simplex,
                    boolean isInit,
                    int numEval);
    }

    /**
     * Register a callback.
     *
     * @param cb Callback.
     */
    public void addObserver(Observer cb) {
        Objects.requireNonNull(cb, "Callback");
        callbacks.add(cb);
    }

    /** {@inheritDoc} */
    @Override
    protected PointValuePair doOptimize() {
        checkParameters();

        // Indirect call to "computeObjectiveValue" in order to update the
        // evaluations counter.
        final MultivariateFunction evalFunc = this::computeObjectiveValue;

        final boolean isMinim = getGoalType() == GoalType.MINIMIZE;
        final Comparator<PointValuePair> comparator = (o1, o2) -> {
            final double v1 = o1.getValue();
            final double v2 = o2.getValue();
            return isMinim ? Double.compare(v1, v2) : Double.compare(v2, v1);
        };

        // Start points for additional search.
        final List<PointValuePair> bestList = new ArrayList<>();

        Simplex currentSimplex = initialSimplex.translate(getStartPoint()).evaluate(evalFunc, comparator);
        notifyObservers(currentSimplex, true);
        double temperature = Double.NaN; // Only used with simulated annealing.
        Simplex previousSimplex = null;

        if (simulatedAnnealing != null) {
            temperature =
                temperature(currentSimplex.get(0),
                            currentSimplex.get(currentSimplex.getDimension()),
                            simulatedAnnealing.getStartProbability());
        }

        while (true) {
            if (previousSimplex != null) { // Skip check at first iteration.
                if (hasConverged(previousSimplex, currentSimplex)) {
                    return currentSimplex.get(0);
                }
            }

            // We still need to search.
            previousSimplex = currentSimplex;

            if (simulatedAnnealing != null) {
                // Update current temperature.
                temperature =
                    simulatedAnnealing.getCoolingSchedule().apply(temperature,
                                                                  currentSimplex);

                final double endTemperature =
                    temperature(currentSimplex.get(0),
                                currentSimplex.get(currentSimplex.getDimension()),
                                simulatedAnnealing.getEndProbability());

                if (temperature < endTemperature) {
                    break;
                }

                final UnaryOperator<Simplex> update =
                    updateRule.create(evalFunc,
                                      comparator,
                                      simulatedAnnealing.metropolis(temperature));

                for (int i = 0; i < simulatedAnnealing.getEpochDuration(); i++) {
                    // Simplex is transformed (and observers are notified).
                    currentSimplex = applyUpdate(update,
                                                 currentSimplex,
                                                 evalFunc,
                                                 comparator);
                }
            } else {
                // No simulated annealing.
                final UnaryOperator<Simplex> update =
                    updateRule.create(evalFunc, comparator, null);

                // Simplex is transformed (and observers are notified).
                currentSimplex = applyUpdate(update,
                                             currentSimplex,
                                             evalFunc,
                                             comparator);
            }

            if (additionalSearch != 0) {
                // In "bestList", we must keep track of at least two points
                // in order to be able to compute the new initial simplex for
                // the additional search.
                final int max = Math.max(additionalSearch, 2);

                // Store best points.
                for (int i = 0; i < currentSimplex.getSize(); i++) {
                    keepIfBetter(currentSimplex.get(i),
                                 comparator,
                                 bestList,
                                 max);
                }
            }

            incrementIterationCount();
        }

        // No convergence.

        if (additionalSearch > 0) {
            // Additional optimizations.
            // Reference to counter in the "main" search in order to retrieve
            // the total number of evaluations in the "best list" search.
            final IntSupplier evalCount = () -> getEvaluations();

            return bestListSearch(evalFunc,
                                  comparator,
                                  bestList,
                                  evalCount);
        }

        throw new MathInternalError(); // Should never happen.
    }

    /**
     * Scans the list of (required and optional) optimization data that
     * characterize the problem.
     *
     * @param optData Optimization data.
     * The following data will be looked for:
     * <ul>
     *  <li>{@link Simplex}</li>
     *  <li>{@link Simplex.TransformFactory}</li>
     *  <li>{@link SimulatedAnnealing}</li>
     *  <li>{@link PopulationSize}</li>
     * </ul>
     */
    @Override
    protected void parseOptimizationData(OptimizationData... optData) {
        // Allow base class to register its own data.
        super.parseOptimizationData(optData);

        // The existing values (as set by the previous call) are reused
        // if not provided in the argument list.
        for (OptimizationData data : optData) {
            if (data instanceof Simplex) {
                initialSimplex = (Simplex) data;
            } else if (data instanceof Simplex.TransformFactory) {
                updateRule = (Simplex.TransformFactory) data;
            } else if (data instanceof SimulatedAnnealing) {
                simulatedAnnealing = (SimulatedAnnealing) data;
            } else if (data instanceof PopulationSize) {
                populationSize = ((PopulationSize) data).getPopulationSize();
            }
        }
    }

    /**
     * Detects whether the simplex has shrunk below the user-defined
     * tolerance.
     *
     * @param previous Simplex at previous iteration.
     * @param current Simplex at current iteration.
     * @return {@code true} if convergence is considered achieved.
     */
    private boolean hasConverged(Simplex previous,
                                 Simplex current) {
        final ConvergenceChecker<PointValuePair> checker = getConvergenceChecker();

        for (int i = 0; i < current.getSize(); i++) {
            final PointValuePair prev = previous.get(i);
            final PointValuePair curr = current.get(i);

            if (!checker.converged(getIterations(), prev, curr)) {
                return false;
            }
        }

        return true;
    }

    /**
     * @throws MathUnsupportedOperationException if bounds were passed to the
     * {@link #optimize(OptimizationData[]) optimize} method.
     * @throws NullPointerException if no initial simplex or no transform rule
     * was passed to the {@link #optimize(OptimizationData[]) optimize} method.
     * @throws IllegalArgumentException if {@link #populationSize} is negative.
     */
    private void checkParameters() {
        Objects.requireNonNull(updateRule, "Update rule");
        Objects.requireNonNull(initialSimplex, "Initial simplex");

        if (getLowerBound() != null ||
            getUpperBound() != null) {
            throw new MathUnsupportedOperationException(LocalizedFormats.CONSTRAINT);
        }

        if (populationSize < 0) {
            throw new IllegalArgumentException("Population size");
        }

        additionalSearch = simulatedAnnealing == null ?
            Math.max(0, populationSize) :
            Math.max(1, populationSize);
    }

    /**
     * Computes the temperature as a function of the acceptance probability
     * and the fitness difference between two of the simplex vertices (usually
     * the best and worst points).
     *
     * @param p1 Simplex point.
     * @param p2 Simplex point.
     * @param prob Acceptance probability.
     * @return the temperature.
     */
    private double temperature(PointValuePair p1,
                               PointValuePair p2,
                               double prob) {
        return -Math.abs(p1.getValue() - p2.getValue()) / Math.log(prob);
    }

    /**
     * Stores the given {@code candidate} if its fitness is better than
     * that of the last (assumed to be the worst) point in {@code list}.
     *
     * @param candidate Point to be stored.
     * @param comp Fitness comparator.
     * @param list Starting points (modified in-place).
     * @param max Maximum size of the {@code list}.
     */
    private static void keepIfBetter(PointValuePair candidate,
                                     Comparator<PointValuePair> comp,
                                     List<PointValuePair> list,
                                     int max) {
        final int listSize = list.size();
        final double[] candidatePoint = candidate.getPoint();
        if (listSize == 0) {
            list.add(candidate);
        } else if (listSize < max) {
            // List is not fully populated yet.
            for (PointValuePair p : list) {
                final double[] pPoint = p.getPoint();
                if (Arrays.equals(pPoint, candidatePoint)) {
                    // Point was already stored.
                    return;
                } else {
                    // Store candidate.
                    list.add(candidate);
                    return;
                }
            }
        } else {
            final int last = max - 1;
            if (comp.compare(candidate, list.get(last)) < 0) {
                for (PointValuePair p : list) {
                    final double[] pPoint = p.getPoint();
                    if (Arrays.equals(pPoint, candidatePoint)) {
                        // Point was already stored.
                        return;
                    }
                }

                // Store better candidate and reorder the list.
                list.set(last, candidate);
                Collections.sort(list, comp);
            }
        }
    }

    /**
     * Computes the smallest distance between the given {@code point}
     * and any of the other points in the {@code list}.
     *
     * @param point Point.
     * @param list List.
     * @return the smallest distance.
     */
    private static double shortestDistance(PointValuePair point,
                                           List<PointValuePair> list) {
        double minDist = Double.POSITIVE_INFINITY;

        final double[] p = point.getPoint();
        for (PointValuePair other : list) {
            final double[] pOther = other.getPoint();
            if (!Arrays.equals(p, pOther)) {
                final double dist = MathArrays.distance(p, pOther);
                if (dist < minDist) {
                    minDist = dist;
                }
            }
        }

        return minDist;
    }

    /**
     * Perform additional optimizations.
     *
     * @param evalFunc Objective function.
     * @param comp Fitness comparator.
     * @param bestList Best points encountered during the "main" search.
     * List is assumed to be ordered from best to worst.
     * @param evalCount Evaluation counter.
     * @return the optimum.
     */
    private PointValuePair bestListSearch(MultivariateFunction evalFunc,
                                          Comparator<PointValuePair> comp,
                                          List<PointValuePair> bestList,
                                          IntSupplier evalCount) {
        PointValuePair best = bestList.get(0); // Overall best result.

        // Additional local optimizations using each of the best
        // points visited during the main search.
        for (int i = 0; i < additionalSearch; i++) {
            final PointValuePair start = bestList.get(i);
            // Find shortest distance to the other points.
            final double dist = shortestDistance(start, bestList);
            final double[] init = start.getPoint();
            // Create smaller initial simplex.
            final Simplex simplex = Simplex.equalSidesAlongAxes(init.length,
                                                                SIMPLEX_SIDE_RATIO * dist);

            final PointValuePair r = directSearch(init,
                                                  simplex,
                                                  evalFunc,
                                                  getConvergenceChecker(),
                                                  getGoalType(),
                                                  callbacks,
                                                  evalCount);
            if (comp.compare(r, best) < 0) {
                best = r; // New overall best.
            }
        }

        return best;
    }

    /**
     * @param init Start point.
     * @param simplex Initial simplex.
     * @param eval Objective function.
     * Note: It is assumed that evaluations of this function are
     * incrementing the main counter.
     * @param checker Convergence checker.
     * @param goalType Whether to minimize or maximize the objective function.
     * @param cbList Callbacks.
     * @param evalCount Evaluation counter.
     * @return the optimum.
     */
    private static PointValuePair directSearch(double[] init,
                                               Simplex simplex,
                                               MultivariateFunction eval,
                                               ConvergenceChecker<PointValuePair> checker,
                                               GoalType goalType,
                                               List<Observer> cbList,
                                               final IntSupplier evalCount) {
        final SimplexOptimizer optim = new SimplexOptimizer(checker);

        for (Observer cOrig : cbList) {
            final SimplexOptimizer.Observer cNew = (spx, isInit, numEval) ->
                cOrig.update(spx, isInit, evalCount.getAsInt());

            optim.addObserver(cNew);
        }

        return optim.optimize(MaxEval.unlimited(),
                              new ObjectiveFunction(eval),
                              goalType,
                              new InitialGuess(init),
                              simplex,
                              new NelderMeadTransform());
    }

    /**
     * @param simplex Current simplex.
     * @param isInit Set to {@code true} at the start of a new search
     * (either "main" or "best list"), after the evaluation of the initial
     * simplex's vertices.
     */
    private void notifyObservers(Simplex simplex,
                                 boolean isInit) {
        for (Observer cb : callbacks) {
            cb.update(simplex,
                      isInit,
                      getEvaluations());
        }
    }

    /**
     * Applies the {@code update} to the given {@code simplex} (and notifies
     * observers).
     *
     * @param update Simplex transformation.
     * @param simplex Current simplex.
     * @param eval Objective function.
     * @param comp Fitness comparator.
     * @return the transformed simplex.
     */
    private Simplex applyUpdate(UnaryOperator<Simplex> update,
                                Simplex simplex,
                                MultivariateFunction eval,
                                Comparator<PointValuePair> comp) {
        final Simplex transformed = update.apply(simplex).evaluate(eval, comp);

        notifyObservers(transformed, false);

        return transformed;
    }
}