Implement node presolve

benchmarks/linear_programming/run_mps_files.sh

@@ -79,6 +79,7 @@ Optional Arguments:
     --batch-num        Batch number
     --n-batches        Number of batches
     --log-to-console   Log to console
+    --model-list       File containing a list of models to run
     -h, --help         Show this help message and exit
 
 Examples:
@@ -168,6 +169,11 @@ while [[ $# -gt 0 ]]; do
             LOG_TO_CONSOLE="$2"
             shift 2
             ;;
+        --model-list)
+            echo "MODEL_LIST: $2"
+            MODEL_LIST="$2"
+            shift 2
+            ;;
         *)
             echo "Unknown argument: $1"
             print_help
@@ -194,7 +200,7 @@ PRESOLVE=${PRESOLVE:-true}
 BATCH_NUM=${BATCH_NUM:-0}
 N_BATCHES=${N_BATCHES:-1}
 LOG_TO_CONSOLE=${LOG_TO_CONSOLE:-true}
-
+MODEL_LIST=${MODEL_LIST:-}
 # Determine GPU list
 if [[ -n "$CUDA_VISIBLE_DEVICES" ]]; then
     IFS=',' read -ra GPU_LIST <<< "$CUDA_VISIBLE_DEVICES"
@@ -206,8 +212,49 @@ else
 fi
 GPU_COUNT=${#GPU_LIST[@]}
 
-# Gather all mps files into an array
-mapfile -t mps_files < <(ls "$MPS_DIR"/*.mps)
+# Ensure all entries in MODEL_LIST have .mps extension
+if [[ -n "$MODEL_LIST" && -f "$MODEL_LIST" ]]; then
+    # Create a temporary file to store the updated model list
+    TMP_MODEL_LIST=$(mktemp)
+    while IFS= read -r line || [[ -n "$line" ]]; do
+        # Skip empty lines
+        [[ -z "$line" ]] && continue
+        # If the line does not end with .mps, append it
+        if [[ "$line" != *.mps ]]; then
+            echo "${line}.mps" >> "$TMP_MODEL_LIST"
+        else
+            echo "$line" >> "$TMP_MODEL_LIST"
+        fi
+    done < "$MODEL_LIST"
+    # Overwrite the original MODEL_LIST with the updated one
+    mv "$TMP_MODEL_LIST" "$MODEL_LIST"
+fi
+
+
+# Gather all mps files into an array, either from the model list or from the directory
+if [[ -n "$MODEL_LIST" ]]; then
+    if [[ ! -f "$MODEL_LIST" ]]; then
+        echo "Model list file not found: $MODEL_LIST"
+        exit 1
+    fi
+    mapfile -t mps_files < <(grep -v '^\s*$' "$MODEL_LIST" | sed "s|^|$MPS_DIR/|")
+    # Optionally, check that all files exist
+    missing_files=()
+    for f in "${mps_files[@]}"; do
+        if [[ ! -f "$f" ]]; then
+            missing_files+=("$f")
+        fi
+    done
+    if (( ${#missing_files[@]} > 0 )); then
+        echo "The following files from the model list do not exist in $MPS_DIR:"
+        for f in "${missing_files[@]}"; do
+            echo "  $f"
+        done
+        exit 1
+    fi
+else
+    mapfile -t mps_files < <(ls "$MPS_DIR"/*.mps)
+fi
 
 # Calculate batch size and start/end indices
 batch_size=$(( (${#mps_files[@]} + N_BATCHES - 1) / N_BATCHES ))

cpp/src/dual_simplex/branch_and_bound.cpp

@@ -529,6 +529,8 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
                          std::move(up_child));  // child pointers moved into the tree
   lp_problem_t leaf_problem =
     original_lp;  // Make a copy of the original LP. We will modify its bounds at each leaf
+  csc_matrix_t<i_t, f_t> Arow(1, 1, 1);
+  original_lp.A.transpose(Arow);
   f_t gap            = get_upper_bound() - lower_bound;
   i_t nodes_explored = 0;
   settings.log.printf(
@@ -631,7 +633,11 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
     // Set the correct bounds for the leaf problem
     leaf_problem.lower = original_lp.lower;
     leaf_problem.upper = original_lp.upper;
-    node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper);
+    std::vector<bool> bounds_changed(original_lp.num_cols, false);
+    // Technically, we can get the already strengthened bounds from the node/parent instead of
+    // getting it from the original problem and re-strengthening. But this requires storing
+    // two vectors at each node and potentially cause memory issues
+    node_ptr->get_variable_bounds(leaf_problem.lower, leaf_problem.upper, bounds_changed);
 
     std::vector<variable_status_t>& leaf_vstatus = node_ptr->vstatus;
     lp_solution_t<i_t, f_t> leaf_solution(leaf_problem.num_rows, leaf_problem.num_cols);
@@ -642,28 +648,43 @@ mip_status_t branch_and_bound_t<i_t, f_t>::solve(mip_solution_t<i_t, f_t>& solut
     std::vector<f_t> leaf_edge_norms      = edge_norms;  // = node.steepest_edge_norms;
     simplex_solver_settings_t lp_settings = settings;
     lp_settings.set_log(false);
-    lp_settings.cut_off      = upper_bound + settings.dual_tol;
-    lp_settings.inside_mip   = 2;
-    dual::status_t lp_status = dual_phase2(2,
-                                           0,
-                                           lp_start_time,
-                                           leaf_problem,
-                                           lp_settings,
-                                           leaf_vstatus,
-                                           leaf_solution,
-                                           node_iter,
-                                           leaf_edge_norms);
-    if (lp_status == dual::status_t::NUMERICAL) {
-      settings.log.printf("Numerical issue node %d. Resolving from scratch.\n", nodes_explored);
-      lp_status_t second_status = solve_linear_program_advanced(
-        leaf_problem, lp_start_time, lp_settings, leaf_solution, leaf_vstatus, leaf_edge_norms);
-      lp_status = convert_lp_status_to_dual_status(second_status);
+    lp_settings.cut_off    = upper_bound + settings.dual_tol;
+    lp_settings.inside_mip = 2;
+
+    // in B&B we only have equality constraints, leave it empty for default
+    std::vector<char> row_sense;
+    bool feasible =
+      bound_strengthening(row_sense, lp_settings, leaf_problem, Arow, var_types, bounds_changed);
+
+    dual::status_t lp_status = dual::status_t::DUAL_UNBOUNDED;
+
+    // If the problem is infeasible after bounds strengthening, we don't need to solve the LP
+    if (feasible) {
+      lp_status = dual_phase2(2,
+                              0,
+                              lp_start_time,
+                              leaf_problem,
+                              lp_settings,
+                              leaf_vstatus,
+                              leaf_solution,
+                              node_iter,
+                              leaf_edge_norms);
+      if (lp_status == dual::status_t::NUMERICAL) {
+        settings.log.printf("Numerical issue node %d. Resolving from scratch.\n", nodes_explored);
+        lp_status_t second_status = solve_linear_program_advanced(
+          leaf_problem, lp_start_time, lp_settings, leaf_solution, leaf_vstatus, leaf_edge_norms);
+        lp_status = convert_lp_status_to_dual_status(second_status);
+      }
     }
     total_lp_solve_time += toc(lp_start_time);
     total_lp_iters += node_iter;
 
     nodes_explored++;
     if (lp_status == dual::status_t::DUAL_UNBOUNDED) {
+      if (!feasible) {
+        settings.log.printf("Infeasible after bounds strengthening. Fathoming node %d.\n",
+                            nodes_explored);
+      }
       node_ptr->lower_bound = inf;
       std::vector<mip_node_t<i_t, f_t>*> stack;
       node_ptr->set_status(node_status_t::INFEASIBLE, stack);

cpp/src/dual_simplex/mip_node.hpp

@@ -79,22 +79,27 @@ class mip_node_t {
     children[1] = nullptr;
   }
 
-  void get_variable_bounds(std::vector<f_t>& lower, std::vector<f_t>& upper) const
+  void get_variable_bounds(std::vector<f_t>& lower,
+                           std::vector<f_t>& upper,
+                           std::vector<bool>& bounds_changed) const
   {
+    std::fill(bounds_changed.begin(), bounds_changed.end(), false);
     // Apply the bounds at the current node
     assert(lower.size() > branch_var);
     assert(upper.size() > branch_var);
     lower[branch_var]                = branch_var_lower;
     upper[branch_var]                = branch_var_upper;
+    bounds_changed[branch_var]       = true;
     mip_node_t<i_t, f_t>* parent_ptr = parent;
     while (parent_ptr != nullptr) {
       if (parent_ptr->node_id == 0) { break; }
       assert(parent_ptr->branch_var >= 0);
       assert(lower.size() > parent_ptr->branch_var);
       assert(upper.size() > parent_ptr->branch_var);
-      lower[parent_ptr->branch_var] = parent_ptr->branch_var_lower;
-      upper[parent_ptr->branch_var] = parent_ptr->branch_var_upper;
-      parent_ptr                    = parent_ptr->parent;
+      lower[parent_ptr->branch_var]          = parent_ptr->branch_var_lower;
+      upper[parent_ptr->branch_var]          = parent_ptr->branch_var_upper;
+      bounds_changed[parent_ptr->branch_var] = true;
+      parent_ptr                             = parent_ptr->parent;
     }
   }
 

cpp/src/dual_simplex/phase2.hpp

@@ -38,7 +38,22 @@ enum class status_t {
   CONCURRENT_LIMIT = 6,
   UNSET            = 7
 };
+
+static std::string status_to_string(status_t status)
+{
+  switch (status) {
+    case status_t::OPTIMAL: return "OPTIMAL";
+    case status_t::DUAL_UNBOUNDED: return "DUAL_UNBOUNDED";
+    case status_t::NUMERICAL: return "NUMERICAL";
+    case status_t::CUTOFF: return "CUTOFF";
+    case status_t::TIME_LIMIT: return "TIME_LIMIT";
+    case status_t::ITERATION_LIMIT: return "ITERATION_LIMIT";
+    case status_t::CONCURRENT_LIMIT: return "CONCURRENT_LIMIT";
+    case status_t::UNSET: return "UNSET";
+  }
+  return "UNKNOWN";
 }
+}  // namespace dual
 
 template <typename i_t, typename f_t>
 dual::status_t dual_phase2(i_t phase,