Skip to content

Commit 0f8d2b2

Browse files
SebKrantzSebKrantz
committed
Making this example work.
1 parent 93e83fe commit 0f8d2b2

1 file changed

Lines changed: 73 additions & 65 deletions

File tree

‎examples/paper_example04.jl‎

Lines changed: 73 additions & 65 deletions
Original file line numberDiff line numberDiff line change
@@ -1,126 +1,134 @@
1-
1+
using OptimalTransportNetworks
22
using Random
33
using LinearAlgebra
44
using Plots
55

66
# Initialize parameters
7-
param = init_parameters(TolKappa=1e-5, K=100, LaborMobility=false)
7+
param = init_parameters(K = 100, labor_mobility = false)
88
width, height = 13, 13
99

1010
# Create graph
11-
param, g = create_graph(param, width, height, Type="map")
11+
param, g = create_graph(param, width, height, type = "map")
1212

1313
# Set fundamentals
1414
Random.seed!(5)
15-
param.N = 1
16-
param.Zjn = fill(param.minpop, g.J) # matrix of productivity
17-
param.Hj = ones(g.J) # matrix of housing
18-
param.Lj = zeros(g.J) # matrix of population
15+
minpop = minimum(param[:Lj])
16+
param[:Zjn] = fill(minpop, g[:J], 1) # matrix of productivity
17+
param[:Lj] = fill(1e-6, g[:J]) # matrix of population
1918

2019
Ni = find_node(g, 5, 5) # center
21-
param.Zjn[Ni] = 1 # more productive node
22-
param.Lj[Ni] = 1 # more productive node
20+
param[:Zjn][Ni] = 1 # more productive node
21+
param[:Lj][Ni] = 1 # more productive node
2322

2423
nb_cities = 20 # draw a number of random cities in space
2524
for i in 1:nb_cities-1
2625
newdraw = false
2726
while newdraw == false
28-
j = round(Int, 1 + rand() * (g.J - 1))
29-
if param.Lj[j] <= param.minpop
27+
j = round(Int, 1 + rand() * (g[:J] - 1))
28+
if param[:Lj][j] <= minpop
3029
newdraw = true
31-
param.Lj[j] = 1
30+
param[:Lj][j] = 1
3231
end
3332
end
3433
end
3534

36-
param.hj = param.Hj ./ param.Lj
37-
param.hj[param.Lj .== 0] = 1
35+
param[:hj] = param[:Hj] ./ param[:Lj]
36+
param[:hj][param[:Lj] .== 1e-6] .= 1
3837

3938
# Draw geography
40-
z = zeros(g.J) # altitude of each node
41-
obstacles = []
39+
z = zeros(g[:J]) # altitude of each node
4240

43-
geography = Dict("z" => z, "obstacles" => obstacles)
44-
g = apply_geography(g, geography)
41+
geographies = Dict()
42+
geographies[:blank] = (z = z, obstacles = nothing)
43+
g = apply_geography(g, geographies[:blank])
4544

4645
param0 = param # store initial params
4746
g0 = g # store initial graph
4847

4948
# Blank geography
50-
results = Dict(1 => optimal_network(param, g))
49+
results = Dict(:blank => optimal_network(param, g))
5150

5251
# Mountain
5352
mountain_size = 0.75
5453
mountain_height = 1
5554
mount_x = 10
5655
mount_y = 10
57-
geography[2] = geography[1]
58-
geography[2]["z"] = mountain_height * exp.(-((g.x .- mount_x).^2 .+ (g.y .- mount_y).^2) / (2 * mountain_size^2))
59-
60-
g = apply_geography(g, geography[2], AlphaUp_i=10, AlphaDown_i=10)
61-
results[2] = optimal_network(param, g)
56+
geographies[:mountain] = (z = mountain_height * exp.(-((g[:x] .- mount_x).^2 .+ (g[:y] .- mount_y).^2) / (2 * mountain_size^2)),
57+
obstacles = nothing)
58+
g = apply_geography(g, geographies[:mountain], alpha_up_i = 10, alpha_down_i = 10)
59+
results[:mountain] = optimal_network(param, g)
6260

6361
# Adding river and access by land
6462
g = g0
65-
geography[3] = geography[2]
66-
geography[3]["obstacles"] = [6 + (1-1)*width 6 + (2-1)*width;
67-
6 + (2-1)*width 6 + (3-1)*width;
68-
6 + (3-1)*width 7 + (4-1)*width;
69-
7 + (4-1)*width 8 + (5-1)*width;
70-
8 + (5-1)*width 9 + (5-1)*width;
71-
11 + (5-1)*width 12 + (5-1)*width;
72-
12 + (5-1)*width 13 + (5-1)*width]
73-
74-
g = apply_geography(g, geography[3], AcrossObstacleDelta_i=Inf, AlphaUp_i=10, AlphaDown_i=10)
75-
results[3] = optimal_network(param, g)
63+
geographies[:river] = (z = geographies[:mountain].z,
64+
obstacles = [6 + (1-1)*width 6 + (2-1)*width;
65+
6 + (2-1)*width 6 + (3-1)*width;
66+
6 + (3-1)*width 7 + (4-1)*width;
67+
7 + (4-1)*width 8 + (5-1)*width;
68+
8 + (5-1)*width 9 + (5-1)*width;
69+
11 + (5-1)*width 12 + (5-1)*width;
70+
12 + (5-1)*width 13 + (5-1)*width])
71+
72+
g = apply_geography(g, geographies[:river], across_abstacle_delta_i = Inf,
73+
alpha_up_i = 10, alpha_down_i = 10)
74+
results[:river] = optimal_network(param, g)
7675

7776
# Reinit and put another river and bridges
7877
g = g0
79-
geography[4] = geography[1]
80-
geography[4]["z"] = mountain_height * exp.(-((g.x .- mount_x).^2 .+ (g.y .- mount_y).^2) / (2 * mountain_size^2))
81-
82-
geography[4]["obstacles"] = [6 + (1-1)*width 6 + (2-1)*width;
78+
geographies[:bridges] = (z = mountain_height * exp.(-((g[:x] .- mount_x).^2 .+ (g[:y] .- mount_y).^2) / (2 * mountain_size^2)),
79+
obstacles = [6 + (1-1)*width 6 + (2-1)*width;
8380
6 + (2-1)*width 6 + (3-1)*width;
8481
6 + (3-1)*width 7 + (4-1)*width;
8582
7 + (4-1)*width 8 + (5-1)*width;
8683
8 + (5-1)*width 9 + (5-1)*width;
8784
9 + (5-1)*width 10 + (5-1)*width;
8885
10 + (5-1)*width 11 + (5-1)*width;
8986
11 + (5-1)*width 12 + (5-1)*width;
90-
12 + (5-1)*width 13 + (5-1)*width]
87+
12 + (5-1)*width 13 + (5-1)*width])
9188

92-
g = apply_geography(g, geography[4], AlphaUp_i=10, AlphaDown_i=10, AcrossObstacleDelta_i=2, AlongObstacleDelta_i=Inf)
93-
results[4] = optimal_network(param, g)
89+
g = apply_geography(g, geographies[:bridges], alpha_up_i = 10, alpha_down_i = 10,
90+
across_abstacle_delta_i = 2, along_obstacle_delta_i = Inf)
91+
results[:bridges] = optimal_network(param, g)
9492

9593
# Allowing for water transport
9694
g = g0
97-
geography[5] = geography[1]
98-
geography[5]["z"] = mountain_height * exp.(-((g.x .- mount_x).^2 .+ (g.y .- mount_y).^2) / (2 * mountain_size^2))
99-
geography[5]["obstacles"] = [6 + (1-1)*width 6 + (2-1)*width;
100-
6 + (2-1)*width 6 + (3-1)*width;
101-
6 + (3-1)*width 7 + (4-1)*width;
102-
7 + (4-1)*width 8 + (5-1)*width;
103-
8 + (5-1)*width 9 + (5-1)*width;
104-
9 + (5-1)*width 10 + (5-1)*width;
105-
10 + (5-1)*width 11 + (5-1)*width;
106-
11 + (5-1)*width 12 + (5-1)*width;
107-
12 + (5-1)*width 13 + (5-1)*width]
108-
109-
g = apply_geography(g, geography[5], AlphaUp_i=10, AlphaDown_i=10, AcrossObstacleDelta_i=2, AlongObstacleDelta_i=0.5)
110-
results[5] = optimal_network(param, g)
95+
geographies[:water_transport] = (z = mountain_height * exp.(-((g[:x] .- mount_x).^2 .+ (g[:y] .- mount_y).^2) / (2 * mountain_size^2)),
96+
obstacles = [6 + (1-1)*width 6 + (2-1)*width;
97+
6 + (2-1)*width 6 + (3-1)*width;
98+
6 + (3-1)*width 7 + (4-1)*width;
99+
7 + (4-1)*width 8 + (5-1)*width;
100+
8 + (5-1)*width 9 + (5-1)*width;
101+
9 + (5-1)*width 10 + (5-1)*width;
102+
10 + (5-1)*width 11 + (5-1)*width;
103+
11 + (5-1)*width 12 + (5-1)*width;
104+
12 + (5-1)*width 13 + (5-1)*width])
105+
106+
g = apply_geography(g, geographies[:water_transport], alpha_up_i = 10, alpha_down_i = 10,
107+
across_abstacle_delta_i = 2, along_obstacle_delta_i = 0.5)
108+
results[:water_transport] = optimal_network(param, g)
111109

112110
# Increasing returns to transport
113-
param.gamma = 2
114-
geography[6] = geography[5]
115-
results[6] = optimal_network(param, g)
111+
param[:gamma] = 2
112+
geographies[:increasing_returns] = geographies[:water_transport]
113+
results[:increasing_returns] = optimal_network(param, g)
116114

117115
# Plot results
118-
s = ["geography_blank", "geography_mountain", "geography_river", "geography_bridges", "geography_water_transport", "geography_increasing_returns"]
116+
simulation = ["blank", "mountain", "river", "bridges", "water_transport", "increasing_returns"]
119117
obstacles = ["off", "off", "on", "on", "on", "on"]
120-
121-
for j in 1:length(s)
122-
fig = plot_graph(param, g, results[j].Ijk, Geography=true, GeographyStruct=geography[j], Sizes=3*results[j].cj .* (param.Lj .> param.minpop) / maximum(results[j].cj), Mesh=true, MeshTransparency=0.2, Obstacles=obstacles[j], Shades=param.Zjn, MinEdgeThickness=1.5)
123-
savefig(fig, "$(s[j]).png")
118+
plots = Vector{Any}(undef, length(simulation))
119+
120+
i = 0
121+
for s in simulation
122+
i += 1
123+
plots[i] = plot_graph(g, results[Symbol(s)][:Ijk],
124+
geography = geographies[Symbol(s)], obstacles = obstacles[i] == "on",
125+
mesh = true, mesh_transparency = 0.2,
126+
node_sizes = results[Symbol(s)][:cj] .* (param[:Lj] .> minpop),
127+
node_shades = param[:Zjn],
128+
edge_min_thickness = 1.5)
129+
title!(plots[i], "Geography $(s)")
124130
end
125131

126-
# Please note that this translation assumes that the functions `init_parameters`, `create_graph`, `find_node`, `apply_geography`, `optimal_network`, and `plot_graph` have been previously defined in Julia with the same functionality as in the original Matlab code.
132+
# Combine plots
133+
final_plot = plot(plots..., layout = (2, 3), size = (3*400, 2*400))
134+
display(final_plot)

0 commit comments

Comments
 (0)