From 465d6e29daf7757c6e1a9e5acea2cd9fbe4b1668 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Fri, 19 Jun 2026 14:59:42 +0200 Subject: [PATCH] added routine for converting from standard form to range form Signed-off-by: Nicolas L. Guidotti --- cpp/src/dual_simplex/presolve.cpp | 149 ++++++++++++++++++++ cpp/src/dual_simplex/presolve.hpp | 13 ++ cpp/tests/dual_simplex/unit_tests/solve.cpp | 120 ++++++++++++++++ 3 files changed, 282 insertions(+) diff --git a/cpp/src/dual_simplex/presolve.cpp b/cpp/src/dual_simplex/presolve.cpp index d438fd133b..34cb1747c1 100644 --- a/cpp/src/dual_simplex/presolve.cpp +++ b/cpp/src/dual_simplex/presolve.cpp @@ -676,6 +676,147 @@ i_t add_artifical_variables(lp_problem_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 +void convert_simplex_problem(const lp_problem_t& simplex_problem, + const std::vector& var_types, + const simplex_solver_settings_t& settings, + const dualize_info_t& dualize_info, + const std::vector& new_slacks, + user_problem_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& simplex_A = simplex_problem.A; + assert((var_types.empty() || static_cast(var_types.size()) == n) && + "var_types must span the full simplex problem (structural + slack columns)"); + + const i_t num_slacks = static_cast(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 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(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& 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 void convert_user_problem(const user_problem_t& user_problem, const simplex_solver_settings_t& settings, @@ -1892,6 +2033,14 @@ template void convert_user_problem( std::vector& new_slacks, dualize_info_t& dualize_info); +template void convert_simplex_problem( + const lp_problem_t& simplex_problem, + const std::vector& var_types, + const simplex_solver_settings_t& settings, + const dualize_info_t& dualize_info, + const std::vector& new_slacks, + user_problem_t& user_problem); + template void convert_user_lp_with_guess( const user_problem_t& user_problem, const lp_solution_t& initial_solution, diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp index 22e578047f..96894c73db 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -205,6 +205,19 @@ void convert_user_problem(const user_problem_t& user_problem, std::vector& new_slacks, dualize_info_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 +void convert_simplex_problem(const lp_problem_t& simplex_problem, + const std::vector& var_types, + const simplex_solver_settings_t& settings, + const dualize_info_t& dualize_info, + const std::vector& new_slacks, + user_problem_t& user_problem); + template void convert_user_problem_with_guess(const user_problem_t& user_problem, const std::vector& guess, diff --git a/cpp/tests/dual_simplex/unit_tests/solve.cpp b/cpp/tests/dual_simplex/unit_tests/solve.cpp index 9bf8360ed1..06ac2f0d34 100644 --- a/cpp/tests/dual_simplex/unit_tests/solve.cpp +++ b/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 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 settings; + dual_simplex::lp_problem_t simplex_problem( + &handle, original.num_rows, original.num_cols, original.A.col_start[original.A.n]); + std::vector new_slacks; + dual_simplex::dualize_info_t 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(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 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 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{'L', 'E', 'L', 'E'})); + EXPECT_EQ(recovered.rhs, (std::vector{8.0, 4.0, -3.0, 1.0})); + EXPECT_EQ(recovered.num_range_rows, 1); + EXPECT_EQ(recovered.range_rows, (std::vector{3})); + EXPECT_EQ(recovered.range_value, (std::vector{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 simplex_again( + &handle, recovered.num_rows, recovered.num_cols, recovered.A.col_start[recovered.A.n]); + std::vector new_slacks_again; + dual_simplex::dualize_info_t 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