STSP Subtour Elimination with half-done Callback Interface

stsp_subtour_elimination.jl

using CPLEX
using JuMP
using Test

# TSPLIB could be more clear on how to compute the distances.
function euclidean_distance(x1, y1, x2, y2)
  round(√((x1 - x2)^2 + (y1 - y2)^2))
end

function get_instance()
  # See: https://people.sc.fsu.edu/~jburkardt/datasets/tsp/tsp.html
  # The optimal value is 33523 but we get 33522 because of rounding.
  att48 = [
    (6734, 1453), 
    (2233,   10),
    (5530, 1424),
    ( 401,  841),
    (3082, 1644),
    (7608, 4458),
    (7573, 3716),
    (7265, 1268),
    (6898, 1885),
    (1112, 2049),
    (5468, 2606),
    (5989, 2873),
    (4706, 2674),
    (4612, 2035),
    (6347, 2683),
    (6107,  669),
    (7611, 5184),
    (7462, 3590),
    (7732, 4723),
    (5900, 3561),
    (4483, 3369),
    (6101, 1110),
    (5199, 2182),
    (1633, 2809),
    (4307, 2322),
    ( 675, 1006),
    (7555, 4819),
    (7541, 3981),
    (3177,  756),
    (7352, 4506),
    (7545, 2801),
    (3245, 3305),
    (6426, 3173),
    (4608, 1198),
    (  23, 2216),
    (7248, 3779),
    (7762, 4595),
    (7392, 2244),
    (3484, 2829),
    (6271, 2135),
    (4985,  140),
    (1916, 1569),
    (7280, 4899),
    (7509, 3239),
    (  10, 2676),
    (6807, 2993),
    (5185, 3258),
    (3023, 1942)
  ]
  n = length(att48)
  d = zeros(n, n)
  for i = 1:n, j = (i+1):n
    d[i, j] = euclidean_distance(
      att48[i][1], att48[i][2], att48[j][1], att48[j][2]
    )
  end

  return n, d
end

# Based on: http://examples.gurobi.com/traveling-salesman-problem/
# n: number of nodes/vertexes
# edges: the current set of edges in the to-be-accepted-or-rejected
#   solution; each edge needs to pertain to exacly one cycle, there may
#   be one or more cycles.
# return: an array of all cycles, each cycle is an array of the vertexes
#   in the order they are visited (starting/finishing in an arbitrary vertex
#   of the cycle)
function tours(
  n :: N,
  edges :: Vector{Tuple{N, N}}
) :: Vector{Vector{N}} where {N}
  # visited: each node is processed a single time; visited[i] is true if
  # node i has been already processed (assigned to a cycle).
  visited = fill(false, n)
  # cycles: the list of tours that is returned at the end of the method.
  cycles = Vector{Vector{N}}()
  # linked: for every node, an array of the other two nodes that connect
  # to them. Ex.: if there is a cycle 1 -> 2 -> 3 -> 4 -> 1, then linked[3]
  # will have a vector [2, 4] (i.e., the neighbors of 3).
  linked = [Vector{N}() for i in 1:n]
  # The loop below just initizalize 'linked' to be as described above.
  for (a, b) in edges
    push!(linked[b], a)
    push!(linked[a], b)
  end

  # This outer loop only finishes when every node of the graph was already
  # visited. Note that this is not when all the nodes of the first cycle
  # are visited, but all nodes in the whole graph (all cycles).
  while true
    # Gets the first node that was not yet visited (from the whole graph).
    current = findfirst(i -> !visited[i], 1:n)
    # Create a cycle with this node.
    thiscycle = [current]
    # The inner loop builds a single cycle and then finishes. It is started
    # one time for each cycle, and iterate the number of nodes in that cycle.
    while true
      # The current node is always already inserted in the cycle, and now that
      # their neighbors will be checked it can be marked as visited to never
      # be considered again.
      visited[current] = true
      # From the two neighbors of the current link we get only the not visited.
      unvisited_neighbors = [v for v in linked[current] if !visited[v]]
      # If all neighbor were visited, then this means we have already iterated
      # the whole cycle and ended up in the first node visited (the one found
      # by findfirst before this loop). This loop may be ended.
      isempty(unvisited_neighbors) && break
      # If there is yet a unvisited neighbor, then we can keep walking the
      # cycle by going to its direction (i.e., having it as the next current
      # node). We may also already register this new current node as
      # pertaining to the same cycle.
      current = first(unvisited_neighbors)
      push!(thiscycle, current)
    end
    # Just after the loop 'thiscycle' has a new complete cycle (built in the
    # inner loop), we add it to the complete list of cycles (to be returned at
    # the end of the method).
    push!(cycles, thiscycle)
    # If all nodes were visited, then the outer loop ends. This could be
    # the outer while condition.
    sum(visited) == n && break
  end
  # Finally the list of all cycles in the graph is returned.
  return cycles
end

# Example: [1, 5, 2, 4] -> [(1, 5), (2, 5), (2, 4), (1, 4)]
# Note: the first value of a tuple is always the smallest of the two.
function tour2edges(tour :: Vector{T}) :: Vector{Tuple{T, T}} where {T}
  lv = last(tour)
  edges = Vector{Tuple{Int, Int}}()
  for v in tour
    edge = (lv < v ? (lv, v) : (v, lv))
    push!(edges, edge)
    lv = v
  end
  
  return edges
end

function stsp_subtour_elimination_example()
    # For now, direct_model is needed. For the next versions you can use
    # Model instead. It does not make much difference as you will see.
    model = direct_model(CPLEX.Optimizer())

    # We need to use single-thread to avoid crashes. Julia is not yet
    # completely thread-safe in the interaction with C code.
    MOI.set(backend(model), MOI.NumberOfThreads(), 1)

    n, d = get_instance()

    # Create the model (except the lazy constraints).
    @variable(model, x[i = 1:n, j = (i+1):n], Bin) # if edge is used or not
    for i = 1:n
      # For every vertex there should be exactly two edges touching it.
      @constraint(model,
        sum(x[j, i] for j = 1:(i - 1)) + sum(x[i, j] for j = (i+1):n) == 2
      )
    end
    # Minimize the distance.
    @objective(model, Min, sum(d[i, j] * x[i, j] for i = 1:n, j = (i+1):n))

    # The callback that will implement lazy cuts. It is a function defined
    # inside a function so it will capture the variables like 'model' and 'x'
    # and be able to use them.
    function my_callback(cb_context::CPLEX.CallbackContext, context_id::Clong)
      # Do not run if the callback was called in a moment different from
      # "just found an integer solution lets check if its valid".
      (context_id != CPLEX.CPX_CALLBACKCONTEXT_CANDIDATE ||
        CPLEX.cbcandidateispoint(cb_context) == 0) && return
      # If this is not called, the the MOI.get below do not get the var values.
      CPLEX.callbackgetcandidatepoint(backend(model), cb_context, context_id)
      # Set of edges in solution.
      edges = Vector{Tuple{Int, Int}}()

      # Extract the set of edges in the solution.
      for i = 1:n, j = (i+1):n
        # If using JuMP master branch:
        # xval = callback_value(cb_data, x[i, j])

        # If using latest JuMP release:
        #@show backend(model).variable_info
        xval = MOI.get(
          backend(model),
          MOI.CallbackVariablePrimal(cb_context),
          index(x[i, j])
        )
        # If the edge is used add it to the edge set.
        xval > 0.5 && push!(edges, (i, j))
      end

      # From the solution edges to one or more tours.
      ts = tours(n, edges)
      # If it is just one tour, then the solution is valid
      # (not many disconnected subtours).
      length(ts) == 1 && return

      # Now some dark magic is needed.
      for tour in ts
        tedges = tour2edges(tour)
        tvars = [x[i, j] for (i, j) in tedges]
        # To get the CPLEX variable column index from the JuMP variables,
        # the line below is needed.
        ind = Cint[backend(model).variable_info[index(var)].column - 1 for var in tvars]
        # This is the CPLEX function for rejecting a solution and giving a
        # lazy constraint that invalidates it. It is documented in:
        # https://www.ibm.com/support/knowledgecenter/SSSA5P_12.9.0/ilog.odms.cplex.help/refcallablelibrary/cpxapi/callbackrejectcandidate.html
        CPLEX.cbrejectcandidate(
          cb_context, # one of the callback parameters
          Cint(1), # rcnt: number of constraints to add
          Cint(length(tvars)), # nzcnt: number of nonzeros of added constraints
          Cdouble[length(tvars) - 1], # rhs: righthand sides for the constraints
          Cchar['L'], # sense: sense for the constraints
          Cint[0], # rmatbeg: indexes of the two arrays below where
          # each constraint start (the two arrays below may have a
          # variable number of elements for each constraint)
          ind, # rmatind: indexes of nonzero columns
          Cdouble[repeat([1.0], length(ind))...] # rmatval: coefficients of nonzero columns
        )
      end
      return
    end

    # This function annex the callback to the model.
    MOI.set(model, CPLEX.CallbackFunction(), my_callback)

    optimize!(model)

    @test termination_status(model) == MOI.OPTIMAL
    @test primal_status(model) == MOI.FEASIBLE_POINT

    # Show the objective value (from the solver), the objective value
    # (recomputed from the variables), the edges of the solution and
    # the tour it makes.
    @show objective_value(model)
    edges = Vector{Tuple{Int, Int}}()
    iobj = 0.0
    for i = 1:n, j = (i+1):n
      if value(x[i, j]) > 0.5
        push!(edges, (i, j))
        iobj += d[i, j]
      end
    end
    @show iobj
    @show edges
    stours = tours(n, edges)
    @test length(stours) == 1
    solution = first(stours)
    @show solution
end

stsp_subtour_elimination_example()