Make graph largely independent of the parameters structure.

Author: SebKrantz

src/main/create_graph.jl

@@ -2,38 +2,35 @@
 # using LinearAlgebra
 
 """
-    create_graph(param, w = 11, h = 11; type = "map", kwargs...) -> Dict, Dict
+    create_graph(param, w = 11, h = 11; type = "map", N = 1, kwargs...) -> Dict
 
 Initialize the underlying graph, population and productivity parameters.
 
 # Arguments
-- `param::Dict`: Dict that contains the model parameters
-- `w::Int64=11`: Number of nodes along the width of the underlying graph (integer)  
-- `h::Int64=11`: Number of nodes along the height of the underlying graph (integer, odd if triangle)
+- `param::Dict`: Dict that contains the model parameters (only needed for checks)
+- `w::Int64=11`: Number of nodes along the width of the underlying graph if type != "custom" (integer)  
+- `h::Int64=11`: Number of nodes along the height of the underlying graph if type != "custom" (integer, odd if triangle)
 
 # Keyword Arguments
 - `type::String="map"`: Either "map", "square", "triangle", or "custom" 
-- `omega::Vector{Float64}`: Vector of Pareto weights for each node or region in partial mobility case (default ones(J or nregions))
-- `Zjn::Matrix{Float64}`: J x N matrix of producties per node (j = 1:J) and good (n = 1:N) (default ones(J, N))
 - `adjacency::BitMatrix`: J x J Adjacency matrix (only used for custom network)
 - `x::Vector{Float64}`: x coordinate (longitude) of each node (only used for custom network)
 - `y::Vector{Float64}`: y coordinate (latitude) of each node (only used for custom network)
+- `omega::Vector{Float64}`: Vector of Pareto weights for each node or region in partial mobility case (default ones(J or nregions))
+- `m::Vector{Float64}=ones(N)`: Vector of weights Nx1 in the cross congestion cost function
+- `Zjn::Matrix{Float64}`: J x N matrix of producties per node (j = 1:J) and good (n = 1:N) (default ones(J, N))
 - `Lj::Vector{Float64}`: Vector of populations in each node (j = 1:J) (only for fixed labour case)
 - `Hj::Vector{Float64}`: Vector of immobile good in each node (j = 1:J) (e.g. housing, default ones(J))
 - `Lr::Vector{Float64}`: Vector of populations in each region (r = 1:nregions) (only for partial mobility)
 - `region::Vector{Int64}`: Vector indicating region of each location (only for partial mobility)
 
-# Notes
-- `create_graph()` will overwrite any parameters `Zjn`, `Lj`, `Hj`, `omega`, `Lr`, already set in the `param` Dict with the default values. So these values should be set inside `create_graph()`, or after `create_graph()` has been called.
-
 # Examples
 ```julia
-param, graph = create_graph(init_parameters())
+graph = create_graph()
 ```
 """
 function create_graph(param, w = 11, h = 11; type = "map", kwargs...)
 
-    param = namedtuple_to_dict(param)
     options = retrieve_options_create_graph(param, w, h, type; kwargs...)
 
     if type == "map"
@@ -46,26 +43,26 @@ function create_graph(param, w = 11, h = 11; type = "map", kwargs...)
         graph = create_custom(options[:adjacency], options[:x], options[:y])
     end
 
-    param[:J] = graph[:J]
-    param[:Zjn] = haskey(options, :Zjn) ? options[:Zjn] : ones(graph[:J], param[:N])
-    param[:Hj] = haskey(options, :Hj) ? options[:Hj] : ones(graph[:J])
+    graph[:N] = param[:N] # Ensure graph has all information
+    graph[:Zjn] = haskey(options, :Zjn) ? options[:Zjn] : ones(graph[:J], graph[:N])
+    graph[:Hj] = haskey(options, :Hj) ? options[:Hj] : ones(graph[:J])
 
     if param[:mobility] == false
-        param[:Lj] = haskey(options, :Lj) ? options[:Lj] : ones(graph[:J]) / graph[:J]
-        param[:hj] = param[:Hj] ./ param[:Lj]
-        param[:hj][param[:Lj] .== 0] .= 1
-        param[:omegaj] = options[:omega]
+        graph[:Lj] = haskey(options, :Lj) ? options[:Lj] : ones(graph[:J]) / graph[:J]
+        graph[:hj] = graph[:Hj] ./ graph[:Lj]
+        graph[:hj][graph[:Lj] .== 0] .= 1
+        graph[:omegaj] = options[:omega]
     elseif param[:mobility] == 0.5
         graph[:region] = options[:region]
-        param[:omegar] = options[:omega]
-        param[:Lr] = options[:Lr]
-        param[:nregions] = options[:nregions]
+        graph[:omegar] = options[:omega]
+        graph[:Lr] = options[:Lr]
+        graph[:nregions] = options[:nregions]
     end
 
     # Running checks on population / productivity / housing parameters
-    check_graph_param(param)
+    check_graph_param(graph, param)
 
-    return param, graph
+    return graph
 end
 
 function isadjacency(M)

src/main/init_parameters.jl

@@ -1,35 +1,30 @@
 
-# ==============================================================
-#  OPTIMAL TRANSPORT NETWORKS IN SPATIAL EQUILIBRIUM
-#  by P. Fajgelbaum, E. Schaal, D. Henricot, C. Mantovani 2017-19
-# ================================================= version 1.0.4
 
 """
     init_parameters(; kwargs...) -> Dict
 
-Returns a `param` dict with the model parameters.
+Returns a `param` dict with the model parameters. These are independent of the graph structure and dimensions.
 
 # Keyword Arguments
 - `alpha::Float64=0.5`: Cobb-Douglas coefficient on final good c^alpha * h^(1-alpha)
 - `beta::Float64=1`: Parameter governing congestion in transport cost
 - `gamma::Float64=1`: Elasticity of transport cost relative to infrastructure
-- `K::Float64=1`: Amount of concrete/asphalt
 - `sigma::Float64=5`: Elasticity of substitution across goods (CES)
 - `rho::Float64=2`: Curvature in utility (c^alpha * h^(1-alpha))^(1-rho)/(1-rho)
 - `a::Float64=0.8`: Curvature of the production function L^alpha
-- `m::Vector{Float64}=ones(N)`: Vector of weights Nx1 in the cross congestion cost function
-- `N::Int64=1`: Number of goods
 - `nu::Float64=1`: Elasticity of substitution b/w goods in transport costs if cross-good congestion
+- `N::Int64=1`: Number of goods traded in the economy (used for checks in `create_graph()`)
+- `K::Float64=1`: Amount of concrete/asphalt
 - `labor_mobility::Any=false`: Switch for labor mobility (true/false or 'partial')
 - `cross_good_congestion::Bool=false`: Switch for cross-good congestion
 - `annealing::Bool=true`: Switch for the use of annealing at the end of iterations (only if gamma > beta)
 - `verbose::Bool=true`: Switch to turn on/off text output (from Ipopt or other optimizers)
 - `duality::Bool=false`: Switch to turn on/off duality whenever available
 - `warm_start::Bool=true`: Use the previous solution as a warm start for the next iteration
-- `kappa_min::Float64=1e-5`: Minimum value for road capacities κ
+- `kappa_min::Float64=1e-5`: Minimum value for road capacities K
 - `min_iter::Int64=20`: Minimum number of iterations
 - `max_iter::Int64=200`: Maximum number of iterations
-- `tol::Float64=1e-5`: Tolerance for convergence of road capacities κ
+- `tol::Float64=1e-5`: Tolerance for convergence of road capacities K
 
 # Optional Parameters
 - `optimizer = Ipopt.Optimizer`: Optimizer to be used
@@ -43,7 +38,7 @@ Returns a `param` dict with the model parameters.
 param = init_parameters(K = 10, labor_mobility = true)
 ```
 """
-function init_parameters(; alpha = 0.5, beta = 1, gamma = 1, K = 1, sigma = 5, rho = 2, a = 0.8, N = 1, m = ones(N), nu = 1, 
+function init_parameters(; alpha = 0.5, beta = 1, gamma = 1, sigma = 5, rho = 2, a = 0.8, nu = 1, K = 1,
                          labor_mobility = false, cross_good_congestion=false, annealing=true, 
                          verbose = true, duality = false, warm_start = true, 
                          kappa_min = 1e-5, min_iter = 20, max_iter = 200, tol = 1e-5, kwargs...)
@@ -65,8 +60,6 @@ function init_parameters(; alpha = 0.5, beta = 1, gamma = 1, K = 1, sigma = 5, r
     param[:sigma] = sigma
     param[:rho] = rho
     param[:a] = a
-    param[:m] = m
-    param[:N] = N
     param[:nu] = nu
     if labor_mobility === "partial" || labor_mobility === 0.5
         param[:mobility] = 0.5
@@ -98,57 +91,61 @@ function init_parameters(; alpha = 0.5, beta = 1, gamma = 1, K = 1, sigma = 5, r
     param[:F] = (L, a) -> L^a
     param[:Fprime] = (L, a) -> a * L^(a - 1)
 
-    # CHECK CONSISTENCY WITH ENDOWMENTS/PRODUCTIVITY (if applicable)
-    check_graph_param(param)
+    # # CHECK CONSISTENCY WITH ENDOWMENTS/PRODUCTIVITY (if applicable)
+    # check_graph_param(param)
 
     return param
 end
 
-function check_graph_param(param)
-
-    ## These are graph parameters that should not be in the parameters dict
-
-    # if haskey(param, :x) && length(param[:x]) != param[:J]
-    #     @warn "x does not have the right length J = $(param[:J])."
-    # end
-
-    # if haskey(param, :y) && length(param[:y]) != param[:J]
-    #     @warn "y does not have the right length J = $(param[:J])."
-    # end
-
-    # if haskey(param, :adjacency) && size(param[:adjacency]) != (param[:J], param[:J])
-    #     @warn "adjacency matrix does not have the right dimensions J x J."
-    # end
-
-    # if haskey(param, :region) && length(param[:region]) != param[:J]
-    #     @warn "region does not have the right length J = $(param[:J])."
-    # end
+# Check if the parameters are consistent with the graph structure
+function check_graph_param(graph, param)
 
-    if haskey(param, :omegaj) && length(param[:omegaj]) != param[:J]
-        @warn "omegaj does not have the right length J = $(param[:J])."
+    if haskey(param, :omegaj) && !haskey(graph, :omegaj)
+        graph[:omegaj] = param[:omegaj]
+    end
+    if haskey(graph, :omegaj) && length(graph[:omegaj]) != graph[:J]
+        @warn "omegaj does not have the right length J = $(graph[:J])."
     end
 
-    if haskey(param, :omegar) && length(param[:omegar]) != param[:nregions]
-        @warn "omegar does not have the right length nregions = $(param[:nregions])."
+    if haskey(param, :omegar) && !haskey(graph, :omegar)
+        graph[:omegar] = param[:omegar]
+    end
+    if haskey(graph, :omegar) && length(graph[:omegar]) != graph[:nregions]
+        @warn "omegar does not have the right length nregions = $(graph[:nregions])."
     end
 
-    if haskey(param, :Lr) && length(param[:Lr]) != param[:nregions]
-        @warn "Lr does not have the right length nregions = $(param[:nregions])."
+    if haskey(param, :Lr) && !haskey(graph, :Lr)
+        graph[:Lr] = param[:Lr]
+    end
+    if haskey(graph, :Lr) && length(graph[:Lr]) != graph[:nregions]
+        @warn "Lr does not have the right length nregions = $(graph[:nregions])."
     end
 
-    if haskey(param, :Lj) && param[:mobility] == 0 && length(param[:Lj]) != param[:J]
-        @warn "Lj does not have the right length J = $(param[:J])."
+    if haskey(param, :Lj) && !haskey(graph, :Lj) && param[:mobility] == 0
+        graph[:Lj] = param[:Lj]
+    end
+    if haskey(graph, :Lj) && param[:mobility] == 0 && length(graph[:Lj]) != graph[:J]
+        @warn "Lj does not have the right length J = $(graph[:J])."
     end
 
-    if haskey(param, :Zjn) && size(param[:Zjn]) != (param[:J], param[:N])
-        @warn "Zjn does not have the right size J ($(param[:J])) x N ($(param[:N]))."
+    if haskey(param, :Zjn) && !haskey(graph, :Zjn)
+        graph[:Zjn] = param[:Zjn]
+    end
+    if haskey(graph, :Zjn) && size(graph[:Zjn]) != (graph[:J], graph[:N])
+        @warn "Zjn does not have the right size J ($(graph[:J])) x N ($(graph[:N]))."
     end
 
-    if haskey(param, :Hj) && length(param[:Hj]) != param[:J]
-        @warn "Hj does not have the right length J = $(param[:J])."
+    if haskey(param, :Hj) && !haskey(graph, :Hj)
+        graph[:Hj] = param[:Hj]
+    end
+    if haskey(graph, :Hj) && length(graph[:Hj]) != graph[:J]
+        @warn "Hj does not have the right length J = $(graph[:J])."
     end
 
-    if haskey(param, :hj) && length(param[:hj]) != param[:J]
-        @warn "hj does not have the right length J = $(param[:J])."
+    if haskey(param, :hj) && !haskey(graph, :hj)
+        graph[:hj] = param[:hj]
+    end
+    if haskey(graph, :hj) && length(graph[:hj]) != graph[:J]
+        @warn "hj does not have the right length J = $(graph[:J])."
     end
 end