Travelling salesman with a directional constraint
<p>I am trying to order an array of 3D coordinates by their order along a path. A sample:</p>

<pre><code>points = np.array([[ 0.81127451,  0.22794118,  0.52009804],
                   [ 0.62986425,  0.4546003 ,  0.12971342],
                   [ 0.50666667,  0.41137255,  0.65215686],
                   [ 0.79526144,  0.58186275,  0.04738562],
                   [ 0.55163399,  0.49803922,  0.24117647],
                   [ 0.47385621,  0.64084967,  0.10653595]])
</code></pre>

<p>The points are in random order, but there is always a single path through them. I am finding the path with an adapted travelling salesman problem (TSP) soution, using the <a href="http://www.akira.ruc.dk/~keld/research/LKH/">LKH solver</a> (Helsgaun 2009). It involves two modifications:</p>

<ul>
<li>Add a point at or near the origin. This finds the best starting point in every instance I've tackled so far. This was my idea, I have no other basis for it.</li>
<li>Add a point at a distance of zero from every point. This makes the solver find a route to the other end of the path. This idea was from <a href="http://stackoverflow.com/questions/6733999/what-is-the-problem-name-for-traveling-salesman-problemtsp-without-considering">this SO question</a>.</li>
</ul>

<p>Note that the TSP does not involve <em>position</em>, only the <em>distances</em> between nodes. So the solver does 'know' (or care) that I'm working in 3D. I just make a distance matrix like so:</p>

<pre><code>import numpy as np
from scipy.spatial.distance import pdist, squareform

# Add a point near the origin.
points = np.vstack([[[0.25, 0, 0.5]], points])
dists = squareform(pdist(points, 'euclidean'))

# Normalize to int16s because the solver likes it.
d = 32767 * dists / np.sqrt(3)
d = d.astype(np.int16)

# Add a point that is zero units from every other point.
row, col = d.shape
d = np.insert(d, row, 0, axis=0)
d = np.insert(d, col, 0, axis=1)
</code></pre>

<p>I pass this to <a href="https://github.com/kwinkunks/pytsp">my fork of <code>pytsp</code></a>, which passes it to the LKH solver. And everything is fine... except when the path crosses itself. TSP solutions cannot have closed loops, so I always get the open loop shown on the right here:</p>

<p><a href="http://i.stack.imgur.com/c4lDk.png"><img src="http://i.stack.imgur.com/c4lDk.png" alt="Travelling salesman paths"></a></p>

<p><em>Note that this is an analogous 2D version of my situation. Note also that the points are imperfectly aligned, even along the 'straight' bits.</em></p>

<p><strong>So my question is: how can I help the solver preserve the direction of the path whenever possible?</strong> I've got two ill-formed ideas, but so far been unable to implement anything:</p>

<ul>
<li>Use another metric instead of L2. But I don't think this can work, because at a given junction, there's nothing inherently different about the 'wrong' point. Its wrongness depends on the previous point. And we don't know yet which is the previous point (that's what we're trying to figure out). So I think this is no good.</li>
<li>Evaluate the local colinearity of every set of three points (e.g. using the determinant of every triple). Modulate the local '3D slope' (not sure what I mean) by this colinearity coefficient. Give every point another dimension expressing this local alignment. Now the norm will reflect local alignment and (hopefully) roughly colinear things will join up.</li>
</ul>

<p>I have put these files on Dropbox:</p>

<ul>
<li><a href="https://www.dropbox.com/s/qx5fcytojra07yg/raw_data.npy?dl=0">Raw data NumPy file</a> </li>
<li><a href="https://www.dropbox.com/s/63nun41m4vbkhzl/ordered_data.npy?dl=0">Ordered data NumPy file</a></li>
</ul>

<p>Thank you for reading; any ideas appreciated.</p>

<h3>Reference</h3>

<p>K. Helsgaun, General k-opt submoves for the Lin-Kernighan TSP heuristic. Mathematical Programming Computation, 2009, <a href="http://dx.doi.org/10.1007/s12532-009-0004-6">doi: 10.1007/s12532-009-0004-6</a>.</p>

<p>Judging by the documentation on pytsp, the distance matrix doesn't have to be symmetric. This means that you could modify the L2 norm to incorporate information on a preferred direction into that matrix. Say you have a preferred direction for some pairs of points (i,j), then for each of these point you could divide <code>dists[i,j]</code> by <code>(1+a)</code> and multiply <code>dists[j,i]</code> by <code>(1+a)</code> to make that direction more favourable. This means that if your algorithm is sure to find the global optimum, you can force it to satisfy your preferred direction by taking <code>a</code> is sufficiently large.</p>

<p>Also, I'm not sure it's impossible to have closed loops in a solution where the distance matrix is taken from 3D data. It seems to me that the 'no closed loops' is a result (of the triangle inequality) specific to 2D.</p>

python numpy linear-algebra graph-algorithm traveling-salesman