001package net.sf.cpsolver.ifs.algorithms;
002
003import java.text.DecimalFormat;
004
005import net.sf.cpsolver.ifs.heuristics.NeighbourSelection;
006import net.sf.cpsolver.ifs.model.Model;
007import net.sf.cpsolver.ifs.model.Neighbour;
008import net.sf.cpsolver.ifs.model.Value;
009import net.sf.cpsolver.ifs.model.Variable;
010import net.sf.cpsolver.ifs.solution.Solution;
011import net.sf.cpsolver.ifs.util.DataProperties;
012import net.sf.cpsolver.ifs.util.JProf;
013import net.sf.cpsolver.ifs.util.ToolBox;
014
015/**
016 * Simulated annealing. In each iteration, one of the given neighbourhoods is selected first,
017 * then a neighbour is generated and it is accepted with probability
018 * {@link SimulatedAnnealing#prob(double)}. The search is guided by the
019 * temperature, which starts at <i>SimulatedAnnealing.InitialTemperature</i>.
020 * After each <i>SimulatedAnnealing.TemperatureLength</i> iterations, the
021 * temperature is reduced by <i>SimulatedAnnealing.CoolingRate</i>. If there was
022 * no improvement in the past <i>SimulatedAnnealing.ReheatLengthCoef *
023 * SimulatedAnnealing.TemperatureLength</i> iterations, the temperature is
024 * increased by <i>SimulatedAnnealing.ReheatRate</i>. If there was no
025 * improvement in the past <i>SimulatedAnnealing.RestoreBestLengthCoef *
026 * SimulatedAnnealing.TemperatureLength</i> iterations, the best ever found
027 * solution is restored. <br>
028 * <br>
029 * If <i>SimulatedAnnealing.StochasticHC</i> is true, the acceptance probability
030 * is computed using stochastic hill climbing criterion, i.e.,
031 * <code>1.0 / (1.0 + Math.exp(value/temperature))</code>, otherwise it is
032 * cumputed using simlated annealing criterion, i.e.,
033 * <code>(value<=0.0?1.0:Math.exp(-value/temperature))</code>. If
034 * <i>SimulatedAnnealing.RelativeAcceptance</i> neighbour value
035 * {@link Neighbour#value()} is taken as the value of the selected
036 * neighbour (difference between the new and the current solution, if the
037 * neighbour is accepted), otherwise the value is computed as the difference
038 * between the value of the current solution if the neighbour is accepted and
039 * the best ever found solution. <br>
040 * <br>
041 * Custom neighbours can be set using SimulatedAnnealing.Neighbours property that should
042 * contain semicolon separated list of {@link NeighbourSelection}. By default, 
043 * each neighbour selection is selected with the same probability (each has 1 point in
044 * a roulette wheel selection). It can be changed by adding &nbsp;@n at the end
045 * of the name of the class, for example:<br>
046 * <code>
047 * SimulatedAnnealing.Neighbours=net.sf.cpsolver.ifs.algorithms.neighbourhoods.RandomMove;net.sf.cpsolver.ifs.algorithms.neighbourhoods.RandomSwapMove@0.1
048 * </code>
049 * <br>
050 * Selector RandomSwapMove is 10&times; less probable to be selected than other selectors.
051 * When SimulatedAnnealing.Random is true, all selectors are selected with the same probability, ignoring these weights.
052 * <br><br>
053 * When SimulatedAnnealing.Update is true, {@link NeighbourSelector#update(Neighbour, long)} is called 
054 * after each iteration (on the selector that was used) and roulette wheel selection 
055 * that is using {@link NeighbourSelector#getPoints()} is used to pick a selector in each iteration. 
056 * See {@link NeighbourSelector} for more details. 
057 * <br>
058 * 
059 * @version IFS 1.2 (Iterative Forward Search)<br>
060 *          Copyright (C) 2014 Tomáš Müller<br>
061 *          <a href="mailto:muller@unitime.org">muller@unitime.org</a><br>
062 *          <a href="http://muller.unitime.org">http://muller.unitime.org</a><br>
063 * <br>
064 *          This library is free software; you can redistribute it and/or modify
065 *          it under the terms of the GNU Lesser General Public License as
066 *          published by the Free Software Foundation; either version 3 of the
067 *          License, or (at your option) any later version. <br>
068 * <br>
069 *          This library is distributed in the hope that it will be useful, but
070 *          WITHOUT ANY WARRANTY; without even the implied warranty of
071 *          MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
072 *          Lesser General Public License for more details. <br>
073 * <br>
074 *          You should have received a copy of the GNU Lesser General Public
075 *          License along with this library; if not see
076 *          <a href='http://www.gnu.org/licenses/'>http://www.gnu.org/licenses/</a>.
077 */
078public class SimulatedAnnealing<V extends Variable<V, T>, T extends Value<V, T>> extends NeighbourSearch<V,T> {
079    private DecimalFormat iDF5 = new DecimalFormat("0.00000");
080    private DecimalFormat iDF10 = new DecimalFormat("0.0000000000");
081    private double iInitialTemperature = 1.5;
082    private double iCoolingRate = 0.95;
083    private double iReheatRate = -1;
084    private long iTemperatureLength = 25000;
085    private long iReheatLength = 0;
086    private long iRestoreBestLength = 0;
087    private double iTemperature = 0.0;
088    private double iReheatLengthCoef = 5.0;
089    private double iRestoreBestLengthCoef = -1;
090    private long iLastImprovingIter = 0;
091    private long iLastReheatIter = 0;
092    private long iLastCoolingIter = 0;
093    private int iAcceptIter[] = new int[] { 0, 0, 0 };
094    private boolean iStochasticHC = false;
095    private int iMoves = 0;
096    private double iAbsValue = 0;
097    private double iBestValue = 0;
098
099    private boolean iRelativeAcceptance = true;
100
101    /**
102     * Constructor. Following problem properties are considered:
103     * <ul>
104     * <li>SimulatedAnnealing.InitialTemperature ... initial temperature (default 1.5)
105     * <li>SimulatedAnnealing.TemperatureLength ... temperature length (number of iterations between temperature decrements, default 25000)
106     * <li>SimulatedAnnealing.CoolingRate ... temperature cooling rate (default 0.95)
107     * <li>SimulatedAnnealing.ReheatLengthCoef ... temperature re-heat length coefficient (multiple of temperature length, default 5)
108     * <li>SimulatedAnnealing.ReheatRate ... temperature re-heating rate (default (1/coolingRate)^(reheatLengthCoef*1.7))
109     * <li>SimulatedAnnealing.RestoreBestLengthCoef ... restore best length coefficient (multiple of re-heat length, default reheatLengthCoef^2)
110     * <li>SimulatedAnnealing.StochasticHC ... true for stochastic search acceptance criterion, false for simulated annealing acceptance (default false)
111     * <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)
112     * <li>SimulatedAnnealing.Neighbours ... semicolon separated list of classes implementing {@link NeighbourSelection}
113     * <li>SimulatedAnnealing.AdditionalNeighbours ... semicolon separated list of classes implementing {@link NeighbourSelection}
114     * <li>SimulatedAnnealing.Random ... when true, a neighbour selector is selected randomly
115     * <li>SimulatedAnnealing.Update ... when true, a neighbour selector is selected using {@link NeighbourSelector#getPoints()} weights (roulette wheel selection)
116     * </ul>
117     * 
118     * @param properties
119     *            problem properties
120     */
121    public SimulatedAnnealing(DataProperties properties) {
122        super(properties);
123        iInitialTemperature = properties.getPropertyDouble(getParameterBaseName() + ".InitialTemperature", iInitialTemperature);
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        if (iReheatRate < 0)
132            iReheatRate = Math.pow(1 / iCoolingRate, iReheatLengthCoef * 1.7);
133        if (iRestoreBestLengthCoef < 0)
134            iRestoreBestLengthCoef = iReheatLengthCoef * iReheatLengthCoef;
135    }
136    
137    /** Setup the temperature */
138    @Override
139    protected void activate(Solution<V, T> solution) {
140        super.activate(solution);
141        iTemperature = iInitialTemperature;
142        iReheatLength = Math.round(iReheatLengthCoef * iTemperatureLength);
143        iRestoreBestLength = Math.round(iRestoreBestLengthCoef * iTemperatureLength);
144        iLastImprovingIter = iIter;
145    }
146
147    /**
148     * Cool temperature
149     */
150    protected void cool(Solution<V, T> solution) {
151        iTemperature *= iCoolingRate;
152        iLog.info("Iter=" + iIter / 1000 + "k, NonImpIter=" + iDF2.format((iIter - iLastImprovingIter) / 1000.0)
153                + "k, Speed=" + iDF2.format(1000.0 * iIter / (JProf.currentTimeMillis() - iT0)) + " it/s");
154        iLog.info("Temperature decreased to " + iDF5.format(iTemperature) + " " + "(#moves=" + iMoves + ", rms(value)="
155                + iDF2.format(Math.sqrt(iAbsValue / iMoves)) + ", " + "accept=-"
156                + iDF2.format(100.0 * iAcceptIter[0] / iMoves) + "/"
157                + iDF2.format(100.0 * iAcceptIter[1] / iMoves) + "/+"
158                + iDF2.format(100.0 * iAcceptIter[2] / iMoves) + "%, "
159                + (prob(-1) < 1.0 ? "p(-1)=" + iDF2.format(100.0 * prob(-1)) + "%, " : "") + 
160                "p(+0.1)=" + iDF2.format(100.0 * prob(0.1)) + "%, " +
161                "p(+1)=" + iDF2.format(100.0 * prob(1)) + "%, " +
162                "p(+10)=" + iDF5.format(100.0 * prob(10)) + "%)");
163        logNeibourStatus();
164        iLastCoolingIter = iIter;
165        iAcceptIter = new int[] { 0, 0, 0 };
166        iMoves = 0;
167        iAbsValue = 0;
168    }
169
170    /**
171     * Reheat temperature
172     */
173    protected void reheat(Solution<V, T> solution) {
174        iTemperature *= iReheatRate;
175        iLog.info("Iter=" + iIter / 1000 + "k, NonImpIter=" + iDF2.format((iIter - iLastImprovingIter) / 1000.0)
176                + "k, Speed=" + iDF2.format(1000.0 * iIter / (JProf.currentTimeMillis() - iT0)) + " it/s");
177        iLog.info("Temperature increased to " + iDF5.format(iTemperature) + " "
178                + (prob(-1) < 1.0 ? "p(-1)=" + iDF2.format(100.0 * prob(-1)) + "%, " : "") + "p(+1)="
179                + iDF2.format(100.0 * prob(1)) + "%, " + "p(+10)=" + iDF5.format(100.0 * prob(10)) + "%, " + "p(+100)="
180                + iDF10.format(100.0 * prob(100)) + "%)");
181        logNeibourStatus();
182        iLastReheatIter = iIter;
183        iProgress.setPhase("Simulated Annealing [" + iDF2.format(iTemperature) + "]...");
184    }
185
186    /**
187     * restore best ever found solution
188     */
189    protected void restoreBest(Solution<V, T> solution) {
190        iLog.info("Best restored");
191        iLastImprovingIter = iIter;
192    }
193
194    /**
195     * Neighbour acceptance probability
196     * 
197     * @param value
198     *            absolute or relative value of the proposed change (neighbour)
199     * @return probability of acceptance of a change (neighbour)
200     */
201    protected double prob(double value) {
202        if (iStochasticHC)
203            return 1.0 / (1.0 + Math.exp(value / iTemperature));
204        else
205            return (value <= 0.0 ? 1.0 : Math.exp(-value / iTemperature));
206    }
207
208    /**
209     * True if the given neighbour is to be be accepted
210     * 
211     * @param solution
212     *            current solution
213     * @param neighbour
214     *            proposed move
215     * @return true if generated random number is below
216     *         {@link SimulatedAnnealing#prob(double)}
217     */
218    @Override
219    protected boolean accept(Model<V, T> model, Neighbour<V, T> neighbour, double value, boolean lazy) {
220        iMoves ++;
221        iAbsValue += value * value;
222        double prob = prob(iRelativeAcceptance ? value : (lazy ? model.getTotalValue() : value + model.getTotalValue()) - iBestValue);
223        if (prob >= 1.0 || ToolBox.random() < prob) {
224            iAcceptIter[value < 0.0 ? 0 : value > 0.0 ? 2 : 1]++;
225            return true;
226        }
227        return false;
228    }
229
230    /**
231     * Increment iteration counter, cool/reheat/restoreBest if necessary
232     */
233    @Override
234    protected void incIteration(Solution<V, T> solution) {
235        super.incIteration(solution);
236        iIter++;
237        if (iIter > iLastImprovingIter + iRestoreBestLength)
238            restoreBest(solution);
239        if (iIter > Math.max(iLastReheatIter, iLastImprovingIter) + iReheatLength)
240            reheat(solution);
241        if (iIter > iLastCoolingIter + iTemperatureLength)
242            cool(solution);
243        iProgress.setProgress(Math.round(100.0 * (iIter - Math.max(iLastReheatIter, iLastImprovingIter)) / iReheatLength));
244    }
245
246    /**
247     * Memorize the iteration when the last best solution was found.
248     */
249    @Override
250    public void bestSaved(Solution<V, T> solution) {
251        if (Math.abs(iBestValue - solution.getBestValue()) >= 1.0) {
252            iLastImprovingIter = iIter;
253            iBestValue = solution.getBestValue();
254        }
255        iLastImprovingIter = iIter;
256    }
257
258    @Override
259    public String getParameterBaseName() { return "SimulatedAnnealing"; }
260}