diff --git a/.github/workflows/clang-format.yml b/.github/workflows/clang-format.yml index b621d680e7..8e368317c8 100644 --- a/.github/workflows/clang-format.yml +++ b/.github/workflows/clang-format.yml @@ -10,7 +10,8 @@ jobs: - uses: actions/checkout@v4 - uses: DoozyX/clang-format-lint-action@v0.18 with: - source: 'app/ src/Highs.h ./src/lp_data ./src/mip ./src/model ./src/simplex ./src/presolve ./src/simplex ./src/util ./src/test' - #./src/test ./interfaces' + source: + 'app/ src/Highs.h ./src/lp_data ./src/mip ./src/model ./src/simplex ./src/presolve ./src/simplex ./src/util ./src/test ./src/qpsolver' + # ./src/test ./interfaces' extensions: 'h,cpp,c' clangFormatVersion: 18 diff --git a/src/qpsolver/a_asm.cpp b/src/qpsolver/a_asm.cpp index 821058619e..4fb9706c4e 100644 --- a/src/qpsolver/a_asm.cpp +++ b/src/qpsolver/a_asm.cpp @@ -1,8 +1,12 @@ #include "qpsolver/a_asm.hpp" + #include "qpsolver/quass.hpp" #include "util/HighsCDouble.h" -QpAsmStatus solveqp_actual(Instance& instance, Settings& settings, QpHotstartInformation& startinfo, Statistics& stats, QpModelStatus& status, QpSolution& solution, HighsTimer& qp_timer) { +QpAsmStatus solveqp_actual(Instance& instance, Settings& settings, + QpHotstartInformation& startinfo, Statistics& stats, + QpModelStatus& status, QpSolution& solution, + HighsTimer& qp_timer) { Runtime rt(instance, stats); rt.settings = settings; Quass quass(rt); @@ -21,54 +25,55 @@ QpAsmStatus solveqp_actual(Instance& instance, Settings& settings, QpHotstartInf solution.dualcon = rt.dualcon; return QpAsmStatus::kOk; - } std::string qpBasisStatusToString(const BasisStatus qp_basis_status) { switch (qp_basis_status) { - case BasisStatus::kInactive: - return "Inactive"; - case BasisStatus::kActiveAtLower: - return "Active at lower bound"; - case BasisStatus::kActiveAtUpper: - return "Active at upper bound"; - case BasisStatus::kInactiveInBasis: - return "Inactive in basis"; - default: - return "Unidentified QP basis status"; + case BasisStatus::kInactive: + return "Inactive"; + case BasisStatus::kActiveAtLower: + return "Active at lower bound"; + case BasisStatus::kActiveAtUpper: + return "Active at upper bound"; + case BasisStatus::kInactiveInBasis: + return "Inactive in basis"; + default: + return "Unidentified QP basis status"; } } std::string qpModelStatusToString(const QpModelStatus qp_model_status) { switch (qp_model_status) { - case QpModelStatus::kNotset: - return "Not set"; - case QpModelStatus::kUndetermined: - return "Undertermined"; - case QpModelStatus::kOptimal: - return "Optimal"; - case QpModelStatus::kUnbounded: - return "Unbounded"; - case QpModelStatus::kInfeasible: - return "Infeasible"; - case QpModelStatus::kIterationLimit: - return "Iteration limit"; - case QpModelStatus::kTimeLimit: - return "Time ;limit"; - case QpModelStatus::kLargeNullspace: - return "Large nullspace"; - case QpModelStatus::kError: - return "Error"; - default: - return "Unidentified QP model status"; + case QpModelStatus::kNotset: + return "Not set"; + case QpModelStatus::kUndetermined: + return "Undertermined"; + case QpModelStatus::kOptimal: + return "Optimal"; + case QpModelStatus::kUnbounded: + return "Unbounded"; + case QpModelStatus::kInfeasible: + return "Infeasible"; + case QpModelStatus::kIterationLimit: + return "Iteration limit"; + case QpModelStatus::kTimeLimit: + return "Time ;limit"; + case QpModelStatus::kLargeNullspace: + return "Large nullspace"; + case QpModelStatus::kError: + return "Error"; + default: + return "Unidentified QP model status"; } } -void assessQpPrimalFeasibility(const Instance& instance, const double primal_feasibility_tolerance, - const std::vector& var_value, const std::vector& con_value, - HighsInt& num_var_infeasibilities, double& max_var_infeasibility, double& sum_var_infeasibilities, - HighsInt& num_con_infeasibilities, double& max_con_infeasibility, double& sum_con_infeasibilities, - double& max_con_residual, double& sum_con_residuals) { +void assessQpPrimalFeasibility( + const Instance& instance, const double primal_feasibility_tolerance, + const std::vector& var_value, const std::vector& con_value, + HighsInt& num_var_infeasibilities, double& max_var_infeasibility, + double& sum_var_infeasibilities, HighsInt& num_con_infeasibilities, + double& max_con_infeasibility, double& sum_con_infeasibilities, + double& max_con_residual, double& sum_con_residuals) { num_var_infeasibilities = 0; max_var_infeasibility = 0; sum_var_infeasibilities = 0; @@ -91,14 +96,16 @@ void assessQpPrimalFeasibility(const Instance& instance, const double primal_fea var_infeasibility = primal - upper; } if (var_infeasibility > 0) { - if (var_infeasibility > primal_feasibility_tolerance) - num_var_infeasibilities++; + if (var_infeasibility > primal_feasibility_tolerance) + num_var_infeasibilities++; max_var_infeasibility = - std::max(var_infeasibility, max_var_infeasibility); + std::max(var_infeasibility, max_var_infeasibility); sum_var_infeasibilities += var_infeasibility; } - for (HighsInt iEl = instance.A.mat.start[iVar]; iEl < instance.A.mat.start[iVar+1]; iEl++) { - con_value_quad[instance.A.mat.index[iEl]] += primal * instance.A.mat.value[iEl]; + for (HighsInt iEl = instance.A.mat.start[iVar]; + iEl < instance.A.mat.start[iVar + 1]; iEl++) { + con_value_quad[instance.A.mat.index[iEl]] += + primal * instance.A.mat.value[iEl]; } } for (HighsInt iCon = 0; iCon < instance.num_con; iCon++) { @@ -112,15 +119,14 @@ void assessQpPrimalFeasibility(const Instance& instance, const double primal_fea con_infeasibility = primal - upper; } if (con_infeasibility > 0) { - if (con_infeasibility > primal_feasibility_tolerance) - num_con_infeasibilities++; + if (con_infeasibility > primal_feasibility_tolerance) + num_con_infeasibilities++; max_con_infeasibility = - std::max(con_infeasibility, max_con_infeasibility); + std::max(con_infeasibility, max_con_infeasibility); sum_con_infeasibilities += con_infeasibility; } - double con_residual = std::fabs(primal-double(con_value_quad[iCon])); + double con_residual = std::fabs(primal - double(con_value_quad[iCon])); max_con_residual = std::max(con_residual, max_con_residual); sum_con_residuals += con_residual; } } - diff --git a/src/qpsolver/a_asm.hpp b/src/qpsolver/a_asm.hpp index 5d80757bc7..7cb9938840 100644 --- a/src/qpsolver/a_asm.hpp +++ b/src/qpsolver/a_asm.hpp @@ -2,9 +2,9 @@ #define __SRC_LIB_QPSOLVER_ASM_HPP__ #include "qpsolver/instance.hpp" -#include "qpsolver/statistics.hpp" #include "qpsolver/qpconst.hpp" #include "qpsolver/settings.hpp" +#include "qpsolver/statistics.hpp" #include "util/HighsTimer.h" enum class QpAsmStatus { @@ -24,7 +24,8 @@ struct QpSolution { std::vector status_var; std::vector status_con; - QpSolution(Instance& instance) : primal(QpVector(instance.num_var)), + QpSolution(Instance& instance) + : primal(QpVector(instance.num_var)), rowactivity(QpVector(instance.num_con)), dualvar(instance.num_var), dualcon(instance.num_con), @@ -32,7 +33,6 @@ struct QpSolution { status_con(instance.num_con) {} }; - struct QpHotstartInformation { std::vector active; std::vector inactive; @@ -44,24 +44,27 @@ struct QpHotstartInformation { : primal(QpVector(num_var)), rowact(QpVector(num_row)) {} }; -// the purpose of this is the pure algorithmic solution of a QP instance with given hotstart information. -// scenarios: -// 1) start from a given phase1 solution +// the purpose of this is the pure algorithmic solution of a QP instance with +// given hotstart information. scenarios: 1) start from a given phase1 solution // 2) start from a user-given hotstart solution -// 3) start from a qp solution that was attained from a scaled instance and cleanup -// 4) start from a qp solution that was attained from a perturbed instance and cleanup -// 5) start from a qp solution and cleanup after recomputing basis and reduced hessian factorization +// 3) start from a qp solution that was attained from a scaled instance and +// cleanup 4) start from a qp solution that was attained from a perturbed +// instance and cleanup 5) start from a qp solution and cleanup after +// recomputing basis and reduced hessian factorization std::string qpBasisStatusToString(const BasisStatus qp_basis_status); std::string qpModelStatusToString(const QpModelStatus qp_model_status); -void assessQpPrimalFeasibility(const Instance& instance, const double primal_feasibility_tolerance, - const std::vector& var_value, const std::vector& con_value, - HighsInt& num_var_infeasibilities, double& max_var_infeasibility, double& sum_var_infeasibilities, - HighsInt& num_con_infeasibilities, double& max_con_infeasibility, double& sum_con_infeasibilities, - double& max_con_residual, double& sum_con_residuals); - -QpAsmStatus solveqp_actual(Instance& instance, Settings& settings, QpHotstartInformation& startinfo, Statistics& stats, QpModelStatus& status, QpSolution& solution, HighsTimer& qp_timer); - +void assessQpPrimalFeasibility( + const Instance& instance, const double primal_feasibility_tolerance, + const std::vector& var_value, const std::vector& con_value, + HighsInt& num_var_infeasibilities, double& max_var_infeasibility, + double& sum_var_infeasibilities, HighsInt& num_con_infeasibilities, + double& max_con_infeasibility, double& sum_con_infeasibilities, + double& max_con_residual, double& sum_con_residuals); +QpAsmStatus solveqp_actual(Instance& instance, Settings& settings, + QpHotstartInformation& startinfo, Statistics& stats, + QpModelStatus& status, QpSolution& solution, + HighsTimer& qp_timer); #endif diff --git a/src/qpsolver/a_quass.cpp b/src/qpsolver/a_quass.cpp index b62d22ebc9..b639991eda 100644 --- a/src/qpsolver/a_quass.cpp +++ b/src/qpsolver/a_quass.cpp @@ -1,56 +1,55 @@ #include "qpsolver/a_quass.hpp" -#include "qpsolver/a_asm.hpp" -#include "qpsolver/feasibility_highs.hpp" +#include "qpsolver/a_asm.hpp" #include "qpsolver/feasibility_bounded.hpp" +#include "qpsolver/feasibility_highs.hpp" -static QpAsmStatus quass2highs(Instance& instance, - Settings& settings, - Statistics& stats, - QpModelStatus& qp_model_status, - QpSolution& qp_solution, - HighsModelStatus& highs_model_status, - HighsBasis& highs_basis, - HighsSolution& highs_solution) { +static QpAsmStatus quass2highs(Instance& instance, Settings& settings, + Statistics& stats, + QpModelStatus& qp_model_status, + QpSolution& qp_solution, + HighsModelStatus& highs_model_status, + HighsBasis& highs_basis, + HighsSolution& highs_solution) { settings.qp_model_status_log.fire(qp_model_status); QpAsmStatus qp_asm_return_status = QpAsmStatus::kError; switch (qp_model_status) { - case QpModelStatus::kOptimal: - highs_model_status = HighsModelStatus::kOptimal; - qp_asm_return_status = QpAsmStatus::kOk; - break; - case QpModelStatus::kUnbounded: - highs_model_status = HighsModelStatus::kUnbounded; - qp_asm_return_status = QpAsmStatus::kOk; - break; - case QpModelStatus::kInfeasible: - highs_model_status = HighsModelStatus::kInfeasible; - qp_asm_return_status = QpAsmStatus::kOk; - break; - case QpModelStatus::kIterationLimit: - highs_model_status = HighsModelStatus::kIterationLimit; - qp_asm_return_status = QpAsmStatus::kWarning; - break; - case QpModelStatus::kTimeLimit: - highs_model_status = HighsModelStatus::kTimeLimit; - qp_asm_return_status = QpAsmStatus::kWarning; - break; - case QpModelStatus::kUndetermined: - highs_model_status = HighsModelStatus::kSolveError; - qp_asm_return_status = QpAsmStatus::kError; - return QpAsmStatus::kError; - case QpModelStatus::kLargeNullspace: - highs_model_status = HighsModelStatus::kSolveError; - return QpAsmStatus::kError; - case QpModelStatus::kError: - highs_model_status = HighsModelStatus::kSolveError; - return QpAsmStatus::kError; - case QpModelStatus::kNotset: - highs_model_status = HighsModelStatus::kNotset; - return QpAsmStatus::kError; - default: - highs_model_status = HighsModelStatus::kNotset; - return QpAsmStatus::kError; + case QpModelStatus::kOptimal: + highs_model_status = HighsModelStatus::kOptimal; + qp_asm_return_status = QpAsmStatus::kOk; + break; + case QpModelStatus::kUnbounded: + highs_model_status = HighsModelStatus::kUnbounded; + qp_asm_return_status = QpAsmStatus::kOk; + break; + case QpModelStatus::kInfeasible: + highs_model_status = HighsModelStatus::kInfeasible; + qp_asm_return_status = QpAsmStatus::kOk; + break; + case QpModelStatus::kIterationLimit: + highs_model_status = HighsModelStatus::kIterationLimit; + qp_asm_return_status = QpAsmStatus::kWarning; + break; + case QpModelStatus::kTimeLimit: + highs_model_status = HighsModelStatus::kTimeLimit; + qp_asm_return_status = QpAsmStatus::kWarning; + break; + case QpModelStatus::kUndetermined: + highs_model_status = HighsModelStatus::kSolveError; + qp_asm_return_status = QpAsmStatus::kError; + return QpAsmStatus::kError; + case QpModelStatus::kLargeNullspace: + highs_model_status = HighsModelStatus::kSolveError; + return QpAsmStatus::kError; + case QpModelStatus::kError: + highs_model_status = HighsModelStatus::kSolveError; + return QpAsmStatus::kError; + case QpModelStatus::kNotset: + highs_model_status = HighsModelStatus::kNotset; + return QpAsmStatus::kError; + default: + highs_model_status = HighsModelStatus::kNotset; + return QpAsmStatus::kError; } assert(qp_asm_return_status != QpAsmStatus::kError); @@ -59,7 +58,8 @@ static QpAsmStatus quass2highs(Instance& instance, highs_solution.col_dual.resize(instance.num_var); for (HighsInt iCol = 0; iCol < instance.num_var; iCol++) { highs_solution.col_value[iCol] = qp_solution.primal.value[iCol]; - highs_solution.col_dual[iCol] = instance.sense * qp_solution.dualvar.value[iCol]; + highs_solution.col_dual[iCol] = + instance.sense * qp_solution.dualvar.value[iCol]; } // extract constraint activity highs_solution.row_value.resize(instance.num_con); @@ -67,7 +67,8 @@ static QpAsmStatus quass2highs(Instance& instance, // Negate the vector and Hessian for (HighsInt iRow = 0; iRow < instance.num_con; iRow++) { highs_solution.row_value[iRow] = qp_solution.rowactivity.value[iRow]; - highs_solution.row_dual[iRow] = instance.sense * qp_solution.dualcon.value[iRow]; + highs_solution.row_dual[iRow] = + instance.sense * qp_solution.dualcon.value[iRow]; } highs_solution.value_valid = true; highs_solution.dual_valid = true; @@ -76,9 +77,12 @@ static QpAsmStatus quass2highs(Instance& instance, highs_basis.col_status.resize(instance.num_var); highs_basis.row_status.resize(instance.num_con); - const bool debug_report = false;// instance.num_var + instance.num_con < 100; + const bool debug_report = + false; // instance.num_var + instance.num_con < 100; for (HighsInt i = 0; i < instance.num_var; i++) { - if (debug_report) printf("Column %2d: status %s\n", int(i), qpBasisStatusToString(qp_solution.status_var[i]).c_str()); + if (debug_report) + printf("Column %2d: status %s\n", int(i), + qpBasisStatusToString(qp_solution.status_var[i]).c_str()); if (qp_solution.status_var[i] == BasisStatus::kActiveAtLower) { highs_basis.col_status[i] = HighsBasisStatus::kLower; } else if (qp_solution.status_var[i] == BasisStatus::kActiveAtUpper) { @@ -92,7 +96,9 @@ static QpAsmStatus quass2highs(Instance& instance, } for (HighsInt i = 0; i < instance.num_con; i++) { - if (debug_report) printf("Row %2d: status %s\n", int(i), qpBasisStatusToString(qp_solution.status_con[i]).c_str()); + if (debug_report) + printf("Row %2d: status %s\n", int(i), + qpBasisStatusToString(qp_solution.status_con[i]).c_str()); if (qp_solution.status_con[i] == BasisStatus::kActiveAtLower) { highs_basis.row_status[i] = HighsBasisStatus::kLower; } else if (qp_solution.status_con[i] == BasisStatus::kActiveAtUpper) { @@ -109,14 +115,10 @@ static QpAsmStatus quass2highs(Instance& instance, return qp_asm_return_status; } -QpAsmStatus solveqp(Instance& instance, - Settings& settings, - Statistics& stats, - HighsModelStatus& highs_model_status, - HighsBasis& highs_basis, - HighsSolution& highs_solution, - HighsTimer& qp_timer) { - +QpAsmStatus solveqp(Instance& instance, Settings& settings, Statistics& stats, + HighsModelStatus& highs_model_status, + HighsBasis& highs_basis, HighsSolution& highs_solution, + HighsTimer& qp_timer) { QpModelStatus qp_model_status = QpModelStatus::kUndetermined; QpSolution qp_solution(instance); @@ -128,12 +130,11 @@ QpAsmStatus solveqp(Instance& instance, // perturb instance, store perturbance information // regularize - for (HighsInt i=0; i active, buffer_row_ep(rt.instance.num_var) { buffer_vec2hvec.setup(rt.instance.num_var); - for (HighsInt i=0; i active, active_constraint_index.push_back(active[i]); basisstatus[active_constraint_index[i]] = status[i]; } - for (size_t i = 0; i< inactive.size(); i++) { + for (size_t i = 0; i < inactive.size(); i++) { non_active_constraint_index.push_back(inactive[i]); basisstatus[non_active_constraint_index[i]] = BasisStatus::kInactiveInBasis; } @@ -38,14 +39,15 @@ void Basis::build() { updatessinceinvert = 0; - baseindex.resize(active_constraint_index.size() + non_active_constraint_index.size()); + baseindex.resize(active_constraint_index.size() + + non_active_constraint_index.size()); constraintindexinbasisfactor.clear(); basisfactor = HFactor(); constraintindexinbasisfactor.assign(Atran.num_row + Atran.num_col, -1); - assert((HighsInt)(non_active_constraint_index.size() + active_constraint_index.size()) == - Atran.num_row); + assert((HighsInt)(non_active_constraint_index.size() + + active_constraint_index.size()) == Atran.num_row); HighsInt counter = 0; for (HighsInt i : non_active_constraint_index) { @@ -66,12 +68,12 @@ void Basis::build() { Atran.value.resize(1); } basisfactor.setup(Atran.num_col, Atran.num_row, Atran.start.data(), - Atran.index.data(), Atran.value.data(), - baseindex.data()); + Atran.index.data(), Atran.value.data(), baseindex.data()); basisfactor.build(); for (size_t i = 0; - i < active_constraint_index.size() + non_active_constraint_index.size(); i++) { + i < active_constraint_index.size() + non_active_constraint_index.size(); + i++) { constraintindexinbasisfactor[baseindex[i]] = i; } } @@ -83,13 +85,14 @@ void Basis::rebuild() { constraintindexinbasisfactor.clear(); constraintindexinbasisfactor.assign(Atran.num_row + Atran.num_col, -1); - assert((HighsInt)(non_active_constraint_index.size() + active_constraint_index.size()) == - Atran.num_row); + assert((HighsInt)(non_active_constraint_index.size() + + active_constraint_index.size()) == Atran.num_row); basisfactor.build(); for (size_t i = 0; - i < active_constraint_index.size() + non_active_constraint_index.size(); i++) { + i < active_constraint_index.size() + non_active_constraint_index.size(); + i++) { constraintindexinbasisfactor[baseindex[i]] = i; } reinversion_hint = false; @@ -112,7 +115,7 @@ void Basis::report() { // // Remaining qp_num_con indices may be degenerate, otherwise they // are off their bounds. They are analogous to primal simplex basic - // variables, in that their values are solved for. + // variables, in that their values are solved for. // // Hence the correspondence between the QP basis and a HiGHS // (simplex) basis is as follows @@ -127,7 +130,7 @@ void Basis::report() { // // BasisStatus::kInactiveInBasis: HighsBasisStatus::kNonbasic // - // + // const HighsInt qp_num_var = Atran.num_row; const HighsInt qp_num_con = Atran.num_col; const HighsInt num_active_in_basis = active_constraint_index.size(); @@ -144,39 +147,39 @@ void Basis::report() { for (HighsInt i = 0; i < qp_num_var; i++) { switch (basisstatus[qp_num_con + i]) { - case BasisStatus::kInactive: - num_var_inactive++; - continue; - case BasisStatus::kActiveAtLower: - num_var_active_at_lower++; - continue; - case BasisStatus::kActiveAtUpper: - num_var_active_at_upper++; - continue; - case BasisStatus::kInactiveInBasis: - num_var_inactive_in_basis++; - continue; - default: - assert(111==123); + case BasisStatus::kInactive: + num_var_inactive++; + continue; + case BasisStatus::kActiveAtLower: + num_var_active_at_lower++; + continue; + case BasisStatus::kActiveAtUpper: + num_var_active_at_upper++; + continue; + case BasisStatus::kInactiveInBasis: + num_var_inactive_in_basis++; + continue; + default: + assert(111 == 123); } } for (HighsInt i = 0; i < qp_num_con; i++) { switch (basisstatus[i]) { - case BasisStatus::kInactive: - num_con_inactive++; - continue; - case BasisStatus::kActiveAtLower: - num_con_active_at_lower++; - continue; - case BasisStatus::kActiveAtUpper: - num_con_active_at_upper++; - continue; - case BasisStatus::kInactiveInBasis: - num_con_inactive_in_basis++; - continue; - default: - assert(111==123); + case BasisStatus::kInactive: + num_con_inactive++; + continue; + case BasisStatus::kActiveAtLower: + num_con_active_at_lower++; + continue; + case BasisStatus::kActiveAtUpper: + num_con_active_at_upper++; + continue; + case BasisStatus::kInactiveInBasis: + num_con_inactive_in_basis++; + continue; + default: + assert(111 == 123); } } @@ -185,38 +188,38 @@ void Basis::report() { printf("basis: "); for (HighsInt a_ : active_constraint_index) { if (a_ < qp_num_con) { - printf("c%-3d ", int(a_)); + printf("c%-3d ", int(a_)); } else { - printf("v%-3d ", int(a_-qp_num_con)); + printf("v%-3d ", int(a_ - qp_num_con)); } - } + } printf(" - "); for (HighsInt n_ : non_active_constraint_index) { if (n_ < qp_num_con) { - printf("c%-3d ", int(n_)); + printf("c%-3d ", int(n_)); } else { - printf("v%-3d ", int(n_-qp_num_con)); + printf("v%-3d ", int(n_ - qp_num_con)); } - } + } printf("\n"); } - printf("Basis::report: QP(%6d [inact %6d; act %6d], %6d)", - int(qp_num_var), int(num_inactive_in_basis), int(num_active_in_basis), - int(qp_num_con)); - printf(" (inact / lo / up / basis) for var (%6d / %6d / %6d / %6d) and con (%6d / %6d / %6d / %6d)\n", - int(num_var_inactive), - int(num_var_active_at_lower), - int(num_var_active_at_upper), - int(num_var_inactive_in_basis), - int(num_con_inactive), - int(num_con_active_at_lower), - int(num_con_active_at_upper), - int(num_con_inactive_in_basis)); + printf("Basis::report: QP(%6d [inact %6d; act %6d], %6d)", int(qp_num_var), + int(num_inactive_in_basis), int(num_active_in_basis), int(qp_num_con)); + printf( + " (inact / lo / up / basis) for var (%6d / %6d / %6d / %6d) and con (%6d " + "/ %6d / %6d / %6d)\n", + int(num_var_inactive), int(num_var_active_at_lower), + int(num_var_active_at_upper), int(num_var_inactive_in_basis), + int(num_con_inactive), int(num_con_active_at_lower), + int(num_con_active_at_upper), int(num_con_inactive_in_basis)); assert(qp_num_var == num_inactive_in_basis + num_active_in_basis); assert(qp_num_con == num_var_inactive + num_con_inactive); - assert(num_inactive_in_basis == num_var_inactive_in_basis + num_con_inactive_in_basis); - assert(num_active_in_basis == num_var_active_at_lower + num_var_active_at_upper + num_con_active_at_lower + num_con_active_at_upper); + assert(num_inactive_in_basis == + num_var_inactive_in_basis + num_con_inactive_in_basis); + assert(num_active_in_basis == + num_var_active_at_lower + num_var_active_at_upper + + num_con_active_at_lower + num_con_active_at_upper); } // move that constraint into V section basis (will correspond to Nullspace @@ -229,7 +232,8 @@ void Basis::deactivate(HighsInt conid) { non_active_constraint_index.push_back(conid); } -QpSolverStatus Basis::activate(const Settings& settings, HighsInt conid, BasisStatus newstatus, +QpSolverStatus Basis::activate(const Settings& settings, HighsInt conid, + BasisStatus newstatus, HighsInt nonactivetoremove, Pricing* pricing) { // printf("activ %" HIGHSINT_FORMAT "\n", conid); if (!contains(active_constraint_index, (HighsInt)conid)) { @@ -257,8 +261,8 @@ QpSolverStatus Basis::activate(const Settings& settings, HighsInt conid, BasisSt return QpSolverStatus::OK; } -void Basis::updatebasis(const Settings& settings, HighsInt newactivecon, HighsInt droppedcon, - Pricing* pricing) { +void Basis::updatebasis(const Settings& settings, HighsInt newactivecon, + HighsInt droppedcon, Pricing* pricing) { if (newactivecon == droppedcon) { return; } @@ -283,7 +287,8 @@ void Basis::updatebasis(const Settings& settings, HighsInt newactivecon, HighsIn basisfactor.update(&col_aq, &row_ep, &row_out, &hint); updatessinceinvert++; - if (updatessinceinvert >= settings.reinvertfrequency || hint != kHintNotChanged) { + if (updatessinceinvert >= settings.reinvertfrequency || + hint != kHintNotChanged) { reinversion_hint = true; } // since basis changed, buffered values are no longer valid @@ -292,7 +297,7 @@ void Basis::updatebasis(const Settings& settings, HighsInt newactivecon, HighsIn } QpVector& Basis::btran(const QpVector& rhs, QpVector& target, bool buffer, - HighsInt p) { + HighsInt p) { HVector rhs_hvec = vec2hvec(rhs); basisfactor.btranCall(rhs_hvec, 1.0); if (buffer) { @@ -325,7 +330,7 @@ QpVector Basis::btran(const QpVector& rhs, bool buffer, HighsInt p) { } QpVector& Basis::ftran(const QpVector& rhs, QpVector& target, bool buffer, - HighsInt q) { + HighsInt q) { HVector rhs_hvec = vec2hvec(rhs); basisfactor.ftranCall(rhs_hvec, 1.0); if (buffer) { @@ -394,7 +399,7 @@ QpVector Basis::recomputex(const Instance& inst) { } QpVector& Basis::Ztprod(const QpVector& rhs, QpVector& target, bool buffer, - HighsInt q) { + HighsInt q) { ftran(rhs, Ztprod_res, buffer, q); target.reset(); diff --git a/src/qpsolver/basis.hpp b/src/qpsolver/basis.hpp index c92d0aba40..0bdaca6376 100644 --- a/src/qpsolver/basis.hpp +++ b/src/qpsolver/basis.hpp @@ -85,9 +85,8 @@ class Basis { HVector col_aq; bool reinversion_hint = false; - public: - + public: Basis(Runtime& rt, std::vector active, std::vector atlower, std::vector inactive); @@ -99,7 +98,9 @@ class Basis { HighsInt getnumactive() const { return active_constraint_index.size(); }; - HighsInt getnuminactive() const { return non_active_constraint_index.size(); }; + HighsInt getnuminactive() const { + return non_active_constraint_index.size(); + }; const std::vector& getactive() const { return active_constraint_index; @@ -121,28 +122,29 @@ class Basis { // Nullspace from now on) void deactivate(HighsInt conid); - QpSolverStatus activate(const Settings& settings, HighsInt conid, BasisStatus atlower, - HighsInt nonactivetoremove, Pricing* pricing); + QpSolverStatus activate(const Settings& settings, HighsInt conid, + BasisStatus atlower, HighsInt nonactivetoremove, + Pricing* pricing); - void updatebasis(const Settings& settings, HighsInt newactivecon, HighsInt droppedcon, - Pricing* pricing); + void updatebasis(const Settings& settings, HighsInt newactivecon, + HighsInt droppedcon, Pricing* pricing); QpVector btran(const QpVector& rhs, bool buffer = false, HighsInt p = -1); QpVector ftran(const QpVector& rhs, bool buffer = false, HighsInt q = -1); QpVector& btran(const QpVector& rhs, QpVector& target, bool buffer = false, - HighsInt p = -1); + HighsInt p = -1); QpVector& ftran(const QpVector& rhs, QpVector& target, bool buffer = false, - HighsInt q = -1); + HighsInt q = -1); QpVector recomputex(const Instance& inst); void write(std::string filename); QpVector& Ztprod(const QpVector& rhs, QpVector& target, bool buffer = false, - HighsInt q = -1); + HighsInt q = -1); QpVector& Zprod(const QpVector& rhs, QpVector& target); }; diff --git a/src/qpsolver/crashsolution.hpp b/src/qpsolver/crashsolution.hpp index a2bc61eff3..31ac7b8f47 100644 --- a/src/qpsolver/crashsolution.hpp +++ b/src/qpsolver/crashsolution.hpp @@ -2,11 +2,12 @@ #define __SRC_LIB_CRASHSOLUTION_HPP__ #include + #include "runtime.hpp" -inline -bool isfreevar(Instance& instance, HighsInt idx) { - return instance.var_lo[idx] == -std::numeric_limits::infinity() && instance.var_up[idx] == std::numeric_limits::infinity(); +inline bool isfreevar(Instance& instance, HighsInt idx) { + return instance.var_lo[idx] == -std::numeric_limits::infinity() && + instance.var_up[idx] == std::numeric_limits::infinity(); } #endif diff --git a/src/qpsolver/dantzigpricing.hpp b/src/qpsolver/dantzigpricing.hpp index 44d06542f8..89b109e5c1 100644 --- a/src/qpsolver/dantzigpricing.hpp +++ b/src/qpsolver/dantzigpricing.hpp @@ -52,7 +52,7 @@ class DantzigPricing : public Pricing { public: DantzigPricing(Runtime& rt, Basis& bas, ReducedCosts& rc) - : runtime(rt), basis(bas), redcosts(rc){}; + : runtime(rt), basis(bas), redcosts(rc) {}; HighsInt price(const QpVector& x, const QpVector& gradient) { HighsInt minidx = chooseconstrainttodrop(redcosts.getReducedCosts()); diff --git a/src/qpsolver/devexharrispricing.hpp b/src/qpsolver/devexharrispricing.hpp index e86b6e4291..cbcd6dc476 100644 --- a/src/qpsolver/devexharrispricing.hpp +++ b/src/qpsolver/devexharrispricing.hpp @@ -54,7 +54,7 @@ class DevexHarrisPricing : public Pricing { public: DevexHarrisPricing(Runtime& rt, Basis& bas, ReducedCosts& rc) : runtime(rt), - basis(bas), + basis(bas), redcosts(rc), weights(std::vector(rt.instance.num_var, 1.0)) {}; diff --git a/src/qpsolver/devexpricing.hpp b/src/qpsolver/devexpricing.hpp index ba34ac04e5..a88f3ecea5 100644 --- a/src/qpsolver/devexpricing.hpp +++ b/src/qpsolver/devexpricing.hpp @@ -58,7 +58,7 @@ class DevexPricing : public Pricing { : runtime(rt), basis(bas), redcosts(rc), - weights(std::vector(rt.instance.num_var, 1.0)){}; + weights(std::vector(rt.instance.num_var, 1.0)) {}; // B lambda = g // lambda = inv(B)g diff --git a/src/qpsolver/factor.hpp b/src/qpsolver/factor.hpp index 9a48288aed..1510357498 100644 --- a/src/qpsolver/factor.hpp +++ b/src/qpsolver/factor.hpp @@ -36,7 +36,7 @@ class CholeskyFactor { HighsInt min_k_max = min(new_k_max, current_k_max); for (HighsInt i = 0; i < min_k_max; i++) { for (HighsInt j = 0; j < min_k_max; j++) { - assert(i * (new_k_max) + j < l_size); + assert(i * (new_k_max) + j < l_size); L[i * (new_k_max) + j] = L_old[i * current_k_max + j]; } } @@ -52,7 +52,6 @@ class CholeskyFactor { L.resize(current_k_max * current_k_max); } - void recompute() { std::vector> orig; HighsInt dim_ns = basis.getinactive().size(); @@ -97,7 +96,8 @@ class CholeskyFactor { uptodate = true; } - QpSolverStatus expand(const QpVector& yp, QpVector& gyp, QpVector& l, QpVector& m) { + QpSolverStatus expand(const QpVector& yp, QpVector& gyp, QpVector& l, + QpVector& m) { if (!uptodate) { return QpSolverStatus::OK; } @@ -124,14 +124,15 @@ class CholeskyFactor { // b*b -a*a = mu // k(b-a) = 1 // b + a = k*mu - // Commented out unreachable code + // Commented out unreachable code // const double tolerance = 0.001; // // double beta = max(tolerance, sqrt(m.norm2() / L[0] + fabs(mu))); // double k = 1 / (beta + sqrt(beta * beta - mu)); // double alpha = k * mu - beta; // - // printf("k = %d, alpha = %lf, beta = %lf, k = %lf\n", (int)current_k, alpha, + // printf("k = %d, alpha = %lf, beta = %lf, k = %lf\n", + // (int)current_k, alpha, // beta, k); // // a.clear(); diff --git a/src/qpsolver/feasibility_bounded.hpp b/src/qpsolver/feasibility_bounded.hpp index 4a92c813e1..5b8c0de534 100644 --- a/src/qpsolver/feasibility_bounded.hpp +++ b/src/qpsolver/feasibility_bounded.hpp @@ -5,13 +5,13 @@ #include "qpsolver/a_asm.hpp" #include "qpsolver/crashsolution.hpp" -static void computeStartingPointBounded(Instance& instance, - Settings& settings, - Statistics& stats, - QpModelStatus& modelstatus, - QpHotstartInformation& result, - HighsTimer& timer) { - // compute initial feasible point for problems with bounds only (no general linear constraints) +static void computeStartingPointBounded(Instance& instance, Settings& settings, + Statistics& stats, + QpModelStatus& modelstatus, + QpHotstartInformation& result, + HighsTimer& timer) { + // compute initial feasible point for problems with bounds only (no general + // linear constraints) // compute Qx + c = 0 --> x = Q^-1c std::vector L; @@ -19,16 +19,18 @@ static void computeStartingPointBounded(Instance& instance, // compute cholesky factorization of Q for (size_t col = 0; col < (size_t)instance.num_var; col++) { - for (size_t idx = instance.Q.mat.start[col]; idx < (size_t)instance.Q.mat.start[col+1]; idx++) { + for (size_t idx = instance.Q.mat.start[col]; + idx < (size_t)instance.Q.mat.start[col + 1]; idx++) { double sum = 0; - size_t row = instance.Q.mat.index[idx]; + size_t row = instance.Q.mat.index[idx]; if (row == col) { for (size_t k = 0; k < row; k++) sum += L[k * instance.num_var + row] * L[k * instance.num_var + row]; L[row * instance.num_var + row] = sqrt(instance.Q.mat.value[idx] - sum); } else { for (size_t k = 0; k < row; k++) - sum += (L[k * instance.num_var + col] * L[k * instance.num_var + row]); + sum += + (L[k * instance.num_var + col] * L[k * instance.num_var + row]); L[row * instance.num_var + col] = (instance.Q.mat.value[idx] - sum) / L[row * instance.num_var + row]; } @@ -37,7 +39,7 @@ static void computeStartingPointBounded(Instance& instance, // solve for c QpVector res = -instance.c; - for (HighsInt r = 0; r initialinactive; std::vector atlower; - for (int i=0; i 0.5/settings.hessianregularizationfactor - && instance.var_up[i] == std::numeric_limits::infinity() - && instance.c.value[i] < 0.0) { - modelstatus = QpModelStatus::kUnbounded; - return; - } else if (res.value[i] < 0.5/settings.hessianregularizationfactor - && instance.var_lo[i] == std::numeric_limits::infinity() - && instance.c.value[i] > 0.0) { - modelstatus = QpModelStatus::kUnbounded; - return; + for (int i = 0; i < instance.num_var; i++) { + if (res.value[i] > 0.5 / settings.hessianregularizationfactor && + instance.var_up[i] == std::numeric_limits::infinity() && + instance.c.value[i] < 0.0) { + modelstatus = QpModelStatus::kUnbounded; + return; + } else if (res.value[i] < 0.5 / settings.hessianregularizationfactor && + instance.var_lo[i] == std::numeric_limits::infinity() && + instance.c.value[i] > 0.0) { + modelstatus = QpModelStatus::kUnbounded; + return; } else if (res.value[i] <= instance.var_lo[i]) { - res.value[i] = instance.var_lo[i]; - initialactive.push_back(i + instance.num_con); - atlower.push_back(BasisStatus::kActiveAtLower); + res.value[i] = instance.var_lo[i]; + initialactive.push_back(i + instance.num_con); + atlower.push_back(BasisStatus::kActiveAtLower); } else if (res.value[i] >= instance.var_up[i]) { - res.value[i] = instance.var_up[i]; - initialactive.push_back(i + instance.num_con); - atlower.push_back(BasisStatus::kActiveAtUpper); + res.value[i] = instance.var_up[i]; + initialactive.push_back(i + instance.num_con); + atlower.push_back(BasisStatus::kActiveAtUpper); } else { - initialinactive.push_back(i + instance.num_con); + initialinactive.push_back(i + instance.num_con); } if (fabs(res.value[i]) > 1e-4) { x0.value[i] = res.value[i]; diff --git a/src/qpsolver/feasibility_highs.hpp b/src/qpsolver/feasibility_highs.hpp index 8594be00df..5c95a2d751 100644 --- a/src/qpsolver/feasibility_highs.hpp +++ b/src/qpsolver/feasibility_highs.hpp @@ -5,15 +5,11 @@ #include "qpsolver/a_asm.hpp" #include "qpsolver/crashsolution.hpp" -static void computeStartingPointHighs(Instance& instance, - Settings& settings, - Statistics& stats, - QpModelStatus& modelstatus, - QpHotstartInformation& result, - HighsModelStatus& highs_model_status, - HighsBasis& highs_basis, - HighsSolution& highs_solution, - HighsTimer& timer) { +static void computeStartingPointHighs( + Instance& instance, Settings& settings, Statistics& stats, + QpModelStatus& modelstatus, QpHotstartInformation& result, + HighsModelStatus& highs_model_status, HighsBasis& highs_basis, + HighsSolution& highs_solution, HighsTimer& timer) { bool have_starting_point = false; const bool debug_report = false; if (highs_solution.value_valid) { @@ -28,24 +24,26 @@ static void computeStartingPointHighs(Instance& instance, double sum_con_infeasibilities = 0; double max_con_residual = 0; double sum_con_residuals = 0; - - assessQpPrimalFeasibility(instance, primal_feasibility_tolerance, - highs_solution.col_value, highs_solution.row_value, - num_var_infeasibilities, max_var_infeasibility, sum_var_infeasibilities, - num_con_infeasibilities, max_con_infeasibility, sum_con_infeasibilities, - max_con_residual, sum_con_residuals); - if (debug_report) printf("computeStartingPointHighs highs_solution has (num / max / sum) " - "var (%d / %g / %g) and " - "con (%d / %g / %g) infeasibilities " - "with (max = %g; sum = %g) residuals\n", - int(num_var_infeasibilities), max_var_infeasibility, sum_var_infeasibilities, - int(num_con_infeasibilities), max_con_infeasibility, sum_con_infeasibilities, - max_con_residual, sum_con_residuals); - have_starting_point = - num_var_infeasibilities == 0 && - num_con_infeasibilities == 0 && - highs_basis.valid; + assessQpPrimalFeasibility( + instance, primal_feasibility_tolerance, highs_solution.col_value, + highs_solution.row_value, num_var_infeasibilities, + max_var_infeasibility, sum_var_infeasibilities, num_con_infeasibilities, + max_con_infeasibility, sum_con_infeasibilities, max_con_residual, + sum_con_residuals); + + if (debug_report) + printf( + "computeStartingPointHighs highs_solution has (num / max / sum) " + "var (%d / %g / %g) and " + "con (%d / %g / %g) infeasibilities " + "with (max = %g; sum = %g) residuals\n", + int(num_var_infeasibilities), max_var_infeasibility, + sum_var_infeasibilities, int(num_con_infeasibilities), + max_con_infeasibility, sum_con_infeasibilities, max_con_residual, + sum_con_residuals); + have_starting_point = num_var_infeasibilities == 0 && + num_con_infeasibilities == 0 && highs_basis.valid; } // compute initial feasible point HighsBasis use_basis; @@ -61,10 +59,11 @@ static void computeStartingPointHighs(Instance& instance, highs.setOptionValue("presolve", "on"); - const double use_time_limit = settings.time_limit < kHighsInf ? - settings.time_limit - timer.readRunHighsClock() : - kHighsInf; - + const double use_time_limit = + settings.time_limit < kHighsInf + ? settings.time_limit - timer.readRunHighsClock() + : kHighsInf; + highs.setOptionValue("time_limit", use_time_limit); HighsLp lp; @@ -84,11 +83,11 @@ static void computeStartingPointHighs(Instance& instance, // create artificial bounds for free variables: false by default assert(!settings.phase1boundfreevars); if (settings.phase1boundfreevars) { - for (HighsInt i=0; i initial_active; std::vector initial_inactive; std::vector initial_status; - const HighsInt num_highs_basis_status = HighsInt(HighsBasisStatus::kNonbasic)+1; + const HighsInt num_highs_basis_status = + HighsInt(HighsBasisStatus::kNonbasic) + 1; std::vector debug_row_status_count; debug_row_status_count.assign(num_highs_basis_status, 0); for (HighsInt i = 0; i < HighsInt(use_basis.row_status.size()); i++) { @@ -181,7 +182,7 @@ static void computeStartingPointHighs(Instance& instance, // Shouldn't happen, since free rows are basic in a logical // basis and remain basic, or are removed by presolve and // restored as basic in postsolve - assert(111==222); + assert(111 == 222); // That said, a free row that is nonbasic in the Highs basis // must be counted as inactive in the QP basis for accounting // purposes @@ -215,7 +216,7 @@ static void computeStartingPointHighs(Instance& instance, initial_active.push_back(instance.num_con + i); initial_status.push_back(BasisStatus::kActiveAtLower); } - + } else if (status == HighsBasisStatus::kUpper) { if (isfreevar(instance, i)) { initial_inactive.push_back(instance.num_con + i); @@ -223,7 +224,7 @@ static void computeStartingPointHighs(Instance& instance, initial_active.push_back(instance.num_con + i); initial_status.push_back(BasisStatus::kActiveAtUpper); } - + } else if (status == HighsBasisStatus::kZero) { initial_inactive.push_back(instance.num_con + i); } else if (status != HighsBasisStatus::kBasic) { @@ -236,29 +237,29 @@ static void computeStartingPointHighs(Instance& instance, if (debug_report) { printf("QP solver initial basis: (Lo / Bs / Up / Ze / Nb) for cols ("); - for (HighsInt k = 0; k < num_highs_basis_status; k++) - printf("%s%d", k==0 ? "" : " / ", int(debug_col_status_count[k])); + for (HighsInt k = 0; k < num_highs_basis_status; k++) + printf("%s%d", k == 0 ? "" : " / ", int(debug_col_status_count[k])); printf(") and rows ("); - for (HighsInt k = 0; k < num_highs_basis_status; k++) - printf("%s%d", k==0 ? "" : " / ", int(debug_row_status_count[k])); + for (HighsInt k = 0; k < num_highs_basis_status; k++) + printf("%s%d", k == 0 ? "" : " / ", int(debug_row_status_count[k])); printf(")\n"); } assert((HighsInt)(initial_active.size() + initial_inactive.size()) == - instance.num_var); + instance.num_var); if (!have_starting_point) { // When starting from a feasible basis, there will generally be // inactive variables in the basis that aren't free for (HighsInt ia : initial_inactive) { if (ia < instance.num_con) { - // printf("free row %d\n", (int)ia); - assert(instance.con_lo[ia] == -kHighsInf); - assert(instance.con_up[ia] == kHighsInf); + // printf("free row %d\n", (int)ia); + assert(instance.con_lo[ia] == -kHighsInf); + assert(instance.con_up[ia] == kHighsInf); } else { - // printf("free col %d\n", (int)ia); - assert(instance.var_lo[ia - instance.num_con] == -kHighsInf); - assert(instance.var_up[ia - instance.num_con] == kHighsInf); + // printf("free col %d\n", (int)ia); + assert(instance.var_lo[ia - instance.num_con] == -kHighsInf); + assert(instance.var_up[ia - instance.num_con] == kHighsInf); } } } diff --git a/src/qpsolver/gradient.hpp b/src/qpsolver/gradient.hpp index c0ec9d0807..f62366a2f4 100644 --- a/src/qpsolver/gradient.hpp +++ b/src/qpsolver/gradient.hpp @@ -1,8 +1,8 @@ #ifndef __SRC_LIB_GRADIENT_HPP__ #define __SRC_LIB_GRADIENT_HPP__ -#include "runtime.hpp" #include "qpvector.hpp" +#include "runtime.hpp" class Gradient { Runtime& runtime; diff --git a/src/qpsolver/instance.hpp b/src/qpsolver/instance.hpp index 060667307a..0ff9c32376 100644 --- a/src/qpsolver/instance.hpp +++ b/src/qpsolver/instance.hpp @@ -12,7 +12,7 @@ struct SumNum { }; struct Instance { - HighsInt sense = 1; // Minimization + HighsInt sense = 1; // Minimization HighsInt num_var = 0; HighsInt num_con = 0; double offset = 0; diff --git a/src/qpsolver/pricing.hpp b/src/qpsolver/pricing.hpp index 884fab2e56..3a1857b946 100644 --- a/src/qpsolver/pricing.hpp +++ b/src/qpsolver/pricing.hpp @@ -6,8 +6,8 @@ class Pricing { public: virtual HighsInt price(const QpVector& x, const QpVector& gradient) = 0; - virtual void update_weights(const QpVector& aq, const QpVector& ep, HighsInt p, - HighsInt q) = 0; + virtual void update_weights(const QpVector& aq, const QpVector& ep, + HighsInt p, HighsInt q) = 0; virtual void recompute() = 0; virtual ~Pricing() {} }; diff --git a/src/qpsolver/qpconst.hpp b/src/qpsolver/qpconst.hpp index 633ef2de37..3d1ad9304a 100644 --- a/src/qpsolver/qpconst.hpp +++ b/src/qpsolver/qpconst.hpp @@ -5,7 +5,7 @@ enum class QpSolverStatus { OK, NOTPOSITIVDEFINITE, DEGENERATE }; enum class QpModelStatus { - kNotset, // 0 + kNotset, // 0 kUndetermined, kOptimal, kUnbounded, @@ -23,5 +23,4 @@ enum class BasisStatus { kInactiveInBasis }; - #endif diff --git a/src/qpsolver/quass.cpp b/src/qpsolver/quass.cpp index 2b0f825e69..2f7cd40648 100644 --- a/src/qpsolver/quass.cpp +++ b/src/qpsolver/quass.cpp @@ -1,10 +1,12 @@ #include "qpsolver/quass.hpp" +#include + #include #include -#include #include "Highs.h" +#include "lp_data/HighsAnalysis.h" #include "qpsolver/basis.hpp" #include "qpsolver/crashsolution.hpp" #include "qpsolver/dantzigpricing.hpp" @@ -13,18 +15,18 @@ #include "qpsolver/factor.hpp" #include "qpsolver/gradient.hpp" #include "qpsolver/instance.hpp" -#include "lp_data/HighsAnalysis.h" +#include "qpsolver/perturbation.hpp" #include "qpsolver/ratiotest.hpp" #include "qpsolver/reducedcosts.hpp" #include "qpsolver/reducedgradient.hpp" +#include "qpsolver/scaling.hpp" #include "qpsolver/snippets.hpp" #include "qpsolver/steepestedgepricing.hpp" -#include "qpsolver/scaling.hpp" -#include "qpsolver/perturbation.hpp" Quass::Quass(Runtime& rt) : runtime(rt) {} -static void loginformation(Runtime& rt, Basis& basis, CholeskyFactor& factor, HighsTimer& timer) { +static void loginformation(Runtime& rt, Basis& basis, CholeskyFactor& factor, + HighsTimer& timer) { rt.statistics.iteration.push_back(rt.statistics.num_iterations); rt.statistics.nullspacedimension.push_back(rt.instance.num_var - basis.getnumactive()); @@ -38,7 +40,8 @@ static void loginformation(Runtime& rt, Basis& basis, CholeskyFactor& factor, Hi rt.statistics.density_nullspace.push_back(0.0); } -static void tidyup(QpVector& p, QpVector& rowmove, Basis& basis, Runtime& runtime) { +static void tidyup(QpVector& p, QpVector& rowmove, Basis& basis, + Runtime& runtime) { for (unsigned acon : basis.getactive()) { if ((HighsInt)acon >= runtime.instance.num_con) { p.value[acon - runtime.instance.num_con] = 0.0; @@ -49,7 +52,7 @@ static void tidyup(QpVector& p, QpVector& rowmove, Basis& basis, Runtime& runtim } static void computerowmove(Runtime& runtime, Basis& basis, QpVector& p, - QpVector& rowmove) { + QpVector& rowmove) { runtime.instance.A.mat_vec(p, rowmove); return; // rowmove.reset(); @@ -77,9 +80,10 @@ static void computerowmove(Runtime& runtime, Basis& basis, QpVector& p, // VECTOR static QpVector& computesearchdirection_minor(Runtime& rt, Basis& bas, - CholeskyFactor& cf, - ReducedGradient& redgrad, QpVector& p) { - QpVector g2 = -redgrad.get(); // TODO PERF: buffer QpVector + CholeskyFactor& cf, + ReducedGradient& redgrad, + QpVector& p) { + QpVector g2 = -redgrad.get(); // TODO PERF: buffer QpVector g2.sanitize(); cf.solve(g2); @@ -89,11 +93,10 @@ static QpVector& computesearchdirection_minor(Runtime& rt, Basis& bas, } // VECTOR -static QpVector& computesearchdirection_major(Runtime& runtime, Basis& basis, - CholeskyFactor& factor, const QpVector& yp, - Gradient& gradient, QpVector& gyp, QpVector& l, - QpVector& m, QpVector& p) { - QpVector yyp = yp; // TODO PERF: buffer QpVector +static QpVector& computesearchdirection_major( + Runtime& runtime, Basis& basis, CholeskyFactor& factor, const QpVector& yp, + Gradient& gradient, QpVector& gyp, QpVector& l, QpVector& m, QpVector& p) { + QpVector yyp = yp; // TODO PERF: buffer QpVector // if (gradient.getGradient().dot(yp) > 0.0) { // yyp.scale(-1.0); // } @@ -102,7 +105,7 @@ static QpVector& computesearchdirection_major(Runtime& runtime, Basis& basis, basis.Ztprod(gyp, m); l = m; factor.solveL(l); - QpVector v = l; // TODO PERF: buffer QpVector + QpVector v = l; // TODO PERF: buffer QpVector factor.solveLT(v); basis.Zprod(v, p); if (gradient.getGradient().dot(yyp) < 0.0) { @@ -118,7 +121,8 @@ static QpVector& computesearchdirection_major(Runtime& runtime, Basis& basis, } static double computemaxsteplength(Runtime& runtime, const QpVector& p, - Gradient& gradient, QpVector& buffer_Qp, bool& zcd) { + Gradient& gradient, QpVector& buffer_Qp, + bool& zcd) { double denominator = p * runtime.instance.Q.mat_vec(p, buffer_Qp); if (fabs(denominator) > runtime.settings.pQp_zero_threshold) { double numerator = -(p * gradient.getGradient()); @@ -133,9 +137,9 @@ static double computemaxsteplength(Runtime& runtime, const QpVector& p, } } -static QpSolverStatus reduce(Runtime& rt, Basis& basis, const HighsInt newactivecon, - QpVector& buffer_d, HighsInt& maxabsd, - HighsInt& constrainttodrop) { +static QpSolverStatus reduce(Runtime& rt, Basis& basis, + const HighsInt newactivecon, QpVector& buffer_d, + HighsInt& maxabsd, HighsInt& constrainttodrop) { HighsInt idx = indexof(basis.getinactive(), newactivecon); if (idx != -1) { maxabsd = idx; @@ -146,7 +150,8 @@ static QpSolverStatus reduce(Runtime& rt, Basis& basis, const HighsInt newactive } // TODO: this operation is inefficient. - QpVector aq = rt.instance.A.t().extractcol(newactivecon); // TODO PERF: buffer QpVector + QpVector aq = + rt.instance.A.t().extractcol(newactivecon); // TODO PERF: buffer QpVector basis.Ztprod(aq, buffer_d, true, newactivecon); maxabsd = 0; @@ -169,7 +174,7 @@ static QpSolverStatus reduce(Runtime& rt, Basis& basis, const HighsInt newactive } static std::unique_ptr getPricing(Runtime& runtime, Basis& basis, - ReducedCosts& redcosts) { + ReducedCosts& redcosts) { switch (runtime.settings.pricing) { case PricingStrategy::SteepestEdge: return std::unique_ptr( @@ -186,7 +191,7 @@ static std::unique_ptr getPricing(Runtime& runtime, Basis& basis, static void regularize(Runtime& rt) { if (!rt.settings.hessianregularization) { - return; + return; } // add small diagonal to hessian for (HighsInt i = 0; i < rt.instance.num_var; i++) { @@ -271,18 +276,20 @@ static bool check_reinvert_due(Basis& basis) { return basis.getreinversionhint(); } -static void reinvert(Basis& basis, CholeskyFactor& factor, Gradient& grad, ReducedCosts& rc, ReducedGradient& rg, std::unique_ptr& pricing) { +static void reinvert(Basis& basis, CholeskyFactor& factor, Gradient& grad, + ReducedCosts& rc, ReducedGradient& rg, + std::unique_ptr& pricing) { basis.rebuild(); factor.recompute(); grad.recompute(); rc.recompute(); rg.recompute(); - //pricing->recompute(); + // pricing->recompute(); } -void Quass::solve(const QpVector& x0, const QpVector& ra, Basis& b0, HighsTimer& timer) { - - //feenableexcept(FE_ALL_EXCEPT & ~FE_INEXACT & ~FE_UNDERFLOW); +void Quass::solve(const QpVector& x0, const QpVector& ra, Basis& b0, + HighsTimer& timer) { + // feenableexcept(FE_ALL_EXCEPT & ~FE_INEXACT & ~FE_UNDERFLOW); runtime.statistics.time_start = std::chrono::high_resolution_clock::now(); Basis& basis = b0; @@ -316,18 +323,17 @@ void Quass::solve(const QpVector& x0, const QpVector& ra, Basis& b0, HighsTimer& runtime.relaxed_for_ratiotest = ratiotest_relax_instance(runtime); - - HighsInt last_logging_iteration = runtime.statistics.num_iterations-1; + HighsInt last_logging_iteration = runtime.statistics.num_iterations - 1; double last_logging_time = 0; double logging_time_interval = 10; - + const HighsInt current_num_active = basis.getnumactive(); bool atfsep = current_num_active == runtime.instance.num_var; while (true) { // check iteration limit if (runtime.statistics.num_iterations >= runtime.settings.iteration_limit) { runtime.status = QpModelStatus::kIterationLimit; - break; + break; } // check time limit @@ -336,9 +342,9 @@ void Quass::solve(const QpVector& x0, const QpVector& ra, Basis& b0, HighsTimer& break; } - if (basis.getnuminactive() > runtime.settings.nullspace_limit) { - runtime.settings.nullspace_limit_log.fire(runtime.settings.nullspace_limit); + runtime.settings.nullspace_limit_log.fire( + runtime.settings.nullspace_limit); runtime.status = QpModelStatus::kLargeNullspace; return; } @@ -346,22 +352,22 @@ void Quass::solve(const QpVector& x0, const QpVector& ra, Basis& b0, HighsTimer& // LOGGING double run_time = timer.readRunHighsClock(); if ((runtime.statistics.num_iterations % - runtime.settings.reportingfequency == - 0 || - run_time-last_logging_time > logging_time_interval) && - runtime.statistics.num_iterations > last_logging_iteration) { + runtime.settings.reportingfequency == + 0 || + run_time - last_logging_time > logging_time_interval) && + runtime.statistics.num_iterations > last_logging_iteration) { bool log_report = true; - if (runtime.statistics.num_iterations > 10*runtime.settings.reportingfequency) { - runtime.settings.reportingfequency *= 10; - log_report = false; + if (runtime.statistics.num_iterations > + 10 * runtime.settings.reportingfequency) { + runtime.settings.reportingfequency *= 10; + log_report = false; } - if (run_time > 10*logging_time_interval) - logging_time_interval *= 2.0; + if (run_time > 10 * logging_time_interval) logging_time_interval *= 2.0; if (log_report) { - last_logging_time = run_time; - last_logging_iteration = runtime.statistics.num_iterations; - loginformation(runtime, basis, factor, timer); - runtime.settings.iteration_log.fire(runtime.statistics); + last_logging_time = run_time; + last_logging_iteration = runtime.statistics.num_iterations; + loginformation(runtime, basis, factor, timer); + runtime.settings.iteration_log.fire(runtime.statistics); } } @@ -419,7 +425,9 @@ void Quass::solve(const QpVector& x0, const QpVector& ra, Basis& b0, HighsTimer& runtime.instance.Q.mat_vec(p, buffer_Qp); } if (p.norm2() < runtime.settings.pnorm_zero_threshold || - maxsteplength == 0.0 || (false && fabs(gradient.getGradient().dot(p)) < runtime.settings.improvement_zero_threshold)) { + maxsteplength == 0.0 || + (false && fabs(gradient.getGradient().dot(p)) < + runtime.settings.improvement_zero_threshold)) { atfsep = true; } else { // Now perform a real iteration @@ -456,8 +464,7 @@ void Quass::solve(const QpVector& x0, const QpVector& ra, Basis& b0, HighsTimer& atfsep = false; } } else { - if (stepres.alpha == - std::numeric_limits::infinity()) { + if (stepres.alpha == std::numeric_limits::infinity()) { // unbounded runtime.status = QpModelStatus::kUnbounded; return; @@ -495,18 +502,20 @@ void Quass::solve(const QpVector& x0, const QpVector& ra, Basis& b0, HighsTimer& runtime.dualcon.resparsify(); runtime.dualvar.resparsify(); - //QpVector actual_dual_var(runtime.instance.num_var); - //QpVector actual_dual_con(runtime.instance.num_con); - //compute_actual_duals(runtime, basis, redcosts.getReducedCosts(), actual_dual_con, actual_dual_var); - //printf("max primal violation = %.20lf\n", compute_primal_violation(runtime)); - //printf("max dual violation = %.20lf\n", compute_dual_violation(runtime.instance, runtime.primal, actual_dual_con, actual_dual_var)); + // QpVector actual_dual_var(runtime.instance.num_var); + // QpVector actual_dual_con(runtime.instance.num_con); + // compute_actual_duals(runtime, basis, redcosts.getReducedCosts(), + // actual_dual_con, actual_dual_var); printf("max primal violation = + // %.20lf\n", compute_primal_violation(runtime)); printf("max dual violation + // = %.20lf\n", compute_dual_violation(runtime.instance, runtime.primal, + // actual_dual_con, actual_dual_var)); // extract basis status - for (HighsInt i=0; i::infinity()) { - bound -= runtime.settings.ratiotest_d; - } - } - - for (double& bound : relaxed_instance.con_up) { - if (bound != std::numeric_limits::infinity()) { - bound += runtime.settings.ratiotest_d; - } - } - - for (double& bound : relaxed_instance.var_lo) { - if (bound != -std::numeric_limits::infinity()) { - bound -= runtime.settings.ratiotest_d; - } - } - - for (double& bound : relaxed_instance.var_up) { - if (bound != std::numeric_limits::infinity()) { - bound += runtime.settings.ratiotest_d; - } - } + for (double& bound : relaxed_instance.con_lo) { + if (bound != -std::numeric_limits::infinity()) { + bound -= runtime.settings.ratiotest_d; + } + } + + for (double& bound : relaxed_instance.con_up) { + if (bound != std::numeric_limits::infinity()) { + bound += runtime.settings.ratiotest_d; + } + } + + for (double& bound : relaxed_instance.var_lo) { + if (bound != -std::numeric_limits::infinity()) { + bound -= runtime.settings.ratiotest_d; + } + } + + for (double& bound : relaxed_instance.var_up) { + if (bound != std::numeric_limits::infinity()) { + bound += runtime.settings.ratiotest_d; + } + } return relaxed_instance; } @@ -131,7 +132,7 @@ RatiotestResult ratiotest(Runtime& runtime, const QpVector& p, alphastart); case RatiotestStrategy::TwoPass: default: // to fix -Wreturn-type warning - return ratiotest_twopass(runtime, p, rowmove, runtime.relaxed_for_ratiotest, - alphastart); + return ratiotest_twopass(runtime, p, rowmove, + runtime.relaxed_for_ratiotest, alphastart); } } diff --git a/src/qpsolver/reducedcosts.hpp b/src/qpsolver/reducedcosts.hpp index e31e9f7ad2..2b00adcd14 100644 --- a/src/qpsolver/reducedcosts.hpp +++ b/src/qpsolver/reducedcosts.hpp @@ -3,8 +3,8 @@ #include "qpsolver/basis.hpp" #include "qpsolver/gradient.hpp" -#include "qpsolver/runtime.hpp" #include "qpsolver/qpvector.hpp" +#include "qpsolver/runtime.hpp" class ReducedCosts { Basis& basis; diff --git a/src/qpsolver/reducedgradient.hpp b/src/qpsolver/reducedgradient.hpp index b124e29c4b..97f90b4354 100644 --- a/src/qpsolver/reducedgradient.hpp +++ b/src/qpsolver/reducedgradient.hpp @@ -2,8 +2,8 @@ #define __SRC_LIB_REDUCEDGRADIENT_HPP__ #include "qpsolver/basis.hpp" -#include "qpsolver/runtime.hpp" #include "qpsolver/qpvector.hpp" +#include "qpsolver/runtime.hpp" class ReducedGradient { QpVector rg; diff --git a/src/qpsolver/runtime.hpp b/src/qpsolver/runtime.hpp index 5c9c199a9a..2d045bf0e6 100644 --- a/src/qpsolver/runtime.hpp +++ b/src/qpsolver/runtime.hpp @@ -1,11 +1,11 @@ #ifndef __SRC_LIB_RUNTIME_HPP__ #define __SRC_LIB_RUNTIME_HPP__ -#include "util/HighsTimer.h" #include "instance.hpp" +#include "qpsolver/qpconst.hpp" #include "settings.hpp" #include "statistics.hpp" -#include "qpsolver/qpconst.hpp" +#include "util/HighsTimer.h" struct Runtime { Instance instance; diff --git a/src/qpsolver/scaling.cpp b/src/qpsolver/scaling.cpp index fdb1c8ea4e..e0c4463367 100644 --- a/src/qpsolver/scaling.cpp +++ b/src/qpsolver/scaling.cpp @@ -3,15 +3,13 @@ #include #include -static -double largestpoweroftwo(double value) { +static double largestpoweroftwo(double value) { double l = log2(value); HighsInt il = (HighsInt)l; return powf(1.0, il); } -static -void scale_rows(Runtime& rt) { +static void scale_rows(Runtime& rt) { if (!rt.settings.rowscaling) { return; } @@ -52,8 +50,7 @@ void scale_rows(Runtime& rt) { } } -static -void scale_cols(Runtime& rt) { +static void scale_cols(Runtime& rt) { if (!rt.settings.varscaling) { return; } diff --git a/src/qpsolver/settings.hpp b/src/qpsolver/settings.hpp index 37d9d12e24..18441e5164 100644 --- a/src/qpsolver/settings.hpp +++ b/src/qpsolver/settings.hpp @@ -18,14 +18,26 @@ struct Settings { PricingStrategy pricing = PricingStrategy::Devex; - double pnorm_zero_threshold = 1e-11; // if ||p|| < this threshold, p is determined to not be an improving search direction - double improvement_zero_threshold = 1e-4; // if p^t gradient < this threshold, p is determined to not be an improving search direction - double d_zero_threshold = 1e-12; // minimal value for pivot, will declare degeneracy if no larger pivot is found - double lambda_zero_threshold = 1e-9; // used for pricing / optimality checking - double pQp_zero_threshold = 1e-7; // if p'Qp < this, p is determined to not have curvature, a simplex-like iteration is performed. + double pnorm_zero_threshold = + 1e-11; // if ||p|| < this threshold, p is determined to not be an + // improving search direction + double improvement_zero_threshold = + 1e-4; // if p^t gradient < this threshold, p is determined to not be an + // improving search direction + double d_zero_threshold = 1e-12; // minimal value for pivot, will declare + // degeneracy if no larger pivot is found + double lambda_zero_threshold = + 1e-9; // used for pricing / optimality checking + double pQp_zero_threshold = + 1e-7; // if p'Qp < this, p is determined to not have curvature, a + // simplex-like iteration is performed. - bool hessianregularization = false; // if true, a small multiple of the identity matrix will be added to the Hessian - double hessianregularizationfactor = 1e-7; // multiple of identity matrix added to hessian in case of regularization + bool hessianregularization = + false; // if true, a small multiple of the identity matrix will be added + // to the Hessian + double hessianregularizationfactor = + 1e-7; // multiple of identity matrix added to hessian in case of + // regularization Phase1Strategy phase1strategy = Phase1Strategy::HIGHS; bool phase1movefreevarsbasic = false; @@ -37,7 +49,7 @@ struct Settings { Eventhandler nullspace_limit_log; HighsInt nullspace_limit = 4000; - + HighsInt reinvertfrequency = 1000; HighsInt gradientrecomputefrequency = 100; HighsInt reducedgradientrecomputefrequency = diff --git a/src/qpsolver/steepestedgepricing.hpp b/src/qpsolver/steepestedgepricing.hpp index 0a216394ee..95beecb363 100644 --- a/src/qpsolver/steepestedgepricing.hpp +++ b/src/qpsolver/steepestedgepricing.hpp @@ -57,22 +57,21 @@ class SteepestEdgePricing : public Pricing { basis(bas), redcosts(rc), weights(std::vector(rt.instance.num_var, 1.0)) { - compute_exact_weights(); - }; + compute_exact_weights(); + }; HighsInt price(const QpVector& x, const QpVector& gradient) { HighsInt minidx = chooseconstrainttodrop(redcosts.getReducedCosts()); return minidx; } - double compute_exact_weight(HighsInt i) { QpVector y_i = basis.btran(QpVector::unit(runtime.instance.num_var, i)); return y_i.dot(y_i); } void compute_exact_weights() { - for (int i=0; i 1e-2) { - //printf("weights[%d] = %lf, should be %lf\n", i, weight_is, weight_comp); + if (fabs(weight_comp - weight_is) > 1e-2) { + // printf("weights[%d] = %lf, should be %lf\n", i, weight_is, + // weight_comp); return false; } return true; } bool check_weights() { - std::vector correct_weights; std::vector incorrect_weights; - for (int i=0; i= 1e-2) { - // printf("old weight[p] discrepancy: updated = %lf, computed=%lf\n", old_weight_p_updated, old_weight_p_computed); - //} + // if (fabs(old_weight_p_computed - old_weight_p_updated) >= 1e-2) { + // printf("old weight[p] discrepancy: updated = %lf, computed=%lf\n", + // old_weight_p_updated, old_weight_p_computed); + // } double weight_p = old_weight_p_computed; - double t_p = aq.value[rowindex_p]; for (HighsInt i = 0; i < runtime.instance.num_var; i++) { if (i != rowindex_p) { double t_i = aq.value[i]; - weights[i] = weights[i] - - 2 * (t_i / t_p) * delta.value[i] - + ((t_i * t_i) / (t_p * t_p)) * weight_p; - //printf("weights[%d] = %lf\n", i, weights[i]); + weights[i] = weights[i] - 2 * (t_i / t_p) * delta.value[i] + + ((t_i * t_i) / (t_p * t_p)) * weight_p; + // printf("weights[%d] = %lf\n", i, weights[i]); } } - //QpVector new_ep = basis.btran(QpVector::unit(runtime.instance.num_var, rowindex_p)); - //double computed_weight = new_ep.dot(new_ep); + // QpVector new_ep = basis.btran(QpVector::unit(runtime.instance.num_var, + // rowindex_p)); double computed_weight = new_ep.dot(new_ep); double new_weight_p_updated = weight_p / (t_p * t_p); - - //if (fabs(updated_weight - computed_weight) > 1e-4) { - // printf("updated weight %lf vs computed weight %lf. aq[p] = %lf\n", updated_weight, computed_weight, t_p); - // printf("old weight = %lf, aq[p] = %lf, ^2 = %lf, new weight = %lf\n", weight_p, t_p, t_p*t_p, updated_weight); - //} + + // if (fabs(updated_weight - computed_weight) > 1e-4) { + // printf("updated weight %lf vs computed weight %lf. aq[p] = %lf\n", + // updated_weight, computed_weight, t_p); printf("old weight = %lf, aq[p] + // = %lf, ^2 = %lf, new weight = %lf\n", weight_p, t_p, t_p*t_p, + // updated_weight); + // } weights[rowindex_p] = new_weight_p_updated; } };