Skip to content

Fix incorrect infeasible list #694

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
chris-maes wants to merge 1 commit into
base: release/25.12
Choose a base branch
from
Open

Fix incorrect infeasible list #694

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 13 additions & 1 deletion cpp/src/dual_simplex/basis_solves.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -613,6 +613,8 @@ i_t factorize_basis(const csc_matrix_t<i_t, f_t>& A,
template <typename i_t, typename f_t>
i_t basis_repair(const csc_matrix_t<i_t, f_t>& A,
const simplex_solver_settings_t<i_t, f_t>& settings,
const std::vector<f_t>& lower,
const std::vector<f_t>& upper,
const std::vector<i_t>& deficient,
const std::vector<i_t>& slacks_needed,
std::vector<i_t>& basis_list,
Expand Down Expand Up @@ -658,7 +660,15 @@ i_t basis_repair(const csc_matrix_t<i_t, f_t>& A,
nonbasic_list[nonbasic_map[replace_j]] = bad_j;
vstatus[replace_j] = variable_status_t::BASIC;
// This is the main issue. What value should bad_j take on.
vstatus[bad_j] = variable_status_t::NONBASIC_FREE;
if (lower[bad_j] == -inf && upper[bad_j] == inf) {
vstatus[bad_j] = variable_status_t::NONBASIC_FREE;
} else if (lower[bad_j] > -inf) {
vstatus[bad_j] = variable_status_t::NONBASIC_LOWER;
} else if (upper[bad_j] < inf) {
vstatus[bad_j] = variable_status_t::NONBASIC_UPPER;
} else {
assert(1 == 0);
}
}

return 0;
Expand Down Expand Up @@ -849,6 +859,8 @@ template int factorize_basis<int>(const csc_matrix_t<int, double>& A,

template int basis_repair<int, double>(const csc_matrix_t<int, double>& A,
const simplex_solver_settings_t<int, double>& settings,
const std::vector<double>& lower,
const std::vector<double>& upper,
const std::vector<int>& deficient,
const std::vector<int>& slacks_needed,
std::vector<int>& basis_list,
Expand Down
2 changes: 2 additions & 0 deletions cpp/src/dual_simplex/basis_solves.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,8 @@ i_t factorize_basis(const csc_matrix_t<i_t, f_t>& A,
template <typename i_t, typename f_t>
i_t basis_repair(const csc_matrix_t<i_t, f_t>& A,
const simplex_solver_settings_t<i_t, f_t>& settings,
const std::vector<f_t>& lower,
const std::vector<f_t>& upper,
const std::vector<i_t>& deficient,
const std::vector<i_t>& slacks_needed,
std::vector<i_t>& basis_list,
Expand Down
4 changes: 3 additions & 1 deletion cpp/src/dual_simplex/basis_updates.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2046,6 +2046,8 @@ template <typename i_t, typename f_t>
int basis_update_mpf_t<i_t, f_t>::refactor_basis(
const csc_matrix_t<i_t, f_t>& A,
const simplex_solver_settings_t<i_t, f_t>& settings,
const std::vector<f_t>& lower,
const std::vector<f_t>& upper,
std::vector<i_t>& basic_list,
std::vector<i_t>& nonbasic_list,
std::vector<variable_status_t>& vstatus)
Expand All @@ -2066,7 +2068,7 @@ int basis_update_mpf_t<i_t, f_t>::refactor_basis(
deficient,
slacks_needed) == -1) {
settings.log.debug("Initial factorization failed\n");
basis_repair(A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
basis_repair(A, settings, lower, upper, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);

#ifdef CHECK_BASIS_REPAIR
const i_t m = A.m;
Expand Down
2 changes: 2 additions & 0 deletions cpp/src/dual_simplex/basis_updates.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -373,6 +373,8 @@ class basis_update_mpf_t {
// Compute L*U = A(p, basic_list)
int refactor_basis(const csc_matrix_t<i_t, f_t>& A,
const simplex_solver_settings_t<i_t, f_t>& settings,
const std::vector<f_t>& lower,
const std::vector<f_t>& upper,
std::vector<i_t>& basic_list,
std::vector<i_t>& nonbasic_list,
std::vector<variable_status_t>& vstatus);
Expand Down
6 changes: 3 additions & 3 deletions cpp/src/dual_simplex/crossover.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -786,7 +786,7 @@ i_t primal_push(const lp_problem_t<i_t, f_t>& lp,
if (rank != m) {
settings.log.debug("Failed to factorize basis. rank %d m %d\n", rank, m);
basis_repair(
lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
lp.A, settings, lp.lower, lp.upper, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
if (factorize_basis(
lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) == -1) {
settings.log.printf("Failed to factorize basis after repair. rank %d m %d\n", rank, m);
Expand Down Expand Up @@ -1132,7 +1132,7 @@ crossover_status_t crossover(const lp_problem_t<i_t, f_t>& lp,
rank = factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed);
if (rank != m) {
settings.log.debug("Failed to factorize basis. rank %d m %d\n", rank, m);
basis_repair(lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
basis_repair(lp.A, settings, lp.lower, lp.upper, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
if (factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) ==
-1) {
settings.log.printf("Failed to factorize basis after repair. rank %d m %d\n", rank, m);
Expand Down Expand Up @@ -1323,7 +1323,7 @@ crossover_status_t crossover(const lp_problem_t<i_t, f_t>& lp,
factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed);
if (rank != m) {
settings.log.debug("Failed to factorize basis. rank %d m %d\n", rank, m);
basis_repair(lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
basis_repair(lp.A, settings, lp.lower, lp.upper, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
if (factorize_basis(
lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) == -1) {
settings.log.printf("Failed to factorize basis after repair. rank %d m %d\n", rank, m);
Expand Down
45 changes: 25 additions & 20 deletions cpp/src/dual_simplex/phase2.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -623,14 +623,17 @@ f_t compute_initial_primal_infeasibilities(const lp_problem_t<i_t, f_t>& lp,
const std::vector<i_t>& basic_list,
const std::vector<f_t>& x,
std::vector<f_t>& squared_infeasibilities,
std::vector<i_t>& infeasibility_indices)
std::vector<i_t>& infeasibility_indices,
f_t& primal_inf)
{
const i_t m = lp.num_rows;
const i_t n = lp.num_cols;
squared_infeasibilities.resize(n, 0.0);
squared_infeasibilities.resize(n);
std::fill(squared_infeasibilities.begin(), squared_infeasibilities.end(), 0.0);
infeasibility_indices.reserve(n);
infeasibility_indices.clear();
f_t primal_inf = 0.0;
f_t primal_inf_squared = 0.0;
primal_inf = 0.0;
for (i_t k = 0; k < m; ++k) {
const i_t j = basic_list[k];
const f_t lower_infeas = lp.lower[j] - x[j];
Expand All @@ -640,10 +643,11 @@ f_t compute_initial_primal_infeasibilities(const lp_problem_t<i_t, f_t>& lp,
const f_t square_infeas = infeas * infeas;
squared_infeasibilities[j] = square_infeas;
infeasibility_indices.push_back(j);
primal_inf += square_infeas;
primal_inf_squared += square_infeas;
primal_inf += infeas;
}
}
return primal_inf;
return primal_inf_squared;
}

template <typename i_t, typename f_t>
Expand Down Expand Up @@ -2241,7 +2245,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
assert(superbasic_list.size() == 0);
assert(nonbasic_list.size() == n - m);

if (ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus) > 0) {
if (ft.refactor_basis(lp.A, settings, lp.lower, lp.upper, basic_list, nonbasic_list, vstatus) > 0) {
return dual::status_t::NUMERICAL;
}

Expand All @@ -2268,7 +2272,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,

#ifdef COMPUTE_DUAL_RESIDUAL
std::vector<f_t> dual_res1;
compute_dual_residual(lp.A, objective, y, z, dual_res1);
phase2::compute_dual_residual(lp.A, objective, y, z, dual_res1);
f_t dual_res_norm = vector_norm_inf<i_t, f_t>(dual_res1);
if (dual_res_norm > settings.tight_tol) {
settings.log.printf("|| A'*y + z - c || %e\n", dual_res_norm);
Expand Down Expand Up @@ -2357,8 +2361,9 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
std::vector<uint8_t> bounded_variables(n, 0);
phase2::compute_bounded_info(lp.lower, lp.upper, bounded_variables);

f_t primal_infeasibility = phase2::compute_initial_primal_infeasibilities(
lp, settings, basic_list, x, squared_infeasibilities, infeasibility_indices);
f_t primal_infeasibility;
f_t primal_infeasibility_squared = phase2::compute_initial_primal_infeasibilities(
lp, settings, basic_list, x, squared_infeasibilities, infeasibility_indices, primal_infeasibility);

#ifdef CHECK_BASIC_INFEASIBILITIES
phase2::check_basic_infeasibilities(basic_list, basic_mark, infeasibility_indices, 0);
Expand Down Expand Up @@ -2557,8 +2562,8 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
phase2::compute_primal_solution_from_basis(
lp, ft, basic_list, nonbasic_list, vstatus, unperturbed_x);
x = unperturbed_x;
primal_infeasibility = phase2::compute_initial_primal_infeasibilities(
lp, settings, basic_list, x, squared_infeasibilities, infeasibility_indices);
primal_infeasibility_squared = phase2::compute_initial_primal_infeasibilities(
lp, settings, basic_list, x, squared_infeasibilities, infeasibility_indices, primal_infeasibility);
settings.log.printf("Updated primal infeasibility: %e\n", primal_infeasibility);

objective = lp.objective;
Expand Down Expand Up @@ -2594,8 +2599,8 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
phase2::compute_primal_solution_from_basis(
lp, ft, basic_list, nonbasic_list, vstatus, unperturbed_x);
x = unperturbed_x;
primal_infeasibility = phase2::compute_initial_primal_infeasibilities(
lp, settings, basic_list, x, squared_infeasibilities, infeasibility_indices);
primal_infeasibility_squared = phase2::compute_initial_primal_infeasibilities(
lp, settings, basic_list, x, squared_infeasibilities, infeasibility_indices, primal_infeasibility);

const f_t orig_dual_infeas = phase2::dual_infeasibility(
lp, settings, vstatus, z, settings.tight_tol, settings.dual_tol);
Expand Down Expand Up @@ -2810,7 +2815,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
delta_xB_0_sparse.i,
squared_infeasibilities,
infeasibility_indices,
primal_infeasibility);
primal_infeasibility_squared);
// Update primal infeasibilities due to changes in basic variables
// from the leaving and entering variables
phase2::update_primal_infeasibilities(lp,
Expand All @@ -2822,7 +2827,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
scaled_delta_xB_sparse.i,
squared_infeasibilities,
infeasibility_indices,
primal_infeasibility);
primal_infeasibility_squared);
// Update the entering variable
phase2::update_single_primal_infeasibility(lp.lower,
lp.upper,
Expand Down Expand Up @@ -2883,14 +2888,14 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
#endif
if (should_refactor) {
bool should_recompute_x = false;
if (ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus) > 0) {
if (ft.refactor_basis(lp.A, settings, lp.lower, lp.upper, basic_list, nonbasic_list, vstatus) > 0) {
should_recompute_x = true;
settings.log.printf("Failed to factorize basis. Iteration %d\n", iter);
if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; }
i_t count = 0;
i_t deficient_size;
while ((deficient_size =
ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus)) > 0) {
ft.refactor_basis(lp.A, settings, lp.lower, lp.upper, basic_list, nonbasic_list, vstatus)) > 0) {
settings.log.printf("Failed to repair basis. Iteration %d. %d deficient columns.\n",
iter,
static_cast<int>(deficient_size));
Expand All @@ -2912,8 +2917,8 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
lp, ft, basic_list, nonbasic_list, vstatus, unperturbed_x);
x = unperturbed_x;
}
phase2::compute_initial_primal_infeasibilities(
lp, settings, basic_list, x, squared_infeasibilities, infeasibility_indices);
primal_infeasibility_squared = phase2::compute_initial_primal_infeasibilities(
lp, settings, basic_list, x, squared_infeasibilities, infeasibility_indices, primal_infeasibility);
}
#ifdef CHECK_BASIC_INFEASIBILITIES
phase2::check_basic_infeasibilities(basic_list, basic_mark, infeasibility_indices, 7);
Expand Down Expand Up @@ -2951,7 +2956,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
iter,
compute_user_objective(lp, obj),
infeasibility_indices.size(),
primal_infeasibility,
primal_infeasibility_squared,
sum_perturb,
now);
}
Expand Down
2 changes: 1 addition & 1 deletion cpp/src/dual_simplex/primal.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -298,7 +298,7 @@ primal::status_t primal_phase2(i_t phase,
factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed);
if (rank != m) {
settings.log.debug("Failed to factorize basis. rank %d m %d\n", rank, m);
basis_repair(lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
basis_repair(lp.A, settings, lp.lower, lp.upper, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
if (factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) ==
-1) {
settings.log.printf("Failed to factorize basis after repair. rank %d m %d\n", rank, m);
Expand Down
Loading