001    package net.sf.cpsolver.exam.heuristics;
002    
003    import java.text.DecimalFormat;
004    
005    import net.sf.cpsolver.exam.neighbours.ExamRandomMove;
006    import net.sf.cpsolver.exam.neighbours.ExamRoomMove;
007    import net.sf.cpsolver.exam.neighbours.ExamSimpleNeighbour;
008    import net.sf.cpsolver.exam.neighbours.ExamTimeMove;
009    import net.sf.cpsolver.ifs.heuristics.NeighbourSelection;
010    import net.sf.cpsolver.ifs.model.Neighbour;
011    import net.sf.cpsolver.ifs.solution.Solution;
012    import net.sf.cpsolver.ifs.solution.SolutionListener;
013    import net.sf.cpsolver.ifs.solver.Solver;
014    import net.sf.cpsolver.ifs.util.DataProperties;
015    import net.sf.cpsolver.ifs.util.Progress;
016    import net.sf.cpsolver.ifs.util.ToolBox;
017    
018    
019    import org.apache.log4j.Logger;
020    
021    /**
022     * Simulated annealing. In each iteration, one of the following three neighbourhoods is selected first
023     * <ul>
024     * <li>random move ({@link ExamRandomMove})
025     * <li>period swap ({@link ExamTimeMove})
026     * <li>room swap ({@link ExamRoomMove})
027     * </ul>,
028     * then a neighbour is generated and it is accepted with probability {@link ExamSimulatedAnnealing#prob(double)}.
029     * The search is guided by the temperature, which starts at <i>SimulatedAnnealing.InitialTemperature</i>. 
030     * After each <i>SimulatedAnnealing.TemperatureLength</i> iterations, the temperature is reduced 
031     * by <i>SimulatedAnnealing.CoolingRate</i>. If there was no improvement in the past 
032     * <i>SimulatedAnnealing.ReheatLengthCoef * SimulatedAnnealing.TemperatureLength</i> iterations, 
033     * the temperature is increased by <i>SimulatedAnnealing.ReheatRate</i>.
034     * If there was no improvement in the past 
035     * <i>SimulatedAnnealing.RestoreBestLengthCoef * SimulatedAnnealing.TemperatureLength</i> iterations,
036     * the best ever found solution is restored.
037     * <br><br>
038     * If <i>SimulatedAnnealing.StochasticHC</i> is true, the acceptance probability is computed using
039     * stochastic hill climbing criterion, i.e., <code>1.0 / (1.0 + Math.exp(value/temperature))</code>,
040     * otherwise it is cumputed using simlated annealing criterion, i.e.,
041     * <code>(value<=0.0?1.0:Math.exp(-value/temperature))</code>.
042     * If <i>SimulatedAnnealing.RelativeAcceptance</i> neighbour value {@link ExamSimpleNeighbour#value()} is taken
043     * as the value of the selected neighbour (difference between the new and the current solution, if the
044     * neighbour is accepted), otherwise the value is computed as the difference
045     * between the value of the current solution if the neighbour is accepted and the best ever found solution.
046     * <br><br>
047     * 
048     * @version
049     * ExamTT 1.1 (Examination Timetabling)<br>
050     * Copyright (C) 2008 Tomáš Müller<br>
051     * <a href="mailto:muller@unitime.org">muller@unitime.org</a><br>
052     * Lazenska 391, 76314 Zlin, Czech Republic<br>
053     * <br>
054     * This library is free software; you can redistribute it and/or
055     * modify it under the terms of the GNU Lesser General Public
056     * License as published by the Free Software Foundation; either
057     * version 2.1 of the License, or (at your option) any later version.
058     * <br><br>
059     * This library is distributed in the hope that it will be useful,
060     * but WITHOUT ANY WARRANTY; without even the implied warranty of
061     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
062     * Lesser General Public License for more details.
063     * <br><br>
064     * You should have received a copy of the GNU Lesser General Public
065     * License along with this library; if not, write to the Free Software
066     * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
067     */
068    public class ExamSimulatedAnnealing implements NeighbourSelection, SolutionListener {
069        private static Logger sLog = Logger.getLogger(ExamSimulatedAnnealing.class);
070        private static DecimalFormat sDF2 = new DecimalFormat("0.00");
071        private static DecimalFormat sDF5 = new DecimalFormat("0.00000");
072        private static DecimalFormat sDF10 = new DecimalFormat("0.0000000000");
073        private double iInitialTemperature = 1.5;
074        private double iCoolingRate = 0.95;
075        private double iReheatRate = -1;   
076        private long iTemperatureLength = 250000;
077        private long iReheatLength = 0;
078        private long iRestoreBestLength = 0;
079        private double iTemperature = 0.0;
080        private double iReheatLengthCoef = 5.0;
081        private double iRestoreBestLengthCoef = -1;
082        private long iIter = 0;
083        private long iLastImprovingIter = 0;
084        private long iLastReheatIter = 0;
085        private long iLastCoolingIter = 0;
086        private int iAcceptIter[] = new int[] {0,0,0};
087        private boolean iStochasticHC = false;
088        private int iMoves = 0;
089        private double iAbsValue = 0;
090        private long iT0 = -1;
091        private double iBestValue = 0;
092        private Progress iProgress = null;
093        private boolean iActive;
094        
095        private NeighbourSelection[] iNeighbours = null;
096        
097        private boolean iRelativeAcceptance = true;
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         * </ul>
111         * @param properties problem properties
112         */
113        public ExamSimulatedAnnealing(DataProperties properties) {
114            iInitialTemperature = properties.getPropertyDouble("SimulatedAnnealing.InitialTemperature", iInitialTemperature);
115            iReheatRate = properties.getPropertyDouble("SimulatedAnnealing.ReheatRate", iReheatRate);
116            iCoolingRate = properties.getPropertyDouble("SimulatedAnnealing.CoolingRate", iCoolingRate);
117            iRelativeAcceptance = properties.getPropertyBoolean("SimulatedAnnealing.RelativeAcceptance", iRelativeAcceptance);
118            iStochasticHC = properties.getPropertyBoolean("SimulatedAnnealing.StochasticHC", iStochasticHC);
119            iTemperatureLength = properties.getPropertyLong("SimulatedAnnealing.TemperatureLength", iTemperatureLength);
120            iReheatLengthCoef = properties.getPropertyDouble("SimulatedAnnealing.ReheatLengthCoef", iReheatLengthCoef);
121            iRestoreBestLengthCoef = properties.getPropertyDouble("SimulatedAnnealing.RestoreBestLengthCoef", iRestoreBestLengthCoef);
122            if (iReheatRate<0) iReheatRate = Math.pow(1/iCoolingRate,iReheatLengthCoef*1.7);
123            if (iRestoreBestLengthCoef<0) iRestoreBestLengthCoef = iReheatLengthCoef * iReheatLengthCoef;
124            iNeighbours = new NeighbourSelection[] {
125                    new ExamRandomMove(properties),
126                    new ExamRoomMove(properties),
127                    new ExamTimeMove(properties)
128            };
129        }
130        
131        /**
132         * Initialization
133         */
134        public void init(Solver solver) {
135            iTemperature = iInitialTemperature;
136            iReheatLength = Math.round(iReheatLengthCoef*iTemperatureLength);
137            iRestoreBestLength = Math.round(iRestoreBestLengthCoef*iTemperatureLength);
138            solver.currentSolution().addSolutionListener(this);
139            for (int i=0;i<iNeighbours.length;i++)
140                iNeighbours[i].init(solver);
141            solver.setUpdateProgress(false);
142            iProgress = Progress.getInstance(solver.currentSolution().getModel());
143            iActive = false;
144        }
145        
146        /**
147         * Cool temperature
148         */
149        protected void cool(Solution solution) {
150            iTemperature *= iCoolingRate;
151            sLog.info("Iter="+iIter/1000+"k, NonImpIter="+sDF2.format((iIter-iLastImprovingIter)/1000.0)+"k, Speed="+sDF2.format(1000.0*iIter/(System.currentTimeMillis()-iT0))+" it/s");
152            sLog.info("Temperature decreased to "+sDF5.format(iTemperature)+" " +
153                    "(#moves="+iMoves+", rms(value)="+sDF2.format(Math.sqrt(iAbsValue/iMoves))+", "+
154                    "accept=-"+sDF2.format(100.0*iAcceptIter[0]/iTemperatureLength)+"/"+sDF2.format(100.0*iAcceptIter[1]/iTemperatureLength)+"/+"+sDF2.format(100.0*iAcceptIter[2]/iTemperatureLength)+"%, " +
155                    (prob(-1)<1.0?"p(-1)="+sDF2.format(100.0*prob(-1))+"%, ":"")+
156                    "p(+1)="+sDF2.format(100.0*prob(1))+"%, "+
157                    "p(+10)="+sDF5.format(100.0*prob(10))+"%)");
158            iLastCoolingIter=iIter;
159            iAcceptIter=new int[] {0,0,0};
160            iMoves = 0; iAbsValue = 0;
161        }
162        
163        /**
164         * Reheat temperature
165         */
166        protected void reheat(Solution solution) {
167            iTemperature *= iReheatRate;
168            sLog.info("Iter="+iIter/1000+"k, NonImpIter="+sDF2.format((iIter-iLastImprovingIter)/1000.0)+"k, Speed="+sDF2.format(1000.0*iIter/(System.currentTimeMillis()-iT0))+" it/s");
169            sLog.info("Temperature increased to "+sDF5.format(iTemperature)+" "+
170                    (prob(-1)<1.0?"p(-1)="+sDF2.format(100.0*prob(-1))+"%, ":"")+
171                    "p(+1)="+sDF2.format(100.0*prob(1))+"%, "+
172                    "p(+10)="+sDF5.format(100.0*prob(10))+"%, "+
173                    "p(+100)="+sDF10.format(100.0*prob(100))+"%)");
174            iLastReheatIter=iIter;
175            iProgress.setPhase("Simulated Annealing ["+sDF2.format(iTemperature)+"]...");
176        }
177        
178        /**
179         * restore best ever found solution 
180         */
181        protected void restoreBest(Solution solution) {
182            sLog.info("Best restored");
183            iLastImprovingIter=iIter;
184        }
185        
186        /**
187         * Generate neighbour -- select neighbourhood randomly, select neighbour
188         */
189        public Neighbour genMove(Solution solution) {
190            while (true) {
191                incIter(solution);
192                NeighbourSelection ns = iNeighbours[ToolBox.random(iNeighbours.length)];
193                Neighbour n = ns.selectNeighbour(solution);
194                if (n!=null) return n;
195            }
196        }
197        
198        /**
199         * Neighbour acceptance probability
200         * @param value absolute or relative value of the proposed change (neighbour)
201         * @return probability of acceptance of a change (neighbour)
202         */
203        protected double prob(double value) {
204            if (iStochasticHC)
205                return 1.0 / (1.0 + Math.exp(value/iTemperature));
206            else
207                return (value<=0.0?1.0:Math.exp(-value/iTemperature));
208        }
209        
210        /**
211         * True if the given neighboir is to be be accepted
212         * @param solution current solution
213         * @param neighbour proposed move
214         * @return true if generated random number is below {@link ExamSimulatedAnnealing#prob(double)}
215         */
216        protected boolean accept(Solution solution, Neighbour neighbour) {
217            double value = (iRelativeAcceptance?neighbour.value():solution.getModel().getTotalValue()+neighbour.value()-solution.getBestValue());
218            double prob = prob(value);
219            if (prob>=1.0 || ToolBox.random()<prob) {
220                iAcceptIter[neighbour.value()<0.0?0:neighbour.value()>0.0?2:1]++;
221                return true;
222            }
223            return false;
224        }
225        
226        /**
227         * Increment iteration counter, cool/reheat/restoreBest if necessary
228         */
229        protected void incIter(Solution solution) {
230            if (iT0<0) iT0 = System.currentTimeMillis();
231            iIter++;
232            if (iIter>iLastImprovingIter+iRestoreBestLength) restoreBest(solution);
233            if (iIter>Math.max(iLastReheatIter,iLastImprovingIter)+iReheatLength) reheat(solution);
234            if (iIter>iLastCoolingIter+iTemperatureLength) cool(solution);
235            iProgress.setProgress(Math.round(100.0 * (iIter-Math.max(iLastReheatIter,iLastImprovingIter)) / iReheatLength));
236        }
237        
238        /**
239         * Select neighbour -- generate a move {@link ExamSimulatedAnnealing#genMove(Solution)} until an acceptable neighbour
240         * is found {@link ExamSimulatedAnnealing#accept(Solution, Neighbour)}, keep increasing iteration 
241         * {@link ExamSimulatedAnnealing#incIter(Solution)}. 
242         */
243        public Neighbour selectNeighbour(Solution solution) {
244            if (!iActive) {
245                iProgress.setPhase("Simulated Annealing ["+sDF2.format(iTemperature)+"]...");
246                iActive = true;
247            }
248            Neighbour neighbour = null;
249            while ((neighbour=genMove(solution))!=null) {
250                iMoves++; iAbsValue += neighbour.value() * neighbour.value();
251                if (accept(solution,neighbour)) break;
252            }
253            if (neighbour==null) iActive = false;
254            return (neighbour==null?null:neighbour);
255        }
256        
257        /**
258         * Memorize the iteration when the last best solution was found.
259         */
260        public void bestSaved(Solution solution) {
261            if (Math.abs(iBestValue-solution.getBestValue())>=1.0) {
262                iLastImprovingIter = iIter;
263                iBestValue = solution.getBestValue();
264            }
265            iLastImprovingIter = iIter;
266        }
267        public void solutionUpdated(Solution solution) {}
268        public void getInfo(Solution solution, java.util.Dictionary info) {}
269        public void getInfo(Solution solution, java.util.Dictionary info, java.util.Vector variables) {}
270        public void bestCleared(Solution solution) {}
271        public void bestRestored(Solution solution){}    
272        
273    
274    
275    }