Merge pull request #8 from SebKrantz/main

Update.

Author: SebKrantz

Project.toml

@@ -1,29 +1,34 @@
 name = "OptimalTransportNetworks"
 uuid = "e2b46e68-897f-4e4e-ba36-a93c9789fd96"
 authors = ["Sebastian Krantz <sebastian.krantz@graduateinstitute.ch>"]
-version = "0.1.4"
+version = "0.1.7"
 
 [deps]
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-julia = "1.8.5"
-Dierckx = "0.5.0"
-Ipopt = "1.4.0"
-JuMP = "1.20.0"
-MathOptInterface = "1.0.0"
-Plots = "1.19.0"
-Random = "1.0.0"
-SparseArrays = "1.0.0"
-Statistics = "1.0.0"
-LinearAlgebra = "1.0.0"
+julia = "1.8.5, 1.9"
+Ipopt = "1.4.0, 2"
+JuMP = "1.20.0, 2"
+LinearAlgebra = "1.0.0, 2"
+NearestNeighbors = "0.4.10, 1"
+Plots = "1.19.0, 2"
+Dierckx = "0.5.0, 1"
+Random = "1.0.0, 2"
+SparseArrays = "1.0.0, 2"
+StaticArrays = "1.1.0, 2"
+Statistics = "1.0.0, 2"
+
+
+
 
 

README.md

@@ -1,13 +1,19 @@
 # OptimalTransportNetworks.jl
+
+[![](https://img.shields.io/badge/docs-stable-blue.svg)](https://sebkrantz.github.io/OptimalTransportNetworks.jl/stable)
+[![](https://img.shields.io/badge/docs-dev-blue.svg)](https://sebkrantz.github.io/OptimalTransportNetworks.jl/dev)
+[![Julia version](https://juliahub.com/docs/General/OptimalTransportNetworks/stable/version.svg)](https://juliahub.com/ui/Packages/General/OptimalTransportNetworks)
+
 **Optimal Transport Networks in Spatial Equilibrium - in Julia**
 
+
 Modern Julia ([JuMP](https://github.com/jump-dev/JuMP.jl)) translation of the [MATLAB OptimalTransportNetworkToolbox (v1.0.4b)](https://github.com/SebKrantz/OptimalTransportNetworkToolbox) implementing the quantitative spatial economic model of:
 
 Fajgelbaum, P. D., & Schaal, E. (2020). Optimal transport networks in spatial equilibrium. *Econometrica, 88*(4), 1411-1452.
 
 The model/software uses duality principles to optimize over the space of networks, nesting an optimal flows problem and a neoclasical general-equilibrium trade model into a global network design problem to derive the optimal (welfare maximizing) transport network (extension) from any primitive set of economic fundamantals [population per location, productivity per location for each of *N* traded goods, endowment of a non-traded good, and (optionally) a pre-existing transport network]. 
 
-For more information about the model see [this folder](https://github.com/SebKrantz/OptimalTransportNetworks.jl/tree/main/misc/paper_materials) and the [MATLAB User Guide](https://raw.githubusercontent.com/SebKrantz/OptimalTransportNetworkToolbox/main/docs/User%20Guide.pdf). 
+For more information about the model see [this folder](https://github.com/SebKrantz/OptimalTransportNetworkToolbox/tree/main/docs/paper_materials) and the [MATLAB User Guide](https://raw.githubusercontent.com/SebKrantz/OptimalTransportNetworkToolbox/main/docs/User%20Guide.pdf). 
 
 The model is the first of its kind and a pathbreaking contribution towards the welfare maximizing planning of transport infrastructure. Its creation has been funded by the European Union through an [ERC Research Grant](https://cordis.europa.eu/project/id/804095). The author of this Julia library has no personal connections to the authors, but has used their Matlab library for research purposes and belives that it deserves an accessible open-source implementation. Community efforts to further improve the code are welcome. In particular, there is a probabilistic extenstion to solving the model using MCMC methods which may be more suitable for large networks, implemented in:
 

__init__.jl

@@ -8,7 +8,7 @@ Fajgelbaum, P. D., & Schaal, E. (2020). Optimal transport networks in spatial eq
 
 The model/software uses duality principles to optimize over the space of networks, nesting an optimal flows problem and a neoclasical general-equilibrium trade model into a global network design problem to derive the optimal (welfare maximizing) transport network (extension) from any primitive set of economic fundamantals [population per location, productivity per location for each of *N* traded goods, endowment of a non-traded good, and (optionally) a pre-existing transport network]. 
 
-For more information about the model see [this folder](https://github.com/SebKrantz/OptimalTransportNetworks.jl/tree/main/misc/paper_materials) and the [MATLAB User Guide](https://raw.githubusercontent.com/SebKrantz/OptimalTransportNetworkToolbox/main/docs/User%20Guide.pdf). 
+For more information about the model see [this folder](https://github.com/SebKrantz/OptimalTransportNetworkToolbox/tree/main/docs/paper_materials) and the [MATLAB User Guide](https://raw.githubusercontent.com/SebKrantz/OptimalTransportNetworkToolbox/main/docs/User%20Guide.pdf). 
 
 The model is the first of its kind and a pathbreaking contribution towards the welfare maximizing planning of transport infrastructure. Its creation has been funded by the European Union through an [ERC Research Grant](https://cordis.europa.eu/project/id/804095). The author of this Julia library has no personal connections to the authors, but has used their Matlab library for research purposes and belives that it deserves an accessible open-source implementation. Community efforts to further improve the code are welcome. In particular, there is a probabilistic extenstion to solving the model using MCMC methods which may be more suitable for large networks, implemented in:
 

docs/src/index.md

@@ -16,8 +16,8 @@ param, graph = create_graph(param, w, h, type = "map")
 # ----------------
 # Draw populations
 
-param[:Zjn] = ones(param[:J], 1) .* 1e-3  # matrix of productivity (not 0 to avoid numerical glitches)
-param[:Lj] = ones(param[:J]) .* 1e-3  # matrix of population
+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
 
 Ni = find_node(graph, ceil(w/2), ceil(h/2))  # center
 param[:Zjn][Ni] = 1  # more productive node
@@ -29,15 +29,15 @@ for i in 1:ncities-1
     newdraw = false
     while newdraw == false
         j = round(Int, 1 + rand() * (param[:J] - 1))
-        if param[:Lj][j] <= 1e-3
+        if param[:Lj][j] <= 1e-6
             newdraw = true
             param[:Lj][j] = 1
         end
     end
 end
 
 param[:hj] = param[:Hj] ./ param[:Lj]
-param[:hj][param[:Lj] .== 1e-3] .= 1  # catch errors in places with infinite housing per capita
+param[:hj][param[:Lj] .== 1e-6] .= 1  # catch errors in places with infinite housing per capita
 
 
 # --------------
@@ -68,6 +68,14 @@ geography = (z = z, obstacles = obstacles)
 # to changes in elevation in building costs (symmetric up/down)
 graph = apply_geography(graph, geography, alpha_up_i = 10, alpha_down_i = 10)
 
+# Plot endowments          
+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_color = :seismic,
+            edge_min_thickness = 1, edge_max_thickness = 4)
+
 # =======================
 # COMPUTE OPTIMAL NETWORK
 # =======================
@@ -78,24 +86,12 @@ graph = apply_geography(graph, geography, alpha_up_i = 10, alpha_down_i = 10)
 # PLOT RESULTS
 # ============
 
-sizes = 2 .* results[:cj] .* (param[:Lj] .> 1e-3) / maximum(results[:cj])
-shades = results[:cj] .* (param[:Lj] .> 1e-3) / maximum(results[:cj])
+sizes = 2 .* results[:cj] .* (param[:Lj] .> 1e-6) / maximum(results[:cj])
+shades = results[:cj] .* (param[:Lj] .> 1e-6) / maximum(results[:cj])
 
-plot_graph(graph, results[:Ijk], 
+plot_graph(graph, results[:Ijk], aspect_ratio = 3/4,
            geography = geography, obstacles = true,
            mesh = true, mesh_transparency = 0.2, 
            node_sizes = sizes, node_shades = shades, 
            edge_min_thickness = 1, edge_max_thickness = 4)
 
-
-# sizes = param[:Lj] .+ 1
-# shades = param[:Zjn]
-           
-# plot_graph(graph, # results[:Ijk], edges = false,
-#             geography = geography, obstacles = true,
-#             mesh = true, mesh_transparency = 0.2, 
-#             node_sizes = sizes, node_shades = shades, 
-#             node_color = :seismic,
-#             edge_min_thickness = 1, edge_max_thickness = 4)
-           
-                     

examples/example04.jl

@@ -1,49 +1,44 @@
 using OptimalTransportNetworks
 using Random
-using LinearAlgebra
 using Plots
 
 # Initialize parameters
 param = init_parameters(K = 100, labor_mobility = false)
 width, height = 13, 13
 
 # Create graph
-param, g = create_graph(param, width, height, type = "map")
+param, g0 = create_graph(param, width, height, type = "map")
 
 # Set fundamentals
-Random.seed!(5)
-minpop = minimum(param[:Lj])
-param[:Zjn] = fill(minpop, g[:J], 1) # matrix of productivity
-param[:Lj] = fill(1e-6, g[:J]) # matrix of population
+param[:Zjn] = ones(param[:J], 1) .* 1e-7  # matrix of productivity (not 0 to avoid numerical glitches)
+param[:Lj] = ones(param[:J]) .* 1e-7  # matrix of population
 
-Ni = find_node(g, 5, 5) # center
+Ni = find_node(g0, ceil(width/2), ceil(height/2)) # center
 param[:Zjn][Ni] = 1 # more productive node
 param[:Lj][Ni] = 1 # more productive node
 
+Random.seed!(10) # use 5, 8, 10 or 11
 nb_cities = 20 # draw a number of random cities in space
 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() * (param[:J] - 1))
+        if param[:Lj][j] <= 1e-7
             newdraw = true
             param[:Lj][j] = 1
         end
     end
 end
 
 param[:hj] = param[:Hj] ./ param[:Lj]
-param[:hj][param[:Lj] .== 1e-6] .= 1
+param[:hj][param[:Lj] .== 1e-7] .= 1  # catch errors in places with infinite housing per capita
 
 # Draw geography
-z = zeros(g[:J]) # altitude of each node
-
+z = zeros(param[:J]) # altitude of each node
 geographies = Dict()
 geographies[:blank] = (z = z, obstacles = nothing)
-g = apply_geography(g, geographies[:blank])
 
-param0 = param # store initial params
-g0 = g # store initial graph
+g = apply_geography(g0, geographies[:blank])
 
 # Blank geography
 results = Dict(:blank => optimal_network(param, g))
@@ -53,13 +48,12 @@ mountain_size = 0.75
 mountain_height = 1
 mount_x = 10
 mount_y = 10
-geographies[:mountain] = (z = mountain_height * exp.(-((g[:x] .- mount_x).^2 .+ (g[:y] .- mount_y).^2) / (2 * mountain_size^2)),
+geographies[:mountain] = (z = mountain_height * exp.(-((g0[:x] .- mount_x).^2 .+ (g0[:y] .- mount_y).^2) / (2 * mountain_size^2)),
                           obstacles = nothing)
-g = apply_geography(g, geographies[:mountain], alpha_up_i = 10, alpha_down_i = 10)
+g = apply_geography(g0, geographies[:mountain], alpha_up_i = 10, alpha_down_i = 10)
 results[:mountain] = optimal_network(param, g)
 
 # Adding river and access by land
-g = g0
 geographies[:river] = (z = geographies[:mountain].z, 
                        obstacles = [6 + (1-1)*width 6 + (2-1)*width;
                           6 + (2-1)*width 6 + (3-1)*width;
@@ -69,13 +63,11 @@ geographies[:river] = (z = geographies[:mountain].z,
                           11 + (5-1)*width 12 + (5-1)*width;
                           12 + (5-1)*width 13 + (5-1)*width])
 
-g = apply_geography(g, geographies[:river], across_abstacle_delta_i = Inf, 
-                    alpha_up_i = 10, alpha_down_i = 10)
+g = apply_geography(g0, geographies[:river], alpha_up_i = 10, alpha_down_i = 10)
 results[:river] = optimal_network(param, g)
 
 # Reinit and put another river and bridges
-g = g0
-geographies[:bridges] = (z = mountain_height * exp.(-((g[:x] .- mount_x).^2 .+ (g[:y] .- mount_y).^2) / (2 * mountain_size^2)),
+geographies[:bridges] = (z = mountain_height * exp.(-((g0[:x] .- mount_x).^2 .+ (g0[:y] .- mount_y).^2) / (2 * mountain_size^2)),
                          obstacles = [6 + (1-1)*width 6 + (2-1)*width;
                              6 + (2-1)*width 6 + (3-1)*width;
                              6 + (3-1)*width 7 + (4-1)*width;
@@ -86,13 +78,12 @@ geographies[:bridges] = (z = mountain_height * exp.(-((g[:x] .- mount_x).^2 .+ (
                              11 + (5-1)*width 12 + (5-1)*width;
                              12 + (5-1)*width 13 + (5-1)*width])
 
-g = apply_geography(g, geographies[:bridges], alpha_up_i = 10, alpha_down_i = 10, 
-                    across_abstacle_delta_i = 2, along_obstacle_delta_i = Inf)
+g = apply_geography(g0, geographies[:bridges], alpha_up_i = 10, alpha_down_i = 10, 
+                    across_obstacle_delta_i = 2)
 results[:bridges] = optimal_network(param, g)
 
 # Allowing for water transport
-g = g0
-geographies[:water_transport] = (z = mountain_height * exp.(-((g[:x] .- mount_x).^2 .+ (g[:y] .- mount_y).^2) / (2 * mountain_size^2)),
+geographies[:water_transport] = (z = mountain_height * exp.(-((g0[:x] .- mount_x).^2 .+ (g0[:y] .- mount_y).^2) / (2 * mountain_size^2)),
                                  obstacles = [6 + (1-1)*width 6 + (2-1)*width;
                                               6 + (2-1)*width 6 + (3-1)*width;
                                               6 + (3-1)*width 7 + (4-1)*width;
@@ -103,8 +94,10 @@ geographies[:water_transport] = (z = mountain_height * exp.(-((g[:x] .- mount_x)
                                               11 + (5-1)*width 12 + (5-1)*width;
                                               12 + (5-1)*width 13 + (5-1)*width])
 
-g = apply_geography(g, geographies[:water_transport], alpha_up_i = 10, alpha_down_i = 10, 
-                    across_abstacle_delta_i = 2, along_obstacle_delta_i = 0.5)
+g = apply_geography(g0, geographies[:water_transport], alpha_up_i = 10, alpha_down_i = 10, 
+                    across_abstacle_delta_i = 2, 
+                    along_obstacle_delta_i = 0.5,
+                    along_abstacle_delta_tau = 1)
 results[:water_transport] = optimal_network(param, g)
 
 # Increasing returns to transport
@@ -120,12 +113,12 @@ plots = Vector{Any}(undef, length(simulation))
 i = 0
 for s in simulation
     i += 1
-    plots[i] = plot_graph(g, results[Symbol(s)][:Ijk], 
+    plots[i] = plot_graph(g0, results[Symbol(s)][:Ijk], 
                           geography = geographies[Symbol(s)], obstacles = obstacles[i] == "on",
                           mesh = true, mesh_transparency = 0.2,
-                          node_sizes = results[Symbol(s)][:cj] .* (param[:Lj] .> minpop), 
-                          node_shades = param[:Zjn], 
-                          edge_min_thickness = 1.5)
+                          node_sizes = results[Symbol(s)][:cj] .* (param[:Lj] .> 1e-7), 
+                          node_shades = param[:Zjn], node_color = :seismic,
+                          edge_min_thickness = 1, edge_max_thickness = 4)
     title!(plots[i], "Geography $(s)")
 end