Skip to content

Commit 87135b6

Browse files
Opt-MuccaOpt-Mucca
authored
Merge pull request #2869 from fwesselm/dualFixingExtension3
Extension of dual fixing for single equations
2 parents b174cc9 + 35c5bf6 commit 87135b6

File tree

2 files changed

+175
-46
lines changed

2 files changed

+175
-46
lines changed

‎highs/presolve/HPresolve.cpp‎

Lines changed: 172 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -4588,30 +4588,29 @@ HPresolve::Result HPresolve::detectDominatedCol(
45884588

45894589
void HPresolve::computeLocks(
45904590
HighsInt col, bool considerObjective,
4591-
std::function<bool(HighsInt, bool)> lockCallback) const {
4591+
std::function<bool(HighsInt, bool, bool)> lockCallback) const {
45924592
// if the callback returns true, we stop examining the locks
45934593
if (considerObjective) {
45944594
// consider objective function
45954595
if (model->col_cost_[col] < 0) {
45964596
// downlock
4597-
if (lockCallback(-1, false)) return;
4597+
if (lockCallback(-1, true, false)) return;
45984598
} else if (model->col_cost_[col] > 0) {
45994599
// uplock
4600-
if (lockCallback(-1, true)) return;
4600+
if (lockCallback(-1, false, true)) return;
46014601
}
46024602
}
46034603

46044604
// check coefficients
46054605
for (const auto& nz : getColumnVector(col)) {
46064606
// implied lower bound -> downlock
4607-
if (yieldsImpliedLowerBound(nz.index(), nz.value()) &&
4608-
lockCallback(nz.index(), false))
4609-
break;
4610-
4607+
bool hasDownLock = yieldsImpliedLowerBound(nz.index(), nz.value());
46114608
// implied upper bound -> uplock
4612-
if (yieldsImpliedUpperBound(nz.index(), nz.value()) &&
4613-
lockCallback(nz.index(), true))
4614-
break;
4609+
bool hasUpLock = yieldsImpliedUpperBound(nz.index(), nz.value());
4610+
// callback
4611+
if (hasDownLock || hasUpLock) {
4612+
if (lockCallback(nz.index(), hasDownLock, hasUpLock)) break;
4613+
}
46154614
}
46164615
}
46174616

@@ -4685,6 +4684,124 @@ HPresolve::Result HPresolve::dualFixing(HighsPostsolveStack& postsolve_stack,
46854684
return Result::kOk;
46864685
};
46874686

4687+
// lambda that checks column for single equation handling
4688+
auto checkColumn = [&](HighsInt col, HighsInt colDirection, HighsInt row) {
4689+
if (colDirection * model->col_cost_[col] < 0) return false;
4690+
for (const auto& colNz : getColumnVector(col)) {
4691+
// skip equation we are inspecting; coefficient is positive (see sign
4692+
// adjustment performed by caller)
4693+
if (colNz.index() == row) continue;
4694+
// skip redundant rows
4695+
if (isRedundant(colNz.index())) continue;
4696+
// another ranged row blocks variable in both directions -> cannot be a
4697+
// candidate for fixing
4698+
if (isRanged(colNz.index())) return false;
4699+
// compute row direction
4700+
HighsInt rowDirection =
4701+
model->row_lower_[colNz.index()] == -kHighsInf &&
4702+
model->row_upper_[colNz.index()] != kHighsInf
4703+
? 1
4704+
: -1;
4705+
// coefficient must have the correct sign
4706+
if (colDirection * rowDirection * colNz.value() < 0) return false;
4707+
}
4708+
return true;
4709+
};
4710+
4711+
// lambda that computes activities for single equation handling
4712+
auto computeActivity = [&](HighsInt col, double val, HighsCDouble& activity,
4713+
bool& activityFinite, HighsInt direction) {
4714+
// direction = 1 -> update upper bound on activity (supremum)
4715+
// direction = -1 -> update lower bound on activity (infimum)
4716+
double bound =
4717+
direction * val > 0 ? model->col_upper_[col] : model->col_lower_[col];
4718+
activityFinite = activityFinite && std::abs(bound) != kHighsInf;
4719+
if (activityFinite) activity += static_cast<HighsCDouble>(val) * bound;
4720+
};
4721+
4722+
// lambda for fixing variables
4723+
auto fixCols = [&](const std::vector<std::pair<HighsInt, double>>& vars,
4724+
HighsInt direction) {
4725+
for (const auto& elm : vars) {
4726+
if (direction * elm.second > 0)
4727+
HPRESOLVE_CHECKED_CALL(fixColToUpper(postsolve_stack, elm.first));
4728+
else
4729+
HPRESOLVE_CHECKED_CALL(fixColToLower(postsolve_stack, elm.first));
4730+
}
4731+
return Result::kOk;
4732+
};
4733+
4734+
// lambda for handling single equations
4735+
auto handleSingleEquation = [&](HighsInt row) {
4736+
assert(isEquation(row));
4737+
std::vector<std::pair<HighsInt, double>> sPlus;
4738+
std::vector<std::pair<HighsInt, double>> sMinus;
4739+
sPlus.reserve(rowsize[row]);
4740+
sMinus.reserve(rowsize[row]);
4741+
std::vector<HighsInt> sMark(model->num_col_, 0);
4742+
for (const auto& rowNz : getRowVector(row)) {
4743+
// flip column direction if coefficient is negative
4744+
// (equivalent to negating the column's coefficients)
4745+
HighsInt colDirection = std::copysign(HighsInt{1}, rowNz.value());
4746+
if (checkColumn(rowNz.index(), colDirection, row)) {
4747+
// store in S_+
4748+
sPlus.push_back(std::make_pair(rowNz.index(), rowNz.value()));
4749+
sMark[rowNz.index()] = 1;
4750+
} else if (checkColumn(rowNz.index(), -colDirection, row)) {
4751+
// store in S_-
4752+
sMinus.push_back(std::make_pair(rowNz.index(), rowNz.value()));
4753+
sMark[rowNz.index()] = -1;
4754+
}
4755+
}
4756+
// return if both sets are empty
4757+
if (sPlus.empty() && sMinus.empty()) return Result::kOk;
4758+
// compute activities
4759+
HighsCDouble activityTPlus = 0.0;
4760+
HighsCDouble activityTMinus = 0.0;
4761+
HighsCDouble activitySCPlus = 0.0;
4762+
HighsCDouble activitySCMinus = 0.0;
4763+
bool activityTPlusFinite = true;
4764+
bool activityTMinusFinite = true;
4765+
bool activitySCPlusFinite = true;
4766+
bool activitySCMinusFinite = true;
4767+
for (const auto& rowNz : getRowVector(row)) {
4768+
if (sMark[rowNz.index()] <= 0 ||
4769+
model->integrality_[rowNz.index()] != HighsVarType::kContinuous)
4770+
// T_+
4771+
computeActivity(rowNz.index(), rowNz.value(), activityTPlus,
4772+
activityTPlusFinite, HighsInt{1});
4773+
else
4774+
// S^C_+
4775+
computeActivity(rowNz.index(), rowNz.value(), activitySCPlus,
4776+
activitySCPlusFinite, HighsInt{-1});
4777+
if (sMark[rowNz.index()] >= 0 ||
4778+
model->integrality_[rowNz.index()] != HighsVarType::kContinuous)
4779+
// T_-
4780+
computeActivity(rowNz.index(), rowNz.value(), activityTMinus,
4781+
activityTMinusFinite, HighsInt{-1});
4782+
else
4783+
// S^C_-
4784+
computeActivity(rowNz.index(), rowNz.value(), activitySCMinus,
4785+
activitySCMinusFinite, HighsInt{1});
4786+
// break if activities are not finite
4787+
if ((sMinus.empty() || !activityTPlusFinite || !activitySCPlusFinite) &&
4788+
(sPlus.empty() || !activityTMinusFinite || !activitySCMinusFinite))
4789+
break;
4790+
}
4791+
// fix variables
4792+
if (!sMinus.empty() && activityTPlusFinite && activitySCPlusFinite &&
4793+
activityTPlus + activitySCPlus <=
4794+
model->row_lower_[row] + primal_feastol)
4795+
// fix all variables in S_-
4796+
HPRESOLVE_CHECKED_CALL(fixCols(sMinus, HighsInt{1}));
4797+
else if (!sPlus.empty() && activityTMinusFinite && activitySCMinusFinite &&
4798+
activityTMinus + activitySCMinus >=
4799+
model->row_lower_[row] - primal_feastol)
4800+
// fix all variables in S_+
4801+
HPRESOLVE_CHECKED_CALL(fixCols(sPlus, HighsInt{-1}));
4802+
return Result::kOk;
4803+
};
4804+
46884805
// lambda for computing tighter bounds
46894806
auto hasTighterBound = [&](HighsInt col, HighsInt direction,
46904807
double currentBound, double& newBound) {
@@ -4719,12 +4836,13 @@ HPresolve::Result HPresolve::dualFixing(HighsPostsolveStack& postsolve_stack,
47194836
HighsInt numUpLocks = 0;
47204837
HighsInt downLockRow = -1;
47214838
HighsInt upLockRow = -1;
4722-
auto lockCallback = [&](HighsInt row, bool isUpLock) {
4839+
auto lockCallback = [&](HighsInt row, bool hasDownLock, bool hasUpLock) {
47234840
// count locks and remember row index
4724-
if (isUpLock) {
4841+
if (hasUpLock) {
47254842
numUpLocks++;
47264843
upLockRow = row;
4727-
} else {
4844+
}
4845+
if (hasDownLock) {
47284846
numDownLocks++;
47294847
downLockRow = row;
47304848
}
@@ -4743,20 +4861,34 @@ HPresolve::Result HPresolve::dualFixing(HighsPostsolveStack& postsolve_stack,
47434861
else
47444862
HPRESOLVE_CHECKED_CALL(fixColToUpper(postsolve_stack, col));
47454863
} else {
4746-
if (mipsolver != nullptr && model->col_lower_[col] != -kHighsInf &&
4747-
model->col_upper_[col] != kHighsInf) {
4748-
// try substitution
4749-
if (numDownLocks == 1 && downLockRow != -1) {
4750-
HPRESOLVE_CHECKED_CALL(substituteCol(col, downLockRow, HighsInt{1},
4751-
model->col_upper_[col],
4752-
model->col_lower_[col]));
4753-
if (colDeleted[col]) return Result::kOk;
4754-
}
4755-
if (numUpLocks == 1 && upLockRow != -1) {
4756-
HPRESOLVE_CHECKED_CALL(substituteCol(col, upLockRow, HighsInt{-1},
4757-
model->col_lower_[col],
4758-
model->col_upper_[col]));
4864+
bool hasSingleDownLock = numDownLocks == 1 && downLockRow != -1;
4865+
bool hasSingleUpLock = numUpLocks == 1 && upLockRow != -1;
4866+
if (hasSingleDownLock || hasSingleUpLock) {
4867+
HighsInt equationRow =
4868+
hasSingleDownLock && isEquation(downLockRow)
4869+
? downLockRow
4870+
: (hasSingleUpLock && isEquation(upLockRow) ? upLockRow : -1);
4871+
if (equationRow != -1) {
4872+
// see section 6.1 "Extension of dual fixing for single equations",
4873+
// Achterberg et al., Presolve Reductions in Mixed Integer
4874+
// Programming, INFORMS Journal on Computing 32(2):473-506.
4875+
HPRESOLVE_CHECKED_CALL(handleSingleEquation(equationRow));
47594876
if (colDeleted[col]) return Result::kOk;
4877+
} else if (mipsolver != nullptr && model->col_lower_[col] != -kHighsInf &&
4878+
model->col_upper_[col] != kHighsInf) {
4879+
// try substitution
4880+
if (hasSingleDownLock) {
4881+
HPRESOLVE_CHECKED_CALL(substituteCol(col, downLockRow, HighsInt{1},
4882+
model->col_upper_[col],
4883+
model->col_lower_[col]));
4884+
if (colDeleted[col]) return Result::kOk;
4885+
}
4886+
if (hasSingleUpLock) {
4887+
HPRESOLVE_CHECKED_CALL(substituteCol(col, upLockRow, HighsInt{-1},
4888+
model->col_lower_[col],
4889+
model->col_upper_[col]));
4890+
if (colDeleted[col]) return Result::kOk;
4891+
}
47604892
}
47614893
}
47624894
// try to strengthen bounds
@@ -5073,10 +5205,11 @@ HPresolve::Result HPresolve::enumerateSolutions(
50735205
if ((std::get<4>(rows[i]) & std::get<4>(rows[ii])) == 0) continue;
50745206
// get indices of binary variables in the row
50755207
getBinaryRow(r2, binvars2, numnzs2);
5076-
// check if there is too much overlap
5208+
// check if there is too much overlap (compute overlap coefficient)
50775209
numComparisons++;
50785210
size_t overlap = computeRowOverlap(binvars, binvars2, numnzs, numnzs2);
5079-
if ((200 * overlap) / (numnzs + numnzs2) > maxPercentageRowOverlap) {
5211+
if ((100 * overlap) / std::min(numnzs, numnzs2) >
5212+
maxPercentageRowOverlap) {
50805213
// mark row for removal
50815214
numRowsRemoved++;
50825215
r2 = -1;
@@ -5121,17 +5254,19 @@ HPresolve::Result HPresolve::enumerateSolutions(
51215254
std::vector<double> col_lower(domain.col_lower_);
51225255
std::vector<double> col_upper(domain.col_upper_);
51235256

5257+
// lambda for finding a variable to branch on
5258+
auto findBranchVar = [&](size_t numVars) {
5259+
// find variable for branching
5260+
for (size_t i = 0; i < numVars; i++)
5261+
if (!domain.isFixed(vars[i])) return vars[i];
5262+
return HighsInt{-1};
5263+
};
5264+
51245265
// lambda for branching (just performs initial lower branch)
51255266
auto doBranch = [&](size_t numVars, HighsInt& numBranches) {
51265267
// find variable for branching
5127-
HighsInt branchvar = -1;
5128-
for (size_t i = 0; i < numVars; i++) {
5129-
if (!domain.isFixed(vars[i])) {
5130-
branchvar = vars[i];
5131-
break;
5132-
}
5133-
}
5134-
if (branchvar < 0) return;
5268+
HighsInt branchvar = findBranchVar(numVars);
5269+
assert(branchvar >= 0);
51355270

51365271
// branch downwards
51375272
branches[++numBranches] = {domain.getDomainChangeStack().size(),
@@ -5164,13 +5299,6 @@ HPresolve::Result HPresolve::enumerateSolutions(
51645299
return (numBranches >= 0);
51655300
};
51665301

5167-
// lambda for checking if a solution was found
5168-
auto solutionFound = [&](size_t numVars) {
5169-
for (size_t i = 0; i < numVars; i++)
5170-
if (!domain.isFixed(vars[i])) return false;
5171-
return true;
5172-
};
5173-
51745302
// lambda for checking whether the values of two binary variables are
51755303
// identical in all feasible solutions
51765304
auto identicalVars = [&](size_t numSolutions, size_t index1, size_t index2) {
@@ -5281,7 +5409,7 @@ HPresolve::Result HPresolve::enumerateSolutions(
52815409
while (true) {
52825410
bool backtrack = domain.infeasible();
52835411
if (!backtrack) {
5284-
backtrack = solutionFound(numVars);
5412+
backtrack = findBranchVar(numVars) < 0;
52855413
if (backtrack)
52865414
handleSolution(numVars, numSolutions, numWorstCaseBounds,
52875415
minNumActiveCols, maxNumActiveCols);

‎highs/presolve/HPresolve.h‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -404,8 +404,9 @@ class HPresolve {
404404
Result detectDominatedCol(HighsPostsolveStack& postsolve_stack, HighsInt col,
405405
bool handleSingletonRows = true);
406406

407-
void computeLocks(HighsInt col, bool considerObjective,
408-
std::function<bool(HighsInt, bool)> lockCallback) const;
407+
void computeLocks(
408+
HighsInt col, bool considerObjective,
409+
std::function<bool(HighsInt, bool, bool)> lockCallback) const;
409410

410411
Result dualFixing(HighsPostsolveStack& postsolve_stack, HighsInt col);
411412

0 commit comments

Comments
 (0)