diff --git a/include/pyoptinterface/copt_model.hpp b/include/pyoptinterface/copt_model.hpp index 3196bd8..3af52ec 100644 --- a/include/pyoptinterface/copt_model.hpp +++ b/include/pyoptinterface/copt_model.hpp @@ -13,68 +13,74 @@ extern "C" int COPT_SearchParamAttr(copt_prob *prob, const char *name, int *p_type); } -#define APILIST \ - B(COPT_GetRetcodeMsg); \ - B(COPT_CreateProb); \ - B(COPT_DeleteProb); \ - B(COPT_WriteMps); \ - B(COPT_WriteLp); \ - B(COPT_WriteCbf); \ - B(COPT_WriteBin); \ - B(COPT_WriteBasis); \ - B(COPT_WriteSol); \ - B(COPT_WriteMst); \ - B(COPT_WriteParam); \ - B(COPT_AddCol); \ - B(COPT_DelCols); \ - B(COPT_AddRow); \ - B(COPT_AddQConstr); \ - B(COPT_AddSOSs); \ - B(COPT_AddCones); \ - B(COPT_DelRows); \ - B(COPT_DelQConstrs); \ - B(COPT_DelSOSs); \ - B(COPT_DelCones); \ - B(COPT_DelQuadObj); \ - B(COPT_ReplaceColObj); \ - B(COPT_SetObjConst); \ - B(COPT_SetObjSense); \ - B(COPT_SetQuadObj); \ - B(COPT_Solve); \ - B(COPT_SearchParamAttr); \ - B(COPT_SetIntParam); \ - B(COPT_SetDblParam); \ - B(COPT_GetIntParam); \ - B(COPT_GetDblParam); \ - B(COPT_GetIntAttr); \ - B(COPT_GetDblAttr); \ - B(COPT_GetColInfo); \ - B(COPT_GetColName); \ - B(COPT_SetColNames); \ - B(COPT_GetColType); \ - B(COPT_SetColType); \ - B(COPT_SetColLower); \ - B(COPT_SetColUpper); \ - B(COPT_GetRowInfo); \ - B(COPT_GetQConstrInfo); \ - B(COPT_GetRowName); \ - B(COPT_GetQConstrName); \ - B(COPT_SetRowNames); \ - B(COPT_SetQConstrNames); \ - B(COPT_AddMipStart); \ - B(COPT_GetQConstrRhs); \ - B(COPT_SetRowLower); \ - B(COPT_SetRowUpper); \ - B(COPT_SetQConstrRhs); \ - B(COPT_GetElem); \ - B(COPT_SetElem); \ - B(COPT_SetColObj); \ - B(COPT_GetBanner); \ - B(COPT_CreateEnv); \ - B(COPT_CreateEnvWithConfig); \ - B(COPT_DeleteEnv); \ - B(COPT_CreateEnvConfig); \ - B(COPT_DeleteEnvConfig); \ +#define APILIST \ + B(COPT_GetRetcodeMsg); \ + B(COPT_CreateProb); \ + B(COPT_DeleteProb); \ + B(COPT_WriteMps); \ + B(COPT_WriteLp); \ + B(COPT_WriteCbf); \ + B(COPT_WriteBin); \ + B(COPT_WriteBasis); \ + B(COPT_WriteSol); \ + B(COPT_WriteMst); \ + B(COPT_WriteParam); \ + B(COPT_AddCol); \ + B(COPT_DelCols); \ + B(COPT_AddRow); \ + B(COPT_AddQConstr); \ + B(COPT_AddSOSs); \ + B(COPT_AddCones); \ + B(COPT_DelRows); \ + B(COPT_DelQConstrs); \ + B(COPT_DelSOSs); \ + B(COPT_DelCones); \ + B(COPT_DelQuadObj); \ + B(COPT_ReplaceColObj); \ + B(COPT_SetObjConst); \ + B(COPT_SetObjSense); \ + B(COPT_SetQuadObj); \ + B(COPT_Solve); \ + B(COPT_SearchParamAttr); \ + B(COPT_SetIntParam); \ + B(COPT_SetDblParam); \ + B(COPT_GetIntParam); \ + B(COPT_GetDblParam); \ + B(COPT_GetIntAttr); \ + B(COPT_GetDblAttr); \ + B(COPT_GetColInfo); \ + B(COPT_GetColName); \ + B(COPT_SetColNames); \ + B(COPT_GetColType); \ + B(COPT_SetColType); \ + B(COPT_SetColLower); \ + B(COPT_SetColUpper); \ + B(COPT_GetRowInfo); \ + B(COPT_GetQConstrInfo); \ + B(COPT_GetRowName); \ + B(COPT_GetQConstrName); \ + B(COPT_SetRowNames); \ + B(COPT_SetQConstrNames); \ + B(COPT_AddMipStart); \ + B(COPT_GetQConstrRhs); \ + B(COPT_SetRowLower); \ + B(COPT_SetRowUpper); \ + B(COPT_SetQConstrRhs); \ + B(COPT_GetElem); \ + B(COPT_SetElem); \ + B(COPT_SetColObj); \ + B(COPT_GetBanner); \ + B(COPT_SetCallback); \ + B(COPT_GetCallbackInfo); \ + B(COPT_AddCallbackSolution); \ + B(COPT_AddCallbackLazyConstr); \ + B(COPT_AddCallbackUserCut); \ + B(COPT_Interrupt); \ + B(COPT_CreateEnv); \ + B(COPT_CreateEnvWithConfig); \ + B(COPT_DeleteEnv); \ + B(COPT_CreateEnvConfig); \ + B(COPT_DeleteEnvConfig); \ B(COPT_SetEnvConfig); namespace copt @@ -125,6 +131,27 @@ struct COPTfreemodelT }; }; +class COPTModel; +using COPTCallback = std::function; + +struct COPTCallbackUserdata +{ + void *model = nullptr; + COPTCallback callback; + int n_variables = 0; + int where = 0; + // store result of cbget + bool cb_get_mipsol_called = false; + std::vector mipsol; + bool cb_get_mipnoderel_called = false; + std::vector mipnoderel; + bool cb_get_mipincumbent_called = false; + std::vector mipincumbent; + // Cache for cbsolution + bool cb_solution_called = false; + std::vector heuristic_solution; +}; + class COPTModel { public: @@ -229,6 +256,32 @@ class COPTModel int _constraint_index(const ConstraintIndex &constraint); int _checked_constraint_index(const ConstraintIndex &constraint); + // Callback + void set_callback(const COPTCallback &callback, int cbctx); + + // For callback + bool has_callback = false; + void *m_cbdata = nullptr; + COPTCallbackUserdata m_callback_userdata; + + int cb_get_info_int(const std::string &what); + double cb_get_info_double(const std::string &what); + void cb_get_info_doublearray(const std::string &what); + + double cb_get_solution(const VariableIndex &variable); + double cb_get_relaxation(const VariableIndex &variable); + double cb_get_incumbent(const VariableIndex &variable); + void cb_set_solution(const VariableIndex &variable, double value); + double cb_submit_solution(); + + void cb_exit(); + + void cb_add_lazy_constraint(const ScalarAffineFunction &function, ConstraintSense sense, + CoeffT rhs); + void cb_add_lazy_constraint(const ExprBuilder &function, ConstraintSense sense, CoeffT rhs); + void cb_add_user_cut(const ScalarAffineFunction &function, ConstraintSense sense, CoeffT rhs); + void cb_add_user_cut(const ExprBuilder &function, ConstraintSense sense, CoeffT rhs); + private: MonotoneIndexer m_variable_index; diff --git a/include/pyoptinterface/gurobi_model.hpp b/include/pyoptinterface/gurobi_model.hpp index 61ffd67..f98c02a 100644 --- a/include/pyoptinterface/gurobi_model.hpp +++ b/include/pyoptinterface/gurobi_model.hpp @@ -266,10 +266,12 @@ class GurobiModel void cb_get_info_doublearray(int what); double cb_get_solution(const VariableIndex &variable); - double cb_get_noderel(const VariableIndex &variable); + double cb_get_relaxation(const VariableIndex &variable); void cb_set_solution(const VariableIndex &variable, double value); double cb_submit_solution(); + void cb_exit(); + void cb_add_lazy_constraint(const ScalarAffineFunction &function, ConstraintSense sense, CoeffT rhs); void cb_add_lazy_constraint(const ExprBuilder &function, ConstraintSense sense, CoeffT rhs); diff --git a/lib/copt_model.cpp b/lib/copt_model.cpp index 9ad0a91..9fcb2d0 100644 --- a/lib/copt_model.cpp +++ b/lib/copt_model.cpp @@ -538,6 +538,11 @@ void COPTModel::set_objective(const ExprBuilder &function, ObjectiveSense sense) void COPTModel::optimize() { + if (has_callback) + { + // Store the number of variables for the callback + m_callback_userdata.n_variables = get_raw_attribute_int("Cols"); + } int error = copt::COPT_Solve(m_model.get()); check_error(error); } @@ -978,3 +983,187 @@ void COPTEnvConfig::set(const char *param_name, const char *value) int error = copt::COPT_SetEnvConfig(m_config, param_name, value); check_error(error); } + +// Callback +int RealCOPTCallbackFunction(copt_prob *prob, void *cbdata, int cbctx, void *userdata) +{ + auto real_userdata = static_cast(userdata); + auto model = static_cast(real_userdata->model); + auto &callback = real_userdata->callback; + + model->m_cbdata = cbdata; + model->m_callback_userdata.where = cbctx; + model->m_callback_userdata.cb_get_mipsol_called = false; + model->m_callback_userdata.cb_get_mipnoderel_called = false; + model->m_callback_userdata.cb_get_mipincumbent_called = false; + model->m_callback_userdata.cb_solution_called = false; + callback(model, cbctx); + + return COPT_RETCODE_OK; +} + +void COPTModel::set_callback(const COPTCallback &callback, int cbctx) +{ + m_callback_userdata.model = this; + m_callback_userdata.callback = callback; + + int error = copt::COPT_SetCallback(m_model.get(), RealCOPTCallbackFunction, cbctx, + &m_callback_userdata); + check_error(error); + + has_callback = true; +} + +int COPTModel::cb_get_info_int(const std::string &what) +{ + int retval; + int error = copt::COPT_GetCallbackInfo(m_cbdata, what.c_str(), &retval); + check_error(error); + return retval; +} + +double COPTModel::cb_get_info_double(const std::string &what) +{ + double retval; + int error = copt::COPT_GetCallbackInfo(m_cbdata, what.c_str(), &retval); + check_error(error); + return retval; +} + +void COPTModel::cb_get_info_doublearray(const std::string &what) +{ + int n_vars = m_callback_userdata.n_variables; + double *val = nullptr; + if (what == COPT_CBINFO_MIPCANDIDATE) + { + m_callback_userdata.mipsol.resize(n_vars); + val = m_callback_userdata.mipsol.data(); + } + else if (what == COPT_CBINFO_RELAXSOLUTION) + { + m_callback_userdata.mipnoderel.resize(n_vars); + val = m_callback_userdata.mipnoderel.data(); + } + else if (what == COPT_CBINFO_INCUMBENT) + { + m_callback_userdata.mipincumbent.resize(n_vars); + val = m_callback_userdata.mipincumbent.data(); + } + else + { + throw std::runtime_error("Invalid what for cb_get_info_doublearray"); + } + int error = copt::COPT_GetCallbackInfo(m_cbdata, what.c_str(), val); + check_error(error); +} + +double COPTModel::cb_get_solution(const VariableIndex &variable) +{ + auto &userdata = m_callback_userdata; + if (!userdata.cb_get_mipsol_called) + { + cb_get_info_doublearray(COPT_CBINFO_MIPCANDIDATE); + userdata.cb_get_mipsol_called = true; + } + auto index = _variable_index(variable); + return userdata.mipsol[index]; +} + +double COPTModel::cb_get_relaxation(const VariableIndex &variable) +{ + auto &userdata = m_callback_userdata; + if (!userdata.cb_get_mipnoderel_called) + { + cb_get_info_doublearray(COPT_CBINFO_RELAXSOLUTION); + userdata.cb_get_mipnoderel_called = true; + } + auto index = _variable_index(variable); + return userdata.mipnoderel[index]; +} + +double COPTModel::cb_get_incumbent(const VariableIndex &variable) +{ + auto &userdata = m_callback_userdata; + if (!userdata.cb_get_mipincumbent_called) + { + cb_get_info_doublearray(COPT_CBINFO_INCUMBENT); + userdata.cb_get_mipincumbent_called = true; + } + auto index = _variable_index(variable); + return userdata.mipincumbent[index]; +} + +void COPTModel::cb_set_solution(const VariableIndex &variable, double value) +{ + auto &userdata = m_callback_userdata; + if (!userdata.cb_solution_called) + { + userdata.heuristic_solution.resize(userdata.n_variables, COPT_UNDEFINED); + userdata.cb_solution_called = true; + } + userdata.heuristic_solution[_variable_index(variable)] = value; +} + +double COPTModel::cb_submit_solution() +{ + if (!m_callback_userdata.cb_solution_called) + { + throw std::runtime_error("No solution is set in the callback!"); + } + double obj; + int error = copt::COPT_AddCallbackSolution(m_cbdata, + m_callback_userdata.heuristic_solution.data(), &obj); + check_error(error); + return obj; +} + +void COPTModel::cb_exit() +{ + int error = copt::COPT_Interrupt(m_model.get()); + check_error(error); +} + +void COPTModel::cb_add_lazy_constraint(const ScalarAffineFunction &function, ConstraintSense sense, + CoeffT rhs) +{ + AffineFunctionPtrForm ptr_form; + ptr_form.make(this, function); + + int numnz = ptr_form.numnz; + int *cind = ptr_form.index; + double *cval = ptr_form.value; + char g_sense = copt_con_sense(sense); + double g_rhs = rhs - function.constant.value_or(0.0); + + int error = copt::COPT_AddCallbackLazyConstr(m_cbdata, numnz, cind, cval, g_sense, g_rhs); + check_error(error); +} + +void COPTModel::cb_add_lazy_constraint(const ExprBuilder &function, ConstraintSense sense, + CoeffT rhs) +{ + ScalarAffineFunction f(function); + cb_add_lazy_constraint(f, sense, rhs); +} + +void COPTModel::cb_add_user_cut(const ScalarAffineFunction &function, ConstraintSense sense, + CoeffT rhs) +{ + AffineFunctionPtrForm ptr_form; + ptr_form.make(this, function); + + int numnz = ptr_form.numnz; + int *cind = ptr_form.index; + double *cval = ptr_form.value; + char g_sense = copt_con_sense(sense); + double g_rhs = rhs - function.constant.value_or(0.0); + + int error = copt::COPT_AddCallbackUserCut(m_cbdata, numnz, cind, cval, g_sense, g_rhs); + check_error(error); +} + +void COPTModel::cb_add_user_cut(const ExprBuilder &function, ConstraintSense sense, CoeffT rhs) +{ + ScalarAffineFunction f(function); + cb_add_user_cut(f, sense, rhs); +} diff --git a/lib/copt_model_ext.cpp b/lib/copt_model_ext.cpp index d81815a..503ac03 100644 --- a/lib/copt_model_ext.cpp +++ b/lib/copt_model_ext.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include "pyoptinterface/copt_model.hpp" @@ -108,11 +109,39 @@ NB_MODULE(copt_model_ext, m) nb::overload_cast(&COPTModelMixin::set_objective_as_constant), nb::arg("expr"), nb::arg("sense") = ObjectiveSense::Minimize) + .def("cb_add_lazy_constraint", + nb::overload_cast( + &COPTModelMixin::cb_add_lazy_constraint), + nb::arg("expr"), nb::arg("sense"), nb::arg("rhs")) + .def("cb_add_lazy_constraint", + nb::overload_cast( + &COPTModelMixin::cb_add_lazy_constraint), + nb::arg("expr"), nb::arg("sense"), nb::arg("rhs")) + .def("cb_add_user_cut", + nb::overload_cast( + &COPTModelMixin::cb_add_user_cut), + nb::arg("expr"), nb::arg("sense"), nb::arg("rhs")) + .def("cb_add_user_cut", + nb::overload_cast( + &COPTModelMixin::cb_add_user_cut), + nb::arg("expr"), nb::arg("sense"), nb::arg("rhs")) + + .def("optimize", &COPTModelMixin::optimize, nb::call_guard()) + // clang-format off - BIND_F(optimize) BIND_F(version_string) BIND_F(get_raw_model) + BIND_F(set_callback) + BIND_F(cb_get_info_int) + BIND_F(cb_get_info_double) + BIND_F(cb_get_solution) + BIND_F(cb_get_relaxation) + BIND_F(cb_get_incumbent) + BIND_F(cb_set_solution) + BIND_F(cb_submit_solution) + BIND_F(cb_exit) + BIND_F(raw_parameter_attribute_type) BIND_F(set_raw_parameter_int) diff --git a/lib/copt_model_ext_constants.cpp b/lib/copt_model_ext_constants.cpp index f568eaf..09713ee 100644 --- a/lib/copt_model_ext_constants.cpp +++ b/lib/copt_model_ext_constants.cpp @@ -51,7 +51,18 @@ void bind_copt_constants(nb::module_ &m) COPT.attr("IMPRECISE") = COPT_LPSTATUS_IMPRECISE; COPT.attr("INTERRUPTED") = COPT_MIPSTATUS_INTERRUPTED; + COPT.attr("CBCONTEXT_INCUMBENT") = COPT_CBCONTEXT_INCUMBENT; COPT.attr("CBCONTEXT_MIPRELAX") = COPT_CBCONTEXT_MIPRELAX; COPT.attr("CBCONTEXT_MIPSOL") = COPT_CBCONTEXT_MIPSOL; COPT.attr("CBCONTEXT_MIPNODE") = COPT_CBCONTEXT_MIPNODE; + + COPT.attr("CBINFO_BESTOBJ") = COPT_CBINFO_BESTOBJ; + COPT.attr("CBINFO_BESTBND") = COPT_CBINFO_BESTBND; + COPT.attr("CBINFO_HASINCUMBENT") = COPT_CBINFO_HASINCUMBENT; + COPT.attr("CBINFO_INCUMBENT") = COPT_CBINFO_INCUMBENT; + COPT.attr("CBINFO_MIPCANDIDATE") = COPT_CBINFO_MIPCANDIDATE; + COPT.attr("CBINFO_MIPCANDOBJ") = COPT_CBINFO_MIPCANDOBJ; + COPT.attr("CBINFO_RELAXSOLUTION") = COPT_CBINFO_RELAXSOLUTION; + COPT.attr("CBINFO_RELAXSOLOBJ") = COPT_CBINFO_RELAXSOLOBJ; + COPT.attr("CBINFO_NODESTATUS") = COPT_CBINFO_NODESTATUS; } diff --git a/lib/gurobi_model.cpp b/lib/gurobi_model.cpp index 745c2c8..c320a5d 100644 --- a/lib/gurobi_model.cpp +++ b/lib/gurobi_model.cpp @@ -1108,7 +1108,7 @@ double GurobiModel::cb_get_solution(const VariableIndex &variable) return userdata.mipsol[index]; } -double GurobiModel::cb_get_noderel(const VariableIndex &variable) +double GurobiModel::cb_get_relaxation(const VariableIndex &variable) { auto &userdata = m_callback_userdata; if (!userdata.cb_get_mipnoderel_called) @@ -1160,6 +1160,11 @@ void GurobiModel::cb_add_lazy_constraint(const ScalarAffineFunction &function, check_error(error); } +void GurobiModel::cb_exit() +{ + gurobi::GRBterminate(m_model.get()); +} + void GurobiModel::cb_add_lazy_constraint(const ExprBuilder &function, ConstraintSense sense, CoeffT rhs) { diff --git a/lib/gurobi_model_ext.cpp b/lib/gurobi_model_ext.cpp index eea6d70..133eef7 100644 --- a/lib/gurobi_model_ext.cpp +++ b/lib/gurobi_model_ext.cpp @@ -140,8 +140,9 @@ NB_MODULE(gurobi_model_ext, m) &GurobiModelMixin::cb_add_user_cut), nb::arg("expr"), nb::arg("sense"), nb::arg("rhs")) + .def("optimize", &GurobiModelMixin::optimize, nb::call_guard()) + // clang-format off - BIND_F(optimize) BIND_F(update) BIND_F(version_string) BIND_F(get_raw_model) @@ -150,9 +151,10 @@ NB_MODULE(gurobi_model_ext, m) BIND_F(cb_get_info_int) BIND_F(cb_get_info_double) BIND_F(cb_get_solution) - BIND_F(cb_get_noderel) + BIND_F(cb_get_relaxation) BIND_F(cb_set_solution) BIND_F(cb_submit_solution) + BIND_F(cb_exit) BIND_F(raw_parameter_type) BIND_F(set_raw_parameter_int) diff --git a/lib/highs_model_ext.cpp b/lib/highs_model_ext.cpp index 7595028..25d2889 100644 --- a/lib/highs_model_ext.cpp +++ b/lib/highs_model_ext.cpp @@ -102,8 +102,9 @@ NB_MODULE(highs_model_ext, m) nb::overload_cast(&HighsModelMixin::set_objective_as_constant), nb::arg("expr"), nb::arg("sense") = ObjectiveSense::Minimize) + .def("optimize", &HighsModelMixin::optimize, nb::call_guard()) + // clang-format off - BIND_F(optimize) BIND_F(version_string) BIND_F(get_raw_model) diff --git a/lib/mosek_model_ext.cpp b/lib/mosek_model_ext.cpp index 463cb31..c5c0030 100644 --- a/lib/mosek_model_ext.cpp +++ b/lib/mosek_model_ext.cpp @@ -104,8 +104,9 @@ NB_MODULE(mosek_model_ext, m) nb::overload_cast(&MOSEKModelMixin::set_objective_as_constant), nb::arg("expr"), nb::arg("sense") = ObjectiveSense::Minimize) + .def("optimize", &MOSEKModelMixin::optimize, nb::call_guard()) + // clang-format off - BIND_F(optimize) BIND_F(version_string) BIND_F(get_raw_model) diff --git a/tests/callback_tsp.py b/tests/tsp_cb.py similarity index 53% rename from tests/callback_tsp.py rename to tests/tsp_cb.py index f801bd4..4fe3d2f 100644 --- a/tests/callback_tsp.py +++ b/tests/tsp_cb.py @@ -1,17 +1,10 @@ # This file is adapted from the examples/python/tsp.py in Gurobi installation. -# We use this file to ensure our callback implementation is correct and the result is compared with gurobipy +# We use this file to ensure our callback implementation is correct and the result is compared with gurobipy/coptpy # this test is currently run manually # Copyright 2024, Gurobi Optimization, LLC -# Solve a traveling salesman problem on a randomly generated set of points -# using lazy constraints. The base MIP model only includes 'degree-2' -# constraints, requiring each node to have exactly two incident edges. -# Solutions to this model may contain subtours - tours that don't visit every -# city. The lazy constraint callback adds new constraints to cut them off. - -import sys -import logging +from typing import List, Tuple import math import random import time @@ -19,18 +12,16 @@ from itertools import combinations import pyoptinterface as poi -from pyoptinterface import gurobi +from pyoptinterface import gurobi, copt import gurobipy as gp from gurobipy import GRB +import coptpy as cp +from coptpy import COPT -def shortest_subtour(edges): - """Given a list of edges, return the shortest subtour (as a list of nodes) - found by following those edges. It is assumed there is exactly one 'in' - edge and one 'out' edge for every node represented in the edge list.""" - # Create a mapping from each node to its neighbours +def shortest_subtour(edges: List[Tuple[int, int]]) -> List[int]: node_neighbors = defaultdict(list) for i, j in edges: node_neighbors[i].append(j) @@ -55,59 +46,16 @@ def shortest_subtour(edges): return shortest -class TSPCallback: - """Callback class implementing lazy constraints for the TSP. At MIPSOL - callbacks, solutions are checked for subtours and subtour elimination - constraints are added if needed.""" - +class GurobiTSPCallback: def __init__(self, nodes, x): self.nodes = nodes self.x = x - def run_poi(self, model, where): - """Callback entry point: call lazy constraints routine when new - solutions are found. Stop the optimization if there is an exception in - user code.""" - if where == GRB.Callback.MIPSOL: - try: - self.eliminate_subtours_poi(model) - except Exception: - logging.exception("Exception occurred in MIPSOL callback") - exit() - - def eliminate_subtours_poi(self, model): - """Extract the current solution, check for subtours, and formulate lazy - constraints to cut off the current solution if subtours are found. - Assumes we are at MIPSOL.""" - edges = [] - for (i, j), xij in self.x.items(): - v = model.cb_get_solution(xij) - if v > 0.5: - edges.append((i, j)) - tour = shortest_subtour(edges) - if len(tour) < len(self.nodes): - # add subtour elimination constraint for every pair of cities in tour - model.cb_add_lazy_constraint( - poi.quicksum(self.x[i, j] for i, j in combinations(tour, 2)), - poi.Leq, - len(tour) - 1, - ) - def __call__(self, model, where): - """Callback entry point: call lazy constraints routine when new - solutions are found. Stop the optimization if there is an exception in - user code.""" if where == GRB.Callback.MIPSOL: - try: - self.eliminate_subtours_gurobipy(model) - except Exception: - logging.exception("Exception occurred in MIPSOL callback") - exit() + self.eliminate_subtours_gurobipy(model) def eliminate_subtours_gurobipy(self, model): - """Extract the current solution, check for subtours, and formulate lazy - constraints to cut off the current solution if subtours are found. - Assumes we are at MIPSOL.""" values = model.cbGetSolution(self.x) edges = [(i, j) for (i, j), v in values.items() if v > 0.5] tour = shortest_subtour(edges) @@ -119,7 +67,7 @@ def eliminate_subtours_gurobipy(self, model): ) -def solve_tsp_poi(nodes, distances): +def solve_tsp_gurobipy(nodes, distances): """ Solve a dense symmetric TSP using the following base formulation: @@ -130,26 +78,75 @@ def solve_tsp_poi(nodes, distances): and subtours eliminated using lazy constraints. """ - m = gurobi.Model() + m = gp.Model() + + x = m.addVars(distances.keys(), obj=distances, vtype=GRB.BINARY, name="e") + x.update({(j, i): v for (i, j), v in x.items()}) + + # Create degree 2 constraints + for i in nodes: + m.addConstr(gp.quicksum(x[i, j] for j in nodes if i != j) == 2) + + m.Params.OutputFlag = 0 + m.Params.LazyConstraints = 1 + cb = GurobiTSPCallback(nodes, x) + m.optimize(cb) + + edges = [(i, j) for (i, j), v in x.items() if v.X > 0.5] + tour = shortest_subtour(edges) + assert set(tour) == set(nodes) + + return tour, m.ObjVal + + +class POITSPCallback: + def __init__(self, nodes, x): + self.nodes = nodes + self.x = x + + def run_gurobi(self, model, where): + if where == GRB.Callback.MIPSOL: + self.eliminate_subtours_poi(model) + + def run_copt(self, model, where): + if where == COPT.CBCONTEXT_MIPSOL: + self.eliminate_subtours_poi(model) + + def eliminate_subtours_poi(self, model): + edges = [] + for (i, j), xij in self.x.items(): + v = model.cb_get_solution(xij) + if v > 0.5: + edges.append((i, j)) + tour = shortest_subtour(edges) + if len(tour) < len(self.nodes): + # add subtour elimination constraint for every pair of cities in tour + model.cb_add_lazy_constraint( + poi.quicksum(self.x[i, j] for i, j in combinations(tour, 2)), + poi.Leq, + len(tour) - 1, + ) + - # Create variables, and add symmetric keys to the resulting dictionary - # 'x', such that (i, j) and (j, i) refer to the same variable. +def solve_tsp_poi(f, nodes, distances): + m = f() x = m.add_variables(distances.keys(), domain=poi.VariableDomain.Binary, name="e") m.set_objective(poi.quicksum(distances[k] * x[k] for k in distances)) for i, j in distances.keys(): x[j, i] = x[i, j] - # Create degree 2 constraints for i in nodes: m.add_linear_constraint( poi.quicksum(x[i, j] for j in nodes if i != j), poi.Eq, 2 ) - # Optimize model using lazy constraints to eliminate subtours - m.set_raw_parameter("OutputFlag", 0) - m.set_raw_parameter("LazyConstraints", 1) - cb = TSPCallback(nodes, x) - m.set_callback(cb.run_poi) + m.set_model_attribute(poi.ModelAttribute.Silent, True) + cb = POITSPCallback(nodes, x) + if isinstance(m, gurobi.Model): + m.set_raw_parameter("LazyConstraints", 1) + m.set_callback(cb.run_gurobi) + elif isinstance(m, copt.Model): + m.set_callback(cb.run_copt, COPT.CBCONTEXT_MIPSOL) m.optimize() # Extract the solution as a tour @@ -162,39 +159,52 @@ def solve_tsp_poi(nodes, distances): return tour, objval -def solve_tsp_gurobipy(nodes, distances): - """ - Solve a dense symmetric TSP using the following base formulation: +class COPTTSPCallback(cp.CallbackBase): + def __init__(self, nodes, x): + super().__init__() + self.nodes = nodes + self.x = x - min sum_ij d_ij x_ij - s.t. sum_j x_ij == 2 forall i in V - x_ij binary forall (i,j) in E + def callback(self): + if self.where() == COPT.CBCONTEXT_MIPSOL: + self.eliminate_subtours_coptpy() - and subtours eliminated using lazy constraints. - """ + def eliminate_subtours_coptpy(self): + values = self.getSolution(self.x) + edges = [(i, j) for (i, j), v in values.items() if v > 0.5] + tour = shortest_subtour(edges) + if len(tour) < len(self.nodes): + # add subtour elimination constraint for every pair of cities in tour + self.addLazyConstr( + cp.quicksum(self.x[i, j] for i, j in combinations(tour, 2)) + <= len(tour) - 1 + ) - with gp.Env() as env, gp.Model(env=env) as m: - # Create variables, and add symmetric keys to the resulting dictionary - # 'x', such that (i, j) and (j, i) refer to the same variable. - x = m.addVars(distances.keys(), obj=distances, vtype=GRB.BINARY, name="e") - x.update({(j, i): v for (i, j), v in x.items()}) - # Create degree 2 constraints - for i in nodes: - m.addConstr(gp.quicksum(x[i, j] for j in nodes if i != j) == 2) +def solve_tsp_coptpy(nodes, distances): + env = cp.Envr() + m = env.createModel("TSP Callback Example") - # Optimize model using lazy constraints to eliminate subtours - m.Params.OutputFlag = 0 - m.Params.LazyConstraints = 1 - cb = TSPCallback(nodes, x) - m.optimize(cb) + x = m.addVars(distances.keys(), vtype=COPT.BINARY, nameprefix="e") + for (i, j), v in x.items(): + v.setInfo(COPT.Info.Obj, distances[i, j]) + for i, j in distances.keys(): + x[j, i] = x[i, j] - # Extract the solution as a tour - edges = [(i, j) for (i, j), v in x.items() if v.X > 0.5] - tour = shortest_subtour(edges) - assert set(tour) == set(nodes) + # Create degree 2 constraints + for i in nodes: + m.addConstr(cp.quicksum(x[i, j] for j in nodes if i != j) == 2) - return tour, m.ObjVal + m.Param.Logging = 0 + cb = COPTTSPCallback(nodes, x) + m.setCallback(cb, COPT.CBCONTEXT_MIPSOL) + m.solve() + + edges = [(i, j) for (i, j), v in x.items() if v.x > 0.5] + tour = shortest_subtour(edges) + assert set(tour) == set(nodes) + + return tour, m.objval def create_map(npoints, seed): @@ -212,10 +222,8 @@ def create_map(npoints, seed): return nodes, distances -if __name__ == "__main__": - seed = 87651234 - - for npoints in range(10, 80, 10): +def test_gurobi(npoints_series, seed): + for npoints in npoints_series: nodes, distances = create_map(npoints, seed) print(f"npoints = {npoints}") @@ -226,9 +234,40 @@ def create_map(npoints, seed): print(f"\t gurobipy time: {t1 - t0:g} seconds") t0 = time.time() - tour2, cost2 = solve_tsp_poi(nodes, distances) + f = gurobi.Model + tour2, cost2 = solve_tsp_poi(f, nodes, distances) t1 = time.time() print(f"\t poi time: {t1 - t0:g} seconds") assert tour1 == tour2 assert abs(cost1 - cost2) < 1e-6 + + +def test_copt(npoints_series, seed): + for npoints in npoints_series: + nodes, distances = create_map(npoints, seed) + + print(f"npoints = {npoints}") + + t0 = time.time() + tour1, cost1 = solve_tsp_coptpy(nodes, distances) + t1 = time.time() + print(f"\t coptpy time: {t1 - t0:g} seconds") + + t0 = time.time() + f = copt.Model + tour2, cost2 = solve_tsp_poi(f, nodes, distances) + t1 = time.time() + print(f"\t poi time: {t1 - t0:g} seconds") + + assert tour1 == tour2 + assert abs(cost1 - cost2) < 1e-6 + + +if __name__ == "__main__": + seed = 987651234 + + X = range(10, 90, 20) + + test_copt(X, seed) + test_gurobi(X, seed)