SimulatedAnnealing.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;

import java.util.function.BiFunction;
import java.util.function.DoublePredicate;
import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.math4.legacy.stat.descriptive.moment.StandardDeviation;
import org.apache.commons.math4.legacy.optim.OptimizationData;
import org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv.Simplex;

/**
 * Simulated annealing setup.
 *
 * @since 4.0
 */
public class SimulatedAnnealing implements OptimizationData {
    /** Number of iterations at fixed temperature. */
    private final int epochDuration;
    /** Initial acceptance probability. */
    private final double startProbability;
    /** Final acceptance probability. */
    private final double endProbability;
    /** Cooling function. */
    private final CoolingSchedule coolingSchedule;
    /** RNG. */
    private final UniformRandomProvider rng;

    /**
     * @param epoch Number of iterations performed at fixed temperature.
     * @param startProb Initial acceptance probablility.
     * @param endProb Final acceptance probablility.
     * @param cooling Computes the temperature as a function of the initial
     * temperature and the epoch.
     * It is called for computing a new temperature after each cycle of
     * {@code epoch} iterations.
     * Simulated annealing <em>assumes</em> that the function decreases
     * monotically wrt the epoch (cf. {@link CoolingSchedule#decreasingExponential(double)
     * provided implementation}).
     * @param random Random number generator.
     * @throws IllegalArgumentException if {@code epoch < 1} or
     * {@code startProb} or {@code endProb} is outside the {@code [0, 1]}
     * interval.
     */
    public SimulatedAnnealing(int epoch,
                              double startProb,
                              double endProb,
                              CoolingSchedule cooling,
                              UniformRandomProvider random) {
        if (epoch < 1) {
            throw new IllegalArgumentException("Epoch out of range: " +
                                               epoch);
        }
        if (startProb < 0 ||
            startProb > 1) {
            throw new IllegalArgumentException("Initial acceptance probability out of range: " +
                                               startProb);
        }
        if (endProb < 0 ||
            endProb > 1) {
            throw new IllegalArgumentException("Final acceptance probability out of range: " +
                                               endProb);
        }
        if (endProb >= startProb) {
            throw new IllegalArgumentException("Final probability larger than initial probability");
        }

        epochDuration = epoch;
        startProbability = startProb;
        endProbability = endProb;
        coolingSchedule = cooling;
        rng = random;
    }

    /**
     * @return the epoch duration.
     */
    public int getEpochDuration() {
        return epochDuration;
    }

    /**
     * @return the acceptance probability at the beginning of the SA process.
     */
    public double getStartProbability() {
        return startProbability;
    }

    /**
     * @return the acceptance probability at the end of the SA process.
     */
    public double getEndProbability() {
        return endProbability;
    }

    /**
     * @return the cooling schedule.
     */
    public CoolingSchedule getCoolingSchedule() {
        return coolingSchedule;
    }

    /**
     * Specifies the cooling schedule.
     * It computes the current temperature as a function of two arguments:
     * <ol>
     *  <li>the previous temperature,</li>
     *  <li>the current simplex.</li>
     * </ol>
     */
    public interface CoolingSchedule extends BiFunction<Double, Simplex, Double> {
        /**
         * Power-law cooling scheme:
         * \[
         *   T_i = T_0 * f^i
         * \], where \( i \) is the current iteration.
         * <p>
         * Note: Simplex argument (of the returned function) is not used.
         *
         * @param f Factor by which the temperature is decreased.
         * @return the cooling schedule.
         */
        static CoolingSchedule decreasingExponential(final double f) {
            if (f <= 0 ||
                f >= 1) {
                throw new IllegalArgumentException("Factor out of range: " + f);
            }

            return (previousTemperature, simplex) -> f * previousTemperature;
        }

        /**
         * Aarst and van Laarhoven (1985) scheme:
         * \[
         *   T_{i + 1} = \frac{T_{i}}{1 + \frac{T_i \ln(1 + \delta)}{3 \sigma}}
         * \]
         * <p>
         * The simplex argument is used to compute the standard deviation
         * (\(\sigma\)) of all the vertices' objective function value.
         *
         * @param delta Trajectory parameter. Values smaller than 1 entail slow
         * convergence; values larger than 1 entail convergence to local optimum.
         * @return the cooling schedule.
         */
        static CoolingSchedule aarstAndVanLaarhoven(final double delta) {
            if (delta <= 0) {
                throw new IllegalArgumentException("Trajectory parameter out of range: " +
                                                   delta);
            }

            return (previousTemperature, simplex) -> {
                // Standard deviation of the values of the objective function.
                final StandardDeviation stddev = new StandardDeviation();
                for (int i = 0; i < simplex.getSize(); i++) {
                    stddev.increment(simplex.get(i).getValue());
                }
                final double sigma = stddev.getResult();

                final double a = previousTemperature * Math.log(1 + delta);
                final double b = 3 * sigma;
                return previousTemperature / (1 + a / b);
            };
        }
    }

    /**
     * Factory for the Metropolis check for accepting a worse state.
     * \( e^{-|\Delta E| / T} \geq \mathrm{rand}(0, 1) \).
     *
     * <p>
     * It is assumed that this check is performed <em>after</em> ensuring
     * that the alternate state is <em>worse</em> than the current best state.
     *
     * @param temperature Current temperature.
     * @return the acceptance test.
     */
    public DoublePredicate metropolis(final double temperature) {
        // Absolute value takes care of both minimization and maximization cases.
        return deltaE -> Math.exp(Math.abs(deltaE) / temperature) >= rng.nextDouble();
    }
}