Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP: Make BQPD matrix-free #114

Closed
wants to merge 25 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
25 commits
Select commit Hold shift + click to select a range
f7acb00
BQPD: replaced kktalphac with alphac + removed unused functions in un…
cvanaret Nov 27, 2024
7da07df
Moved gdotx function (Hessian-vector product) from Fortran file to BQ…
cvanaret Nov 27, 2024
5c3fcc7
[CI] macOS build: try to remove the install of Fortran compiler
cvanaret Nov 27, 2024
d3da1c5
[CI] bad idea: restored the macOS build
cvanaret Nov 27, 2024
6785634
Added compute_hessian_vector_product() in Model and OptimizationProblem
cvanaret Nov 28, 2024
72f1e53
Forgot EvaluationErrors.hpp
cvanaret Nov 28, 2024
f651263
Commented out BQPD QP functional test
cvanaret Nov 29, 2024
e5db730
Temporarily commented out all BQPD tests
cvanaret Nov 29, 2024
0784d0f
Added LagrangeNewtonSubproblem object (wraps OptimizationProblem and …
cvanaret Nov 29, 2024
b37be2a
Temporarily commented out the HiGHS functional tests
cvanaret Nov 30, 2024
8a72953
HiGHS was not linked, forgot to do the changes to HiGHSSolver. It now…
cvanaret Nov 30, 2024
8e873d8
BQPD: replaced kktalphac with alphac + removed unused functions in un…
cvanaret Nov 27, 2024
a68f652
Moved gdotx function (Hessian-vector product) from Fortran file to BQ…
cvanaret Nov 27, 2024
c761b0a
[CI] macOS build: try to remove the install of Fortran compiler
cvanaret Nov 27, 2024
eedd0d1
[CI] bad idea: restored the macOS build
cvanaret Nov 27, 2024
7037434
Added compute_hessian_vector_product() in Model and OptimizationProblem
cvanaret Nov 28, 2024
c2cb165
Forgot EvaluationErrors.hpp
cvanaret Nov 28, 2024
ec08b39
Commented out BQPD QP functional test
cvanaret Nov 29, 2024
2c7a326
Temporarily commented out all BQPD tests
cvanaret Nov 29, 2024
50ff733
Added LagrangeNewtonSubproblem object (wraps OptimizationProblem and …
cvanaret Nov 29, 2024
cbde375
Temporarily commented out the HiGHS functional tests
cvanaret Nov 30, 2024
8cfad9b
HiGHS was not linked, forgot to do the changes to HiGHSSolver. It now…
cvanaret Nov 30, 2024
1f53646
Merge branch 'bqpd_matrix_free' of https://github.com/cvanaret/Uno in…
cvanaret Dec 2, 2024
36f62bc
Finally figured out how to link the Fortran compiler when linking sta…
cvanaret Dec 2, 2024
de23554
Merge branch 'main' into bqpd_matrix_free
cvanaret Dec 2, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 5 additions & 4 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -49,9 +49,10 @@ file(GLOB UNO_SOURCE_FILES
uno/ingredients/globalization_strategies/switching_methods/filter_methods/filters/*.cpp
uno/ingredients/globalization_strategies/switching_methods/funnel_methods/*.cpp
uno/ingredients/hessian_models/*.cpp
uno/ingredients/subproblems/*.cpp
uno/ingredients/subproblems/inequality_constrained_methods/*.cpp
uno/ingredients/subproblems/interior_point_methods/*.cpp
uno/ingredients/local_models/*.cpp
uno/ingredients/inequality_handling_methods/*.cpp
uno/ingredients/inequality_handling_methods/inequality_constrained_methods/*.cpp
uno/ingredients/inequality_handling_methods/interior_point_methods/*.cpp
uno/model/*.cpp
uno/optimization/*.cpp
uno/options/*.cpp
Expand Down Expand Up @@ -142,7 +143,7 @@ find_library(BQPD bqpd)
if(NOT BQPD)
message(WARNING "Optional library BQPD was not found.")
else()
list(APPEND UNO_SOURCE_FILES uno/solvers/BQPD/BQPDSolver.cpp uno/solvers/BQPD/wdotd.f)
list(APPEND UNO_SOURCE_FILES uno/solvers/BQPD/BQPDSolver.cpp)
list(APPEND TESTS_UNO_SOURCE_FILES unotest/functional_tests/BQPDSolverTests.cpp)
link_to_uno(bqpd ${BQPD})
endif()
Expand Down
232 changes: 129 additions & 103 deletions bindings/AMPL/AMPLModel.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -87,34 +87,6 @@ namespace uno {
ASL_free(&this->asl);
}

void AMPLModel::generate_variables() {
for (size_t variable_index: Range(this->number_variables)) {
this->variable_lower_bounds[variable_index] = (this->asl->i.LUv_ != nullptr) ? this->asl->i.LUv_[2*variable_index] : -INF<double>;
this->variable_upper_bounds[variable_index] = (this->asl->i.LUv_ != nullptr) ? this->asl->i.LUv_[2*variable_index + 1] : INF<double>;
if (this->variable_lower_bounds[variable_index] == this->variable_upper_bounds[variable_index]) {
WARNING << "Variable x" << variable_index << " has identical bounds\n";
this->fixed_variables.emplace_back(variable_index);
}
}
AMPLModel::determine_bounds_types(this->variable_lower_bounds, this->variable_upper_bounds, this->variable_status);
// figure out the bounded variables
for (size_t variable_index: Range(this->number_variables)) {
const BoundType status = this->get_variable_bound_type(variable_index);
if (status == BOUNDED_LOWER || status == BOUNDED_BOTH_SIDES) {
this->lower_bounded_variables.emplace_back(variable_index);
if (status == BOUNDED_LOWER) {
this->single_lower_bounded_variables.emplace_back(variable_index);
}
}
if (status == BOUNDED_UPPER || status == BOUNDED_BOTH_SIDES) {
this->upper_bounded_variables.emplace_back(variable_index);
if (status == BOUNDED_UPPER) {
this->single_upper_bounded_variables.emplace_back(variable_index);
}
}
}
}

double AMPLModel::evaluate_objective(const Vector<double>& x) const {
fint error_flag = 0;
double result = this->objective_sign * (*(this->asl)->p.Objval)(this->asl, 0, const_cast<double*>(x.data()), &error_flag);
Expand All @@ -129,8 +101,8 @@ namespace uno {
fint error_flag = 0;
// prevent ASL to crash by catching all evaluation errors
Jmp_buf err_jmp_uno;
asl->i.err_jmp_ = &err_jmp_uno;
asl->i.err_jmp1_ = &err_jmp_uno;
this->asl->i.err_jmp_ = &err_jmp_uno;
this->asl->i.err_jmp1_ = &err_jmp_uno;
if (setjmp(err_jmp_uno.jb)) {
error_flag = 1;
}
Expand All @@ -140,14 +112,15 @@ namespace uno {
throw GradientEvaluationError();
}

// create the Uno sparse vector
ograd* asl_variables_tmp = this->asl->i.Ograd_[0];
while (asl_variables_tmp != nullptr) {
const size_t variable_index = static_cast<size_t>(asl_variables_tmp->varno);
// fill the result
gradient.clear();
ograd* sparse_vector = this->asl->i.Ograd_[0];
while (sparse_vector != nullptr) {
const size_t variable_index = static_cast<size_t>(sparse_vector->varno);
// scale by the objective sign
const double partial_derivative = this->objective_sign*this->asl_gradient[variable_index];
gradient.insert(variable_index, partial_derivative);
asl_variables_tmp = asl_variables_tmp->next;
sparse_vector = sparse_vector->next;
}
}

Expand All @@ -162,7 +135,7 @@ namespace uno {
}
*/

void AMPLModel::evaluate_constraints(const Vector<double>& x, std::vector<double>& constraints) const {
void AMPLModel::evaluate_constraints(const Vector<double>& x, Vector<double>& constraints) const {
fint error_flag = 0;
(*(this->asl)->p.Conval)(this->asl, const_cast<double*>(x.data()), constraints.data(), &error_flag);
if (0 < error_flag) {
Expand All @@ -180,14 +153,14 @@ namespace uno {
throw GradientEvaluationError();
}

// construct the Uno sparse vector
// fill the result
gradient.clear();
cgrad* asl_variables_tmp = this->asl->i.Cgrad_[constraint_index];
cgrad* sparse_vector = this->asl->i.Cgrad_[constraint_index];
size_t sparse_asl_index = 0;
while (asl_variables_tmp != nullptr) {
const size_t variable_index = static_cast<size_t>(asl_variables_tmp->varno);
while (sparse_vector != nullptr) {
const size_t variable_index = static_cast<size_t>(sparse_vector->varno);
gradient.insert(variable_index, this->asl_gradient[sparse_asl_index]);
asl_variables_tmp = asl_variables_tmp->next;
sparse_vector = sparse_vector->next;
sparse_asl_index++;
}
}
Expand All @@ -199,68 +172,6 @@ namespace uno {
}
}

void AMPLModel::compute_lagrangian_hessian_sparsity() {
// compute the maximum number of nonzero elements, provided that all multipliers are non-zero
// int (*Sphset) (ASL*, SputInfo**, int nobj, int ow, int y, int uptri);
const int objective_number = -1;
const int upper_triangular = 1;
this->number_asl_hessian_nonzeros = static_cast<size_t>((*(this->asl)->p.Sphset)(this->asl, nullptr, objective_number, 1, 1, upper_triangular));
this->asl_hessian.reserve(this->number_asl_hessian_nonzeros);

// sparsity pattern
[[maybe_unused]] const fint* asl_column_start = this->asl->i.sputinfo_->hcolstarts;
// check that the column pointers are sorted in increasing order
assert(in_increasing_order(asl_column_start, this->number_variables + 1) && "AMPLModel::evaluate_lagrangian_hessian: column starts are not ordered");
}

size_t AMPLModel::number_objective_gradient_nonzeros() const {
return static_cast<size_t>(this->asl->i.nzo_);
}

size_t AMPLModel::number_jacobian_nonzeros() const {
return static_cast<size_t>(this->asl->i.nzc_);
}

size_t AMPLModel::number_hessian_nonzeros() const {
return this->number_asl_hessian_nonzeros;
}

const Collection<size_t>& AMPLModel::get_equality_constraints() const {
return this->equality_constraints_collection;
}

const Collection<size_t>& AMPLModel::get_inequality_constraints() const {
return this->inequality_constraints_collection;
}

const Collection<size_t>& AMPLModel::get_linear_constraints() const {
return this->linear_constraints_collection;
}

const SparseVector<size_t>& AMPLModel::get_slacks() const {
return this->slacks;
}

const Collection<size_t>& AMPLModel::get_single_lower_bounded_variables() const {
return this->single_lower_bounded_variables_collection;
}

const Collection<size_t>& AMPLModel::get_single_upper_bounded_variables() const {
return this->single_upper_bounded_variables_collection;
}

const Vector<size_t>& AMPLModel::get_fixed_variables() const {
return this->fixed_variables;
}

const Collection<size_t>& AMPLModel::get_lower_bounded_variables() const {
return this->lower_bounded_variables_collection;
}

const Collection<size_t>& AMPLModel::get_upper_bounded_variables() const {
return this->upper_bounded_variables_collection;
}

void AMPLModel::evaluate_lagrangian_hessian(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
SymmetricMatrix<size_t, double>& hessian) const {
assert(hessian.capacity() >= this->number_asl_hessian_nonzeros);
Expand Down Expand Up @@ -297,6 +208,31 @@ namespace uno {
this->asl->i.x_known = 0;
}

void AMPLModel::compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
Vector<double>& result) const {
fint error_flag = 0;
// prevent ASL to crash by catching all evaluation errors
Jmp_buf err_jmp_uno;
this->asl->i.err_jmp_ = &err_jmp_uno;
this->asl->i.err_jmp1_ = &err_jmp_uno;
if (setjmp(err_jmp_uno.jb)) {
error_flag = 1;
}

// scale by the objective sign
objective_multiplier *= this->objective_sign;
const int objective_number = -1;
// flip the signs of the multipliers: in AMPL, the Lagrangian is f + lambda.g, while Uno uses f - lambda.g
this->multipliers_with_flipped_sign = -multipliers;

// compute the Hessian-vector product
(this->asl->p.Hvcomp)(this->asl, result.data(), const_cast<double*>(x.data()), objective_number, &objective_multiplier,
const_cast<double*>(this->multipliers_with_flipped_sign.data()));
if (error_flag) {
throw HessianEvaluationError();
}
}

double AMPLModel::variable_lower_bound(size_t variable_index) const {
return this->variable_lower_bounds[variable_index];
}
Expand All @@ -309,6 +245,30 @@ namespace uno {
return this->variable_status[variable_index];
}

const Collection<size_t>& AMPLModel::get_lower_bounded_variables() const {
return this->lower_bounded_variables_collection;
}

const Collection<size_t>& AMPLModel::get_upper_bounded_variables() const {
return this->upper_bounded_variables_collection;
}

const SparseVector<size_t>& AMPLModel::get_slacks() const {
return this->slacks;
}

const Collection<size_t>& AMPLModel::get_single_lower_bounded_variables() const {
return this->single_lower_bounded_variables_collection;
}

const Collection<size_t>& AMPLModel::get_single_upper_bounded_variables() const {
return this->single_upper_bounded_variables_collection;
}

const Vector<size_t>& AMPLModel::get_fixed_variables() const {
return this->fixed_variables;
}

double AMPLModel::constraint_lower_bound(size_t constraint_index) const {
return this->constraint_lower_bounds[constraint_index];
}
Expand All @@ -325,6 +285,18 @@ namespace uno {
return this->constraint_status[constraint_index];
}

const Collection<size_t>& AMPLModel::get_equality_constraints() const {
return this->equality_constraints_collection;
}

const Collection<size_t>& AMPLModel::get_inequality_constraints() const {
return this->inequality_constraints_collection;
}

const Collection<size_t>& AMPLModel::get_linear_constraints() const {
return this->linear_constraints_collection;
}

// initial primal point
void AMPLModel::initial_primal_point(Vector<double>& x) const {
assert(x.size() >= this->number_variables);
Expand Down Expand Up @@ -379,6 +351,46 @@ namespace uno {
}
}

size_t AMPLModel::number_objective_gradient_nonzeros() const {
return static_cast<size_t>(this->asl->i.nzo_);
}

size_t AMPLModel::number_jacobian_nonzeros() const {
return static_cast<size_t>(this->asl->i.nzc_);
}

size_t AMPLModel::number_hessian_nonzeros() const {
return this->number_asl_hessian_nonzeros;
}

void AMPLModel::generate_variables() {
for (size_t variable_index: Range(this->number_variables)) {
this->variable_lower_bounds[variable_index] = (this->asl->i.LUv_ != nullptr) ? this->asl->i.LUv_[2*variable_index] : -INF<double>;
this->variable_upper_bounds[variable_index] = (this->asl->i.LUv_ != nullptr) ? this->asl->i.LUv_[2*variable_index + 1] : INF<double>;
if (this->variable_lower_bounds[variable_index] == this->variable_upper_bounds[variable_index]) {
WARNING << "Variable x" << variable_index << " has identical bounds\n";
this->fixed_variables.emplace_back(variable_index);
}
}
AMPLModel::determine_bounds_types(this->variable_lower_bounds, this->variable_upper_bounds, this->variable_status);
// figure out the bounded variables
for (size_t variable_index: Range(this->number_variables)) {
const BoundType status = this->get_variable_bound_type(variable_index);
if (status == BOUNDED_LOWER || status == BOUNDED_BOTH_SIDES) {
this->lower_bounded_variables.emplace_back(variable_index);
if (status == BOUNDED_LOWER) {
this->single_lower_bounded_variables.emplace_back(variable_index);
}
}
if (status == BOUNDED_UPPER || status == BOUNDED_BOTH_SIDES) {
this->upper_bounded_variables.emplace_back(variable_index);
if (status == BOUNDED_UPPER) {
this->single_upper_bounded_variables.emplace_back(variable_index);
}
}
}
}

void AMPLModel::generate_constraints() {
for (size_t constraint_index: Range(this->number_constraints)) {
this->constraint_lower_bounds[constraint_index] = (this->asl->i.LUrhs_ != nullptr) ? this->asl->i.LUrhs_[2*constraint_index] : -INF<double>;
Expand Down Expand Up @@ -407,6 +419,20 @@ namespace uno {
}
}

void AMPLModel::compute_lagrangian_hessian_sparsity() {
// compute the maximum number of nonzero elements, provided that all multipliers are non-zero
// int (*Sphset) (ASL*, SputInfo**, int nobj, int ow, int y, int uptri);
const int objective_number = -1;
const int upper_triangular = 1;
this->number_asl_hessian_nonzeros = static_cast<size_t>((*(this->asl)->p.Sphset)(this->asl, nullptr, objective_number, 1, 1, upper_triangular));
this->asl_hessian.reserve(this->number_asl_hessian_nonzeros);

// sparsity pattern
[[maybe_unused]] const fint* asl_column_start = this->asl->i.sputinfo_->hcolstarts;
// check that the column pointers are sorted in increasing order
assert(in_increasing_order(asl_column_start, this->number_variables + 1) && "AMPLModel::evaluate_lagrangian_hessian: column starts are not ordered");
}

void AMPLModel::determine_bounds_types(const std::vector<double>& lower_bounds, const std::vector<double>& upper_bounds, std::vector<BoundType>& status) {
assert(lower_bounds.size() == status.size());
assert(upper_bounds.size() == status.size());
Expand Down
4 changes: 3 additions & 1 deletion bindings/AMPL/AMPLModel.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -33,11 +33,13 @@ namespace uno {

[[nodiscard]] double evaluate_objective(const Vector<double>& x) const override;
void evaluate_objective_gradient(const Vector<double>& x, SparseVector<double>& gradient) const override;
void evaluate_constraints(const Vector<double>& x, std::vector<double>& constraints) const override;
void evaluate_constraints(const Vector<double>& x, Vector<double>& constraints) const override;
void evaluate_constraint_gradient(const Vector<double>& x, size_t constraint_index, SparseVector<double>& gradient) const override;
void evaluate_constraint_jacobian(const Vector<double>& x, RectangularMatrix<double>& constraint_jacobian) const override;
void evaluate_lagrangian_hessian(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
SymmetricMatrix<size_t, double>& hessian) const override;
void compute_hessian_vector_product(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
Vector<double>& result) const override;

[[nodiscard]] double variable_lower_bound(size_t variable_index) const override;
[[nodiscard]] double variable_upper_bound(size_t variable_index) const override;
Expand Down
4 changes: 2 additions & 2 deletions uno/Uno.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
#include "ingredients/globalization_mechanisms/GlobalizationMechanism.hpp"
#include "ingredients/globalization_mechanisms/GlobalizationMechanismFactory.hpp"
#include "ingredients/globalization_strategies/GlobalizationStrategyFactory.hpp"
#include "ingredients/subproblems/SubproblemFactory.hpp"
#include "ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.hpp"
#include "linear_algebra/Vector.hpp"
#include "model/Model.hpp"
#include "optimization/Iterate.hpp"
Expand Down Expand Up @@ -156,9 +156,9 @@ namespace uno {
void Uno::print_available_strategies() {
std::cout << "Available strategies:\n";
std::cout << "- Constraint relaxation strategies: " << join(ConstraintRelaxationStrategyFactory::available_strategies(), ", ") << '\n';
std::cout << "- Inequality-handling methods: " << join(InequalityHandlingMethodFactory::available_strategies(), ", ") << '\n';
std::cout << "- Globalization mechanisms: " << join(GlobalizationMechanismFactory::available_strategies(), ", ") << '\n';
std::cout << "- Globalization strategies: " << join(GlobalizationStrategyFactory::available_strategies(), ", ") << '\n';
std::cout << "- Subproblems: " << join(SubproblemFactory::available_strategies(), ", ") << '\n';
std::cout << "- QP solvers: " << join(QPSolverFactory::available_solvers(), ", ") << '\n';
std::cout << "- LP solvers: " << join(LPSolverFactory::available_solvers(), ", ") << '\n';
std::cout << "- Linear solvers: " << join(SymmetricIndefiniteLinearSolverFactory::available_solvers(), ", ") << '\n';
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,8 @@
#include "ingredients/globalization_strategies/GlobalizationStrategy.hpp"
#include "ingredients/globalization_strategies/GlobalizationStrategyFactory.hpp"
#include "optimization/Direction.hpp"
#include "ingredients/subproblems/Subproblem.hpp"
#include "ingredients/subproblems/SubproblemFactory.hpp"
#include "ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp"
#include "ingredients/inequality_handling_methods/InequalityHandlingMethodFactory.hpp"
#include "linear_algebra/SymmetricMatrix.hpp"
#include "model/Model.hpp"
#include "optimization/Iterate.hpp"
Expand All @@ -22,8 +22,8 @@ namespace uno {
size_t number_objective_gradient_nonzeros, size_t number_jacobian_nonzeros, size_t number_hessian_nonzeros, const Options& options):
model(model),
globalization_strategy(GlobalizationStrategyFactory::create(options.get_string("globalization_strategy"), options)),
subproblem(SubproblemFactory::create(number_variables, number_constraints, number_objective_gradient_nonzeros, number_jacobian_nonzeros,
number_hessian_nonzeros, options)),
subproblem(InequalityHandlingMethodFactory::create(number_variables, number_constraints, number_objective_gradient_nonzeros,
number_jacobian_nonzeros, number_hessian_nonzeros, options)),
progress_norm(norm_from_string(options.get_string("progress_norm"))),
residual_norm(norm_from_string(options.get_string("residual_norm"))),
residual_scaling_threshold(options.get_double("residual_scaling_threshold")),
Expand Down
Loading
Loading