Skip to content

Commit

Permalink
HiGHS was not linked, forgot to do the changes to HiGHSSolver. It now…
Browse files Browse the repository at this point in the history 
… takes a subproblem and evaluates the functions locally
  • Loading branch information
@cvanaret
cvanaret committed Nov 30, 2024
1 parent b37be2a commit 8a72953
Show file tree
Hide file tree
Showing 4 changed files with 80 additions and 75 deletions.
32 changes: 17 additions & 15 deletions uno/solvers/BQPD/BQPDSolver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<double>& initial_point, Direction& direction,
Expand All @@ -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) {
Expand Down Expand Up @@ -151,8 +151,8 @@ namespace uno {
WSC.ll += sizeof(intptr_t);
}

void BQPDSolver::solve_subproblem(const Vector<double>& initial_point, Direction& direction,
const WarmstartInformation& warmstart_information, const LagrangeNewtonSubproblem& subproblem) {
void BQPDSolver::solve_subproblem(const LagrangeNewtonSubproblem& subproblem, const Vector<double>& initial_point, Direction& direction,
const WarmstartInformation& warmstart_information) {
direction.primals = initial_point;
const int n = static_cast<int>(subproblem.number_variables);
const int m = static_cast<int>(subproblem.number_constraints);
Expand Down Expand Up @@ -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++) {
Expand All @@ -317,13 +313,19 @@ void hessian_vector_product(int *n, const double x[], const double /*ws*/[], con
std::copy(reinterpret_cast<const char *>(lws), reinterpret_cast<const char *>(lws) + sizeof(intptr_t), reinterpret_cast<char *>(&pointer_to_subproblem));
const uno::LagrangeNewtonSubproblem* subproblem = reinterpret_cast<const uno::LagrangeNewtonSubproblem*>(pointer_to_subproblem);

uno::Vector<double> my_x(*n);
for (size_t variable_index: uno::Range(*n)) {
my_x[variable_index] = x[variable_index];
assert(static_cast<size_t>(*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<double>
// TODO improve that
uno::Vector<double> primals(subproblem->number_variables);
for (size_t variable_index: uno::Range(subproblem->number_variables)) {
primals[variable_index] = x[variable_index];
}
uno::Vector<double> result(*n);
subproblem->compute_hessian_vector_product(my_x, result);
for (size_t variable_index: uno::Range(*n)) {
uno::Vector<double> 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];
}
}
4 changes: 2 additions & 2 deletions uno/solvers/BQPD/BQPDSolver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,8 +84,8 @@ namespace uno {
RectangularMatrix<double> constraint_jacobian; /*!< Sparse Jacobian of the constraints */

void set_up_subproblem(const WarmstartInformation& warmstart_information, const LagrangeNewtonSubproblem& subproblem);
void solve_subproblem(const Vector<double>& initial_point, Direction& direction, const WarmstartInformation& warmstart_information,
const LagrangeNewtonSubproblem& subproblem);
void solve_subproblem(const LagrangeNewtonSubproblem& subproblem, const Vector<double>& 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<double>& linear_objective,
Expand Down
98 changes: 50 additions & 48 deletions uno/solvers/HiGHS/HiGHSSolver.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,15 +1,18 @@
#include <cassert>
#include "HiGHSSolver.hpp"
#include "ingredients/local_models/LagrangeNewtonSubproblem.hpp"
#include "linear_algebra/RectangularMatrix.hpp"
#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
Expand All @@ -29,66 +32,78 @@ namespace uno {
this->highs_solver.setOptionValue("output_flag", "false");
}

void HiGHSSolver::build_linear_subproblem(size_t number_variables, size_t number_constraints, const std::vector<double>& variables_lower_bounds,
const std::vector<double>& variables_upper_bounds, const std::vector<double>& constraints_lower_bounds,
const std::vector<double>& constraints_upper_bounds, const SparseVector<double>& linear_objective,
const RectangularMatrix<double>& constraint_jacobian) {
this->model.lp_.num_col_ = static_cast<HighsInt>(number_variables);
this->model.lp_.num_row_ = static_cast<HighsInt>(number_constraints);
void HiGHSSolver::solve_LP(const LagrangeNewtonSubproblem& subproblem, const Vector<double>& /*initial_point*/, Direction& direction,
const WarmstartInformation& warmstart_information) {
this->set_up_subproblem(subproblem, warmstart_information);
this->solve_subproblem(subproblem, direction);
}

// build the LP in the HiGHS format
void HiGHSSolver::set_up_subproblem(const LagrangeNewtonSubproblem& subproblem, const WarmstartInformation& warmstart_information) {
this->model.lp_.num_col_ = static_cast<HighsInt>(subproblem.number_variables);
this->model.lp_.num_row_ = static_cast<HighsInt>(subproblem.number_constraints);

// variable bounds
for (size_t variable_index = 0; variable_index < 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";
}
}
}

void HiGHSSolver::solve_subproblem(Direction& direction, size_t number_variables, size_t number_constraints) {
void HiGHSSolver::solve_subproblem(const LagrangeNewtonSubproblem& subproblem, Direction& direction) {
// solve the LP
HighsStatus return_status = this->highs_solver.passModel(this->model);
assert(return_status == HighsStatus::kOk);
Expand Down Expand Up @@ -116,7 +131,7 @@ namespace uno {
direction.status = SubproblemStatus::OPTIMAL;
const HighsSolution& solution = this->highs_solver.getSolution();
// read the primal solution and bound dual solution
for (size_t variable_index = 0; variable_index < number_variables; variable_index++) {
for (size_t variable_index: Range(subproblem.number_variables)) {
direction.primals[variable_index] = solution.col_value[variable_index];
const double bound_multiplier = solution.col_dual[variable_index];
if (0. < bound_multiplier) {
Expand All @@ -127,23 +142,10 @@ namespace uno {
}
}
// read the dual solution
for (size_t constraint_index = 0; constraint_index < number_constraints; constraint_index++) {
for (size_t constraint_index: Range(subproblem.number_constraints)) {
direction.multipliers.constraints[constraint_index] = solution.row_dual[constraint_index];
}
const HighsInfo& info = this->highs_solver.getInfo();
direction.subproblem_objective = info.objective_function_value;
}

void HiGHSSolver::solve_LP(size_t number_variables, size_t number_constraints, const std::vector<double>& variables_lower_bounds,
const std::vector<double>& variables_upper_bounds, const std::vector<double>& constraints_lower_bounds,
const std::vector<double>& constraints_upper_bounds, const SparseVector<double>& linear_objective,
const RectangularMatrix<double>& constraint_jacobian, const Vector<double>& /*initial_point*/, Direction& direction,
const WarmstartInformation& /*warmstart_information*/, const OptimizationProblem& problem) {
// build the LP in the HiGHS format
this->build_linear_subproblem(number_variables, number_constraints, variables_lower_bounds, variables_upper_bounds, constraints_lower_bounds,
constraints_upper_bounds, linear_objective, constraint_jacobian);

// solve the LP
this->solve_subproblem(direction, number_variables, number_constraints);
}
} // namespace
21 changes: 11 additions & 10 deletions uno/solvers/HiGHS/HiGHSSolver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -13,22 +16,20 @@ namespace uno {
HiGHSSolver(size_t number_variables, size_t number_constraints, size_t number_jacobian_nonzeros, size_t number_hessian_nonzeros,
const Options& options);

void solve_LP(size_t number_variables, size_t number_constraints, const std::vector<double>& variables_lower_bounds,
const std::vector<double>& variables_upper_bounds, const std::vector<double>& constraints_lower_bounds,
const std::vector<double>& constraints_upper_bounds, const SparseVector<double>& linear_objective,
const RectangularMatrix<double>& constraint_jacobian, const Vector<double>& initial_point, Direction& direction,
const WarmstartInformation& warmstart_information, const OptimizationProblem& problem) override;
void solve_LP(const LagrangeNewtonSubproblem& subproblem, const Vector<double>& initial_point, Direction& direction,
const WarmstartInformation& warmstart_information) override;

protected:
HighsModel model;
Highs highs_solver;
const bool print_subproblem;

void build_linear_subproblem(size_t number_variables, size_t number_constraints, const std::vector<double>& variables_lower_bounds,
const std::vector<double>& variables_upper_bounds, const std::vector<double>& constraints_lower_bounds,
const std::vector<double>& constraints_upper_bounds, const SparseVector<double>& linear_objective,
const RectangularMatrix<double>& constraint_jacobian, const OptimizationProblem& problem);
void solve_subproblem(Direction& direction, size_t number_variables, size_t number_constraints);
SparseVector<double> linear_objective;
Vector<double> constraints;
RectangularMatrix<double> constraint_jacobian;

void set_up_subproblem(const LagrangeNewtonSubproblem& subproblem, const WarmstartInformation& warmstart_information);
void solve_subproblem(const LagrangeNewtonSubproblem& subproblem, Direction& direction);
};
} // namespace

Expand Down

0 comments on commit 8a72953

Please sign in to comment.