Simulated annealing: Difference between revisions

Content added Content deleted
(added Perl 6 programming solution)
(→‎{{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>#!/usr/bin/env 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 ($ab, $cd) = ($j/100).floor, ($j%100).floor;
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


sub Es(@path) { # 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 $d
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 ($su, $sv) = s[u], s[v];
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
return exp( -ΔE / T(k, kmax, kT))
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;
say "kT =", kT;
(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: 0 T: 1.0 Es: 492.6876644465769
k: 200000 T: 0.8 Es: 166.1
k: 100000 T: 0.9 Es: 182.6232539415397
k: 300000 T: 0.7 Es: 174.7
k: 200000 T: 0.8 Es: 178.1973980739033
k: 400000 T: 0.6 Es: 146.9
k: 300000 T: 0.7 Es: 166.7332678035799
k: 500000 T: 0.5 Es: 140.2
k: 400000 T: 0.6 Es: 142.9907760779216
k: 600000 T: 0.4 Es: 127.5
k: 500000 T: 0.5 Es: 141.8232565352380
k: 700000 T: 0.3 Es: 115.9
k: 600000 T: 0.4 Es: 131.9340245844944
k: 800000 T: 0.2 Es: 111.9
k: 700000 T: 0.3 Es: 118.1083199652142
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) : 115.0995529529687
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 1 10 20 21 11 12 2 3 4 5 14 13 23 22 31 41 32 33 44
0 10 20 30 40 50 60 84 85 86 96 97 87 88 98 99 89 79 78 77
34 24 25 35 36 37 27 17 8 7 6 16 15 26 45 46 57 58 48 47
67 68 69 59 58 57 56 66 76 95 94 93 92 91 90 80 70 81 82 83
38 28 18 19 9 29 39 49 59 69 79 89 99 98 97 88 78 68 67 77
73 72 71 62 63 64 74 75 65 55 54 53 52 61 51 41 31 21 22 32
87 96 86 76 65 54 55 56 66 74 75 85 95 94 84 93 92 91 90 80
42 43 44 45 46 35 34 24 23 33 25 15 16 26 36 47 37 27 17 18
70 71 81 82 83 73 72 62 51 52 42 43 53 64 63 61 60 50 40 30
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}}==
Retrieved from "https://rosettacode.org/wiki/Simulated_annealing"