Merge pull request #52 from harshangrjn/cheeger_update

Cheeger update

Author: harshangrjn

CHANGELOG.md

@@ -1,6 +1,10 @@
 LaplacianOpt.jl Change Log
 =========================
 
+### v0.7.1
+- Clean up of `round_zeros_ones!`, `laplacian_matrix`, `_violated_eigen_vector`, `weighted_adjacency_matrix`, `edge_combinations`, `optimal_graph_edges` functions
+- Added support for cheeger constant of graph based on set cardinality (updated MILP formulation) and set volume definitions
+
 ### v0.7.0
 - Solver logging option added in `optimizers.jl`
 - Added option `minors_on_augment_edges` to include principal minor cuts only corresponding to vertices with augmentable edges - helps in reducing run times

Project.toml

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

README.md

@@ -37,13 +37,13 @@ Please report any issues via the Github **[issue tracker](https://github.com/har
 This work was supported by Los Alamos National Laboratory (LANL)'s LDRD Early Career Research Award (20190590ECR) and [LANL-TAMU's collaborative research project](https://nationallabsoffice.tamus.edu/the-texas-am-university-system-and-los-alamos-national-laboratory-partner-to-design-robust-networks/) grant. The primary developer of this package is [Harsha Nagarajan](http://harshanagarajan.com) ([@harshangrjn](https://github.com/harshangrjn)). 
 
 ## Citing LaplacianOpt
-If you find LaplacianOpt.jl useful in your work, we request you to cite the following papers [\[link-1\]](https://doi.org/10.1109/ECC.2015.7330770) [\[link-2\]](https://arxiv.org/abs/2304.08571): 
+If you find LaplacianOpt.jl useful in your work, we request you to cite the following papers [\[link-1\]](https://doi.org/10.1109/ECC.2015.7330770) [\[link-2\]](https://doi.org/10.1109/TCNS.2024.3431408): 
 ```bibtex
 @article{LOpt_TCNS2024,
-  title={Optimal Robust Network Design: Formulations and Algorithms for Maximizing Algebraic Connectivity},
+  title={Optimal robust network design: Formulations and algorithms for maximizing algebraic connectivity},
   author={Somisetty, Neelkamal and Nagarajan, Harsha and Darbha, Swaroop},
-  journal={IEEE Transactions on Control of Network Systems (accepted)},
-  url = {https://arxiv.org/abs/2304.08571},
+  journal={IEEE Transactions on Control of Network Systems},
+  url = {https://doi.org/10.1109/TCNS.2024.3431408},
   year={2024},
   publisher={IEEE}
 }

docs/src/index.md

@@ -34,13 +34,13 @@ Pkg.test("LaplacianOpt")
 ```
 
 ## Citing LaplacianOpt
-If you find LaplacianOpt.jl useful in your work, we request you to cite the following papers [\[link-1\]](https://doi.org/10.1109/ECC.2015.7330770) [\[link-2\]](https://arxiv.org/abs/2304.08571): 
+If you find LaplacianOpt.jl useful in your work, we request you to cite the following papers [\[link-1\]](https://doi.org/10.1109/ECC.2015.7330770) [\[link-2\]](https://doi.org/10.1109/TCNS.2024.3431408): 
 ```bibtex
 @article{LOpt_TCNS2024,
-  title={Optimal Robust Network Design: Formulations and Algorithms for Maximizing Algebraic Connectivity},
+  title={Optimal robust network design: Formulations and algorithms for maximizing algebraic connectivity},
   author={Somisetty, Neelkamal and Nagarajan, Harsha and Darbha, Swaroop},
-  journal={IEEE Transactions on Control of Network Systems (accepted)},
-  url = {https://arxiv.org/abs/2304.08571},
+  journal={IEEE Transactions on Control of Network Systems},
+  url = {https://doi.org/10.1109/TCNS.2024.3431408},
   year={2024},
   publisher={IEEE}
 }

examples/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 CPLEX = "a076750e-1247-5638-91d2-ce28b192dca0"
 GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
+Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 Gurobi = "2e9cd046-0924-5485-92f1-d5272153d98b"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"

examples/optimizers.jl

@@ -1,13 +1,15 @@
 #=====================================#
 # MIP solvers (commercial, but fast)  #
 #=====================================#
-function get_gurobi(; solver_log = true)
+function get_gurobi(; solver_log = true, time_limit = 1E5)
     GRB_ENV = Gurobi.Env()
     return optimizer_with_attributes(
         () -> Gurobi.Optimizer(GRB_ENV),
         MOI.Silent() => !solver_log,
         # "MIPFocus" => 3, # Focus on optimality over feasibility 
         "Presolve" => 1,
+        "TimeLimit" => time_limit,
+        # "NonConvex" => 2
     )
 end
 

examples/run_examples.jl

@@ -1,15 +1,16 @@
 import LaplacianOpt as LOpt
 using JuMP
 using CPLEX
+using Gurobi
 # using GLPK
 
 include("optimizers.jl")
 
-#-----------------------------------#
-#            MIP solver             #
-# (> Cplex 22.1 performs the best)  #
-#-----------------------------------#
-lopt_optimizer = get_cplex(solver_log = true)
+#------------------------------------#
+#            MIP solver              #
+# (Gurobi 11.0.1+ performs the best) #
+#------------------------------------#
+lopt_optimizer = get_gurobi(solver_log = true)
 
 #-------------------------------------#
 #      User-defined input graphs      #

src/constraints.jl

@@ -272,11 +272,19 @@ function constraint_cheeger_cuts(
     lom::LaplacianOptModel,
     optimizer::MOI.OptimizerWithAttributes,
 )
+    ac_lower_bound = floor(lom.options.best_lower_bound, digits = 3)
+
+    if isapprox(ac_lower_bound, 0, atol = 1E-4)
+        Memento.info(
+            _LOGGER,
+            "Cheeger cuts may be ineffective: best λ₂ lower bound is approximately 0.",
+        )
+        return
+    end
+
     result_dict =
         LOpt.cheeger_constant(z_val .* lom.data["adjacency_augment_graph"], optimizer)
 
-    ac_lower_bound = floor(lom.options.best_lower_bound, digits = 3)
-
     if !(isapprox(result_dict["cheeger_constant"], 0, atol = 1E-5)) && (
         result_dict["cheeger_constant"] <
         (ac_lower_bound * lom.options.cheeger_cuts_factor)

src/data.jl

@@ -129,63 +129,52 @@ function parse_file(file_path::String)
 end
 
 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"]
     adjacency_augment_graph = data_dict["adjacency_augment_graph"]
 
-    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
-            end
-            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,
-            ))
-                Memento.error("Adjacency matrix of the augment graph has to be symmetric")
-            end
-            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))
-                num_edges_to_augment += 1
-            end
+    num_edges_existing, num_edges_to_augment = 0, 0
+
+    for i in 1:num_nodes-1, j in i+1:num_nodes
+        adjacency_base_graph[i, j] =
+            isapprox(abs(adjacency_base_graph[i, j]), 0, atol = 1E-6) ? 0 :
+            adjacency_base_graph[i, j]
+        adjacency_augment_graph[i, j] =
+            isapprox(abs(adjacency_augment_graph[i, j]), 0, atol = 1E-6) ? 0 :
+            adjacency_augment_graph[i, j]
+
+        if !isapprox(adjacency_base_graph[i, j], adjacency_base_graph[j, i], atol = 1E-5)
+            Memento.error(_LOGGER, "Adjacency matrix of the base graph must be symmetric")
+        end
+        if !isapprox(
+            adjacency_augment_graph[i, j],
+            adjacency_augment_graph[j, i],
+            atol = 1E-5,
+        )
+            Memento.error(
+                _LOGGER,
+                "Adjacency matrix of the augment graph must be symmetric",
+            )
         end
-    end
 
-    # 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)
+        num_edges_existing +=
+            !isapprox(adjacency_base_graph[i, j], 0, atol = 1E-6) ? 1 : 0
+        num_edges_to_augment +=
+            !isapprox(adjacency_augment_graph[i, j], 0, atol = 1E-6) ? 1 : 0
     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["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["is_base_graph_connected"] = is_base_graph_connected
+    is_base_graph_connected =
+        (num_edges_existing > 0) &&
+        !isapprox(abs(LOpt.algebraic_connectivity(adjacency_base_graph)), 0, atol = 1E-6)
 
-    return data_dict_new
+    return Dict{String,Any}(
+        "num_nodes" => num_nodes,
+        "adjacency_base_graph" => adjacency_base_graph,
+        "adjacency_augment_graph" => adjacency_augment_graph,
+        "num_edges_existing" => num_edges_existing,
+        "num_edges_to_augment" => num_edges_to_augment,
+        "is_base_graph_connected" => is_base_graph_connected,
+    )
 end
 
 """
@@ -213,57 +202,56 @@ 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})
     num_nodes = data["num_nodes"]
+    num_edges_existing = data["num_edges_existing"]
+    num_edges_to_augment = data["num_edges_to_augment"]
+    augment_budget = data["augment_budget"]
 
-    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
+    if num_edges_to_augment == 0
+        Memento.error(_LOGGER, "At least one edge is required for augmentation.")
+    elseif num_edges_existing == num_nodes * (num_nodes - 1) / 2
         Memento.error(
             _LOGGER,
-            "Input graph is already a complete graph; augmentation is unnecessary",
+            "The graph is already complete; augmentation is unnecessary.",
         )
-    elseif (data["num_edges_existing"] == 0) && (data["augment_budget"] < (num_nodes - 1))
+    elseif num_edges_existing == 0 && augment_budget < (num_nodes - 1)
         Memento.error(
             _LOGGER,
-            "Detected trivial solutions with disconnected graphs. `augment_budget` may be insufficient.",
+            "Graph is disconnected; `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")
+    elseif !isinteger(augment_budget) || augment_budget < 0
+        Memento.error(_LOGGER, "Augmentation budget must be a positive integer.")
     end
 
-    # Detect mutually exclusive sets of edges between base and augment graph
+    # Check mutually exclusive edge sets
     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",
+            "Edge sets of base and augment graphs must be mutually exclusive.",
         )
     end
 
-    # Detect free vertices 
+    # Check for free vertices
     A = data["adjacency_augment_graph"] + data["adjacency_base_graph"]
     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.",
+                "Free vertices detected, leading to disconnected graphs.",
             )
         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"])
+    # Check Hamiltonian cycle feasibility
+    if data["graph_type"] == "hamiltonian_cycle" && num_edges_existing == 0
+        if !(num_edges_to_augment >= num_nodes && augment_budget == num_nodes)
             Memento.error(
                 _LOGGER,
-                "Detected infeasibility due to the number of augmentation edges incompatible for a hamiltonian cycle",
+                "Infeasibility due to incompatible augmentation for Hamiltonian cycle.",
             )
         end
     end