Tengo una matriz de tamaño nudoso 4D (98,359,256,269) que quiero umbral. En este momento, tengo dos listas separadas que mantienen las coordenadas de las primeras 2 dimensiones y las últimas 2 dimensiones. (mag_ang para las 2 primeras dimensiones e índices para las 2 últimas).

Tamaño de los índices: (61821,2)

Tamaño de mag_ang: (35182,2)

Actualmente, mi código se ve así:

inner_points = []

for k in indices:
    x = k[0]
    y = k[1]
    for i,ctr in enumerate(mag_ang):
        mag = ctr[0]
        ang = ctr[1]
        if X[mag][ang][x][y] > 10:
            inner_points.append((y,x))

Este código funciona pero es bastante lento y me pregunto si hay alguna forma más pitónica / más rápida de hacer esto.

0
sabrinazuraimi 6 oct. 2019 a las 10:30

4 respuestas

La mejor respuesta

(EDITAR: se agregó un segundo método alternativo)

Utilice numpy indexación de múltiples arreglos:

import time

import numpy as np

n_mag, n_ang, n_x, n_y = 10, 12, 5, 6
shape = n_mag, n_ang, n_x, n_y
X = np.random.random_sample(shape) * 20

nb_indices = 100 # 61821
indices = np.c_[np.random.randint(0, n_x, nb_indices), np.random.randint(0, n_y, nb_indices)]

nb_mag_ang = 50 # 35182
mag_ang = np.c_[np.random.randint(0, n_mag, nb_mag_ang), np.random.randint(0, n_ang, nb_mag_ang)]

# original method
inner_points = []
start = time.time()
for x, y in indices:
    for mag, ang in mag_ang:
        if X[mag][ang][x][y] > 10:
            inner_points.append((y, x))
end = time.time()
print(end - start)

# faster method 1:
inner_points_faster1 = []
start = time.time()
for x, y in indices:
    if np.any(X[mag_ang[:, 0], mag_ang[:, 1], x, y] > 10):
        inner_points_faster1.append((y, x))
end = time.time()
print(end - start)

# faster method 2:
start = time.time()
# note: depending on the real size of mag_ang and indices, you may wish to do this the other way round ?
found = X[:, :, indices[:, 0], indices[:, 1]][mag_ang[:, 0], mag_ang[:, 1], :] > 10
# 'found' shape is (nb_mag_ang x nb_indices)
assert found.shape == (nb_mag_ang, nb_indices)
matching_indices_mask = found.any(axis=0)
inner_points_faster2 = indices[matching_indices_mask, :]
end = time.time()
print(end - start)

# finally assert equality of findings
inner_points = np.unique(np.array(inner_points))
inner_points_faster1 = np.unique(np.array(inner_points_faster1))
inner_points_faster2 = np.unique(inner_points_faster2)
assert np.array_equal(inner_points, inner_points_faster1)
assert np.array_equal(inner_points, inner_points_faster2)

Rendimientos

0.04685807228088379
0.0
0.0

(por supuesto, si aumenta la forma, el tiempo no será cero para el segundo y el tercero)

Nota final: aquí utilizo "único" al final, pero tal vez sería conveniente hacerlo por adelantado para las matrices indices y mag_ang (excepto si está seguro de que ya son únicas)

1
smarie 6 oct. 2019 a las 19:20

Aquí hay algunas formas vecorizadas de aprovechar broadcasting -

thresh = 10
mask = X[mag_ang[:,0],mag_ang[:,1],indices[:,0,None],indices[:,1,None]]>thresh
r = np.where(mask)[0]
inner_points_out = indices[r][:,::-1]

Para matrices más grandes, podemos comparar primero y luego indexar para obtener la máscara:

mask = (X>thresh)[mag_ang[:,0],mag_ang[:,1],indices[:,0,None],indices[:,1,None]]

Si solo está interesado en las coordenadas únicas de indices, use la máscara directamente -

inner_points_out = indices[mask.any(1)][:,::-1]

Para matrices grandes, también podemos aprovechar los núcleos múltiples con el módulo numexpr.

Por lo tanto, primero importa el módulo -

import numexpr as ne

Luego, reemplace (X>thresh) con ne.evaluate('X>thresh') en los cálculos enumerados anteriormente.

1
Divakar 6 oct. 2019 a las 08:58

Use numpy directamente. Si indices y mag_ang son matrices vacías de dos columnas cada una para la coordenada apropiada:

(x, y), (mag, ang) = indices.T, mag_ang.T
index_matrix = np.meshgrid(mag, ang, x, y).T.reshape(-1,4)
inner_mag, inner_ang, inner_x, inner_y = np.where(X[index_matrix] > 10)

Ahora las variables inner... tienen matrices para cada coordenada. Para obtener una lista única de pares, puede comprimir los inner_y y inner_x.

1
kabanus 6 oct. 2019 a las 08:25

Utilice np.where

inner = np.where(X > 10)
a, b, x, y = zip(*inner)
inner_points = np.vstack([y, x]).T
0
meTchaikovsky 6 oct. 2019 a las 08:08
58255155