Store the fixed variables in the model. If there are fixed variables,…

… fail in PrimalDualInteriorPointSubproblem for the moment. TO-DO: move the fixed bounds to the set of general constraints

bindings/AMPL/AMPLModel.cpp

@@ -13,7 +13,6 @@
 #include "tools/Infinity.hpp"
 #include "options/Options.hpp"
 #include "symbolic/Concatenation.hpp"
-#include "symbolic/ScalarMultiple.hpp"
 #include "symbolic/UnaryNegation.hpp"
 #include "Uno.hpp"
 
@@ -27,7 +26,6 @@ namespace uno {
       int n_discrete = asl->i.nbv_ + asl->i.niv_ + asl->i.nlvbi_ + asl->i.nlvci_ + asl->i.nlvoi_;
       if (0 < n_discrete) {
          throw std::runtime_error("Error: " + std::to_string(n_discrete) + " variables are discrete, which Uno cannot handle");
-         // asl->i.need_nl_ = 0;
       }
 
       // preallocate initial primal and dual solutions
@@ -56,13 +54,14 @@ namespace uno {
          variable_status(this->number_variables),
          constraint_type(this->number_constraints),
          constraint_status(this->number_constraints),
-         multipliers_with_flipped_sign(this->number_variables, this->number_constraints),
+         multipliers_with_flipped_sign(this->number_constraints),
          equality_constraints_collection(this->equality_constraints),
          inequality_constraints_collection(this->inequality_constraints),
          lower_bounded_variables_collection(this->lower_bounded_variables),
          upper_bounded_variables_collection(this->upper_bounded_variables),
          single_lower_bounded_variables_collection(this->single_lower_bounded_variables),
-         single_upper_bounded_variables_collection(this->single_upper_bounded_variables) {
+         single_upper_bounded_variables_collection(this->single_upper_bounded_variables),
+         fixed_variables_collection(this->fixed_variables) {
       // evaluate the constraint Jacobian in sparse mode
       this->asl->i.congrd_mode = 1;
 
@@ -71,6 +70,7 @@ namespace uno {
       this->upper_bounded_variables.reserve(this->number_variables);
       this->single_lower_bounded_variables.reserve(this->number_variables);
       this->single_upper_bounded_variables.reserve(this->number_variables);
+      this->fixed_variables.reserve(this->number_variables);
       this->generate_variables();
 
       // constraints
@@ -93,6 +93,7 @@ namespace uno {
          this->variable_upper_bounds[variable_index] = (this->asl->i.LUv_ != nullptr) ? this->asl->i.LUv_[2*variable_index + 1] : INF<double>;
          if (this->variable_lower_bounds[variable_index] == this->variable_upper_bounds[variable_index]) {
             WARNING << "Variable x" << variable_index << " has identical bounds\n";
+            this->fixed_variables.emplace_back(variable_index);
          }
       }
       AMPLModel::determine_bounds_types(this->variable_lower_bounds, this->variable_upper_bounds, this->variable_status);
@@ -283,16 +284,16 @@ namespace uno {
       // evaluate the Hessian: store the matrix in a preallocated array this->asl_hessian
       const int objective_number = -1;
       // flip the signs of the multipliers: in AMPL, the Lagrangian is f + lambda.g, while Uno uses f - lambda.g
-      this->multipliers_with_flipped_sign.constraints = -multipliers;
+      this->multipliers_with_flipped_sign = -multipliers;
       if (this->fixed_hessian_sparsity) {
          (*(this->asl)->p.Sphes)(this->asl, nullptr, const_cast<double*>(this->asl_hessian.data()), objective_number, &objective_multiplier,
-               const_cast<double*>(this->multipliers_with_flipped_sign.constraints.data()));
+               const_cast<double*>(this->multipliers_with_flipped_sign.data()));
       }
       else {
          double* objective_multiplier_pointer = (objective_multiplier != 0.) ? &objective_multiplier : nullptr;
          const bool all_zeros_multipliers = are_all_zeros(multipliers);
          (*(this->asl)->p.Sphes)(this->asl, nullptr, const_cast<double*>(this->asl_hessian.data()), objective_number, objective_multiplier_pointer,
-               all_zeros_multipliers ? nullptr : const_cast<double*>(this->multipliers_with_flipped_sign.constraints.data()));
+               all_zeros_multipliers ? nullptr : const_cast<double*>(this->multipliers_with_flipped_sign.data()));
       }
 
       // generate the sparsity pattern in the right sparse format
@@ -376,24 +377,25 @@ namespace uno {
             this->asl->p.solve_code_ = 500;
          }
 
-         // flip the signs of the multipliers if we maximize
-         this->multipliers_with_flipped_sign.constraints = this->objective_sign * iterate.multipliers.constraints;
-         this->multipliers_with_flipped_sign.lower_bounds = this->objective_sign * iterate.multipliers.lower_bounds;
-         this->multipliers_with_flipped_sign.upper_bounds = this->objective_sign * iterate.multipliers.upper_bounds;
+         // flip the signs of the multipliers and the objective if we maximize
+         iterate.multipliers.constraints *= this->objective_sign;
+         iterate.multipliers.lower_bounds *= this->objective_sign;
+         iterate.multipliers.upper_bounds *= this->objective_sign;
+         iterate.evaluations.objective *= this->objective_sign;
 
          // include the bound duals in the .sol file, using suffixes
          SufDecl lower_bound_suffix{const_cast<char*>("lower_bound_duals"), nullptr, ASL_Sufkind_var | ASL_Sufkind_real, 0};
          SufDecl upper_bound_suffix{const_cast<char*>("upper_bound_duals"), nullptr, ASL_Sufkind_var | ASL_Sufkind_real, 0};
          std::array<SufDecl, 2> suffixes{lower_bound_suffix, upper_bound_suffix};
-         suf_declare_ASL(this->asl, suffixes.data(), 2);
-         suf_rput_ASL(this->asl, "lower_bound_duals", ASL_Sufkind_var, this->multipliers_with_flipped_sign.lower_bounds.data());
-         suf_rput_ASL(this->asl, "upper_bound_duals", ASL_Sufkind_var, this->multipliers_with_flipped_sign.upper_bounds.data());
+         suf_declare_ASL(this->asl, suffixes.data(), suffixes.size());
+         suf_rput_ASL(this->asl, "lower_bound_duals", ASL_Sufkind_var, iterate.multipliers.lower_bounds.data());
+         suf_rput_ASL(this->asl, "upper_bound_duals", ASL_Sufkind_var, iterate.multipliers.upper_bounds.data());
 
          Option_Info option_info{};
          option_info.wantsol = 9; // write the solution without printing the message to stdout
          std::string message = "Uno ";
          message.append(Uno::current_version()).append(": ").append(status_to_message(termination_status));
-         write_sol_ASL(this->asl, message.data(), iterate.primals.data(), this->multipliers_with_flipped_sign.constraints.data(), &option_info);
+         write_sol_ASL(this->asl, message.data(), iterate.primals.data(), iterate.multipliers.constraints.data(), &option_info);
       }
    }
 

bindings/AMPL/AMPLModel.hpp

@@ -47,6 +47,7 @@ namespace uno {
       [[nodiscard]] const SparseVector<size_t>& get_slacks() const override;
       [[nodiscard]] const Collection<size_t>& get_single_lower_bounded_variables() const override;
       [[nodiscard]] const Collection<size_t>& get_single_upper_bounded_variables() const override;
+      [[nodiscard]] const Collection<size_t>& get_fixed_variables() const override { return this->fixed_variables_collection; }
 
       [[nodiscard]] double constraint_lower_bound(size_t constraint_index) const override;
       [[nodiscard]] double constraint_upper_bound(size_t constraint_index) const override;
@@ -84,7 +85,7 @@ namespace uno {
       std::vector<BoundType> constraint_status; /*!< Status of the constraints (EQUAL_BOUNDS, BOUNDED_LOWER, BOUNDED_UPPER, BOUNDED_BOTH_SIDES,
     * UNBOUNDED) */
       std::vector<size_t> linear_constraints;
-      mutable Multipliers multipliers_with_flipped_sign;
+      mutable Vector<double> multipliers_with_flipped_sign;
 
       // lists of variables and constraints + corresponding collection objects
       std::vector<size_t> equality_constraints{};
@@ -100,6 +101,8 @@ namespace uno {
       CollectionAdapter<std::vector<size_t>&> single_lower_bounded_variables_collection;
       std::vector<size_t> single_upper_bounded_variables{}; // indices of the single upper-bounded variables
       CollectionAdapter<std::vector<size_t>&> single_upper_bounded_variables_collection;
+      std::vector<size_t> fixed_variables;
+      CollectionAdapter<std::vector<size_t>&> fixed_variables_collection;
 
       void generate_variables();
       void generate_constraints();

bindings/AMPL/uno_ampl.cpp

@@ -50,8 +50,7 @@ namespace uno {
          Uno uno = Uno(*globalization_mechanism, options);
 
          // solve the instance
-         Result result = uno.solve(*model, initial_iterate, options);
-         uno.print_optimization_summary(result);
+         uno.solve(*model, initial_iterate, options);
          // std::cout << "memory_allocation_amount = " << memory_allocation_amount << '\n';
       }
       catch (std::exception& exception) {

uno/Uno.cpp

@@ -29,58 +29,58 @@ namespace uno {
 
    Level Logger::level = INFO;
 
-   Result Uno::solve(const Model& model, Iterate& current_iterate, const Options& options) {
+   void Uno::solve(const Model& model, Iterate& current_iterate, const Options& options) {
       Timer timer{};
       Statistics statistics = Uno::create_statistics(model, options);
 
-      // use the initial primal-dual point to initialize the strategies and generate the initial iterate
-      this->initialize(statistics, current_iterate, options);
-      current_iterate.status = TerminationStatus::NOT_OPTIMAL;
-      // allocate the trial iterate once and for all here
-      Iterate trial_iterate(current_iterate);
-
-      size_t major_iterations = 0;
       try {
-         bool termination = false;
-         // check for termination
-         while (not termination) {
-            major_iterations++;
+         // use the initial primal-dual point to initialize the strategies and generate the initial iterate
+         this->initialize(statistics, current_iterate, options);
+         current_iterate.status = TerminationStatus::NOT_OPTIMAL;
+         // allocate the trial iterate once and for all here
+         Iterate trial_iterate(current_iterate);
+
+         size_t major_iterations = 0;
+         try {
+            bool termination = false;
+            // check for termination
+            while (not termination) {
+               major_iterations++;
+               statistics.start_new_line();
+               statistics.set("iter", major_iterations);
+               DEBUG << "### Outer iteration " << major_iterations << '\n';
+
+               // compute an acceptable iterate by solving a subproblem at the current point
+               this->globalization_mechanism.compute_next_iterate(statistics, model, current_iterate, trial_iterate);
+               termination = this->termination_criteria(trial_iterate.status, major_iterations, timer.get_duration());
+               // the trial iterate becomes the current iterate for the next iteration
+               std::swap(current_iterate, trial_iterate);
+            }
+         }
+         catch (std::exception& exception) {
             statistics.start_new_line();
-            statistics.set("iter", major_iterations);
-            DEBUG << "### Outer iteration " << major_iterations << '\n';
-
-            // compute an acceptable iterate by solving a subproblem at the current point
-            this->globalization_mechanism.compute_next_iterate(statistics, model, current_iterate, trial_iterate);
-            termination = this->termination_criteria(trial_iterate.status, major_iterations, timer.get_duration());
-            // the trial iterate becomes the current iterate for the next iteration
-            std::swap(current_iterate, trial_iterate);
+            statistics.set("status", exception.what());
+            if (Logger::level == INFO) statistics.print_current_line();
+            DEBUG << exception.what() << '\n';
          }
+         if (Logger::level == INFO) statistics.print_footer();
+
+         Uno::postprocess_iterate(model, current_iterate, current_iterate.status);
+         Result result = this->create_result(model, current_iterate, major_iterations, timer);
+         this->print_optimization_summary(result);
       }
-      catch (std::exception& exception) {
-         statistics.start_new_line();
-         statistics.set("status", exception.what());
-         if (Logger::level == INFO) statistics.print_current_line();
-         DEBUG << exception.what() << '\n';
+      catch (const std::exception& e) {
+         DISCRETE  << "An error occurred at the initial iterate: " << e.what()  << '\n';
       }
-      if (Logger::level == INFO) statistics.print_footer();
-
-      Uno::postprocess_iterate(model, current_iterate, current_iterate.status);
-      return this->create_result(model, current_iterate, major_iterations, timer);
    }
 
    void Uno::initialize(Statistics& statistics, Iterate& current_iterate, const Options& options) {
-      try {
-         statistics.start_new_line();
-         statistics.set("iter", 0);
-         statistics.set("status", "initial point");
-         this->globalization_mechanism.initialize(statistics, current_iterate, options);
-         options.print_used();
-         if (Logger::level == INFO) statistics.print_current_line();
-      }
-      catch (const std::exception& e) {
-         DISCRETE << RED << "An error occurred at the initial iterate: " << e.what() << RESET << '\n';
-         throw;
-      }
+      statistics.start_new_line();
+      statistics.set("iter", 0);
+      statistics.set("status", "initial point");
+      this->globalization_mechanism.initialize(statistics, current_iterate, options);
+      options.print_used();
+      if (Logger::level == INFO) statistics.print_current_line();
    }
 
    Statistics Uno::create_statistics(const Model& model, const Options& options) {

uno/Uno.hpp

@@ -27,7 +27,7 @@ namespace uno {
    public:
       Uno(GlobalizationMechanism& globalization_mechanism, const Options& options);
 
-      [[nodiscard]] Result solve(const Model& model, Iterate& initial_iterate, const Options& options);
+      void solve(const Model& model, Iterate& initial_iterate, const Options& options);
 
       static std::string current_version();
       static void print_available_strategies();

uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp

@@ -58,6 +58,9 @@ namespace uno {
       if (problem.has_inequality_constraints()) {
          throw std::runtime_error("The problem has inequality constraints. Create an instance of HomogeneousEqualityConstrainedModel");
       }
+      if (problem.has_fixed_variables()) {
+         throw std::runtime_error("The problem has fixed variables. Move them to the set of general constraints.");
+      }
 
       // TODO: enforce linear constraints at initial point
       //if (options.get_bool("enforce_linear_constraints")) {

uno/linear_algebra/Vector.hpp

@@ -89,6 +89,10 @@ namespace uno {
          }
       }
 
+      void operator*=(ElementType factor) {
+         this->scale(factor);
+      }
+
       ElementType* data() { return this->vector.data(); }
       const ElementType* data() const { return this->vector.data(); }
 

uno/model/BoundRelaxedModel.hpp

@@ -40,6 +40,7 @@ namespace uno {
       [[nodiscard]] const SparseVector<size_t>& get_slacks() const override { return this->model->get_slacks(); }
       [[nodiscard]] const Collection<size_t>& get_single_lower_bounded_variables() const override { return this->model->get_single_lower_bounded_variables(); }
       [[nodiscard]] const Collection<size_t>& get_single_upper_bounded_variables() const override { return this->model->get_single_upper_bounded_variables(); }
+      [[nodiscard]] const Collection<size_t>& get_fixed_variables() const override { return this->model->get_fixed_variables(); }
 
       [[nodiscard]] double constraint_lower_bound(size_t constraint_index) const override { return this->model->constraint_lower_bound(constraint_index); }
       [[nodiscard]] double constraint_upper_bound(size_t constraint_index) const override { return this->model->constraint_upper_bound(constraint_index); }

uno/model/HomogeneousEqualityConstrainedModel.hpp

@@ -45,6 +45,7 @@ namespace uno {
       [[nodiscard]] const Collection<size_t>& get_equality_constraints() const override { return this->equality_constraints; }
       [[nodiscard]] const Collection<size_t>& get_inequality_constraints() const override { return this->inequality_constraints; }
       [[nodiscard]] const std::vector<size_t>& get_linear_constraints() const override { return this->model->get_linear_constraints(); }
+      [[nodiscard]] const Collection<size_t>& get_fixed_variables() const override { return this->model->get_fixed_variables(); }
 
       void initial_primal_point(Vector<double>& x) const override;
       void initial_dual_point(Vector<double>& multipliers) const override { this->model->initial_dual_point(multipliers); }

uno/model/Model.hpp

@@ -64,6 +64,7 @@ namespace uno {
       [[nodiscard]] virtual const SparseVector<size_t>& get_slacks() const = 0;
       [[nodiscard]] virtual const Collection<size_t>& get_single_lower_bounded_variables() const = 0;
       [[nodiscard]] virtual const Collection<size_t>& get_single_upper_bounded_variables() const = 0;
+      [[nodiscard]] virtual const Collection<size_t>& get_fixed_variables() const = 0;
 
       [[nodiscard]] virtual double constraint_lower_bound(size_t constraint_index) const = 0;
       [[nodiscard]] virtual double constraint_upper_bound(size_t constraint_index) const = 0;

uno/model/ScaledModel.hpp

@@ -30,6 +30,7 @@ namespace uno {
       [[nodiscard]] const SparseVector<size_t>& get_slacks() const override { return this->model->get_slacks(); }
       [[nodiscard]] const Collection<size_t>& get_single_lower_bounded_variables() const override { return this->model->get_single_lower_bounded_variables(); }
       [[nodiscard]] const Collection<size_t>& get_single_upper_bounded_variables() const override { return this->model->get_single_upper_bounded_variables(); }
+      [[nodiscard]] const Collection<size_t>& get_fixed_variables() const override { return this->model->get_fixed_variables(); }
 
       [[nodiscard]] double constraint_lower_bound(size_t constraint_index) const override;
       [[nodiscard]] double constraint_upper_bound(size_t constraint_index) const override;

uno/reformulation/OptimizationProblem.hpp

@@ -34,6 +34,7 @@ namespace uno {
 
       [[nodiscard]] bool is_constrained() const;
       [[nodiscard]] bool has_inequality_constraints() const;
+      [[nodiscard]] bool has_fixed_variables() const;
 
       // function evaluations
       [[nodiscard]] virtual double get_objective_multiplier() const = 0;
@@ -75,6 +76,10 @@ namespace uno {
       return (not this->model.get_inequality_constraints().is_empty());
    }
 
+   inline bool OptimizationProblem::has_fixed_variables() const {
+      return (not this->model.get_fixed_variables().is_empty());
+   }
+
    inline size_t OptimizationProblem::get_number_original_variables() const {
       return this->model.number_variables;
    }