Add IIS

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) {

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,

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<double> col_lower;
     std::vector<double> 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<HighsInt> ecol_index;
     std::vector<double> 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<bool> 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",