Integrated knapsack solver: now debug its use

Author: jajhall

highs/lp_data/HighsLp.cpp

@@ -37,7 +37,7 @@ bool HighsLp::isBip() const {
   return true;
 }
 
-bool HighsLp::isKnapsack(double& knapsack_rhs) const {
+bool HighsLp::isKnapsack(HighsInt& knapsack_rhs) const {
   // Has to kave one constraint
   if (this->num_row_ != 1) return false;
   // Has to be a binary integer programming problem

highs/lp_data/HighsLp.h

@@ -66,7 +66,7 @@ class HighsLp {
   bool equalScaling(const HighsLp& lp) const;
   bool isMip() const;
   bool isBip() const;
-  bool isKnapsack(double& knapsack_rhs) const;
+  bool isKnapsack(HighsInt& knapsack_rhs) const;
   bool hasSemiVariables() const;
   bool hasInfiniteCost(const double infinite_cost) const;
   bool hasMods() const;

highs/mip/HighsMipSolver.cpp

@@ -90,11 +90,16 @@ void HighsMipSolver::run() {
 
   // Determine whether this is a knapsack problem and, if so, at least
   // update the data on knapsack sub-MIPs
-  double knapsack_rhs = 0;
+  HighsInt knapsack_rhs = 0;
   if (orig_model_->isKnapsack(knapsack_rhs)) {
     mipdata_->knapsack_data_.num_problem++;
     mipdata_->knapsack_data_.sum_variables += orig_model_->num_col_;
     mipdata_->knapsack_data_.sum_rhs += knapsack_rhs;
+    // Solve as a knapsack
+    HighsStatus call_status = mipdata_->heuristics.solveKnapsack();
+    printf("HighsMipSolver: Knapsack solver return status is %d\n", int(call_status));
+    cleanupSolve();
+    return;
   }
 
   analysis_.mipTimerStart(kMipClockRunPresolve);

highs/mip/HighsPrimalHeuristics.cpp

@@ -1646,3 +1646,110 @@ void HighsPrimalHeuristics::flushStatistics() {
   mipsolver.mipdata_->total_lp_iterations += lp_iterations;
   lp_iterations = 0;
 }
+
+double knapsackRecurrence(const std::vector<HighsInt>& weights,
+			  const std::vector<double>& values,
+			  const HighsInt num_col,
+			  const double rhs,
+			  std::vector<std::vector<double>> &dp,
+			  std::vector<std::vector<bool>> &use) {
+    if (num_col == 0 || rhs == 0)
+        return 0;  // Base case
+
+    if (dp[num_col][rhs] != -1) return dp[num_col][rhs];  // Check if result is already computed
+
+    // Exclude the item
+    double exclude = knapsackRecurrence(weights, values, num_col - 1, rhs, dp, use);
+
+    // Include the item (if it fits in the knapsack)
+    double include = 0;
+    if (weights[num_col - 1] <= rhs)
+      include = values[num_col - 1] + knapsackRecurrence(weights, values, num_col - 1, rhs - weights[num_col - 1], dp, use);
+
+    use[num_col][rhs] = include > exclude;
+    dp[num_col][rhs] = use[num_col][rhs] ? include : exclude;  // Store the result
+    
+    return dp[num_col][rhs];
+
+}
+HighsStatus HighsPrimalHeuristics::solveKnapsack() {
+  HighsInt knapsack_rhs_;
+  assert(mipsolver.orig_model_->isKnapsack(knapsack_rhs_));
+  const HighsInt knapsack_rhs = knapsack_rhs_;
+  
+  HighsLp lp = *(mipsolver.orig_model_);
+  //  const HighsLp& lp = mipsolver.mipdata_->orig_model_;
+  const bool upper = lp.row_upper_[0] < kHighsInf;
+  const HighsInt sign = upper ? 1 : -1;
+  if (knapsack_rhs < 0) {
+    mipsolver.modelstatus_ = HighsModelStatus::kInfeasible;
+    return HighsStatus::kOk;
+  } else if (knapsack_rhs == 0) {
+    // Trivial knapsack with zero solution
+    mipsolver.solution_.assign(lp.num_col_, 0);
+    mipsolver.solution_objective_ = lp.offset_;
+    mipsolver.modelstatus_ = HighsModelStatus::kOptimal;
+    return HighsStatus::kOk;
+  }
+  std::vector<HighsInt> weights(lp.num_col_, 0);
+  for (HighsInt iEl = 0; iEl < lp.a_matrix_.numNz(); iEl++) {
+    HighsInt iCol = lp.a_matrix_.index_[iEl];
+    weights[iCol] = HighsInt(sign * lp.a_matrix_.value_[iEl]);
+  }
+  HighsInt sense = HighsInt(lp.sense_);
+  std::vector<double> values;
+  for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++)
+    values.push_back(sense*lp.col_cost_[iCol]);
+  std::vector<std::vector<double>> dp(lp.num_col_ + 1, std::vector<double>(knapsack_rhs + 1, -1));
+  std::vector<std::vector<bool>> use(lp.num_col_ + 1, std::vector<bool>(knapsack_rhs + 1, -1));
+  double knapsack_optimal_objective_value = knapsackRecurrence(weights, values, lp.num_col_, knapsack_rhs, dp, use);
+  // Deduce the solution
+  std::vector<HighsInt> knapsack_solution(lp.num_col_, 0);
+  HighsInt knapsack_solution_rhs = knapsack_rhs;
+  for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) {
+    if (use[iCol][knapsack_solution_rhs]) {
+      knapsack_solution[iCol] = 1;
+      knapsack_solution_rhs -= weights[iCol];
+    }
+  }
+  const double row_violation = std::max(0, -knapsack_solution_rhs);
+  const double rel_row_violation = row_violation / (1.0 * knapsack_rhs);
+
+  if (rel_row_violation > 1e-12) {
+    highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kError,
+		 "HighsPrimalHeuristics::solveKnapsack() Relative capacity violation of %g\n", rel_row_violation);
+   mipsolver.modelstatus_ = HighsModelStatus::kSolveError;
+   return HighsStatus::kError;
+  }
+  const double solution_objective = lp.offset_ + sense * knapsack_optimal_objective_value;
+
+  // Compute the objective directly
+  double check_objective = lp.offset_;
+  for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) 
+    check_objective += knapsack_solution[iCol]*lp.col_cost_[iCol];
+
+  double abs_dl_solution_objective = std::fabs(mipsolver.solution_objective_ - check_objective);
+  double rel_dl_solution_objective = abs_dl_solution_objective / (1.0 + std::fabs(mipsolver.solution_objective_));
+  if (rel_dl_solution_objective > 1e-12) {
+    highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kError,
+		 "HighsPrimalHeuristics::solveKnapsack() Relative optimal objective value mismatch of %g\n", rel_dl_solution_objective);
+   mipsolver.modelstatus_ = HighsModelStatus::kSolveError;
+   return HighsStatus::kError;
+  }
+  // Copy in the solution and compute the objective directly
+  mipsolver.solution_.resize(lp.num_col_);
+  for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) 
+    mipsolver.solution_[iCol] = knapsack_solution[iCol];
+  mipsolver.solution_objective_ = solution_objective;
+
+  mipsolver.bound_violation_ = 0;
+  mipsolver.integrality_violation_ = 0;
+  mipsolver.row_violation_ = 0;
+  mipsolver.dual_bound_ = mipsolver.solution_objective_;
+  mipsolver.primal_bound_ = mipsolver.solution_objective_;
+  mipsolver.gap_ = 0;
+
+  mipsolver.modelstatus_ = HighsModelStatus::kOptimal;
+  return HighsStatus::kOk;				      
+}
+

highs/mip/HighsPrimalHeuristics.h

@@ -69,6 +69,8 @@ class HighsPrimalHeuristics {
   void shifting(const std::vector<double>& relaxationsol);
 
   void ziRound(const std::vector<double>& relaxationsol);
+
+  HighsStatus solveKnapsack();
 };
 
 #endif