Расчет векторизации в матрице с взаимозависимыми значениями

Question

Я отслеживаю несколько дискретных временных рядов при нескольких временных разрешениях, в результате чего получается матрица SxRxB, где S - это количество временных рядов, R - количество различных разрешений, а B - буфер, т.е. сколько значений в каждой серии. помнит. Каждая серия дискретна и использует ограниченный диапазон натуральных чисел для представления своих значений. Я назову эти «символы» здесь.

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

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

Заранее спасибо!

def supports_loop(data, num_series, resolutions, buffer_size, vocab_size):
    # For small test matrices we can calculate the complete matrix without problems
    indices = []
    indices.append(xrange(num_series))
    indices.append(xrange(vocab_size))
    indices.append(xrange(num_series))
    indices.append(xrange(vocab_size))
    indices.append(xrange(resolutions))

    # This is huge! :/
    # dimensions:
    #   series and value for which we calculate,
    #   series and value which precedes that measurement,
    #   resolution
    ratios = np.full((num_series, vocab_size, num_series, vocab_size, resolutions), 0.0)

    for idx in itertools.product(*indices):
        s0, v0 = idx[0],idx[1]  # the series and symbol for which we calculate
        s1, v1 = idx[2],idx[3]  # the series and symbol which should precede the we're calculating for
        res = idx[4]

        # Find the positions where s0==v0
        found0 = np.where(data[s0, res, :] == v0)[0]
        if found0.size == 0:
            continue
        #print('found {}={} at {}'.format(s0, v0, found0))

        # Check how often s1==v1 right before s0==v0
        candidates = (s1, res, (found0 - 1 + buffer_size) % buffer_size)
        found01 = np.count_nonzero(data[candidates] == v1)
        if found01 == 0:
            continue

        print('found {}={} following {}={} at {}'.format(s0, v0, s1, v1, found01))
        # total01 = number of positions where either s0 or s1 is defined (i.e. >=0)
        total01 = len(np.argwhere((data[s0, res, :] >= 0) & (data[s1, res, :] >= 0)))
        ratio = (float(found01) / total01) if total01 > 0 else 0.0
        ratios[idx] = ratio

    return ratios


def stackoverflow_example(fnc):
    data = np.array([
        [[0, 0, 1],  # series 0, resolution 0
         [1, 3, 2]], # series 0, resolution 1

        [[2, 1, 2],  # series 1, resolution 0
         [3, 3, 3]], # series 1, resoltuion 1
    ])

    num_series = data.shape[0]
    resolutions = data.shape[1]
    buffer_size = data.shape[2]
    vocab_size = np.max(data)+1

    ratios = fnc(data, num_series, resolutions, buffer_size, vocab_size)
    coordinates = np.argwhere(ratios > 0.0)
    nz_values = ratios[ratios > 0.0]
    print(np.hstack((coordinates, nz_values[:,None])))
    print('0/0 precedes 0/0 in 1 out of 3 cases: {}'.format(np.isclose(ratios[0,0,0,0,0], 1.0/3.0)))
    print('1/2 precedes 0/0 in 2 out of 3 cases: {}'.format(np.isclose(ratios[0,0,1,2,0], 2.0/3.0)))

Ожидаемый результат (21 пара, 5 столбцов для координат с последующим найденным количеством):

[[0 0 0 0 0 1]
 [0 0 0 1 0 1]
 [0 0 1 2 0 2]
 [0 1 0 0 0 1]
 [0 1 0 2 1 1]
 [0 1 1 1 0 1]
 [0 1 1 3 1 1]
 [0 2 0 3 1 1]
 [0 2 1 3 1 1]
 [0 3 0 1 1 1]
 [0 3 1 3 1 1]
 [1 1 0 0 0 1]
 [1 1 1 2 0 1]
 [1 2 0 0 0 1]
 [1 2 0 1 0 1]
 [1 2 1 1 0 1]
 [1 2 1 2 0 1]
 [1 3 0 1 1 1]
 [1 3 0 2 1 1]
 [1 3 0 3 1 1]
 [1 3 1 3 1 3]]

В приведенном выше примере 0 в серии 0 следует за 2 в серии 1 в двух из трех случаев (поскольку буферы имеют круглую форму), поэтому отношение в [0, 0, 1, 2, 0] будет ~ 0,6666. Также серия 0, значение 0 следует за собой в одном из трех случаев, поэтому отношение при [0, 0, 0, 0, 0] будет ~ 0,3333. Есть и другие, которые> 0,0 тоже.

---

Я тестирую каждый ответ на двух наборах данных: крошечный (как показано выше) и более реалистичный (100 серий, 5 разрешений, 10 значений в серии, 50 символов).

Результаты

Answer        Time (tiny)     Time (huge)     All pairs found (tiny=21)
-----------------------------------------------------------------------
Baseline      ~1ms            ~675s (!)       Yes
Saedeas       ~0.13ms         ~1.4ms          No (!)
Saedeas2      ~0.20ms         ~4.0ms          Yes, +cross resolutions
Elliot_1      ~0.70ms         ~100s (!)       Yes
Elliot_2      ~1ms            ~21s (!)        Yes
Kuppern_1     ~0.39ms         ~2.4s (!)       Yes
Kuppern_2     ~0.18ms         ~28ms           Yes
Kuppern_3     ~0.19ms         ~24ms           Yes
David         ~0.21ms         ~27ms           Yes

Saedeas 2nd подход - явный победитель! Большое спасибо всем вам:)

Answer 1

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

Внесенные мною изменения были направлены на удаление конкатенации, которая происходит в конце функции run(). Это создавало новый массив и не нужно. Вместо этого мы создаем массив полного размера в начале и просто не используем последнюю строку до конца.

Другое изменение, которое я сделал, - то, что тайлинг single был немного неэффективным. Я заменил это на немного более быстрый код.

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

код ниже;

import numpy as np
import logging
import sys
import time
import itertools
import timeit


logging.basicConfig(stream=sys.stdout,
                    level=logging.DEBUG,
                    format='%(message)s')


def run():
    series = 2
    resolutions = 2
    buffer_len = 3
    symbols = range(50)

    #mat = np.random.choice(symbols, size=(series, resolutions, buffer_len))

    mat = np.array([
            [[0, 0, 1],  # series 0, resolution 0
             [1, 3, 2]],  # series 0, resolution 1
            [[2, 1, 2],  # series 1, resolution 0
             [3, 3, 3]],  # series 1, resoltuion 1
            # [[4, 5, 6, 10],
            #  [7, 8, 9, 11]],
        ])

    # logging.debug("Original:")
    # logging.debug(mat)

    start = time.time()
    index_mat = np.indices((series, resolutions, buffer_len))

    # This loop shifts all series but the one being looked at, and zips the
    # element being looked at with every other member of that row
    cross_pairs = np.empty((series, resolutions, buffer_len, series, 2), int)
    #cross_pairs = []
    right_shift_indices = [index_mat[0], index_mat[1], (index_mat[2] - 1) % buffer_len]

    for i in range(series):
        right_shift_indices[2][i] = (right_shift_indices[2][i] + 1) % buffer_len


        # create a new matrix from the modified indices
        mat_shifted = mat[right_shift_indices]
        mat_shifted_t = mat_shifted.T.reshape(-1, series)
        single = mat_shifted_t[:, i]

        #print np.tile(single,(series-1,1)).T
        #print single.reshape(-1,1).repeat(series-1,1)
        #print single.repeat(series-1).reshape(-1,series-1)

        mat_shifted_t = np.delete(mat_shifted_t, i, axis=1)

        #cross_pairs[i,:,:,:-1] = (np.dstack((np.tile(single, (mat_shifted_t.shape[1], 1)).T, mat_shifted_t))).reshape(resolutions, buffer_len, (series-1), 2, order='F')
        #cross_pairs[i,:,:,:-1] = (np.dstack((single.reshape(-1,1).repeat(series-1,1), mat_shifted_t))).reshape(resolutions, buffer_len, (series-1), 2, order='F')
        cross_pairs[i,:,:,:-1] = np.dstack((single.repeat(series-1).reshape(-1,series-1), mat_shifted_t)).reshape(resolutions, buffer_len, (series-1), 2, order='F')

        right_shift_indices[2][i] = (right_shift_indices[2][i] - 1) % buffer_len
        #cross_pairs.extend([zip(itertools.repeat(x[i]), np.append(x[:i], x[i+1:])) for x in mat_shifted_t])

    #consecutive_pairs = np.empty((series, resolutions, buffer_len, 2, 2), int)
    #print "1", consecutive_pairs.shape
    # tedious code to put this stuff in the right shape
    in_series_zips = np.stack([mat[:, :, :-1], mat[:, :, 1:]], axis=3)
    circular_in_series_zips = np.stack([mat[:, :, -1], mat[:, :, 0]], axis=2)
    # This creates the final array.
    # Index 0 is the preceding series
    # Index 1 is the resolution
    # Index 2 is the location in the buffer
    # Index 3 is for the first n-1 elements, the following series, and for the last element
    #         it's the next element of the Index 0 series
    # Index 4 is the index into the two element pair
    cross_pairs[:,:,:-1,-1] = in_series_zips
    cross_pairs[:,:,-1,-1] = circular_in_series_zips

    end = time.time()
    #logging.debug("Pairs encountered:")
    #logging.debug(pairs)
    logging.info("Elapsed: {}".format(end - start))

if __name__ == '__main__':
    run()

Answer 2

Для начала, вы делаете себе небольшую медвежью услугу, не вкладывая явно циклы for. Вы повторяете много усилий и ничего не экономите с точки зрения памяти. Когда цикл является вложенным, вы можете переместить некоторые вычисления с одного уровня на другой и выяснить, какие внутренние циклы можно векторизовать.

def supports_5_loop(data, num_series, resolutions, buffer_size, vocab_size):
    ratios = np.full((num_series, vocab_size, num_series, vocab_size, resolutions), 0.0)
    for res in xrange(resolutions):
        for s0 in xrange(num_series):
            # Find the positions where s0==v0
            for v0 in np.unique(data[s0, res]):
                # only need to find indices once for each series and value
                found0 = np.where(data[s0, res, :] == v0)[0]
                for s1 in xrange(num_series):
                    # Check how often s1==v1 right before s0==v0
                    candidates = (s1, res, (found0 - 1 + buffer_size) % buffer_size)
                    total01 = np.logical_or(data[s0, res, :] >= 0, data[s1, res, :] >= 0).sum()
                    # can skip inner loops if there are no candidates
                    if total01 == 0:
                        continue
                    for v1 in xrange(vocab_size):
                        found01 = np.count_nonzero(data[candidates] == v1)
                        if found01 == 0:
                            continue

                        ratio = (float(found01) / total01)
                        ratios[(s0, v0, s1, v1, res)] = ratio

    return ratios

В моменты времени вы увидите, что большая часть ускорения происходит из-за неповторяющихся усилий.

Как только вы создали вложенную структуру, вы можете начать смотреть на векторизацию и другие оптимизации.

def supports_4_loop(data, num_series, resolutions, buffer_size, vocab_size):
    # For small test matrices we can calculate the complete matrix without problems

    # This is huge! :/
    # dimensions:
    #   series and value for which we calculate,
    #   series and value which precedes that measurement,
    #   resolution
    ratios = np.full((num_series, vocab_size, num_series, vocab_size, resolutions), 0.0)

    for res in xrange(resolutions):
        for s0 in xrange(num_series):
            # find the counts where either s0 or s1 are present
            total01 = np.logical_or(data[s0, res] >= 0,
                                    data[:, res] >= 0).sum(axis=1)
            s1s = np.where(total01)[0]
            # Find the positions where s0==v0
            v0s, counts = np.unique(data[s0, res], return_counts=True)
            # sorting before searching will show gains as the datasets
            # get larger
            indarr = np.argsort(data[s0, res])
            i0 = 0
            for v0, count in itertools.izip(v0s, counts):
                found0 = indarr[i0:i0+count]
                i0 += count
                for s1 in s1s:
                    candidates = data[(s1, res, (found0 - 1) % buffer_size)]
                    # can replace the innermost loop with numpy functions
                    v1s, counts = np.unique(candidates, return_counts=True)
                    ratios[s0, v0, s1, v1s, res] = counts / total01[s1]

    return ratios

К сожалению, я мог по-настоящему векторизовать только самый внутренний цикл, и это только увеличило ускорение на 10%. Вне самого внутреннего цикла вы не можете гарантировать, что все векторы имеют одинаковый размер, поэтому вы не можете построить массив.

In [121]: (np.all(supports_loop(data, num_series, resolutions, buffer_size, vocab_size) == supports_5_loop(data, num_series, resolutions, buffer_size, vocab_size)))
Out[121]: True

In [122]: (np.all(supports_loop(data, num_series, resolutions, buffer_size, vocab_size) == supports_4_loop(data, num_series, resolutions, buffer_size, vocab_size)))
Out[122]: True
In [123]: %timeit(supports_loop(data, num_series, resolutions, buffer_size, vocab_size))
2.29 ms ± 73.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [124]: %timeit(supports_5_loop(data, num_series, resolutions, buffer_size, vocab_size))
949 µs ± 5.37 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [125]: %timeit(supports_4_loop(data, num_series, resolutions, buffer_size, vocab_size))
843 µs ± 3.21 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Answer 3

Хитрость, которая делает этот векторизованным, состоит в том, чтобы создать массив comb[i] = buffer1[i]+buffer2[i-1]*voc_size для каждой пары серий. Каждая комбинация получает уникальное значение в массиве. И можно найти комбинацию, выполнив v1[i] = comb[i] % voc_size, v2[i] = comb[i]//voc_size. Пока число серий не очень велико (я думаю, <10000), нет смысла делать какие-либо дальнейшие векторизации. </p>

def support_vectorized(data, num_series, resolutions, buffer_size, vocab_size):
    ratios = np.zeros((num_series, vocab_size, num_series, vocab_size, resolutions))
    prev = np.roll(data, 1, axis=2)  # Get previous values
    prev *= vocab_size  # To separate prev from data
    for i, series in enumerate(data):
        for j, prev_series in enumerate(prev):
            comb = series + prev_series
            for k, buffer in enumerate(comb):
                idx, counts = np.unique(buffer, return_counts=True)
                v = idx % vocab_size    
                v2 = idx // vocab_size
                ratios[i, v, j, v2, k] = counts/buffer_size
    return ratios

Если S или R велики, возможна полная векторизация, но при этом используется много памяти:

def row_unique(comb):
    comb.sort(axis=-1)
    changes = np.concatenate((
        np.ones((comb.shape[0], comb.shape[1], comb.shape[2], 1), dtype="bool"),
        comb[:, :,:, 1:] != comb[:, :, :, :-1]), axis=-1)
    vals = comb[changes]
    idxs = np.nonzero(changes)
    tmp = np.hstack((idxs[-1], 0))
    counts = np.where(tmp[1:], np.diff(tmp), comb.shape[-1]-tmp[:-1])
    return idxs, vals, counts


def supports_full_vectorized(data, num_series, resolutions, buffer_size, vocab_size):
    ratios = np.zeros((num_series, vocab_size, num_series, vocab_size, resolutions))
    prev = np.roll(data, 1, axis=2)*vocab_size
    comb = data + prev[:, None]  # Create every combination
    idxs, vals, counts = row_unique(comb)  # Get unique values and counts for each row
    ratios[idxs[1], vals % vocab_size, idxs[0], vals // vocab_size, idxs[2]] = counts/buffer_size
    return ratios

Однако для S=100 это медленнее, чем предыдущее решение. Промежуточная точка - держать цикл for над сериями, чтобы уменьшить использование памяти:

def row_unique2(comb):
    comb.sort(axis=-1)
    changes = np.concatenate((
        np.ones((comb.shape[0], comb.shape[1], 1), dtype="bool"),
        comb[:, :, 1:] != comb[:, :, :-1]), axis=-1)
    vals = comb[changes]
    idxs = np.nonzero(changes)
    tmp = np.hstack((idxs[-1], 0))
    counts = np.where(tmp[1:], np.diff(tmp), comb.shape[-1]-tmp[:-1])
    return idxs, vals, counts


def supports_half_vectorized(data, num_series, resolutions, buffer_size, vocab_size):
    prev = np.roll(data, 1, axis=2)*vocab_size
    ratios = np.zeros((num_series, vocab_size, num_series, vocab_size, resolutions))
    for i, series in enumerate(data):
        comb = series + prev
        idxs, vals, counts = row_unique2(comb)
        ratios[i, vals % vocab_size, idxs[0], vals // vocab_size, idxs[1]] = counts/buffer_size
    return ratios

Время выполнения для различных решений показывает, что support_half_vectorized является самым быстрым

In [41]: S, R, B, voc_size = (100, 5, 1000, 29)

In [42]: data = np.random.randint(voc_size, size=S*R*B).reshape((S, R, B))

In [43]: %timeit support_vectorized(data, S, R, B, voc_size)
1 loop, best of 3: 4.84 s per loop

In [44]: %timeit supports_full_vectorized(data, S, R, B, voc_size)
1 loop, best of 3: 5.3 s per loop

In [45]: %timeit supports_half_vectorized(data, S, R, B, voc_size)
1 loop, best of 3: 4.36 s per loop

In [46]: %timeit supports_4_loop(data, S, R, B, voc_size)
1 loop, best of 3: 36.7 s per loop

Answer 4

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

import numpy as np
import time
from collections import Counter

series = 2
resolutions = 2
buffer_len = 3
symbols = range(3)

#mat = np.random.choice(symbols, size=(series, resolutions, buffer_len)).astype('uint8')

mat = np.array([
        [[0, 0, 1],  # series 0, resolution 0
         [1, 3, 2]],  # series 0, resolution 1
        [[2, 1, 2],  # series 1, resolution 0
         [3, 3, 3]],  # series 1, resoltuion 1
    ])

start = time.time()
index_mat = np.indices(mat.shape)

right_shift_indices = np.roll(index_mat, -1, axis=3)
mat_shifted = mat[right_shift_indices[0], right_shift_indices[1], right_shift_indices[2]]

# These construct all the pairs directly
first_series = np.repeat(range(series), series*resolutions*buffer_len)
second_series = np.tile(np.repeat(range(series), resolutions*buffer_len), series)
res_loop = np.tile(np.repeat(range(resolutions), buffer_len), series*series)
mat_unroll = np.repeat(mat, series, axis=0)
shift_unroll = np.tile(mat_shifted, series)

# Constructs the pairs
pairs = zip(np.ravel(first_series),
            np.ravel(second_series),
            np.ravel(res_loop),
            np.ravel(mat_unroll),
            np.ravel(shift_unroll))

pair_time = time.time() - start
results = Counter(pairs)
end = time.time() - start

print("Mat: {}").format(mat)
print("Pairs: {}").format(results)
print("Number of Pairs: {}".format(len(pairs)))
print("Pair time is: {}".format(pair_time))
print("Count time is: {}".format(end-pair_time))
print("Total time is: {}".format(end))

Основная идея состояла в том, чтобы циклически сдвигать каждый буфер на соответствующую величину в зависимости от того, каким временным рядом он был (я думаю, именно этим занимался ваш текущий код). Затем я могу сгенерировать все пары символов, просто сжав списки со смещением на 1 вдоль оси серии.

Пример вывода:

Mat: [[[0 0 1]
  [1 3 2]]

 [[2 1 2]
  [3 3 3]]]
Pairs: Counter({(1, 1, 1, 3, 3): 3, (1, 0, 0, 2, 0): 2, (0, 0, 0, 0, 0): 1, (1, 1, 0, 2, 2): 1, (1, 1, 0, 2, 1): 1, (0, 1, 0, 0, 2): 1, (1, 0, 1, 3, 3): 1, (0, 0, 1, 1, 3): 1, (0, 0, 1, 3, 2): 1, (1, 0, 0, 1, 1): 1, (0, 1, 0, 0, 1): 1, (0, 1, 1, 2, 3): 1, (0, 1, 0, 1, 2): 1, (1, 1, 0, 1, 2): 1, (0, 1, 1, 3, 3): 1, (1, 0, 1, 3, 2): 1, (0, 0, 0, 0, 1): 1, (0, 1, 1, 1, 3): 1, (0, 0, 1, 2, 1): 1, (0, 0, 0, 1, 0): 1, (1, 0, 1, 3, 1): 1})
Number of Pairs: 24
Pair time is: 0.000135183334351
Count time is: 5.10215759277e-05
Total time is: 0.000186204910278

Редактировать: Истинная последняя попытка. Полностью векторизовано.