Ускорение вложенных циклов с NumPy

Question

Я пытаюсь решить проблему динамического программирования, и я придумал простой алгоритм на основе цикла, который заполняет двумерный массив на основе ряда if операторов, таких как:

s = # some string of size n
opt = numpy.zeros(shape=(n, n))

for j in range(0, n):
    for i in range(j, -1, -1):
        if j - i == 0:
            opt[i, j] = 1
        elif j - i == 1:
            opt[i, j] = 2 if s[i] == s[j] else 1
        elif s[i] == s[j] and opt[i + 1, j - 1] == (j - 1) - (i + 1) + 1:
            opt[i, j] = 2 + opt[i + 1, j - 1]
        else:
            opt[i, j] = max(opt[i + 1, j], opt[i, j - 1], opt[i + 1, j - 1])

К сожалению, этот код очень медленный для больших значений N. Я обнаружил, что гораздо лучше использовать встроенные функции, такие как numpy.where и numpy.fill, для заполнения значений массива, а не для циклов,но я изо всех сил пытаюсь найти какие-либо примеры, которые объясняют, как эти функции (или другие оптимизированные numpy методы) можно заставить работать с серией if операторов, как это делает мой алгоритм.Как можно переписать приведенный выше код со встроенными библиотеками numpy, чтобы сделать его лучше оптимизированным для Python?

Answer 1

Я не думаю, что np.where и np.fill могут решить вашу проблему.np.where используется для возврата элементов массива numpy, которые удовлетворяют определенному условию, но в вашем случае это условие NOT для VALUES массива numpy, но для значенийиз переменных i и j.

Для вашего конкретного вопроса я бы порекомендовал использовать Cython для оптимизации вашего кода специально для больших значений N. Cython - это, по сути, интерфейс между Python и C. Красота Cython заключается в том, что онпозволяет вам сохранить ваш синтаксис Python, но оптимизировать его, используя структуры Си.Он позволяет вам определять типы переменных в стиле C, чтобы ускорить ваши вычисления.Например, определение i и j как целых чисел с использованием Cython значительно ускорит процесс, поскольку типы i и j проверяются на каждой итерации цикла.

Кроме того, Cython позволит вам определять классические, быстрые, 2Dмассивы с использованием C. Затем вы можете использовать указатели для быстрого доступа к элементу этого двумерного массива вместо использования массивов.В вашем случае opt будет тем 2D-массивом.

Answer 2

Вот векторизованное решение.

Создает диагональные представления в выходном массиве, которые позволяют нам накапливать в диагональном направлении.

Пошаговое объяснение:

• оцените s [i] == s [j] в диагональном виде.

• оставьте только те, которые подключены к основной или первой субдиагональнойсерией истин в верхнем правом нижнем левом направлении

• замените все истины на 2 с, кроме главной диагонали, которая вместо этого получает 1 с;возьмите накопленную сумму в нижнем левом направлении в верхнем правом направлении

• наконец, возьмите совокупный максимум в нижнем восходящем и левом правом направлении

Как этоне совсем очевидно, что это то же самое, что и зацикленный код, который я тестировал на довольно многих примерах (с использованием функции stresstest ниже), и это кажется правильным.И примерно в 7 раз быстрее для строк средней длины (1–100 символов).

import numpy as np

def loopy(s):
    n = len(s)
    opt = np.zeros(shape=(n, n), dtype=int)
    for j in range(0, n):
        for i in range(j, -1, -1):
            if j - i == 0:
                opt[i, j] = 1
            elif j - i == 1:
                opt[i, j] = 2 if s[i] == s[j] else 1
            elif s[i] == s[j] and opt[i + 1, j - 1] == (j - 1) - (i + 1) + 1:
                opt[i, j] = 2 + opt[i + 1, j - 1]
            else:
                opt[i, j] = max(opt[i + 1, j], opt[i, j - 1], opt[i + 1, j - 1])
    return opt

def vect(s):
    n = len(s)
    h = (n+1) // 2
    s = np.array([s, s]).view('U1').ravel()
    opt = np.zeros((n+2*h-1, n+2*h-1), int)
    y, x = opt.strides
    hh = np.lib.stride_tricks.as_strided(opt[h-1:, h-1:], (2, h, n), (x, x-y, x+y))
    p, o, c = np.ogrid[:2, :h, :n]
    hh[...] = 2 * np.logical_and.accumulate(s[c+o+p] == s[c-o], axis=1)
    np.einsum('ii->i', opt)[...] = 1
    hh[...] = hh.cumsum(axis=1)
    opt = np.maximum.accumulate(opt[-h-1:None if h == 1 else h-2:-1, h-1:-h], axis=0)[::-1]
    return np.maximum.accumulate(opt, axis=1)

def stresstest(n=100):
    from string import ascii_lowercase
    import random
    from timeit import timeit
    Tv, Tl = 0, 0
    for i in range(n):
        s = ''.join(random.choices(ascii_lowercase[:random.randint(2, 26)], k=random.randint(1, 100)))
        print(s, end=' ')
        assert np.all(vect(s) == loopy(s))
        Tv += timeit(lambda: vect(s), number=10)
        Tl += timeit(lambda: loopy(s), number=10)
    print()
    print(f"total time loopy {Tl}, vect {Tv}")

Демонстрация:

>>> stresstest(20)
caccbbdbcfbfdcacebbecffacabeddcfdededeeafaebeaeedaaedaabebfacbdd fckjhrmupcqmihlohjog dffffgalbdbhkjigladhgdjaaagelddehahbbhejkibdgjhlkbcihiejdgidljfalfhlaglcgcih eacdebdcfcdcccaacfccefbccbced agglljlhfj mvwlkedblhvwbsmvtbjpqhgbaolnceqpgkhfivtbkwgbvujskkoklgforocj jljiqlidcdolcpmbfdqbdpjjjhbklcqmnmkfckkch ohsxiviwanuafkjocpexjmdiwlcmtcbagksodasdriieikvxphksedajwrbpee mcwdxsoghnuvxglhxcxxrezcdkahpijgujqqrqaideyhepfmrgxndhyifg omhppjaenjprnd roubpjfjbiafulerejpdniniuljqpouimsfukudndgtjggtbcjbchhfcdhrgf krutrwnttvqdemuwqwidvntpvptjqmekjctvbbetrvehsgxqfsjhoivdvwonvjd adiccabdbifigeigdfaieecceciaghadiaigibehdaichfibeaggcgdciahfegefigghgebhddciaei llobdegpmebejvotsr rtnsevatjvuowmquaulfmgiwsophuvlablslbwrpnhtekmpphsenarhrptgbjvlseeqstewjgfhopqwgmcbcihljeguv gcjlfihmfjbkdmimjknamfbahiccbhnceiahbnhghnlleimmieglgbfjbnmemdgddndhinncegnmgmfmgahhhjkg nhbnfhp cyjcygpaaeotcpwfhnumcfveq snyefmeuyjhcglyluezrx hcjhejhdaejchedbce 
total time loopy 0.2523909523151815, vect 0.03500175685621798

Answer 3

Ваши операторы if и левые части ваших операторов присваивания содержат ссылки на массив, который вы изменяете в цикле.Это означает, что не будет общего способа преобразовать ваш цикл в операции с массивами.Таким образом, вы застряли с каким-то циклом for.

Если бы вместо этого у вас был более простой цикл:

for j in range(0, n):
    for i in range(j, -1, -1):
        if j - i == 0:
            opt[i, j] = 1
        elif j - i == 1:
            opt[i, j] = 2
        elif s[i] == s[j]:
            opt[i, j] = 3
        else:
            opt[i, j] = 4

, вы могли бы создавать логические массивы (используя некоторые вещание ), которые представляют три ваших условия:

import numpy as np

# get arrays i and j that represent the row and column indices
i,j = np.ogrid[:n, :n]
# construct an array with the characters from s
sarr = np.fromiter(s, dtype='U1').reshape(1, -1)

cond1 = i==j             # result will be a bool arr with True wherever row index equals column index
cond2 = j==i+1           # result will be a bool arr with True wherever col index equals (row index + 1)
cond3 = sarr==sarr.T     # result will be a bool arr with True wherever s[i]==s[j]

Затем вы можете использовать numpy.select для построения желаемого opt:

opt = np.select([cond1, cond2, cond3], [1, 2, 3], default=4)

Для n=5и s='abbca', это даст:

array([[1, 2, 4, 4, 3],
       [4, 1, 2, 4, 4],
       [4, 3, 1, 2, 4],
       [4, 4, 4, 1, 2],
       [3, 4, 4, 4, 1]])