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}