diff --git a/cpp/src/dual_simplex/basis_solves.cpp b/cpp/src/dual_simplex/basis_solves.cpp index 2ee5e1d01f..3b26b2a8b9 100644 --- a/cpp/src/dual_simplex/basis_solves.cpp +++ b/cpp/src/dual_simplex/basis_solves.cpp @@ -619,12 +619,13 @@ i_t basis_repair(const csc_matrix_t& A, const std::vector& slacks_needed, std::vector& basis_list, std::vector& nonbasic_list, + std::vector& superbasic_list, std::vector& vstatus) { const i_t m = A.m; const i_t n = A.n; assert(basis_list.size() == m); - assert(nonbasic_list.size() == n - m); + assert(nonbasic_list.size() + superbasic_list.size() == n - m); // Create slack_map std::vector slack_map(m); // slack_map[i] = j if column j is e_i @@ -649,6 +650,13 @@ i_t basis_repair(const csc_matrix_t& A, for (i_t k = 0; k < num_nonbasic; ++k) { nonbasic_map[nonbasic_list[k]] = k; } + // Create a superbasic_map + std::vector superbasic_map( + n, -1); // superbasic_map[j] = p if superbasic[p] = j, -1 if j is basic/nonbasic + const i_t num_superbasic = superbasic_list.size(); + for (i_t k = 0; k < num_superbasic; ++k) { + superbasic_map[superbasic_list[k]] = k; + } const i_t columns_to_replace = deficient.size(); for (i_t k = 0; k < columns_to_replace; ++k) { @@ -656,19 +664,25 @@ i_t basis_repair(const csc_matrix_t& A, const i_t replace_i = slacks_needed[k]; const i_t replace_j = slack_map[replace_i]; basis_list[deficient[k]] = replace_j; - assert(nonbasic_map[replace_j] != -1); - 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. - 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; + if (nonbasic_map[replace_j] != -1) { + nonbasic_list[nonbasic_map[replace_j]] = bad_j; + // This is the main issue. What value should bad_j take on. + 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); + } + } else if (superbasic_map[replace_j] != -1) { + superbasic_list[superbasic_map[replace_j]] = bad_j; + vstatus[bad_j] = variable_status_t::SUPERBASIC; } else { - assert(1 == 0); + assert(nonbasic_map[replace_j] != -1 || superbasic_map[replace_j] != -1); } + vstatus[replace_j] = variable_status_t::BASIC; } return 0; @@ -865,6 +879,7 @@ template int basis_repair(const csc_matrix_t& A, const std::vector& slacks_needed, std::vector& basis_list, std::vector& nonbasic_list, + std::vector& superbasic_list, std::vector& vstatus); template int form_b(const csc_matrix_t& A, diff --git a/cpp/src/dual_simplex/basis_solves.hpp b/cpp/src/dual_simplex/basis_solves.hpp index 295bedccdd..59b4725e42 100644 --- a/cpp/src/dual_simplex/basis_solves.hpp +++ b/cpp/src/dual_simplex/basis_solves.hpp @@ -48,6 +48,7 @@ i_t basis_repair(const csc_matrix_t& A, const std::vector& slacks_needed, std::vector& basis_list, std::vector& nonbasic_list, + std::vector& superbasic_list, std::vector& vstatus); // Form the basis matrix B = A(:, basic_list) diff --git a/cpp/src/dual_simplex/basis_updates.cpp b/cpp/src/dual_simplex/basis_updates.cpp index 1d94f41c7f..dd262622c2 100644 --- a/cpp/src/dual_simplex/basis_updates.cpp +++ b/cpp/src/dual_simplex/basis_updates.cpp @@ -2265,6 +2265,7 @@ int basis_update_mpf_t::refactor_basis( { std::vector deficient; std::vector slacks_needed; + std::vector superbasic_list; // Empty superbasic list if (L0_.m != A.m) { resize(A.m); } std::vector q; @@ -2281,8 +2282,16 @@ int basis_update_mpf_t::refactor_basis( if (status == CONCURRENT_HALT_RETURN) { return CONCURRENT_HALT_RETURN; } if (status == -1) { settings.log.debug("Initial factorization failed\n"); - basis_repair( - A, settings, lower, upper, deficient, slacks_needed, basic_list, nonbasic_list, vstatus); + basis_repair(A, + settings, + lower, + upper, + deficient, + slacks_needed, + basic_list, + nonbasic_list, + superbasic_list, + vstatus); #ifdef CHECK_BASIS_REPAIR const i_t m = A.m; diff --git a/cpp/src/dual_simplex/branch_and_bound.cpp b/cpp/src/dual_simplex/branch_and_bound.cpp index 7f6e2c1921..d354159046 100644 --- a/cpp/src/dual_simplex/branch_and_bound.cpp +++ b/cpp/src/dual_simplex/branch_and_bound.cpp @@ -1991,7 +1991,9 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_relax_soln_, iter, edge_norms_); - f_t dual_phase2_time = toc(dual_phase2_start_time); + exploration_stats_.total_lp_iters += iter; + root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + f_t dual_phase2_time = toc(dual_phase2_start_time); if (dual_phase2_time > 1.0) { settings_.log.debug("Dual phase2 time %.2f seconds\n", dual_phase2_time); } @@ -2002,11 +2004,27 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } if (cut_status != dual::status_t::OPTIMAL) { - settings_.log.printf("Cut status %s\n", dual::status_to_string(cut_status).c_str()); - return mip_status_t::NUMERICAL; + settings_.log.printf("Numerical issue at root node. Resolving from scratch\n"); + lp_status_t scratch_status = + solve_linear_program_with_advanced_basis(original_lp_, + exploration_stats_.start_time, + lp_settings, + root_relax_soln_, + basis_update, + basic_list, + nonbasic_list, + root_vstatus_, + edge_norms_); + if (scratch_status == lp_status_t::OPTIMAL) { + // We recovered + cut_status = convert_lp_status_to_dual_status(scratch_status); + exploration_stats_.total_lp_iters += root_relax_soln_.iterations; + root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + } else { + settings_.log.printf("Cut status %s\n", dual::status_to_string(cut_status).c_str()); + return mip_status_t::NUMERICAL; + } } - exploration_stats_.total_lp_iters += root_relax_soln_.iterations; - root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); local_lower_bounds_.assign(settings_.num_bfs_workers, root_objective_); diff --git a/cpp/src/dual_simplex/crossover.cpp b/cpp/src/dual_simplex/crossover.cpp index 81d5ec1e6d..afc8c66744 100644 --- a/cpp/src/dual_simplex/crossover.cpp +++ b/cpp/src/dual_simplex/crossover.cpp @@ -516,6 +516,7 @@ i_t dual_push(const lp_problem_t& lp, slacks_needed, basic_list, nonbasic_list, + superbasic_list, vstatus); rank = factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed); @@ -585,6 +586,55 @@ f_t primal_residual(const lp_problem_t& lp, const lp_solution_t(primal_residual); } +template +void find_primal_superbasic_variables(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const lp_solution_t& initial_solution, + lp_solution_t& solution, + std::vector& vstatus, + std::vector& nonbasic_list, + std::vector& superbasic_list) +{ + const i_t n = lp.num_cols; + const f_t fixed_tolerance = settings.fixed_tol; + constexpr f_t basis_threshold = 1e-6; + nonbasic_list.clear(); + superbasic_list.clear(); + + for (i_t j = 0; j < n; ++j) { + if (vstatus[j] != variable_status_t::BASIC) { + const f_t lower_infeas = lp.lower[j] - initial_solution.x[j]; + const f_t lower_bound_slack = initial_solution.x[j] - lp.lower[j]; + const f_t upper_infeas = initial_solution.x[j] - lp.upper[j]; + const f_t upper_bound_slack = lp.upper[j] - initial_solution.x[j]; + if (std::abs(lp.lower[j] - lp.upper[j]) < fixed_tolerance) { + vstatus[j] = variable_status_t::NONBASIC_FIXED; + nonbasic_list.push_back(j); + } else if (lower_infeas > 0 && lp.lower[j] > -inf) { + vstatus[j] = variable_status_t::NONBASIC_LOWER; + solution.x[j] = lp.lower[j]; + nonbasic_list.push_back(j); + } else if (upper_infeas > 0 && lp.upper[j] < inf) { + vstatus[j] = variable_status_t::NONBASIC_UPPER; + solution.x[j] = lp.upper[j]; + nonbasic_list.push_back(j); + } else if (lower_bound_slack < basis_threshold && lp.lower[j] > -inf) { + vstatus[j] = variable_status_t::NONBASIC_LOWER; + nonbasic_list.push_back(j); + } else if (upper_bound_slack < basis_threshold && lp.upper[j] < inf) { + vstatus[j] = variable_status_t::NONBASIC_UPPER; + nonbasic_list.push_back(j); + } else if (lp.lower[j] == -inf && lp.upper[j] == inf) { + vstatus[j] = variable_status_t::NONBASIC_FREE; + nonbasic_list.push_back(j); + } else { + vstatus[j] = variable_status_t::SUPERBASIC; + superbasic_list.push_back(j); + } + } + } +} + template f_t primal_ratio_test(const lp_problem_t& lp, const simplex_solver_settings_t& settings, @@ -784,6 +834,7 @@ i_t primal_push(const lp_problem_t& lp, } basic_list[basic_leaving_index] = s; nonbasic_list.push_back(leaving_index); + superbasic_list.pop_back(); // Remove superbasic variable // Refactor or Update bool should_refactor = ft.num_updates() > 100; @@ -820,7 +871,11 @@ i_t primal_push(const lp_problem_t& lp, slacks_needed, basic_list, nonbasic_list, + superbasic_list, vstatus); + // We need to be careful. As basis_repair may have changed the superbasic list + find_primal_superbasic_variables( + lp, settings, solution, solution, vstatus, nonbasic_list, superbasic_list); rank = factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed); if (rank == CONCURRENT_HALT_RETURN) { @@ -844,11 +899,9 @@ i_t primal_push(const lp_problem_t& lp, vstatus[s] = variable_status_t::NONBASIC_UPPER; nonbasic_list.push_back(s); } + superbasic_list.pop_back(); // Remove superbasic variable } - // Remove superbasic variable - superbasic_list.pop_back(); - num_pushes++; if (num_pushes % settings.iteration_log_frequency == 0 || toc(last_print_time) > 10.0 || superbasic_list.size() == 0) { @@ -1178,6 +1231,7 @@ crossover_status_t crossover(const lp_problem_t& lp, slacks_needed, basic_list, nonbasic_list, + superbasic_list, vstatus); rank = factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed); if (rank == CONCURRENT_HALT_RETURN) { @@ -1210,41 +1264,8 @@ crossover_status_t crossover(const lp_problem_t& lp, settings.log.debug("nonbasic list size %ld n - m %d\n", nonbasic_list.size(), n - m); print_crossover_info(lp, settings, vstatus, solution, "Dual push complete"); - nonbasic_list.clear(); - superbasic_list.clear(); - - for (i_t j = 0; j < n; ++j) { - if (vstatus[j] != variable_status_t::BASIC) { - const f_t lower_infeas = lp.lower[j] - initial_solution.x[j]; - const f_t lower_bound_slack = initial_solution.x[j] - lp.lower[j]; - const f_t upper_infeas = initial_solution.x[j] - lp.upper[j]; - const f_t upper_bound_slack = lp.upper[j] - initial_solution.x[j]; - if (std::abs(lp.lower[j] - lp.upper[j]) < fixed_tolerance) { - vstatus[j] = variable_status_t::NONBASIC_FIXED; - nonbasic_list.push_back(j); - } else if (lower_infeas > 0 && lp.lower[j] > -inf) { - vstatus[j] = variable_status_t::NONBASIC_LOWER; - solution.x[j] = lp.lower[j]; - nonbasic_list.push_back(j); - } else if (upper_infeas > 0 && lp.upper[j] < inf) { - vstatus[j] = variable_status_t::NONBASIC_UPPER; - solution.x[j] = lp.upper[j]; - nonbasic_list.push_back(j); - } else if (lower_bound_slack < basis_threshold && lp.lower[j] > -inf) { - vstatus[j] = variable_status_t::NONBASIC_LOWER; - nonbasic_list.push_back(j); - } else if (upper_bound_slack < basis_threshold && lp.upper[j] < inf) { - vstatus[j] = variable_status_t::NONBASIC_UPPER; - nonbasic_list.push_back(j); - } else if (lp.lower[j] == -inf && lp.upper[j] == inf) { - vstatus[j] = variable_status_t::NONBASIC_FREE; - nonbasic_list.push_back(j); - } else { - vstatus[j] = variable_status_t::SUPERBASIC; - superbasic_list.push_back(j); - } - } - } + find_primal_superbasic_variables( + lp, settings, initial_solution, solution, vstatus, nonbasic_list, superbasic_list); if (superbasic_list.size() > 0) { std::vector save_x = solution.x; @@ -1381,6 +1402,7 @@ crossover_status_t crossover(const lp_problem_t& lp, slacks_needed, basic_list, nonbasic_list, + superbasic_list, vstatus); rank = factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed); diff --git a/cpp/src/dual_simplex/primal.cpp b/cpp/src/dual_simplex/primal.cpp index 38cddc0e24..98f5f4193b 100644 --- a/cpp/src/dual_simplex/primal.cpp +++ b/cpp/src/dual_simplex/primal.cpp @@ -308,6 +308,7 @@ primal::status_t primal_phase2(i_t phase, slacks_needed, basic_list, nonbasic_list, + superbasic_list, vstatus); rank = factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed); if (rank == CONCURRENT_HALT_RETURN) {