Merge branch 'main' into ma27

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<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 {
@@ -255,26 +260,17 @@ namespace uno {
    const Collection<size_t>& AMPLModel::get_upper_bounded_variables() const {
       return this->upper_bounded_variables_collection;
    }
-
-   bool are_all_zeros(const Vector<double>& multipliers) {
-      return std::all_of(multipliers.begin(), multipliers.end(), [](double xj) {
-         return xj == 0.;
-      });
-   }
-
+
    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);
+
       // register the vector of variables
       (*(this->asl)->p.Xknown)(this->asl, const_cast<double*>(x.data()), nullptr);
 
-      // scale by the objective sign
-      objective_multiplier *= this->objective_sign;
-
-      // compute the number of nonzeros
-      assert(hessian.capacity() >= this->number_asl_hessian_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;
       (*(this->asl)->p.Sphes)(this->asl, nullptr, const_cast<double*>(this->asl_hessian.data()), objective_number, &objective_multiplier,

bindings/AMPL/AMPLModel.hpp

@@ -107,7 +107,7 @@ namespace uno {
       void generate_variables();
       void generate_constraints();
 
-      void set_number_hessian_nonzeros();
+      void compute_lagrangian_hessian_sparsity();
       static void determine_bounds_types(const std::vector<double>& lower_bounds, const std::vector<double>& upper_bounds, std::vector<BoundType>& status);
    };