Numpy: создать массив элементов в одном массиве, которые находятся в диапазоне вокруг элемента другого массива

Question

это моя проблема:

Дано:
массива А имеет такую ​​форму: 12000,3 и массив В этом shape: 150,

Первый столбец A содержит значения времени, а B также содержит значения времени, отобранные с другой скоростью, поэтому они точно не совпадают.

Задача: Создать массив C с формой 150,3, который содержит строки массива A, где первый столбец находится в временном окне вокруг одной из временных точек в массиве B. Временное окно определяется временем «до» и времени «после»

Решение:

Он работает с 1D списки, используя список понимание, таких как: C = [е для е в а, если е> (B - раньше) и е < (B + after)]

Но попробуйте это с помощью массивов не работал.

Я начал простую по целочисленным массивам, пытаясь использовать логическую индексацию, но уже это не удается. Либо я получаю все False, либо сообщение об ошибке.

A = np.array([1,2,3,4,5,6]) 
B = np.array([1,3,5]) 

C = A[A in B] 
C = A[A in B.any] 
C = A[A == B] 

np.select делает то же самое.

В идеале это должно выглядеть следующим образом:

C = A[A > (B.any-before) and A < (B.any+after)] 

Большое спасибо за вашу помощь!

Answer 1

Хитрость здесь заключается в том, чтобы передавать массив A так, чтобы разные A и B могли транслироваться вместе, а затем уменьшали исходный массив, чтобы вернуть измерение B. Также стоит отметить, что умножение двух 1 вместе дает 1, а добавление 1 и нулей также дает 1. Кроме того, логическое значение любого положительного целого равно True.

A = A.reshape(A.size,1) 
B = np.ravel(B) 
bool_array = np.sum((A > B - before) * (A < B + after), axis=1).astype('bool') 
C = A[bool_array] 
A = np.ravel(A) 

Answer 2

Я использую это сейчас. Я должен петля через короткий массив, но он по-прежнему гораздо быстрее, чем список понимание:

C_all = [] 
for i in range(len(B): 
    #A and B must be np.arrays 
    bool_array = ((A > B - before) * (A < B + after), axis=1).astype('bool') 
    C = B[bool_array] 
    C_all.append(C) 

Answer 3

Это выглядит как проблема недавно я столкнулся (при выполнении слоистых передискретизаций, если кто-то заботится). Решение, которое я закончил использовать, было примерно таким:

EDIT: в зависимости от размера вашего массива, может быть гораздо быстрее (и более pythonic) использовать трансляцию, поэтому я включил это и небольшое испытание для сравнения два.

import scipy 
import timeit 
import itertools 

def AinBiterative(A,B): 
    M = B.shape[0] 
    N = A.shape[0] 

    which = scipy.arange(M)*(N-1)/M 
    # which will contain which elements from A are selected by B, this is a first guess 
    shift = (B >= A[:-1][which,0]).astype(int) - (B <= A[1:][which,0]).astype(int) 
    while not (shift == 0).all(): 
     which+= shift 
     shift = (B >= A[:-1][which,0]).astype(int) - (B <= A[1:][which,0]).astype(int) 

    return A[which,:] 

def AinBbroadcast(A,B): 
    comp = A[:,0][scipy.newaxis,:] < B[:,scipy.newaxis] 
    return A[comp.sum(axis=-1)-1,:] 

def wrapper(func, *args, **kwargs): 
    def wrapped(): 
     return func(*args, **kwargs) 
    return wrapped 

Ms = scipy.power(10, range(1,5)) 
Ns = scipy.power(10, range(1,6)) 

for M, N in itertools.product(Ms, Ns): 
    A = scipy.int_(N*scipy.randn(N,3)) 
    A[:,0] = scipy.arange(N, dtype=int) 
    B = scipy.rand(M)*(N-1) 
    B = scipy.sort(B) 

    wrapped = wrapper(AinBiterative, A, B) 
    iterative_time = timeit.timeit(wrapped, number=10) 

    wrapped = wrapper(AinBbroadcast, A, B) 
    broadcast_time = timeit.timeit(wrapped, number=10) 

    print '{:<6d}, {:<6d}: {:1.2e} vs {:1.2e}'.format(M, N, iterative_time, broadcast_time) 

я предположил, что все B вы ищете на самом деле в пределах диапазона, натянутой на А (отсюда N-1) и что А и В являются монотонными. Не обязательно быть целым, ни одно время не должно быть равномерно распределенным.

Answer 4

Вы можете использовать numpy's broadcasting, чтобы выполнить это.

>>> A = np.random.randn(12000, 3) 
>>> B = np.random.randn(150) 
>>> Ad = A[:, 0] 

Здесь Ad вектор, содержащий от времени A. Numpy имеет очень мощное радиовещание,

Ad[:, None] - B[None, :] 

будет возвращать 12000 x 150 массив, сравнивая каждый A каждые B. Это очень полезно и быстро, если у вас достаточно памяти для хранения этого массива. 12000 x 150 довольно маленький, так что это прекрасно, но имейте в виду, что размеры растут много.

Таким образом, используя вышеупомянутую идею вещания, вы можете транслировать свои правила на всем пространстве продукта:

>>> before = 0.001 
>>> after = 0.001 
>>> rule1 = (Ad[:, None] > B[None, :] - before) # 12000 x 150 
>>> rule2 = (Ad[:, None] < B[None, :] + after) # 12000 x 150 
>>> mask = rule1 & rule2      # 12000 x 150 
>>> valid = mask.any(axis=1) 

Обратите внимание, что rule1 и rule2 плотные сравнения (невероятно быстро в NumPy), что сравнить выше правила между каждым элементом Ad и каждым элементом B.

Действительный здесь проверяет, находится ли каждая запись A внутри не менее 1 окно B. Обратите внимание: вы можете использовать 1-строчную маску без создания временных rule1 и rule2, просто хотите сделать код более удобочитаемым.

>>> C = A[valid] 
>>> C.shape, A.shape 
(943, 3), (12000, 3) 

Мы отфильтровали 943 элемента из 12000!

---

Функцию из него:

def window_filter(A, B, before=0.001, after=0.001): 
    """Returns indexes from A that are within a window of B""" 
    rule1 = (A[:, None] > B[None, :] - before) 
    rule2 = (A[:, None] < B[None, :] + after) 
    mask = rule1 & rule2 
    valid = mask.any(axis=1) 
    return valid 

>>> C = A[window_filter(A[:, 0], B)] 
>>> C.shape, A.shape 
(943, 3), (12000, 3) 

А также довольно быстро (4ms):

>>> %timeit window_filter(A[:, 0], B) 
100 loops, best of 3: 4.53 ms per loop 

---

Кроме того, если вы хотите дополнительно указать более ограничивающего, например elements from A that are inside **only 1** window of B, вы можете заменить valid как

valid = (mask.sum(axis=1) == 1)