Merge branch 'dev'

dev/main.cpp

@@ -17,59 +17,17 @@
 #include "idol/optimizers/branch-and-bound/branching-rules/factories/PseudoCost.h"
 #include "idol/optimizers/branch-and-bound/node-selection-rules/factories/BestEstimate.h"
 #include "idol/optimizers/branch-and-bound/branching-rules/factories/BracnhingWithPriority.h"
+#include "idol/optimizers/wrappers/Mosek/Mosek.h"
 
 using namespace idol;
 
 int main(int t_argc, char** t_argv) {
 
     Env env;
 
-    // Read instance
-    const auto instance = Problems::FLP::read_instance_1991_Cornuejols_et_al("/home/henri/Research/idol/examples/facility.data.txt");
-    const unsigned int n_customers = instance.n_customers();
-    const unsigned int n_facilities = instance.n_facilities();
-
-    // Make model
-
     Model model(env);
 
-    auto x = model.add_vars(Dim<1>(n_facilities), 0., 1., Binary, "x");
-    auto y = model.add_vars(Dim<2>(n_facilities, n_customers), 0., 1., Continuous, "y");
-
-    for (unsigned int i = 0 ; i < n_facilities ; ++i) {
-        model.add_ctr(idol_Sum(j, Range(n_customers), instance.demand(j) * y[i][j]) <= instance.capacity(i) * x[i]);
-    }
-
-    for (unsigned int j = 0 ; j < n_customers ; ++j) {
-        model.add_ctr(idol_Sum(i, Range(n_facilities), y[i][j]) == 1);
-    }
-
-    model.set_obj_expr(idol_Sum(i, Range(n_facilities),
-                                instance.fixed_cost(i) * x[i]
-                                + idol_Sum(j, Range(n_customers),
-                                           instance.per_unit_transportation_cost(i, j) *
-                                           instance.demand(j) *
-                                           y[i][j]
-                                )
-                       )
-    );
-
-    // Set backend options
-    model.use(
-            BranchAndBound()
-                    .with_node_optimizer(HiGHS::ContinuousRelaxation())
-                    .with_branching_rule(
-                        BranchingWithPriority()
-                                        .add_branching_rule(MostInfeasible(x.begin(), x.end()))
-                                        .add_branching_rule(MostInfeasible(y.at(5).begin(), y.at(5).end()))
-                    )
-                    .with_node_selection_rule(BestEstimate())
-                    .with_log_level(Trace, Blue)
-    );
-
-    model.optimize();
-
-    std::cout << save_primal(model) << std::endl;
+    model.use(Mosek());
 
     return 0;
 }

examples/assignment.example.cpp

@@ -64,7 +64,7 @@ int main(int t_argc, const char** t_argv) {
                                       .with_log_level(Info, Yellow)
                                       .with_dual_price_smoothing_stabilization(DantzigWolfe::Neame(.3))
                                       .with_infeasibility_strategy(DantzigWolfe::FarkasPricing())
-                                      .with_hard_branching(false)
+                                      .with_hard_branching(true)
                       )
                       .with_subtree_depth(0)
                       .with_branching_rule(MostInfeasible())

lib/include/idol/optimizers/dantzig-wolfe/DantzigWolfeDecomposition.h

@@ -26,6 +26,7 @@ class idol::DantzigWolfeDecomposition : public OptimizerFactoryWithDefaultParame
     std::unique_ptr<DantzigWolfe::LoggerFactory> m_logger_factory;
     std::optional<unsigned int> m_max_parallel_sub_problems;
     std::optional<bool> m_use_hard_branching;
+    std::optional<bool> m_use_infeasible_column_removal_when_branching;
     std::optional<DantzigWolfe::SubProblem> m_default_sub_problem_spec;
     Map<unsigned int, DantzigWolfe::SubProblem> m_sub_problem_specs;
 
@@ -50,6 +51,8 @@ class idol::DantzigWolfeDecomposition : public OptimizerFactoryWithDefaultParame
 
     DantzigWolfeDecomposition& with_hard_branching(bool t_value);
 
+    DantzigWolfeDecomposition& with_infeasible_columns_removal(bool t_value);
+
     DantzigWolfeDecomposition& with_max_parallel_sub_problems(unsigned int t_n_sub_problems);
 
     DantzigWolfeDecomposition& with_logger(const DantzigWolfe::LoggerFactory& t_logger);

lib/include/idol/optimizers/dantzig-wolfe/Formulation.h

@@ -38,6 +38,7 @@ class idol::DantzigWolfe::Formulation {
     void dispatch_linking_constraint(const Ctr& t_original_ctr, const Row& t_row, CtrType t_type);
     std::pair<Expr<Var, Var>, std::vector<Constant>> decompose_expression(const LinExpr<Var> &t_linear, const QuadExpr<Var, Var>& t_quadratic);
     void dispatch_objective_function(const Model& t_original_formulation);
+    bool is_feasible(const Solution::Primal& t_primal, unsigned int t_sub_problem_id);
 
     void apply_sub_problem_bound_on_master(bool t_is_lb, const idol::Var &t_var, unsigned int t_sub_problem_id, double t_value);
     LinExpr<Var> reformulate_sub_problem_variable(const Var &t_var, unsigned int t_sub_problem_id);
@@ -75,9 +76,9 @@ class idol::DantzigWolfe::Formulation {
 
     double get_original_space_var_primal(const Var& t_var, const Solution::Primal& t_master_primal) const;
 
-    void update_var_lb(const Var& t_var, double t_lb, bool t_hard);
+    void update_var_lb(const Var& t_var, double t_lb, bool t_hard, bool t_remove_infeasible_columns);
 
-    void update_var_ub(const Var& t_var, double t_ub, bool t_hard);
+    void update_var_ub(const Var& t_var, double t_ub, bool t_hard, bool t_remove_infeasible_columns);
 
     void remove_column_if(unsigned int t_sub_problem_id, const std::function<bool(const Var &, const Solution::Primal &)> &t_indicator_for_removal);
 
@@ -94,6 +95,8 @@ class idol::DantzigWolfe::Formulation {
     void remove(const Ctr& t_ctr);
 
     unsigned int get_n_present_generators() const;
+
+    void load_columns_from_pool();
 };
 
 #endif //IDOL_FORMULATION_H

lib/include/idol/optimizers/dantzig-wolfe/Optimizers_DantzigWolfeDecomposition.h

@@ -24,14 +24,16 @@ class idol::Optimizers::DantzigWolfeDecomposition : public Algorithm {
     std::unique_ptr<DantzigWolfe::LoggerFactory::Strategy> m_logger;
     std::vector<DantzigWolfe::SubProblem> m_sub_problem_specifications;
     unsigned int m_max_parallel_pricing;
-    bool m_use_hard_branching = false;
+    bool m_use_hard_branching;
+    bool m_remove_infeasible_columns;
 public:
     DantzigWolfeDecomposition(const Model& t_model,
                               idol::DantzigWolfe::Formulation&& t_formulation,
                               const OptimizerFactory& t_master_optimizer_factory,
                               const DantzigWolfe::DualPriceSmoothingStabilization& t_stabilization,
                               unsigned int t_max_parallel_pricing,
                               bool t_use_hard_branching,
+                              bool t_remove_infeasible_columns,
                               std::vector<DantzigWolfe::SubProblem>&& t_sub_problem_specifications,
                               const DantzigWolfe::InfeasibilityStrategyFactory& t_strategy,
                               const DantzigWolfe::LoggerFactory& t_logger_factory);

lib/src/optimizers/dantzig-wolfe/DantzigWolfeDecomposition.cpp

@@ -24,6 +24,7 @@ idol::DantzigWolfeDecomposition::DantzigWolfeDecomposition(const DantzigWolfeDec
       m_dual_price_smoothing_stabilization(t_src.m_dual_price_smoothing_stabilization ? t_src.m_dual_price_smoothing_stabilization->clone() : nullptr),
       m_max_parallel_sub_problems(t_src.m_max_parallel_sub_problems),
       m_use_hard_branching(t_src.m_use_hard_branching),
+      m_use_infeasible_column_removal_when_branching(t_src.m_use_infeasible_column_removal_when_branching),
       m_infeasibility_strategy(t_src.m_infeasibility_strategy ? t_src.m_infeasibility_strategy->clone() : nullptr),
       m_default_sub_problem_spec(t_src.m_default_sub_problem_spec),
       m_sub_problem_specs(t_src.m_sub_problem_specs),
@@ -57,12 +58,16 @@ idol::Optimizer *idol::DantzigWolfeDecomposition::operator()(const Model &t_mode
         logger_factory = std::make_unique<DantzigWolfe::Loggers::Debug>();
     }
 
+    const bool use_hard_branching = m_use_hard_branching.has_value() && m_use_hard_branching.value();
+    const bool remove_infeasible_column = m_use_infeasible_column_removal_when_branching.has_value() ? m_use_infeasible_column_removal_when_branching.value() : use_hard_branching;
+
     return new Optimizers::DantzigWolfeDecomposition(t_model,
                                                      std::move(dantzig_wolfe_formulation),
                                                      *m_master_optimizer_factory,
                                                      m_dual_price_smoothing_stabilization ? *m_dual_price_smoothing_stabilization : *dual_price_smoothing,
                                                      m_max_parallel_sub_problems.has_value() ? m_max_parallel_sub_problems.value() : 1,
-                                                     m_use_hard_branching.has_value() && m_use_hard_branching.value(),
+                                                     use_hard_branching,
+                                                     remove_infeasible_column,
                                                      std::move(sub_problems_specifications),
                                                      m_infeasibility_strategy ? *m_infeasibility_strategy : *default_strategy,
                                                      *logger_factory
@@ -221,3 +226,14 @@ idol::DantzigWolfeDecomposition::with_logger(const idol::DantzigWolfe::LoggerFac
     return *this;
 }
 
+idol::DantzigWolfeDecomposition &idol::DantzigWolfeDecomposition::with_infeasible_columns_removal(bool t_value) {
+
+    if (m_use_infeasible_column_removal_when_branching.has_value()) {
+        throw Exception("Infeasible column removal has already been configured.");
+    }
+
+    m_use_infeasible_column_removal_when_branching = t_value;
+
+    return *this;
+}
+

lib/src/optimizers/dantzig-wolfe/Formulation.cpp

@@ -416,7 +416,7 @@ double idol::DantzigWolfe::Formulation::get_original_space_var_primal(const idol
 
 }
 
-void idol::DantzigWolfe::Formulation::update_var_lb(const idol::Var &t_var, double t_lb, bool t_hard) {
+void idol::DantzigWolfe::Formulation::update_var_lb(const idol::Var &t_var, double t_lb, bool t_hard, bool t_remove_infeasible_columns) {
 
     const unsigned int sub_problem_id = t_var.get(m_decomposition_by_var);
 
@@ -425,10 +425,11 @@ void idol::DantzigWolfe::Formulation::update_var_lb(const idol::Var &t_var, doub
         return;
     }
 
-    remove_column_if(sub_problem_id, [&](const Var& t_object, const Solution::Primal& t_generator) {
-        //return true;
-        return t_generator.get(t_var) < t_lb;
-    });
+    if (t_remove_infeasible_columns) {
+        remove_column_if(sub_problem_id, [&](const Var &t_object, const Solution::Primal &t_generator) {
+            return Row(t_var, t_lb).is_violated(t_generator, GreaterOrEqual, Tolerance::Feasibility);
+        });
+    }
 
     if (t_hard) {
         sub_problem(sub_problem_id).set_var_lb(t_var, t_lb);
@@ -439,7 +440,7 @@ void idol::DantzigWolfe::Formulation::update_var_lb(const idol::Var &t_var, doub
 
 }
 
-void idol::DantzigWolfe::Formulation::update_var_ub(const idol::Var &t_var, double t_ub, bool t_hard) {
+void idol::DantzigWolfe::Formulation::update_var_ub(const idol::Var &t_var, double t_ub, bool t_hard, bool t_remove_infeasible_columns) {
 
     const unsigned int sub_problem_id = t_var.get(m_decomposition_by_var);
 
@@ -448,10 +449,11 @@ void idol::DantzigWolfe::Formulation::update_var_ub(const idol::Var &t_var, doub
         return;
     }
 
-    remove_column_if(sub_problem_id, [&](const Var& t_object, const Solution::Primal& t_generator) {
-        // return true;
-        return t_generator.get(t_var) > t_ub;
-    });
+    if (t_remove_infeasible_columns) {
+        remove_column_if(sub_problem_id, [&](const Var &t_object, const Solution::Primal &t_generator) {
+            return Row(t_var, t_ub).is_violated(t_generator, LessOrEqual, Tolerance::Feasibility);
+        });
+    }
 
     if (t_hard) {
         sub_problem(sub_problem_id).set_var_ub(t_var, t_ub);
@@ -668,3 +670,64 @@ unsigned int idol::DantzigWolfe::Formulation::get_n_present_generators() const {
 
     return result;
 }
+
+void idol::DantzigWolfe::Formulation::load_columns_from_pool() {
+
+    const unsigned int n_sub_problems = m_sub_problems.size();
+
+    for (unsigned int sub_problem_id = 0 ; sub_problem_id < n_sub_problems ; ++sub_problem_id) {
+
+        for (const auto& [var, generator] : m_pools[sub_problem_id].values()) {
+
+            if (!m_master.has(var) && is_feasible(generator, sub_problem_id)) {
+
+                m_master.add(var, TempVar(
+                    0,
+                    Inf,
+                    Continuous,
+                    m_generation_patterns[sub_problem_id].fix(generator)
+                ));
+                m_present_generators[sub_problem_id].emplace_back(var, generator);
+
+            }
+
+        }
+
+    }
+
+}
+
+bool
+idol::DantzigWolfe::Formulation::is_feasible(const idol::Solution::Primal &t_primal, unsigned int t_sub_problem_id) {
+
+    const auto& model = m_sub_problems[t_sub_problem_id];
+
+    // Check bounds
+    for (const auto& var : model.vars()) {
+
+        const double lb = model.get_var_lb(var);
+        const double ub = model.get_var_ub(var);
+
+        if (
+                Row(var, lb).is_violated(t_primal, GreaterOrEqual, Tolerance::Feasibility)
+                || Row(var, ub).is_violated(t_primal, LessOrEqual, Tolerance::Feasibility)
+            ) {
+            return false;
+        }
+
+    }
+
+    // Check constraints
+    for (const auto& ctr : model.ctrs()) {
+
+        const auto& row = model.get_ctr_row(ctr);
+        const auto& type = model.get_ctr_type(ctr);
+
+        if (row.is_violated(t_primal, type, Tolerance::Feasibility)) {
+            return false;
+        }
+
+    }
+
+    return true;
+}

lib/src/optimizers/dantzig-wolfe/Optimizers_DantzigWolfeDecomposition.cpp

@@ -9,6 +9,7 @@ idol::Optimizers::DantzigWolfeDecomposition::DantzigWolfeDecomposition(const Mod
                                                                        const DantzigWolfe::DualPriceSmoothingStabilization& t_stabilization,
                                                                        unsigned int t_max_parallel_pricing,
                                                                        bool t_use_hard_branching,
+                                                                       bool t_remove_infeasible_columns,
                                                                        std::vector<idol::DantzigWolfe::SubProblem>&& t_sub_problem_specifications,
                                                                        const idol::DantzigWolfe::InfeasibilityStrategyFactory& t_strategy,
                                                                        const DantzigWolfe::LoggerFactory& t_logger_factory)
@@ -17,6 +18,7 @@ idol::Optimizers::DantzigWolfeDecomposition::DantzigWolfeDecomposition(const Mod
       m_master_optimizer_factory(t_master_optimizer_factory.clone()),
       m_max_parallel_pricing(t_max_parallel_pricing),
       m_use_hard_branching(t_use_hard_branching),
+      m_remove_infeasible_columns(t_remove_infeasible_columns),
       m_sub_problem_specifications(std::move(t_sub_problem_specifications)),
       m_stabilization(t_stabilization()),
       m_strategy(t_strategy()),
@@ -44,6 +46,8 @@ void idol::Optimizers::DantzigWolfeDecomposition::hook_before_optimize() {
 
 void idol::Optimizers::DantzigWolfeDecomposition::hook_optimize() {
 
+    m_formulation.load_columns_from_pool();
+
     m_strategy->execute(*this);
 
     set_status(m_strategy->status());
@@ -173,11 +177,11 @@ void idol::Optimizers::DantzigWolfeDecomposition::update_var_type(const idol::Va
 }
 
 void idol::Optimizers::DantzigWolfeDecomposition::update_var_lb(const idol::Var &t_var) {
-    m_formulation.update_var_lb(t_var, parent().get_var_lb(t_var), m_use_hard_branching);
+    m_formulation.update_var_lb(t_var, parent().get_var_lb(t_var), m_use_hard_branching, m_remove_infeasible_columns);
 }
 
 void idol::Optimizers::DantzigWolfeDecomposition::update_var_ub(const idol::Var &t_var) {
-    m_formulation.update_var_ub(t_var, parent().get_var_ub(t_var), m_use_hard_branching);
+    m_formulation.update_var_ub(t_var, parent().get_var_ub(t_var), m_use_hard_branching, m_remove_infeasible_columns);
 }
 
 void idol::Optimizers::DantzigWolfeDecomposition::update_var_obj(const idol::Var &t_var) {

lib/src/optimizers/wrappers/Mosek/Optimizers_Mosek.cpp

@@ -127,7 +127,7 @@ idol::MosekVar idol::Optimizers::Mosek::hook_add(const Var &t_var, bool t_add_co
 
     MosekVar result;
 
-    result.variable = m_model->variable(t_var.name(), 1, mosek::fusion::Domain::unbounded());
+    result.variable = m_model->variable(/* t_var.name(), */ 1, mosek::fusion::Domain::unbounded());
 
     const double lb = parent().get_var_lb(t_var);
     const double ub = parent().get_var_ub(t_var);
@@ -225,11 +225,11 @@ idol::MosekCtr idol::Optimizers::Mosek::hook_add(const Ctr &t_ctr) {
                     );
         }
 
-        result.constraint = m_model->constraint(t_ctr.name(), std::move(expression), mosek::fusion::Domain::inRotatedQCone());
+        result.constraint = m_model->constraint(/* t_ctr.name(), */ std::move(expression), mosek::fusion::Domain::inRotatedQCone());
 
         return result;
 #else
-      throw Exception("idol/modeling/ quadratic expressions with Mosek requires Eigen, please set the USE_EIGEN cmake option to YES to use this feature.");
+      throw Exception("Modeling quadratic expressions with Mosek requires Eigen, please set the USE_EIGEN cmake option to YES to use this feature.");
 #endif
 
     }
@@ -241,13 +241,13 @@ idol::MosekCtr idol::Optimizers::Mosek::hook_add(const Ctr &t_ctr) {
     // Set constraint type
     switch (type) {
         case LessOrEqual:
-            result.constraint = m_model->constraint(t_ctr.name(), std::move(expr), mosek::fusion::Domain::lessThan(0.));
+            result.constraint = m_model->constraint(/* t_ctr.name(), */ std::move(expr), mosek::fusion::Domain::lessThan(0.));
         break;
         case GreaterOrEqual:
-            result.constraint = m_model->constraint(t_ctr.name(), std::move(expr), mosek::fusion::Domain::greaterThan(0.));
+            result.constraint = m_model->constraint(/* t_ctr.name(), */ std::move(expr), mosek::fusion::Domain::greaterThan(0.));
             break;
         case Equal:
-            result.constraint = m_model->constraint(t_ctr.name(), std::move(expr), mosek::fusion::Domain::equalsTo(0.));
+            result.constraint = m_model->constraint(/* t_ctr.name(), */ std::move(expr), mosek::fusion::Domain::equalsTo(0.));
             break;
         default: throw Exception("Enum out of bounds.");
     }