diff --git a/highs/lp_data/HighsOptions.h b/highs/lp_data/HighsOptions.h index 063918be135..06298d71445 100644 --- a/highs/lp_data/HighsOptions.h +++ b/highs/lp_data/HighsOptions.h @@ -477,6 +477,8 @@ struct HighsOptionsStruct { bool less_infeasible_DSE_check; bool less_infeasible_DSE_choose_row; bool use_original_HFactor_logic; + bool presolve_dfprobing; + bool presolve_gdf; // bool allow_pdlp_cleanup; bool run_centring; HighsInt max_centring_steps; @@ -1763,6 +1765,17 @@ class HighsOptions : public HighsOptionsStruct { advanced, ¢ring_ratio_tolerance, 0, 100, kHighsInf); records.push_back(record_double); + record_bool = new OptionRecordBool( + "presolve_dfprobing", + "Use the dual fixing aumgented probing technique in presolve", advanced, + &presolve_dfprobing, true); + records.push_back(record_bool); + + record_bool = new OptionRecordBool( + "presolve_gdf", "Use the generalized dual fixing technique in presolve", + advanced, &presolve_gdf, true); + records.push_back(record_bool); + // Set up the log_options aliases log_options.clear(); log_options.log_stream = diff --git a/highs/mip/HighsDomain.cpp b/highs/mip/HighsDomain.cpp index 7bbcc3ef174..ec7c698c27a 100644 --- a/highs/mip/HighsDomain.cpp +++ b/highs/mip/HighsDomain.cpp @@ -75,6 +75,7 @@ HighsDomain::HighsDomain(HighsMipSolver& mipsolver) : mipsolver(&mipsolver) { changedcols_.reserve(mipsolver.numCol()); infeasible_reason = Reason::unspecified(); infeasible_ = false; + dfprobingPropagation.domain = this; } void HighsDomain::addCutpool(HighsCutPool& cutpool) { @@ -637,6 +638,590 @@ void HighsDomain::CutpoolPropagation::updateActivityUbChange( } } +HighsDomain::DualfixingProbingPropagation::DualfixingProbingPropagation( + const DualfixingProbingPropagation& other) + : redundantPropagateFlag_(other.redundantPropagateFlag_), + redundantPropagateVec_(other.redundantPropagateVec_), + zeroCostVarsDirection_(other.zeroCostVarsDirection_), + zeroCostFixedVariables_(other.zeroCostFixedVariables_), + colLowerLockOriginal_(other.colLowerLockOriginal_), + colUpperLockOriginal_(other.colUpperLockOriginal_), + colLowerLockReduced_(other.colLowerLockReduced_), + colUpperLockReduced_(other.colUpperLockReduced_), + candidatesVec_(other.candidatesVec_), + candidatesFlag_(other.candidatesFlag_), + lockNeedClear_(other.lockNeedClear_), + gdfCandidatesVec_(other.gdfCandidatesVec_), + gdfCandidatesFlag_(other.gdfCandidatesFlag_), + gdfLbReachable0_(other.gdfLbReachable0_), + gdfLbReachable1_(other.gdfLbReachable1_), + gdfUbReachable0_(other.gdfUbReachable0_), + gdfUbReachable1_(other.gdfUbReachable1_) { + ; +} + +void HighsDomain::DualfixingProbingPropagation::recomputeLocks() { + mipsolver = domain->mipsolver; + redundantPropagateFlag_.assign(2 * mipsolver->numRow(), false); + redundantPropagateVec_.clear(); + redundantPropagateVec_.reserve(2 * mipsolver->numRow()); + zeroCostVarsDirection_.assign(2 * mipsolver->numCol(), + FIXDIRECTION_NOT_DECIDED); + zeroCostFixedVariables_.clear(); + zeroCostFixedVariables_.reserve(2 * mipsolver->numCol()); + + startZeroCostFixing_ = false; + previousSize_ = 0; + + colLowerLockOriginal_.assign(mipsolver->numCol(), 0); + colUpperLockOriginal_.assign(mipsolver->numCol(), 0); + colLowerLockReduced_.assign(mipsolver->numCol(), 0); + colUpperLockReduced_.assign(mipsolver->numCol(), 0); + + candidatesVec_.clear(); + candidatesVec_.reserve(mipsolver->numCol()); + candidatesFlag_.assign(mipsolver->numCol(), false); + lockNeedClear_.clear(); + lockNeedClear_.reserve(mipsolver->numCol()); + + gdfCandidatesVec_.clear(); + gdfCandidatesVec_.reserve(mipsolver->numCol()); + gdfCandidatesFlag_.assign(mipsolver->numCol(), false); + + gdfLbReachable0_.assign(mipsolver->numCol(), 0); + gdfLbReachable1_.assign(mipsolver->numCol(), 0); + gdfUbReachable0_.assign(mipsolver->numCol(), 0); + gdfUbReachable1_.assign(mipsolver->numCol(), 0); + + // compute the original locks for each variable + const auto model = mipsolver->model_; + for (HighsInt iCol = 0; iCol < model->a_matrix_.num_col_; iCol++) { + for (HighsInt k = model->a_matrix_.start_[iCol]; + k < model->a_matrix_.start_[iCol + 1]; k++) { + const HighsInt iRow = model->a_matrix_.index_[k]; + const double iValue = model->a_matrix_.value_[k]; + const double lhs = model->row_lower_[iRow], rhs = model->row_upper_[iRow]; + if ((iValue > 0 && rhs != kHighsInf) || (iValue < 0 && lhs != -kHighsInf)) + colUpperLockOriginal_[iCol]++; + if ((iValue > 0 && lhs != -kHighsInf) || (iValue < 0 && rhs != kHighsInf)) + colLowerLockOriginal_[iCol]++; + } + } +} + +void HighsDomain::DualfixingProbingPropagation::updateRhsRedundant( + HighsInt row) { + if (!isEnabled()) return; + + if (domain->activitymaxinf_[row] != 0 || + redundantPropagateFlag_[2 * row + 1] || + mipsolver->model_->row_upper_[row] == kHighsInf) + return; + + if (domain->getMaxActivity(row) <= + mipsolver->model_->row_upper_[row] + mipsolver->mipdata_->feastol) { + redundantPropagateVec_.push_back(2 * row + 1); + redundantPropagateFlag_[2 * row + 1] = 1; + } +} + +void HighsDomain::DualfixingProbingPropagation::updateLhsRedundant( + HighsInt row) { + if (!isEnabled()) return; + + if (domain->activitymininf_[row] != 0 || redundantPropagateFlag_[2 * row] || + mipsolver->model_->row_lower_[row] == -kHighsInf) + return; + + if (domain->getMinActivity(row) >= + mipsolver->model_->row_lower_[row] - mipsolver->mipdata_->feastol) { + redundantPropagateVec_.push_back(2 * row); + redundantPropagateFlag_[2 * row] = 1; + } +} + +void HighsDomain::DualfixingProbingPropagation::propagate() { + // The boolean variable ``startZeroCostFixing_'' is used to flag if we allow + // variables with zero cost can be fixed in domain propagation. The process of + // domain propagtion in probing is executed in two phases: + // Phase 1: Apply classic domain propagation, and additionally fix + // variables with non-zero objective coefficients using dual fixing Phase + // 2: Apply classic domain propagation, and additionally fix variables + // (including those with zero objective coefficients) using dual fixing + // In Phase 1, ``startZeroCostFixing_'' is set to be ``false'' to exclude + // variable with zero objective coefficients. In Phase 2, + // ``startZeroCostFixing_'' is set to be ``true''. Note that + // (1) For all the bound changes in Phase 1, reductions deduced from them + // are valid for all optimal solutions; (2) For the bound changes in Phase + // 2, reductions deduced from them can only be used to derive global valid + // reductions (i.e., variable fixing, global bound tightening, and variable + // substitution). + if (!isEnabled()) return; + + assert(candidatesVec_.empty()); + vector domainchangeDFProbing; + + // tool lambda functions + auto addToCandidate = [&](HighsInt k) { + if (candidatesFlag_[k]) + return; + else { + candidatesVec_.push_back(k); + candidatesFlag_[k] = true; + } + }; + + // debug functions to check locks + auto checkVariableLowerLock = [&](HighsInt iCol) { + auto model = mipsolver->model_; + if (ableToFixToLb(iCol)) { + for (HighsInt k = model->a_matrix_.start_[iCol]; + k < model->a_matrix_.start_[iCol + 1]; k++) { + const HighsInt iRow = model->a_matrix_.index_[k]; + const double iValue = model->a_matrix_.value_[k]; + const double blower = model->row_lower_[iRow], + bupper = model->row_upper_[iRow]; + const bool lhsOk = iValue > 0 && domain->getMinActivity(iRow) >= + blower - domain->feastol(); + const bool rhsOk = iValue < 0 && domain->getMaxActivity(iRow) <= + bupper + domain->feastol(); + if (!lhsOk && !rhsOk) { + std::cout << "Lower lock: variable " << iCol << " at row = " << iRow + << " coef = " << iValue << " not redundant at constraint " + << iRow << ", minact = " << domain->getMinActivity(iRow) + << ", maxact = " << domain->getMaxActivity(iRow) + << " lhs = " << blower << " rhs = " << bupper << std::endl; + } + } + } + }; + + auto checkVariableUpperLock = [&](HighsInt iCol) { + auto model = mipsolver->model_; + if (ableToFixToUb(iCol)) { + for (HighsInt k = model->a_matrix_.start_[iCol]; + k < model->a_matrix_.start_[iCol + 1]; k++) { + const HighsInt iRow = model->a_matrix_.index_[k]; + const double iValue = model->a_matrix_.value_[k]; + const double blower = model->row_lower_[iRow], + bupper = model->row_upper_[iRow]; + const bool lhsOk = iValue < 0 && domain->getMinActivity(iRow) >= + blower - domain->feastol(); + const bool rhsOk = iValue > 0 && domain->getMaxActivity(iRow) <= + bupper + domain->feastol(); + if (!lhsOk && !rhsOk) { + std::cout << "Upper lock: variable " << iCol << " at row = " << iRow + << " coef = " << iValue << " not redundant at constraint " + << iRow << ", minact = " << domain->getMinActivity(iRow) + << ", maxact = " << domain->getMaxActivity(iRow) + << " lhs = " << blower << " rhs = " << bupper << std::endl; + } + } + } + }; + + auto addFixLower = [&](int iCol) { + HighsDomainChange* thisbchg = new HighsDomainChange; + thisbchg->column = iCol; + thisbchg->boundtype = HighsBoundType::kUpper; + thisbchg->boundval = domain->col_lower_[iCol]; + domainchangeDFProbing.push_back(thisbchg); + }; + + auto addFixUpper = [&](int iCol) { + HighsDomainChange* thisbchg = new HighsDomainChange; + thisbchg->column = iCol; + thisbchg->boundtype = HighsBoundType::kLower; + thisbchg->boundval = domain->col_upper_[iCol]; + domainchangeDFProbing.push_back(thisbchg); + }; + + // only record - we do not actually fix them now as their objective + // coefficients are zero + auto collectFixLower = [&](int iCol) { + zeroCostFixedVariables_.emplace_back(iCol, FIXDIRECTION_LOWER_BOUND); + }; + + auto collectFixUpper = [&](int iCol) { + zeroCostFixedVariables_.emplace_back(iCol, FIXDIRECTION_UPPER_BOUND); + }; + + // exit if no new redundant constraints are found + HighsInt maxLockLeft = redundantPropagateVec_.size() - previousSize_; + if (maxLockLeft == 0) return; + + for (; previousSize_ < redundantPropagateVec_.size(); + ++previousSize_, --maxLockLeft) { + const HighsInt i = redundantPropagateVec_[previousSize_]; + const HighsInt iRow = i / 2; + assert(iRow < mipsolver->numRow()); + + if (i % 2 == 0) { // lower redundant + HighsInt rstart = mipsolver->mipdata_->ARstart_[iRow]; + HighsInt rend = mipsolver->mipdata_->ARstart_[iRow + 1]; + for (auto k = rstart; k < rend; ++k) { + const HighsInt iCol = mipsolver->mipdata_->ARindex_[k]; + if (domain->isFixed(iCol)) continue; + const double iValue = mipsolver->mipdata_->ARvalue_[k]; + const double cost = mipsolver->model_->col_cost_[iCol]; + + // do not insert to candidates if the lock is not reduced enough + bool lowerNoInsert = colLowerLockReduced_[iCol] + maxLockLeft < + colLowerLockOriginal_[iCol]; + bool upperNoInsert = colUpperLockReduced_[iCol] + maxLockLeft < + colUpperLockOriginal_[iCol]; + + if (iValue > 0 && + cost >= mipsolver->options_mip_->dual_feasibility_tolerance) { + lockNeedClear_.insert(iCol); + colLowerLockReduced_[iCol]++; + lowerNoInsert = + lowerNoInsert && colLowerLockReduced_[iCol] + maxLockLeft < + colLowerLockOriginal_[iCol]; + } else if (iValue < 0 && + cost <= + mipsolver->options_mip_->dual_feasibility_tolerance) { + lockNeedClear_.insert(iCol); + colUpperLockReduced_[iCol]++; + upperNoInsert = + upperNoInsert && colUpperLockReduced_[iCol] + maxLockLeft < + colUpperLockOriginal_[iCol]; + } + + if (!lowerNoInsert || !upperNoInsert) addToCandidate(iCol); + } + } else { // upper redundant + HighsInt rstart = mipsolver->mipdata_->ARstart_[iRow]; + HighsInt rend = mipsolver->mipdata_->ARstart_[iRow + 1]; + for (auto k = rstart; k < rend; k++) { + const HighsInt iCol = mipsolver->mipdata_->ARindex_[k]; + if (domain->isFixed(iCol)) continue; + const double iValue = mipsolver->mipdata_->ARvalue_[k]; + const double cost = mipsolver->model_->col_cost_[iCol]; + + // do not insert to candidates if the lock is not reduced enough + bool lowerNoInsert = colLowerLockReduced_[iCol] + maxLockLeft < + colLowerLockOriginal_[iCol]; + bool upperNoInsert = colUpperLockReduced_[iCol] + maxLockLeft < + colUpperLockOriginal_[iCol]; + + if (iValue < 0 && + cost >= mipsolver->options_mip_->dual_feasibility_tolerance) { + lockNeedClear_.insert(iCol); + colLowerLockReduced_[iCol]++; + lowerNoInsert = + lowerNoInsert && colLowerLockReduced_[iCol] + maxLockLeft < + colLowerLockOriginal_[iCol]; + } else if (iValue > 0 && + cost <= + mipsolver->options_mip_->dual_feasibility_tolerance) { + lockNeedClear_.insert(iCol); + colUpperLockReduced_[iCol]++; + upperNoInsert = + upperNoInsert && colUpperLockReduced_[iCol] + maxLockLeft < + colUpperLockOriginal_[iCol]; + } + + if (!lowerNoInsert || !upperNoInsert) addToCandidate(iCol); + } + } + } + + for (auto iCol : candidatesVec_) { + if (domain->isFixed(iCol)) continue; + const bool canBeFixedToLower = + colLowerLockReduced_[iCol] == colLowerLockOriginal_[iCol]; + const bool canBeFixedToUpper = + colUpperLockReduced_[iCol] == colUpperLockOriginal_[iCol]; + if (!canBeFixedToLower && !canBeFixedToUpper) continue; + + if (fabs(mipsolver->model_->col_cost_[iCol]) <= + mipsolver->options_mip_->dual_feasibility_tolerance) { + if (startZeroCostFixing_) { + // not fixed before + if (zeroCostVarsDirection_[iCol] == FIXDIRECTION_NOT_DECIDED) { + // both directions are ok - depending on cost (no tolerance) + if (canBeFixedToLower && canBeFixedToUpper) { + if (mipsolver->model_->col_cost_[iCol] >= 0) { + addFixLower(iCol); + zeroCostVarsDirection_[iCol] = FIXDIRECTION_LOWER_BOUND; + } else { + addFixUpper(iCol); + zeroCostVarsDirection_[iCol] = FIXDIRECTION_UPPER_BOUND; + } + } + // fix depending on the direction + else if (canBeFixedToLower) { + addFixLower(iCol); + zeroCostVarsDirection_[iCol] = FIXDIRECTION_LOWER_BOUND; + } else if (canBeFixedToUpper) { + addFixUpper(iCol); + zeroCostVarsDirection_[iCol] = FIXDIRECTION_UPPER_BOUND; + } + } + // fix to lb + else if (zeroCostVarsDirection_[iCol] == FIXDIRECTION_LOWER_BOUND && + canBeFixedToLower) + addFixLower(iCol); + // fix to ub + else if (zeroCostVarsDirection_[iCol] == FIXDIRECTION_UPPER_BOUND && + canBeFixedToUpper) + addFixUpper(iCol); + + continue; + } + // do not perfrom zero cost variable fixing, just collect them and choose + // directions + else { + // not fixed before + if (zeroCostVarsDirection_[iCol] == FIXDIRECTION_NOT_DECIDED) { + // both directions are ok - depending on cost (no tolerance) + if (canBeFixedToLower && canBeFixedToUpper) { + if (mipsolver->model_->col_cost_[iCol] >= 0) { + collectFixLower(iCol); + zeroCostVarsDirection_[iCol] = FIXDIRECTION_LOWER_BOUND; + } else { + collectFixUpper(iCol); + zeroCostVarsDirection_[iCol] = FIXDIRECTION_UPPER_BOUND; + } + } else if (canBeFixedToLower) { // fix to lower and set its direction + collectFixLower(iCol); + zeroCostVarsDirection_[iCol] = FIXDIRECTION_LOWER_BOUND; + } else if (canBeFixedToUpper) { + collectFixUpper(iCol); + zeroCostVarsDirection_[iCol] = FIXDIRECTION_UPPER_BOUND; + } + } else if (zeroCostVarsDirection_[iCol] == FIXDIRECTION_UPPER_BOUND && + canBeFixedToUpper) { // fix to upper + collectFixUpper(iCol); + } else if (zeroCostVarsDirection_[iCol] == FIXDIRECTION_LOWER_BOUND && + canBeFixedToLower) { // fix to lower + collectFixLower(iCol); + } + // we have collected this column + continue; + } + } + + if (mipsolver->model_->col_cost_[iCol] >= + mipsolver->options_mip_->dual_feasibility_tolerance) { + if (canBeFixedToLower) { + // checkVariableLowerLock(iCol); + addFixLower(iCol); + continue; + } + } + + if (mipsolver->model_->col_cost_[iCol] <= + mipsolver->options_mip_->dual_feasibility_tolerance) { + if (canBeFixedToUpper) { + // checkVariableUpperLock(iCol); + addFixUpper(iCol); + continue; + } + } + } + + // clear candidate info + for (const auto x : candidatesVec_) { + candidatesFlag_[x] = false; + } + candidatesVec_.clear(); + + // change bound + size_t j = 0; + for (; j != domainchangeDFProbing.size() && !domain->infeasible_; ++j) { + domain->changeBound(*domainchangeDFProbing[j], Reason::unspecified()); + delete domainchangeDFProbing[j]; + } + + // clear the remaining domain changes if infeasible + for (j++; j < domainchangeDFProbing.size(); ++j) { + assert(domain->infeasible_); + delete domainchangeDFProbing[j]; + } + + // record the current number of redundant constraints + previousSize_ = redundantPropagateVec_.size(); +} + +void HighsDomain::DualfixingProbingPropagation::updateGDFInfo( + HighsInt probing_variable, bool val) { + // tool lambda functions + auto addToCandidate = [&](HighsInt k) { + if (gdfCandidatesFlag_[k]) + return; + else { + gdfCandidatesVec_.push_back(k); + gdfCandidatesFlag_[k] = true; + } + }; + + // only redundant constraints are useful in GDF + for (const auto x : redundantPropagateVec_) { + const HighsInt iRow = x / 2; + const bool isRhs = x % 2; + HighsInt rstart = mipsolver->mipdata_->ARstart_[iRow]; + HighsInt rend = mipsolver->mipdata_->ARstart_[iRow + 1]; + + for (auto k = rstart; k < rend; k++) { + const HighsInt iCol = mipsolver->mipdata_->ARindex_[k]; + const double iValue = mipsolver->mipdata_->ARvalue_[k]; + const double cost = mipsolver->model_->col_cost_[iCol]; + bool considered = false; + if (mipsolver->model_->col_lower_[iCol] == + mipsolver->model_->col_upper_[iCol] || + mipsolver->mipdata_->implications.colsubstituted[iCol]) + continue; + + if (iValue > 0) { + if (isRhs) { // consider upper bound reachable + const double globalUb = mipsolver->model_->col_upper_[iCol]; + const double probingUb = domain->col_upper_[iCol]; + if (!ableToFixToUb(iCol) || domain->getMaxActivity(iRow) == kHighsInf) + continue; + const bool upper_bound_reachable = + domain->getMaxActivity(iRow) + iValue * (globalUb - probingUb) <= + mipsolver->model_->row_upper_[iRow] + domain->feastol(); + if (upper_bound_reachable) { + considered = true; + // special treat if the current variable is the probing variable + if (iCol == probing_variable && val == 0) { + gdfUbReachable0_[iCol]++; + gdfUbReachable1_[iCol]++; + } else { + if (val == 0) gdfUbReachable0_[iCol]++; + if (val == 1) gdfUbReachable1_[iCol]++; + } + } + } else { // consider lower bound reachable + const double globalLb = mipsolver->model_->col_lower_[iCol]; + const double probingLb = domain->col_lower_[iCol]; + if (!ableToFixToLb(iCol) || + domain->getMinActivity(iRow) == -kHighsInf) + continue; + const bool lower_bound_reachable = + domain->getMinActivity(iRow) + iValue * (globalLb - probingLb) >= + mipsolver->model_->row_lower_[iRow] - domain->feastol(); + if (lower_bound_reachable) { + considered = true; + // special treat if the current variable is the probing variable + if (iCol == probing_variable && val == 1) { + gdfLbReachable0_[iCol]++; + gdfLbReachable1_[iCol]++; + } else { + if (val == 0) gdfLbReachable0_[iCol]++; + if (val == 1) gdfLbReachable1_[iCol]++; + } + } + } + } + + else { + if (isRhs) { // consider lower bound reachable + const double globalLb = mipsolver->model_->col_lower_[iCol]; + const double probingLb = domain->col_lower_[iCol]; + if (!ableToFixToLb(iCol) || domain->getMaxActivity(iRow) == kHighsInf) + continue; + const bool lower_bound_reachable = + domain->getMaxActivity(iRow) + iValue * (globalLb - probingLb) <= + mipsolver->model_->row_upper_[iRow] + domain->feastol(); + if (lower_bound_reachable) { + considered = true; + // special treat if the current variable is the probing variable + if (iCol == probing_variable && val == 1) { + gdfLbReachable0_[iCol]++; + gdfLbReachable1_[iCol]++; + } else { + if (val == 0) gdfLbReachable0_[iCol]++; + if (val == 1) gdfLbReachable1_[iCol]++; + } + } + } else { // consider upper bound reachable + const double globalUb = mipsolver->model_->col_upper_[iCol]; + const double probingUb = domain->col_upper_[iCol]; + if (!ableToFixToUb(iCol) || + domain->getMinActivity(iRow) == -kHighsInf) + continue; + const bool upper_bound_reachable = + domain->getMinActivity(iRow) + iValue * (globalUb - probingUb) >= + mipsolver->model_->row_lower_[iRow] - domain->feastol(); + if (upper_bound_reachable) { + considered = true; + // special treat if the current variable is the probing variable + if (iCol == probing_variable && val == 0) { + gdfUbReachable0_[iCol]++; + gdfUbReachable1_[iCol]++; + } else { + if (val == 0) gdfUbReachable0_[iCol]++; + if (val == 1) gdfUbReachable1_[iCol]++; + } + } + } + } + + if (considered) addToCandidate(iCol); + } + } +} + +HighsInt HighsDomain::DualfixingProbingPropagation::processGDFFixing() { + std::vector gdfFixingStack_; + + // derive global fixings from the GDF information + for (const auto iCol : gdfCandidatesVec_) { + const HighsInt lowerLock = colLowerLockOriginal_[iCol]; + const HighsInt upperLock = colUpperLockOriginal_[iCol]; + if (ableToFixToLb(iCol) && lowerLock > 0 && + gdfLbReachable0_[iCol] == lowerLock && + gdfLbReachable1_[iCol] == lowerLock) { + HighsDomainChange* thisbchg = new HighsDomainChange; + thisbchg->column = iCol; + thisbchg->boundtype = HighsBoundType::kUpper; + thisbchg->boundval = domain->col_lower_[iCol]; + gdfFixingStack_.push_back(thisbchg); + } + // a variable cannot be fixed to lb and ub simultaneously + else if (ableToFixToUb(iCol) && upperLock > 0 && + gdfUbReachable0_[iCol] == upperLock && + gdfUbReachable1_[iCol] == upperLock) { + HighsDomainChange* thisbchg = new HighsDomainChange; + thisbchg->column = iCol; + thisbchg->boundtype = HighsBoundType::kLower; + thisbchg->boundval = domain->col_upper_[iCol]; + gdfFixingStack_.push_back(thisbchg); + } + } + + // apply bound change + size_t j = 0; + for (; j != gdfFixingStack_.size() && !domain->infeasible_; ++j) { + domain->changeBound(*gdfFixingStack_[j], Reason::unspecified()); + delete gdfFixingStack_[j]; + } + + // clear the remaining domain changes if infeasible + for (; j < gdfFixingStack_.size(); ++j) { + assert(domain->infeasible_); + delete gdfFixingStack_[j]; + } + + gdfFixingStack_.clear(); + + return (HighsInt)j; +} + +void HighsDomain::DualfixingProbingPropagation::clearGDFInfo() { + for (const auto x : gdfCandidatesVec_) { + gdfLbReachable0_[x] = 0; + gdfLbReachable1_[x] = 0; + gdfUbReachable0_[x] = 0; + gdfUbReachable1_[x] = 0; + gdfCandidatesFlag_[x] = false; + } + gdfCandidatesVec_.clear(); +} + namespace highs { template <> struct RbTreeTraits< @@ -1554,12 +2139,18 @@ void HighsDomain::updateActivityLbChange(HighsInt col, double oldbound, assert(tmpinf == activitymininf_[mip->a_matrix_.index_[i]]); } #endif - + // If dfprobingPropagation.isZeroObjFixingEnabled() is true, + // then we cannot record redundant rows for lifting, as this bound change + // could disregarded. if (recordRedundantRows_ && + !dfprobingPropagation.isZeroObjFixingEnabled() && mip->row_lower_[mip->a_matrix_.index_[i]] != -kHighsInf && mip->row_upper_[mip->a_matrix_.index_[i]] == kHighsInf) updateRedundantRows(mip->a_matrix_.index_[i]); + if (newbound >= oldbound + mipsolver->mipdata_->feastol) + dfprobingPropagation.updateLhsRedundant(mip->a_matrix_.index_[i]); + if (deltamin <= 0) { updateThresholdLbChange(col, newbound, mip->a_matrix_.value_[i], capacityThreshold_[mip->a_matrix_.index_[i]]); @@ -1603,12 +2194,18 @@ void HighsDomain::updateActivityLbChange(HighsInt col, double oldbound, assert(tmpinf == activitymaxinf_[mip->a_matrix_.index_[i]]); } #endif - + // If dfprobingPropagation.isZeroObjFixingEnabled() is true, + // then we cannot record redundant rows for lifting, as this bound change + // could disregarded. if (recordRedundantRows_ && + !dfprobingPropagation.isZeroObjFixingEnabled() && mip->row_lower_[mip->a_matrix_.index_[i]] == -kHighsInf && mip->row_upper_[mip->a_matrix_.index_[i]] != kHighsInf) updateRedundantRows(mip->a_matrix_.index_[i]); + if (newbound >= oldbound + mipsolver->mipdata_->feastol) + dfprobingPropagation.updateRhsRedundant(mip->a_matrix_.index_[i]); + if (deltamax >= 0) { updateThresholdLbChange(col, newbound, mip->a_matrix_.value_[i], capacityThreshold_[mip->a_matrix_.index_[i]]); @@ -1721,12 +2318,18 @@ void HighsDomain::updateActivityUbChange(HighsInt col, double oldbound, assert(tmpinf == activitymaxinf_[mip->a_matrix_.index_[i]]); } #endif - + // If dfprobingPropagation.isZeroObjFixingEnabled() is true, + // then we cannot record redundant rows for lifting, as this bound change + // could disregarded. if (recordRedundantRows_ && + !dfprobingPropagation.isZeroObjFixingEnabled() && mip->row_lower_[mip->a_matrix_.index_[i]] == -kHighsInf && mip->row_upper_[mip->a_matrix_.index_[i]] != kHighsInf) updateRedundantRows(mip->a_matrix_.index_[i]); + if (newbound <= oldbound - mipsolver->mipdata_->feastol) + dfprobingPropagation.updateRhsRedundant(mip->a_matrix_.index_[i]); + if (deltamax >= 0) { updateThresholdUbChange(col, newbound, mip->a_matrix_.value_[i], capacityThreshold_[mip->a_matrix_.index_[i]]); @@ -1773,12 +2376,18 @@ void HighsDomain::updateActivityUbChange(HighsInt col, double oldbound, assert(tmpinf == activitymininf_[mip->a_matrix_.index_[i]]); } #endif - + // If dfprobingPropagation.isZeroObjFixingEnabled() is true, + // then we cannot record redundant rows for lifting, as this bound change + // could disregarded. if (recordRedundantRows_ && + !dfprobingPropagation.isZeroObjFixingEnabled() && mip->row_lower_[mip->a_matrix_.index_[i]] != -kHighsInf && mip->row_upper_[mip->a_matrix_.index_[i]] == kHighsInf) updateRedundantRows(mip->a_matrix_.index_[i]); + if (newbound <= oldbound - mipsolver->mipdata_->feastol) + dfprobingPropagation.updateLhsRedundant(mip->a_matrix_.index_[i]); + if (deltamin <= 0) { updateThresholdUbChange(col, newbound, mip->a_matrix_.value_[i], capacityThreshold_[mip->a_matrix_.index_[i]]); @@ -2372,6 +2981,10 @@ bool HighsDomain::propagate() { if (!conflictprop.propagateConflictInds_.empty()) return true; } + if (!infeasible_ && dfprobingPropagation.isActive() && + mipsolver->options_mip_->presolve_dfprobing) + return true; + return false; }; @@ -2547,6 +3160,19 @@ bool HighsDomain::propagate() { propagateinds.clear(); } } + + if (!infeasible_ && dfprobingPropagation.isActive() && + mipsolver->options_mip_->presolve_dfprobing) { + // std::cout << "Activated by nRedundantIndices = " << + // dfprobingPropagation.redundantPropagateVec_.size() << std::endl; + dfprobingPropagation.propagate(); + if (!havePropagationRows() && + !dfprobingPropagation.isZeroObjFixingEnabled()) { + dfprobingPropagation.enableZeroObjFixing(); + dfprobingPropagation.setZeroCostFixingPosition(domchgstack_.size()); + dfprobingPropagation.propagate(); + } + } } return true; diff --git a/highs/mip/HighsDomain.h b/highs/mip/HighsDomain.h index 92747afd954..978ade04852 100644 --- a/highs/mip/HighsDomain.h +++ b/highs/mip/HighsDomain.h @@ -12,6 +12,7 @@ #include #include #include +#include #include #include "HighsPseudocost.h" @@ -235,6 +236,128 @@ class HighsDomain { void propagateConflict(HighsInt conflict); }; + struct DualfixingProbingPropagation { + HighsDomain* domain; + HighsMipSolver* mipsolver; + + // row lower and upper, length = 2 * rownum + std::vector redundantPropagateFlag_; + std::vector redundantPropagateVec_; + + // For zero-cost variables, we need to know which direction we can fix them + enum DFPROBING_FIX_DIRECTION { + FIXDIRECTION_NOT_DECIDED = 0, + FIXDIRECTION_LOWER_BOUND, + FIXDIRECTION_UPPER_BOUND, + }; + std::vector zeroCostVarsDirection_; + std::vector> zeroCostFixedVariables_; + + // Flag and position in the domchgstack of the first zero-cost variable that + // can be fixed to its lower or upper bound. + bool startZeroCostFixing_ = false; + size_t zeroCostStartPos_; + + bool enabled_ = false; + size_t previousSize_; + + // Original lower and upper locks, and the reduced locks after propagation. + std::vector colLowerLockOriginal_; + std::vector colUpperLockOriginal_; + std::vector colLowerLockReduced_; + std::vector colUpperLockReduced_; + + // temporary buffers for DFProbing + std::vector candidatesVec_; + std::vector candidatesFlag_; + std::unordered_set lockNeedClear_; + + // temporary buffers for GDF + std::vector gdfCandidatesVec_; + std::vector gdfCandidatesFlag_; + + // GDF reachable-row counts, indexed by column id. For each + // variable touched during probing, we only need to know how many + // rows make this variable lower/upper bound reachable. + std::vector gdfLbReachable0_; + std::vector gdfLbReachable1_; + std::vector gdfUbReachable0_; + std::vector gdfUbReachable1_; + + void enablePropagator() { enabled_ = true; } + + void disablePropagator() { enabled_ = false; } + + bool isEnabled() { return enabled_; } + + // active only when new redundant rows are found. + bool isActive() { + return enabled_ && redundantPropagateVec_.size() > previousSize_; + } + + // mark the position when the first zero-cost variable can be fixed to its + // lower or upper bound. + void setZeroCostFixingPosition(HighsInt v) { zeroCostStartPos_ = v; } + + size_t getZeroCostFixingPosition() { return zeroCostStartPos_; } + + void enableZeroObjFixing() { startZeroCostFixing_ = true; } + + void disableZeroObjFixing() { startZeroCostFixing_ = false; } + + bool isZeroObjFixingEnabled() { return startZeroCostFixing_; } + + bool ableToFixToLb(int col) { + return mipsolver->model_->col_cost_[col] >= + -mipsolver->options_mip_->dual_feasibility_tolerance && + mipsolver->model_->col_lower_[col] > -kHighsInf; + } + + bool ableToFixToUb(int col) { + return mipsolver->model_->col_cost_[col] <= + mipsolver->options_mip_->dual_feasibility_tolerance && + mipsolver->model_->col_upper_[col] < kHighsInf; + } + + // remove redundant information + void clearRedundantInfo() { + previousSize_ = 0; + if (!redundantPropagateVec_.empty()) { // clear buffers + for (const auto x : redundantPropagateVec_) + redundantPropagateFlag_[x] = false; + + redundantPropagateVec_.clear(); + } + + for (size_t i = 0; i < redundantPropagateFlag_.size(); ++i) + assert(!redundantPropagateFlag_[i]); + + zeroCostFixedVariables_.clear(); + + for (const auto x : lockNeedClear_) + colLowerLockReduced_[x] = colUpperLockReduced_[x] = 0; + lockNeedClear_.clear(); + } + + DualfixingProbingPropagation() { ; }; + + DualfixingProbingPropagation(HighsDomain* domain) : domain(domain) {}; + + DualfixingProbingPropagation(const DualfixingProbingPropagation& other); + + ~DualfixingProbingPropagation() { ; }; + + void recomputeLocks(); + void updateRhsRedundant(HighsInt row); + void updateLhsRedundant(HighsInt row); + void propagate(); + + // functionalities for GDF + void updateGDFInfo(HighsInt probing_variable, bool val); + HighsInt processGDFFixing(); + void clearGDFInfo(); + }; + private: struct ObjectivePropagation { HighsDomain* domain = nullptr; @@ -320,6 +443,7 @@ class HighsDomain { private: std::deque cutpoolpropagation; std::deque conflictPoolPropagation; + DualfixingProbingPropagation dfprobingPropagation; bool infeasible_ = false; Reason infeasible_reason; @@ -351,6 +475,8 @@ class HighsDomain { std::vector col_lower_; std::vector col_upper_; + bool inProbing_ = false; + HighsDomain(HighsMipSolver& mipsolver); HighsDomain(const HighsDomain& other) @@ -370,6 +496,7 @@ class HighsDomain { mipsolver(other.mipsolver), cutpoolpropagation(other.cutpoolpropagation), conflictPoolPropagation(other.conflictPoolPropagation), + dfprobingPropagation(other.dfprobingPropagation), infeasible_(other.infeasible_), infeasible_reason(other.infeasible_reason), infeasible_pos(other.infeasible_pos), @@ -383,6 +510,7 @@ class HighsDomain { for (ConflictPoolPropagation& conflictprop : conflictPoolPropagation) conflictprop.domain = this; if (objProp_.domain) objProp_.domain = this; + dfprobingPropagation.domain = this; } HighsDomain& operator=(const HighsDomain& other) { @@ -414,6 +542,7 @@ class HighsDomain { for (ConflictPoolPropagation& conflictprop : conflictPoolPropagation) conflictprop.domain = this; if (objProp_.domain) objProp_.domain = this; + dfprobingPropagation.domain = this; return *this; } @@ -686,6 +815,10 @@ class HighsDomain { void setRecordRedundantRows(bool val) { recordRedundantRows_ = val; }; bool isRedundantRow(HighsInt row) const; + + DualfixingProbingPropagation& getDfProbingPropagation() { + return dfprobingPropagation; + } }; #endif diff --git a/highs/mip/HighsImplications.cpp b/highs/mip/HighsImplications.cpp index b04cf9a9e1c..787c911c12e 100644 --- a/highs/mip/HighsImplications.cpp +++ b/highs/mip/HighsImplications.cpp @@ -27,6 +27,17 @@ bool HighsImplications::computeImplications(HighsInt col, bool val) { const auto& domchgreason = globaldomain.getDomainChangeReason(); size_t changedend = globaldomain.getChangedCols().size(); + // get two flags + const bool useDFProbing = + globaldomain.inProbing_ && mipsolver.options_mip_->presolve_dfprobing; + const bool useGDF = + globaldomain.inProbing_ && mipsolver.options_mip_->presolve_gdf; + // record redundant rows if any of the two flags is true + if (useDFProbing || useGDF) { + globaldomain.getDfProbingPropagation().clearRedundantInfo(); + globaldomain.getDfProbingPropagation().enablePropagator(); + } + HighsInt stackimplicstart = domchgstack.size() + 1; HighsInt numImplications = -stackimplicstart; if (val) @@ -53,6 +64,8 @@ bool HighsImplications::computeImplications(HighsInt col, bool val) { auto isInfeasible = [&](HighsInt col, bool val) { if (!globaldomain.infeasible()) return false; + if (globaldomain.inProbing_) + globaldomain.getDfProbingPropagation().disablePropagator(); storeLiftingOpportunities(col, val); doBacktrack(changedend); cliquetable.vertexInfeasible(globaldomain, col, val); @@ -62,6 +75,8 @@ bool HighsImplications::computeImplications(HighsInt col, bool val) { if (isInfeasible(col, val)) return true; globaldomain.propagate(); + if (useDFProbing || useGDF) + globaldomain.getDfProbingPropagation().disablePropagator(); if (isInfeasible(col, val)) return true; @@ -73,23 +88,54 @@ bool HighsImplications::computeImplications(HighsInt col, bool val) { std::vector implics; implics.reserve(numImplications); + // data structure to cache implications for non-binary variables + std::vector implics_tentative; + implics_tentative.reserve(numImplications); + HighsInt numEntries = mipsolver.mipdata_->cliquetable.getNumEntries(); HighsInt maxEntries = 100000 + mipsolver.numNonzero(); + const HighsInt tentativeStart = + useDFProbing + ? globaldomain.getDfProbingPropagation().getZeroCostFixingPosition() + : kHighsIInf32; + if (useDFProbing) + implics_tentative.assign(domchgstack.begin() + stackimplicstart, + domchgstack.begin() + stackimplicend); + for (HighsInt i = stackimplicstart; i < stackimplicend; ++i) { if (domchgreason[i].type == HighsDomain::Reason::kCliqueTable && ((domchgreason[i].index >> 1) == col || numEntries >= maxEntries)) continue; + if (i >= tentativeStart) // record tentative implications + continue; + implics.push_back(domchgstack[i]); } // inform caller about lifting opportunities storeLiftingOpportunities(col, val); + // update information to derive generalized dual fixings + if (useGDF) globaldomain.getDfProbingPropagation().updateGDFInfo(col, val); + // backtrack doBacktrack(changedend); + if (!implics_tentative.empty()) { + // add the implications of binary variables to the clique table + auto binstart_tmp = + std::partition(implics_tentative.begin(), implics_tentative.end(), + [&](const HighsDomainChange& a) { + return !globaldomain.isBinary(a.column); + }); + // store the tentative bound changes of binary variables separately + for (auto i = binstart_tmp; i != implics_tentative.end(); ++i) + recordTentativeCliques(val, *i); + implics_tentative.erase(binstart_tmp, implics_tentative.end()); + } + // add the implications of binary variables to the clique table auto binstart = std::partition(implics.begin(), implics.end(), [&](const HighsDomainChange& a) { @@ -147,6 +193,10 @@ bool HighsImplications::computeImplications(HighsInt col, bool val) { implications[loc].implics = std::move(implics); this->numImplications += implications[loc].implics.size(); } + if (!implics_tentative.empty()) { + pdqsort(implics_tentative.begin(), implics_tentative.end()); + implications[loc].implics_tentative = std::move(implics_tentative); + } return false; } @@ -300,6 +350,18 @@ bool HighsImplications::runProbing(HighsInt col, HighsInt& numReductions) { if (globaldomain.isBinary(col) && !implicationsCached(col, 1) && !implicationsCached(col, 0) && mipsolver.mipdata_->cliquetable.getSubstitution(col) == nullptr) { + const bool useDFProbing = + globaldomain.inProbing_ && mipsolver.options_mip_->presolve_dfprobing; + const bool useGDF = + globaldomain.inProbing_ && mipsolver.options_mip_->presolve_gdf; + // setup for dfprobingPropagation + if (useDFProbing) { + clearTentativeClique(); + globaldomain.getDfProbingPropagation().setZeroCostFixingPosition( + kHighsIInf32); + } + if (useGDF) globaldomain.getDfProbingPropagation().clearGDFInfo(); + bool infeasible = computeImplications(col, 1); if (globaldomain.infeasible()) return true; if (infeasible) return true; @@ -312,11 +374,98 @@ bool HighsImplications::runProbing(HighsInt col, HighsInt& numReductions) { if (mipsolver.mipdata_->cliquetable.getSubstitution(col) != nullptr) return true; + if (useDFProbing && !binaryInvolvedInds_.empty()) { + HighsCliqueTable& cliquetable = mipsolver.mipdata_->cliquetable; + HighsCliqueTable::CliqueVar clique[2]; + // Loop over binary variables that are tighened at least once + for (auto k : binaryInvolvedInds_) { + // Skip non-binary variables (being fixed now) or those can be + // substituted by other binary variables + if (!globaldomain.isBinary(k) || colsubstituted[k]) continue; + // Return if infeasible + if (globaldomain.infeasible()) return true; + // Get the information how x[k] is fixed in probing on x[col] = 0 and + // x[col] = 1 For the meaning of ``data'', please see lines 71-89 in + // HighsImplications.h + uint8_t data = binaryInvolvedFlags_[k]; + if (data == 0) // flag for no reduction + continue; + + if (data == + binaryFixType::kGlobalLower) { // x[k] is fixed at 0 under both + // x[col] = 0 and x[col] = 1 + // fix x[k] = 0 by adding two cliques (i.e., these two cliques should + // be added in computeImplications() to derive global reductions) + clique[0] = HighsCliqueTable::CliqueVar(col, 0); + clique[1] = HighsCliqueTable::CliqueVar(k, 1); + cliquetable.addClique(mipsolver, &clique[0], 2); + clique[0] = HighsCliqueTable::CliqueVar(col, 1); + clique[1] = HighsCliqueTable::CliqueVar(k, 1); + cliquetable.addClique(mipsolver, &clique[0], 2); + data = 0; + } else if (data == binaryFixType::kGlobalUpper) { // x[k] is fixed at 1 + // under both x[col] + // = 0 and x[col] = 1 + // fix x[k] = 1 by adding two cliques (i.e., these two cliques should + // be added in computeImplications() to derive global reductions) + clique[0] = HighsCliqueTable::CliqueVar(col, 0); + clique[1] = HighsCliqueTable::CliqueVar(k, 0); + cliquetable.addClique(mipsolver, &clique[0], 2); + clique[0] = HighsCliqueTable::CliqueVar(col, 1); + clique[1] = HighsCliqueTable::CliqueVar(k, 0); + cliquetable.addClique(mipsolver, &clique[0], 2); + data = 0; + } else if (data == + binaryFixType:: + kSubstituteComplement) { // x[k] is fixed at 0 under + // x[col] = 1, and is fixed at + // 1 under x[col] = 0; this + // makes x[col] + x[k] = 1 + // Adding two cliques (i.e., these two cliques should be added in + // computeImplications() to derive global reductions) + clique[0] = HighsCliqueTable::CliqueVar(col, 1); + clique[1] = HighsCliqueTable::CliqueVar(k, 1); + cliquetable.addClique(mipsolver, &clique[0], 2); + clique[0] = HighsCliqueTable::CliqueVar(col, 0); + clique[1] = HighsCliqueTable::CliqueVar(k, 0); + cliquetable.addClique(mipsolver, &clique[0], 2); + data = 0; + } else if (data == + binaryFixType:: + kSubstituteEqual) { // x[k] is fixed at 0 under x[col] = + // 0, and is fixed at 1 under x[col] + // = 1; this makes x[col] = x[k] + // Adding two cliques (i.e., these two cliques should be added in + // computeImplications() to derive global reductions) + clique[0] = HighsCliqueTable::CliqueVar(col, 1); + clique[1] = HighsCliqueTable::CliqueVar(k, 0); + cliquetable.addClique(mipsolver, &clique[0], 2); + clique[0] = HighsCliqueTable::CliqueVar(col, 0); + clique[1] = HighsCliqueTable::CliqueVar(k, 1); + cliquetable.addClique(mipsolver, &clique[0], 2); + data = 0; + } + } + + // clear the tentative bound changes for binary variables obtained from + // probing on x[col] + clearTentativeClique(); + } + // analyze implications + // also include the bound changes of non-binary variables here, to derive + // tighter global bounds and variable substitutions + const bool haveTentativeImplics_zero = + !implications[2 * col].implics_tentative.empty(); + const bool haveTentativeImplics_one = + !implications[2 * col + 1].implics_tentative.empty(); + const std::vector& implicsdown = - getImplications(col, 0, infeasible); + haveTentativeImplics_zero ? getImplications_tentative(col, 0) + : getImplications(col, 0, infeasible); const std::vector& implicsup = - getImplications(col, 1, infeasible); + haveTentativeImplics_one ? getImplications_tentative(col, 1) + : getImplications(col, 1, infeasible); HighsInt nimplicsdown = implicsdown.size(); HighsInt nimplicsup = implicsup.size(); HighsInt u = 0; @@ -382,6 +531,19 @@ bool HighsImplications::runProbing(HighsInt col, HighsInt& numReductions) { } } + // clear tentative implications + if (haveTentativeImplics_zero) + implications[2 * col].implics_tentative.clear(); + if (haveTentativeImplics_one) + implications[2 * col + 1].implics_tentative.clear(); + + if (useGDF) { + // fix variables using generalized dual fixing + HighsInt nfix = globaldomain.getDfProbingPropagation().processGDFFixing(); + // propagate if necessary + if (nfix > 0) globaldomain.propagate(); + } + return true; } diff --git a/highs/mip/HighsImplications.h b/highs/mip/HighsImplications.h index 50fa737395e..9cdaf6b2bb3 100644 --- a/highs/mip/HighsImplications.h +++ b/highs/mip/HighsImplications.h @@ -25,6 +25,13 @@ class HighsImplications { struct Implics { std::vector implics; + /* The "tentative" implications: + A implication of type x_j \ge (\ell^1_j - \ell^0_j) x_k + \ell^0_j is + called "tentative", if (1) c_j = 0 (2) x_j is fixed by applying dual + fixing in probing These implications can only be used to perform globally + valid reductions. Therefore, special treatment is required. + */ + std::vector implics_tentative; bool computed = false; }; std::vector implications; @@ -57,6 +64,30 @@ class HighsImplications { const HighsMipSolver& mipsolver; std::vector substitutions; std::vector colsubstituted; + + // vector used to derive global reductions from dfprobing + std::vector binaryInvolvedInds_; + enum binaryFixType { + kNoReduction = 0b0000, + kGlobalLower = 0b1010, + kGlobalUpper = 0b0101, + kSubstituteComplement = 0b1001, + kSubstituteEqual = 0b0110, + }; + /* + Possible values for binaryInvolvedFlags_ + 0 (0000, kNoReduction): Not involved + 2 (0010): fixed to 0 in second side probing + 1 (0001): fixed to 1 in second side probing + 8 (1000): fixed to 0 in first side probing + 4 (0100): fixed to 1 in first side probing + 10(1010, kGlobalLower): fixed to 0 in both side probing (global fixing!) + 5 (0101, kGlobalUpper): fixed to 1 in both side probing (global fixing!) + 9 (1001, kSubstituteComplement): substitutation type 1 --- x1 + x2 = 1 + 6 (0110, kSubstituteEqual): substitutation type 2 --- x1 = x2 + */ + std::vector binaryInvolvedFlags_; + HighsImplications(const HighsMipSolver& mipsolver) : mipsolver(mipsolver) { HighsInt numcol = mipsolver.numCol(); implications.resize(2 * static_cast(numcol)); @@ -67,6 +98,9 @@ class HighsImplications { numImplications = 0; numVarBounds = 0; maxVarBounds = calcMaxVarBounds(numcol); + + binaryInvolvedInds_.reserve(numcol); + binaryInvolvedFlags_.assign(numcol, 0b0000); } std::function @@ -92,6 +126,8 @@ class HighsImplications { maxVarBounds = calcMaxVarBounds(numcol); nextCleanupCall = mipsolver.numNonzero(); + binaryInvolvedInds_.reserve(numcol); + binaryInvolvedFlags_.assign(numcol, 0b0000); } constexpr static int64_t calcMaxVarBounds(HighsInt numcol) { @@ -115,6 +151,12 @@ class HighsImplications { return implications[loc].implics; } + const std::vector& getImplications_tentative(HighsInt col, + bool val) { + HighsInt loc = 2 * col + val; + return implications[loc].implics_tentative; + } + bool implicationsCached(HighsInt col, bool val) { HighsInt loc = 2 * col + val; return implications[loc].computed; @@ -192,6 +234,72 @@ class HighsImplications { bool& infeasible, bool allowBoundChanges = true) const; void applyImplications(HighsDomain& domain, HighsInt col, HighsInt val); + + // collect tentative binary implications + void recordTentativeCliques(bool val, const HighsDomainChange& bchg) { + const int iCol = bchg.column; + if (val == 0) { // probing x_k = 0 + if (bchg.boundtype == HighsBoundType::kLower) { // fixed to 1 + if (!isFixedTo1(val, iCol)) { + if (binaryInvolvedFlags_[iCol] == 0) + binaryInvolvedInds_.push_back(iCol); + binaryInvolvedFlags_[iCol] += 0b0001; // 0001 + } + } else { // fixed to 0 + if (!isFixedTo0(val, iCol)) { + if (binaryInvolvedFlags_[iCol] == 0) + binaryInvolvedInds_.push_back(iCol); + binaryInvolvedFlags_[iCol] += 0b0010; // 0010 + } + } + } else { // probing x_k = 1 + if (bchg.boundtype == HighsBoundType::kLower) { // fixed to 1 + if (!isFixedTo1(val, iCol)) { + if (binaryInvolvedFlags_[iCol] == 0) + binaryInvolvedInds_.push_back(iCol); + binaryInvolvedFlags_[iCol] += 0b0100; // 0100 + } + } else { // fixed to 0 + if (!isFixedTo0(val, iCol)) { + if (binaryInvolvedFlags_[iCol] == 0) + binaryInvolvedInds_.push_back(iCol); + binaryInvolvedFlags_[iCol] += 0b1000; // 1000 + } + } + } + } + // clear tentative binary implications + void clearTentativeClique() { + for (auto iCol : binaryInvolvedInds_) + binaryInvolvedFlags_[iCol] = binaryFixType::kNoReduction; + binaryInvolvedInds_.clear(); + } + // tools for recordTentativeCliques + bool isFixedTo0(bool val, HighsInt iCol) { + if (binaryInvolvedFlags_[iCol] == 0) return false; + + uint8_t mask; + if (val == 0) { // probing at x = 0, last two digits + mask = 1 << (1); + return (binaryInvolvedFlags_[iCol] & mask) != 0; + } else { // probing at x = 1, first two digits + mask = 1 << (3); + return (binaryInvolvedFlags_[iCol] & mask) != 0; + } + } + // tools for recordTentativeCliques + bool isFixedTo1(bool val, HighsInt iCol) { + if (binaryInvolvedFlags_[iCol] == 0) return false; + + uint8_t mask; + if (val == 0) { // probing at x = 0, last two digits + mask = 1; + return (binaryInvolvedFlags_[iCol] & mask) != 0; + } else { // probint at x = 1, first two digits + mask = 1 << (2); + return (binaryInvolvedFlags_[iCol] & mask) != 0; + } + } }; #endif diff --git a/highs/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp index 90a3e0c8bab..93201952049 100644 --- a/highs/presolve/HPresolve.cpp +++ b/highs/presolve/HPresolve.cpp @@ -1863,6 +1863,10 @@ HPresolve::Result HPresolve::runProbing(HighsPostsolveStack& postsolve_stack) { } }; + // setup for dfprobing and gdf + if (options->presolve_dfprobing || options->presolve_gdf) + domain.getDfProbingPropagation().recomputeLocks(); + for (const auto& binvar : binaries) { // Count the binaries considered iBin++; @@ -1928,7 +1932,9 @@ HPresolve::Result HPresolve::runProbing(HighsPostsolveStack& postsolve_stack) { HighsInt numBoundChgs = 0; HighsInt numNewCliques = -cliquetable.numCliques(); + domain.inProbing_ = true; const bool probing_result = implications.runProbing(i, numBoundChgs); + domain.inProbing_ = false; if (!probing_result) continue; probingContingent += numBoundChgs; numNewCliques += cliquetable.numCliques();