Constructing A 3d Cube Of Points From A List

I have a list pts containing N points (Python floats). I wish to construct a NumPy array of dimension N*N*N*3 such that the array is equivalent to: for i in xrange(0, N): for

Solution 1:

I think that the following should work:

pts = np.array(pts)  #Skip if pts is a numpy array already
lp = len(pts)
arr = np.zeros((lp,lp,lp,3))
arr[:,:,:,0] = pts[:,None,None]  #None is the same as np.newaxis
arr[:,:,:,1] = pts[None,:,None]
arr[:,:,:,2] = pts[None,None,:]

A quick test:

import numpy as np
import timeit

defmeth1(pts):
   pts = np.array(pts)  #Skip if pts is a numpy array already
   lp = len(pts)
   arr = np.zeros((lp,lp,lp,3))
   arr[:,:,:,0] = pts[:,None,None]  #None is the same as np.newaxis
   arr[:,:,:,1] = pts[None,:,None]
   arr[:,:,:,2] = pts[None,None,:]
   return arr

defmeth2(pts):
   lp = len(pts)
   N = lp
   arr = np.zeros((lp,lp,lp,3))
   for i in xrange(0, N):
      for j in xrange(0, N):
         for k in xrange(0, N):
            arr[i,j,k,0] = pts[i]
            arr[i,j,k,1] = pts[j]
            arr[i,j,k,2] = pts[k]

   return arr

pts = range(10)
a1 = meth1(pts)
a2 = meth2(pts)

print np.all(a1 == a2)

NREPEAT = 10000print timeit.timeit('meth1(pts)','from __main__ import meth1,pts',number=NREPEAT)
print timeit.timeit('meth2(pts)','from __main__ import meth2,pts',number=NREPEAT)

results in:

True0.873255968094#my way11.4249279499#original

So this new method is an order of magnitude faster as well.

Solution 2:

import numpy as np
N = 10
pts = xrange(0,N)
l = [ [ [ [ pts[i],pts[j],pts[k] ]  for k in xrange(0,N) ] for j in xrange(0,N) ] for i in xrange(0,N) ]
x = np.array(l, np.int32)
print x.shape # (10,10,10,3)

Solution 3:

This can be done in two lines:

defmeth3(pts):
    arrs = np.broadcast_arrays(*np.ix_(pts, pts, pts))
    return np.concatenate([a[...,None] for a in arrs], axis=3)

However, this method is not as fast as mgilson's answer, because concatenate is annoyingly slow. A generalized version of his answer performs roughly as well, though, and can generate the result you want (i.e. an n-dimensional cartesian product contained within an n-dimensional grid) for any set of arrays.

defmeth4(arrs):     # or meth4(*arrs) for a simplified interface
    arr = np.empty([len(a) for a in arrs] + [len(arrs)])
    for i, a inenumerate(np.ix_(*arrs)):
        arr[...,i] = a
    return arr

This accepts any sequence of sequences, as long as it can be converted into a sequence of numpy arrays:

>>> meth4([[0, 1], [2, 3]])
array([[[ 0.,  2.],
        [ 0.,  3.]],

       [[ 1.,  2.],
        [ 1.,  3.]]])

And the cost of this generality isn't too high -- it's only twice as slow for small pts arrays:

>>>(meth4([pts, pts, pts]) == meth1(pts)).all()
True
>>>%timeit meth4([pts, pts, pts])
10000 loops, best of 3: 27.4 us per loop
>>>%timeit meth1(pts)
100000 loops, best of 3: 13.1 us per loop

And it's actually a bit faster for larger ones (although the speed gain is probably due to my use of empty instead of zeros):

>>>pts = np.linspace(0, 1, 100)>>>%timeit meth4([pts, pts, pts])
100 loops, best of 3: 13.4 ms per loop
>>>%timeit meth1(pts)
100 loops, best of 3: 16.7 ms per loop