[WIP] Refactoring of the reformulation chain

.github/julia/runtests_uno_ipopt_ssids.jl

@@ -22,7 +22,7 @@ function Optimizer(options)
     return AmplNLWriter.Optimizer(Uno_jll.amplexe, options)
 end
 
-Optimizer_Uno_ipopt() = Optimizer(["logger=INFO", "preset=ipopt", "linear_solver=SSIDS", "unbounded_objective_threshold=-1e15"])
+Optimizer_Uno_ipopt() = Optimizer(["logger=SILENT", "preset=ipopt", "linear_solver=SSIDS", "unbounded_objective_threshold=-1e15"])
 
 # This testset runs https://github.com/jump-dev/MINLPTests.jl
 @testset "MINLPTests" begin

interfaces/Fortran/example_uno.f90

@@ -125,6 +125,8 @@ program example_uno
     success = uno_set_lagrangian_hessian(model, number_hessian_nonzeros, hessian_triangular_part, hessian_row_indices, &
                                          hessian_column_indices, lagrangian_hessian)
     success = uno_set_lagrangian_sign_convention(model, lagrangian_sign_convention)
+    ! The Hessian model was overwritten. Set it again
+    success = uno_set_solver_string_option(solver, "hessian_model", hessian_model)
     call uno_optimize(solver, model)
 
     !---------------------------------------------------

uno/Uno.cpp

@@ -2,7 +2,6 @@
 // Licensed under the MIT license. See LICENSE file in the project directory for details.
 
 #include "Uno.hpp"
-#include "ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.hpp"
 #include "ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategyFactory.hpp"
 #include "ingredients/globalization_mechanisms/GlobalizationMechanismFactory.hpp"
 #include "ingredients/globalization_strategies/GlobalizationStrategyFactory.hpp"
@@ -78,25 +77,24 @@ namespace uno {
       model.initial_primal_point(current_iterate.primals);
       model.initial_dual_point(current_iterate.multipliers.constraints);
       model.reset_number_evaluations();
-      OptimizationStatus optimization_status = OptimizationStatus::SUCCESS;
       EvaluationCache evaluation_cache{model};
-
       Statistics statistics = Uno::create_statistics(model);
       WarmstartInformation warmstart_information{};
       warmstart_information.whole_problem_changed();
-      const size_t max_iterations = options.get_unsigned_int("max_iterations"); // maximum number of iterations
-      const double time_limit = options.get_double("time_limit"); // CPU time limit (can be inf)
 
+      OptimizationStatus optimization_status = OptimizationStatus::SUCCESS;
       try { // catch errors at startup/initial iterate
-         // pick the ingredients based on the user-defined options
-         Uno::pick_ingredients(model, options);
+         // set the ingredients based on the user-defined options
+         this->globalization_mechanism = GlobalizationMechanismFactory::create(model, options);
          // use the initial primal-dual point to initialize the strategies and generate the initial iterate
          this->initialize(statistics, model, current_iterate, options, evaluation_cache);
          // allocate the trial iterate once and for all here
          Iterate trial_iterate(current_iterate);
 
-         try {
+         try { // catch errors during the optimization process
             bool termination = false;
+            const size_t max_iterations = options.get_unsigned_int("max_iterations"); // maximum number of iterations
+            const double time_limit = options.get_double("time_limit"); // CPU time limit (can be inf)
             if (max_iterations == 0) {
                termination = true;
                optimization_status = OptimizationStatus::ITERATION_LIMIT;
@@ -112,7 +110,6 @@ namespace uno {
                warmstart_information.iterate_changed();
                this->globalization_mechanism->compute_next_iterate(statistics, model, current_iterate, trial_iterate,
                   this->direction, evaluation_cache, warmstart_information, user_callbacks);
-               GlobalizationMechanism::set_dual_residuals_statistics(statistics, trial_iterate);
                const bool user_termination = user_callbacks.user_termination(trial_iterate.primals, trial_iterate.multipliers,
                   trial_iterate.objective_multiplier, trial_iterate.progress.infeasibility, trial_iterate.residuals.stationarity,
                   trial_iterate.residuals.complementarity);
@@ -165,27 +162,20 @@ namespace uno {
       std::cout << "- Presets: filtersqp, ipopt\n";
    }
 
-   void Uno::pick_ingredients(const Model& model, Options& options) {
-      this->globalization_mechanism = GlobalizationMechanismFactory::create(model, options);
-   }
-
    void Uno::initialize(Statistics& statistics, const Model& model, Iterate& current_iterate, const Options& options,
          EvaluationCache& evaluation_cache) {
       statistics.start_new_line();
       statistics.set("Major", 0);
       statistics.set("Status", "initial point");
 
       model.project_onto_variable_bounds(current_iterate.primals);
-      this->globalization_mechanism->initialize(statistics, current_iterate, this->direction, evaluation_cache);
-      GlobalizationMechanism::set_primal_statistics(statistics, model, current_iterate, evaluation_cache.current_evaluations);
-      GlobalizationMechanism::set_dual_residuals_statistics(statistics, current_iterate);
+      this->globalization_mechanism->initialize(statistics, model, current_iterate, this->direction, evaluation_cache, options);
 
       options.print_non_default();
       if (Logger::level == INFO) {
          statistics.print_header();
          statistics.print_current_line();
       }
-      current_iterate.status = SolutionStatus::NOT_OPTIMAL;
    }
 
    Statistics Uno::create_statistics(const Model& model) {

uno/Uno.hpp

@@ -34,7 +34,6 @@ namespace uno {
       std::unique_ptr<GlobalizationMechanism> globalization_mechanism{};
       Direction direction{};
 
-      void pick_ingredients(const Model& model, Options& options);
       void initialize(Statistics& statistics, const Model& model, Iterate& current_iterate, const Options& options,
          EvaluationCache& evaluation_cache);
       [[nodiscard]] static Statistics create_statistics(const Model& model);

uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.cpp

@@ -3,17 +3,22 @@
 
 #include "ConstraintRelaxationStrategy.hpp"
 #include "ingredients/globalization_strategies/GlobalizationStrategy.hpp"
-#include "ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp"
+#include "ingredients/subproblem/Subproblem.hpp"
 #include "linear_algebra/VectorView.hpp"
 #include "model/Model.hpp"
 #include "optimization/Direction.hpp"
+#include "optimization/EvaluationCache.hpp"
 #include "optimization/Iterate.hpp"
 #include "optimization/Multipliers.hpp"
 #include "optimization/OptimizationProblem.hpp"
 #include "options/Options.hpp"
+#include "tools/Logger.hpp"
+#include "tools/Statistics.hpp"
+#include "tools/UserCallbacks.hpp"
 
 namespace uno {
    ConstraintRelaxationStrategy::ConstraintRelaxationStrategy(const Options& options):
+         progress_norm(norm_from_string(options.get_string("progress_norm"))),
          residual_norm(norm_from_string(options.get_string("residual_norm"))),
          residual_scaling_threshold(options.get_double("residual_scaling_threshold")),
          primal_tolerance(options.get_double("primal_tolerance")),
@@ -37,6 +42,61 @@ namespace uno {
    }
    */
 
+   size_t ConstraintRelaxationStrategy::get_number_subproblems_solved() const {
+      return this->number_subproblems_solved;
+   }
+
+   // protected member functions
+
+   void ConstraintRelaxationStrategy::evaluate_progress_measures(const OptimizationProblem& problem, Iterate& iterate,
+         Evaluations& evaluations) const {
+      problem.set_infeasibility_measure(iterate, evaluations, this->progress_norm);
+      problem.set_objective_measure(iterate, evaluations);
+      problem.set_auxiliary_measure(iterate);
+   }
+
+   bool ConstraintRelaxationStrategy::is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy,
+         const Subproblem& subproblem, const SolverWorkspace& solver_workspace, Iterate& current_iterate, Iterate& trial_iterate,
+         const Direction& direction, double step_length, EvaluationCache& evaluation_cache, UserCallbacks& user_callbacks) const {
+      subproblem.problem.postprocess_iterate(trial_iterate);
+      const double objective_multiplier = subproblem.problem.get_objective_multiplier();
+
+      // evaluate progress measures
+      trial_iterate.objective_multiplier = objective_multiplier;
+      //if (this->subproblem_definition_changed) {
+         //DEBUG << "The subproblem definition changed, the globalization strategy is reset and the auxiliary measure is recomputed\n";
+         //globalization_strategy.reset();
+         subproblem.problem.set_auxiliary_measure(current_iterate);
+         //this->subproblem_definition_changed = false;
+      //}
+      this->evaluate_progress_measures(subproblem.problem, trial_iterate, evaluation_cache.trial_evaluations);
+
+      bool accept_iterate = false;
+      if (direction.norm == 0.) {
+         DEBUG << "Zero step acceptable\n";
+         evaluation_cache.trial_evaluations.objective = subproblem.problem.model.evaluate_objective(trial_iterate.primals);
+         accept_iterate = true;
+         statistics.set("Status", "0 primal step");
+      }
+      else {
+         // determine acceptance wrt the globalization strategy
+         const ProgressMeasures predicted_reductions = subproblem.compute_predicted_reductions(direction, step_length,
+            this->progress_norm, evaluation_cache.current_evaluations, solver_workspace);
+         accept_iterate = globalization_strategy.is_iterate_acceptable(statistics, current_iterate.progress, trial_iterate.progress,
+            predicted_reductions, objective_multiplier);
+         // check that the derivatives exist at the accepted trial iterate (an exception is thrown upon evaluation failure)
+         if (accept_iterate) {
+            evaluation_cache.trial_evaluations.evaluate_objective_gradient(subproblem.problem.model, trial_iterate.primals);
+            evaluation_cache.trial_evaluations.evaluate_jacobian(subproblem.problem.model, trial_iterate.primals);
+         }
+      }
+      if (accept_iterate) {
+         user_callbacks.notify_acceptable_iterate(trial_iterate.primals, trial_iterate.multipliers, objective_multiplier,
+            trial_iterate.progress.infeasibility, trial_iterate.residuals.stationarity, trial_iterate.residuals.complementarity);
+      }
+      return accept_iterate;
+   }
+
    // stationarity errors:
    // - for KKT conditions: with standard multipliers and current objective multiplier
    // - for FJ conditions: with standard multipliers and 0 objective multiplier

uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.hpp

@@ -15,11 +15,14 @@ namespace uno {
    class Direction;
    class EvaluationCache;
    class Evaluations;
+   class GlobalizationStrategy;
    class Model;
    class Multipliers;
    class OptimizationProblem;
    class Options;
+   class SolverWorkspace;
    class Statistics;
+   class Subproblem;
    class UserCallbacks;
    class WarmstartInformation;
 
@@ -29,24 +32,25 @@ namespace uno {
       virtual ~ConstraintRelaxationStrategy();
 
       virtual void initialize(Statistics& statistics, Iterate& initial_iterate, Direction& direction, bool uses_trust_region,
-         EvaluationCache& evaluation_cache) = 0;
+         EvaluationCache& evaluation_cache, const Options& options) = 0;
 
       // direction computation
       virtual void compute_feasible_direction(Statistics& statistics, Iterate& current_iterate, Direction& direction,
          double trust_region_radius, Evaluations& current_evaluations, WarmstartInformation& warmstart_information) = 0;
       [[nodiscard]] virtual bool solving_feasibility_problem() const = 0;
       virtual void switch_to_feasibility_problem(Statistics& statistics, Iterate& current_iterate, Evaluations& current_evaluations,
-         bool uses_trust_region, WarmstartInformation& warmstart_information) = 0;
+         WarmstartInformation& warmstart_information) = 0;
 
       // trial iterate acceptance
       [[nodiscard]] virtual bool is_iterate_acceptable(Statistics& statistics, const Model& model, Iterate& current_iterate,
          Iterate& trial_iterate, const Direction& direction, double step_length, EvaluationCache& evaluation_cache,
          WarmstartInformation& warmstart_information, UserCallbacks& user_callbacks) = 0;
 
       [[nodiscard]] virtual std::string get_name() const = 0;
-      [[nodiscard]] virtual size_t get_number_subproblems_solved() const = 0;
+      [[nodiscard]] size_t get_number_subproblems_solved() const;
 
    protected:
+      const Norm progress_norm;
       const Norm residual_norm;
       const double residual_scaling_threshold;
       const double primal_tolerance;
@@ -56,7 +60,12 @@ namespace uno {
       size_t loose_tolerance_consecutive_iterations{0};
       const size_t loose_tolerance_iteration_threshold;
       const double unbounded_objective_threshold;
+      size_t number_subproblems_solved{0};
 
+      void evaluate_progress_measures(const OptimizationProblem& problem, Iterate& iterate, Evaluations& evaluations) const;
+      bool is_iterate_acceptable(Statistics& statistics, GlobalizationStrategy& globalization_strategy, const Subproblem& subproblem,
+         const SolverWorkspace& solver_workspace, Iterate& current_iterate, Iterate& trial_iterate,
+         const Direction& direction, double step_length, EvaluationCache& evaluation_cache, UserCallbacks& user_callbacks) const;
       void compute_residuals(const OptimizationProblem& problem, Iterate& iterate, Evaluations& evaluations) const;
       [[nodiscard]] double compute_stationarity_scaling(const Model& model, const Multipliers& multipliers) const;
       [[nodiscard]] double compute_complementarity_scaling(const Model& model, const Multipliers& multipliers) const;