In [1]: import numpy as np
In [2]: a = np.array([[2,0],[3,0],[3,1],[5,0],[5,1],[5,2]])
In [3]: b = np.zeros((6,3), dtype='int32')
In [4]: b[a[:,0], a[:,1]] = 10
In [5]: b
Out[5]:
array([[ 0, 0, 0],
[ 0, 0, 0],
[10, 0, 0],
[10, 10, 0],
[ 0, 0, 0],
[10, 10, 10]])
Why it works:
If you index b with two numpy arrays in an assignment,
b[x, y] = z
then think of NumPy as moving simultaneously over each element of x and each element of y and each element of z (let's call them xval, yval and zval), and assigning to b[xval, yval] the value zval. When z is a constant, "moving over z just returns the same value each time.
That's what we want, with x being the first column of a and y being the second column of a. Thus, choose x = a[:, 0], and y = a[:, 1].
b[a[:,0], a[:,1]] = 10
Why b[a] = 10 does not work
When you write b[a], think of NumPy as creating a new array by moving over each element of a, (let's call each one idx) and placing in the new array the value of b[idx] at the location of idx in a.
idx is a value in a. So it is an int32. b is of shape (6,3), so b[idx] is a row of b of shape (3,). For example, when idx is
In [37]: a[1,1]
Out[37]: 0
b[a[1,1]] is
In [38]: b[a[1,1]]
Out[38]: array([0, 0, 0])
So
In [33]: b[a].shape
Out[33]: (6, 2, 3)
So let's repeat: NumPy is creating a new array by moving over each element of a and placing in the new array the value of b[idx] at the location of idx in a. As idx moves over a, an array of shape (6,2) would be created. But since b[idx] is itself of shape (3,), at each location in the (6,2)-shaped array, a (3,)-shaped value is being placed. The result is an array of shape (6,2,3).
Now, when you make an assignment like
b[a] = 10
a temporary array of shape (6,2,3) with values b[a] is created, then the assignment is performed. Since 10 is a constant, this assignment places the value 10 at each location in the (6,2,3)-shaped array.
Then the values from the temporary array are reassigned back to b.
See reference to docs. Thus the values in the (6,2,3)-shaped array are copied back to the (6,3)-shaped b array. Values overwrite each other. But the main point is you do not obtain the assignments you desire.
In [1]: import numpy as np
In [2]: a = np.array([[2,0],[3,0],[3,1],[5,0],[5,1],[5,2]])
In [3]: b = np.zeros((6,3), dtype='int32')
In [4]: b[a[:,0], a[:,1]] = 10
In [5]: b
Out[5]:
array([[ 0, 0, 0],
[ 0, 0, 0],
[10, 0, 0],
[10, 10, 0],
[ 0, 0, 0],
[10, 10, 10]])
Why it works:
If you index b with two numpy arrays in an assignment,
b[x, y] = z
then think of NumPy as moving simultaneously over each element of x and each element of y and each element of z (let's call them xval, yval and zval), and assigning to b[xval, yval] the value zval. When z is a constant, "moving over z just returns the same value each time.
That's what we want, with x being the first column of a and y being the second column of a. Thus, choose x = a[:, 0], and y = a[:, 1].
b[a[:,0], a[:,1]] = 10
Why b[a] = 10 does not work
When you write b[a], think of NumPy as creating a new array by moving over each element of a, (let's call each one idx) and placing in the new array the value of b[idx] at the location of idx in a.
idx is a value in a. So it is an int32. b is of shape (6,3), so b[idx] is a row of b of shape (3,). For example, when idx is
In [37]: a[1,1]
Out[37]: 0
b[a[1,1]] is
In [38]: b[a[1,1]]
Out[38]: array([0, 0, 0])
So
In [33]: b[a].shape
Out[33]: (6, 2, 3)
So let's repeat: NumPy is creating a new array by moving over each element of a and placing in the new array the value of b[idx] at the location of idx in a. As idx moves over a, an array of shape (6,2) would be created. But since b[idx] is itself of shape (3,), at each location in the (6,2)-shaped array, a (3,)-shaped value is being placed. The result is an array of shape (6,2,3).
Now, when you make an assignment like
b[a] = 10
a temporary array of shape (6,2,3) with values b[a] is created, then the assignment is performed. Since 10 is a constant, this assignment places the value 10 at each location in the (6,2,3)-shaped array.
Then the values from the temporary array are reassigned back to b.
See reference to docs. Thus the values in the (6,2,3)-shaped array are copied back to the (6,3)-shaped b array. Values overwrite each other. But the main point is you do not obtain the assignments you desire.
TL;DR: Use advanced indexing: b[*a.T] = 10
You can also transpose the index array a, convert the result into a tuple and index the array b and assign a value. Converting the index array into a tuple (or unpacking it inside a []) ensures that multidimensional indexing works as expected. This is assignment by advanced indexing.
a = np.array([[2, 0], [3, 0], [3, 1], [5, 0], [5, 1], [5, 2]])
b = np.zeros((6,3), dtype ='int32')
b[*a.T] = 10
# or
b[tuple(a.T)] = 10
# or
b[(*a.T,)] = 10
# or
b[(*a.T.tolist(),)] = 10
All of them produce the expected output of
array([[ 0, 0, 0],
[ 0, 0, 0],
[10, 0, 0],
[10, 10, 0],
[ 0, 0, 0],
[10, 10, 10]])
I don't know of an automatic way to do it, but if
a = [[1,2],[3,4],[5,6]]
and you want to find the location of 3, you can do:
x = [x for x in a if 3 in x][0]
print 'The index is (%d,%d)'%(a.index(x),x.index(3))
The output is:
The index is (1,0)
For two dimensional list; you can iterate over rows and using .index function for looking for item:
def find(l, elem):
for row, i in enumerate(l):
try:
column = i.index(elem)
except ValueError:
continue
return row, column
return -1
tl = [[1,2,3],[4,5,6],[7,8,9]]
print(find(tl, 6)) # (1,2)
print(find(tl, 1)) # (0,0)
print(find(tl, 9)) # (2,2)
print(find(tl, 12)) # -1
Another option is multidimensional list-of-locations indexing:
import numpy as np
ncol = 10 # 10 in your case
nrow = 500 # 500 in your case
# just creating some test data:
x = np.arange(ncol*nrow).reshape(nrow,ncol)
y = (ncol * np.random.random_sample((nrow, 1))).astype(int)
print(x)
print(y)
print(x[np.arange(nrow),y.T].T)
The syntax is explained here. You basically need an array of indices for each dimension. In the first dimension this is simply [0,...,500] in your case and the second dimension is your y-array. We need to transpose it (.T), because it has to have the same shape as the first and the output array. The second transposition is not really needed, but gives you the shape you want.
EDIT:
The question of performance came up and I tried the three methods mentioned so far. You'll need line_profiler to run the following with
kernprof -l -v tmp.py
where tmp.py is:
import numpy as np
@profile
def calc(x,y):
z = np.arange(nrow)
a = x[z,y.T].T # mine, with the suggested speed up
b = x[:,y].diagonal().T # Christoph Terasa
c = np.array([i[j] for i, j in zip(x, y)]) # tobias_k
return (a,b,c)
ncol = 5 # 10 in your case
nrow = 10 # 500 in your case
x = np.arange(ncol*nrow).reshape(nrow,ncol)
y = (ncol * np.random.random_sample((nrow, 1))).astype(int)
a, b, c = calc(x,y)
print(a==b)
print(b==c)
The output for my python 2.7.6:
Line # Hits Time Per Hit % Time Line Contents
==============================================================
3 @profile
4 def calc(x,y):
5 1 4 4.0 0.1 z = np.arange(nrow)
6 1 35 35.0 0.8 a = x[z,y.T].T
7 1 3409 3409.0 76.7 b = x[:,y].diagonal().T
8 501 995 2.0 22.4 c = np.array([i[j] for i, j in zip(x, y)])
9
10 1 1 1.0 0.0 return (a,b,c)
Where %Time or Time are the relevant columns. I don't know how to profile memory consumption, someone else would have to do that. For now it looks like my solution is the fastest for the requested dimensions.
While not really intuitive from a syntactic perspective
X[:,Y].diagonal()[0]
will give you the values you're looking for. The fancy indexing selects from each row of X all values in Y, and diagonal selects only those at the indexes where i == j. The indexing with [0] at the end just flattens the 2d array.
We can compute the linear indices and then use np.take -
np.take(lut, np.ravel_multi_index(arr.T, lut.shape)).T
If you are open to alternatives, we can reshape the indices array to 2D, convert to tuples, index into the data array with it, to give us 1D, which could be reshaped back to 2D -
lut[tuple(arr.reshape(-1,arr.shape[-1]).T)].reshape(arr.shape[:2])
Sample run -
In [49]: lut = np.random.randint(11,99,(13,13,13))
In [50]: arr = np.arange(12).reshape([2,2,3])
In [51]: lut[ arr[:,:,0],arr[:,:,1],arr[:,:,2] ] # Original approach
Out[51]:
array([[41, 21],
[94, 22]])
In [52]: np.take(lut, np.ravel_multi_index(arr.T, lut.shape)).T
Out[52]:
array([[41, 21],
[94, 22]])
In [53]: lut[tuple(arr.reshape(-1,arr.shape[-1]).T)].reshape(arr.shape[:2])
Out[53]:
array([[41, 21],
[94, 22]])
We can avoid the double transposing for the np.take approach, like so -
In [55]: np.take(lut, np.ravel_multi_index(arr.transpose(2,0,1), lut.shape))
Out[55]:
array([[41, 21],
[94, 22]])
Generalize to multi-dimensional arrays of generic dimensions
This could be generalized to ndarrays of generic no. of dims, like so -
np.take(lut, np.ravel_multi_index(np.rollaxis(arr,-1,0), lut.shape))
The tuple-based approach should work without any change.
Here's a sample run for the same -
In [95]: lut = np.random.randint(11,99,(13,13,13,13))
In [96]: arr = np.random.randint(0,13,(2,3,4,4))
In [97]: lut[ arr[:,:,:,0] , arr[:,:,:,1],arr[:,:,:,2],arr[:,:,:,3] ]
Out[97]:
array([[[95, 11, 40, 75],
[38, 82, 11, 38],
[30, 53, 69, 21]],
[[61, 74, 33, 94],
[90, 35, 89, 72],
[52, 64, 85, 22]]])
In [98]: np.take(lut, np.ravel_multi_index(np.rollaxis(arr,-1,0), lut.shape))
Out[98]:
array([[[95, 11, 40, 75],
[38, 82, 11, 38],
[30, 53, 69, 21]],
[[61, 74, 33, 94],
[90, 35, 89, 72],
[52, 64, 85, 22]]])
I did not try in 3-dimensions. But in 2-dimensions I get the result I want with using numpy.take :
np.take(np.take(T,ix,axis=0), iy,axis=1 )
Perhaps you can expand that to 3-dimensions.
As an example I can adress with two 1-dim arrays for the indices ix and iy the 2-dimensional stencil for the discrete Laplace equation,
ฮT = T[ix-1,iy] + T[ix+1, iy] + T[ix,iy-1] + T[ix,iy+1] - 4*T[ix,iy]
Introducing for a leaner writing:
def q(ฮฆ,kx,ky):
return np.take(np.take(ฮฆ,kx,axis=0), ky,axis=1 )
then I can run the following python code with numpy.take:
nx = 6; ny= 10
T = np.arange(nx*ny).reshape(nx, ny)
ix = np.linspace(1,nx-2,nx-2,dtype=int)
iy = np.linspace(1,ny-2,ny-2,dtype=int)
ฮT = q(T,ix-1,iy) + q(T,ix+1,iy) + q(T,ix,iy-1) + q(T,ix,iy+1) - 4.0 * q(T,ix,iy)
The broadcasting is more like:
In [280]: y[np.array([0,2,4])[...,None], np.array([1,2])]
Out[280]:
array([[ 1, 2],
[15, 16],
[29, 30]])
I added a dimension to [0,2,4] making it 2d. broadcast_arrays can be used to see what the broadcasted arrays look like:
In [281]: np.broadcast_arrays(np.array([0,2,4])[...,None], np.array([1,2]))
Out[281]:
[array([[0, 0],
[2, 2],
[4, 4]]),
array([[1, 2],
[1, 2],
[1, 2]])]
np.broadcast_arrays([[0],[2],[4]], [1,2]) is the samething without the array wrappers. np.meshgrid([0,2,4], [1,2], indexing='ij') is another way of producing these indexing arrays.
(the lists produced by meshgrid or broadcast_arrays could be used as the argument for y[_].)
So it's right to say [1,2] is broadcast with the index array, but it omits the bit about adjusting dimensions.
A little earlier they have this example:
y[np.array([0,2,4])]
which is equivalent to y[np.array([0,2,4]), :]. It picks 3 rows, and all items from them. The 1:3 case can be thought of as an extension of this, picking 3 rows, and then 2 columns.
y[[0,2,4],:][:,1:3]
This might be a better way of thinking about the indexing if broadcasting is too confusing.
There's another docs page that might handle this better
http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html
In this docs, basic indexing involves slices and integers
y[:,1:3], y[1,:], y[1, 1:3]
Advanced indexing involves an array (or list)
y[[0,2,4],:]
This produces the same result as y[::2,:], except the list case produces a copy, the slice (basic) a view.
y[[0,2,4], [1,2,3]] is a case of pure advance index array indexing, the result is 3 items, ones at (0,1), (2,2), and (4,3).
y[[0,2,4], 1:3] is a case that this docs calls Combining advanced and basic indexing, 'advanced' from `[0,2,4]', basic from '1:3'.
http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html#combining-advanced-and-basic-indexing
Looking at a more complex index array might add some insight.
In [222]: i=[[0,2],[1,4]]
Used with another list, it is 'pure' advanced, and the result is broadcasted:
In [224]: y[i, [1,2]]
Out[224]:
array([[ 1, 16],
[ 8, 30]])
The index arrays are:
In [234]: np.broadcast_arrays(i, [1,2])
Out[234]:
[array([[0, 2],
[1, 4]]),
array([[1, 2],
[1, 2]])]
The [1,2] list is just expanded to a (2,2) array.
Using it with a slice is an example of this mixed advanced/basic, and the result is 3d (2,2,2).
In [223]: y[i, 1:3]
Out[223]:
array([[[ 1, 2],
[15, 16]],
[[ 8, 9],
[29, 30]]])
The equivalent with broadcasting is
y[np.array(i)[...,None], [1,2]]
y[data,beginIndex:endIndex]
import numpy as np
y = np.arange(35).reshape(5,7)
print(y)
[[ 0 1 2 3 4 5 6]
[ 7 8 9 10 11 12 13]
[14 15 16 17 18 19 20]
[21 22 23 24 25 26 27]
[28 29 30 31 32 33 34]]
print(y[np.array([0,2,4]),1:3])
[[ 1 2]
[15 16]
[29 30]]