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 Dec 2, 2024
1 parent cbde375 commit 8cfad9b
Show file tree
Hide file tree
Showing 4 changed files with 63 additions and 48 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
68 changes: 37 additions & 31 deletions uno/solvers/HiGHS/HiGHSSolver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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,
Expand Down Expand Up @@ -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<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 < 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
7 changes: 7 additions & 0 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 @@ -21,6 +24,10 @@ namespace uno {
Highs highs_solver;
const bool print_subproblem;

SparseVector<double> linear_objective;
Vector<double> constraints;
RectangularMatrix<double> constraint_jacobian;

void build_linear_subproblem(const LagrangeNewtonSubproblem& subproblem, const WarmstartInformation& warmstart_information);
void solve_subproblem(const LagrangeNewtonSubproblem& subproblem, Direction& direction);
};
Expand Down

0 comments on commit 8cfad9b

Please sign in to comment.