Conversion from the standard form to range form

cpp/src/dual_simplex/presolve.cpp

@@ -676,6 +676,147 @@ i_t add_artifical_variables(lp_problem_t<i_t, f_t>& problem,
   return 0;
 }
 
+// Inverse of convert_user_problem (LP/MIP only -- no quadratic or SOC handling):
+// take a problem in simplex standard form
+//   minimize   c^T x
+//   subject to A x = b
+//              l <= x <= u
+// and express it as a user_problem in range (row-bounded) form
+//   minimize   c^T x
+//   subject to l_row <= A x <= u_row
+//              l <= x <= u
+template <typename i_t, typename f_t>
+void convert_simplex_problem(const lp_problem_t<i_t, f_t>& simplex_problem,
+                             const std::vector<variable_type_t>& var_types,
+                             const simplex_solver_settings_t<i_t, f_t>& settings,
+                             const dualize_info_t<i_t, f_t>& dualize_info,
+                             const std::vector<i_t>& new_slacks,
+                             user_problem_t<i_t, f_t>& user_problem)
+{
+  constexpr bool verbose = false;
+  if (verbose) {
+    settings.log.printf("Converting simplex problem with %d rows and %d columns and %d nonzeros\n",
+                        simplex_problem.num_rows,
+                        simplex_problem.num_cols,
+                        simplex_problem.A.col_start[simplex_problem.num_cols]);
+  }
+
+  // The simplex_problem is expected to be the primal in standard form. The
+  // dualize path of convert_user_problem hands the dual to the solver and
+  // stashes the primal in dualize_info.primal_problem, so a faithful inverse
+  // of a dualized problem must start from that primal instead.
+  assert(!dualize_info.solving_dual &&
+         "convert_simplex_problem expects the primal problem; reconstruct the primal from "
+         "dualize_info.primal_problem before calling");
+
+  const i_t m                             = simplex_problem.num_rows;
+  const i_t n                             = simplex_problem.num_cols;
+  const csc_matrix_t<i_t, f_t>& simplex_A = simplex_problem.A;
+  assert((var_types.empty() || static_cast<i_t>(var_types.size()) == n) &&
+         "var_types must span the full simplex problem (structural + slack columns)");
+
+  const i_t num_slacks = static_cast<i_t>(new_slacks.size());
+  const i_t new_n      = n - num_slacks;
+  const i_t new_nnz    = simplex_A.col_start[n] - num_slacks;
+
+  // Reset the destination so we can populate user_problem in place.
+  user_problem.num_rows = m;
+  user_problem.num_cols = new_n;
+  user_problem.range_rows.clear();
+  user_problem.range_value.clear();
+  user_problem.objective.clear();
+  user_problem.lower.clear();
+  user_problem.upper.clear();
+  user_problem.var_types.clear();
+  user_problem.objective.reserve(new_n);
+  user_problem.lower.reserve(new_n);
+  user_problem.upper.reserve(new_n);
+  if (!var_types.empty()) { user_problem.var_types.reserve(new_n); }
+
+  // Rows with no associated slack stay equalities at their rhs; the slack loop
+  // below overrides the rows that had one.
+  user_problem.row_sense.assign(m, 'E');
+  user_problem.rhs = simplex_problem.rhs;
+
+  // Recover each row's constraint from its slack. The forward equality
+  // a_i^T x_struct + c * s = b_i (with l_s <= s <= u_s) gives
+  // a_i^T x_struct in [min(b - c l_s, b - c u_s), max(b - c l_s, b - c u_s)],
+  // which maps back to row_sense / rhs / range. Flag the slack columns so they
+  // can be dropped from the matrix below.
+  std::vector<bool> is_slack(n, false);
+  for (i_t s : new_slacks) {
+    assert(s >= 0 && s < n && "slack index out of range");
+    is_slack[s] = true;
+
+    const i_t col_start = simplex_A.col_start[s];
+    const i_t col_end   = simplex_A.col_start[s + 1];
+    assert(col_end - col_start == 1 && "slack/artificial column must have a single nonzero");
+    const i_t row = simplex_A.i[col_start];
+    const f_t c   = simplex_A.x[col_start];
+    const f_t b   = simplex_problem.rhs[row];
+    const f_t t1  = b - c * simplex_problem.lower[s];
+    const f_t t2  = b - c * simplex_problem.upper[s];
+    const f_t lo  = std::min(t1, t2);
+    const f_t hi  = std::max(t1, t2);
+    if (lo == hi) {
+      user_problem.row_sense[row] = 'E';
+      user_problem.rhs[row]       = lo;
+    } else if (lo == -inf) {
+      user_problem.row_sense[row] = 'L';
+      user_problem.rhs[row]       = hi;
+    } else if (hi == inf) {
+      user_problem.row_sense[row] = 'G';
+      user_problem.rhs[row]       = lo;
+    } else {
+      user_problem.row_sense[row] = 'E';
+      user_problem.rhs[row]       = lo;
+      user_problem.range_rows.push_back(row);
+      user_problem.range_value.push_back(hi - lo);
+    }
+  }
+  user_problem.num_range_rows = static_cast<i_t>(user_problem.range_rows.size());
+
+  // Rebuild the constraint matrix and column data in place, dropping the
+  // slack/artificial columns (each contributes exactly one nonzero). var_types
+  // is filtered to the kept (structural) columns in the same pass.
+  csc_matrix_t<i_t, f_t>& user_A = user_problem.A;
+  user_A.m                       = m;
+  user_A.n                       = new_n;
+  user_A.nz_max                  = new_nnz;
+  user_A.col_start.assign(new_n + 1, 0);
+  user_A.i.assign(new_nnz, 0);
+  user_A.x.assign(new_nnz, 0);
+
+  i_t nz    = 0;
+  i_t new_j = 0;
+  for (i_t j = 0; j < n; ++j) {
+    if (is_slack[j]) { continue; }
+    user_A.col_start[new_j] = nz;
+    for (i_t p = simplex_A.col_start[j]; p < simplex_A.col_start[j + 1]; ++p) {
+      user_A.i[nz] = simplex_A.i[p];
+      user_A.x[nz] = simplex_A.x[p];
+      ++nz;
+    }
+    user_problem.objective.push_back(simplex_problem.objective[j]);
+    user_problem.lower.push_back(simplex_problem.lower[j]);
+    user_problem.upper.push_back(simplex_problem.upper[j]);
+    if (!var_types.empty()) { user_problem.var_types.push_back(var_types[j]); }
+    ++new_j;
+  }
+  user_A.col_start[new_n] = nz;
+  assert(new_j == new_n);
+  assert(nz == new_nnz);
+
+  user_problem.obj_scale             = simplex_problem.obj_scale;
+  user_problem.obj_constant          = simplex_problem.obj_constant;
+  user_problem.objective_is_integral = simplex_problem.objective_is_integral;
+  user_problem.objective_step        = simplex_problem.objective_step;
+
+  // LP/MIP only: no quadratic objective and no second-order cones.
+  user_problem.cone_var_start = 0;
+  user_problem.second_order_cone_dims.clear();
+}
+
 template <typename i_t, typename f_t>
 void convert_user_problem(const user_problem_t<i_t, f_t>& user_problem,
                           const simplex_solver_settings_t<i_t, f_t>& settings,
@@ -1892,6 +2033,14 @@ template void convert_user_problem<int, double>(
   std::vector<int>& new_slacks,
   dualize_info_t<int, double>& dualize_info);
 
+template void convert_simplex_problem<int, double>(
+  const lp_problem_t<int, double>& simplex_problem,
+  const std::vector<variable_type_t>& var_types,
+  const simplex_solver_settings_t<int, double>& settings,
+  const dualize_info_t<int, double>& dualize_info,
+  const std::vector<int>& new_slacks,
+  user_problem_t<int, double>& user_problem);
+
 template void convert_user_lp_with_guess<int, double>(
   const user_problem_t<int, double>& user_problem,
   const lp_solution_t<int, double>& initial_solution,

cpp/src/dual_simplex/presolve.hpp

@@ -205,6 +205,19 @@ void convert_user_problem(const user_problem_t<i_t, f_t>& user_problem,
                           std::vector<i_t>& new_slacks,
                           dualize_info_t<i_t, f_t>& dualize_info);
 
+// Inverse of convert_user_problem (LP/MIP only): reconstruct a range-form
+// user_problem from a simplex standard-form lp_problem, dropping the
+// slack/artificial columns recorded in new_slacks. var_types spans the full
+// simplex problem (structural + slack columns) and is filtered to the
+// structural columns.
+template <typename i_t, typename f_t>
+void convert_simplex_problem(const lp_problem_t<i_t, f_t>& simplex_problem,
+                             const std::vector<variable_type_t>& var_types,
+                             const simplex_solver_settings_t<i_t, f_t>& settings,
+                             const dualize_info_t<i_t, f_t>& dualize_info,
+                             const std::vector<i_t>& new_slacks,
+                             user_problem_t<i_t, f_t>& user_problem);
+
 template <typename i_t, typename f_t>
 void convert_user_problem_with_guess(const user_problem_t<i_t, f_t>& user_problem,
                                      const std::vector<f_t>& guess,

cpp/tests/dual_simplex/unit_tests/solve.cpp

@@ -331,4 +331,124 @@ TEST(dual_simplex, dual_variable_greater_than)
   EXPECT_NEAR(solution.z[1], 0.0, 1e-6);
 }
 
+// Round-trip a MIP through convert_user_problem (range form -> simplex standard
+// form, appending one slack/artificial column per row) and
+// convert_simplex_problem (the inverse: drop the slacks and recover the row
+// bounds). The problem exercises every row type: '<=', '==', '>=', and a range
+// row.
+//
+// '>=' rows are folded into '<=' rows by negating their coefficients/rhs in the
+// forward pass, so the inverse recovers them as the equivalent negated '<=' row
+// rather than the original '>='. The recovered problem is therefore feasibly
+// equivalent but not textually identical. We assert two things:
+//   1. the directly-predictable recovered fields (sizes, row_sense, rhs, range,
+//      objective, bounds, var_types), and
+//   2. the round-trip invariant: re-running the forward conversion on the
+//      recovered problem reproduces the original standard-form problem exactly.
+TEST(dual_simplex, convert_simplex_problem_mip_round_trip)
+{
+  cuopt::init_logger_t log("", true);
+  namespace dual_simplex = cuopt::linear_programming::dual_simplex;
+  raft::handle_t handle{};
+
+  // minimize  x0 + 2 x1 + 3 x2
+  // subject to 2 x0 + x1        <= 8   (row 0, 'L')
+  //              x0        + x2   = 4   (row 1, 'E')
+  //                   x1   + x2  >= 3   (row 2, 'G')
+  //            1 <=  x0 + x1     <= 7   (row 3, range row stored as 'E' + range)
+  //            x0 integer in [0, 10]
+  //            x1 continuous in [0, inf)
+  //            x2 integer in [0, 5]
+  dual_simplex::user_problem_t<int, double> original(&handle);
+  constexpr int m  = 4;
+  constexpr int n  = 3;
+  constexpr int nz = 8;
+
+  original.num_rows = m;
+  original.num_cols = n;
+  original.objective.assign({1.0, 2.0, 3.0});
+  original.A.m      = m;
+  original.A.n      = n;
+  original.A.nz_max = nz;
+  original.A.reallocate(nz);
+  // column 0 (x0): rows {0, 1, 3} -> {2, 1, 1}
+  // column 1 (x1): rows {0, 2, 3} -> {1, 1, 1}
+  // column 2 (x2): rows {1, 2}    -> {1, 1}
+  original.A.col_start.assign({0, 3, 6, 8});
+  original.A.i = {0, 1, 3, 0, 2, 3, 1, 2};
+  original.A.x = {2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
+  original.rhs.assign({8.0, 4.0, 3.0, 1.0});
+  original.row_sense.assign({'L', 'E', 'G', 'E'});
+  original.lower.assign({0.0, 0.0, 0.0});
+  original.upper.assign({10.0, dual_simplex::inf, 5.0});
+  // row 3 is the range row 1 <= a^T x <= 7: stored as an 'E' range row, which
+  // convert_range_rows reads as [rhs, rhs + range] = [1, 7].
+  original.range_rows.assign({3});
+  original.range_value.assign({6.0});
+  original.num_range_rows = 1;
+  original.obj_constant   = 0.0;
+  original.var_types      = {dual_simplex::variable_type_t::INTEGER,
+                             dual_simplex::variable_type_t::CONTINUOUS,
+                             dual_simplex::variable_type_t::INTEGER};
+
+  // Forward: range form -> simplex standard form.
+  dual_simplex::simplex_solver_settings_t<int, double> settings;
+  dual_simplex::lp_problem_t<int, double> simplex_problem(
+    &handle, original.num_rows, original.num_cols, original.A.col_start[original.A.n]);
+  std::vector<int> new_slacks;
+  dual_simplex::dualize_info_t<int, double> dualize_info;
+  dual_simplex::convert_user_problem(original, settings, simplex_problem, new_slacks, dualize_info);
+
+  // Each row gets exactly one slack/artificial column, appended after the
+  // structural columns.
+  EXPECT_EQ(new_slacks.size(), static_cast<size_t>(m));
+  EXPECT_EQ(simplex_problem.num_cols, n + m);
+
+  // var_types spans the full simplex problem; the appended columns are
+  // continuous (mirrors full_variable_types in branch_and_bound).
+  std::vector<dual_simplex::variable_type_t> var_types = original.var_types;
+  var_types.resize(simplex_problem.num_cols, dual_simplex::variable_type_t::CONTINUOUS);
+
+  // Inverse: simplex standard form -> range form.
+  dual_simplex::user_problem_t<int, double> recovered(&handle);
+  dual_simplex::convert_simplex_problem(
+    simplex_problem, var_types, settings, dualize_info, new_slacks, recovered);
+
+  // (1) Directly-predictable recovered fields.
+  EXPECT_EQ(recovered.num_rows, m);
+  EXPECT_EQ(recovered.num_cols, n);
+  // row 0 '<=' -> 'L' rhs 8; row 1 '==' -> 'E' rhs 4;
+  // row 2 '>=' -> recovered as the negated 'L' (-x1 - x2 <= -3);
+  // row 3 range [1, 7] -> canonical 'E' range row: rhs 1 with range 6.
+  EXPECT_EQ(recovered.row_sense, (std::vector<char>{'L', 'E', 'L', 'E'}));
+  EXPECT_EQ(recovered.rhs, (std::vector<double>{8.0, 4.0, -3.0, 1.0}));
+  EXPECT_EQ(recovered.num_range_rows, 1);
+  EXPECT_EQ(recovered.range_rows, (std::vector<int>{3}));
+  EXPECT_EQ(recovered.range_value, (std::vector<double>{6.0}));
+  // Column data (objective / bounds / types) carries over unchanged.
+  EXPECT_EQ(recovered.objective, original.objective);
+  EXPECT_EQ(recovered.lower, original.lower);
+  EXPECT_EQ(recovered.upper, original.upper);
+  EXPECT_EQ(recovered.var_types, original.var_types);
+
+  // (2) Round-trip invariant: converting the recovered problem forward again
+  // must reproduce the original standard-form problem exactly.
+  dual_simplex::lp_problem_t<int, double> simplex_again(
+    &handle, recovered.num_rows, recovered.num_cols, recovered.A.col_start[recovered.A.n]);
+  std::vector<int> new_slacks_again;
+  dual_simplex::dualize_info_t<int, double> dualize_info_again;
+  dual_simplex::convert_user_problem(
+    recovered, settings, simplex_again, new_slacks_again, dualize_info_again);
+
+  EXPECT_EQ(simplex_again.num_rows, simplex_problem.num_rows);
+  EXPECT_EQ(simplex_again.num_cols, simplex_problem.num_cols);
+  EXPECT_EQ(simplex_again.A.col_start, simplex_problem.A.col_start);
+  EXPECT_EQ(simplex_again.A.i, simplex_problem.A.i);
+  EXPECT_EQ(simplex_again.A.x, simplex_problem.A.x);
+  EXPECT_EQ(simplex_again.rhs, simplex_problem.rhs);
+  EXPECT_EQ(simplex_again.lower, simplex_problem.lower);
+  EXPECT_EQ(simplex_again.upper, simplex_problem.upper);
+  EXPECT_EQ(new_slacks_again, new_slacks);
+}
+
 }  // namespace cuopt::linear_programming::dual_simplex::test