Now extracting primal phase 1 duals for sensitivity test

ERGO-Code · Jun 18, 2024 · a463332 · a463332
1 parent f0e03cd
commit a463332
Show file tree

Hide file tree

Showing 5 changed files with 84 additions and 17 deletions.
diff --git 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<double>& 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
@@ -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<double> 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
@@ -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
@@ -210,6 +210,7 @@ class HEkk {
   double debug_max_relative_dual_steepest_edge_weight_error = 0;
 
   std::vector<HighsSimplexBadBasisChangeRecord> bad_basis_change_;
+  std::vector<double> primal_phase1_dual_;
 
  private:
   bool isUnconstrainedLp();

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