diff --git a/uno/solvers/BQPD/BQPDSolver.cpp b/uno/solvers/BQPD/BQPDSolver.cpp index e9b6104e..6192f184 100644 --- a/uno/solvers/BQPD/BQPDSolver.cpp +++ b/uno/solvers/BQPD/BQPDSolver.cpp @@ -80,7 +80,7 @@ namespace uno { if (warmstart_information.objective_changed || warmstart_information.constraints_changed || warmstart_information.problem_changed) { this->save_hessian(subproblem); } - this->solve_subproblem(initial_point, direction, warmstart_information, subproblem); + this->solve_subproblem(subproblem, initial_point, direction, warmstart_information); } void BQPDSolver::solve_LP(const LagrangeNewtonSubproblem& subproblem, const Vector& initial_point, Direction& direction, @@ -89,7 +89,7 @@ namespace uno { DEBUG << "LP:\n"; } this->set_up_subproblem(warmstart_information, subproblem); - this->solve_subproblem(initial_point, direction, warmstart_information, subproblem); + this->solve_subproblem(subproblem, initial_point, direction, warmstart_information); } void BQPDSolver::set_up_subproblem(const WarmstartInformation& warmstart_information, const LagrangeNewtonSubproblem& subproblem) { @@ -151,8 +151,8 @@ namespace uno { WSC.ll += sizeof(intptr_t); } - void BQPDSolver::solve_subproblem(const Vector& initial_point, Direction& direction, - const WarmstartInformation& warmstart_information, const LagrangeNewtonSubproblem& subproblem) { + void BQPDSolver::solve_subproblem(const LagrangeNewtonSubproblem& subproblem, const Vector& initial_point, Direction& direction, + const WarmstartInformation& warmstart_information) { direction.primals = initial_point; const int n = static_cast(subproblem.number_variables); const int m = static_cast(subproblem.number_constraints); @@ -293,11 +293,7 @@ namespace uno { } } // namespace -void hessian_vector_product(int *n, const double x[], const double /*ws*/[], const int lws[], double v[]) { - for (int i = 0; i < *n; i++) { - v[i] = 0.; - } - +void hessian_vector_product([[maybe_unused]] int *n, const double x[], const double /*ws*/[], const int lws[], double v[]) { /* int footer_start = lws[0]; for (int i = 0; i < *n; i++) { @@ -317,13 +313,19 @@ void hessian_vector_product(int *n, const double x[], const double /*ws*/[], con std::copy(reinterpret_cast(lws), reinterpret_cast(lws) + sizeof(intptr_t), reinterpret_cast(&pointer_to_subproblem)); const uno::LagrangeNewtonSubproblem* subproblem = reinterpret_cast(pointer_to_subproblem); - uno::Vector my_x(*n); - for (size_t variable_index: uno::Range(*n)) { - my_x[variable_index] = x[variable_index]; + assert(static_cast(*n) == subproblem->number_variables && "BQPD and the subproblem do not have the same number of variables"); + for (size_t i = 0; i < subproblem->number_variables; i++) { + v[i] = 0.; + } + // convert x[] and v[] into Vector + // TODO improve that + uno::Vector primals(subproblem->number_variables); + for (size_t variable_index: uno::Range(subproblem->number_variables)) { + primals[variable_index] = x[variable_index]; } - uno::Vector result(*n); - subproblem->compute_hessian_vector_product(my_x, result); - for (size_t variable_index: uno::Range(*n)) { + uno::Vector result(subproblem->number_variables); + subproblem->compute_hessian_vector_product(primals, result); + for (size_t variable_index: uno::Range(subproblem->number_variables)) { v[variable_index] = result[variable_index]; } } \ No newline at end of file diff --git a/uno/solvers/BQPD/BQPDSolver.hpp b/uno/solvers/BQPD/BQPDSolver.hpp index 2d4d8644..a4db6336 100644 --- a/uno/solvers/BQPD/BQPDSolver.hpp +++ b/uno/solvers/BQPD/BQPDSolver.hpp @@ -84,8 +84,8 @@ namespace uno { RectangularMatrix constraint_jacobian; /*!< Sparse Jacobian of the constraints */ void set_up_subproblem(const WarmstartInformation& warmstart_information, const LagrangeNewtonSubproblem& subproblem); - void solve_subproblem(const Vector& initial_point, Direction& direction, const WarmstartInformation& warmstart_information, - const LagrangeNewtonSubproblem& subproblem); + void solve_subproblem(const LagrangeNewtonSubproblem& subproblem, const Vector& initial_point, Direction& direction, + const WarmstartInformation& warmstart_information); void categorize_constraints(size_t number_variables, Direction& direction); void save_hessian(const LagrangeNewtonSubproblem& subproblem); void save_gradients_to_local_format(size_t number_constraints, const SparseVector& linear_objective, diff --git a/uno/solvers/HiGHS/HiGHSSolver.cpp b/uno/solvers/HiGHS/HiGHSSolver.cpp index 5b43a4fb..5dd11331 100644 --- a/uno/solvers/HiGHS/HiGHSSolver.cpp +++ b/uno/solvers/HiGHS/HiGHSSolver.cpp @@ -5,12 +5,14 @@ #include "linear_algebra/SparseVector.hpp" #include "linear_algebra/Vector.hpp" #include "optimization/Direction.hpp" +#include "optimization/WarmstartInformation.hpp" #include "options/Options.hpp" #include "symbolic/VectorView.hpp" namespace uno { HiGHSSolver::HiGHSSolver(size_t number_variables, size_t number_constraints, size_t number_jacobian_nonzeros, size_t /*number_hessian_nonzeros*/, - const Options& options): LPSolver(), print_subproblem(options.get_bool("print_subproblem")) { + const Options& options): LPSolver(), print_subproblem(options.get_bool("print_subproblem")), + linear_objective(number_variables), constraints(number_constraints), constraint_jacobian(number_constraints, number_variables) { this->model.lp_.sense_ = ObjSense::kMinimize; this->model.lp_.offset_ = 0.; // the linear part of the objective is a dense vector @@ -37,8 +39,6 @@ namespace uno { } void HiGHSSolver::solve_subproblem(const LagrangeNewtonSubproblem& subproblem, Direction& direction) { - - /* HighsStatus passModel(const HighsInt num_col, const HighsInt num_row, const HighsInt num_nz, const HighsInt a_format, @@ -95,61 +95,67 @@ namespace uno { direction.subproblem_objective = info.objective_function_value; } + // build the LP in the HiGHS format void HiGHSSolver::build_linear_subproblem(const LagrangeNewtonSubproblem& subproblem, const WarmstartInformation& warmstart_information) { - /* this->model.lp_.num_col_ = static_cast(subproblem.number_variables); this->model.lp_.num_row_ = static_cast(subproblem.number_constraints); - // variable bounds - for (size_t variable_index = 0; variable_index < subproblem.number_variables; variable_index++) { - this->model.lp_.col_lower_[variable_index] = variables_lower_bounds[variable_index]; - this->model.lp_.col_upper_[variable_index] = variables_upper_bounds[variable_index]; - // reset the linear part of the objective - this->model.lp_.col_cost_[variable_index] = 0.; + if (warmstart_information.variable_bounds_changed) { + subproblem.set_direction_bounds(this->model.lp_.col_lower_, this->model.lp_.col_upper_); } // linear part of the objective - for (const auto [variable_index, value]: linear_objective) { - this->model.lp_.col_cost_[variable_index] = value; + if (warmstart_information.objective_changed) { + for (size_t variable_index: Range(subproblem.number_variables)) { + this->model.lp_.col_cost_[variable_index] = 0.; + } + subproblem.evaluate_objective_gradient(this->linear_objective); + for (const auto [variable_index, value]: this->linear_objective) { + this->model.lp_.col_cost_[variable_index] = value; + } } // constraint bounds - for (size_t constraint_index = 0; constraint_index < number_constraints; constraint_index++) { - this->model.lp_.row_lower_[constraint_index] = constraints_lower_bounds[constraint_index]; - this->model.lp_.row_upper_[constraint_index] = constraints_upper_bounds[constraint_index]; + if (warmstart_information.constraints_changed) { + subproblem.evaluate_constraints(this->constraints); + } + if (warmstart_information.constraint_bounds_changed) { + subproblem.set_constraint_bounds(this->constraints, this->model.lp_.row_lower_, this->model.lp_.row_upper_); } // constraint matrix - this->model.lp_.a_matrix_.value_.clear(); - this->model.lp_.a_matrix_.index_.clear(); - this->model.lp_.a_matrix_.start_.clear(); - - size_t number_nonzeros = 0; - this->model.lp_.a_matrix_.start_.emplace_back(number_nonzeros); - for (size_t constraint_index = 0; constraint_index < number_constraints; constraint_index++) { - for (const auto [variable_index, value]: constraint_jacobian[constraint_index]) { - this->model.lp_.a_matrix_.value_.emplace_back(value); - this->model.lp_.a_matrix_.index_.emplace_back(variable_index); - number_nonzeros++; - } + if (warmstart_information.constraints_changed) { + // TODO evaluate directly into this->model.lp_.a_matrix_ + subproblem.evaluate_constraint_jacobian(this->constraint_jacobian); + this->model.lp_.a_matrix_.value_.clear(); + this->model.lp_.a_matrix_.index_.clear(); + this->model.lp_.a_matrix_.start_.clear(); + size_t number_nonzeros = 0; this->model.lp_.a_matrix_.start_.emplace_back(number_nonzeros); + for (size_t constraint_index = 0; constraint_index < subproblem.number_constraints; constraint_index++) { + for (const auto [variable_index, value]: this->constraint_jacobian[constraint_index]) { + this->model.lp_.a_matrix_.value_.emplace_back(value); + this->model.lp_.a_matrix_.index_.emplace_back(variable_index); + number_nonzeros++; + } + this->model.lp_.a_matrix_.start_.emplace_back(number_nonzeros); + } } if (this->print_subproblem) { - DEBUG << "Linear objective part: "; print_vector(DEBUG, view(this->model.lp_.col_cost_, 0, number_variables)); + DEBUG << "Linear objective part: "; print_vector(DEBUG, view(this->model.lp_.col_cost_, 0, subproblem.number_variables)); DEBUG << "Jacobian:\n"; DEBUG << "J = "; print_vector(DEBUG, this->model.lp_.a_matrix_.value_); DEBUG << "with column start: "; print_vector(DEBUG, this->model.lp_.a_matrix_.start_); DEBUG << "and row index: "; print_vector(DEBUG, this->model.lp_.a_matrix_.index_); - for (size_t variable_index = 0; variable_index < number_variables; variable_index++) { + for (size_t variable_index = 0; variable_index < subproblem.number_variables; variable_index++) { DEBUG << "d" << variable_index << " in [" << this->model.lp_.col_lower_[variable_index] << ", " << this->model.lp_.col_upper_[variable_index] << "]\n"; } - for (size_t constraint_index = 0; constraint_index < number_constraints; constraint_index++) { + for (size_t constraint_index = 0; constraint_index < subproblem.number_constraints; constraint_index++) { DEBUG << "linearized c" << constraint_index << " in [" << this->model.lp_.row_lower_[constraint_index] << ", " << this->model.lp_.row_upper_[constraint_index]<< "]\n"; } } - */ } } // namespace \ No newline at end of file diff --git a/uno/solvers/HiGHS/HiGHSSolver.hpp b/uno/solvers/HiGHS/HiGHSSolver.hpp index 940080b7..7222515d 100644 --- a/uno/solvers/HiGHS/HiGHSSolver.hpp +++ b/uno/solvers/HiGHS/HiGHSSolver.hpp @@ -3,6 +3,9 @@ #include "solvers/LPSolver.hpp" #include "Highs.h" +#include "linear_algebra/RectangularMatrix.hpp" +#include "linear_algebra/SparseVector.hpp" +#include "linear_algebra/Vector.hpp" namespace uno { // forward declaration @@ -21,6 +24,10 @@ namespace uno { Highs highs_solver; const bool print_subproblem; + SparseVector linear_objective; + Vector constraints; + RectangularMatrix constraint_jacobian; + void build_linear_subproblem(const LagrangeNewtonSubproblem& subproblem, const WarmstartInformation& warmstart_information); void solve_subproblem(const LagrangeNewtonSubproblem& subproblem, Direction& direction); };