    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}