In [1]:
d1(x1, x2) = length(x1) == length(x2) ? sum(abs.(x1 .- x2)) :
error("Inputs not same length.")
D1(X₁, X₂) = [d1(i, j) for i in eachrow(X₁), j in eachrow(X₂)]
P = [50 150 220 295 420]'
r, f = 1, 1
w = r * f
k = [150, 200, 150, 150, 200]
C = w * D1(P, P)
Out[1]:
5×5 Matrix{Int64}:
0 100 170 245 370
100 0 70 145 270
170 70 0 75 200
245 145 75 0 125
370 270 200 125 0
In [2]:
# Find first NF to add: 3
N = 1:size(C, 1)
hcat(N, k, sum(C, dims=2), k + sum(C, dims=2))
Out[2]:
5×4 Matrix{Int64}:
1 150 885 1035
2 200 585 785
3 150 515 665
4 150 590 740
5 200 965 1165
In [3]:
y = [3]
setdiff(N, y)
Out[3]:
4-element Vector{Int64}:
1
2
4
5
In [4]:
# Find next NF to add: 1
for i in setdiff(N, y)
y′ = vcat(y, i)
[println(" ",r, ":",k[r],C[r, :]) for r in y′]
println(" ",minimum(C[y′, :], dims=1))
TC′ = sum(k[y′]) + sum(minimum(C[y′, :], dims=1))
println(i, " : ", sum(k[y′]), " + ", sum(minimum(C[y′, :], dims=1))," = ", TC′)
end
3:150[170, 70, 0, 75, 200]
1:150[0, 100, 170, 245, 370]
[0 70 0 75 200]
1 : 300 + 345 = 645
3:150[170, 70, 0, 75, 200]
2:200[100, 0, 70, 145, 270]
[100 0 0 75 200]
2 : 350 + 375 = 725
3:150[170, 70, 0, 75, 200]
4:150[245, 145, 75, 0, 125]
[170 70 0 0 125]
4 : 300 + 365 = 665
3:150[170, 70, 0, 75, 200]
5:200[370, 270, 200, 125, 0]
[170 70 0 75 0]
5 : 350 + 315 = 665
In [5]:
y = vcat(y, 1)
Out[5]:
2-element Vector{Int64}:
3
1
In [6]:
# Find next NF to add: none reduce TC, so STOP
for i in setdiff(N, y)
y′ = vcat(y, i)
[println(" ",r, ":",k[r],C[r, :]) for r in y′]
println(" ",minimum(C[y′, :], dims=1))
TC′ = sum(k[y′]) + sum(minimum(C[y′, :], dims=1))
println(i, " : ", sum(k[y′]), " + ", sum(minimum(C[y′, :], dims=1))," = ", TC′)
end
3:150[170, 70, 0, 75, 200]
1:150[0, 100, 170, 245, 370]
2:200[100, 0, 70, 145, 270]
[0 0 0 75 200]
2 : 500 + 275 = 775
3:150[170, 70, 0, 75, 200]
1:150[0, 100, 170, 245, 370]
4:150[245, 145, 75, 0, 125]
[0 70 0 0 125]
4 : 450 + 195 = 645
3:150[170, 70, 0, 75, 200]
1:150[0, 100, 170, 245, 370]
5:200[370, 270, 200, 125, 0]
[0 70 0 75 0]
5 : 500 + 145 = 645
In [7]:
# Write as a general procedure, given k and C as inputs
function ufladd(k, C)
fTC(y) = sum(k[y]) + sum(minimum(C[y, :], dims=1))
y = Int[]
TCᵒ, done = Inf, false
while !done
TC, i = Inf, nothing # Stops if y = all NF
for i′ = setdiff(1:size(C, 1), y) # since i′ = []
TC′ = fTC(vcat(y, i′))
if TC′ < TC
TC, i = TC′, i′
end
end
if TC < TCᵒ # TC = Inf if y = all NF
TCᵒ, y = TC, push!(y, i)
else
done = true
end
end
return y, TCᵒ
end
ufladd(k, C)
Out[7]:
([3, 1], 645)
UFL DROP Construction Procedure¶
Using same data as UFL ADD
In [8]:
# Calc TC with no dropped NFs
y = 1:size(C, 1)
TC = sum(k[y]) + sum(minimum(C[y, :], dims=1))
Out[8]:
850
In [9]:
# Find first NF to drop: 2
for i in y
y′ = setdiff(y, i)
[println(" ",r, ":",k[r],C[r, :]) for r in y′]
println(" ",minimum(C[y′, :], dims=1))
TC′ = sum(k[y′]) + sum(minimum(C[y′, :], dims=1))
println(i, " : ", sum(k[y′]), " + ", sum(minimum(C[y′, :], dims=1))," = ", TC′)
end
2:200[100, 0, 70, 145, 270]
3:150[170, 70, 0, 75, 200]
4:150[245, 145, 75, 0, 125]
5:200[370, 270, 200, 125, 0]
[100 0 0 0 0]
1 : 700 + 100 = 800
1:150[0, 100, 170, 245, 370]
3:150[170, 70, 0, 75, 200]
4:150[245, 145, 75, 0, 125]
5:200[370, 270, 200, 125, 0]
[0 70 0 0 0]
2 : 650 + 70 = 720
1:150[0, 100, 170, 245, 370]
2:200[100, 0, 70, 145, 270]
4:150[245, 145, 75, 0, 125]
5:200[370, 270, 200, 125, 0]
[0 0 70 0 0]
3 : 700 + 70 = 770
1:150[0, 100, 170, 245, 370]
2:200[100, 0, 70, 145, 270]
3:150[170, 70, 0, 75, 200]
5:200[370, 270, 200, 125, 0]
[0 0 0 75 0]
4 : 700 + 75 = 775
1:150[0, 100, 170, 245, 370]
2:200[100, 0, 70, 145, 270]
3:150[170, 70, 0, 75, 200]
4:150[245, 145, 75, 0, 125]
[0 0 0 0 125]
5 : 650 + 125 = 775
In [10]:
y = setdiff(y, 2)
Out[10]:
4-element Vector{Int64}:
1
3
4
5
In [11]:
# Find next NF to drop: 4 or 5, pick 4
for i in y
y′ = setdiff(y, i)
[println(" ",r, ":",k[r],C[r, :]) for r in y′]
println(" ",minimum(C[y′, :], dims=1))
TC′ = sum(k[y′]) + sum(minimum(C[y′, :], dims=1))
println(i, " : ", sum(k[y′]), " + ", sum(minimum(C[y′, :], dims=1))," = ", TC′)
end
3:150[170, 70, 0, 75, 200]
4:150[245, 145, 75, 0, 125]
5:200[370, 270, 200, 125, 0]
[170 70 0 0 0]
1 : 500 + 240 = 740
1:150[0, 100, 170, 245, 370]
4:150[245, 145, 75, 0, 125]
5:200[370, 270, 200, 125, 0]
[0 100 75 0 0]
3 : 500 + 175 = 675
1:150[0, 100, 170, 245, 370]
3:150[170, 70, 0, 75, 200]
5:200[370, 270, 200, 125, 0]
[0 70 0 75 0]
4 : 500 + 145 = 645
1:150[0, 100, 170, 245, 370]
3:150[170, 70, 0, 75, 200]
4:150[245, 145, 75, 0, 125]
[0 70 0 0 125]
5 : 450 + 195 = 645
In [12]:
y = setdiff(y, 4)
Out[12]:
3-element Vector{Int64}:
1
3
5
In [13]:
# Find next NF to add: none reduce TC, so STOP
for i in y
y′ = setdiff(y, i)
[println(" ",r, ":",k[r],C[r, :]) for r in y′]
println(" ",minimum(C[y′, :], dims=1))
TC′ = sum(k[y′]) + sum(minimum(C[y′, :], dims=1))
println(i, " : ", sum(k[y′]), " + ", sum(minimum(C[y′, :], dims=1))," = ", TC′)
end
3:150[170, 70, 0, 75, 200]
5:200[370, 270, 200, 125, 0]
[170 70 0 75 0]
1 : 350 + 315 = 665
1:150[0, 100, 170, 245, 370]
5:200[370, 270, 200, 125, 0]
[0 100 170 125 0]
3 : 350 + 395 = 745
1:150[0, 100, 170, 245, 370]
3:150[170, 70, 0, 75, 200]
[0 70 0 75 200]
5 : 300 + 345 = 645
Download file ufl_heuristics.jl posted along with this notebook, and then use include it to input all of its functions for use:
In [14]:
include("ufl_heuristics.jl")
ufldrop(k, C)
Out[14]:
([1, 3, 5], 645)
UFL Exchange Improvement Procedure¶
In [15]:
# First use UFLADD to construct a solution `y` that can then be improved
y, TC = ufladd(k, C)
uflxchg(k, C, y)
Out[15]:
([4, 1], 600)
UFL Hybrid Procedure¶
In [16]:
ufl(k, C)
Add: 645 Xchg: 600 Add: 600 Drop: 600
Out[16]:
([4, 1], 600)
p-Median Construction Procedure¶
In [17]:
p = 2
pmedian(p, C)
Out[17]:
([2, 5], 295)
2. Ex: NC & SC 10K+ Cities¶
In [18]:
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₂)]
function lonlat2name(X, D, name, st)
for i = 1:size(D, 1)
d = D[i, :]
println("NF", i, ": ", round(minimum(d)), " miles from ",
name[argmin(d)], ", ", st[argmin(d)])
end
end
Out[18]:
lonlat2name (generic function with 1 method)
In [19]:
using Logjam.DataTools, DataFrames
include("ufl_heuristics.jl")
df = filter!(r -> any(i -> r.STFIP == i, st2fips.([:NC,:SC])) &&
r.POP > 10_000, usplace())
select!(df, :NAME, :ST, :LAT, :LON, :POP)
P = hcat(df.LON, df.LAT)
w = df.POP
D = Dgc(P, P)
C = w[:]' .* D
Out[19]:
155×155 Matrix{Float64}:
0.0 9.52936e5 1.26403e7 1.93727e6 … 2.69232e6 1.87705e6
1.14833e6 0.0 1.92483e7 3.15753e6 4.13654e6 2.61596e6
2.19587e6 2.77484e6 0.0 1.27228e6 1.05608e6 2.38361e6
1.66736e6 2.25519e6 6.30335e6 0.0 2.13336e6 2.73321e6
1.05173e6 8.59114e5 1.6674e7 2.34719e6 3.88762e6 2.9804e6
1.3553e6 4.99219e5 1.98684e7 3.0905e6 … 4.45741e6 3.11635e6
1.24345e6 6.26805e5 1.87838e7 2.8265e6 4.27439e6 3.09141e6
6.16389e5 1.44203e6 9.47312e6 1.59848e6 1.88681e6 1.49064e6
4.12052e5 1.29487e6 1.02697e7 1.54549e6 2.19322e6 1.74832e6
1.38335e6 6.76204e5 1.95452e7 2.96233e6 4.45698e6 3.22952e6
3.7931e6 2.29626e6 3.3756e7 5.77613e6 … 7.64915e6 5.52991e6
1.17327e6 1.71006e5 1.93483e7 3.24807e6 4.09446e6 2.4837e6
9.29875e5 1.71426e6 7.57604e6 1.37237e6 1.49366e6 1.51704e6
⋮ ⋱
2.24114e6 2.76844e6 5.1776e6 2.09199e6 2.53149e5 1.68267e6
1.89325e6 2.16672e6 1.28539e7 3.04903e6 1.95162e6 1.52803e5
2.10641e6 1.40063e6 2.23322e7 4.30376e6 … 4.38312e6 2.09763e6
1.71327e6 2.01857e6 1.26076e7 2.91609e6 1.94286e6 73629.6
1.73462e6 2.04375e6 1.25174e7 2.91474e6 1.91444e6 79207.6
2.04826e6 2.59485e6 5.7289e6 2.05966e6 2.46599e5 1.45274e6
2.19798e6 1.48578e6 2.26852e7 4.39909e6 4.44418e6 2.12091e6
2.67704e6 2.33399e6 2.11185e7 4.53311e6 … 3.84864e6 1.50507e6
2.03388e6 2.60771e6 4.47688e6 1.85018e6 51639.7 1.63983e6
8.49583e5 1.60037e6 8.91089e6 1.70392e6 1.61541e6 1.25948e6
2.05941e6 2.62572e6 4.65012e6 1.89601e6 0.0 1.62772e6
1.771e6 2.04819e6 1.29458e7 2.99623e6 2.00773e6 0.0
In [20]:
# Check broadcasting:
w[1] * D[2, 1], w[2] * D[1, 2]
Out[20]:
(1.148331477137756e6, 952936.2233599344)
In [21]:
y, TC = ufladd(0, C, p = 2)
y, TC = uflxchg(0, C, y)
Out[21]:
([6, 8], 3.935223149761722e8)
In [22]:
y, TC = pmedian(2, C)
Out[22]:
([6, 8], 3.935223149761722e8)
In [23]:
lonlat2name(P[y, :], Dgc(P[y, :], hcat(df.LON, df.LAT)), df.NAME, df.ST)
NF1: 0.0 miles from Cary, NC NF2: 0.0 miles from Charlotte, NC
3. Ex: Best Retail Warehouse Locations in U.S.¶
In [24]:
using Logjam.DataTools, DataFrames
include("ufl_heuristics.jl")
df = filter!(r -> r.ISCUS == true && r.POP > 0, uszcta3())
select!(df, :ZCTA3, :LAT, :LON, :POP)
P = hcat(df.LON, df.LAT)
w = df.POP
D = Dgc(P, P)
C = w[:]' .* D
Out[24]:
882×882 Matrix{Float64}:
0.0 1.27581e6 4.80354e6 2.18512e6 … 9.2409e8 4.91286e7
3.51217e6 0.0 5.41746e6 2.77588e6 9.257e8 4.92168e7
1.77188e7 7.25901e6 0.0 3.21568e6 9.08958e8 4.8294e7
1.24675e7 5.75327e6 4.974e6 0.0 9.21586e8 4.89931e7
2.25424e7 8.77433e6 9.50444e6 2.97611e6 9.3507e8 4.97385e7
1.87714e7 6.75069e6 9.60909e6 3.61903e6 … 9.38799e8 4.99417e7
1.85466e7 6.78733e6 9.43232e6 3.40518e6 9.3798e8 4.98968e7
2.73191e7 9.9899e6 1.16637e7 4.61751e6 9.44071e8 5.0234e7
3.57961e7 1.33727e7 1.33056e7 5.41801e6 9.45415e8 5.0311e7
4.18636e7 1.54333e7 1.51848e7 6.66314e6 9.51715e8 5.06587e7
3.45029e7 1.23299e7 1.3883e7 6.19699e6 … 9.51841e8 5.06621e7
3.6616e7 1.33181e7 1.4167e7 6.16469e6 9.5125e8 5.0631e7
3.60784e7 1.31345e7 1.40081e7 6.05778e6 9.50721e8 5.06018e7
⋮ ⋱ ⋮
1.13975e9 4.14776e8 3.04181e8 1.99091e8 … 8.26229e7 5.98634e6
1.15063e9 4.18704e8 3.07132e8 2.01052e8 7.39824e7 5.77462e6
1.15008e9 4.18497e8 3.06983e8 2.00969e8 7.04601e7 5.62251e6
1.16344e9 4.23344e8 3.10605e8 2.03319e8 7.87306e7 6.13573e6
1.16173e9 4.22683e8 3.10148e8 2.03096e8 7.00666e7 5.72333e6
1.0916e9 3.97255e8 2.91128e8 1.90713e8 … 3.955e7 3.32021e6
1.10906e9 4.0357e8 2.95868e8 1.93833e8 2.8577e7 3.37339e6
1.03629e9 3.7716e8 2.76135e8 1.8103e8 5.01937e7 2.03863e6
1.04366e9 3.79836e8 2.78131e8 1.82319e8 4.632e7 2.07917e6
1.0363e9 3.77165e8 2.76138e8 1.81032e8 5.02776e7 2.04599e6
1.07985e9 3.92941e8 2.87952e8 1.88747e8 … 0.0 1.86568e6
1.04058e9 3.78672e8 2.77307e8 1.81874e8 3.38165e7 0.0
In [25]:
y, TC = pmedian(9, C)
Out[25]:
([834, 281, 873, 686, 312, 729, 559, 397, 60], 5.710476191147255e10)
In [26]:
df2 = filter!(r -> r.ISCUS == true && r.POP > 100_000, usplace())
lonlat2name(P[y, :], Dgc(P[y, :], hcat(df2.LON, df2.LAT)), df2.NAME, df2.ST)
NF1: 14.0 miles from Lancaster, CA NF2: 40.0 miles from Athens-Clarke County unified government (balance), GA NF3: 4.0 miles from Tacoma, WA NF4: 58.0 miles from Tyler, TX NF5: 14.0 miles from Lakeland, FL NF6: 5.0 miles from Denver, CO NF7: 17.0 miles from Rockford, IL NF8: 43.0 miles from Columbus, OH NF9: 3.0 miles from Newark, NJ
Using Integrated Formulation from previous notebook:
