Merge pull request #15 from SebKrantz/main

Update devel

NEWS.md

@@ -0,0 +1,13 @@
+# 0.1.6
+## Breaking Changes
+* All spatial parameters, including `Lj`, `Lr`, `Hj`, `hj`, `Zjn`, `omegaj`, `omegar`, and `region` are now stored in the `graph` structure created by `create_graph()`. `create_graph()` therefore only returns the `graph` structure, instead of both the (updated) parameters and the graph. Converesely, `init_parameters()` only contains parameters that are independent of the particular geography defined by `create_graph()`.
+
+## Improvements
+* Minor improvements to Simulated Annealing.
+* Better spline options for plotting frictions surface (geography).
+* More faithful translation of `apply_geography()`. 
+
+
+# 0.1.5
+
+* Removed the MATLAB toolbox and corresponding documentation (PDF files) from the repo to decrease size. A new repo was created for the MATLAB toolbox at [SebKrantz/OptimalTransportNetworkToolbox](https://github.com/SebKrantz/OptimalTransportNetworkToolbox). This repo and especially the `docs` folder continue to be very useful for Julia users, but are no longer part of the Julia library. 

README.md

@@ -42,15 +42,15 @@ import MathOptInterface as MOI
 import MathOptSymbolicAD
 
 param[:model_attr] = Dict(:backend => (MOI.AutomaticDifferentiationBackend(), 
-                                        MathOptSymbolicAD.DefaultBackend())) 
+                                       MathOptSymbolicAD.DefaultBackend())) 
                                 # Or:  MathOptSymbolicAD.ThreadedBackend()
 ```
 
 * It is recommended to use Coin-HSL linear solvers for [Ipopt](https://github.com/jump-dev/Ipopt.jl) to speed up computations. In my opinion the simplest way to use them is to get a (free for academics) license and download the binaries [here](https://licences.stfc.ac.uk/product/coin-hsl), extract them somewhere, and then set the `hsllib` (place here the path to where you extracted `libhsl.dylib`, it may also be called `libcoinhsl.dylib`, in which case you may have to rename it to `libhsl.dylib`) and `linear_solver` options as follows:
 
 ```julia
 param[:optimizer_attr] = Dict(:hsllib => "/usr/local/lib/libhsl.dylib", # Adjust path
-                                :linear_solver => "ma57") # Use ma57, ma86 or ma97
+                              :linear_solver => "ma57") # Use ma57, ma86 or ma97
 ```
 
 The [Ipopt.jl README](https://github.com/jump-dev/Ipopt.jl?tab=readme-ov-file#linear-solvers) suggests to use the larger LibHSL package for which there exists a Julia module and proceed similarly. In addition, users may try an [optimized BLAS](https://github.com/jump-dev/Ipopt.jl?tab=readme-ov-file#blas-and-lapack) and see if it yields significant performance gains (and let me know if it does). 

docs/src/api.md

@@ -16,17 +16,28 @@ Order   = [:function]
 
 # Exported Functions
 
-```@autodocs
-Modules = [OptimalTransportNetworks]
-Private = false
-Order = [:function]
+```@docs
+init_parameters
+create_graph
+apply_geography
+optimal_network
+annealing
+plot_graph
+find_node
+add_node
+remove_node
 ```
 
 # Documented Internal Functions
 
-```@autodocs
-Modules = [OptimalTransportNetworks]
-Public = false
-Private = true
-Order = [:function]
-```
+```@docs
+OptimalTransportNetworks.dict_to_namedtuple
+OptimalTransportNetworks.namedtuple_to_dict
+OptimalTransportNetworks.create_map
+OptimalTransportNetworks.create_square
+OptimalTransportNetworks.create_triangle
+OptimalTransportNetworks.create_custom
+OptimalTransportNetworks.represent_edges
+OptimalTransportNetworks.create_auxdata
+OptimalTransportNetworks.get_model
+```

docs/src/index.md

@@ -37,15 +37,15 @@ import MathOptInterface as MOI
 import MathOptSymbolicAD
 
 param[:model_attr] = Dict(:backend => (MOI.AutomaticDifferentiationBackend(), 
-                                        MathOptSymbolicAD.DefaultBackend())) 
+                                       MathOptSymbolicAD.DefaultBackend())) 
                                 # Or:  MathOptSymbolicAD.ThreadedBackend()
 ```
 
 * It is recommended to use Coin-HSL linear solvers for [Ipopt](https://github.com/jump-dev/Ipopt.jl) to speed up computations. In my opinion the simplest way to use them is to get a (free for academics) license and download the binaries [here](https://licences.stfc.ac.uk/product/coin-hsl), extract them somewhere, and then set the `hsllib` (place here the path to where you extracted `libhsl.dylib`, it may also be called `libcoinhsl.dylib`, in which case you may have to rename it to `libhsl.dylib`) and `linear_solver` options as follows:
 
 ```julia
 param[:optimizer_attr] = Dict(:hsllib => "/usr/local/lib/libhsl.dylib", # Adjust path
-                                :linear_solver => "ma57") # Use ma57, ma86 or ma97
+                              :linear_solver => "ma57") # Use ma57, ma86 or ma97
 ```
 
 The [Ipopt.jl README](https://github.com/jump-dev/Ipopt.jl?tab=readme-ov-file#linear-solvers) suggests to use the larger LibHSL package for which there exists a Julia module and proceed similarly. In addition, users may try an [optimized BLAS](https://github.com/jump-dev/Ipopt.jl?tab=readme-ov-file#blas-and-lapack) and see if it yields significant performance gains (and let me know if it does). 

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

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(param, 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))])