Update examples.

Author: SebKrantz

examples/example01.jl

@@ -17,13 +17,13 @@ param = init_parameters(labor_mobility = true, K = 10, gamma = 1, beta = 1, verb
 # ------------
 # Init network
 
-param, graph = create_graph(param, 11, 11, type = "map") # create a map network of 11x11 nodes located in [0,10]x[0,10]
+graph = create_graph(param, 11, 11, type = "map") # create a map network of 11x11 nodes located in [0,10]x[0,10]
 # note: by default, productivity and population are equalized everywhere
 
 # Customize graph
-param[:Zjn] = fill(0.1, param[:J], param[:N]) # set most places to low productivity
+graph[:Zjn] = fill(0.1, graph[:J], param[:N]) # set most places to low productivity
 Ni = find_node(graph, 6, 6) # Find index of the central node at (6,6)
-param[:Zjn][Ni, :] .= 1 # central node more productive
+graph[:Zjn][Ni, :] .= 1 # central node more productive
 
 # ==========
 # RESOLUTION

examples/example02.jl

@@ -13,31 +13,31 @@ param = init_parameters(labor_mobility = false, K = 100, gamma = 1, beta = 1, N
 # Init network
 
 w, h = 13, 13
-param, graph = create_graph(param, w, h, type = "triangle") # create a triangular network of 21x21
+graph = create_graph(param, w, h, type = "triangle") # create a triangular network of 21x21
 
 ncities = 20 # number of random cities
-param[:Lj] .= 0 # set population to minimum everywhere
-param[:Zjn] .= 0.01 # set low productivity everywhere
+graph[:Lj] .= 0 # set population to minimum everywhere
+graph[:Zjn] .= 0.01 # set low productivity everywhere
 
 Ni = find_node(graph, ceil(Int, w/2), ceil(Int, h/2)) # Find the central node
-param[:Zjn][Ni] = 1 # central node is more productive
-param[:Lj][Ni] = 1 # cities are equally populated
+graph[:Zjn][Ni] = 1 # central node is more productive
+graph[:Lj][Ni] = 1 # cities are equally populated
 
 # Draw the rest of the cities randomly
 Random.seed!(5) # reinit random number generator
 for i in 1:ncities
     newdraw = false
     while !newdraw
         j = round(Int, 1 + rand() * (graph[:J] - 1))
-        if param[:Lj][j] != 1 / (ncities + 1) # make sure node j is unpopulated           
+        if graph[:Lj][j] != 1 / (ncities + 1) # make sure node j is unpopulated           
             newdraw = true
-            param[:Lj][j] = 1
+            graph[:Lj][j] = 1
         end
     end
 end
 
 # For Ipopt: population cannot be zero!
-param[:Lj][param[:Lj] .== 0] .= 1e-6
+graph[:Lj][graph[:Lj] .== 0] .= 1e-6
 
 # ==========
 # RESOLUTION
@@ -50,7 +50,7 @@ param[:Lj][param[:Lj] .== 0] .= 1e-6
     # second, in the nonconvex case (gamma>beta)
     param[:gamma] = 2 # change only gamma, keep other parameters
     res[2] = optimal_network(param, graph) # optimize by iterating on FOCs
-    res[3] = annealing(param, graph, res[2][:Ijk]) # improve with annealing, starting from previous result
+    res[3] = annealing(graph, res[2][:Ijk]) # improve with annealing, starting from previous result
 end
 
 # ==========

examples/example03.jl

@@ -15,21 +15,21 @@ param = init_parameters(labor_mobility = true, K = 100, gamma = 2, beta = 1, N =
 # ------------
 # Init network
 
-param, graph = create_graph(param, 7, 7, type = "triangle") # create a triangular network of 21x21
+graph = create_graph(param, 7, 7, type = "triangle") # create a triangular network of 21x21
 
-param[:Zjn][:, 1:param[:N]-1] .= 0 # default locations cannot produce goods 1-10
-param[:Zjn][:, param[:N]] .= 1 # but they can all produce good 11 (agricultural)
+graph[:Zjn][:, 1:param[:N]-1] .= 0 # default locations cannot produce goods 1-10
+graph[:Zjn][:, param[:N]] .= 1 # but they can all produce good 11 (agricultural)
 
 # Draw the cities randomly
 Random.seed!(5) # reinit random number generator
 for i in 1:param[:N]-1
     newdraw = false
     while newdraw == false
         j = round(Int, 1 + rand() * (graph[:J] - 1))
-        if any(param[:Zjn][j, 1:param[:N]-1] .> 0) == false # make sure node j does not produce any differentiated good
+        if any(graph[:Zjn][j, 1:param[:N]-1] .> 0) == false # make sure node j does not produce any differentiated good
             newdraw = true
-            param[:Zjn][j, 1:param[:N]] .= 0
-            param[:Zjn][j, i] = 1
+            graph[:Zjn][j, 1:param[:N]] .= 0
+            graph[:Zjn][j, i] = 1
         end
     end
 end

examples/example04.jl

@@ -11,33 +11,33 @@ import Random
 param = init_parameters(K = 100, labor_mobility = false,
                         N = 1, gamma = 1, beta = 1, duality = false)
 w = 13; h = 13
-param, graph = create_graph(param, w, h, type = "map")
+graph = create_graph(param, w, h, type = "map")
 
 # ----------------
 # Draw populations
 
-param[:Zjn] = ones(param[:J], 1) .* 1e-6  # matrix of productivity (not 0 to avoid numerical glitches)
-param[:Lj] = ones(param[:J]) .* 1e-6  # matrix of population
+graph[:Zjn] = ones(graph[:J], 1) .* 1e-6  # matrix of productivity (not 0 to avoid numerical glitches)
+graph[:Lj] = ones(graph[:J]) .* 1e-6  # matrix of population
 
 Ni = find_node(graph, ceil(w/2), ceil(h/2))  # center
-param[:Zjn][Ni] = 1  # more productive node
-param[:Lj][Ni] = 1  # more productive node
+graph[:Zjn][Ni] = 1  # more productive node
+graph[:Lj][Ni] = 1  # more productive node
 
 Random.seed!(5)
 ncities = 20  # draw a number of random cities in space
 for i in 1:ncities-1
     newdraw = false
     while newdraw == false
-        j = round(Int, 1 + rand() * (param[:J] - 1))
-        if param[:Lj][j] <= 1e-6
+        j = round(Int, 1 + rand() * (graph[:J] - 1))
+        if graph[:Lj][j] <= 1e-6
             newdraw = true
-            param[:Lj][j] = 1
+            graph[:Lj][j] = 1
         end
     end
 end
 
-param[:hj] = param[:Hj] ./ param[:Lj]
-param[:hj][param[:Lj] .== 1e-6] .= 1  # catch errors in places with infinite housing per capita
+graph[:hj] = graph[:Hj] ./ graph[:Lj]
+graph[:hj][graph[:Lj] .== 1e-6] .= 1  # catch errors in places with infinite housing per capita
 
 
 # --------------
@@ -72,7 +72,7 @@ graph = apply_geography(graph, geography, alpha_up_i = 10, alpha_down_i = 10)
 plot_graph(graph, aspect_ratio = 3/4, 
             geography = geography, obstacles = true,
             mesh = true, mesh_transparency = 0.2, 
-            node_sizes = param[:Lj], node_shades = param[:Zjn], 
+            node_sizes = graph[:Lj], node_shades = graph[:Zjn], 
             node_color = :seismic,
             edge_min_thickness = 1, edge_max_thickness = 4)
 
@@ -86,8 +86,8 @@ plot_graph(graph, aspect_ratio = 3/4,
 # PLOT RESULTS
 # ============
 
-sizes = 2 .* results[:cj] .* (param[:Lj] .> 1e-6) / maximum(results[:cj])
-shades = results[:cj] .* (param[:Lj] .> 1e-6) / maximum(results[:cj])
+sizes = 2 .* results[:cj] .* (graph[:Lj] .> 1e-6) / maximum(results[:cj])
+shades = results[:cj] .* (graph[:Lj] .> 1e-6) / maximum(results[:cj])
 
 plot_graph(graph, results[:Ijk], aspect_ratio = 3/4,
            geography = geography, obstacles = true,

examples/paper_example01.jl

@@ -10,41 +10,41 @@ using Plots
 
 K = [1, 100]
 param = init_parameters()
-param, g = create_graph(param, 9, 9, type = "map")
+graph = create_graph(param, 9, 9, type = "map")
 
 #  Define fundamentals
-param[:Zjn] *= 0.1; # matrix of productivity
-Ni = find_node(g, 5, 5); # center
-param[:Zjn][Ni, :] .= 1; # more productive node
+graph[:Zjn] *= 0.1; # matrix of productivity
+Ni = find_node(graph, 5, 5); # center
+graph[:Zjn][Ni, :] .= 1; # more productive node
 
 # Plot the mesh with population 
-plot_graph(g, edges = false, mesh = true, node_sizes = param[:Lj])
+plot_graph(graph, edges = false, mesh = true, node_sizes = graph[:Lj])
 
 # Plot the mesh with productivity
-plot_graph(g, edges = false, mesh = true, node_sizes = param[:Zjn])
+plot_graph(graph, edges = false, mesh = true, node_sizes = graph[:Zjn])
 
 # Compute networks
 results = [] 
 for i in 1:length(K) 
     param[:K] = K[i]
-    push!(results, optimal_network(param, g))
+    push!(results, optimal_network(param, graph))
 end
 
 # Plot networks
 
 plots = [] # Initialize an empty array to hold the subplots
 for i in 1:length(K) 
     shades = sizes = results[i][:Cj] / maximum(results[i][:Cj])
-    p = plot_graph(g, results[i][:Ijk], node_sizes = sizes, node_shades = shades, node_sizes_scale = 40)
+    p = plot_graph(graph, results[i][:Ijk], node_sizes = sizes, node_shades = shades, node_sizes_scale = 40)
     title!(p, "(a) Transport Network (I_{jk})")
     push!(plots, p) 
-    p = plot_graph(g, results[1][:Qjkn][:, :, 1], edge_color = :brown, arrows = true, node_sizes = sizes, node_shades = shades, node_sizes_scale = 40)
+    p = plot_graph(graph, results[1][:Qjkn][:, :, 1], edge_color = :brown, arrows = true, node_sizes = sizes, node_shades = shades, node_sizes_scale = 40)
     title!(p, "(b) Shipping (Q_{jk})")
     push!(plots, p) 
-    p = plot_graph(g, results[i][:Ijk], map = results[i][:Pjn] / maximum(results[i][:Pjn]), node_sizes = sizes, node_shades = shades, node_sizes_scale = 40)
+    p = plot_graph(graph, results[i][:Ijk], map = results[i][:Pjn] / maximum(results[i][:Pjn]), node_sizes = sizes, node_shades = shades, node_sizes_scale = 40)
     title!(p, "(c) Prices (P_{j})")
     push!(plots, p) 
-    p = plot_graph(g, results[i][:Ijk], map = results[i][:cj] / maximum(results[i][:cj]), node_sizes = sizes, node_shades = shades, node_sizes_scale = 40)
+    p = plot_graph(graph, results[i][:Ijk], map = results[i][:cj] / maximum(results[i][:cj]), node_sizes = sizes, node_shades = shades, node_sizes_scale = 40)
     title!(p, "(d) Consumption (c_{j})")
     push!(plots, p) 
 end

examples/paper_example02.jl

@@ -11,42 +11,42 @@ height = 9
 nb_cities = 20
 
 param = init_parameters(K = 100, annealing = false)
-param, g = create_graph(param, width, height, type = "triangle") # case with random cities, one good
+graph = create_graph(param, width, height, type = "triangle") # case with random cities, one good
 
 # set fundamentals
 
 Random.seed!(5)
-minpop = minimum(param[:Lj])
-param[:Zjn] = fill(minpop, g[:J], 1)  # matrix of productivity
+minpop = minimum(graph[:Lj])
+graph[:Zjn] = fill(minpop, graph[:J], 1)  # matrix of productivity
 
-Ni = find_node(g, ceil(width/2), ceil(height/2)) # find center
-param[:Zjn][Ni] = 1 # more productive node
-param[:Lj][Ni] = 1  # more productive node
+Ni = find_node(graph, ceil(width/2), ceil(height/2)) # find center
+graph[:Zjn][Ni] = 1 # more productive node
+graph[:Lj][Ni] = 1  # more productive node
 
 for i in 1:nb_cities-1
     newdraw = false
     while newdraw == false
-        j = round(Int, 1 + rand() * (g[:J] - 1))
-        if param[:Lj][j] <= minpop
+        j = round(Int, 1 + rand() * (graph[:J] - 1))
+        if graph[:Lj][j] <= minpop
             newdraw = true
-            param[:Lj][j] = 1
+            graph[:Lj][j] = 1
         end
     end
 end
 
-param[:hj] = param[:Hj] ./ param[:Lj]
-param[:hj][param[:Lj] .== minpop] .= 1
+graph[:hj] = graph[:Hj] ./ graph[:Lj]
+graph[:hj][graph[:Lj] .== minpop] .= 1
 
 # Convex case
 results = Vector{Any}(undef, 3)
-results[1] = optimal_network(param, g)
+results[1] = optimal_network(param, graph)
 
 # Nonconvex - no annealing
 param[:gamma] = 2
-results[2] = optimal_network(param, g)
+results[2] = optimal_network(param, graph)
 
 # Nonconvex - annealing
-results[3] = annealing(param, g, results[2][:Ijk], perturbation_method = "rebranching")
+results[3] = annealing(graph, results[2][:Ijk], perturbation_method = "rebranching")
 
 welfare_increase = (results[3][:welfare] / results[2][:welfare]) ^ (1 / (param[:alpha] * (1 - param[:rho]))) # compute welfare increase in consumption equivalent
 
@@ -55,18 +55,18 @@ welfare_increase = (results[3][:welfare] / results[2][:welfare]) ^ (1 / (param[:
 
 plots = Vector{Any}(undef, 6) # Initialize an empty array to hold the subplots
 
-plots[1] = plot_graph(g, results[1][:Ijk], node_sizes = results[1][:Lj], node_shades = results[1][:Cj], node_sizes_scale = 40)
+plots[1] = plot_graph(graph, results[1][:Ijk], node_sizes = results[1][:Lj], node_shades = results[1][:Cj], node_sizes_scale = 40)
 title!(plots[1], "Convex Network (I_{jk})")
-plots[2] = plot_graph(g, results[1][:Qjkn][:, :, 1], edge_color = :brown, arrows = true, node_sizes = results[1][:Lj], node_shades = results[1][:Cj], node_sizes_scale = 40)
+plots[2] = plot_graph(graph, results[1][:Qjkn][:, :, 1], edge_color = :brown, arrows = true, node_sizes = results[1][:Lj], node_shades = results[1][:Cj], node_sizes_scale = 40)
 title!(plots[2], "Convex Shipping (Q_{jk})")
 
-plots[3] = plot_graph(g, results[2][:Ijk], node_sizes = results[2][:Lj], node_shades = results[2][:Cj], node_sizes_scale = 40)
+plots[3] = plot_graph(graph, results[2][:Ijk], node_sizes = results[2][:Lj], node_shades = results[2][:Cj], node_sizes_scale = 40)
 title!(plots[3], "Nonconvex Network (I_{jk})")
-plots[4] = plot_graph(g, results[2][:Qjkn][:, :, 1], edge_color = :brown, arrows = true, node_sizes = results[2][:Lj], node_shades = results[2][:Cj], node_sizes_scale = 40)
+plots[4] = plot_graph(graph, results[2][:Qjkn][:, :, 1], edge_color = :brown, arrows = true, node_sizes = results[2][:Lj], node_shades = results[2][:Cj], node_sizes_scale = 40)
 title!(plots[4], "Nonconvex Shipping (Q_{jk})")
 
 plots[5] = plots[3]
-plots[6] = plot_graph(g, results[3][:Ijk], node_sizes = results[3][:Lj], node_shades = results[3][:Cj], node_sizes_scale = 40)
+plots[6] = plot_graph(graph, results[3][:Ijk], node_sizes = results[3][:Lj], node_shades = results[3][:Cj], node_sizes_scale = 40)
 title!(plots[6], "Nonconvex Network (I_{jk})")
 annotate!(plots[6], [(0.5, 1.04, text("With Annealing. Welfare increase: $(round((welfare_increase-1)*100, digits = 2))%", :black, :center, 10))])
 

examples/paper_example03.jl

@@ -14,41 +14,41 @@ height = 9
 # Init and Solve network
 
 param = init_parameters(K = 10000, labor_mobility = true, annealing = false)
-param, g = create_graph(param, width, height, type = "triangle")
+graph = create_graph(param, width, height, type = "triangle")
 
 # set fundamentals
 
 param[:N] = 1 + ngoods
-param[:Zjn] = zeros(g[:J], param[:N]) # matrix of productivity
+graph[:Zjn] = zeros(graph[:J], param[:N]) # matrix of productivity
 
-param[:Zjn][:,1] .= 1 # the entire countryside produces the homogenenous good 1
+graph[:Zjn][:,1] .= 1 # the entire countryside produces the homogenenous good 1
 
 if ngoods > 0
-    Ni = find_node(g, 5, 5) # center
-    param[:Zjn][Ni, 2] = 1 # central node
-    param[:Zjn][Ni, 1] = 0 # central node
+    Ni = find_node(graph, 5, 5) # center
+    graph[:Zjn][Ni, 2] = 1 # central node
+    graph[:Zjn][Ni, 1] = 0 # central node
 
     Random.seed!(5)
     for i in 2:ngoods
         newdraw = false
         while newdraw == false
-            j = round(Int, 1 + rand() * (g[:J] - 1))
-            if param[:Zjn][j, 1] > 0
+            j = round(Int, 1 + rand() * (graph[:J] - 1))
+            if graph[:Zjn][j, 1] > 0
                 newdraw = true
-                param[:Zjn][j, i+1] = 1
-                param[:Zjn][j, 1] = 0
+                graph[:Zjn][j, i+1] = 1
+                graph[:Zjn][j, 1] = 0
             end
         end
     end
 end
 
 # Convex case
 results = Array{Any}(undef, 2)
-results[1] = optimal_network(param, g)
+results[1] = optimal_network(param, graph)
 
 # Nonconvex
 param[:gamma] = 2
-results[2] = optimal_network(param, g, I0 = results[1][:Ijk])
+results[2] = optimal_network(graph, I0 = results[1][:Ijk])
 
 # Plot results
 
@@ -59,13 +59,13 @@ for j in 1:2
     # Initialize an empty array to hold the subplots
     plots = Vector{Any}(undef, (1 + param[:N])) 
     # Plot network
-    plots[1] = plot_graph(g, results[j][:Ijk], node_shades = results[j][:Lj], node_sizes = results[j][:Lj], node_sizes_scale = 40)
+    plots[1] = plot_graph(graph, results[j][:Ijk], node_shades = results[j][:Lj], node_sizes = results[j][:Lj], node_sizes_scale = 40)
     title!(plots[1], "(a) Transport Network")
     # Plot goods flows
     for i in 1:param[:N]
-        plots[i+1] = plot_graph(g, results[j][:Qjkn][:, :, i], edge_color = :brown, arrows = true, arrow_style = "thin", 
+        plots[i+1] = plot_graph(graph, results[j][:Qjkn][:, :, i], edge_color = :brown, arrows = true, arrow_style = "thin", 
                                 node_sizes = results[j][:Yjn][:, i], node_sizes_scale = 40,
-                                node_shades = param[:Zjn][:, i])
+                                node_shades = graph[:Zjn][:, i])
         title!(plots[i+1], string('(', Char(96 + i + 1), ')', " Flows Good ", i))
     end
     # Combine plots