From 1df27e21b55c3a354f9504acf2e262a1bd176ea6 Mon Sep 17 00:00:00 2001 From: jajhall Date: Mon, 14 Oct 2024 20:38:04 +0100 Subject: [PATCH 1/7] Recovered useful code from fix-1966 --- src/lp_data/HighsOptions.h | 14 ++++++++++++++ src/util/HighsSparseMatrix.cpp | 26 ++++++++++++++++++++++++++ src/util/HighsSparseMatrix.h | 3 +++ 3 files changed, 43 insertions(+) diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h index 208bcd1802..bf853a4839 100644 --- a/src/lp_data/HighsOptions.h +++ b/src/lp_data/HighsOptions.h @@ -389,6 +389,8 @@ struct HighsOptionsStruct { bool less_infeasible_DSE_choose_row; bool use_original_HFactor_logic; bool run_centring; + double primal_residual_tolerance; + double dual_residual_tolerance; HighsInt max_centring_steps; double centring_ratio_tolerance; @@ -524,6 +526,8 @@ struct HighsOptionsStruct { run_centring(false), max_centring_steps(0), centring_ratio_tolerance(0.0), + primal_residual_tolerance(0.0), + dual_residual_tolerance(0.0), icrash(false), icrash_dualize(false), icrash_strategy(""), @@ -1378,6 +1382,16 @@ class HighsOptions : public HighsOptionsStruct { advanced, ¢ring_ratio_tolerance, 0, 100, kHighsInf); records.push_back(record_double); + record_double = new OptionRecordDouble( + "primal_residual_tolerance", "Primal residual tolerance", advanced, + &primal_residual_tolerance, 1e-10, 1e-7, kHighsInf); + records.push_back(record_double); + + record_double = new OptionRecordDouble( + "dual_residual_tolerance", "Dual residual tolerance", advanced, + &dual_residual_tolerance, 1e-10, 1e-7, kHighsInf); + records.push_back(record_double); + // Set up the log_options aliases log_options.clear(); log_options.log_stream = diff --git a/src/util/HighsSparseMatrix.cpp b/src/util/HighsSparseMatrix.cpp index 16648e1506..a493d2e7ad 100644 --- a/src/util/HighsSparseMatrix.cpp +++ b/src/util/HighsSparseMatrix.cpp @@ -1174,6 +1174,32 @@ void HighsSparseMatrix::productQuad(vector& result, } } +void HighsSparseMatrix::productTransposeQuad( + vector& result, const vector& row, + const HighsInt debug_report) const { + assert(this->formatOk()); + assert((int)row.size() >= this->num_row_); + result.assign(this->num_col_, 0.0); + if (this->isColwise()) { + for (HighsInt iCol = 0; iCol < this->num_col_; iCol++) { + HighsCDouble value = 0.0; + for (HighsInt iEl = this->start_[iCol]; iEl < this->start_[iCol + 1]; + iEl++) + value += row[this->index_[iEl]] * this->value_[iEl]; + result[iCol] = double(value); + } + } else { + std::vector value(this->num_col_, 0); + for (HighsInt iRow = 0; iRow < this->num_row_; iRow++) { + for (HighsInt iEl = this->start_[iRow]; iEl < this->start_[iRow + 1]; + iEl++) + value[this->index_[iEl]] += row[iRow] * this->value_[iEl]; + } + for (HighsInt iCol = 0; iCol < this->num_col_; iCol++) + result[iCol] = double(value[iCol]); + } +} + void HighsSparseMatrix::productTransposeQuad( vector& result_value, vector& result_index, const HVector& column, const HighsInt debug_report) const { diff --git a/src/util/HighsSparseMatrix.h b/src/util/HighsSparseMatrix.h index 22fef9ea40..7333a37e64 100644 --- a/src/util/HighsSparseMatrix.h +++ b/src/util/HighsSparseMatrix.h @@ -94,6 +94,9 @@ class HighsSparseMatrix { void productTranspose(vector& result, const vector& x) const; void productQuad(vector& result, const vector& x, const HighsInt debug_report = kDebugReportOff) const; + void productTransposeQuad( + vector& result_value, const vector& x, + const HighsInt debug_report = kDebugReportOff) const; void productTransposeQuad( vector& result_value, vector& result_index, const HVector& x, const HighsInt debug_report = kDebugReportOff) const; From 30ecb6ff9084baf03d2032e76781ec5d9dd9a10a Mon Sep 17 00:00:00 2001 From: jajhall Date: Mon, 14 Oct 2024 20:45:14 +0100 Subject: [PATCH 2/7] And the unit tests from fix-1966 --- check/TestIpm.cpp | 29 +++++++++++++++++++++++++++++ check/TestLpSolvers.cpp | 29 +++++++++++++++++++++++++++++ 2 files changed, 58 insertions(+) diff --git a/check/TestIpm.cpp b/check/TestIpm.cpp index 911d01a67a..e9201c89bb 100644 --- a/check/TestIpm.cpp +++ b/check/TestIpm.cpp @@ -81,3 +81,32 @@ TEST_CASE("test-analytic-centre-box", "[highs_ipm]") { REQUIRE(solution_norm < 1e-6); if (dev_run) printf("Analytic centre solution norm is %g\n", solution_norm); } + +TEST_CASE("test-1966", "[highs_ipm]") { + Highs highs; + // highs.setOptionValue("output_flag", dev_run); + HighsLp lp; + lp.num_col_ = 2; + lp.num_row_ = 2; + lp.col_cost_ = {2 - 1}; + lp.col_lower_ = {0, 0}; + lp.col_upper_ = {kHighsInf, kHighsInf}; + lp.row_lower_ = {-kHighsInf, 2}; + lp.row_upper_ = {1, kHighsInf}; + lp.a_matrix_.start_ = {0, 2, 4}; + lp.a_matrix_.index_ = {0, 1, 0, 1}; + lp.a_matrix_.value_ = {1, -1, 1, -1}; + highs.passModel(lp); + highs.setOptionValue("solver", kIpmString); + highs.setOptionValue("presolve", kHighsOffString); + highs.run(); + HighsInfo info = highs.getInfo(); + printf("Num primal infeasibilities = %d\n", + int(info.num_primal_infeasibilities)); + printf("Max primal infeasibilities = %g\n", info.max_primal_infeasibility); + printf("Sum primal infeasibilities = %g\n", info.sum_primal_infeasibilities); + printf("Num dual infeasibilities = %d\n", + int(info.num_dual_infeasibilities)); + printf("Max dual infeasibilities = %g\n", info.max_dual_infeasibility); + printf("Sum dual infeasibilities = %g\n", info.sum_dual_infeasibilities); +} diff --git a/check/TestLpSolvers.cpp b/check/TestLpSolvers.cpp index 5202a1dee5..cedf9d6a83 100644 --- a/check/TestLpSolvers.cpp +++ b/check/TestLpSolvers.cpp @@ -438,3 +438,32 @@ TEST_CASE("dual-objective-upper-bound", "[highs_lp_solver]") { if (dev_run) printf("\nOptimal objective value error = %g\n", error); REQUIRE(error < 1e-10); } + +TEST_CASE("blending-lp", "[highs_lp_solver]") { + Highs highs; + // highs.setOptionValue("output_flag", dev_run); + HighsLp lp; + lp.num_col_ = 2; + lp.num_row_ = 2; + lp.col_cost_ = {-8, -10}; + lp.col_lower_ = {0, 0}; + lp.col_upper_ = {kHighsInf, kHighsInf}; + lp.row_lower_ = {-kHighsInf, -kHighsInf}; + lp.row_upper_ = {80, 120}; + lp.a_matrix_.start_ = {0, 2, 4}; + lp.a_matrix_.index_ = {0, 1, 0, 1}; + lp.a_matrix_.value_ = {1, 1, 2, 4}; + highs.passModel(lp); + highs.setOptionValue("solver", kIpmString); + highs.setOptionValue("presolve", kHighsOffString); + highs.run(); + HighsInfo info = highs.getInfo(); + printf("Num primal infeasibilities = %d\n", + int(info.num_primal_infeasibilities)); + printf("Max primal infeasibilities = %g\n", info.max_primal_infeasibility); + printf("Sum primal infeasibilities = %g\n", info.sum_primal_infeasibilities); + printf("Num dual infeasibilities = %d\n", + int(info.num_dual_infeasibilities)); + printf("Max dual infeasibilities = %g\n", info.max_dual_infeasibility); + printf("Sum dual infeasibilities = %g\n", info.sum_dual_infeasibilities); +} From 7ccf279e882b706f486040401f564c43a85275a0 Mon Sep 17 00:00:00 2001 From: jajhall Date: Mon, 14 Oct 2024 20:46:17 +0100 Subject: [PATCH 3/7] Formatted --- check/TestNames.cpp | 3 +-- src/lp_data/HighsOptions.h | 2 +- src/pdlp/CupdlpWrapper.cpp | 14 ++++++-------- src/util/HFactor.h | 2 +- 4 files changed, 9 insertions(+), 12 deletions(-) diff --git a/check/TestNames.cpp b/check/TestNames.cpp index f0d11fc4a8..4aceee0ca6 100644 --- a/check/TestNames.cpp +++ b/check/TestNames.cpp @@ -120,13 +120,12 @@ TEST_CASE("highs-names", "[highs_names]") { } TEST_CASE("highs-model-name", "[model_names]") { - Highs highs; const HighsLp& lp = highs.getLp(); std::string name = lp.model_name_; REQUIRE(name == ""); - + highs.passModelName("new_name"); name = lp.model_name_; REQUIRE(name == "new_name"); diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h index bf853a4839..9b38657551 100644 --- a/src/lp_data/HighsOptions.h +++ b/src/lp_data/HighsOptions.h @@ -559,7 +559,7 @@ struct HighsOptionsStruct { #endif mip_improving_solution_save(false), mip_improving_solution_report_sparse(false), - mip_improving_solution_file("") {}; + mip_improving_solution_file(""){}; }; // For now, but later change so HiGHS properties are string based so that new diff --git a/src/pdlp/CupdlpWrapper.cpp b/src/pdlp/CupdlpWrapper.cpp index d41f14e595..31fa31e1e8 100644 --- a/src/pdlp/CupdlpWrapper.cpp +++ b/src/pdlp/CupdlpWrapper.cpp @@ -207,7 +207,7 @@ HighsStatus solveLpCupdlp(const HighsOptions& options, HighsTimer& timer, #if CUPDLP_DEBUG analysePdlpSolution(options, lp, highs_solution); #endif - + free(cost); free(lower); free(upper); @@ -225,7 +225,7 @@ HighsStatus solveLpCupdlp(const HighsOptions& options, HighsTimer& timer, free(prob->lower); free(prob->upper); free(prob->rhs); - + free(prob->hasLower); free(prob->hasUpper); @@ -233,7 +233,7 @@ HighsStatus solveLpCupdlp(const HighsOptions& options, HighsTimer& timer, free(prob->data->csr_matrix->rowMatIdx); free(prob->data->csr_matrix->rowMatElem); free(prob->data->csr_matrix); - + free(prob->data->csc_matrix->colMatBeg); free(prob->data->csc_matrix->colMatIdx); free(prob->data->csc_matrix->colMatElem); @@ -246,13 +246,11 @@ HighsStatus solveLpCupdlp(const HighsOptions& options, HighsTimer& timer, free(csc_cpu->colMatBeg); free(csc_cpu->colMatIdx); free(csc_cpu->colMatElem); - + free(csc_cpu); - if (scaling->rowScale != nullptr) - free(scaling->rowScale); - if (scaling->colScale != nullptr) - free(scaling->colScale); + if (scaling->rowScale != nullptr) free(scaling->rowScale); + if (scaling->colScale != nullptr) free(scaling->colScale); free(scaling); return HighsStatus::kOk; diff --git a/src/util/HFactor.h b/src/util/HFactor.h index 03cb8d8b1d..2f75260d8a 100644 --- a/src/util/HFactor.h +++ b/src/util/HFactor.h @@ -132,7 +132,7 @@ class HFactor { build_timer_(nullptr), nwork(0), u_merit_x(0), - u_total_x(0) {}; + u_total_x(0){}; /** * @brief Copy problem size and pointers of constraint matrix, and set From 40b35a33b286de2be20b781fa5c5b0a09821a5b7 Mon Sep 17 00:00:00 2001 From: Julian Hall Date: Tue, 15 Oct 2024 16:59:32 +0100 Subject: [PATCH 4/7] Now checking residual errors after IPM and PDLP, and correcting row values and column duals accordingly --- check/TestIpm.cpp | 45 ++++++++++++----- check/TestLpSolvers.cpp | 22 +++++---- src/lp_data/Highs.cpp | 4 -- src/lp_data/HighsSolution.cpp | 92 +++++++++++++++++++++++++++++++++-- src/lp_data/HighsSolution.h | 3 ++ src/lp_data/HighsSolve.cpp | 8 +++ 6 files changed, 144 insertions(+), 30 deletions(-) diff --git a/check/TestIpm.cpp b/check/TestIpm.cpp index e9201c89bb..ff60bbd403 100644 --- a/check/TestIpm.cpp +++ b/check/TestIpm.cpp @@ -84,11 +84,12 @@ TEST_CASE("test-analytic-centre-box", "[highs_ipm]") { TEST_CASE("test-1966", "[highs_ipm]") { Highs highs; - // highs.setOptionValue("output_flag", dev_run); + highs.setOptionValue("output_flag", dev_run); + const HighsInfo& info = highs.getInfo(); HighsLp lp; lp.num_col_ = 2; lp.num_row_ = 2; - lp.col_cost_ = {2 - 1}; + lp.col_cost_ = {2, -1}; lp.col_lower_ = {0, 0}; lp.col_upper_ = {kHighsInf, kHighsInf}; lp.row_lower_ = {-kHighsInf, 2}; @@ -99,14 +100,36 @@ TEST_CASE("test-1966", "[highs_ipm]") { highs.passModel(lp); highs.setOptionValue("solver", kIpmString); highs.setOptionValue("presolve", kHighsOffString); + + if (dev_run) printf("\nWith default residual tolerances\n"); + highs.run(); + if (dev_run) { + highs.writeSolution("", kSolutionStylePretty); + printf("Num primal infeasibilities = %d\n", + int(info.num_primal_infeasibilities)); + printf("Max primal infeasibility = %g\n", info.max_primal_infeasibility); + printf("Sum primal infeasibilities = %g\n", info.sum_primal_infeasibilities); + printf("Num dual infeasibilities = %d\n", + int(info.num_dual_infeasibilities)); + printf("Max dual infeasibility = %g\n", info.max_dual_infeasibility); + printf("Sum dual infeasibilities = %g\n", info.sum_dual_infeasibilities); + } + highs.clearSolver(); + + if (dev_run) + printf("\nWith infinite residual tolerances\n"); + highs.setOptionValue("primal_residual_tolerance", 1e30); + highs.setOptionValue("dual_residual_tolerance", 1e30); highs.run(); - HighsInfo info = highs.getInfo(); - printf("Num primal infeasibilities = %d\n", - int(info.num_primal_infeasibilities)); - printf("Max primal infeasibilities = %g\n", info.max_primal_infeasibility); - printf("Sum primal infeasibilities = %g\n", info.sum_primal_infeasibilities); - printf("Num dual infeasibilities = %d\n", - int(info.num_dual_infeasibilities)); - printf("Max dual infeasibilities = %g\n", info.max_dual_infeasibility); - printf("Sum dual infeasibilities = %g\n", info.sum_dual_infeasibilities); + if (dev_run) { + highs.writeSolution("", kSolutionStylePretty); + printf("Num primal infeasibilities = %d\n", + int(info.num_primal_infeasibilities)); + printf("Max primal infeasibility = %g\n", info.max_primal_infeasibility); + printf("Sum primal infeasibilities = %g\n", info.sum_primal_infeasibilities); + printf("Num dual infeasibilities = %d\n", + int(info.num_dual_infeasibilities)); + printf("Max dual infeasibility = %g\n", info.max_dual_infeasibility); + printf("Sum dual infeasibilities = %g\n", info.sum_dual_infeasibilities); + } } diff --git a/check/TestLpSolvers.cpp b/check/TestLpSolvers.cpp index cedf9d6a83..e4f8d4bc97 100644 --- a/check/TestLpSolvers.cpp +++ b/check/TestLpSolvers.cpp @@ -439,9 +439,9 @@ TEST_CASE("dual-objective-upper-bound", "[highs_lp_solver]") { REQUIRE(error < 1e-10); } -TEST_CASE("blending-lp", "[highs_lp_solver]") { +TEST_CASE("blending-lp-ipm", "[highs_lp_solver]") { Highs highs; - // highs.setOptionValue("output_flag", dev_run); + highs.setOptionValue("output_flag", dev_run); HighsLp lp; lp.num_col_ = 2; lp.num_row_ = 2; @@ -458,12 +458,14 @@ TEST_CASE("blending-lp", "[highs_lp_solver]") { highs.setOptionValue("presolve", kHighsOffString); highs.run(); HighsInfo info = highs.getInfo(); - printf("Num primal infeasibilities = %d\n", - int(info.num_primal_infeasibilities)); - printf("Max primal infeasibilities = %g\n", info.max_primal_infeasibility); - printf("Sum primal infeasibilities = %g\n", info.sum_primal_infeasibilities); - printf("Num dual infeasibilities = %d\n", - int(info.num_dual_infeasibilities)); - printf("Max dual infeasibilities = %g\n", info.max_dual_infeasibility); - printf("Sum dual infeasibilities = %g\n", info.sum_dual_infeasibilities); + if (dev_run) { + printf("Num primal infeasibilities = %d\n", + int(info.num_primal_infeasibilities)); + printf("Max primal infeasibility = %g\n", info.max_primal_infeasibility); + printf("Sum primal infeasibilities = %g\n", info.sum_primal_infeasibilities); + printf("Num dual infeasibilities = %d\n", + int(info.num_dual_infeasibilities)); + printf("Max dual infeasibility = %g\n", info.max_dual_infeasibility); + printf("Sum dual infeasibilities = %g\n", info.sum_dual_infeasibilities); + } } diff --git a/src/lp_data/Highs.cpp b/src/lp_data/Highs.cpp index c420839a5d..172df328c4 100644 --- a/src/lp_data/Highs.cpp +++ b/src/lp_data/Highs.cpp @@ -1959,10 +1959,6 @@ HighsStatus Highs::getReducedColumn(const HighsInt col, double* col_vector, HighsStatus Highs::getKappa(double& kappa, const bool exact, const bool report) { - printf( - "Highs::getKappa basis_.valid = %d, ekk_instance_.status_.has_invert = " - "%d\n", - int(basis_.valid), int(ekk_instance_.status_.has_invert)); if (!ekk_instance_.status_.has_invert) return invertRequirementError("getBasisInverseRow"); kappa = ekk_instance_.computeBasisCondition(this->model_.lp_, exact, report); diff --git a/src/lp_data/HighsSolution.cpp b/src/lp_data/HighsSolution.cpp index 29e40168a0..6da9ba4ab4 100644 --- a/src/lp_data/HighsSolution.cpp +++ b/src/lp_data/HighsSolution.cpp @@ -282,7 +282,7 @@ void getKktFailures(const HighsOptions& options, const HighsLp& lp, relative_primal_infeasibility, dual_infeasibility, value_residual); if (!status_value_ok) highsLogUser(options.log_options, HighsLogType::kError, - "getKktFailures: %s %d status-value error: [%g; %g; %g] has " + "getKktFailures: %s %d status-value error: [%23.18g; %23.18g; %23.18g] has " "residual %g\n", iVar < lp.num_col_ ? "Column" : "Row ", iVar < lp.num_col_ ? int(iVar) : int(iVar - lp.num_col_), @@ -531,14 +531,22 @@ bool getVariableKktFailures(const double primal_feasibility_tolerance, // If the variable is basic, then consider it not to be at a bound // so that any dual value yields an infeasibility value if (*status_pointer == HighsBasisStatus::kBasic) at_a_bound = false; + // With very large values, accuracy is lost in adding/subtracting + // the feasibility tolerance from the bounds, so skip if this may + // occur + // // Check that kLower and kUpper are consistent with value and // bounds - for debugging QP basis errors if (*status_pointer == HighsBasisStatus::kLower) { - status_value_ok = value >= lower - primal_feasibility_tolerance && - value <= lower + primal_feasibility_tolerance; + if (std::fabs(lower) / primal_feasibility_tolerance < 1e-16) + status_value_ok = + value >= lower - primal_feasibility_tolerance && + value <= lower + primal_feasibility_tolerance; } else if (*status_pointer == HighsBasisStatus::kUpper) { - status_value_ok = value >= upper - primal_feasibility_tolerance && - value <= upper + primal_feasibility_tolerance; + if (std::fabs(upper) / primal_feasibility_tolerance < 1e-16) + status_value_ok = + value >= upper - primal_feasibility_tolerance && + value <= upper + primal_feasibility_tolerance; } } if (at_a_bound) { @@ -1369,6 +1377,80 @@ void accommodateAlienBasis(HighsLpSolverObject& solver_object) { assert(num_basic_variables == num_row); } +void correctResiduals(HighsLpSolverObject& solver_object) { + HighsOptions& options = solver_object.options_; + HighsSolution& solution = solver_object.solution_; + const bool& have_primal_solution = solution.value_valid; + const bool& have_dual_solution = solution.dual_valid; + assert(have_primal_solution); + assert(have_dual_solution); + HighsLp& lp = solver_object.lp_; + std::vector check_row_value; + std::vector check_col_dual; + lp.a_matrix_.productQuad(check_row_value, solution.col_value); + if (have_dual_solution) { + lp.a_matrix_.productTransposeQuad(check_col_dual, solution.row_dual); + for (HighsInt iCol=0; iCol < lp.num_col_; iCol++) + check_col_dual[iCol] -= lp.col_cost_[iCol]; + } + double norm_primal_residual = 0; + HighsInt num_primal_correction = 0; + double max_primal_correction = 0; + double sum_primal_correction = 0; + double use_primal_residual_tolerance = options.primal_residual_tolerance; + + double norm_dual_residual = 0; + HighsInt num_dual_correction = 0; + double max_dual_correction = 0; + double sum_dual_correction = 0; + double use_dual_residual_tolerance = options.dual_residual_tolerance; + + for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++) { + double primal_residual = check_row_value[iRow] - solution.row_value[iRow]; + double abs_primal_residual = std::fabs(primal_residual); + // printf("correctResiduals: Row %d; check_row_value = %11.4g; row_value = %11.4g; primal_residual = %11.4g\n", + // int(iRow), check_row_value[iRow], solution.row_value[iRow], abs_primal_residual); + if (abs_primal_residual > use_primal_residual_tolerance) { + solution.row_value[iRow] += primal_residual; + double check_primal_residual = std::fabs(check_row_value[iRow] - solution.row_value[iRow]); + assert(check_primal_residual < 1e-10); + num_primal_correction++; + max_primal_correction = std::max(abs_primal_residual, max_primal_correction); + sum_primal_correction += abs_primal_residual; + } + norm_primal_residual = std::max(abs_primal_residual, norm_primal_residual); + } + if (have_dual_solution) { + for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) { + double dual_residual = check_col_dual[iCol] + solution.col_dual[iCol]; + double abs_dual_residual = std::fabs(dual_residual); + // printf("correctResiduals: Col %d; check_col_dual = %11.4g; col_dual = %11.4g; col_cost_ = %11.4g; dual_residual = %11.4g\n", + // int(iCol), check_col_dual[iCol], solution.col_dual[iCol], lp.col_cost_[iCol], abs_dual_residual); + if (abs_dual_residual > use_dual_residual_tolerance) { + solution.col_dual[iCol] -= dual_residual; + double check_dual_residual = std::fabs(check_col_dual[iCol] + solution.col_dual[iCol]); + const bool ok_dual_correction = check_dual_residual < 1e-10; + // if (!ok_dual_correction) { + // printf("correctResiduals: Col %d; check_col_dual = %11.4g; col_dual = %11.4g; check_dual_residual = %g\n", + // int(iCol), check_col_dual[iCol], solution.col_dual[iCol], check_dual_residual); + // } + assert(ok_dual_correction); + num_dual_correction++; + max_dual_correction = std::max(abs_dual_residual, max_dual_correction); + sum_dual_correction += abs_dual_residual; + } + norm_dual_residual = std::max(abs_dual_residual, norm_dual_residual); + } + } + + if (num_primal_correction > 0 || num_dual_correction > 0) + highsLogUser(options.log_options, HighsLogType::kWarning, + "LP solver yields |primal_residual| = %g; |dual_residual| = %g with num / max / sum primal (%d / %g / %g) and dual (%d / %g / %g) corrections\n", + norm_primal_residual, norm_dual_residual, + int(num_primal_correction), max_primal_correction, sum_primal_correction, + int(num_dual_correction), max_dual_correction, sum_dual_correction); +} + void resetModelStatusAndHighsInfo(HighsLpSolverObject& solver_object) { solver_object.model_status_ = HighsModelStatus::kNotset; solver_object.highs_info_.objective_function_value = 0; diff --git a/src/lp_data/HighsSolution.h b/src/lp_data/HighsSolution.h index f869b6b0ee..1efbedda52 100644 --- a/src/lp_data/HighsSolution.h +++ b/src/lp_data/HighsSolution.h @@ -123,6 +123,8 @@ HighsStatus formSimplexLpBasisAndFactor( void accommodateAlienBasis(HighsLpSolverObject& solver_object); +void correctResiduals(HighsLpSolverObject& solver_object); + void resetModelStatusAndHighsInfo(HighsLpSolverObject& solver_object); void resetModelStatusAndHighsInfo(HighsModelStatus& model_status, HighsInfo& highs_info); @@ -134,4 +136,5 @@ bool isDualSolutionRightSize(const HighsLp& lp, const HighsSolution& solution); bool isSolutionRightSize(const HighsLp& lp, const HighsSolution& solution); bool isBasisRightSize(const HighsLp& lp, const HighsBasis& basis); + #endif // LP_DATA_HIGHSSOLUTION_H_ diff --git a/src/lp_data/HighsSolve.cpp b/src/lp_data/HighsSolve.cpp index c03837ac01..9e11c91888 100644 --- a/src/lp_data/HighsSolve.cpp +++ b/src/lp_data/HighsSolve.cpp @@ -72,6 +72,10 @@ HighsStatus solveLp(HighsLpSolverObject& solver_object, const string message) { return_status, "solveLpCupdlp"); } if (return_status == HighsStatus::kError) return return_status; + // IPM (and PDLP?) can claim optimality with large primal and/or + // dual residual errors, so correct any residual errors that + // exceed the tolerance + correctResiduals(solver_object); // Non-error return requires a primal solution assert(solver_object.solution_.value_valid); // Get the objective and any KKT failures @@ -79,6 +83,10 @@ HighsStatus solveLp(HighsLpSolverObject& solver_object, const string message) { solver_object.lp_.objectiveValue(solver_object.solution_.col_value); getLpKktFailures(options, solver_object.lp_, solver_object.solution_, solver_object.basis_, solver_object.highs_info_); + if (solver_object.model_status_ == HighsModelStatus::kOptimal && + solver_object.highs_info_.num_primal_infeasibilities > 0 || + solver_object.highs_info_.num_dual_infeasibilities) + solver_object.model_status_ = HighsModelStatus::kUnknown; if (options.solver == kIpmString || options.run_centring) { // Setting the IPM-specific values of (highs_)info_ has been done in // solveLpIpx From 1210410ed55d4965c9e60608b6d487f7082ed8e4 Mon Sep 17 00:00:00 2001 From: JAJHall Date: Tue, 15 Oct 2024 22:06:44 +0100 Subject: [PATCH 5/7] WIP --- check/TestCallbacks.cpp | 2 +- check/TestIpm.cpp | 2 +- check/TestPdlp.cpp | 2 +- check/TestSpecialLps.cpp | 2 +- src/lp_data/HighsSolve.cpp | 9 +++++---- 5 files changed, 9 insertions(+), 8 deletions(-) diff --git a/check/TestCallbacks.cpp b/check/TestCallbacks.cpp index d97bc0d271..e6de40de4d 100644 --- a/check/TestCallbacks.cpp +++ b/check/TestCallbacks.cpp @@ -6,7 +6,7 @@ #include "catch.hpp" #include "lp_data/HighsCallback.h" -const bool dev_run = false; +const bool dev_run = true; const double egout_optimal_objective = 568.1007; const double egout_objective_target = 610; diff --git a/check/TestIpm.cpp b/check/TestIpm.cpp index ff60bbd403..651cb113ed 100644 --- a/check/TestIpm.cpp +++ b/check/TestIpm.cpp @@ -5,7 +5,7 @@ // I use dev_run to switch on/off printing and logging used for // development of the unit test -const bool dev_run = false; +const bool dev_run = true; const double inf = kHighsInf; TEST_CASE("test-analytic-centre", "[highs_ipm]") { diff --git a/check/TestPdlp.cpp b/check/TestPdlp.cpp index 0514e72dc5..93e75d2c46 100644 --- a/check/TestPdlp.cpp +++ b/check/TestPdlp.cpp @@ -3,7 +3,7 @@ #include "SpecialLps.h" #include "catch.hpp" -const bool dev_run = false; +const bool dev_run = true; const double double_equal_tolerance = 1e-3; TEST_CASE("pdlp-distillation-lp", "[pdlp]") { diff --git a/check/TestSpecialLps.cpp b/check/TestSpecialLps.cpp index 66d99ed304..c39470f356 100644 --- a/check/TestSpecialLps.cpp +++ b/check/TestSpecialLps.cpp @@ -3,7 +3,7 @@ #include "SpecialLps.h" #include "catch.hpp" -const bool dev_run = false; +const bool dev_run = true; void solve(Highs& highs, std::string presolve, std::string solver, const HighsModelStatus require_model_status, diff --git a/src/lp_data/HighsSolve.cpp b/src/lp_data/HighsSolve.cpp index 9e11c91888..536ca04748 100644 --- a/src/lp_data/HighsSolve.cpp +++ b/src/lp_data/HighsSolve.cpp @@ -72,10 +72,11 @@ HighsStatus solveLp(HighsLpSolverObject& solver_object, const string message) { return_status, "solveLpCupdlp"); } if (return_status == HighsStatus::kError) return return_status; - // IPM (and PDLP?) can claim optimality with large primal and/or - // dual residual errors, so correct any residual errors that - // exceed the tolerance - correctResiduals(solver_object); + if (solver_object.model_status_ == HighsModelStatus::kOptimal) + // IPM (and PDLP?) can claim optimality with large primal and/or + // dual residual errors, so correct any residual errors that + // exceed the tolerance + correctResiduals(solver_object); // Non-error return requires a primal solution assert(solver_object.solution_.value_valid); // Get the objective and any KKT failures From 42e520b053b9051d2307808cd43b800b4b108b11 Mon Sep 17 00:00:00 2001 From: JAJHall Date: Tue, 15 Oct 2024 22:40:20 +0100 Subject: [PATCH 6/7] Formatted --- check/TestCallbacks.cpp | 2 +- check/TestIpm.cpp | 19 +++++++------ check/TestLpSolvers.cpp | 7 +++-- check/TestPdlp.cpp | 2 +- check/TestSpecialLps.cpp | 2 +- src/lp_data/HighsOptions.h | 2 +- src/lp_data/HighsSolution.cpp | 52 +++++++++++++++-------------------- src/lp_data/HighsSolution.h | 1 - src/lp_data/HighsSolve.cpp | 6 ++-- 9 files changed, 43 insertions(+), 50 deletions(-) diff --git a/check/TestCallbacks.cpp b/check/TestCallbacks.cpp index e6de40de4d..d97bc0d271 100644 --- a/check/TestCallbacks.cpp +++ b/check/TestCallbacks.cpp @@ -6,7 +6,7 @@ #include "catch.hpp" #include "lp_data/HighsCallback.h" -const bool dev_run = true; +const bool dev_run = false; const double egout_optimal_objective = 568.1007; const double egout_objective_target = 610; diff --git a/check/TestIpm.cpp b/check/TestIpm.cpp index 651cb113ed..cb2d1a3419 100644 --- a/check/TestIpm.cpp +++ b/check/TestIpm.cpp @@ -5,7 +5,7 @@ // I use dev_run to switch on/off printing and logging used for // development of the unit test -const bool dev_run = true; +const bool dev_run = false; const double inf = kHighsInf; TEST_CASE("test-analytic-centre", "[highs_ipm]") { @@ -106,29 +106,30 @@ TEST_CASE("test-1966", "[highs_ipm]") { if (dev_run) { highs.writeSolution("", kSolutionStylePretty); printf("Num primal infeasibilities = %d\n", - int(info.num_primal_infeasibilities)); + int(info.num_primal_infeasibilities)); printf("Max primal infeasibility = %g\n", info.max_primal_infeasibility); - printf("Sum primal infeasibilities = %g\n", info.sum_primal_infeasibilities); + printf("Sum primal infeasibilities = %g\n", + info.sum_primal_infeasibilities); printf("Num dual infeasibilities = %d\n", - int(info.num_dual_infeasibilities)); + int(info.num_dual_infeasibilities)); printf("Max dual infeasibility = %g\n", info.max_dual_infeasibility); printf("Sum dual infeasibilities = %g\n", info.sum_dual_infeasibilities); } highs.clearSolver(); - if (dev_run) - printf("\nWith infinite residual tolerances\n"); + if (dev_run) printf("\nWith infinite residual tolerances\n"); highs.setOptionValue("primal_residual_tolerance", 1e30); highs.setOptionValue("dual_residual_tolerance", 1e30); highs.run(); if (dev_run) { highs.writeSolution("", kSolutionStylePretty); printf("Num primal infeasibilities = %d\n", - int(info.num_primal_infeasibilities)); + int(info.num_primal_infeasibilities)); printf("Max primal infeasibility = %g\n", info.max_primal_infeasibility); - printf("Sum primal infeasibilities = %g\n", info.sum_primal_infeasibilities); + printf("Sum primal infeasibilities = %g\n", + info.sum_primal_infeasibilities); printf("Num dual infeasibilities = %d\n", - int(info.num_dual_infeasibilities)); + int(info.num_dual_infeasibilities)); printf("Max dual infeasibility = %g\n", info.max_dual_infeasibility); printf("Sum dual infeasibilities = %g\n", info.sum_dual_infeasibilities); } diff --git a/check/TestLpSolvers.cpp b/check/TestLpSolvers.cpp index e4f8d4bc97..572a478cf2 100644 --- a/check/TestLpSolvers.cpp +++ b/check/TestLpSolvers.cpp @@ -460,11 +460,12 @@ TEST_CASE("blending-lp-ipm", "[highs_lp_solver]") { HighsInfo info = highs.getInfo(); if (dev_run) { printf("Num primal infeasibilities = %d\n", - int(info.num_primal_infeasibilities)); + int(info.num_primal_infeasibilities)); printf("Max primal infeasibility = %g\n", info.max_primal_infeasibility); - printf("Sum primal infeasibilities = %g\n", info.sum_primal_infeasibilities); + printf("Sum primal infeasibilities = %g\n", + info.sum_primal_infeasibilities); printf("Num dual infeasibilities = %d\n", - int(info.num_dual_infeasibilities)); + int(info.num_dual_infeasibilities)); printf("Max dual infeasibility = %g\n", info.max_dual_infeasibility); printf("Sum dual infeasibilities = %g\n", info.sum_dual_infeasibilities); } diff --git a/check/TestPdlp.cpp b/check/TestPdlp.cpp index 93e75d2c46..0514e72dc5 100644 --- a/check/TestPdlp.cpp +++ b/check/TestPdlp.cpp @@ -3,7 +3,7 @@ #include "SpecialLps.h" #include "catch.hpp" -const bool dev_run = true; +const bool dev_run = false; const double double_equal_tolerance = 1e-3; TEST_CASE("pdlp-distillation-lp", "[pdlp]") { diff --git a/check/TestSpecialLps.cpp b/check/TestSpecialLps.cpp index c39470f356..66d99ed304 100644 --- a/check/TestSpecialLps.cpp +++ b/check/TestSpecialLps.cpp @@ -3,7 +3,7 @@ #include "SpecialLps.h" #include "catch.hpp" -const bool dev_run = true; +const bool dev_run = false; void solve(Highs& highs, std::string presolve, std::string solver, const HighsModelStatus require_model_status, diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h index 9b38657551..bf853a4839 100644 --- a/src/lp_data/HighsOptions.h +++ b/src/lp_data/HighsOptions.h @@ -559,7 +559,7 @@ struct HighsOptionsStruct { #endif mip_improving_solution_save(false), mip_improving_solution_report_sparse(false), - mip_improving_solution_file(""){}; + mip_improving_solution_file("") {}; }; // For now, but later change so HiGHS properties are string based so that new diff --git a/src/lp_data/HighsSolution.cpp b/src/lp_data/HighsSolution.cpp index 6da9ba4ab4..27cd16153e 100644 --- a/src/lp_data/HighsSolution.cpp +++ b/src/lp_data/HighsSolution.cpp @@ -282,7 +282,8 @@ void getKktFailures(const HighsOptions& options, const HighsLp& lp, relative_primal_infeasibility, dual_infeasibility, value_residual); if (!status_value_ok) highsLogUser(options.log_options, HighsLogType::kError, - "getKktFailures: %s %d status-value error: [%23.18g; %23.18g; %23.18g] has " + "getKktFailures: %s %d status-value error: [%23.18g; " + "%23.18g; %23.18g] has " "residual %g\n", iVar < lp.num_col_ ? "Column" : "Row ", iVar < lp.num_col_ ? int(iVar) : int(iVar - lp.num_col_), @@ -539,14 +540,12 @@ bool getVariableKktFailures(const double primal_feasibility_tolerance, // bounds - for debugging QP basis errors if (*status_pointer == HighsBasisStatus::kLower) { if (std::fabs(lower) / primal_feasibility_tolerance < 1e-16) - status_value_ok = - value >= lower - primal_feasibility_tolerance && - value <= lower + primal_feasibility_tolerance; + status_value_ok = value >= lower - primal_feasibility_tolerance && + value <= lower + primal_feasibility_tolerance; } else if (*status_pointer == HighsBasisStatus::kUpper) { if (std::fabs(upper) / primal_feasibility_tolerance < 1e-16) - status_value_ok = - value >= upper - primal_feasibility_tolerance && - value <= upper + primal_feasibility_tolerance; + status_value_ok = value >= upper - primal_feasibility_tolerance && + value <= upper + primal_feasibility_tolerance; } } if (at_a_bound) { @@ -1390,7 +1389,7 @@ void correctResiduals(HighsLpSolverObject& solver_object) { lp.a_matrix_.productQuad(check_row_value, solution.col_value); if (have_dual_solution) { lp.a_matrix_.productTransposeQuad(check_col_dual, solution.row_dual); - for (HighsInt iCol=0; iCol < lp.num_col_; iCol++) + for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) check_col_dual[iCol] -= lp.col_cost_[iCol]; } double norm_primal_residual = 0; @@ -1408,14 +1407,14 @@ void correctResiduals(HighsLpSolverObject& solver_object) { for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++) { double primal_residual = check_row_value[iRow] - solution.row_value[iRow]; double abs_primal_residual = std::fabs(primal_residual); - // printf("correctResiduals: Row %d; check_row_value = %11.4g; row_value = %11.4g; primal_residual = %11.4g\n", - // int(iRow), check_row_value[iRow], solution.row_value[iRow], abs_primal_residual); if (abs_primal_residual > use_primal_residual_tolerance) { solution.row_value[iRow] += primal_residual; - double check_primal_residual = std::fabs(check_row_value[iRow] - solution.row_value[iRow]); + double check_primal_residual = + std::fabs(check_row_value[iRow] - solution.row_value[iRow]); assert(check_primal_residual < 1e-10); num_primal_correction++; - max_primal_correction = std::max(abs_primal_residual, max_primal_correction); + max_primal_correction = + std::max(abs_primal_residual, max_primal_correction); sum_primal_correction += abs_primal_residual; } norm_primal_residual = std::max(abs_primal_residual, norm_primal_residual); @@ -1424,31 +1423,24 @@ void correctResiduals(HighsLpSolverObject& solver_object) { for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) { double dual_residual = check_col_dual[iCol] + solution.col_dual[iCol]; double abs_dual_residual = std::fabs(dual_residual); - // printf("correctResiduals: Col %d; check_col_dual = %11.4g; col_dual = %11.4g; col_cost_ = %11.4g; dual_residual = %11.4g\n", - // int(iCol), check_col_dual[iCol], solution.col_dual[iCol], lp.col_cost_[iCol], abs_dual_residual); if (abs_dual_residual > use_dual_residual_tolerance) { - solution.col_dual[iCol] -= dual_residual; - double check_dual_residual = std::fabs(check_col_dual[iCol] + solution.col_dual[iCol]); - const bool ok_dual_correction = check_dual_residual < 1e-10; - // if (!ok_dual_correction) { - // printf("correctResiduals: Col %d; check_col_dual = %11.4g; col_dual = %11.4g; check_dual_residual = %g\n", - // int(iCol), check_col_dual[iCol], solution.col_dual[iCol], check_dual_residual); - // } - assert(ok_dual_correction); - num_dual_correction++; - max_dual_correction = std::max(abs_dual_residual, max_dual_correction); - sum_dual_correction += abs_dual_residual; + solution.col_dual[iCol] -= dual_residual; + num_dual_correction++; + max_dual_correction = std::max(abs_dual_residual, max_dual_correction); + sum_dual_correction += abs_dual_residual; } norm_dual_residual = std::max(abs_dual_residual, norm_dual_residual); } } if (num_primal_correction > 0 || num_dual_correction > 0) - highsLogUser(options.log_options, HighsLogType::kWarning, - "LP solver yields |primal_residual| = %g; |dual_residual| = %g with num / max / sum primal (%d / %g / %g) and dual (%d / %g / %g) corrections\n", - norm_primal_residual, norm_dual_residual, - int(num_primal_correction), max_primal_correction, sum_primal_correction, - int(num_dual_correction), max_dual_correction, sum_dual_correction); + highsLogUser( + options.log_options, HighsLogType::kWarning, + "LP solver residuals: primal = %g; dual = %g yield num/max/sum primal " + "(%d/%g/%g) and dual (%d/%g/%g) corrections\n", + norm_primal_residual, norm_dual_residual, int(num_primal_correction), + max_primal_correction, sum_primal_correction, int(num_dual_correction), + max_dual_correction, sum_dual_correction); } void resetModelStatusAndHighsInfo(HighsLpSolverObject& solver_object) { diff --git a/src/lp_data/HighsSolution.h b/src/lp_data/HighsSolution.h index 1efbedda52..36342274eb 100644 --- a/src/lp_data/HighsSolution.h +++ b/src/lp_data/HighsSolution.h @@ -136,5 +136,4 @@ bool isDualSolutionRightSize(const HighsLp& lp, const HighsSolution& solution); bool isSolutionRightSize(const HighsLp& lp, const HighsSolution& solution); bool isBasisRightSize(const HighsLp& lp, const HighsBasis& basis); - #endif // LP_DATA_HIGHSSOLUTION_H_ diff --git a/src/lp_data/HighsSolve.cpp b/src/lp_data/HighsSolve.cpp index 536ca04748..9c013b4f2e 100644 --- a/src/lp_data/HighsSolve.cpp +++ b/src/lp_data/HighsSolve.cpp @@ -85,9 +85,9 @@ HighsStatus solveLp(HighsLpSolverObject& solver_object, const string message) { getLpKktFailures(options, solver_object.lp_, solver_object.solution_, solver_object.basis_, solver_object.highs_info_); if (solver_object.model_status_ == HighsModelStatus::kOptimal && - solver_object.highs_info_.num_primal_infeasibilities > 0 || - solver_object.highs_info_.num_dual_infeasibilities) - solver_object.model_status_ = HighsModelStatus::kUnknown; + (solver_object.highs_info_.num_primal_infeasibilities > 0 || + solver_object.highs_info_.num_dual_infeasibilities)) + solver_object.model_status_ = HighsModelStatus::kUnknown; if (options.solver == kIpmString || options.run_centring) { // Setting the IPM-specific values of (highs_)info_ has been done in // solveLpIpx From 7efad63936ca6e1ec4600118917b81e4f6eb2a6c Mon Sep 17 00:00:00 2001 From: JAJHall Date: Tue, 15 Oct 2024 22:48:18 +0100 Subject: [PATCH 7/7] Correct formatting? --- src/util/HFactor.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/util/HFactor.h b/src/util/HFactor.h index 2f75260d8a..03cb8d8b1d 100644 --- a/src/util/HFactor.h +++ b/src/util/HFactor.h @@ -132,7 +132,7 @@ class HFactor { build_timer_(nullptr), nwork(0), u_merit_x(0), - u_total_x(0){}; + u_total_x(0) {}; /** * @brief Copy problem size and pointers of constraint matrix, and set