Add cuts to MIP solver

cpp/src/dual_simplex/basis_updates.cpp

@@ -1108,6 +1108,151 @@ i_t basis_update_t<i_t, f_t>::lower_triangular_multiply(const csc_matrix_t<i_t,
   return new_nz;
 }
 
+// Start of middle product form: basis_update_mpf_t
+
+template <typename i_t, typename f_t>
+i_t basis_update_mpf_t<i_t, f_t>::append_cuts(const csr_matrix_t<i_t, f_t>& cuts_basic)
+{
+  const i_t m = L0_.m;
+
+  // Solve for U^T W^T = C_B^T
+  // We do this one row at a time of C_B
+  csc_matrix_t<i_t, f_t> WT(m, cuts_basic.m, 0);
+
+  i_t WT_nz = 0;
+  for (i_t k = 0; k < cuts_basic.m; k++) {
+    sparse_vector_t<i_t, f_t> rhs(cuts_basic, k);
+    u_transpose_solve(rhs);
+    WT.col_start[k] = WT_nz;
+    for (i_t q = 0; q < rhs.i.size(); q++) {
+      WT.i.push_back(rhs.i[q]);
+      WT.x.push_back(rhs.x[q]);
+      WT_nz++;
+    }
+  }
+  WT.col_start[cuts_basic.m] = WT_nz;
+
+  csc_matrix_t<i_t, f_t> V(cuts_basic.m, m, 0);
+  if (num_updates_ > 0) {
+    // W = V T_0 ... T_{num_updates_ - 1}
+    // or V = W T_{num_updates_ - 1}^{-1} ... T_0^{-1}
+    // or V^T = T_0^{-T} ... T_{num_updates_ - 1}^{-T} W^T
+    // We can compute V^T column by column so that we have
+    // V^T(:, h) = T_0^{-T} ... T_{num_updates_ - 1}^{-T} W^T(:, h)
+    // or
+    // V(h, :) = T_0^{-T} ... T_{num_updates_ - 1}^{-T} W^T(:, h)
+
+    csr_matrix_t<i_t, f_t> V_row(cuts_basic.m, m, 0);
+    i_t V_nz           = 0;
+    const f_t zero_tol = 1e-13;
+    for (i_t h = 0; h < cuts_basic.m; h++) {
+      sparse_vector_t rhs(WT, h);
+      scatter_into_workspace(rhs);
+      i_t nz = rhs.i.size();
+      for (i_t k = num_updates_ - 1; k >= 0; --k) {
+        // T_k^{-T} = ( I - v u^T/(1 + u^T v))
+        // T_k^{-T} * b = b - v * (u^T * b) / (1 + u^T * v) = b - theta * v, theta = u^T b / mu
+
+        const i_t u_col = 2 * k;
+        const i_t v_col = 2 * k + 1;
+        const f_t mu    = mu_values_[k];
+
+        // dot = u^T * b
+        f_t dot         = dot_product(u_col, xi_workspace_, x_workspace_);
+        const f_t theta = dot / mu;
+        if (std::abs(theta) > zero_tol) {
+          add_sparse_column(S_, v_col, -theta, xi_workspace_, nz, x_workspace_);
+        }
+      }
+      gather_into_sparse_vector(nz, rhs);
+      V_row.row_start[h] = V_nz;
+      for (i_t q = 0; q < rhs.i.size(); q++) {
+        V_row.j.push_back(rhs.i[q]);
+        V_row.x.push_back(rhs.x[q]);
+        V_nz++;
+      }
+    }
+    V_row.row_start[cuts_basic.m] = V_nz;
+
+    V_row.to_compressed_col(V);
+  } else {
+    // W = V
+    WT.transpose(V);
+  }
+
+  // Extend u_i, v_i for i = 0, ..., num_updates_ - 1
+  S_.m += cuts_basic.m;
+
+  // Adjust L and U
+  // L = [ L0  0 ]
+  //     [ V   I ]
+
+  i_t V_nz = V.col_start[m];
+  i_t L_nz = L0_.col_start[m];
+  csc_matrix_t<i_t, f_t> new_L(m + cuts_basic.m, m + cuts_basic.m, L_nz + V_nz + cuts_basic.m);
+  L_nz = 0;
+  for (i_t j = 0; j < m; ++j) {
+    new_L.col_start[j]  = L_nz;
+    const i_t col_start = L0_.col_start[j];
+    const i_t col_end   = L0_.col_start[j + 1];
+    for (i_t p = col_start; p < col_end; ++p) {
+      new_L.i[L_nz] = L0_.i[p];
+      new_L.x[L_nz] = L0_.x[p];
+      L_nz++;
+    }
+    const i_t V_col_start = V.col_start[j];
+    const i_t V_col_end   = V.col_start[j + 1];
+    for (i_t p = V_col_start; p < V_col_end; ++p) {
+      new_L.i[L_nz] = V.i[p] + m;
+      new_L.x[L_nz] = V.x[p];
+      L_nz++;
+    }
+  }
+  for (i_t j = m; j < m + cuts_basic.m; ++j) {
+    new_L.col_start[j] = L_nz;
+    new_L.i[L_nz]      = j;
+    new_L.x[L_nz]      = 1.0;
+    L_nz++;
+  }
+  new_L.col_start[m + cuts_basic.m] = L_nz;
+
+  L0_ = new_L;
+
+  // Adjust U
+  // U = [ U0 0 ]
+  //     [ 0  I ]
+
+  i_t U_nz = U0_.col_start[m];
+  U0_.col_start.resize(m + cuts_basic.m + 1);
+  U0_.i.resize(U_nz + cuts_basic.m);
+  U0_.x.resize(U_nz + cuts_basic.m);
+  for (i_t k = m; k < m + cuts_basic.m; ++k) {
+    U0_.col_start[k] = U_nz;
+    U0_.i[U_nz]      = k;
+    U0_.x[U_nz]      = 1.0;
+    U_nz++;
+  }
+  U0_.col_start[m + cuts_basic.m] = U_nz;
+  U0_.n                           = m + cuts_basic.m;
+  U0_.m                           = m + cuts_basic.m;
+
+  compute_transposes();
+
+  // Adjust row_permutation_ and inverse_row_permutation_
+  row_permutation_.resize(m + cuts_basic.m);
+  inverse_row_permutation_.resize(m + cuts_basic.m);
+  for (i_t k = m; k < m + cuts_basic.m; ++k) {
+    row_permutation_[k] = k;
+  }
+  inverse_permutation(row_permutation_, inverse_row_permutation_);
+
+  // Adjust workspace sizes
+  xi_workspace_.resize(2 * (m + cuts_basic.m), 0);
+  x_workspace_.resize(m + cuts_basic.m, 0.0);
+
+  return 0;
+}
+
 template <typename i_t, typename f_t>
 void basis_update_mpf_t<i_t, f_t>::gather_into_sparse_vector(i_t nz,
                                                              sparse_vector_t<i_t, f_t>& out) const

cpp/src/dual_simplex/basis_updates.hpp

@@ -291,6 +291,8 @@ class basis_update_mpf_t {
     reset_stats();
   }
 
+  i_t append_cuts(const csr_matrix_t<i_t, f_t>& cuts_basic);
+
   f_t estimate_solution_density(f_t rhs_nz, f_t sum, i_t& num_calls, bool& use_hypersparse) const
   {
     num_calls++;

cpp/src/dual_simplex/solve.cpp

@@ -293,6 +293,156 @@ lp_status_t solve_linear_program_with_advanced_basis(
   return lp_status;
 }
 
+template <typename i_t, typename f_t>
+lp_status_t solve_linear_program_with_cuts(const f_t start_time,
+                                           const simplex_solver_settings_t<i_t, f_t>& settings,
+                                           const csr_matrix_t<i_t, f_t>& cuts,
+                                           const std::vector<f_t>& cut_rhs,
+                                           lp_problem_t<i_t, f_t>& lp,
+                                           lp_solution_t<i_t, f_t>& solution,
+                                           basis_update_mpf_t<i_t, f_t>& basis_update,
+                                           std::vector<i_t>& basic_list,
+                                           std::vector<i_t>& nonbasic_list,
+                                           std::vector<variable_status_t>& vstatus,
+                                           std::vector<f_t>& edge_norms)
+{
+  // Given a set of cuts: C*x <= d that are currently violated
+  // by the current solution x* (i.e. C*x* > d), this function
+  // adds the cuts into the LP and solves again.
+
+  const i_t p = cuts.m;
+  if (cut_rhs.size() != static_cast<size_t>(p)) {
+    settings.log.printf("cut_rhs must have the same number of rows as cuts\n");
+    return lp_status_t::NUMERICAL_ISSUES;
+  }
+  settings.log.printf("Number of cuts %d\n", p);
+  settings.log.printf("Original lp rows %d\n", lp.num_rows);
+  settings.log.printf("Original lp cols %d\n", lp.num_cols);
+
+  csr_matrix_t<i_t, f_t> new_A_row(lp.num_rows, lp.num_cols, 1);
+  lp.A.to_compressed_row(new_A_row);
+
+  new_A_row.append_rows(cuts);
+
+  csc_matrix_t<i_t, f_t> new_A_col(lp.num_rows + p, lp.num_cols, 1);
+  new_A_row.to_compressed_col(new_A_col);
+
+  // Add in slacks variables for the new rows
+  lp.lower.resize(lp.num_cols + p);
+  lp.upper.resize(lp.num_cols + p);
+  lp.objective.resize(lp.num_cols + p);
+  i_t nz = new_A_col.col_start[lp.num_cols];
+  new_A_col.col_start.resize(lp.num_cols + p + 1);
+  new_A_col.i.resize(nz + p);
+  new_A_col.x.resize(nz + p);
+  i_t k = lp.num_rows;
+  for (i_t j = lp.num_cols; j < lp.num_cols + p; j++) {
+    new_A_col.col_start[j] = nz;
+    new_A_col.i[nz]        = k++;
+    new_A_col.x[nz]        = 1.0;
+    nz++;
+    lp.lower[j]     = 0.0;
+    lp.upper[j]     = inf;
+    lp.objective[j] = 0.0;
+  }
+  new_A_col.col_start[lp.num_cols + p] = nz;
+  new_A_col.n                          = lp.num_cols + p;
+
+  lp.A         = new_A_col;
+  i_t old_rows = lp.num_rows;
+  lp.num_rows += p;
+  i_t old_cols = lp.num_cols;
+  lp.num_cols += p;
+
+
+  lp.rhs.resize(lp.num_rows);
+  for (i_t k = old_rows; k < old_rows + p; k++) {
+    const i_t h = k - old_rows;
+    lp.rhs[k]   = cut_rhs[h];
+  }
+
+  // Construct C_B = C(:, basic_list)
+  std::vector<i_t> C_col_degree(p, 0);
+  i_t cuts_nz = cuts.row_start[p];
+  for (i_t q = 0; q < cuts_nz; q++) {
+    const i_t j = cuts.j[q];
+    C_col_degree[j]++;
+  }
+
+  std::vector<i_t> in_basis(old_cols, 0);
+  const i_t num_basic = static_cast<i_t>(basic_list.size());
+  i_t C_B_nz          = 0;
+  for (i_t k = 0; k < num_basic; k++) {
+    const i_t j = basic_list[k];
+    in_basis[j] = 1;
+    C_B_nz += C_col_degree[j];
+  }
+
+  csr_matrix_t<i_t, f_t> C_B(num_basic, num_basic, C_B_nz);
+  nz = 0;
+  for (i_t i = 0; i < p; i++) {
+    C_B.row_start[i]    = nz;
+    const i_t row_start = cuts.row_start[i];
+    const i_t row_end   = cuts.row_start[i + 1];
+    for (i_t q = row_start; q < row_end; q++) {
+      const i_t j = cuts.j[q];
+      if (in_basis[j] == 0) { continue; }
+      C_B.j[nz] = j;
+      C_B.x[nz] = cuts.x[q];
+      nz++;
+    }
+  }
+  C_B.row_start[p] = nz;
+  settings.log.printf("predicted nz %d actual nz %d\n", C_B_nz, nz);
+  if (nz != C_B_nz) { return lp_status_t::NUMERICAL_ISSUES; }
+
+  // Adjust the basis update to include the new cuts
+  basis_update.append_cuts(C_B);
+
+  // Adjust the vstatus
+  vstatus.resize(lp.num_cols);
+  for (i_t j = old_cols; j < lp.num_cols; j++) {
+    vstatus[j] = variable_status_t::BASIC;
+  }
+
+  basic_list.resize(lp.num_rows, 0);
+  i_t h = old_cols;
+  for (i_t j = old_rows; j < lp.num_rows; j++) {
+    basic_list[j] = h++;
+  }
+  // Adjust the solution
+  solution.x.resize(lp.num_cols, 0.0);
+  solution.y.resize(lp.num_rows, 0.0);
+  solution.z.resize(lp.num_cols, 0.0);
+
+  // For now just clear the edge norms
+  edge_norms.clear();
+  i_t iter              = 0;
+  dual::status_t status = dual_phase2_with_advanced_basis(2,
+                                                          0,
+                                                          false,
+                                                          start_time,
+                                                          lp,
+                                                          settings,
+                                                          vstatus,
+                                                          basis_update,
+                                                          basic_list,
+                                                          nonbasic_list,
+                                                          solution,
+                                                          iter,
+                                                          edge_norms);
+
+  lp_status_t lp_status;
+  if (status == dual::status_t::OPTIMAL) { lp_status = lp_status_t::OPTIMAL; }
+  if (status == dual::status_t::DUAL_UNBOUNDED) { lp_status = lp_status_t::INFEASIBLE; }
+  if (status == dual::status_t::TIME_LIMIT) { lp_status = lp_status_t::TIME_LIMIT; }
+  if (status == dual::status_t::ITERATION_LIMIT) { lp_status = lp_status_t::ITERATION_LIMIT; }
+  if (status == dual::status_t::CONCURRENT_LIMIT) { lp_status = lp_status_t::CONCURRENT_LIMIT; }
+  if (status == dual::status_t::NUMERICAL) { lp_status = lp_status_t::NUMERICAL_ISSUES; }
+  if (status == dual::status_t::CUTOFF) { lp_status = lp_status_t::CUTOFF; }
+  return lp_status;
+}
+
 template <typename i_t, typename f_t>
 lp_status_t solve_linear_program_with_barrier(const user_problem_t<i_t, f_t>& user_problem,
                                               const simplex_solver_settings_t<i_t, f_t>& settings,
@@ -661,6 +811,19 @@ template lp_status_t solve_linear_program_with_advanced_basis(
   std::vector<variable_status_t>& vstatus,
   std::vector<double>& edge_norms);
 
+template lp_status_t solve_linear_program_with_cuts(
+  const double start_time,
+  const simplex_solver_settings_t<int, double>& settings,
+  const csr_matrix_t<int, double>& cuts,
+  const std::vector<double>& cut_rhs,
+  lp_problem_t<int, double>& lp,
+  lp_solution_t<int, double>& solution,
+  basis_update_mpf_t<int, double>& basis_update,
+  std::vector<int>& basic_list,
+  std::vector<int>& nonbasic_list,
+  std::vector<variable_status_t>& vstatus,
+  std::vector<double>& edge_norms);
+
 template lp_status_t solve_linear_program_with_barrier(
   const user_problem_t<int, double>& user_problem,
   const simplex_solver_settings_t<int, double>& settings,

cpp/src/dual_simplex/solve.hpp

@@ -61,6 +61,19 @@ lp_status_t solve_linear_program_with_advanced_basis(
   std::vector<variable_status_t>& vstatus,
   std::vector<f_t>& edge_norms);
 
+template <typename i_t, typename f_t>
+lp_status_t solve_linear_program_with_cuts(const f_t start_time,
+                                           const simplex_solver_settings_t<i_t, f_t>& settings,
+                                           const csr_matrix_t<i_t, f_t>& cuts,
+                                           const std::vector<f_t>& cut_rhs,
+                                           lp_problem_t<i_t, f_t>& lp,
+                                           lp_solution_t<i_t, f_t>& solution,
+                                           basis_update_mpf_t<i_t, f_t>& basis_update,
+                                           std::vector<i_t>& basic_list,
+                                           std::vector<i_t>& nonbasic_list,
+                                           std::vector<variable_status_t>& vstatus,
+                                           std::vector<f_t>& edge_norms);
+
 template <typename i_t, typename f_t>
 lp_status_t solve_linear_program_with_barrier(const user_problem_t<i_t, f_t>& user_problem,
                                               const simplex_solver_settings_t<i_t, f_t>& settings,