Skip to content

Commit

Permalink
Now using reciprocal filter
Browse files Browse the repository at this point in the history
  • Loading branch information
@jajhall
jajhall committed Jun 14, 2024
1 parent d68df22 commit 392bc07
Show file tree
Hide file tree
Showing 4 changed files with 144 additions and 53 deletions.
14 changes: 10 additions & 4 deletions check/TestIis.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -175,7 +175,11 @@ TEST_CASE("lp-get-iis-galenet", "[iis]") {

TEST_CASE("lp-get-iis-avgas", "[iis]") {
std::string model = "avgas";
testMps(model, kIisStrategyFromLpRowPriority, HighsModelStatus::kOptimal);
// For the whole LP calculation the elasticity filter only
// identified feasibility, so the model status is not set
testMps(model, kIisStrategyFromLpRowPriority, HighsModelStatus::kNotset);
// For the ray calculation the model is solved, so its status is
// known
testMps(model, kIisStrategyFromRayRowPriority, HighsModelStatus::kOptimal);
}

Expand Down Expand Up @@ -263,8 +267,6 @@ void testIis(const std::string& model, const HighsIis& iis) {
"%s\n",
int(iisCol), int(iCol), iis.iisBoundStatusToString(iis_bound).c_str(),
highs.modelStatusToString(model_status).c_str());
// if (model_status != HighsModelStatus::kOptimal)
// highs.writeSolution("", kSolutionStylePretty);
REQUIRE(model_status == HighsModelStatus::kOptimal);
REQUIRE(highs.changeColBounds(iCol, lower, upper) == HighsStatus::kOk);
}
Expand Down Expand Up @@ -309,7 +311,11 @@ void testMps(std::string& model, const HighsInt iis_strategy,
highs.setOptionValue("output_flag", dev_run);

REQUIRE(highs.readModel(model_file) == HighsStatus::kOk);
REQUIRE(highs.run() == HighsStatus::kOk);
if (iis_strategy == kIisStrategyFromRayRowPriority ||
iis_strategy == kIisStrategyFromRayColPriority) {
// For a ray strategy, solve the LP first
REQUIRE(highs.run() == HighsStatus::kOk);
}
highs.setOptionValue("iis_strategy", iis_strategy);
HighsIis iis;
REQUIRE(highs.getIis(iis) == HighsStatus::kOk);
Expand Down
52 changes: 34 additions & 18 deletions src/lp_data/HighsIis.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,6 @@
* @brief Class-independent utilities for HiGHS
*/

#include "lp_data/HighsIis.h"

#include "Highs.h"

void HighsIis::invalidate() {
Expand All @@ -23,6 +21,7 @@ void HighsIis::invalidate() {
this->row_index_.clear();
this->col_bound_.clear();
this->row_bound_.clear();
this->info_.clear();
}

std::string HighsIis::iisBoundStatusToString(HighsInt bound_status) const {
Expand Down Expand Up @@ -216,7 +215,6 @@ HighsStatus HighsIis::getData(const HighsLp& lp, const HighsOptions& options,

HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options,
const HighsBasis* basis) {
this->invalidate();
const HighsLogOptions& log_options = options.log_options;
const bool row_priority =
options.iis_strategy == kIisStrategyFromRayRowPriority ||
Expand All @@ -225,6 +223,7 @@ HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options,
for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) this->addCol(iCol);
for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++) this->addRow(iRow);
Highs highs;
const HighsInfo& info = highs.getInfo();
highs.setOptionValue("output_flag", kIisDevReport);
highs.setOptionValue("presolve", kHighsOffString);
const HighsLp& incumbent_lp = highs.getLp();
Expand All @@ -239,9 +238,24 @@ HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options,
assert(run_status == HighsStatus::kOk);
// Solve the LP
if (basis) highs.setBasis(*basis);
run_status = highs.run();
assert(run_status == HighsStatus::kOk);

// Lambda for gathering data when solving an LP
auto solveLp = [&]() -> HighsStatus {
HighsIisInfo iis_info;
iis_info.simplex_time = -highs.getRunTime();
iis_info.simplex_iterations = -info.simplex_iteration_count;
run_status = highs.run();
assert(run_status == HighsStatus::kOk);
if (run_status != HighsStatus::kOk) return run_status;
iis_info.simplex_time += highs.getRunTime();
iis_info.simplex_iterations += info.simplex_iteration_count;
this->info_.push_back(iis_info);
return run_status;
};

run_status = solveLp();
if (run_status != HighsStatus::kOk) return run_status;

assert(highs.getModelStatus() == HighsModelStatus::kInfeasible);

const bool use_sensitivity_filter = false;
Expand All @@ -252,7 +266,8 @@ HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options,
//
highs.setOptionValue("output_flag", true);
// Solve the LP
run_status = highs.run();
run_status = solveLp();
if (run_status != HighsStatus::kOk) return run_status;
highs.writeSolution("", kSolutionStylePretty);
highs.setOptionValue("output_flag", output_flag);
}
Expand All @@ -275,15 +290,15 @@ HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options,
double upper = row_deletion ? lp.row_upper_[iX] : lp.col_upper_[iX];
// Record whether the upper bound has been dropped due to the lower bound
// being kept
bool drop_upper = false;
if (lower > -kHighsInf) {
// Drop the lower bound temporarily
run_status = row_deletion
? highs.changeRowBounds(iX, -kHighsInf, upper)
: highs.changeColBounds(iX, -kHighsInf, upper);
assert(run_status == HighsStatus::kOk);
run_status = highs.run();
assert(run_status == HighsStatus::kOk);
// Solve the LP
run_status = solveLp();
if (run_status != HighsStatus::kOk) return run_status;
HighsModelStatus model_status = highs.getModelStatus();
if (model_status == HighsModelStatus::kOptimal) {
// Now feasible, so restore the lower bound
Expand All @@ -292,18 +307,21 @@ HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options,
assert(run_status == HighsStatus::kOk);
// If the lower bound must be kept, then any finite upper bound
// must be dropped
const bool apply_reciprocal_rule = false;
const bool apply_reciprocal_rule = true;
if (apply_reciprocal_rule) {
if (upper < kHighsInf) {
// Drop the upper bound permanently
upper = kHighsInf;
run_status = row_deletion
? highs.changeRowBounds(iX, lower, kHighsInf)
? highs.changeRowBounds(iX, lower, upper)
: highs.changeColBounds(iX, lower, upper);
assert(run_status == HighsStatus::kOk);
drop_upper = true;
}
// continue;
assert(upper >= kHighsInf);
// Since upper = kHighsInf, allow the loop to run so that
// bound status is set as if upper were set to kHighsInf
// by relaxing it and finding that the LP was still
// infeasible
}
} else {
// Bound can be dropped permanently
Expand All @@ -316,12 +334,10 @@ HighsStatus HighsIis::compute(const HighsLp& lp, const HighsOptions& options,
run_status = row_deletion ? highs.changeRowBounds(iX, lower, kHighsInf)
: highs.changeColBounds(iX, lower, kHighsInf);
assert(run_status == HighsStatus::kOk);
run_status = highs.run();
assert(run_status == HighsStatus::kOk);
// Solve the LP
run_status = solveLp();
if (run_status != HighsStatus::kOk) return run_status;
HighsModelStatus model_status = highs.getModelStatus();
// If the upper bound has been dropped due to the reciprical
// rule, the LP must be infeasible
if (drop_upper) assert(model_status == HighsModelStatus::kInfeasible);
if (model_status == HighsModelStatus::kOptimal) {
// Now feasible, so restore the upper bound
run_status = row_deletion ? highs.changeRowBounds(iX, lower, upper)
Expand Down
6 changes: 6 additions & 0 deletions src/lp_data/HighsIis.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,11 @@ enum IisBoundStatus {
kIisBoundStatusBoxed // 4
};

struct HighsIisInfo {
double simplex_time = 0;
HighsInt simplex_iterations = 0;
};

class HighsIis {
public:
HighsIis() {}
Expand Down Expand Up @@ -54,6 +59,7 @@ class HighsIis {
std::vector<HighsInt> row_index_;
std::vector<HighsInt> col_bound_;
std::vector<HighsInt> row_bound_;
std::vector<HighsIisInfo> info_;
};

#endif // LP_DATA_HIGHSIIS_H_
125 changes: 94 additions & 31 deletions src/lp_data/HighsInterface.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1573,9 +1573,17 @@ HighsStatus Highs::getIisInterface() {
// detected in presolve, so solve without presolve
std::string presolve = options_.presolve;
options_.presolve = kHighsOffString;
HighsStatus return_status = this->run();

HighsIisInfo iis_info;
iis_info.simplex_time = -this->getRunTime();
iis_info.simplex_iterations = -info_.simplex_iteration_count;
HighsStatus run_status = this->run();
options_.presolve = presolve;
if (return_status != HighsStatus::kOk) return return_status;
if (run_status != HighsStatus::kOk) return run_status;
iis_info.simplex_time += this->getRunTime();
iis_info.simplex_iterations += -info_.simplex_iteration_count;
this->iis_.info_.push_back(iis_info);

// Model should remain infeasible!
if (this->model_status_ != HighsModelStatus::kInfeasible) {
highsLogUser(
Expand All @@ -1600,7 +1608,7 @@ HighsStatus Highs::getIisInterface() {
std::vector<double> rhs;
HighsInt iRow = ekk_instance_.info_.dual_ray_row_;
rhs.assign(num_row, 0);
rhs[iRow] = 1; // ekk_instance_.info_.dual_ray_sign_;
rhs[iRow] = 1;
std::vector<double> dual_ray_value(num_row);
HighsInt* dual_ray_num_nz = 0;
basisSolveInterface(rhs, dual_ray_value.data(), dual_ray_num_nz, NULL,
Expand All @@ -1622,12 +1630,45 @@ HighsStatus Highs::getIisInterface() {
assert(check_lp_before.equalButForScalingAndNames(check_lp_after));
if (return_status != HighsStatus::kOk) return return_status;
}
HighsStatus return_status =
this->iis_.getData(lp, options_, basis_, infeasible_row_subset);
if (return_status == HighsStatus::kOk) {
// Existence of non-empty IIS => infeasibility
if (this->iis_.col_index_.size() > 0 || this->iis_.row_index_.size() > 0)
this->model_status_ = HighsModelStatus::kInfeasible;
HighsStatus return_status = HighsStatus::kOk;
if (infeasible_row_subset.size() == 0) {
// No subset of infeasible rows, so model is feasible
this->iis_.valid_ = true;
} else {
return_status =
this->iis_.getData(lp, options_, basis_, infeasible_row_subset);
if (return_status == HighsStatus::kOk) {
// Existence of non-empty IIS => infeasibility
if (this->iis_.col_index_.size() > 0 || this->iis_.row_index_.size() > 0)
this->model_status_ = HighsModelStatus::kInfeasible;
}
// Analyse the LP solution data
const HighsInt num_lp_solved = this->iis_.info_.size();
double min_time = kHighsInf;
double sum_time = 0;
double max_time = 0;
HighsInt min_iterations = kHighsIInf;
HighsInt sum_iterations = 0;
HighsInt max_iterations = 0;
for (HighsInt iX = 0; iX < num_lp_solved; iX++) {
double time = this->iis_.info_[iX].simplex_time;
HighsInt iterations = this->iis_.info_[iX].simplex_iterations;
min_time = std::min(time, min_time);
sum_time += time;
max_time = std::max(time, max_time);
min_iterations = std::min(iterations, min_iterations);
sum_iterations += iterations;
max_iterations = std::max(iterations, max_iterations);
}
highsLogUser(options_.log_options, HighsLogType::kInfo,
" (min / average / max) iteration count (%6d / %6.2g / % 6d)"
" and time (%6.2f / %6.2f / % 6.2f) \n",
int(this->iis_.col_index_.size()),
int(this->iis_.row_index_.size()), int(num_lp_solved),
int(min_iterations),
num_lp_solved > 0 ? (1.0 * sum_iterations) / num_lp_solved : 0,
int(max_iterations), min_time,
num_lp_solved > 0 ? sum_time / num_lp_solved : 0, max_time);
}
return return_status;
}
Expand Down Expand Up @@ -1843,15 +1884,31 @@ HighsStatus Highs::computeInfeasibleRows(
this->writeModel("");
this->setOptionValue("output_flag", output_flag);
}
run_status = this->run();
assert(run_status == HighsStatus::kOk);
// Lambda for gathering data when solving an LP
auto solveLp = [&]() -> HighsStatus {
HighsIisInfo iis_info;
iis_info.simplex_time = -this->getRunTime();
iis_info.simplex_iterations = -info_.simplex_iteration_count;
run_status = this->run();
assert(run_status == HighsStatus::kOk);
if (run_status != HighsStatus::kOk) return run_status;
iis_info.simplex_time += this->getRunTime();
iis_info.simplex_iterations += info_.simplex_iteration_count;
this->iis_.info_.push_back(iis_info);
return run_status;
};

run_status = solveLp();
if (run_status != HighsStatus::kOk) return run_status;
if (kIisDevReport) this->writeSolution("", kSolutionStylePretty);
HighsModelStatus model_status = this->getModelStatus();
assert(model_status == HighsModelStatus::kOptimal);
// Model status should be optimal, unless model is unbounded
assert(this->model_status_ == HighsModelStatus::kOptimal ||
this->model_status_ == HighsModelStatus::kUnbounded);

const HighsSolution& solution = this->getSolution();
// Now fix e-variables that are positive and re-solve until e-LP is infeasible
HighsInt loop_k = 0;
bool feasible_model = false;
for (;;) {
if (kIisDevReport)
printf("\nElasticity filter pass %d\n==============\n", int(loop_k));
Expand Down Expand Up @@ -1889,14 +1946,17 @@ HighsStatus Highs::computeInfeasibleRows(
num_fixed++;
}
}
assert(num_fixed > 0);
run_status = this->run();
assert(run_status == HighsStatus::kOk);
if (num_fixed == 0) {
// No elastic variables were positive, so problem is feasible
feasible_model = true;
break;
}
HighsStatus run_status = solveLp();
if (run_status != HighsStatus::kOk) return run_status;
if (kIisDevReport) this->writeSolution("", kSolutionStylePretty);
HighsModelStatus model_status = this->getModelStatus();
if (model_status == HighsModelStatus::kInfeasible) break;
loop_k++;
if (loop_k > 10) assert(1666 == 1999);
}

infeasible_row_subset.clear();
Expand All @@ -1917,26 +1977,31 @@ HighsStatus Highs::computeInfeasibleRows(
HighsInt iRow = row_of_ecol[eCol];
if (lp.col_upper_[row_ecol_offset + eCol] == 0) {
num_enforced_row_ecol++;
infeasible_row_subset.push_back(iRow);
if (kIisDevReport)
printf(
"Row e-col %2d (column %2d) corresponds to row %2d with bound "
"%g "
"and is enforced\n",
"%g and is enforced\n",
int(eCol), int(row_ecol_offset + eCol), int(iRow),
bound_of_row_of_ecol[eCol]);
infeasible_row_subset.push_back(iRow);
}
}
if (kIisDevReport) {

if (feasible_model)
assert(num_enforced_col_ecol == 0 && num_enforced_row_ecol == 0);

highsLogUser(
options_.log_options, HighsLogType::kInfo,
"Elasticity filter after %d passes enforces bounds on %d cols and %d "
"rows\n",
int(loop_k), int(num_enforced_col_ecol), int(num_enforced_row_ecol));

if (kIisDevReport)
printf(
"\nElasticity filter after %d passes enforces bounds on %d cols and %d "
"rows\n",
int(loop_k), int(num_enforced_col_ecol), int(num_enforced_row_ecol));
printf(
"Highs::computeInfeasibleRows: Before clearing additional rows and "
"columns - model status is %s\n",
this->modelStatusToString(this->model_status_).c_str());
}

// Delete any additional rows and columns, and restore the original
// column costs and bounds
run_status = this->deleteRows(original_num_row, lp.num_row_ - 1);
Expand All @@ -1957,11 +2022,9 @@ HighsStatus Highs::computeInfeasibleRows(
assert(lp.num_col_ == original_num_col);
assert(lp.num_row_ == original_num_row);

if (kIisDevReport)
printf(
"Highs::computeInfeasibleRows: After clearing additional rows and "
"columns - model status is %s\n",
this->modelStatusToString(this->model_status_).c_str());
// If the model is feasible, then the status of model is not known
if (feasible_model) this->model_status_ = HighsModelStatus::kNotset;

return HighsStatus::kOk;
}

Expand Down

0 comments on commit 392bc07

Please sign in to comment.