SimulatedAnnealing.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math4.legacy.optim.nonlinear.scalar;

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

  24. /**
  25.  * Simulated annealing setup.
  26.  *
  27.  * @since 4.0
  28.  */
  29. public class SimulatedAnnealing implements OptimizationData {
  30.     /** Number of iterations at fixed temperature. */
  31.     private final int epochDuration;
  32.     /** Initial acceptance probability. */
  33.     private final double startProbability;
  34.     /** Final acceptance probability. */
  35.     private final double endProbability;
  36.     /** Cooling function. */
  37.     private final CoolingSchedule coolingSchedule;
  38.     /** RNG. */
  39.     private final UniformRandomProvider rng;

  40.     /**
  41.      * @param epoch Number of iterations performed at fixed temperature.
  42.      * @param startProb Initial acceptance probablility.
  43.      * @param endProb Final acceptance probablility.
  44.      * @param cooling Computes the temperature as a function of the initial
  45.      * temperature and the epoch.
  46.      * It is called for computing a new temperature after each cycle of
  47.      * {@code epoch} iterations.
  48.      * Simulated annealing <em>assumes</em> that the function decreases
  49.      * monotically wrt the epoch (cf. {@link CoolingSchedule#decreasingExponential(double)
  50.      * provided implementation}).
  51.      * @param random Random number generator.
  52.      * @throws IllegalArgumentException if {@code epoch < 1} or
  53.      * {@code startProb} or {@code endProb} is outside the {@code [0, 1]}
  54.      * interval.
  55.      */
  56.     public SimulatedAnnealing(int epoch,
  57.                               double startProb,
  58.                               double endProb,
  59.                               CoolingSchedule cooling,
  60.                               UniformRandomProvider random) {
  61.         if (epoch < 1) {
  62.             throw new IllegalArgumentException("Epoch out of range: " +
  63.                                                epoch);
  64.         }
  65.         if (startProb < 0 ||
  66.             startProb > 1) {
  67.             throw new IllegalArgumentException("Initial acceptance probability out of range: " +
  68.                                                startProb);
  69.         }
  70.         if (endProb < 0 ||
  71.             endProb > 1) {
  72.             throw new IllegalArgumentException("Final acceptance probability out of range: " +
  73.                                                endProb);
  74.         }
  75.         if (endProb >= startProb) {
  76.             throw new IllegalArgumentException("Final probability larger than initial probability");
  77.         }

  78.         epochDuration = epoch;
  79.         startProbability = startProb;
  80.         endProbability = endProb;
  81.         coolingSchedule = cooling;
  82.         rng = random;
  83.     }

  84.     /**
  85.      * @return the epoch duration.
  86.      */
  87.     public int getEpochDuration() {
  88.         return epochDuration;
  89.     }

  90.     /**
  91.      * @return the acceptance probability at the beginning of the SA process.
  92.      */
  93.     public double getStartProbability() {
  94.         return startProbability;
  95.     }

  96.     /**
  97.      * @return the acceptance probability at the end of the SA process.
  98.      */
  99.     public double getEndProbability() {
  100.         return endProbability;
  101.     }

  102.     /**
  103.      * @return the cooling schedule.
  104.      */
  105.     public CoolingSchedule getCoolingSchedule() {
  106.         return coolingSchedule;
  107.     }

  108.     /**
  109.      * Specifies the cooling schedule.
  110.      * It computes the current temperature as a function of two arguments:
  111.      * <ol>
  112.      *  <li>the previous temperature,</li>
  113.      *  <li>the current simplex.</li>
  114.      * </ol>
  115.      */
  116.     public interface CoolingSchedule extends BiFunction<Double, Simplex, Double> {
  117.         /**
  118.          * Power-law cooling scheme:
  119.          * \[
  120.          *   T_i = T_0 * f^i
  121.          * \], where \( i \) is the current iteration.
  122.          * <p>
  123.          * Note: Simplex argument (of the returned function) is not used.
  124.          *
  125.          * @param f Factor by which the temperature is decreased.
  126.          * @return the cooling schedule.
  127.          */
  128.         static CoolingSchedule decreasingExponential(final double f) {
  129.             if (f <= 0 ||
  130.                 f >= 1) {
  131.                 throw new IllegalArgumentException("Factor out of range: " + f);
  132.             }

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

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

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

  160.                 final double a = previousTemperature * Math.log(1 + delta);
  161.                 final double b = 3 * sigma;
  162.                 return previousTemperature / (1 + a / b);
  163.             };
  164.         }
  165.     }

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