diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 0000000..693fd7c --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1,8 @@ +# Configuration file for JuliaFormatter.jl + # For more information, see: https://domluna.github.io/JuliaFormatter.jl/stable/config/ + + always_for_in = true + always_use_return = true + margin = 90 + remove_extra_newlines = true + short_to_long_function_def = true \ No newline at end of file diff --git a/.github/workflows/format_check.yml b/.github/workflows/format_check.yml new file mode 100644 index 0000000..8509ce4 --- /dev/null +++ b/.github/workflows/format_check.yml @@ -0,0 +1,31 @@ +name: format-check + on: + push: + branches: + - master + - release-* + pull_request: + types: [opened, synchronize, reopened] + jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: julia-actions/setup-julia@latest + with: + version: '1' + - uses: actions/checkout@v1 + - name: Format check + shell: julia --color=yes {0} + run: | + using Pkg + Pkg.add(PackageSpec(name="JuliaFormatter", version="1")) + using JuliaFormatter + format("src/MOI_wrapper", verbose=true) + format("test", verbose=true) + out = String(read(Cmd(`git diff`))) + if isempty(out) + exit(0) + end + @error "Some files have not been formatted !!!" + write(stdout, out) + exit(1) \ No newline at end of file diff --git a/CHANGELOG.md b/CHANGELOG.md index 7873e96..fa4e23c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,10 @@ LaplacianOpt.jl Change Log ========================= +### v0.3.1 +- Minor change in `src/variables.jl` +- Added JuliaFormatter.toml and formatting workfow + ### v0.3.0 - Includes adjacency of base and augments graphs in results dictionary - Constraints added to support cycle graphs with max algberaic connectivity using `hamiltonian_cycle` in `graph_type` diff --git a/Project.toml b/Project.toml index 8ca7f15..8cf2cf6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LaplacianOpt" uuid = "bb20392f-64fb-4001-92e8-14b3aedd5a9e" authors = ["Harsha Nagarajan"] -version = "0.3.0" +version = "0.3.1" [deps] DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" diff --git a/docs/make.jl b/docs/make.jl index fcd359a..a9d6841 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -5,27 +5,32 @@ using DocumenterTools: Themes for w in ("light", "dark") header = read(joinpath(@__DIR__, "src/assets/themes/style.scss"), String) theme = read(joinpath(@__DIR__, "src/assets/themes/$(w)defs.scss"), String) - write(joinpath(@__DIR__, "src/assets/themes/$(w).scss"), header*"\n"*theme) + write(joinpath(@__DIR__, "src/assets/themes/$(w).scss"), header * "\n" * theme) end -Themes.compile(joinpath(@__DIR__, "src/assets/themes/light.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-light.css")) -Themes.compile(joinpath(@__DIR__, "src/assets/themes/dark.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-dark.css")) +Themes.compile( + joinpath(@__DIR__, "src/assets/themes/light.scss"), + joinpath(@__DIR__, "src/assets/themes/documenter-light.css"), +) +Themes.compile( + joinpath(@__DIR__, "src/assets/themes/dark.scss"), + joinpath(@__DIR__, "src/assets/themes/documenter-dark.css"), +) makedocs( modules = [LaplacianOpt], sitename = "LaplacianOpt", - format = Documenter.HTML(mathengine = Documenter.MathJax(), - prettyurls = get(ENV, "CI", nothing) == "true"), + format = Documenter.HTML( + mathengine = Documenter.MathJax(), + prettyurls = get(ENV, "CI", nothing) == "true", + ), strict = true, authors = "Harsha Nagarajan", pages = [ - "Introduction" => "index.md", - "Quick Start guide" => "quickguide.md", + "Introduction" => "index.md", + "Quick Start guide" => "quickguide.md", "Function References" => "function_references.md", ], ) -deploydocs( - repo = "github.com/harshangrjn/LaplacianOpt.jl.git", - push_preview = true -) +deploydocs(repo = "github.com/harshangrjn/LaplacianOpt.jl.git", push_preview = true) diff --git a/examples/optimizer.jl b/examples/optimizer.jl index 849deb0..52d54d6 100644 --- a/examples/optimizer.jl +++ b/examples/optimizer.jl @@ -3,18 +3,22 @@ #=====================================# function get_gurobi() - return JuMP.optimizer_with_attributes(Gurobi.Optimizer, - MOI.Silent() => false, - # "MIPFocus" => 3, # Focus on optimality over feasibility - "Presolve" => 1) + return JuMP.optimizer_with_attributes( + Gurobi.Optimizer, + MOI.Silent() => false, + # "MIPFocus" => 3, # Focus on optimality over feasibility + "Presolve" => 1, + ) end function get_cplex() - return JuMP.optimizer_with_attributes(CPLEX.Optimizer, - MOI.Silent() => false, - # "CPX_PARAM_EPGAP" => 1E-4, - # "CPX_PARAM_MIPEMPHASIS" => 2 # Focus on optimality over feasibility - "CPX_PARAM_PREIND" => 1) + return JuMP.optimizer_with_attributes( + CPLEX.Optimizer, + MOI.Silent() => false, + # "CPX_PARAM_EPGAP" => 1E-4, + # "CPX_PARAM_MIPEMPHASIS" => 2 # Focus on optimality over feasibility + "CPX_PARAM_PREIND" => 1, + ) end #======================================# @@ -22,13 +26,11 @@ end #======================================# function get_cbc() - return JuMP.optimizer_with_attributes(Cbc.Optimizer, - MOI.Silent() => false) + return JuMP.optimizer_with_attributes(Cbc.Optimizer, MOI.Silent() => false) end function get_glpk() - return JuMP.optimizer_with_attributes(GLPK.Optimizer, - MOI.Silent() => false) + return JuMP.optimizer_with_attributes(GLPK.Optimizer, MOI.Silent() => false) end #========================================================# @@ -36,19 +38,22 @@ end #========================================================# function get_ipopt() - return JuMP.optimizer_with_attributes(Ipopt.Optimizer, - MOI.Silent() => true, - "sb" => "yes", - "max_iter" => Int(1E4)) + return JuMP.optimizer_with_attributes( + Ipopt.Optimizer, + MOI.Silent() => true, + "sb" => "yes", + "max_iter" => Int(1E4), + ) end - #=================================================================# # Local mixed-integer nonlinear programming solver (open-source) # #=================================================================# function get_juniper() - return JuMP.optimizer_with_attributes(Juniper.Optimizer, - MOI.Silent() => false, - "mip_solver" => get_gurobi(), - "nl_solver" => get_ipopt()) -end \ No newline at end of file + return JuMP.optimizer_with_attributes( + Juniper.Optimizer, + MOI.Silent() => false, + "mip_solver" => get_gurobi(), + "nl_solver" => get_ipopt(), + ) +end diff --git a/examples/run_examples.jl b/examples/run_examples.jl index 1ec3057..8adccd0 100644 --- a/examples/run_examples.jl +++ b/examples/run_examples.jl @@ -18,14 +18,15 @@ lopt_optimizer = get_cplex() #= Option I: Let the package parse a JSON file and obtain the data dictionary (data_dict) Sample JSON format: Check examples/instances -=# +=# function data_I() num_nodes = 8 - instance = 1 + instance = 1 # Data format has to be as given in this JSON file - file_path = joinpath(@__DIR__, "instances/$(num_nodes)_nodes/$(num_nodes)_$(instance).json") + file_path = + joinpath(@__DIR__, "instances/$(num_nodes)_nodes/$(num_nodes)_$(instance).json") data_dict = LOpt.parse_file(file_path) - augment_budget = (num_nodes-1) # spanning tree constraint + augment_budget = (num_nodes - 1) # spanning tree constraint return data_dict, augment_budget end @@ -34,7 +35,7 @@ end num_nodes, adjacency_base_graph and adjacency_augment_graph =# function data_II() - data_dict = Dict{String, Any}() + data_dict = Dict{String,Any}() data_dict["num_nodes"] = 4 data_dict["adjacency_base_graph"] = [0 2 0 0; 2 0 3 0; 0 3 0 4; 0 0 4 0] data_dict["adjacency_augment_graph"] = [0 0 4 8; 0 0 0 7; 4 0 0 0; 8 7 0 0] @@ -47,23 +48,25 @@ end #-------------------------------# data_dict, augment_budget = data_II() -params = Dict{String, Any}( - "data_dict" => data_dict, - "augment_budget" => augment_budget, - "eigen_cuts_full" => true, +params = Dict{String,Any}( + "data_dict" => data_dict, + "augment_budget" => augment_budget, + "eigen_cuts_full" => true, "soc_linearized_cuts" => false, - "eigen_cuts_2minors" => false, - "eigen_cuts_3minors" => false, - "topology_flow_cuts" => true, + "eigen_cuts_2minors" => false, + "eigen_cuts_3minors" => false, + "topology_flow_cuts" => true, # "time_limit" => 3600, - ) +) #----------------------------------------------------------------# # Optimization model and visualize solution (optional) # # Graph plots can be located inside `examples/plots` folder # #----------------------------------------------------------------# -result = LOpt.run_LOpt(params, - lopt_optimizer; - visualize_solution = false, # Make this true to plot the graph solution - visualizing_tool = "tikz", # "graphviz" is another option - display_edge_weights = false) \ No newline at end of file +result = LOpt.run_LOpt( + params, + lopt_optimizer; + visualize_solution = false, # Make this true to plot the graph solution + visualizing_tool = "tikz", # "graphviz" is another option + display_edge_weights = false, +) diff --git a/src/LaplacianOpt.jl b/src/LaplacianOpt.jl index c6139bb..59e0749 100644 --- a/src/LaplacianOpt.jl +++ b/src/LaplacianOpt.jl @@ -9,10 +9,10 @@ import Graphs import TikzGraphs import TikzPictures -const MOI = MathOptInterface -const LA = LinearAlgebra +const MOI = MathOptInterface +const LA = LinearAlgebra const LOpt = LaplacianOpt - + # Create our module level logger (this will get precompiled) const _LOGGER = Memento.getlogger(@__MODULE__) @@ -22,13 +22,16 @@ __init__() = Memento.register(_LOGGER) "Suppresses information and warning messages output by LaplacianOpt, for fine grained control use the Memento package" function silence() - Memento.info(_LOGGER, "Suppressing information and warning messages for the rest of this session. Use the Memento package for more fine-grained control of logging.") - Memento.setlevel!(Memento.getlogger(LaplacianOpt), "error") + Memento.info( + _LOGGER, + "Suppressing information and warning messages for the rest of this session. Use the Memento package for more fine-grained control of logging.", + ) + return Memento.setlevel!(Memento.getlogger(LaplacianOpt), "error") end "allows the user to set the logging level without the need to add Memento" function logger_config!(level) - Memento.config!(Memento.getlogger("LaplacianOpt"), level) + return Memento.config!(Memento.getlogger("LaplacianOpt"), level) end include("data.jl") diff --git a/src/constraints.jl b/src/constraints.jl index 01035a3..59b784f 100644 --- a/src/constraints.jl +++ b/src/constraints.jl @@ -3,80 +3,107 @@ #------------------------------------------------------# function constraint_build_W_var_matrix(lom::LaplacianOptModel) - num_nodes = lom.data["num_nodes"] - num_edges_existing = lom.data["num_edges_existing"] - adjacency_base_graph = lom.data["adjacency_base_graph"] + num_nodes = lom.data["num_nodes"] + num_edges_existing = lom.data["num_edges_existing"] + adjacency_base_graph = lom.data["adjacency_base_graph"] adjacency_augment_graph = lom.data["adjacency_augment_graph"] adjacency_full_graph = adjacency_augment_graph (num_edges_existing > 0) && (adjacency_full_graph += adjacency_base_graph) - + # Diagonal entries - JuMP.@constraint(lom.model, [i=1:num_nodes], lom.variables[:W_var][i,i] == sum(adjacency_full_graph[i,:] .* lom.variables[:z_var][i,:]) - lom.variables[:γ_var]*(num_nodes-1)/num_nodes) + JuMP.@constraint( + lom.model, + [i = 1:num_nodes], + lom.variables[:W_var][i, i] == + sum(adjacency_full_graph[i, :] .* lom.variables[:z_var][i, :]) - + lom.variables[:γ_var] * (num_nodes - 1) / num_nodes + ) # Off-diagonal entries - JuMP.@constraint(lom.model, [i=1:(num_nodes-1), j=(i+1):num_nodes], lom.variables[:W_var][i,j] == -(adjacency_full_graph[i,j] * lom.variables[:z_var][i,j]) + lom.variables[:γ_var]/num_nodes) + JuMP.@constraint( + lom.model, + [i = 1:(num_nodes-1), j = (i+1):num_nodes], + lom.variables[:W_var][i, j] == + -(adjacency_full_graph[i, j] * lom.variables[:z_var][i, j]) + + lom.variables[:γ_var] / num_nodes + ) return end function constraint_single_vertex_cutset(lom::LaplacianOptModel) - num_nodes = lom.data["num_nodes"] - num_edges_existing = lom.data["num_edges_existing"] - adjacency_base_graph = lom.data["adjacency_base_graph"] + num_nodes = lom.data["num_nodes"] + num_edges_existing = lom.data["num_edges_existing"] + adjacency_base_graph = lom.data["adjacency_base_graph"] adjacency_augment_graph = lom.data["adjacency_augment_graph"] - graph_type = lom.data["graph_type"] + graph_type = lom.data["graph_type"] adjacency_full_graph = adjacency_augment_graph (num_edges_existing > 0) && (adjacency_full_graph += adjacency_base_graph) if graph_type == "hamiltonian_cycle" - for i = 1:num_nodes - if sum(adjacency_full_graph[i,:] .> 0) >= 2 - JuMP.@constraint(lom.model, sum(lom.variables[:z_var][i,j] for j in setdiff((1:num_nodes), i)) == 2) + for i in 1:num_nodes + if sum(adjacency_full_graph[i, :] .> 0) >= 2 + JuMP.@constraint( + lom.model, + sum(lom.variables[:z_var][i, j] for j in setdiff((1:num_nodes), i)) == + 2 + ) end end else # for any other graph_type - for i = 1:num_nodes - if sum(adjacency_full_graph[i,:] .> 0) >= 1 - JuMP.@constraint(lom.model, sum(lom.variables[:z_var][i,j] for j in setdiff((1:num_nodes), i)) >= 1) + for i in 1:num_nodes + if sum(adjacency_full_graph[i, :] .> 0) >= 1 + JuMP.@constraint( + lom.model, + sum(lom.variables[:z_var][i, j] for j in setdiff((1:num_nodes), i)) >= + 1 + ) end end end - + return end function constraint_augment_edges_budget(lom::LaplacianOptModel) - num_nodes = lom.data["num_nodes"] - augment_budget = lom.data["augment_budget"] + num_nodes = lom.data["num_nodes"] + augment_budget = lom.data["augment_budget"] adjacency_augment_graph = lom.data["adjacency_augment_graph"] - JuMP.@constraint(lom.model, sum(lom.variables[:z_var][i,j] for i=1:(num_nodes-1), j=(i+1):num_nodes - if adjacency_augment_graph[i,j] > 0) == augment_budget) - + JuMP.@constraint( + lom.model, + sum( + lom.variables[:z_var][i, j] for + i = 1:(num_nodes-1), j = (i+1):num_nodes if adjacency_augment_graph[i, j] > 0 + ) == augment_budget + ) + return end function constraint_lazycallback_wrapper(lom::LaplacianOptModel; max_span_tree = false) - function constraint_lazycallback_cuts(cb_cuts) - status = JuMP.callback_node_status(cb_cuts, lom.model) # if status == MOI.CALLBACK_NODE_STATUS_FRACTIONAL - + # if lom.data["optimizer"] != "glpk" # Memento.info(_LOGGER, "Callback status: solution is fractional") # end if status == MOI.CALLBACK_NODE_STATUS_UNKNOWN - Memento.error(_LOGGER, "Callback status: unknown - solution may not be integer feasible") + Memento.error( + _LOGGER, + "Callback status: unknown - solution may not be integer feasible", + ) - else status == MOI.CALLBACK_NODE_STATUS_INTEGER + else + status == MOI.CALLBACK_NODE_STATUS_INTEGER - if !max_span_tree - if lom.data["eigen_cuts_full"] + if !max_span_tree + if lom.data["eigen_cuts_full"] W_val = JuMP.callback_value.(Ref(cb_cuts), lom.variables[:W_var]) LOpt.constraint_eigen_cuts_on_full_matrix(W_val, cb_cuts, lom) end @@ -97,44 +124,57 @@ function constraint_lazycallback_wrapper(lom::LaplacianOptModel; max_span_tree = end end - if (!(lom.data["is_base_graph_connected"]) && (lom.data["topology_flow_cuts"])) || (max_span_tree) - z_val = abs.(JuMP.callback_value.(Ref(cb_cuts), lom.variables[:z_var])) + if ( + !(lom.data["is_base_graph_connected"]) && + (lom.data["topology_flow_cuts"]) + ) || (max_span_tree) + z_val = abs.(JuMP.callback_value.(Ref(cb_cuts), lom.variables[:z_var])) LOpt.get_rounded_zeros_and_ones!(z_val) LOpt.constraint_topology_flow_cuts(z_val, cb_cuts, lom) end - end + end end if lom.data["lazy_callback_status"] MOI.set(lom.model, MOI.LazyConstraintCallback(), constraint_lazycallback_cuts) end - end -function constraint_soc_cuts_on_2minors(W_val::Matrix{<:Number}, cb_cuts, lom::LaplacianOptModel) - +function constraint_soc_cuts_on_2minors( + W_val::Matrix{<:Number}, + cb_cuts, + lom::LaplacianOptModel, +) num_nodes = lom.data["num_nodes"] - for i = 1:(num_nodes-1) - for j = (i+1):num_nodes - W_val_minor = W_val[i,j]^2 - W_val[i,i] * W_val[j,j] - if W_val_minor >= lom.data["tol_zero"] - coefficient_vec = [-(W_val[i,j])^2 - 2 * W_val[i,j] * W_val[i,i]] - - variable_vec = [lom.variables[:W_var][i,i], lom.variables[:W_var][i,j]] - + for i in 1:(num_nodes-1) + for j in (i+1):num_nodes + W_val_minor = W_val[i, j]^2 - W_val[i, i] * W_val[j, j] + if W_val_minor >= lom.data["tol_zero"] + coefficient_vec = [ + -(W_val[i, j])^2 + 2 * W_val[i, j] * W_val[i, i] + ] + + variable_vec = [lom.variables[:W_var][i, i], lom.variables[:W_var][i, j]] + # Taylor's approximation of rotated SOC constraint (x^2 <= y*z) - if !isapprox(W_val[i,i], 0, atol=1E-6) - con = JuMP.@build_constraint(coefficient_vec[1]*variable_vec[1] + coefficient_vec[2]*variable_vec[2] - <= W_val[i,i]^2 * lom.variables[:W_var][j,j]) - - MOI.submit(lom.model, MOI.LazyConstraint(cb_cuts), con) + if !isapprox(W_val[i, i], 0, atol = 1E-6) + con = JuMP.@build_constraint( + coefficient_vec[1] * variable_vec[1] + + coefficient_vec[2] * variable_vec[2] <= + W_val[i, i]^2 * lom.variables[:W_var][j, j] + ) + + MOI.submit(lom.model, MOI.LazyConstraint(cb_cuts), con) end - + if lom.data["lazycuts_logging"] - Memento.info(_LOGGER, "Polyhedral relaxation cuts: 2x2 minor linearized SOC") - end + Memento.info( + _LOGGER, + "Polyhedral relaxation cuts: 2x2 minor linearized SOC", + ) + end end end end @@ -142,35 +182,44 @@ function constraint_soc_cuts_on_2minors(W_val::Matrix{<:Number}, cb_cuts, lom::L return end -function constraint_eigen_cuts_on_full_matrix(W_val::Matrix{<:Number}, cb_cuts, lom::LaplacianOptModel) - +function constraint_eigen_cuts_on_full_matrix( + W_val::Matrix{<:Number}, + cb_cuts, + lom::LaplacianOptModel, +) LOpt._add_eigen_cut_lazy(W_val, cb_cuts, lom, collect(1:lom.data["num_nodes"])) return end -function constraint_eigen_cuts_on_2minors(W_val::Matrix{<:Number}, cb_cuts, lom::LaplacianOptModel) - +function constraint_eigen_cuts_on_2minors( + W_val::Matrix{<:Number}, + cb_cuts, + lom::LaplacianOptModel, +) num_nodes = lom.data["num_nodes"] - + if num_nodes >= 3 - for i = 1:(num_nodes-1) - for j = (i+1):num_nodes - LOpt._add_eigen_cut_lazy(W_val, cb_cuts, lom, [i,j]) + for i in 1:(num_nodes-1) + for j in (i+1):num_nodes + LOpt._add_eigen_cut_lazy(W_val, cb_cuts, lom, [i, j]) end end end return end -function constraint_eigen_cuts_on_3minors(W_val::Matrix{<:Number}, cb_cuts, lom::LaplacianOptModel) - +function constraint_eigen_cuts_on_3minors( + W_val::Matrix{<:Number}, + cb_cuts, + lom::LaplacianOptModel, +) num_nodes = lom.data["num_nodes"] if num_nodes >= 4 - for i = 1:(num_nodes-2) - for j = (i+1):num_nodes - for k = (j+1):num_nodes - LOpt._add_eigen_cut_lazy(W_val, cb_cuts, lom, [i,j,k]) + for i in 1:(num_nodes-2) + for j in (i+1):num_nodes + for k in (j+1):num_nodes + LOpt._add_eigen_cut_lazy(W_val, cb_cuts, lom, [i, j, k]) end end end @@ -178,13 +227,21 @@ function constraint_eigen_cuts_on_3minors(W_val::Matrix{<:Number}, cb_cuts, lom: return end -function _add_eigen_cut_lazy(W_val::Matrix{<:Number}, cb_cuts, lom::LaplacianOptModel, idx::Vector{Int64}) +function _add_eigen_cut_lazy( + W_val::Matrix{<:Number}, + cb_cuts, + lom::LaplacianOptModel, + idx::Vector{Int64}, +) W_val_minor = W_val[idx, idx] W_var_minor = lom.variables[:W_var][idx, idx] - violated_eigen_vec = LOpt._violated_eigen_vector(W_val_minor, tol = lom.data["tol_psd"]) - if violated_eigen_vec !== nothing - con = JuMP.@build_constraint(violated_eigen_vec' * W_var_minor * violated_eigen_vec >= 0) + violated_eigen_vec = + LOpt._violated_eigen_vector(W_val_minor, tol = lom.data["tol_psd"]) + if violated_eigen_vec !== nothing + con = JuMP.@build_constraint( + violated_eigen_vec' * W_var_minor * violated_eigen_vec >= 0 + ) MOI.submit(lom.model, MOI.LazyConstraint(cb_cuts), con) if lom.data["lazycuts_logging"] @@ -197,88 +254,130 @@ function _add_eigen_cut_lazy(W_val::Matrix{<:Number}, cb_cuts, lom::LaplacianOpt end end end - end -function constraint_topology_flow_cuts(z_val::Matrix{<:Number}, cb_cuts, lom::LaplacianOptModel) - - num_nodes = lom.data["num_nodes"] +function constraint_topology_flow_cuts( + z_val::Matrix{<:Number}, + cb_cuts, + lom::LaplacianOptModel, +) + num_nodes = lom.data["num_nodes"] adjacency_augment_graph = lom.data["adjacency_augment_graph"] - graph_type = lom.data["graph_type"] - num_edges_existing = lom.data["num_edges_existing"] - augment_budget = lom.data["augment_budget"] - cc_lazy = Graphs.connected_components(Graphs.SimpleGraph(abs.(z_val))) - - is_spanning_tree = false - if (num_edges_existing == 0) && (augment_budget == (num_nodes-1)) + graph_type = lom.data["graph_type"] + num_edges_existing = lom.data["num_edges_existing"] + augment_budget = lom.data["augment_budget"] + cc_lazy = Graphs.connected_components(Graphs.SimpleGraph(abs.(z_val))) + + is_spanning_tree = false + if (num_edges_existing == 0) && (augment_budget == (num_nodes - 1)) is_spanning_tree = true end max_cc = 5 # increase this to any greater integer value and it works if length(cc_lazy) > max_cc - Memento.info(_LOGGER, "Polyhedral relaxation: flow cuts not added for integer solutions with $(length(cc_lazy)) connected components") - elseif length(cc_lazy) == 1 + Memento.info( + _LOGGER, + "Polyhedral relaxation: flow cuts not added for integer solutions with $(length(cc_lazy)) connected components", + ) + elseif length(cc_lazy) == 1 return end - min_cc_size = minimum(length.(cc_lazy)) - min_cc_loc = argmin(length.(cc_lazy)) - cc_size_threshold = ceil(lom.data["num_nodes"]/4) - + min_cc_size = minimum(length.(cc_lazy)) + min_cc_loc = argmin(length.(cc_lazy)) + cc_size_threshold = ceil(lom.data["num_nodes"] / 4) + # Subtour elimination (for cycle or tree graphs) - if (2 <= min_cc_size <= cc_size_threshold) && ((graph_type == "hamiltonian_cycle") || (is_spanning_tree)) - for k = 1:(length(cc_lazy)) + if (2 <= min_cc_size <= cc_size_threshold) && + ((graph_type == "hamiltonian_cycle") || (is_spanning_tree)) + for k in 1:(length(cc_lazy)) cc_lazy_size = length(cc_lazy[k]) if cc_lazy_size <= cc_size_threshold - cc = cc_lazy[k] - con = JuMP.@build_constraint(sum(lom.variables[:z_var][cc[i],cc[j]] for i=1:(cc_lazy_size-1), j=(i+1):cc_lazy_size - if !(isapprox(adjacency_augment_graph[i,j], 0, atol=1E-6))) <= (cc_lazy_size - 1)) + cc = cc_lazy[k] + con = JuMP.@build_constraint( + sum( + lom.variables[:z_var][cc[i], cc[j]] for + i = 1:(cc_lazy_size-1), j = (i+1):cc_lazy_size if + !(isapprox(adjacency_augment_graph[i, j], 0, atol = 1E-6)) + ) <= (cc_lazy_size - 1) + ) MOI.submit(lom.model, MOI.LazyConstraint(cb_cuts), con) end end - - # Flow cuts for connected components - elseif length(cc_lazy) in 2:max_cc - num_edges_cutset = 1 + # Flow cuts for connected components + elseif length(cc_lazy) in 2:max_cc + num_edges_cutset = 1 if graph_type == "hamiltonian_cycle" num_edges_cutset = 2 end - for k = 1:(length(cc_lazy)-1) + for k in 1:(length(cc_lazy)-1) # From cutset cutset_f = cc_lazy[k] # To cutset ii = setdiff(1:length(cc_lazy), k) cutset_t = reduce(vcat, cc_lazy[ii]) - if LOpt._is_flow_cut_valid(cutset_f, cutset_t, num_edges_cutset, adjacency_augment_graph) - con = JuMP.@build_constraint(sum(lom.variables[:z_var][i,j] for i in cutset_f, j in cutset_t - if !(isapprox(adjacency_augment_graph[i,j], 0, atol=1E-6))) >= num_edges_cutset) + if LOpt._is_flow_cut_valid( + cutset_f, + cutset_t, + num_edges_cutset, + adjacency_augment_graph, + ) + con = JuMP.@build_constraint( + sum( + lom.variables[:z_var][i, j] for + i in cutset_f, j in cutset_t if + !(isapprox(adjacency_augment_graph[i, j], 0, atol = 1E-6)) + ) >= num_edges_cutset + ) MOI.submit(lom.model, MOI.LazyConstraint(cb_cuts), con) end end if lom.data["lazycuts_logging"] - Memento.info(_LOGGER, "Polyhedral relaxation cuts: graph connectivity (#cc = $(length(cc_lazy)))") - end + Memento.info( + _LOGGER, + "Polyhedral relaxation cuts: graph connectivity (#cc = $(length(cc_lazy)))", + ) + end end - return + return end function constraint_topology_multi_commodity_flow(lom::LaplacianOptModel) - num_nodes = lom.data["num_nodes"] - + # Flow balance on source node - JuMP.@constraint(lom.model, [k=1:(num_nodes-1)], sum(lom.variables[:flow_var][1,2:num_nodes,k]) - sum(lom.variables[:flow_var][2:num_nodes,1,k]) == 1) + JuMP.@constraint( + lom.model, + [k = 1:(num_nodes-1)], + sum(lom.variables[:flow_var][1, 2:num_nodes, k]) - + sum(lom.variables[:flow_var][2:num_nodes, 1, k]) == 1 + ) # Flow balance on target node - JuMP.@constraint(lom.model, [k=1:(num_nodes-1)], sum(lom.variables[:flow_var][:,k+1,k]) - sum(lom.variables[:flow_var][k+1,:,k]) == 1) #(k-th commodity corresponds to k+1th node) + JuMP.@constraint( + lom.model, + [k = 1:(num_nodes-1)], + sum(lom.variables[:flow_var][:, k+1, k]) - + sum(lom.variables[:flow_var][k+1, :, k]) == 1 + ) #(k-th commodity corresponds to k+1th node) # Flow balance on remaining nodes - JuMP.@constraint(lom.model, [k=1:(num_nodes-1), v=(2:num_nodes); k!=v-1], sum(lom.variables[:flow_var][:,v,k]) - sum(lom.variables[:flow_var][v,:,k]) == 0) + JuMP.@constraint( + lom.model, + [k = 1:(num_nodes-1), v = (2:num_nodes); k != v - 1], + sum(lom.variables[:flow_var][:, v, k]) - sum(lom.variables[:flow_var][v, :, k]) == + 0 + ) # Flow <= capacity - JuMP.@constraint(lom.model, [k=1:(num_nodes-1), i=1:num_nodes, j=1:num_nodes, kdash=1:(num_nodes-1)], lom.variables[:flow_var][i,j,k] + lom.variables[:flow_var][j,i,kdash]' <= lom.variables[:z_var][i,j]) - + JuMP.@constraint( + lom.model, + [k = 1:(num_nodes-1), i = 1:num_nodes, j = 1:num_nodes, kdash = 1:(num_nodes-1)], + lom.variables[:flow_var][i, j, k] + lom.variables[:flow_var][j, i, kdash]' <= + lom.variables[:z_var][i, j] + ) + return -end \ No newline at end of file +end diff --git a/src/data.jl b/src/data.jl index b54a559..e4b20eb 100644 --- a/src/data.jl +++ b/src/data.jl @@ -4,11 +4,10 @@ Given the input params, this function preprocesses the data, catches any error and missing data, and outputs a detailed dictionary which forms the basis for building an optimization model. """ -function get_data(params::Dict{String, Any}) - +function get_data(params::Dict{String,Any}) if "data_dict" in keys(params) data_dict = params["data_dict"] - else + else Memento.error(_LOGGER, "Data dictionary has to be specified in input params") end @@ -21,15 +20,27 @@ function get_data(params::Dict{String, Any}) augment_budget = params["augment_budget"] else augment_budget = data_dict["num_edges_to_augment"] - Memento.info(_LOGGER, "Setting edge augmentation budget to $(data_dict["num_edges_to_augment"])") + Memento.info( + _LOGGER, + "Setting edge augmentation budget to $(data_dict["num_edges_to_augment"])", + ) end else augment_budget = data_dict["num_edges_to_augment"] - Memento.info(_LOGGER, "Setting edge augmentation budget to $(data_dict["num_edges_to_augment"])") - end - - Memento.info(_LOGGER, "Number of edges in the base graph: $(data_dict["num_edges_existing"])") - Memento.info(_LOGGER, "Number of candidate edges to augment: $(data_dict["num_edges_to_augment"])") + Memento.info( + _LOGGER, + "Setting edge augmentation budget to $(data_dict["num_edges_to_augment"])", + ) + end + + Memento.info( + _LOGGER, + "Number of edges in the base graph: $(data_dict["num_edges_existing"])", + ) + Memento.info( + _LOGGER, + "Number of candidate edges to augment: $(data_dict["num_edges_to_augment"])", + ) Memento.info(_LOGGER, "Augment budget: $(augment_budget)") # Solution type @@ -120,7 +131,10 @@ function get_data(params::Dict{String, Any}) end if topology_multi_commodity && data_dict["is_base_graph_connected"] topology_multi_commodity = false - Memento.info(_LOGGER, "Deactivating topology multi commodity formulation as the base graph is connected") + Memento.info( + _LOGGER, + "Deactivating topology multi commodity formulation as the base graph is connected", + ) end if topology_multi_commodity Memento.info(_LOGGER, "Applying topology multi commodity constraints") @@ -132,7 +146,10 @@ function get_data(params::Dict{String, Any}) end if topology_flow_cuts && data_dict["is_base_graph_connected"] topology_flow_cuts = false - Memento.info(_LOGGER, "Deactivating topology flow cuts as the base graph is connected") + Memento.info( + _LOGGER, + "Deactivating topology flow cuts as the base graph is connected", + ) elseif topology_multi_commodity topology_flow_cuts = false end @@ -140,7 +157,11 @@ function get_data(params::Dict{String, Any}) Memento.info(_LOGGER, "Applying topology flow cuts") end - if eigen_cuts_full || topology_flow_cuts || soc_linearized_cuts || eigen_cuts_2minors || eigen_cuts_3minors + if eigen_cuts_full || + topology_flow_cuts || + soc_linearized_cuts || + eigen_cuts_2minors || + eigen_cuts_3minors lazy_callback_status = true end @@ -158,27 +179,29 @@ function get_data(params::Dict{String, Any}) time_limit = 10800 end - data = Dict{String, Any}("num_nodes" => num_nodes, - "num_edges_existing" => data_dict["num_edges_existing"], - "num_edges_to_augment" => data_dict["num_edges_to_augment"], - "augment_budget" => augment_budget, - "adjacency_base_graph" => data_dict["adjacency_base_graph"], - "adjacency_augment_graph" => data_dict["adjacency_augment_graph"], - "is_base_graph_connected" => data_dict["is_base_graph_connected"], - "solution_type" => solution_type, - "graph_type" => graph_type, - "tol_zero" => tol_zero, - "tol_psd" => tol_psd, - "eigen_cuts_full" => eigen_cuts_full, - "soc_linearized_cuts" => soc_linearized_cuts, - "eigen_cuts_2minors" => eigen_cuts_2minors, - "eigen_cuts_3minors" => eigen_cuts_3minors, - "lazycuts_logging" => lazycuts_logging, - "topology_flow_cuts" => topology_flow_cuts, - "topology_multi_commodity" => topology_multi_commodity, - "lazy_callback_status" => lazy_callback_status, - "relax_integrality" => relax_integrality, - "time_limit" => time_limit) + data = Dict{String,Any}( + "num_nodes" => num_nodes, + "num_edges_existing" => data_dict["num_edges_existing"], + "num_edges_to_augment" => data_dict["num_edges_to_augment"], + "augment_budget" => augment_budget, + "adjacency_base_graph" => data_dict["adjacency_base_graph"], + "adjacency_augment_graph" => data_dict["adjacency_augment_graph"], + "is_base_graph_connected" => data_dict["is_base_graph_connected"], + "solution_type" => solution_type, + "graph_type" => graph_type, + "tol_zero" => tol_zero, + "tol_psd" => tol_psd, + "eigen_cuts_full" => eigen_cuts_full, + "soc_linearized_cuts" => soc_linearized_cuts, + "eigen_cuts_2minors" => eigen_cuts_2minors, + "eigen_cuts_3minors" => eigen_cuts_3minors, + "lazycuts_logging" => lazycuts_logging, + "topology_flow_cuts" => topology_flow_cuts, + "topology_multi_commodity" => topology_multi_commodity, + "lazy_callback_status" => lazy_callback_status, + "relax_integrality" => relax_integrality, + "time_limit" => time_limit, + ) # Optimizer if "optimizer" in keys(params) @@ -192,17 +215,17 @@ end function parse_file(file_path::String) data_dict = JSON.parsefile(file_path) - + if "num_nodes" in keys(data_dict) num_nodes = data_dict["num_nodes"] else Memento.error(_LOGGER, "Number of nodes is missing in the input data file") end - adjacency_base_graph = Matrix{Float64}(zeros(num_nodes, num_nodes)) + adjacency_base_graph = Matrix{Float64}(zeros(num_nodes, num_nodes)) adjacency_augment_graph = Matrix{Float64}(zeros(num_nodes, num_nodes)) - for k=1:length(data_dict["edges_existing"]) + for k in 1:length(data_dict["edges_existing"]) i, j = data_dict["edges_existing"][k][1] w_ij = data_dict["edges_existing"][k][2] @@ -211,12 +234,12 @@ function parse_file(file_path::String) end LOpt._catch_data_input_error(num_nodes, i, j, w_ij) - - adjacency_base_graph[i,j] = w_ij - adjacency_base_graph[j,i] = w_ij + + adjacency_base_graph[i, j] = w_ij + adjacency_base_graph[j, i] = w_ij end - for k=1:length(data_dict["edges_to_augment"]) + for k in 1:length(data_dict["edges_to_augment"]) i, j = data_dict["edges_to_augment"][k][1] w_ij = data_dict["edges_to_augment"][k][2] @@ -225,44 +248,54 @@ function parse_file(file_path::String) end LOpt._catch_data_input_error(num_nodes, i, j, w_ij) - - adjacency_augment_graph[i,j] = w_ij - adjacency_augment_graph[j,i] = w_ij + + adjacency_augment_graph[i, j] = w_ij + adjacency_augment_graph[j, i] = w_ij end - data_dict_new = Dict{String, Any}("num_nodes" => num_nodes, - "adjacency_base_graph" => adjacency_base_graph, - "adjacency_augment_graph" => adjacency_augment_graph) - + data_dict_new = Dict{String,Any}( + "num_nodes" => num_nodes, + "adjacency_base_graph" => adjacency_base_graph, + "adjacency_augment_graph" => adjacency_augment_graph, + ) + return data_dict_new end -function _pre_process_data(data_dict::Dict{String, Any}) - num_edges_existing = 0 - num_edges_to_augment = 0 +function _pre_process_data(data_dict::Dict{String,Any}) + num_edges_existing = 0 + num_edges_to_augment = 0 - num_nodes = data_dict["num_nodes"] - adjacency_base_graph = data_dict["adjacency_base_graph"] + num_nodes = data_dict["num_nodes"] + adjacency_base_graph = data_dict["adjacency_base_graph"] adjacency_augment_graph = data_dict["adjacency_augment_graph"] - for i=1:(num_nodes-1) - for j=(i+1):num_nodes - if isapprox(abs(adjacency_base_graph[i,j]), 0, atol = 1E-6) - adjacency_base_graph[i,j] = 0 + for i in 1:(num_nodes-1) + for j in (i+1):num_nodes + if isapprox(abs(adjacency_base_graph[i, j]), 0, atol = 1E-6) + adjacency_base_graph[i, j] = 0 end - if isapprox(abs(adjacency_augment_graph[i,j]), 0, atol = 1E-6) - adjacency_augment_graph[i,j] = 0 + if isapprox(abs(adjacency_augment_graph[i, j]), 0, atol = 1E-6) + adjacency_augment_graph[i, j] = 0 end - if !(isapprox(adjacency_base_graph[i,j], adjacency_base_graph[j,i], atol = 1E-5)) + if !(isapprox( + adjacency_base_graph[i, j], + adjacency_base_graph[j, i], + atol = 1E-5, + )) Memento.error("Adjacency matrix of the base graph has to be symmetric") end - if !(isapprox(adjacency_augment_graph[i,j], adjacency_augment_graph[j,i], atol = 1E-5)) + if !(isapprox( + adjacency_augment_graph[i, j], + adjacency_augment_graph[j, i], + atol = 1E-5, + )) Memento.error("Adjacency matrix of the augment graph has to be symmetric") end - if !(isapprox(adjacency_base_graph[i,j], 0, atol = 1E-6)) + if !(isapprox(adjacency_base_graph[i, j], 0, atol = 1E-6)) num_edges_existing += 1 end - if !(isapprox(adjacency_augment_graph[i,j], 0, atol = 1E-6)) + if !(isapprox(adjacency_augment_graph[i, j], 0, atol = 1E-6)) num_edges_to_augment += 1 end end @@ -271,17 +304,21 @@ function _pre_process_data(data_dict::Dict{String, Any}) # Base graph connectivity is_base_graph_connected = false if num_edges_existing > 0 - !(isapprox(abs(LOpt.algebraic_connectivity(adjacency_base_graph)), 0, atol=1E-6)) && (is_base_graph_connected = true) + !(isapprox( + abs(LOpt.algebraic_connectivity(adjacency_base_graph)), + 0, + atol = 1E-6, + )) && (is_base_graph_connected = true) end - data_dict_new = Dict{String, Any}() - data_dict_new["num_nodes"] = num_nodes - data_dict_new["adjacency_base_graph"] = adjacency_base_graph + data_dict_new = Dict{String,Any}() + data_dict_new["num_nodes"] = num_nodes + data_dict_new["adjacency_base_graph"] = adjacency_base_graph data_dict_new["adjacency_augment_graph"] = adjacency_augment_graph - data_dict_new["num_edges_existing"] = num_edges_existing - data_dict_new["num_edges_to_augment"] = num_edges_to_augment + data_dict_new["num_edges_existing"] = num_edges_existing + data_dict_new["num_edges_to_augment"] = num_edges_to_augment data_dict_new["is_base_graph_connected"] = is_base_graph_connected - + return data_dict_new end @@ -293,10 +330,13 @@ any input data error in the JSON file. """ function _catch_data_input_error(num_nodes::Int64, i::Int64, j::Int64, w_ij::Number) if (i > num_nodes) || (j > num_nodes) - Memento.error(_LOGGER, "Node pair ($i,$j) does not match with total number of nodes, $num_nodes") + Memento.error( + _LOGGER, + "Node pair ($i,$j) does not match with total number of nodes, $num_nodes", + ) end - if !(isapprox(abs(w_ij), 0, atol=1E-6)) & (w_ij < 0) + if !(isapprox(abs(w_ij), 0, atol = 1E-6)) & (w_ij < 0) Memento.error(_LOGGER, "Graphs with negative weights are not supported") end end @@ -307,39 +347,58 @@ end Given the pre-processed data dictionary, this function detects any infeasibility before building the optimization model. """ -function _detect_infeasbility_in_data(data::Dict{String, Any}) - +function _detect_infeasbility_in_data(data::Dict{String,Any}) num_nodes = data["num_nodes"] if data["num_edges_to_augment"] == 0 - Memento.error(_LOGGER, "At least one edge has to be specified to be able to augment to the base graph") - # Assuming undirected graph - elseif data["num_edges_existing"] == num_nodes*(num_nodes-1)/2 - Memento.error(_LOGGER, "Input graph is already a complete graph; augmentation is unnecessary") - elseif (data["num_edges_existing"] == 0) && (data["augment_budget"] < (num_nodes-1)) - Memento.error(_LOGGER, "Detected trivial solutions with disconnected graphs. `augment_budget` may be insufficient.") + Memento.error( + _LOGGER, + "At least one edge has to be specified to be able to augment to the base graph", + ) + # Assuming undirected graph + elseif data["num_edges_existing"] == num_nodes * (num_nodes - 1) / 2 + Memento.error( + _LOGGER, + "Input graph is already a complete graph; augmentation is unnecessary", + ) + elseif (data["num_edges_existing"] == 0) && (data["augment_budget"] < (num_nodes - 1)) + Memento.error( + _LOGGER, + "Detected trivial solutions with disconnected graphs. `augment_budget` may be insufficient.", + ) elseif !isinteger(data["augment_budget"]) || (data["augment_budget"] < -1E-6) Memento.error(_LOGGER, "Edge augmentation budget has to be a positive integer") - end + end # Detect mutually exclusive sets of edges between base and augment graph - if maximum((data["adjacency_augment_graph"] .> 0) + (data["adjacency_base_graph"] .> 0)) > 1 - Memento.error(_LOGGER, "Edge sets of base and augment graphs have to be mutually exclusive") - end + if maximum( + (data["adjacency_augment_graph"] .> 0) + (data["adjacency_base_graph"] .> 0), + ) > 1 + Memento.error( + _LOGGER, + "Edge sets of base and augment graphs have to be mutually exclusive", + ) + end # Detect free vertices A = data["adjacency_augment_graph"] + data["adjacency_base_graph"] - for i=1:num_nodes - if isapprox(sum(A[i,:]), 0, atol=1E-6) - Memento.error(_LOGGER, "Detected trivial solutions with disconnected graphs due to free vertices.") + for i in 1:num_nodes + if isapprox(sum(A[i, :]), 0, atol = 1E-6) + Memento.error( + _LOGGER, + "Detected trivial solutions with disconnected graphs due to free vertices.", + ) end end # Detect tour infeasibility - if (data["graph_type"] == "hamiltonian_cycle") && (data["num_edges_existing"] == 0) - if !(data["num_edges_to_augment"] >= data["num_nodes"]) || !(data["augment_budget"] == data["num_nodes"]) - Memento.error(_LOGGER, "Detected infeasibility due to the number of augmentation edges incompatible for a hamiltonian cycle") + if (data["graph_type"] == "hamiltonian_cycle") && (data["num_edges_existing"] == 0) + if !(data["num_edges_to_augment"] >= data["num_nodes"]) || + !(data["augment_budget"] == data["num_nodes"]) + Memento.error( + _LOGGER, + "Detected infeasibility due to the number of augmentation edges incompatible for a hamiltonian cycle", + ) end end - -end \ No newline at end of file +end diff --git a/src/log.jl b/src/log.jl index a260da4..8103a2e 100644 --- a/src/log.jl +++ b/src/log.jl @@ -1,110 +1,132 @@ -function visualize_solution(results::Dict{String, Any}, data::Dict{String, Any}; - visualizing_tool = "tikz", - plot_file_format = "pdf", - display_edge_weights = true) - - num_edges_existing = data["num_edges_existing"] - adjacency_base_graph = data["adjacency_base_graph"] +function visualize_solution( + results::Dict{String,Any}, + data::Dict{String,Any}; + visualizing_tool = "tikz", + plot_file_format = "pdf", + display_edge_weights = true, +) + num_edges_existing = data["num_edges_existing"] + adjacency_base_graph = data["adjacency_base_graph"] adjacency_augment_graph = data["adjacency_augment_graph"] - adjacency_full_graph = adjacency_augment_graph + adjacency_full_graph = adjacency_augment_graph (num_edges_existing > 0) && (adjacency_full_graph += adjacency_base_graph) if results["primal_status"] != MOI.FEASIBLE_POINT - Memento.error(_LOGGER, "Non-feasible primal status. Graph solution may not be exact") + Memento.error( + _LOGGER, + "Non-feasible primal status. Graph solution may not be exact", + ) end if !data["relax_integrality"] Memento.info(_LOGGER, "Plotting the graph of integral solution") - + if visualizing_tool == "tikz" - LOpt.plot_tikzgraph(results["solution"]["z_var"] .* adjacency_full_graph, - plot_file_format = plot_file_format, - display_edge_weights = display_edge_weights) + LOpt.plot_tikzgraph( + results["solution"]["z_var"] .* adjacency_full_graph, + plot_file_format = plot_file_format, + display_edge_weights = display_edge_weights, + ) elseif visualizing_tool == "graphviz" - LOpt.plot_graphviz(results["solution"]["z_var"] .* adjacency_full_graph, - display_edge_weights = display_edge_weights) + LOpt.plot_graphviz( + results["solution"]["z_var"] .* adjacency_full_graph, + display_edge_weights = display_edge_weights, + ) end - else - Memento.info(_LOGGER, "Cannot plot as the obtained solutions are non-integral; fractional values can be found in the results dictionary") + else + Memento.info( + _LOGGER, + "Cannot plot as the obtained solutions are non-integral; fractional values can be found in the results dictionary", + ) return end end -function plot_tikzgraph(adjacency_matrix::Matrix{<:Number}; - plot_file_format = "pdf", - display_edge_weights = false) - - num_nodes = size(adjacency_matrix)[1] +function plot_tikzgraph( + adjacency_matrix::Matrix{<:Number}; + plot_file_format = "pdf", + display_edge_weights = false, +) + num_nodes = size(adjacency_matrix)[1] solution_graph = Graphs.SimpleGraph(num_nodes) - edge_labels = Dict{Tuple{Int64, Int64}, String}() - - for i=1:(num_nodes-1) - for j=(i+1):num_nodes + edge_labels = Dict{Tuple{Int64,Int64},String}() - if !isapprox(abs(adjacency_matrix[i,j]), 0, atol=1E-6) + for i in 1:(num_nodes-1) + for j in (i+1):num_nodes + if !isapprox(abs(adjacency_matrix[i, j]), 0, atol = 1E-6) Graphs.add_edge!(solution_graph, i, j) - edge_labels[(i,j)] = string(ceil(adjacency_matrix[i,j], digits=2)) + edge_labels[(i, j)] = string(ceil(adjacency_matrix[i, j], digits = 2)) end - end end if display_edge_weights # Plot with edge weights - graphs do not look great - t = TikzGraphs.plot(solution_graph, - edge_labels=edge_labels, - node_style="draw, rounded corners, fill=blue!10", - TikzGraphs.Layouts.SpringElectrical()) + t = TikzGraphs.plot( + solution_graph, + edge_labels = edge_labels, + node_style = "draw, rounded corners, fill=blue!10", + TikzGraphs.Layouts.SpringElectrical(), + ) else - t = TikzGraphs.plot(solution_graph, - node_style="draw, rounded corners, fill=blue!10", - TikzGraphs.Layouts.SpringElectrical()) + t = TikzGraphs.plot( + solution_graph, + node_style = "draw, rounded corners, fill=blue!10", + TikzGraphs.Layouts.SpringElectrical(), + ) end - file_path = joinpath(dirname(pathof(LaplacianOpt)),"..", "examples/plots/plot_$(num_nodes)") + file_path = + joinpath(dirname(pathof(LaplacianOpt)), "..", "examples/plots/plot_$(num_nodes)") if plot_file_format == "pdf" TikzPictures.save(TikzPictures.PDF(file_path), t) elseif plot_file_format == "tex" TikzPictures.save(TikzPictures.TEX(file_path), t) end - end -function plot_graphviz(adjacency_matrix::Matrix{<:Number}; - display_edge_weights = true) - +function plot_graphviz(adjacency_matrix::Matrix{<:Number}; display_edge_weights = true) num_nodes = size(adjacency_matrix)[1] - file_path = joinpath(dirname(pathof(LaplacianOpt)),"..", "examples/plots/plot_$(num_nodes).dot") + file_path = joinpath( + dirname(pathof(LaplacianOpt)), + "..", + "examples/plots/plot_$(num_nodes).dot", + ) open(file_path, "w") do file - write(file, "graph G { \n") write(file, "layout=neato; \n") write(file, "size=\"10,5\"; \n") # write(file, "rankdir=LR; \n") - write(file, "node [fontname=\"Helvetica\", fontsize=20, shape = circle, width=0.4, fixedsize=true, style=\"filled\", fillcolor=\"0.650 0.200 1.000\"]; \n") - - for i=1:(num_nodes-1) - for j=(i+1):num_nodes - if !isapprox(abs(adjacency_matrix[i,j]), 0, atol=1E-6) - if display_edge_weights - w_ij = string(ceil(adjacency_matrix[i,j], digits=3)) - write(file, "$i -- $j [label = \"$w_ij\", fontsize=9, fontname=\"Helvetica\"]; \n") - else + write( + file, + "node [fontname=\"Helvetica\", fontsize=20, shape = circle, width=0.4, fixedsize=true, style=\"filled\", fillcolor=\"0.650 0.200 1.000\"]; \n", + ) + + for i in 1:(num_nodes-1) + for j in (i+1):num_nodes + if !isapprox(abs(adjacency_matrix[i, j]), 0, atol = 1E-6) + if display_edge_weights + w_ij = string(ceil(adjacency_matrix[i, j], digits = 3)) + write( + file, + "$i -- $j [label = \"$w_ij\", fontsize=9, fontname=\"Helvetica\"]; \n", + ) + else write(file, "$i -- $j [fontsize=9, fontname=\"Helvetica\"]; \n") end end end end - write(file, "}") + return write(file, "}") end # Further to convert the .dot in to a pdf file, use the following command in the terminal # execute the dot_to_pdf.sh script file in the examples/plots folder -end \ No newline at end of file +end diff --git a/src/lopt_model.jl b/src/lopt_model.jl index 3eb948d..8509d49 100644 --- a/src/lopt_model.jl +++ b/src/lopt_model.jl @@ -1,8 +1,7 @@ -function build_LOModel(data::Dict{String, Any}) - +function build_LOModel(data::Dict{String,Any}) lom = LaplacianOptModel(data) LOpt.variable_LOModel(lom) - + if lom.data["solution_type"] == "exact" LOpt.constraint_LOModel(lom) end @@ -13,7 +12,6 @@ function build_LOModel(data::Dict{String, Any}) end function variable_LOModel(lom::LaplacianOptModel) - LOpt.variable_lifted_W_matrix(lom) LOpt.variable_edge_onoff(lom) lom.data["topology_multi_commodity"] && LOpt.variable_multi_commodity_flow(lom) @@ -23,13 +21,13 @@ function variable_LOModel(lom::LaplacianOptModel) end function constraint_LOModel(lom::LaplacianOptModel) - LOpt.constraint_build_W_var_matrix(lom) LOpt.constraint_single_vertex_cutset(lom) LOpt.constraint_augment_edges_budget(lom) - lom.data["topology_multi_commodity"] && LOpt.constraint_topology_multi_commodity_flow(lom) + lom.data["topology_multi_commodity"] && + LOpt.constraint_topology_multi_commodity_flow(lom) LOpt.constraint_lazycallback_wrapper(lom) - + return end @@ -39,7 +37,7 @@ function objective_LOModel(lom::LaplacianOptModel) return end -function optimize_LOModel!(lom::LaplacianOptModel; optimizer=nothing) +function optimize_LOModel!(lom::LaplacianOptModel; optimizer = nothing) if lom.data["relax_integrality"] JuMP.relax_integrality(lom.model) end @@ -50,72 +48,87 @@ function optimize_LOModel!(lom::LaplacianOptModel; optimizer=nothing) if lom.model.moi_backend.state == MOI.Utilities.NO_OPTIMIZER JuMP.set_optimizer(lom.model, optimizer) else - Memento.warn(_LOGGER, "Model already contains optimizer, cannot use optimizer specified in `optimize_LOModel!`") + Memento.warn( + _LOGGER, + "Model already contains optimizer, cannot use optimizer specified in `optimize_LOModel!`", + ) end end - if JuMP.mode(lom.model) != JuMP.DIRECT && lom.model.moi_backend.state == MOI.Utilities.NO_OPTIMIZER - Memento.error(_LOGGER, "No optimizer specified in `optimize_LOModel!` or the given JuMP model.") + if JuMP.mode(lom.model) != JuMP.DIRECT && + lom.model.moi_backend.state == MOI.Utilities.NO_OPTIMIZER + Memento.error( + _LOGGER, + "No optimizer specified in `optimize_LOModel!` or the given JuMP model.", + ) end - + start_time = time() _, solve_time, solve_bytes_alloc, sec_in_gc = @timed JuMP.optimize!(lom.model) try solve_time = JuMP.solve_time(lom.model) catch - Memento.warn(_LOGGER, "The given optimizer does not provide the SolveTime() attribute, falling back on @timed. This is not a rigorous timing value."); + Memento.warn( + _LOGGER, + "The given optimizer does not provide the SolveTime() attribute, falling back on @timed. This is not a rigorous timing value.", + ) end - + Memento.debug(_LOGGER, "JuMP model optimize time: $(time() - start_time)") - lom.result = LOpt.build_LOModel_result(lom, solve_time) + lom.result = LOpt.build_LOModel_result(lom, solve_time) return lom.result end -function run_LOpt(params::Dict{String, Any}, - lom_optimizer::MOI.OptimizerWithAttributes; - visualize_solution = false, - visualizing_tool = "tikz", - display_edge_weights = false) - - data = LOpt.get_data(params) - model_lopt = LOpt.build_LOModel(data) +function run_LOpt( + params::Dict{String,Any}, + lom_optimizer::MOI.OptimizerWithAttributes; + visualize_solution = false, + visualizing_tool = "tikz", + display_edge_weights = false, +) + data = LOpt.get_data(params) + model_lopt = LOpt.build_LOModel(data) result_lopt = LOpt.optimize_LOModel!(model_lopt, optimizer = lom_optimizer) if visualize_solution - LOpt.visualize_solution(result_lopt, - data, - visualizing_tool = visualizing_tool, - display_edge_weights = display_edge_weights) + LOpt.visualize_solution( + result_lopt, + data, + visualizing_tool = visualizing_tool, + display_edge_weights = display_edge_weights, + ) end return result_lopt end -function run_MaxSpanTree(params::Dict{String, Any}, - lom_optimizer::MOI.OptimizerWithAttributes; - visualize_solution = false, - visualizing_tool = "tikz", - display_edge_weights = false, - lazy_callback = false) - - data = LOpt.get_data(params) - model_mst = LOpt.build_MaxSpanTree_model(data, lazy_callback) +function run_MaxSpanTree( + params::Dict{String,Any}, + lom_optimizer::MOI.OptimizerWithAttributes; + visualize_solution = false, + visualizing_tool = "tikz", + display_edge_weights = false, + lazy_callback = false, +) + data = LOpt.get_data(params) + model_mst = LOpt.build_MaxSpanTree_model(data, lazy_callback) result_mst = LOpt.optimize_LOModel!(model_mst, optimizer = lom_optimizer) if visualize_solution - LOpt.visualize_solution(result_mst, - data, - visualizing_tool = visualizing_tool, - display_edge_weights = display_edge_weights) + LOpt.visualize_solution( + result_mst, + data, + visualizing_tool = visualizing_tool, + display_edge_weights = display_edge_weights, + ) end return result_mst end -function build_MaxSpanTree_model(data::Dict{String, Any}, lazy_callback::Bool) - +function build_MaxSpanTree_model(data::Dict{String,Any}, lazy_callback::Bool) m_mst = LaplacianOptModel(data) LOpt.variable_MaxSpanTree_model(m_mst, lazy_callback) LOpt.constraint_MaxSpanTree_model(m_mst, lazy_callback) @@ -125,16 +138,14 @@ function build_MaxSpanTree_model(data::Dict{String, Any}, lazy_callback::Bool) end function objective_MaxSpanTree_model(lom::LaplacianOptModel) - LOpt.objective_maximize_spanning_tree_cost(lom) return end function variable_MaxSpanTree_model(lom::LaplacianOptModel, lazy_callback::Bool) - LOpt.variable_edge_onoff(lom) - + if !lazy_callback LOpt.variable_multi_commodity_flow(lom) end @@ -143,15 +154,14 @@ function variable_MaxSpanTree_model(lom::LaplacianOptModel, lazy_callback::Bool) end function constraint_MaxSpanTree_model(lom::LaplacianOptModel, lazy_callback::Bool) - LOpt.constraint_single_vertex_cutset(lom) LOpt.constraint_augment_edges_budget(lom) - + if lazy_callback LOpt.constraint_lazycallback_wrapper(lom, max_span_tree = true) else LOpt.constraint_topology_multi_commodity_flow(lom) end - + return end diff --git a/src/objective.jl b/src/objective.jl index 3f8a2e9..af3d77a 100644 --- a/src/objective.jl +++ b/src/objective.jl @@ -7,9 +7,15 @@ function objective_maximize_algebraic_connectivity(lom::LaplacianOptModel) end function objective_maximize_spanning_tree_cost(lom::LaplacianOptModel) - num_nodes = lom.data["num_nodes"] + num_nodes = lom.data["num_nodes"] adjacency_augment_graph = lom.data["adjacency_augment_graph"] - JuMP.@objective(lom.model, Max, sum(adjacency_augment_graph[i,j] * lom.variables[:z_var][i,j] - for i=1:(num_nodes-1), j=(i+1):num_nodes if adjacency_augment_graph[i,j] > 0)) -end \ No newline at end of file + JuMP.@objective( + lom.model, + Max, + sum( + adjacency_augment_graph[i, j] * lom.variables[:z_var][i, j] for + i = 1:(num_nodes-1), j = (i+1):num_nodes if adjacency_augment_graph[i, j] > 0 + ) + ) +end diff --git a/src/relaxations.jl b/src/relaxations.jl index 8ef2ed7..bf213bf 100644 --- a/src/relaxations.jl +++ b/src/relaxations.jl @@ -5,16 +5,15 @@ Computes the valid domain of a given JuMP variable taking into account bounds and the varaible's implicit bounds (e.g. binary). """ function variable_domain(var::JuMP.VariableRef) + lb = -Inf + (JuMP.has_lower_bound(var)) && (lb = JuMP.lower_bound(var)) + (JuMP.is_binary(var)) && (lb = max(lb, 0.0)) - lb = -Inf - (JuMP.has_lower_bound(var)) && (lb = JuMP.lower_bound(var)) - (JuMP.is_binary(var)) && (lb = max(lb, 0.0)) + ub = Inf + (JuMP.has_upper_bound(var)) && (ub = JuMP.upper_bound(var)) + (JuMP.is_binary(var)) && (ub = min(ub, 1.0)) - ub = Inf - (JuMP.has_upper_bound(var)) && (ub = JuMP.upper_bound(var)) - (JuMP.is_binary(var)) && (ub = min(ub, 1.0)) - - return (lower_bound=lb, upper_bound=ub) + return (lower_bound = lb, upper_bound = ub) end """ @@ -29,42 +28,38 @@ z <= JuMP.lower_bound(x)*y + JuMP.upper_bound(y)*x - JuMP.lower_bound(x)*JuMP.up z <= JuMP.upper_bound(x)*y + JuMP.lower_bound(y)*x - JuMP.upper_bound(x)*JuMP.lower_bound(y) ``` """ -function relaxation_bilinear(m::JuMP.Model, - xy::JuMP.VariableRef, - x::JuMP.VariableRef, - y::JuMP.VariableRef) - +function relaxation_bilinear( + m::JuMP.Model, + xy::JuMP.VariableRef, + x::JuMP.VariableRef, + y::JuMP.VariableRef, +) lb_x, ub_x = LOpt.variable_domain(x) lb_y, ub_y = LOpt.variable_domain(y) - + if (lb_x == 0) && (ub_x == 1) && ((lb_y != 0) || (ub_y != 1)) - - JuMP.@constraint(m, xy >= lb_y*x) - JuMP.@constraint(m, xy >= y + ub_y*x - ub_y) - JuMP.@constraint(m, xy <= ub_y*x) - JuMP.@constraint(m, xy <= y + lb_y*x - lb_y) + JuMP.@constraint(m, xy >= lb_y * x) + JuMP.@constraint(m, xy >= y + ub_y * x - ub_y) + JuMP.@constraint(m, xy <= ub_y * x) + JuMP.@constraint(m, xy <= y + lb_y * x - lb_y) elseif (lb_y == 0) && (ub_y == 1) && ((lb_x != 0) || (ub_x != 1)) - - JuMP.@constraint(m, xy >= lb_x*y) - JuMP.@constraint(m, xy >= ub_x*y + x - ub_x) - JuMP.@constraint(m, xy <= lb_x*y + x - lb_x) - JuMP.@constraint(m, xy <= ub_x*y) + JuMP.@constraint(m, xy >= lb_x * y) + JuMP.@constraint(m, xy >= ub_x * y + x - ub_x) + JuMP.@constraint(m, xy <= lb_x * y + x - lb_x) + JuMP.@constraint(m, xy <= ub_x * y) elseif (lb_x == 0) && (ub_x == 1) && (lb_y == 0) && (ub_y == 1) - - JuMP.@constraint(m, xy <= x) - JuMP.@constraint(m, xy <= y) - JuMP.@constraint(m, xy >= x + y - 1) + JuMP.@constraint(m, xy <= x) + JuMP.@constraint(m, xy <= y) + JuMP.@constraint(m, xy >= x + y - 1) else - - JuMP.@constraint(m, xy >= lb_x*y + lb_y*x - lb_x*lb_y) - JuMP.@constraint(m, xy >= ub_x*y + ub_y*x - ub_x*ub_y) - JuMP.@constraint(m, xy <= lb_x*y + ub_y*x - lb_x*ub_y) - JuMP.@constraint(m, xy <= ub_x*y + lb_y*x - ub_x*lb_y) - + JuMP.@constraint(m, xy >= lb_x * y + lb_y * x - lb_x * lb_y) + JuMP.@constraint(m, xy >= ub_x * y + ub_y * x - ub_x * ub_y) + JuMP.@constraint(m, xy <= lb_x * y + ub_y * x - lb_x * ub_y) + JuMP.@constraint(m, xy <= ub_x * y + lb_y * x - ub_x * lb_y) end return m - end \ No newline at end of file +end diff --git a/src/solution.jl b/src/solution.jl index b1ff91a..664c312 100644 --- a/src/solution.jl +++ b/src/solution.jl @@ -4,7 +4,10 @@ function build_LOModel_result(lom::LaplacianOptModel, solve_time::Number) try result_count = JuMP.result_count(lom.model) catch - Memento.warn(_LOGGER, "the given optimizer does not provide the ResultCount() attribute, assuming the solver returned a solution which may be incorrect."); + Memento.warn( + _LOGGER, + "the given optimizer does not provide the ResultCount() attribute, assuming the solver returned a solution which may be incorrect.", + ) end solution = Dict{String,Any}() @@ -12,18 +15,21 @@ function build_LOModel_result(lom::LaplacianOptModel, solve_time::Number) if result_count > 0 solution = LOpt.build_LOModel_solution(lom) else - Memento.warn(_LOGGER, "LaplacianOpt model has no results - solution cannot be built") + Memento.warn( + _LOGGER, + "LaplacianOpt model has no results - solution cannot be built", + ) end result = Dict{String,Any}( - "optimizer" => JuMP.solver_name(lom.model), - "termination_status" => JuMP.termination_status(lom.model), - "primal_status" => JuMP.primal_status(lom.model), - "objective" => get_objective_value(lom.model), - "objective_ub" => get_objective_bound(lom.model), - "solve_time" => solve_time, - "solution" => solution, - "adjacency_base_graph" => lom.data["adjacency_base_graph"], + "optimizer" => JuMP.solver_name(lom.model), + "termination_status" => JuMP.termination_status(lom.model), + "primal_status" => JuMP.primal_status(lom.model), + "objective" => get_objective_value(lom.model), + "objective_ub" => get_objective_bound(lom.model), + "solve_time" => solve_time, + "solution" => solution, + "adjacency_base_graph" => lom.data["adjacency_base_graph"], "adjacency_augment_graph" => lom.data["adjacency_augment_graph"], ) @@ -32,22 +38,22 @@ end "" function get_objective_value(model::JuMP.Model) - obj_val = NaN try obj_val = JuMP.objective_value(model) catch - Memento.warn(_LOGGER, "Objective value is unbounded. Problem may be infeasible or not constrained properly"); + Memento.warn( + _LOGGER, + "Objective value is unbounded. Problem may be infeasible or not constrained properly", + ) end return obj_val end - "" function get_objective_bound(model::JuMP.Model) - obj_lb = -Inf try @@ -59,12 +65,11 @@ function get_objective_bound(model::JuMP.Model) end function build_LOModel_solution(lom::LaplacianOptModel) - solution = Dict{String,Any}() - + for i in keys(lom.variables) solution[String(i)] = JuMP.value.(lom.variables[i]) end return solution -end \ No newline at end of file +end diff --git a/src/types.jl b/src/types.jl index 28d7d1d..74a0c2f 100644 --- a/src/types.jl +++ b/src/types.jl @@ -6,24 +6,22 @@ export LaplacianOptModel The composite mutable struct, `LaplacianOptModel`, holds dictionaries for input data, abstract JuMP model for optimization, variable references and result from solving the JuMP model. """ -mutable struct LaplacianOptModel - data :: Dict{String,Any} - model :: JuMP.Model - variables :: Dict{Symbol,Any} - result :: Dict{String,Any} +mutable struct LaplacianOptModel + data::Dict{String,Any} + model::JuMP.Model + variables::Dict{Symbol,Any} + result::Dict{String,Any} "Contructor for struct `LaplacianOptModel`" function LaplacianOptModel(data::Dict{String,Any}) - data = data - model = JuMP.Model() + model = JuMP.Model() variables = Dict{Symbol,Any}() result = Dict{String,Any}() lom = new(data, model, variables, result) return lom end - end """ @@ -33,11 +31,11 @@ The composite mutable struct, `GraphData`, holds matrices of adjacency matrix, l fiedler vector and algebraic connectivity. """ mutable struct GraphData - adjacency ::Array{<:Number} - laplacian ::Array{<:Number} - fiedler ::Vector{<:Number} - ac ::Number - + adjacency::Array{<:Number} + laplacian::Array{<:Number} + fiedler::Vector{<:Number} + ac::Number + function GraphData(adjacency::Array{<:Number}) adj_matrix = adjacency laplacian = LOpt.laplacian_matrix(adjacency) @@ -47,6 +45,5 @@ mutable struct GraphData graph_data = new(adj_matrix, laplacian, fiedler, ac) return graph_data - end end diff --git a/src/utility.jl b/src/utility.jl index 25cddfe..4042b43 100644 --- a/src/utility.jl +++ b/src/utility.jl @@ -4,19 +4,18 @@ Given a vector of numbers, this function updates the input vector by rounding the values closest to 0, 1 and -1. """ function get_rounded_zeros_and_ones!(v::Array{<:Number}; tol = 1E-6) - for i in findall(abs.(v) .<= tol) v[i] = 0 end - for i in findall(((1-tol) .<= abs.(v) .<= (1+tol))) - if v[i] > 0 + for i in findall(((1 - tol) .<= abs.(v) .<= (1 + tol))) + if v[i] > 0 v[i] = 1 - elseif v[i] < 0 + elseif v[i] < 0 v[i] = -1 end end -end +end """ optimal_graph_edges(adjacency_matrix::Array{<:Number}) @@ -25,25 +24,24 @@ Returns a vector of tuples of edges corresponding to an input adjacency matrix of the graph. """ function optimal_graph_edges(adjacency_matrix::Array{<:Number}) + edges = Vector{Tuple{Int64,Int64}}() - edges = Vector{Tuple{Int64, Int64}}() - num_nodes = size(adjacency_matrix)[1] if !LA.issymmetric(adjacency_matrix) Memento.error(_LOGGER, "Input adjacency matrix is asymmetric") end - for i=1:num_nodes - if !isapprox(adjacency_matrix[i,i], 0) + for i in 1:num_nodes + if !isapprox(adjacency_matrix[i, i], 0) Memento.error(_LOGGER, "Input adjacency matrix cannot have self loops") end end - for i=1:(num_nodes-1) - for j=(i+1):num_nodes - if !isapprox(adjacency_matrix[i,j], 0) - push!(edges, (i,j)) + for i in 1:(num_nodes-1) + for j in (i+1):num_nodes + if !isapprox(adjacency_matrix[i, j], 0) + push!(edges, (i, j)) end end end @@ -58,7 +56,6 @@ Given a weighted adjacency matrix as an input, this function returns the weighted Laplacian matrix of the graph. """ function laplacian_matrix(adjacency_matrix::Array{<:Number}) - if !LA.issymmetric(adjacency_matrix) Memento.error(_LOGGER, "Input adjacency matrix is asymmetric") end @@ -67,26 +64,30 @@ function laplacian_matrix(adjacency_matrix::Array{<:Number}) laplacian_matrix = zeros(Float64, num_nodes, num_nodes) - for i = 1:num_nodes - for j = i:num_nodes - - if i == j - if !isapprox(adjacency_matrix[i,j], 0) - Memento.error(_LOGGER, "Input adjacency matrix cannot have self loops") + for i in 1:num_nodes + for j in i:num_nodes + if i == j + if !isapprox(adjacency_matrix[i, j], 0) + Memento.error( + _LOGGER, + "Input adjacency matrix cannot have self loops", + ) end - laplacian_matrix[i,j] = sum(adjacency_matrix[i,:]) - else - if adjacency_matrix[i,j] <= -1E-6 - Memento.error(_LOGGER, "Input adjacency matrix cannot have negative weights") + laplacian_matrix[i, j] = sum(adjacency_matrix[i, :]) + else + if adjacency_matrix[i, j] <= -1E-6 + Memento.error( + _LOGGER, + "Input adjacency matrix cannot have negative weights", + ) end - if !isapprox(adjacency_matrix[i,j], 0) - laplacian_matrix[i,j] = -adjacency_matrix[i,j] - laplacian_matrix[j,i] = -adjacency_matrix[i,j] + if !isapprox(adjacency_matrix[i, j], 0) + laplacian_matrix[i, j] = -adjacency_matrix[i, j] + laplacian_matrix[j, i] = -adjacency_matrix[i, j] end end - end end @@ -100,18 +101,16 @@ Returns the Fiedler vector or the eigenvector corresponding to the second smallest eigenvalue of the Laplacian matrix for an input weighted adjacency matrix of the graph. """ function fiedler_vector(adjacency_matrix::Array{<:Number}) - L_mat = LOpt.laplacian_matrix(adjacency_matrix) - ac = LOpt.algebraic_connectivity(adjacency_matrix) + ac = LOpt.algebraic_connectivity(adjacency_matrix) - fiedler = LA.eigvecs(L_mat)[:,2] + fiedler = LA.eigvecs(L_mat)[:, 2] - if !isapprox(fiedler' * L_mat * fiedler, ac, atol=1E-6) + if !isapprox(fiedler' * L_mat * fiedler, ac, atol = 1E-6) Memento.error(_LOGGER, "Evaluated eigenvector is not the Fiedler vector") end - - return fiedler + return fiedler end """ @@ -121,11 +120,10 @@ Returns the algebraic connectivity or the second smallest eigenvalue of the Laplacian matrix, for an input weighted adjacency matrix of the graph. """ function algebraic_connectivity(adjacency_matrix::Array{<:Number}) - L_mat = LOpt.laplacian_matrix(adjacency_matrix) - ac = sort(LA.eigvals(L_mat))[2] - + ac = sort(LA.eigvals(L_mat))[2] + return ac end @@ -135,15 +133,16 @@ end Given from and to sets of vertices of connected components of a graph, this function returns a boolean if the input `num_edges` number of candidate edge exist to augment or not, between these two sets of vertices. """ -function _is_flow_cut_valid(cutset_f::Vector{Int64}, - cutset_t::Vector{Int64}, - num_edges::Int64, - adjacency::Array{<:Number}) - - k = 0 +function _is_flow_cut_valid( + cutset_f::Vector{Int64}, + cutset_t::Vector{Int64}, + num_edges::Int64, + adjacency::Array{<:Number}, +) + k = 0 for i in cutset_f for j in cutset_t - if !(isapprox(adjacency[i,j], 0, atol = 1E-6)) + if !(isapprox(adjacency[i, j], 0, atol = 1E-6)) k += 1 if k == num_edges return true @@ -169,15 +168,17 @@ function _violated_eigen_vector(W::Array{<:Number}; tol = 1E-6) end if LA.eigmin(W) <= -tol - if W_eigvals[1] == LA.eigmin(W) - violated_eigen_vec = LA.eigvecs(W)[:,1] + violated_eigen_vec = LA.eigvecs(W)[:, 1] LOpt.get_rounded_zeros_and_ones!(violated_eigen_vec) return violated_eigen_vec else - Memento.warn(_LOGGER, "Eigen cut corresponding to the negative eigenvalue could not be evaluated") + Memento.warn( + _LOGGER, + "Eigen cut corresponding to the negative eigenvalue could not be evaluated", + ) end else return end -end \ No newline at end of file +end diff --git a/src/variables.jl b/src/variables.jl index cb34610..b4545c8 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -4,12 +4,12 @@ function variable_lifted_W_matrix(lom::LaplacianOptModel) num_nodes = lom.data["num_nodes"] - adjacency_base_graph = lom.data["adjacency_base_graph"] - lom.variables[:W_var] = JuMP.@variable(lom.model, W_var[1:num_nodes, 1:num_nodes], Symmetric) + lom.variables[:W_var] = + JuMP.@variable(lom.model, W_var[1:num_nodes, 1:num_nodes], Symmetric) - for i = 1:num_nodes - JuMP.set_lower_bound(W_var[i,i], 0) + for i in 1:num_nodes + JuMP.set_lower_bound(W_var[i, i], 0) end return @@ -20,25 +20,26 @@ function variable_edge_onoff(lom::LaplacianOptModel) adjacency_base_graph = lom.data["adjacency_base_graph"] adjacency_augment_graph = lom.data["adjacency_augment_graph"] - lom.variables[:z_var] = JuMP.@variable(lom.model, z_var[1:num_nodes, 1:num_nodes], Bin, Symmetric) - - for i=1:num_nodes - for j=i:num_nodes + lom.variables[:z_var] = + JuMP.@variable(lom.model, z_var[1:num_nodes, 1:num_nodes], Bin, Symmetric) + + for i in 1:num_nodes + for j in i:num_nodes # No self-loop at nodes if i == j - JuMP.fix(z_var[i,i], 0) + JuMP.fix(z_var[i, i], 0) continue end - if isapprox(adjacency_augment_graph[i,j], 0, atol = 1E-6) + if isapprox(adjacency_augment_graph[i, j], 0, atol = 1E-6) # (i,j) in base graph - if !isapprox(adjacency_base_graph[i,j], 0, atol = 1E-6) - JuMP.fix(z_var[i,j], 1) - JuMP.fix(z_var[j,i], 1) + if !isapprox(adjacency_base_graph[i, j], 0, atol = 1E-6) + JuMP.fix(z_var[i, j], 1) + JuMP.fix(z_var[j, i], 1) else - # (i,j) neither in base not augment graph - JuMP.fix(z_var[i,j], 0) - JuMP.fix(z_var[j,i], 0) + # (i,j) neither in base not augment graph + JuMP.fix(z_var[i, j], 0) + JuMP.fix(z_var[j, i], 0) end end end @@ -48,7 +49,6 @@ function variable_edge_onoff(lom::LaplacianOptModel) end function variable_algebraic_connectivity(lom::LaplacianOptModel) - lom.variables[:γ_var] = JuMP.@variable(lom.model, γ_var >= 0) return @@ -57,7 +57,10 @@ end function variable_multi_commodity_flow(lom::LaplacianOptModel) num_nodes = lom.data["num_nodes"] - lom.variables[:flow_var] = JuMP.@variable(lom.model, 0 <= flow_var[1:num_nodes, 1:num_nodes, 1:(num_nodes-1)] <= 1) + lom.variables[:flow_var] = JuMP.@variable( + lom.model, + 0 <= flow_var[1:num_nodes, 1:num_nodes, 1:(num_nodes-1)] <= 1 + ) return -end \ No newline at end of file +end diff --git a/test/lo_model_tests.jl b/test/lo_model_tests.jl index 81959d7..8e09818 100644 --- a/test/lo_model_tests.jl +++ b/test/lo_model_tests.jl @@ -1,45 +1,48 @@ @testset "Algebraic Connectivity: Optimal solution tests" begin - num_nodes = 5 - instance = 1 - file_path = joinpath(@__DIR__, "..","examples/instances/$(num_nodes)_nodes/$(num_nodes)_$(instance).json") + instance = 1 + file_path = joinpath( + @__DIR__, + "..", + "examples/instances/$(num_nodes)_nodes/$(num_nodes)_$(instance).json", + ) data_dict = LOpt.parse_file(file_path) - params = Dict{String, Any}( + params = Dict{String,Any}( "data_dict" => data_dict, - "augment_budget" => (num_nodes-1), + "augment_budget" => (num_nodes - 1), "solution_type" => "exact", "presolve" => true, "optimizer_log" => true, - "relax_integrality" => false + "relax_integrality" => false, ) result_lo = LaplacianOpt.run_LOpt(params, glpk_optimizer) @test result_lo["termination_status"] == MOI.OPTIMAL @test result_lo["primal_status"] == MOI.FEASIBLE_POINT - @test isapprox(result_lo["objective"], 11.1592, atol=1E-4) - @test isapprox(result_lo["solution"]["z_var"][1,4], 1.0) - @test isapprox(result_lo["solution"]["z_var"][4,1], 1.0) - @test isapprox(result_lo["solution"]["z_var"][2,4], 1.0) - @test isapprox(result_lo["solution"]["z_var"][4,2], 1.0) - @test isapprox(result_lo["solution"]["z_var"][3,4], 1.0) - @test isapprox(result_lo["solution"]["z_var"][4,3], 1.0) - @test isapprox(result_lo["solution"]["z_var"][3,5], 1.0) - @test isapprox(result_lo["solution"]["z_var"][5,3], 1.0) + @test isapprox(result_lo["objective"], 11.1592, atol = 1E-4) + @test isapprox(result_lo["solution"]["z_var"][1, 4], 1.0) + @test isapprox(result_lo["solution"]["z_var"][4, 1], 1.0) + @test isapprox(result_lo["solution"]["z_var"][2, 4], 1.0) + @test isapprox(result_lo["solution"]["z_var"][4, 2], 1.0) + @test isapprox(result_lo["solution"]["z_var"][3, 4], 1.0) + @test isapprox(result_lo["solution"]["z_var"][4, 3], 1.0) + @test isapprox(result_lo["solution"]["z_var"][3, 5], 1.0) + @test isapprox(result_lo["solution"]["z_var"][5, 3], 1.0) @test LA.issymmetric(result_lo["solution"]["W_var"]) - + L_val = result_lo["solution"]["W_var"] - for i=1:num_nodes - L_val[i,i] += result_lo["solution"]["γ_var"] * (num_nodes - 1)/num_nodes + for i in 1:num_nodes + L_val[i, i] += result_lo["solution"]["γ_var"] * (num_nodes - 1) / num_nodes end - for i=1:(num_nodes - 1) - for j=(i+1):num_nodes - L_val[i,j] -= result_lo["solution"]["γ_var"] / num_nodes - L_val[j,i] = L_val[i,j] + for i in 1:(num_nodes-1) + for j in (i+1):num_nodes + L_val[i, j] -= result_lo["solution"]["γ_var"] / num_nodes + L_val[j, i] = L_val[i, j] end end @@ -48,181 +51,201 @@ end @testset "Max Span Tree: Optimal solution tests" begin - num_nodes = 5 - instance = 1 - file_path = joinpath(@__DIR__, "..","examples/instances/$(num_nodes)_nodes/$(num_nodes)_$(instance).json") + instance = 1 + file_path = joinpath( + @__DIR__, + "..", + "examples/instances/$(num_nodes)_nodes/$(num_nodes)_$(instance).json", + ) data_dict = LOpt.parse_file(file_path) - params = Dict{String, Any}( - "data_dict" => data_dict, - "augment_budget" => (num_nodes-1) - ) + params = + Dict{String,Any}("data_dict" => data_dict, "augment_budget" => (num_nodes - 1)) result_mst = LaplacianOpt.run_MaxSpanTree(params, glpk_optimizer) @test result_mst["termination_status"] == MOI.OPTIMAL @test result_mst["primal_status"] == MOI.FEASIBLE_POINT - @test isapprox(result_mst["objective"], 113.7588, atol=1E-4) - @test isapprox(result_mst["solution"]["z_var"][1,2], 1.0) - @test isapprox(result_mst["solution"]["z_var"][2,1], 1.0) - @test isapprox(result_mst["solution"]["z_var"][2,5], 1.0) - @test isapprox(result_mst["solution"]["z_var"][5,2], 1.0) - @test isapprox(result_mst["solution"]["z_var"][3,4], 1.0) - @test isapprox(result_mst["solution"]["z_var"][4,3], 1.0) - @test isapprox(result_mst["solution"]["z_var"][3,5], 1.0) - @test isapprox(result_mst["solution"]["z_var"][5,3], 1.0) - + @test isapprox(result_mst["objective"], 113.7588, atol = 1E-4) + @test isapprox(result_mst["solution"]["z_var"][1, 2], 1.0) + @test isapprox(result_mst["solution"]["z_var"][2, 1], 1.0) + @test isapprox(result_mst["solution"]["z_var"][2, 5], 1.0) + @test isapprox(result_mst["solution"]["z_var"][5, 2], 1.0) + @test isapprox(result_mst["solution"]["z_var"][3, 4], 1.0) + @test isapprox(result_mst["solution"]["z_var"][4, 3], 1.0) + @test isapprox(result_mst["solution"]["z_var"][3, 5], 1.0) + @test isapprox(result_mst["solution"]["z_var"][5, 3], 1.0) end @testset "Max Span Tree: Lazy callback tests" begin - num_nodes = 15 - instance = 1 - file_path = joinpath(@__DIR__, "..","examples/instances/$(num_nodes)_nodes/$(num_nodes)_$(instance).json") + instance = 1 + file_path = joinpath( + @__DIR__, + "..", + "examples/instances/$(num_nodes)_nodes/$(num_nodes)_$(instance).json", + ) data_dict = LOpt.parse_file(file_path) - params = Dict{String, Any}( - "data_dict" => data_dict, - "augment_budget" => (num_nodes-1) - ) + params = + Dict{String,Any}("data_dict" => data_dict, "augment_budget" => (num_nodes - 1)) + + result_mst = + LaplacianOpt.run_MaxSpanTree(params, glpk_optimizer, lazy_callback = true) - result_mst = LaplacianOpt.run_MaxSpanTree(params, glpk_optimizer, lazy_callback=true) - @test result_mst["termination_status"] == MOI.OPTIMAL @test result_mst["primal_status"] == MOI.FEASIBLE_POINT - @test isapprox(result_mst["objective"], 3625.89701005, atol=1E-6) - @test isapprox(result_mst["solution"]["z_var"][1,7], 1.0) - @test isapprox(result_mst["solution"]["z_var"][2,4], 1.0) - @test isapprox(result_mst["solution"]["z_var"][3,4], 1.0) - @test isapprox(result_mst["solution"]["z_var"][4,14], 1.0) - @test isapprox(result_mst["solution"]["z_var"][4,15], 1.0) - @test isapprox(result_mst["solution"]["z_var"][5,14], 1.0) - @test isapprox(result_mst["solution"]["z_var"][6,15], 1.0) - @test isapprox(result_mst["solution"]["z_var"][7,12], 1.0) - @test isapprox(result_mst["solution"]["z_var"][7,15], 1.0) - @test isapprox(result_mst["solution"]["z_var"][8,13], 1.0) - @test isapprox(result_mst["solution"]["z_var"][8,14], 1.0) - @test isapprox(result_mst["solution"]["z_var"][9,10], 1.0) - @test isapprox(result_mst["solution"]["z_var"][9,13], 1.0) - @test isapprox(result_mst["solution"]["z_var"][11,12], 1.0) + @test isapprox(result_mst["objective"], 3625.89701005, atol = 1E-6) + @test isapprox(result_mst["solution"]["z_var"][1, 7], 1.0) + @test isapprox(result_mst["solution"]["z_var"][2, 4], 1.0) + @test isapprox(result_mst["solution"]["z_var"][3, 4], 1.0) + @test isapprox(result_mst["solution"]["z_var"][4, 14], 1.0) + @test isapprox(result_mst["solution"]["z_var"][4, 15], 1.0) + @test isapprox(result_mst["solution"]["z_var"][5, 14], 1.0) + @test isapprox(result_mst["solution"]["z_var"][6, 15], 1.0) + @test isapprox(result_mst["solution"]["z_var"][7, 12], 1.0) + @test isapprox(result_mst["solution"]["z_var"][7, 15], 1.0) + @test isapprox(result_mst["solution"]["z_var"][8, 13], 1.0) + @test isapprox(result_mst["solution"]["z_var"][8, 14], 1.0) + @test isapprox(result_mst["solution"]["z_var"][9, 10], 1.0) + @test isapprox(result_mst["solution"]["z_var"][9, 13], 1.0) + @test isapprox(result_mst["solution"]["z_var"][11, 12], 1.0) # Test multi commodity flow - result_mst_1 = LaplacianOpt.run_MaxSpanTree(params, glpk_optimizer, lazy_callback=false) + result_mst_1 = + LaplacianOpt.run_MaxSpanTree(params, glpk_optimizer, lazy_callback = false) @test result_mst_1["termination_status"] == MOI.OPTIMAL @test result_mst_1["primal_status"] == MOI.FEASIBLE_POINT - @test isapprox(result_mst_1["objective"], 3625.89701005, atol=1E-6) - + @test isapprox(result_mst_1["objective"], 3625.89701005, atol = 1E-6) end @testset "SOC relaxations - 1: constraint_soc_cuts_on_2minors" begin - num_nodes = 5 - instance = 5 - file_path = joinpath(@__DIR__, "..","examples/instances/$(num_nodes)_nodes/$(num_nodes)_$(instance).json") + instance = 5 + file_path = joinpath( + @__DIR__, + "..", + "examples/instances/$(num_nodes)_nodes/$(num_nodes)_$(instance).json", + ) data_dict = LOpt.parse_file(file_path) - params_1 = Dict{String, Any}( + params_1 = Dict{String,Any}( "data_dict" => data_dict, - "augment_budget" => (num_nodes-1), + "augment_budget" => (num_nodes - 1), "eigen_cuts_full" => false, "soc_linearized_cuts" => true, - "topology_flow_cuts" => true + "topology_flow_cuts" => true, ) result_1 = LaplacianOpt.run_LOpt(params_1, glpk_optimizer) - + @test result_1["termination_status"] == MOI.OPTIMAL @test result_1["primal_status"] == MOI.FEASIBLE_POINT - @test isapprox(result_1["objective"], 15.992792245, atol=1E-6) - @test isapprox(result_1["solution"]["z_var"][1,5], 1.0) - @test isapprox(result_1["solution"]["z_var"][2,5], 1.0) - @test isapprox(result_1["solution"]["z_var"][3,4], 1.0) - @test isapprox(result_1["solution"]["z_var"][3,5], 1.0) - @test isapprox(result_1["solution"]["z_var"][4,3], 1.0) - - params_2 = Dict{String, Any}( + @test isapprox(result_1["objective"], 15.992792245, atol = 1E-6) + @test isapprox(result_1["solution"]["z_var"][1, 5], 1.0) + @test isapprox(result_1["solution"]["z_var"][2, 5], 1.0) + @test isapprox(result_1["solution"]["z_var"][3, 4], 1.0) + @test isapprox(result_1["solution"]["z_var"][3, 5], 1.0) + @test isapprox(result_1["solution"]["z_var"][4, 3], 1.0) + + params_2 = Dict{String,Any}( "data_dict" => data_dict, - "augment_budget" => (num_nodes-1), + "augment_budget" => (num_nodes - 1), "eigen_cuts_full" => true, - "soc_linearized_cuts" => true + "soc_linearized_cuts" => true, ) result_2 = LaplacianOpt.run_LOpt(params_2, glpk_optimizer) - + @test result_2["termination_status"] == MOI.OPTIMAL @test result_2["primal_status"] == MOI.FEASIBLE_POINT @test result_2["objective"] <= result_1["objective"] - @test isapprox(result_2["objective"], 13.065372790, atol=1E-6) - @test isapprox(result_1["solution"]["z_var"], result_2["solution"]["z_var"], atol = 1E-6) + @test isapprox(result_2["objective"], 13.065372790, atol = 1E-6) + @test isapprox( + result_1["solution"]["z_var"], + result_2["solution"]["z_var"], + atol = 1E-6, + ) end @testset "Minor-based eigen relaxations: 2x2 and 3x3" begin - num_nodes = 5 - instance = 3 - file_path = joinpath(@__DIR__, "..","examples/instances/$(num_nodes)_nodes/$(num_nodes)_$(instance).json") + instance = 3 + file_path = joinpath( + @__DIR__, + "..", + "examples/instances/$(num_nodes)_nodes/$(num_nodes)_$(instance).json", + ) data_dict = LOpt.parse_file(file_path) - params_1 = Dict{String, Any}( + params_1 = Dict{String,Any}( "data_dict" => data_dict, - "augment_budget" => (num_nodes-1), + "augment_budget" => (num_nodes - 1), "eigen_cuts_full" => true, - "eigen_cuts_2minors" => true, - "eigen_cuts_3minors" => true, - "topology_multi_commodity" => true + "eigen_cuts_2minors" => true, + "eigen_cuts_3minors" => true, + "topology_multi_commodity" => true, ) result_1 = LaplacianOpt.run_LOpt(params_1, glpk_optimizer) - + @test result_1["termination_status"] == MOI.OPTIMAL @test result_1["primal_status"] == MOI.FEASIBLE_POINT - @test isapprox(result_1["objective"], 8.356739455, atol=1E-6) - @test isapprox(result_1["solution"]["z_var"][1,4], 1.0) - @test isapprox(result_1["solution"]["z_var"][2,4], 1.0) - @test isapprox(result_1["solution"]["z_var"][2,5], 1.0) - @test isapprox(result_1["solution"]["z_var"][3,4], 1.0) + @test isapprox(result_1["objective"], 8.356739455, atol = 1E-6) + @test isapprox(result_1["solution"]["z_var"][1, 4], 1.0) + @test isapprox(result_1["solution"]["z_var"][2, 4], 1.0) + @test isapprox(result_1["solution"]["z_var"][2, 5], 1.0) + @test isapprox(result_1["solution"]["z_var"][3, 4], 1.0) end @testset "Test base graph with existing edges" begin function data_I() - data_dict = Dict{String, Any}() + data_dict = Dict{String,Any}() data_dict["num_nodes"] = 4 - data_dict["adjacency_base_graph"] = [0 2 0 0 - 2 0 3 0 - 0 3 0 4 - 0 0 4 0] - data_dict["adjacency_augment_graph"] = [0 0 4 8 - 0 0 0 7 - 4 0 0 0 - 8 7 0 0] + data_dict["adjacency_base_graph"] = [ + 0 2 0 0 + 2 0 3 0 + 0 3 0 4 + 0 0 4 0 + ] + data_dict["adjacency_augment_graph"] = [ + 0 0 4 8 + 0 0 0 7 + 4 0 0 0 + 8 7 0 0 + ] augment_budget = 2 return data_dict, augment_budget end data_dict, augment_budget = data_I() - params = Dict{String, Any}( + params = Dict{String,Any}( "data_dict" => data_dict, "augment_budget" => augment_budget, "eigen_cuts_full" => true, "soc_linearized_cuts" => true, - "eigen_cuts_2minors" => true, - "eigen_cuts_3minors" => true, + "eigen_cuts_2minors" => true, + "eigen_cuts_3minors" => true, "topology_multi_commodity" => true, ) result = LaplacianOpt.run_LOpt(params, glpk_optimizer) - + @test result["termination_status"] == MOI.OPTIMAL @test result["primal_status"] == MOI.FEASIBLE_POINT @test isapprox(result["objective"], 7.8551986404, atol = 1E-6) - @test isapprox(result["solution"]["z_var"][1,4], 1.0) - @test isapprox(result["solution"]["z_var"][2,4], 1.0) + @test isapprox(result["solution"]["z_var"][1, 4], 1.0) + @test isapprox(result["solution"]["z_var"][2, 4], 1.0) function data_II() num_nodes = 5 - instance = 1 - file_path = joinpath(@__DIR__, "..", "examples/instances/$(num_nodes)_nodes/$(num_nodes)_$(instance)_test.json") + instance = 1 + file_path = joinpath( + @__DIR__, + "..", + "examples/instances/$(num_nodes)_nodes/$(num_nodes)_$(instance)_test.json", + ) data_dict = LOpt.parse_file(file_path) augment_budget = 3 return data_dict, augment_budget @@ -230,56 +253,59 @@ end data_dict, augment_budget = data_II() - params = Dict{String, Any}( - "data_dict" => data_dict, - "augment_budget" => augment_budget - ) + params = + Dict{String,Any}("data_dict" => data_dict, "augment_budget" => augment_budget) result = LaplacianOpt.run_LOpt(params, glpk_optimizer) - + @test result["termination_status"] == MOI.OPTIMAL @test result["primal_status"] == MOI.FEASIBLE_POINT @test isapprox(result["objective"], 34.404789725, atol = 1E-6) - @test isapprox(result["solution"]["z_var"][1,4], 1.0) - @test isapprox(result["solution"]["z_var"][2,4], 1.0) - @test isapprox(result["solution"]["z_var"][2,5], 1.0) + @test isapprox(result["solution"]["z_var"][1, 4], 1.0) + @test isapprox(result["solution"]["z_var"][2, 4], 1.0) + @test isapprox(result["solution"]["z_var"][2, 5], 1.0) end @testset "Test hamiltonian cycle problem" begin function data_II() - data_dict = Dict{String, Any}() + data_dict = Dict{String,Any}() data_dict["num_nodes"] = 5 - data_dict["adjacency_base_graph"] = zeros(5,5) - data_dict["adjacency_augment_graph"] = [0 10 14 12 15; 10 0 9 6 7; 14 9 0 8 9; 12 6 8 0 3; 15 7 9 3 0] + data_dict["adjacency_base_graph"] = zeros(5, 5) + data_dict["adjacency_augment_graph"] = + [0 10 14 12 15; 10 0 9 6 7; 14 9 0 8 9; 12 6 8 0 3; 15 7 9 3 0] augment_budget = 5 return data_dict, augment_budget end data_dict, augment_budget = data_II() - params = Dict{String, Any}( + params = Dict{String,Any}( "data_dict" => data_dict, "augment_budget" => augment_budget, "eigen_cuts_full" => true, - "graph_type" => "hamiltonian_cycle" + "graph_type" => "hamiltonian_cycle", ) result = LaplacianOpt.run_LOpt(params, glpk_optimizer) - + @test result["termination_status"] == MOI.OPTIMAL @test result["primal_status"] == MOI.FEASIBLE_POINT @test isapprox(result["objective"], 11.53910488161, atol = 1E-6) - @test isapprox(result["solution"]["z_var"][1,3], 1.0) - @test isapprox(result["solution"]["z_var"][1,5], 1.0) - @test isapprox(result["solution"]["z_var"][2,4], 1.0) - @test isapprox(result["solution"]["z_var"][2,5], 1.0) - @test isapprox(result["solution"]["z_var"][3,4], 1.0) + @test isapprox(result["solution"]["z_var"][1, 3], 1.0) + @test isapprox(result["solution"]["z_var"][1, 5], 1.0) + @test isapprox(result["solution"]["z_var"][2, 4], 1.0) + @test isapprox(result["solution"]["z_var"][2, 5], 1.0) + @test isapprox(result["solution"]["z_var"][3, 4], 1.0) end @testset "Test topology flow cuts" begin function data_I() num_nodes = 8 - instance = 1 - file_path = joinpath(@__DIR__, "..", "examples/instances/$(num_nodes)_nodes/$(num_nodes)_$(instance).json") + instance = 1 + file_path = joinpath( + @__DIR__, + "..", + "examples/instances/$(num_nodes)_nodes/$(num_nodes)_$(instance).json", + ) data_dict = LOpt.parse_file(file_path) augment_budget = 7 return data_dict, augment_budget @@ -287,16 +313,16 @@ end data_dict, augment_budget = data_I() - params = Dict{String, Any}( - "data_dict" => data_dict, - "augment_budget" => augment_budget, - "eigen_cuts_full" => false, - "eigen_cuts_2minors" => true, - "topology_flow_cuts" => true, - "time_limit" => 1.5, + params = Dict{String,Any}( + "data_dict" => data_dict, + "augment_budget" => augment_budget, + "eigen_cuts_full" => false, + "eigen_cuts_2minors" => true, + "topology_flow_cuts" => true, + "time_limit" => 1.5, ) result = LaplacianOpt.run_LOpt(params, glpk_optimizer) @test result["termination_status"] == MOI.TIME_LIMIT @test result["primal_status"] == MOI.FEASIBLE_POINT -end \ No newline at end of file +end diff --git a/test/log_tests.jl b/test/log_tests.jl index 1f21b64..1c77a0a 100644 --- a/test/log_tests.jl +++ b/test/log_tests.jl @@ -1,31 +1,50 @@ @testset "visualize_solution tests" begin - num_nodes = 5 - instance = 1 - file_path = joinpath(@__DIR__, "..","examples/instances/$(num_nodes)_nodes/$(num_nodes)_$(instance).json") + instance = 1 + file_path = joinpath( + @__DIR__, + "..", + "examples/instances/$(num_nodes)_nodes/$(num_nodes)_$(instance).json", + ) data_dict = LOpt.parse_file(file_path) - params = Dict{String, Any}( + params = Dict{String,Any}( "data_dict" => data_dict, - "augment_budget" => (num_nodes-1), + "augment_budget" => (num_nodes - 1), "tol_zero" => 1E-4, "tol_psd" => 1E-3, "eigen_cuts_full" => true, "topology_flow_cuts" => true, - "lazycuts_logging" => true + "lazycuts_logging" => true, ) data = LOpt.get_data(params) - model_lopt = LOpt.build_LOModel(data) + model_lopt = LOpt.build_LOModel(data) result_lopt = LOpt.optimize_LOModel!(model_lopt, optimizer = glpk_optimizer) - LOpt.visualize_solution(result_lopt, data, visualizing_tool = "tikz", plot_file_format = "tex") - LOpt.visualize_solution(result_lopt, data, visualizing_tool = "tikz", plot_file_format = "tex", display_edge_weights=true) + LOpt.visualize_solution( + result_lopt, + data, + visualizing_tool = "tikz", + plot_file_format = "tex", + ) + LOpt.visualize_solution( + result_lopt, + data, + visualizing_tool = "tikz", + plot_file_format = "tex", + display_edge_weights = true, + ) LOpt.visualize_solution(result_lopt, data, visualizing_tool = "graphviz") - LOpt.visualize_solution(result_lopt, data, visualizing_tool = "graphviz", display_edge_weights=false) + LOpt.visualize_solution( + result_lopt, + data, + visualizing_tool = "graphviz", + display_edge_weights = false, + ) @test result_lopt["termination_status"] == MOI.OPTIMAL @test result_lopt["primal_status"] == MOI.FEASIBLE_POINT # Nothing more to test for coverage since the visualizaiton depends on exterior packages. So, I believe the dot files generated are accurate. -end \ No newline at end of file +end diff --git a/test/relaxations_test.jl b/test/relaxations_test.jl index 151d360..97f628f 100644 --- a/test/relaxations_test.jl +++ b/test/relaxations_test.jl @@ -1,11 +1,11 @@ @testset "relaxation_bilinear tests" begin LOpt.silence() - + m = JuMP.Model(glpk_optimizer) LB = [-1, -2.5, 0, -3, 0] UB = [2.5, 0, 1, 3.2, 1] - JuMP.@variable(m, LB[i] <= x[i=1:5] <= UB[i]) + JuMP.@variable(m, LB[i] <= x[i = 1:5] <= UB[i]) JuMP.@variable(m, z[1:6]) LOpt.relaxation_bilinear(m, z[1], x[1], x[2]) @@ -18,8 +18,8 @@ JuMP.@constraint(m, sum(x) >= 3.3) JuMP.@constraint(m, sum(x) <= 3.8) - coefficients = [0.4852, -0.7795, 0.1788, 0.7778, 0.3418, -0.3407] - JuMP.set_objective_function(m, sum(coefficients[i] * z[i] for i=1:length(z))) + coefficients = [0.4852, -0.7795, 0.1788, 0.7778, 0.3418, -0.3407] + JuMP.set_objective_function(m, sum(coefficients[i] * z[i] for i in 1:length(z))) obj_vals = Vector{Float64}() @@ -29,8 +29,7 @@ @test JuMP.termination_status(m) == MOI.OPTIMAL push!(obj_vals, JuMP.objective_value(m)) end - - @test obj_vals[1] <= -9.922644 # global optimum value - @test obj_vals[2] >= 4.908516 # global optimum value -end \ No newline at end of file + @test obj_vals[1] <= -9.922644 # global optimum value + @test obj_vals[2] >= 4.908516 # global optimum value +end diff --git a/test/runtests.jl b/test/runtests.jl index 8d7be2d..a8d1330 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,8 +8,8 @@ import GLPK import MathOptInterface const LOpt = LaplacianOpt -const LA = LinearAlgebra -const MOI = MathOptInterface +const LA = LinearAlgebra +const MOI = MathOptInterface # Suppress warnings during testing LOpt.logger_config!("error") @@ -17,10 +17,8 @@ LOpt.logger_config!("error") glpk_optimizer = JuMP.optimizer_with_attributes(GLPK.Optimizer, MOI.Silent() => true) @testset "LaplacianOpt" begin - include("utility_tests.jl") include("lo_model_tests.jl") include("relaxations_test.jl") include("log_tests.jl") - -end \ No newline at end of file +end diff --git a/test/utility_tests.jl b/test/utility_tests.jl index de0bb95..0c418f8 100644 --- a/test/utility_tests.jl +++ b/test/utility_tests.jl @@ -3,47 +3,42 @@ A2 = [0 -0.99999; 0.999999 2.000000001; 1.0000001 -1.000001; -0.000001 0.0000001] LOpt.get_rounded_zeros_and_ones!(A1, tol = 1E-5) - @test A1[1,2] == -1 - @test A1[2,1] == 1 - @test A1[2,2] > 2 - @test A1[3,1] == 1 - @test A1[3,2] == -1 - @test A1[4,1] == 0 - @test A1[4,2] == 0 + @test A1[1, 2] == -1 + @test A1[2, 1] == 1 + @test A1[2, 2] > 2 + @test A1[3, 1] == 1 + @test A1[3, 2] == -1 + @test A1[4, 1] == 0 + @test A1[4, 2] == 0 LOpt.get_rounded_zeros_and_ones!(A2, tol = 1E-6) - @test A2[1,2] > -1 - + @test A2[1, 2] > -1 end @testset "utility tests: optimal_graph_edges" begin - z_val = [0.0 1 1; 1 0 0; 1 0 0] edges = LOpt.optimal_graph_edges(z_val) - @test edges[1] == (1,2) - @test edges[2] == (1,3) - + @test edges[1] == (1, 2) + @test edges[2] == (1, 3) end @testset "utility tests: algebraic_connectivity and fiedler_vector" begin - adj_mat_1 = [0 1 2; 1 0 3; 2 3 0.0] graph_data = LOpt.GraphData(adj_mat_1) - @test isapprox(graph_data.ac, 4.267949192, atol=1E-6) + @test isapprox(graph_data.ac, 4.267949192, atol = 1E-6) v = [0.7886751345948128, -0.5773502691896261, -0.21132486540518727] - @test isapprox(v, graph_data.fiedler, atol=1E-6) + @test isapprox(v, graph_data.fiedler, atol = 1E-6) adj_mat_2 = [0 1 2; 1 0 0; 2 0 0.0] graph_data = LOpt.GraphData(adj_mat_2) - @test isapprox(graph_data.ac, 1.26794919243, atol=1E-6) - v = [ -0.2113248654051879, 0.7886751345948118, -0.5773502691896268] - @test isapprox(v, graph_data.fiedler, atol=1E-6) + @test isapprox(graph_data.ac, 1.26794919243, atol = 1E-6) + v = [-0.2113248654051879, 0.7886751345948118, -0.5773502691896268] + @test isapprox(v, graph_data.fiedler, atol = 1E-6) # Disconnected graph adj_mat_3 = [0 1 1 0; 1 0 0 0; 1 0 0 0; 0 0 0 0.0] graph_data = LOpt.GraphData(adj_mat_3) - @test isapprox(graph_data.ac, 0, atol=1E-6) - v = [ 0.5773502691896258, 0.5773502691896273, 0.5773502691896241, 0.0] - @test isapprox(v, graph_data.fiedler, atol=1E-6) - -end \ No newline at end of file + @test isapprox(graph_data.ac, 0, atol = 1E-6) + v = [0.5773502691896258, 0.5773502691896273, 0.5773502691896241, 0.0] + @test isapprox(v, graph_data.fiedler, atol = 1E-6) +end