From ed17c2ed20cb851dc91ba1d229ee7d1852f844f5 Mon Sep 17 00:00:00 2001 From: Mark Turner Date: Fri, 6 Jun 2025 14:18:29 +0200 Subject: [PATCH 01/19] Add strongcg draft --- highs/mip/HighsCutGeneration.cpp | 53 ++++++++++++++++++++++++++++++-- 1 file changed, 51 insertions(+), 2 deletions(-) diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index f518b7d8009..711e5805fb3 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/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 = true; double maxabsdelta = 0.0; constexpr double maxCMirScale = 1e6; constexpr double f0min = 0.005; @@ -545,6 +548,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, } else { updateViolationAndNorm(i, vals[i], continuouscontribution, continuoussqrnorm); + if (vals[i] <= 0) strongcg = false; } } @@ -687,6 +691,42 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, } } + // If inequality has no neg coef on cont vars: Calc strongcg cut + 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); + // TODO: Should stronger safeguards be used, i.e. multiply delta by 2, if f0 is too large? + double k = fast_ceil(1 / f0) - 1; + + double contscale = scale * oneoveroneminusf0; + double sqrnorm = contscale * contscale * continuoussqrnorm; + double viol = contscale * continuouscontribution - downrhs; + + for (HighsInt j : integerinds) { + double scalaj = vals[j] * scale; + double downaj = fast_floor(scalaj + kHighsTiny); + double fj = scalaj - downaj; + if (fj <= f0) { + double aj = downaj; + updateViolationAndNorm(j, aj, viol, sqrnorm); + } + else { + double pj = fast_ceil(k * (fj - f0) * oneoveroneminusf0); + double aj = downaj + (pj / (k + 1)); + updateViolationAndNorm(j, aj, viol, sqrnorm); + } + } + double efficacy = viol / sqrt(sqrnorm); + // Use the strong CG cut instead of the CMIR if efficacy is larger + if (efficacy < bestefficacy + feastol) { + strongcg = false; + } + } + HighsCDouble scale = 1.0 / HighsCDouble(bestdelta); HighsCDouble scalrhs = rhs * scale; double downrhs = floor(double(scalrhs)); @@ -694,6 +734,8 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, HighsCDouble f0 = scalrhs - downrhs; HighsCDouble oneoveroneminusf0 = 1.0 / (1.0 - f0); + double k = strongcg ? ceil(static_cast(1 / f0)) : 0; + rhs = downrhs * bestdelta; integralSupport = true; integralCoefficients = false; @@ -711,8 +753,15 @@ 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) { + HighsCDouble pj = ceil(k * (fj - f0) * oneoveroneminusf0); + aj += pj / (k + 1); + } + else { + aj += (fj - f0) * oneoveroneminusf0; + } + } vals[j] = double(aj * bestdelta); } } From 213d1f71ef8823edc1bc37e5df2c7b2daea10159 Mon Sep 17 00:00:00 2001 From: Mark Turner Date: Fri, 6 Jun 2025 14:34:05 +0200 Subject: [PATCH 02/19] Add missing 1 --- highs/mip/HighsCutGeneration.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index 711e5805fb3..06a3060933e 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/highs/mip/HighsCutGeneration.cpp @@ -734,7 +734,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, HighsCDouble f0 = scalrhs - downrhs; HighsCDouble oneoveroneminusf0 = 1.0 / (1.0 - f0); - double k = strongcg ? ceil(static_cast(1 / f0)) : 0; + double k = strongcg ? ceil(static_cast(1 / f0)) - 1 : 0; rhs = downrhs * bestdelta; integralSupport = true; From 0eccd2c78bbc3c335dc5141e7aab4b0e1ef86626 Mon Sep 17 00:00:00 2001 From: Mark Turner Date: Fri, 6 Jun 2025 16:28:47 +0200 Subject: [PATCH 03/19] Correct viol for strongcg. --- highs/mip/HighsCutGeneration.cpp | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index 06a3060933e..9fbe4e49025 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/highs/mip/HighsCutGeneration.cpp @@ -545,10 +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); - if (vals[i] <= 0) strongcg = false; + strongcg = false; } } @@ -702,9 +704,9 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, // TODO: Should stronger safeguards be used, i.e. multiply delta by 2, if f0 is too large? double k = fast_ceil(1 / f0) - 1; - double contscale = scale * oneoveroneminusf0; - double sqrnorm = contscale * contscale * continuoussqrnorm; - double viol = contscale * continuouscontribution - downrhs; + // 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; @@ -715,16 +717,19 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, updateViolationAndNorm(j, aj, viol, sqrnorm); } else { - double pj = fast_ceil(k * (fj - f0) * oneoveroneminusf0); + double pj = fast_ceil(k * (fj - f0) * oneoveroneminusf0 - kHighsTiny); double aj = downaj + (pj / (k + 1)); updateViolationAndNorm(j, aj, viol, sqrnorm); } } double efficacy = viol / sqrt(sqrnorm); // Use the strong CG cut instead of the CMIR if efficacy is larger - if (efficacy < bestefficacy + feastol) { + if (efficacy < bestefficacy + 10 * feastol) { strongcg = false; } + else { + bestefficacy = efficacy; + } } HighsCDouble scale = 1.0 / HighsCDouble(bestdelta); From 7ea536423e9f12bd12178e288513558aaa7be778 Mon Sep 17 00:00:00 2001 From: Mark Turner Date: Fri, 6 Jun 2025 19:35:15 +0200 Subject: [PATCH 04/19] Add some tolerance to constraint improvement check --- highs/mip/HighsCutGeneration.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index 9fbe4e49025..8efdfaa7185 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/highs/mip/HighsCutGeneration.cpp @@ -712,7 +712,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, double scalaj = vals[j] * scale; double downaj = fast_floor(scalaj + kHighsTiny); double fj = scalaj - downaj; - if (fj <= f0) { + if (fj <= f0 + 0.001) { double aj = downaj; updateViolationAndNorm(j, aj, viol, sqrnorm); } @@ -760,7 +760,8 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, HighsCDouble aj = downaj; if (fj > f0) { if (strongcg) { - HighsCDouble pj = ceil(k * (fj - f0) * oneoveroneminusf0); + HighsCDouble pj = ceil(k * (fj - f0) * oneoveroneminusf0 - + kHighsTiny); aj += pj / (k + 1); } else { From 2a2fdae0b8f95b81411512b000a7aefb996bcc55 Mon Sep 17 00:00:00 2001 From: Mark Turner Date: Fri, 6 Jun 2025 20:22:54 +0200 Subject: [PATCH 05/19] Add tolerance to actual cut --- highs/mip/HighsCutGeneration.cpp | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index 8efdfaa7185..a1bd7f8c9a6 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/highs/mip/HighsCutGeneration.cpp @@ -760,9 +760,11 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, HighsCDouble aj = downaj; if (fj > f0) { if (strongcg) { - HighsCDouble pj = ceil(k * (fj - f0) * oneoveroneminusf0 - - kHighsTiny); - aj += pj / (k + 1); + if (fj - f0 > feastol) { + HighsCDouble pj = ceil(k * (fj - f0) * oneoveroneminusf0 - + kHighsTiny); + aj += pj / (k + 1); + } } else { aj += (fj - f0) * oneoveroneminusf0; @@ -783,8 +785,15 @@ 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 { + // We round conservatively for scoring => actual strongcg cut can be stronger + assert(bestefficacy - checkefficacy < 0.001 || + fabs(checkefficacy) >= 1e30); + } } #endif From 6d418599a18d0d6afd582545ea67bee71e9e9772 Mon Sep 17 00:00:00 2001 From: Mark Turner Date: Tue, 10 Jun 2025 10:20:47 +0200 Subject: [PATCH 06/19] Update strongcg logic --- highs/mip/HighsCutGeneration.cpp | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index a1bd7f8c9a6..4867f0fecf8 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/highs/mip/HighsCutGeneration.cpp @@ -523,7 +523,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, integerinds.clear(); integerinds.reserve(rowlen); - bool strongcg = true; + bool strongcg = !onlyInitialCMIRScale; double maxabsdelta = 0.0; constexpr double maxCMirScale = 1e6; constexpr double f0min = 0.005; @@ -693,7 +693,10 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, } } - // If inequality has no neg coef on cont vars: Calc strongcg cut + // 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; @@ -702,7 +705,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, double f0 = scalrhs - downrhs; double oneoveroneminusf0 = 1.0 / (1.0 - f0); // TODO: Should stronger safeguards be used, i.e. multiply delta by 2, if f0 is too large? - double k = fast_ceil(1 / f0) - 1; + double k = fast_ceil((1 / f0) - kHighsTiny) - 1; // All coefficients of continuous variables are 0 in strong CG cut double sqrnorm = 0; @@ -712,7 +715,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, double scalaj = vals[j] * scale; double downaj = fast_floor(scalaj + kHighsTiny); double fj = scalaj - downaj; - if (fj <= f0 + 0.001) { + if (fj <= f0 + 0.0001) { double aj = downaj; updateViolationAndNorm(j, aj, viol, sqrnorm); } @@ -760,9 +763,8 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, HighsCDouble aj = downaj; if (fj > f0) { if (strongcg) { - if (fj - f0 > feastol) { - HighsCDouble pj = ceil(k * (fj - f0) * oneoveroneminusf0 - - kHighsTiny); + if (fj - f0 > epsilon) { + HighsCDouble pj = ceil(k * (fj - f0) * oneoveroneminusf0); aj += pj / (k + 1); } } @@ -790,7 +792,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, fabs(checkefficacy) >= 1e30); } else { - // We round conservatively for scoring => actual strongcg cut can be stronger + // Rounded conservatively for scoring => strongcg cut can be stronger assert(bestefficacy - checkefficacy < 0.001 || fabs(checkefficacy) >= 1e30); } From a8b844415b10bcadcdc3df8d8f6ef7617a789808 Mon Sep 17 00:00:00 2001 From: Mark Turner Date: Tue, 10 Jun 2025 12:06:17 +0200 Subject: [PATCH 07/19] Add explicit tolerances not tied to feastol --- highs/mip/HighsCutGeneration.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index 4867f0fecf8..8ff9e74adb8 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/highs/mip/HighsCutGeneration.cpp @@ -704,7 +704,6 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, double downrhs = fast_floor(scalrhs); double f0 = scalrhs - downrhs; double oneoveroneminusf0 = 1.0 / (1.0 - f0); - // TODO: Should stronger safeguards be used, i.e. multiply delta by 2, if f0 is too large? double k = fast_ceil((1 / f0) - kHighsTiny) - 1; // All coefficients of continuous variables are 0 in strong CG cut @@ -715,7 +714,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, double scalaj = vals[j] * scale; double downaj = fast_floor(scalaj + kHighsTiny); double fj = scalaj - downaj; - if (fj <= f0 + 0.0001) { + if (fj <= f0 + 1e-5) { double aj = downaj; updateViolationAndNorm(j, aj, viol, sqrnorm); } @@ -727,7 +726,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, } double efficacy = viol / sqrt(sqrnorm); // Use the strong CG cut instead of the CMIR if efficacy is larger - if (efficacy < bestefficacy + 10 * feastol) { + if (efficacy < bestefficacy + 1e-5) { strongcg = false; } else { From 0ba53e0b924311f374e89b0602e10b459e15b8f1 Mon Sep 17 00:00:00 2001 From: Mark Turner Date: Tue, 10 Jun 2025 12:24:50 +0200 Subject: [PATCH 08/19] Format else statements --- highs/mip/HighsCutGeneration.cpp | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index 8ff9e74adb8..4a93e3128a1 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/highs/mip/HighsCutGeneration.cpp @@ -545,8 +545,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, maxabsdelta = max(maxabsdelta, delta); deltas.push_back(delta); } - } - else if (vals[i] < 0) { + } else if (vals[i] < 0) { // Positive coefficient values later set to 0 so have no contribution updateViolationAndNorm(i, vals[i], continuouscontribution, continuoussqrnorm); @@ -717,8 +716,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, if (fj <= f0 + 1e-5) { double aj = downaj; updateViolationAndNorm(j, aj, viol, sqrnorm); - } - else { + } else { double pj = fast_ceil(k * (fj - f0) * oneoveroneminusf0 - kHighsTiny); double aj = downaj + (pj / (k + 1)); updateViolationAndNorm(j, aj, viol, sqrnorm); @@ -728,8 +726,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, // Use the strong CG cut instead of the CMIR if efficacy is larger if (efficacy < bestefficacy + 1e-5) { strongcg = false; - } - else { + } else { bestefficacy = efficacy; } } @@ -766,8 +763,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, HighsCDouble pj = ceil(k * (fj - f0) * oneoveroneminusf0); aj += pj / (k + 1); } - } - else { + } else { aj += (fj - f0) * oneoveroneminusf0; } } @@ -789,8 +785,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, if (!strongcg) { assert(fabs(checkefficacy - bestefficacy) < 0.001 || fabs(checkefficacy) >= 1e30); - } - else { + } else { // Rounded conservatively for scoring => strongcg cut can be stronger assert(bestefficacy - checkefficacy < 0.001 || fabs(checkefficacy) >= 1e30); From 2d7900a298d84341e90b6310fc393dceab4d574c Mon Sep 17 00:00:00 2001 From: Mark Turner Date: Tue, 10 Jun 2025 15:48:00 +0200 Subject: [PATCH 09/19] Catch case where all coefs are rounded 0 --- highs/mip/HighsCutGeneration.cpp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index 4a93e3128a1..2b6cac04a58 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/highs/mip/HighsCutGeneration.cpp @@ -722,12 +722,16 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, updateViolationAndNorm(j, aj, viol, sqrnorm); } } - double efficacy = viol / sqrt(sqrnorm); - // Use the strong CG cut instead of the CMIR if efficacy is larger - if (efficacy < bestefficacy + 1e-5) { + if (sqrnorm <= kHighsTiny) { strongcg = false; } else { - bestefficacy = efficacy; + double efficacy = viol / sqrt(sqrnorm); + // Use the strong CG cut instead of the CMIR if efficacy is larger + if (efficacy < bestefficacy + 1e-5) { + strongcg = false; + } else { + bestefficacy = efficacy; + } } } From b2ffd809d6f9244614fb4e8898e1b0ab080209ef Mon Sep 17 00:00:00 2001 From: Mark Turner Date: Wed, 11 Jun 2025 10:59:35 +0200 Subject: [PATCH 10/19] Change KHighsTiny to epsilon in approximate cut calc --- highs/mip/HighsCutGeneration.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index 2b6cac04a58..ebea505f20b 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/highs/mip/HighsCutGeneration.cpp @@ -703,7 +703,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, double downrhs = fast_floor(scalrhs); double f0 = scalrhs - downrhs; double oneoveroneminusf0 = 1.0 / (1.0 - f0); - double k = fast_ceil((1 / f0) - kHighsTiny) - 1; + double k = fast_ceil((1 / f0) - epsilon) - 1; // All coefficients of continuous variables are 0 in strong CG cut double sqrnorm = 0; From a0248492453ae913a48d02dfd5d24958e0d42c20 Mon Sep 17 00:00:00 2001 From: Mark Turner Date: Wed, 11 Jun 2025 11:38:10 +0200 Subject: [PATCH 11/19] Added additional numerical safeguard --- highs/mip/HighsCutGeneration.cpp | 52 ++++++++++++++++++-------------- 1 file changed, 29 insertions(+), 23 deletions(-) diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index ebea505f20b..b9f5a5cb475 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/highs/mip/HighsCutGeneration.cpp @@ -703,34 +703,40 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, double downrhs = fast_floor(scalrhs); double f0 = scalrhs - downrhs; double oneoveroneminusf0 = 1.0 / (1.0 - f0); + // Skip numerically troublesome cuts double k = fast_ceil((1 / f0) - epsilon) - 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; - if (fj <= f0 + 1e-5) { - double aj = downaj; - updateViolationAndNorm(j, aj, viol, sqrnorm); - } else { - double pj = fast_ceil(k * (fj - f0) * oneoveroneminusf0 - kHighsTiny); - double aj = downaj + (pj / (k + 1)); - updateViolationAndNorm(j, aj, viol, sqrnorm); - } - } - if (sqrnorm <= kHighsTiny) { + double checkk = fast_ceil((1 / f0) + epsilon) - 1; + if (fabs(k - checkk) >= 0.5) { strongcg = false; } else { - double efficacy = viol / sqrt(sqrnorm); - // Use the strong CG cut instead of the CMIR if efficacy is larger - if (efficacy < bestefficacy + 1e-5) { + // 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; + if (fj <= f0 + 1e-5) { + double aj = downaj; + updateViolationAndNorm(j, aj, viol, sqrnorm); + } else { + double pj = fast_ceil( + k * (fj - f0) * oneoveroneminusf0 - kHighsTiny); + double aj = downaj + (pj / (k + 1)); + updateViolationAndNorm(j, aj, viol, sqrnorm); + } + } + if (sqrnorm <= kHighsTiny) { strongcg = false; } else { - bestefficacy = efficacy; + double efficacy = viol / sqrt(sqrnorm); + // Use the strong CG cut instead of the CMIR if efficacy is larger + if (efficacy < bestefficacy + 1e-5) { + strongcg = false; + } else { + bestefficacy = efficacy; + } } } } From ee50b06fd1ff0d9a9d8a98876f24a07b7f18fa60 Mon Sep 17 00:00:00 2001 From: Mark Turner Date: Wed, 11 Jun 2025 12:44:46 +0200 Subject: [PATCH 12/19] Relax approx tolerances. Enforce tolerance for ceil --- highs/mip/HighsCutGeneration.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index b9f5a5cb475..8866e3b19a3 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/highs/mip/HighsCutGeneration.cpp @@ -717,12 +717,12 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, double scalaj = vals[j] * scale; double downaj = fast_floor(scalaj + kHighsTiny); double fj = scalaj - downaj; - if (fj <= f0 + 1e-5) { + if (fj <= f0 + 1e-6) { double aj = downaj; updateViolationAndNorm(j, aj, viol, sqrnorm); } else { - double pj = fast_ceil( - k * (fj - f0) * oneoveroneminusf0 - kHighsTiny); + double pj = + fast_ceil(k * (fj - f0) * oneoveroneminusf0 - (10 * epsilon)); double aj = downaj + (pj / (k + 1)); updateViolationAndNorm(j, aj, viol, sqrnorm); } @@ -732,7 +732,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, } else { double efficacy = viol / sqrt(sqrnorm); // Use the strong CG cut instead of the CMIR if efficacy is larger - if (efficacy < bestefficacy + 1e-5) { + if (efficacy < bestefficacy + 1e-6) { strongcg = false; } else { bestefficacy = efficacy; @@ -770,7 +770,8 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, if (fj > f0) { if (strongcg) { if (fj - f0 > epsilon) { - HighsCDouble pj = ceil(k * (fj - f0) * oneoveroneminusf0); + HighsCDouble pj = + ceil(k * (fj - f0) * oneoveroneminusf0 - epsilon); aj += pj / (k + 1); } } else { From 9903204b6d39cffef2991ad383310ba615947f3f Mon Sep 17 00:00:00 2001 From: Mark Turner Date: Wed, 11 Jun 2025 14:25:52 +0200 Subject: [PATCH 13/19] Make clang format happy --- highs/mip/HighsCutGeneration.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index 8866e3b19a3..080f303327d 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/highs/mip/HighsCutGeneration.cpp @@ -706,7 +706,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, // Skip numerically troublesome cuts double k = fast_ceil((1 / f0) - epsilon) - 1; double checkk = fast_ceil((1 / f0) + epsilon) - 1; - if (fabs(k - checkk) >= 0.5) { + if (checkk - k > 0.5) { strongcg = false; } else { // All coefficients of continuous variables are 0 in strong CG cut @@ -722,7 +722,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, updateViolationAndNorm(j, aj, viol, sqrnorm); } else { double pj = - fast_ceil(k * (fj - f0) * oneoveroneminusf0 - (10 * epsilon)); + fast_ceil(k * (fj - f0) * oneoveroneminusf0 - (10 * epsilon)); double aj = downaj + (pj / (k + 1)); updateViolationAndNorm(j, aj, viol, sqrnorm); } @@ -770,8 +770,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, if (fj > f0) { if (strongcg) { if (fj - f0 > epsilon) { - HighsCDouble pj = - ceil(k * (fj - f0) * oneoveroneminusf0 - epsilon); + HighsCDouble pj = ceil(k * (fj - f0) * oneoveroneminusf0 - epsilon); aj += pj / (k + 1); } } else { From 8c54b3e8f0cde875229c7122b411072b2346a0ec Mon Sep 17 00:00:00 2001 From: Mark Turner Date: Tue, 17 Jun 2025 12:57:38 +0200 Subject: [PATCH 14/19] Calc strongcg for conflict + modk --- highs/mip/HighsCutGeneration.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index 080f303327d..33465cb9e67 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/highs/mip/HighsCutGeneration.cpp @@ -523,7 +523,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, integerinds.clear(); integerinds.reserve(rowlen); - bool strongcg = !onlyInitialCMIRScale; + bool strongcg = true; double maxabsdelta = 0.0; constexpr double maxCMirScale = 1e6; constexpr double f0min = 0.005; From bf28f6a48a1f8f2b0df224ddc9251d5eda5e963b Mon Sep 17 00:00:00 2001 From: Mark Turner Date: Mon, 4 Aug 2025 11:52:04 +0200 Subject: [PATCH 15/19] Update tolerances --- highs/mip/HighsCutGeneration.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index 33465cb9e67..2b31a323a61 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/highs/mip/HighsCutGeneration.cpp @@ -704,9 +704,9 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, double f0 = scalrhs - downrhs; double oneoveroneminusf0 = 1.0 / (1.0 - f0); // Skip numerically troublesome cuts - double k = fast_ceil((1 / f0) - epsilon) - 1; - double checkk = fast_ceil((1 / f0) + epsilon) - 1; - if (checkk - k > 0.5) { + double oneoverf0 = 1 / f0; + double k = fast_ceil(oneoverf0) - 1; + if (oneoverf0 - k < 1e-3 || oneoverf0 - k > 1 - 1e-3) { strongcg = false; } else { // All coefficients of continuous variables are 0 in strong CG cut @@ -717,12 +717,12 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, double scalaj = vals[j] * scale; double downaj = fast_floor(scalaj + kHighsTiny); double fj = scalaj - downaj; - if (fj <= f0 + 1e-6) { + if (fj <= f0 + feastol) { double aj = downaj; updateViolationAndNorm(j, aj, viol, sqrnorm); } else { double pj = - fast_ceil(k * (fj - f0) * oneoveroneminusf0 - (10 * epsilon)); + fast_ceil(k * (fj - f0) * oneoveroneminusf0 - feastol); double aj = downaj + (pj / (k + 1)); updateViolationAndNorm(j, aj, viol, sqrnorm); } @@ -732,7 +732,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, } else { double efficacy = viol / sqrt(sqrnorm); // Use the strong CG cut instead of the CMIR if efficacy is larger - if (efficacy < bestefficacy + 1e-6) { + if (efficacy < bestefficacy + epsilon) { strongcg = false; } else { bestefficacy = efficacy; From 775b04acdd5abe291dad0640d7607b9b095e0c88 Mon Sep 17 00:00:00 2001 From: Mark Turner Date: Tue, 14 Oct 2025 14:10:31 +0200 Subject: [PATCH 16/19] Make linter happy --- highs/mip/HighsCutGeneration.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index 2b31a323a61..d26777cfb42 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/highs/mip/HighsCutGeneration.cpp @@ -721,8 +721,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, double aj = downaj; updateViolationAndNorm(j, aj, viol, sqrnorm); } else { - double pj = - fast_ceil(k * (fj - f0) * oneoveroneminusf0 - feastol); + double pj = fast_ceil(k * (fj - f0) * oneoveroneminusf0 - feastol); double aj = downaj + (pj / (k + 1)); updateViolationAndNorm(j, aj, viol, sqrnorm); } From b89e5ce59dffe85d3c206f443321f215b149ee21 Mon Sep 17 00:00:00 2001 From: Mark Turner Date: Tue, 25 Nov 2025 12:05:57 +0100 Subject: [PATCH 17/19] Make strongCG more conersvative --- highs/mip/HighsCutGeneration.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index d26777cfb42..36e8c008fce 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/highs/mip/HighsCutGeneration.cpp @@ -523,7 +523,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, integerinds.clear(); integerinds.reserve(rowlen); - bool strongcg = true; + bool strongcg = !onlyInitialCMIRScale; double maxabsdelta = 0.0; constexpr double maxCMirScale = 1e6; constexpr double f0min = 0.005; @@ -721,7 +721,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, double aj = downaj; updateViolationAndNorm(j, aj, viol, sqrnorm); } else { - double pj = fast_ceil(k * (fj - f0) * oneoveroneminusf0 - feastol); + double pj = fast_ceil(k * (fj - f0) * oneoveroneminusf0 - 1e-4); double aj = downaj + (pj / (k + 1)); updateViolationAndNorm(j, aj, viol, sqrnorm); } @@ -769,7 +769,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, if (fj > f0) { if (strongcg) { if (fj - f0 > epsilon) { - HighsCDouble pj = ceil(k * (fj - f0) * oneoveroneminusf0 - epsilon); + HighsCDouble pj = ceil(k * (fj - f0) * oneoveroneminusf0 - feastol); aj += pj / (k + 1); } } else { From 13a5edfbad074d27c6525a9260db713148da5b97 Mon Sep 17 00:00:00 2001 From: Mark Turner Date: Thu, 18 Dec 2025 14:10:16 +0100 Subject: [PATCH 18/19] Add an additional complementation attempt for strongcg --- highs/mip/HighsCutGeneration.cpp | 74 +++++++++++++++++++++++++++----- 1 file changed, 64 insertions(+), 10 deletions(-) diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index 36e8c008fce..4d1c9cd190d 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/highs/mip/HighsCutGeneration.cpp @@ -549,6 +549,7 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, // 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; } } @@ -704,11 +705,10 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, double f0 = scalrhs - downrhs; double oneoveroneminusf0 = 1.0 / (1.0 - f0); // Skip numerically troublesome cuts - double oneoverf0 = 1 / f0; - double k = fast_ceil(oneoverf0) - 1; - if (oneoverf0 - k < 1e-3 || oneoverf0 - k > 1 - 1e-3) { + 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; @@ -717,20 +717,18 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, double scalaj = vals[j] * scale; double downaj = fast_floor(scalaj + kHighsTiny); double fj = scalaj - downaj; - if (fj <= f0 + feastol) { - double aj = downaj; - updateViolationAndNorm(j, aj, viol, sqrnorm); - } else { + double aj = downaj; + if (fj >= f0 + 10 * feastol) { double pj = fast_ceil(k * (fj - f0) * oneoveroneminusf0 - 1e-4); - double aj = downaj + (pj / (k + 1)); - updateViolationAndNorm(j, aj, viol, sqrnorm); + aj += pj / (k + 1); } + updateViolationAndNorm(j, aj, viol, sqrnorm); } if (sqrnorm <= kHighsTiny) { strongcg = false; } else { double efficacy = viol / sqrt(sqrnorm); - // Use the strong CG cut instead of the CMIR if efficacy is larger + // Use the strongCG cut instead of the CMIR if efficacy is larger if (efficacy < bestefficacy + epsilon) { strongcg = false; } else { @@ -740,6 +738,62 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, } } + 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)); From 70837feb70b5c181261e08a4e1fdaf0220751356 Mon Sep 17 00:00:00 2001 From: Mark Turner Date: Wed, 4 Mar 2026 11:00:37 +0100 Subject: [PATCH 19/19] Remove attempt at complementing strongcg --- highs/mip/HighsCutGeneration.cpp | 56 -------------------------------- 1 file changed, 56 deletions(-) diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index 4d1c9cd190d..a4a906c6885 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/highs/mip/HighsCutGeneration.cpp @@ -738,62 +738,6 @@ bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, } } - 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));