Simulated annealing: Difference between revisions
Content deleted Content added
added Perl 6 programming solution |
SqrtNegInf (talk | contribs) →{{header|Perl 6}}: city count configurable, code more idiomatic, even fewer sigils |
||
Line 646: | Line 646: | ||
=={{header|Perl 6}}== |
=={{header|Perl 6}}== |
||
{{trans|Go}} |
{{trans|Go}} |
||
<lang perl6># |
<lang perl6># simulation setup |
||
my \cities = 100; # number of cities |
|||
my \kmax = 1e6; # iterations to run |
|||
my \kT = 1; # initial 'temperature' |
|||
die 'City count must be a perfect square.' if cities.sqrt != cities.sqrt.Int; |
|||
my $dists = calcDists; |
|||
my \dirs = <1 -1 10 -10 9 11 -11 -9>; # all 8 neighbors |
|||
# locations of (up to) 8 neighbors, with grid size derived from number of cities |
|||
sub calcDists { # distances |
|||
my \gs = cities.sqrt; |
|||
loop (my @dists, my $j = 0; $j < 10000; $j++) { |
|||
my \neighbors = [1, -1, gs, -gs, gs-1, gs+1, -(gs+1), -(gs-1)]; |
|||
my ($a, $b) = ($ab/10).floor, ($ab.Int % 10).floor; |
|||
# matrix of distances between cities |
|||
my ($c, $d) = ($cd/10).floor, ($cd.Int % 10); |
|||
my \D = [;]; |
|||
@dists[$j] = sqrt(($a-$c)² + ($b-$d)²) |
|||
for 0 ..^ cities² -> \j { |
|||
} |
|||
my (\ab, \cd) = (j/cities, j%cities)».Int; |
|||
return @dists |
|||
my (\a, \b, \c, \d) = (ab/gs, ab%gs, cd/gs, cd%gs)».Int; |
|||
D[ab;cd] = sqrt (a - c)² + (b - d)² |
|||
} |
} |
||
sub T(\k, \kmax, \kT) { (1 - k/kmax) × kT } # temperature function, decreases to 0 |
|||
sub dist(\ci, \cj) { return @$dists[cj×100+ci] } # index into lookup table |
|||
sub P(\ΔE, \k, \kmax, \kT) { exp( -ΔE / T(k, kmax, kT)) } # probability to move from s to s_next |
|||
sub Es(\path) { sum (D[ path[$_]; path[$_+1] ] for 0 ..^ +path-1) } # energy at s, to be minimized |
|||
# variation of E, from state s to state s_next |
|||
sub delta-E(\s, \u, \v) { |
|||
loop (my $d = 0,my $i = 0; $i < +@path-1; $i++ ) { |
|||
my (\a, \b, \c, \d) = D[s[u-1];s[u]], D[s[u+1];s[u]], D[s[v-1];s[v]], D[s[v+1];s[v]]; |
|||
$d += dist @path[$i], @path[$i+1] |
|||
my (\na, \nb, \nc, \nd) = D[s[u-1];s[v]], D[s[u+1];s[v]], D[s[v-1];s[u]], D[s[v+1];s[u]]; |
|||
} |
|||
return |
if v == u+1 { return (na + nd) - (a + d) } |
||
elsif u == v+1 { return (nc + nb) - (c + b) } |
|||
else { return (na + nb + nc + nd) - (a + b + c + d) } |
|||
} |
} |
||
my \s = my @s = flat 0, (1 ..^ cities).pick(*), 0; |
|||
sub T(\k, \kmax, \kT) { # temperature function, decreases to 0 |
|||
return (1 - k/kmax) * kT |
|||
} |
|||
# E(s0), initial state |
|||
sub dE(\s, \u, \v) { # variation of E, from state s to state s_next |
|||
my |
my \E-min-global = my \E-min = my $s0 = Es(s); |
||
# old |
|||
my ($a, $b, $c, $d) = dist(s[u-1], $su), dist(s[u+1], $su), dist(s[v-1], $sv), dist(s[v+1], $sv); |
|||
# new |
|||
my ($na, $nb, $nc, $nd) = dist(s[u-1], $sv), dist(s[u+1], $sv), dist(s[v-1], $su), dist(s[v+1], $su); |
|||
if v == u+1 { |
|||
return ($na + $nd) - ($a + $d) |
|||
} elsif u == v+1 { |
|||
return ($nc + $nb) - ($c + $b) |
|||
} else { |
|||
return ($na + $nb + $nc + $nd) - ($a + $b + $c + $d) |
|||
} |
|||
} |
|||
for 0 ..^ kmax -> \k { |
|||
sub P(\ΔE, \k, \kmax, \kT) { # probability to move from s to s_next |
|||
printf "k:%8u T:%4.1f Es: %3.1f\n" , k, T(k, kmax, kT), Es(s) |
|||
if k % (kmax/10) == 0; |
|||
} |
|||
# valid candidate cities (exist, adjacent) |
|||
sub PrintPath(\p) { |
|||
my \cv = neighbors[(^8).roll] + s[ my \u = 1 + (^(cities-1)).roll ]; |
|||
say "Path: "; |
|||
next if cv ≤ 0 or cv ≥ cities or D[s[u];cv] > sqrt(2); |
|||
loop (my $i = 0; $i < +p; $i++) { |
|||
if $i > 0 and $i%20 == 0 { say " "; } |
|||
print " ", p[$i] |
|||
} |
|||
say " "; |
|||
} |
|||
my \v = s[cv]; |
|||
sub sa(\kmax, \kT) { |
|||
my \ΔE = delta-E(s, u, v); |
|||
srand(12345); |
|||
if ΔE < 0 or P(ΔE, k, kmax, kT) ≥ rand { # always move if negative |
|||
my @PathRecord = my @s = flat 0, (1..99).pick(99), 0; |
|||
(s[u], s[v]) = (s[v], s[u]); |
|||
E-min += ΔE; |
|||
say "E(s0) : ", Es(@s); # random starter |
|||
E-min-global = E-min if E-min < E-min-global; |
|||
my $EminRecord = my $Emin = Es(@s); # E0 |
|||
} |
|||
loop (my $k = 0; $k < kmax; $k++ ) { |
|||
if ($k%(kmax/10)) == 0 { |
|||
sprintf("k:%8u T:%4.1f Es: %8.13f",$k,T($k,kmax,kT),Es(@s)).put |
|||
} |
|||
my $u = 1 + (^99).roll; # city index 1 to 99 |
|||
my $cv = @s[$u] + dirs[(^8).roll]; # city number |
|||
next if $cv ≤ 0 or $cv ≥ 100 ; # bogus city |
|||
next if dist(@s[$u], $cv) > 5 ; # check true neighbor (eg 0 9) |
|||
my $v = @s[$cv]; # city index |
|||
my $ΔE = dE(@s, $u, $v); |
|||
if $ΔE < 0 or P($ΔE,$k,kmax,kT) ≥ 1.rand { #always move if negative |
|||
(@s[$u], @s[$v]) = (@s[$v], @s[$u]); |
|||
$Emin += $ΔE; |
|||
if $Emin < $EminRecord { |
|||
$EminRecord = $Emin; |
|||
@PathRecord = @s |
|||
} |
|||
} |
|||
} |
|||
say "\nE(s_final) : ", $Emin; |
|||
PrintPath @s; |
|||
say "Global optium : ",$EminRecord; |
|||
PrintPath @PathRecord; |
|||
} |
} |
||
say "\nE(s_final): " ~ E-min-global.fmt('%.1f'); |
|||
sa(1e6, 1)</lang> |
|||
say "Path:\n" ~ s».fmt('%2d').rotor(20,:partial).join: "\n";</lang> |
|||
{{out}} |
{{out}} |
||
<pre>k: 0 T: 1.0 Es: 522.0 |
|||
<pre>kT =1 |
|||
k: 100000 T: 0.9 Es: 185.3 |
|||
E(s0) : 492.6876644465769 |
|||
k: |
k: 200000 T: 0.8 Es: 166.1 |
||
k: |
k: 300000 T: 0.7 Es: 174.7 |
||
k: |
k: 400000 T: 0.6 Es: 146.9 |
||
k: |
k: 500000 T: 0.5 Es: 140.2 |
||
k: |
k: 600000 T: 0.4 Es: 127.5 |
||
k: |
k: 700000 T: 0.3 Es: 115.9 |
||
k: |
k: 800000 T: 0.2 Es: 111.9 |
||
k: |
k: 900000 T: 0.1 Es: 109.4 |
||
k: 800000 T: 0.2 Es: 113.0995529529691 |
|||
k: 900000 T: 0.1 Es: 115.0995529529691 |
|||
E(s_final) |
E(s_final): 109.4 |
||
Path: |
|||
0 1 10 20 21 11 12 2 3 4 15 14 13 23 22 31 41 32 33 44 |
|||
34 24 25 35 36 37 27 17 8 7 6 5 16 26 45 46 56 58 48 47 |
|||
38 28 18 19 9 29 39 49 59 69 79 89 99 98 97 88 78 68 67 77 |
|||
87 96 86 76 65 54 55 57 66 74 75 84 85 95 94 93 92 91 90 80 |
|||
70 71 81 82 83 73 72 61 51 52 42 43 53 64 63 62 60 50 40 30 |
|||
0 |
|||
Global optium : 113.0995529529687 |
|||
Path: |
Path: |
||
0 |
0 10 20 30 40 50 60 84 85 86 96 97 87 88 98 99 89 79 78 77 |
||
67 68 69 59 58 57 56 66 76 95 94 93 92 91 90 80 70 81 82 83 |
|||
73 72 71 62 63 64 74 75 65 55 54 53 52 61 51 41 31 21 22 32 |
|||
42 43 44 45 46 35 34 24 23 33 25 15 16 26 36 47 37 27 17 18 |
|||
28 38 48 49 39 29 19 9 8 7 6 5 4 14 13 12 11 2 3 1 |
|||
0 |
0</pre> |
||
</pre> |
|||
=={{header|Phix}}== |
=={{header|Phix}}== |