Adds StrongCG cuts

highs/mip/HighsCutGeneration.cpp

@@ -500,6 +500,8 @@ bool HighsCutGeneration::separateLiftedMixedIntegerCover() {
   return true;
 }
 
+static double fast_ceil(double x) { return (int64_t)x + (x > (int64_t)x); }
+
 static double fast_floor(double x) { return (int64_t)x - (x < (int64_t)x); }
 
 bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy,
@@ -521,6 +523,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy,
   integerinds.clear();
   integerinds.reserve(rowlen);
 
+  bool strongcg = !onlyInitialCMIRScale;
   double maxabsdelta = 0.0;
   constexpr double maxCMirScale = 1e6;
   constexpr double f0min = 0.005;
@@ -542,9 +545,12 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy,
         maxabsdelta = max(maxabsdelta, delta);
         deltas.push_back(delta);
       }
-    } else {
+    } else if (vals[i] < 0) {
+      // Positive coefficient values later set to 0 so have no contribution
       updateViolationAndNorm(i, vals[i], continuouscontribution,
                              continuoussqrnorm);
+      // StrongCG cannot be computed when negative coefficients for cont exist
+      strongcg = false;
     }
   }
 
@@ -687,13 +693,116 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy,
     }
   }
 
+  // Calc strongcg cut as potential alternative to CMIR
+  if (bestefficacy <= 0 || rowlen == 0) {
+    strongcg = false;
+  }
+  if (strongcg) {
+    double delta = bestdelta;
+    double scale = 1.0 / delta;
+    double scalrhs = double(rhs) * scale;
+    double downrhs = fast_floor(scalrhs);
+    double f0 = scalrhs - downrhs;
+    double oneoveroneminusf0 = 1.0 / (1.0 - f0);
+    // Skip numerically troublesome cuts
+    if (fractionality(1 / f0) < 1e-3) {
+      strongcg = false;
+    } else {
+      double k = fast_ceil(1 / f0) - 1;
+      // All coefficients of continuous variables are 0 in strong CG cut
+      double sqrnorm = 0;
+      double viol = -downrhs;
+
+      for (HighsInt j : integerinds) {
+        double scalaj = vals[j] * scale;
+        double downaj = fast_floor(scalaj + kHighsTiny);
+        double fj = scalaj - downaj;
+        double aj = downaj;
+        if (fj >= f0 + 10 * feastol) {
+          double pj = fast_ceil(k * (fj - f0) * oneoveroneminusf0 - 1e-4);
+          aj += pj / (k + 1);
+        }
+        updateViolationAndNorm(j, aj, viol, sqrnorm);
+      }
+      if (sqrnorm <= kHighsTiny) {
+        strongcg = false;
+      } else {
+        double efficacy = viol / sqrt(sqrnorm);
+        // Use the strongCG cut instead of the CMIR if efficacy is larger
+        if (efficacy < bestefficacy + epsilon) {
+          strongcg = false;
+        } else {
+          bestefficacy = efficacy;
+        }
+      }
+    }
+  }
+
+  if (strongcg) {
+    // try to flip complementation of integers to increase efficacy
+    double delta = bestdelta;
+    double scale = 1.0 / delta;
+    for (HighsInt i : integerinds) {
+      if (upper[i] == kHighsInf) continue;
+      if (solval[i] <= feastol) continue;
+
+      flipComplementation(i);
+
+      double scalrhs = double(rhs) * scale;
+      double downrhs = fast_floor(scalrhs);
+
+      double f0 = scalrhs - downrhs;
+      if (f0 < f0min || f0 > f0max) {
+        flipComplementation(i);
+        continue;
+      }
+
+      double oneoveroneminusf0 = 1.0 / (1.0 - f0);
+      if (oneoveroneminusf0 > maxCMirScale) {
+        flipComplementation(i);
+        continue;
+      }
+
+      if (fractionality(1 / f0) < 1e-3) {
+        flipComplementation(i);
+        continue;
+      }
+
+      double k = fast_ceil(1 / f0) - 1;
+
+      double sqrnorm = 0;
+      double viol = -downrhs;
+
+      for (HighsInt j : integerinds) {
+        double scalaj = vals[j] * scale;
+        double downaj = fast_floor(scalaj + kHighsTiny);
+        double fj = scalaj - downaj;
+        double aj = downaj;
+        if (fj - f0 >= 10 * feastol) {
+          double pj = fast_ceil(k * (fj - f0) * oneoveroneminusf0 - 1e-4);
+          aj += (pj / (k + 1));
+        }
+        updateViolationAndNorm(j, aj, viol, sqrnorm);
+      }
+
+      double efficacy = viol / sqrt(sqrnorm);
+      if (efficacy > bestefficacy) {
+        bestefficacy = efficacy;
+      } else {
+        flipComplementation(i);
+      }
+    }
+  }
+
   HighsCDouble scale = 1.0 / HighsCDouble(bestdelta);
   HighsCDouble scalrhs = rhs * scale;
   double downrhs = floor(double(scalrhs));
 
   HighsCDouble f0 = scalrhs - downrhs;
   HighsCDouble oneoveroneminusf0 = 1.0 / (1.0 - f0);
 
+  double k = strongcg ? ceil(static_cast<double>(1 / f0)) - 1 : 0;
+
   rhs = downrhs * bestdelta;
   integralSupport = true;
   integralCoefficients = false;
@@ -711,8 +820,16 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy,
       double downaj = floor(double(scalaj + kHighsTiny));
       HighsCDouble fj = scalaj - downaj;
       HighsCDouble aj = downaj;
-      if (fj > f0) aj += (fj - f0) * oneoveroneminusf0;
-
+      if (fj > f0) {
+        if (strongcg) {
+          if (fj - f0 > epsilon) {
+            HighsCDouble pj = ceil(k * (fj - f0) * oneoveroneminusf0 - feastol);
+            aj += pj / (k + 1);
+          }
+        } else {
+          aj += (fj - f0) * oneoveroneminusf0;
+        }
+      }
       vals[j] = double(aj * bestdelta);
     }
   }
@@ -728,8 +845,14 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy,
     }
     double checkefficacy = checkviol / sqrt(checknorm);
     // the efficacy can become infinite if the cut 0 <= -1 is derived
-    assert(fabs(checkefficacy - bestefficacy) < 0.001 ||
-           fabs(checkefficacy) >= 1e30);
+    if (!strongcg) {
+      assert(fabs(checkefficacy - bestefficacy) < 0.001 ||
+             fabs(checkefficacy) >= 1e30);
+    } else {
+      // Rounded conservatively for scoring => strongcg cut can be stronger
+      assert(bestefficacy - checkefficacy < 0.001 ||
+             fabs(checkefficacy) >= 1e30);
+    }
   }
 #endif