This solution uses Google's OR-Tools. Install ortools via python -m pip install --upgrade --user ortools.
from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp
Generate points:
import numpy as np
points = np.random.RandomState(42).rand(100,2)
Calcualte distance matrix as in question:
z = np.array([[complex(*c) for c in points]]) # notice the [[ ... ]]
distance_matrix_pre = abs(z.T-z)
Since the docs say:
Note: Since the routing solver does all computations with integers, the distance callback must return an integer distance for
any two locations. If any of the entries of data['distance_matrix']
are not integers, you need to round either the matrix entries, or the
return values of the callback, to integers. See Scaling the distance
matrix for an example that shows how to avoid problems caused by
rounding error.
I am going to multiply by a 10000 every entry in our distance matrix and round down to an int. If want to use more decimal places, multiply by a larger number.
distance_matrix = np.floor(distance_matrix_pre*10000).astype(int).tolist()
Keep following Google's tutorial:
# Create the routing index manager.
manager = pywrapcp.RoutingIndexManager(len(points),1,0)
# Create Routing Model.
routing = pywrapcp.RoutingModel(manager)
def distance_callback(from_index, to_index):
"""Returns the distance between the two nodes."""
# Convert from routing variable Index to distance matrix NodeIndex.
from_node = manager.IndexToNode(from_index)
to_node = manager.IndexToNode(to_index)
return distance_matrix[from_node][to_node]
transit_callback_index = routing.RegisterTransitCallback(distance_callback)
# Define cost of each arc.
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)
The output depends on initial path-selection strategy used (similar to how a numerical solution can depend on inital guess). Here we specify initial strategy:
# Setting first solution heuristic.
search_parameters = pywrapcp.DefaultRoutingSearchParameters()
search_parameters.first_solution_strategy = (
routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)
In Google's words:
The code sets the first solution strategy to PATH_CHEAPEST_ARC,
which creates an initial route for the solver by repeatedly adding
edges with the least weight that don't lead to a previously visited
node (other than the depot). For other options, see First solution
strategy.
Solve the problem with initial strategy defined above (ie with PATH_CHEAPEST_ARC):
# Solve the problem.
solution = routing.SolveWithParameters(search_parameters)
Get a route:
if solution:
index = routing.Start(0)
route = [manager.IndexToNode(index)]
while not routing.IsEnd(index):
index = solution.Value(routing.NextVar(index))
route.append(manager.IndexToNode(index))
Plot result:
plt.scatter(points[:,0],points[:,1])
for a, b in zip(route[:-1],route[1:]):
x = points[[a,b]].T[0]
y = points[[a,b]].T[1]
plt.plot(x, y,c='r',zorder=-1)
plt.gca().set_aspect('equal')
The list of other initial strategies which can be used (source):
initial_strategy_list = \
['AUTOMATIC','PATH_CHEAPEST_ARC','PATH_MOST_CONSTRAINED_ARC','EVALUATOR_STRATEGY','SAVINGS','SWEEP','CHRISTOFIDES','ALL_UNPERFORMED','BEST_INSERTION','PARALLEL_CHEAPEST_INSERTION','LOCAL_CHEAPEST_INSERTION','GLOBAL_CHEAPEST_ARC','LOCAL_CHEAPEST_ARC','FIRST_UNBOUND_MIN_VALUE']
If want to increase confidence in the solution, could iterate over all of them and find the best one.
