Implement GUB branching

docs/reference/environments.rst

@@ -16,6 +16,11 @@ Branching
 .. autoclass:: ecole.environment.Branching
 .. autoclass:: ecole.dynamics.BranchingDynamics
 
+BranchingGUB
+^^^^^^^^^^^^
+.. autoclass:: ecole.environment.BranchingGUB
+.. autoclass:: ecole.dynamics.BranchingGUBDynamics
+
 Configuring
 ^^^^^^^^^^^
 .. autoclass:: ecole.environment.Configuring

libecole/CMakeLists.txt

@@ -35,6 +35,7 @@ add_library(
 	src/observation/pseudocosts.cpp
 
 	src/dynamics/branching.cpp
+	src/dynamics/branching-gub.cpp
 	src/dynamics/configuring.cpp
 )
 set_target_properties(ecole-lib PROPERTIES OUTPUT_NAME ecole)

libecole/include/ecole/dynamics/branching-gub.hpp

@@ -0,0 +1,24 @@
+#pragma once
+
+#include <cstddef>
+#include <optional>
+
+#include <nonstd/span.hpp>
+#include <xtensor/xtensor.hpp>
+
+#include "ecole/dynamics/dynamics.hpp"
+
+namespace ecole::dynamics {
+
+class BranchingGUBDynamics :
+	public EnvironmentDynamics<nonstd::span<std::size_t const>, std::optional<xt::xtensor<std::size_t, 1>>> {
+public:
+	using Action = nonstd::span<std::size_t const>;
+	using ActionSet = std::optional<xt::xtensor<std::size_t, 1>>;
+
+	std::tuple<bool, ActionSet> reset_dynamics(scip::Model& model) override;
+
+	std::tuple<bool, ActionSet> step_dynamics(scip::Model& model, Action const& var_idx) override;
+};
+
+}  // namespace ecole::dynamics

libecole/include/ecole/dynamics/branching.hpp

@@ -19,7 +19,7 @@ class BranchingDynamics : public EnvironmentDynamics<std::size_t, std::optional<
 
 	std::tuple<bool, ActionSet> reset_dynamics(scip::Model& model) override;
 
-	std::tuple<bool, ActionSet> step_dynamics(scip::Model& model, std::size_t const& action) override;
+	std::tuple<bool, ActionSet> step_dynamics(scip::Model& model, std::size_t const& var_idx) override;
 };
 
 }  // namespace ecole::dynamics

libecole/include/ecole/environment/branching-gub.hpp

@@ -0,0 +1,18 @@
+#pragma once
+
+#include "ecole/dynamics/branching-gub.hpp"
+#include "ecole/environment/environment.hpp"
+#include "ecole/information/nothing.hpp"
+#include "ecole/observation/nodebipartite.hpp"
+#include "ecole/reward/isdone.hpp"
+
+namespace ecole::environment {
+
+template <
+	typename ObservationFunction = observation::NodeBipartite,
+	typename RewardFunction = reward::IsDone,
+	typename InformationFunction = information::Nothing>
+using BranchingGUB =
+	Environment<dynamics::BranchingGUBDynamics, ObservationFunction, RewardFunction, InformationFunction>;
+
+}  // namespace ecole::environment

libecole/include/ecole/scip/model.hpp

@@ -124,7 +124,7 @@ class Model {
 	[[nodiscard]] bool is_solved() const noexcept;
 
 	void solve_iter();
-	void solve_iter_branch(SCIP_VAR* var);
+	void solve_iter_branch(SCIP_RESULT result);
 	void solve_iter_stop();
 	[[nodiscard]] bool solve_iter_is_done();
 

libecole/include/ecole/scip/scimpl.hpp

@@ -2,6 +2,7 @@
 
 #include <memory>
 
+#include <nonstd/span.hpp>
 #include <scip/scip.h>
 
 #include "ecole/utility/reverse-control.hpp"
@@ -22,7 +23,7 @@ class Scimpl {
 	Scimpl copy_orig();
 
 	void solve_iter();
-	void solve_iter_branch(SCIP_VAR* var);
+	void solve_iter_branch(SCIP_RESULT result);
 	void solve_iter_stop();
 	bool solve_iter_is_done();
 

libecole/src/dynamics/branching-gub.cpp

@@ -0,0 +1,198 @@
+/***********************************************************
+ *  Implementation of SCIPbranchGUB (aim to merge in SCIP) *
+ ***********************************************************/
+
+#include <cassert>
+
+#include <fmt/format.h>
+
+#include <scip/cons_linear.h>
+#include <scip/debug.h>
+#include <scip/scip.h>
+#include <scip/scip_branch.h>
+#include <scip/scip_numerics.h>
+#include <scip/struct_scip.h>
+#include <scip/struct_var.h>
+
+namespace {
+
+SCIP_RETCODE SCIPbranchGUB_add_child(
+	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=*/FALSE,
+		/*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
+SCIPbranchGUB(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;
+	}
+
+	/* Check individual variables and compute the sum of their LP relaxation value */
+	SCIP_Real pseudo_sol_sum = 0.;
+	SCIP_VAR const* const* const vars_end = vars + nvars;
+	for (SCIP_VAR** var_iter = vars; var_iter < vars_end; ++var_iter) {
+		/* Validate that variables are integer, not fixed */
+		assert(SCIPvarIsActive(*var_iter));
+		assert(SCIPvarGetProbindex(*var_iter) >= 0);
+		if (SCIPvarGetType(*var_iter) == SCIP_VARTYPE_CONTINUOUS) {
+			SCIPerrorMessage("cannot branch on constraint containing continuous variable <%s>\n", SCIPvarGetName(*var_iter));
+			return SCIP_INVALIDDATA;
+		}
+		if (SCIPisEQ(scip, SCIPvarGetLbLocal(*var_iter), SCIPvarGetUbLocal(*var_iter))) {
+			SCIPerrorMessage(
+				"cannot branch on constraint containing variable <%s> with fixed domain [%.15g,%.15g]\n",
+				SCIPvarGetName(*var_iter),
+				SCIPvarGetLbLocal(*var_iter),
+				SCIPvarGetUbLocal(*var_iter));
+			return SCIP_INVALIDDATA;
+		}
+
+		SCIP_Real val = SCIPvarGetSol(*var_iter, SCIPhasCurrentNodeLP(scip));
+
+		/* avoid branching on infinite values in pseudo solution */
+		if (SCIPisInfinity(scip, -val) || SCIPisInfinity(scip, val)) {
+			SCIPerrorMessage("cannot branch on variables containing infinite values");
+			return SCIP_INVALIDDATA;
+		}
+		pseudo_sol_sum += val;
+	}
+
+	SCIP_Real const downbound = SCIPfeasFloor(scip, pseudo_sol_sum);
+	SCIP_Real const upbound = SCIPfeasCeil(scip, pseudo_sol_sum);
+	SCIP_Real const estimate = SCIPnodeGetLowerbound(SCIPgetCurrentNode(scip));
+	if (SCIPisEQ(scip, downbound, upbound)) {
+		SCIPerrorMessage("cannot branch on a variables whose sum of LP solution value is integer");
+		return SCIP_INVALIDDATA;
+	}
+
+	SCIP_Real const inf = SCIPinfinity(scip);
+	SCIP_Real* ones = nullptr;
+	SCIP_Retcode retcode = SCIP_OKAY;
+	SCIP_CALL(SCIPallocBufferArray(scip, &ones, nvars));
+	for (int i = 0; i < nvars; ++i) {
+		ones[i] = 1.;
+	}
+
+	SCIP_CALL_TERMINATE(
+		retcode, SCIPbranchGUB_add_child(scip, 1, estimate, vars, ones, nvars, -inf, downbound, downchild), TERM);
+	SCIP_CALL_TERMINATE(
+		retcode, SCIPbranchGUB_add_child(scip, 1, estimate, vars, ones, nvars, upbound, inf, upchild), TERM);
+
+TERM:
+	SCIPfreeBufferArray(scip, &ones);
+	return retcode;
+}
+
+}  // namespace
+
+/********************************************
+ *  Implementation of BranchingGUBDynamics  *
+ ********************************************/
+
+#include <algorithm>
+#include <stdexcept>
+#include <vector>
+
+#include <xtensor/xtensor.hpp>
+
+#include "ecole/dynamics/branching-gub.hpp"
+#include "ecole/scip/model.hpp"
+#include "ecole/scip/utils.hpp"
+
+namespace ecole::dynamics {
+
+namespace {
+
+std::optional<xt::xtensor<std::size_t, 1>> action_set(scip::Model const& model) {
+	if (model.get_stage() != SCIP_STAGE_SOLVING) {
+		return {};
+	}
+	auto const branch_cands = model.lp_branch_cands();
+	auto branch_cols = xt::xtensor<std::size_t, 1>::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 BranchingGUBDynamics::reset_dynamics(scip::Model& model) -> std::tuple<bool, ActionSet> {
+	model.solve_iter();
+	if (model.solve_iter_is_done()) {
+		return {true, {}};
+	}
+	return {false, action_set(model)};
+}
+
+auto BranchingGUBDynamics::step_dynamics(scip::Model& model, Action const& var_idx) -> std::tuple<bool, ActionSet> {
+	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_idx.begin(), var_idx.end(), is_out_of_bounds)) {
+		throw std::invalid_argument{"Branching index is larger than the number of columns."};
+	}
+
+	// Get variables associated with indices
+	auto vars = std::vector<SCIP_VAR*>(var_idx.size());
+	auto const idx_to_var = [lp_cols](auto idx) { return SCIPcolGetVar(lp_cols[idx]); };
+	std::transform(var_idx.begin(), var_idx.end(), vars.begin(), idx_to_var);
+
+	scip::call(SCIPbranchGUB, model.get_scip_ptr(), vars.data(), static_cast<int>(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

libecole/src/dynamics/branching.cpp

@@ -1,13 +1,10 @@
 #include <algorithm>
-#include <memory>
-#include <utility>
+#include <stdexcept>
 
 #include <xtensor/xtensor.hpp>
 
 #include "ecole/dynamics/branching.hpp"
-#include "ecole/exception.hpp"
 #include "ecole/scip/model.hpp"
-
 #include "ecole/scip/utils.hpp"
 
 namespace ecole::dynamics {
@@ -22,11 +19,8 @@ std::optional<xt::xtensor<std::size_t, 1>> action_set(scip::Model const& model,
 	}
 	auto const branch_cands = pseudo ? model.pseudo_branch_cands() : model.lp_branch_cands();
 	auto branch_cols = xt::xtensor<std::size_t, 1>::from_shape({branch_cands.size()});
-	std::transform(  //
-		branch_cands.begin(),
-		branch_cands.end(),
-		branch_cols.begin(),
-		[](auto const var) { return SCIPcolGetLPPos(SCIPvarGetCol(var)); });
+	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;
@@ -36,25 +30,25 @@ std::optional<xt::xtensor<std::size_t, 1>> action_set(scip::Model const& model,
 
 auto BranchingDynamics::reset_dynamics(scip::Model& model) -> std::tuple<bool, ActionSet> {
 	model.solve_iter();
-	auto const done = model.solve_iter_is_done();
-	if (done) {
-		return {done, {}};
+	if (model.solve_iter_is_done()) {
+		return {true, {}};
 	}
-	return {done, action_set(model, pseudo_candidates)};
+	return {false, action_set(model, pseudo_candidates)};
 }
 
-auto BranchingDynamics::step_dynamics(scip::Model& model, std::size_t const& action) -> std::tuple<bool, ActionSet> {
+auto BranchingDynamics::step_dynamics(scip::Model& model, std::size_t const& var_idx) -> std::tuple<bool, ActionSet> {
 	auto const lp_cols = model.lp_columns();
-	if (action >= lp_cols.size()) {
-		throw Exception{"Branching index is larger than the number of columns."};
+	if (var_idx >= lp_cols.size()) {
+		throw std::invalid_argument{"Branching index is larger than the number of columns."};
 	}
-	model.solve_iter_branch(SCIPcolGetVar(lp_cols[action]));
+	auto* const var = SCIPcolGetVar(lp_cols[var_idx]);
+	scip::call(SCIPbranchVar, model.get_scip_ptr(), var, nullptr, nullptr, nullptr);
+	model.solve_iter_branch(SCIP_BRANCHED);
 
-	auto const done = model.solve_iter_is_done();
-	if (done) {
-		return {done, {}};
+	if (model.solve_iter_is_done()) {
+		return {true, {}};
 	}
-	return {done, action_set(model, pseudo_candidates)};
+	return {false, action_set(model, pseudo_candidates)};
 }
 
 }  // namespace ecole::dynamics

libecole/src/scip/model.cpp

@@ -3,8 +3,8 @@
 #include <cstddef>
 #include <cstdlib>
 #include <cstring>
-#include <exception>
 #include <iterator>
+#include <stdexcept>
 #include <string>
 
 #include <fmt/format.h>
@@ -201,8 +201,8 @@ void Model::solve_iter() {
 	scimpl->solve_iter();
 }
 
-void Model::solve_iter_branch(SCIP_VAR* var) {
-	scimpl->solve_iter_branch(var);
+void Model::solve_iter_branch(SCIP_RESULT result) {
+	scimpl->solve_iter_branch(result);
 }
 
 void Model::solve_iter_stop() {