[Numpy-discussion] Interleaved Arrays and

[Neil Martinsen-Burrell]

On 2009-06-16 16:05 , Robert wrote:
> Neil Martinsen-Burrell wrote:
>> On 06/16/2009 02:18 PM, Robert wrote:
>>>    >>>   n = 10
>>>    >>>   xx = np.ones(n)
>>>    >>>   yy = np.arange(n)
>>>    >>>   aa = np.column_stack((xx,yy))
>>>    >>>   bb = np.column_stack((xx+1,yy))
>>>    >>>   aa
>>> array([[ 1.,  0.],
>>>           [ 1.,  1.],
>>>           [ 1.,  2.],
>>>           [ 1.,  3.],
>>>           [ 1.,  4.],
>>>           [ 1.,  5.],
>>>           [ 1.,  6.],
>>>           [ 1.,  7.],
>>>           [ 1.,  8.],
>>>           [ 1.,  9.]])
>>>    >>>   bb
>>> array([[ 2.,  0.],
>>>           [ 2.,  1.],
>>>           [ 2.,  2.],
>>>           [ 2.,  3.],
>>>           [ 2.,  4.],
>>>           [ 2.,  5.],
>>>           [ 2.,  6.],
>>>           [ 2.,  7.],
>>>           [ 2.,  8.],
>>>           [ 2.,  9.]])
>>>    >>>   np.column_stack((aa,bb))
>>> array([[ 1.,  0.,  2.,  0.],
>>>           [ 1.,  1.,  2.,  1.],
>>>           [ 1.,  2.,  2.,  2.],
>>>           [ 1.,  3.,  2.,  3.],
>>>           [ 1.,  4.,  2.,  4.],
>>>           [ 1.,  5.,  2.,  5.],
>>>           [ 1.,  6.,  2.,  6.],
>>>           [ 1.,  7.,  2.,  7.],
>>>           [ 1.,  8.,  2.,  8.],
>>>           [ 1.,  9.,  2.,  9.]])
>>>    >>>   cc = _
>>>    >>>   cc.reshape((n*2,2))
>>> array([[ 1.,  0.],
>>>           [ 2.,  0.],
>>>           [ 1.,  1.],
>>>           [ 2.,  1.],
>>>           [ 1.,  2.],
>>>           [ 2.,  2.],
>>>           [ 1.,  3.],
>>>           [ 2.,  3.],
>>>           [ 1.,  4.],
>>>           [ 2.,  4.],
>>>           [ 1.,  5.],
>>>           [ 2.,  5.],
>>>           [ 1.,  6.],
>>>           [ 2.,  6.],
>>>           [ 1.,  7.],
>>>           [ 2.,  7.],
>>>           [ 1.,  8.],
>>>           [ 2.,  8.],
>>>           [ 1.,  9.],
>>>           [ 2.,  9.]])
>>>    >>>
>>>
>>>
>>> However I feel too, there is a intuitive abbrev function like
>>> 'interleave' or so missing in numpy shape_base or so.
>>
>> Using fancy indexing, you can set strided portions of an array equal to
>> another array.  So::
>>
>> In [2]: aa = np.empty((10,2))
>>
>> In [3]: aa[:, 0] = 1
>>
>> In [4]: aa[:,1] = np.arange(10)
>>
>> In [5]: bb = np.empty((10,2))
>>
>> In [6]: bb[:,0] = 2
>>
>> In [7]: bb[:,1] = aa[:,1] # this works
>>
>> In [8]: cc = np.empty((20,2))
>>
>> In [9]: cc[::2,:] = aa
>>
>> In [10]: cc[1::2,:] = bb
>>
>> In [11]: cc
>> Out[11]:
>> array([[ 1.,  0.],
>>          [ 2.,  0.],
>>          [ 1.,  1.],
>>          [ 2.,  1.],
>>          [ 1.,  2.],
>>          [ 2.,  2.],
>>          [ 1.,  3.],
>>          [ 2.,  3.],
>>          [ 1.,  4.],
>>          [ 2.,  4.],
>>          [ 1.,  5.],
>>          [ 2.,  5.],
>>          [ 1.,  6.],
>>          [ 2.,  6.],
>>          [ 1.,  7.],
>>          [ 2.,  7.],
>>          [ 1.,  8.],
>>          [ 2.,  8.],
>>          [ 1.,  9.],
>>          [ 2.,  9.]])
>>
>> Using this syntax, interleave could be a one-liner.
>>
>> -Neil
>
> that method of 'filling an empty with a pattern' was mentioned in
> the other (general) interleaving question. It requires however a
> lot of particular numbers and :'s in the code, and requires even
> more statements which can hardly be written in functional style -
> in one line?. The other approach is more jount, free of fancy
> indexing assignments.

jount?  I think that assigning to a strided index is very clear, but 
that is a difference of opinion.  All of the calls to np.empty are the 
equivalent of the column_stack's in your example.  I think that 
operations on segments of arrays are fundamental to an array-processing 
language such as NumPy.  Using ";" you can put as many of those 
statements as you would like one line. :)

> The general interleaving should work efficiently in one like this:
>
> np.column_stack/concatenate((r,g,b,....), axis=...).reshape(..)
>
> But as all this is not intuitive, something like this should be in
> numpy perhaps? :
>
> def interleave( tup_arrays, axis = None )

Here is a minimally tested implementation.  If anyone really wants this 
for numpy, I'll gladly add comments and tests.  I couldn't figure out 
how to automatically find the greatest dtype, so I added an argument to 
specify, otherwise it uses the type of the first array.

def interleave(arrays, axis=0, dtype=None):
     assert len(arrays) > 0
     first = arrays[0]
     assert all([arr.shape == first.shape for arr in arrays])
     new_shape = list(first.shape)
     new_shape[axis] *= len(arrays)
     if dtype is None:
         new_dtype = first.dtype
     else:
         new_dtype = dtype
     interleaved = np.empty(new_shape, new_dtype)
     axis_slice = [slice(None, None, None)]*axis + \
                  [slice(0,None,len(arrays))] + [Ellipsis]
     for i, arr in enumerate(arrays):
         axis_slice[axis] = slice(i, None, len(arrays))
         interleaved[tuple(axis_slice)] = arr
     return interleaved