diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 6ce9a4f4d0..5ff757a92f 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -1947,6 +1947,12 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut exploration_stats_.nodes_explored = 0; original_lp_.A.to_compressed_row(Arow_); + settings_.log.printf("Reduced cost strengthening enabled: %d\n", + settings_.reduced_cost_strengthening); + + variable_bounds_t variable_bounds( + original_lp_, settings_, var_types_, Arow_, new_slacks_); + if (guess_.size() != 0) { raft::common::nvtx::range scope_guess("BB::check_initial_guess"); std::vector crushed_guess; @@ -2117,8 +2123,10 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut var_types_, basis_update, root_relax_soln_.x, + root_relax_soln_.y, basic_list, - nonbasic_list); + nonbasic_list, + variable_bounds); f_t cut_generation_time = toc(cut_start_time); if (cut_generation_time > 1.0) { settings_.log.debug("Cut generation time %.2f seconds\n", cut_generation_time); @@ -2135,7 +2143,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut i_t num_cuts = cut_pool.get_best_cuts(cuts_to_add, cut_rhs, cut_types); if (num_cuts == 0) { break; } cut_info.record_cut_types(cut_types); -#ifdef PRINT_CUT_POOL_TYPES +#if 1 cut_pool.print_cutpool_types(); print_cut_types("In LP ", cut_types, settings_); printf("Cut pool size: %d\n", cut_pool.pool_size()); @@ -2179,6 +2187,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_vstatus_, edge_norms_); var_types_.resize(original_lp_.num_cols, variable_type_t::CONTINUOUS); + variable_bounds.resize(original_lp_.num_cols); mutex_original_lp_.unlock(); f_t add_cuts_time = toc(add_cuts_start_time); if (add_cuts_time > 1.0) { @@ -2227,6 +2236,9 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } if (!feasible) { settings_.log.printf("Bound strengthening detected infeasibility\n"); +#ifdef WRITE_BOUND_STRENGTHENING_INFEASIBLE_MPS + original_lp_.write_mps("bound_strengthening_infeasible.mps"); +#endif return mip_status_t::INFEASIBLE; } @@ -2248,7 +2260,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut iter, edge_norms_); 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); @@ -2278,9 +2289,13 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut 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()); +#ifdef WRITE_CUT_INFEASIBLE_MPS + original_lp_.write_mps("cut_infeasible.mps"); +#endif return mip_status_t::NUMERICAL; } } + root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); f_t remove_cuts_start_time = tic(); mutex_original_lp_.lock(); @@ -2299,6 +2314,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut basic_list, nonbasic_list, basis_update); + variable_bounds.resize(original_lp_.num_cols); mutex_original_lp_.unlock(); f_t remove_cuts_time = toc(remove_cuts_start_time); if (remove_cuts_time > 1.0) { @@ -2327,8 +2343,9 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut f_t change_in_objective = root_objective_ - last_objective; const f_t factor = settings_.cut_change_threshold; const f_t min_objective = 1e-3; - if (change_in_objective <= factor * std::max(min_objective, std::abs(root_relax_objective))) { - settings_.log.debug( + if (0 && + change_in_objective <= factor * std::max(min_objective, std::abs(root_relax_objective))) { + settings_.log.printf( "Change in objective %.16e is less than 1e-3 of root relax objective %.16e\n", change_in_objective, root_relax_objective); diff --git a/cpp/src/cuts/cuts.cpp b/cpp/src/cuts/cuts.cpp index cc2555611a..0a9ed5a19f 100644 --- a/cpp/src/cuts/cuts.cpp +++ b/cpp/src/cuts/cuts.cpp @@ -12,17 +12,18 @@ #include +#include +#include + namespace cuopt::linear_programming::dual_simplex { template -void cut_pool_t::add_cut(cut_type_t cut_type, - const sparse_vector_t& cut, - f_t rhs) +void cut_pool_t::add_cut(cut_type_t cut_type, const inequality_t& cut) { // TODO: Need to deduplicate cuts and only add if the cut is not already in the pool - for (i_t p = 0; p < cut.i.size(); p++) { - const i_t j = cut.i[p]; + for (i_t p = 0; p < cut.size(); p++) { + const i_t j = cut.index(p); if (j >= original_vars_) { settings_.log.printf( "Cut has variable %d that is greater than original_vars_ %d\n", j, original_vars_); @@ -30,14 +31,14 @@ void cut_pool_t::add_cut(cut_type_t cut_type, } } - sparse_vector_t cut_squeezed; + inequality_t cut_squeezed; cut.squeeze(cut_squeezed); - if (cut_squeezed.i.size() == 0) { + if (cut_squeezed.size() == 0) { settings_.log.printf("Cut has no coefficients\n"); return; } - cut_storage_.append_row(cut_squeezed); - rhs_storage_.push_back(rhs); + cut_storage_.append_row(cut_squeezed.vector); + rhs_storage_.push_back(cut_squeezed.rhs); cut_type_.push_back(cut_type); cut_age_.push_back(0); } @@ -99,7 +100,6 @@ f_t cut_pool_t::cut_orthogonality(i_t i, i_t j) template void cut_pool_t::score_cuts(std::vector& x_relax) { - const f_t min_cut_distance = 1e-4; cut_distances_.resize(cut_storage_.m, 0.0); cut_norms_.resize(cut_storage_.m, 0.0); @@ -107,7 +107,7 @@ void cut_pool_t::score_cuts(std::vector& x_relax) for (i_t i = 0; i < cut_storage_.m; i++) { f_t violation; f_t cut_dist = cut_distance(i, x_relax, violation, cut_norms_[i]); - cut_distances_[i] = cut_dist <= min_cut_distance ? 0.0 : cut_dist; + cut_distances_[i] = cut_dist <= min_cut_distance_ ? 0.0 : cut_dist; if (verbose) { settings_.log.printf("Cut %d type %d distance %+e violation %+e cut_norm %e\n", i, @@ -138,7 +138,7 @@ void cut_pool_t::score_cuts(std::vector& x_relax) const i_t i = sorted_indices.back(); sorted_indices.pop_back(); - if (cut_distances_[i] <= min_cut_distance) { break; } + if (cut_distances_[i] <= min_cut_distance_) { break; } f_t cut_ortho = 1.0; const i_t best_cuts_size = best_cuts_.size(); @@ -171,6 +171,7 @@ i_t cut_pool_t::get_best_cuts(csr_matrix_t& best_cuts, best_cut_types.reserve(scored_cuts_); for (i_t i : best_cuts_) { + if (cut_distances_[i] <= min_cut_distance_) { continue; } sparse_vector_t cut(cut_storage_, i); cut.negate(); best_cuts.append_row(cut); @@ -180,7 +181,7 @@ i_t cut_pool_t::get_best_cuts(csr_matrix_t& best_cuts, age_cuts(); - return static_cast(best_cuts_.size()); + return static_cast(best_rhs.size()); } template @@ -274,8 +275,7 @@ i_t knapsack_generation_t::generate_knapsack_cuts( const std::vector& var_types, const std::vector& xstar, i_t knapsack_row, - sparse_vector_t& cut, - f_t& cut_rhs) + inequality_t& cut) { const bool verbose = false; // Get the row associated with the knapsack constraint @@ -346,44 +346,31 @@ i_t knapsack_generation_t::generate_knapsack_cuts( if (solution[k] == 0.0) { cover_size++; } } - cut.i.clear(); - cut.x.clear(); - cut.i.reserve(cover_size); - cut.x.reserve(cover_size); + cut.reserve(cover_size); + cut.clear(); h = 0; for (i_t k = 0; k < knapsack_inequality.i.size(); k++) { const i_t j = knapsack_inequality.i[k]; if (!is_slack_[j]) { - if (solution[h] == 0.0) { - cut.i.push_back(j); - cut.x.push_back(-1.0); - } + if (solution[h] == 0.0) { cut.push_back(j, -1.0); } h++; } } - cut_rhs = -cover_size + 1; + cut.rhs = -cover_size + 1; cut.sort(); // The cut is in the form: - sum_{j in cover} x_j >= -cover_size + 1 // Which is equivalent to: sum_{j in cover} x_j <= cover_size - 1 // Verify the cut is violated - f_t dot = cut.dot(xstar); - f_t violation = dot - cut_rhs; + f_t dot = cut.vector.dot(xstar); + f_t violation = dot - cut.rhs; if (verbose) { settings.log.printf("Knapsack cut %d violation %e < 0\n", knapsack_row, violation); } if (violation >= -tol) { return -1; } - -#ifdef PRINT_KNAPSACK_CUT - settings.log.printf("knapsack cut (cover %d): \n", cover_size); - for (i_t k = 0; k < cut.i.size(); k++) { - settings.log.printf("x%d coeff %g value %g\n", cut.i[k], -cut.x[k], xstar[cut.i[k]]); - } - settings.log.printf("cut_rhs %g\n", -cut_rhs); -#endif return 0; } @@ -560,8 +547,10 @@ void cut_generation_t::generate_cuts(const lp_problem_t& lp, const std::vector& var_types, basis_update_mpf_t& basis_update, const std::vector& xstar, + const std::vector& ystar, const std::vector& basic_list, - const std::vector& nonbasic_list) + const std::vector& nonbasic_list, + variable_bounds_t& variable_bounds) { // Generate Gomory and CG Cuts if (settings.mixed_integer_gomory_cuts != 0 || settings.strong_chvatal_gomory_cuts != 0) { @@ -587,7 +576,7 @@ void cut_generation_t::generate_cuts(const lp_problem_t& lp, // Generate MIR and CG cuts if (settings.mir_cuts != 0 || settings.strong_chvatal_gomory_cuts != 0) { f_t cut_start_time = tic(); - generate_mir_cuts(lp, settings, Arow, new_slacks, var_types, xstar); + generate_mir_cuts(lp, settings, Arow, new_slacks, var_types, xstar, ystar, variable_bounds); f_t cut_generation_time = toc(cut_start_time); if (cut_generation_time > 1.0) { settings.log.debug("MIR and CG cut generation time %.2f seconds\n", cut_generation_time); @@ -606,11 +595,10 @@ void cut_generation_t::generate_knapsack_cuts( { if (knapsack_generation_.num_knapsack_constraints() > 0) { for (i_t knapsack_row : knapsack_generation_.get_knapsack_constraints()) { - sparse_vector_t cut(lp.num_cols, 0); - f_t cut_rhs; + inequality_t cut(lp.num_cols); i_t knapsack_status = knapsack_generation_.generate_knapsack_cuts( - lp, settings, Arow, new_slacks, var_types, xstar, knapsack_row, cut, cut_rhs); - if (knapsack_status == 0) { cut_pool_.add_cut(cut_type_t::KNAPSACK, cut, cut_rhs); } + lp, settings, Arow, new_slacks, var_types, xstar, knapsack_row, cut); + if (knapsack_status == 0) { cut_pool_.add_cut(cut_type_t::KNAPSACK, cut); } } } } @@ -622,127 +610,96 @@ void cut_generation_t::generate_mir_cuts( csr_matrix_t& Arow, const std::vector& new_slacks, const std::vector& var_types, - const std::vector& xstar) + const std::vector& xstar, + const std::vector& ystar, + variable_bounds_t& variable_bounds) { - f_t mir_start_time = tic(); - mixed_integer_rounding_cut_t mir(lp, settings, new_slacks, xstar); + f_t mir_start_time = tic(); + constexpr bool verbose = false; + complemented_mixed_integer_rounding_cut_t complemented_mir(lp, settings, new_slacks); strong_cg_cut_t cg(lp, var_types, xstar); - std::vector slack_map(lp.num_rows, -1); - for (i_t slack : new_slacks) { - const i_t col_start = lp.A.col_start[slack]; - const i_t col_end = lp.A.col_start[slack + 1]; - const i_t col_len = col_end - col_start; - assert(col_len == 1); - const i_t i = lp.A.i[col_start]; - slack_map[i] = slack; - } + std::vector scores; + complemented_mir.compute_initial_scores_for_rows(lp, settings, Arow, xstar, ystar, scores); - // Compute initial scores for all rows - std::vector score(lp.num_rows, 0.0); + // Push all the scores onto the priority queue + std::priority_queue> score_queue; for (i_t i = 0; i < lp.num_rows; i++) { - const i_t row_start = Arow.row_start[i]; - const i_t row_end = Arow.row_start[i + 1]; - - const i_t row_nz = row_end - row_start; - i_t num_integer_in_row = 0; - i_t num_continuous_in_row = 0; - for (i_t p = row_start; p < row_end; p++) { - const i_t j = Arow.j[p]; - if (var_types[j] == variable_type_t::INTEGER) { - num_integer_in_row++; - } else { - num_continuous_in_row++; - } - } - - if (num_integer_in_row == 0) { - score[i] = 0.0; - - } else { - f_t nz_score = lp.num_cols - row_nz; - - const i_t slack = slack_map[i]; - assert(slack >= 0); - const f_t slack_value = xstar[slack]; - - f_t slack_score = -std::log10(1e-16 + std::abs(slack_value)); - - const f_t nz_weight = 1.0; - const f_t slack_weight = 1.0; - const f_t integer_weight = 1.0; - - score[i] = - nz_weight * nz_score + slack_weight * slack_score + integer_weight * num_integer_in_row; - } + score_queue.push(std::make_pair(scores[i], i)); } - // Sort the rows by score - std::vector sorted_indices; - best_score_last_permutation(score, sorted_indices); - // These data structures are used to track the rows that have been aggregated // The invariant is that aggregated_rows is empty and aggregated_mark is all zeros // at the beginning of each iteration of the for loop below std::vector aggregated_rows; std::vector aggregated_mark(lp.num_rows, 0); - const i_t max_cuts = std::min(lp.num_rows, 1000); - f_t work_estimate = 0.0; - for (i_t h = 0; h < max_cuts; h++) { - // Get the row with the highest score - const i_t i = sorted_indices.back(); - sorted_indices.pop_back(); - const f_t max_score = score[i]; + // Transform the relaxation solution + std::vector transformed_xstar; + complemented_mir.bound_substitution(lp, variable_bounds, var_types, xstar, transformed_xstar); - const i_t row_nz = Arow.row_start[i + 1] - Arow.row_start[i]; - const i_t slack = slack_map[i]; + const i_t max_cuts = std::min(lp.num_rows, 100000); + f_t work_estimate = 0.0; + i_t num_cuts = 0; + while (num_cuts < max_cuts && !score_queue.empty()) { + // Get the row with the highest score from the queue + auto [max_score, i] = score_queue.top(); + score_queue.pop(); + // skip stale score entries + if (max_score != scores[i]) { continue; } + + // Add the current row to the aggregated set + aggregated_mark[i] = 1; + aggregated_rows.push_back(i); + + const i_t row_nz = Arow.row_length(i); + const i_t slack = complemented_mir.slack_cols(i); const f_t slack_value = xstar[slack]; if (max_score <= 0.0) { break; } if (work_estimate > 2e9) { break; } - sparse_vector_t inequality(Arow, i); - work_estimate += inequality.i.size(); - f_t inequality_rhs = lp.rhs[i]; + inequality_t inequality(Arow, i, lp.rhs[i]); + work_estimate += inequality.size(); + const bool generate_cg_cut = settings.strong_chvatal_gomory_cuts != 0; - f_t fractional_part_rhs = fractional_part(inequality_rhs); + f_t fractional_part_rhs = fractional_part(inequality.rhs); if (generate_cg_cut && fractional_part_rhs > 1e-6 && fractional_part_rhs < (1 - 1e-6)) { // Try to generate a CG cut - sparse_vector_t cg_inequality = inequality; - f_t cg_inequality_rhs = inequality_rhs; - if (fractional_part(inequality_rhs) < 0.5) { + + inequality_t cg_inequality = inequality; + if (fractional_part(inequality.rhs) < 0.5) { // Multiply by -1 to force the fractional part to be greater than 0.5 - cg_inequality_rhs *= -1; cg_inequality.negate(); } - sparse_vector_t cg_cut(lp.num_cols, 0); - f_t cg_cut_rhs; - i_t cg_status = cg.generate_strong_cg_cut( - lp, settings, var_types, cg_inequality, cg_inequality_rhs, xstar, cg_cut, cg_cut_rhs); - if (cg_status == 0) { cut_pool_.add_cut(cut_type_t::CHVATAL_GOMORY, cg_cut, cg_cut_rhs); } + inequality_t cg_cut; + i_t cg_status = + cg.generate_strong_cg_cut(lp, settings, var_types, cg_inequality, xstar, cg_cut); + if (cg_status == 0) { cut_pool_.add_cut(cut_type_t::CHVATAL_GOMORY, cg_cut); } } + if (settings.mir_cuts == 0) { continue; } + // Remove the slack from the equality to get an inequality - work_estimate += inequality.i.size(); + work_estimate += inequality.size(); i_t negate_inequality = 1; - for (i_t k = 0; k < inequality.i.size(); k++) { - const i_t j = inequality.i[k]; + for (i_t k = 0; k < inequality.size(); k++) { + const i_t j = inequality.index(k); if (j == slack) { - if (inequality.x[k] != 1.0) { - if (inequality.x[k] == -1.0 && lp.lower[j] >= 0.0) { + if (inequality.coeff(k) != 1.0) { + if (inequality.coeff(k) == -1.0 && lp.lower[j] >= 0.0) { negate_inequality = 0; } else { settings.log.debug("Bad slack %d in inequality: aj %e lo %e up %e\n", j, - inequality.x[k], + inequality.coeff(k), lp.lower[j], lp.upper[j]); negate_inequality = -1; break; } } - inequality.x[k] = 0.0; + inequality.vector.x[k] = 0.0; } } @@ -751,112 +708,55 @@ void cut_generation_t::generate_mir_cuts( if (negate_inequality) { // inequaility'*x <= inequality_rhs // But for MIR we need: inequality'*x >= inequality_rhs - inequality_rhs *= -1; inequality.negate(); - work_estimate += inequality.i.size(); + work_estimate += inequality.size(); } // We should now have: inequality'*x >= inequality_rhs - // Transform the relaxation solution - std::vector transformed_xstar; - mir.relaxation_to_nonnegative(lp, xstar, transformed_xstar); - work_estimate += transformed_xstar.size(); + for (i_t k = 0; k < inequality.size(); k++) { + const i_t j = inequality.index(k); + if (var_types[j] == variable_type_t::INTEGER) { + if (transformed_xstar[j] > complemented_mir.new_upper(j) / 2.0) { + settings.log.printf("!!!!!! j %d transformed x_j %e new_upper_j/2.0 %e\n", + j, + transformed_xstar[j], + complemented_mir.new_upper(j) / 2.0); + } + } + } - sparse_vector_t cut(lp.num_cols, 0); - f_t cut_rhs; bool add_cut = false; i_t num_aggregated = 0; const i_t max_aggregated = 6; + f_t min_abs_multiplier = 1.0; + f_t max_abs_multiplier = 1.0; work_estimate += lp.num_cols; while (!add_cut && num_aggregated < max_aggregated) { - sparse_vector_t transformed_inequality; + inequality_t transformed_inequality; inequality.squeeze(transformed_inequality); - f_t transformed_rhs = inequality_rhs; - work_estimate += transformed_inequality.i.size(); - - mir.to_nonnegative(lp, transformed_inequality, transformed_rhs); - work_estimate += transformed_inequality.i.size(); - std::vector> transformed_cuts; - std::vector transformed_cut_rhs; - std::vector transformed_violations; - - // Generate cut for delta = 1 - { - sparse_vector_t cut_1(lp.num_cols, 0); - f_t cut_1_rhs; - mir.generate_cut_nonnegative( - transformed_inequality, transformed_rhs, var_types, cut_1, cut_1_rhs); - f_t cut_1_violation = mir.compute_violation(cut_1, cut_1_rhs, transformed_xstar); - if (cut_1_violation > 1e-6) { - transformed_cuts.push_back(cut_1); - transformed_cut_rhs.push_back(cut_1_rhs); - transformed_violations.push_back(cut_1_violation); - } - work_estimate += transformed_inequality.i.size(); - } - - // Generate a cut for delta = max { |a_j|, j in I} - { - f_t max_coeff = 0.0; - for (i_t k = 0; k < transformed_inequality.i.size(); k++) { - const i_t j = transformed_inequality.i[k]; - if (var_types[j] == variable_type_t::INTEGER) { - const f_t abs_aj = std::abs(transformed_inequality.x[k]); - if (abs_aj > max_coeff) { max_coeff = abs_aj; } - } - } - work_estimate += transformed_inequality.i.size(); - - if (max_coeff > 1e-6 && max_coeff != 1.0) { - sparse_vector_t scaled_inequality = transformed_inequality; - const i_t nz = transformed_inequality.i.size(); - for (i_t k = 0; k < nz; k++) { - scaled_inequality.x[k] /= max_coeff; - } - const f_t scaled_rhs = transformed_rhs / max_coeff; - sparse_vector_t cut_2(lp.num_cols, 0); - f_t cut_2_rhs; - mir.generate_cut_nonnegative(scaled_inequality, scaled_rhs, var_types, cut_2, cut_2_rhs); - f_t cut_2_violation = mir.compute_violation(cut_2, cut_2_rhs, transformed_xstar); - if (cut_2_violation > 1e-6) { - transformed_cuts.push_back(cut_2); - transformed_cut_rhs.push_back(cut_2_rhs); - transformed_violations.push_back(cut_2_violation); - } - work_estimate += 5 * transformed_inequality.i.size(); - } - } - - if (!transformed_violations.empty()) { - std::vector permuted(transformed_violations.size()); - std::iota(permuted.begin(), permuted.end(), 0); - std::sort(permuted.begin(), permuted.end(), [&](i_t i, i_t j) { - return transformed_violations[i] > transformed_violations[j]; - }); - work_estimate += transformed_violations.size() * std::log2(transformed_violations.size()); - // Get the biggest violation - const i_t best_index = permuted[0]; - f_t max_viol = transformed_violations[best_index]; - cut = transformed_cuts[best_index]; - cut_rhs = transformed_cut_rhs[best_index]; - - if (max_viol > 1e-6) { - // TODO: Divide by 1/2*violation, 1/4*violation, 1/8*violation - // Transform back to the original variables - mir.to_original(lp, cut, cut_rhs); - mir.remove_small_coefficients(lp.lower, lp.upper, cut, cut_rhs); - mir.substitute_slacks(lp, Arow, cut, cut_rhs); - f_t viol = mir.compute_violation(cut, cut_rhs, xstar); - work_estimate += 10 * cut.i.size(); - add_cut = true; - } + work_estimate += transformed_inequality.size(); + + complemented_mir.transform_inequality(variable_bounds, var_types, transformed_inequality); + work_estimate += transformed_inequality.size(); + + inequality_t cut; + bool cut_found = complemented_mir.cut_generation_heuristic( + transformed_inequality, var_types, transformed_xstar, cut, work_estimate); + // Note cut is in the transformed variables + + if (cut_found) { + // Transform back to the original variables + complemented_mir.untransform_inequality(variable_bounds, var_types, cut); + complemented_mir.remove_small_coefficients(lp.lower, lp.upper, cut); + complemented_mir.substitute_slacks(lp, Arow, cut); + f_t viol = complemented_mir.compute_violation(cut, xstar); + work_estimate += 10 * cut.size(); + if (viol > 1e-6) { add_cut = true; } } if (add_cut) { - if (settings.mir_cuts != 0) { - cut_pool_.add_cut(cut_type_t::MIXED_INTEGER_ROUNDING, cut, cut_rhs); - } + if (settings.mir_cuts != 0) { cut_pool_.add_cut(cut_type_t::MIXED_INTEGER_ROUNDING, cut); } break; } else { // Perform aggregation to try and find a cut @@ -865,24 +765,26 @@ void cut_generation_t::generate_mir_cuts( i_t num_continuous = 0; f_t max_off_bound = 0.0; i_t max_off_bound_var = -1; - for (i_t p = 0; p < inequality.i.size(); p++) { - const i_t j = inequality.i[p]; + for (i_t p = 0; p < inequality.size(); p++) { + const i_t j = inequality.index(p); + const f_t aj = inequality.coeff(p); + if (aj == 0.0) { continue; } if (var_types[j] == variable_type_t::CONTINUOUS) { num_continuous++; - const f_t off_lower = lp.lower[j] > -inf ? xstar[j] - lp.lower[j] : std::abs(xstar[j]); - const f_t off_upper = lp.upper[j] < inf ? lp.upper[j] - xstar[j] : std::abs(xstar[j]); - const f_t off_bound = std::max(off_lower, off_upper); - const i_t col_start = lp.A.col_start[j]; - const i_t col_end = lp.A.col_start[j + 1]; - const i_t col_len = col_end - col_start; + const f_t lb_star_j = complemented_mir.get_lb_star(j); + const f_t ub_star_j = complemented_mir.get_ub_star(j); + const f_t off_lower = lb_star_j > -inf ? xstar[j] - lb_star_j : std::abs(xstar[j]); + const f_t off_upper = ub_star_j < inf ? ub_star_j - xstar[j] : std::abs(xstar[j]); + const f_t off_bound = std::min(off_lower, off_upper); + const i_t col_len = lp.A.col_length(j); if (off_bound > max_off_bound && col_len > 1) { max_off_bound = off_bound; max_off_bound_var = j; } } } - work_estimate += 10 * inequality.i.size(); + work_estimate += 10 * inequality.size(); if (num_continuous == 0 || max_off_bound < 1e-6) { break; } @@ -890,8 +792,8 @@ void cut_generation_t::generate_mir_cuts( if (max_off_bound_var >= 0) { const i_t col_start = lp.A.col_start[max_off_bound_var]; const i_t col_end = lp.A.col_start[max_off_bound_var + 1]; - const i_t col_len = col_end - col_start; - const i_t max_potential_rows = 10; + const i_t col_len = lp.A.col_length(max_off_bound_var); + const i_t max_potential_rows = col_len; if (col_len > 1) { std::vector potential_rows; potential_rows.reserve(col_len); @@ -901,35 +803,48 @@ void cut_generation_t::generate_mir_cuts( const i_t i = lp.A.i[q]; const f_t val = lp.A.x[q]; // Can't use rows that have already been aggregated - if (std::abs(val) > threshold && aggregated_mark[i] == 0) { - potential_rows.push_back(i); - } + if (std::abs(val) > threshold && !aggregated_mark[i]) { potential_rows.push_back(i); } if (potential_rows.size() >= max_potential_rows) { break; } } work_estimate += 5 * (col_end - col_start); - if (!potential_rows.empty()) { - std::sort(potential_rows.begin(), potential_rows.end(), [&](i_t a, i_t b) { - return score[a] > score[b]; - }); - work_estimate += 10 * std::log2(10); - - const i_t pivot_row = potential_rows[0]; - - sparse_vector_t pivot_row_inequality(Arow, pivot_row); - f_t pivot_row_rhs = lp.rhs[pivot_row]; - work_estimate += pivot_row_inequality.i.size(); - mir.combine_rows(lp, - Arow, - max_off_bound_var, - pivot_row_inequality, - pivot_row_rhs, - inequality, - inequality_rhs); + bool did_aggregate = false; + while (!potential_rows.empty()) { + const i_t pivot_row = + *std::max_element(potential_rows.begin(), potential_rows.end(), [&](i_t a, i_t b) { + return scores[a] < scores[b]; + }); + work_estimate += potential_rows.size(); + + inequality_t pivot_row_inequality(Arow, pivot_row, lp.rhs[pivot_row]); + work_estimate += pivot_row_inequality.size(); + // Save inequality before combine_rows mutates it, so we can restore on rejection + inequality_t saved_inequality = inequality; + f_t multiplier = complemented_mir.combine_rows( + lp, Arow, max_off_bound_var, pivot_row_inequality, inequality); + if (max_abs_multiplier / std::abs(multiplier) > 10000 || + std::abs(multiplier) / min_abs_multiplier > 10000) { + printf("Multiplier %e is too small/large max %e min %e\n", + multiplier, + max_abs_multiplier, + min_abs_multiplier); + inequality = saved_inequality; + // Erase the pivot row from the potential rows + potential_rows.erase( + std::remove(potential_rows.begin(), potential_rows.end(), pivot_row), + potential_rows.end()); + continue; + } + max_abs_multiplier = std::max(max_abs_multiplier, std::abs(multiplier)); + min_abs_multiplier = std::min(min_abs_multiplier, std::abs(multiplier)); aggregated_rows.push_back(pivot_row); aggregated_mark[pivot_row] = 1; - work_estimate += inequality.i.size() + pivot_row_inequality.i.size(); - } else { + work_estimate += inequality.size() + pivot_row_inequality.size(); + did_aggregate = true; + break; + } + + if (!did_aggregate) { // No potential rows to aggregate break; } @@ -942,32 +857,26 @@ void cut_generation_t::generate_mir_cuts( if (add_cut) { // We were successful in generating a cut. - // Set the score of the aggregated rows to zero + // Set the score of the aggregated rows to a lower value for (i_t row : aggregated_rows) { - score[row] = 0.0; + scores[row] = 0.99 * scores[row]; + score_queue.push(std::make_pair(scores[row], row)); } + work_estimate += aggregated_rows.size() * std::log2(score_queue.size()); } // Clear the aggregated mark + work_estimate += 2 * aggregated_rows.size(); for (i_t row : aggregated_rows) { aggregated_mark[row] = 0; } - work_estimate += 2 * aggregated_rows.size(); // Clear the aggregated rows aggregated_rows.clear(); // Set the score of the current row to zero - score[i] = 0.0; - - // Re-sort the rows by score - // It's possible this could be made more efficient by storing the rows in a data structure - // that allows us to: - // 1. Get the row with the best score - // 2. Get the row with a nonzero in column j that has the best score - // 3. Remove the rows that have been aggregated - // 4. Remove the current row - best_score_last_permutation(score, sorted_indices); - work_estimate += score.size() * std::log2(score.size()); + scores[i] = 0.0; + score_queue.push(std::make_pair(scores[i], i)); + work_estimate += std::log2(std::max(1, static_cast(score_queue.size()))); } } @@ -984,100 +893,118 @@ void cut_generation_t::generate_gomory_cuts( const std::vector& nonbasic_list) { tableau_equality_t tableau(lp, basis_update, nonbasic_list); - mixed_integer_rounding_cut_t mir(lp, settings, new_slacks, xstar); + mixed_integer_gomory_cut_t gomory_cut; + complemented_mixed_integer_rounding_cut_t complemented_mir(lp, settings, new_slacks); + simplex_solver_settings_t variable_settings = settings; + variable_settings.sub_mip = 1; + variable_bounds_t variable_bounds(lp, variable_settings, var_types, Arow, new_slacks); strong_cg_cut_t cg(lp, var_types, xstar); + std::vector transformed_xstar; + complemented_mir.bound_substitution(lp, variable_bounds, var_types, xstar, transformed_xstar); for (i_t i = 0; i < lp.num_rows; i++) { - sparse_vector_t inequality(lp.num_cols, 0); - f_t inequality_rhs; + inequality_t inequality(lp.num_cols); const i_t j = basic_list[i]; if (var_types[j] != variable_type_t::INTEGER) { continue; } const f_t x_j = xstar[j]; - if (std::abs(x_j - std::round(x_j)) < settings.integer_tol) { continue; } - i_t tableau_status = tableau.generate_base_equality(lp, - settings, - Arow, - var_types, - basis_update, - xstar, - basic_list, - nonbasic_list, - i, - inequality, - inequality_rhs); + if (fractional_part(x_j) < 0.05 || fractional_part(x_j) > 0.95) { continue; } + + i_t tableau_status = tableau.generate_base_equality( + lp, settings, Arow, var_types, basis_update, xstar, basic_list, nonbasic_list, i, inequality); if (tableau_status == 0) { // Generate a CG cut const bool generate_cg_cut = settings.strong_chvatal_gomory_cuts != 0; if (generate_cg_cut) { // Try to generate a CG cut - sparse_vector_t cg_inequality = inequality; - f_t cg_inequality_rhs = inequality_rhs; - if (fractional_part(inequality_rhs) < 0.5) { + inequality_t cg_inequality = inequality; + if (fractional_part(inequality.rhs) < 0.5) { // Multiply by -1 to force the fractional part to be greater than 0.5 - cg_inequality_rhs *= -1; cg_inequality.negate(); } - sparse_vector_t cg_cut(lp.num_cols, 0); - f_t cg_cut_rhs; - i_t cg_status = cg.generate_strong_cg_cut( - lp, settings, var_types, cg_inequality, cg_inequality_rhs, xstar, cg_cut, cg_cut_rhs); - if (cg_status == 0) { cut_pool_.add_cut(cut_type_t::CHVATAL_GOMORY, cg_cut, cg_cut_rhs); } + inequality_t cg_cut(lp.num_cols); + i_t cg_status = + cg.generate_strong_cg_cut(lp, settings, var_types, cg_inequality, xstar, cg_cut); + if (cg_status == 0) { cut_pool_.add_cut(cut_type_t::CHVATAL_GOMORY, cg_cut); } } if (settings.mixed_integer_gomory_cuts == 0) { continue; } - // Given the base inequality, generate a MIR cut - sparse_vector_t cut_A(lp.num_cols, 0); - f_t cut_A_rhs; - i_t mir_status = mir.generate_cut( - inequality, inequality_rhs, lp.upper, lp.lower, var_types, cut_A, cut_A_rhs); - bool A_valid = false; - f_t cut_A_distance = 0.0; - if (mir_status == 0) { - if (cut_A.i.size() == 0) { continue; } - mir.substitute_slacks(lp, Arow, cut_A, cut_A_rhs); - if (cut_A.i.size() == 0) { + // Transform the inequality + inequality_t transformed_inequality = inequality; + complemented_mir.transform_inequality(variable_bounds, var_types, transformed_inequality); + + // Generate a MIR cut from the transformed inequality + inequality_t cut_A_float(lp.num_cols); + bool cut_ok = complemented_mir.generate_cut_nonnegative_maintain_indicies( + transformed_inequality, var_types, cut_A_float); + + // Transform the cut back to the original variables + complemented_mir.untransform_inequality(variable_bounds, var_types, cut_A_float); + complemented_mir.remove_small_coefficients(lp.lower, lp.upper, cut_A_float); + + inequality_t cut_A(lp.num_cols); + if (cut_ok) { cut_ok = gomory_cut.rational_coefficients(var_types, cut_A_float, cut_A); } + + // See if the inequality is violated by the original relaxation solution + f_t cut_A_violation = complemented_mir.compute_violation(cut_A, xstar); + bool A_valid = false; + f_t cut_A_distance = 0.0; + if (cut_ok && cut_A_violation > 1e-6) { + if (cut_A.size() == 0) { continue; } + complemented_mir.substitute_slacks(lp, Arow, cut_A); + complemented_mir.remove_small_coefficients(lp.lower, lp.upper, cut_A); + if (cut_A.size() == 0) { A_valid = false; } else { // Check that the cut is violated - f_t dot = cut_A.dot(xstar); - f_t cut_norm = cut_A.norm2_squared(); - if (dot >= cut_A_rhs) { continue; } - cut_A_distance = (cut_A_rhs - dot) / std::sqrt(cut_norm); + f_t dot = cut_A.vector.dot(xstar); + f_t cut_norm = cut_A.vector.norm2_squared(); + if (dot >= cut_A.rhs) { continue; } + cut_A_distance = (cut_A.rhs - dot) / std::sqrt(cut_norm); A_valid = true; } } // Negate the base inequality inequality.negate(); - inequality_rhs *= -1; - - sparse_vector_t cut_B(lp.num_cols, 0); - f_t cut_B_rhs; - - mir_status = mir.generate_cut( - inequality, inequality_rhs, lp.upper, lp.lower, var_types, cut_B, cut_B_rhs); - bool B_valid = false; - f_t cut_B_distance = 0.0; - if (mir_status == 0) { - if (cut_B.i.size() == 0) { continue; } - mir.substitute_slacks(lp, Arow, cut_B, cut_B_rhs); - if (cut_B.i.size() == 0) { + + inequality_t cut_B_float(lp.num_cols); + + transformed_inequality = inequality; + complemented_mir.transform_inequality(variable_bounds, var_types, transformed_inequality); + + cut_ok = complemented_mir.generate_cut_nonnegative_maintain_indicies( + transformed_inequality, var_types, cut_B_float); + // Transform the cut back to the original variables + complemented_mir.untransform_inequality(variable_bounds, var_types, cut_B_float); + complemented_mir.remove_small_coefficients(lp.lower, lp.upper, cut_B_float); + + inequality_t cut_B(lp.num_cols); + if (cut_ok) { cut_ok = gomory_cut.rational_coefficients(var_types, cut_B_float, cut_B); } + + bool B_valid = false; + f_t cut_B_distance = 0.0; + f_t cut_B_violation = complemented_mir.compute_violation(cut_B, xstar); + if (cut_ok && cut_B_violation > 1e-6) { + if (cut_B.size() == 0) { continue; } + complemented_mir.substitute_slacks(lp, Arow, cut_B); + complemented_mir.remove_small_coefficients(lp.lower, lp.upper, cut_B); + if (cut_B.size() == 0) { B_valid = false; } else { // Check that the cut is violated - f_t dot = cut_B.dot(xstar); - f_t cut_norm = cut_B.norm2_squared(); - if (dot >= cut_B_rhs) { continue; } - cut_B_distance = (cut_B_rhs - dot) / std::sqrt(cut_norm); + f_t dot = cut_B.vector.dot(xstar); + f_t cut_norm = cut_B.vector.norm2_squared(); + if (dot >= cut_B.rhs) { continue; } + cut_B_distance = (cut_B.rhs - dot) / std::sqrt(cut_norm); B_valid = true; } } if ((cut_A_distance > cut_B_distance) && A_valid) { - cut_pool_.add_cut(cut_type_t::MIXED_INTEGER_GOMORY, cut_A, cut_A_rhs); + cut_pool_.add_cut(cut_type_t::MIXED_INTEGER_GOMORY, cut_A); } else if (B_valid) { - cut_pool_.add_cut(cut_type_t::MIXED_INTEGER_GOMORY, cut_B, cut_B_rhs); + cut_pool_.add_cut(cut_type_t::MIXED_INTEGER_GOMORY, cut_B); } } } @@ -1094,8 +1021,7 @@ i_t tableau_equality_t::generate_base_equality( const std::vector& basic_list, const std::vector& nonbasic_list, i_t i, - sparse_vector_t& inequality, - f_t& inequality_rhs) + inequality_t& inequality) { // Let's look for Gomory cuts const i_t j = basic_list[i]; @@ -1233,504 +1159,1155 @@ i_t tableau_equality_t::generate_base_equality( settings_.log.printf("b_bar[%d] = %e\n", i, b_bar[i]); #endif - inequality = a_bar; - inequality_rhs = b_bar_[i]; + inequality.vector = a_bar; + inequality.rhs = b_bar_[i]; return 0; } template -mixed_integer_rounding_cut_t::mixed_integer_rounding_cut_t( - const lp_problem_t& lp, - const simplex_solver_settings_t& settings, - const std::vector& new_slacks, - const std::vector& xstar) - : num_vars_(lp.num_cols), - settings_(settings), - x_workspace_(num_vars_, 0.0), - x_mark_(num_vars_, 0), - has_lower_(num_vars_, 0), - has_upper_(num_vars_, 0), - is_slack_(num_vars_, 0), - slack_rows_(num_vars_, 0), - bound_info_(num_vars_, 0) +bool mixed_integer_gomory_cut_t::rational_approximation(f_t x, + int64_t max_denominator, + int64_t& numerator, + int64_t& denominator) { + int64_t a, p0 = 0, q0 = 1, p1 = 1, q1 = 0; + f_t val = x; + bool negative = false; + + if (x < 0) { + negative = true; + val = -val; + } + + while (1) { + a = (int64_t)std::floor(val); + if (a < 0 || a > INT64_MAX) { return false; } // Protect against overflow + int64_t p2 = a * p1 + p0; + int64_t q2 = a * q1 + q0; + if (q2 > max_denominator) { break; } + p0 = p1; + q0 = q1; + p1 = p2; + q1 = q2; + + f_t rem = val - a; + if (rem < 1e-14) { break; } + val = 1.0 / rem; + } + + numerator = negative ? -p1 : p1; + denominator = q1; + + f_t approx = static_cast(numerator) / static_cast(denominator); + f_t err = std::abs(approx - x); + return err <= 1e-14; +} + +template +bool mixed_integer_gomory_cut_t::rational_coefficients( + const std::vector& var_types, + const inequality_t& input_inequality, + inequality_t& rational_inequality) +{ + rational_inequality = input_inequality; + + std::vector numerators; + std::vector denominators; + std::vector indices; + for (i_t k = 0; k < input_inequality.size(); k++) { + const i_t j = rational_inequality.index(k); + const f_t x = rational_inequality.coeff(k); + if (var_types[j] == variable_type_t::INTEGER) { + int64_t numerator, denominator; + if (!rational_approximation(x, static_cast(1000), numerator, denominator)) { + return false; + } + numerators.push_back(numerator); + denominators.push_back(denominator); + indices.push_back(k); + rational_inequality.vector.x[k] = static_cast(numerator) / static_cast(denominator); + } + } + + int64_t gcd_numerators = gcd(numerators); + int64_t lcm_denominators = lcm(denominators); + + f_t scalar = static_cast(lcm_denominators) / static_cast(gcd_numerators); + if (scalar < 0) { return false; } + if (std::abs(scalar) > 1000) { return false; } + + rational_inequality.scale(scalar); + + return true; +} + +template +int64_t mixed_integer_gomory_cut_t::gcd(const std::vector& integers) +{ + if (integers.empty()) { return 0; } + + int64_t result = integers[0]; + for (size_t i = 1; i < integers.size(); ++i) { + result = std::gcd(result, integers[i]); + } + return result; +} + +template +int64_t mixed_integer_gomory_cut_t::lcm(const std::vector& integers) +{ + if (integers.empty()) { return 0; } + int64_t result = + std::reduce(std::next(integers.begin()), integers.end(), integers[0], [](int64_t a, int64_t b) { + return std::lcm(a, b); + }); + return result; +} + +template +variable_bounds_t::variable_bounds_t(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& var_types, + const csr_matrix_t& Arow, + const std::vector& new_slacks) + : upper_offsets(lp.num_cols + 1, 0), + lower_offsets(lp.num_cols + 1, 0), + upper_activities_(lp.num_rows, 0.0), + lower_activities_(lp.num_rows, 0.0), + num_pos_inf_(lp.num_rows, 0), + num_neg_inf_(lp.num_rows, 0) +{ + if (settings.sub_mip) { + return; // Don't compute the variable upper/lower bounds inside sub-MIP + } + f_t start_time = tic(); + + // Construct the slack map + slack_map_.resize(lp.num_rows, -1); + std::vector slack_coeff(lp.num_rows, 0.0); for (i_t j : new_slacks) { - is_slack_[j] = 1; const i_t col_start = lp.A.col_start[j]; - const i_t i = lp.A.i[col_start]; - slack_rows_[j] = i; - assert(std::abs(lp.A.x[col_start]) == 1.0); + const i_t col_end = lp.A.col_start[j + 1]; + const i_t col_len = col_end - col_start; + assert(col_len == 1); + const i_t i = lp.A.i[col_start]; + slack_map_[i] = j; + slack_coeff[i] = lp.A.x[col_start]; } - needs_complement_ = false; - for (i_t j = 0; j < num_vars_; j++) { - if (lp.lower[j] < 0) { - settings_.log.debug("Variable %d has negative lower bound %e\n", j, lp.lower[j]); + // The constraints are in the form: + // sum_j a_j x_j + sigma * slack = beta + + std::vector num_integer_in_row(lp.num_rows, 0); + // Compute the upper activities of the constraints + for (i_t i = 0; i < lp.num_rows; i++) { + const i_t row_start = Arow.row_start[i]; + const i_t row_end = Arow.row_start[i + 1]; + const i_t slack_index = slack_map_[i]; + f_t activity = 0.0; + for (i_t p = row_start; p < row_end; p++) { + const i_t j = Arow.j[p]; + if (j == slack_index) { continue; } + const f_t aj = Arow.x[p]; + const f_t uj = lp.upper[j]; + const f_t lj = lp.lower[j]; + + if (aj > 0.0) { + if (uj < inf) { + activity += aj * uj; + } else { + num_pos_inf_[i]++; + } + } else { // a_j < 0.0 + if (lj > -inf) { + activity += aj * lj; + } else { + num_pos_inf_[i]++; + } + } + + if (var_types[j] == variable_type_t::INTEGER) { num_integer_in_row[i]++; } } - const f_t uj = lp.upper[j]; - const f_t lj = lp.lower[j]; - const f_t xstar_j = xstar[j]; - if (uj < inf) { - if (uj - xstar_j <= xstar_j - lj) { - has_upper_[j] = 1; - bound_info_[j] = 1; - needs_complement_ = true; - } else if (lj != 0.0) { - has_lower_[j] = 1; - bound_info_[j] = -1; - needs_complement_ = true; + upper_activities_[i] = activity; + } + + // Compute the lower activities of the constraints + for (i_t i = 0; i < lp.num_rows; i++) { + const i_t row_start = Arow.row_start[i]; + const i_t row_end = Arow.row_start[i + 1]; + const i_t slack_index = slack_map_[i]; + f_t activity = 0.0; + for (i_t p = row_start; p < row_end; p++) { + const i_t j = Arow.j[p]; + if (j == slack_index) { continue; } + const f_t aj = Arow.x[p]; + const f_t uj = lp.upper[j]; + const f_t lj = lp.lower[j]; + if (aj > 0.0) { + if (lj > -inf) { + activity += aj * lj; + } else { + num_neg_inf_[i]++; + } + } else { // a_j < 0.0 + if (uj < inf) { + activity += aj * uj; + } else { + num_neg_inf_[i]++; + } } - continue; } + lower_activities_[i] = activity; + } - if (lj > -inf && lj != 0.0) { - has_lower_[j] = 1; - bound_info_[j] = -1; - needs_complement_ = true; + // Now go through all continuous variables and use the activiites to get upper variable bounds + i_t upper_edges = 0; + for (i_t j = 0; j < lp.num_cols; j++) { + upper_offsets[j] = upper_edges; + if (var_types[j] != variable_type_t::CONTINUOUS) { continue; } + const i_t col_start = lp.A.col_start[j]; + const i_t col_end = lp.A.col_start[j + 1]; + for (i_t p = col_start; p < col_end; p++) { + const i_t i = lp.A.i[p]; + if (num_integer_in_row[i] < 1) { continue; } + if (num_neg_inf_[i] > 2 && num_pos_inf_[i] > 2) { continue; } + const i_t row_start = Arow.row_start[i]; + const i_t row_end = Arow.row_start[i + 1]; + const i_t row_len = row_end - row_start; + if (row_len < 2) { continue; } + const f_t a_ij = lp.A.x[p]; + const f_t slack_lower = lp.lower[slack_map_[i]]; + const f_t slack_upper = lp.upper[slack_map_[i]]; + const f_t slack_coeff_i = slack_coeff[i]; + const f_t sigma_slack_lower = slack_coeff_i == 1.0 ? slack_lower : -slack_upper; + const f_t sigma_slack_upper = slack_coeff_i == 1.0 ? slack_upper : -slack_lower; + + if (sigma_slack_lower > -inf) { + const f_t beta = lp.rhs[i] - sigma_slack_lower; + // sum_k a_ik x_k <= beta + + // If we have too many variables in the row that would cause the activity to be infinite, + // we cannot derive an variable bound + if (a_ij > 0.0 && num_neg_inf_[i] <= 2) { + const f_t lower_activity_j = lower_activity(lp.lower[j], lp.upper[j], a_ij); + + // This is inefficient if num_neg_inf_[i] > 0 + // If num_neg_inf_[i] == 1 and var_types[s] != INTEGER, we can't derive a bound + // If num_neg_inf_[i] == 2 and var_types[s ^ j] != INTEGER, we can't derive a bound + // If num_neg_inf_[i] == 2 and var_types[s ^ j] == INTEGER, and lower_activity_j != -inf, + // we can't derive a bound + for (i_t q = row_start; q < row_end; q++) { + const i_t l = Arow.j[q]; + if (var_types[l] == variable_type_t::CONTINUOUS) { continue; } + // sum_{k != l, k != j} a_ik x_k + a_ij x_j + a_il x_l <= beta + // a_ij x_j <= -a_il x_l + beta - sum_{k != l, k != j} a_ik x_k + const f_t a_il = Arow.x[q]; + const f_t lower_activity_l = lower_activity(lp.lower[l], lp.upper[l], a_il); + const f_t sum = adjusted_lower_activity( + lower_activities_[i], num_neg_inf_[i], lower_activity_j, lower_activity_l); + if (sum > -inf) { + // We have a valid variable upper bound + // x_j <= -a_il/a_ij * x_l + beta/a_ij - 1/a_ij * sum_{k != l, k != j} a_ik * + // bound(x_k) + upper_variables.push_back(l); + upper_weights.push_back(-a_il / a_ij); + upper_biases.push_back(beta / a_ij - (1.0 / a_ij) * sum); + upper_edges++; + } + } + } + } + + if (sigma_slack_upper < inf) { + const f_t beta = lp.rhs[i] - sigma_slack_upper; + // sum_k a_ik x_k >= beta + + // If we have too many variables in the row that would cause the activity to be infinite, + // we cannot derive an variable bound + if (a_ij < 0.0 && num_pos_inf_[i] <= 2) { + const f_t upper_activity_j = upper_activity(lp.lower[j], lp.upper[j], a_ij); + + for (i_t q = row_start; q < row_end; q++) { + const i_t l = Arow.j[q]; + if (var_types[l] == variable_type_t::CONTINUOUS) { continue; } + // sum_{k != l, k != j} a_ik x_k + a_ij x_j + a_il x_l >= beta + // a_ij x_j >= -a_il x_l + beta - sum_{k != l, k != j} a_ik x_k + const f_t a_il = Arow.x[q]; + const f_t upper_activity_l = upper_activity(lp.lower[l], lp.upper[l], a_il); + const f_t sum = adjusted_upper_activity( + upper_activities_[i], num_pos_inf_[i], upper_activity_j, upper_activity_l); + if (sum < inf) { + // We have a valid variable upper bound + // x_j <= -a_il/a_ij * x_l + beta/a_ij - 1/a_ij * sum_{k != l, k != j} a_ik * + // bound(x_k) + upper_variables.push_back(l); + upper_weights.push_back(-a_il / a_ij); + upper_biases.push_back(beta / a_ij - (1.0 / a_ij) * sum); + upper_edges++; + } + } + } + } + } + } + upper_offsets[lp.num_cols] = upper_edges; + settings.log.printf("%d variable upper bounds in %.2f seconds\n", upper_edges, toc(start_time)); + + // Now go through all continuous variables and use the activiites to get lower variable bounds + i_t lower_edges = 0; + for (i_t j = 0; j < lp.num_cols; j++) { + lower_offsets[j] = lower_edges; + if (var_types[j] != variable_type_t::CONTINUOUS) { continue; } + const i_t col_start = lp.A.col_start[j]; + const i_t col_end = lp.A.col_start[j + 1]; + for (i_t p = col_start; p < col_end; p++) { + const i_t i = lp.A.i[p]; + if (num_integer_in_row[i] < 1) { continue; } + const i_t row_start = Arow.row_start[i]; + const i_t row_end = Arow.row_start[i + 1]; + const i_t row_len = row_end - row_start; + if (row_len < 2) { continue; } + const f_t a_ij = lp.A.x[p]; + const f_t slack_lower = lp.lower[slack_map_[i]]; + const f_t slack_upper = lp.upper[slack_map_[i]]; + const f_t slack_coeff_i = slack_coeff[i]; + const f_t sigma_slack_lower = slack_coeff_i == 1.0 ? slack_lower : -slack_upper; + const f_t sigma_slack_upper = slack_coeff_i == 1.0 ? slack_upper : -slack_lower; + + if (sigma_slack_lower > -inf) { + const f_t beta = lp.rhs[i] - sigma_slack_lower; + // sum_k a_ik x_k <= beta + + // If we have too many variables in the row that would cause the activity to be infinite, + // we cannot derive a variable bound + if (a_ij < 0.0 && num_neg_inf_[i] <= 2) { + const f_t lower_activity_j = lower_activity(lp.lower[j], lp.upper[j], a_ij); + + for (i_t q = row_start; q < row_end; q++) { + const i_t l = Arow.j[q]; + if (var_types[l] == variable_type_t::CONTINUOUS) { continue; } + // sum_{k != l, k != j} a_ik x_k + a_ij x_j + a_il x_l <= beta + // a_ij x_j <= -a_il x_l + beta - sum_{k != l, k != j} a_ik x_k + // x_j >= -a_il/a_ij * x_l + beta/a_ij - 1/a_ij * sum_{k != l, k != j} a_ik * bound(x_k) + const f_t a_il = Arow.x[q]; + const f_t lower_activity_l = lower_activity(lp.lower[l], lp.upper[l], a_il); + const f_t sum = adjusted_lower_activity( + lower_activities_[i], num_neg_inf_[i], lower_activity_j, lower_activity_l); + if (sum > -inf) { + // We have a valid variable lower bound + // x_j >= -a_il/a_ij * x_l + beta/a_ij - 1/a_ij * sum_{k != l, k != j} a_ik * + // bound(x_k) + lower_variables.push_back(l); + lower_weights.push_back(-a_il / a_ij); + lower_biases.push_back(beta / a_ij - (1.0 / a_ij) * sum); + lower_edges++; + } + } + } + } + + if (sigma_slack_upper < inf) { + const f_t beta = lp.rhs[i] - sigma_slack_upper; + // sum_k a_ik x_k >= beta + + // If we have too many variables in the row that would cause the activity to be infinite, + // we cannot derive a variable bound + if (a_ij > 0.0 && num_pos_inf_[i] <= 2) { + const f_t upper_activity_j = upper_activity(lp.lower[j], lp.upper[j], a_ij); + + for (i_t q = row_start; q < row_end; q++) { + const i_t l = Arow.j[q]; + if (var_types[l] == variable_type_t::CONTINUOUS) { continue; } + // sum_{k != l, k != j} a_ik x_k + a_ij x_j + a_il x_l >= beta + // a_ij x_j >= -a_il x_l + beta - sum_{k != l, k != j} a_ik x_k + const f_t a_il = Arow.x[q]; + const f_t upper_activity_l = upper_activity(lp.lower[l], lp.upper[l], a_il); + const f_t sum = adjusted_upper_activity( + upper_activities_[i], num_pos_inf_[i], upper_activity_j, upper_activity_l); + if (sum < inf) { + // We have a valid variable lower bound + // x_j >= -a_il/a_ij * x_l + beta/a_ij - 1/a_ij * sum_{k != l, k != j} a_ik * + // bound(x_k) + lower_variables.push_back(l); + lower_weights.push_back(-a_il / a_ij); + lower_biases.push_back(beta / a_ij - (1.0 / a_ij) * sum); + lower_edges++; + } + } + } + } } } + lower_offsets[lp.num_cols] = lower_edges; + settings.log.printf("%d variable lower bounds in %.2f seconds\n", lower_edges, toc(start_time)); } template -void mixed_integer_rounding_cut_t::to_nonnegative(const lp_problem_t& lp, - sparse_vector_t& inequality, - f_t& rhs) +complemented_mixed_integer_rounding_cut_t::complemented_mixed_integer_rounding_cut_t( + const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& new_slacks) + : is_slack_(lp.num_cols, 0), + slack_rows_(lp.num_cols, -1), + slack_cols_(lp.num_rows, -1), + lb_variable_(lp.num_cols, -1), + lb_star_(lp.num_cols, 0.0), + ub_variable_(lp.num_cols, -1), + ub_star_(lp.num_cols, 0.0), + transformed_upper_(lp.num_cols, inf), + bound_changed_(lp.num_cols, 0), + scratch_pad_(lp.num_cols) { - const i_t nz = inequality.i.size(); - for (i_t k = 0; k < nz; k++) { - const i_t j = inequality.i[k]; - const f_t aj = inequality.x[k]; - if (bound_info_[j] == -1) { - // v_j = x_j - l_j, v_j >= 0 - // x_j = v_j + l_j - // sum_{k != j} a_k x_j + a_j x_j <= beta - // sum_{k != j} a_k x_j + a_j (v_j + l_j) <= beta - // sum_{k != j} a_k x_j + a_j v_j <= beta - a_j l_j - const f_t lj = lp.lower[j]; - rhs -= aj * lj; - } else if (bound_info_[j] == 1) { - // w_j = u_j - x_j, w_j >= 0 - // x_j = u_j - w_j - // sum_{k != j} a_k x_k + a_j x_j <= beta - // sum_{k != j} a_k x_k + a_j (u_j - w_j) <= beta - // sum_{k != j} a_k x_k - a_j w_j <= beta - a_j u_j - const f_t uj = lp.upper[j]; - inequality.x[k] *= -1.0; - rhs -= aj * uj; - } + for (i_t j : new_slacks) { + is_slack_[j] = 1; + const i_t col_start = lp.A.col_start[j]; + const i_t i = lp.A.i[col_start]; + slack_rows_[j] = i; + slack_cols_[i] = j; + assert(std::abs(lp.A.x[col_start]) == 1.0); } } template -void mixed_integer_rounding_cut_t::relaxation_to_nonnegative( +void complemented_mixed_integer_rounding_cut_t::compute_initial_scores_for_rows( const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const csr_matrix_t& Arow, const std::vector& xstar, - std::vector& xstar_nonnegative) + const std::vector& ystar, + std::vector& scores) { - xstar_nonnegative = xstar; - const i_t n = lp.num_cols; - for (i_t j = 0; j < n; ++j) { - if (bound_info_[j] == -1) { - // v_j = x_j - l_j - const f_t lj = lp.lower[j]; - xstar_nonnegative[j] -= lj; - } else if (bound_info_[j] == 1) { - // w_j = u_j - x_j - const f_t uj = lp.upper[j]; - xstar_nonnegative[j] = uj - xstar_nonnegative[j]; + const bool verbose = false; + const i_t n = lp.num_cols; + const f_t obj_norm = vector_norm2(lp.objective); + const f_t obj_denom = std::max(1.0, obj_norm); + + // Compute initial scores for all rows + scores.resize(lp.num_rows, 0.0); + for (i_t i = 0; i < lp.num_rows; i++) { + const i_t row_start = Arow.row_start[i]; + const i_t row_end = Arow.row_start[i + 1]; + + const i_t row_nz = row_end - row_start; + f_t row_norm = 0.0; + for (i_t p = row_start; p < row_end; p++) { + const f_t a_j = Arow.x[p]; + row_norm += a_j * a_j; + } + row_norm = std::sqrt(row_norm); + + const f_t density = static_cast(row_nz) / static_cast(n); + const f_t dual = std::abs(ystar[i]); + + const i_t slack = slack_cols_[i]; + assert(slack >= 0); + const f_t slack_value = std::max(xstar[slack], 0.0); + const f_t slack_denom = std::max(0.1, std::sqrt(row_norm)); + + const f_t nz_weight = 0.0001; + const f_t dual_weight = 1.0; + const f_t slack_weight = 0.001; + + scores[i] = nz_weight * (1.0 - density) + dual_weight * std::max(dual / obj_denom, 0.0001) + + slack_weight * (1.0 - slack_value / slack_denom); + + if (verbose) { + settings.log.printf("Scores[%d] = %e density %.2f dual %e slack %e\n", + i, + scores[i], + density, + dual, + slack_value); } } } template -void mixed_integer_rounding_cut_t::to_original(const lp_problem_t& lp, - sparse_vector_t& inequality, - f_t& rhs) +bool complemented_mixed_integer_rounding_cut_t::cut_generation_heuristic( + const inequality_t& transformed_inequality, + const std::vector& var_types, + const std::vector& transformed_xstar, + inequality_t& transformed_cut, + f_t& work_estimate) { - const i_t nz = inequality.i.size(); - for (i_t k = 0; k < nz; k++) { - const i_t j = inequality.i[k]; - const f_t dj = inequality.x[k]; - if (bound_info_[j] == -1) { - // v_j = x_j - l_j, v_j >= 0 - // sum_{k != j} d_k x_k + d_j v_j >= beta - // sum_{k != j} d_k x_k + d_j (x_j - l_j) >= beta - // sum_{k != j} d_k x_k + d_j x_j >= beta + d_j l_j - const f_t lj = lp.lower[j]; - rhs += dj * lj; - } else if (bound_info_[j] == 1) { - // w_j = u_j - x_j, w_j >= 0 - // sum_{k != j} d_k x_k + d_j w_j >= beta - // sum_{k != j} d_k x_k + d_j (u_j - x_j) >= beta - // sum_{k != j} d_k x_k - d_j x_j >= beta - d_j u_j - const f_t uj = lp.upper[j]; - inequality.x[k] *= -1.0; - rhs -= dj * uj; + std::vector deltas_to_try; + deltas_to_try.reserve(transformed_inequality.size()); + deltas_to_try.push_back(1.0); + work_estimate += transformed_inequality.size(); + i_t num_integers = 0; + f_t max_coeff = 0.0; + for (i_t k = 0; k < transformed_inequality.size(); k++) { + const i_t j = transformed_inequality.index(k); + const f_t abs_aj = std::abs(transformed_inequality.coeff(k)); + if (var_types[j] == variable_type_t::INTEGER) { + num_integers++; + max_coeff = std::max(max_coeff, abs_aj); + const f_t x_j = transformed_xstar[j]; + const f_t new_upper_j = new_upper(j); + const f_t dist_upper = new_upper_j - x_j; + const f_t dist_lower = x_j; + const bool between_bounds = x_j > 1e-6 && (new_upper_j == inf || dist_upper > 0.0); + if (between_bounds) { deltas_to_try.push_back(abs_aj); } + } + } + if (max_coeff > 1e-6 && max_coeff != 1.0) { + deltas_to_try.push_back(max_coeff); + deltas_to_try.push_back(max_coeff + 1.0); + } + + std::vector complemented_indices; + complemented_indices.reserve(num_integers); + std::vector distance_from_midpoint; + distance_from_midpoint.reserve(num_integers); + std::vector integer_indices; + integer_indices.reserve(num_integers); + for (i_t k = 0; k < transformed_inequality.size(); k++) { + const i_t j = transformed_inequality.index(k); + if (var_types[j] == variable_type_t::INTEGER && new_upper(j) < inf) { + const f_t x_j = transformed_xstar[j]; + const f_t new_upper_j = new_upper(j); + if (x_j > 1e-6 && new_upper_j < inf) { + const f_t midpoint_j = new_upper_j / 2.0; + distance_from_midpoint.push_back(x_j - midpoint_j); + integer_indices.push_back(k); + } + } + } + + std::vector perm(integer_indices.size()); + best_score_first_permutation(distance_from_midpoint, perm); + work_estimate += + integer_indices.size() > 0 ? integer_indices.size() * std::log2(integer_indices.size()) : 0; + + bool cut_found = false; + + inequality_t complemented_inequality = transformed_inequality; + work_estimate += 4 * transformed_inequality.size(); + + f_t delta = 0.0; + f_t best_violation = 0.0; + + // First try without any complementation + for (const f_t tmp_delta : deltas_to_try) { + bool cut_ok = scale_uncomplement_and_generate_cut(var_types, + transformed_xstar, + complemented_indices, + complemented_inequality, + tmp_delta, + transformed_cut, + work_estimate); + if (!cut_ok) { continue; } + // Check if the cut is violated + best_violation = compute_violation(transformed_cut, transformed_xstar); + work_estimate += 4 * transformed_cut.size(); + if (best_violation > 1e-6) { + cut_found = true; + delta = tmp_delta; + break; } } + + if (!cut_found) { + // Complement an integer variable + for (const i_t idx : perm) { + const i_t l = integer_indices[idx]; + const i_t j = complemented_inequality.index(l); + // We have an integer variable x_j <= b_j + // We create a new variable xbar_j such that + // x_j + xbar_j = b_j + // x_j = b_j - xbar_j, xbar_j = b_j - x_j + // + // The inequality + // sum_{k != j} a_k x_k + a_j x_j >= beta + // becomes + // sum_{k != j} a_k x_k + a_j (b_j - xbar_j) >= beta + // sum_{k != j} a_k x_k - a_j xbar_j >= beta - a_j b_j + const f_t b_j = new_upper(j); + const f_t a_j = complemented_inequality.coeff(l); + + complemented_inequality.vector.x[l] = -a_j; + complemented_inequality.rhs -= a_j * b_j; + complemented_indices.push_back(l); + + for (const f_t tmp_delta : deltas_to_try) { + bool cut_ok = scale_uncomplement_and_generate_cut(var_types, + transformed_xstar, + complemented_indices, + complemented_inequality, + tmp_delta, + transformed_cut, + work_estimate); + if (!cut_ok) { continue; } + // Check if the cut is violated + best_violation = compute_violation(transformed_cut, transformed_xstar); + work_estimate += 4 * transformed_cut.size(); + if (best_violation > 1e-6) { + cut_found = true; + delta = tmp_delta; + break; + } + } + if (cut_found) { break; } + } + } + + if (!cut_found) { return false; } + + // We have found a cut. Now try to improve the violation by scaling the cut by 1/2, 1/4, 1/8, etc. + std::vector scaled_deltas_to_try = {delta / 2.0, delta / 4.0, delta / 8.0}; + for (const f_t tmp_delta : scaled_deltas_to_try) { + inequality_t tmp_cut_delta; + bool cut_ok = scale_uncomplement_and_generate_cut(var_types, + transformed_xstar, + complemented_indices, + complemented_inequality, + tmp_delta, + tmp_cut_delta, + work_estimate); + if (!cut_ok) { continue; } + + // Check if the cut is violated + f_t violation = compute_violation(tmp_cut_delta, transformed_xstar); + work_estimate += 4 * tmp_cut_delta.size(); + if (violation > best_violation) { + best_violation = violation; + transformed_cut = tmp_cut_delta; + delta = tmp_delta; + } + } + + std::vector best_complemented_indices = complemented_indices; + work_estimate += 2 * best_complemented_indices.size(); + + // Try to improve the violation by complementing integer variables + complemented_inequality = transformed_inequality; + work_estimate += 4 * transformed_inequality.size(); + complemented_indices.clear(); + for (const i_t idx : perm) { + const i_t l = integer_indices[idx]; + const i_t j = complemented_inequality.index(l); + // We have an integer variable x_j <= b_j + // We create a new variable xbar_j such that + // x_j + xbar_j = b_j + // x_j = b_j - xbar_j, xbar_j = b_j - x_j + // + // The inequality + // sum_{k != j} a_k x_k + a_j x_j >= beta + // becomes + // sum_{k != j} a_k x_k + a_j (b_j - xbar_j) >= beta + // sum_{k != j} a_k x_k - a_j xbar_j >= beta - a_j b_j + const f_t b_j = new_upper(j); + const f_t a_j = complemented_inequality.coeff(l); + + complemented_inequality.vector.x[l] = -a_j; + complemented_inequality.rhs -= a_j * b_j; + complemented_indices.push_back(l); + + inequality_t tmp_cut_delta; + + bool cut_ok = scale_uncomplement_and_generate_cut(var_types, + transformed_xstar, + complemented_indices, + complemented_inequality, + delta, + tmp_cut_delta, + work_estimate); + if (!cut_ok) { continue; } + // Check if the cut is violated + f_t violation = compute_violation(tmp_cut_delta, transformed_xstar); + work_estimate += 4 * tmp_cut_delta.size(); + if (violation > best_violation) { + best_violation = violation; + best_complemented_indices = complemented_indices; + transformed_cut = tmp_cut_delta; + } + } + + return best_violation > 1e-6; +} + +template +bool complemented_mixed_integer_rounding_cut_t::scale_uncomplement_and_generate_cut( + const std::vector& var_types, + const std::vector& transformed_xstar, + const std::vector& complemented_indices, + const inequality_t& complemented_inequality, + f_t delta, + inequality_t& cut_delta, + f_t& work_estimate) +{ + inequality_t scaled_inequality = complemented_inequality; + if (delta != 1.0) { scaled_inequality.scale(1.0 / delta); } + bool cut_ok = generate_cut_nonnegative_maintain_indicies(scaled_inequality, var_types, cut_delta); + if (!cut_ok) { return false; } + work_estimate += 4 * scaled_inequality.size(); + + // Now we need to transform the complemented variables back + for (i_t h = 0; h < complemented_indices.size(); h++) { + const i_t l = complemented_indices[h]; + const i_t j = complemented_inequality.index(l); + // Our cut is of the form + // sum_{k != j} d_k x_k + d_j xbar_j >= alpha + // we have that xbar_j = b_j - x_j + // So + // sum_{k != j} d_k x_k + d_j (b_j - x_j) >= alpha + // Or + // sum_{k != j} d_k x_k - d_j x_j >= alpha - d_j b_j + + const f_t b_j = new_upper(j); + const f_t d_j = cut_delta.coeff(l); + cut_delta.vector.x[l] = -d_j; + cut_delta.rhs -= d_j * b_j; + } + work_estimate += 5 * complemented_indices.size(); + return true; } template -void mixed_integer_rounding_cut_t::remove_small_coefficients( +void complemented_mixed_integer_rounding_cut_t::remove_small_coefficients( const std::vector& lower_bounds, const std::vector& upper_bounds, - sparse_vector_t& cut, - f_t& cut_rhs) + inequality_t& cut) { - const i_t nz = cut.i.size(); + const i_t nz = cut.size(); i_t removed = 0; - for (i_t k = 0; k < cut.i.size(); k++) { - const i_t j = cut.i[k]; + for (i_t k = 0; k < cut.size(); k++) { + const i_t j = cut.index(k); // Check for small coefficients - const f_t aj = cut.x[k]; + const f_t aj = cut.coeff(k); if (std::abs(aj) < 1e-6) { if (aj >= 0.0 && upper_bounds[j] < inf) { // Move this to the right-hand side - cut_rhs -= aj * upper_bounds[j]; - cut.x[k] = 0.0; + cut.rhs -= aj * upper_bounds[j]; + cut.vector.x[k] = 0.0; removed++; } else if (aj <= 0.0 && lower_bounds[j] > -inf) { - cut_rhs += aj * lower_bounds[j]; - cut.x[k] = 0.0; + cut.rhs -= aj * lower_bounds[j]; + cut.vector.x[k] = 0.0; removed++; continue; } else { + // We need to keep the coefficient } } } if (removed > 0) { - sparse_vector_t new_cut(cut.n, 0); + inequality_t new_cut(cut.vector.n); cut.squeeze(new_cut); cut = new_cut; } } template -i_t mixed_integer_rounding_cut_t::generate_cut_nonnegative( - const sparse_vector_t& a, - f_t beta, +void complemented_mixed_integer_rounding_cut_t::bound_substitution( + const lp_problem_t& lp, + const variable_bounds_t& variable_bounds, const std::vector& var_types, - sparse_vector_t& cut, - f_t& cut_rhs) + const std::vector& xstar, + std::vector& transformed_xstar) { - auto f = [](f_t q_1, f_t q_2) -> f_t { - f_t q_1_hat = q_1 - std::floor(q_1); - f_t q_2_hat = q_2 - std::floor(q_2); - return std::min(q_1_hat, q_2_hat) + q_2_hat * std::floor(q_1); - }; - - auto h = [](f_t q) -> f_t { return std::max(q, 0.0); }; + transformed_xstar.resize(lp.num_cols); + // Perform bound substitution for continuous variables + for (i_t j = 0; j < lp.num_cols; j++) { + if (var_types[j] != variable_type_t::CONTINUOUS) { continue; } + // Step 1: Decide whether to use variable or simple bounds + const f_t uj = lp.upper[j]; + const f_t lj = lp.lower[j]; + const f_t xstar_j = xstar[j]; - std::vector cut_indices; - cut_indices.reserve(a.i.size()); - f_t R = (beta - std::floor(beta)) * std::ceil(beta); - - for (i_t k = 0; k < a.i.size(); k++) { - const i_t jj = a.i[k]; - f_t aj = a.x[k]; - if (var_types[jj] == variable_type_t::INTEGER) { - x_workspace_[jj] += f(aj, beta); - if (!x_mark_[jj] && x_workspace_[jj] != 0.0) { - x_mark_[jj] = 1; - cut_indices.push_back(jj); - } - } else { - x_workspace_[jj] += h(aj); - if (!x_mark_[jj] && x_workspace_[jj] != 0.0) { - x_mark_[jj] = 1; - cut_indices.push_back(jj); + // Set lb and lb_star to the simple lower bound + lb_variable_[j] = -1; + lb_star_[j] = lj; + + // Set ub and ub_star to the simple upper bound + ub_variable_[j] = -1; + ub_star_[j] = uj; + + // Check the variable lower bound and update lb and lb_star + // if these yield a tighter bound + const i_t lower_variable_start = variable_bounds.lower_offsets[j]; + const i_t lower_variable_end = variable_bounds.lower_offsets[j + 1]; + for (i_t p = lower_variable_start; p < lower_variable_end; p++) { + const i_t i = variable_bounds.lower_variables[p]; + const f_t gamma = variable_bounds.lower_weights[p]; + const f_t alpha = variable_bounds.lower_biases[p]; + // x_j >= gamma * x_i + alpha + + const f_t xstar_i = xstar[i]; + const f_t val = gamma * xstar_i + alpha; + if (val > lb_star_[j]) { + lb_variable_[j] = p; + lb_star_[j] = val; } } - } - - cut.i.reserve(cut_indices.size()); - cut.x.reserve(cut_indices.size()); - cut.i.clear(); - cut.x.clear(); - for (i_t k = 0; k < cut_indices.size(); k++) { - const i_t j = cut_indices[k]; - cut.i.push_back(j); - cut.x.push_back(x_workspace_[j]); - } - // Clear the workspace - for (i_t jj : cut_indices) { - x_workspace_[jj] = 0.0; - x_mark_[jj] = 0; - } - -#ifdef CHECK_WORKSPACE - for (i_t j = 0; j < x_workspace_.size(); j++) { - if (x_workspace_[j] != 0.0) { - printf("After generate_cut: Dirty x_workspace_[%d] = %e\n", j, x_workspace_[j]); - assert(x_workspace_[j] == 0.0); + // Check the variable upper bound and update ub and ub_star + // if these yield a tighter bound + const i_t upper_variable_start = variable_bounds.upper_offsets[j]; + const i_t upper_variable_end = variable_bounds.upper_offsets[j + 1]; + for (i_t p = upper_variable_start; p < upper_variable_end; p++) { + const i_t i = variable_bounds.upper_variables[p]; + const f_t gamma = variable_bounds.upper_weights[p]; + const f_t alpha = variable_bounds.upper_biases[p]; + // x_j <= gamma * x_i + alpha + + const f_t xstar_i = xstar[i]; + const f_t val = gamma * xstar_i + alpha; + if (val < ub_star_[j]) { + ub_variable_[j] = p; + ub_star_[j] = val; + } } - if (x_mark_[j] != 0) { - printf("After generate_cut: Dirty x_mark_[%d] = %d\n", j, x_mark_[j]); - assert(x_mark_[j] == 0); + + // Step 2: Decide to use the lower or upper bound + const bool has_finite_lower_bound = lb_star_[j] > -inf; + const bool has_finite_upper_bound = ub_star_[j] < inf; + if (!has_finite_lower_bound && !has_finite_upper_bound) { + transformed_xstar[j] = xstar_j; + transformed_upper_[j] = inf; + bound_changed_[j] = 0; + continue; } - } -#endif + if (has_finite_lower_bound && + (!has_finite_upper_bound || (xstar_j - lb_star_[j] <= ub_star_[j] - xstar_j))) { + // Use the lower bound + // lb_star_j <= x_j <= ub_star_j + // v_j = x_j - lb_star_j, + // 0 <= v_j <= ub_star_j - lb_star_j + transformed_upper_[j] = ub_star_[j] - lb_star_[j]; + transformed_xstar[j] = xstar_j - lb_star_[j]; + bound_changed_[j] = (lb_star_[j] == 0.0) ? 0 : -1; + } else if (has_finite_upper_bound) { + // Use the upper bound + // lb_star_j <= x_j <= ub_star_j + // x_j + w_j = ub_star_j, + // w_j = ub_star_j - x_j, + // x_j = ub_star_j - w_j + // 0 <= w_j <= ub_star_j - lb_star_j + transformed_upper_[j] = ub_star_[j] - lb_star_[j]; + transformed_xstar[j] = ub_star_[j] - xstar_j; + bound_changed_[j] = 1; + } + } + + // Perform bound substitution for the integer variables + for (i_t j = 0; j < lp.num_cols; j++) { + if (var_types[j] != variable_type_t::INTEGER) { continue; } + const f_t uj = lp.upper[j]; + const f_t lj = lp.lower[j]; + const f_t xstar_j = xstar[j]; - // The new cut is: g'*x >= R - // But we want to have it in the form h'*x <= b - cut.sort(); + lb_star_[j] = lj; + ub_star_[j] = uj; + transformed_xstar[j] = xstar_j; + transformed_upper_[j] = uj; + bound_changed_[j] = 0; - cut_rhs = R; + if (uj < inf) { + if (uj - xstar_j <= xstar_j - lj) { + // Use the upper bound + // lj <= x_j <= uj + // x_j + w_j = uj, + // w_j = uj - x_j, + // x_j = uj - w_j + // 0 <= w_j <= uj - lj + transformed_upper_[j] = uj - lj; + transformed_xstar[j] = uj - xstar_j; + bound_changed_[j] = 1; + } else if (lj != 0.0) { + // Use the lower bound + // lj <= x_j <= uj + // v_j = x_j - lj, + // 0 <= v_j <= uj - lj + transformed_upper_[j] = uj - lj; + transformed_xstar[j] = xstar_j - lj; + bound_changed_[j] = -1; + } + continue; + } -#ifdef CHECK_REPEATED_INDICES - // Check for repeated indicies - std::vector check(num_vars_, 0); - for (i_t p = 0; p < cut.i.size(); p++) { - if (check[cut.i[p]] != 0) { - printf("repeated index in generated cut\n"); - assert(check[cut.i[p]] == 0); + if (lj > -inf && lj != 0.0) { + // Use the lower bound + // lj <= x_j <= uj + // v_j = x_j - lj, + // 0 <= v_j <= uj - lj + transformed_upper_[j] = uj - lj; + transformed_xstar[j] = xstar_j - lj; + bound_changed_[j] = -1; } - check[cut.i[p]] = 1; } -#endif - - if (cut.i.size() == 0) { return -1; } - - return 0; } template -i_t mixed_integer_rounding_cut_t::generate_cut( - const sparse_vector_t& a, - f_t beta, - const std::vector& upper_bounds, - const std::vector& lower_bounds, - const std::vector& var_types, - sparse_vector_t& cut, - f_t& cut_rhs) +void complemented_mixed_integer_rounding_cut_t::transform_inequality( + const variable_bounds_t& variable_bounds, + const std::vector& var_type, + inequality_t& inequality) { -#ifdef CHECK_WORKSPACE - for (i_t j = 0; j < x_workspace_.size(); j++) { - if (x_workspace_[j] != 0.0) { - printf("Before generate_cut: Dirty x_workspace_[%d] = %e\n", j, x_workspace_[j]); - printf("num_vars_ %d\n", num_vars_); - printf("x_workspace_.size() %ld\n", x_workspace_.size()); - assert(x_workspace_[j] == 0.0); - } - if (x_mark_[j] != 0) { - printf("Before generate_cut: Dirty x_mark_[%d] = %d\n", j, x_mark_[j]); - assert(x_mark_[j] == 0); + const i_t nz = inequality.size(); + for (i_t k = 0; k < nz; k++) { + const i_t j = inequality.index(k); + const f_t aj = inequality.coeff(k); + if (var_type[j] != variable_type_t::CONTINUOUS) { + scratch_pad_.add_to_pad(j, aj); + continue; } - } -#endif - - auto f = [](f_t q_1, f_t q_2) -> f_t { - f_t q_1_hat = q_1 - std::floor(q_1); - f_t q_2_hat = q_2 - std::floor(q_2); - return std::min(q_1_hat, q_2_hat) + q_2_hat * std::floor(q_1); - }; - - auto h = [](f_t q) -> f_t { return std::max(q, 0.0); }; - - std::vector cut_indices; - cut_indices.reserve(a.i.size()); - f_t R; - if (!needs_complement_) { - R = (beta - std::floor(beta)) * std::ceil(beta); - - for (i_t k = 0; k < a.i.size(); k++) { - const i_t jj = a.i[k]; - f_t aj = a.x[k]; - if (var_types[jj] == variable_type_t::INTEGER) { - x_workspace_[jj] += f(aj, beta); - if (!x_mark_[jj] && x_workspace_[jj] != 0.0) { - x_mark_[jj] = 1; - cut_indices.push_back(jj); - } + if (bound_changed_[j] == -1) { + if (lb_variable_[j] == -1) { + // v_j = x_j - l_j, v_j >= 0 + // x_j = v_j + l_j + // sum_{k != j} a_k x_k + a_j x_j >= beta + // sum_{k != j} a_k x_k + a_j (v_j + l_j) >= beta + // sum_{k != j} a_k x_k + a_j v_j >= beta - a_j l_j + const f_t lj = lb_star_[j]; + inequality.rhs -= aj * lj; + scratch_pad_.add_to_pad(j, aj); } else { - x_workspace_[jj] += h(aj); - if (!x_mark_[jj] && x_workspace_[jj] != 0.0) { - x_mark_[jj] = 1; - cut_indices.push_back(jj); - } - } - } - } else { - // Compute r - f_t r = beta; - for (i_t k = 0; k < a.i.size(); k++) { - const i_t jj = a.i[k]; - if (has_upper_[jj]) { - const f_t uj = upper_bounds[jj]; - r -= uj * a.x[k]; - continue; - } - if (has_lower_[jj]) { - const f_t lj = lower_bounds[jj]; - r -= lj * a.x[k]; + // v_j = x_j - lb*_j, v_j >= 0 + // x_j = v_j + lb*_j + // lb*_j = gamma * x_i + alpha + // x_j = v_j + gamma * x_i + alpha + // sum_{k != j} a_k x_k + a_j x_j >= beta + // sum_{k != j} a_k x_k + a_j (v_j + gamma * x_i + alpha) >= beta + // sum_{k != j} a_k x_k + a_j v_j + a_j * gamma * x_i >= beta - a_j alpha + const i_t p = lb_variable_[j]; + const f_t alpha = variable_bounds.lower_biases[p]; + const f_t gamma = variable_bounds.lower_weights[p]; + const i_t i = variable_bounds.lower_variables[p]; + inequality.rhs -= aj * alpha; + scratch_pad_.add_to_pad(j, aj); + scratch_pad_.add_to_pad(i, aj * gamma); } - } - - // Compute R - R = std::ceil(r) * (r - std::floor(r)); - for (i_t k = 0; k < a.i.size(); k++) { - const i_t jj = a.i[k]; - const f_t aj = a.x[k]; - if (has_upper_[jj]) { - const f_t uj = upper_bounds[jj]; - if (var_types[jj] == variable_type_t::INTEGER) { - R -= f(-aj, r) * uj; - } else { - R -= h(-aj) * uj; - } - } else if (has_lower_[jj]) { - const f_t lj = lower_bounds[jj]; - if (var_types[jj] == variable_type_t::INTEGER) { - R += f(aj, r) * lj; - } else { - R += h(aj) * lj; - } - } - } - - // Compute the cut coefficients - for (i_t k = 0; k < a.i.size(); k++) { - const i_t jj = a.i[k]; - const f_t aj = a.x[k]; - if (has_upper_[jj]) { - if (var_types[jj] == variable_type_t::INTEGER) { - // Upper intersect I - x_workspace_[jj] -= f(-aj, r); - if (!x_mark_[jj] && x_workspace_[jj] != 0.0) { - x_mark_[jj] = 1; - cut_indices.push_back(jj); - } - } else { - // Upper intersect C - f_t h_j = h(-aj); - if (h_j != 0.0) { - x_workspace_[jj] -= h_j; - if (!x_mark_[jj]) { - x_mark_[jj] = 1; - cut_indices.push_back(jj); - } - } - } - } else if (var_types[jj] == variable_type_t::INTEGER) { - // I \ Upper - x_workspace_[jj] += f(aj, r); - if (!x_mark_[jj] && x_workspace_[jj] != 0.0) { - x_mark_[jj] = 1; - cut_indices.push_back(jj); - } + } else if (bound_changed_[j] == 1) { + if (ub_variable_[j] == -1) { + // w_j = u_j - x_j, w_j >= 0 + // x_j = u_j - w_j + // sum_{k != j} a_k x_k + a_j x_j >= beta + // sum_{k != j} a_k x_k + a_j (u_j - w_j) >= beta + // sum_{k != j} a_k x_k - a_j w_j >= beta - a_j u_j + const f_t uj = ub_star_[j]; + inequality.rhs -= aj * uj; + scratch_pad_.add_to_pad(j, -aj); } else { - // C \ Upper - f_t h_j = h(aj); - if (h_j != 0.0) { - x_workspace_[jj] += h_j; - if (!x_mark_[jj]) { - x_mark_[jj] = 1; - cut_indices.push_back(jj); - } - } + // w_j = ub*_j - x_j, w_j >= 0 + // x_j = ub*_j - w_j + // ub*_j = gamma * x_i + alpha + // x_j = gamma * x_i + alpha - w_j + // sum_{k != j} a_k x_k + a_j x_j >= beta + // sum_{k != j} a_k x_k + a_j (ub*_j - w_j) >= beta + // sum_{k != j} a_k x_k + a_j (gamma * x_i + alpha) - a_j w_j >= beta + // sum_{k != j} a_k x_k + a_j gamma * x_i - a_j w_j >= beta - a_j alpha + const i_t p = ub_variable_[j]; + const f_t alpha = variable_bounds.upper_biases[p]; + const f_t gamma = variable_bounds.upper_weights[p]; + const i_t i = variable_bounds.upper_variables[p]; + inequality.rhs -= aj * alpha; + scratch_pad_.add_to_pad(j, -aj); + scratch_pad_.add_to_pad(i, aj * gamma); } + } else if (bound_changed_[j] == 0) { + scratch_pad_.add_to_pad(j, aj); } } + scratch_pad_.get_pad(inequality.vector.i, inequality.vector.x); + // At this point we have converted all the continuous variables to be nonnegative + // Note that since continuous variables had VUB or VLB, they modified + // the integer variables. - cut.i.reserve(cut_indices.size()); - cut.x.reserve(cut_indices.size()); - cut.i.clear(); - cut.x.clear(); - for (i_t k = 0; k < cut_indices.size(); k++) { - const i_t jj = cut_indices[k]; + // We clear the scratch pad. As it is no longer needed. + scratch_pad_.clear_pad(); - // Check for small coefficients - const f_t aj = x_workspace_[jj]; - if (std::abs(aj) < 1e-6) { - if (aj >= 0.0 && upper_bounds[jj] < inf) { - // Move this to the right-hand side - R -= aj * upper_bounds[jj]; - continue; - } else if (aj <= 0.0 && lower_bounds[jj] > -inf) { - R += aj * lower_bounds[jj]; - continue; - } else { - } + // We now convert all the integer variables to be nonnegative + const i_t nz_after = inequality.size(); + for (i_t k = 0; k < nz_after; k++) { + const i_t j = inequality.index(k); + if (var_type[j] != variable_type_t::INTEGER) { continue; } + const f_t aj = inequality.coeff(k); + if (bound_changed_[j] == -1) { + // v_j = x_j - l_j, v_j >= 0 + // x_j = v_j + l_j + // sum_{k != j} a_k x_k + a_j x_j >= beta + // sum_{k != j} a_k x_k + a_j (v_j + l_j) >= beta + // sum_{k != j} a_k x_k + a_j v_j >= beta - a_j l_j + const f_t lj = lb_star_[j]; + inequality.rhs -= aj * lj; + } else if (bound_changed_[j] == 1) { + // w_j = u_j - x_j, w_j >= 0 + // x_j = u_j - w_j + // sum_{k != j} a_k x_j + a_j x_j >= beta + // sum_{k != j} a_k x_j + a_j (u_j - w_j) >= beta + // sum_{k != j} a_k x_j - a_j w_j >= beta - a_j u_j + const f_t uj = ub_star_[j]; + inequality.rhs -= aj * uj; + inequality.vector.x[k] *= -1.0; } - cut.i.push_back(jj); - cut.x.push_back(x_workspace_[jj]); } +} - // Clear the workspace - for (i_t jj : cut_indices) { - x_workspace_[jj] = 0.0; - x_mark_[jj] = 0; +template +void complemented_mixed_integer_rounding_cut_t::untransform_inequality( + const variable_bounds_t& variable_bounds, + const std::vector& var_type, + inequality_t& inequality) +{ + // First convert all the integers variables back to their original form: l_j <= x_j <= u_j + const i_t nz = inequality.size(); + for (i_t k = 0; k < nz; k++) { + const i_t j = inequality.index(k); + if (var_type[j] != variable_type_t::INTEGER) { continue; } + const f_t dj = inequality.coeff(k); + if (bound_changed_[j] == -1) { + // v_j = x_j - l_j, v_j >= 0 + // sum_{k != j} d_k x_k + d_j v_j >= beta + // sum_{k != j} d_k x_k + d_j (x_j - l_j) >= beta + // sum_{k != j} d_k x_k + d_j x_j >= beta + d_j l_j + const f_t lj = lb_star_[j]; + inequality.rhs += dj * lj; + } else if (bound_changed_[j] == 1) { + // w_j = u_j - x_j, w_j >= 0 + // sum_{k != j} d_k x_k + d_j w_j >= beta + // sum_{k != j} d_k x_k + d_j (u_j - x_j) >= beta + // sum_{k != j} d_k x_k - d_j x_j >= beta - d_j u_j + const f_t uj = ub_star_[j]; + inequality.rhs -= dj * uj; + inequality.vector.x[k] *= -1.0; + } } - -#ifdef CHECK_WORKSPACE - for (i_t j = 0; j < x_workspace_.size(); j++) { - if (x_workspace_[j] != 0.0) { - printf("After generate_cut: Dirty x_workspace_[%d] = %e\n", j, x_workspace_[j]); - assert(x_workspace_[j] == 0.0); + // Then undo the VUB/VLB substitions and bring continuous variables back to their original form + for (i_t k = 0; k < nz; k++) { + const i_t j = inequality.index(k); + const f_t dj = inequality.coeff(k); + if (var_type[j] != variable_type_t::CONTINUOUS) { + scratch_pad_.add_to_pad(j, dj); + continue; } - if (x_mark_[j] != 0) { - printf("After generate_cut: Dirty x_mark_[%d] = %d\n", j, x_mark_[j]); - assert(x_mark_[j] == 0); + if (bound_changed_[j] == -1) { + if (lb_variable_[j] == -1) { + // v_j = x_j - l_j, v_j >= 0 + // sum_{k != j} d_k x_k + d_j v_j >= beta + // sum_{k != j} d_k x_k + d_j (x_j - l_j) >= beta + // sum_{k != j} d_k x_k + d_j x_j >= beta + d_j l_j + const f_t lj = lb_star_[j]; + inequality.rhs += dj * lj; + scratch_pad_.add_to_pad(j, dj); + } else { + // v_j = x_j - lb*_j, v_j >= 0 + // lb*_j = gamma * x_i + alpha + // v_j = x_j - gamma * x_i - alpha + // sum_{k != j} d_k x_k + d_j v_j >= beta + // sum_{k != j} d_k x_k + d_j (x_j - gamma * x_i - alpha) >= beta + // sum_{k != j} d_k x_k + d_j x_j - d_j * gamma * x_i >= beta + d_j alpha + const i_t p = lb_variable_[j]; + const f_t alpha = variable_bounds.lower_biases[p]; + const f_t gamma = variable_bounds.lower_weights[p]; + const i_t i = variable_bounds.lower_variables[p]; + inequality.rhs += dj * alpha; + scratch_pad_.add_to_pad(j, dj); + scratch_pad_.add_to_pad(i, -dj * gamma); + } + } else if (bound_changed_[j] == 1) { + if (ub_variable_[j] == -1) { + // w_j = u_j - x_j, w_j >= 0 + // sum_{k != j} d_k x_k + d_j w_j >= beta + // sum_{k != j} d_k x_k + d_j (u_j - x_j) >= beta + // sum_{k != j} d_k x_k - d_j x_j >= beta - d_j u_j + const f_t uj = ub_star_[j]; + inequality.rhs -= dj * uj; + scratch_pad_.add_to_pad(j, -dj); + } else { + // w_j = ub*_j - x_j, w_j >= 0 + // ub*_j = gamma * x_i + alpha + // w_j = gamma * x_i + alpha - x_j + // sum_{k != j} d_k x_k + d_j w_j >= beta + // sum_{k != j} d_k x_k + d_j (gamma * x_i + alpha - x_j) >= beta + // sum_{k != j} d_k x_k + d_j gamma * x_i - d_j x_j >= beta - d_j alpha + const i_t p = ub_variable_[j]; + const f_t alpha = variable_bounds.upper_biases[p]; + const f_t gamma = variable_bounds.upper_weights[p]; + const i_t i = variable_bounds.upper_variables[p]; + inequality.rhs -= dj * alpha; + scratch_pad_.add_to_pad(j, -dj); + scratch_pad_.add_to_pad(i, dj * gamma); + } + } else { + scratch_pad_.add_to_pad(j, dj); } } -#endif - // The new cut is: g'*x >= R - // But we want to have it in the form h'*x <= b - cut.sort(); + scratch_pad_.get_pad(inequality.vector.i, inequality.vector.x); + scratch_pad_.clear_pad(); +} - cut_rhs = R; +template +bool complemented_mixed_integer_rounding_cut_t:: + generate_cut_nonnegative_maintain_indicies(const inequality_t& inequality, + const std::vector& var_types, + inequality_t& cut) +{ + auto f = [](f_t q_1, f_t q_2) -> f_t { + f_t q_1_hat = q_1 - std::floor(q_1); + f_t q_2_hat = q_2 - std::floor(q_2); + return std::min(q_1_hat, q_2_hat) + q_2_hat * std::floor(q_1); + }; + + auto h = [](f_t q) -> f_t { return std::max(q, 0.0); }; + + cut.vector = inequality.vector; + const f_t beta = inequality.rhs; + const f_t f_beta = fractional_part(beta); + cut.rhs = f_beta * std::ceil(beta); + if (f_beta < 0.05 || f_beta > 0.95) { return false; } -#ifdef CHECK_REPEATED_INDICES - // Check for repeated indicies - std::vector check(num_vars_, 0); - for (i_t p = 0; p < cut.i.size(); p++) { - if (check[cut.i[p]] != 0) { - printf("repeated index in generated cut\n"); - assert(check[cut.i[p]] == 0); + for (i_t k = 0; k < inequality.size(); k++) { + const i_t j = inequality.index(k); + f_t aj = inequality.coeff(k); + if (var_types[j] == variable_type_t::INTEGER) { + cut.vector.x[k] = f(aj, beta); + } else { + cut.vector.x[k] = h(aj); + } + if (cut.vector.x[k] != cut.vector.x[k]) { + printf("cut.x[%d] %e != cut.x[%d] %e. aj %e beta %e var type %d\n", + k, + cut.vector.x[k], + k, + cut.vector.x[k], + aj, + beta, + static_cast(var_types[j])); + exit(1); } - check[cut.i[p]] = 1; } -#endif - if (cut.i.size() == 0) { return -1; } + return true; +} - return 0; +template +f_t complemented_mixed_integer_rounding_cut_t::compute_violation( + const inequality_t& cut, const std::vector& xstar) +{ + f_t dot = cut.vector.dot(xstar); + f_t cut_violation = cut.rhs - dot; + return cut_violation; } template -void mixed_integer_rounding_cut_t::substitute_slacks(const lp_problem_t& lp, - csr_matrix_t& Arow, - sparse_vector_t& cut, - f_t& cut_rhs) +void complemented_mixed_integer_rounding_cut_t::substitute_slacks( + const lp_problem_t& lp, csr_matrix_t& Arow, inequality_t& cut) { // Remove slacks from the cut // So that the cut is only over the original variables bool found_slack = false; i_t cut_nz = 0; std::vector cut_indices; - cut_indices.reserve(cut.i.size()); + cut_indices.reserve(cut.size()); -#ifdef CHECK_WORKSPACE - for (i_t j = 0; j < x_workspace_.size(); j++) { - if (x_workspace_[j] != 0.0) { - printf("Begin Dirty x_workspace_[%d] = %e\n", j, x_workspace_[j]); - assert(x_workspace_[j] == 0.0); - } - if (x_mark_[j] != 0) { - printf("Begin Dirty x_mark_[%d] = %d\n", j, x_mark_[j]); - assert(x_mark_[j] == 0); - } - } -#endif - - for (i_t k = 0; k < cut.i.size(); k++) { - const i_t j = cut.i[k]; - const f_t cj = cut.x[k]; + for (i_t k = 0; k < cut.size(); k++) { + const i_t j = cut.index(k); + const f_t cj = cut.coeff(k); if (is_slack_[j]) { found_slack = true; const i_t slack_start = lp.A.col_start[j]; @@ -1766,198 +2343,97 @@ void mixed_integer_rounding_cut_t::substitute_slacks(const lp_problem_ // sum_{k != j} C(k) * x_k + sum_{h != j} -C(j)/alpha * A(i, h) * x_h >= cut_rhs - C(j)/alpha // * rhs_i const i_t i = slack_rows_[j]; - cut_rhs -= cj * lp.rhs[i] / alpha; + cut.rhs -= cj * lp.rhs[i] / alpha; const i_t row_start = Arow.row_start[i]; const i_t row_end = Arow.row_start[i + 1]; for (i_t q = row_start; q < row_end; q++) { const i_t h = Arow.j[q]; if (h != j) { const f_t aih = Arow.x[q]; - x_workspace_[h] -= cj * aih / alpha; - if (!x_mark_[h]) { - x_mark_[h] = 1; - cut_indices.push_back(h); - cut_nz++; - } + scratch_pad_.add_to_pad(h, -cj * aih / alpha); } else { const f_t aij = Arow.x[q]; if (std::abs(aij) != 1.0) { - settings_.log.printf( - "Slack row %d has non-unit coefficient %e for variable %d\n", i, aij, j); + printf("Slack row %d has non-unit coefficient %e for variable %d\n", i, aij, j); assert(std::abs(aij) == 1.0); } } } } else { - x_workspace_[j] += cj; - if (!x_mark_[j]) { - x_mark_[j] = 1; - cut_indices.push_back(j); - cut_nz++; - } + scratch_pad_.add_to_pad(j, cj); } } if (found_slack) { - cut.i.reserve(cut_nz); - cut.x.reserve(cut_nz); - cut.i.clear(); - cut.x.clear(); - - for (i_t k = 0; k < cut_nz; k++) { - const i_t j = cut_indices[k]; - - // Check for small coefficients - const f_t aj = x_workspace_[j]; - if (std::abs(aj) < 1e-6) { - if (aj >= 0.0 && lp.upper[j] < inf) { - // Move this to the right-hand side - cut_rhs -= aj * lp.upper[j]; - continue; - } else if (aj <= 0.0 && lp.lower[j] > -inf) { - cut_rhs += aj * lp.lower[j]; - continue; - } else { - } - } - - cut.i.push_back(j); - cut.x.push_back(x_workspace_[j]); - } + scratch_pad_.get_pad(cut.vector.i, cut.vector.x); // Sort the cut cut.sort(); } // Clear the workspace - for (i_t jj : cut_indices) { - x_workspace_[jj] = 0.0; - x_mark_[jj] = 0; - } - -#ifdef CHECK_WORKSPACE - for (i_t j = 0; j < x_workspace_.size(); j++) { - if (x_workspace_[j] != 0.0) { - printf("End Dirty x_workspace_[%d] = %e\n", j, x_workspace_[j]); - assert(x_workspace_[j] == 0.0); - } - if (x_mark_[j] != 0) { - printf("End Dirty x_mark_[%d] = %d\n", j, x_mark_[j]); - assert(x_mark_[j] == 0); - } - } -#endif + scratch_pad_.clear_pad(); } template -f_t mixed_integer_rounding_cut_t::compute_violation(const sparse_vector_t& cut, - f_t cut_rhs, - const std::vector& xstar) -{ - f_t dot = cut.dot(xstar); - f_t cut_violation = cut_rhs - dot; - return cut_violation; -} - -template -void mixed_integer_rounding_cut_t::combine_rows( +f_t complemented_mixed_integer_rounding_cut_t::combine_rows( const lp_problem_t& lp, csr_matrix_t& Arow, i_t xj, - const sparse_vector_t& pivot_row, - f_t pivot_row_rhs, - sparse_vector_t& inequality, - f_t& inequality_rhs) + const inequality_t& pivot_row, + inequality_t& inequality) { -#ifdef CHECK_WORKSPACE - for (i_t k = 0; k < x_workspace_.size(); k++) { - if (x_workspace_[k] != 0.0) { - printf("Dirty x_workspace_[%d] = %e\n", k, x_workspace_[k]); - assert(x_workspace_[k] == 0.0); - } - if (x_mark_[k] != 0) { - printf("Dirty x_mark_[%d] = %d\n", k, x_mark_[k]); - assert(x_mark_[k] == 0); - } - } -#endif - - indices_.clear(); - indices_.reserve(pivot_row.i.size() + inequality.i.size()); - // Find the coefficient associated with variable xj in the pivot row f_t a_l_j = 0.0; - for (i_t k = 0; k < pivot_row.i.size(); k++) { - const i_t j = pivot_row.i[k]; + for (i_t k = 0; k < pivot_row.size(); k++) { + const i_t j = pivot_row.index(k); if (j == xj) { - a_l_j = pivot_row.x[k]; + a_l_j = pivot_row.coeff(k); break; } } - if (a_l_j == 0) { return; } + if (a_l_j == 0) { + printf("Pivot row has no coefficient for variable %d\n", xj); + return 0.0; + } f_t a_i_j = 0.0; - i_t nz = 0; // Store the inequality in the workspace // and save the coefficient associated with variable xj - for (i_t k = 0; k < inequality.i.size(); k++) { - const i_t j = inequality.i[k]; + for (i_t k = 0; k < inequality.size(); k++) { + const i_t j = inequality.index(k); if (j != xj) { - x_workspace_[j] = inequality.x[k]; - x_mark_[j] = 1; - indices_.push_back(j); - nz++; + scratch_pad_.add_to_pad(j, inequality.coeff(k)); } else { - a_i_j = inequality.x[k]; + a_i_j = inequality.coeff(k); } } + if (a_i_j == 0.0) { + printf("Inequality has zero coefficient for variable %d\n", xj); + scratch_pad_.clear_pad(); + return 0.0; + } f_t pivot_value = a_i_j / a_l_j; // Adjust the rhs of the inequality - inequality_rhs -= pivot_value * pivot_row_rhs; + inequality.rhs -= pivot_value * pivot_row.rhs; // Adjust the coefficients of the inequality // based on the nonzeros in the pivot row - for (i_t k = 0; k < pivot_row.i.size(); k++) { - const i_t j = pivot_row.i[k]; - if (j != xj) { - x_workspace_[j] -= pivot_value * pivot_row.x[k]; - if (!x_mark_[j]) { - x_mark_[j] = 1; - indices_.push_back(j); - nz++; - } - } + for (i_t k = 0; k < pivot_row.size(); k++) { + const i_t j = pivot_row.index(k); + if (j != xj) { scratch_pad_.add_to_pad(j, -pivot_value * pivot_row.coeff(k)); } } // Store the new inequality - inequality.i.resize(nz); - inequality.x.resize(nz); - for (i_t k = 0; k < nz; k++) { - inequality.i[k] = indices_[k]; - inequality.x[k] = x_workspace_[indices_[k]]; - } - -#ifdef CHECK_REPEATED_INDICES - // Check for repeated indices - std::vector check(num_vars_, 0); - for (i_t k = 0; k < inequality.i.size(); k++) { - if (check[inequality.i[k]] == 1) { - printf("repeated index\n"); - assert(check[inequality.i[k]] == 0); - } - check[inequality.i[k]] = 1; - } -#endif + scratch_pad_.get_pad(inequality.vector.i, inequality.vector.x); // Clear the workspace - for (i_t j : indices_) { - x_workspace_[j] = 0.0; - x_mark_[j] = 0; - } - indices_.clear(); + scratch_pad_.clear_pad(); + + return -pivot_value; } template @@ -1993,15 +2469,14 @@ i_t strong_cg_cut_t::remove_continuous_variables_integers_nonnegative( const lp_problem_t& lp, const simplex_solver_settings_t& settings, const std::vector& var_types, - sparse_vector_t& inequality, - f_t& inequality_rhs) + inequality_t& inequality) { const bool verbose = false; // Count the number of continuous variables in the inequality i_t num_continuous = 0; - const i_t nz = inequality.i.size(); + const i_t nz = inequality.size(); for (i_t k = 0; k < nz; k++) { - const i_t j = inequality.i[k]; + const i_t j = inequality.index(k); if (var_types[j] == variable_type_t::CONTINUOUS) { num_continuous++; } } @@ -2009,10 +2484,10 @@ i_t strong_cg_cut_t::remove_continuous_variables_integers_nonnegative( // We assume the inequality is of the form sum_j a_j x_j <= rhs for (i_t k = 0; k < nz; k++) { - const i_t j = inequality.i[k]; + const i_t j = inequality.index(k); const f_t l_j = lp.lower[j]; const f_t u_j = lp.upper[j]; - const f_t a_j = inequality.x[k]; + const f_t a_j = inequality.coeff(k); if (var_types[j] == variable_type_t::CONTINUOUS) { if (a_j == 0.0) { continue; } @@ -2022,13 +2497,13 @@ i_t strong_cg_cut_t::remove_continuous_variables_integers_nonnegative( // sum_{k != j} a_k x_k + a_j x_j <= rhs // sum_{k != j} a_k x_k + a_j (v_j + l_j) <= rhs // sum_{k != j} a_k x_k + a_j v_j <= rhs - a_j l_j - inequality_rhs -= a_j * l_j; + inequality.rhs -= a_j * l_j; transformed_variables_[j] = -1; // We now have a_j * v_j with a_j, v_j >= 0 // So we have sum_{k != j} a_k x_k <= sum_{k != j} a_k x_k + a_j v_j <= rhs - a_j l_j // So we can now drop the continuous variable v_j - inequality.x[k] = 0.0; + inequality.vector.x[k] = 0.0; } else if (a_j < 0.0 && u_j < inf) { // w_j = u_j - x_j >= 0 @@ -2036,13 +2511,13 @@ i_t strong_cg_cut_t::remove_continuous_variables_integers_nonnegative( // sum_{k != j} a_k x_k + a_j x_j <= rhs // sum_{k != j} a_k x_k + a_j (u_j - w_j) <= rhs // sum_{k != j} a_k x_k - a_j w_j <= rhs - a_j u_j - inequality_rhs -= a_j * u_j; + inequality.rhs -= a_j * u_j; transformed_variables_[j] = 1; // We now have a_j * w_j with a_j, w_j >= 0 // So we have sum_{k != j} a_k x_k <= sum_{k != j} a_k x_k + a_j w_j <= rhs - a_j u_j // So we can now drop the continuous variable w_j - inequality.x[k] = 0.0; + inequality.vector.x[k] = 0.0; } else { // We can't keep the coefficient of the continuous variable positive // This means we can't eliminate the continuous variable @@ -2058,7 +2533,7 @@ i_t strong_cg_cut_t::remove_continuous_variables_integers_nonnegative( // sum_{k != j} a_k x_k + a_j x_j <= rhs // sum_{k != j} a_k x_k + a_j (v_j + l_j) <= rhs // sum_{k != j} a_k x_k + a_j v_j <= rhs - a_j l_j - inequality_rhs -= a_j * l_j; + inequality.rhs -= a_j * l_j; } else if (transformed_variables_[j] == 1) { // We are closer to the finite upper bound // w_j = u_j - x_j >= 0 @@ -2066,44 +2541,43 @@ i_t strong_cg_cut_t::remove_continuous_variables_integers_nonnegative( // sum_{k != j} a_k x_k + a_j x_j <= rhs // sum_{k != j} a_k x_k + a_j (u_j - w_j) <= rhs // sum_{k != j} a_k x_k - a_j w_j <= rhs - a_j u_j - inequality_rhs -= a_j * u_j; - inequality.x[k] *= -1.0; + inequality.rhs -= a_j * u_j; + inequality.vector.x[k] *= -1.0; } } } // Squeeze out the zero coefficents - sparse_vector_t new_inequality(inequality.n, 0); - inequality.squeeze(new_inequality); - inequality = new_inequality; + sparse_vector_t new_inequality_vector(inequality.vector.n, 0); + inequality.vector.squeeze(new_inequality_vector); + inequality.vector = new_inequality_vector; return 0; } template void strong_cg_cut_t::to_original_integer_variables(const lp_problem_t& lp, - sparse_vector_t& cut, - f_t& cut_rhs) + inequality_t& cut) { // We expect a cut of the form sum_j a_j y_j <= rhs // where y_j >= 0 is a transformed variable // We need to convert it back into a cut on the original variables - for (i_t k = 0; k < cut.i.size(); k++) { - const i_t j = cut.i[k]; - const f_t a_j = cut.x[k]; + for (i_t k = 0; k < cut.size(); k++) { + const i_t j = cut.index(k); + const f_t a_j = cut.coeff(k); if (transformed_variables_[j] == -1) { // sum_{k != j} a_k x_k + a_j v_j <= rhs // v_j = x_j - l_j >= 0, // sum_{k != j} a_k x_k + a_j (x_j - l_j) <= rhs // sum_{k != j} a_k x_k + a_j x_j <= rhs + a_j l_j - cut_rhs += a_j * lp.lower[j]; + cut.rhs += a_j * lp.lower[j]; } else if (transformed_variables_[j] == 1) { // sum_{k != j} a_k x_k + a_j w_j <= rhs // w_j = u_j - x_j >= 0 // sum_{k != j} a_k x_k + a_j (u_j - x_j) <= rhs // sum_{k != j} a_k x_k - a_j x_j <= rhs - a_j u_j - cut_rhs -= a_j * lp.upper[j]; - cut.x[k] *= -1.0; + cut.rhs -= a_j * lp.upper[j]; + cut.vector.x[k] *= -1.0; } } } @@ -2112,44 +2586,37 @@ template i_t strong_cg_cut_t::generate_strong_cg_cut_integer_only( const simplex_solver_settings_t& settings, const std::vector& var_types, - const sparse_vector_t& inequality, - f_t inequality_rhs, - sparse_vector_t& cut, - f_t& cut_rhs) + const inequality_t& inequality, + inequality_t& cut) { // We expect an inequality of the form sum_j a_j x_j <= rhs // where all the variables x_j are integer and nonnegative // We then apply the CG cut: // sum_j floor(a_j) x_j <= floor(rhs) - cut.i.reserve(inequality.i.size()); - cut.x.reserve(inequality.i.size()); - cut.i.clear(); - cut.x.clear(); + cut.reserve(inequality.size()); + cut.clear(); - f_t a_0 = inequality_rhs; + f_t a_0 = inequality.rhs; f_t f_a_0 = fractional_part(a_0); if (f_a_0 == 0.0) { // f(a_0) == 0.0 so we do a weak CG cut - cut.i.reserve(inequality.i.size()); - cut.x.reserve(inequality.i.size()); - cut.i.clear(); - cut.x.clear(); - for (i_t k = 0; k < inequality.i.size(); k++) { - const i_t j = inequality.i[k]; - const f_t a_j = inequality.x[k]; + cut.reserve(inequality.size()); + cut.clear(); + for (i_t k = 0; k < inequality.size(); k++) { + const i_t j = inequality.index(k); + const f_t a_j = inequality.coeff(k); if (var_types[j] == variable_type_t::INTEGER) { - cut.i.push_back(j); - cut.x.push_back(std::floor(a_j)); + cut.push_back(j, std::floor(a_j)); } else { return -1; } } - cut_rhs = std::floor(inequality_rhs); + cut.rhs = std::floor(inequality.rhs); } else { return generate_strong_cg_cut_helper( - inequality.i, inequality.x, inequality_rhs, var_types, cut, cut_rhs); + inequality.vector.i, inequality.vector.x, inequality.rhs, var_types, cut); } return 0; } @@ -2160,8 +2627,7 @@ i_t strong_cg_cut_t::generate_strong_cg_cut_helper( const std::vector& coefficients, f_t rhs, const std::vector& var_types, - sparse_vector_t& cut, - f_t& cut_rhs) + inequality_t& cut) { const bool verbose = false; const i_t nz = indicies.size(); @@ -2182,10 +2648,8 @@ i_t strong_cg_cut_t::generate_strong_cg_cut_helper( f_t upper = 1.0 / static_cast(k); if (verbose) { printf("f_a_0 %e lower %e upper %e alpha %e\n", f_a_0, lower, upper, alpha); } if (f_a_0 >= lower && f_a_0 < upper) { - cut.i.reserve(nz); - cut.x.reserve(nz); - cut.i.clear(); - cut.x.clear(); + cut.reserve(nz); + cut.clear(); for (i_t q = 0; q < nz; q++) { const i_t j = indicies[q]; const f_t a_j = coefficients[q]; @@ -2193,8 +2657,7 @@ i_t strong_cg_cut_t::generate_strong_cg_cut_helper( const f_t f_a_j = fractional_part(a_j); const f_t tol = 1e-4; if (f_a_j <= f_a_0 + tol) { - cut.i.push_back(j); - cut.x.push_back((k + 1.0) * std::floor(a_j)); + cut.push_back(j, (k + 1.0) * std::floor(a_j)); if (verbose) { printf("j %d a_j %e f_a_j %e k %d\n", j, a_j, f_a_j, k); } } else { // Find p such that p <= k * f(a_j) < p + 1 @@ -2203,11 +2666,9 @@ i_t strong_cg_cut_t::generate_strong_cg_cut_helper( const f_t rhs_j = f_a_0 + static_cast(p) / static_cast(k) * alpha; const i_t coeff = (k + 1) * static_cast(std::floor(a_j)) + p; if (f_a_j > rhs_j + tol) { - cut.i.push_back(j); - cut.x.push_back(static_cast(coeff + 1)); + cut.push_back(j, static_cast(coeff + 1)); } else { - cut.i.push_back(j); - cut.x.push_back(static_cast(coeff)); + cut.push_back(j, static_cast(coeff)); } } } else { @@ -2218,11 +2679,11 @@ i_t strong_cg_cut_t::generate_strong_cg_cut_helper( if (verbose) { printf("Error: k %d lower %e f(a_0) %e upper %e\n", k, lower, f_a_0, upper); } return -1; } - cut_rhs = (k + 1.0) * std::floor(rhs); + cut.rhs = (k + 1.0) * std::floor(rhs); if (verbose) { - printf("Generated strong CG cut: k %d f_a_0 %e cut_rhs %e\n", k, f_a_0, cut_rhs); - for (i_t q = 0; q < cut.i.size(); q++) { - if (cut.x[q] != 0.0) { printf("%.16e x%d ", cut.x[q], cut.i[q]); } + printf("Generated strong CG cut: k %d f_a_0 %e cut_rhs %e\n", k, f_a_0, cut.rhs); + for (i_t q = 0; q < cut.size(); q++) { + if (cut.vector.x[q] != 0.0) { printf("%.16e x%d ", cut.vector.x[q], cut.vector.i[q]); } } printf("\n"); printf("Original inequality rhs %e nz %ld\n", rhs, coefficients.size()); @@ -2239,11 +2700,9 @@ i_t strong_cg_cut_t::generate_strong_cg_cut( const lp_problem_t& lp, const simplex_solver_settings_t& settings, const std::vector& var_types, - const sparse_vector_t& inequality, - const f_t inequality_rhs, + const inequality_t& inequality, const std::vector& xstar, - sparse_vector_t& cut, - f_t& cut_rhs) + inequality_t& cut) { #ifdef PRINT_INEQUALITY_INFO for (i_t k = 0; k < inequality.i.size(); k++) { @@ -2258,36 +2717,33 @@ i_t strong_cg_cut_t::generate_strong_cg_cut( // and transform integer variables to be nonnegative // Copy the inequality since remove continuous variables will modify it - sparse_vector_t cg_inequality = inequality; - f_t cg_inequality_rhs = inequality_rhs; - i_t status = remove_continuous_variables_integers_nonnegative( - lp, settings, var_types, cg_inequality, cg_inequality_rhs); + inequality_t cg_inequality = inequality; + i_t status = + remove_continuous_variables_integers_nonnegative(lp, settings, var_types, cg_inequality); if (status != 0) { // Try negating the equality and see if that helps cg_inequality = inequality; cg_inequality.negate(); - cg_inequality_rhs = -inequality_rhs; - status = remove_continuous_variables_integers_nonnegative( - lp, settings, var_types, cg_inequality, cg_inequality_rhs); + status = + remove_continuous_variables_integers_nonnegative(lp, settings, var_types, cg_inequality); } if (status == 0) { // We have an inequality with no continuous variables // Generate a CG cut - status = generate_strong_cg_cut_integer_only( - settings, var_types, cg_inequality, cg_inequality_rhs, cut, cut_rhs); + status = generate_strong_cg_cut_integer_only(settings, var_types, cg_inequality, cut); if (status != 0) { return -1; } // Convert the CG cut back to the original variables - to_original_integer_variables(lp, cut, cut_rhs); + to_original_integer_variables(lp, cut); // Check for violation - f_t dot = cut.dot(xstar); + f_t dot = cut.vector.dot(xstar); // If the cut is violated we will have: sum_j a_j xstar_j > rhs - f_t violation = dot - cut_rhs; + f_t violation = dot - cut.rhs; const f_t min_violation_threshold = 1e-6; if (violation > min_violation_threshold) { // Note that no slacks are currently present. Since slacks are currently treated as @@ -2296,7 +2752,6 @@ i_t strong_cg_cut_t::generate_strong_cg_cut( // The CG cut is in the form: sum_j a_j x_j <= rhs // The cut pool wants the cut in the form: sum_j a_j x_j >= rhs cut.negate(); - cut_rhs *= -1.0; return 0; } } @@ -2784,7 +3239,8 @@ template class cut_pool_t; template class cut_generation_t; template class knapsack_generation_t; template class tableau_equality_t; -template class mixed_integer_rounding_cut_t; +template class complemented_mixed_integer_rounding_cut_t; +template class variable_bounds_t; template int add_cuts(const simplex_solver_settings_t& settings, const csr_matrix_t& cuts, diff --git a/cpp/src/cuts/cuts.hpp b/cpp/src/cuts/cuts.hpp index 4f55e96e4d..46f74193a6 100644 --- a/cpp/src/cuts/cuts.hpp +++ b/cpp/src/cuts/cuts.hpp @@ -32,6 +32,59 @@ enum cut_type_t : int8_t { MAX_CUT_TYPE = 4 }; +template +struct inequality_t { + inequality_t() : vector(), rhs(0.0) {} + inequality_t(i_t num_cols) : vector(num_cols, 0), rhs(0.0) {} + inequality_t(csr_matrix_t& A, i_t row, f_t rhs_value) : vector(A, row), rhs(rhs_value) + { + } + sparse_vector_t vector; + f_t rhs; + + void push_back(i_t j, f_t x) + { + vector.i.push_back(j); + vector.x.push_back(x); + } + void clear() + { + vector.i.clear(); + vector.x.clear(); + } + void reserve(size_t n) + { + vector.i.reserve(n); + vector.x.reserve(n); + } + size_t size() const { return vector.i.size(); } + i_t index(i_t k) const { return vector.i[k]; } + f_t coeff(i_t k) const { return vector.x[k]; } + void negate() + { + vector.negate(); + rhs *= -1.0; + } + void sort() { vector.sort(); } + void squeeze(inequality_t& out) const + { + vector.squeeze(out.vector); + out.rhs = rhs; + } + void scale(f_t factor) + { + vector.scale(factor); + rhs *= factor; + } + void print() const + { + for (i_t k = 0; k < size(); k++) { + printf("%g x%d ", coeff(k), index(k)); + } + printf("\nrhs %g\n", rhs); + } +}; + template struct cut_info_t { bool has_cuts() const @@ -72,10 +125,8 @@ void print_cut_types(const std::string& prefix, cut_info.record_cut_types(cut_types); settings.log.printf("%s: ", prefix.c_str()); for (i_t i = 0; i < MAX_CUT_TYPE; i++) { - settings.log.printf("%s cuts: %d ", cut_info.cut_type_names[i], cut_info.num_cuts[i]); - if (i < MAX_CUT_TYPE - 1) { settings.log.printf(", "); } + settings.log.printf("%s cuts: %d\n", cut_info.cut_type_names[i], cut_info.num_cuts[i]); } - settings.log.printf("\n"); } template @@ -136,7 +187,7 @@ class cut_pool_t { // Add a cut in the form: cut'*x >= rhs. // We expect that the cut is violated by the current relaxation xstar // cut'*xstart < rhs - void add_cut(cut_type_t cut_type, const sparse_vector_t& cut, f_t rhs); + void add_cut(cut_type_t cut_type, const inequality_t& cut); void score_cuts(std::vector& x_relax); @@ -172,6 +223,7 @@ class cut_pool_t { std::vector cut_orthogonality_; std::vector cut_scores_; std::vector best_cuts_; + const f_t min_cut_distance_{1e-4}; }; template @@ -190,8 +242,7 @@ class knapsack_generation_t { const std::vector& var_types, const std::vector& xstar, i_t knapsack_row, - sparse_vector_t& cut, - f_t& cut_rhs); + inequality_t& cut); i_t num_knapsack_constraints() const { return knapsack_constraints_.size(); } const std::vector& get_knapsack_constraints() const { return knapsack_constraints_; } @@ -214,10 +265,13 @@ class knapsack_generation_t { const simplex_solver_settings_t& settings_; }; -// Forward declaration +// Forward declarations template class mixed_integer_rounding_cut_t; +template +class variable_bounds_t; + template class cut_generation_t { public: @@ -238,8 +292,10 @@ class cut_generation_t { const std::vector& var_types, basis_update_mpf_t& basis_update, const std::vector& xstar, + const std::vector& ystar, const std::vector& basic_list, - const std::vector& nonbasic_list); + const std::vector& nonbasic_list, + variable_bounds_t& variable_bounds); private: // Generate all mixed integer gomory cuts @@ -259,7 +315,9 @@ class cut_generation_t { csr_matrix_t& Arow, const std::vector& new_slacks, const std::vector& var_types, - const std::vector& xstar); + const std::vector& xstar, + const std::vector& ystar, + variable_bounds_t& variable_bounds); // Generate all knapsack cuts void generate_knapsack_cuts(const lp_problem_t& lp, @@ -273,6 +331,77 @@ class cut_generation_t { knapsack_generation_t knapsack_generation_; }; +template +class scratch_pad_t { + public: + scratch_pad_t(i_t num_vars) : workspace_(num_vars, 0.0), mark_(num_vars, 0) + { + indices_.reserve(num_vars); + } + + // O(1) to add a value to the pad + void add_to_pad(i_t j, f_t value) + { + workspace_[j] += value; + if (!mark_[j]) { + mark_[j] = 1; + indices_.push_back(j); + } + } + + // O(nz) to clear the pad + void clear_pad() + { + for (i_t j : indices_) { + workspace_[j] = 0.0; + mark_[j] = 0; + } + indices_.clear(); + } + + // O(nz) to get the pad + void get_pad(std::vector& indices, std::vector& values) + { + indices.reserve(indices_.size()); + values.reserve(indices_.size()); + indices.clear(); + values.clear(); + const i_t nz = indices_.size(); + for (i_t k = 0; k < nz; k++) { + const i_t j = indices_[k]; + const f_t val = workspace_[j]; + if (val != 0.0) { + indices.push_back(j); + values.push_back(val); + } + } + } + + private: + std::vector workspace_; + std::vector mark_; + std::vector indices_; +}; + +template +class mixed_integer_gomory_cut_t { + public: + mixed_integer_gomory_cut_t() {} + + bool rational_coefficients(const std::vector& var_types, + const inequality_t& input_inequality, + inequality_t& rational_inequality); + + private: + bool rational_approximation(f_t x, + int64_t max_denominator, + int64_t& numerator, + int64_t& denominator); + + int64_t gcd(const std::vector& integers); + int64_t lcm(const std::vector& integers); +}; + template class tableau_equality_t { public: @@ -301,8 +430,7 @@ class tableau_equality_t { const std::vector& basic_list, const std::vector& nonbasic_list, i_t i, - sparse_vector_t& inequality, - f_t& inequality_rhs); + inequality_t& inequality); private: std::vector b_bar_; @@ -313,93 +441,227 @@ class tableau_equality_t { }; template -class mixed_integer_rounding_cut_t { +class variable_bounds_t { public: - mixed_integer_rounding_cut_t(const lp_problem_t& lp, - const simplex_solver_settings_t& settings, - const std::vector& new_slacks, - const std::vector& xstar); + variable_bounds_t(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& var_types, + const csr_matrix_t& Arow, + const std::vector& new_slacks); + + std::vector upper_offsets; + std::vector upper_variables; + std::vector upper_weights; + std::vector upper_biases; + + std::vector lower_offsets; + std::vector lower_variables; + std::vector lower_weights; + std::vector lower_biases; + + void resize(i_t new_num_cols) + { + const i_t current_upper_nz = upper_offsets.back(); + upper_offsets.resize(new_num_cols + 1, current_upper_nz); + const i_t current_lower_nz = lower_offsets.back(); + lower_offsets.resize(new_num_cols + 1, current_lower_nz); + } + + private: + f_t lower_activity(f_t lower_bound, f_t upper_bound, f_t coefficient) + { + return (coefficient > 0.0 ? lower_bound : upper_bound) * coefficient; + } + + f_t upper_activity(f_t lower_bound, f_t upper_bound, f_t coefficient) + { + return (coefficient > 0.0 ? upper_bound : lower_bound) * coefficient; + } + + // Returns the lower activity adjusted for the number of lower inf variables + // adjusted_lower_activity = { activity - lower_activity_i - lower_activity_j, if num_lower_inf = + // 0 + // { activity - lower_activity_i , if num_lower_inf = + // 1, lower_activity_j = -inf { activity - lower_activity_j , if + // num_lower_inf = 1, lower_activity_i != -inf { activity , if + // num_lower_inf = 2, lower_activity_i = lower_activity_j = -inf { -inf + // , if num_lower_inf > 2 + f_t adjusted_lower_activity(f_t activity, + i_t num_lower_inf, + f_t lower_activity_i, + f_t lower_activity_j) + { + if (num_lower_inf == 0) { + return activity - lower_activity_i - lower_activity_j; + } else if (num_lower_inf == 1 && lower_activity_j == -inf) { + return activity - lower_activity_i; + } else if (num_lower_inf == 1 && lower_activity_i == -inf) { + return activity - lower_activity_j; + } else if (num_lower_inf == 2 && lower_activity_i == -inf && lower_activity_j == -inf) { + return activity; + } else { + return -inf; + } + } + + // Returns the upper activity adjusted for the number of upper inf variables + // adjusted_upper_activity = { activity - upper_activity_i - upper_activity_j, if num_upper_inf = + // 0 + // { activity - upper_activity_i , if num_upper_inf = + // 1, upper_activity_j = inf { activity - upper_activity_j , if + // num_upper_inf = 1, upper_activity_i != inf { activity , if + // num_upper_inf = 2, upper_activity_i = upper_activity_j = inf { inf , + // if num_upper_inf > 2 + f_t adjusted_upper_activity(f_t activity, + i_t num_upper_inf, + f_t upper_activity_i, + f_t upper_activity_j) + { + if (num_upper_inf == 0) { + return activity - upper_activity_i - upper_activity_j; + } else if (num_upper_inf == 1 && upper_activity_j == inf) { + return activity - upper_activity_i; + } else if (num_upper_inf == 1 && upper_activity_i == inf) { + return activity - upper_activity_j; + } else if (num_upper_inf == 2 && upper_activity_i == inf && upper_activity_j == inf) { + return activity; + } else { + return inf; + } + } + + std::vector upper_activities_; + std::vector num_pos_inf_; + std::vector lower_activities_; + std::vector num_neg_inf_; - // Convert an inequality of the form: sum_j a_j x_j >= beta + std::vector slack_map_; +}; + +template +class complemented_mixed_integer_rounding_cut_t { + public: + complemented_mixed_integer_rounding_cut_t(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const std::vector& new_slacks); + + void compute_initial_scores_for_rows(const lp_problem_t& lp, + const simplex_solver_settings_t& settings, + const csr_matrix_t& Arow, + const std::vector& xstar, + const std::vector& ystar, + std::vector& score); + + // Perform bound substitution for the continuous variables using simple bounds + // and variable bounds. And bound substitution for the integer variables + // using simple bounds. + void bound_substitution(const lp_problem_t& lp, + const variable_bounds_t& variable_bounds, + const std::vector& var_types, + const std::vector& xstar, + std::vector& transformed_xstar); + + // Converts an inequality of the form: sum_j a_j x_j >= beta // with l_j <= x_j <= u_j into the form: // sum_{j not in L union U} d_j x_j + sum_{j in L} d_j v_j // + sum_{j in U} d_j w_j >= delta, // where v_j = x_j - l_j for j in L - // and w_j = u_j - x_j for j in Us - void to_nonnegative(const lp_problem_t& lp, - sparse_vector_t& inequality, - f_t& rhs); - - void relaxation_to_nonnegative(const lp_problem_t& lp, - const std::vector& xstar, - std::vector& xstar_nonnegative); + // and w_j = u_j - x_j for j in U + void transform_inequality(const variable_bounds_t& variable_bounds, + const std::vector& var_type, + inequality_t& inequality); - // Convert an inequality of the form: + // Converts an inequality of the form: // sum_{j not in L union U} d_j x_j + sum_{j in L} d_j v_j - // + sum_{j in U} d_j w_j >= delta + // + sum_{j in U} d_j w_j >= delta, // where v_j = x_j - l_j for j in L // and w_j = u_j - x_j for j in U - // back to an inequality on the original variables - // sum_j a_j x_j >= beta - void to_original(const lp_problem_t& lp, - sparse_vector_t& inequality, - f_t& rhs); + // back to the form: sum_j a_j x_j >= beta + // with l_j <= x_j <= u_j + void untransform_inequality(const variable_bounds_t& variable_bounds, + const std::vector& var_type, + inequality_t& inequality); + + bool cut_generation_heuristic(const inequality_t& transformed_inequality, + const std::vector& var_types, + const std::vector& transformed_xstar, + inequality_t& transformed_cut, + f_t& work_estimate); + + bool scale_uncomplement_and_generate_cut(const std::vector& var_types, + const std::vector& transformed_xstar, + const std::vector& complemented_indices, + const inequality_t& complemented_inequality, + f_t delta, + inequality_t& cut_delta, + f_t& work_estimate); + + // This routine takes an inequality and generates the MIR cut + bool generate_cut_nonnegative_maintain_indicies(const inequality_t& inequality, + const std::vector& var_types, + inequality_t& cut); + + f_t compute_violation(const inequality_t& cut, const std::vector& xstar); + + f_t new_upper(i_t j) const { return transformed_upper_[j]; } // Given a cut of the form sum_j d_j x_j >= beta // with l_j <= x_j <= u_j, try to remove coefficients d_j // with | d_j | < epsilon void remove_small_coefficients(const std::vector& lower_bounds, const std::vector& upper_bounds, - sparse_vector_t& cut, - f_t& cut_rhs); - - // Given an inequality sum_j a_j x_j >= beta, x_j >= 0, x_j in Z, j in I - // generate an MIR cut of the form sum_j d_j x_j >= delta - i_t generate_cut_nonnegative(const sparse_vector_t& a, - f_t beta, - const std::vector& var_types, - sparse_vector_t& cut, - f_t& cut_rhs); - - f_t compute_violation(const sparse_vector_t& cut, - f_t cut_rhs, - const std::vector& xstar); - - i_t generate_cut(const sparse_vector_t& a, - f_t beta, - const std::vector& upper_bounds, - const std::vector& lower_bounds, - const std::vector& var_types, - sparse_vector_t& cut, - f_t& cut_rhs); + inequality_t& cut); void substitute_slacks(const lp_problem_t& lp, csr_matrix_t& Arow, - sparse_vector_t& cut, - f_t& cut_rhs); + inequality_t& cut); // Combine the pivot row with the inequality to eliminate the variable j // The new inequality is returned in inequality and inequality_rhs - void combine_rows(const lp_problem_t& lp, - csr_matrix_t& Arow, - i_t j, - const sparse_vector_t& pivot_row, - f_t pivot_row_rhs, - sparse_vector_t& inequality, - f_t& inequality_rhs); + // The multiplier for the pivot row is returned + f_t combine_rows(const lp_problem_t& lp, + csr_matrix_t& Arow, + i_t j, + const inequality_t& pivot_row, + inequality_t& inequality); + + const f_t get_lb_star(i_t j) const { return lb_star_[j]; } + const f_t get_ub_star(i_t j) const { return ub_star_[j]; } + + const i_t slack_rows(i_t j) const { return slack_rows_[j]; } + const i_t slack_cols(i_t i) const { return slack_cols_[i]; } + + bool scale_and_generate_mir_cut(const std::vector& var_types, + const std::vector& transformed_xstar, + const inequality_t& inequality, + f_t divisor, + std::vector>& cuts, + std::vector& violations, + std::vector& deltas); + + bool check_violation_and_add_cut(const inequality_t& inequality, + const std::vector& xstar, + f_t divisor, + std::vector>& cuts, + std::vector& violations, + std::vector& deltas); private: - i_t num_vars_; - const simplex_solver_settings_t& settings_; - std::vector x_workspace_; - std::vector x_mark_; - std::vector has_lower_; - std::vector has_upper_; std::vector is_slack_; - std::vector slack_rows_; - std::vector indices_; - std::vector bound_info_; - bool needs_complement_; + std::vector + slack_rows_; // slack_rows_[j] = i, if variable j is slack for row i, -1 is sentinal value + std::vector + slack_cols_; // slack_cols_[i] = j, if variable j is slack for row i -1 is sentinal value + + std::vector lb_variable_; + std::vector lb_star_; + std::vector ub_variable_; + std::vector ub_star_; + + std::vector bound_changed_; + std::vector transformed_upper_; + + scratch_pad_t scratch_pad_; }; template @@ -412,37 +674,29 @@ class strong_cg_cut_t { i_t generate_strong_cg_cut(const lp_problem_t& lp, const simplex_solver_settings_t& settings, const std::vector& var_types, - const sparse_vector_t& inequality, - const f_t inequality_rhs, + const inequality_t& inequality, const std::vector& xstar, - sparse_vector_t& cut, - f_t& cut_rhs); + inequality_t& cut); i_t remove_continuous_variables_integers_nonnegative( const lp_problem_t& lp, const simplex_solver_settings_t& settings, const std::vector& var_types, - sparse_vector_t& inequality, - f_t& inequality_rhs); + inequality_t& inequality); - void to_original_integer_variables(const lp_problem_t& lp, - sparse_vector_t& cut, - f_t& cut_rhs); + void to_original_integer_variables(const lp_problem_t& lp, inequality_t& cut); i_t generate_strong_cg_cut_integer_only(const simplex_solver_settings_t& settings, const std::vector& var_types, - const sparse_vector_t& inequality, - f_t inequality_rhs, - sparse_vector_t& cut, - f_t& cut_rhs); + const inequality_t& inequality, + inequality_t& cut); private: i_t generate_strong_cg_cut_helper(const std::vector& indicies, const std::vector& coefficients, f_t rhs, const std::vector& var_types, - sparse_vector_t& cut, - f_t& cut_rhs); + inequality_t& cut); std::vector transformed_variables_; }; diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index 426d9a7535..0794ed1681 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -412,19 +412,23 @@ void compute_delta_z(const csc_matrix_t& A_transpose, delta_z[leaving_index] = direction; #ifdef CHECK_CHANGE_IN_REDUCED_COST - delta_y_sparse.to_dense(delta_y); + const i_t m = A_transpose.n; + const i_t n = A_transpose.m; + std::vector delta_y_dense(m); + delta_y.to_dense(delta_y_dense); std::vector delta_z_check(n); std::vector delta_z_mark_check(n, 0); std::vector delta_z_indices_check; phase2::compute_reduced_cost_update(lp, basic_list, nonbasic_list, - delta_y, + delta_y_dense, leaving_index, direction, delta_z_mark_check, delta_z_indices_check, - delta_z_check); + delta_z_check, + work_estimate); f_t error_check = 0.0; for (i_t k = 0; k < n; ++k) { const f_t diff = std::abs(delta_z[k] - delta_z_check[k]); @@ -1726,6 +1730,7 @@ i_t compute_delta_x(const lp_problem_t& lp, const std::vector& basic_list, const std::vector& delta_x_flip, const sparse_vector_t& rhs_sparse, + const std::vector& delta_z, const std::vector& x, sparse_vector_t& utilde_sparse, sparse_vector_t& scaled_delta_xB_sparse, @@ -1782,6 +1787,9 @@ i_t compute_delta_x(const lp_problem_t& lp, scaled_delta_xB_sparse.negate(); work_estimate += 2 * scaled_delta_xB_sparse.i.size() + scaled_delta_xB.size(); scale = -scaled_delta_xB[basic_leaving_index]; + } else if (delta_z[entering_index] != 0.0) { + printf("Using delta_z for entering index %d %e\n", entering_index, delta_z[entering_index]); + scale = -delta_z[entering_index]; } else { return -1; } @@ -2029,8 +2037,8 @@ void check_primal_infeasibilities(const lp_problem_t& lp, const simplex_solver_settings_t& settings, const std::vector& basic_list, const std::vector& x, - const ins_vector& squared_infeasibilities, - const ins_vector& infeasibility_indices) + const std::vector& squared_infeasibilities, + const std::vector& infeasibility_indices) { const i_t m = basic_list.size(); for (i_t k = 0; k < m; ++k) { @@ -2054,14 +2062,30 @@ void check_primal_infeasibilities(const lp_problem_t& lp, } } if (!found) { settings.log.printf("Infeasibility index not found %d\n", j); } + } else { + bool found = false; + i_t h; + for (h = 0; h < infeasibility_indices.size(); ++h) { + if (infeasibility_indices[h] == j) { + found = true; + break; + } + } + if (found) { + settings.log.printf("Incorrect infeasible index %d/%d infeas %e sq %e\n", + j, + h, + infeas, + squared_infeasibilities[j]); + } } } } template void check_basic_infeasibilities(const std::vector& basic_list, - const ins_vector& basic_mark, - const ins_vector& infeasibility_indices, + const std::vector& basic_mark, + const std::vector& infeasibility_indices, i_t info) { for (i_t k = 0; k < infeasibility_indices.size(); ++k) { @@ -2104,8 +2128,8 @@ template void check_basis_mark(const simplex_solver_settings_t& settings, const std::vector& basic_list, const std::vector& nonbasic_list, - const ins_vector& basic_mark, - const ins_vector& nonbasic_mark) + const std::vector& basic_mark, + const std::vector& nonbasic_mark) { const i_t m = basic_list.size(); const i_t n = basic_mark.size(); @@ -2925,7 +2949,9 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, #ifdef COMPUTE_DUAL_RESIDUAL std::vector dual_residual; std::vector zeros(n, 0.0); - phase2::compute_dual_residual(lp.A, zeros, delta_y, delta_z, dual_residual); + std::vector delta_y_dense(m); + delta_y_sparse.to_dense(delta_y_dense); + phase2::compute_dual_residual(lp.A, zeros, delta_y_dense, delta_z, dual_residual); // || A'*delta_y + delta_z ||_inf f_t dual_residual_norm = vector_norm_inf(dual_residual); settings.log.printf( @@ -3182,6 +3208,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, timers.vector_time += timers.stop_timer(); #ifdef COMPUTE_DUAL_RESIDUAL + std::vector dual_res1; phase2::compute_dual_residual(lp.A, objective, y, z, dual_res1); f_t dual_res_norm = vector_norm_inf(dual_res1); if (dual_res_norm > settings.dual_tol) { @@ -3241,6 +3268,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, basic_list, delta_x_flip, rhs_sparse, + delta_z, x, utilde_sparse, scaled_delta_xB_sparse, @@ -3406,7 +3434,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, if (should_refactor) { PHASE2_NVTX_RANGE("DualSimplex::refactorization"); num_refactors++; - bool should_recompute_x = false; + bool should_recompute_x = true; // Need for numerically difficult problems like cbs-cta i_t refactor_status = ft.refactor_basis( lp.A, settings, lp.lower, lp.upper, start_time, basic_list, nonbasic_list, vstatus); if (refactor_status == CONCURRENT_HALT_RETURN) { return dual::status_t::CONCURRENT_LIMIT; } diff --git a/cpp/src/dual_simplex/presolve.hpp b/cpp/src/dual_simplex/presolve.hpp index 17e6176e3b..a068ed04ab 100644 --- a/cpp/src/dual_simplex/presolve.hpp +++ b/cpp/src/dual_simplex/presolve.hpp @@ -13,6 +13,13 @@ #include #include +#include +#include +#include +#include +#include +#include + namespace cuopt::linear_programming::dual_simplex { template @@ -42,6 +49,113 @@ struct lp_problem_t { f_t obj_constant; f_t obj_scale; // 1.0 for min, -1.0 for max bool objective_is_integral{false}; + + void write_problem(const std::string& path) const + { + FILE* fid = fopen(path.c_str(), "w"); + if (fid) { + fwrite(&num_rows, sizeof(i_t), 1, fid); + fwrite(&num_cols, sizeof(i_t), 1, fid); + fwrite(&obj_constant, sizeof(f_t), 1, fid); + fwrite(&obj_scale, sizeof(f_t), 1, fid); + i_t is_integral = objective_is_integral ? 1 : 0; + fwrite(&is_integral, sizeof(i_t), 1, fid); + fwrite(objective.data(), sizeof(f_t), num_cols, fid); + fwrite(rhs.data(), sizeof(f_t), num_rows, fid); + fwrite(lower.data(), sizeof(f_t), num_cols, fid); + fwrite(upper.data(), sizeof(f_t), num_cols, fid); + fwrite(A.col_start.data(), sizeof(i_t), A.col_start.size(), fid); + fwrite(A.i.data(), sizeof(i_t), A.i.size(), fid); + fwrite(A.x.data(), sizeof(f_t), A.x.size(), fid); + fclose(fid); + } + } + + void read_problem(const std::string& path) + { + FILE* fid = fopen(path.c_str(), "r"); + if (fid) { + fread(&num_rows, sizeof(i_t), 1, fid); + fread(&num_cols, sizeof(i_t), 1, fid); + fread(&obj_constant, sizeof(f_t), 1, fid); + fread(&obj_scale, sizeof(f_t), 1, fid); + i_t is_integral; + fread(&is_integral, sizeof(i_t), 1, fid); + objective_is_integral = is_integral == 1; + objective.resize(num_cols); + fread(objective.data(), sizeof(f_t), num_cols, fid); + rhs.resize(num_rows); + fread(rhs.data(), sizeof(f_t), num_rows, fid); + lower.resize(num_cols); + fread(lower.data(), sizeof(f_t), num_cols, fid); + upper.resize(num_cols); + fread(upper.data(), sizeof(f_t), num_cols, fid); + A.n = num_cols; + A.m = num_rows; + A.col_start.resize(num_cols + 1); + fread(A.col_start.data(), sizeof(i_t), num_cols + 1, fid); + A.i.resize(A.col_start[num_cols]); + fread(A.i.data(), sizeof(i_t), A.i.size(), fid); + A.x.resize(A.i.size()); + fread(A.x.data(), sizeof(f_t), A.x.size(), fid); + fclose(fid); + } + } + + void write_mps(const std::string& path) const + { + std::ofstream mps_file(path); + if (!mps_file.is_open()) { + printf("Failed to open file %s\n", path.c_str()); + return; + } + mps_file << std::setprecision(std::numeric_limits::max_digits10); + mps_file << "NAME " << "cuopt_lp_problem_t" << "\n"; + mps_file << "ROWS\n"; + mps_file << " N OBJ\n"; + for (i_t i = 0; i < num_rows; i++) { + mps_file << " E R" << i << "\n"; + } + mps_file << "COLUMNS\n"; + for (i_t j = 0; j < num_cols; j++) { + const i_t col_start = A.col_start[j]; + const i_t col_end = A.col_start[j + 1]; + mps_file << " " << "C" << j << " OBJ " << objective[j] << "\n"; + for (i_t k = col_start; k < col_end; k++) { + const i_t i = A.i[k]; + const f_t x = A.x[k]; + std::string col_name = "C" + std::to_string(j); + std::string row_name = "R" + std::to_string(i); + mps_file << " " << col_name << " " << row_name << " " << x << "\n"; + } + } + mps_file << "RHS\n"; + for (i_t i = 0; i < num_rows; i++) { + mps_file << " RHS1 R" << i << " " << rhs[i] << "\n"; + } + + mps_file << "BOUNDS\n"; + for (i_t j = 0; j < num_cols; j++) { + const f_t lb = lower[j]; + const f_t ub = upper[j]; + std::string col_name = "C" + std::to_string(j); + if (lb == -std::numeric_limits::infinity() && + ub == std::numeric_limits::infinity()) { + mps_file << " FR BOUND1 " << col_name << "\n"; + } else { + if (lb == -std::numeric_limits::infinity()) { + mps_file << " MI BOUND1 " << col_name << "\n"; + } else { + mps_file << " LO BOUND1 " << col_name << " " << lb << "\n"; + } + if (ub != std::numeric_limits::infinity()) { + mps_file << " UP BOUND1 " << col_name << " " << ub << "\n"; + } + } + } + mps_file << "ENDATA\n"; + mps_file.close(); + } }; template diff --git a/cpp/src/dual_simplex/sparse_matrix.hpp b/cpp/src/dual_simplex/sparse_matrix.hpp index 56e3ca82c6..a78e0a4044 100644 --- a/cpp/src/dual_simplex/sparse_matrix.hpp +++ b/cpp/src/dual_simplex/sparse_matrix.hpp @@ -48,7 +48,12 @@ class csc_matrix_t { // Adjust to i and x vectors for a new number of nonzeros void reallocate(i_t new_nz); + // Get the number of nonzeros in the matrix i_t nnz() const { return col_start[n]; } + + // Get the number of nonzeros in column j + i_t col_length(i_t j) const { return col_start[j + 1] - col_start[j]; } + // Convert the CSC matrix to a CSR matrix i_t to_compressed_row( cuopt::linear_programming::dual_simplex::csr_matrix_t& Arow) const; @@ -145,6 +150,9 @@ class csr_matrix_t { { } + // Get the number of nonzeros in row i + i_t row_length(i_t i) const { return row_start[i + 1] - row_start[i]; } + // Convert the CSR matrix to CSC i_t to_compressed_col(csc_matrix_t& Acol) const; diff --git a/cpp/src/dual_simplex/sparse_vector.cpp b/cpp/src/dual_simplex/sparse_vector.cpp index 0d34d4c390..1f3798e7dd 100644 --- a/cpp/src/dual_simplex/sparse_vector.cpp +++ b/cpp/src/dual_simplex/sparse_vector.cpp @@ -240,6 +240,15 @@ void sparse_vector_t::negate() } } +template +void sparse_vector_t::scale(f_t factor) +{ + const i_t nz = x.size(); + for (i_t k = 0; k < nz; ++k) { + x[k] *= factor; + } +} + template f_t sparse_vector_t::find_coefficient(i_t index) const { diff --git a/cpp/src/dual_simplex/sparse_vector.hpp b/cpp/src/dual_simplex/sparse_vector.hpp index 6ea642ee07..d9ca540d18 100644 --- a/cpp/src/dual_simplex/sparse_vector.hpp +++ b/cpp/src/dual_simplex/sparse_vector.hpp @@ -48,7 +48,11 @@ class sparse_vector_t { void sort(); // compute the squared 2-norm of the sparse vector f_t norm2_squared() const; + // negate the coefficients in the sparse vector void negate(); + // scale the coefficients in the sparse vector by a factor + void scale(f_t factor); + // find the coefficient of a given index f_t find_coefficient(i_t index) const; void clear() diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index 235d4500d2..2170c2437d 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -224,6 +224,10 @@ solution_t mip_solver_t::run_solver() branch_and_bound_settings.cut_min_orthogonality = context.settings.cut_min_orthogonality; branch_and_bound_settings.mip_batch_pdlp_strong_branching = context.settings.mip_batch_pdlp_strong_branching; + branch_and_bound_settings.reduced_cost_strengthening = + context.settings.reduced_cost_strengthening == -1 + ? 2 + : context.settings.reduced_cost_strengthening; if (context.settings.num_cpu_threads < 0) { branch_and_bound_settings.num_threads = std::max(1, omp_get_max_threads() - 1);