Preprocessing: fixed enforce_linear_constraints() in case the primal …

…vector is longer than the direction
cvanaret · Nov 12, 2024 · d57ee9c · d57ee9c
1 parent a58a6e5
commit d57ee9c
Show file tree

Hide file tree

Showing 4 changed files with 22 additions and 11 deletions.
diff --git a/uno/preprocessing/Preprocessing.cpp b/uno/preprocessing/Preprocessing.cpp
@@ -85,13 +85,13 @@ namespace uno {
       return infeasible_linear_constraints;
    }
 
-   void Preprocessing::enforce_linear_constraints(const Model& model, Vector<double>& x, Multipliers& multipliers, QPSolver& qp_solver) {
+   void Preprocessing::enforce_linear_constraints(const Model& model, Vector<double>& primals, Multipliers& multipliers, QPSolver& qp_solver) {
       const auto& linear_constraints = model.get_linear_constraints();
       INFO << "\nPreprocessing phase: the problem has " << linear_constraints.size() << " linear constraints\n";
       if (not linear_constraints.empty()) {
          // evaluate the constraints
          std::vector<double> constraints(model.number_constraints);
-         model.evaluate_constraints(x, constraints);
+         model.evaluate_constraints(primals, constraints);
          const size_t infeasible_linear_constraints = count_infeasible_linear_constraints(model, constraints);
          INFO << "There are " << infeasible_linear_constraints << " infeasible linear constraints at the initial point\n";
          if (0 < infeasible_linear_constraints) {
@@ -105,15 +105,15 @@ namespace uno {
             RectangularMatrix<double> constraint_jacobian(linear_constraints.size(), model.number_variables);
             size_t linear_constraint_index = 0;
             for (size_t constraint_index: linear_constraints) {
-               model.evaluate_constraint_gradient(x, constraint_index, constraint_jacobian[linear_constraint_index]);
+               model.evaluate_constraint_gradient(primals, constraint_index, constraint_jacobian[linear_constraint_index]);
                linear_constraint_index++;
             }
             // variable bounds
             std::vector<double> variables_lower_bounds(model.number_variables);
             std::vector<double> variables_upper_bounds(model.number_variables);
             for (size_t variable_index: Range(model.number_variables)) {
-               variables_lower_bounds[variable_index] = model.variable_lower_bound(variable_index) - x[variable_index];
-               variables_upper_bounds[variable_index] = model.variable_upper_bound(variable_index) - x[variable_index];
+               variables_lower_bounds[variable_index] = model.variable_lower_bound(variable_index) - primals[variable_index];
+               variables_upper_bounds[variable_index] = model.variable_upper_bound(variable_index) - primals[variable_index];
             }
             // constraint bounds
             std::vector<double> constraints_lower_bounds(linear_constraints.size());
@@ -137,15 +137,17 @@ namespace uno {
             }
 
             // take the step
-            x += direction.primals;
-            multipliers.lower_bounds += direction.multipliers.lower_bounds;
-            multipliers.upper_bounds += direction.multipliers.upper_bounds;
+            for (size_t variable_index: Range(model.number_variables)) {
+               primals[variable_index] += direction.primals[variable_index];
+               multipliers.lower_bounds[variable_index] += direction.multipliers.lower_bounds[variable_index];
+               multipliers.upper_bounds[variable_index] += direction.multipliers.upper_bounds[variable_index];
+            }
             linear_constraint_index = 0;
             for (size_t constraint_index: linear_constraints) {
                multipliers.constraints[constraint_index] += direction.multipliers.constraints[linear_constraint_index];
                linear_constraint_index++;
             }
-            DEBUG3 << "Linear feasible initial point: " << x << '\n';
+            DEBUG3 << "Linear feasible initial point: " << view(primals, 0, model.number_variables) << '\n';
          }
       }
    }

diff --git a/uno/preprocessing/Preprocessing.hpp b/uno/preprocessing/Preprocessing.hpp
@@ -25,7 +25,7 @@ namespace uno {
       static void compute_least_square_multipliers(const Model& model, SymmetricMatrix<size_t, double>& matrix, Vector<double>& rhs,
             DirectSymmetricIndefiniteLinearSolver<size_t, double>& linear_solver, Iterate& current_iterate, Vector<double>& multipliers,
             double multiplier_max_norm);
-      static void enforce_linear_constraints(const Model& model, Vector<double>& x, Multipliers& multipliers, QPSolver& qp_solver);
+      static void enforce_linear_constraints(const Model& model, Vector<double>& primals, Multipliers& multipliers, QPSolver& qp_solver);
    };
 } // namespace
 

diff --git a/uno/solvers/BQPD/BQPDSolver.cpp b/uno/solvers/BQPD/BQPDSolver.cpp
@@ -169,7 +169,7 @@ namespace uno {
       if (warmstart_information.problem_changed) {
          mode = BQPDMode::ACTIVE_SET_EQUALITIES;
       }
-         // if only the variable bounds changed, reuse the active set estimate and the Jacobian information
+      // if only the variable bounds changed, reuse the active set estimate and the Jacobian information
       else if (warmstart_information.variable_bounds_changed && not warmstart_information.objective_changed &&
                not warmstart_information.constraints_changed && not warmstart_information.constraint_bounds_changed) {
          mode = BQPDMode::UNCHANGED_ACTIVE_SET_AND_JACOBIAN;

diff --git a/uno/symbolic/VectorView.hpp b/uno/symbolic/VectorView.hpp
@@ -58,6 +58,15 @@ namespace uno {
    VectorView<Expression> view(Expression&& expression, size_t start, size_t end) {
       return {std::forward<Expression>(expression), start, end};
    }
+
+   template <typename Expression>
+   std::ostream& operator<<(std::ostream& stream, const VectorView<Expression>& view) {
+      for (size_t index: Range(view.size())) {
+         stream << view[index] << " ";
+      }
+      stream << '\n';
+      return stream;
+   }
 } // namespace
 
 #endif //UNO_VECTORVIEW_H