Merge pull request #41 from harshangrjn/cycle_graph_support

Cycle graph added

CHANGELOG.md

@@ -1,6 +1,14 @@
 LaplacianOpt.jl Change Log
 =========================
 
+### 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`
+- `_is_flow_cut_valid` can handle variable number of edges to be verified in cutset
+- Added support for subtour elimination constraints
+- Added option for `time_limit` in params (default = 10800 s)
+- Clean up in `data.jl` for logging
+
 ### v0.2.1
 - Update in `log.jl` for handling close to zero integral solutions before rounding
 - Update in `lopt_model.jl` to handle displaying edge wts in plotting 

Project.toml

@@ -1,7 +1,7 @@
 name = "LaplacianOpt"
 uuid = "bb20392f-64fb-4001-92e8-14b3aedd5a9e"
 authors = ["Harsha Nagarajan"]
-version = "0.2.1"
+version = "0.3.0"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"

examples/run_examples.jl

@@ -1,7 +1,7 @@
 import LaplacianOpt as LOpt
 using JuMP
-# using Gurobi
 using CPLEX
+# using Gurobi
 # using GLPK
 
 include("optimizer.jl")
@@ -29,7 +29,7 @@ function data_I()
     return data_dict, augment_budget
 end
 
-#=
+#=  
  Option II: Directly input the data dictionary (data_dict) with 
             num_nodes, adjacency_base_graph and adjacency_augment_graph
 =#
@@ -38,7 +38,7 @@ function data_II()
     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]
-    augment_budget = 2
+    augment_budget = 3
     return data_dict, augment_budget
 end
 
@@ -53,7 +53,9 @@ params = Dict{String, Any}(
     "eigen_cuts_full"     => true,
     "soc_linearized_cuts" => false,
     "eigen_cuts_2minors"  => false,
-    "eigen_cuts_3minors"  => false
+    "eigen_cuts_3minors"  => false,
+    "topology_flow_cuts"  => true,
+    # "time_limit"          => 3600,
     )
 
 #----------------------------------------------------------------#
@@ -62,7 +64,6 @@ params = Dict{String, Any}(
 #----------------------------------------------------------------#
 result = LOpt.run_LOpt(params, 
                        lopt_optimizer;
-                       # Make this true to plot the graph solution
-                       visualize_solution   = false,   
+                       visualize_solution   = false,  # Make this true to plot the graph solution
                        visualizing_tool     = "tikz", # "graphviz" is another option
                        display_edge_weights = false)

src/constraints.jl

@@ -15,7 +15,7 @@ function constraint_build_W_var_matrix(lom::LaplacianOptModel)
     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
@@ -25,13 +25,22 @@ function constraint_single_vertex_cutset(lom::LaplacianOptModel)
     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"]
 
     adjacency_full_graph = adjacency_augment_graph
     (num_edges_existing > 0) && (adjacency_full_graph += adjacency_base_graph)
 
-    for i = 1:num_nodes
-        if sum(adjacency_full_graph[i,:]) > 0
-            JuMP.@constraint(lom.model, sum(lom.variables[:z_var][i,j] for j in setdiff((1:num_nodes), i)) >= 1)
+    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)
+            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)
+            end
         end
     end
 
@@ -162,7 +171,7 @@ function constraint_eigen_cuts_on_3minors(W_val::Matrix{<:Number}, cb_cuts, lom:
             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])
-                end     
+                end
             end
         end
     end
@@ -193,27 +202,58 @@ end
 
 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"]
-
-    cc_lazy = Graphs.connected_components(Graphs.SimpleGraph(abs.(z_val)))
+    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)))
 
-    max_cc  = 5 # increase this to any greater integer value and it works
+    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 
+        return
     end
 
-    if length(cc_lazy) in 2:max_cc
-
-        for k = 1:(length(cc_lazy)-1)
+    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))
+            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))
+                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 
+        if graph_type == "hamiltonian_cycle"
+            num_edges_cutset = 2
+        end
+
+        for k = 1:(length(cc_lazy)-1)   
             # From cutset  
-            cutset_f = cc_lazy[k]        
+            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, adjacency_augment_graph)
+            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))) >= 1)
+                                            if !(isapprox(adjacency_augment_graph[i,j], 0, atol=1E-6))) >= num_edges_cutset)
 
                 MOI.submit(lom.model, MOI.LazyConstraint(cb_cuts), con)
             end

src/data.jl

@@ -39,14 +39,25 @@ function get_data(params::Dict{String, Any})
         solution_type = "exact"
     end
 
-    # Relax Integrality 
+    # Relax integrality 
     if "relax_integrality" in keys(params)
         relax_integrality = params["relax_integrality"]
     else
         # default value
         relax_integrality = false
     end
 
+    # Graph type
+    if "graph_type" in keys(params)
+        graph_type = params["graph_type"]
+        if !(graph_type in ["any", "hamiltonian_cycle"])
+            Memento.error(_LOGGER, "Invalid graph type, $graph_type, in input params")
+        end
+    else
+        # default value
+        graph_type = "any"
+    end
+
     # Tolerance to verify zero values
     if "tol_zero" in keys(params)
         tol_zero = params["tol_zero"]
@@ -65,57 +76,71 @@ function get_data(params::Dict{String, Any})
 
     if "eigen_cuts_full" in keys(params)
         eigen_cuts_full = params["eigen_cuts_full"]
-        if eigen_cuts_full
-            Memento.info(_LOGGER, "Applying full-sized eigen cuts")
-        end
     else
         #default value
-        Memento.info(_LOGGER, "Applying full-sized eigen cuts")
         eigen_cuts_full = true
     end
+    if eigen_cuts_full
+        Memento.info(_LOGGER, "Applying full-sized eigen cuts")
+    end
 
     if "soc_linearized_cuts" in keys(params)
         soc_linearized_cuts = params["soc_linearized_cuts"]
-        if soc_linearized_cuts
-            Memento.info(_LOGGER, "Applying linearized SOC cuts (2x2 minors)")
-        end
     else
         #default value
         soc_linearized_cuts = false
     end
+    if soc_linearized_cuts
+        Memento.info(_LOGGER, "Applying linearized SOC cuts (2x2 minors)")
+    end
 
     if "eigen_cuts_2minors" in keys(params)
         eigen_cuts_2minors = params["eigen_cuts_2minors"]
-        if eigen_cuts_2minors
-            Memento.info(_LOGGER, "Applying eigen cuts (2x2 minors)")
-        end
     else
         #default value
         eigen_cuts_2minors = false
     end
+    if eigen_cuts_2minors
+        Memento.info(_LOGGER, "Applying eigen cuts (2x2 minors)")
+    end
 
     if "eigen_cuts_3minors" in keys(params)
         eigen_cuts_3minors = params["eigen_cuts_3minors"]
-        if eigen_cuts_3minors
-            Memento.info(_LOGGER, "Applying eigen cuts (3x3 minors)")
-        end
     else
         #default value
         eigen_cuts_3minors = false
     end
+    if eigen_cuts_3minors
+        Memento.info(_LOGGER, "Applying eigen cuts (3x3 minors)")
+    end
 
+    topology_multi_commodity = false # default
+    if "topology_multi_commodity" in keys(params)
+        topology_multi_commodity = params["topology_multi_commodity"]
+    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")
+    end
+    if topology_multi_commodity
+        Memento.info(_LOGGER, "Applying topology multi commodity constraints")
+    end
+
+    topology_flow_cuts = true # default
     if "topology_flow_cuts" in keys(params)
         topology_flow_cuts = params["topology_flow_cuts"]
-    elseif data_dict["is_base_graph_connected"]
+    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")
+    elseif topology_multi_commodity
         topology_flow_cuts = false
-    else
-        #default value
+    end
+    if topology_flow_cuts
         Memento.info(_LOGGER, "Applying topology flow cuts")
-        topology_flow_cuts = true
     end
 
-    if eigen_cuts_full || topology_flow_cuts || soc_linearized_cuts
+    if eigen_cuts_full || topology_flow_cuts || soc_linearized_cuts || eigen_cuts_2minors || eigen_cuts_3minors
         lazy_callback_status = true
     end
 
@@ -126,24 +151,34 @@ function get_data(params::Dict{String, Any})
         lazycuts_logging = false
     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,
-                             "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,
-                             "lazy_callback_status"    => lazy_callback_status,
-                             "relax_integrality"       => relax_integrality)
+    if "time_limit" in keys(params)
+        time_limit = params["time_limit"]
+    else
+        #default value (in seconds)
+        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)
 
     # Optimizer
     if "optimizer" in keys(params)
@@ -299,4 +334,12 @@ function _detect_infeasbility_in_data(data::Dict{String, Any})
             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")
+        end
+    end
+
 end

src/lopt_model.jl

@@ -16,6 +16,7 @@ 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)
     LOpt.variable_algebraic_connectivity(lom)
 
     return
@@ -26,6 +27,7 @@ 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)
     LOpt.constraint_lazycallback_wrapper(lom)
 
     return
@@ -42,6 +44,8 @@ function optimize_LOModel!(lom::LaplacianOptModel; optimizer=nothing)
         JuMP.relax_integrality(lom.model)
     end
 
+    JuMP.set_time_limit_sec(lom.model, lom.data["time_limit"])
+
     if JuMP.mode(lom.model) != JuMP.DIRECT && optimizer !== nothing
         if lom.model.moi_backend.state == MOI.Utilities.NO_OPTIMIZER
             JuMP.set_optimizer(lom.model, optimizer)