Skip to content

Commit

Permalink
Merge pull request #127 from cvanaret/hessian
Browse files Browse the repository at this point in the history 
Move Hessian evaluation after first derivatives
  • Loading branch information
@cvanaret
cvanaret authored Dec 5, 2024
2 parents 3a92daa + 51ba604 commit 57393ad
Show file tree
Hide file tree
Showing 2 changed files with 97 additions and 97 deletions.
186 changes: 93 additions & 93 deletions bindings/AMPL/AMPLModel.cpp
View file
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 Down Expand Up @@ -199,74 +171,12 @@ 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,
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);
//(*(this->asl)->p.Xknown)(this->asl, const_cast<double*>(x.data()), nullptr);

// evaluate the Hessian: store the matrix in a preallocated array this->asl_hessian
const int objective_number = -1;
Expand Down Expand Up @@ -294,7 +204,7 @@ namespace uno {
hessian.finalize_column(column_index);
}
// unregister the vector of variables
this->asl->i.x_known = 0;
//this->asl->i.x_known = 0;
}

double AMPLModel::variable_lower_bound(size_t variable_index) const {
Expand All @@ -309,6 +219,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 +259,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 +325,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 +393,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
8 changes: 4 additions & 4 deletions uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,10 +45,6 @@ namespace uno {

void QPSubproblem::evaluate_functions(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate,
const Multipliers& current_multipliers, const WarmstartInformation& warmstart_information) {
// Lagrangian Hessian
if (warmstart_information.objective_changed || warmstart_information.constraints_changed) {
this->hessian_model->evaluate(statistics, problem, current_iterate.primals, current_multipliers.constraints);
}
// objective gradient, constraints and constraint Jacobian
if (warmstart_information.objective_changed) {
problem.evaluate_objective_gradient(current_iterate, this->objective_gradient);
Expand All @@ -57,6 +53,10 @@ namespace uno {
problem.evaluate_constraints(current_iterate, this->constraints);
problem.evaluate_constraint_jacobian(current_iterate, this->constraint_jacobian);
}
// Lagrangian Hessian
if (warmstart_information.objective_changed || warmstart_information.constraints_changed) {
this->hessian_model->evaluate(statistics, problem, current_iterate.primals, current_multipliers.constraints);
}
}

void QPSubproblem::solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, const Multipliers& current_multipliers,
Expand Down

0 comments on commit 57393ad

Please sign in to comment.