add reduced cost fixing

examples/CMakeLists.txt

@@ -8,4 +8,4 @@ function(ADD_EXAMPLE NAME)
 endfunction()
 
 add_subdirectory(mixed-integer)
-add_subdirectory(bilevel)
+add_subdirectory(bilevel)

examples/mixed-integer/facility.example.cpp

@@ -11,6 +11,8 @@
 #include "idol/mixed-integer/optimizers/wrappers/GLPK/GLPK.h"
 #include "idol/mixed-integer/optimizers/wrappers/Gurobi/Gurobi.h"
 #include "idol/mixed-integer/optimizers/branch-and-bound/node-selection-rules/factories/BestBound.h"
+#include "idol/mixed-integer/optimizers/branch-and-bound/branching-rules/factories/MostInfeasible.h"
+#include "idol/mixed-integer/optimizers/callbacks/ReducedCostFixing.h"
 
 using namespace idol;
 
@@ -57,10 +59,10 @@ int main(int t_argc, const char** t_argv) {
     // Set backend options
     model.use(
             BranchAndBound()
-                    .with_node_optimizer(GLPK::ContinuousRelaxation())
-                    //.add_callback(Cuts::KnapsackCover())
+                    .with_node_optimizer(Gurobi::ContinuousRelaxation())
                     .with_branching_rule(PseudoCost())
                     .with_node_selection_rule(BestEstimate())
+                    .add_callback(ReducedCostFixing())
                     .with_logs(true)
     );
 

lib/CMakeLists.txt

@@ -227,6 +227,7 @@ add_library(idol STATIC
         include/idol/general/utils/GenerationPattern.h
         include/idol/mixed-integer/modeling/models/KKT.h
         src/mixed-integer/modeling/models/KKT.cpp
+        include/idol/mixed-integer/optimizers/callbacks/ReducedCostFixing.h
 )
 
 find_package(OpenMP REQUIRED)

lib/include/idol/mixed-integer/optimizers/branch-and-bound/BranchAndBound.h

@@ -36,7 +36,6 @@ class idol::BranchAndBound : public OptimizerFactoryWithDefaultParameters<Branch
 
     std::optional<unsigned int> m_subtree_depth;
     std::optional<unsigned int> m_log_frequency;
-    std::optional<bool> m_scaling;
 public:
     /**
      * This type is used to exploit [SFINAE](https://en.wikipedia.org/wiki/Substitution_failure_is_not_an_error)
@@ -175,8 +174,6 @@ class idol::BranchAndBound : public OptimizerFactoryWithDefaultParameters<Branch
 
     BranchAndBound<NodeT>& with_logger(const Logs::BranchAndBound::Factory<NodeT>& t_log_factory);
 
-    BranchAndBound<NodeT>& with_scaling(bool t_value);
-
     /**
      * Adds a callback which will be called by the optimizer.
      *
@@ -206,6 +203,7 @@ class idol::BranchAndBound : public OptimizerFactoryWithDefaultParameters<Branch
      * @return the optimizer factory itself
      */
     BranchAndBound<NodeT>& add_callback(const CallbackFactory& t_callback);
+
 };
 
 template<class NodeT>
@@ -226,18 +224,6 @@ idol::BranchAndBound<NodeT>::with_logger(const typename idol::Logs::BranchAndBou
     return *this;
 }
 
-template<class NodeT>
-idol::BranchAndBound<NodeT> &idol::BranchAndBound<NodeT>::with_scaling(bool t_value) {
-
-    if (m_scaling.has_value()) {
-        throw Exception("Scaling has already been configured.");
-    }
-
-    m_scaling = t_value;
-
-    return *this;
-}
-
 template<class NodeT>
 idol::BranchAndBound<NodeT> &idol::BranchAndBound<NodeT>::operator+=(const idol::OptimizerFactory &t_node_optimizer) {
     return with_node_optimizer(t_node_optimizer);
@@ -328,7 +314,6 @@ idol::BranchAndBound<NodeT>::BranchAndBound(const BranchAndBound &t_rhs)
           m_branching_rule_factory(t_rhs.m_branching_rule_factory ? t_rhs.m_branching_rule_factory->clone() : nullptr),
           m_node_selection_rule_factory(t_rhs.m_node_selection_rule_factory ? t_rhs.m_node_selection_rule_factory->clone() : nullptr),
           m_subtree_depth(t_rhs.m_subtree_depth),
-          m_scaling(t_rhs.m_scaling),
           m_logger_factory(t_rhs.m_logger_factory ? t_rhs.m_logger_factory->clone() : nullptr) {
 
     for (auto& cb : t_rhs.m_callbacks) {
@@ -352,20 +337,17 @@ idol::Optimizer *idol::BranchAndBound<NodeT>::operator()(const Model &t_model) c
         throw Exception("No node selection rule has been given, please call BranchAndBound::with_node_selection_rule to configure.");
     }
 
-    auto* callback_interface = new BranchAndBoundCallbackI<NodeT>();
-
     std::unique_ptr<Logs::BranchAndBound::Factory<NodeT>> default_logger_factory;
     if (!m_logger_factory) {
         default_logger_factory = std::make_unique<Logs::BranchAndBound::Info<NodeT>>();
     }
 
     auto* result = new Optimizers::BranchAndBound<NodeT>(t_model,
-                                     *m_relaxation_optimizer_factory,
-                                     *m_branching_rule_factory,
-                                     *m_node_selection_rule_factory,
-                                     callback_interface,
-                                     m_scaling.has_value() && m_scaling.value(),
-                                     m_logger_factory ? *m_logger_factory : *default_logger_factory);
+                                                         *m_relaxation_optimizer_factory,
+                                                         *m_branching_rule_factory,
+                                                         *m_node_selection_rule_factory,
+                                                         new BranchAndBoundCallbackI<NodeT>(),
+                                                         m_logger_factory ? *m_logger_factory : *default_logger_factory);
 
     this->handle_default_parameters(result);
 

lib/include/idol/mixed-integer/optimizers/branch-and-bound/Optimizers_BranchAndBound.h

@@ -106,7 +106,6 @@ class idol::Optimizers::BranchAndBound : public Algorithm {
                               const BranchingRuleFactory<NodeInfoT>& t_branching_rule_factory,
                               const NodeSelectionRuleFactory<NodeInfoT>& t_node_selection_rule_factory,
                               AbstractBranchAndBoundCallbackI<NodeInfoT>* t_callback,
-                              bool t_scaling,
                               const Logs::BranchAndBound::Factory<NodeInfoT>& t_logger_factory);
 
     [[nodiscard]] std::string name() const override { return "Branch-and-Bound"; }
@@ -424,7 +423,6 @@ idol::Optimizers::BranchAndBound<NodeInfoT>::BranchAndBound(const Model &t_model
                                                       const BranchingRuleFactory<NodeInfoT>& t_branching_rule_factory,
                                                       const NodeSelectionRuleFactory<NodeInfoT>& t_node_selection_rule_factory,
                                                       AbstractBranchAndBoundCallbackI<NodeInfoT>* t_callback,
-                                                      bool t_scaling,
                                                       const Logs::BranchAndBound::Factory<NodeInfoT>& t_logger_factory)
     : Algorithm(t_model),
       m_relaxation_optimizer_factory(t_node_optimizer.clone()),
@@ -452,7 +450,7 @@ void idol::Optimizers::BranchAndBound<NodeInfoT>::create_relaxations() {
         auto* relaxation = original_model.clone();
         relaxation->use(*m_relaxation_optimizer_factory);
         m_relaxations.emplace_back(relaxation);
-        m_node_updators.emplace_back(dynamic_cast<NodeUpdator<NodeInfoT>*>(NodeInfoT::create_updator(*relaxation)));
+        m_node_updators.emplace_back(dynamic_cast<NodeUpdator<NodeInfoT>*>(NodeInfoT::create_updator(parent(), *relaxation)));
     }
 
 }
@@ -463,6 +461,7 @@ idol::Node<NodeInfoT> idol::Optimizers::BranchAndBound<NodeInfoT>::create_root_n
     auto root_node = Node<NodeInfoT>::create_root_node();
     assert(root_node.id() == 0);
     ++m_n_created_nodes;
+
     return root_node;
 
 }
@@ -637,14 +636,27 @@ template<class NodeInfoT>
 void idol::Optimizers::BranchAndBound<NodeInfoT>::solve(TreeNode& t_node,
                                                         unsigned int t_relaxation_id) const {
 
-    m_node_updators[t_relaxation_id]->prepare(t_node);
+    auto& node_updator = *m_node_updators[t_relaxation_id];
+    auto& relaxation = *m_relaxations[t_relaxation_id];
+
+    relaxation.optimizer().set_param_best_bound_stop(std::min(get_best_obj(), get_param_best_bound_stop()));
+    relaxation.optimizer().set_param_time_limit(get_remaining_time());
 
-    m_relaxations[t_relaxation_id]->optimizer().set_param_best_bound_stop(std::min(get_best_obj(), get_param_best_bound_stop()));
-    m_relaxations[t_relaxation_id]->optimizer().set_param_time_limit(get_remaining_time());
+    std::cout << "Preparing node " << t_node.id() << std::endl;
+    node_updator.prepare(t_node);
+    std::cout << "END" << std::endl;
 
-    m_relaxations[t_relaxation_id]->optimize();
+    for (const auto& var : parent().vars()) {
+        const double lb = relaxation.get_var_lb(var);
+        const double ub = relaxation.get_var_ub(var);
+        if (lb > ub + Tolerance::Integer) {
+            std::cout << "Inconsistent bounds for variable " << var << ": " << lb << " > " << ub << std::endl;
+        }
+    }
 
-    t_node.info().save(parent(), *m_relaxations[t_relaxation_id]);
+    relaxation.optimize();
+
+    t_node.info().save(parent(), relaxation);
 
     m_branching_rule->on_node_solved(t_node);
 
@@ -696,15 +708,7 @@ void idol::Optimizers::BranchAndBound<NodeInfoT>::analyze(const BranchAndBound::
 
     if (status == Fail || status == Loaded) {
 
-        /*
-        if (m_n_postponed_nodes < m_max_postponed_nodes) {
-            *t_explore_children_flag = true;
-            std::cout << "Postponing Node " << t_node.id() << " since returned status is " << status << "." << std::endl;
-            idol_Log(Trace, "Postponing Node " << t_node.id() << " since returned status is " << status << ".");
-            ++m_n_postponed_nodes;
-            return;
-        }
-         */
+        std::cout << "HERE?" << std::endl;
 
         set_status(Fail);
         set_reason(NotSpecified);
@@ -746,7 +750,7 @@ void idol::Optimizers::BranchAndBound<NodeInfoT>::analyze(const BranchAndBound::
         m_root_node_best_obj = get_best_obj();
     }
 
-    if (side_effects.n_added_lazy_cuts > 0 || side_effects.n_added_user_cuts > 0) {
+    if (side_effects.n_added_lazy_cuts > 0 || side_effects.n_added_user_cuts > 0 || side_effects.n_added_local_variable_branching > 0) {
         *t_reoptimize_flag = true;
         return;
     }

lib/include/idol/mixed-integer/optimizers/branch-and-bound/callbacks/AbstractBranchAndBoundCallbackI.h

@@ -28,6 +28,7 @@ namespace idol {
 struct idol::SideEffectRegistry {
     unsigned int n_added_lazy_cuts = 0;
     unsigned int n_added_user_cuts = 0;
+    unsigned int n_added_local_variable_branching = 0;
 };
 
 template<class NodeInfoT>

lib/include/idol/mixed-integer/optimizers/branch-and-bound/callbacks/BranchAndBoundCallback.h

@@ -41,6 +41,8 @@ class idol::BranchAndBoundCallbackI : public AbstractBranchAndBoundCallbackI<Nod
 
     void add_lazy_cut(const TempCtr &t_cut);
 
+    void add_local_variable_branching(const idol::Var &t_var, CtrType t_type, double t_rhs);
+
     [[nodiscard]] const Node<NodeInfoT>& node() const;
 
     [[nodiscard]] const Model& relaxation() const;
@@ -74,6 +76,12 @@ class idol::BranchAndBoundCallbackI : public AbstractBranchAndBoundCallbackI<Nod
     void log_after_termination() override;
 };
 
+template<class NodeInfoT>
+void idol::BranchAndBoundCallbackI<NodeInfoT>::add_local_variable_branching(const idol::Var &t_var, CtrType t_type, double t_rhs) {
+    m_node->info().add_branching_variable(t_var, t_type, t_rhs);
+    m_registry->n_added_local_variable_branching++;
+}
+
 template<class NodeInfoT>
 void idol::BranchAndBoundCallbackI<NodeInfoT>::terminate() {
     m_parent->terminate();
@@ -150,6 +158,8 @@ class idol::BranchAndBoundCallback {
      */
     void add_lazy_cut(const TempCtr& t_cut);
 
+    void add_local_variable_branching(const Var &t_var, CtrType t_type, double t_rhs);
+
     /**
      * Returns the node which is currently explored
      * @return the node which is currently explored
@@ -206,6 +216,12 @@ class idol::BranchAndBoundCallback {
     friend class BranchAndBoundCallbackI<NodeInfoT>;
 };
 
+template<class NodeInfoT>
+void idol::BranchAndBoundCallback<NodeInfoT>::add_local_variable_branching(const Var &t_var, CtrType t_type, double t_rhs) {
+    throw_if_no_interface();
+    m_interface->add_local_variable_branching(t_var, t_type, t_rhs);
+}
+
 template<class NodeInfoT>
 double idol::BranchAndBoundCallback<NodeInfoT>::best_obj() const {
     throw_if_no_interface();

lib/include/idol/mixed-integer/optimizers/branch-and-bound/nodes/DefaultNodeInfo.h

@@ -15,7 +15,6 @@ namespace idol {
     namespace Optimizers {
         template<class NodeInfoT> class BranchAndBound;
     }
-
 }
 
 class idol::DefaultNodeInfo {
@@ -50,7 +49,8 @@ class idol::DefaultNodeInfo {
 
     [[nodiscard]] auto constraint_branching_decisions() const { return ConstIteratorForward(m_constraint_branching_decisions); }
 
-    static DefaultNodeUpdator<DefaultNodeInfo>* create_updator(Model& t_relaxation);
+    static DefaultNodeUpdator<DefaultNodeInfo>* create_updator(const Model& t_src_model, Model& t_relaxation);
+
 private:
     PrimalPoint m_primal_solution;
     std::optional<double> m_sum_of_infeasibilities;

lib/include/idol/mixed-integer/optimizers/branch-and-bound/nodes/DefaultNodeUpdator.h

@@ -20,6 +20,7 @@ namespace idol {
 template<class NodeInfoT>
 class idol::DefaultNodeUpdator : public NodeUpdator<NodeInfoT> {
 
+    const Model& m_parent;
     Model& m_relaxation;
 
     Map<Var, double> m_changed_lower_bounds; // variable -> old_bound
@@ -37,7 +38,7 @@ class idol::DefaultNodeUpdator : public NodeUpdator<NodeInfoT> {
 
     void apply_constraint_branching_decisions(const Node<NodeInfoT> &t_node);
 public:
-    explicit DefaultNodeUpdator(Model& t_relaxation);
+    explicit DefaultNodeUpdator(const Model& t_parent, Model& t_relaxation);
 
     DefaultNodeUpdator(const DefaultNodeUpdator&) = delete;
 
@@ -71,7 +72,8 @@ idol::DefaultNodeUpdator<NodeInfoT>::apply_constraint_branching_decisions(const
 }
 
 template<class NodeT>
-idol::DefaultNodeUpdator<NodeT>::DefaultNodeUpdator(idol::Model &t_relaxation) : m_relaxation(t_relaxation) {
+idol::DefaultNodeUpdator<NodeT>::DefaultNodeUpdator(const Model& t_parent, idol::Model &t_relaxation)
+    : m_parent(t_parent), m_relaxation(t_relaxation) {
 
 }
 
@@ -117,12 +119,10 @@ void idol::DefaultNodeUpdator<NodeInfoT>::apply_variable_branching_decisions(con
                                                                                 Map<Var, double> &t_changed_upper_bounds
                                                                                 ) {
 
-    if (t_node.level() == 0) {
-        return;
-    }
-
     for (const auto& branching_decision : t_node.info().variable_branching_decisions()) {
 
+        std::cout << branching_decision.variable << " " << branching_decision.type << " " << branching_decision.bound << std::endl;
+
         switch (branching_decision.type) {
             case LessOrEqual:
                 apply_variable_branching_decision(branching_decision, m_changed_upper_bounds, t_changed_upper_bounds);
@@ -136,6 +136,10 @@ void idol::DefaultNodeUpdator<NodeInfoT>::apply_variable_branching_decisions(con
 
     }
 
+    if (t_node.level() == 0) {
+        return;
+    }
+
     apply_variable_branching_decisions(t_node.parent(), t_changed_lower_bounds, t_changed_upper_bounds);
 
 }

lib/include/idol/mixed-integer/optimizers/callbacks/ReducedCostFixing.h

@@ -0,0 +1,59 @@
+//
+// Created by henri on 24.11.24.
+//
+
+#ifndef IDOL_REDUCEDCOSTFIXING_H
+#define IDOL_REDUCEDCOSTFIXING_H
+
+#include "idol/mixed-integer/optimizers/branch-and-bound/callbacks/BranchAndBoundCallback.h"
+
+namespace idol {
+    template<class> class ReducedCostFixing;
+}
+
+template<class NodeInfoT = idol::DefaultNodeInfo>
+class idol::ReducedCostFixing : public BranchAndBoundCallbackFactory<NodeInfoT> {
+public:
+    class Strategy : public BranchAndBoundCallback<NodeInfoT> {
+    public:
+    protected:
+        void operator()(CallbackEvent t_event) override {
+
+            if (t_event != InvalidSolution) {
+                return;
+            }
+
+            const auto& original_model = this->original_model();
+            const auto& relaxation = this->relaxation();
+            const double best_obj = this->best_obj();
+            const double current_obj = this->relaxation().get_best_obj();
+
+            for (const auto &var : original_model.vars()) {
+
+                double reduced_cost = relaxation.get_var_reduced_cost(var);
+
+                if (current_obj + reduced_cost > best_obj + Tolerance::MIPAbsoluteGap) {
+                    const double relaxation_lb = original_model.get_var_lb(var);
+                    this->add_local_variable_branching(var, LessOrEqual, relaxation_lb);
+                }
+
+                if (current_obj - reduced_cost > best_obj + Tolerance::MIPAbsoluteGap) {
+                    const double relaxation_ub = original_model.get_var_ub(var);
+                    this->add_local_variable_branching(var, GreaterOrEqual, relaxation_ub);
+                }
+
+            }
+
+        }
+    };
+
+    BranchAndBoundCallback<NodeInfoT> *operator()() override {
+        return new Strategy();
+    }
+
+    BranchAndBoundCallbackFactory<NodeInfoT> *clone() const override {
+        return new ReducedCostFixing();
+    }
+};
+
+#endif //IDOL_REDUCEDCOSTFIXING_H