diff --git a/src/lp_data/HighsIis.cpp b/src/lp_data/HighsIis.cpp index 9f244415ca..50fef4e969 100644 --- a/src/lp_data/HighsIis.cpp +++ b/src/lp_data/HighsIis.cpp @@ -211,7 +211,7 @@ HighsStatus HighsIis::getData(const HighsLp& lp, const HighsOptions& options, this->col_index_[iCol] = from_col[this->col_index_[iCol]]; for (HighsInt iRow = 0; iRow < HighsInt(this->row_index_.size()); iRow++) this->row_index_[iRow] = from_row[this->row_index_[iRow]]; - if (kIisDevReport) this->report("On exit", lp); + if (kIisDevReportVerbose) this->report("On exit", lp); return HighsStatus::kOk; } @@ -226,7 +226,7 @@ HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options, for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++) this->addRow(iRow); Highs highs; const HighsInfo& info = highs.getInfo(); - highs.setOptionValue("output_flag", kIisDevReport); + highs.setOptionValue("output_flag", kIisDevReportVerbose); highs.setOptionValue("presolve", kHighsOffString); const HighsLp& incumbent_lp = highs.getLp(); const HighsBasis& incumbent_basis = highs.getBasis(); @@ -248,12 +248,36 @@ HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options, HighsInt iX = -1; bool drop_lower = false; + HighsInt num_run_call = 0; + const HighsInt check_run_call = 154; // kHighsIInf; + // Lambda for gathering data when solving an LP auto solveLp = [&]() -> HighsStatus { HighsIisInfo iis_info; iis_info.simplex_time = -highs.getRunTime(); iis_info.simplex_iterations = -info.simplex_iteration_count; + bool output_flag; + HighsInt log_dev_level; + HighsInt highs_analysis_level; + highs.getOptionValue("output_flag", output_flag); + highs.getOptionValue("log_dev_level", log_dev_level); + highs.getOptionValue("highs_analysis_level", highs_analysis_level); + + num_run_call++; + if (num_run_call == check_run_call) { + highs.setOptionValue("output_flag", true); + highs.setOptionValue("log_dev_level", 3); + highs.setOptionValue("highs_analysis_level", 4); + highs.writeModel("form.mps"); + } run_status = highs.run(); + highs.setOptionValue("output_flag", output_flag); + highs.setOptionValue("log_dev_level", log_dev_level); + highs.setOptionValue("highs_analysis_level", highs_analysis_level); + if (run_status != HighsStatus::kOk) { + printf("HighsIis::compute highs.run() %d returns status %s\n", + int(num_run_call), highsStatusToString(run_status).c_str()); + } assert(run_status == HighsStatus::kOk); if (run_status != HighsStatus::kOk) return run_status; HighsModelStatus model_status = highs.getModelStatus(); @@ -262,7 +286,6 @@ HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options, 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; @@ -360,6 +383,8 @@ HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options, // Pass twice: rows before columns, or columns before rows, according to // row_priority + std::string check_type = "Col"; + HighsInt check_ix = 81; for (HighsInt k = 0; k < 2; k++) { row_deletion = (row_priority && k == 0) || (!row_priority && k == 1); std::string type = row_deletion ? "Row" : "Col"; @@ -375,6 +400,10 @@ HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options, double upper = row_deletion ? lp.row_upper_[iX] : lp.col_upper_[iX]; // Record whether the upper bound has been dropped due to the lower bound // being kept + if (check_type == type && check_ix == iX) { + printf("CheckType %s, index %d, will be num_run_call = %d\n", + check_type.c_str(), int(check_ix), int(num_run_call + 1)); + } if (lower > -kHighsInf) { // Drop the lower bound temporarily bool drop_lower = true; @@ -489,7 +518,7 @@ HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options, } if (k == 1) continue; // End of first pass: look to simplify second pass - if (kIisDevReport) this->report("End of deletion", incumbent_lp); + if (kIisDevReportVerbose) this->report("End of deletion", incumbent_lp); if (row_deletion) { // Mark empty columns as dropped for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) { @@ -511,9 +540,9 @@ HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options, } } } - if (kIisDevReport) this->report("End of pass 1", incumbent_lp); + if (kIisDevReportVerbose) this->report("End of pass 1", incumbent_lp); } - if (kIisDevReport) this->report("End of pass 2", incumbent_lp); + if (kIisDevReportVerbose) this->report("End of pass 2", incumbent_lp); HighsInt iss_num_col = 0; for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) { if (this->col_bound_[iCol] != kIisBoundStatusDropped) { diff --git a/src/lp_data/HighsIis.h b/src/lp_data/HighsIis.h index cb01ea3500..65fdb1fa5b 100644 --- a/src/lp_data/HighsIis.h +++ b/src/lp_data/HighsIis.h @@ -16,7 +16,8 @@ #include "lp_data/HighsLp.h" -const bool kIisDevReport = false; +const bool kIisDevReportBrief = true; +const bool kIisDevReportVerbose = false; enum IisBoundStatus { kIisBoundStatusDropped = -1, diff --git a/src/lp_data/HighsInterface.cpp b/src/lp_data/HighsInterface.cpp index 6176c39bf1..3eded75ac8 100644 --- a/src/lp_data/HighsInterface.cpp +++ b/src/lp_data/HighsInterface.cpp @@ -1821,13 +1821,15 @@ HighsStatus Highs::elasticityFilter( const bool has_elastic_rows = has_local_rhs_penalty || has_global_elastic_rhs; assert(has_elastic_columns || has_elastic_rows); + const bool has_col_names = lp.col_names_.size() > 0; + const bool has_row_names = lp.row_names_.size() > 0; + const HighsInt col_ecol_offset = lp.num_col_; if (has_elastic_columns) { // Accumulate bounds to be used for columns std::vector col_lower; std::vector col_upper; // When defining names, need to know the column number - const bool has_col_names = lp.col_names_.size() > 0; erow_start.push_back(0); for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) { const double lower = lp.col_lower_[iCol]; @@ -1901,7 +1903,7 @@ HighsStatus Highs::elasticityFilter( HighsInt num_new_col = col_of_ecol.size(); HighsInt num_new_row = erow_start.size() - 1; HighsInt num_new_nz = erow_start[num_new_row]; - if (kIisDevReport) + if (kIisDevReportBrief) printf( "Elasticity filter: For columns there are %d variables and %d " "constraints\n", @@ -1954,7 +1956,6 @@ HighsStatus Highs::elasticityFilter( std::vector ecol_index; std::vector ecol_value; ecol_start.push_back(0); - const bool has_row_names = lp.row_names_.size() > 0; for (HighsInt iRow = 0; iRow < original_num_row; iRow++) { // Get the penalty for violating the bounds on this row const double penalty = @@ -2044,7 +2045,7 @@ HighsStatus Highs::elasticityFilter( original_num_row, original_col_cost, original_col_lower, original_col_upper, original_integrality); - if (kIisDevReport) this->writeSolution("", kSolutionStylePretty); + if (kIisDevReportVerbose) this->writeSolution("", kSolutionStylePretty); // Model status should be optimal, unless model is unbounded assert(this->model_status_ == HighsModelStatus::kOptimal || this->model_status_ == HighsModelStatus::kUnbounded); @@ -2058,49 +2059,61 @@ HighsStatus Highs::elasticityFilter( // Now fix e-variables that are positive and re-solve until e-LP is infeasible HighsInt loop_k = 0; bool feasible_model = false; + // Use strict zero solution value rather than within tolerance since + // cplex2 is almost feasible, and only requires an elastic variable + // value of 8.87022e-10 to be feasible. + const double use_primal_tolerance = 0; for (;;) { - if (kIisDevReport) + if (kIisDevReportBrief) printf("\nElasticity filter pass %d\n==============\n", int(loop_k)); + // An elastic variable can be fixed at zero, but have positive + // value (within the tolerance) if basic, so allowing it to be + // re-fixed can cause an infinite loop. Happens with cplex2 when + // fixed elastic variable is basic at 8.87022e-10 HighsInt num_fixed = 0; if (has_elastic_columns) { for (HighsInt eCol = 0; eCol < col_of_ecol.size(); eCol++) { - HighsInt iCol = col_of_ecol[eCol]; - if (solution.col_value[col_ecol_offset + eCol] > - this->options_.primal_feasibility_tolerance) { - if (kIisDevReport) + const HighsInt iCol = col_ecol_offset + eCol; + if (lp.col_upper_[iCol] == 0) continue; + const HighsInt original_col = col_of_ecol[eCol]; + if (solution.col_value[iCol] > use_primal_tolerance) { + if (kIisDevReportBrief) printf( "E-col %2d (column %2d) corresponds to column %2d with bound " "%g " "and has solution value %g\n", - int(eCol), int(col_ecol_offset + eCol), int(iCol), - bound_of_col_of_ecol[eCol], - solution.col_value[col_ecol_offset + eCol]); - this->changeColBounds(col_ecol_offset + eCol, 0, 0); + int(eCol), int(iCol), int(original_col), + bound_of_col_of_ecol[eCol], solution.col_value[iCol]); + this->changeColBounds(iCol, 0, 0); num_fixed++; } } } if (has_elastic_rows) { for (HighsInt eCol = 0; eCol < row_of_ecol.size(); eCol++) { - HighsInt iRow = row_of_ecol[eCol]; - if (solution.col_value[row_ecol_offset + eCol] > - this->options_.primal_feasibility_tolerance) { - if (kIisDevReport) + const HighsInt iCol = row_ecol_offset + eCol; + if (lp.col_upper_[iCol] == 0) continue; + const HighsInt original_row = row_of_ecol[eCol]; + if (solution.col_value[iCol] > use_primal_tolerance) { + if (kIisDevReportBrief) printf( "E-row %2d (column %2d) corresponds to row %2d with bound " "%g " "and has solution value %g\n", - int(eCol), int(row_ecol_offset + eCol), int(iRow), - bound_of_row_of_ecol[eCol], - solution.col_value[row_ecol_offset + eCol]); - this->changeColBounds(row_ecol_offset + eCol, 0, 0); + int(eCol), int(iCol), int(original_row), + bound_of_row_of_ecol[eCol], solution.col_value[iCol]); + this->changeColBounds(iCol, 0, 0); num_fixed++; } } } if (num_fixed == 0) { - // No elastic variables were positive, so problem is feasible - feasible_model = true; + // No new elastic variables were fixed, so break. If this was + // the first pass, then the original problem is feasible. If the + // original problem is feasible, its model status cannot be set + // so, according to the truth of feasible_model, this will be + // done in elasticityFilterReturn. + feasible_model = loop_k == 0; break; } HighsStatus run_status = solveLp(); @@ -2109,7 +2122,7 @@ HighsStatus Highs::elasticityFilter( original_num_col, original_num_row, original_col_cost, original_col_lower, original_col_upper, original_integrality); - if (kIisDevReport) this->writeSolution("", kSolutionStylePretty); + if (kIisDevReportVerbose) this->writeSolution("", kSolutionStylePretty); HighsModelStatus model_status = this->getModelStatus(); if (model_status == HighsModelStatus::kInfeasible) break; loop_k++; @@ -2120,29 +2133,39 @@ HighsStatus Highs::elasticityFilter( HighsInt num_enforced_row_ecol = 0; if (has_elastic_columns) { for (HighsInt eCol = 0; eCol < col_of_ecol.size(); eCol++) { - HighsInt iCol = col_of_ecol[eCol]; - if (lp.col_upper_[col_ecol_offset + eCol] == 0) { + const HighsInt original_col = col_of_ecol[eCol]; + const HighsInt iCol = col_ecol_offset + eCol; + if (lp.col_upper_[iCol] == 0) { num_enforced_col_ecol++; printf( - "Col e-col %2d (column %2d) corresponds to column %2d with bound " - "%g " + "Col e-col %4d (column %4d) corresponds to column %4d (%8s) with " + "bound " + "%11.4g " "and is enforced\n", - int(eCol), int(col_ecol_offset + eCol), int(iCol), + int(eCol), int(iCol), int(original_col), + has_col_names ? lp.col_names_[original_col].c_str() : "", bound_of_col_of_ecol[eCol]); } } } if (has_elastic_rows) { + std::vector in_infeasible_row_subset; + in_infeasible_row_subset.assign(original_num_row, false); for (HighsInt eCol = 0; eCol < row_of_ecol.size(); eCol++) { - HighsInt iRow = row_of_ecol[eCol]; - if (lp.col_upper_[row_ecol_offset + eCol] == 0) { + const HighsInt original_row = row_of_ecol[eCol]; + assert(original_row < original_num_row); + const HighsInt iCol = row_ecol_offset + eCol; + if (lp.col_upper_[iCol] == 0) { + assert(!in_infeasible_row_subset[original_row]); num_enforced_row_ecol++; - infeasible_row_subset.push_back(iRow); - if (kIisDevReport) + infeasible_row_subset.push_back(original_row); + if (kIisDevReportBrief) printf( - "Row e-col %2d (column %2d) corresponds to row %2d with bound " - "%g and is enforced\n", - int(eCol), int(row_ecol_offset + eCol), int(iRow), + "Row e-col %4d (column %4d) corresponds to row %4d (%8s) with " + "bound %11.4g " + "and is enforced\n", + int(eCol), int(iCol), int(original_row), + has_row_names ? lp.row_names_[original_row].c_str() : "", bound_of_row_of_ecol[eCol]); } } @@ -2156,7 +2179,7 @@ HighsStatus Highs::elasticityFilter( "rows\n", int(loop_k), int(num_enforced_col_ecol), int(num_enforced_row_ecol)); - if (kIisDevReport) + if (kIisDevReportBrief) printf( "\nElasticity filter after %d passes enforces bounds on %d cols and %d " "rows\n",