diff --git a/docs/reference/environments.rst b/docs/reference/environments.rst index 3fa6416c0..dce79ccda 100644 --- a/docs/reference/environments.rst +++ b/docs/reference/environments.rst @@ -16,6 +16,11 @@ Branching .. autoclass:: ecole.environment.Branching .. autoclass:: ecole.dynamics.BranchingDynamics +BranchingSum +^^^^^^^^^^^^ +.. autoclass:: ecole.environment.BranchingSum +.. autoclass:: ecole.dynamics.BranchingSumDynamics + Configuring ^^^^^^^^^^^ .. autoclass:: ecole.environment.Configuring diff --git a/libecole/CMakeLists.txt b/libecole/CMakeLists.txt index a757ac042..088da2e76 100644 --- a/libecole/CMakeLists.txt +++ b/libecole/CMakeLists.txt @@ -37,6 +37,7 @@ add_library( src/observation/pseudocosts.cpp src/dynamics/branching.cpp + src/dynamics/branching-sum.cpp src/dynamics/configuring.cpp src/dynamics/primal-search.cpp ) diff --git a/libecole/include/ecole/dynamics/branching-sum.hpp b/libecole/include/ecole/dynamics/branching-sum.hpp new file mode 100644 index 000000000..c9602ec44 --- /dev/null +++ b/libecole/include/ecole/dynamics/branching-sum.hpp @@ -0,0 +1,26 @@ +#pragma once + +#include +#include + +#include +#include + +#include "ecole/dynamics/dynamics.hpp" +#include "ecole/export.hpp" + +namespace ecole::dynamics { + +class ECOLE_EXPORT BranchingSumDynamics : + public EnvironmentDynamics, std::optional>> { +public: + using Action = nonstd::span; + using ActionSet = std::optional>; + + ECOLE_EXPORT auto reset_dynamics(scip::Model& model) -> std::tuple override; + + ECOLE_EXPORT auto step_dynamics(scip::Model& model, Action const& var_indices) + -> std::tuple override; +}; + +} // namespace ecole::dynamics diff --git a/libecole/include/ecole/environment/branching-sum.hpp b/libecole/include/ecole/environment/branching-sum.hpp new file mode 100644 index 000000000..88c208708 --- /dev/null +++ b/libecole/include/ecole/environment/branching-sum.hpp @@ -0,0 +1,18 @@ +#pragma once + +#include "ecole/dynamics/branching-sum.hpp" +#include "ecole/environment/environment.hpp" +#include "ecole/information/nothing.hpp" +#include "ecole/observation/node-bipartite.hpp" +#include "ecole/reward/is-done.hpp" + +namespace ecole::environment { + +template < + typename ObservationFunction = observation::NodeBipartite, + typename RewardFunction = reward::IsDone, + typename InformationFunction = information::Nothing> +using BranchingSum = + Environment; + +} // namespace ecole::environment diff --git a/libecole/src/dynamics/branching-sum.cpp b/libecole/src/dynamics/branching-sum.cpp new file mode 100644 index 000000000..24262b82a --- /dev/null +++ b/libecole/src/dynamics/branching-sum.cpp @@ -0,0 +1,223 @@ +/*********************************************************** + * Implementation of SCIPbranchSum (aim to merge in SCIP) * + ***********************************************************/ + +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace { + +/** Validate that variables are integer, not fixed, and their LP solution value finite. */ +SCIP_RETCODE SCIPbranchSum_ValidateVar(SCIP* scip, SCIP_VAR* var, SCIP_Real var_sol) { + assert(SCIPvarIsActive(var)); + assert(SCIPvarGetProbindex(var) >= 0); + if (SCIPvarGetType(var) == SCIP_VARTYPE_CONTINUOUS) { + SCIPerrorMessage("cannot branch on constraint containing continuous variable <%s>\n", SCIPvarGetName(var)); + return SCIP_INVALIDDATA; + } + if (SCIPisEQ(scip, SCIPvarGetLbLocal(var), SCIPvarGetUbLocal(var))) { + SCIPerrorMessage( + "cannot branch on constraint containing variable <%s> with fixed domain [%.15g,%.15g]\n", + SCIPvarGetName(var), + SCIPvarGetLbLocal(var), + SCIPvarGetUbLocal(var)); + return SCIP_INVALIDDATA; + } + + if (SCIPisInfinity(scip, -var_sol) || SCIPisInfinity(scip, var_sol)) { + SCIPerrorMessage("cannot branch on variables containing infinite values"); + return SCIP_INVALIDDATA; + } + + return SCIP_OKAY; +} + +SCIP_RETCODE SCIPbranchSum_AddChild( + SCIP* scip, + SCIP_Real priority, + SCIP_Real estimate, + SCIP_VAR** vars, + SCIP_Real* ones, + int nvars, + SCIP_Real lhs, + SCIP_Real rhs, + SCIP_NODE** node_out) { + SCIP_NODE* node = nullptr; + SCIP_CALL(SCIPcreateChild(scip, &node, priority, estimate)); + auto name = fmt::format("branching-{}", SCIPnodeGetNumber(node)); + SCIP_CONS* cons = nullptr; + SCIP_CALL(SCIPcreateConsLinear( + scip, + &cons, + name.c_str(), + nvars, + vars, + ones, + lhs, + rhs, + /*initial=*/TRUE, + /*separate=*/TRUE, + /*enforce=*/TRUE, + /*check=*/FALSE, + /*propagate=*/TRUE, + /*local=*/TRUE, + /*modifiable=*/FALSE, + /*dynamic=*/FALSE, + /*removable=*/FALSE, + /*stickingatnode=*/TRUE)); + SCIP_RETCODE retcode = SCIP_OKAY; + SCIP_CALL_TERMINATE(retcode, SCIPaddConsNode(scip, node, cons, nullptr), TERM); + if (node_out != nullptr) { + *node_out = node; + } + +TERM: + SCIP_CALL(SCIPreleaseCons(scip, &cons)); + return retcode; +} + +SCIP_RETCODE +SCIPbranchSum(SCIP* scip, SCIP_VAR** vars, int nvars, SCIP_NODE** downchild, SCIP_NODE** upchild) { + /* Check input parameters */ + assert(scip != nullptr); + if (SCIPgetStage(scip) != SCIP_STAGE_SOLVING) { + SCIPerrorMessage("cannot branch when not solving\n"); + return SCIP_INVALIDCALL; + } + if (nvars <= 0) { + SCIPerrorMessage("cannot branch on empty variable set\n"); + return SCIP_INVALIDDATA; + } + + /* Compute the node estimate and priority but taking the average of the variables'. + * Compute the sum of the LP solution value of the variables. */ + SCIP_Real pseudo_sol_sum = 0.; + SCIP_Real estimate_down = 0.; + SCIP_Real estimate_up = 0.; + SCIP_Real priority_down = 0.; + SCIP_Real priority_up = 0.; + { + SCIP_VAR const* const* const vars_end = vars + nvars; + for (SCIP_VAR** var_iter = vars; var_iter < vars_end; ++var_iter) { + SCIP_Real const sol = SCIPvarGetSol(*var_iter, SCIPhasCurrentNodeLP(scip)); + SCIP_Real const sol_floor = SCIPfeasFloor(scip, sol); + SCIP_Real const sol_ceil = SCIPfeasCeil(scip, sol); + SCIPbranchSum_ValidateVar(scip, *var_iter, sol); + estimate_down += SCIPcalcChildEstimate(scip, *var_iter, sol_floor); + estimate_up += SCIPcalcChildEstimate(scip, *var_iter, sol_ceil); + priority_down += SCIPcalcNodeselPriority(scip, *var_iter, SCIP_BRANCHDIR_DOWNWARDS, sol_floor); + priority_up += SCIPcalcNodeselPriority(scip, *var_iter, SCIP_BRANCHDIR_UPWARDS, sol_ceil); + pseudo_sol_sum += sol; + } + estimate_down /= static_cast(nvars); + estimate_up /= static_cast(nvars); + priority_down /= static_cast(nvars); + priority_up /= static_cast(nvars); + } + + /* If the sum of the LP solution values is integral, we cannot branch without portentially excluding feasible + * solutions */ + SCIP_Real const downbound = SCIPfeasFloor(scip, pseudo_sol_sum); + SCIP_Real const upbound = SCIPfeasCeil(scip, pseudo_sol_sum); + if (SCIPisEQ(scip, downbound, upbound)) { + SCIPerrorMessage("cannot branch on a variables whose sum of LP solution value is integer"); + return SCIP_INVALIDDATA; + } + + SCIP_Real* ones = nullptr; + SCIP_CALL(SCIPallocBufferArray(scip, &ones, nvars)); + for (int i = 0; i < nvars; ++i) { + ones[i] = 1.; + } + + SCIP_Retcode retcode = SCIP_OKAY; + SCIP_Real const inf = SCIPinfinity(scip); + SCIP_CALL_TERMINATE( + retcode, + SCIPbranchSum_AddChild(scip, priority_down, estimate_down, vars, ones, nvars, -inf, downbound, downchild), + TERM); + SCIP_CALL_TERMINATE( + retcode, SCIPbranchSum_AddChild(scip, priority_up, estimate_up, vars, ones, nvars, upbound, inf, upchild), TERM); + +TERM: + SCIPfreeBufferArray(scip, &ones); + return retcode; +} + +} // namespace + +/******************************************** + * Implementation of BranchingSumDynamics * + ********************************************/ + +#include +#include +#include + +#include + +#include "ecole/dynamics/branching-sum.hpp" +#include "ecole/scip/model.hpp" +#include "ecole/scip/utils.hpp" + +namespace ecole::dynamics { + +namespace { + +std::optional> action_set(scip::Model const& model) { + if (model.stage() != SCIP_STAGE_SOLVING) { + return {}; + } + auto const branch_cands = model.lp_branch_cands(); + auto branch_cols = xt::xtensor::from_shape({branch_cands.size()}); + auto const var_to_idx = [](auto const var) { return SCIPcolGetLPPos(SCIPvarGetCol(var)); }; + std::transform(branch_cands.begin(), branch_cands.end(), branch_cols.begin(), var_to_idx); + + assert(branch_cols.size() > 0); + return branch_cols; +} + +} // namespace + +auto BranchingSumDynamics::reset_dynamics(scip::Model& model) -> std::tuple { + model.solve_iter_start_branch(); + if (model.solve_iter_is_done()) { + return {true, {}}; + } + return {false, action_set(model)}; +} + +auto BranchingSumDynamics::step_dynamics(scip::Model& model, Action const& var_indices) -> std::tuple { + auto const lp_cols = model.lp_columns(); + + // Check that input indices are within range + auto const is_out_of_bounds = [size = lp_cols.size()](auto idx) { return idx >= size; }; + if (std::any_of(var_indices.begin(), var_indices.end(), is_out_of_bounds)) { + throw std::invalid_argument{"Branching index is larger than the number of columns."}; + } + + { // Get variables associated with indices and branch + auto vars = std::vector(var_indices.size()); + auto const idx_to_var = [lp_cols](auto idx) { return SCIPcolGetVar(lp_cols[idx]); }; + std::transform(var_indices.begin(), var_indices.end(), vars.begin(), idx_to_var); + + scip::call(SCIPbranchSum, model.get_scip_ptr(), vars.data(), static_cast(vars.size()), nullptr, nullptr); + model.solve_iter_branch(SCIP_BRANCHED); + } + + if (model.solve_iter_is_done()) { + return {true, {}}; + } + return {false, action_set(model)}; +} + +} // namespace ecole::dynamics diff --git a/libecole/tests/CMakeLists.txt b/libecole/tests/CMakeLists.txt index 2f4c20188..c881f8ebd 100644 --- a/libecole/tests/CMakeLists.txt +++ b/libecole/tests/CMakeLists.txt @@ -47,6 +47,7 @@ add_executable( src/observation/test-hutter-2011.cpp src/dynamics/test-branching.cpp + src/dynamics/test-branching-sum.cpp src/dynamics/test-configuring.cpp src/dynamics/test-primal-search.cpp diff --git a/libecole/tests/src/dynamics/test-branching-sum.cpp b/libecole/tests/src/dynamics/test-branching-sum.cpp new file mode 100644 index 000000000..beba1de78 --- /dev/null +++ b/libecole/tests/src/dynamics/test-branching-sum.cpp @@ -0,0 +1,129 @@ +#include +#include +#include +#include +#include + +#include +#include +#include + +#include "ecole/dynamics/branching-sum.hpp" +#include "ecole/exception.hpp" +#include "ecole/random.hpp" +#include "ecole/traits.hpp" + +#include "conftest.hpp" +#include "dynamics/unit-tests.hpp" + +using namespace ecole; + +TEST_CASE("BranchingSumDynamics unit tests with single branching", "[unit][dynamics]") { + auto const policy = [](auto const& action_set, auto const& /*model*/) { return std::array{action_set.value()[0]}; }; + dynamics::unit_tests(dynamics::BranchingSumDynamics{}, policy); +} + +/** + * A policy for multivariable branching. + * + * Try randomly to find two variables whose LP solution value sum is not integer. + * As it may not always exists, stop otherwise and branch on a single variable. + */ +struct MultiBranchingPolicy { + std::size_t n_multi = 0; + RandomEngine rng{0}; // NOLINT(cert-msc32-c, cert-msc51-cpp) We want reproducibility in tests + + auto operator()(trait::action_set_of_t const& action_set, scip::Model& model) + -> std::vector { + auto const& as = action_set.value(); + auto choice = std::uniform_int_distribution{0, as.size() - 1}; + auto indices = std::array{as[choice(rng)], as[choice(rng)]}; + auto const is_lp_sum_integral = [&model](auto idx1, auto idx2) { + auto const idx_to_var = [&](auto idx) { return SCIPcolGetVar(model.lp_columns()[idx]); }; + auto const is_lp = SCIPhasCurrentNodeLP(model.get_scip_ptr()); + auto const sum = SCIPvarGetSol(idx_to_var(idx1), is_lp) + SCIPvarGetSol(idx_to_var(idx2), is_lp); + return SCIPfeasFloor(model.get_scip_ptr(), sum) == SCIPfeasCeil(model.get_scip_ptr(), sum); + }; + auto constexpr n_trials = 10; + for (std::size_t n = 0; n < n_trials; ++n) { + if ((indices[0] != indices[1]) && (!is_lp_sum_integral(indices[0], indices[1]))) { + n_multi++; + return {indices[0], indices[1]}; + } + indices = {as[choice(rng)], as[choice(rng)]}; + } + return {indices[0]}; + } +}; + +TEST_CASE("BranchingSumDynamics unit tests with multi branching", "[unit][dynamics]") { + auto policy = MultiBranchingPolicy{}; + dynamics::unit_tests(dynamics::BranchingSumDynamics{}, policy); +} + +TEST_CASE("BranchingSumDynamics return valid action set", "[dynamics][slow]") { + auto dyn = dynamics::BranchingSumDynamics{}; + auto model = get_model(); + + SECTION("Return valid action set") { + auto const [done, action_set] = dyn.reset_dynamics(model); + REQUIRE(action_set.has_value()); + auto const& branch_cands = action_set.value(); + REQUIRE(branch_cands.size() > 0); + REQUIRE(branch_cands.size() <= model.lp_columns().size()); + REQUIRE(xt::all(branch_cands >= 0)); + REQUIRE(xt::all(branch_cands < model.lp_columns().size())); + REQUIRE(xt::unique(branch_cands).size() == branch_cands.size()); + } +} + +/** Compute the final dual bound of a Model solving a copy. */ +auto final_dual_bound(scip::Model const& model) { + auto model_copy = model.copy_orig(); + model_copy.solve(); + return model_copy.dual_bound(); +} + +TEST_CASE("BranchingSumDynamics can solve instance", "[dynamics][slow]") { + auto dyn = dynamics::BranchingSumDynamics{}; + auto model = get_model(); + auto const dual_bound = final_dual_bound(model); + + SECTION("Solve instance with multiple variables") { + auto policy = MultiBranchingPolicy{}; + auto [done, action_set] = dyn.reset_dynamics(model); + while (!done) { + REQUIRE(action_set.has_value()); + auto const action = policy(action_set, model); + std::tie(done, action_set) = dyn.step_dynamics(model, action); + } + REQUIRE(model.is_solved()); + REQUIRE(model.dual_bound() == Approx(dual_bound)); // Floating point approximation + REQUIRE(model.primal_bound() == Approx(dual_bound)); // Floating point approximation + } + + SECTION("Solve instance with single variables") { + auto [done, action_set] = dyn.reset_dynamics(model); + while (!done) { + REQUIRE(action_set.has_value()); + auto const action = action_set.value()[0]; + std::tie(done, action_set) = dyn.step_dynamics(model, {&action, 1}); + } + REQUIRE(model.is_solved()); + REQUIRE(model.dual_bound() == Approx(dual_bound)); // Floating point approximation + REQUIRE(model.primal_bound() == Approx(dual_bound)); // Floating point approximation + } +} + +TEST_CASE("BranchingSumDynamics handles invalid inputs", "[dynamics]") { + auto dyn = dynamics::BranchingSumDynamics{}; + auto model = get_model(); + + SECTION("Throw on invalid branching variable") { + auto const [done, action_set] = dyn.reset_dynamics(model); + REQUIRE_FALSE(done); + REQUIRE(action_set.has_value()); + auto const action = model.lp_columns().size() + 1; + REQUIRE_THROWS_AS(dyn.step_dynamics(model, {&action, 1}), std::invalid_argument); + } +} diff --git a/python/ecole/src/ecole/core/dynamics.cpp b/python/ecole/src/ecole/core/dynamics.cpp index d2cdde7b9..42222967f 100644 --- a/python/ecole/src/ecole/core/dynamics.cpp +++ b/python/ecole/src/ecole/core/dynamics.cpp @@ -4,6 +4,7 @@ #include #include +#include "ecole/dynamics/branching-sum.hpp" #include "ecole/dynamics/branching.hpp" #include "ecole/dynamics/configuring.hpp" #include "ecole/dynamics/primal-search.hpp" @@ -132,6 +133,90 @@ void bind_submodule(pybind11::module_ const& m) { )"); } + { + using idx_t = typename BranchingSumDynamics::Action::value_type; + dynamics_class{m, "BranchingSumDynamics", R"( + Sum of variables branching Dynamics. + + Based on a SCIP `branching callback `_ + with maximal priority and no depth limit. + The dynamics give the control back to the user every time the callback would be called. + The user receives as an action set the list of branching candidates, and is expected to select + a subset of them to branch on their sum. + + .. warning:: + The function used to perform branching is provided by Ecole and has not been extensively tested on + a large variety of problem instances. + )"} + .def_reset_dynamics(R"( + Start solving up to first branching node. + + Start solving with SCIP defaults (``SCIPsolve``) and give back control to the user on the + first branching decision. + Users can inherit from this dynamics to change the defaults settings such as presolving + and cutting planes. + + Parameters + ---------- + model: + The state of the Markov Decision Process. Passed by the environment. + + Returns + ------- + done: + Whether the instance is solved. + This can happen without branching, for instance if the instance is solved during presolving. + action_set: + List of branching candidates (``SCIPgetPseudoBranchCands``). + )") + .def_set_dynamics_random_state(R"( + Set seeds on the :py:class:`~ecole.scip.Model`. + + Set seed parameters, including permutation, LP, and shift. + + Parameters + ---------- + model: + The state of the Markov Decision Process. Passed by the environment. + random_engine: + The source of randomness. Passed by the environment. + )") + .def( + "step_dynamics", + [](BranchingSumDynamics& self, scip::Model& model, Numpy const& action) { + auto const vars = nonstd::span{action.data(), static_cast(action.size())}; + auto const release = py::gil_scoped_release{}; + return self.step_dynamics(model, vars); + }, + py::arg("model"), + py::arg("action"), + R"( + Branch and resume solving until next branching. + + Branching is done on the sum of the given variables using their LP or pseudo solution value. + To make a valid branching, the sum of the LP or pseudo solution value of the given variables + must be non integer. + In the opposite case, an exception will be thrown, + The control is given back to the user on the next branching decision or when done. + + Parameters + ---------- + model: + The state of the Markov Decision Process. Passed by the environment. + action: + A subset of the variables given in the action set. + Not all subsets are valid (see above). + + Returns + ------- + done: + Whether the instance is solved. + action_set: + List of branching candidates (``SCIPgetPseudoBranchCands``). + )") + .def(py::init<>()); + } + { dynamics_class{m, "ConfiguringDynamics", R"( Setting solving parameters Dynamics. diff --git a/python/ecole/src/ecole/environment.py b/python/ecole/src/ecole/environment.py index 4929423c3..6bb819b1f 100644 --- a/python/ecole/src/ecole/environment.py +++ b/python/ecole/src/ecole/environment.py @@ -213,6 +213,11 @@ class Branching(Environment): __DefaultObservationFunction__ = ecole.observation.NodeBipartite +class BranchingSum(Environment): + __Dynamics__ = ecole.dynamics.BranchingSumDynamics + __DefaultObservationFunction__ = ecole.observation.NodeBipartite + + class Configuring(Environment): __Dynamics__ = ecole.dynamics.ConfiguringDynamics diff --git a/python/ecole/tests/test_dynamics.py b/python/ecole/tests/test_dynamics.py index c9ffa452c..26edb7bdc 100644 --- a/python/ecole/tests/test_dynamics.py +++ b/python/ecole/tests/test_dynamics.py @@ -79,6 +79,32 @@ def setup_method(self, method): self.dynamics = ecole.dynamics.BranchingDynamics(True) +class TestBranchingSum_List(DynamicsUnitTests): + @staticmethod + def assert_action_set(action_set): + assert isinstance(action_set, np.ndarray) + assert action_set.ndim == 1 + assert action_set.size > 0 + assert action_set.dtype == np.uint64 + + @staticmethod + def policy(action_set): + return [action_set[0]] + + @staticmethod + def bad_policy(action_set): + return [1 << 31] + + def setup_method(self, method): + self.dynamics = ecole.dynamics.BranchingSumDynamics() + + +class TestBranchingSum_Numpy(TestBranchingSum_List): + @staticmethod + def policy(action_set): + return np.array([action_set[0]]) + + class TestConfiguring(DynamicsUnitTests): @staticmethod def assert_action_set(action_set):