From de602f0bef58f60fe6c2434de36461573968ec68 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Wed, 6 Nov 2024 15:06:49 +0100 Subject: [PATCH] Refactoring: split all the classes into hpp and cpp --- CMakeLists.txt | 2 + .../ConstraintRelaxationStrategy.cpp | 2 +- .../FeasibilityRestoration.cpp | 2 +- .../FeasibilityRestoration.hpp | 1 + .../l1Relaxation.cpp | 4 +- .../l1Relaxation.hpp | 2 + .../GlobalizationMechanism.hpp | 2 +- .../l1MeritFunction.cpp | 2 +- .../hessian_models/HessianModelFactory.cpp | 27 ++ .../hessian_models/HessianModelFactory.hpp | 25 +- .../InequalityConstrainedMethod.cpp | 3 +- .../LPSubproblem.cpp | 4 +- .../LPSubproblem.hpp | 1 + .../QPSubproblem.cpp | 4 +- .../QPSubproblem.hpp | 1 + .../PrimalDualInteriorPointSubproblem.cpp | 11 +- uno/model/BoundRelaxedModel.cpp | 25 ++ uno/model/BoundRelaxedModel.hpp | 23 +- uno/model/FixedBoundsConstraintsModel.cpp | 219 ++++++++++++++ uno/model/FixedBoundsConstraintsModel.hpp | 229 +++----------- .../HomogeneousEqualityConstrainedModel.cpp | 214 +++++++++++++ .../HomogeneousEqualityConstrainedModel.hpp | 178 ++--------- uno/model/ScaledModel.cpp | 180 +++++++++++ uno/model/ScaledModel.hpp | 140 ++------- .../Direction.cpp | 0 .../Direction.hpp | 7 +- uno/optimization/Multipliers.cpp | 32 ++ uno/optimization/Multipliers.hpp | 26 -- uno/optimization/WarmstartInformation.cpp | 47 +++ uno/optimization/WarmstartInformation.hpp | 40 --- uno/preprocessing/Preprocessing.cpp | 2 +- uno/reformulation/ElasticVariables.cpp | 46 +++ uno/reformulation/ElasticVariables.hpp | 43 +-- uno/reformulation/OptimalityProblem.cpp | 90 ++++++ uno/reformulation/OptimalityProblem.hpp | 84 ------ uno/reformulation/OptimizationProblem.cpp | 33 ++ uno/reformulation/OptimizationProblem.hpp | 27 -- uno/reformulation/l1RelaxedProblem.cpp | 281 ++++++++++++++++++ uno/reformulation/l1RelaxedProblem.hpp | 275 ----------------- uno/solvers/BQPD/BQPDSolver.cpp | 16 +- uno/solvers/HiGHS/HiGHSSolver.cpp | 2 +- uno/solvers/LPSolverFactory.cpp | 56 ++++ uno/solvers/LPSolverFactory.hpp | 54 +--- uno/solvers/QPSolver.hpp | 5 +- uno/solvers/QPSolverFactory.cpp | 46 +++ uno/solvers/QPSolverFactory.hpp | 42 +-- ...SymmetricIndefiniteLinearSolverFactory.cpp | 74 +++++ ...SymmetricIndefiniteLinearSolverFactory.hpp | 71 +---- uno/tools/Logger.cpp | 33 ++ uno/tools/Logger.hpp | 51 +--- unotest/HiGHSSolverTests.cpp | 2 +- 51 files changed, 1569 insertions(+), 1217 deletions(-) create mode 100644 uno/ingredients/hessian_models/HessianModelFactory.cpp create mode 100644 uno/model/BoundRelaxedModel.cpp create mode 100644 uno/model/FixedBoundsConstraintsModel.cpp create mode 100644 uno/model/HomogeneousEqualityConstrainedModel.cpp create mode 100644 uno/model/ScaledModel.cpp rename uno/{ingredients/subproblems => optimization}/Direction.cpp (100%) rename uno/{ingredients/subproblems => optimization}/Direction.hpp (86%) create mode 100644 uno/optimization/Multipliers.cpp create mode 100644 uno/optimization/WarmstartInformation.cpp create mode 100644 uno/reformulation/ElasticVariables.cpp create mode 100644 uno/reformulation/OptimalityProblem.cpp create mode 100644 uno/reformulation/OptimizationProblem.cpp create mode 100644 uno/reformulation/l1RelaxedProblem.cpp create mode 100644 uno/solvers/LPSolverFactory.cpp create mode 100644 uno/solvers/QPSolverFactory.cpp create mode 100644 uno/solvers/SymmetricIndefiniteLinearSolverFactory.cpp create mode 100644 uno/tools/Logger.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 7704ea44..274563a7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -56,6 +56,8 @@ file(GLOB UNO_SOURCE_FILES uno/optimization/*.cpp uno/options/*.cpp uno/preprocessing/*.cpp + uno/reformulation/*.cpp + uno/solvers/*.cpp uno/tools/*.cpp ) diff --git a/uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.cpp b/uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.cpp index 38545965..e8b5853d 100644 --- a/uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.cpp @@ -4,7 +4,7 @@ #include "ConstraintRelaxationStrategy.hpp" #include "ingredients/globalization_strategies/GlobalizationStrategy.hpp" #include "ingredients/globalization_strategies/GlobalizationStrategyFactory.hpp" -#include "ingredients/subproblems/Direction.hpp" +#include "optimization/Direction.hpp" #include "ingredients/subproblems/Subproblem.hpp" #include "ingredients/subproblems/SubproblemFactory.hpp" #include "linear_algebra/SymmetricMatrix.hpp" diff --git a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp index a66c095a..4b9c4a42 100644 --- a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp @@ -4,7 +4,7 @@ #include #include "FeasibilityRestoration.hpp" #include "ingredients/globalization_strategies/GlobalizationStrategy.hpp" -#include "ingredients/subproblems/Direction.hpp" +#include "optimization/Direction.hpp" #include "ingredients/subproblems/Subproblem.hpp" #include "linear_algebra/SymmetricIndefiniteLinearSystem.hpp" #include "model/Model.hpp" diff --git a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.hpp b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.hpp index ff73052f..b4dbf130 100644 --- a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.hpp +++ b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.hpp @@ -6,6 +6,7 @@ #include #include "ConstraintRelaxationStrategy.hpp" +#include "ingredients/globalization_strategies/ProgressMeasures.hpp" #include "reformulation/OptimalityProblem.hpp" #include "reformulation/l1RelaxedProblem.hpp" diff --git a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp index ef19172d..e5364103 100644 --- a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp @@ -4,7 +4,7 @@ #include #include "l1Relaxation.hpp" #include "ingredients/globalization_strategies/GlobalizationStrategy.hpp" -#include "ingredients/subproblems/Direction.hpp" +#include "optimization/Direction.hpp" #include "ingredients/subproblems/Subproblem.hpp" #include "optimization/Iterate.hpp" #include "optimization/WarmstartInformation.hpp" @@ -166,7 +166,7 @@ namespace uno { DEBUG << "Further aggressively decrease the penalty parameter to " << this->penalty_parameter << '\n'; } else { - DEBUG << RED << "l1Relaxation: all multipliers are almost 0. The penalty parameter won't be decreased" << RESET << '\n'; + DEBUG << "l1Relaxation: all multipliers are almost 0. The penalty parameter won't be decreased\n"; } } diff --git a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp index bcbdb5f0..759986f3 100644 --- a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp +++ b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp @@ -6,6 +6,8 @@ #include #include "ConstraintRelaxationStrategy.hpp" +#include "ingredients/globalization_strategies/ProgressMeasures.hpp" +#include "optimization/Multipliers.hpp" #include "reformulation/l1RelaxedProblem.hpp" namespace uno { diff --git a/uno/ingredients/globalization_mechanisms/GlobalizationMechanism.hpp b/uno/ingredients/globalization_mechanisms/GlobalizationMechanism.hpp index 02b2f915..0fb8af8e 100644 --- a/uno/ingredients/globalization_mechanisms/GlobalizationMechanism.hpp +++ b/uno/ingredients/globalization_mechanisms/GlobalizationMechanism.hpp @@ -4,7 +4,7 @@ #ifndef UNO_GLOBALIZATIONMECHANISM_H #define UNO_GLOBALIZATIONMECHANISM_H -#include "ingredients/subproblems/Direction.hpp" +#include "optimization/Direction.hpp" namespace uno { // forward declarations diff --git a/uno/ingredients/globalization_strategies/l1MeritFunction.cpp b/uno/ingredients/globalization_strategies/l1MeritFunction.cpp index 508e8fea..126dcdd5 100644 --- a/uno/ingredients/globalization_strategies/l1MeritFunction.cpp +++ b/uno/ingredients/globalization_strategies/l1MeritFunction.cpp @@ -30,7 +30,7 @@ namespace uno { double constrained_predicted_reduction = l1MeritFunction::constrained_merit_function(predicted_reduction, objective_multiplier); DEBUG << "Constrained predicted reduction: " << constrained_predicted_reduction << '\n'; if (constrained_predicted_reduction <= 0.) { - WARNING << YELLOW << "The direction is not a descent direction for the merit function. You should decrease the penalty parameter.\n" << RESET; + WARNING << "The direction is not a descent direction for the merit function. You should decrease the penalty parameter.\n"; } // compute current exact penalty const double current_merit_value = l1MeritFunction::constrained_merit_function(current_progress, objective_multiplier); diff --git a/uno/ingredients/hessian_models/HessianModelFactory.cpp b/uno/ingredients/hessian_models/HessianModelFactory.cpp new file mode 100644 index 00000000..a7ee8cde --- /dev/null +++ b/uno/ingredients/hessian_models/HessianModelFactory.cpp @@ -0,0 +1,27 @@ +// Copyright (c) 2018-2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "HessianModelFactory.hpp" +#include "HessianModel.hpp" +#include "ConvexifiedHessian.hpp" +#include "ExactHessian.hpp" +#include "ZeroHessian.hpp" +#include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp" + +namespace uno { + std::unique_ptr HessianModelFactory::create(const std::string& hessian_model, size_t dimension, size_t maximum_number_nonzeros, + bool convexify, const Options& options) { + if (hessian_model == "exact") { + if (convexify) { + return std::make_unique(dimension, maximum_number_nonzeros + dimension, options); + } + else { + return std::make_unique(dimension, maximum_number_nonzeros, options); + } + } + else if (hessian_model == "zero") { + return std::make_unique(dimension, options); + } + throw std::invalid_argument("Hessian model " + hessian_model + " does not exist"); + } +} // namespace \ No newline at end of file diff --git a/uno/ingredients/hessian_models/HessianModelFactory.hpp b/uno/ingredients/hessian_models/HessianModelFactory.hpp index 8a134b9a..04024525 100644 --- a/uno/ingredients/hessian_models/HessianModelFactory.hpp +++ b/uno/ingredients/hessian_models/HessianModelFactory.hpp @@ -6,34 +6,17 @@ #include #include -#include "HessianModel.hpp" -#include "ConvexifiedHessian.hpp" -#include "ExactHessian.hpp" -#include "ZeroHessian.hpp" -#include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp" namespace uno { + // forward declarations + class HessianModel; + class Options; + class HessianModelFactory { public: static std::unique_ptr create(const std::string& hessian_model, size_t dimension, size_t maximum_number_nonzeros, bool convexify, const Options& options); }; - - inline std::unique_ptr HessianModelFactory::create(const std::string& hessian_model, size_t dimension, size_t maximum_number_nonzeros, - bool convexify, const Options& options) { - if (hessian_model == "exact") { - if (convexify) { - return std::make_unique(dimension, maximum_number_nonzeros + dimension, options); - } - else { - return std::make_unique(dimension, maximum_number_nonzeros, options); - } - } - else if (hessian_model == "zero") { - return std::make_unique(dimension, options); - } - throw std::invalid_argument("Hessian model " + hessian_model + " does not exist"); - } } // namespace #endif // UNO_HESSIANMODELFACTORY_H, \ No newline at end of file diff --git a/uno/ingredients/subproblems/inequality_constrained_methods/InequalityConstrainedMethod.cpp b/uno/ingredients/subproblems/inequality_constrained_methods/InequalityConstrainedMethod.cpp index d2fc14a2..7091167e 100644 --- a/uno/ingredients/subproblems/inequality_constrained_methods/InequalityConstrainedMethod.cpp +++ b/uno/ingredients/subproblems/inequality_constrained_methods/InequalityConstrainedMethod.cpp @@ -2,7 +2,8 @@ // Licensed under the MIT license. See LICENSE file in the project directory for details. #include "InequalityConstrainedMethod.hpp" -#include "ingredients/subproblems/Direction.hpp" +#include "optimization/Direction.hpp" +#include "optimization/Iterate.hpp" #include "linear_algebra/Vector.hpp" #include "reformulation/l1RelaxedProblem.hpp" #include "options/Options.hpp" diff --git a/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.cpp b/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.cpp index 87648434..7d3c9298 100644 --- a/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.cpp +++ b/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.cpp @@ -2,7 +2,7 @@ // Licensed under the MIT license. See LICENSE file in the project directory for details. #include "LPSubproblem.hpp" -#include "ingredients/subproblems/Direction.hpp" +#include "optimization/Direction.hpp" #include "optimization/WarmstartInformation.hpp" #include "reformulation/OptimizationProblem.hpp" #include "solvers/LPSolver.hpp" @@ -18,6 +18,8 @@ namespace uno { zero_hessian(SymmetricMatrix::zero(number_variables)) { } + LPSubproblem::~LPSubproblem() { } + void LPSubproblem::generate_initial_iterate(const OptimizationProblem& /*problem*/, Iterate& /*initial_iterate*/) { } diff --git a/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.hpp b/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.hpp index 40472df7..b0670c0e 100644 --- a/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.hpp +++ b/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.hpp @@ -16,6 +16,7 @@ namespace uno { public: LPSubproblem(size_t number_variables, size_t number_constraints, size_t number_objective_gradient_nonzeros, size_t number_jacobian_nonzeros, const Options& options); + ~LPSubproblem(); void generate_initial_iterate(const OptimizationProblem& problem, Iterate& initial_iterate) override; void solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, const Multipliers& current_multipliers, diff --git a/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.cpp b/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.cpp index e3f3b1d1..b4e81449 100644 --- a/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.cpp +++ b/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.cpp @@ -2,7 +2,7 @@ // Licensed under the MIT license. See LICENSE file in the project directory for details. #include "QPSubproblem.hpp" -#include "ingredients/subproblems/Direction.hpp" +#include "optimization/Direction.hpp" #include "ingredients/hessian_models/HessianModelFactory.hpp" #include "linear_algebra/SymmetricMatrix.hpp" #include "optimization/Iterate.hpp" @@ -29,6 +29,8 @@ namespace uno { options)) { } + QPSubproblem::~QPSubproblem() { } + void QPSubproblem::initialize_statistics(Statistics& statistics, const Options& options) { if (this->use_regularization) { statistics.add_column("regularization", Statistics::double_width, options.get_int("statistics_regularization_column_order")); diff --git a/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.hpp b/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.hpp index bf1bc2f1..b0eea453 100644 --- a/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.hpp +++ b/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.hpp @@ -15,6 +15,7 @@ namespace uno { public: QPSubproblem(size_t number_variables, size_t number_constraints, size_t number_objective_gradient_nonzeros, size_t number_jacobian_nonzeros, size_t number_hessian_nonzeros, const Options& options); + ~QPSubproblem(); void initialize_statistics(Statistics& statistics, const Options& options) override; void generate_initial_iterate(const OptimizationProblem& problem, Iterate& initial_iterate) override; diff --git a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp index 1eaf841f..abec25bc 100644 --- a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp +++ b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp @@ -3,7 +3,8 @@ #include #include "PrimalDualInteriorPointSubproblem.hpp" -#include "ingredients/subproblems/Direction.hpp" +#include "optimization/Direction.hpp" +#include "optimization/Iterate.hpp" #include "ingredients/hessian_models/HessianModelFactory.hpp" #include "linear_algebra/SparseStorageFactory.hpp" #include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp" @@ -49,12 +50,12 @@ namespace uno { l1_constraint_violation_coefficient(options.get_double("l1_constraint_violation_coefficient")) { } - inline void PrimalDualInteriorPointSubproblem::initialize_statistics(Statistics& statistics, const Options& options) { + void PrimalDualInteriorPointSubproblem::initialize_statistics(Statistics& statistics, const Options& options) { statistics.add_column("regularization", Statistics::double_width - 1, options.get_int("statistics_regularization_column_order")); statistics.add_column("barrier param.", Statistics::double_width - 1, options.get_int("statistics_barrier_parameter_column_order")); } - inline void PrimalDualInteriorPointSubproblem::generate_initial_iterate(const OptimizationProblem& problem, Iterate& initial_iterate) { + void PrimalDualInteriorPointSubproblem::generate_initial_iterate(const OptimizationProblem& problem, Iterate& initial_iterate) { if (problem.has_inequality_constraints()) { throw std::runtime_error("The problem has inequality constraints. Create an instance of HomogeneousEqualityConstrainedModel"); } @@ -497,7 +498,7 @@ namespace uno { } } else { - WARNING << YELLOW << "Barrier subproblem: the bounds are in the wrong order in the lower bound multiplier reset" << RESET << '\n'; + WARNING << "Barrier subproblem: the bounds are in the wrong order in the lower bound multiplier reset\n"; } } @@ -515,7 +516,7 @@ namespace uno { } } else { - WARNING << YELLOW << "Barrier subproblem: the bounds are in the wrong order in the upper bound multiplier reset" << RESET << '\n'; + WARNING << "Barrier subproblem: the bounds are in the wrong order in the upper bound multiplier reset\n"; } } } diff --git a/uno/model/BoundRelaxedModel.cpp b/uno/model/BoundRelaxedModel.cpp new file mode 100644 index 00000000..c64e7324 --- /dev/null +++ b/uno/model/BoundRelaxedModel.cpp @@ -0,0 +1,25 @@ +// Copyright (c) 2018-2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "BoundRelaxedModel.hpp" +#include "optimization/Iterate.hpp" +#include "options/Options.hpp" + +namespace uno { + BoundRelaxedModel::BoundRelaxedModel(std::unique_ptr original_model, const Options& options): + Model(original_model->name + " -> bounds relaxed", original_model->number_variables, original_model->number_constraints, + original_model->objective_sign), + model(std::move(original_model)), + relaxation_factor(options.get_double("tolerance")) { + } + + double BoundRelaxedModel::variable_lower_bound(size_t variable_index) const { + const double lower_bound = this->model->variable_lower_bound(variable_index); + return lower_bound - this->relaxation_factor * std::max(1., std::abs(lower_bound)); + } + + double BoundRelaxedModel::variable_upper_bound(size_t variable_index) const { + const double upper_bound = this->model->variable_upper_bound(variable_index); + return upper_bound + this->relaxation_factor * std::max(1., std::abs(upper_bound)); + } +} // namespace \ No newline at end of file diff --git a/uno/model/BoundRelaxedModel.hpp b/uno/model/BoundRelaxedModel.hpp index 36dbc4f6..0595e95c 100644 --- a/uno/model/BoundRelaxedModel.hpp +++ b/uno/model/BoundRelaxedModel.hpp @@ -4,11 +4,13 @@ #ifndef UNO_BOUNDRELAXEDMODEL_H #define UNO_BOUNDRELAXEDMODEL_H +#include #include "Model.hpp" -#include "optimization/Iterate.hpp" -#include "options/Options.hpp" namespace uno { + // forward declaration + class Options; + class BoundRelaxedModel: public Model { public: BoundRelaxedModel(std::unique_ptr original_model, const Options& options); @@ -64,23 +66,6 @@ namespace uno { const std::unique_ptr model{}; const double relaxation_factor; }; - - inline BoundRelaxedModel::BoundRelaxedModel(std::unique_ptr original_model, const Options& options): - Model(original_model->name + " -> bounds relaxed", original_model->number_variables, original_model->number_constraints, - original_model->objective_sign), - model(std::move(original_model)), - relaxation_factor(options.get_double("tolerance")) { - } - - inline double BoundRelaxedModel::variable_lower_bound(size_t variable_index) const { - const double lower_bound = this->model->variable_lower_bound(variable_index); - return lower_bound - this->relaxation_factor * std::max(1., std::abs(lower_bound)); - } - - inline double BoundRelaxedModel::variable_upper_bound(size_t variable_index) const { - const double upper_bound = this->model->variable_upper_bound(variable_index); - return upper_bound + this->relaxation_factor * std::max(1., std::abs(upper_bound)); - } } // namespace #endif // UNO_BOUNDRELAXEDMODEL_H \ No newline at end of file diff --git a/uno/model/FixedBoundsConstraintsModel.cpp b/uno/model/FixedBoundsConstraintsModel.cpp new file mode 100644 index 00000000..96a6d08a --- /dev/null +++ b/uno/model/FixedBoundsConstraintsModel.cpp @@ -0,0 +1,219 @@ +// Copyright (c) 2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "FixedBoundsConstraintsModel.hpp" +#include "optimization/Iterate.hpp" + +namespace uno { + FixedBoundsConstraintsModel::FixedBoundsConstraintsModel(std::unique_ptr original_model, const Options& /*options*/) : + Model(original_model->name + " -> no fixed bounds", original_model->number_variables, + // move the fixed variables to the set of general constraints + original_model->number_constraints + original_model->get_fixed_variables().size(), + original_model->objective_sign), + model(std::move(original_model)), + fixed_variables(), // empty vector + lower_bounded_variables_collection(this->lower_bounded_variables), + upper_bounded_variables_collection(this->upper_bounded_variables), + equality_constraints(concatenate(this->model->get_equality_constraints(), Range(this->model->number_constraints, this->number_constraints))), + linear_constraints(concatenate(this->model->get_linear_constraints(), Range(this->model->number_constraints, this->number_constraints))) { + this->lower_bounded_variables.reserve(this->model->get_lower_bounded_variables().size()); + this->upper_bounded_variables.reserve(this->model->get_upper_bounded_variables().size()); + + for (size_t variable_index: Range(this->model->number_variables)) { + if (is_finite(this->variable_lower_bound(variable_index))) { + this->lower_bounded_variables.emplace_back(variable_index); + } + if (is_finite(this->variable_upper_bound(variable_index))) { + this->upper_bounded_variables.emplace_back(variable_index); + } + } + } + + double FixedBoundsConstraintsModel::evaluate_objective(const Vector& x) const { + return this->model->evaluate_objective(x); + } + + void FixedBoundsConstraintsModel::evaluate_objective_gradient(const Vector& x, SparseVector& gradient) const { + this->model->evaluate_objective_gradient(x, gradient); + } + + void FixedBoundsConstraintsModel::evaluate_constraints(const Vector& x, std::vector& constraints) const { + this->model->evaluate_constraints(x, constraints); + // add the fixed variables + size_t current_constraint = this->model->number_constraints; + for (size_t fixed_variable_index: this->model->get_fixed_variables()) { + constraints[current_constraint] = x[fixed_variable_index]; + current_constraint++; + } + } + + void FixedBoundsConstraintsModel::evaluate_constraint_gradient(const Vector& x, size_t constraint_index, SparseVector& gradient) const { + if (constraint_index < this->model->number_constraints) { + // original constraint + this->model->evaluate_constraint_gradient(x, constraint_index, gradient); + } + else { + // fixed variable + const size_t variable_index = this->model->get_fixed_variables()[constraint_index - this->model->number_constraints]; + gradient.insert(variable_index, 1.); + } + } + + void FixedBoundsConstraintsModel::evaluate_constraint_jacobian(const Vector& x, RectangularMatrix& constraint_jacobian) const { + this->model->evaluate_constraint_jacobian(x, constraint_jacobian); + // add the fixed variables + size_t current_constraint = this->model->number_constraints; + for (size_t fixed_variable_index: this->model->get_fixed_variables()) { + constraint_jacobian[current_constraint].insert(fixed_variable_index, 1.); + current_constraint++; + } + } + + void FixedBoundsConstraintsModel::evaluate_lagrangian_hessian(const Vector& x, double objective_multiplier, const Vector& multipliers, + SymmetricMatrix& hessian) const { + this->model->evaluate_lagrangian_hessian(x, objective_multiplier, multipliers, hessian); + } + + double FixedBoundsConstraintsModel::variable_lower_bound(size_t variable_index) const { + if (this->model->variable_lower_bound(variable_index) == this->model->variable_upper_bound(variable_index)) { + // remove bounds of fixed variables + return -INF; + } + return this->model->variable_lower_bound(variable_index); + } + + double FixedBoundsConstraintsModel::variable_upper_bound(size_t variable_index) const { + if (this->model->variable_lower_bound(variable_index) == this->model->variable_upper_bound(variable_index)) { +// remove bounds of fixed variables + return INF; + } + return this->model->variable_upper_bound(variable_index); + } + + BoundType FixedBoundsConstraintsModel::get_variable_bound_type(size_t variable_index) const { + if (this->model->variable_lower_bound(variable_index) == this->model->variable_upper_bound(variable_index)) { +// fixed variable: remove the bounds + return BoundType::UNBOUNDED; + } + else { + return this->model->get_variable_bound_type(variable_index); + } + } + + const Collection& FixedBoundsConstraintsModel::get_lower_bounded_variables() const { + return this->lower_bounded_variables_collection; + } + + const Collection& FixedBoundsConstraintsModel::get_upper_bounded_variables() const { + return this->upper_bounded_variables_collection; + } + + const SparseVector& FixedBoundsConstraintsModel::get_slacks() const { + return this->model->get_slacks(); + } + const Collection& FixedBoundsConstraintsModel::get_single_lower_bounded_variables() const { + return this->model->get_single_lower_bounded_variables(); + } + const Collection& FixedBoundsConstraintsModel::get_single_upper_bounded_variables() const { + return this->model->get_single_upper_bounded_variables(); + } + const Vector& FixedBoundsConstraintsModel::get_fixed_variables() const { + return this->fixed_variables; + } + + double FixedBoundsConstraintsModel::constraint_lower_bound(size_t constraint_index) const { + if (constraint_index < this->model->number_constraints) { +// original constraint + return this->model->constraint_lower_bound(constraint_index); + } + else { +// fixed variable + const size_t variable_index = this->model->get_fixed_variables()[constraint_index - this->model->number_constraints]; + return this->model->variable_lower_bound(variable_index); + } + } + + double FixedBoundsConstraintsModel::constraint_upper_bound(size_t constraint_index) const { + if (constraint_index < this->model->number_constraints) { +// original constraint + return this->model->constraint_upper_bound(constraint_index); + } + else { +// fixed variable + const size_t variable_index = this->model->get_fixed_variables()[constraint_index - this->model->number_constraints]; + return this->model->variable_lower_bound(variable_index); + } + } + + FunctionType FixedBoundsConstraintsModel::get_constraint_type(size_t constraint_index) const { + if (constraint_index < this->model->number_constraints) { +// original constraint + return this->model->get_constraint_type(constraint_index); + } + else { +// fixed variables + return FunctionType::LINEAR; + } + } + + BoundType FixedBoundsConstraintsModel::get_constraint_bound_type(size_t constraint_index) const { + if (constraint_index < this->model->number_constraints) { +// original constraint + return this->model->get_constraint_bound_type(constraint_index); + } + else { +// fixed variables + return BoundType::EQUAL_BOUNDS; + } + } + + const Collection& FixedBoundsConstraintsModel::get_equality_constraints() const { + return this->equality_constraints; + } + const Collection& FixedBoundsConstraintsModel::get_inequality_constraints() const { + return this->model->get_inequality_constraints(); + } + const Collection& FixedBoundsConstraintsModel::get_linear_constraints() const { + return this->linear_constraints; + } + + void FixedBoundsConstraintsModel::initial_primal_point(Vector& x) const { + this->model->initial_primal_point(x); +// set the fixed variables + for (size_t variable_index: this->model->get_fixed_variables()) { + x[variable_index] = this->model->variable_lower_bound(variable_index); + } + } + + void FixedBoundsConstraintsModel::initial_dual_point(Vector& multipliers) const { + this->model->initial_dual_point(multipliers); + } + + void FixedBoundsConstraintsModel::postprocess_solution(Iterate& iterate, TerminationStatus termination_status) const { +// move the multipliers back from the general constraints to the bound constraints + size_t current_constraint = this->model->number_constraints; + for (size_t variable_index: this->model->get_fixed_variables()) { + const double constraint_multiplier = iterate.multipliers.constraints[current_constraint]; + if (0. < constraint_multiplier) { + iterate.multipliers.lower_bounds[variable_index] = constraint_multiplier; + } + else { + iterate.multipliers.upper_bounds[variable_index] = constraint_multiplier; + } + current_constraint++; + } + this->model->postprocess_solution(iterate, termination_status); + } + + size_t FixedBoundsConstraintsModel::number_objective_gradient_nonzeros() const { + return this->model->number_objective_gradient_nonzeros(); + } + + size_t FixedBoundsConstraintsModel::number_jacobian_nonzeros() const { + return this->model->number_jacobian_nonzeros() + this->model->get_fixed_variables().size(); + } + + size_t FixedBoundsConstraintsModel::number_hessian_nonzeros() const { + return this->model->number_hessian_nonzeros(); + } +} // namespace \ No newline at end of file diff --git a/uno/model/FixedBoundsConstraintsModel.hpp b/uno/model/FixedBoundsConstraintsModel.hpp index 3a5e7376..f5707008 100644 --- a/uno/model/FixedBoundsConstraintsModel.hpp +++ b/uno/model/FixedBoundsConstraintsModel.hpp @@ -4,14 +4,17 @@ #ifndef UNO_FIXEDBOUNDSCONSTRAINTSMODEL_H #define UNO_FIXEDBOUNDSCONSTRAINTSMODEL_H +#include #include "Model.hpp" -#include "optimization/Iterate.hpp" -#include "options/Options.hpp" +#include "linear_algebra/Vector.hpp" #include "symbolic/CollectionAdapter.hpp" #include "symbolic/Concatenation.hpp" #include "symbolic/Range.hpp" namespace uno { + // forward declaration + class Options; + // move the fixed variables to the set of general constraints (required in barrier methods) // for instance, the constraint x[2] == 1 is not interpreted as 1 <= x[2] <= 1, but as the // linear constraint 0 x[1] + 1 x[2] + ... + 0 x[n] == 1 @@ -19,172 +22,40 @@ namespace uno { public: FixedBoundsConstraintsModel(std::unique_ptr original_model, const Options& options); - [[nodiscard]] double evaluate_objective(const Vector& x) const override { return this->model->evaluate_objective(x); } - - void evaluate_objective_gradient(const Vector& x, SparseVector& gradient) const override { - this->model->evaluate_objective_gradient(x, gradient); - } - - void evaluate_constraints(const Vector& x, std::vector& constraints) const override { - this->model->evaluate_constraints(x, constraints); - // add the fixed variables - size_t current_constraint = this->model->number_constraints; - for (size_t fixed_variable_index: this->model->get_fixed_variables()) { - constraints[current_constraint] = x[fixed_variable_index]; - current_constraint++; - } - } - - void evaluate_constraint_gradient(const Vector& x, size_t constraint_index, SparseVector& gradient) const override { - if (constraint_index < this->model->number_constraints) { - // original constraint - this->model->evaluate_constraint_gradient(x, constraint_index, gradient); - } - else { - // fixed variable - const size_t variable_index = this->model->get_fixed_variables()[constraint_index - this->model->number_constraints]; - gradient.insert(variable_index, 1.); - } - } - - void evaluate_constraint_jacobian(const Vector& x, RectangularMatrix& constraint_jacobian) const override { - this->model->evaluate_constraint_jacobian(x, constraint_jacobian); - // add the fixed variables - size_t current_constraint = this->model->number_constraints; - for (size_t fixed_variable_index: this->model->get_fixed_variables()) { - constraint_jacobian[current_constraint].insert(fixed_variable_index, 1.); - current_constraint++; - } - } - + [[nodiscard]] double evaluate_objective(const Vector& x) const override; + void evaluate_objective_gradient(const Vector& x, SparseVector& gradient) const override; + void evaluate_constraints(const Vector& x, std::vector& constraints) const override; + void evaluate_constraint_gradient(const Vector& x, size_t constraint_index, SparseVector& gradient) const override; + void evaluate_constraint_jacobian(const Vector& x, RectangularMatrix& constraint_jacobian) const override; void evaluate_lagrangian_hessian(const Vector& x, double objective_multiplier, const Vector& multipliers, - SymmetricMatrix& hessian) const override { - this->model->evaluate_lagrangian_hessian(x, objective_multiplier, multipliers, hessian); - } - - // only these two functions are redefined - [[nodiscard]] double variable_lower_bound(size_t variable_index) const override { - if (this->model->variable_lower_bound(variable_index) == this->model->variable_upper_bound(variable_index)) { - // remove bounds of fixed variables - return -INF; - } - return this->model->variable_lower_bound(variable_index); - } - - [[nodiscard]] double variable_upper_bound(size_t variable_index) const override { - if (this->model->variable_lower_bound(variable_index) == this->model->variable_upper_bound(variable_index)) { - // remove bounds of fixed variables - return INF; - } - return this->model->variable_upper_bound(variable_index); - } - - [[nodiscard]] BoundType get_variable_bound_type(size_t variable_index) const override { - if (this->model->variable_lower_bound(variable_index) == this->model->variable_upper_bound(variable_index)) { - // fixed variable: remove the bounds - return BoundType::UNBOUNDED; - } - else { - return this->model->get_variable_bound_type(variable_index); - } - } - - [[nodiscard]] const Collection& get_lower_bounded_variables() const override { - return this->lower_bounded_variables_collection; - } - - [[nodiscard]] const Collection& get_upper_bounded_variables() const override { - return this->upper_bounded_variables_collection; - } - - [[nodiscard]] const SparseVector& get_slacks() const override { return this->model->get_slacks(); } - [[nodiscard]] const Collection& get_single_lower_bounded_variables() const override { return this->model->get_single_lower_bounded_variables(); } - [[nodiscard]] const Collection& get_single_upper_bounded_variables() const override { return this->model->get_single_upper_bounded_variables(); } - [[nodiscard]] const Vector& get_fixed_variables() const override { return this->fixed_variables; } - - [[nodiscard]] double constraint_lower_bound(size_t constraint_index) const override { - if (constraint_index < this->model->number_constraints) { - // original constraint - return this->model->constraint_lower_bound(constraint_index); - } - else { - // fixed variable - const size_t variable_index = this->model->get_fixed_variables()[constraint_index - this->model->number_constraints]; - return this->model->variable_lower_bound(variable_index); - } - } - - [[nodiscard]] double constraint_upper_bound(size_t constraint_index) const override { - if (constraint_index < this->model->number_constraints) { - // original constraint - return this->model->constraint_upper_bound(constraint_index); - } - else { - // fixed variable - const size_t variable_index = this->model->get_fixed_variables()[constraint_index - this->model->number_constraints]; - return this->model->variable_lower_bound(variable_index); - } - } - - [[nodiscard]] FunctionType get_constraint_type(size_t constraint_index) const override { - if (constraint_index < this->model->number_constraints) { - // original constraint - return this->model->get_constraint_type(constraint_index); - } - else { - // fixed variables - return FunctionType::LINEAR; - } - } - - [[nodiscard]] BoundType get_constraint_bound_type(size_t constraint_index) const override { - if (constraint_index < this->model->number_constraints) { - // original constraint - return this->model->get_constraint_bound_type(constraint_index); - } - else { - // fixed variables - return BoundType::EQUAL_BOUNDS; - } - } - - [[nodiscard]] const Collection& get_equality_constraints() const override { return this->equality_constraints; } - [[nodiscard]] const Collection& get_inequality_constraints() const override { return this->model->get_inequality_constraints(); } - [[nodiscard]] const Collection& get_linear_constraints() const override { return this->linear_constraints; } - - void initial_primal_point(Vector& x) const override { - this->model->initial_primal_point(x); - // set the fixed variables - for (size_t variable_index: this->model->get_fixed_variables()) { - x[variable_index] = this->model->variable_lower_bound(variable_index); - } - } - - void initial_dual_point(Vector& multipliers) const override { this->model->initial_dual_point(multipliers); } - - void postprocess_solution(Iterate& iterate, TerminationStatus termination_status) const override { - // move the multipliers back from the general constraints to the bound constraints - size_t current_constraint = this->model->number_constraints; - for (size_t variable_index: this->model->get_fixed_variables()) { - const double constraint_multiplier = iterate.multipliers.constraints[current_constraint]; - if (0. < constraint_multiplier) { - iterate.multipliers.lower_bounds[variable_index] = constraint_multiplier; - } - else { - iterate.multipliers.upper_bounds[variable_index] = constraint_multiplier; - } - current_constraint++; - } - this->model->postprocess_solution(iterate, termination_status); - } - - [[nodiscard]] size_t number_objective_gradient_nonzeros() const override { return this->model->number_objective_gradient_nonzeros(); } - - [[nodiscard]] size_t number_jacobian_nonzeros() const override { - return this->model->number_jacobian_nonzeros() + this->model->get_fixed_variables().size(); - } - - [[nodiscard]] size_t number_hessian_nonzeros() const override { return this->model->number_hessian_nonzeros(); } + SymmetricMatrix& hessian) const override; + + [[nodiscard]] double variable_lower_bound(size_t variable_index) const override; + [[nodiscard]] double variable_upper_bound(size_t variable_index) const override; + [[nodiscard]] BoundType get_variable_bound_type(size_t variable_index) const override; + [[nodiscard]] const Collection& get_lower_bounded_variables() const override; + [[nodiscard]] const Collection& get_upper_bounded_variables() const override; + [[nodiscard]] const SparseVector& get_slacks() const override; + [[nodiscard]] const Collection& get_single_lower_bounded_variables() const override; + [[nodiscard]] const Collection& get_single_upper_bounded_variables() const override; + [[nodiscard]] const Vector& get_fixed_variables() const override; + + [[nodiscard]] double constraint_lower_bound(size_t constraint_index) const override; + [[nodiscard]] double constraint_upper_bound(size_t constraint_index) const override; + [[nodiscard]] FunctionType get_constraint_type(size_t constraint_index) const override; + [[nodiscard]] BoundType get_constraint_bound_type(size_t constraint_index) const override; + [[nodiscard]] const Collection& get_equality_constraints() const override; + [[nodiscard]] const Collection& get_inequality_constraints() const override; + [[nodiscard]] const Collection& get_linear_constraints() const override; + + void initial_primal_point(Vector& x) const override; + void initial_dual_point(Vector& multipliers) const override; + + void postprocess_solution(Iterate& iterate, TerminationStatus termination_status) const override; + + [[nodiscard]] size_t number_objective_gradient_nonzeros() const override; + [[nodiscard]] size_t number_jacobian_nonzeros() const override; + [[nodiscard]] size_t number_hessian_nonzeros() const override; private: const std::unique_ptr model{}; @@ -196,30 +67,6 @@ namespace uno { Concatenation&, ForwardRange> equality_constraints; Concatenation&, ForwardRange> linear_constraints; }; - - inline FixedBoundsConstraintsModel::FixedBoundsConstraintsModel(std::unique_ptr original_model, const Options& /*options*/): - Model(original_model->name + " -> no fixed bounds", original_model->number_variables, - // move the fixed variables to the set of general constraints - original_model->number_constraints + original_model->get_fixed_variables().size(), - original_model->objective_sign), - model(std::move(original_model)), - fixed_variables(), // empty vector - lower_bounded_variables_collection(this->lower_bounded_variables), - upper_bounded_variables_collection(this->upper_bounded_variables), - equality_constraints(concatenate(this->model->get_equality_constraints(), Range(this->model->number_constraints, this->number_constraints))), - linear_constraints(concatenate(this->model->get_linear_constraints(), Range(this->model->number_constraints, this->number_constraints))) { - this->lower_bounded_variables.reserve(this->model->get_lower_bounded_variables().size()); - this->upper_bounded_variables.reserve(this->model->get_upper_bounded_variables().size()); - - for (size_t variable_index: Range(this->model->number_variables)) { - if (is_finite(this->variable_lower_bound(variable_index))) { - this->lower_bounded_variables.emplace_back(variable_index); - } - if (is_finite(this->variable_upper_bound(variable_index))) { - this->upper_bounded_variables.emplace_back(variable_index); - } - } - } } // namespace #endif // UNO_FIXEDBOUNDSCONSTRAINTSMODEL_H \ No newline at end of file diff --git a/uno/model/HomogeneousEqualityConstrainedModel.cpp b/uno/model/HomogeneousEqualityConstrainedModel.cpp new file mode 100644 index 00000000..48d44940 --- /dev/null +++ b/uno/model/HomogeneousEqualityConstrainedModel.cpp @@ -0,0 +1,214 @@ +// Copyright (c) 2018-2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "HomogeneousEqualityConstrainedModel.hpp" +#include "linear_algebra/SymmetricMatrix.hpp" +#include "optimization/Iterate.hpp" +#include "symbolic/Concatenation.hpp" +#include "symbolic/CollectionAdapter.hpp" +#include "symbolic/Range.hpp" +#include "tools/Infinity.hpp" + +namespace uno { + // Transform the problem into an equality-constrained problem with constraints c(x) = 0. This implies: + // - inequality constraints get a slack + // - equality constraints are shifted by their RHS + HomogeneousEqualityConstrainedModel::HomogeneousEqualityConstrainedModel(std::unique_ptr original_model): + Model(original_model->name + " -> equality constrained", original_model->number_variables + + original_model->get_inequality_constraints().size(), original_model->number_constraints, original_model->objective_sign), + // transfer ownership of the pointer + model(std::move(original_model)), + constraint_index_of_inequality_index(this->model->get_inequality_constraints().size()), + slack_index_of_constraint_index(this->model->number_constraints), + // all constraints are equality constraints + equality_constraints(Range(this->number_constraints)), + inequality_constraints(Range(0)), + slacks(this->model->get_inequality_constraints().size()), + lower_bounded_slacks(this->slacks.size()), + upper_bounded_slacks(this->slacks.size()), + lower_bounded_variables(concatenate(this->model->get_lower_bounded_variables(), adapt(this->lower_bounded_slacks))), + upper_bounded_variables(concatenate(this->model->get_upper_bounded_variables(), adapt(this->upper_bounded_slacks))), + single_lower_bounded_variables(concatenate(this->model->get_single_lower_bounded_variables(), adapt(this->single_lower_bounded_slacks))), + single_upper_bounded_variables(concatenate(this->model->get_single_upper_bounded_variables(), adapt(this->single_upper_bounded_slacks))){ + // register the inequality constraint of each slack + size_t inequality_index = 0; + for (const size_t constraint_index: this->model->get_inequality_constraints()) { + const size_t slack_variable_index = this->model->number_variables + inequality_index; + this->constraint_index_of_inequality_index[inequality_index] = constraint_index; + this->slack_index_of_constraint_index[constraint_index] = slack_variable_index; + this->slacks.insert(constraint_index, slack_variable_index); + if (is_finite(this->model->constraint_lower_bound(constraint_index))) { + this->lower_bounded_slacks.emplace_back(slack_variable_index); + if (not is_finite(this->model->constraint_upper_bound(constraint_index))) { + this->single_lower_bounded_slacks.emplace_back(slack_variable_index); + } + } + if (is_finite(this->model->constraint_upper_bound(constraint_index))) { + this->upper_bounded_slacks.emplace_back(slack_variable_index); + if (not is_finite(this->model->constraint_lower_bound(constraint_index))) { + this->single_upper_bounded_slacks.emplace_back(slack_variable_index); + } + } + inequality_index++; + } + } + + double HomogeneousEqualityConstrainedModel::evaluate_objective(const Vector& x) const { + return this->model->evaluate_objective(x); + } + + void HomogeneousEqualityConstrainedModel::evaluate_objective_gradient(const Vector& x, SparseVector& gradient) const { + this->model->evaluate_objective_gradient(x, gradient); + } + + void HomogeneousEqualityConstrainedModel::evaluate_constraints(const Vector& x, std::vector& constraints) const { + this->model->evaluate_constraints(x, constraints); + // inequality constraints: add the slacks + for (const auto [constraint_index, slack_index]: this->get_slacks()) { + constraints[constraint_index] -= x[slack_index]; + } + + // equality constraints: make sure they are homogeneous (c(x) = 0) + for (const size_t constraint_index: this->model->get_equality_constraints()) { + constraints[constraint_index] -= this->model->constraint_lower_bound(constraint_index); + } + } + + void HomogeneousEqualityConstrainedModel::evaluate_constraint_gradient(const Vector& x, size_t constraint_index, + SparseVector& gradient) const { + this->model->evaluate_constraint_gradient(x, constraint_index, gradient); + // if the original constraint is an inequality, add the slack contribution + if (this->model->get_constraint_bound_type(constraint_index) != EQUAL_BOUNDS) { + const size_t slack_variable_index = this->slack_index_of_constraint_index[constraint_index]; + gradient.insert(slack_variable_index, -1.); + } + } + + void HomogeneousEqualityConstrainedModel::evaluate_constraint_jacobian(const Vector& x, RectangularMatrix& constraint_jacobian) const { + this->model->evaluate_constraint_jacobian(x, constraint_jacobian); + // add the slack contributions + for (const auto [constraint_index, slack_index]: this->get_slacks()) { + constraint_jacobian[constraint_index].insert(slack_index, -1.); + } + } + + void HomogeneousEqualityConstrainedModel::evaluate_lagrangian_hessian(const Vector& x, double objective_multiplier, + const Vector& multipliers, SymmetricMatrix& hessian) const { + this->model->evaluate_lagrangian_hessian(x, objective_multiplier, multipliers, hessian); + // extend the dimension of the Hessian by finalizing the remaining columns (note: the slacks do not enter the Hessian) + for (size_t constraint_index: Range(this->model->number_variables, this->number_variables)) { + hessian.finalize_column(constraint_index); + } + } + + double HomogeneousEqualityConstrainedModel::variable_lower_bound(size_t variable_index) const { + if (variable_index < this->model->number_variables) { // original variable + return this->model->variable_lower_bound(variable_index); + } + else { // slack variable + const size_t slack_index = variable_index - this->model->number_variables; + const size_t constraint_index = this->constraint_index_of_inequality_index[slack_index]; + return this->model->constraint_lower_bound(constraint_index); + } + } + + double HomogeneousEqualityConstrainedModel::variable_upper_bound(size_t variable_index) const { + if (variable_index < this->model->number_variables) { // original variable + return this->model->variable_upper_bound(variable_index); + } + else { // slack variable + const size_t inequality_index = variable_index - this->model->number_variables; + const size_t constraint_index = this->constraint_index_of_inequality_index[inequality_index]; + return this->model->constraint_upper_bound(constraint_index); + } + } + + BoundType HomogeneousEqualityConstrainedModel::get_variable_bound_type(size_t variable_index) const { + if (variable_index < this->model->number_variables) { // original variable + return this->model->get_variable_bound_type(variable_index); + } + else { // slack variable + const size_t inequality_index = variable_index - this->model->number_variables; + const size_t constraint_index = this->constraint_index_of_inequality_index[inequality_index]; + return this->model->get_constraint_bound_type(constraint_index); + } + } + + const Collection& HomogeneousEqualityConstrainedModel::get_lower_bounded_variables() const { + return this->lower_bounded_variables; + } + + const Collection& HomogeneousEqualityConstrainedModel::get_upper_bounded_variables() const { + return this->upper_bounded_variables; + } + + const SparseVector& HomogeneousEqualityConstrainedModel::get_slacks() const { + return this->slacks; + } + + const Collection& HomogeneousEqualityConstrainedModel::get_single_lower_bounded_variables() const { + return this->single_lower_bounded_variables; + } + + const Collection& HomogeneousEqualityConstrainedModel::get_single_upper_bounded_variables() const { + return this->single_upper_bounded_variables; + } + + double HomogeneousEqualityConstrainedModel::constraint_lower_bound(size_t /*constraint_index*/) const { + return 0.; } // c(x) = 0 + + double HomogeneousEqualityConstrainedModel::constraint_upper_bound(size_t /*constraint_index*/) const { + return 0.; } + + FunctionType HomogeneousEqualityConstrainedModel::get_constraint_type(size_t constraint_index) const { + return this->model->get_constraint_type(constraint_index); } + + BoundType HomogeneousEqualityConstrainedModel::get_constraint_bound_type(size_t /*constraint_index*/) const { + return EQUAL_BOUNDS; + } + + const Collection& HomogeneousEqualityConstrainedModel::get_equality_constraints() const { + return this->equality_constraints; + } + + const Collection& HomogeneousEqualityConstrainedModel::get_inequality_constraints() const { + return this->inequality_constraints; + } + + const Collection& HomogeneousEqualityConstrainedModel::get_linear_constraints() const { + return this->model->get_linear_constraints(); + } + + const Vector& HomogeneousEqualityConstrainedModel::get_fixed_variables() const { + return this->model->get_fixed_variables(); + } + + void HomogeneousEqualityConstrainedModel::initial_primal_point(Vector& x) const { + this->model->initial_primal_point(x); + // set the slacks + for (const auto [_, slack_index]: this->get_slacks()) { + x[slack_index] = 0.; + } + } + + void HomogeneousEqualityConstrainedModel::initial_dual_point(Vector& multipliers) const { + this->model->initial_dual_point(multipliers); + } + + void HomogeneousEqualityConstrainedModel::postprocess_solution(Iterate& iterate, TerminationStatus termination_status) const { + // discard the slacks + iterate.number_variables = this->model->number_variables; + this->model->postprocess_solution(iterate, termination_status); + } + + size_t HomogeneousEqualityConstrainedModel::number_objective_gradient_nonzeros() const { + return this->model->number_objective_gradient_nonzeros(); } + + size_t HomogeneousEqualityConstrainedModel::number_jacobian_nonzeros() const { + return this->model->number_jacobian_nonzeros() + this->slacks.size(); + } + + size_t HomogeneousEqualityConstrainedModel::number_hessian_nonzeros() const { + return this->model->number_hessian_nonzeros(); + } +} // namespace \ No newline at end of file diff --git a/uno/model/HomogeneousEqualityConstrainedModel.hpp b/uno/model/HomogeneousEqualityConstrainedModel.hpp index 7b056459..1250c9db 100644 --- a/uno/model/HomogeneousEqualityConstrainedModel.hpp +++ b/uno/model/HomogeneousEqualityConstrainedModel.hpp @@ -4,13 +4,11 @@ #ifndef UNO_HOMOGENEOUSEQUALITYCONSTRAINEDMODEL_H #define UNO_HOMOGENEOUSEQUALITYCONSTRAINEDMODEL_H +#include #include "Model.hpp" -#include "linear_algebra/SymmetricMatrix.hpp" -#include "optimization/Iterate.hpp" -#include "symbolic/Concatenation.hpp" +#include "linear_algebra/SparseVector.hpp" #include "symbolic/CollectionAdapter.hpp" -#include "symbolic/Range.hpp" -#include "tools/Infinity.hpp" +#include "symbolic/Concatenation.hpp" namespace uno { // generate an equality-constrained model by: @@ -21,8 +19,8 @@ namespace uno { public: explicit HomogeneousEqualityConstrainedModel(std::unique_ptr original_model); - [[nodiscard]] double evaluate_objective(const Vector& x) const override { return this->model->evaluate_objective(x); } - void evaluate_objective_gradient(const Vector& x, SparseVector& gradient) const override { this->model->evaluate_objective_gradient(x, gradient); } + [[nodiscard]] double evaluate_objective(const Vector& x) const override; + void evaluate_objective_gradient(const Vector& x, SparseVector& gradient) const override; void evaluate_constraints(const Vector& x, std::vector& constraints) const override; void evaluate_constraint_gradient(const Vector& x, size_t constraint_index, SparseVector& gradient) const override; void evaluate_constraint_jacobian(const Vector& x, RectangularMatrix& constraint_jacobian) const override; @@ -32,28 +30,28 @@ namespace uno { [[nodiscard]] double variable_lower_bound(size_t variable_index) const override; [[nodiscard]] double variable_upper_bound(size_t variable_index) const override; [[nodiscard]] BoundType get_variable_bound_type(size_t variable_index) const override; - [[nodiscard]] const Collection& get_lower_bounded_variables() const override { return this->lower_bounded_variables; } - [[nodiscard]] const Collection& get_upper_bounded_variables() const override { return this->upper_bounded_variables; } - [[nodiscard]] const SparseVector& get_slacks() const override { return this->slacks; } - [[nodiscard]] const Collection& get_single_lower_bounded_variables() const override { return this->single_lower_bounded_variables; } - [[nodiscard]] const Collection& get_single_upper_bounded_variables() const override { return this->single_upper_bounded_variables; } - - [[nodiscard]] double constraint_lower_bound(size_t /*constraint_index*/) const override { return 0.; } // c(x) = 0 - [[nodiscard]] double constraint_upper_bound(size_t /*constraint_index*/) const override { return 0.; } - [[nodiscard]] FunctionType get_constraint_type(size_t constraint_index) const override { return this->model->get_constraint_type(constraint_index); } - [[nodiscard]] BoundType get_constraint_bound_type(size_t /*constraint_index*/) const override { return EQUAL_BOUNDS; } - [[nodiscard]] const Collection& get_equality_constraints() const override { return this->equality_constraints; } - [[nodiscard]] const Collection& get_inequality_constraints() const override { return this->inequality_constraints; } - [[nodiscard]] const Collection& get_linear_constraints() const override { return this->model->get_linear_constraints(); } - [[nodiscard]] const Vector& get_fixed_variables() const override { return this->model->get_fixed_variables(); } + [[nodiscard]] const Collection& get_lower_bounded_variables() const override; + [[nodiscard]] const Collection& get_upper_bounded_variables() const override; + [[nodiscard]] const SparseVector& get_slacks() const override; + [[nodiscard]] const Collection& get_single_lower_bounded_variables() const override; + [[nodiscard]] const Collection& get_single_upper_bounded_variables() const override; + + [[nodiscard]] double constraint_lower_bound(size_t /*constraint_index*/) const override; + [[nodiscard]] double constraint_upper_bound(size_t /*constraint_index*/) const override; + [[nodiscard]] FunctionType get_constraint_type(size_t constraint_index) const override; + [[nodiscard]] BoundType get_constraint_bound_type(size_t /*constraint_index*/) const override; + [[nodiscard]] const Collection& get_equality_constraints() const override; + [[nodiscard]] const Collection& get_inequality_constraints() const override; + [[nodiscard]] const Collection& get_linear_constraints() const override; + [[nodiscard]] const Vector& get_fixed_variables() const override; void initial_primal_point(Vector& x) const override; - void initial_dual_point(Vector& multipliers) const override { this->model->initial_dual_point(multipliers); } + void initial_dual_point(Vector& multipliers) const override; void postprocess_solution(Iterate& iterate, TerminationStatus termination_status) const override; - [[nodiscard]] size_t number_objective_gradient_nonzeros() const override { return this->model->number_objective_gradient_nonzeros(); } - [[nodiscard]] size_t number_jacobian_nonzeros() const override { return this->model->number_jacobian_nonzeros() + this->slacks.size(); } - [[nodiscard]] size_t number_hessian_nonzeros() const override { return this->model->number_hessian_nonzeros(); } + [[nodiscard]] size_t number_objective_gradient_nonzeros() const override; + [[nodiscard]] size_t number_jacobian_nonzeros() const override; + [[nodiscard]] size_t number_hessian_nonzeros() const override; protected: const std::unique_ptr model{}; @@ -72,136 +70,6 @@ namespace uno { Concatenation&, CollectionAdapter>> single_lower_bounded_variables; Concatenation&, CollectionAdapter>> single_upper_bounded_variables; }; - - // Transform the problem into an equality-constrained problem with constraints c(x) = 0. This implies: - // - inequality constraints get a slack - // - equality constraints are shifted by their RHS - inline HomogeneousEqualityConstrainedModel::HomogeneousEqualityConstrainedModel(std::unique_ptr original_model): - Model(original_model->name + " -> equality constrained", original_model->number_variables + - original_model->get_inequality_constraints().size(), original_model->number_constraints, original_model->objective_sign), - // transfer ownership of the pointer - model(std::move(original_model)), - constraint_index_of_inequality_index(this->model->get_inequality_constraints().size()), - slack_index_of_constraint_index(this->model->number_constraints), - // all constraints are equality constraints - equality_constraints(Range(this->number_constraints)), - inequality_constraints(Range(0)), - slacks(this->model->get_inequality_constraints().size()), - lower_bounded_slacks(this->slacks.size()), - upper_bounded_slacks(this->slacks.size()), - lower_bounded_variables(concatenate(this->model->get_lower_bounded_variables(), adapt(this->lower_bounded_slacks))), - upper_bounded_variables(concatenate(this->model->get_upper_bounded_variables(), adapt(this->upper_bounded_slacks))), - single_lower_bounded_variables(concatenate(this->model->get_single_lower_bounded_variables(), adapt(this->single_lower_bounded_slacks))), - single_upper_bounded_variables(concatenate(this->model->get_single_upper_bounded_variables(), adapt(this->single_upper_bounded_slacks))){ - // register the inequality constraint of each slack - size_t inequality_index = 0; - for (const size_t constraint_index: this->model->get_inequality_constraints()) { - const size_t slack_variable_index = this->model->number_variables + inequality_index; - this->constraint_index_of_inequality_index[inequality_index] = constraint_index; - this->slack_index_of_constraint_index[constraint_index] = slack_variable_index; - this->slacks.insert(constraint_index, slack_variable_index); - if (is_finite(this->model->constraint_lower_bound(constraint_index))) { - this->lower_bounded_slacks.emplace_back(slack_variable_index); - if (not is_finite(this->model->constraint_upper_bound(constraint_index))) { - this->single_lower_bounded_slacks.emplace_back(slack_variable_index); - } - } - if (is_finite(this->model->constraint_upper_bound(constraint_index))) { - this->upper_bounded_slacks.emplace_back(slack_variable_index); - if (not is_finite(this->model->constraint_lower_bound(constraint_index))) { - this->single_upper_bounded_slacks.emplace_back(slack_variable_index); - } - } - inequality_index++; - } - } - - inline void HomogeneousEqualityConstrainedModel::evaluate_constraints(const Vector& x, std::vector& constraints) const { - this->model->evaluate_constraints(x, constraints); - // inequality constraints: add the slacks - for (const auto [constraint_index, slack_index]: this->get_slacks()) { - constraints[constraint_index] -= x[slack_index]; - } - - // equality constraints: make sure they are homogeneous (c(x) = 0) - for (const size_t constraint_index: this->model->get_equality_constraints()) { - constraints[constraint_index] -= this->model->constraint_lower_bound(constraint_index); - } - } - - inline void HomogeneousEqualityConstrainedModel::evaluate_constraint_gradient(const Vector& x, size_t constraint_index, - SparseVector& gradient) const { - this->model->evaluate_constraint_gradient(x, constraint_index, gradient); - // if the original constraint is an inequality, add the slack contribution - if (this->model->get_constraint_bound_type(constraint_index) != EQUAL_BOUNDS) { - const size_t slack_variable_index = this->slack_index_of_constraint_index[constraint_index]; - gradient.insert(slack_variable_index, -1.); - } - } - - inline void HomogeneousEqualityConstrainedModel::evaluate_constraint_jacobian(const Vector& x, RectangularMatrix& constraint_jacobian) const { - this->model->evaluate_constraint_jacobian(x, constraint_jacobian); - // add the slack contributions - for (const auto [constraint_index, slack_index]: this->get_slacks()) { - constraint_jacobian[constraint_index].insert(slack_index, -1.); - } - } - - inline void HomogeneousEqualityConstrainedModel::evaluate_lagrangian_hessian(const Vector& x, double objective_multiplier, - const Vector& multipliers, SymmetricMatrix& hessian) const { - this->model->evaluate_lagrangian_hessian(x, objective_multiplier, multipliers, hessian); - // extend the dimension of the Hessian by finalizing the remaining columns (note: the slacks do not enter the Hessian) - for (size_t constraint_index: Range(this->model->number_variables, this->number_variables)) { - hessian.finalize_column(constraint_index); - } - } - - inline double HomogeneousEqualityConstrainedModel::variable_lower_bound(size_t variable_index) const { - if (variable_index < this->model->number_variables) { // original variable - return this->model->variable_lower_bound(variable_index); - } - else { // slack variable - const size_t slack_index = variable_index - this->model->number_variables; - const size_t constraint_index = this->constraint_index_of_inequality_index[slack_index]; - return this->model->constraint_lower_bound(constraint_index); - } - } - - inline double HomogeneousEqualityConstrainedModel::variable_upper_bound(size_t variable_index) const { - if (variable_index < this->model->number_variables) { // original variable - return this->model->variable_upper_bound(variable_index); - } - else { // slack variable - const size_t inequality_index = variable_index - this->model->number_variables; - const size_t constraint_index = this->constraint_index_of_inequality_index[inequality_index]; - return this->model->constraint_upper_bound(constraint_index); - } - } - - inline BoundType HomogeneousEqualityConstrainedModel::get_variable_bound_type(size_t variable_index) const { - if (variable_index < this->model->number_variables) { // original variable - return this->model->get_variable_bound_type(variable_index); - } - else { // slack variable - const size_t inequality_index = variable_index - this->model->number_variables; - const size_t constraint_index = this->constraint_index_of_inequality_index[inequality_index]; - return this->model->get_constraint_bound_type(constraint_index); - } - } - - inline void HomogeneousEqualityConstrainedModel::initial_primal_point(Vector& x) const { - this->model->initial_primal_point(x); - // set the slacks - for (const auto [_, slack_index]: this->get_slacks()) { - x[slack_index] = 0.; - } - } - - inline void HomogeneousEqualityConstrainedModel::postprocess_solution(Iterate& iterate, TerminationStatus termination_status) const { - // discard the slacks - iterate.number_variables = this->model->number_variables; - this->model->postprocess_solution(iterate, termination_status); - } } // namespace #endif // UNO_HOMOGENEOUSEQUALITYCONSTRAINEDMODEL_H diff --git a/uno/model/ScaledModel.cpp b/uno/model/ScaledModel.cpp new file mode 100644 index 00000000..9e6c38ed --- /dev/null +++ b/uno/model/ScaledModel.cpp @@ -0,0 +1,180 @@ +// Copyright (c) 2018-2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "ScaledModel.hpp" +#include "Model.hpp" +#include "optimization/Iterate.hpp" +#include "options/Options.hpp" + +namespace uno { + ScaledModel::ScaledModel(std::unique_ptr original_model, Iterate& initial_iterate, const Options& options): + Model(original_model->name + " -> scaled", original_model->number_variables, original_model->number_constraints, + original_model->objective_sign), + model(std::move(original_model)), + scaling(this->model->number_constraints, options.get_double("function_scaling_threshold")) { + if (options.get_bool("scale_functions")) { + // evaluate the gradients at the current point + initial_iterate.evaluate_objective_gradient(*this->model); + initial_iterate.evaluate_constraint_jacobian(*this->model); + this->scaling.compute(initial_iterate.evaluations.objective_gradient, initial_iterate.evaluations.constraint_jacobian); + // scale the gradients + scale(initial_iterate.evaluations.objective_gradient, this->scaling.get_objective_scaling()); + for (size_t constraint_index: Range(this->model->number_constraints)) { + scale(initial_iterate.evaluations.constraint_jacobian[constraint_index], this->scaling.get_constraint_scaling(constraint_index)); + } + // since the definition of the constraints changed, reset the evaluation flags + initial_iterate.is_objective_gradient_computed = false; + initial_iterate.is_constraint_jacobian_computed = false; + } + // check the scaling factors + assert(0 < this->scaling.get_objective_scaling() && "Objective scaling failed."); + for ([[maybe_unused]] size_t constraint_index: Range(this->number_constraints)) { + assert(0 < this->scaling.get_constraint_scaling(constraint_index) && "Constraint scaling failed."); + } + } + + double ScaledModel::evaluate_objective(const Vector& x) const { + const double objective = this->model->evaluate_objective(x); + return this->scaling.get_objective_scaling()*objective; + } + + void ScaledModel::evaluate_objective_gradient(const Vector& x, SparseVector& gradient) const { + this->model->evaluate_objective_gradient(x, gradient); + scale(gradient, this->scaling.get_objective_scaling()); + } + + void ScaledModel::evaluate_constraints(const Vector& x, std::vector& constraints) const { + this->model->evaluate_constraints(x, constraints); + for (size_t constraint_index: Range(this->number_constraints)) { + constraints[constraint_index] *= this->scaling.get_constraint_scaling(constraint_index); + } + } + + void ScaledModel::evaluate_constraint_gradient(const Vector& x, size_t constraint_index, SparseVector& gradient) const { + this->model->evaluate_constraint_gradient(x, constraint_index, gradient); + scale(gradient, this->scaling.get_constraint_scaling(constraint_index)); + } + + void ScaledModel::evaluate_constraint_jacobian(const Vector& x, RectangularMatrix& constraint_jacobian) const { + this->model->evaluate_constraint_jacobian(x, constraint_jacobian); + for (size_t constraint_index: Range(this->number_constraints)) { + scale(constraint_jacobian[constraint_index], this->scaling.get_constraint_scaling(constraint_index)); + } + } + + void ScaledModel::evaluate_lagrangian_hessian(const Vector& x, double objective_multiplier, const Vector& multipliers, + SymmetricMatrix& hessian) const { + // scale the objective and constraint multipliers + const double scaled_objective_multiplier = objective_multiplier*this->scaling.get_objective_scaling(); + // TODO preallocate this vector + static Vector scaled_multipliers(this->number_constraints); + for (size_t constraint_index: Range(this->number_constraints)) { + scaled_multipliers[constraint_index] = this->scaling.get_constraint_scaling(constraint_index) * multipliers[constraint_index]; + } + this->model->evaluate_lagrangian_hessian(x, scaled_objective_multiplier, scaled_multipliers, hessian); + } + + double ScaledModel::variable_lower_bound(size_t variable_index) const { + return this->model->variable_lower_bound(variable_index); + } + + double ScaledModel::variable_upper_bound(size_t variable_index) const { + return this->model->variable_upper_bound(variable_index); + } + + BoundType ScaledModel::get_variable_bound_type(size_t variable_index) const { + return this->model->get_variable_bound_type(variable_index); + } + + const Collection& ScaledModel::get_lower_bounded_variables() const { + return this->model->get_lower_bounded_variables(); + } + + const Collection& ScaledModel::get_upper_bounded_variables() const { + return this->model->get_upper_bounded_variables(); + } + + const SparseVector& ScaledModel::get_slacks() const { + return this->model->get_slacks(); + } + + const Collection& ScaledModel::get_single_lower_bounded_variables() const { + return this->model->get_single_lower_bounded_variables(); + } + + const Collection& ScaledModel::get_single_upper_bounded_variables() const { + return this->model->get_single_upper_bounded_variables(); + } + + const Vector& ScaledModel::get_fixed_variables() const { + return this->model->get_fixed_variables(); + } + + double ScaledModel::constraint_lower_bound(size_t constraint_index) const { + return this->scaling.get_constraint_scaling(constraint_index) * this->model->constraint_lower_bound(constraint_index); + } + + double ScaledModel::constraint_upper_bound(size_t constraint_index) const { + return this->scaling.get_constraint_scaling(constraint_index) * this->model->constraint_upper_bound(constraint_index); + } + + FunctionType ScaledModel::get_constraint_type(size_t constraint_index) const { + return this->model->get_constraint_type(constraint_index); + } + + BoundType ScaledModel::get_constraint_bound_type(size_t constraint_index) const { + return this->model->get_constraint_bound_type(constraint_index); + } + + const Collection& ScaledModel::get_equality_constraints() const { + return this->model->get_equality_constraints(); + } + + const Collection& ScaledModel::get_inequality_constraints() const { + return this->model->get_inequality_constraints(); + } + + const Collection& ScaledModel::get_linear_constraints() const { + return this->model->get_linear_constraints(); + } + + void ScaledModel::initial_primal_point(Vector& x) const { + this->model->initial_primal_point(x); + } + + void ScaledModel::initial_dual_point(Vector& multipliers) const { + this->model->initial_dual_point(multipliers); + } + + void ScaledModel::postprocess_solution(Iterate& iterate, TerminationStatus termination_status) const { + // unscale the objective value + if (iterate.is_objective_computed) { + iterate.evaluations.objective /= this->scaling.get_objective_scaling(); + } + + // unscale the constraint multipliers + for (size_t constraint_index: Range(iterate.number_constraints)) { + iterate.multipliers.constraints[constraint_index] *= this->scaling.get_constraint_scaling(constraint_index) / this->scaling.get_objective_scaling(); + } + + // unscale the bound multipliers + for (size_t variable_index: Range(iterate.number_variables)) { + iterate.multipliers.lower_bounds[variable_index] /= this->scaling.get_objective_scaling(); + iterate.multipliers.upper_bounds[variable_index] /= this->scaling.get_objective_scaling(); + } + this->model->postprocess_solution(iterate, termination_status); + } + + size_t ScaledModel::number_objective_gradient_nonzeros() const { + return this->model->number_objective_gradient_nonzeros(); + } + + size_t ScaledModel::number_jacobian_nonzeros() const { + return this->model->number_jacobian_nonzeros(); + } + + size_t ScaledModel::number_hessian_nonzeros() const { + return this->model->number_hessian_nonzeros(); + } +} // namespace + diff --git a/uno/model/ScaledModel.hpp b/uno/model/ScaledModel.hpp index 56ff390a..c3a8199e 100644 --- a/uno/model/ScaledModel.hpp +++ b/uno/model/ScaledModel.hpp @@ -4,12 +4,14 @@ #ifndef UNO_SCALEDMODEL_H #define UNO_SCALEDMODEL_H +#include #include "Model.hpp" -#include "optimization/Iterate.hpp" #include "preprocessing/Scaling.hpp" -#include "options/Options.hpp" namespace uno { + // forward declaration + class Options; + class ScaledModel: public Model { public: ScaledModel(std::unique_ptr original_model, Iterate& initial_iterate, const Options& options); @@ -22,130 +24,36 @@ namespace uno { void evaluate_lagrangian_hessian(const Vector& x, double objective_multiplier, const Vector& multipliers, SymmetricMatrix& hessian) const override; - [[nodiscard]] double variable_lower_bound(size_t variable_index) const override { return this->model->variable_lower_bound(variable_index); } - [[nodiscard]] double variable_upper_bound(size_t variable_index) const override { return this->model->variable_upper_bound(variable_index); } - [[nodiscard]] BoundType get_variable_bound_type(size_t variable_index) const override { return this->model->get_variable_bound_type(variable_index); } - [[nodiscard]] const Collection& get_lower_bounded_variables() const override { return this->model->get_lower_bounded_variables(); } - [[nodiscard]] const Collection& get_upper_bounded_variables() const override { return this->model->get_upper_bounded_variables(); } - [[nodiscard]] const SparseVector& get_slacks() const override { return this->model->get_slacks(); } - [[nodiscard]] const Collection& get_single_lower_bounded_variables() const override { return this->model->get_single_lower_bounded_variables(); } - [[nodiscard]] const Collection& get_single_upper_bounded_variables() const override { return this->model->get_single_upper_bounded_variables(); } - [[nodiscard]] const Vector& get_fixed_variables() const override { return this->model->get_fixed_variables(); } + [[nodiscard]] double variable_lower_bound(size_t variable_index) const override; + [[nodiscard]] double variable_upper_bound(size_t variable_index) const override; + [[nodiscard]] BoundType get_variable_bound_type(size_t variable_index) const override; + [[nodiscard]] const Collection& get_lower_bounded_variables() const override; + [[nodiscard]] const Collection& get_upper_bounded_variables() const override; + [[nodiscard]] const SparseVector& get_slacks() const override; + [[nodiscard]] const Collection& get_single_lower_bounded_variables() const override; + [[nodiscard]] const Collection& get_single_upper_bounded_variables() const override; + [[nodiscard]] const Vector& get_fixed_variables() const override; [[nodiscard]] double constraint_lower_bound(size_t constraint_index) const override; [[nodiscard]] double constraint_upper_bound(size_t constraint_index) const override; - [[nodiscard]] FunctionType get_constraint_type(size_t constraint_index) const override { return this->model->get_constraint_type(constraint_index); } - [[nodiscard]] BoundType get_constraint_bound_type(size_t constraint_index) const override { return this->model->get_constraint_bound_type(constraint_index); } - [[nodiscard]] const Collection& get_equality_constraints() const override { return this->model->get_equality_constraints(); } - [[nodiscard]] const Collection& get_inequality_constraints() const override { return this->model->get_inequality_constraints(); } - [[nodiscard]] const Collection& get_linear_constraints() const override { return this->model->get_linear_constraints(); } - - void initial_primal_point(Vector& x) const override { this->model->initial_primal_point(x); } - void initial_dual_point(Vector& multipliers) const override { this->model->initial_dual_point(multipliers); } + [[nodiscard]] FunctionType get_constraint_type(size_t constraint_index) const override; + [[nodiscard]] BoundType get_constraint_bound_type(size_t constraint_index) const override; + [[nodiscard]] const Collection& get_equality_constraints() const override; + [[nodiscard]] const Collection& get_inequality_constraints() const override; + [[nodiscard]] const Collection& get_linear_constraints() const override; + + void initial_primal_point(Vector& x) const override; + void initial_dual_point(Vector& multipliers) const override; void postprocess_solution(Iterate& iterate, TerminationStatus termination_status) const override; - [[nodiscard]] size_t number_objective_gradient_nonzeros() const override { return this->model->number_objective_gradient_nonzeros(); } - [[nodiscard]] size_t number_jacobian_nonzeros() const override { return this->model->number_jacobian_nonzeros(); } - [[nodiscard]] size_t number_hessian_nonzeros() const override { return this->model->number_hessian_nonzeros(); } + [[nodiscard]] size_t number_objective_gradient_nonzeros() const override; + [[nodiscard]] size_t number_jacobian_nonzeros() const override; + [[nodiscard]] size_t number_hessian_nonzeros() const override; private: const std::unique_ptr model{}; Scaling scaling; }; - - inline ScaledModel::ScaledModel(std::unique_ptr original_model, Iterate& initial_iterate, const Options& options): - Model(original_model->name + " -> scaled", original_model->number_variables, original_model->number_constraints, - original_model->objective_sign), - model(std::move(original_model)), - scaling(this->model->number_constraints, options.get_double("function_scaling_threshold")) { - if (options.get_bool("scale_functions")) { - // evaluate the gradients at the current point - initial_iterate.evaluate_objective_gradient(*this->model); - initial_iterate.evaluate_constraint_jacobian(*this->model); - this->scaling.compute(initial_iterate.evaluations.objective_gradient, initial_iterate.evaluations.constraint_jacobian); - // scale the gradients - scale(initial_iterate.evaluations.objective_gradient, this->scaling.get_objective_scaling()); - for (size_t constraint_index: Range(this->model->number_constraints)) { - scale(initial_iterate.evaluations.constraint_jacobian[constraint_index], this->scaling.get_constraint_scaling(constraint_index)); - } - // since the definition of the constraints changed, reset the evaluation flags - initial_iterate.is_objective_gradient_computed = false; - initial_iterate.is_constraint_jacobian_computed = false; - } - // check the scaling factors - assert(0 < this->scaling.get_objective_scaling() && "Objective scaling failed."); - for ([[maybe_unused]] size_t constraint_index: Range(this->number_constraints)) { - assert(0 < this->scaling.get_constraint_scaling(constraint_index) && "Constraint scaling failed."); - } - } - - inline double ScaledModel::evaluate_objective(const Vector& x) const { - const double objective = this->model->evaluate_objective(x); - return this->scaling.get_objective_scaling()*objective; - } - - inline void ScaledModel::evaluate_objective_gradient(const Vector& x, SparseVector& gradient) const { - this->model->evaluate_objective_gradient(x, gradient); - scale(gradient, this->scaling.get_objective_scaling()); - } - - inline void ScaledModel::evaluate_constraints(const Vector& x, std::vector& constraints) const { - this->model->evaluate_constraints(x, constraints); - for (size_t constraint_index: Range(this->number_constraints)) { - constraints[constraint_index] *= this->scaling.get_constraint_scaling(constraint_index); - } - } - - inline void ScaledModel::evaluate_constraint_gradient(const Vector& x, size_t constraint_index, SparseVector& gradient) const { - this->model->evaluate_constraint_gradient(x, constraint_index, gradient); - scale(gradient, this->scaling.get_constraint_scaling(constraint_index)); - } - - inline void ScaledModel::evaluate_constraint_jacobian(const Vector& x, RectangularMatrix& constraint_jacobian) const { - this->model->evaluate_constraint_jacobian(x, constraint_jacobian); - for (size_t constraint_index: Range(this->number_constraints)) { - scale(constraint_jacobian[constraint_index], this->scaling.get_constraint_scaling(constraint_index)); - } - } - - inline void ScaledModel::evaluate_lagrangian_hessian(const Vector& x, double objective_multiplier, const Vector& multipliers, - SymmetricMatrix& hessian) const { - // scale the objective and constraint multipliers - const double scaled_objective_multiplier = objective_multiplier*this->scaling.get_objective_scaling(); - // TODO preallocate this vector - static Vector scaled_multipliers(this->number_constraints); - for (size_t constraint_index: Range(this->number_constraints)) { - scaled_multipliers[constraint_index] = this->scaling.get_constraint_scaling(constraint_index) * multipliers[constraint_index]; - } - this->model->evaluate_lagrangian_hessian(x, scaled_objective_multiplier, scaled_multipliers, hessian); - } - - inline double ScaledModel::constraint_lower_bound(size_t constraint_index) const { - return this->scaling.get_constraint_scaling(constraint_index) * this->model->constraint_lower_bound(constraint_index); - } - - inline double ScaledModel::constraint_upper_bound(size_t constraint_index) const { - return this->scaling.get_constraint_scaling(constraint_index) * this->model->constraint_upper_bound(constraint_index); - } - - inline void ScaledModel::postprocess_solution(Iterate& iterate, TerminationStatus termination_status) const { - // unscale the objective value - if (iterate.is_objective_computed) { - iterate.evaluations.objective /= this->scaling.get_objective_scaling(); - } - - // unscale the constraint multipliers - for (size_t constraint_index: Range(iterate.number_constraints)) { - iterate.multipliers.constraints[constraint_index] *= this->scaling.get_constraint_scaling(constraint_index) / this->scaling.get_objective_scaling(); - } - - // unscale the bound multipliers - for (size_t variable_index: Range(iterate.number_variables)) { - iterate.multipliers.lower_bounds[variable_index] /= this->scaling.get_objective_scaling(); - iterate.multipliers.upper_bounds[variable_index] /= this->scaling.get_objective_scaling(); - } - this->model->postprocess_solution(iterate, termination_status); - } } // namespace #endif // UNO_SCALEDMODEL_H \ No newline at end of file diff --git a/uno/ingredients/subproblems/Direction.cpp b/uno/optimization/Direction.cpp similarity index 100% rename from uno/ingredients/subproblems/Direction.cpp rename to uno/optimization/Direction.cpp diff --git a/uno/ingredients/subproblems/Direction.hpp b/uno/optimization/Direction.hpp similarity index 86% rename from uno/ingredients/subproblems/Direction.hpp rename to uno/optimization/Direction.hpp index afe2d915..c00d36d4 100644 --- a/uno/ingredients/subproblems/Direction.hpp +++ b/uno/optimization/Direction.hpp @@ -7,17 +7,12 @@ #include #include #include -#include "SubproblemStatus.hpp" +#include "ingredients/subproblems/SubproblemStatus.hpp" #include "linear_algebra/Vector.hpp" #include "optimization/Multipliers.hpp" #include "tools/Infinity.hpp" namespace uno { - /*! \struct ConstraintActivity - * \brief Constraints at lower or upper bound at the optimum solution - * - * Description of the active or infeasible constraints: at lower or upper bound at the optimum solution - */ struct ActiveConstraints { std::vector at_lower_bound; /*!< List of constraint indices at their lower bound */ std::vector at_upper_bound; /*!< List of constraint indices at their upper bound */ diff --git a/uno/optimization/Multipliers.cpp b/uno/optimization/Multipliers.cpp new file mode 100644 index 00000000..118768b7 --- /dev/null +++ b/uno/optimization/Multipliers.cpp @@ -0,0 +1,32 @@ +// Copyright (c) 2018-2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "Multipliers.hpp" + +namespace uno { + Multipliers::Multipliers(size_t number_variables, size_t number_constraints) : lower_bounds(number_variables), + upper_bounds(number_variables), constraints(number_constraints) { + } + + void Multipliers::reset() { + this->constraints.fill(0.); + this->lower_bounds.fill(0.); + this->upper_bounds.fill(0.); + } + + bool Multipliers::not_all_zero(size_t number_variables, double tolerance) const { + // constraint multipliers + for (double multiplier_j: this->constraints) { + if (tolerance < std::abs(multiplier_j)) { + return true; + } + } + // bound multipliers + for (size_t variable_index: Range(number_variables)) { + if (tolerance < std::abs(this->lower_bounds[variable_index] + this->upper_bounds[variable_index])) { + return true; + } + } + return false; + } +} // namespace \ No newline at end of file diff --git a/uno/optimization/Multipliers.hpp b/uno/optimization/Multipliers.hpp index 4d3cb464..9780fdf7 100644 --- a/uno/optimization/Multipliers.hpp +++ b/uno/optimization/Multipliers.hpp @@ -22,32 +22,6 @@ namespace uno { void reset(); [[nodiscard]] bool not_all_zero(size_t number_variables, double tolerance) const; }; - - inline Multipliers::Multipliers(size_t number_variables, size_t number_constraints) : lower_bounds(number_variables), - upper_bounds(number_variables), constraints(number_constraints) { - } - - inline void Multipliers::reset() { - this->constraints.fill(0.); - this->lower_bounds.fill(0.); - this->upper_bounds.fill(0.); - } - - inline bool Multipliers::not_all_zero(size_t number_variables, double tolerance) const { - // constraint multipliers - for (double multiplier_j: this->constraints) { - if (tolerance < std::abs(multiplier_j)) { - return true; - } - } - // bound multipliers - for (size_t variable_index: Range(number_variables)) { - if (tolerance < std::abs(this->lower_bounds[variable_index] + this->upper_bounds[variable_index])) { - return true; - } - } - return false; - } } // namespace #endif // UNO_MULTIPLIERS_H diff --git a/uno/optimization/WarmstartInformation.cpp b/uno/optimization/WarmstartInformation.cpp new file mode 100644 index 00000000..b51e43c4 --- /dev/null +++ b/uno/optimization/WarmstartInformation.cpp @@ -0,0 +1,47 @@ +// Copyright (c) 2018-2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include +#include "WarmstartInformation.hpp" + +namespace uno { + void WarmstartInformation::display() const { + std::cout << "Objective: " << std::boolalpha << this->objective_changed << '\n'; + std::cout << "Constraints: " << std::boolalpha << this->constraints_changed << '\n'; + std::cout << "Constraint bounds: " << std::boolalpha << this->constraint_bounds_changed << '\n'; + std::cout << "Variable bounds: " << std::boolalpha << this->variable_bounds_changed << '\n'; + std::cout << "Problem: " << std::boolalpha << this->problem_changed << '\n'; + } + + void WarmstartInformation::set_cold_start() { + this->objective_changed = true; + this->constraints_changed = true; + this->constraint_bounds_changed = true; + this->variable_bounds_changed = true; + this->problem_changed = true; + } + + void WarmstartInformation::set_hot_start() { + this->objective_changed = true; + this->constraints_changed = true; + this->constraint_bounds_changed = true; + this->variable_bounds_changed = true; + this->problem_changed = false; + } + + void WarmstartInformation::only_objective_changed() { + this->objective_changed = true; + this->constraints_changed = false; + this->constraint_bounds_changed = false; + this->variable_bounds_changed = false; + this->problem_changed = false; + } + + void WarmstartInformation::only_variable_bounds_changed() { + this->objective_changed = false; + this->constraints_changed = false; + this->constraint_bounds_changed = false; + this->variable_bounds_changed = true; + this->problem_changed = false; + } +} // namespace \ No newline at end of file diff --git a/uno/optimization/WarmstartInformation.hpp b/uno/optimization/WarmstartInformation.hpp index ca2f891e..5158df76 100644 --- a/uno/optimization/WarmstartInformation.hpp +++ b/uno/optimization/WarmstartInformation.hpp @@ -18,46 +18,6 @@ namespace uno { void only_objective_changed(); void only_variable_bounds_changed(); }; - - inline void WarmstartInformation::display() const { - std::cout << "Objective: " << std::boolalpha << this->objective_changed << '\n'; - std::cout << "Constraints: " << std::boolalpha << this->constraints_changed << '\n'; - std::cout << "Constraint bounds: " << std::boolalpha << this->constraint_bounds_changed << '\n'; - std::cout << "Variable bounds: " << std::boolalpha << this->variable_bounds_changed << '\n'; - std::cout << "Problem: " << std::boolalpha << this->problem_changed << '\n'; - } - - inline void WarmstartInformation::set_cold_start() { - this->objective_changed = true; - this->constraints_changed = true; - this->constraint_bounds_changed = true; - this->variable_bounds_changed = true; - this->problem_changed = true; - } - - inline void WarmstartInformation::set_hot_start() { - this->objective_changed = true; - this->constraints_changed = true; - this->constraint_bounds_changed = true; - this->variable_bounds_changed = true; - this->problem_changed = false; - } - - inline void WarmstartInformation::only_objective_changed() { - this->objective_changed = true; - this->constraints_changed = false; - this->constraint_bounds_changed = false; - this->variable_bounds_changed = false; - this->problem_changed = false; - } - - inline void WarmstartInformation::only_variable_bounds_changed() { - this->objective_changed = false; - this->constraints_changed = false; - this->constraint_bounds_changed = false; - this->variable_bounds_changed = true; - this->problem_changed = false; - } } // namespace #endif // UNO_WARMSTARTINFORMATION_H diff --git a/uno/preprocessing/Preprocessing.cpp b/uno/preprocessing/Preprocessing.cpp index 4c211e93..6ae7be16 100644 --- a/uno/preprocessing/Preprocessing.cpp +++ b/uno/preprocessing/Preprocessing.cpp @@ -2,7 +2,7 @@ // Licensed under the MIT license. See LICENSE file in the project directory for details. #include "Preprocessing.hpp" -#include "ingredients/subproblems/Direction.hpp" +#include "optimization/Direction.hpp" #include "linear_algebra/SymmetricMatrix.hpp" #include "linear_algebra/CSCSparseStorage.hpp" #include "linear_algebra/RectangularMatrix.hpp" diff --git a/uno/reformulation/ElasticVariables.cpp b/uno/reformulation/ElasticVariables.cpp new file mode 100644 index 00000000..72b50b08 --- /dev/null +++ b/uno/reformulation/ElasticVariables.cpp @@ -0,0 +1,46 @@ +// Copyright (c) 2018-2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "ElasticVariables.hpp" +#include "model/Model.hpp" +#include "tools/Infinity.hpp" + +namespace uno { + ElasticVariables::ElasticVariables(size_t number_positive_variables, size_t number_negative_variables): + positive(number_positive_variables), negative(number_negative_variables) { } + + ElasticVariables ElasticVariables::generate(const Model& model) { + const ElasticVariablesSizes sizes = ElasticVariables::count(model); + ElasticVariables elastic_variables(sizes.positive, sizes.negative); + + // generate elastic variables to relax the constraints + size_t elastic_index = model.number_variables; + for (size_t constraint_index: Range(model.number_constraints)) { + if (is_finite(model.constraint_upper_bound(constraint_index))) { + // nonnegative variable p that captures the positive part of the constraint violation + elastic_variables.positive.insert(constraint_index, elastic_index); + elastic_index++; + } + if (is_finite(model.constraint_lower_bound(constraint_index))) { + // nonpositive variable n that captures the negative part of the constraint violation + elastic_variables.negative.insert(constraint_index, elastic_index); + elastic_index++; + } + } + return elastic_variables; + } + + ElasticVariablesSizes ElasticVariables::count(const Model& model) { + ElasticVariablesSizes number_elastic_variables; + // if the subproblem uses slack variables, the bounds of the constraints are [0, 0] + for (size_t constraint_index: Range(model.number_constraints)) { + if (is_finite(model.constraint_upper_bound(constraint_index))) { + number_elastic_variables.positive++; + } + if (is_finite(model.constraint_lower_bound(constraint_index))) { + number_elastic_variables.negative++; + } + } + return number_elastic_variables; + } +} // namespace \ No newline at end of file diff --git a/uno/reformulation/ElasticVariables.hpp b/uno/reformulation/ElasticVariables.hpp index 1fd31963..7c1c7c15 100644 --- a/uno/reformulation/ElasticVariables.hpp +++ b/uno/reformulation/ElasticVariables.hpp @@ -4,9 +4,12 @@ #ifndef UNO_ELASTICVARIABLES_H #define UNO_ELASTICVARIABLES_H -#include "tools/Infinity.hpp" +#include "linear_algebra/SparseVector.hpp" namespace uno { + // forward declaration + class Model; + struct ElasticVariablesSizes { size_t positive{0}; size_t negative{0}; @@ -25,44 +28,6 @@ namespace uno { protected: static ElasticVariablesSizes count(const Model& model); }; - - inline ElasticVariables::ElasticVariables(size_t number_positive_variables, size_t number_negative_variables): - positive(number_positive_variables), negative(number_negative_variables) { } - - inline ElasticVariables ElasticVariables::generate(const Model& model) { - const ElasticVariablesSizes sizes = ElasticVariables::count(model); - ElasticVariables elastic_variables(sizes.positive, sizes.negative); - - // generate elastic variables to relax the constraints - size_t elastic_index = model.number_variables; - for (size_t constraint_index: Range(model.number_constraints)) { - if (is_finite(model.constraint_upper_bound(constraint_index))) { - // nonnegative variable p that captures the positive part of the constraint violation - elastic_variables.positive.insert(constraint_index, elastic_index); - elastic_index++; - } - if (is_finite(model.constraint_lower_bound(constraint_index))) { - // nonpositive variable n that captures the negative part of the constraint violation - elastic_variables.negative.insert(constraint_index, elastic_index); - elastic_index++; - } - } - return elastic_variables; - } - - inline ElasticVariablesSizes ElasticVariables::count(const Model& model) { - ElasticVariablesSizes number_elastic_variables; - // if the subproblem uses slack variables, the bounds of the constraints are [0, 0] - for (size_t constraint_index: Range(model.number_constraints)) { - if (is_finite(model.constraint_upper_bound(constraint_index))) { - number_elastic_variables.positive++; - } - if (is_finite(model.constraint_lower_bound(constraint_index))) { - number_elastic_variables.negative++; - } - } - return number_elastic_variables; - } } // namespace #endif // UNO_ELASTICVARIABLES_H diff --git a/uno/reformulation/OptimalityProblem.cpp b/uno/reformulation/OptimalityProblem.cpp new file mode 100644 index 00000000..29aa94bb --- /dev/null +++ b/uno/reformulation/OptimalityProblem.cpp @@ -0,0 +1,90 @@ +// Copyright (c) 2018-2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "OptimalityProblem.hpp" +#include "optimization/Iterate.hpp" +#include "optimization/LagrangianGradient.hpp" +#include "symbolic/Expression.hpp" + +namespace uno { + OptimalityProblem::OptimalityProblem(const Model& model): OptimizationProblem(model, model.number_variables, model.number_constraints) { + } + + void OptimalityProblem::evaluate_objective_gradient(Iterate& iterate, SparseVector& objective_gradient) const { + iterate.evaluate_objective_gradient(this->model); + // TODO change this + objective_gradient = iterate.evaluations.objective_gradient; + } + + void OptimalityProblem::evaluate_constraints(Iterate& iterate, std::vector& constraints) const { + iterate.evaluate_constraints(this->model); + constraints = iterate.evaluations.constraints; + } + + void OptimalityProblem::evaluate_constraint_jacobian(Iterate& iterate, RectangularMatrix& constraint_jacobian) const { + iterate.evaluate_constraint_jacobian(this->model); + // TODO change this + constraint_jacobian = iterate.evaluations.constraint_jacobian; + } + + void OptimalityProblem::evaluate_lagrangian_hessian(const Vector& x, const Vector& multipliers, + SymmetricMatrix& hessian) const { + this->model.evaluate_lagrangian_hessian(x, this->get_objective_multiplier(), multipliers, hessian); + } + + // Lagrangian gradient split in two parts: objective contribution and constraints' contribution + void OptimalityProblem::evaluate_lagrangian_gradient(LagrangianGradient& lagrangian_gradient, Iterate& iterate, + const Multipliers& multipliers) const { + lagrangian_gradient.objective_contribution.fill(0.); + lagrangian_gradient.constraints_contribution.fill(0.); + + // objective gradient + for (auto [variable_index, derivative]: iterate.evaluations.objective_gradient) { + lagrangian_gradient.objective_contribution[variable_index] += derivative; + } + + // constraints + for (size_t constraint_index: Range(this->number_constraints)) { + if (multipliers.constraints[constraint_index] != 0.) { + for (auto [variable_index, derivative]: iterate.evaluations.constraint_jacobian[constraint_index]) { + lagrangian_gradient.constraints_contribution[variable_index] -= multipliers.constraints[constraint_index] * derivative; + } + } + } + + // bound constraints of original variables + for (size_t variable_index: Range(this->number_variables)) { + lagrangian_gradient.constraints_contribution[variable_index] -= (multipliers.lower_bounds[variable_index] + + multipliers.upper_bounds[variable_index]); + } + } + + double OptimalityProblem::complementarity_error(const Vector& primals, const std::vector& constraints, + const Multipliers& multipliers, double shift_value, Norm residual_norm) const { + // bound constraints + const Range variables_range = Range(this->model.number_variables); + const VectorExpression variable_complementarity{variables_range, [&](size_t variable_index) { + if (0. < multipliers.lower_bounds[variable_index]) { + return multipliers.lower_bounds[variable_index] * (primals[variable_index] - this->model.variable_lower_bound(variable_index)) - shift_value; + } + if (multipliers.upper_bounds[variable_index] < 0.) { + return multipliers.upper_bounds[variable_index] * (primals[variable_index] - this->model.variable_upper_bound(variable_index)) - shift_value; + } + return 0.; + }}; + + // inequality constraints + const VectorExpression constraint_complementarity{this->model.get_inequality_constraints(), [&](size_t constraint_index) { + if (0. < multipliers.constraints[constraint_index]) { // lower bound + return multipliers.constraints[constraint_index] * (constraints[constraint_index] - this->model.constraint_lower_bound(constraint_index)) - + shift_value; + } + else if (multipliers.constraints[constraint_index] < 0.) { // upper bound + return multipliers.constraints[constraint_index] * (constraints[constraint_index] - this->model.constraint_upper_bound(constraint_index)) - + shift_value; + } + return 0.; + }}; + return norm(residual_norm, variable_complementarity, constraint_complementarity); + } +} // namespace \ No newline at end of file diff --git a/uno/reformulation/OptimalityProblem.hpp b/uno/reformulation/OptimalityProblem.hpp index 1040d48f..9f488dae 100644 --- a/uno/reformulation/OptimalityProblem.hpp +++ b/uno/reformulation/OptimalityProblem.hpp @@ -5,9 +5,6 @@ #define UNO_OPTIMALITYPROBLEM_H #include "OptimizationProblem.hpp" -#include "optimization/Iterate.hpp" -#include "optimization/LagrangianGradient.hpp" -#include "symbolic/Expression.hpp" namespace uno { class OptimalityProblem: public OptimizationProblem { @@ -38,87 +35,6 @@ namespace uno { [[nodiscard]] double complementarity_error(const Vector& primals, const std::vector& constraints, const Multipliers& multipliers, double shift_value, Norm residual_norm) const override; }; - - inline OptimalityProblem::OptimalityProblem(const Model& model): OptimizationProblem(model, model.number_variables, model.number_constraints) { - } - - inline void OptimalityProblem::evaluate_objective_gradient(Iterate& iterate, SparseVector& objective_gradient) const { - iterate.evaluate_objective_gradient(this->model); - // TODO change this - objective_gradient = iterate.evaluations.objective_gradient; - } - - inline void OptimalityProblem::evaluate_constraints(Iterate& iterate, std::vector& constraints) const { - iterate.evaluate_constraints(this->model); - constraints = iterate.evaluations.constraints; - } - - inline void OptimalityProblem::evaluate_constraint_jacobian(Iterate& iterate, RectangularMatrix& constraint_jacobian) const { - iterate.evaluate_constraint_jacobian(this->model); - // TODO change this - constraint_jacobian = iterate.evaluations.constraint_jacobian; - } - - inline void OptimalityProblem::evaluate_lagrangian_hessian(const Vector& x, const Vector& multipliers, - SymmetricMatrix& hessian) const { - this->model.evaluate_lagrangian_hessian(x, this->get_objective_multiplier(), multipliers, hessian); - } - - // Lagrangian gradient split in two parts: objective contribution and constraints' contribution - inline void OptimalityProblem::evaluate_lagrangian_gradient(LagrangianGradient& lagrangian_gradient, Iterate& iterate, - const Multipliers& multipliers) const { - lagrangian_gradient.objective_contribution.fill(0.); - lagrangian_gradient.constraints_contribution.fill(0.); - - // objective gradient - for (auto [variable_index, derivative]: iterate.evaluations.objective_gradient) { - lagrangian_gradient.objective_contribution[variable_index] += derivative; - } - - // constraints - for (size_t constraint_index: Range(this->number_constraints)) { - if (multipliers.constraints[constraint_index] != 0.) { - for (auto [variable_index, derivative]: iterate.evaluations.constraint_jacobian[constraint_index]) { - lagrangian_gradient.constraints_contribution[variable_index] -= multipliers.constraints[constraint_index] * derivative; - } - } - } - - // bound constraints of original variables - for (size_t variable_index: Range(this->number_variables)) { - lagrangian_gradient.constraints_contribution[variable_index] -= (multipliers.lower_bounds[variable_index] + - multipliers.upper_bounds[variable_index]); - } - } - - inline double OptimalityProblem::complementarity_error(const Vector& primals, const std::vector& constraints, - const Multipliers& multipliers, double shift_value, Norm residual_norm) const { - // bound constraints - const Range variables_range = Range(this->model.number_variables); - const VectorExpression variable_complementarity{variables_range, [&](size_t variable_index) { - if (0. < multipliers.lower_bounds[variable_index]) { - return multipliers.lower_bounds[variable_index] * (primals[variable_index] - this->model.variable_lower_bound(variable_index)) - shift_value; - } - if (multipliers.upper_bounds[variable_index] < 0.) { - return multipliers.upper_bounds[variable_index] * (primals[variable_index] - this->model.variable_upper_bound(variable_index)) - shift_value; - } - return 0.; - }}; - - // inequality constraints - const VectorExpression constraint_complementarity{this->model.get_inequality_constraints(), [&](size_t constraint_index) { - if (0. < multipliers.constraints[constraint_index]) { // lower bound - return multipliers.constraints[constraint_index] * (constraints[constraint_index] - this->model.constraint_lower_bound(constraint_index)) - - shift_value; - } - else if (multipliers.constraints[constraint_index] < 0.) { // upper bound - return multipliers.constraints[constraint_index] * (constraints[constraint_index] - this->model.constraint_upper_bound(constraint_index)) - - shift_value; - } - return 0.; - }}; - return norm(residual_norm, variable_complementarity, constraint_complementarity); - } } // namespace #endif // UNO_OPTIMALITYPROBLEM_H diff --git a/uno/reformulation/OptimizationProblem.cpp b/uno/reformulation/OptimizationProblem.cpp new file mode 100644 index 00000000..210addd6 --- /dev/null +++ b/uno/reformulation/OptimizationProblem.cpp @@ -0,0 +1,33 @@ +// Copyright (c) 2018-2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "OptimizationProblem.hpp" + +namespace uno { + OptimizationProblem::OptimizationProblem(const Model& model, size_t number_variables, size_t number_constraints): + model(model), number_variables(number_variables), number_constraints(number_constraints) { + } + + bool OptimizationProblem::is_constrained() const { + return (0 < this->number_constraints); + } + + bool OptimizationProblem::has_inequality_constraints() const { + return (not this->model.get_inequality_constraints().empty()); + } + + bool OptimizationProblem::has_fixed_variables() const { + return (not this->model.get_fixed_variables().empty()); + } + + size_t OptimizationProblem::get_number_original_variables() const { + return this->model.number_variables; + } + + double OptimizationProblem::stationarity_error(const LagrangianGradient& lagrangian_gradient, double objective_multiplier, + Norm residual_norm) { + // norm of the scaled Lagrangian gradient + const auto scaled_lagrangian = objective_multiplier * lagrangian_gradient.objective_contribution + lagrangian_gradient.constraints_contribution; + return norm(residual_norm, scaled_lagrangian); + } +} // namespace \ No newline at end of file diff --git a/uno/reformulation/OptimizationProblem.hpp b/uno/reformulation/OptimizationProblem.hpp index 74e19e6d..cb01c884 100644 --- a/uno/reformulation/OptimizationProblem.hpp +++ b/uno/reformulation/OptimizationProblem.hpp @@ -63,33 +63,6 @@ namespace uno { [[nodiscard]] virtual double complementarity_error(const Vector& primals, const std::vector& constraints, const Multipliers& multipliers, double shift_value, Norm residual_norm) const = 0; }; - - inline OptimizationProblem::OptimizationProblem(const Model& model, size_t number_variables, size_t number_constraints): - model(model), number_variables(number_variables), number_constraints(number_constraints) { - } - - inline bool OptimizationProblem::is_constrained() const { - return (0 < this->number_constraints); - } - - inline bool OptimizationProblem::has_inequality_constraints() const { - return (not this->model.get_inequality_constraints().empty()); - } - - inline bool OptimizationProblem::has_fixed_variables() const { - return (not this->model.get_fixed_variables().empty()); - } - - inline size_t OptimizationProblem::get_number_original_variables() const { - return this->model.number_variables; - } - - inline double OptimizationProblem::stationarity_error(const LagrangianGradient& lagrangian_gradient, double objective_multiplier, - Norm residual_norm) { - // norm of the scaled Lagrangian gradient - const auto scaled_lagrangian = objective_multiplier * lagrangian_gradient.objective_contribution + lagrangian_gradient.constraints_contribution; - return norm(residual_norm, scaled_lagrangian); - } } // namespace #endif // UNO_OPTIMIZATIONPROBLEM_H \ No newline at end of file diff --git a/uno/reformulation/l1RelaxedProblem.cpp b/uno/reformulation/l1RelaxedProblem.cpp new file mode 100644 index 00000000..06f98e2c --- /dev/null +++ b/uno/reformulation/l1RelaxedProblem.cpp @@ -0,0 +1,281 @@ +// Copyright (c) 2018-2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "l1RelaxedProblem.hpp" +#include "linear_algebra/SymmetricMatrix.hpp" +#include "model/Model.hpp" +#include "optimization/Iterate.hpp" +#include "optimization/LagrangianGradient.hpp" +#include "symbolic/Expression.hpp" +#include "symbolic/VectorExpression.hpp" +#include "symbolic/Concatenation.hpp" +#include "tools/Infinity.hpp" + +namespace uno { + l1RelaxedProblem::l1RelaxedProblem(const Model& model, double objective_multiplier, double constraint_violation_coefficient, + double proximal_coefficient, double const* proximal_center): + // call delegating constructor + l1RelaxedProblem(model, ElasticVariables::generate(model), objective_multiplier, constraint_violation_coefficient, proximal_coefficient, + proximal_center) { + } + + // private delegating constructor + l1RelaxedProblem::l1RelaxedProblem(const Model& model, ElasticVariables&& elastic_variables, double objective_multiplier, + double constraint_violation_coefficient, double proximal_coefficient, double const* proximal_center): + OptimizationProblem(model, model.number_variables + elastic_variables.size(), model.number_constraints), + objective_multiplier(objective_multiplier), + constraint_violation_coefficient(constraint_violation_coefficient), + proximal_coefficient(proximal_coefficient), + proximal_center(proximal_center), + elastic_variables(std::forward(elastic_variables)), + // lower bounded variables are the model variables + the elastic variables + lower_bounded_variables(concatenate(this->model.get_lower_bounded_variables(), Range(model.number_variables, + model.number_variables + this->elastic_variables.size()))), + single_lower_bounded_variables(concatenate(this->model.get_single_lower_bounded_variables(), + Range(model.number_variables, model.number_variables + this->elastic_variables.size()))) { + } + + double l1RelaxedProblem::get_objective_multiplier() const { + return this->objective_multiplier; + } + + void l1RelaxedProblem::evaluate_objective_gradient(Iterate& iterate, SparseVector& objective_gradient) const { + // scale nabla f(x) by rho + if (this->objective_multiplier != 0.) { + iterate.evaluate_objective_gradient(this->model); + // TODO change this + objective_gradient = iterate.evaluations.objective_gradient; + scale(objective_gradient, this->objective_multiplier); + } + else { + objective_gradient.clear(); + } + + // constraint violation (through elastic variables) contribution + for (const auto [_, elastic_index]: this->elastic_variables.positive) { + objective_gradient.insert(elastic_index, this->constraint_violation_coefficient); + } + for (const auto [_, elastic_index]: this->elastic_variables.negative) { + objective_gradient.insert(elastic_index, this->constraint_violation_coefficient); + } + + // proximal contribution + if (this->proximal_center != nullptr && this->proximal_coefficient != 0.) { + for (size_t variable_index: Range(this->model.number_variables)) { + const double scaling = std::min(1., 1./std::abs(this->proximal_center[variable_index])); + const double proximal_term = this->proximal_coefficient * scaling * scaling * (iterate.primals[variable_index] - this->proximal_center[variable_index]); + objective_gradient.insert(variable_index, proximal_term); + } + } + } + + void l1RelaxedProblem::evaluate_constraints(Iterate& iterate, std::vector& constraints) const { + iterate.evaluate_constraints(this->model); + constraints = iterate.evaluations.constraints; + // add the contribution of the elastics + for (const auto [constraint_index, elastic_index]: this->elastic_variables.positive) { + constraints[constraint_index] -= iterate.primals[elastic_index]; + } + for (const auto [constraint_index, elastic_index]: this->elastic_variables.negative) { + constraints[constraint_index] += iterate.primals[elastic_index]; + } + } + + void l1RelaxedProblem::evaluate_constraint_jacobian(Iterate& iterate, RectangularMatrix& constraint_jacobian) const { + iterate.evaluate_constraint_jacobian(this->model); + // TODO change this + constraint_jacobian = iterate.evaluations.constraint_jacobian; + // add the contribution of the elastics + for (const auto [constraint_index, elastic_index]: this->elastic_variables.positive) { + constraint_jacobian[constraint_index].insert(elastic_index, -1.); + } + for (const auto [constraint_index, elastic_index]: this->elastic_variables.negative) { + constraint_jacobian[constraint_index].insert(elastic_index, 1.); + } + } + + void l1RelaxedProblem::evaluate_lagrangian_hessian(const Vector& x, const Vector& multipliers, + SymmetricMatrix& hessian) const { + this->model.evaluate_lagrangian_hessian(x, this->objective_multiplier, multipliers, hessian); + + // proximal contribution + if (this->proximal_center != nullptr && this->proximal_coefficient != 0.) { + for (size_t variable_index: Range(this->model.number_variables)) { + const double scaling = std::min(1., 1./std::abs(this->proximal_center[variable_index])); + const double proximal_term = this->proximal_coefficient * scaling * scaling; + hessian.insert(proximal_term, variable_index, variable_index); + } + } + + // extend the dimension of the Hessian by finalizing the remaining columns (note: the elastics do not enter the Hessian) + for (size_t constraint_index: Range(this->model.number_variables, this->number_variables)) { + hessian.finalize_column(constraint_index); + } + } + + // Lagrangian gradient split in two parts: objective contribution and constraints' contribution + void l1RelaxedProblem::evaluate_lagrangian_gradient(LagrangianGradient& lagrangian_gradient, Iterate& iterate, + const Multipliers& multipliers) const { + lagrangian_gradient.objective_contribution.fill(0.); + lagrangian_gradient.constraints_contribution.fill(0.); + + // objective gradient + for (auto [variable_index, derivative]: iterate.evaluations.objective_gradient) { + lagrangian_gradient.objective_contribution[variable_index] += derivative; + } + + // constraints + for (size_t constraint_index: Range(this->number_constraints)) { + if (multipliers.constraints[constraint_index] != 0.) { + for (auto [variable_index, derivative]: iterate.evaluations.constraint_jacobian[constraint_index]) { + lagrangian_gradient.constraints_contribution[variable_index] -= multipliers.constraints[constraint_index] * derivative; + } + } + } + + // bound constraints of original variables + for (size_t variable_index: Range(this->model.number_variables)) { + lagrangian_gradient.constraints_contribution[variable_index] -= (multipliers.lower_bounds[variable_index] + + multipliers.upper_bounds[variable_index]); + } + + // elastic variables + for (const auto [constraint_index, elastic_index]: this->elastic_variables.positive) { + lagrangian_gradient.constraints_contribution[elastic_index] += this->constraint_violation_coefficient + + multipliers.constraints[constraint_index] - multipliers.lower_bounds[elastic_index]; + } + for (const auto [constraint_index, elastic_index]: this->elastic_variables.negative) { + lagrangian_gradient.constraints_contribution[elastic_index] += this->constraint_violation_coefficient - + multipliers.constraints[constraint_index] - multipliers.lower_bounds[elastic_index]; + } + + // proximal contribution + if (this->proximal_center != nullptr && this->proximal_coefficient != 0.) { + for (size_t variable_index: Range(this->model.number_variables)) { + const double scaling = std::min(1., 1./std::abs(this->proximal_center[variable_index])); + const double proximal_term = this->proximal_coefficient * scaling * scaling; + lagrangian_gradient.constraints_contribution[variable_index] += proximal_term * (iterate.primals[variable_index] - + this->proximal_center[variable_index]); + } + } + } + + // complementary slackness error: expression for violated constraints depends on the definition of the relaxed problem + double l1RelaxedProblem::complementarity_error(const Vector& primals, const std::vector& constraints, + const Multipliers& multipliers, double shift_value, Norm residual_norm) const { + // bound constraints + const Range variables_range = Range(this->number_variables); + const VectorExpression bounds_complementarity{variables_range, [&](size_t variable_index) { + if (0. < multipliers.lower_bounds[variable_index]) { + return multipliers.lower_bounds[variable_index] * (primals[variable_index] - this->variable_lower_bound(variable_index)) - shift_value; + } + if (multipliers.upper_bounds[variable_index] < 0.) { + return multipliers.upper_bounds[variable_index] * (primals[variable_index] - this->variable_upper_bound(variable_index)) - shift_value; + } + return 0.; + }}; + + // general constraints + // TODO use the values of the relaxed constraints + const Range constraints_range = Range(this->number_constraints); + const VectorExpression constraints_complementarity{constraints_range, [&](size_t constraint_index) { + if (this->model.get_constraint_bound_type(constraint_index) != EQUAL_BOUNDS) { + if (0. < multipliers.constraints[constraint_index]) { // lower bound + return multipliers.constraints[constraint_index] * (constraints[constraint_index] - this->constraint_lower_bound(constraint_index)) - shift_value; + } + else if (multipliers.constraints[constraint_index] < 0.) { // upper bound + return multipliers.constraints[constraint_index] * (constraints[constraint_index] - this->constraint_upper_bound(constraint_index)) - shift_value; + } + } + return 0.; + }}; + return norm(residual_norm, bounds_complementarity); + } + + double l1RelaxedProblem::variable_lower_bound(size_t variable_index) const { + if (variable_index < this->model.number_variables) { // model variable + return this->model.variable_lower_bound(variable_index); + } + else { // elastic variable in [0, +inf[ + return 0.; + } + } + + double l1RelaxedProblem::variable_upper_bound(size_t variable_index) const { + if (variable_index < this->model.number_variables) { // model variable + return this->model.variable_upper_bound(variable_index); + } + else { // elastic variable in [0, +inf[ + return INF; + } + } + + double l1RelaxedProblem::constraint_lower_bound(size_t constraint_index) const { + return this->model.constraint_lower_bound(constraint_index); + } + + double l1RelaxedProblem::constraint_upper_bound(size_t constraint_index) const { + return this->model.constraint_upper_bound(constraint_index); + } + + const Collection& l1RelaxedProblem::get_lower_bounded_variables() const { + return this->lower_bounded_variables; + } + + const Collection& l1RelaxedProblem::get_upper_bounded_variables() const { + // same set as the model + return this->model.get_upper_bounded_variables(); + } + + const Collection& l1RelaxedProblem::get_single_lower_bounded_variables() const { + return this->single_lower_bounded_variables; + } + + const Collection& l1RelaxedProblem::get_single_upper_bounded_variables() const { + // same set as the model + return this->model.get_single_upper_bounded_variables(); + } + + size_t l1RelaxedProblem::number_objective_gradient_nonzeros() const { + // elastic contribution + size_t number_nonzeros = this->elastic_variables.size(); + + // objective contribution + if (this->objective_multiplier != 0.) { + number_nonzeros += this->model.number_objective_gradient_nonzeros(); + } + return number_nonzeros; + } + + size_t l1RelaxedProblem::number_jacobian_nonzeros() const { + return this->model.number_jacobian_nonzeros() + this->elastic_variables.size(); + } + + size_t l1RelaxedProblem::number_hessian_nonzeros() const { + return this->model.number_hessian_nonzeros(); + } + + void l1RelaxedProblem::set_objective_multiplier(double new_objective_multiplier) { + assert(0. <= new_objective_multiplier && "The objective multiplier should be non-negative"); + this->objective_multiplier = new_objective_multiplier; + } + + void l1RelaxedProblem::set_proximal_multiplier(double new_proximal_coefficient) { + this->proximal_coefficient = new_proximal_coefficient; + } + + void l1RelaxedProblem::set_proximal_center(double const* new_proximal_center) { + this->proximal_center = new_proximal_center; + } + + void l1RelaxedProblem::set_elastic_variable_values(Iterate& iterate, const std::function& + elastic_setting_function) const { + iterate.set_number_variables(this->number_variables); + for (const auto [constraint_index, elastic_index]: this->elastic_variables.positive) { + elastic_setting_function(iterate, constraint_index, elastic_index, -1.); + } + for (const auto [constraint_index, elastic_index]: this->elastic_variables.negative) { + elastic_setting_function(iterate, constraint_index, elastic_index, 1.); + } + } +} // namespace \ No newline at end of file diff --git a/uno/reformulation/l1RelaxedProblem.hpp b/uno/reformulation/l1RelaxedProblem.hpp index 4a4cf069..e5ab02ed 100644 --- a/uno/reformulation/l1RelaxedProblem.hpp +++ b/uno/reformulation/l1RelaxedProblem.hpp @@ -6,15 +6,7 @@ #include "OptimizationProblem.hpp" #include "ElasticVariables.hpp" -#include "linear_algebra/SymmetricMatrix.hpp" -#include "model/Model.hpp" -#include "optimization/Iterate.hpp" -#include "optimization/LagrangianGradient.hpp" -#include "symbolic/Expression.hpp" -#include "symbolic/VectorExpression.hpp" #include "symbolic/Concatenation.hpp" -#include "symbolic/Range.hpp" -#include "tools/Infinity.hpp" namespace uno { class l1RelaxedProblem: public OptimizationProblem { @@ -66,273 +58,6 @@ namespace uno { l1RelaxedProblem(const Model& model, ElasticVariables&& elastic_variables, double objective_multiplier, double constraint_violation_coefficient, double proximal_coefficient, double const* proximal_center); }; - - inline l1RelaxedProblem::l1RelaxedProblem(const Model& model, double objective_multiplier, double constraint_violation_coefficient, - double proximal_coefficient, double const* proximal_center): - // call delegating constructor - l1RelaxedProblem(model, ElasticVariables::generate(model), objective_multiplier, constraint_violation_coefficient, proximal_coefficient, - proximal_center) { - } - - // private delegating constructor - inline l1RelaxedProblem::l1RelaxedProblem(const Model& model, ElasticVariables&& elastic_variables, double objective_multiplier, - double constraint_violation_coefficient, double proximal_coefficient, double const* proximal_center): - OptimizationProblem(model, model.number_variables + elastic_variables.size(), model.number_constraints), - objective_multiplier(objective_multiplier), - constraint_violation_coefficient(constraint_violation_coefficient), - proximal_coefficient(proximal_coefficient), - proximal_center(proximal_center), - elastic_variables(std::forward(elastic_variables)), - // lower bounded variables are the model variables + the elastic variables - lower_bounded_variables(concatenate(this->model.get_lower_bounded_variables(), Range(model.number_variables, - model.number_variables + this->elastic_variables.size()))), - single_lower_bounded_variables(concatenate(this->model.get_single_lower_bounded_variables(), - Range(model.number_variables, model.number_variables + this->elastic_variables.size()))) { - } - - inline double l1RelaxedProblem::get_objective_multiplier() const { - return this->objective_multiplier; - } - - inline void l1RelaxedProblem::evaluate_objective_gradient(Iterate& iterate, SparseVector& objective_gradient) const { - // scale nabla f(x) by rho - if (this->objective_multiplier != 0.) { - iterate.evaluate_objective_gradient(this->model); - // TODO change this - objective_gradient = iterate.evaluations.objective_gradient; - scale(objective_gradient, this->objective_multiplier); - } - else { - objective_gradient.clear(); - } - - // constraint violation (through elastic variables) contribution - for (const auto [_, elastic_index]: this->elastic_variables.positive) { - objective_gradient.insert(elastic_index, this->constraint_violation_coefficient); - } - for (const auto [_, elastic_index]: this->elastic_variables.negative) { - objective_gradient.insert(elastic_index, this->constraint_violation_coefficient); - } - - // proximal contribution - if (this->proximal_center != nullptr && this->proximal_coefficient != 0.) { - for (size_t variable_index: Range(this->model.number_variables)) { - const double scaling = std::min(1., 1./std::abs(this->proximal_center[variable_index])); - const double proximal_term = this->proximal_coefficient * scaling * scaling * (iterate.primals[variable_index] - this->proximal_center[variable_index]); - objective_gradient.insert(variable_index, proximal_term); - } - } - } - - inline void l1RelaxedProblem::evaluate_constraints(Iterate& iterate, std::vector& constraints) const { - iterate.evaluate_constraints(this->model); - constraints = iterate.evaluations.constraints; - // add the contribution of the elastics - for (const auto [constraint_index, elastic_index]: this->elastic_variables.positive) { - constraints[constraint_index] -= iterate.primals[elastic_index]; - } - for (const auto [constraint_index, elastic_index]: this->elastic_variables.negative) { - constraints[constraint_index] += iterate.primals[elastic_index]; - } - } - - inline void l1RelaxedProblem::evaluate_constraint_jacobian(Iterate& iterate, RectangularMatrix& constraint_jacobian) const { - iterate.evaluate_constraint_jacobian(this->model); - // TODO change this - constraint_jacobian = iterate.evaluations.constraint_jacobian; - // add the contribution of the elastics - for (const auto [constraint_index, elastic_index]: this->elastic_variables.positive) { - constraint_jacobian[constraint_index].insert(elastic_index, -1.); - } - for (const auto [constraint_index, elastic_index]: this->elastic_variables.negative) { - constraint_jacobian[constraint_index].insert(elastic_index, 1.); - } - } - - inline void l1RelaxedProblem::evaluate_lagrangian_hessian(const Vector& x, const Vector& multipliers, - SymmetricMatrix& hessian) const { - this->model.evaluate_lagrangian_hessian(x, this->objective_multiplier, multipliers, hessian); - - // proximal contribution - if (this->proximal_center != nullptr && this->proximal_coefficient != 0.) { - for (size_t variable_index: Range(this->model.number_variables)) { - const double scaling = std::min(1., 1./std::abs(this->proximal_center[variable_index])); - const double proximal_term = this->proximal_coefficient * scaling * scaling; - hessian.insert(proximal_term, variable_index, variable_index); - } - } - - // extend the dimension of the Hessian by finalizing the remaining columns (note: the elastics do not enter the Hessian) - for (size_t constraint_index: Range(this->model.number_variables, this->number_variables)) { - hessian.finalize_column(constraint_index); - } - } - - // Lagrangian gradient split in two parts: objective contribution and constraints' contribution - inline void l1RelaxedProblem::evaluate_lagrangian_gradient(LagrangianGradient& lagrangian_gradient, Iterate& iterate, - const Multipliers& multipliers) const { - lagrangian_gradient.objective_contribution.fill(0.); - lagrangian_gradient.constraints_contribution.fill(0.); - - // objective gradient - for (auto [variable_index, derivative]: iterate.evaluations.objective_gradient) { - lagrangian_gradient.objective_contribution[variable_index] += derivative; - } - - // constraints - for (size_t constraint_index: Range(this->number_constraints)) { - if (multipliers.constraints[constraint_index] != 0.) { - for (auto [variable_index, derivative]: iterate.evaluations.constraint_jacobian[constraint_index]) { - lagrangian_gradient.constraints_contribution[variable_index] -= multipliers.constraints[constraint_index] * derivative; - } - } - } - - // bound constraints of original variables - for (size_t variable_index: Range(this->model.number_variables)) { - lagrangian_gradient.constraints_contribution[variable_index] -= (multipliers.lower_bounds[variable_index] + - multipliers.upper_bounds[variable_index]); - } - - // elastic variables - for (const auto [constraint_index, elastic_index]: this->elastic_variables.positive) { - lagrangian_gradient.constraints_contribution[elastic_index] += this->constraint_violation_coefficient + - multipliers.constraints[constraint_index] - multipliers.lower_bounds[elastic_index]; - } - for (const auto [constraint_index, elastic_index]: this->elastic_variables.negative) { - lagrangian_gradient.constraints_contribution[elastic_index] += this->constraint_violation_coefficient - - multipliers.constraints[constraint_index] - multipliers.lower_bounds[elastic_index]; - } - - // proximal contribution - if (this->proximal_center != nullptr && this->proximal_coefficient != 0.) { - for (size_t variable_index: Range(this->model.number_variables)) { - const double scaling = std::min(1., 1./std::abs(this->proximal_center[variable_index])); - const double proximal_term = this->proximal_coefficient * scaling * scaling; - lagrangian_gradient.constraints_contribution[variable_index] += proximal_term * (iterate.primals[variable_index] - - this->proximal_center[variable_index]); - } - } - } - - // complementary slackness error: expression for violated constraints depends on the definition of the relaxed problem - inline double l1RelaxedProblem::complementarity_error(const Vector& primals, const std::vector& constraints, - const Multipliers& multipliers, double shift_value, Norm residual_norm) const { - // bound constraints - const Range variables_range = Range(this->number_variables); - const VectorExpression bounds_complementarity{variables_range, [&](size_t variable_index) { - if (0. < multipliers.lower_bounds[variable_index]) { - return multipliers.lower_bounds[variable_index] * (primals[variable_index] - this->variable_lower_bound(variable_index)) - shift_value; - } - if (multipliers.upper_bounds[variable_index] < 0.) { - return multipliers.upper_bounds[variable_index] * (primals[variable_index] - this->variable_upper_bound(variable_index)) - shift_value; - } - return 0.; - }}; - - // general constraints - // TODO use the values of the relaxed constraints - const Range constraints_range = Range(this->number_constraints); - const VectorExpression constraints_complementarity{constraints_range, [&](size_t constraint_index) { - if (this->model.get_constraint_bound_type(constraint_index) != EQUAL_BOUNDS) { - if (0. < multipliers.constraints[constraint_index]) { // lower bound - return multipliers.constraints[constraint_index] * (constraints[constraint_index] - this->constraint_lower_bound(constraint_index)) - shift_value; - } - else if (multipliers.constraints[constraint_index] < 0.) { // upper bound - return multipliers.constraints[constraint_index] * (constraints[constraint_index] - this->constraint_upper_bound(constraint_index)) - shift_value; - } - } - return 0.; - }}; - return norm(residual_norm, bounds_complementarity); - } - - inline double l1RelaxedProblem::variable_lower_bound(size_t variable_index) const { - if (variable_index < this->model.number_variables) { // model variable - return this->model.variable_lower_bound(variable_index); - } - else { // elastic variable in [0, +inf[ - return 0.; - } - } - - inline double l1RelaxedProblem::variable_upper_bound(size_t variable_index) const { - if (variable_index < this->model.number_variables) { // model variable - return this->model.variable_upper_bound(variable_index); - } - else { // elastic variable in [0, +inf[ - return INF; - } - } - - inline double l1RelaxedProblem::constraint_lower_bound(size_t constraint_index) const { - return this->model.constraint_lower_bound(constraint_index); - } - - inline double l1RelaxedProblem::constraint_upper_bound(size_t constraint_index) const { - return this->model.constraint_upper_bound(constraint_index); - } - - inline const Collection& l1RelaxedProblem::get_lower_bounded_variables() const { - return this->lower_bounded_variables; - } - - inline const Collection& l1RelaxedProblem::get_upper_bounded_variables() const { - // same set as the model - return this->model.get_upper_bounded_variables(); - } - - inline const Collection& l1RelaxedProblem::get_single_lower_bounded_variables() const { - return this->single_lower_bounded_variables; - } - - inline const Collection& l1RelaxedProblem::get_single_upper_bounded_variables() const { - // same set as the model - return this->model.get_single_upper_bounded_variables(); - } - - inline size_t l1RelaxedProblem::number_objective_gradient_nonzeros() const { - // elastic contribution - size_t number_nonzeros = this->elastic_variables.size(); - - // objective contribution - if (this->objective_multiplier != 0.) { - number_nonzeros += this->model.number_objective_gradient_nonzeros(); - } - return number_nonzeros; - } - - inline size_t l1RelaxedProblem::number_jacobian_nonzeros() const { - return this->model.number_jacobian_nonzeros() + this->elastic_variables.size(); - } - - inline size_t l1RelaxedProblem::number_hessian_nonzeros() const { - return this->model.number_hessian_nonzeros(); - } - - inline void l1RelaxedProblem::set_objective_multiplier(double new_objective_multiplier) { - assert(0. <= new_objective_multiplier && "The objective multiplier should be non-negative"); - this->objective_multiplier = new_objective_multiplier; - } - - inline void l1RelaxedProblem::set_proximal_multiplier(double new_proximal_coefficient) { - this->proximal_coefficient = new_proximal_coefficient; - } - - inline void l1RelaxedProblem::set_proximal_center(double const* new_proximal_center) { - this->proximal_center = new_proximal_center; - } - - inline void l1RelaxedProblem::set_elastic_variable_values(Iterate& iterate, const std::function& - elastic_setting_function) const { - iterate.set_number_variables(this->number_variables); - for (const auto [constraint_index, elastic_index]: this->elastic_variables.positive) { - elastic_setting_function(iterate, constraint_index, elastic_index, -1.); - } - for (const auto [constraint_index, elastic_index]: this->elastic_variables.negative) { - elastic_setting_function(iterate, constraint_index, elastic_index, 1.); - } - } } // namespace #endif // UNO_L1RELAXEDPROBLEM_H diff --git a/uno/solvers/BQPD/BQPDSolver.cpp b/uno/solvers/BQPD/BQPDSolver.cpp index f42edfdb..6d54e4d2 100644 --- a/uno/solvers/BQPD/BQPDSolver.cpp +++ b/uno/solvers/BQPD/BQPDSolver.cpp @@ -4,7 +4,7 @@ #include #include #include "BQPDSolver.hpp" -#include "ingredients/subproblems/Direction.hpp" +#include "optimization/Direction.hpp" #include "linear_algebra/RectangularMatrix.hpp" #include "linear_algebra/SymmetricMatrix.hpp" #include "linear_algebra/Vector.hpp" @@ -285,28 +285,28 @@ namespace uno { case BQPDStatus::UNBOUNDED_PROBLEM: return SubproblemStatus::UNBOUNDED_PROBLEM; case BQPDStatus::BOUND_INCONSISTENCY: - DEBUG << YELLOW << "BQPD error: bound inconsistency\n" << RESET; + DEBUG << "BQPD error: bound inconsistency\n"; return SubproblemStatus::ERROR; case BQPDStatus::INFEASIBLE: return SubproblemStatus::INFEASIBLE; // errors case BQPDStatus::INCORRECT_PARAMETER: - DEBUG << YELLOW << "BQPD error: incorrect parameter\n" << RESET; + DEBUG << "BQPD error: incorrect parameter\n"; return SubproblemStatus::ERROR; case BQPDStatus::LP_INSUFFICIENT_SPACE: - DEBUG << YELLOW << "BQPD error: LP insufficient space\n" << RESET; + DEBUG << "BQPD error: LP insufficient space\n"; return SubproblemStatus::ERROR; case BQPDStatus::HESSIAN_INSUFFICIENT_SPACE: - DEBUG << YELLOW << "BQPD kmax too small, continue anyway\n" << RESET; + DEBUG << "BQPD kmax too small, continue anyway\n"; return SubproblemStatus::ERROR; case BQPDStatus::SPARSE_INSUFFICIENT_SPACE: - DEBUG << YELLOW << "BQPD error: sparse insufficient space\n" << RESET; + DEBUG << "BQPD error: sparse insufficient space\n"; return SubproblemStatus::ERROR; case BQPDStatus::MAX_RESTARTS_REACHED: - DEBUG << YELLOW << "BQPD max restarts reached\n" << RESET; + DEBUG << "BQPD max restarts reached\n"; return SubproblemStatus::ERROR; case BQPDStatus::UNDEFINED: - DEBUG << YELLOW << "BQPD error: undefined\n" << RESET; + DEBUG << "BQPD error: undefined\n"; return SubproblemStatus::ERROR; } throw std::invalid_argument("The BQPD ifail is not consistent with the Uno status values"); diff --git a/uno/solvers/HiGHS/HiGHSSolver.cpp b/uno/solvers/HiGHS/HiGHSSolver.cpp index a6c769c0..e0d29943 100644 --- a/uno/solvers/HiGHS/HiGHSSolver.cpp +++ b/uno/solvers/HiGHS/HiGHSSolver.cpp @@ -3,7 +3,7 @@ #include "linear_algebra/RectangularMatrix.hpp" #include "linear_algebra/SparseVector.hpp" #include "linear_algebra/Vector.hpp" -#include "ingredients/subproblems/Direction.hpp" +#include "optimization/Direction.hpp" #include "options/Options.hpp" #include "symbolic/VectorView.hpp" diff --git a/uno/solvers/LPSolverFactory.cpp b/uno/solvers/LPSolverFactory.cpp new file mode 100644 index 00000000..528b85a6 --- /dev/null +++ b/uno/solvers/LPSolverFactory.cpp @@ -0,0 +1,56 @@ +// Copyright (c) 2018-2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include +#include "LPSolverFactory.hpp" +#include "linear_algebra/Vector.hpp" +#include "options/Options.hpp" +#include "solvers/LPSolver.hpp" + +#ifdef HAS_BQPD +#include "solvers/BQPD/BQPDSolver.hpp" +#endif +#ifdef HAS_HIGHS +#include "solvers/HiGHS/HiGHSSolver.hpp" +#endif + +namespace uno { + std::unique_ptr LPSolverFactory::create([[maybe_unused]] size_t number_variables, [[maybe_unused]] size_t number_constraints, + [[maybe_unused]] size_t number_objective_gradient_nonzeros, [[maybe_unused]] size_t number_jacobian_nonzeros, + [[maybe_unused]] const Options& options) { + try { + [[maybe_unused]] const std::string& LP_solver_name = options.get_string("LP_solver"); +#ifdef HAS_BQPD + if (LP_solver_name == "BQPD") { + return std::make_unique(number_variables, number_constraints, number_objective_gradient_nonzeros, + number_jacobian_nonzeros, 0, BQPDProblemType::LP, options); + } +#endif +#ifdef HAS_HIGHS + if (LP_solver_name == "HiGHS") { + return std::make_unique(number_variables, number_constraints, number_jacobian_nonzeros, 0, options); + } +#endif + std::string message = "The LP solver "; + message.append(LP_solver_name).append(" is unknown").append("\n").append("The following values are available: ") + .append(join(LPSolverFactory::available_solvers(), ", ")); + throw std::invalid_argument(message); + } + catch (const std::out_of_range& exception) { + std::string message = exception.what(); + message.append("\n").append("The following values are available: ").append(join(LPSolverFactory::available_solvers(), ", ")); + throw std::out_of_range(message); + } + } + + std::vector LPSolverFactory::available_solvers() { + std::vector solvers{}; +#ifdef HAS_BQPD + solvers.emplace_back("BQPD"); +#endif +#ifdef HAS_HIGHS + solvers.emplace_back("HiGHS"); +#endif + return solvers; + } +} // namespace \ No newline at end of file diff --git a/uno/solvers/LPSolverFactory.hpp b/uno/solvers/LPSolverFactory.hpp index 8613ec27..df0aa144 100644 --- a/uno/solvers/LPSolverFactory.hpp +++ b/uno/solvers/LPSolverFactory.hpp @@ -5,60 +5,20 @@ #define UNO_LPSOLVERFACTORY_H #include -#include -#include "linear_algebra/Vector.hpp" -#include "options/Options.hpp" -#include "solvers/LPSolver.hpp" - -#ifdef HAS_BQPD -#include "solvers/BQPD/BQPDSolver.hpp" -#endif -#ifdef HAS_HIGHS -#include "solvers/HiGHS/HiGHSSolver.hpp" -#endif +#include namespace uno { + // forward declarations + class LPSolver; + class Options; + class LPSolverFactory { public: static std::unique_ptr create([[maybe_unused]] size_t number_variables, [[maybe_unused]] size_t number_constraints, [[maybe_unused]] size_t number_objective_gradient_nonzeros, [[maybe_unused]] size_t number_jacobian_nonzeros, - [[maybe_unused]] const Options& options) { - try { - [[maybe_unused]] const std::string& LP_solver_name = options.get_string("LP_solver"); -#ifdef HAS_BQPD - if (LP_solver_name == "BQPD") { - return std::make_unique(number_variables, number_constraints, number_objective_gradient_nonzeros, - number_jacobian_nonzeros, 0, BQPDProblemType::LP, options); - } -#endif -#ifdef HAS_HIGHS - if (LP_solver_name == "HiGHS") { - return std::make_unique(number_variables, number_constraints, number_jacobian_nonzeros, 0, options); - } -#endif - std::string message = "The LP solver "; - message.append(LP_solver_name).append(" is unknown").append("\n").append("The following values are available: ") - .append(join(LPSolverFactory::available_solvers(), ", ")); - throw std::invalid_argument(message); - } - catch (const std::out_of_range& exception) { - std::string message = exception.what(); - message.append("\n").append("The following values are available: ").append(join(LPSolverFactory::available_solvers(), ", ")); - throw std::out_of_range(message); - } - } + [[maybe_unused]] const Options& options); - // return the list of available LP solvers - static std::vector available_solvers() { - std::vector solvers{}; -#ifdef HAS_BQPD - solvers.emplace_back("BQPD"); -#endif -#ifdef HAS_HIGHS - solvers.emplace_back("HiGHS"); -#endif - return solvers; - } + static std::vector available_solvers(); }; } // namespace diff --git a/uno/solvers/QPSolver.hpp b/uno/solvers/QPSolver.hpp index 80c9132b..fb1e93c0 100644 --- a/uno/solvers/QPSolver.hpp +++ b/uno/solvers/QPSolver.hpp @@ -24,7 +24,7 @@ namespace uno { */ class QPSolver : public LPSolver { public: - QPSolver(); + QPSolver(): LPSolver() { } ~QPSolver() override = default; void solve_LP(size_t number_variables, size_t number_constraints, const std::vector& variables_lower_bounds, @@ -39,9 +39,6 @@ namespace uno { const RectangularMatrix& constraint_jacobian, const SymmetricMatrix& hessian, const Vector& initial_point, Direction& direction, const WarmstartInformation& warmstart_information) = 0; }; - - inline QPSolver::QPSolver(): LPSolver() { - } } // namespace #endif // UNO_QPSOLVER_H \ No newline at end of file diff --git a/uno/solvers/QPSolverFactory.cpp b/uno/solvers/QPSolverFactory.cpp new file mode 100644 index 00000000..ea71c71c --- /dev/null +++ b/uno/solvers/QPSolverFactory.cpp @@ -0,0 +1,46 @@ +// Copyright (c) 2018-2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include +#include "QPSolverFactory.hpp" +#include "linear_algebra/Vector.hpp" +#include "options/Options.hpp" +#include "solvers/QPSolver.hpp" + +#ifdef HAS_BQPD +#include "solvers/BQPD/BQPDSolver.hpp" +#endif + +namespace uno { + std::unique_ptr QPSolverFactory::create([[maybe_unused]] size_t number_variables, [[maybe_unused]] size_t number_constraints, + [[maybe_unused]] size_t number_objective_gradient_nonzeros, [[maybe_unused]] size_t number_jacobian_nonzeros, + [[maybe_unused]] size_t number_hessian_nonzeros, [[maybe_unused]] const Options& options) { + try { + [[maybe_unused]] const std::string& QP_solver_name = options.get_string("QP_solver"); +#ifdef HAS_BQPD + if (QP_solver_name == "BQPD") { + return std::make_unique(number_variables, number_constraints, number_objective_gradient_nonzeros, number_jacobian_nonzeros, + number_hessian_nonzeros, BQPDProblemType::QP, options); + } +#endif + std::string message = "The QP solver "; + message.append(QP_solver_name).append(" is unknown").append("\n").append("The following values are available: ") + .append(join(QPSolverFactory::available_solvers(), ", ")); + throw std::invalid_argument(message); + } + catch (const std::out_of_range& exception) { + std::string message = exception.what(); + message.append("\n").append("The following values are available: ").append(join(QPSolverFactory::available_solvers(), ", ")); + throw std::out_of_range(message); + } + } + + // return the list of available QP solvers + std::vector QPSolverFactory::available_solvers() { + std::vector solvers{}; +#ifdef HAS_BQPD + solvers.emplace_back("BQPD"); +#endif + return solvers; + } +} // namespace diff --git a/uno/solvers/QPSolverFactory.hpp b/uno/solvers/QPSolverFactory.hpp index 64d4e7c6..7d141369 100644 --- a/uno/solvers/QPSolverFactory.hpp +++ b/uno/solvers/QPSolverFactory.hpp @@ -5,50 +5,22 @@ #define UNO_QPSOLVERFACTORY_H #include -#include -#include "linear_algebra/Vector.hpp" -#include "options/Options.hpp" -#include "solvers/QPSolver.hpp" - -#ifdef HAS_BQPD -#include "solvers/BQPD/BQPDSolver.hpp" -#endif +#include namespace uno { + // forward declarations + class Options; + class QPSolver; + class QPSolverFactory { public: // create a QP solver static std::unique_ptr create([[maybe_unused]] size_t number_variables, [[maybe_unused]] size_t number_constraints, [[maybe_unused]] size_t number_objective_gradient_nonzeros, [[maybe_unused]] size_t number_jacobian_nonzeros, - [[maybe_unused]] size_t number_hessian_nonzeros, [[maybe_unused]] const Options& options) { - try { - [[maybe_unused]] const std::string& QP_solver_name = options.get_string("QP_solver"); -#ifdef HAS_BQPD - if (QP_solver_name == "BQPD") { - return std::make_unique(number_variables, number_constraints, number_objective_gradient_nonzeros, number_jacobian_nonzeros, - number_hessian_nonzeros, BQPDProblemType::QP, options); - } -#endif - std::string message = "The QP solver "; - message.append(QP_solver_name).append(" is unknown").append("\n").append("The following values are available: ") - .append(join(QPSolverFactory::available_solvers(), ", ")); - throw std::invalid_argument(message); - } - catch (const std::out_of_range& exception) { - std::string message = exception.what(); - message.append("\n").append("The following values are available: ").append(join(QPSolverFactory::available_solvers(), ", ")); - throw std::out_of_range(message); - } - } + [[maybe_unused]] size_t number_hessian_nonzeros, [[maybe_unused]] const Options& options); // return the list of available QP solvers - static std::vector available_solvers() { - std::vector solvers{}; -#ifdef HAS_BQPD - solvers.emplace_back("BQPD"); -#endif - return solvers; - } + static std::vector available_solvers(); }; } // namespace diff --git a/uno/solvers/SymmetricIndefiniteLinearSolverFactory.cpp b/uno/solvers/SymmetricIndefiniteLinearSolverFactory.cpp new file mode 100644 index 00000000..edf79213 --- /dev/null +++ b/uno/solvers/SymmetricIndefiniteLinearSolverFactory.cpp @@ -0,0 +1,74 @@ +// Copyright (c) 2018-2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include +#include +#include "SymmetricIndefiniteLinearSolverFactory.hpp" +#include "DirectSymmetricIndefiniteLinearSolver.hpp" +#include "linear_algebra/Vector.hpp" +#include "options/Options.hpp" + +#if defined(HAS_HSL) || defined(HAS_MA57) +#include "solvers/MA57/MA57Solver.hpp" +#endif + +#ifdef HAS_HSL +namespace uno { + extern "C" { + bool LIBHSL_isfunctional(); + } +} +#endif + +#ifdef HAS_MUMPS +#include "solvers/MUMPS/MUMPSSolver.hpp" +#endif + +namespace uno { + std::unique_ptr> SymmetricIndefiniteLinearSolverFactory::create([[maybe_unused]] size_t dimension, + [[maybe_unused]] size_t number_nonzeros, const Options& options) { + try { + [[maybe_unused]] const std::string& linear_solver_name = options.get_string("linear_solver"); +#if defined(HAS_HSL) || defined(HAS_MA57) + if (linear_solver_name == "MA57" +#ifdef HAS_HSL + && LIBHSL_isfunctional() +#endif + ) { + return std::make_unique(dimension, number_nonzeros); + } +#endif +#ifdef HAS_MUMPS + if (linear_solver_name == "MUMPS") { + return std::make_unique(dimension, number_nonzeros); + } +#endif + std::string message = "The linear solver "; + message.append(linear_solver_name).append(" is unknown").append("\n").append("The following values are available: ") + .append(join(SymmetricIndefiniteLinearSolverFactory::available_solvers(), ", ")); + throw std::invalid_argument(message); + } + catch (const std::out_of_range& exception) { + std::string message = exception.what(); + message.append("\n").append("The following values are available: ") + .append(join(SymmetricIndefiniteLinearSolverFactory::available_solvers(), ", ")); + throw std::out_of_range(message); + } + } + + // return the list of available solvers + std::vector SymmetricIndefiniteLinearSolverFactory::available_solvers() { + std::vector solvers{}; +#ifdef HAS_HSL + if (LIBHSL_isfunctional()) { + solvers.emplace_back("MA57"); + } +#elif defined(HAS_MA57) + solvers.emplace_back("MA57"); +#endif +#ifdef HAS_MUMPS + solvers.emplace_back("MUMPS"); +#endif + return solvers; + } +} // namespace \ No newline at end of file diff --git a/uno/solvers/SymmetricIndefiniteLinearSolverFactory.hpp b/uno/solvers/SymmetricIndefiniteLinearSolverFactory.hpp index 0b285ed1..17542609 100644 --- a/uno/solvers/SymmetricIndefiniteLinearSolverFactory.hpp +++ b/uno/solvers/SymmetricIndefiniteLinearSolverFactory.hpp @@ -5,78 +5,21 @@ #define UNO_LINEARSOLVERFACTORY_H #include -#include -#include -#include "DirectSymmetricIndefiniteLinearSolver.hpp" -#include "linear_algebra/Vector.hpp" -#include "options/Options.hpp" -#include "tools/Logger.hpp" +#include -#if defined(HAS_MA57) || defined(HAS_HSL) -#include "solvers/MA57/MA57Solver.hpp" -#endif - -#ifdef HAS_HSL namespace uno { - extern "C" { - bool LIBHSL_isfunctional(); - } -} -#endif - -#ifdef HAS_MUMPS -#include "solvers/MUMPS/MUMPSSolver.hpp" -#endif + // forward declarations + template + class DirectSymmetricIndefiniteLinearSolver; + class Options; -namespace uno { class SymmetricIndefiniteLinearSolverFactory { public: static std::unique_ptr> create([[maybe_unused]] size_t dimension, - [[maybe_unused]] size_t number_nonzeros, const Options& options) { - try { - [[maybe_unused]] const std::string& linear_solver_name = options.get_string("linear_solver"); -#if defined(HAS_HSL) - if (linear_solver_name == "MA57" && LIBHSL_isfunctional()) { - return std::make_unique(dimension, number_nonzeros); - } -#elif defined(HAS_MA57) - if (linear_solver_name == "MA57") { - return std::make_unique(dimension, number_nonzeros); - } -#endif -#ifdef HAS_MUMPS - if (linear_solver_name == "MUMPS") { - return std::make_unique(dimension, number_nonzeros); - } -#endif - std::string message = "The linear solver "; - message.append(linear_solver_name).append(" is unknown").append("\n").append("The following values are available: ") - .append(join(SymmetricIndefiniteLinearSolverFactory::available_solvers(), ", ")); - throw std::invalid_argument(message); - } - catch (const std::out_of_range& exception) { - std::string message = exception.what(); - message.append("\n").append("The following values are available: ") - .append(join(SymmetricIndefiniteLinearSolverFactory::available_solvers(), ", ")); - throw std::out_of_range(message); - } - } + [[maybe_unused]] size_t number_nonzeros, const Options& options); // return the list of available solvers - static std::vector available_solvers() { - std::vector solvers{}; -#ifdef HAS_HSL - if (LIBHSL_isfunctional()) { - solvers.emplace_back("MA57"); - } -#elif defined(HAS_MA57) - solvers.emplace_back("MA57"); -#endif -#ifdef HAS_MUMPS - solvers.emplace_back("MUMPS"); -#endif - return solvers; - } + static std::vector available_solvers(); }; } // namespace diff --git a/uno/tools/Logger.cpp b/uno/tools/Logger.cpp new file mode 100644 index 00000000..307e7929 --- /dev/null +++ b/uno/tools/Logger.cpp @@ -0,0 +1,33 @@ +// Copyright (c) 2018-2024 Charlie Vanaret +// Licensed under the MIT license. See LICENSE file in the project directory for details. + +#include "Logger.hpp" + +namespace uno { + void Logger::set_logger(const std::string& logger_level) { + if (logger_level == "SILENT") { + Logger::level = SILENT; + } + else if (logger_level == "DISCRETE") { + Logger::level = DISCRETE; + } + else if (logger_level == "WARNING") { + Logger::level = WARNING; + } + else if (logger_level == "INFO") { + Logger::level = INFO; + } + else if (logger_level == "DEBUG") { + Logger::level = DEBUG; + } + else if (logger_level == "DEBUG2") { + Logger::level = DEBUG2; + } + else if (logger_level == "DEBUG3") { + Logger::level = DEBUG3; + } + else { + throw std::out_of_range("The logger level " + logger_level + " was not found"); + } + } +} // namespace \ No newline at end of file diff --git a/uno/tools/Logger.hpp b/uno/tools/Logger.hpp index 00d34f2f..bc473694 100644 --- a/uno/tools/Logger.hpp +++ b/uno/tools/Logger.hpp @@ -4,18 +4,10 @@ #ifndef UNO_LOGGER_H #define UNO_LOGGER_H +#include #include namespace uno { - #define RED "\x1B[31m" - // #define GREEN "\x1B[32m" - #define YELLOW "\x1B[33m" - // #define BLUE "\x1B[34m" - // #define MAGENTA "\x1B[35m" - // #define CYAN "\x1B[36m" - // #define WHITE "\x1B[37m" - #define RESET "\x1B[0m" - enum Level { SILENT = 0, DISCRETE, WARNING, INFO, DEBUG, DEBUG2, DEBUG3 }; @@ -28,45 +20,18 @@ namespace uno { template const Level& operator<<(const Level& level, T& element) { - if (level <= Logger::level) { - std::cout << element; - } - return level; + if (level <= Logger::level) { + std::cout << element; + } + return level; } template const Level& operator<<(const Level& level, const T& element) { - if (level <= Logger::level) { - std::cout << element; - } - return level; - } - - inline void Logger::set_logger(const std::string& logger_level) { - if (logger_level == "SILENT") { - Logger::level = SILENT; - } - else if (logger_level == "DISCRETE") { - Logger::level = DISCRETE; - } - else if (logger_level == "WARNING") { - Logger::level = WARNING; - } - else if (logger_level == "INFO") { - Logger::level = INFO; - } - else if (logger_level == "DEBUG") { - Logger::level = DEBUG; - } - else if (logger_level == "DEBUG2") { - Logger::level = DEBUG2; - } - else if (logger_level == "DEBUG3") { - Logger::level = DEBUG3; - } - else { - throw std::out_of_range("The logger level " + logger_level + " was not found"); + if (level <= Logger::level) { + std::cout << element; } + return level; } } // namespace diff --git a/unotest/HiGHSSolverTests.cpp b/unotest/HiGHSSolverTests.cpp index d21304d5..d61a7213 100644 --- a/unotest/HiGHSSolverTests.cpp +++ b/unotest/HiGHSSolverTests.cpp @@ -2,7 +2,7 @@ // Licensed under the MIT license. See LICENSE file in the project directory for details. #include -#include "ingredients/subproblems/Direction.hpp" +#include "optimization/Direction.hpp" #include "linear_algebra/RectangularMatrix.hpp" #include "linear_algebra/SparseVector.hpp" #include "optimization/WarmstartInformation.hpp"