Skip to content

Commit 3a03e04

Browse files
Opt-MuccaOpt-Mucca
committed
Refactor. Add comments
1 parent c9a7977 commit 3a03e04

File tree

3 files changed

+91
-58
lines changed

3 files changed

+91
-58
lines changed

‎highs/mip/HighsMipSolverData.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -49,7 +49,7 @@ std::string HighsMipSolverData::solutionSourceToString(
4949
if (code) return "L";
5050
return "Sub-MIP";
5151
} else if (solution_source == kSolutionSourceLocalMip) {
52-
if (code) return "M";
52+
if (code) return "O";
5353
return "Local MIP";
5454
} else if (solution_source == kSolutionSourceEmptyMip) {
5555
if (code) return "P";

‎highs/mip/HighsMipSolverData.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -50,7 +50,7 @@ enum MipSolutionSource : int {
5050
kSolutionSourceShifting, // I
5151
kSolutionSourceFeasibilityJump, // J
5252
kSolutionSourceSubMip, // L
53-
kSolutionSourceLocalMip, // M
53+
kSolutionSourceLocalMip, // O
5454
kSolutionSourceEmptyMip, // P
5555
kSolutionSourceRandomizedRounding, // R
5656
kSolutionSourceSolveLp, // S

‎highs/mip/HighsPrimalHeuristics.cpp‎

Lines changed: 89 additions & 56 deletions
Original file line numberDiff line numberDiff line change
@@ -1599,7 +1599,7 @@ bool HighsPrimalHeuristics::localMip() {
15991599
HighsSparseMatrix a_matrix_row;
16001600
a_matrix_row.createRowwise(a_matrix);
16011601

1602-
// Create global objects used in lambda expressions
1602+
// Initialise values and reserve vectors
16031603
HighsInt num_violated_rows_sample = 250;
16041604
HighsInt num_satisfied_rows_sample = 20;
16051605
HighsInt iters = 0;
@@ -1662,20 +1662,22 @@ bool HighsPrimalHeuristics::localMip() {
16621662
return false;
16631663
};
16641664

1665-
struct Violated {
1665+
struct RowStatus {
16661666
std::vector<HighsInt> viol;
16671667
std::vector<HighsInt> viol_index;
16681668

1669-
explicit Violated(const HighsInt n) : viol(n, -1) { viol_index.reserve(n); }
1669+
explicit RowStatus(const HighsInt n) : viol(n, -1) {
1670+
viol_index.reserve(n);
1671+
}
16701672
};
16711673

1672-
auto add_entry = [&](Violated& v, HighsInt r) {
1674+
auto add_entry = [&](RowStatus& v, HighsInt r) {
16731675
assert(v.viol[r] == -1);
16741676
v.viol[r] = static_cast<HighsInt>(v.viol_index.size());
16751677
v.viol_index.push_back(r);
16761678
};
16771679

1678-
auto remove_entry = [&](Violated& v, HighsInt r) {
1680+
auto remove_entry = [&](RowStatus& v, HighsInt r) {
16791681
assert(v.viol[r] != -1);
16801682
v.viol_index[v.viol[r]] = v.viol_index.back();
16811683
v.viol[v.viol_index.back()] = v.viol[r];
@@ -1688,11 +1690,11 @@ bool HighsPrimalHeuristics::localMip() {
16881690
HighsInt momentum = 1;
16891691
std::vector<HighsInt> allow_neg_delta(mipsolver.numCol());
16901692
std::vector<HighsInt> allow_pos_delta(mipsolver.numCol());
1691-
Violated viol(mipsolver.numRow());
1692-
Violated satisfied(mipsolver.numRow());
1693+
RowStatus viol(mipsolver.numRow());
1694+
RowStatus satisfied(mipsolver.numRow());
16931695
std::vector<HighsCDouble> activities(mipsolver.numRow());
16941696

1695-
auto calc_activites = [&](bool reset_violated) {
1697+
auto calc_activities = [&](bool reset_violated) {
16961698
std::fill(activities.begin(), activities.end(), 0.0);
16971699
if (reset_violated) {
16981700
viol.viol_index.clear();
@@ -1751,8 +1753,10 @@ bool HighsPrimalHeuristics::localMip() {
17511753
// printf("Applying move col %d with delta %g\n", c, delta);
17521754
};
17531755

1754-
auto one_opt_calc_delta = [&](HighsInt c) -> double {
1755-
// Only check lower bound delta if positive column cost, likewise for neg.
1756+
// Calculate how much a column can be shifted while maintaining feasibility
1757+
// of all rows containing given column.
1758+
auto one_opt_calc_col_delta = [&](HighsInt c) -> double {
1759+
// Only check lower bound delta if positive column cost and vice versa.
17561760
const HighsInt dir = mipsolver.colCost(c) > 0 ? -1 : 1;
17571761
bool integral = mipsolver.isColIntegral(c);
17581762
double delta = dir == -1 ? sol[c] - globaldom.col_lower_[c]
@@ -1785,11 +1789,14 @@ bool HighsPrimalHeuristics::localMip() {
17851789
};
17861790

17871791
std::vector<double> one_opt_deltas(mipsolver.numCol());
1792+
1793+
// Find single column shift with the best objective change that maintains
1794+
// feasibility of all rows containing the column
17881795
auto one_opt = [&](bool recalc) -> bool {
17891796
std::vector<double> one_opt_scores(mipsolver.numCol());
17901797
if (recalc) {
17911798
for (const HighsInt c : cols_with_obj) {
1792-
one_opt_deltas[c] = one_opt_calc_delta(c);
1799+
one_opt_deltas[c] = one_opt_calc_col_delta(c);
17931800
}
17941801
}
17951802
HighsInt best_col = -1;
@@ -1803,7 +1810,7 @@ bool HighsPrimalHeuristics::localMip() {
18031810
}
18041811
if (best_col == -1) return false;
18051812
apply_move(best_col, one_opt_deltas[best_col]);
1806-
one_opt_deltas[best_col] = one_opt_calc_delta(best_col);
1813+
one_opt_deltas[best_col] = one_opt_calc_col_delta(best_col);
18071814
HighsInt start = a_matrix.start_[best_col];
18081815
HighsInt end = a_matrix.start_[best_col + 1];
18091816
for (HighsInt i = start; i < end; ++i) {
@@ -1813,50 +1820,56 @@ bool HighsPrimalHeuristics::localMip() {
18131820
for (HighsInt j = rstart; j < rend; ++j) {
18141821
HighsInt c = a_matrix_row.index_[j];
18151822
if (mipsolver.colCost(c) != 0) {
1816-
one_opt_deltas[c] = one_opt_calc_delta(c);
1823+
one_opt_deltas[c] = one_opt_calc_col_delta(c);
18171824
}
18181825
}
18191826
}
18201827
return true;
18211828
};
18221829

1823-
auto feas_one_opt = [&](HighsInt c, HighsInt r, double coef,
1824-
std::vector<std::pair<HighsInt, double>>& moves) {
1825-
// Get direction column can be moved to maximise constraint feasibility
1826-
double row_lower = mipsolver.model_->row_lower_[r];
1827-
double row_upper = mipsolver.model_->row_upper_[r];
1828-
bool lower_inf = row_lower == -kHighsInf;
1829-
bool upper_inf = row_upper == kHighsInf;
1830-
double slack_lower = static_cast<double>(activities[r]) - row_lower;
1831-
double slack_upper = row_upper - static_cast<double>(activities[r]);
1832-
HighsInt increase_activity;
1833-
if (!upper_inf && activities[r] > row_upper + feastol) {
1834-
increase_activity = -1;
1835-
} else if (!lower_inf && activities[r] < row_lower - feastol) {
1836-
increase_activity = 1;
1837-
} else {
1838-
return;
1839-
}
1840-
double delta = increase_activity == 1 ? std::abs(slack_lower / coef)
1841-
: std::abs(slack_upper / coef);
1842-
HighsInt dir = coef < 0 ? -increase_activity : increase_activity;
1843-
if (mipsolver.isColIntegral(c)) {
1844-
delta = std::ceil(delta - feastol);
1845-
}
1846-
delta = dir == -1
1830+
// Find shift of given column that respects column bounds and maximises
1831+
// feasibility of given row.
1832+
auto calc_col_max_row_feas_delta =
1833+
[&](HighsInt c, HighsInt r, double coef,
1834+
std::vector<std::pair<HighsInt, double>>& moves) {
1835+
// Get direction column can be moved to maximise row feasibility
1836+
double row_lower = mipsolver.model_->row_lower_[r];
1837+
double row_upper = mipsolver.model_->row_upper_[r];
1838+
bool lower_inf = row_lower == -kHighsInf;
1839+
bool upper_inf = row_upper == kHighsInf;
1840+
double slack_lower = static_cast<double>(activities[r]) - row_lower;
1841+
double slack_upper = row_upper - static_cast<double>(activities[r]);
1842+
HighsInt increase_activity;
1843+
if (!upper_inf && activities[r] > row_upper + feastol) {
1844+
increase_activity = -1;
1845+
} else if (!lower_inf && activities[r] < row_lower - feastol) {
1846+
increase_activity = 1;
1847+
} else {
1848+
return;
1849+
}
1850+
double delta = increase_activity == 1 ? std::abs(slack_lower / coef)
1851+
: std::abs(slack_upper / coef);
1852+
HighsInt dir = coef < 0 ? -increase_activity : increase_activity;
1853+
if (mipsolver.isColIntegral(c)) {
1854+
delta = std::ceil(delta - feastol);
1855+
}
1856+
delta =
1857+
dir == -1
18471858
? std::min(std::abs(sol[c] - globaldom.col_lower_[c]), delta)
18481859
: std::min(std::abs(globaldom.col_upper_[c] - sol[c]), delta);
1849-
if (mipsolver.isColIntegral(c)) {
1850-
delta = std::floor(delta + feastol);
1851-
}
1852-
if (delta < feastol) return;
1853-
if (!tabu(c, dir * delta)) return;
1854-
// TODO: Need to add some safety net for unbounded variables
1855-
assert(sol[c] + delta * dir > globaldom.col_lower_[c] - feastol);
1856-
assert(sol[c] + delta * dir < globaldom.col_upper_[c] + feastol);
1857-
moves.emplace_back(c, dir * delta);
1858-
};
1860+
// Don't try and shift unbounded columns
1861+
if (delta == kHighsInf) return;
1862+
if (mipsolver.isColIntegral(c)) {
1863+
delta = std::floor(delta + feastol);
1864+
}
1865+
if (delta < feastol) return;
1866+
if (!tabu(c, dir * delta)) return;
1867+
assert(sol[c] + delta * dir > globaldom.col_lower_[c] - feastol);
1868+
assert(sol[c] + delta * dir < globaldom.col_upper_[c] + feastol);
1869+
moves.emplace_back(c, dir * delta);
1870+
};
18591871

1872+
// Sample row indices
18601873
auto sample_indices = [&](std::vector<HighsInt>& indices,
18611874
HighsInt max_samples,
18621875
HighsInt& n) -> std::vector<HighsInt>& {
@@ -1889,6 +1902,7 @@ bool HighsPrimalHeuristics::localMip() {
18891902
return sampled_indices;
18901903
};
18911904

1905+
// Sample operations (how much to move a column)
18921906
auto sample_operations = [&](std::vector<std::pair<HighsInt, double>>& moves,
18931907
HighsInt max_samples, HighsInt& n) {
18941908
HighsInt available = static_cast<HighsInt>(moves.size());
@@ -1903,6 +1917,7 @@ bool HighsPrimalHeuristics::localMip() {
19031917
}
19041918
};
19051919

1920+
// Propose column shift that improves the objective. Ignore feasibility.
19061921
auto explore_breakthrough =
19071922
[&](std::vector<std::pair<HighsInt, double>>& moves) {
19081923
for (HighsInt c : cols_with_obj) {
@@ -1912,15 +1927,15 @@ bool HighsPrimalHeuristics::localMip() {
19121927
std::max(static_cast<double>(obj - bestobj), feastol);
19131928
HighsInt dir = obj_coef < 0 ? 1 : -1;
19141929
double delta = obj_change / std::abs(obj_coef);
1915-
// TODO: THis feastol should be a negative.... Why does it work better?
1930+
// TODO: THis feastol should be a negative.... Why does it work
1931+
// better?
19161932
if (mipsolver.isColIntegral(c)) {
19171933
delta = std::ceil(delta - kHighsTiny);
19181934
}
19191935
delta =
19201936
dir == -1
19211937
? std::min(delta, std::abs(sol[c] - globaldom.col_lower_[c]))
19221938
: std::min(delta, std::abs(globaldom.col_upper_[c] - sol[c]));
1923-
// TODO: Need to add this in case integral variables have continuous bounds?
19241939
if (mipsolver.isColIntegral(c)) {
19251940
delta = std::floor(delta + feastol);
19261941
}
@@ -1931,6 +1946,8 @@ bool HighsPrimalHeuristics::localMip() {
19311946
}
19321947
};
19331948

1949+
// Iterate over sample violated rows and compute candidate moves
1950+
// Also compute breakthrough moves
19341951
auto explore_violated = [&](std::vector<std::pair<HighsInt, double>>& moves,
19351952
HighsInt& n) {
19361953
if (!viol.viol_index.empty()) {
@@ -1943,7 +1960,7 @@ bool HighsPrimalHeuristics::localMip() {
19431960
for (HighsInt i = start; i < end; ++i) {
19441961
HighsInt c = a_matrix_row.index_[i];
19451962
double val = a_matrix_row.value_[i];
1946-
feas_one_opt(c, r, val, moves);
1963+
calc_col_max_row_feas_delta(c, r, val, moves);
19471964
}
19481965
}
19491966
}
@@ -1953,6 +1970,7 @@ bool HighsPrimalHeuristics::localMip() {
19531970
sample_operations(moves, num_violated_rows_sample, n);
19541971
};
19551972

1973+
// Iterate over sample satisfied rows and compute candidate moves
19561974
auto explore_satisfied = [&](std::vector<std::pair<HighsInt, double>>& moves,
19571975
HighsInt& n) {
19581976
if (!satisfied.viol_index.empty()) {
@@ -1965,13 +1983,24 @@ bool HighsPrimalHeuristics::localMip() {
19651983
for (HighsInt i = start; i < end; ++i) {
19661984
HighsInt c = a_matrix_row.index_[i];
19671985
double val = a_matrix_row.value_[i];
1968-
feas_one_opt(c, r, val, moves);
1986+
calc_col_max_row_feas_delta(c, r, val, moves);
19691987
}
19701988
}
19711989
}
19721990
sample_operations(moves, num_satisfied_rows_sample, n);
19731991
};
19741992

1993+
// Score candidate moves. Given a move (c, delta), the score components are:
1994+
// + obj_weight if objective improves, - obj_weight otherwise
1995+
// For each row:
1996+
// + momentum * weights[row] if row feasibility improves
1997+
// - weights[row] if row feasibility worsens
1998+
// an additional multiplier of two is used if the row becomes feasible or
1999+
// becomes infeasible
2000+
// The candidate selected in order of priority:
2001+
// - Largest positive score
2002+
// - Last candidate with 0 score that improved feasibility TODO: Work on this
2003+
// - Either no candidate or with some small chance a random "bad" candidate
19752004
auto score_moves = [&](std::vector<std::pair<HighsInt, double>>& moves,
19762005
HighsInt n) {
19772006
HighsInt best_i = -1;
@@ -2015,10 +2044,10 @@ bool HighsPrimalHeuristics::localMip() {
20152044
} else if (was_satisfied && !now_satisfied) {
20162045
score -= 2 * weights[r];
20172046
} else if (!now_satisfied && !was_satisfied) {
2018-
if (old_violation > new_violation) {
2047+
if (old_violation > new_violation - kHighsTiny) {
20192048
score += momentum * weights[r];
20202049
feas_improvements++;
2021-
} else {
2050+
} else if (new_violation > old_violation - kHighsTiny) {
20222051
score -= weights[r];
20232052
}
20242053
}
@@ -2064,10 +2093,12 @@ bool HighsPrimalHeuristics::localMip() {
20642093

20652094
auto terminate = [&]() { return iters > 5000; };
20662095

2067-
calc_activites(false);
2096+
calc_activities(false);
20682097
bool last_iter_feas = false;
20692098
std::vector<std::pair<HighsInt, double>> candidate_moves;
2070-
candidate_moves.reserve(5000);
2099+
candidate_moves.reserve(std::max(
2100+
static_cast<HighsInt>(cols_with_obj.size()) + num_violated_rows_sample,
2101+
num_satisfied_rows_sample));
20712102
// Main solving loop
20722103
while (!terminate()) {
20732104
// No rows are violated, so check the solution
@@ -2079,7 +2110,7 @@ bool HighsPrimalHeuristics::localMip() {
20792110
momentum = 1;
20802111
recalc_objective();
20812112
// If solution is improving then store it
2082-
// TODO: Should we undo the transformation here??
2113+
// TODO: Undo the transformation here??
20832114
if (obj + feastol < bestobj) {
20842115
found_feas_before = true;
20852116
save_sol();
@@ -2103,6 +2134,7 @@ bool HighsPrimalHeuristics::localMip() {
21032134
for (HighsInt r : viol.viol_index) {
21042135
weights[r]++;
21052136
}
2137+
// Increase or reset momentum depending if candidate is found
21062138
if (candidate == -1) {
21072139
momentum++;
21082140
} else {
@@ -2113,6 +2145,7 @@ bool HighsPrimalHeuristics::localMip() {
21132145
// viol.viol_index.size(), found_feas_before,
21142146
// static_cast<double>(bestobj));
21152147
}
2148+
// Try and add incumbent to storage
21162149
if (found_feas_before) {
21172150
recalc_objective();
21182151
mipsolver.mipdata_->addIncumbent(intsol, static_cast<double>(obj),

0 commit comments

Comments
 (0)