Numpy: поиск количества различных значений из ассоциаций посредством биннинга

Необходимое условие

Этот вопрос является продолжением этого сообщения. Итак, часть введения проблемы будет похожа на этот пост.

Проблема

Допустим, result — это двумерный массив, а values — это одномерный массив. values содержит некоторые значения, связанные с каждым элементом в result. Отображение элемента из values в result сохраняется в x_mapping и y_mapping. Позиция в result может быть связана с разными значениями. Пара (x,y) из x_mapping и y_mapping связана с results[-y,x]. Мне нужно найти уникальное количество значений, сгруппированных по ассоциациям.

Пример для лучшего пояснения.

result массив:

[[ 0.,  0.],
[ 0.,  0.],
[ 0.,  0.],
[ 0.,  0.]]

values массив:

[ 1.,  2.,  1.,  1.,  5.,  6.,  7.,  1.]

Примечание. Здесь массивы result и values имеют одинаковое количество элементов. Но это может быть не так. Между размерами вообще нет никакой связи.

x_mapping и y_mapping имеют отображения из 1D values в 2D result. Размеры x_mapping, y_mapping и values будут одинаковыми.

x_mapping - [0, 1, 0, 0, 0, 0, 0, 0]

y_mapping - [0, 3, 2, 2, 0, 3, 2, 0]

Здесь 1-е значение (значения [0]), 5-е значение (значения [4]) и 8-е значение (значения [7]) имеют x как 0 и y как 0 (x_mapping [0] и y_mapping [0]) и, следовательно, связаны с результатом [0, 0]. Если мы вычислим количество различных значений из этой группы (1,5,1), в результате мы получим 2. @WarrenWeckesser Давайте посмотрим, как пара [1, 3] (x, y) из x_mapping и y_mapping способствует results. Поскольку существует только одно значение, т.е. 2, связанное с этой конкретной группой, results[-3,1] будет иметь единицу, так как количество различных значений, связанных с этой ячейкой, равно единице.

Другой пример. Давайте вычислим значение results[-1,1]. Из сопоставлений, поскольку с ячейкой не связано никакого значения, значение results[-1,1] будет равно нулю.

Точно так же позиция [-2, 0] в results будет иметь значение 2.

Обратите внимание, что если связи вообще нет, то значение по умолчанию для result будет равно нулю.

result после вычисления,

[[ 2.,  0.],
[ 1.,  1.],
[ 2.,  0.],
[ 0.,  0.]]

Текущее рабочее решение

Используя ответ от @Divakar, я смог найти рабочее решение.

x_mapping = np.array([0, 1, 0, 0, 0, 0, 0, 0])
y_mapping = np.array([0, 3, 2, 2, 0, 3, 2, 0])
values = np.array([ 1.,  2.,  1.,  1.,  5.,  6.,  7.,  1.], dtype=np.float32)
result = np.zeros([4, 2], dtype=np.float32) 

m,n = result.shape
out_dtype = result.dtype
lidx = ((-y_mapping)%m)*n + x_mapping

sidx = lidx.argsort()
idx = lidx[sidx]
val = values[sidx]

m_idx = np.flatnonzero(np.r_[True,idx[:-1] != idx[1:]])
unq_ids = idx[m_idx]

r_res = np.zeros(m_idx.size, dtype=np.float32)
for i in range(0, m_idx.shape[0]):
    _next = None
    arr = None
    if i == m_idx.shape[0]-1:
        _next = val.shape[0]
    else:
        _next = m_idx[i+1]
    _start = m_idx[i]

    if _start >= _next:
        arr = val[_start]
    else:
        arr = val[_start:_next]
    r_res[i] = np.unique(arr).size
result.flat[unq_ids] = r_res

Вопрос

Теперь приведенное выше решение требует 15 мс для работы со значениями 19943. Я ищу способ вычислить результат быстрее. Есть ли более эффективный способ сделать это?

Примечание

Я использую Numpy версии 1.14.3 с Python 3.5.2.

Правки

Благодаря @WarrenWeckesser, указав, что я не объяснил, как элемент в results связан с (x,y) из сопоставлений. Я обновил сообщение и добавил примеры для ясности.


python numpy arrays
person tpk    schedule 28.11.2018    source источник
comment
У меня возникли проблемы с согласованием вашего описания того, как вы вычислили result[0,0], с остальными значениями в result (которые генерируются кодом, который, как вы говорите, работает). Например, в массивах x_mapping и y_mapping пара (x, y) [1, 3] встречается один раз. Насколько я понимаю, это индексы столбцов и строк в result. Так почему же result[3, 1] не равно 1? И в вычисленном result у вас есть result[1, 0] = 1 и result[1, 1] = 1, но ни одна из пар (x, y) [0, 1] и [1, 1] не встречается в массивах отображения.   -  person Warren Weckesser    schedule 28.11.2018
comment
@WarrenWeckesser, спасибо, что указали на это. Прошу прощения за то, что не добавил подробностей о том, как пара (x,y) связана с элементами в results. Каждая пара (x,y) связана с results[-y,x]. Я обновил сообщение и добавил примеры для ясности. Спасибо.   -  person tpk    schedule 28.11.2018

Ответы (1)


arrow_upward
1
arrow_downward

Вот одно решение

import numpy as np

x_mapping = np.array([0, 1, 0, 0, 0, 0, 0, 0])
y_mapping = np.array([0, 3, 2, 2, 0, 3, 2, 0])
values = np.array([ 1.,  2.,  1.,  1.,  5.,  6.,  7.,  1.], dtype=np.float32)
result = np.zeros([4, 2], dtype=np.float32)

# Get flat indices
idx_mapping = np.ravel_multi_index((-y_mapping, x_mapping), result.shape, mode='wrap')
# Sort flat indices and reorders values accordingly
reorder = np.argsort(idx_mapping)
idx_mapping = idx_mapping[reorder]
values = values[reorder]
# Get unique values
val_uniq = np.unique(values)
# Find where each unique value appears
val_uniq_hit = values[:, np.newaxis] == val_uniq
# Find reduction indices (slices with the same flat index)
reduce_idx = np.concatenate([[0], np.nonzero(np.diff(idx_mapping))[0] + 1])
# Reduce slices
reduced = np.logical_or.reduceat(val_uniq_hit, reduce_idx)
# Count distinct values on each slice
counts = np.count_nonzero(reduced, axis=1)
# Put counts in result
result.flat[idx_mapping[reduce_idx]] = counts

print(result)
# [[2. 0.]
#  [1. 1.]
#  [2. 0.]
#  [0. 0.]]

Этот метод требует больше памяти (O(len(values) * len(np.unique(values)))), но небольшой тест по сравнению с вашим исходным решением показывает значительное ускорение (хотя это зависит от фактического размера проблемы):

import numpy as np

np.random.seed(100)
result = np.zeros([400, 200], dtype=np.float32)
values = np.random.randint(100, size=(20000,)).astype(np.float32)
x_mapping = np.random.randint(result.shape[1], size=values.shape)
y_mapping = np.random.randint(result.shape[0], size=values.shape)

res1 = solution_orig(x_mapping, y_mapping, values, result)
res2 = solution(x_mapping, y_mapping, values, result)
print(np.allclose(res1, res2))
# True

# Original solution
%timeit solution_orig(x_mapping, y_mapping, values, result)
# 76.2 ms ± 623 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

# This solution
%timeit solution(x_mapping, y_mapping, values, result)
# 13.8 ms ± 51.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Полный код тестовых функций:

import numpy as np

def solution(x_mapping, y_mapping, values, result):
    result = np.array(result)
    idx_mapping = np.ravel_multi_index((-y_mapping, x_mapping), result.shape, mode='wrap')
    reorder = np.argsort(idx_mapping)
    idx_mapping = idx_mapping[reorder]
    values = values[reorder]
    val_uniq = np.unique(values)
    val_uniq_hit = values[:, np.newaxis] == val_uniq
    reduce_idx = np.concatenate([[0], np.nonzero(np.diff(idx_mapping))[0] + 1])
    reduced = np.logical_or.reduceat(val_uniq_hit, reduce_idx)
    counts = np.count_nonzero(reduced, axis=1)
    result.flat[idx_mapping[reduce_idx]] = counts
    return result

def solution_orig(x_mapping, y_mapping, values, result):
    result = np.array(result)
    m,n = result.shape
    out_dtype = result.dtype
    lidx = ((-y_mapping)%m)*n + x_mapping

    sidx = lidx.argsort()
    idx = lidx[sidx]
    val = values[sidx]

    m_idx = np.flatnonzero(np.r_[True,idx[:-1] != idx[1:]])
    unq_ids = idx[m_idx]

    r_res = np.zeros(m_idx.size, dtype=np.float32)
    for i in range(0, m_idx.shape[0]):
        _next = None
        arr = None
        if i == m_idx.shape[0]-1:
            _next = val.shape[0]
        else:
            _next = m_idx[i+1]
        _start = m_idx[i]

        if _start >= _next:
            arr = val[_start]
        else:
            arr = val[_start:_next]
        r_res[i] = np.unique(arr).size
    result.flat[unq_ids] = r_res
    return result
person jdehesa    schedule 28.11.2018
comment
Спасибо за ответы. Я изменил существующее решение с вашей логикой использования np.logical_or.reduceat. Это намного быстрее. Спасибо. - person tpk; 29.11.2018