001package org.cpsolver.ifs.algorithms;
002
003import java.text.DecimalFormat;
004
005import org.cpsolver.ifs.assignment.Assignment;
006import org.cpsolver.ifs.heuristics.NeighbourSelection;
007import org.cpsolver.ifs.model.Model;
008import org.cpsolver.ifs.model.Neighbour;
009import org.cpsolver.ifs.model.Value;
010import org.cpsolver.ifs.model.Variable;
011import org.cpsolver.ifs.solution.Solution;
012import org.cpsolver.ifs.util.DataProperties;
013import org.cpsolver.ifs.util.JProf;
014import org.cpsolver.ifs.util.ToolBox;
015
016
017/**
018 * Simulated annealing. In each iteration, one of the given neighbourhoods is selected first,
019 * then a neighbour is generated and it is accepted with probability
020 * {@link SimulatedAnnealingContext#prob(double)}. The search is guided by the
021 * temperature, which starts at <i>SimulatedAnnealing.InitialTemperature</i>.
022 * After each <i>SimulatedAnnealing.TemperatureLength</i> iterations, the
023 * temperature is reduced by <i>SimulatedAnnealing.CoolingRate</i>. If there was
024 * no improvement in the past <i>SimulatedAnnealing.ReheatLengthCoef *
025 * SimulatedAnnealing.TemperatureLength</i> iterations, the temperature is
026 * increased by <i>SimulatedAnnealing.ReheatRate</i>. If there was no
027 * improvement in the past <i>SimulatedAnnealing.RestoreBestLengthCoef *
028 * SimulatedAnnealing.TemperatureLength</i> iterations, the best ever found
029 * solution is restored. <br>
030 * <br>
031 * If <i>SimulatedAnnealing.StochasticHC</i> is true, the acceptance probability
032 * is computed using stochastic hill climbing criterion, i.e.,
033 * <code>1.0 / (1.0 + Math.exp(value/temperature))</code>, otherwise it is
034 * cumputed using simlated annealing criterion, i.e.,
035 * <code>(value&lt;=0.0?1.0:Math.exp(-value/temperature))</code>. If
036 * <i>SimulatedAnnealing.RelativeAcceptance</i> neighbour value
037 * {@link Neighbour#value(Assignment)} is taken as the value of the selected
038 * neighbour (difference between the new and the current solution, if the
039 * neighbour is accepted), otherwise the value is computed as the difference
040 * between the value of the current solution if the neighbour is accepted and
041 * the best ever found solution. <br>
042 * <br>
043 * Custom neighbours can be set using SimulatedAnnealing.Neighbours property that should
044 * contain semicolon separated list of {@link NeighbourSelection}. By default, 
045 * each neighbour selection is selected with the same probability (each has 1 point in
046 * a roulette wheel selection). It can be changed by adding &nbsp;@n at the end
047 * of the name of the class, for example:
048 * <pre><code>
049 * SimulatedAnnealing.Neighbours=org.cpsolver.ifs.algorithms.neighbourhoods.RandomMove;org.cpsolver.ifs.algorithms.neighbourhoods.RandomSwapMove@0.1
050 * </code></pre>
051 * Selector RandomSwapMove is 10&times; less probable to be selected than other selectors.
052 * When SimulatedAnnealing.Random is true, all selectors are selected with the same probability, ignoring these weights.
053 * <br><br>
054 * When SimulatedAnnealing.Update is true, {@link NeighbourSelector#update(Assignment, Neighbour, long)} is called 
055 * after each iteration (on the selector that was used) and roulette wheel selection 
056 * that is using {@link NeighbourSelector#getPoints()} is used to pick a selector in each iteration. 
057 * See {@link NeighbourSelector} for more details. 
058 * <br>
059 * 
060 * @version IFS 1.3 (Iterative Forward Search)<br>
061 *          Copyright (C) 2014 Tomáš Müller<br>
062 *          <a href="mailto:muller@unitime.org">muller@unitime.org</a><br>
063 *          <a href="http://muller.unitime.org">http://muller.unitime.org</a><br>
064 * <br>
065 *          This library is free software; you can redistribute it and/or modify
066 *          it under the terms of the GNU Lesser General Public License as
067 *          published by the Free Software Foundation; either version 3 of the
068 *          License, or (at your option) any later version. <br>
069 * <br>
070 *          This library is distributed in the hope that it will be useful, but
071 *          WITHOUT ANY WARRANTY; without even the implied warranty of
072 *          MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
073 *          Lesser General Public License for more details. <br>
074 * <br>
075 *          You should have received a copy of the GNU Lesser General Public
076 *          License along with this library; if not see
077 *          <a href='http://www.gnu.org/licenses/'>http://www.gnu.org/licenses/</a>.
078 * @param <V> Variable
079 * @param <T> Value
080 */
081public class SimulatedAnnealing<V extends Variable<V, T>, T extends Value<V, T>> extends NeighbourSearch<V,T> {
082    private DecimalFormat iDF5 = new DecimalFormat("0.00000");
083    private DecimalFormat iDF10 = new DecimalFormat("0.0000000000");
084    private double iInitialTemperature = -1;
085    private double iMaximalTemperature = 1.5;
086    private double iMinimalTemperature = 0.01;
087    private double iCoolingRate = 0.95;
088    private double iReheatRate = -1;
089    private long iTemperatureLength = 25000;
090    private double iReheatLengthCoef = 5.0;
091    private double iRestoreBestLengthCoef = -1;
092    private boolean iStochasticHC = false;
093    private boolean iRelativeAcceptance = true;
094    private Double[] iCoolingRateAdjusts = null;
095    private int iTrainingValues = 10000;
096    private double iTrainingProbability = 0.00001;
097    private double iTimeBetweenCooldowns = 10.0;
098
099    /**
100     * Constructor. Following problem properties are considered:
101     * <ul>
102     * <li>SimulatedAnnealing.InitialTemperature ... initial temperature (default 1.5)
103     * <li>SimulatedAnnealing.TemperatureLength ... temperature length (number of iterations between temperature decrements, default 25000)
104     * <li>SimulatedAnnealing.CoolingRate ... temperature cooling rate (default 0.95)
105     * <li>SimulatedAnnealing.ReheatLengthCoef ... temperature re-heat length coefficient (multiple of temperature length, default 5)
106     * <li>SimulatedAnnealing.ReheatRate ... temperature re-heating rate (default (1/coolingRate)^(reheatLengthCoef*1.7))
107     * <li>SimulatedAnnealing.RestoreBestLengthCoef ... restore best length coefficient (multiple of re-heat length, default reheatLengthCoef^2)
108     * <li>SimulatedAnnealing.StochasticHC ... true for stochastic search acceptance criterion, false for simulated annealing acceptance (default false)
109     * <li>SimulatedAnnealing.RelativeAcceptance ... true for relative acceptance (different between the new and the current solutions, if the neighbour is accepted), false for absolute acceptance (difference between the new and the best solutions, if the neighbour is accepted)
110     * <li>SimulatedAnnealing.Neighbours ... semicolon separated list of classes implementing {@link NeighbourSelection}
111     * <li>SimulatedAnnealing.AdditionalNeighbours ... semicolon separated list of classes implementing {@link NeighbourSelection}
112     * <li>SimulatedAnnealing.Random ... when true, a neighbour selector is selected randomly
113     * <li>SimulatedAnnealing.Update ... when true, a neighbour selector is selected using {@link NeighbourSelector#getPoints()} weights (roulette wheel selection)
114     * </ul>
115     * 
116     * @param properties
117     *            problem properties
118     */
119    public SimulatedAnnealing(DataProperties properties) {
120        super(properties);
121        iInitialTemperature = properties.getPropertyDouble(getParameterBaseName() + ".InitialTemperature", iInitialTemperature);
122        iMaximalTemperature = properties.getPropertyDouble(getParameterBaseName() + ".MaximalTemperature", (iInitialTemperature <= 0.0 ? -1.0 : iMaximalTemperature));
123        iMinimalTemperature = properties.getPropertyDouble(getParameterBaseName() + ".MinimalTemperature", (iInitialTemperature <= 0.0 ? -1.0 : 0.0));
124        iReheatRate = properties.getPropertyDouble(getParameterBaseName() + ".ReheatRate", iReheatRate);
125        iCoolingRate = properties.getPropertyDouble(getParameterBaseName() + ".CoolingRate", iCoolingRate);
126        iRelativeAcceptance = properties.getPropertyBoolean(getParameterBaseName() + ".RelativeAcceptance", iRelativeAcceptance);
127        iStochasticHC = properties.getPropertyBoolean(getParameterBaseName() + ".StochasticHC", iStochasticHC);
128        iTemperatureLength = properties.getPropertyLong(getParameterBaseName() + ".TemperatureLength", iTemperatureLength);
129        iReheatLengthCoef = properties.getPropertyDouble(getParameterBaseName() + ".ReheatLengthCoef", iReheatLengthCoef);
130        iRestoreBestLengthCoef = properties.getPropertyDouble(getParameterBaseName() + ".RestoreBestLengthCoef", iRestoreBestLengthCoef);
131        iCoolingRateAdjusts = properties.getPropertyDoubleArry(getParameterBaseName() + ".CoolingRateAdjustments", null);
132        iTrainingValues = properties.getPropertyInt(getParameterBaseName() + ".TrainingValues", iTrainingValues);
133        iTrainingProbability = properties.getPropertyDouble(getParameterBaseName() + ".TrainingProbability", iTrainingProbability);
134        iTimeBetweenCooldowns = properties.getPropertyDouble(getParameterBaseName() + ".TimeBetweenCooldowns", iTimeBetweenCooldowns);
135        if (iReheatRate < 0)
136            iReheatRate = Math.pow(1 / iCoolingRate, iReheatLengthCoef * 1.7);
137        if (iRestoreBestLengthCoef < 0)
138            iRestoreBestLengthCoef = iReheatLengthCoef * iReheatLengthCoef;
139    }
140    
141    @Override
142    public String getParameterBaseName() { return "SimulatedAnnealing"; }
143    
144    @Override
145    public NeighbourSearchContext createAssignmentContext(Assignment<V, T> assignment) {
146        return new SimulatedAnnealingContext();
147    }
148    
149    public class SimulatedAnnealingContext extends NeighbourSearchContext {
150        private double iTemperature = 0.0;
151        private int iMoves = 0;
152        private double iAbsValue = 0;
153        private double iBestValue = 0;
154        private long iLastImprovingIter = -1;
155        private long iLastBestIter = -1;
156        private long iLastReheatIter = 0;
157        private long iLastCoolingIter = 0;
158        private int iAcceptIter[] = new int[] { 0, 0, 0 };
159        private long iReheatLength = 0;
160        private long iRestoreBestLength = 0;
161        private int iTrainingIterations = 0;
162        private double iTrainingTotal = 0.0;
163
164        /** Setup the temperature */
165        @Override
166        protected void activate(Solution<V, T> solution) {
167            super.activate(solution);
168            iTrainingTotal = 0.0; iTrainingIterations = 0;
169            iTemperature = iInitialTemperature;
170            iReheatLength = Math.round(iReheatLengthCoef * iTemperatureLength);
171            iRestoreBestLength = Math.round(iRestoreBestLengthCoef * iTemperatureLength);
172            iLastImprovingIter = -1;
173            iLastBestIter = -1;
174        }
175        
176        protected double getCoolingRate(int idx) {
177            if (idx < 0 || iCoolingRateAdjusts == null || idx >= iCoolingRateAdjusts.length || iCoolingRateAdjusts[idx] == null) return iCoolingRate;
178            return iCoolingRate * iCoolingRateAdjusts[idx];
179        }
180        
181        /**
182         * Set initial temperature based on the training period
183         * @param solution current solution
184         */
185        protected void train(Solution<V, T> solution) {
186            double value = iTrainingTotal / iTrainingIterations;
187            if (iStochasticHC) {
188                iInitialTemperature = value / Math.log(1.0 / 0.01 * iTrainingProbability - 1.0);
189                if (iMaximalTemperature <= 0.0)
190                    iMaximalTemperature = value / Math.log(1.0 / iTrainingProbability - 1.0);
191                if (iMinimalTemperature < 0.0)
192                    iMinimalTemperature = 0.1 / Math.log(1.0 / iTrainingProbability - 1.0);
193            } else {
194                iInitialTemperature = - value / Math.log(0.01 * iTrainingProbability);
195                if (iMaximalTemperature <= 0.0)
196                    iMaximalTemperature = - value / Math.log(iTrainingProbability);
197                if (iMinimalTemperature < 0.0)
198                    iMinimalTemperature = - 0.1 / Math.log(iTrainingProbability);
199            }
200            iTemperature = iInitialTemperature;
201            info("Iter=" + iIter / 1000 + (iLastImprovingIter < 0 ? "" : "k, NonImpIter=" + iDF2.format((iIter - iLastImprovingIter) / 1000.0))
202                    + "k, Speed=" + iDF2.format(1000.0 * iIter / (JProf.currentTimeMillis() - iT0)) + " it/s, " +
203                    "Value=" + iDF2.format(solution.getModel().getTotalValue(solution.getAssignment())) +
204                    ", Best=" + iDF2.format(solution.getBestValue()) +
205                    " (" + iDF2.format(100.0 * solution.getModel().getTotalValue(solution.getAssignment()) / solution.getBestValue()) + " %)");
206            info("Temperature set to " + iDF5.format(iTemperature) + " " + "(" + 
207                    "p(+0.1)=" + iDF2.format(100.0 * prob(0.1)) + "%, " +
208                    "p(+1)=" + iDF2.format(100.0 * prob(1)) + "%, " +
209                    "p(+10)=" + iDF5.format(100.0 * prob(10)) + "%, " +
210                    "p(+" + iDF2.format(value) + ")=" + iDF5.format(100.0 * prob(value)) + "%)");
211            info("Maximal temperature set to " + iDF5.format(iMaximalTemperature) + " " + "(" + 
212                    "p(+0.1)=" + iDF2.format(100.0 * prob(0.1, iMaximalTemperature)) + "%, " +
213                    "p(+1)=" + iDF2.format(100.0 * prob(1, iMaximalTemperature)) + "%, " +
214                    "p(+10)=" + iDF5.format(100.0 * prob(10, iMaximalTemperature)) + "%, " +
215                    "p(+" + iDF2.format(value) + ")=" + iDF5.format(100.0 * prob(value, iMaximalTemperature)) + "%)");
216            info("Minimal temperature set to " + iDF5.format(iMinimalTemperature) + " " + "(" + 
217                    "p(+0.001)=" + iDF2.format(100.0 * prob(0.001, iMinimalTemperature)) + "%, " +
218                    "p(+0.01)=" + iDF2.format(100.0 * prob(0.01, iMinimalTemperature)) + "%, " +
219                    "p(+0.1)=" + iDF5.format(100.0 * prob(0.1, iMinimalTemperature)) + "%)");
220            logNeibourStatus();
221            double speed = ((double)iIter) / (JProf.currentTimeMillis() - iT0); // iterations per ms
222            iTemperatureLength = Math.round(speed * iTimeBetweenCooldowns * 1000.0);
223            iReheatLength = Math.round(iReheatLengthCoef * iTemperatureLength);
224            iRestoreBestLength = Math.round(iRestoreBestLengthCoef * iTemperatureLength);
225            info("Training speed was " + iDF2.format(1000.0 * speed) + " it/s, temperature length adjusted to " + iTemperatureLength + ".");
226            iIter = 0; iT0 = JProf.currentTimeMillis();
227            iLastImprovingIter = -1; iLastBestIter = 0; iBestValue = solution.getBestValue();
228            iAcceptIter = new int[] { 0, 0, 0 };
229            iMoves = 0;
230            iAbsValue = 0;
231        }
232
233        /**
234         * Cool temperature
235         * @param solution current solution
236         */
237        protected void cool(Solution<V, T> solution) {
238            info("Iter=" + iIter / 1000 + (iLastImprovingIter < 0 ? "" : "k, NonImpIter=" + iDF2.format((iIter - iLastImprovingIter) / 1000.0))
239                    + "k, Speed=" + iDF2.format(1000.0 * iIter / (JProf.currentTimeMillis() - iT0)) + " it/s, " +
240                    "Value=" + iDF2.format(solution.getModel().getTotalValue(solution.getAssignment())) +
241                    ", Best=" + iDF2.format(solution.getBestValue()) +
242                    " (" + iDF2.format(100.0 * solution.getModel().getTotalValue(solution.getAssignment()) / solution.getBestValue()) + " %), " +
243                    "Step=" + iDF2.format(1.0 * (iIter - Math.max(iLastReheatIter, iLastImprovingIter)) / iTemperatureLength));
244            if (iLastImprovingIter < 0 || iIter > iLastImprovingIter + iTemperatureLength) {
245                // Do not update the temperature when solution was improved recently 
246                iTemperature *= getCoolingRate(solution.getAssignment().getIndex() - 1);
247                info("Temperature decreased to " + iDF5.format(iTemperature) + " " + "(#moves=" + iMoves + ", rms(value)="
248                        + iDF2.format(Math.sqrt(iAbsValue / iMoves)) + ", " + "accept=-"
249                        + iDF2.format(100.0 * iAcceptIter[0] / iMoves) + "/"
250                        + iDF2.format(100.0 * iAcceptIter[1] / iMoves) + "/+"
251                        + iDF2.format(100.0 * iAcceptIter[2] / iMoves) + "%, "
252                        + (prob(-1) < 1.0 ? "p(-1)=" + iDF2.format(100.0 * prob(-1)) + "%, " : "") + 
253                        "p(+0.1)=" + iDF2.format(100.0 * prob(0.1)) + "%, " +
254                        "p(+1)=" + iDF5.format(100.0 * prob(1)) + "%, " +
255                        "p(+10)=" + iDF10.format(100.0 * prob(10)) + "%)");
256                iAbsValue = 0;
257                iAcceptIter = new int[] { 0, 0, 0 };
258                iMoves = 0;
259            }
260            logNeibourStatus();
261            iLastCoolingIter = iIter;
262        }
263
264        /**
265         * Reheat temperature
266         * @param solution current solution
267         */
268        protected void reheat(Solution<V, T> solution) {
269            iTemperature *= iReheatRate;
270            info("Iter=" + iIter / 1000 + (iLastImprovingIter < 0 ? "" : "k, NonImpIter=" + iDF2.format((iIter - iLastImprovingIter) / 1000.0))
271                    + "k, Speed=" + iDF2.format(1000.0 * iIter / (JProf.currentTimeMillis() - iT0)) + " it/s, " +
272                    "Value=" + iDF2.format(solution.getModel().getTotalValue(solution.getAssignment())) +
273                    ", Best=" + iDF2.format(solution.getBestValue()) +
274                    " (" + iDF2.format(100.0 * solution.getModel().getTotalValue(solution.getAssignment()) / solution.getBestValue()) + " %)");
275            info("Temperature increased to " + iDF5.format(iTemperature) + " "
276                    + (prob(-1) < 1.0 ? "p(-1)=" + iDF2.format(100.0 * prob(-1)) + "%, " : "") + "p(+0.1)="
277                    + iDF2.format(100.0 * prob(0.1)) + "%, " + "p(+1)=" + iDF5.format(10.0 * prob(1)) + "%, " + "p(+10)="
278                    + iDF10.format(100.0 * prob(10)) + "%)");
279            logNeibourStatus();
280            iLastReheatIter = iIter;
281            if (iTemperature > iMaximalTemperature) {
282                // Max temperature reached -> restore best solution and stop re-heating
283                restoreBest(solution);
284                iLastImprovingIter = -1;
285            }
286            iBestValue = solution.getBestValue();
287            setProgressPhase("Simulated Annealing [" + iDF5.format(iTemperature) + "]...");
288        }
289
290        /**
291         * restore best ever found solution
292         * @param solution current solution
293         */
294        protected void restoreBest(Solution<V, T> solution) {
295            info("Best solution restored.");
296            solution.restoreBest();
297            iLastBestIter = iIter;
298        }
299
300        /**
301         * Neighbour acceptance probability
302         * 
303         * @param value
304         *            absolute or relative value of the proposed change (neighbour)
305         * @return probability of acceptance of a change (neighbour)
306         */
307        protected double prob(double value) {
308            return prob(value, iTemperature);
309        }
310        
311        /**
312         * Neighbour acceptance probability
313         * 
314         * @param value
315         *            absolute or relative value of the proposed change (neighbour)
316         * @param temp
317         *      current temperature
318         * @return probability of acceptance of a change (neighbour)
319         */
320        protected double prob(double value, double temp) {
321            if (iStochasticHC)
322                return 1.0 / (1.0 + Math.exp(value / temp));
323            else
324                return (value <= 0.0 ? 1.0 : Math.exp(-value / temp));
325        }
326
327        /**
328         * True if the given neighbour is to be be accepted
329         * 
330         * @param assignment
331         *            current assignment
332         * @param neighbour
333         *            proposed move
334         * @return true if generated random number is below the generated probability
335         */
336        @Override
337        protected boolean accept(Assignment<V, T> assignment, Model<V, T> model, Neighbour<V, T> neighbour, double value, boolean lazy) {
338            iMoves ++;
339            iAbsValue += value * value;
340            double v = (iRelativeAcceptance ? value : (lazy ? model.getTotalValue(assignment) : value + model.getTotalValue(assignment)) - iBestValue);
341            if (iTrainingIterations < iTrainingValues) {
342                if (v <= 0.0) {
343                    iAcceptIter[value < 0.0 ? 0 : value > 0.0 ? 2 : 1]++;
344                    return true;
345                } else {
346                    iTrainingIterations ++; iTrainingTotal += v;
347                }
348                return false;
349            }
350            double prob = prob(v);
351            if (v > 0) {
352                iTrainingIterations ++; iTrainingTotal += v;
353            }
354            if (prob >= 1.0 || ToolBox.random() < prob) {
355                iAcceptIter[value < 0.0 ? 0 : value > 0.0 ? 2 : 1]++;
356                return true;
357            }
358            return false;
359        }
360
361        /**
362         * Increment iteration counter, cool/reheat/restoreBest if necessary
363         */
364        @Override
365        protected void incIteration(Solution<V, T> solution) {
366            super.incIteration(solution);
367            iIter++;
368            if (isMaster(solution)) {
369                if (iInitialTemperature <= 0.0) {
370                    if (iTrainingIterations < iTrainingValues) {
371                        setProgress(Math.round(100.0 * iTrainingIterations / iTrainingValues));
372                        return;
373                    } else {
374                        train(solution);
375                    }
376                }
377                if (iLastBestIter >= 0 && iIter > iLastBestIter + iRestoreBestLength)
378                    restoreBest(solution);
379                if (iTemperature < iMinimalTemperature) { // Minimal temperature reached, start going up again
380                    // If the current solution is more than 0.1% worse, restore best
381                    if ((solution.getModel().getTotalValue(solution.getAssignment()) - solution.getBestValue()) / Math.abs(solution.getBestValue()) >= 0.001)
382                        restoreBest(solution);
383                    // Set last improving iteration so that the re-heating would continue
384                    iLastImprovingIter = iIter;
385                    // Do the first reheat
386                    reheat(solution);
387                } else if (iLastImprovingIter >= 0 && iIter > Math.max(iLastReheatIter, iLastImprovingIter) + iReheatLength) {
388                    reheat(solution);
389                } else if (iIter > iLastCoolingIter + iTemperatureLength) {
390                    cool(solution);
391                }
392                setProgress(Math.round(100.0 * (iIter - Math.max(iLastReheatIter, iLastImprovingIter)) / iReheatLength));
393            }
394        }
395
396        /**
397         * Memorize the iteration when the last best solution was found.
398         */
399        @Override
400        public void bestSaved(Solution<V, T> solution) {
401            super.bestSaved(solution);
402            if (iLastImprovingIter < 0 || Math.abs(iBestValue - solution.getBestValue()) / Math.max(Math.abs(iBestValue), Math.abs(solution.getBestValue())) >= 0.0001) {
403                iLastImprovingIter = iIter;
404                iLastBestIter = iIter;
405                iBestValue = solution.getBestValue();
406            }
407        }
408    }
409
410}