Incorporated integral scaling into mipIsKnapsack

Author: jajhall

check/TestMipSolver.cpp

@@ -1011,9 +1011,11 @@ TEST_CASE("knapsack", "[highs_test_mip_solver]") {
   const double offset = -100;
   std::vector<double> neg_value(num_item);
   std::vector<double> neg_weight(num_item);
+  std::vector<double> half_weight(num_item);
   for (HighsInt iItem = 0; iItem < num_item; iItem++) {
     neg_value[iItem] = -value[iItem];
     neg_weight[iItem] = -weight[iItem];
+    half_weight[iItem] = 0.5 * weight[iItem];
   }
   HighsLp lp;
   lp.sense_ = ObjSense::kMaximize;
@@ -1076,6 +1078,34 @@ TEST_CASE("knapsack", "[highs_test_mip_solver]") {
   REQUIRE(h.passModel(lp) == HighsStatus::kOk);
   REQUIRE(h.run() == HighsStatus::kOk);
   REQUIRE(h.getModelStatus() == HighsModelStatus::kOptimal);
+  REQUIRE(h.getInfo().objective_function_value == required_objective_value);
+
+  lp.a_matrix_.value_ = half_weight;
+  lp.row_upper_ = {0.5 * capacity};
+  REQUIRE(h.passModel(lp) == HighsStatus::kOk);
+  REQUIRE(h.run() == HighsStatus::kOk);
+  REQUIRE(h.getModelStatus() == HighsModelStatus::kOptimal);
+  REQUIRE(h.getInfo().objective_function_value == required_objective_value);
+
+  lp.a_matrix_.value_ = {1.1, 2.25, 3.05, 6.25, 7, 4.125};
+  lp.row_upper_ = {capacity};
+  required_objective_value = -190;
+  REQUIRE(h.passModel(lp) == HighsStatus::kOk);
+  REQUIRE(h.run() == HighsStatus::kOk);
+  REQUIRE(h.getModelStatus() == HighsModelStatus::kOptimal);
+  REQUIRE(h.getInfo().objective_function_value == required_objective_value);
+
+  lp.a_matrix_.value_ = {1.1334, 2.2501, 3.0534, 6.251, 7, 4.125};
+  REQUIRE(h.passModel(lp) == HighsStatus::kOk);
+  REQUIRE(h.run() == HighsStatus::kOk);
+  REQUIRE(h.getModelStatus() == HighsModelStatus::kOptimal);
+  REQUIRE(h.getInfo().objective_function_value == required_objective_value);
+
+  lp.a_matrix_.value_[1] = 2.250001;
+  REQUIRE(h.passModel(lp) == HighsStatus::kOk);
+  REQUIRE(h.run() == HighsStatus::kOk);
+  REQUIRE(h.getModelStatus() == HighsModelStatus::kOptimal);
+  REQUIRE(h.getInfo().objective_function_value == required_objective_value);
 
   h.resetGlobalScheduler(true);
 }

highs/mip/HighsMipSolver.cpp

@@ -91,10 +91,10 @@ void HighsMipSolver::run() {
   // Determine whether this is a knapsack problem and, if so, at least
   // update the data on knapsack sub-MIPs
   HighsInt capacity = 0;
-  if (mipdata_->mipIsKnapsack(capacity)) {
+  if (mipdata_->mipIsKnapsack()) {
     mipdata_->knapsack_data_.num_problem++;
     mipdata_->knapsack_data_.sum_variables += orig_model_->num_col_;
-    mipdata_->knapsack_data_.sum_capacity += capacity;
+    mipdata_->knapsack_data_.sum_capacity += mipdata_->knapsack_capacity_;
     // Solve as a knapsack
     HighsStatus call_status = mipdata_->heuristics.solveMipKnapsack();
     assert(call_status == HighsStatus::kOk);

highs/mip/HighsMipSolverData.cpp

@@ -685,6 +685,8 @@ void HighsMipSolverData::init() {
   optimality_limit = mipsolver.options_mip_->objective_bound;
   primal_dual_integral.initialise();
   knapsack_data_.initialise();
+  knapsack_capacity_ = 0;
+  knapsack_integral_scale_ = 0;
 
   if (mipsolver.options_mip_->mip_report_level == 0)
     dispfreq = 0;
@@ -2651,7 +2653,7 @@ void HighsMipSolverData::callbackUserSolution(
   }
 }
 
-bool HighsMipSolverData::mipIsKnapsack(HighsInt& capacity) {
+bool HighsMipSolverData::mipIsKnapsack(const bool logging) {
   const HighsLp& lp = *(mipsolver.model_);
   // Has to have one constraint
   if (lp.num_row_ != 1) return false;
@@ -2662,18 +2664,24 @@ bool HighsMipSolverData::mipIsKnapsack(HighsInt& capacity) {
     return false;
   const bool upper = lp.row_upper_[0] < kHighsInf;
   const HighsInt constraint_sign = upper ? 1 : -1;
-  // Now check that all the (signed) coefficients are integer and non-negative
-  for (HighsInt iEl = 0; iEl < lp.a_matrix_.numNz(); iEl++) {
-    double coeff = constraint_sign * lp.a_matrix_.value_[iEl];
-    if (coeff < 0) return false;
-    if (fractionality(coeff) > 0) return false;
-  }
+  // Now check that all the (signed) coefficients are non-negative
+  for (HighsInt iEl = 0; iEl < lp.a_matrix_.numNz(); iEl++)
+    if (constraint_sign * lp.a_matrix_.value_[iEl] < 0) return false;
+  this->knapsack_integral_scale_ =
+      HighsIntegers::integralScale(lp.a_matrix_.value_, 1e-6, 1e-6);
+  if (this->knapsack_integral_scale_ == 0) return false;
+  if (this->knapsack_integral_scale_ > 1000000) return false;
+  if (logging)
+    highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kInfo,
+                 "MIP is a knapsack problem with with scale %d\n",
+                 int(this->knapsack_integral_scale_));
   // Capacity must be integer, but OK to round down any fractional
   // values since activity of constraint is integer
   double double_capacity =
       upper ? lp.row_upper_[0] : constraint_sign * lp.row_lower_[0];
+  double_capacity *= this->knapsack_integral_scale_;
   const double capacity_margin = 1e-6;
-  capacity = std::floor(double_capacity + capacity_margin);
+  this->knapsack_capacity_ = std::floor(double_capacity + capacity_margin);
   // Problem is knapsack!
   return true;
 }

highs/mip/HighsMipSolverData.h

@@ -162,6 +162,8 @@ struct HighsMipSolverData {
   HighsDebugSol debugSolution;
 
   HighsKnapsackData knapsack_data_;
+  HighsInt knapsack_capacity_;
+  HighsInt knapsack_integral_scale_;
 
   HighsMipSolverData(HighsMipSolver& mipsolver)
       : mipsolver(mipsolver),
@@ -220,7 +222,9 @@ struct HighsMipSolverData {
         upper_bound(kHighsInf),
         upper_limit(kHighsInf),
         optimality_limit(kHighsInf),
-        debugSolution(mipsolver) {
+        debugSolution(mipsolver),
+        knapsack_capacity_(0),
+        knapsack_integral_scale_(0) {
     domain.addCutpool(cutpool);
     domain.addConflictPool(conflictPool);
   }
@@ -309,7 +313,7 @@ struct HighsMipSolverData {
   void callbackUserSolution(
       const double mipsolver_objective_value,
       const userMipSolutionCallbackOrigin user_solution_callback_origin);
-  bool mipIsKnapsack(HighsInt& capacity);
+  bool mipIsKnapsack(const bool logging = true);
 };
 
 #endif

highs/mip/HighsPrimalHeuristics.cpp

@@ -1750,10 +1750,7 @@ HighsStatus HighsPrimalHeuristics::solveMipKnapsackReturn(
     const HighsStatus& return_status) {
   const HighsLp& lp = *(mipsolver.model_);
   std::stringstream ss;
-  if (!mipsolver.submip) {
-    ss.str(std::string());
-    ss << "MIP is a knapsack problem: ";
-  }
+  if (!mipsolver.submip) ss.str(std::string());
   if (mipsolver.modelstatus_ == HighsModelStatus::kOptimal) {
     // mipsolver.solution_objective_ is the objective value for the
     // original problem - using the offset and ignoring the
@@ -1771,12 +1768,12 @@ HighsStatus HighsPrimalHeuristics::solveMipKnapsackReturn(
     mipsolver.gap_ = 0;
     if (!mipsolver.submip)
       ss << highsFormatToString(
-          "optimal objective by dynamic programming is %g",
+          "Optimal objective by dynamic programming is %g",
           mipsolver.solution_objective_);
   } else {
     if (!mipsolver.submip)
       ss << highsFormatToString(
-          "model status is %s",
+          "Model status is %s",
           utilModelStatusToString(mipsolver.modelstatus_).c_str());
     mipsolver.solution_.clear();
   }
@@ -1789,16 +1786,13 @@ HighsStatus HighsPrimalHeuristics::solveMipKnapsackReturn(
 HighsStatus HighsPrimalHeuristics::solveMipKnapsack() {
   const HighsLp& lp = *(mipsolver.model_);
   const HighsLogOptions& log_options = mipsolver.options_mip_->log_options;
-  HighsInt capacity_;
-  assert(mipsolver.mipdata_->mipIsKnapsack(capacity_));
+  assert(mipsolver.mipdata_->mipIsKnapsack(false));
 
   const bool upper = lp.row_upper_[0] < kHighsInf;
   const HighsInt constraint_sign = upper ? 1 : -1;
-  double double_capacity =
-      upper ? lp.row_upper_[0] : constraint_sign * lp.row_lower_[0];
-  const double capacity_margin = 1e-6;
-  const HighsInt capacity = std::floor(double_capacity + capacity_margin);
-
+  const HighsInt integral_scale = mipsolver.mipdata_->knapsack_integral_scale_;
+  assert(integral_scale > 0);
+  const HighsInt capacity = mipsolver.mipdata_->knapsack_capacity_;
   if (capacity < 0) {
     mipsolver.modelstatus_ = HighsModelStatus::kInfeasible;
     return solveMipKnapsackReturn(HighsStatus::kOk);
@@ -1816,7 +1810,8 @@ HighsStatus HighsPrimalHeuristics::solveMipKnapsack() {
   for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) {
     for (HighsInt iEl = lp.a_matrix_.start_[iCol];
          iEl < lp.a_matrix_.start_[iCol + 1]; iEl++)
-      weight[iCol] = HighsInt(constraint_sign * lp.a_matrix_.value_[iEl]);
+      weight[iCol] = HighsInt(integral_scale *
+                              (constraint_sign * lp.a_matrix_.value_[iEl]));
   }
   HighsInt sense = HighsInt(lp.sense_);
   // Set up the values for the knapsack solver. Since it solves a