Merge pull request #24 from OptimalTransportNetworks/development

Development

Author: SebKrantz

Project.toml

@@ -1,7 +1,7 @@
 name = "OptimalTransportNetworks"
 uuid = "e2b46e68-897f-4e4e-ba36-a93c9789fd96"
 authors = ["Sebastian Krantz <sebastian.krantz@graduateinstitute.ch>"]
-version = "0.1.6"
+version = "0.1.7"
 
 [deps]
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"

src/OptimalTransportNetworks.jl

@@ -75,6 +75,7 @@ module OptimalTransportNetworks
 
 using LinearAlgebra, JuMP, Plots
 using Statistics: mean
+using SparseArrays: sparse, findnz
 # using MathOptInterface: Parameter
 using Dierckx: Spline2D, evaluate
 using NearestNeighbors: KDTree, knn

src/main/annealing.jl

@@ -106,7 +106,7 @@ function annealing(param, graph, I0; kwargs...)
         model, recover_allocation = get_model(auxdata)
     end
 
-    if !options.verbose
+    if !options.verbose && model !== nothing
         set_silent(model)
     end
 
@@ -139,8 +139,12 @@ function annealing(param, graph, I0; kwargs...)
 
     # rng(0) # set the seed of random number generator for replicability
     acceptance_str = ["rejected", "accepted"]
-    all_vars = all_variables_except_kappa_ex(model)
-    start_values = start_value.(all_vars)
+    if model === nothing
+        start_values = []
+    else
+        all_vars = all_variables_except_kappa_ex(model)
+        start_values = start_value.(all_vars)
+    end
     while T > T_min
         accepted = false
 
@@ -161,26 +165,40 @@ function annealing(param, graph, I0; kwargs...)
         end
 
         k = 0
-        set_start_value.(all_vars, start_values) # reset start values
+        # reset start values
+        if model === nothing
+            start_values = []
+        else
+            set_start_value.(all_vars, start_values) 
+        end
         while k <= num_deepening - 1
 
             auxdata = create_auxdata(param, graph, edges, I1)
-            set_parameter_value.(model.obj_dict[:kappa_ex], auxdata.kappa_ex)
-
-            optimize!(model)
-
-            results = recover_allocation(model, auxdata)
-            score = results[:welfare]
-
-            if (!is_solved_and_feasible(model, allow_almost = true) || isnan(score)) && param.verbose # optimization failed
-                println("optimization failed! k=$(k), return flag=$(termination_status(model))")
-                k = num_deepening - 1
-                score = -Inf
-            end
-
-            if param.warm_start
-                vars_solution = value.(all_vars)
-                set_start_value.(all_vars, vars_solution)
+            if model === nothing
+                results, status, start_values = recover_allocation(start_values, auxdata, options.verbose)
+                score = results[:welfare]
+                if (!any(status .== [0, 1]) || isnan(score)) && param.verbose # optimization failed
+                    println("optimization failed! k=$(k), return flag=$(status)")
+                    k = num_deepening - 1
+                    score = -Inf
+                end
+                if !param.warm_start
+                    start_values = []
+                end
+            else
+                set_parameter_value.(model.obj_dict[:kappa_ex], auxdata.kappa_ex)
+                optimize!(model)
+                results = recover_allocation(model, auxdata)
+                score = results[:welfare]
+                if (!is_solved_and_feasible(model, allow_almost = true) || isnan(score)) && param.verbose # optimization failed
+                    println("optimization failed! k=$(k), return flag=$(termination_status(model))")
+                    k = num_deepening - 1
+                    score = -Inf
+                end
+                if param.warm_start
+                    vars_solution = value.(all_vars)
+                    set_start_value.(all_vars, vars_solution)
+                end
             end
 
             if score > best_score
@@ -242,18 +260,31 @@ function annealing(param, graph, I0; kwargs...)
     # Last deepening before returning found optimum
     has_converged = false
     I0 = best_I
-    set_start_value.(all_vars, start_values) # reset start values
+    # reset start values
+    if model === nothing
+        start_values = []
+    else
+        set_start_value.(all_vars, start_values) 
+    end
     while !has_converged && counter < 100
         # Update auxdata
         auxdata = create_auxdata(param, graph, edges, I0)
-        set_parameter_value.(model.obj_dict[:kappa_ex], auxdata.kappa_ex)
 
         # Solve allocation
-        optimize!(model)
-        if !is_solved_and_feasible(model, allow_almost = true)
-            println("Solver returned with error code $(termination_status(model)).")
+        if model === nothing
+            results, status, start_values = recover_allocation(start_values, auxdata, options.verbose)
+            if !any(status .== [0, 1]) 
+                println("Solver returned with error code $(status).")
+            end
+        else
+            set_parameter_value.(model.obj_dict[:kappa_ex], auxdata.kappa_ex)
+            optimize!(model)
+            if !is_solved_and_feasible(model, allow_almost = true)
+                println("Solver returned with error code $(termination_status(model)).")
+            end
+            results = recover_allocation(model, auxdata)
         end
-        results = recover_allocation(model, auxdata)
+
         score = results[:welfare]
 
         # # Fajgelbaum & Schaal always use initial values here. 

src/main/helper.jl

@@ -39,7 +39,8 @@ function namedtuple_to_dict(namedtuple)
     if namedtuple isa Dict
         return namedtuple
     else
-        return Dict(pairs(namedtuple))
+        # return Dict(Symbol(k) => v for (k, v) in pairs(namedtuple))
+        return Dict(pairs(namedtuple)) # Can be of fixed type
     end
 end
 
@@ -243,7 +244,10 @@ function get_model(auxdata)
             recover_allocation = recover_allocation_partial_mobility_cgc
         end
     elseif param.mobility == 0 && param.cong
-        if all(sum(graph.Zjn .> 0, dims = 2) .<= 1) # Armington case
+        if param.beta <= 1 && param.duality # && param.a < 1
+            model = nothing # model_fixed_duality_cgc(optimizer, auxdata)
+            recover_allocation = solve_allocation_by_duality_cgc #recover_allocation_fixed_duality_cgc
+        elseif all(sum(graph.Zjn .> 0, dims = 2) .<= 1) # Armington case
             model = model_fixed_cgc_armington(optimizer, auxdata)
             recover_allocation = recover_allocation_fixed_cgc_armington
         else
@@ -267,12 +271,12 @@ function get_model(auxdata)
             recover_allocation = recover_allocation_partial_mobility    
         end
     elseif param.mobility == 0 && !param.cong
-        if all(sum(graph.Zjn .> 0, dims = 2) .<= 1) # Armington case
+        if param.beta <= 1 && param.duality # && param.a < 1 
+            model = nothing # model_fixed_duality(optimizer, auxdata)
+            recover_allocation = solve_allocation_by_duality # recover_allocation_fixed_duality
+        elseif all(sum(graph.Zjn .> 0, dims = 2) .<= 1) # Armington case
             model = model_fixed_armington(optimizer, auxdata)
             recover_allocation = recover_allocation_fixed_armington
-        elseif param.beta <= 1 && param.a < 1 && param.duality
-            model = model_fixed_duality(optimizer, auxdata)
-            recover_allocation = recover_allocation_fixed_duality
         else
             model = model_fixed(optimizer, auxdata)
             recover_allocation = recover_allocation_fixed
@@ -283,26 +287,27 @@ function get_model(auxdata)
 
     # --------------
     # CUSTOMIZATIONS
-
-    if haskey(param, :optimizer_attr)
-        for (key, value) in param.optimizer_attr
-            set_optimizer_attribute(model, String(key), value)
+    if model !== nothing
+        if haskey(param, :optimizer_attr)
+            for (key, value) in param.optimizer_attr
+                set_optimizer_attribute(model, String(key), value)
+            end
         end
-    end
 
-    if haskey(param, :model_attr) 
-        for value in values(param.model_attr)
-            if !(value isa Tuple)  
-                error("model_attr must be a dict of tuples.")
+        if haskey(param, :model_attr)
+            for value in values(param.model_attr)
+                if !(value isa Tuple)  
+                    error("model_attr must be a dict of tuples.")
+                end
+                set_optimizer_attribute(model, value[1], value[2])
             end
-            set_optimizer_attribute(model, value[1], value[2])
         end
+        # E.g.:
+        # set_attribute(model,
+        #    MOI.AutomaticDifferentiationBackend(),
+        #    MathOptSymbolicAD.DefaultBackend(),
+        # )
     end
-    # E.g.:
-    # set_attribute(model,
-    #    MOI.AutomaticDifferentiationBackend(),
-    #    MathOptSymbolicAD.DefaultBackend(),
-    # )
     return model, recover_allocation
 end
 

src/main/init_parameters.jl

@@ -16,7 +16,7 @@ Returns a `param` dict with the model parameters. These are independent of the g
 - `N::Int64=1`: Number of goods traded in the economy (used for checks in `create_graph()`)
 - `labor_mobility::Any=false`: Switch for labor mobility (true/false or 'partial')
 - `cross_good_congestion::Bool=false`: Switch for cross-good congestion
-- `nu::Float64=1`: Elasticity of substitution b/w goods in transport costs if cross-good congestion
+- `nu::Float64=1.5`: Elasticity of substitution b/w goods in transport costs if cross-good congestion
 - `m::Vector{Float64}=ones(N)`: Vector of weights Nx1 in the cross congestion cost function
 - `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)
@@ -40,7 +40,7 @@ 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, 
-                         labor_mobility = false, cross_good_congestion=false, nu = 1, m = ones(N), 
+                         labor_mobility = false, cross_good_congestion=false, nu = 1.5, m = ones(N), 
                          annealing=true, 
                          verbose = true, duality = false, warm_start = true, 
                          kappa_min = 1e-5, min_iter = 20, max_iter = 200, tol = 1e-5, kwargs...)
@@ -93,7 +93,8 @@ function init_parameters(; alpha = 0.5, beta = 1, gamma = 1, K = 1, sigma = 5, r
     param[:u] = (c, h) -> ((c / alpha)^alpha * (h / (1 - alpha))^(1 - alpha))^(1 - rho) / (1 - rho)
     param[:uprime] = (c, h) -> ((c / alpha)^alpha * (h / (1 - alpha))^(1 - alpha))^-rho * ((c / alpha)^(alpha - 1) * (h / (1 - alpha))^(1 - alpha))
     param[:usecond] = (c, h) -> -rho * ((c / alpha)^alpha * (h / (1 - alpha))^(1 - alpha))^(-rho - 1) * ((c / alpha)^(alpha - 1) * (h / (1 - alpha))^(1 - alpha))^2 + (alpha - 1) / alpha * ((c / alpha)^alpha * (h / (1 - alpha))^(1 - alpha))^-rho * ((c / alpha)^(alpha - 2) * (h / (1 - alpha))^(1 - alpha))
-    param[:uprimeinv] = (x, h) -> alpha * x^(-1 / (1 + alpha * (rho - 1))) * (h / (1 - alpha))^-((1 - alpha) * (rho - 1) / (1 + alpha * (rho - 1)))
+    param[:uprimeinv] = (x, h) -> alpha * x^(-1 / (1 + alpha * (rho - 1))) * (h / (1 - alpha))^((1 - alpha) * (1 - rho) / (1 + alpha * (rho - 1)))
+    param[:uprimeinvprime] = (x, h) -> -alpha / (1 + alpha * (rho - 1)) * x^(-1 / (1 + alpha * (rho - 1)) - 1) * (h / (1 - alpha))^((1 - alpha) * (1 - rho) / (1 + alpha * (rho - 1)))
 
     # Define production function 
     param[:F] = (L, a) -> L^a