version 0.10.5

CMakeLists.txt

@@ -1,5 +1,5 @@
 cmake_minimum_required(VERSION 3.22)
-project(idol VERSION 0.10.5)
+project(idol VERSION 0.10.6)
 
 set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${idol_SOURCE_DIR}/cmake")

bin/method-managers/robust/BBBB.cpp

@@ -3,6 +3,7 @@
 //
 #include "BBBB.h"
 #include "../milp/MILPMethodManager.h"
+#include "idol/bilevel/optimizers/wrappers/MibS/MibS.h"
 #include "idol/robust/optimizers/bilevel-based-branch-and-bound/MaxMinRelaxation.h"
 #include "idol/mixed-integer/modeling/expressions/operations/operators.h"
 #include "idol/mixed-integer/optimizers/branch-and-bound/BranchAndBound.h"
@@ -11,7 +12,82 @@
 #include "idol/mixed-integer/optimizers/callbacks/heuristics/IntegerMaster.h"
 #include "idol/mixed-integer/optimizers/callbacks/heuristics/RENS.h"
 #include "idol/mixed-integer/optimizers/callbacks/heuristics/SimpleRounding.h"
+#include "idol/mixed-integer/optimizers/wrappers/Gurobi/Gurobi.h"
+#include "idol/robust/optimizers/bilevel-based-branch-and-bound/Optimizers_MaxMinRelaxation.h"
+#include "idol/robust/optimizers/column-and-constraint-generation/ColumnAndConstraintGeneration.h"
+#include "idol/robust/optimizers/column-and-constraint-generation/separation/OptimalitySeparation.h"
+#include "idol/robust/optimizers/column-and-constraint-generation/Optimizers_ColumnAndConstraintGeneration.h"
 #include "idol/robust/optimizers/critical-value-column-and-constraint-generation/CriticalValueColumnAndConstraintGeneration.h"
+#include "idol/robust/optimizers/critical-value-column-and-constraint-generation/Optimizers_CriticalValueColumnAndConstraintGeneration.h"
+
+template<class NodeInfoT = idol::DefaultNodeInfo>
+class EvaluateCallback : public idol::BranchAndBoundCallbackFactory<NodeInfoT> {
+public:
+    class Strategy : public idol::BranchAndBoundCallback<NodeInfoT> {
+        const unsigned int m_max_evaluations = 10;
+        const unsigned int m_node_frequency = 20;
+    protected:
+        void operator()(idol::CallbackEvent t_event) override;
+    };
+
+    idol::BranchAndBoundCallback<NodeInfoT>* operator()() override { return new Strategy(); }
+    [[nodiscard]] idol::BranchAndBoundCallbackFactory<NodeInfoT>* clone() const override { return new EvaluateCallback<NodeInfoT>(*this); }
+};
+
+template <class NodeInfoT>
+void EvaluateCallback<NodeInfoT>::Strategy::operator()(idol::CallbackEvent t_event) {
+
+    if (t_event != idol::InvalidSolution) {
+        return;
+    }
+
+    if (this->node_count() % m_node_frequency != 0) {
+        return;
+    }
+
+    const auto& max_min_relaxation = this->relaxation().optimizer().template as<idol::Optimizers::Robust::MaxMinRelaxation>();
+    const auto& cvccg = max_min_relaxation.get_formulation().model.optimizer().template as<idol::Optimizers::Robust::CriticalValueColumnAndConstraintGeneration>();
+    const auto& branching_candidates = cvccg.branching_candidates();
+
+    auto relaxation_fixed = this->relaxation().copy();
+    relaxation_fixed.optimizer().set_param_time_limit(1e-2);
+    relaxation_fixed.optimize();
+
+    unsigned int n_evaluations = 0;
+    for (const auto& uncertainty : cvccg.get_formulation().uncertainties()) {
+
+        std::vector<idol::CVCCG::Formulation::CurrentlyPresentCut> cuts;
+        std::copy(uncertainty.currently_present_cuts().begin(), uncertainty.currently_present_cuts().end(), std::back_inserter(cuts));
+        std::sort(cuts.begin(), cuts.end(), [](const auto& t_a, const auto& t_b) {
+            return t_a.scenario->scenario.objective_value() < t_b.scenario->scenario.objective_value();
+        });
+
+        for (const auto& cut : cuts) {
+
+            for (const auto& var : branching_candidates) {
+                const double val = cut.scenario->scenario.get(var);
+                relaxation_fixed.set_var_lb(var, val);
+                relaxation_fixed.set_var_ub(var, val);
+            }
+            const double remaining_time = this->original_model().optimizer().get_remaining_time();
+            relaxation_fixed.optimizer().set_param_time_limit(remaining_time);
+            relaxation_fixed.optimize();
+
+            const auto status = relaxation_fixed.get_status();
+            if (status == idol::Optimal || status == idol::Feasible) {
+                auto* info = new idol::DefaultNodeInfo();
+                info->set_primal_solution(idol::save_primal(relaxation_fixed));
+                this->submit_heuristic_solution(info);
+            }
+
+            n_evaluations++;
+            if (n_evaluations >= m_max_evaluations) {
+                return;
+            }
+        }
+    }
+
+}
 
 std::string RobustMethods::BBBB::description() const {
     return "Bilevel-based branch-and-bound.";
@@ -51,6 +127,34 @@ void RobustMethods::BBBB::set_optimizer(idol::Model& t_model, const RobustMethod
 
     const auto sub_milp_optimizer = MILPMethodManager::get_sub_milp_optimizer(args);
 
+    for (const auto& var : t_model.vars()) {
+        if (t_model.get_var_type(var) == idol::Continuous) {
+            t_model.set_var_type(var, idol::Integer);
+        }
+    }
+
+    std::list<idol::PrimalPoint> initial_scenarios;
+
+    if (false) {
+
+        auto ccg = idol::Robust::ColumnAndConstraintGeneration(robust_description, bilevel_description);
+        ccg.with_master_optimizer(idol::Gurobi());
+        ccg.add_separation(idol::Robust::CCG::OptimalitySeparation().with_bilevel_optimizer(idol::Bilevel::MibS()));
+        ccg.with_logs(true);
+        ccg.with_iteration_limit(2);
+
+        t_model.use(ccg);
+        t_model.optimize();
+
+        const auto& optimizer = t_model.optimizer().as<idol::Optimizers::Robust::ColumnAndConstraintGeneration>();
+        for (const auto& scenario : optimizer.get_formulation().generated_scenarios()) {
+            initial_scenarios.emplace_back(scenario);
+        }
+
+        std::cout << "Initializing with " << initial_scenarios.size() << " scenarios." << std::endl;
+
+    }
+
     // Branching Candidates
     std::list<idol::Var> branching_candidates;
     for (const auto& var : t_model.vars()) {
@@ -65,10 +169,16 @@ void RobustMethods::BBBB::set_optimizer(idol::Model& t_model, const RobustMethod
     branch_and_bound.with_logs(true);
     branch_and_bound.with_logger(idol::Logs::BranchAndBound::Info().with_frequency_in_seconds(0).with_node_logs(false));
 
+    branch_and_bound.add_callback(EvaluateCallback());
+
     auto max_min_relaxation = idol::Robust::MaxMinRelaxation(robust_description, bilevel_description);
     max_min_relaxation.with_master_optimizer(*sub_milp_optimizer);
     max_min_relaxation.with_deterministic_optimizer(*sub_milp_optimizer);
     max_min_relaxation.with_indicator(false);
 
+    for (const auto& scenario : initial_scenarios) {
+        max_min_relaxation.add_initial_scenario(scenario);
+    }
+
     t_model.use(branch_and_bound + max_min_relaxation);
 }

lib/include/idol/general/utils/SparseVector.h

@@ -120,6 +120,8 @@ class idol::SparseVector {
 
     void clear() { m_map.clear(); }
 
+    std::pair<double, double> range() const;
+
     void reserve(unsigned int t_capacity) {
 #ifdef IDOL_USE_TSL
         m_map.reserve(t_capacity);
@@ -234,6 +236,17 @@ bool idol::SparseVector<IndexT, ValueT>::is_zero(double t_tolerance) const {
     return true;
 }
 
+template <class IndexT, class ValueT>
+std::pair<double, double> idol::SparseVector<IndexT, ValueT>::range() const {
+    double min = Inf, max = -Inf;
+    for (const auto& [var, val] : m_map) {
+        const double abs_val = std::abs(val);
+        min = std::min(min, abs_val);
+        max = std::max(max, abs_val);
+    }
+    return std::make_pair(min, max);
+}
+
 template<class IndexT, class ValueT>
 idol::SparseVector<IndexT, ValueT>&
 idol::SparseVector<IndexT, ValueT>::operator+=(const SparseVector &t_vector) {

lib/include/idol/mixed-integer/modeling/models/Model.h

@@ -303,6 +303,8 @@ class idol::Model {
     void reset_minor_representation();
 
     static Model read_from_file(Env& t_env, const std::string& t_filename);
+
+    void print_statistics(std::ostream& t_os = std::cout) const;
 };
 
 template<class T, unsigned int N>

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

@@ -466,7 +466,7 @@ void idol::Optimizers::BranchAndBound<NodeInfoT>::submit_heuristic_solution(Node
     }
 
     set_as_incumbent(t_node);
-    log_node_after_solve(t_node);
+    //log_node_after_solve(t_node);
 
     //if (m_branching_rule->is_valid(t_node)) {
         // New incumbent by submission
@@ -585,6 +585,10 @@ void idol::Optimizers::BranchAndBound<NodeInfoT>::hook_before_optimize() {
 template<class NodeInfoT>
 void idol::Optimizers::BranchAndBound<NodeInfoT>::hook_optimize() {
 
+    if (get_param_logs()) {
+        parent().print_statistics(std::cout);
+    }
+
     if (!m_presolve.empty()) {
         m_presolved_model.reset(working_model().clone());
         m_presolve.execute(*m_presolved_model);

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

@@ -324,21 +324,27 @@ idol::BranchAndBoundCallbackI<NodeInfoT>::operator()(Optimizers::BranchAndBound<
                                                      Model *t_relaxation) {
     SideEffectRegistry result;
 
-    m_parent = t_parent;
-    m_node = t_current_node;
-    m_relaxation = t_relaxation;
-    m_registry = &result;
+    const bool is_nested = m_parent != nullptr;
+
+    if (!is_nested) {
+        m_parent = t_parent;
+        m_node = t_current_node;
+        m_relaxation = t_relaxation;
+        m_registry = &result;
+    }
 
     for (auto &cb: m_callbacks) {
-        cb->m_interface = this;
+        if (!is_nested) { cb->m_interface = this; }
         cb->operator()(t_event);
-        cb->m_interface = nullptr;
+        if (!is_nested) { cb->m_interface = nullptr; }
     }
 
-    m_parent = nullptr;
-    m_node.reset();
-    m_relaxation = nullptr;
-    m_registry = nullptr;
+    if (!is_nested) {
+        m_parent = nullptr;
+        m_node.reset();
+        m_relaxation = nullptr;
+        m_registry = nullptr;
+    }
 
     return result;
 }

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

@@ -147,8 +147,10 @@ void idol::Logs::BranchAndBound::Info<NodeInfoT>::Strategy::log_node_after_solve
 
     const double total_time = parent.time().count();
 
-    if (!m_root_node_has_been_printed) {
-        log_root_node(t_node);
+    if (!m_root_node_has_been_printed && t_node.id() == 0) {
+        if (t_node.id() == 0) {
+            log_root_node(t_node);
+        }
         return;
     }
 

lib/include/idol/mixed-integer/optimizers/callbacks/cutting-planes/CglCutCallback.h

@@ -65,7 +65,8 @@ class idol::CglCutCallback : public BranchAndBoundCallbackFactory<NodeInfoT> {
 
         std::list<TempCtr> to_idol_cuts(OsiCuts& t_cuts);
         TempCtr to_idol_cut(OsiRowCut& t_cut);
-        std::vector<std::pair<TempCtr, double>> sort_cuts_by_effectiveness(const std::list<TempCtr>& t_cuts);
+        std::vector<std::pair<TempCtr, double>> sort_and_filter_cuts_by_effectiveness(const std::list<TempCtr>& t_cuts);
+        void standard_scaling(std::list<TempCtr>& t_cuts);
     protected:
         NodeCutContext& get_cut_context();
         const NodeCutContext& get_cut_context() const { return const_cast<Strategy*>(this)->get_cut_context(); }
@@ -184,7 +185,8 @@ void idol::CglCutCallback<NodeInfoT>::Strategy::operator()(CallbackEvent t_event
 
         auto osi_cuts = cut_family->generate(*m_osi_solver, be_aggressive ? 100 : 0);
         auto idol_cuts = to_idol_cuts(osi_cuts);
-        auto sorted_cuts = sort_cuts_by_effectiveness(idol_cuts);
+        standard_scaling(idol_cuts);
+        auto sorted_cuts = sort_and_filter_cuts_by_effectiveness(idol_cuts);
 
         for (auto& [cut, effectiveness] : sorted_cuts) {
 
@@ -290,7 +292,7 @@ idol::TempCtr idol::CglCutCallback<NodeInfoT>::Strategy::to_idol_cut(OsiRowCut&
 }
 
 template <class NodeInfoT>
-std::vector<std::pair<idol::TempCtr, double>> idol::CglCutCallback<NodeInfoT>::Strategy::sort_cuts_by_effectiveness(
+std::vector<std::pair<idol::TempCtr, double>> idol::CglCutCallback<NodeInfoT>::Strategy::sort_and_filter_cuts_by_effectiveness(
     const std::list<TempCtr>& t_cuts) {
     std::vector<std::pair<TempCtr, double>> result;
     const auto primal_solution = this->node().info().primal_solution();
@@ -320,6 +322,34 @@ std::vector<std::pair<idol::TempCtr, double>> idol::CglCutCallback<NodeInfoT>::S
     return result;
 }
 
+template <class NodeInfoT>
+void idol::CglCutCallback<NodeInfoT>::Strategy::standard_scaling(std::list<TempCtr>& t_cuts) {
+
+    for (auto& cut : t_cuts) {
+
+        auto& row = cut.lhs();
+        double& rhs = cut.rhs();
+        double infinity_norm = 0;
+        for (const auto& [var, coeff] : row) {
+            infinity_norm = std::max(std::abs(coeff), infinity_norm);
+        }
+        infinity_norm = std::max(infinity_norm, std::abs(rhs));
+
+        if (is_zero(infinity_norm, Tolerance::Sparsity)) {
+            continue;
+        }
+
+        int e = 0;
+        std::frexp(infinity_norm, &e);
+        double closest_power_of_2 = std::ldexp(1.0, e - 1); // scale = 2^(e-1) or 2^e depending on your normalization choice
+        if (closest_power_of_2 != 1.) {
+            row /= closest_power_of_2;
+            rhs /= closest_power_of_2;
+        }
+
+    }
+}
+
 template <class NodeInfoT>
 idol::CglCutCallback<NodeInfoT>::Strategy::NodeCutContext& idol::CglCutCallback<
     NodeInfoT>::Strategy::get_cut_context() {

lib/include/idol/mixed-integer/optimizers/dantzig-wolfe/ColumnGeneration.h

@@ -8,18 +8,20 @@
 #include "Optimizers_DantzigWolfeDecomposition.h"
 
 class idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration {
+
+    enum NumericalPolicy { Default, ColumnPoolCleanUp, NoDualSmoothing, Failure };
+
     DantzigWolfeDecomposition& m_parent;
     double m_best_bound_stop;
 
-    enum NumericalPolicy { Default, ColumnPoolCleanUp, NoDualSmoothing, Failure };
-    const unsigned int m_max_n_iterations_without_generating_column = 1000;
+    const unsigned int m_max_n_iterations_without_generating_column = 200;
     unsigned int m_n_iterations_without_generating_column = 0;
     NumericalPolicy m_numerical_policy = Default;
 
     SolutionStatus m_status = Loaded;
     SolutionReason m_reason = NotSpecified;
     std::optional<PrimalPoint> m_master_primal_solution;
-    std::optional<DualPoint> m_last_master_solution;
+    std::optional<DualPoint> m_master_dual_solution;
     std::vector<DantzigWolfe::SubProblem::PhaseId> m_sub_problems_phases;
     double m_best_obj = -Inf;
     double m_best_bound = +Inf;
@@ -39,6 +41,7 @@ class idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration {
     void solve_sub_problems_in_parallel();
     void analyze_sub_problems();
     void enrich_master();
+    bool check_numerical_stability();
     void pool_clean_up();
 
     void next_numerical_policy();
@@ -47,6 +50,8 @@ class idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration {
     void log_master();
     void log_sub_problems();
     void log_end();
+
+    friend std::ostream& operator<<(std::ostream& t_os, idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration::NumericalPolicy t_numerical_policy);
 public:
     ColumnGeneration(DantzigWolfeDecomposition& t_parent, bool t_use_farkas_for_infeasibility, double t_best_bound_stop);
 
@@ -69,4 +74,5 @@ class idol::Optimizers::DantzigWolfeDecomposition::ColumnGeneration {
     void execute();
 };
 
+
 #endif //IDOL_COLUMNGENERATION_H

lib/include/idol/robust/optimizers/bilevel-based-branch-and-bound/MaxMinRelaxation.h

@@ -19,6 +19,7 @@ class idol::Robust::MaxMinRelaxation : public OptimizerFactoryWithDefaultParamet
     std::unique_ptr<OptimizerFactory> m_master_optimizer_factory;
     std::unique_ptr<OptimizerFactory> m_deterministic_optimizer_factory;
     std::optional<bool> m_use_indicator;
+    std::list<PrimalPoint> m_initial_scenarios;
 protected:
     [[nodiscard]] Optimizer* create(const Model& t_model) const override;
     MaxMinRelaxation(const MaxMinRelaxation& t_src);
@@ -30,6 +31,7 @@ class idol::Robust::MaxMinRelaxation : public OptimizerFactoryWithDefaultParamet
     MaxMinRelaxation& with_master_optimizer(const OptimizerFactory& t_optimizer);
     MaxMinRelaxation& with_deterministic_optimizer(const OptimizerFactory& t_optimizer);
     MaxMinRelaxation& with_indicator(bool t_value);
+    MaxMinRelaxation& add_initial_scenario(PrimalPoint t_scenario);
 };
 
 #endif //IDOL_ROBUST_BBBB_H