Digamos que tengo una matriz Numpy de N = 10 números flotantes aleatorios:
import numpy as np np.random.seed(99) N = 10 arr = np.random.uniform(0., 10., size=(N,)) print(arr) out[1]: [6.72278559 4.88078399 8.25495174 0.31446388 8.08049963 5.6561742 2.97622499 0.46695721 9.90627399 0.06825733] Quiero encontrar todos los pares únicos de números que no sean diferentes entre sí más que una tolerancia tol = 1. (es decir, diferencia absoluta <= 1). Específicamente, quiero obtener todos los pares únicos de índices. Los índices de cada par cercano deben ordenarse y todos los pares cercanos deben ordenarse por el primer índice. Me las arreglé para escribir el siguiente código de trabajo:
def all_close_pairs(arr, tol=1.): res = set() for i, x1 in enumerate(arr): for j, x2 in enumerate(arr): if i == j: continue if np.isclose(x1, x2, rtol=0., atol=tol): res.add(tuple(sorted([i, j]))) res = np.array(list(res)) return res[res[:,0].argsort()] print(all_close_pairs(arr, tol=1.)) out[2]: [[1 5] [2 4] [3 7] [3 9] [7 9]] Sin embargo, en realidad tengo una matriz de N = 1000 números, y mi código se vuelve extremadamente lento debido a los bucles for anidados. Creo que hay formas mucho más eficientes de hacer esto con la vectorización Numpy. ¿Alguien sabe la forma más rápida de hacer esto en Numpy?
Una solución eficiente es ordenar primero los valores de entrada usando index = np.argsort() . Luego, puede generar la matriz ordenada usando arr[index] y luego iterar sobre los valores de cierre en un tiempo casi lineal si el número de pares es pequeño en una matriz contigua rápida. Si el número de pares es grande, entonces la complejidad es cuadrática debido al número cuadrático de pares generados. La complejidad resultante es: O(n log n + m) donde n es el tamaño de la matriz de entrada y m es el número de pares producidos.
Para encontrar valores cercanos entre sí, una forma eficiente es iterar sobre el valor usando Numba . De hecho, si bien podría ser posible en Numpy, es probable que no sea eficiente debido a la cantidad variable de valores que se compararán. Aquí hay una implementación:
import numba as nb @nb.njit('int32[:,::1](float64[::1], float64)') def findCloseValues(arr, tol): res = [] for i in range(arr.size): val = arr[i] # Iterate over the close numbers (only once) for j in range(i+1, arr.size): # Sadly neither np.isclose or np.abs are implemented in Numba so far if max(val, arr[j]) - min(val, arr[j]) >= tol: break res.append((i, j)) if len(res) == 0: # No pairs: we need to help Numpy to know the shape return np.empty((0, 2), dtype=np.int32) return np.array(res, dtype=np.int32) Finalmente, los índices deben actualizarse para hacer referencia a los índices en la matriz no ordenada y no a la ordenada. Puedes hacerlo usando index[result] .
Aquí está el código resultante:
index = arr.argsort() result = findCloseValues(arr[index], 1.0) print(index[result])Aquí está el resultado (el orden no es el mismo que en la pregunta, pero puede ordenarlo si es necesario):
array([[9, 3], [9, 7], [3, 7], [1, 5], [4, 2]]) Si necesita un algoritmo más rápido, puede usar otro formato de salida: para cada valor de entrada puede proporcionar el rango mínimo/máximo de valores cercanos al valor de entrada objetivo. Para encontrar el rango, puede usar una búsqueda binaria (ver: np.searchsorted ) en la matriz ordenada. El algoritmo resultante se ejecuta en O(n log n) . Sin embargo, no puede obtener los índices en la matriz desordenada ya que el rango no sería contiguo.
Estos son los resultados de rendimiento en una entrada aleatoria con 1000 elementos y una tolerancia de 1,0 en mi máquina:
Reference implementation: ~17000 ms (x 1) Angelicos' implementation: 1773 ms (x ~10) Rivers' implementation: 122 ms (x 139) Rchome's implementation: 20 ms (x 850) Chris' implementation: 4.57 ms (x 3720) This implementation: 0.67 ms (x 25373)Un poco tarde pero una solución completamente numpy:
import numpy as np def close_enough( arr, tol = 1 ): result = np.where( np.triu(np.isclose( arr[ :, None ], arr[ None, : ], rtol = 0.0, atol = tol ), 1)) return np.swapaxes( result, 0, 1 )Ampliado para explicar lo que está sucediendo.
def close_enough( arr, tol = 1 ): bool_arr = np.isclose( arr[ :, None ], arr[ None, : ], rtol = 0.0, atol = tol ) # is_close generates a square array after comparing all elements with all elements. bool_arr = np.triu( bool_arr, 1 ) # Keep the upper right triangle, offset by 1 column. ie zero the main diagonal # and all elements below and to the left. result = np.where( bool_arr ) # Return the row and column indices for Trues return np.swapaxes( result, 0, 1 ) # Return the pairs in rows rather than columnsCon N = 1000, arr = una matriz de flotadores
%timeit close_enough( arr, tol = 1 ) 14.1 ms ± 28.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each) In [19]: %timeit all_close_pairs( arr, tol = 1 ) 54.3 ms ± 268 µs per loop (mean ± std. dev. of 7 runs, 10 loops each) (close_enough( arr, tol = 1) == all_close_pairs( arr, tol = 1 )).all() # TrueEsta es una solución con operaciones numpy puras. Parece bastante rápido en mi máquina, pero no sé qué tipo de velocidad estamos buscando.
def all_close_pairs(arr, tol=1.): N = arr.shape[0] # get indices in the array to consider using meshgrid pair_coords = np.array(np.meshgrid(np.arange(N), np.arange(N))).T # filter out pairs so we get indices in increasing order pair_coords = pair_coords[pair_coords[:, :, 0] < pair_coords[:, :, 1]] # compare indices in your array for closeness is_close = np.isclose(arr[pair_coords[:, 0]], arr[pair_coords[:, 1]], rtol=0, atol=tol) return pair_coords[is_close, :]