import ecole.environment
import torch
import random
import numpy as np
import ecole
import os
import matplotlib.pyplot as plt
def main():
seed = 0
# For reproducibility
torch.manual_seed(seed)
random.seed(seed)
np.random.seed(seed)
generator = torch.Generator()
generator.manual_seed(seed)
n_rows = 500
n_cols = 1000
density = 0.1
n_instances = 1
# Define problem instance
gen = ecole.instance.SetCoverGenerator(
n_rows=n_rows,
n_cols=n_cols,
density=density,
)
gen.seed(seed)
# SCIP parameters used in paper
scip_parameters = {
"separating/maxrounds": 0,
"presolving/maxrestarts": 0,
"limits/time": 3600.0, # Time limit in seconds
}
# Define information function
information_function = {
"n_nodes": ecole.reward.NNodes(), # Total number of nodes processed since the previous state
"time": ecole.reward.SolvingTime(), # Number of seconds spent solving the instance since the previous state
"primal_int": ecole.reward.PrimalIntegral(),
"dual_int": ecole.reward.DualIntegral(),
"primal_dual_int": ecole.reward.PrimalDualIntegral()
}
# SCIP environments
sb_env = ecole.environment.Branching(
observation_function=ecole.observation.StrongBranchingScores(),
reward_function=None,
information_function=information_function,
scip_params=scip_parameters
)
sb_env.seed(seed)
pb_env = ecole.environment.Branching(
observation_function=ecole.observation.Pseudocosts(),
reward_function=None,
information_function=information_function,
scip_params=scip_parameters
)
pb_env.seed(seed)
def record(instance, env):
# Solve using benchmark
observation, action_set, _, done, info = env.reset(instance)
dual_int_arr = []
nodes = []
time = []
while not done:
action = action_set[observation[action_set].argmax()]
observation, action_set, _, done, info = env.step(action)
dual_int_arr.append(info['dual_int'])
nodes.append(info['n_nodes'])
time.append(info['time'])
# Get the optimal or best-found solution
scip_model = env.model.as_pyscipopt()
status = scip_model.getStatus()
flag = False
if status in ["optimal"]:
flag = True
obj_value = scip_model.getObjVal() # Optimal solution
else:
obj_value = scip_model.getPrimalbound() # Best solution found within time limit
return dual_int_arr, obj_value, nodes, time, flag
# Record results
for instance_i in range(n_instances):
instance = next(gen)
sb_dual, sb_obj, sb_nodes, sb_time, sb_is_optimal = record(instance, sb_env)
pb_dual, pb_obj, pb_nodes, pb_time, pb_is_optimal = record(instance, pb_env)
sb_time = np.cumsum(sb_time)
pb_time = np.cumsum(pb_time)
print(f"Strong Branching - Time: {sb_time[-1]:.2f}s, Objective: {sb_obj}, Optimal: {sb_is_optimal}")
print(f"Pseudocost Branching - Time: {pb_time[-1]:.2f}s, Objective: {pb_obj}, Optimal: {pb_is_optimal}")
fig, axis = plt.subplots(3, 1) # Larger figure for better readability
# Plot for reward dual integral
axis[0].plot(sb_dual, marker='^', linestyle='-', label='SB reward dual integral', alpha=0.6)
axis[0].plot(pb_dual, marker='o', linestyle='--', label='PB reward dual integral', alpha=0.7)
axis[0].set_title('Dual Integral')
axis[0].legend()
axis[0].grid(True)
# Plot for nodes processed
axis[1].plot(np.cumsum(sb_nodes), marker='^', linestyle='-', label='SB nodes', alpha=0.6)
axis[1].plot(np.cumsum(pb_nodes), marker='o', linestyle='--', label='PB nodes', alpha=0.7)
axis[1].set_title('Nodes Processed')
axis[1].legend()
axis[1].grid(True)
# Plot for solving time
axis[2].plot(sb_time, marker='^', linestyle='-', label='SB time', alpha=0.6)
axis[2].plot(pb_time, marker='o', linestyle='--', label='PB time', alpha=0.7)
axis[2].set_title('Solving Time')
axis[2].legend()
axis[2].grid(True)
# Adjust layout and display
plt.tight_layout()
plt.show()
return
if __name__ == "__main__":
main()