title: Python, Pairwise 'distance', need a fast way to do it
body: <p>For a side project in my PhD, I engaged in the task of modelling some system in Python. Efficiency wise, my program hits a bottleneck in the following problem, which I'll expose in a Minimal Working Example.</p>
  
  <p>I deal with a large number of segments encoded by their 3D beginning and endpoints, so each segment is represented by 6 scalars.</p>
  
  <p>I need to calculate a pairwise minimal intersegment distance. The analytical expression of the minimal distance between two segments is found in this <a href="http://www.geometrictools.com/Documentation/DistanceLine3Line3.pdf">source</a>. To the MWE:</p>
  
  <pre><code>import numpy as np
  N_segments = 1000
  List_of_segments = np.random.rand(N_segments, 6)
  
  Pairwise_minimal_distance_matrix = np.zeros( (N_segments,N_segments) )
  for i in range(N_segments):
      for j in range(i+1,N_segments): 
  
          p0 = List_of_segments[i,0:3] #beginning point of segment i
          p1 = List_of_segments[i,3:6] #end point of segment i
          q0 = List_of_segments[j,0:3] #beginning point of segment j
          q1 = List_of_segments[j,3:6] #end point of segment j
          #for readability, some definitions
          a = np.dot( p1-p0, p1-p0)
          b = np.dot( p1-p0, q1-q0)
          c = np.dot( q1-q0, q1-q0)
          d = np.dot( p1-p0, p0-q0)
          e = np.dot( q1-q0, p0-q0)
          s = (b*e-c*d)/(a*c-b*b)
          t = (a*e-b*d)/(a*c-b*b)
          #the minimal distance between segment i and j
          Pairwise_minimal_distance_matrix[i,j] = sqrt(sum( (p0+(p1-p0)*s-(q0+(q1-q0)*t))**2)) #minimal distance
  </code></pre>
  
  <p>Now, I realize this is extremely inefficient, and this is why I am here. I have looked extensively in how to avoid the loop, but I run into a bit of a problem. Apparently, this sort of calculations is best done with the <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.cdist.html#scipy.spatial.distance.cdist">cdist</a> of python. However, the custom distance functions it can handle have to be binary functions. This is a problem in my case, because my vectors have specifically a length of 6, and have to bit split into their first and last 3 components. I don't think I can translate the distance calculation into a binary function.</p>
  
  <p>Any input is appreciated.</p>