Simulated annealing: Difference between revisions
→{{header|Phix}}
(Added Julia language) |
|||
Line 386:
76, 75, 84, 83, 93, 92, 91, 100, 90, 11
1</pre>
=={{header|Phix}}==
{{trans|zkl}}
Note that the standard builtin exp() suffered occasional overflows, so this uses b_a_exp() from bigatom.e, but
it does make it much slower.
<lang Phix>function hypot(atom a,b) return sqrt(a*a+b*b) end function
function calc_dists()
sequence dists = repeat(0,10000)
for abcd=1 to 10000 do
integer {ab,cd} = {floor(abcd/100),mod(abcd,100)},
{a,b,c,d} = {floor(ab/10),mod(ab,10),
floor(cd/10),mod(cd,10)}
dists[abcd] = hypot(a-c,b-d)
end for
return dists
end function
constant dists = calc_dists()
function dist(integer ci,cj) return dists[cj*100+ci] end function
function Es(sequence path)
atom d = 0
for i=1 to length(path)-1 do
d += dist(path[i],path[i+1])
end for
return d
end function
-- temperature() function
function T(integer k, kmax, kT) return (1-k/kmax)*kT end function
include bigatom.e -- (just for b_a_exp())
-- deltaE = Es_new - Es_old > 0
-- probability to move if deltaE > 0, -->0 when T --> 0 (frozen state)
function P(atom deltaE, integer k, kmax, kT) return b_a_exp(-deltaE/T(k,kmax,kT)) end function
-- deltaE from path ( .. a u b .. c v d ..) to (.. a v b ... c u d ..)
function dE(sequence s, integer u,v)
-- (note that u,v are 0-based, but 1..99 here)
-- integer sum1 = s[u-1], su = s[u], sup1 = s[u+1],
-- svm1 = s[v-1], sv = s[v], svp1 = s[v+1]
integer sum1 = s[u], su = s[u+1], sup1 = s[u+2],
svm1 = s[v], sv = s[v+1], svp1 = s[v+2]
-- old
atom {a,b,c,d}:={dist(sum1,su), dist(su,sup1), dist(svm1,sv), dist(sv,svp1)},
-- new
{na,nb,nc,nd}:={dist(sum1,sv), dist(sv,sup1), dist(svm1,su), dist(su,svp1)}
return iff(v==u+1?(na+nd)-(a+d):
iff(u==v+1?(nc+nb)-(c+b):
(na+nb+nc+nd)-(a+b+c+d)))
end function
-- all 8 neighbours
constant dirs = {1, -1, 10, -10, 9, 11, -11, -9}
procedure sa(integer kmax, kT=10)
sequence s = 0&shuffle(tagset(99))&0
atom Emin:=Es(s) -- E0
printf(1,"E(s0) %f\n",Emin) -- random starter
for k=0 to kmax do
if mod(k,kmax/10)=0 then
printf(1,"k:%,10d T: %8.4f Es: %8.4f\n",{k,T(k,kmax,kT),Es(s)})
end if
integer u = rand(99), -- city index 1 99
cv = s[u+1]+dirs[rand(8)] -- city number
if cv>0 and cv<100 -- not bogus city
and dist(s[u+1],cv)<5 then -- and true neighbour
integer v = s[cv+1] -- city index
atom deltae := dE(s,u,v);
if deltae<0 -- always move if negative
or P(deltae,k,kmax,kT)>=rnd() then
{s[u+1],s[v+1]} = {s[v+1],s[u+1]}
Emin += deltae
end if
end if
end for
printf(1,"E(s_final) %f\n",Emin)
printf(1,"Path:\n")
pp(s,{pp_IntFmt,"%2d",pp_StrFmt,-2})
end procedure
sa(1_000_000,1)</lang>
{{out}}
<pre>
E(s0) 515.164811
k: 0 T: 1.0000 Es: 515.1648
k: 100,000 T: 0.9000 Es: 189.3123
k: 200,000 T: 0.8000 Es: 198.7498
k: 300,000 T: 0.7000 Es: 158.2189
k: 400,000 T: 0.6000 Es: 165.4813
k: 500,000 T: 0.5000 Es: 156.3467
k: 600,000 T: 0.4000 Es: 142.7928
k: 700,000 T: 0.3000 Es: 128.0352
k: 800,000 T: 0.2000 Es: 121.7794
k: 900,000 T: 0.1000 Es: 121.2328
k: 1,000,000 T: 0.0000 Es: 121.1291
E(s_final) 121.129115
Path:
{ 0,10,62,63,64,65,76,75,84,85,95,86,96,97,87,77,67,66,56,46,47,48,49,59,69,
79,89,99,98,88,78,68,58,57,37,38,27,26,36,35,45,55,54,53,52,43,33,23,22,32,
42,41,51,61,60,50,40,30,31,21,20,11,12, 2, 3, 4, 5, 6,17,18,28,39,29,19, 9,
8, 7,16,15,24,44,74,83,93,94,92,91,71,70,90,80,81,82,72,73,34,25,14,13, 1,
0}
</pre>
=={{header|zkl}}==
| |||