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 Tomas Muller<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.5;
085    private double iCoolingRate = 0.95;
086    private double iReheatRate = -1;
087    private long iTemperatureLength = 25000;
088    private double iReheatLengthCoef = 5.0;
089    private double iRestoreBestLengthCoef = -1;
090    private boolean iStochasticHC = false;
091    private boolean iRelativeAcceptance = true;
092    private Double[] iCoolingRateAdjusts = null;
093    private int iTrainingValues = 10000;
094    private double iTrainingProbability = 0.01;
095
096    /**
097     * Constructor. Following problem properties are considered:
098     * <ul>
099     * <li>SimulatedAnnealing.InitialTemperature ... initial temperature (default 1.5)
100     * <li>SimulatedAnnealing.TemperatureLength ... temperature length (number of iterations between temperature decrements, default 25000)
101     * <li>SimulatedAnnealing.CoolingRate ... temperature cooling rate (default 0.95)
102     * <li>SimulatedAnnealing.ReheatLengthCoef ... temperature re-heat length coefficient (multiple of temperature length, default 5)
103     * <li>SimulatedAnnealing.ReheatRate ... temperature re-heating rate (default (1/coolingRate)^(reheatLengthCoef*1.7))
104     * <li>SimulatedAnnealing.RestoreBestLengthCoef ... restore best length coefficient (multiple of re-heat length, default reheatLengthCoef^2)
105     * <li>SimulatedAnnealing.StochasticHC ... true for stochastic search acceptance criterion, false for simulated annealing acceptance (default false)
106     * <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)
107     * <li>SimulatedAnnealing.Neighbours ... semicolon separated list of classes implementing {@link NeighbourSelection}
108     * <li>SimulatedAnnealing.AdditionalNeighbours ... semicolon separated list of classes implementing {@link NeighbourSelection}
109     * <li>SimulatedAnnealing.Random ... when true, a neighbour selector is selected randomly
110     * <li>SimulatedAnnealing.Update ... when true, a neighbour selector is selected using {@link NeighbourSelector#getPoints()} weights (roulette wheel selection)
111     * </ul>
112     * 
113     * @param properties
114     *            problem properties
115     */
116    public SimulatedAnnealing(DataProperties properties) {
117        super(properties);
118        iInitialTemperature = properties.getPropertyDouble(getParameterBaseName() + ".InitialTemperature", iInitialTemperature);
119        iReheatRate = properties.getPropertyDouble(getParameterBaseName() + ".ReheatRate", iReheatRate);
120        iCoolingRate = properties.getPropertyDouble(getParameterBaseName() + ".CoolingRate", iCoolingRate);
121        iRelativeAcceptance = properties.getPropertyBoolean(getParameterBaseName() + ".RelativeAcceptance", iRelativeAcceptance);
122        iStochasticHC = properties.getPropertyBoolean(getParameterBaseName() + ".StochasticHC", iStochasticHC);
123        iTemperatureLength = properties.getPropertyLong(getParameterBaseName() + ".TemperatureLength", iTemperatureLength);
124        iReheatLengthCoef = properties.getPropertyDouble(getParameterBaseName() + ".ReheatLengthCoef", iReheatLengthCoef);
125        iRestoreBestLengthCoef = properties.getPropertyDouble(getParameterBaseName() + ".RestoreBestLengthCoef", iRestoreBestLengthCoef);
126        iCoolingRateAdjusts = properties.getPropertyDoubleArry(getParameterBaseName() + ".CoolingRateAdjustments", null);
127        iTrainingValues = properties.getPropertyInt(getParameterBaseName() + ".TrainingValues", iTrainingValues);
128        iTrainingProbability = properties.getPropertyDouble(getParameterBaseName() + ".TrainingProbability", iTrainingProbability);
129        if (iReheatRate < 0)
130            iReheatRate = Math.pow(1 / iCoolingRate, iReheatLengthCoef * 1.7);
131        if (iRestoreBestLengthCoef < 0)
132            iRestoreBestLengthCoef = iReheatLengthCoef * iReheatLengthCoef;
133    }
134    
135    @Override
136    public String getParameterBaseName() { return "SimulatedAnnealing"; }
137    
138    @Override
139    public NeighbourSearchContext createAssignmentContext(Assignment<V, T> assignment) {
140        return new SimulatedAnnealingContext();
141    }
142    
143    public class SimulatedAnnealingContext extends NeighbourSearchContext {
144        private double iTemperature = 0.0;
145        private int iMoves = 0;
146        private double iAbsValue = 0;
147        private double iBestValue = 0;
148        private long iLastImprovingIter = -1;
149        private long iLastReheatIter = 0;
150        private long iLastCoolingIter = 0;
151        private int iAcceptIter[] = new int[] { 0, 0, 0 };
152        private long iReheatLength = 0;
153        private long iRestoreBestLength = 0;
154        private int iTrainingIterations = 0;
155        private double iTrainingTotal = 0.0;
156
157        /** Setup the temperature */
158        @Override
159        protected void activate(Solution<V, T> solution) {
160            super.activate(solution);
161            iTrainingTotal = 0.0; iTrainingIterations = 0;
162            iTemperature = iInitialTemperature;
163            iReheatLength = Math.round(iReheatLengthCoef * iTemperatureLength);
164            iRestoreBestLength = Math.round(iRestoreBestLengthCoef * iTemperatureLength);
165            iLastImprovingIter = -1;
166        }
167        
168        protected double getCoolingRate(int idx) {
169            if (idx < 0 || iCoolingRateAdjusts == null || idx >= iCoolingRateAdjusts.length || iCoolingRateAdjusts[idx] == null) return iCoolingRate;
170            return iCoolingRate * iCoolingRateAdjusts[idx];
171        }
172        
173        /**
174         * Set initial temperature based on the training period
175         * @param solution current solution
176         */
177        protected void train(Solution<V, T> solution) {
178            double value = iTrainingTotal / iTrainingIterations;
179            if (iStochasticHC) {
180                iInitialTemperature = - value / Math.log(1.0 / iTrainingProbability - 1.0);
181            } else {
182                iInitialTemperature = - value / Math.log(iTrainingProbability);
183            }
184            iTemperature = iInitialTemperature;
185            info("Iter=" + iIter / 1000 + (iLastImprovingIter < 0 ? "" : "k, NonImpIter=" + iDF2.format((iIter - iLastImprovingIter) / 1000.0))
186                    + "k, Speed=" + iDF2.format(1000.0 * iIter / (JProf.currentTimeMillis() - iT0)) + " it/s, " +
187                    "Value=" + iDF2.format(solution.getModel().getTotalValue(solution.getAssignment())) +
188                    ", Best=" + iDF2.format(solution.getBestValue()) +
189                    " (" + iDF2.format(100.0 * solution.getModel().getTotalValue(solution.getAssignment()) / solution.getBestValue()) + " %)");
190            info("Temperature set to " + iDF5.format(iTemperature) + " " + "(" + 
191                    "p(+0.1)=" + iDF2.format(100.0 * prob(0.1)) + "%, " +
192                    "p(+1)=" + iDF2.format(100.0 * prob(1)) + "%, " +
193                    "p(+10)=" + iDF5.format(100.0 * prob(10)) + "%, " +
194                    "p(+" + iDF2.format(value) + ")=" + iDF5.format(100.0 * prob(value)) + "%)");
195            logNeibourStatus();
196            iIter = 0; iLastImprovingIter = -1; iBestValue = solution.getBestValue();
197            iAcceptIter = new int[] { 0, 0, 0 };
198            iMoves = 0;
199            iAbsValue = 0;
200        }
201
202        /**
203         * Cool temperature
204         * @param solution current solution
205         */
206        protected void cool(Solution<V, T> solution) {
207            iTemperature *= getCoolingRate(solution.getAssignment().getIndex() - 1);
208            info("Iter=" + iIter / 1000 + (iLastImprovingIter < 0 ? "" : "k, NonImpIter=" + iDF2.format((iIter - iLastImprovingIter) / 1000.0))
209                    + "k, Speed=" + iDF2.format(1000.0 * iIter / (JProf.currentTimeMillis() - iT0)) + " it/s, " +
210                    "Value=" + iDF2.format(solution.getModel().getTotalValue(solution.getAssignment())) +
211                    ", Best=" + iDF2.format(solution.getBestValue()) +
212                    " (" + iDF2.format(100.0 * solution.getModel().getTotalValue(solution.getAssignment()) / solution.getBestValue()) + " %), " +
213                    "Step=" + iDF2.format(1.0 * (iIter - Math.max(iLastReheatIter, iLastImprovingIter)) / iTemperatureLength));
214            info("Temperature decreased to " + iDF5.format(iTemperature) + " " + "(#moves=" + iMoves + ", rms(value)="
215                    + iDF2.format(Math.sqrt(iAbsValue / iMoves)) + ", " + "accept=-"
216                    + iDF2.format(100.0 * iAcceptIter[0] / iMoves) + "/"
217                    + iDF2.format(100.0 * iAcceptIter[1] / iMoves) + "/+"
218                    + iDF2.format(100.0 * iAcceptIter[2] / iMoves) + "%, "
219                    + (prob(-1) < 1.0 ? "p(-1)=" + iDF2.format(100.0 * prob(-1)) + "%, " : "") + 
220                    "p(+0.1)=" + iDF2.format(100.0 * prob(0.1)) + "%, " +
221                    "p(+1)=" + iDF2.format(100.0 * prob(1)) + "%, " +
222                    "p(+10)=" + iDF5.format(100.0 * prob(10)) + "%)");
223            logNeibourStatus();
224            iLastCoolingIter = iIter;
225            iAcceptIter = new int[] { 0, 0, 0 };
226            iMoves = 0;
227            iAbsValue = 0;
228        }
229
230        /**
231         * Reheat temperature
232         * @param solution current solution
233         */
234        protected void reheat(Solution<V, T> solution) {
235            iTemperature *= iReheatRate;
236            info("Iter=" + iIter / 1000 + (iLastImprovingIter < 0 ? "" : "k, NonImpIter=" + iDF2.format((iIter - iLastImprovingIter) / 1000.0))
237                    + "k, Speed=" + iDF2.format(1000.0 * iIter / (JProf.currentTimeMillis() - iT0)) + " it/s, " +
238                    "Value=" + iDF2.format(solution.getModel().getTotalValue(solution.getAssignment())) +
239                    ", Best=" + iDF2.format(solution.getBestValue()) +
240                    " (" + iDF2.format(100.0 * solution.getModel().getTotalValue(solution.getAssignment()) / solution.getBestValue()) + " %)");
241            info("Temperature increased to " + iDF5.format(iTemperature) + " "
242                    + (prob(-1) < 1.0 ? "p(-1)=" + iDF2.format(100.0 * prob(-1)) + "%, " : "") + "p(+1)="
243                    + iDF2.format(100.0 * prob(1)) + "%, " + "p(+10)=" + iDF5.format(100.0 * prob(10)) + "%, " + "p(+100)="
244                    + iDF10.format(100.0 * prob(100)) + "%)");
245            logNeibourStatus();
246            iLastReheatIter = iIter;
247            iLastImprovingIter = -1; iBestValue = solution.getBestValue();
248            setProgressPhase("Simulated Annealing [" + iDF2.format(iTemperature) + "]...");
249        }
250
251        /**
252         * restore best ever found solution
253         * @param solution current solution
254         */
255        protected void restoreBest(Solution<V, T> solution) {
256            solution.restoreBest();
257            iLastImprovingIter = -1;
258        }
259
260        /**
261         * Neighbour acceptance probability
262         * 
263         * @param value
264         *            absolute or relative value of the proposed change (neighbour)
265         * @return probability of acceptance of a change (neighbour)
266         */
267        protected double prob(double value) {
268            if (iStochasticHC)
269                return 1.0 / (1.0 + Math.exp(value / iTemperature));
270            else
271                return (value <= 0.0 ? 1.0 : Math.exp(-value / iTemperature));
272        }
273
274        /**
275         * True if the given neighbour is to be be accepted
276         * 
277         * @param assignment
278         *            current assignment
279         * @param neighbour
280         *            proposed move
281         * @return true if generated random number is below the generated probability
282         */
283        @Override
284        protected boolean accept(Assignment<V, T> assignment, Model<V, T> model, Neighbour<V, T> neighbour, double value, boolean lazy) {
285            iMoves ++;
286            iAbsValue += value * value;
287            double v = (iRelativeAcceptance ? value : (lazy ? model.getTotalValue(assignment) : value + model.getTotalValue(assignment)) - iBestValue);
288            if (iTrainingIterations < iTrainingValues) {
289                if (v <= 0.0) {
290                    iAcceptIter[value < 0.0 ? 0 : value > 0.0 ? 2 : 1]++;
291                    return true;
292                } else {
293                    iTrainingIterations ++; iTrainingTotal += v;
294                }
295                return false;
296            }
297            double prob = prob(v);
298            if (v > 0) {
299                iTrainingIterations ++; iTrainingTotal += v;
300            }
301            if (prob >= 1.0 || ToolBox.random() < prob) {
302                iAcceptIter[value < 0.0 ? 0 : value > 0.0 ? 2 : 1]++;
303                return true;
304            }
305            return false;
306        }
307
308        /**
309         * Increment iteration counter, cool/reheat/restoreBest if necessary
310         */
311        @Override
312        protected void incIteration(Solution<V, T> solution) {
313            super.incIteration(solution);
314            iIter++;
315            if (isMaster(solution)) {
316                if (iInitialTemperature <= 0.0) {
317                    if (iTrainingIterations < iTrainingValues) {
318                        setProgress(Math.round(100.0 * iTrainingIterations / iTrainingValues));
319                        return;
320                    } else {
321                        train(solution);
322                    }
323                }
324                if (iLastImprovingIter >= 0 && iIter > iLastImprovingIter + iRestoreBestLength)
325                    restoreBest(solution);
326                if (iLastImprovingIter >= 0 && iIter > Math.max(iLastReheatIter, iLastImprovingIter) + iReheatLength)
327                    reheat(solution);
328                if (iIter > iLastCoolingIter + iTemperatureLength)
329                    cool(solution);
330                setProgress(Math.round(100.0 * (iIter - Math.max(iLastReheatIter, iLastImprovingIter)) / iReheatLength));
331            }
332        }
333
334        /**
335         * Memorize the iteration when the last best solution was found.
336         */
337        @Override
338        public void bestSaved(Solution<V, T> solution) {
339            super.bestSaved(solution);
340            if (Math.abs(iBestValue - solution.getBestValue()) >= 1.0) {
341                iLastImprovingIter = iIter;
342                iBestValue = solution.getBestValue();
343            }
344        }
345    }
346
347}