001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math4.legacy.optim.nonlinear.scalar;
018
019import java.util.function.BiFunction;
020import java.util.function.DoublePredicate;
021import org.apache.commons.rng.UniformRandomProvider;
022import org.apache.commons.math4.legacy.stat.descriptive.moment.StandardDeviation;
023import org.apache.commons.math4.legacy.optim.OptimizationData;
024import org.apache.commons.math4.legacy.optim.nonlinear.scalar.noderiv.Simplex;
025
026/**
027 * Simulated annealing setup.
028 *
029 * @since 4.0
030 */
031public class SimulatedAnnealing implements OptimizationData {
032    /** Number of iterations at fixed temperature. */
033    private final int epochDuration;
034    /** Initial acceptance probability. */
035    private final double startProbability;
036    /** Final acceptance probability. */
037    private final double endProbability;
038    /** Cooling function. */
039    private final CoolingSchedule coolingSchedule;
040    /** RNG. */
041    private final UniformRandomProvider rng;
042
043    /**
044     * @param epoch Number of iterations performed at fixed temperature.
045     * @param startProb Initial acceptance probablility.
046     * @param endProb Final acceptance probablility.
047     * @param cooling Computes the temperature as a function of the initial
048     * temperature and the epoch.
049     * It is called for computing a new temperature after each cycle of
050     * {@code epoch} iterations.
051     * Simulated annealing <em>assumes</em> that the function decreases
052     * monotically wrt the epoch (cf. {@link CoolingSchedule#decreasingExponential(double)
053     * provided implementation}).
054     * @param random Random number generator.
055     * @throws IllegalArgumentException if {@code epoch < 1} or
056     * {@code startProb} or {@code endProb} is outside the {@code [0, 1]}
057     * interval.
058     */
059    public SimulatedAnnealing(int epoch,
060                              double startProb,
061                              double endProb,
062                              CoolingSchedule cooling,
063                              UniformRandomProvider random) {
064        if (epoch < 1) {
065            throw new IllegalArgumentException("Epoch out of range: " +
066                                               epoch);
067        }
068        if (startProb < 0 ||
069            startProb > 1) {
070            throw new IllegalArgumentException("Initial acceptance probability out of range: " +
071                                               startProb);
072        }
073        if (endProb < 0 ||
074            endProb > 1) {
075            throw new IllegalArgumentException("Final acceptance probability out of range: " +
076                                               endProb);
077        }
078        if (endProb >= startProb) {
079            throw new IllegalArgumentException("Final probability larger than initial probability");
080        }
081
082        epochDuration = epoch;
083        startProbability = startProb;
084        endProbability = endProb;
085        coolingSchedule = cooling;
086        rng = random;
087    }
088
089    /**
090     * @return the epoch duration.
091     */
092    public int getEpochDuration() {
093        return epochDuration;
094    }
095
096    /**
097     * @return the acceptance probability at the beginning of the SA process.
098     */
099    public double getStartProbability() {
100        return startProbability;
101    }
102
103    /**
104     * @return the acceptance probability at the end of the SA process.
105     */
106    public double getEndProbability() {
107        return endProbability;
108    }
109
110    /**
111     * @return the cooling schedule.
112     */
113    public CoolingSchedule getCoolingSchedule() {
114        return coolingSchedule;
115    }
116
117    /**
118     * Specifies the cooling schedule.
119     * It computes the current temperature as a function of two arguments:
120     * <ol>
121     *  <li>the previous temperature,</li>
122     *  <li>the current simplex.</li>
123     * </ol>
124     */
125    public interface CoolingSchedule extends BiFunction<Double, Simplex, Double> {
126        /**
127         * Power-law cooling scheme:
128         * \[
129         *   T_i = T_0 * f^i
130         * \], where \( i \) is the current iteration.
131         * <p>
132         * Note: Simplex argument (of the returned function) is not used.
133         *
134         * @param f Factor by which the temperature is decreased.
135         * @return the cooling schedule.
136         */
137        static CoolingSchedule decreasingExponential(final double f) {
138            if (f <= 0 ||
139                f >= 1) {
140                throw new IllegalArgumentException("Factor out of range: " + f);
141            }
142
143            return (previousTemperature, simplex) -> f * previousTemperature;
144        }
145
146        /**
147         * Aarst and van Laarhoven (1985) scheme:
148         * \[
149         *   T_{i + 1} = \frac{T_{i}}{1 + \frac{T_i \ln(1 + \delta)}{3 \sigma}}
150         * \]
151         * <p>
152         * The simplex argument is used to compute the standard deviation
153         * (\(\sigma\)) of all the vertices' objective function value.
154         *
155         * @param delta Trajectory parameter. Values smaller than 1 entail slow
156         * convergence; values larger than 1 entail convergence to local optimum.
157         * @return the cooling schedule.
158         */
159        static CoolingSchedule aarstAndVanLaarhoven(final double delta) {
160            if (delta <= 0) {
161                throw new IllegalArgumentException("Trajectory parameter out of range: " +
162                                                   delta);
163            }
164
165            return (previousTemperature, simplex) -> {
166                // Standard deviation of the values of the objective function.
167                final StandardDeviation stddev = new StandardDeviation();
168                for (int i = 0; i < simplex.getSize(); i++) {
169                    stddev.increment(simplex.get(i).getValue());
170                }
171                final double sigma = stddev.getResult();
172
173                final double a = previousTemperature * Math.log(1 + delta);
174                final double b = 3 * sigma;
175                return previousTemperature / (1 + a / b);
176            };
177        }
178    }
179
180    /**
181     * Factory for the Metropolis check for accepting a worse state.
182     * \( e^{-|\Delta E| / T} \geq \mathrm{rand}(0, 1) \).
183     *
184     * <p>
185     * It is assumed that this check is performed <em>after</em> ensuring
186     * that the alternate state is <em>worse</em> than the current best state.
187     *
188     * @param temperature Current temperature.
189     * @return the acceptance test.
190     */
191    public DoublePredicate metropolis(final double temperature) {
192        // Absolute value takes care of both minimization and maximization cases.
193        return deltaE -> Math.exp(Math.abs(deltaE) / temperature) >= rng.nextDouble();
194    }
195}