Index N Dimensional Array with (N-1) D Array

Index n dimensional array with (n-1) d array

Make use of advanced-indexing -

m,n = a.shape[1:]
I,J = np.ogrid[:m,:n]
a_max_values = a[idx, I, J]
b_max_values = b[idx, I, J]

---

For the general case:

def argmax_to_max(arr, argmax, axis):
    """argmax_to_max(arr, arr.argmax(axis), axis) == arr.max(axis)"""
    new_shape = list(arr.shape)
    del new_shape[axis]

    grid = np.ogrid[tuple(map(slice, new_shape))]
    grid.insert(axis, argmax)

    return arr[tuple(grid)]

Quite a bit more awkward than such a natural operation should be, unfortunately.

For indexing a n dim array with a (n-1) dim array, we could simplify it a bit to give us the grid of indices for all axes, like so -

def all_idx(idx, axis):
    grid = np.ogrid[tuple(map(slice, idx.shape))]
    grid.insert(axis, idx)
    return tuple(grid)

Hence, use it to index into input arrays -

axis = 0
a_max_values = a[all_idx(idx, axis=axis)]
b_max_values = b[all_idx(idx, axis=axis)]

Index n dimensional array with (n-1) d array

Make use of advanced-indexing -

m,n = a.shape[1:]
I,J = np.ogrid[:m,:n]
a_max_values = a[idx, I, J]
b_max_values = b[idx, I, J]

---

For the general case:

def argmax_to_max(arr, argmax, axis):
    """argmax_to_max(arr, arr.argmax(axis), axis) == arr.max(axis)"""
    new_shape = list(arr.shape)
    del new_shape[axis]

    grid = np.ogrid[tuple(map(slice, new_shape))]
    grid.insert(axis, argmax)

    return arr[tuple(grid)]

Quite a bit more awkward than such a natural operation should be, unfortunately.

For indexing a n dim array with a (n-1) dim array, we could simplify it a bit to give us the grid of indices for all axes, like so -

def all_idx(idx, axis):
    grid = np.ogrid[tuple(map(slice, idx.shape))]
    grid.insert(axis, idx)
    return tuple(grid)

Hence, use it to index into input arrays -

axis = 0
a_max_values = a[all_idx(idx, axis=axis)]
b_max_values = b[all_idx(idx, axis=axis)]

Slice numpy ndarry of arbitrary dimension to 1d array given a list of indices

Some improvements to the solution above (there may still be a one-liner or a more performant solution):

def make_custom_slice(n, indices):
    s = indices.copy()
    s[n] = slice(None)
    return tuple(s)

for n in range(arr.ndim):
    print(arr[make_custom_slice(n, indices)])

An integer value idx can be used to replace the slice object slice(idx, idx+1). Because most indices are copied over directly, start with a copy of indices rather than building the list from scratch.

When built in this way, the result of arr[make_custom_slice(n, indices) has the expected dimension and np.squeeze is unnecessary.

A better way to access n-d array element with a 1-d index array in C++?

I figured out a way to solve this myself.

The idea is that use void * pointers, we know that every memory cell holds value or an address of a memory cell, so we can directly compute the offset of the target to the base address.

In this case, we use void *p = arr to get the base address of the n-d array, and then loop over the array idx, to calculate the offset.

For arr[10][10][10][10], the offset between arr[0] and arr[1] is 10 * 10 * 10 * sizeof(int), since arr is 4-d, arr[0] and arr[1] is 3-d, so there is 10 * 10 * 10 = 1000 elements between arr[0] and arr[1], after that, we should know that the offset between two void * adjacent addresses is 1 byte, so we should multiply sizeof(int) to get the correct offset, according to this, we finally get the exact address of the memory cell we want to access.

Finally, we have to cast void * pointer to int * pointer and access the address to get the correct int value, that's it!

With void *(not so good)

#include <bits/stdc++.h>
using namespace std;

#define N 10

int main()
{
    int arr[N][N][N][N] = {0};
    int idx[4] = {1, 2, 3, 4};
    arr[1][2][3][4] = 1;
    cout<<"Expected: "<<arr[1][2][3][4]<<" at "<<&arr[1][2][3][4]<<endl;
    cout<<"Got with ****: ";
    cout<<*(*(*(*(arr + idx[0]) + idx[1]) + idx[2]) + idx[3])<<endl;

    void *p = arr;
    for(int i = 0; i < 4; i++)
        p += idx[i] * int(pow(10, 3-i)) * sizeof(int);
    cout<<"Got with void *:";
    cout<<*((int*)p)<<" at "<<p<<endl;

    return 0;
}

Output

Expected: 1 at 0x7fff5e3a3f18
Got with ****: 1
Got with void *:1 at 0x7fff5e3a3f18

Notice:
There is a warning when compiling it, but I choose to ignore it.

test.cpp: In function 'int main()':
test.cpp:23:53: warning: pointer of type 'void *' used in arithmetic [-Wpointer-arith]
         p += idx[i] * int(pow(10, 3-i)) * sizeof(int);

Use char * instead of void *(better)

Since we want to manipulate pointer byte by byte, it would be better to use char * to replace void *.

#include <bits/stdc++.h>
using namespace std;

#define N 10

int main()
{
    int arr[N][N][N][N] = {0};
    int idx[4] = {1, 2, 3, 4};
    arr[1][2][3][4] = 1;
    cout<<"Expected: "<<arr[1][2][3][4]<<" at "<<&arr[1][2][3][4]<<endl;

    char *p = (char *)arr;
    for(int i = 0; i < 4; i++)
        p += idx[i] * int(pow(10, 3-i)) * sizeof(int);
    cout<<"Got with char *:";
    cout<<*((int*)p)<<" at "<<(void *)p<<endl;

    return 0;
}

Output

Expected: 1 at 0x7fff4ffd7f18
Got with char *:1 at 0x7fff4ffd7f18

With int *(In this specific case)

I have been told it's not a good practice for void * used in arithmetic, it would be better to use int *, so I cast arr into int * pointer and also replace pow.

#include <bits/stdc++.h>
using namespace std;

#define N 10

int main()
{
    int arr[N][N][N][N] = {0};
    int idx[4] = {1, 2, 3, 4};
    arr[1][2][3][4] = 1;
    cout<<"Expected: "<<arr[1][2][3][4]<<" at "<<&arr[1][2][3][4]<<endl;
    cout<<"Got with ****: ";
    cout<<*(*(*(*(arr + idx[0]) + idx[1]) + idx[2]) + idx[3])<<endl;
    int *p = (int *)arr;
    int offset = 1e3;
    for(int i = 0; i < 4; i++)
    {
        p += idx[i] * offset;
        offset /= 10;
    }
    cout<<"Got with int *:";
    cout<<*p<<" at "<<p<<endl;

    return 0;
}

Output

Expected: 1 at 0x7fff5eaf9f08
Got with ****: 1
Got with int *:1 at 0x7fff5eaf9f08

How do I calculate a one-dimensional index from a multidimensional array? For example two dimensions to one dimension : y * xTot + x = i. But bigger

int index = d * aTot * bTot * cTot
          + c * aTot * bTot
          + b * ATot
          + a;

a = index % aTot;
index = (index-a)/aTot;
b = index % bTot;
index = (index-b)/bTot;
...

Indices of unique values in n-dimensional array

Here's a vectorized approach, which works for arrays of an arbitrary amount of dimensions. The idea of this solution is to extend the functionality of the return_index method in np.unique, and return an array of arrays, each containing the N-dimensional indices of unique values in a numpy array.

For a more compact solution, I've defined the following function along with some explanations throughout the different steps:

def ndix_unique(x):
    """
    Returns an N-dimensional array of indices
    of the unique values in x
    ----------
    x: np.array
       Array with arbitrary dimensions
    Returns
    -------
    - 1D-array of sorted unique values
    - Array of arrays. Each array contains the indices where a
      given value in x is found
    """
    x_flat = x.ravel()
    ix_flat = np.argsort(x_flat)
    u, ix_u = np.unique(x_flat[ix_flat], return_index=True)
    ix_ndim = np.unravel_index(ix_flat, x.shape)
    ix_ndim = np.c_[ix_ndim] if x.ndim > 1 else ix_flat
    return u, np.split(ix_ndim, ix_u[1:])

---

Checking with the array from the question -

a = np.array([[1, 0, 1],[2, 2, 0]])

vals, ixs = ndix_unique(a)

print(vals)
array([0, 1, 2])

print(ixs)
[array([[0, 1],
        [1, 2]]), 
 array([[0, 0],
        [0, 2]]), 
 array([[1, 0],
        [1, 1]])]

Lets try with this other case:

a = np.array([[1,1,4],[2,2,1],[3,3,1]])

vals, ixs = ndix_unique(a)

print(vals)
array([1, 2, 3, 4])

print(ixs)
array([array([[0, 0],
              [0, 1],
              [1, 2],
              [2, 2]]),
       array([[1, 0],
              [1, 1]]), 
       array([[2, 0],
              [2, 1]]),
       array([[0, 2]])], dtype=object)

For a 1D array:

a = np.array([1,5,4,3,3])

vals, ixs = ndix_unique(a)

print(vals)
array([1, 3, 4, 5])

print(ixs)
array([array([0]), array([3, 4]), array([2]), array([1])], dtype=object)

Finally another example with a 3D ndarray:

a = np.array([[[1,1,2]],[[2,3,4]]])

vals, ixs = ndix_unique(a)

print(vals)
array([1, 2, 3, 4])

print(ixs)
array([array([[0, 0, 0],
              [0, 0, 1]]),
       array([[0, 0, 2],
              [1, 0, 0]]), 
       array([[1, 0, 1]]),
       array([[1, 0, 2]])], dtype=object)

Index n dimensional array with (n-1) d array

Make use of advanced-indexing -

m,n = a.shape[1:]
I,J = np.ogrid[:m,:n]
a_max_values = a[idx, I, J]
b_max_values = b[idx, I, J]

---

For the general case:

def argmax_to_max(arr, argmax, axis):
    """argmax_to_max(arr, arr.argmax(axis), axis) == arr.max(axis)"""
    new_shape = list(arr.shape)
    del new_shape[axis]

    grid = np.ogrid[tuple(map(slice, new_shape))]
    grid.insert(axis, argmax)

    return arr[tuple(grid)]

Quite a bit more awkward than such a natural operation should be, unfortunately.

For indexing a n dim array with a (n-1) dim array, we could simplify it a bit to give us the grid of indices for all axes, like so -

def all_idx(idx, axis):
    grid = np.ogrid[tuple(map(slice, idx.shape))]
    grid.insert(axis, idx)
    return tuple(grid)

Hence, use it to index into input arrays -

axis = 0
a_max_values = a[all_idx(idx, axis=axis)]
b_max_values = b[all_idx(idx, axis=axis)]

How can I select values along an axis of an nD array with an (n-1)D array of indices of that axis?

For the 2 and 1d case, this indexing works:

A[np.arange(J.shape[0]), J]

Which can be applied to more dimensions by reshaping to 2d (and back):

A.reshape(-1, A.shape[-1])[np.arange(np.prod(A.shape[:-1])).reshape(J.shape), J]

For 3d A this works:

A[np.arange(J.shape[0])[:,None], np.arange(J.shape[1])[None,:], J]

where the 1st 2 arange indices broadcast to the same dimension as J.

With functions in lib.index_tricks, this can be expressed as:

A[np.ogrid[0:J.shape[0],0:J.shape[1]]+[J]]
A[np.ogrid[slice(J.shape[0]),slice(J.shape[1])]+[J]]

or for multiple dimensions:

A[np.ix_(*[np.arange(x) for x in J.shape])+(J,)]
A[np.ogrid[[slice(k) for k in J.shape]]+[J]]

---

For small A and J (eg 2*3*4), J.choose(np.rollaxis(A,-1)) is faster. All of the extra time is in preparing the index tuple. np.ix_ is faster than np.ogrid.

np.choose has a size limit. At its upper end it is slower than ix_:

In [610]: Abig=np.arange(31*31).reshape(31,31)
In [611]: Jbig=np.arange(31)
In [612]: Jbig.choose(np.rollaxis(Abig,-1))
Out[612]: 
array([  0,  32,  64,  96, 128, 160, ... 960])

In [613]: timeit Jbig.choose(np.rollaxis(Abig,-1))
10000 loops, best of 3: 73.1 µs per loop
In [614]: timeit Abig[np.ix_(*[np.arange(x) for x in Jbig.shape])+(Jbig,)]
10000 loops, best of 3: 22.7 µs per loop
In [635]: timeit Abig.ravel()[Jbig+Abig.shape[-1]*np.arange(0,np.prod(Jbig.shape)).reshape(Jbig.shape) ]
10000 loops, best of 3: 44.8 µs per loop

I did similar indexing tests at https://stackoverflow.com/a/28007256/901925, and found that flat indexing was faster for much larger arrays (e.g. n0=1000). That's where I learned about the 32 limit for choice.

---

Related Topics