1# annealing.tcl --
2#     Package implementing simulated annealing for minimizing functions
3#     of one or more parameters
4#
5# Copyright (c) 2007 by Arjen Markus <arjenmarkus@users.sourceforge.net>
6#
7# See the file "license.terms" for information on usage and redistribution
8# of this file, and for a DISCLAIMER OF ALL WARRANTIES.
9#
10# RCS: @(#) $Id: annealing.tcl,v 1.4 2008/02/22 13:34:07 arjenmarkus Exp $
11#------------------------------------------------------------------------------
12
13package require Tcl 8.4
14
15# ::simulation::annealing --
16#     Create the namespace
17#
18namespace eval ::simulation::annealing {
19}
20
21
22# getOption --
23#     Return the value of an option
24#
25# Arguments:
26#     option         Name of the option (without -)
27# Result:
28#     The value or an error message if it does not exist
29#
30proc ::simulation::annealing::getOption {option} {
31    variable ann_option
32
33    if { [info exists ann_option(-$option)] } {
34        return $ann_option(-$option)
35    } else {
36        return -code error "No such option: $option"
37    }
38}
39
40
41# setOption --
42#     Set the value of an option
43#
44# Arguments:
45#     option         Name of the option (without -)
46#     value          Value of the option
47# Result:
48#     None
49#
50proc ::simulation::annealing::setOption {option value} {
51    variable ann_option
52
53    set ann_option(-$option) $value
54}
55
56
57# hasOption --
58#     Return whether the given option exists or not
59#
60# Arguments:
61#     option         Name of the option (without -)
62# Result:
63#     1 if it exists, 0 if not
64#
65proc ::simulation::annealing::hasOption {option} {
66    variable ann_option
67
68    if { [info exists ann_option(-$option)] } {
69        return 1
70    } else {
71        return 0
72    }
73}
74
75
76# findMinimum --
77#     Find the (global) minimum of a function using simulated annealing
78#
79# Arguments:
80#     args            Option-value pairs:
81#                     -parameters list - triples defining parameters and ranges
82#                     -function expr   - expression defining the function
83#                     -code body       - body of code to define the function
84#                                        (takes precedence over -function)
85#                                        should set the variable "result"
86#                     -init code       - code to be run at start up
87#                     -final code      - code to be run at the end
88#                     -trials n        - number of trials before reducing the temperature
89#                     -reduce factor   - reduce the temperature by this factor
90#                                        (between 0 and 1)
91#                     -initial-temp t  - initial temperature
92#                     -scale s         - scale of the function (order of
93#                                        magnitude of the values)
94#                     -estimate-scale y/n - estimate the scale (only if -scale not present)
95#                     -verbose y/n     - Turn verbose printing on (1) or off (0)
96#                     -reportfile file - Channel to write verbose output to
97#                     Any others can be used via the getOption procedure
98#                     in the body.
99#
100# Result:
101#     Estimated minimum and the parameters involved:
102#     function value param1 value param2 value ...
103#
104proc ::simulation::annealing::findMinimum {args} {
105    variable ann_option
106
107    #
108    # Handle the options
109    #
110    set ann_option(-parameters)     {}
111    set ann_option(-function)       {}
112    set ann_option(-code)           {}
113    set ann_option(-init)           {}
114    set ann_option(-final)          {}
115    set ann_option(-trials)         300
116    set ann_option(-reduce)         0.95
117    set ann_option(-initial-temp)   1.0
118    set ann_option(-scale)          {}
119    set ann_option(-estimate-scale) 0
120    set ann_option(-verbose)        0
121    set ann_option(-reportfile)     stdout
122
123    foreach {option value} $args {
124        set ann_option($option) $value
125    }
126
127    if { $ann_option(-scale) == {} } {
128        if { ! $ann_option(-estimate-scale) } {
129            set ann_option(-scale) 1.0
130        }
131    }
132
133    if { $ann_option(-code) != {} } {
134        set ann_option(-function) {}
135    }
136
137    if { $ann_option(-code) == {} && $ann_option(-function) == {} } {
138        return -code error "Neither code nor function given! Nothing to optimize"
139    }
140    if { $ann_option(-parameters) == {} } {
141        return -code error "No parameters given! Nothing to optimize"
142    }
143
144    if { $ann_option(-function) != {} } {
145        set ann_option(-code) "set result \[expr {$ann_option(-function)}\]"
146    }
147
148    #
149    # Create the procedure
150    #
151    proc FindMin {} [string map \
152        [list PARAMETERS $ann_option(-parameters) \
153              CODE  $ann_option(-code) \
154              INIT  $ann_option(-init) \
155              FINAL $ann_option(-final)] {
156        #
157        # Give all parameters a value
158        #
159        foreach {_param_ _min_ _max_} {PARAMETERS} {
160            set $_param_ $_min_
161        }
162
163        set _trials_ [getOption trials]
164        set _temperature_ [getOption initial-temp]
165        set _reduce_      [getOption reduce]
166        set _noparams_    [expr {[llength {PARAMETERS}]/3}]
167        set _verbose_     [getOption verbose]
168        set _reportfile_  [getOption reportfile]
169
170        INIT
171
172        #
173        # Estimate the scale
174        #
175        if { [getOption estimate-scale] == 1 } {
176            set _sum_ 0.0
177            for { set _trial_ 0 } { $_trial_ < $_trials_/3 } { incr _trial_ } {
178                set _randp_ [expr {3*int($_noparams_*rand())}]
179                set _param_ [lindex {PARAMETERS} $_randp_]
180                set _min_   [lindex {PARAMETERS} [expr {$_randp_+1}]]
181                set _max_   [lindex {PARAMETERS} [expr {$_randp_+2}]]
182                set _old_param_ [set $_param_]
183                set $_param_ [expr {$_min_ + rand()*($_max_-$_min_)}]
184
185                CODE
186
187                set _sum_  [expr {$_sum_ + abs($result)}]
188            }
189            set _scale_ [expr {3.0*$_sum_/$_trials_}]
190        } else {
191            set _scale_ [getOption scale]
192        }
193        if { $_verbose_ } {
194            puts $_reportfile_ "Scale value: $_scale_"
195        }
196
197        #
198        # Start the outer loop
199        #
200        set _changes_     1
201
202        #
203        # Get the initial value of the function
204        #
205        CODE
206        set _old_result_ $result
207
208        if { $_verbose_ } {
209            puts $_reportfile_ "Result -- Mean of accepted values -- % accepted"
210        }
211
212        while {1} {
213            set _sum_       $_old_result_
214            set _accepted_  1
215            for { set _trial_ 0 } { $_trial_ < $_trials_} { incr _trial_ } {
216                set _randp_ [expr {3*int($_noparams_*rand())}]
217                set _param_ [lindex {PARAMETERS} $_randp_]
218                set _min_   [lindex {PARAMETERS} [expr {$_randp_+1}]]
219                set _max_   [lindex {PARAMETERS} [expr {$_randp_+2}]]
220                set _old_param_ [set $_param_]
221                set $_param_ [expr {$_min_ + rand()*($_max_-$_min_)}]
222
223                CODE
224
225                #
226                # Accept the new solution?
227                #
228                set _rand_  [expr {rand()}]
229                if { log($_rand_) < -($result-$_old_result_)/($_scale_*$_temperature_) } {
230                    incr _changes_
231                    set _old_result_ $result
232                    set _sum_        [expr {$_sum_ + $result}]
233                    incr _accepted_
234                } else {
235                    set $_param_ $_old_param_
236                }
237            }
238
239            if { $_verbose_ } {
240                puts $_reportfile_ \
241                    [format "%.5g -- %.5g -- %.2f %%" $_old_result_ \
242                        [expr {$_sum_/$_accepted_}] [expr {100.0*double($_changes_)/$_trials_}]]
243            }
244
245            set _temperature_ [expr {$_reduce_ * $_temperature_}]
246            if { $_changes_ == 0 } {
247                break
248            } else {
249                set _changes_ 0
250            }
251        }
252
253        set result [list result $_old_result_] ;# Note: we need the last accepted result!
254        foreach {_param_ _min_ _max_} {PARAMETERS} {
255            lappend result $_param_ [set $_param_]
256        }
257
258        FINAL
259
260        return $result
261    }]
262
263    #
264    # Do the actual computation and return the result
265    #
266    return [FindMin]
267}
268
269
270# findCombinatorialMinimum --
271#     Find the (global) minimum of a combinatorial function using simulated annealing
272#
273# Arguments:
274#     args            Option-value pairs:
275#                     -number-params n     - number of (binary) parameters
276#                     -initial-values list - list of parameter values to start with
277#                     -function expr       - expression defining the function
278#                     -code body           - body of code to define the function
279#                                            (takes precedence over -function)
280#                                            should set the variable "result"
281#                                            The values of the solutions
282#                                            are stored as a list in the
283#                                            variable params
284#                     -init code           - code to be run at start up
285#                     -final code          - code to be run at the end
286#                     -trials n            - number of trials before reducing the temperature
287#                     -reduce factor       - reduce the temperature by this factor
288#                                            (between 0 and 1)
289#                     -initial-temp t      - initial temperature
290#                     -scale s             - scale of the function (order of
291#                                            magnitude of the values)
292#                     -estimate-scale y/n  - estimate the scale (only if -scale not present)
293#                     -verbose y/n         - Turn verbose printing on (1) or off (0)
294#                     -reportfile file     - Channel to write verbose output to
295#                     Any others can be used via the getOption procedure
296#                     in the body.
297#
298# Result:
299#     Estimated minimum and the parameters involved:
300#     function value, list of values
301#
302# Note:
303#     The parameters have the values 0 or 1
304#
305#     The stop criterion is that if the result value does not change in
306#     sqrt(trials) then the iteration stops. Experiments with the
307#     example below show that the function to be minimised can show
308#     a very wide minimum due to the parameters being discrete.
309#     sqrt(trials) is just an arbitrary value.
310#
311proc ::simulation::annealing::findCombinatorialMinimum {args} {
312    variable ann_option
313
314    #
315    # Handle the options
316    #
317    set ann_option(-number-params)  {}
318    set ann_option(-initial-values) {}
319    set ann_option(-function)       {}
320    set ann_option(-code)           {}
321    set ann_option(-init)           {}
322    set ann_option(-final)          {}
323    set ann_option(-trials)         300
324    set ann_option(-reduce)         0.95
325    set ann_option(-initial-temp)   1.0
326    set ann_option(-scale)          {}
327    set ann_option(-estimate-scale) 0
328    set ann_option(-verbose)        0
329    set ann_option(-reportfile)     stdout
330
331    foreach {option value} $args {
332        set ann_option($option) $value
333    }
334
335    if { $ann_option(-scale) == {} } {
336        if { ! $ann_option(-estimate-scale) } {
337            set ann_option(-scale) 1.0
338        }
339    }
340
341    if { $ann_option(-code) != {} } {
342        set ann_option(-function) {}
343    }
344
345    if { $ann_option(-code) == {} && $ann_option(-function) == {} } {
346        return -code error "Neither code nor function given! Nothing to optimize"
347    }
348    if { $ann_option(-number-params) == {} } {
349        return -code error "Number of parameters not given! Nothing to optimize"
350    }
351
352    if { $ann_option(-initial-values) == {} } {
353        for { set i 0 } { $i < $ann_option(-number-params) } { incr i } {
354            lappend ann_option(-initial-values) 0
355        }
356    }
357
358    if { $ann_option(-function) != {} } {
359        set ann_option(-code) "set result \[expr {$ann_option(-function)}\]"
360    }
361
362    #
363    # Create the procedure
364    #
365    proc FindCombMin {params} [string map \
366        [list CODE  $ann_option(-code) \
367              INIT  $ann_option(-init) \
368              FINAL $ann_option(-final)] {
369
370        set _trials_      [getOption trials]
371        set _temperature_ [getOption initial-temp]
372        set _reduce_      [getOption reduce]
373        set _noparams_    [llength $params]
374        set _verbose_     [getOption verbose]
375        set _reportfile_  [getOption reportfile]
376
377        INIT
378
379        #
380        # Estimate the scale
381        #
382        if { [getOption estimate-scale] == 1 } {
383            set _sum_ 0.0
384            set _old_params_ $params
385            for { set _trial_ 0 } { $_trial_ < $_trials_/3 } { incr _trial_ } {
386                set _randp_ [expr {int($_noparams_*rand())}]
387                lset params $_randp_ [expr {rand()>0.5? 0 : 1}]
388
389                CODE
390
391                set _sum_  [expr {$_sum_ + abs($result)}]
392            }
393            set _scale_ [expr {3.0*$_sum_/$_trials_}]
394            set params $_old_params_
395        } else {
396            set _scale_ [getOption scale]
397        }
398        if { $_verbose_ } {
399            puts $_reportfile_ "Scale value: $_scale_"
400        }
401
402        #
403        # Start the outer loop
404        #
405        set _changes_     1
406
407        #
408        # Get the initial value of the function
409        #
410        CODE
411        set _old_result_        $result
412        set _result_same_       0
413        set _result_after_loop_ $result
414
415        if { $_verbose_ } {
416            puts $_reportfile_ "Result -- Mean of accepted values -- % accepted"
417        }
418
419        while {1} {
420            set _sum_       $_old_result_
421            set _accepted_  1
422            for { set _trial_ 0 } { $_trial_ < $_trials_} { incr _trial_ } {
423                set _old_params_ $params
424                set _randp_ [expr {int($_noparams_*rand())}]
425                lset params $_randp_ [expr {rand()>0.5? 0 : 1}]
426
427                CODE
428
429                #
430                # Accept the new solution?
431                #
432                set _rand_  [expr {rand()}]
433                if { log($_rand_) < -($result-$_old_result_)/($_scale_*$_temperature_) } {
434                    incr _changes_
435                    set _old_result_ $result
436                    set _sum_        [expr {$_sum_ + $result}]
437                    incr _accepted_
438                } else {
439                    set params $_old_params_
440                }
441            }
442
443            if { $_verbose_ } {
444                puts $_reportfile_ \
445                    [format "%.5g -- %.5g -- %.2f %%" $_old_result_ \
446                        [expr {$_sum_/$_accepted_}] [expr {100.0*double($_changes_)/$_trials_}]]
447            }
448
449            set _temperature_ [expr {$_reduce_ * $_temperature_}]
450            if { $_changes_ == 0 || $_result_same_ > sqrt($_trials_) } {
451                break
452            } else {
453                set _changes_ 0
454            }
455
456            if { $_result_after_loop_ == $_old_result_ } {
457                incr _result_same_
458            } else {
459                set _result_after_loop_ $_old_result_
460            }
461        }
462
463        set result [list result $_old_result_] ;# Note: we need the last accepted result!
464        lappend result solution $params
465
466        FINAL
467
468        return $result
469    }]
470
471    #
472    # Do the actual computation and return the result
473    #
474    return [FindCombMin $ann_option(-initial-values)]
475}
476
477# Announce the package
478#
479package provide simulation::annealing 0.2
480
481# main --
482#     Example
483#
484if { 0 } {
485puts [::simulation::annealing::findMinimum \
486    -trials 300 \
487    -verbose 1 \
488    -parameters {x -5.0 5.0 y -5.0 5.0} \
489    -function {$x*$x+$y*$y+sin(10.0*$x)+4.0*cos(20.0*$y)}]
490
491puts "Constrained:"
492puts [::simulation::annealing::findMinimum \
493    -trials 3000 \
494    -reduce 0.98 \
495    -parameters {x -5.0 5.0 y -5.0 5.0} \
496    -code {
497        if { hypot($x-5.0,$y-5.0) < 4.0 } {
498            set result [expr {$x*$x+$y*$y+sin(10.0*$x)+4.0*cos(20.0*$y)}]
499        } else {
500           set result 1.0e100
501        }
502    }]
503}
504
505#
506# A simple combinatorial problem:
507# We have 100 items and the function is optimal if the first 10
508# values are 1 and the result is 0. Can we find this solution?
509#
510# What if we have 1000 items? Or 10000 items?
511#
512# WARNING:
513# 10000 items take a very long time!
514#
515if { 0 } {
516    proc cost {params} {
517        set cost 0
518        foreach p [lrange $params 0 9] {
519            if { $p == 0 } {
520                incr cost
521            }
522        }
523        foreach p [lrange $params 10 end] {
524            if { $p == 1 } {
525                incr cost
526            }
527        }
528        return $cost
529    }
530
531    foreach n {100 1000 10000} {
532        break
533        puts "Problem size: $n"
534        puts [::simulation::annealing::findCombinatorialMinimum \
535            -trials 300 \
536            -verbose 0 \
537            -number-params $n \
538            -code {set result [cost $params]}]
539    }
540
541    #
542    # Second problem:
543    #     Only the values of the first 10 items are important -
544    #     they should be 1
545    #
546    proc cost2 {params} {
547        set cost 0
548        foreach p [lrange $params 0 9] {
549            if { $p == 0 } {
550                incr cost
551            }
552        }
553        return $cost
554    }
555
556    foreach n {100 1000 10000} {
557        puts "Problem size: $n"
558        puts [::simulation::annealing::findCombinatorialMinimum \
559            -trials 300 \
560            -verbose 0 \
561            -number-params $n \
562            -code {set result [cost2 $params]}]
563    }
564}
565