1. Two-Opt Improvment¶
In [1]:
r = [1:6; 1]
println(r)
for i = 1:length(r)-2
for j = i+2:min(length(r)-1, length(r)+i-3)
println(vcat(r[1:i], reverse(r[i+1:j]), r[j+1:end]))
end
end
[1, 2, 3, 4, 5, 6, 1] [1, 3, 2, 4, 5, 6, 1] [1, 4, 3, 2, 5, 6, 1] [1, 5, 4, 3, 2, 6, 1] [1, 2, 4, 3, 5, 6, 1] [1, 2, 5, 4, 3, 6, 1] [1, 2, 6, 5, 4, 3, 1] [1, 2, 3, 5, 4, 6, 1] [1, 2, 3, 6, 5, 4, 1] [1, 2, 3, 4, 6, 5, 1]
In [2]:
D = [
0 8 6 9 1 5
3 0 1 5 4 2
9 2 0 3 1 1
8 2 1 0 10 6
6 7 10 1 0 10
6 2 5 2 1 0]
segcost(r, C) = map(i -> C[r[i], r[i+1]], 1:length(r) - 1)
@show r
segcost(r, D)
r = [1, 2, 3, 4, 5, 6, 1]
Out[2]:
6-element Vector{Int64}:
8
1
3
10
10
6
In [3]:
rTCh(r) = sum(segcost(r, D))
rTCh(r)
Out[3]:
38
In [4]:
function twoopt_prt(r, rTCh)
rᵒ, TCᵒ = copy(r), rTCh(r)
println("*", TCᵒ, ": ", rᵒ)
done = false
while !done
done = true
for i = 1:length(r)-2
for j = i+2:min(length(r)-1, length(r)+i-3)
r′ = vcat(r[1:i], reverse(r[i+1:j]), r[j+1:end])
TC = rTCh(r′)
if TC < TCᵒ
rᵒ, TCᵒ = r′, TC
println("*", TCᵒ, ": ", rᵒ)
done = false
break
else
println(" ", TC, ": ", r′)
end
end
if done == false
break
end
end
r = copy(rᵒ)
end
return rᵒ, TCᵒ
end
twoopt_prt(r, rTCh)
*38: [1, 2, 3, 4, 5, 6, 1] 39: [1, 3, 2, 4, 5, 6, 1] *32: [1, 4, 3, 2, 5, 6, 1] *31: [1, 3, 4, 2, 5, 6, 1] 32: [1, 4, 3, 2, 5, 6, 1] 31: [1, 2, 4, 3, 5, 6, 1] *21: [1, 5, 2, 4, 3, 6, 1] 21: [1, 2, 5, 4, 3, 6, 1] 32: [1, 4, 2, 5, 3, 6, 1] 31: [1, 3, 4, 2, 5, 6, 1] *12: [1, 5, 4, 2, 3, 6, 1] 34: [1, 4, 5, 2, 3, 6, 1] 40: [1, 2, 4, 5, 3, 6, 1] 39: [1, 3, 2, 4, 5, 6, 1] 21: [1, 5, 2, 4, 3, 6, 1] 30: [1, 5, 3, 2, 4, 6, 1] 31: [1, 5, 6, 3, 2, 4, 1] 13: [1, 5, 4, 3, 2, 6, 1] 18: [1, 5, 4, 6, 3, 2, 1] 20: [1, 5, 4, 2, 6, 3, 1]
Out[4]:
([1, 5, 4, 2, 3, 6, 1], 12)
In [5]:
function twoopt(r, rTCh)
rᵒ, TCᵒ = copy(r), rTCh(r)
done = false
while !done
done = true
for i = 1:length(r)-2
for j = i+2:min(length(r)-1, length(r)+i-3)
r′ = vcat(r[1:i], reverse(r[i+1:j]), r[j+1:end])
TC = rTCh(r′)
if TC < TCᵒ
rᵒ, TCᵒ = r′, TC
done = false
break
end
end
if done == false
break
end
end
r = copy(rᵒ)
end
return rᵒ, TCᵒ
end
Out[5]:
twoopt (generic function with 1 method)
In [6]:
using Random
Random.seed!(2024)
randtour(n) = (tour = randperm(n); push!(tour, tour[1]))
@show r = randtour(6);
r = randtour(6) = [4, 5, 6, 3, 2, 1, 4]
Multi-start Random Tour Construction + Two-Opt Improvement¶
In [7]:
n, nruns = 6, 100
rᵒ, TCᵒ = nothing, Inf
for i = 1:nruns
r, TC = twoopt(randtour(n), rTCh)
if TC < TCᵒ
rᵒ, TCᵒ = r, TC
println(i, ": ", TCᵒ)
end
end
rᵒ, TCᵒ
1: 20 2: 13 5: 12 20: 9
Out[7]:
([2, 1, 5, 4, 3, 6, 2], 9)
2. The Tour of North Carolina¶
Determine the minimum great-circle distance tour that visits the population centers of each North Carolina county.
- Hard problem in theory, since $n = 100 \implies (n - 1)! \approx 9.33 \times 10^{155}$
- Easy in practice, using multi-start random tour construction and two-opt improvement since each final two-opt verification requires $n(n - 3)/2 = 4,850$ tour cost calculations.
In [8]:
using Logjam.DataTools, DataFrames, CSV
df = filter(r -> (r.STFIP == st2fips(:NC)), uscounty())
Out[8]:
100×10 DataFrame
75 rows omitted
| Row | STFIP | COFIP | NAME | ST | LAT | LON | POP | ALAND | AWATER | CBSA |
|---|---|---|---|---|---|---|---|---|---|---|
| Int64 | Int64 | String31 | String3 | Float64 | Float64 | Int64 | Float64 | Float64 | Int64? | |
| 1 | 37 | 1 | Alamance | NC | 36.0715 | -79.4132 | 171415 | 423.452 | 10.788 | 15500 |
| 2 | 37 | 3 | Alexander | NC | 35.8893 | -81.1935 | 36444 | 259.985 | 3.647 | 25860 |
| 3 | 37 | 7 | Anson | NC | 34.9741 | -80.1108 | 22055 | 531.461 | 5.638 | 16740 |
| 4 | 37 | 13 | Beaufort | NC | 35.5216 | -76.9618 | 44652 | 832.736 | 130.112 | 47820 |
| 5 | 37 | 19 | Brunswick | NC | 34.0418 | -78.2166 | 136693 | 850.083 | 199.46 | 48900 |
| 6 | 37 | 21 | Buncombe | NC | 35.5853 | -82.5518 | 269452 | 656.497 | 3.453 | 11700 |
| 7 | 37 | 23 | Burke | NC | 35.7313 | -81.6412 | 87570 | 506.237 | 8.0 | 25860 |
| 8 | 37 | 25 | Cabarrus | NC | 35.3957 | -80.6211 | 225804 | 361.232 | 2.697 | 16740 |
| 9 | 37 | 27 | Caldwell | NC | 35.873 | -81.5025 | 80652 | 471.888 | 2.723 | 25860 |
| 10 | 37 | 29 | Camden | NC | 36.3693 | -76.1977 | 10355 | 240.33 | 69.927 | 47260 |
| 11 | 37 | 31 | Carteret | NC | 34.7441 | -76.8187 | 67686 | 507.604 | 822.807 | 33980 |
| 12 | 37 | 35 | Catawba | NC | 35.6955 | -81.2453 | 160610 | 401.374 | 14.65 | 25860 |
| 13 | 37 | 37 | Chatham | NC | 35.756 | -79.2084 | 76285 | 681.679 | 27.252 | 20500 |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
| 89 | 37 | 113 | Macon | NC | 35.1569 | -83.3773 | 37014 | 515.579 | 4.094 | missing |
| 90 | 37 | 117 | Martin | NC | 35.8303 | -77.0973 | 22031 | 456.412 | 0.288 | missing |
| 91 | 37 | 121 | Mitchell | NC | 35.9601 | -82.114 | 14903 | 221.251 | 0.63 | missing |
| 92 | 37 | 123 | Montgomery | NC | 35.3406 | -79.8929 | 25751 | 491.537 | 9.995 | missing |
| 93 | 37 | 143 | Perquimans | NC | 36.1806 | -76.4229 | 13005 | 247.173 | 81.759 | missing |
| 94 | 37 | 149 | Polk | NC | 35.2584 | -82.1769 | 19328 | 237.687 | 0.758 | missing |
| 95 | 37 | 163 | Sampson | NC | 35.0367 | -78.3976 | 59036 | 945.93 | 1.901 | missing |
| 96 | 37 | 173 | Swain | NC | 35.4378 | -83.412 | 14117 | 527.728 | 12.516 | missing |
| 97 | 37 | 177 | Tyrrell | NC | 35.8953 | -76.2425 | 3245 | 390.777 | 206.399 | missing |
| 98 | 37 | 185 | Warren | NC | 36.4113 | -78.1345 | 18642 | 429.391 | 14.91 | missing |
| 99 | 37 | 187 | Washington | NC | 35.8632 | -76.6517 | 11003 | 346.51 | 75.399 | missing |
| 100 | 37 | 199 | Yancey | NC | 35.9109 | -82.2946 | 18470 | 312.592 | 0.589 | missing |
In [9]:
function dgc(xy₁, xy₂; unit=:mi)
length(xy₁) == length(xy₂) == 2 || error("Inputs must have length 2.")
unit in [:mi, :km] || error("Unit must be :mi or :km")
Δx, Δy = xy₂[1] - xy₁[1], xy₂[2] - xy₁[2]
a = sind(Δy / 2)^2 + cosd(xy₁[2]) * cosd(xy₂[2]) * sind(Δx / 2)^2
2 * asin(min(sqrt(a), 1.0)) * (unit == :mi ? 3958.75 : 6371.00)
end
Dgc(X₁, X₂) = [dgc(i, j) for i in eachrow(X₁), j in eachrow(X₂)]
XY = hcat(df.LON, df.LAT)
D = Dgc(XY, XY)
Out[9]:
100×100 Matrix{Float64}:
0.0 100.331 85.3701 142.532 … 75.0221 155.086 161.461
100.331 0.0 87.8249 238.757 174.423 254.254 61.6489
85.3701 87.8249 0.0 181.659 148.852 204.213 138.922
142.532 238.757 181.659 0.0 89.8832 29.3224 300.331
155.714 211.418 125.619 124.599 163.784 153.917 264.504
178.996 79.0234 144.007 314.218 … 253.407 331.457 26.7238
126.894 27.3643 100.868 263.178 201.387 279.737 38.653
82.2834 46.8642 40.9705 206.115 155.851 225.23 100.473
117.633 17.3402 99.9844 255.919 191.561 271.576 44.4137
180.398 280.722 239.834 72.504 107.757 43.1829 341.608
172.497 258.994 187.315 54.3249 … 136.877 77.893 318.985
105.804 13.7062 81.0749 240.904 180.663 257.726 60.6543
24.6289 111.588 74.1856 127.185 75.1373 143.444 173.201
⋮ ⋱
231.441 132.822 185.153 362.407 306.33 381.324 80.1265
130.598 229.393 179.721 22.6443 70.4439 25.057 291.015
151.131 51.7355 131.704 290.488 … 224.082 305.703 10.6611
57.2275 82.3016 28.1533 165.484 123.14 185.622 140.518
167.045 267.277 223.361 54.6251 96.636 25.3851 328.505
164.99 70.3883 118.408 294.278 240.011 313.309 45.5667
91.4934 167.999 97.0542 87.6399 96.1141 113.651 227.429
228.436 128.383 189.093 362.897 … 302.802 380.611 70.7225
177.677 277.093 226.857 47.9049 111.405 23.014 338.661
75.0221 174.423 148.852 89.8832 0.0 90.9914 234.608
155.086 254.254 204.213 29.3224 90.9914 0.0 315.85
161.461 61.6489 138.922 300.331 234.608 315.85 0.0
In [10]:
using CairoMakie, Logjam.MapTools
Random.seed!(2024)
r = randtour(nrow(df))
@show rTCh(r)
fig, ax = makemap(df.LON, df.LAT; xexpand=0.05)
lines!(ax, df.LON[r], df.LAT[r], color=:red, linewidth=.5)
fig
rTCh(r) = 15704.014343120863
Out[10]:
In [11]:
Random.seed!(2024)
n, nruns = nrow(df), 100
rᵒ, TCᵒ = nothing, Inf
for i = 1:nruns
r, TC = twoopt(randtour(n), rTCh)
if TC < TCᵒ
rᵒ, TCᵒ = r, TC
println(i, ": ", TCᵒ)
end
end
1: 2244.4562545646195 2: 2239.234111221034 5: 2206.2988301209016 8: 2193.59127803764 14: 2171.6384056131983 35: 2150.95477387751 58: 2141.2129756166273
In [12]:
fig, ax = makemap(df.LON, df.LAT; xexpand=0.05)
lines!(ax, df.LON[rᵒ], df.LAT[rᵒ], color=:red, linewidth=.5)
fig
Out[12]:
Expected Tour Distance¶
In [14]:
@show a = sum(df.ALAND) + sum(df.AWATER)
@show TCᴱ = 0.9sqrt(nrow(df)*a)
@show TCᵒ
TCᴱ/TCᵒ
a = sum(df.ALAND) + sum(df.AWATER) = 53818.558999999994 TCᴱ = 0.9 * sqrt(nrow(df) * a) = 2087.8944606947925 TCᵒ = 2141.2129756166273
Out[14]:
0.9750989203180594