    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.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            if (iTemperature > iInitialTemperature)
    248                iLastImprovingIter = -1;
    249            iBestValue = solution.getBestValue();
    250            setProgressPhase("Simulated Annealing [" + iDF5.format(iTemperature) + "]...");
    251        }
    252
    253        /**
    254         * restore best ever found solution
    255         * @param solution current solution
    256         */
    257        protected void restoreBest(Solution<V, T> solution) {
    258            info("Best solution restored.");
    259            solution.restoreBest();
    260            iLastImprovingIter = -1;
    261        }
    262
    263        /**
    264         * Neighbour acceptance probability
    265         *
    266         * @param value
    267         *            absolute or relative value of the proposed change (neighbour)
    268         * @return probability of acceptance of a change (neighbour)
    269         */
    270        protected double prob(double value) {
    271            if (iStochasticHC)
    272                return 1.0 / (1.0 + Math.exp(value / iTemperature));
    273            else
    274                return (value <= 0.0 ? 1.0 : Math.exp(-value / iTemperature));
    275        }
    276
    277        /**
    278         * True if the given neighbour is to be be accepted
    279         *
    280         * @param assignment
    281         *            current assignment
    282         * @param neighbour
    283         *            proposed move
    284         * @return true if generated random number is below the generated probability
    285         */
    286        @Override
    287        protected boolean accept(Assignment<V, T> assignment, Model<V, T> model, Neighbour<V, T> neighbour, double value, boolean lazy) {
    288            iMoves ++;
    289            iAbsValue += value * value;
    290            double v = (iRelativeAcceptance ? value : (lazy ? model.getTotalValue(assignment) : value + model.getTotalValue(assignment)) - iBestValue);
    291            if (iTrainingIterations < iTrainingValues) {
    292                if (v <= 0.0) {
    293                    iAcceptIter[value < 0.0 ? 0 : value > 0.0 ? 2 : 1]++;
    294                    return true;
    295                } else {
    296                    iTrainingIterations ++; iTrainingTotal += v;
    297                }
    298                return false;
    299            }
    300            double prob = prob(v);
    301            if (v > 0) {
    302                iTrainingIterations ++; iTrainingTotal += v;
    303            }
    304            if (prob >= 1.0 || ToolBox.random() < prob) {
    305                iAcceptIter[value < 0.0 ? 0 : value > 0.0 ? 2 : 1]++;
    306                return true;
    307            }
    308            return false;
    309        }
    310
    311        /**
    312         * Increment iteration counter, cool/reheat/restoreBest if necessary
    313         */
    314        @Override
    315        protected void incIteration(Solution<V, T> solution) {
    316            super.incIteration(solution);
    317            iIter++;
    318            if (isMaster(solution)) {
    319                if (iInitialTemperature <= 0.0) {
    320                    if (iTrainingIterations < iTrainingValues) {
    321                        setProgress(Math.round(100.0 * iTrainingIterations / iTrainingValues));
    322                        return;
    323                    } else {
    324                        train(solution);
    325                    }
    326                }
    327                if (iLastImprovingIter >= 0 && iIter > iLastImprovingIter + iRestoreBestLength)
    328                    restoreBest(solution);
    329                if (iLastImprovingIter >= 0 && iIter > Math.max(iLastReheatIter, iLastImprovingIter) + iReheatLength)
    330                    reheat(solution);
    331                if (iIter > iLastCoolingIter + iTemperatureLength)
    332                    cool(solution);
    333                setProgress(Math.round(100.0 * (iIter - Math.max(iLastReheatIter, iLastImprovingIter)) / iReheatLength));
    334            }
    335        }
    336
    337        /**
    338         * Memorize the iteration when the last best solution was found.
    339         */
    340        @Override
    341        public void bestSaved(Solution<V, T> solution) {
    342            super.bestSaved(solution);
    343            if (Math.abs(iBestValue - solution.getBestValue()) >= 1.0) {
    344                iLastImprovingIter = iIter;
    345                iBestValue = solution.getBestValue();
    346            }
    347        }
    348    }
    349
    350}