|
=={{header|Perl 6}}==
{{trans|Go}}
<lang perl6>#!/usr/bin/env perl6simulation 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, my-gs, ($abgs-1, $cd) =gs+1, -($j/100gs+1).floor, -($j%100gs-1).floor];
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)variation {of #E, from energy atstate s, to bestate minimizeds_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]];
}
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 \E-min-global = my ($su,\E-min = my $sv)s0 = Es(s[u], s[v]);
# 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(printf -ΔE"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 (s[u], s[v]) =" (s[v], kTs[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}}
<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: 0200000 T: 10.08 Es: 492166.68766444657691
k: 100000300000 T: 0.97 Es: 182174.62325394153977
k: 200000400000 T: 0.86 Es: 178146.19739807390339
k: 300000500000 T: 0.75 Es: 166140.73326780357992
k: 400000600000 T: 0.64 Es: 142127.99077607792165
k: 500000700000 T: 0.53 Es: 141115.82325653523809
k: 600000800000 T: 0.42 Es: 131111.93402458449449
k: 700000900000 T: 0.31 Es: 118109.10831996521424
k: 800000 T: 0.2 Es: 113.0995529529691
k: 900000 T: 0.1 Es: 115.0995529529691
E(s_final) : 115109.09955295296874
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:
0 1 10 20 2130 40 1150 1260 284 385 486 596 1497 1387 2388 2298 3199 4189 3279 3378 4477
3467 2468 2569 3559 3658 3757 2756 1766 876 795 694 1693 1592 2691 4590 4680 5770 5881 4882 4783
3873 2872 1871 1962 963 2964 3974 4975 5965 6955 7954 8953 9952 9861 9751 8841 7831 6821 6722 7732
8742 9643 8644 7645 6546 5435 5534 5624 6623 7433 7525 8515 9516 9426 8436 9347 9237 9127 9017 8018
7028 7138 8148 8249 8339 7329 7219 62 519 52 428 43 537 64 636 61 605 50 404 3014 13 12 11 2 3 1
0</pre>
</pre>
=={{header|Phix}}==
| 