diff --git a/check/TestCallbacks.cpp b/check/TestCallbacks.cpp index e8b2e919bd..2706c35a92 100644 --- a/check/TestCallbacks.cpp +++ b/check/TestCallbacks.cpp @@ -426,7 +426,7 @@ TEST_CASE("highs-callback-mip-interrupt", "[highs_callback]") { HighsStatus status = highs.run(); REQUIRE(status == HighsStatus::kWarning); REQUIRE(highs.getModelStatus() == HighsModelStatus::kInterrupt); - REQUIRE(highs.getInfo().objective_function_value > egout_optimal_objective); + REQUIRE(highs.getInfo().objective_function_value >= egout_optimal_objective); highs.resetGlobalScheduler(true); } diff --git a/check/TestMipSolver.cpp b/check/TestMipSolver.cpp index b09b601d58..e84e3f52f0 100644 --- a/check/TestMipSolver.cpp +++ b/check/TestMipSolver.cpp @@ -609,7 +609,7 @@ TEST_CASE("MIP-get-saved-solutions", "[highs_test_mip_solver]") { TEST_CASE("MIP-objective-target", "[highs_test_mip_solver]") { const double egout_optimal_objective = 568.1007; - const double egout_objective_target = 610; + const double egout_objective_target = 680; std::string filename = std::string(HIGHS_DIR) + "/check/instances/egout.mps"; Highs highs; highs.setOptionValue("output_flag", dev_run); diff --git a/highs/lp_data/HighsOptions.h b/highs/lp_data/HighsOptions.h index 7b0fe3bdb1..07f72c07cd 100644 --- a/highs/lp_data/HighsOptions.h +++ b/highs/lp_data/HighsOptions.h @@ -480,6 +480,7 @@ struct HighsOptionsStruct { bool mip_heuristic_run_root_reduced_cost; bool mip_heuristic_run_zi_round; bool mip_heuristic_run_shifting; + bool mip_cut_flow_cover; double mip_min_logging_interval; std::string mip_lp_solver; std::string mip_ipm_solver; @@ -638,6 +639,7 @@ struct HighsOptionsStruct { mip_heuristic_run_root_reduced_cost(false), mip_heuristic_run_zi_round(false), mip_heuristic_run_shifting(false), + mip_cut_flow_cover(false), mip_min_logging_interval(0.0), mip_lp_solver(""), mip_ipm_solver(""), @@ -1217,6 +1219,11 @@ class HighsOptions : public HighsOptionsStruct { &mip_heuristic_run_shifting, false); records.push_back(record_bool); + record_bool = + new OptionRecordBool("mip_cut_flow_cover", "Enable flow cover cuts", + advanced, &mip_cut_flow_cover, true); + records.push_back(record_bool); + record_bool = new OptionRecordBool( "mip_allow_cut_separation_at_nodes", "Whether cut separation at nodes is permitted", advanced, diff --git a/highs/mip/HighsCutGeneration.cpp b/highs/mip/HighsCutGeneration.cpp index f518b7d800..6430e49a92 100644 --- a/highs/mip/HighsCutGeneration.cpp +++ b/highs/mip/HighsCutGeneration.cpp @@ -500,6 +500,238 @@ bool HighsCutGeneration::separateLiftedMixedIntegerCover() { return true; } +// Calculates the lifted simple generalized flow cover cut out of +// Gu, Z., Nemhauser, G. L., & Savelsbergh, M. W. (1999). +// Lifted flow cover inequalities for mixed 0-1 integer programs. +// Mathematical Programming, 85(3), 439-467. +bool HighsCutGeneration::separateLiftedFlowCover() { + const double vubEpsilon = 1e-8; + // Compute the lifting data (ld) first + struct LiftingData { + std::vector m; + std::vector M; + HighsInt r = 0; + HighsInt t = 0; + double mp = kHighsInf; + double ml = 0; + HighsCDouble d1 = 0; + }; + LiftingData ld; + ld.m.resize(snfr.numNnzs); + HighsCDouble sumN2mC2LE = 0.0; + HighsCDouble sumC1LE = 0.0; + HighsCDouble sumC2 = 0.0; + + for (HighsInt i = 0; i != snfr.numNnzs; ++i) { + // col is in N- \ C- + if (snfr.flowCoverStatus[i] == -1 && snfr.coef[i] == -1) { + assert(snfr.vubCoef[i] >= 0); + if (snfr.vubCoef[i] > snfr.lambda + vubEpsilon) { + ld.m[ld.r] = snfr.vubCoef[i]; + ++ld.r; + } else { + sumN2mC2LE += snfr.vubCoef[i]; + } + } else if (snfr.flowCoverStatus[i] == 1 && snfr.coef[i] == -1) { + // col is in C- + assert(snfr.vubCoef[i] > 0); + sumC2 += snfr.vubCoef[i]; + } else if (snfr.flowCoverStatus[i] == 1 && snfr.coef[i] == 1) { + // col is in C+ + assert(snfr.vubCoef[i] > 0); + if (snfr.vubCoef[i] > snfr.lambda + vubEpsilon) { + ld.m[ld.r] = snfr.vubCoef[i]; + ++ld.r; + ld.mp = std::min(ld.mp, snfr.vubCoef[i]); + } else { + sumC1LE += snfr.vubCoef[i]; + } + } else { + // col is in N+ \ C+ + assert(snfr.flowCoverStatus[i] == -1 && snfr.coef[i] == 1); + } + } + + if (ld.mp == kHighsInf) return false; + + ld.M.resize(ld.r + 1); + ld.ml = std::min(snfr.lambda, static_cast(sumC1LE + sumN2mC2LE)); + ld.d1 = sumC2 + snfr.rhs; + + pdqsort_branchless(ld.m.begin(), ld.m.begin() + ld.r, + [&](const double a, const double b) { return a > b; }); + + for (HighsInt i = 0; i != ld.r + 1; ++i) { + if (i == 0) { + ld.M[i] = 0; + } else { + ld.M[i] = ld.M[i - 1] + ld.m[i - 1]; + } + } + + for (HighsInt i = ld.r - 1; i != -1; --i) { + if (ld.m[i] >= ld.mp) { + ld.t = i; + break; + } + } + assert(ld.m[ld.t] == ld.mp); + ld.t++; + + auto getAlphaBeta = [&](double vubcoef) { + HighsInt alpha; + HighsCDouble beta{}; + double vubcoefpluslambda = vubcoef + snfr.lambda; + + HighsInt i = 0; + while (i < ld.r && vubcoefpluslambda >= ld.M[i + 1] + vubEpsilon) { + ++i; + } + + if (vubcoef <= ld.M[i] - vubEpsilon) { + assert(ld.M[i] < vubcoefpluslambda); + alpha = 1; + beta = -i * HighsCDouble(snfr.lambda) + ld.M[i]; + } else { + alpha = 0; + beta = 0; + } + return std::make_pair(alpha, beta); + }; + + auto evaluateLiftingFunction = [&](double vubcoef) { + HighsInt i = 0; + while (i < ld.r && vubcoef + snfr.lambda >= ld.M[i + 1] + vubEpsilon) { + ++i; + } + if (i < ld.t) { + HighsCDouble liftedcoef = i * HighsCDouble(snfr.lambda); + if (ld.M[i] < vubcoef + vubEpsilon) { + return static_cast(liftedcoef); + } + assert(i > 0 && ld.M[i] < vubcoef + snfr.lambda - vubEpsilon && + vubcoef <= ld.M[i]); + liftedcoef += vubcoef; + liftedcoef -= ld.M[i]; + return static_cast(liftedcoef); + } else if (i < ld.r) { + HighsCDouble tmp = HighsCDouble(ld.m[i]) - ld.mp - ld.ml + snfr.lambda; + if (tmp < 0) tmp = 0; + tmp += ld.M[i] + ld.ml; + if (tmp < vubcoef + snfr.lambda - vubEpsilon) { + return static_cast(i * HighsCDouble(snfr.lambda)); + } + assert(ld.M[i] <= vubcoef + snfr.lambda + feastol && + snfr.lambda + vubcoef <= + ld.M[i] + ld.ml + feastol + + std::max(0.0, static_cast( + ld.m[i] - + (HighsCDouble(ld.mp) - snfr.lambda) - + ld.ml))); + return static_cast(i * HighsCDouble(snfr.lambda) + vubcoef - + ld.M[i]); + } + assert(i == ld.r && ld.M[i] <= vubcoef + snfr.lambda + vubEpsilon); + return static_cast(ld.r * HighsCDouble(snfr.lambda) + vubcoef - + ld.M[ld.r]); + }; + + // Calculate the lifted simple generalized flow cover cut + // L- = {i \in N- \ C- : m_i > lambda} + // L-- = N- \ (L- union C-) + HighsCDouble tmpRhs = ld.d1; + rowlen = 0; + + auto addCutNonZero = [&](const double& coef, const HighsInt& index, + const bool complement) -> void { + vals[rowlen] = coef; + inds[rowlen] = index; + if (complement) { + tmpRhs -= vals[rowlen]; + vals[rowlen] = -vals[rowlen]; + } + rowlen++; + }; + + for (HighsInt i = 0; i != snfr.numNnzs; ++i) { + // col is in N- \ C- + if (snfr.flowCoverStatus[i] == -1 && snfr.coef[i] == -1) { + if (snfr.vubCoef[i] > snfr.lambda + vubEpsilon) { + if (snfr.origBinCols[i] != -1) { + // col is in L- + addCutNonZero(-snfr.lambda, snfr.origBinCols[i], + snfr.complementation[i]); + } else { + tmpRhs += snfr.lambda; + } + } else { + // col is in L-- + if (snfr.origContCols[i] != -1 && snfr.aggrContCoef[i] != 0) { + addCutNonZero(-snfr.aggrContCoef[i], snfr.origContCols[i], false); + } + if (snfr.origBinCols[i] != -1 && snfr.aggrBinCoef[i] != 0) { + addCutNonZero(-snfr.aggrBinCoef[i], snfr.origBinCols[i], + snfr.complementation[i]); + } + tmpRhs += snfr.aggrConstant[i]; + } + } else if (snfr.flowCoverStatus[i] == 1 && snfr.coef[i] == -1) { + // col is in C- + if (snfr.origBinCols[i] != -1) { + double liftedbincoef = evaluateLiftingFunction(snfr.vubCoef[i]); + if (liftedbincoef != 0) { + addCutNonZero(-liftedbincoef, snfr.origBinCols[i], + snfr.complementation[i]); + tmpRhs -= liftedbincoef; + } + } + } else if (snfr.flowCoverStatus[i] == -1 && snfr.coef[i] == 1) { + // col is in N+ \ C+ + std::pair alphabeta = + getAlphaBeta(snfr.vubCoef[i]); + assert(alphabeta.first == 0 || alphabeta.first == 1); + if (alphabeta.first == 1) { + assert(alphabeta.second > 0); + if (snfr.origContCols[i] != -1 && snfr.aggrContCoef[i] != 0) { + addCutNonZero(snfr.aggrContCoef[i], snfr.origContCols[i], false); + } + HighsCDouble binvarcoef = snfr.aggrBinCoef[i] - alphabeta.second; + if (snfr.origBinCols[i] != -1) { + if (binvarcoef != 0) { + addCutNonZero(static_cast(binvarcoef), snfr.origBinCols[i], + snfr.complementation[i]); + } + } else { + tmpRhs -= binvarcoef; + } + tmpRhs -= snfr.aggrConstant[i]; + } + } else { + // col is in C+ + assert(snfr.flowCoverStatus[i] == 1 && snfr.coef[i] == 1); + HighsCDouble bincoef = snfr.aggrBinCoef[i]; + HighsCDouble constant = snfr.aggrConstant[i]; + if (snfr.origBinCols[i] != -1 && + snfr.vubCoef[i] >= snfr.lambda + vubEpsilon) { + // col is in C++ + constant += HighsCDouble(snfr.vubCoef[i]) - snfr.lambda; + bincoef -= HighsCDouble(snfr.vubCoef[i]) - snfr.lambda; + } + if (snfr.origBinCols[i] != -1 && bincoef != 0) { + addCutNonZero(static_cast(bincoef), snfr.origBinCols[i], + snfr.complementation[i]); + } + if (snfr.origContCols[i] != -1 && snfr.aggrContCoef[i] != 0) { + addCutNonZero(snfr.aggrContCoef[i], snfr.origContCols[i], false); + } + tmpRhs -= constant; + } + } + + rhs = static_cast(tmpRhs); + return true; +} + static double fast_floor(double x) { return (int64_t)x - (x < (int64_t)x); } bool HighsCutGeneration::cmirCutGenerationHeuristic(double minEfficacy, @@ -1066,7 +1298,8 @@ static void checkNumerics(const double* vals, HighsInt len, double rhs) { bool HighsCutGeneration::generateCut(HighsTransformedLp& transLp, std::vector& inds_, std::vector& vals_, double& rhs_, - bool onlyInitialCMIRScale) { + bool onlyInitialCMIRScale, + bool genFlowCover) { #if 0 if (vals_.size() > 1) { std::vector indsCheck_ = inds_; @@ -1110,6 +1343,32 @@ bool HighsCutGeneration::generateCut(HighsTransformedLp& transLp, } #endif + // Copy data to later generate lifted simple generalized flow cover cut + std::vector flowCoverVals; + std::vector flowCoverInds; + double flowCoverRhs = rhs_; + double flowCoverEfficacy = 0; + if (genFlowCover && !lpRelaxation.getMipSolver().submip && + !lpRelaxation.getMipSolver().mipdata_->continuous_cols.empty() && + lpRelaxation.getMipSolver().options_mip_->mip_cut_flow_cover) { + bool hasContinuousBeforePreprocess = false; + for (size_t i = 0; i != inds_.size(); ++i) { + if (lpRelaxation.isColIntegral(inds_[i]) && + std::abs(vals_[i]) > 10 * feastol) { + hasContinuousBeforePreprocess = true; + break; + } + } + if (!hasContinuousBeforePreprocess) { + genFlowCover = false; + } else { + flowCoverVals = vals_; + flowCoverInds = inds_; + } + } else { + genFlowCover = false; + } + bool intsPositive = true; if (!transLp.transform(vals_, upper, solval, inds_, rhs_, intsPositive)) return false; @@ -1144,28 +1403,72 @@ bool HighsCutGeneration::generateCut(HighsTransformedLp& transLp, } // try to generate a cut - if (!tryGenerateCut(inds_, vals_, hasUnboundedInts, hasGeneralInts, - hasContinuous, 10 * feastol, onlyInitialCMIRScale)) - return false; + bool cmirSuccess = + tryGenerateCut(inds_, vals_, hasUnboundedInts, hasGeneralInts, + hasContinuous, 10 * feastol, onlyInitialCMIRScale); - // remove the complementation if exists - removeComplementation(); + if (cmirSuccess) { + // remove the complementation if exists + removeComplementation(); - // remove zeros in place - for (HighsInt i = rowlen - 1; i >= 0; --i) { - if (vals[i] == 0.0) { - --rowlen; - inds[i] = inds[rowlen]; - vals[i] = vals[rowlen]; + // remove zeros in place + for (HighsInt i = rowlen - 1; i >= 0; --i) { + if (vals[i] == 0.0) { + --rowlen; + inds[i] = inds[rowlen]; + vals[i] = vals[rowlen]; + } } + + // transform the cut back into the original space, i.e. remove the bound + // substitution and replace implicit slack variables + rhs_ = (double)rhs; + vals_.resize(rowlen); + inds_.resize(rowlen); + cmirSuccess = transLp.untransform(vals_, inds_, rhs_); } - // transform the cut back into the original space, i.e. remove the bound - // substitution and replace implicit slack variables - rhs_ = (double)rhs; - vals_.resize(rowlen); - inds_.resize(rowlen); - if (!transLp.untransform(vals_, inds_, rhs_)) return false; + const auto& sol = lpRelaxation.getSolution().col_value; + + // Try to generate a lifted simple generalized flow cover cut + bool flowCoverSuccess = + genFlowCover + ? tryGenerateFlowCoverCut(transLp, flowCoverInds, flowCoverVals, + flowCoverRhs, flowCoverEfficacy) + : false; + if (flowCoverSuccess && cmirSuccess) { + HighsInt rowlen_ = inds_.size(); + double viol = -rhs_; + double sqrnorm = 0; + double efficacy = 0; + for (HighsInt i = 0; i != rowlen_; ++i) { + HighsInt col = inds_[i]; + viol += vals_[i] * sol[col]; + if (vals_[i] >= 0 && + sol[col] <= + lpRelaxation.getMipSolver().mipdata_->domain.col_lower_[col] + + lpRelaxation.getMipSolver().mipdata_->feastol) + continue; + if (vals_[i] < 0 && + sol[col] >= + lpRelaxation.getMipSolver().mipdata_->domain.col_upper_[col] - + lpRelaxation.getMipSolver().mipdata_->feastol) + continue; + sqrnorm += vals_[i] * vals_[i]; + } + if (sqrnorm > 0) efficacy = viol / sqrt(sqrnorm); + if (flowCoverEfficacy < efficacy + 10 * feastol) { + flowCoverSuccess = false; + } + } + if (!flowCoverSuccess && !cmirSuccess) return false; + if (flowCoverSuccess) { + rhs_ = flowCoverRhs; + std::swap(vals_, flowCoverVals); + std::swap(inds_, flowCoverInds); + integralSupport = false; + integralCoefficients = false; + } rowlen = inds_.size(); inds = inds_.data(); @@ -1186,7 +1489,6 @@ bool HighsCutGeneration::generateCut(HighsTransformedLp& transLp, // finally determine the violation of the cut in the original space HighsCDouble violation = -rhs_; - const auto& sol = lpRelaxation.getSolution().col_value; for (HighsInt i = 0; i != rowlen; ++i) violation += sol[inds[i]] * vals_[i]; if (violation <= 10 * feastol) return false; @@ -1299,6 +1601,257 @@ bool HighsCutGeneration::generateConflict(HighsDomain& localdomain, return cutindex != -1; } +void HighsCutGeneration::initSNFRelaxation() { + if (static_cast(snfr.coef.size()) < rowlen) { + snfr.origBinCols.resize(rowlen); + snfr.origContCols.resize(rowlen); + snfr.binSolval.resize(rowlen); + snfr.coef.resize(rowlen); + snfr.vubCoef.resize(rowlen); + snfr.aggrConstant.resize(rowlen); + snfr.aggrBinCoef.resize(rowlen); + snfr.aggrContCoef.resize(rowlen); + snfr.complementation.resize(rowlen); + snfr.flowCoverStatus.resize(rowlen); + } + snfr.rhs = 0; + snfr.lambda = 0; + snfr.numNnzs = 0; +} + +bool HighsCutGeneration::preprocessSNFRelaxation() { + // preprocess the inequality before generating a single node flow relaxation. + // 1. Determine the maximal activity to check for trivial redundancy + // 2. Check for presence of unbounded variables (the SNFR construction + // requires finite bounds on all non-zeros) + // 3. Remove coefficients that are below the feasibility tolerance to avoid + // numerical troubles, use bound constraints to cancel them and + // reject base inequalities where that is not possible due to unbounded + // variables + // 4. Don't consider any inequality with too many non-zeros + // 5. Don't consider any inequality with too few continuous cols + + HighsInt maxLen = 100 + 0.05 * (lpRelaxation.numCols()); + if (rowlen > maxLen) return false; + + HighsInt numZeros = 0; + HighsInt numContCols = 0; + double maxact = -feastol; + double maxAbsVal = 0; + HighsInt slackOffset = lpRelaxation.getMipSolver().numCol(); + const HighsDomain& domain = lpRelaxation.getMipSolver().mipdata_->domain; + + auto getLb = [&](HighsInt col) { + return (col < slackOffset ? domain.col_lower_[col] + : lpRelaxation.slackLower(col - slackOffset)); + }; + + auto getUb = [&](HighsInt col) { + return (col < slackOffset ? domain.col_upper_[col] + : lpRelaxation.slackUpper(col - slackOffset)); + }; + + for (HighsInt i = 0; i < rowlen; ++i) { + HighsInt col = inds[i]; + double lb = getLb(col); + double ub = getUb(col); + maxAbsVal = std::max(std::abs(vals[i]), maxAbsVal); + if (lb == -kHighsInf || ub == kHighsInf) { + return false; + } + if (col < slackOffset && !lpRelaxation.isColIntegral(col)) numContCols++; + } + + if (numContCols == 0) return false; + if (maxAbsVal > 1 + feastol || maxAbsVal < 0.5 - feastol) scale(maxAbsVal); + + for (HighsInt i = 0; i != rowlen; ++i) { + HighsInt col = inds[i]; + double lb = getLb(col); + double ub = getUb(col); + + // relax variables with small contributions if possible + if (std::abs(vals[i]) * (ub - lb) <= 10 * feastol) { + if (vals[i] < 0) { + if (ub == kHighsInf) return false; + rhs -= vals[i] * ub; + ++numZeros; + vals[i] = 0.0; + } else if (vals[i] > 0) { + if (lb == -kHighsInf) return false; + rhs -= vals[i] * lb; + ++numZeros; + vals[i] = 0.0; + } + } + + if (vals[i] > 0) { + maxact += vals[i] * ub; + } else if (vals[i] < 0) { + maxact += vals[i] * lb; + } + } + + if (numZeros != 0) { + // remove zeros in place + for (HighsInt i = rowlen - 1; i >= 0; --i) { + if (vals[i] == 0.0) { + --rowlen; + inds[i] = inds[rowlen]; + vals[i] = vals[rowlen]; + if (--numZeros == 0) break; + } + } + } + + return maxact > rhs; +} + +bool HighsCutGeneration::computeFlowCover() { + // Compute the flow cover, i.e., get sets C+ subset N+ and C- subset N- + // with sum_{j in C+} u_j - sum_{j in C-} u_j = b + lambda, lambda > 0 + std::vector items(snfr.numNnzs, -1); + HighsInt nNonFlowCover = 0; + HighsInt nFlowCover = 0; + HighsInt nitems = 0; + HighsCDouble n1itemsWeight = 0; + HighsCDouble flowCoverWeight = 0; + for (HighsInt i = 0; i < snfr.numNnzs; ++i) { + assert(snfr.coef[i] == 1 || snfr.coef[i] == -1); + assert(snfr.binSolval[i] >= -feastol && snfr.binSolval[i] <= 1 + feastol); + // if u_i = 0 put i into N+ \ C+ or N- \ C-, i.e., not the cover + if (abs(snfr.vubCoef[i]) < feastol) { + snfr.flowCoverStatus[i] = -1; + nNonFlowCover++; + } else if (fractionality(snfr.binSolval[i]) > feastol) { + // x_i is fractional -> becomes an item in knapsack (decides if in cover) + items[nitems] = i; + nitems++; + if (snfr.coef[i] == 1) { + n1itemsWeight += snfr.vubCoef[i]; + } + } else if (snfr.coef[i] == 1 && snfr.binSolval[i] < 0.5) { + // i is in N+ and x_i = 0 -> don't put in cover + snfr.flowCoverStatus[i] = -1; + nNonFlowCover++; + } else if (snfr.coef[i] == 1 && snfr.binSolval[i] > 0.5) { + // i is in N+ and x_i = 1 -> put in cover + snfr.flowCoverStatus[i] = 1; + nFlowCover++; + flowCoverWeight += snfr.vubCoef[i]; + } else if (snfr.coef[i] == -1 && snfr.binSolval[i] > 0.5) { + // i is in N- and x_i = 1 -> put in cover + snfr.flowCoverStatus[i] = 1; + nFlowCover++; + flowCoverWeight -= snfr.vubCoef[i]; + } else { + // i is in N- and x_i = 0 -> don't put in cover + assert(snfr.coef[i] == -1 && snfr.binSolval[i] < 0.5); + snfr.flowCoverStatus[i] = -1; + nNonFlowCover++; + } + } + assert(nNonFlowCover + nFlowCover + nitems == snfr.numNnzs); + + double capacity = static_cast(-snfr.rhs + flowCoverWeight + + n1itemsWeight - feastol); + + // There is no flow cover if capacity is less than zero after fixing + if (capacity < feastol) return false; + + // Solve a knapsack greedily to assign items to C+, C-, N+\C+, N-\C- + // z_j is a binary variable deciding whether the item goes into the knapsack + // max sum_{j in N+} ( x_j - 1 ) z_j + sum_{j in N-} x_j + // sum_{j in N+} u'_j z_j + sum_{j in N-} u'_j z_j < -b + sum_{j in N+} u'_j + // z_j in {0,1} for all j in N+ & N- + + double knapsackWeight = 0; + std::vector weights(nitems); + std::vector profitweightratios(nitems); + for (HighsInt i = 0; i < nitems; ++i) { + weights[i] = snfr.vubCoef[items[i]]; + if (snfr.coef[items[i]] == 1) { + profitweightratios[i] = (1 - snfr.binSolval[items[i]]) / weights[i]; + } else { + profitweightratios[i] = snfr.binSolval[items[i]] / weights[i]; + } + } + std::vector perm(nitems); + std::iota(perm.begin(), perm.end(), 0); + pdqsort_branchless(perm.begin(), perm.end(), + [&](const HighsInt a, const HighsInt b) { + return profitweightratios[a] > profitweightratios[b]; + }); + // Greedily add items to knapsack + for (HighsInt i = 0; i < nitems; ++i) { + const HighsInt j = perm[i]; + const HighsInt k = items[j]; + if (knapsackWeight + weights[j] < capacity) { + knapsackWeight += weights[j]; + if (snfr.coef[k] == 1) { + // j in N+ with z_j = 1 => j in N+ \ C+ + snfr.flowCoverStatus[k] = -1; + nNonFlowCover++; + } else { + // j in N- with z_j = 1 => j in C- + snfr.flowCoverStatus[k] = 1; + nFlowCover++; + flowCoverWeight -= snfr.vubCoef[k]; + } + } else { + if (snfr.coef[k] == 1) { + // j in N+ with z_j = 0 => j in C+ + snfr.flowCoverStatus[k] = 1; + nFlowCover++; + flowCoverWeight += snfr.vubCoef[k]; + } else { + // j in N- with z_j = 0 => j in N- \ C- + snfr.flowCoverStatus[k] = -1; + nNonFlowCover++; + } + } + } + + assert(nFlowCover + nNonFlowCover == snfr.numNnzs); + + snfr.lambda = static_cast(flowCoverWeight - snfr.rhs); + if (snfr.lambda < feastol) return false; + return true; +} + +bool HighsCutGeneration::tryGenerateFlowCoverCut(HighsTransformedLp& transLp, + std::vector& inds_, + std::vector& vals_, + double& rhs_, + double& efficacy) { + rowlen = static_cast(inds_.size()); + this->inds = inds_.data(); + this->vals = vals_.data(); + this->rhs = rhs_; + if (!preprocessSNFRelaxation()) return false; + initSNFRelaxation(); + vals_.resize(rowlen); + inds_.resize(rowlen); + rhs_ = static_cast(this->rhs); + if (!transLp.transformSNFRelaxation(inds_, vals_, rhs_, snfr)) return false; + + // Array resized as each continuous col may create a non-zero bin coef + vals_.resize(2 * snfr.numNnzs); + inds_.resize(2 * snfr.numNnzs); + this->inds = inds_.data(); + this->vals = vals_.data(); + this->rhs = rhs_; + + if (!computeFlowCover()) return false; + if (!separateLiftedFlowCover()) return false; + inds_.resize(rowlen); + vals_.resize(rowlen); + rhs_ = static_cast(rhs); + if (!transLp.cleanup(inds_, vals_, rhs_, efficacy)) return false; + if (efficacy < 10 * feastol) return false; + return true; +} + bool HighsCutGeneration::finalizeAndAddCut(std::vector& inds_, std::vector& vals_, double& rhs_) { diff --git a/highs/mip/HighsCutGeneration.h b/highs/mip/HighsCutGeneration.h index cfe01b3a8c..6263044be9 100644 --- a/highs/mip/HighsCutGeneration.h +++ b/highs/mip/HighsCutGeneration.h @@ -65,6 +65,17 @@ class HighsCutGeneration { bool cmirCutGenerationHeuristic(double minEfficacy, bool onlyInitialCMIRScale = false); + bool computeFlowCover(); + + bool separateLiftedFlowCover(); + + bool preprocessSNFRelaxation(); + + bool tryGenerateFlowCoverCut(HighsTransformedLp& transLp, + std::vector& inds_, + std::vector& vals_, double& rhs_, + double& efficacy); + double scale(double val); bool postprocessCut(); @@ -92,7 +103,8 @@ class HighsCutGeneration { /// separates the LP solution for the given single row relaxation bool generateCut(HighsTransformedLp& transLp, std::vector& inds, std::vector& vals, double& rhs, - bool onlyInitialCMIRScale = false); + bool onlyInitialCMIRScale = false, + bool genFlowCover = false); /// generate a conflict from the given proof constraint which cuts of the /// given local domain @@ -103,6 +115,32 @@ class HighsCutGeneration { /// cutpool if it is violated enough bool finalizeAndAddCut(std::vector& inds, std::vector& vals, double& rhs); + + /// Single Node Flow Relaxation for flow cover cuts + struct SNFRelaxation { + HighsInt numNnzs; // |N-| + |N+| + std::vector coef; // (+-1) coefficient of col in SNFR + std::vector vubCoef; // u_j in y'_j <= u_j x_j in SNFR + std::vector binSolval; // lp[x_j], y'_j <= u_j x_j in SNFR + std::vector origBinCols; // orig x_i, y'_j <= u_j x_j in SNFR + std::vector origContCols; // orig y_i used to make y'_j in SNFR + std::vector aggrBinCoef; // coef of x_i in y'_j aggregation + std::vector aggrContCoef; // coef of y_i in y'_j aggrregation + std::vector aggrConstant; // constant shift in y'_j aggregation + std::vector complementation; // was the original bincol complemented + + std::vector + flowCoverStatus; // (+1) in flow cover (-1) notin flow cover + double rhs; // in \sum_{j \in N+} y'_j - \sum_{j \in N-} y'_j <= b + double lambda; // in sum_{j in C+} u_j - sum_{j in C-} u_j = b + lambda + }; + + private: + SNFRelaxation snfr; + void initSNFRelaxation(); + + public: + SNFRelaxation& getSNFRelaxation() { return snfr; } }; #endif diff --git a/highs/mip/HighsModkSeparator.cpp b/highs/mip/HighsModkSeparator.cpp index 19727670bd..2dc344b321 100644 --- a/highs/mip/HighsModkSeparator.cpp +++ b/highs/mip/HighsModkSeparator.cpp @@ -223,7 +223,7 @@ void HighsModkSeparator::separateLpSolution(HighsLpRelaxation& lpRelaxation, lpAggregator.getCurrentAggregation(inds, vals, false); rhs = 0.0; - cutGen.generateCut(transLp, inds, vals, rhs, true); + cutGen.generateCut(transLp, inds, vals, rhs, true, false); if (k != 2) { lpAggregator.clear(); @@ -237,7 +237,7 @@ void HighsModkSeparator::separateLpSolution(HighsLpRelaxation& lpRelaxation, lpAggregator.getCurrentAggregation(inds, vals, true); rhs = 0.0; - cutGen.generateCut(transLp, inds, vals, rhs, true); + cutGen.generateCut(transLp, inds, vals, rhs, true, false); lpAggregator.clear(); }; diff --git a/highs/mip/HighsPathSeparator.cpp b/highs/mip/HighsPathSeparator.cpp index 443a24b98b..7f8f7d204d 100644 --- a/highs/mip/HighsPathSeparator.cpp +++ b/highs/mip/HighsPathSeparator.cpp @@ -164,6 +164,9 @@ void HighsPathSeparator::separateLpSolution(HighsLpRelaxation& lpRelaxation, colOutArcs[col].second = outArcRows.size(); } + const bool genFlowCover = + mip.mipdata_->num_nodes - mip.mipdata_->num_nodes_before_run == 0 && + !mip.submip; HighsCutGeneration cutGen(lpRelaxation, cutpool); std::vector baseRowInds; std::vector baseRowVals; @@ -350,7 +353,8 @@ void HighsPathSeparator::separateLpSolution(HighsLpRelaxation& lpRelaxation, // generate cut double rhs = 0; - success = cutGen.generateCut(transLp, baseRowInds, baseRowVals, rhs); + success = cutGen.generateCut(transLp, baseRowInds, baseRowVals, rhs, + false, genFlowCover); lpAggregator.getCurrentAggregation(baseRowInds, baseRowVals, true); if (!aggregatedPath.empty() || bestOutArcCol != -1 || @@ -359,7 +363,8 @@ void HighsPathSeparator::separateLpSolution(HighsLpRelaxation& lpRelaxation, // generate reverse cut rhs = 0; - success |= cutGen.generateCut(transLp, baseRowInds, baseRowVals, rhs); + success |= cutGen.generateCut(transLp, baseRowInds, baseRowVals, rhs, + false, genFlowCover); if (success || (bestOutArcCol == -1 && bestInArcCol == -1)) break; diff --git a/highs/mip/HighsTableauSeparator.cpp b/highs/mip/HighsTableauSeparator.cpp index 3d4cd0d5a2..577cd94e45 100644 --- a/highs/mip/HighsTableauSeparator.cpp +++ b/highs/mip/HighsTableauSeparator.cpp @@ -188,6 +188,9 @@ void HighsTableauSeparator::separateLpSolution(HighsLpRelaxation& lpRelaxation, HighsInt numCuts = cutpool.getNumCuts(); const double bestScoreFac[] = {0.0025, 0.01}; + const bool genFlowCover = + mip.mipdata_->num_nodes - mip.mipdata_->num_nodes_before_run == 0 && + !mip.submip; for (const auto& fracvar : fractionalBasisvars) { if (cutpool.getNumCuts() - numCuts >= 1000) break; @@ -229,12 +232,12 @@ void HighsTableauSeparator::separateLpSolution(HighsLpRelaxation& lpRelaxation, baseRowInds.size()); double rhs = 0; - cutGen.generateCut(transLp, baseRowInds, baseRowVals, rhs); + cutGen.generateCut(transLp, baseRowInds, baseRowVals, rhs, false, genFlowCover); if (mip.mipdata_->domain.infeasible()) break; lpAggregator.getCurrentAggregation(baseRowInds, baseRowVals, true); rhs = 0; - cutGen.generateCut(transLp, baseRowInds, baseRowVals, rhs); + cutGen.generateCut(transLp, baseRowInds, baseRowVals, rhs, false, genFlowCover); if (mip.mipdata_->domain.infeasible()) break; lpAggregator.clear(); diff --git a/highs/mip/HighsTransformedLp.cpp b/highs/mip/HighsTransformedLp.cpp index 703b377e09..44432e8f93 100644 --- a/highs/mip/HighsTransformedLp.cpp +++ b/highs/mip/HighsTransformedLp.cpp @@ -561,3 +561,516 @@ bool HighsTransformedLp::untransform(std::vector& vals, return true; } + +// Create a single node flow relaxation (SNFR) from an aggregated +// mixed-integer row and find a valid flow cover. +// Turn \sum c_i x_i + \sum a_i y_i <= a_0 (x_i binary, y_i real non-neg) +// into \sum_{j \in N+} y'_j - \sum_{j \in N-} y'_j <= b, where y'_j <= u_j x_j +bool HighsTransformedLp::transformSNFRelaxation( + std::vector& inds, std::vector& vals, double& rhs, + HighsCutGeneration::SNFRelaxation& snfr) { + const HighsSolution& lpSolution = lprelaxation.getLpSolver().getSolution(); + + HighsCDouble tmpSnfrRhs = rhs; + + const HighsMipSolver& mip = lprelaxation.getMipSolver(); + const HighsInt slackOffset = lprelaxation.numCols(); + + HighsInt numNz = inds.size(); + HighsInt numBinCols = 0; + + auto getLb = [&](HighsInt col) { + return (col < slackOffset ? mip.mipdata_->domain.col_lower_[col] + : lprelaxation.slackLower(col - slackOffset)); + }; + + auto getUb = [&](HighsInt col) { + return (col < slackOffset ? mip.mipdata_->domain.col_upper_[col] + : lprelaxation.slackUpper(col - slackOffset)); + }; + + auto colIsBinary = [&](const HighsInt col, const double lb, const double ub) { + // Check local domain too. Global domain change may not yet be reflected + if (lprelaxation.isColIntegral(col) && lb == 0 && ub == 1 && + lprelaxation.colLower(col) >= 0 && lprelaxation.colUpper(col) <= 1) { + return true; + } + return false; + }; + + auto getLpSolution = [&](HighsInt col) { + HighsInt numCols = lprelaxation.numCols(); + if (col < numCols) { + return lpSolution.col_value[col]; + } else { + return lpSolution.row_value[col - numCols]; + } + }; + + auto remove = [&](HighsInt position) { + numNz--; + if (position < numNz - numBinCols - 1) { + std::swap(vals[position], vals[numNz - numBinCols]); + std::swap(vals[numNz - numBinCols], vals[numNz]); + std::swap(inds[position], inds[numNz - numBinCols]); + std::swap(inds[numNz - numBinCols], inds[numNz]); + } else { + inds[position] = inds[numNz]; + vals[position] = vals[numNz]; + } + inds[numNz] = 0; + vals[numNz] = 0; + numNz--; + }; + + auto checkValidityVB = [&](HighsInt bincol, HighsImplications::VarBound vb, + double coef, double origbincoef, double lb, + double ub, bool isVub, bool& complement, + bool& inclbincoef) { + // a_j coefficient of y_j in row. c_j coefficient of x_j in row. + // u_j, l_j, closest simple bounds imposed on y_j. + // y_j <= u'_j x_j + d_j || y_j >= l'_j x_j + d_j + // variable bound can only be used if following properties respected: + // |l'_j| <= 1e6 && |u'_j| <= 1e6 + // if a_j > 0 + // !isVub: (1) u_j <= d_j + // (2) a_j (u_j - d_j) + c_j <= 0 + // (3) a_j l'_j + c_j <= 0 + // isVub: (1) l_j >= d_j + // (2) a_j (l_j - d_j) + c_j >= 0 + // (3) a_j u'_j + c_j >= 0 + // if a_j < 0 + // !isVub: (1) u_j <= d_j + // (2) a_j (u_j - d_j) + c_j >= 0 + // (3) a_j l'_j + c_j >= 0 + // isVub: (1) l_j >= d_j + // (2) a_j (l_j - d_j) + c_j <= 0 + // (3) a_j u'_j + c_j <= 0 + // Note: c_j may change during transform if two columns use same bin col + // complement: Should x_j be replaced by x'_j = (1 - x_j) + // inclbincoef: Should c_j be used and included in the aggregation + if (bincol == -1) return false; + if (abs(vb.coef) >= 1e+6) return false; + complement = false; + inclbincoef = true; + if (isVub && lb < vb.constant) { + if (lb >= vb.constant + vb.coef) { + complement = true; + } else { + return false; + } + } + if (!isVub && ub > vb.constant) { + if (ub <= vb.constant + vb.coef) { + complement = true; + } else { + return false; + } + } + + const double sign = coef >= 0 ? 1 : -1; + double complorigbincoef = complement ? -origbincoef : origbincoef; + double vbconstant = complement ? vb.constant + vb.coef : vb.constant; + double vbcoef = complement ? -vb.coef : vb.coef; + double val1; + double val2; + + if (isVub) { + val1 = sign * (coef * vbcoef + complorigbincoef); + if (val1 > kHighsInf) return false; + if (val1 < 0) { + val1 -= sign * complorigbincoef; + if (val1 < 0) return false; + inclbincoef = false; + } + if (inclbincoef) { + val2 = sign * (coef * (lb - vbconstant) + complorigbincoef); + if (val2 < 0) { + val2 -= sign * complorigbincoef; + if (val2 < 0) return false; + inclbincoef = false; + val1 -= sign * complorigbincoef; + if (val1 < 0) return false; + } + } else { + val2 = sign * (coef * (lb - vbconstant)); + if (val2 < 0) return false; + } + } else { + val1 = sign * (coef * vbcoef + complorigbincoef); + if (-val1 > kHighsInf) return false; + if (val1 > 0) { + val1 -= sign * complorigbincoef; + if (val1 > 0) return false; + inclbincoef = false; + } + if (inclbincoef) { + val2 = sign * (coef * (ub - vbconstant) + complorigbincoef); + if (val2 > 0) { + val2 -= sign * complorigbincoef; + if (val2 > 0) return false; + inclbincoef = false; + val1 -= sign * complorigbincoef; + if (val1 > 0) return false; + } + } else { + val2 = sign * (coef * (ub - vbconstant)); + if (val2 > 0) return false; + } + } + + return true; + }; + + auto addSNFRentry = [&](HighsInt origbincol, HighsInt origcontcol, + double binsolval, HighsInt coef, double vubcoef, + double aggrconstant, double aggrbincoef, + double aggrcontcoef, bool complement) { + assert(binsolval >= -lprelaxation.getMipSolver().mipdata_->feastol && + binsolval <= 1 + lprelaxation.getMipSolver().mipdata_->feastol); + assert(vubcoef >= -1e-10); + snfr.origBinCols[snfr.numNnzs] = origbincol; + snfr.origContCols[snfr.numNnzs] = origcontcol; + snfr.binSolval[snfr.numNnzs] = binsolval; + snfr.coef[snfr.numNnzs] = coef; + snfr.vubCoef[snfr.numNnzs] = std::max(vubcoef, 0.0); + snfr.aggrConstant[snfr.numNnzs] = aggrconstant; + snfr.aggrBinCoef[snfr.numNnzs] = aggrbincoef; + snfr.aggrContCoef[snfr.numNnzs] = aggrcontcoef; + snfr.complementation[snfr.numNnzs] = complement; + snfr.numNnzs++; + }; + + // Place the non-binary columns to the front (all general ints relaxed) + // Use vectorsum to track original row coefficients of binary columns. + HighsInt i = 0; + while (i < numNz - numBinCols) { + HighsInt col = inds[i]; + double lb = getLb(col); + double ub = getUb(col); + if (colIsBinary(col, lb, ub)) { + numBinCols++; + vectorsum.add(col, vals[i]); + std::swap(inds[i], inds[numNz - numBinCols]); + std::swap(vals[i], vals[numNz - numBinCols]); + continue; + } + ++i; + } + + i = 0; + while (i < numNz) { + HighsInt col = inds[i]; + + double lb = getLb(col); + double ub = getUb(col); + + if (lb == -kHighsInf || ub == kHighsInf) { + vectorsum.clear(); + return false; + } + + if (ub - lb < mip.options_mip_->small_matrix_value) { + rhs -= std::min(lb, ub) * vals[i]; + tmpSnfrRhs -= std::min(lb, ub) * vals[i]; + remove(i); + continue; + } + + // We are using vectorsum to keep track of the binary coefficients + if (i >= numNz - numBinCols && vectorsum.getValue(col) == 0) { + ++i; + continue; + } + + // the code below uses the difference between the column upper and lower + // bounds as the upper bound for the slack from the variable upper bound + // constraint (upper[j] = ub - lb) and thus assumes that the variable upper + // bound constraints are tight. this assumption may not be satisfied when + // new bound changes were derived during cut generation and, therefore, we + // tighten the best variable upper bound. + if (bestVub[col].first != -1 && + bestVub[col].second.maxValue() > ub + mip.mipdata_->feastol) { + bool redundant = false; + bool infeasible = false; + mip.mipdata_->implications.cleanupVub(col, bestVub[col].first, + bestVub[col].second, ub, redundant, + infeasible, false); + } + + // the code below uses the difference between the column upper and lower + // bounds as the upper bound for the slack from the variable lower bound + // constraint (upper[j] = ub - lb) and thus assumes that the variable lower + // bound constraints are tight. this assumption may not be satisfied when + // new bound changes were derived during cut generation and, therefore, we + // tighten the best variable lower bound. + if (bestVlb[col].first != -1 && + bestVlb[col].second.minValue() < lb - mip.mipdata_->feastol) { + bool redundant = false; + bool infeasible = false; + mip.mipdata_->implications.cleanupVlb(col, bestVlb[col].first, + bestVlb[col].second, lb, redundant, + infeasible, false); + } + + // Transform entry into the SNFR + if (i >= numNz - numBinCols) { + // Binary columns can be added directly to the SNFR + if (vals[i] >= 0) { + addSNFRentry(col, -1, getLpSolution(col), 1, vals[i], 0, + vectorsum.getValue(col), 0, false); + } else { + addSNFRentry(col, -1, getLpSolution(col), -1, -vals[i], 0, + -vectorsum.getValue(col), 0, false); + } + } else { + // Decide whether to use {simple, variable} {lower, upper} bound + bool complementvlb = false; + bool inclbincolvlb = true; + bool vlbValid = checkValidityVB( + bestVlb[col].first, bestVlb[col].second, vals[i], + bestVlb[col].first == -1 ? 0 : vectorsum.getValue(bestVlb[col].first), + lb, ub, false, complementvlb, inclbincolvlb); + bool complementvub = false; + bool inclbincolvub = true; + bool vubValid = checkValidityVB( + bestVub[col].first, bestVub[col].second, vals[i], + bestVub[col].first == -1 ? 0 : vectorsum.getValue(bestVub[col].first), + lb, ub, true, complementvub, inclbincolvub); + auto boundType = BoundType::kSimpleLb; + if (lbDist[col] < ubDist[col] - mip.mipdata_->feastol && vlbValid) { + boundType = BoundType::kVariableLb; + } else if (ubDist[col] < lbDist[col] - mip.mipdata_->feastol && + vubValid) { + boundType = BoundType::kVariableUb; + } else if (vals[i] > 0 && vlbValid) { + boundType = BoundType::kVariableLb; + } else if (vals[i] < 0 && vubValid) { + boundType = BoundType::kVariableUb; + } else if (vlbValid) { + boundType = BoundType::kVariableLb; + } else if (vubValid) { + boundType = BoundType::kVariableUb; + } else if (lbDist[col] < ubDist[col] - mip.mipdata_->feastol) { + boundType = BoundType::kSimpleLb; + } else if (ubDist[col] < lbDist[col] - mip.mipdata_->feastol) { + boundType = BoundType::kSimpleUb; + } else if (vals[i] > 0) { + boundType = BoundType::kSimpleLb; + } else { + boundType = BoundType::kSimpleUb; + } + + double vbconstant; + double vbcoef; + double aggrvbcoef; + double aggrconstant; + double binsolval; + double bincoef; + HighsInt vbcol; + switch (boundType) { + case BoundType::kSimpleLb: + // Case (1) -> a_j > 0, y'_j -> N- Case (2) -> a_j < 0, y'_j -> N+ + // (1) y'_j = -a_j(y_j - u_j), 0 <= y'_j <= a_j(u_j - l_j)x_j, x_j=1 + // (2) y'_j = a_j(y_j - u_j), 0 <= y'_j <= -a_j(u_j - l_j)x_j, x_j=1 + // rhs -= a_j * u_j + aggrvbcoef = static_cast(vals[i] * (HighsCDouble(ub) - lb)); + aggrconstant = static_cast(HighsCDouble(vals[i]) * ub); + binsolval = std::min( + 1.0, std::max(0.0, (ub - getLpSolution(col)) / (ub - lb))); + if (vals[i] >= 0) { + addSNFRentry(-1, col, binsolval, -1, aggrvbcoef, aggrconstant, 0, + -vals[i], false); + } else { + addSNFRentry(-1, col, binsolval, 1, -aggrvbcoef, -aggrconstant, 0, + vals[i], false); + } + tmpSnfrRhs -= aggrconstant; + break; + case BoundType::kSimpleUb: + // Case (1) -> a_j > 0, y'_j -> N+ Case (2) -> a_j < 0, y'_j ->N- + // (1) y'_j = a_j(y_j - l_j), 0 <= y'_j <= a_j(u_j - l_j)x_j, x_j = 1 + // (2) y'_j = -a_j(y_j - l_j), 0 <= y'_j <= -a_j(u_j - l_j)x_j, + // x_j = 1 + // rhs -= a_j * l_j + aggrvbcoef = static_cast(vals[i] * (HighsCDouble(ub) - lb)); + aggrconstant = static_cast(HighsCDouble(vals[i]) * lb); + binsolval = std::min( + 1.0, std::max(0.0, (getLpSolution(col) - lb) / (ub - lb))); + if (vals[i] >= 0) { + addSNFRentry(-1, col, binsolval, 1, aggrvbcoef, -aggrconstant, 0, + vals[i], false); + } else { + addSNFRentry(-1, col, binsolval, -1, -aggrvbcoef, aggrconstant, 0, + -vals[i], false); + } + tmpSnfrRhs -= aggrconstant; + break; + case BoundType::kVariableLb: + // vlb: l'_j x_j + d_j <= y_j. c_j coef of x_j in row + // if compl: -l'_j x_j + (l'_j + d_j) <= y_j. -c_j coef of x_j in row + // Case (1) -> a_j > 0, y'_j -> N- Case (2) -> a_j < 0, y'_j ->N+ + // (1) y'_j = -(a_j(y_j - d_j) + c_j * x_j), + // 0 <= y'_j <= -(a_j l'_j + c_j)x_j + // (2) y'_j = a_j(y_j - d_j) + c_j * x_j, + // 0 <= y'_j <= (a_j l'_j + c_j)x_j + // rhs -= a_j * d_j + vbcol = bestVlb[col].first; + vbconstant = bestVlb[col].second.constant + + (complementvlb ? bestVlb[col].second.coef : 0); + vbcoef = complementvlb ? -bestVlb[col].second.coef + : bestVlb[col].second.coef; + bincoef = inclbincolvlb ? complementvlb ? -vectorsum.getValue(vbcol) + : vectorsum.getValue(vbcol) + : 0; + binsolval = lpSolution.col_value[vbcol]; + if (complementvlb) binsolval = 1 - binsolval; + aggrvbcoef = + static_cast(HighsCDouble(vals[i]) * vbcoef + bincoef); + aggrconstant = + static_cast(HighsCDouble(vals[i]) * vbconstant); + if (vals[i] >= 0) { + addSNFRentry(vbcol, col, binsolval, -1, -aggrvbcoef, aggrconstant, + -bincoef, -vals[i], complementvlb); + } else { + addSNFRentry(vbcol, col, binsolval, 1, aggrvbcoef, -aggrconstant, + bincoef, vals[i], complementvlb); + } + tmpSnfrRhs -= + aggrconstant + (complementvlb && inclbincolvlb ? -bincoef : 0); + if (inclbincolvlb) vectorsum.values[vbcol] = 0; + break; + case BoundType::kVariableUb: + // vub: y_j <= u'_j x_j + d_j. c_j coef of x_j in row + // if compl: y_j <= -u'_j x_j + (u'_j + d_j). -c_j coef of x_j in row + // Case (1) -> a_j > 0, y'_j -> N+ Case (2) -> a_j < 0, y'_j ->N- + // (1) y'_j = a_j(y_j - d_j) + c_j * x_j), + // 0 <= y'_j <= (a_j u'_j + c_j)x_j + // (2) y'_j = -(a_j(y_j - d_j) + c_j * x_j), + // 0 <= y'_j <= -(a_j u'_j + c_j)x_j + // rhs -= a_j * d_j + vbcol = bestVub[col].first; + vbconstant = bestVub[col].second.constant + + (complementvub ? bestVub[col].second.coef : 0); + vbcoef = complementvub ? -bestVub[col].second.coef + : bestVub[col].second.coef; + bincoef = inclbincolvub ? complementvub ? -vectorsum.getValue(vbcol) + : vectorsum.getValue(vbcol) + : 0; + binsolval = lpSolution.col_value[vbcol]; + if (complementvub) binsolval = 1 - binsolval; + aggrvbcoef = + static_cast(HighsCDouble(vals[i]) * vbcoef + bincoef); + aggrconstant = + static_cast(HighsCDouble(vals[i]) * vbconstant); + if (vals[i] >= 0) { + addSNFRentry(vbcol, col, binsolval, 1, aggrvbcoef, -aggrconstant, + bincoef, vals[i], complementvub); + } else { + addSNFRentry(vbcol, col, binsolval, -1, -aggrvbcoef, aggrconstant, + -bincoef, -vals[i], complementvub); + } + tmpSnfrRhs -= + aggrconstant + (complementvub && inclbincolvub ? -bincoef : 0); + if (inclbincolvub) vectorsum.values[vbcol] = 0; + break; + } + } + // move to next element + i++; + } + + vectorsum.clear(); + snfr.rhs = static_cast(tmpSnfrRhs); + if (numNz == 0) return false; + return true; +} + +// Remove slack, small coefficients, and calculate the efficacy +bool HighsTransformedLp::cleanup(std::vector& inds, + std::vector& vals, double& rhs, + double& efficacy) { + HighsCDouble tmpRhs = rhs; + const HighsMipSolver& mip = lprelaxation.getMipSolver(); + const HighsInt slackOffset = mip.numCol(); + + auto numNz = static_cast(inds.size()); + + for (HighsInt i = 0; i != numNz; ++i) { + if (vals[i] == 0.0) continue; + const HighsInt col = inds[i]; + if (col < slackOffset) { + vectorsum.add(col, vals[i]); + } else { + const HighsInt row = col - slackOffset; + HighsInt rowlen; + const HighsInt* rowinds; + const double* rowvals; + lprelaxation.getRow(row, rowlen, rowinds, rowvals); + + for (HighsInt j = 0; j != rowlen; ++j) + vectorsum.add(rowinds[j], vals[i] * rowvals[j]); + } + } + + bool abort = false; + auto IsZero = [&](HighsInt col, double val) { + assert(col < mip.numCol()); + double absval = std::abs(val); + if (absval <= mip.options_mip_->small_matrix_value) { + return true; + } + + if (absval <= mip.mipdata_->feastol) { + if (val > 0) { + if (mip.mipdata_->domain.col_lower_[col] == -kHighsInf) + abort = true; + else + tmpRhs -= val * mip.mipdata_->domain.col_lower_[col]; + } else { + if (mip.mipdata_->domain.col_upper_[col] == kHighsInf) + abort = true; + else + tmpRhs -= val * mip.mipdata_->domain.col_upper_[col]; + } + return true; + } + return false; + }; + + vectorsum.cleanup(IsZero); + if (abort) { + vectorsum.clear(); + return false; + } + rhs = static_cast(tmpRhs); + + inds = vectorsum.getNonzeros(); + numNz = static_cast(inds.size()); + vals.resize(numNz); + for (HighsInt i = 0; i != numNz; ++i) vals[i] = vectorsum.getValue(inds[i]); + vectorsum.clear(); + + double viol = -rhs; + double sqrnorm = 0; + const std::vector& lpSolution = lprelaxation.getSolution().col_value; + for (HighsInt i = 0; i != numNz; ++i) { + HighsInt col = inds[i]; + viol += vals[i] * lpSolution[col]; + if (vals[i] >= 0 && + lpSolution[col] <= + mip.mipdata_->domain.col_lower_[col] + mip.mipdata_->feastol) + continue; + if (vals[i] < 0 && lpSolution[col] >= mip.mipdata_->domain.col_upper_[col] - + mip.mipdata_->feastol) + continue; + sqrnorm += vals[i] * vals[i]; + } + if (sqrnorm == 0) { + efficacy = 0; + } else { + efficacy = viol / sqrt(sqrnorm); + } + + return true; +} diff --git a/highs/mip/HighsTransformedLp.h b/highs/mip/HighsTransformedLp.h index 5c3d2d01db..944ee9c947 100644 --- a/highs/mip/HighsTransformedLp.h +++ b/highs/mip/HighsTransformedLp.h @@ -16,6 +16,7 @@ #include +#include "HighsCutGeneration.h" #include "lp_data/HConst.h" #include "mip/HighsImplications.h" #include "util/HighsCDouble.h" @@ -58,6 +59,13 @@ class HighsTransformedLp { bool untransform(std::vector& vals, std::vector& inds, double& rhs, bool integral = false); + + bool transformSNFRelaxation(std::vector& inds, + std::vector& vals, double& rhs, + HighsCutGeneration::SNFRelaxation& snfr); + + bool cleanup(std::vector& inds, std::vector& vals, + double& rhs, double& efficacy); }; #endif