From 4e122a110094c19f31ac29e86de6e002af3f25dc Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Sun, 1 Dec 2024 20:57:15 +0100 Subject: [PATCH] Model: removed flag fixed_hessian_sparsity + cleanup --- bindings/AMPL/AMPLModel.cpp | 50 +++++-------------- bindings/AMPL/AMPLModel.hpp | 3 +- .../SymmetricIndefiniteLinearSystem.hpp | 4 +- uno/model/Model.hpp | 3 -- 4 files changed, 16 insertions(+), 44 deletions(-) diff --git a/bindings/AMPL/AMPLModel.cpp b/bindings/AMPL/AMPLModel.cpp index 06f09a2e..f7af1062 100644 --- a/bindings/AMPL/AMPLModel.cpp +++ b/bindings/AMPL/AMPLModel.cpp @@ -79,8 +79,8 @@ namespace uno { this->linear_constraints.reserve(this->number_constraints); this->generate_constraints(); - // compute number of nonzeros in the Lagrangian Hessian - this->set_number_hessian_nonzeros(); + // compute sparsity pattern and number of nonzeros of Lagrangian Hessian + this->compute_lagrangian_hessian_sparsity(); } AMPLModel::~AMPLModel() { @@ -199,13 +199,18 @@ namespace uno { } } - void AMPLModel::set_number_hessian_nonzeros() { + 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((*(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 { @@ -252,49 +257,20 @@ namespace uno { return this->upper_bounded_variables_collection; } - bool are_all_zeros(const Vector& multipliers) { - return std::all_of(multipliers.begin(), multipliers.end(), [](double xj) { - return xj == 0.; - }); - } - - size_t AMPLModel::compute_hessian_number_nonzeros(double objective_multiplier, const Vector& multipliers) const { - // compute the sparsity - const int objective_number = -1; - const int upper_triangular = 1; - const bool all_zeros_multipliers = are_all_zeros(multipliers); - int number_nonzeros = (*(this->asl)->p.Sphset)(this->asl, nullptr, objective_number, (objective_multiplier != 0.), - not all_zeros_multipliers, upper_triangular); - return static_cast(number_nonzeros); - } - void AMPLModel::evaluate_lagrangian_hessian(const Vector& x, double objective_multiplier, const Vector& multipliers, SymmetricMatrix& hessian) const { + assert(hessian.capacity() >= this->number_asl_hessian_nonzeros); + // register the vector of variables (*(this->asl)->p.Xknown)(this->asl, const_cast(x.data()), nullptr); - // scale by the objective sign - objective_multiplier *= this->objective_sign; - - // compute the number of nonzeros - [[maybe_unused]] const size_t number_nonzeros = this->fixed_hessian_sparsity ? this->number_asl_hessian_nonzeros : - this->compute_hessian_number_nonzeros(objective_multiplier, multipliers); - assert(hessian.capacity() >= number_nonzeros); - // evaluate the Hessian: store the matrix in a preallocated array this->asl_hessian const int objective_number = -1; + objective_multiplier *= this->objective_sign; // 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; - if (this->fixed_hessian_sparsity) { - (*(this->asl)->p.Sphes)(this->asl, nullptr, const_cast(this->asl_hessian.data()), objective_number, &objective_multiplier, - const_cast(this->multipliers_with_flipped_sign.data())); - } - else { - double* objective_multiplier_pointer = (objective_multiplier != 0.) ? &objective_multiplier : nullptr; - const bool all_zeros_multipliers = are_all_zeros(multipliers); - (*(this->asl)->p.Sphes)(this->asl, nullptr, const_cast(this->asl_hessian.data()), objective_number, objective_multiplier_pointer, - all_zeros_multipliers ? nullptr : const_cast(this->multipliers_with_flipped_sign.data())); - } + (*(this->asl)->p.Sphes)(this->asl, nullptr, const_cast(this->asl_hessian.data()), objective_number, &objective_multiplier, + const_cast(this->multipliers_with_flipped_sign.data())); // generate the sparsity pattern in the right sparse format const fint* asl_column_start = this->asl->i.sputinfo_->hcolstarts; diff --git a/bindings/AMPL/AMPLModel.hpp b/bindings/AMPL/AMPLModel.hpp index 370d5928..da9aecce 100644 --- a/bindings/AMPL/AMPLModel.hpp +++ b/bindings/AMPL/AMPLModel.hpp @@ -107,8 +107,7 @@ namespace uno { void generate_variables(); void generate_constraints(); - void set_number_hessian_nonzeros(); - [[nodiscard]] size_t compute_hessian_number_nonzeros(double objective_multiplier, const Vector& multipliers) const; + void compute_lagrangian_hessian_sparsity(); static void determine_bounds_types(const std::vector& lower_bounds, const std::vector& upper_bounds, std::vector& status); }; diff --git a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp index 00b4e5f9..d06c2847 100644 --- a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp +++ b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp @@ -89,11 +89,11 @@ namespace uno { } template - void SymmetricIndefiniteLinearSystem::factorize_matrix(const Model& model, + void SymmetricIndefiniteLinearSystem::factorize_matrix(const Model& /*model*/, DirectSymmetricIndefiniteLinearSolver& linear_solver) { // compute the symbolic factorization only when: // the problem has a non-constant augmented system (ie is not an LP or a QP) or it is the first factorization - if (true || this->number_factorizations == 0 || not model.fixed_hessian_sparsity) { + if (true || this->number_factorizations == 0) { linear_solver.do_symbolic_factorization(this->matrix); } linear_solver.do_numerical_factorization(this->matrix); diff --git a/uno/model/Model.hpp b/uno/model/Model.hpp index de5f12af..389c53ea 100644 --- a/uno/model/Model.hpp +++ b/uno/model/Model.hpp @@ -44,9 +44,6 @@ namespace uno { const size_t number_constraints; /*!< Number of constraints */ const double objective_sign; /*!< Sign of the objective function (1: minimization, -1: maximization) */ - // Hessian - const bool fixed_hessian_sparsity{true}; - [[nodiscard]] virtual double evaluate_objective(const Vector& x) const = 0; virtual void evaluate_objective_gradient(const Vector& x, SparseVector& gradient) const = 0; virtual void evaluate_constraints(const Vector& x, std::vector& constraints) const = 0;