diff --git a/src/Highs.h b/src/Highs.h index 6f7152d77a..11053f2282 100644 --- a/src/Highs.h +++ b/src/Highs.h @@ -1217,6 +1217,14 @@ class Highs { HighsStatus getBasisInverseRowSparse(const HighsInt row, HVector& row_ep_buffer); + /** + * @Brief Get the primal simplex phase 1 dual values. Advanced + * method: for HiGHS IIS calculation + */ + const std::vector& getPrimalPhase1Dual() const { + return ekk_instance_.primal_phase1_dual_; + } + // Start of deprecated methods std::string compilationDate() const { return "deprecated"; } diff --git a/src/lp_data/HighsIis.cpp b/src/lp_data/HighsIis.cpp index 3ea491c57d..ec641467d6 100644 --- a/src/lp_data/HighsIis.cpp +++ b/src/lp_data/HighsIis.cpp @@ -227,6 +227,8 @@ HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options, highs.setOptionValue("output_flag", kIisDevReport); highs.setOptionValue("presolve", kHighsOffString); const HighsLp& incumbent_lp = highs.getLp(); + const HighsBasis& incumbent_basis = highs.getBasis(); + const HighsSolution& solution = highs.getSolution(); HighsStatus run_status = highs.passModel(lp); assert(run_status == HighsStatus::kOk); if (basis) highs.setBasis(*basis); @@ -238,6 +240,11 @@ HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options, assert(run_status == HighsStatus::kOk); // Solve the LP if (basis) highs.setBasis(*basis); + const bool use_sensitivity_filter = false; + std::vector primal_phase1_dual; + bool row_deletion = false; + HighsInt iX = -1; + bool drop_lower = false; // Lambda for gathering data when solving an LP auto solveLp = [&]() -> HighsStatus { @@ -247,6 +254,64 @@ HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options, run_status = highs.run(); assert(run_status == HighsStatus::kOk); if (run_status != HighsStatus::kOk) return run_status; + HighsModelStatus model_status = highs.getModelStatus(); + if (use_sensitivity_filter && + model_status == HighsModelStatus::kInfeasible) { + printf("\nHighsIis::compute %s deletion for %d and %s bound\n", + row_deletion ? "Row" : "Col", int(iX), + drop_lower ? "Lower" : "Upper"); + bool output_flag; + highs.getOptionValue("output_flag", output_flag); + highs.setOptionValue("output_flag", true); + HighsInt simplex_strategy; + highs.getOptionValue("simplex_strategy", simplex_strategy); + highs.setOptionValue("simplex_strategy", kSimplexStrategyPrimal); + // Solve the LP + run_status = highs.run(); + if (run_status != HighsStatus::kOk) return run_status; + highs.writeSolution("", kSolutionStylePretty); + primal_phase1_dual = highs.getPrimalPhase1Dual(); + HighsInt num_zero_dual = 0; + for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) { + const HighsBasisStatus status = incumbent_basis.col_status[iCol]; + const double dual = primal_phase1_dual[iCol]; + const double lower = lp.col_lower_[iCol]; + const double upper = lp.col_upper_[iCol]; + const double value = solution.col_value[iCol]; + if (status != HighsBasisStatus::kBasic && + std::fabs(dual) < options.dual_feasibility_tolerance) { + num_zero_dual++; + // Small dual for nonbasic variable + printf( + "HighsIis::compute Column %d [%g, %g, %g] with status %s has " + "dual %g\n", + int(iCol), lower, value, upper, + highs.basisStatusToString(status).c_str(), dual); + // assert(123 == 456); + } + } + for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++) { + const HighsBasisStatus status = incumbent_basis.row_status[iRow]; + const double dual = primal_phase1_dual[lp.num_col_ + iRow]; + const double lower = lp.row_lower_[iRow]; + const double upper = lp.row_upper_[iRow]; + const double value = solution.row_value[iRow]; + if (status != HighsBasisStatus::kBasic && + std::fabs(dual) < options.dual_feasibility_tolerance) { + num_zero_dual++; + // Small dual for nonbasic variable + printf( + "HighsIis::compute Row %d [%g, %g, %g] with status %s has " + "dual %g\n", + int(iRow), lower, value, upper, + highs.basisStatusToString(status).c_str(), dual); + // assert(123 == 456); + } + } + highs.setOptionValue("output_flag", output_flag); + highs.setOptionValue("simplex_strategy", simplex_strategy); + assert(!num_zero_dual); + } iis_info.simplex_time += highs.getRunTime(); iis_info.simplex_iterations += info.simplex_iteration_count; this->info_.push_back(iis_info); @@ -258,29 +323,14 @@ HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options, assert(highs.getModelStatus() == HighsModelStatus::kInfeasible); - const bool use_sensitivity_filter = false; - if (use_sensitivity_filter) { - bool output_flag; - highs.getOptionValue("output_flag", output_flag); - highs.setOptionValue("simplex_strategy", kSimplexStrategyPrimal); - // - highs.setOptionValue("output_flag", true); - // Solve the LP - run_status = solveLp(); - if (run_status != HighsStatus::kOk) return run_status; - highs.writeSolution("", kSolutionStylePretty); - highs.setOptionValue("output_flag", output_flag); - } - // Pass twice: rows before columns, or columns before rows, according to // row_priority for (HighsInt k = 0; k < 2; k++) { - const bool row_deletion = - (row_priority && k == 0) || (!row_priority && k == 1); + row_deletion = (row_priority && k == 0) || (!row_priority && k == 1); std::string type = row_deletion ? "Row" : "Col"; // Perform deletion pass HighsInt num_index = row_deletion ? lp.num_row_ : lp.num_col_; - for (HighsInt iX = 0; iX < num_index; iX++) { + for (iX = 0; iX < num_index; iX++) { const HighsInt ix_status = row_deletion ? this->row_bound_[iX] : this->col_bound_[iX]; if (ix_status == kIisBoundStatusDropped || @@ -292,6 +342,7 @@ HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options, // being kept if (lower > -kHighsInf) { // Drop the lower bound temporarily + bool drop_lower = true; run_status = row_deletion ? highs.changeRowBounds(iX, -kHighsInf, upper) : highs.changeColBounds(iX, -kHighsInf, upper); diff --git a/src/simplex/HEkk.cpp b/src/simplex/HEkk.cpp index b94d0959f6..459bdd7a52 100644 --- a/src/simplex/HEkk.cpp +++ b/src/simplex/HEkk.cpp @@ -161,6 +161,8 @@ void HEkk::clearEkkData() { this->debug_max_relative_dual_steepest_edge_weight_error = 0; clearBadBasisChange(); + + this->primal_phase1_dual_.clear(); } void HEkk::clearEkkDataInfo() { diff --git a/src/simplex/HEkk.h b/src/simplex/HEkk.h index c91aba199a..27e927a67d 100644 --- a/src/simplex/HEkk.h +++ b/src/simplex/HEkk.h @@ -210,6 +210,7 @@ class HEkk { double debug_max_relative_dual_steepest_edge_weight_error = 0; std::vector bad_basis_change_; + std::vector primal_phase1_dual_; private: bool isUnconstrainedLp(); diff --git a/src/simplex/HEkkPrimal.cpp b/src/simplex/HEkkPrimal.cpp index 8c9787cee4..1bf52aa439 100644 --- a/src/simplex/HEkkPrimal.cpp +++ b/src/simplex/HEkkPrimal.cpp @@ -253,6 +253,11 @@ HighsStatus HEkkPrimal::solve(const bool pass_force_phase2) { // LP identified as not having an optimal solution assert(ekk_instance_.model_status_ == HighsModelStatus::kInfeasible || ekk_instance_.model_status_ == HighsModelStatus::kUnbounded); + // If infeasible, save the primal phase 1 dual values before + // they are overwritten with the duals for the original + // objective + if (ekk_instance_.model_status_ == HighsModelStatus::kInfeasible) + ekk_instance_.primal_phase1_dual_ = ekk_instance_.info_.workDual_; break; } if (solve_phase == kSolvePhaseOptimalCleanup) {