Как найти подпоследовательность в строке

def is_subseq(x, y):
    it = iter(y)
    return all(any(c == ch for c in it) for ch in x)

assert is_subseq('india', 'indonesia')
assert is_subseq('oman', 'romania')
assert is_subseq('mali', 'malawi')
assert not is_subseq('mali', 'banana')
assert not is_subseq('ais', 'indonesia')
assert not is_subseq('ca', 'abc')

Also works for any iterables:

assert is_subseq(['i', 'n', 'd', 'i', 'a'],
                 ['i', 'n', 'd', 'o', 'n', 'e', 's', 'i', 'a'])

UPDATE

Stefan Pochmann suggested this.

def is_subseq(x, y):
    it = iter(y)
    return all(c in it for c in x)

Both versions uses iterators; Iterator yields items that was not yielded in previous iteration.

For example:

>>> it = iter([1,2,3,4])
>>> for x in it:
...     print(x)
...     break
...
1
>>> for x in it:  # `1` is yielded in previous iteration. It's not yielded here.
...     print(x)
...
2
3
4

Нахождение максимальной общей подпоследовательности

В настоящей статье я хотел бы сделать обзор популярных алгоритмов для решения задачи нахождения максимальной общей подпоследовательности или LCS (longest common sequense). Так как акцент сделан на обучении, а не на реальном использовании, в качестве языка для реализации выбран Python, что позволит сократить количество кода и сконцентрироваться на основных идеях.

Определения

Последовательность представляет собой упорядоченный набор элементов. Строка — это частный случай последовательности, дальнейшие примеры будут для простоты рассматривать именно строки, но без изменений их можно использовать и для произвольного текста или чего-нибудь последовательного еще.
Пусть имеется последовательность x, состоящая из элементов x₁x₂…x_m и последовательность y, состоящая из элементов y₁y₂…y_n. z — подпоследовательность x в том случае, если существует строго возрастающий набор индексов элементов x, из которых получается z.
Общей подпоследовательностью для x и y считаем такую последовательность z, которая является одновременно подпоследовательностью x и подпоследовательностью y.
Максимальная общая подпоследовательность — это общая подпоследовательность с максимальной длинной. Далее по тексту будем использовать сокращение LCS.
В качестве примера, пусть x=HABRAHABR, y=HARBOUR, в этом случае LCS(x, y)=HARBR. Можно уже переходить непосредственно к алгоритму вычисления LCS, но, хорошо бы понять, для чего нам может это может понадобиться.

Применение на практике

Наиболее частое применение — использование в программах для сравнения файлов, например GNU diff. Имея найденную для двух текстов LCS, составить список элементарных изменений для превращения x в y или обратно задача тривиальная. В качестве бонуса, на основе длины общей подпоследовательности можно задать метрику для определения схожести двух последовательностей. Все, теперь точно можно переходить к делу.

Первый подход или народное творчество

Для начала пара наблюдений:

1. Если для последовательностей x и y нами уже вычислена LCS(x, y)=z, то, LCS для последовательностей полученных из x и y путем добавления одинакового элемента, будет состоять из z и этого добавленного элемента
2. Если мы добавим к последовательностям x и y по одному разному элементу, то для полученных таким образом xa и yb, LCS должна быть большая из двух: LCS(x,yb) или LCS(xa,y)

Этих наблюдений уже достаточно, чтобы реализовать рекурсию:

def LCS_RECURSIVE(x, y):
    if len(x) == 0 or len(y) == 0:
        return []
    if x[-1] == y[-1]:
        return LCS_RECURSIVE(x[:-1], y[:-1]) + [x[-1]]
    else:
        left = LCS_RECURSIVE(x[:-1], y)
        right = LCS_RECURSIVE(x, y[:-1])
        return left if len(left) > len(right) else right

Теперь можно подумать, что с этой реализацией не так. В наихудшем случае, когда между x и y нет одинаковых элементов LCS_RECURSIVE вызовется 2^{len(x)+len(y)} раз, при уровне рекурсии len(x)+len(y). Даже если закрыть глаза на производительность, код:

from uuid import uuid4
x = [uuid4().hex for _ in xrange(500)]
y = [uuid4().hex for _ in xrange(500)]
print LCS_RECURSIVE(x,y)

без дополнительного вызова sys.setrecursionlimit удачно не выполнится. Не самая удачная реализация.

Динамическое программирование или 64 кб хватит всем

Рассматриваемый алгоритм также известен как алгоритм Нидлмана—Вунша (Needleman-Wunsch).
Весь подход сводится к поэтапному заполнению матрицы, где строки представляют собой элементы x, а колонки элементы y. При заполнении действуют два правила, вытекающие из уже сделанных наблюдений:
1. Если элемент x_i равен y_j то в ячейке (i,j) записывается значение ячейки (i-1,j-1) с добавлением единицы
2. Если элемент x_i не равен y_j то в ячейку (i,j) записывается максимум из значений(i-1,j) и (i,j-1).
Заполнение происходит в двойном цикле по i и j с увеличением значений на единицу, таким образом на каждой итерации нужные на этом шаге значения ячеек уже вычислены:

def fill_dyn_matrix(x, y):
    L = [[0]*(len(y)+1) for _ in xrange(len(x)+1)]
    for x_i,x_elem in enumerate(x):
        for y_i,y_elem in enumerate(y):
            if x_elem == y_elem:
                L[x_i][y_i] = L[x_i-1][y_i-1] + 1
            else:
                L[x_i][y_i] = max((L[x_i][y_i-1],L[x_i-1][y_i]))
    return L

Иллюстрация происходящего:

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

def LCS_DYN(x, y):
    L = fill_dyn_matrix(x, y)
    LCS = []
    x_i,y_i = len(x)-1,len(y)-1
    while x_i >= 0 and y_i >= 0:
        if x[x_i] == y[y_i]:
            LCS.append(x[x_i])
            x_i, y_i = x_i-1, y_i-1
        elif L[x_i-1][y_i] > L[x_i][y_i-1]:
            x_i -= 1
        else:
            y_i -= 1
    LCS.reverse()
    return LCS

Сложность алогоритма — O(len(x)*len(y)), такая же оценка по памяти. Таким образом, если я захочу построчно сравнить два файла из 100000 строк, то нужно будет хранить в памяти матрицу из 10¹⁰ элементов. Т.е. реальное использование грозит получением MemoryError почти на ровном месте. Перед тем как перейти к следующему алгоритму заметим, что при заполнении матрицы L на каждой итерации по элементам x нам нужна только строка, полученная на предыщем ходу. Т.е. если бы задача ограничивалась только нахождением длины LCS без необходимости вычисления самой LCS, то можно было бы снизить использование памяти до O(len(y)), сохраняя одновременно только две строки матрицы L.

Борьба за место или Алгоритм Хишберга (Hirschberg)

Идея в основе этого алгоритма проста: если разделить входную последовательность x=x₁x₂…x_m на две произвольные части по любому граничному индексу i на xb=x₁x₂…x_i и xe=x_i+1x_i+2…x_m, то найдется способ аналогичного разбиения y (найдется такой индекс j, разбивающий y на yb=y₁y₂…y_j и ye=y_j+1y_j+2…y_n), такой, что LCS(x,y) = LCS(xb,yb) + LCS(xe,ye).
Для того, чтобы найти это разбиение y предлагается:

1. Заполнить динамическую матрицу L с помощью fill_dyn_matrix для xs и y. L1 приравняем последней строке вычисленной матрицы L
2. Заполнить матрицу L для *xe и *y (обратные последовательности для xe и y, в терминах Python: list(reversed(x)), list(reversed(y))). Приравняем *L2 последнюю строку вычисленной матрицы L, L2 это риверс от *L2
3. Принять j равным индексу, для которого сумма L1[j]+L2[j] максимальна

Это можно представить, как заполнение матрицы L с двух противоположных концов:

Заметим, что раз есть необходимость только в последней строке матрицы L, то при вычислении хранить всю матрицу не нужно. Немного изменив реализацию fill_dyn_matrix получим:

def lcs_length(x, y):
    curr = [0]*(1 + len(y))
    for x_elem in x:
        prev = curr[:]
        for y_i, y_elem in enumerate(y):
            if x_elem == y_elem:
                curr[y_i + 1] = prev[y_i] + 1
            else:
                curr[y_i + 1] = max(curr[y_i], prev[y_i + 1])
    return curr

Теперь непосредственно, о самом алгоритме. Он представляет собой классический divide and conquer: пока в последовательности x есть элементы, мы делим x пополам, находим подходящее разбиение для y и возвращаем сумму рекурсивных вызовов для пар последовательностей (xb,yb) и (xe,ye). Заметим, что в тривиальном случае, если x состоит из одного элемента и встречается в y мы просто возвращаем последовательность из этого единственного элемента x.

def LCS_HIRSHBERG(x, y):
    x_len = len(x)
    if x_len == 0:
        return []
    elif x_len == 1:
        if x[0] in y:
            return [x[0]]
        else:
            return []
    else:
        i = x_len // 2
        xb, xe = x[:i], x[i:]
        L1 = lcs_length(xb, y)
        L2 = reversed(lcs_length(xe[::-1], y[::-1]))
        SUM = (l1 + l2 for l1,l2 in itertools.izip(L1, L2))
        _, j = max((sum_val, sum_i) for sum_i, sum_val in enumerate(SUM))
        yb, ye = y[:j], y[j:]
        return LCS_HIRSHBERG(xb, yb) + LCS_HIRSHBERG(xe, ye)

Теперь требования по памяти у нас O(len(y)), время, необходимое для вычисления, удвоилось, но асимптотически по-прежнему O(len(x)len(y)).

Алгоритм Ханта-Шуманского (Hunt-Szymanski Algorithm)

Первое что нам потребуется это создание хэш таблицы matchlist, которая будет содержать индексы элементов второй последовательностей. Время необходимое для этого O(len(y)). Следующий код на питоне делает это:

y_matchlist = {}
for index, elem in enumerate(y):
    indexes = y_matchlist.setdefault(elem, [])
    indexes.append(index)
    y_matchlist[elem] = indexes

Для последовательности «HARBOUR» хэш будет следующим {‘A’: [1], ‘B’: [3], ‘H’: [0], ‘O’: [4], ‘R’: [2, 6], ‘U’: [5]}.

Далее, перебирая элементы последовательности x, заполняем массив THRESH соответсвующими индексами из заготовленного matchlist, таким образом, что значением k-го элемента THRESH должен быть индекс y_index, при условии что THRESH[k-1] < y_index и y_index < THRESH[k]. Таким образом, в любой момент времени массив THRESH отсортирован и для нахождения подходящего k мы можем использовать бинарный поиск. При обновлении элемента THRESH мы также добавляем соответствующий индексу y_index элемент последовательности к нашей LCS. Внести ясность может следующий код:

x_length, y_length = len(x), len(y)
min_length = min(x_length, y_length)
THRESH  = list(itertools.repeat(y_length,  min_length+1))
LINK_s1 = list(itertools.repeat(None,      min_length+1))
LINK_s2 = list(itertools.repeat(None,      min_length+1))
THRESH[0], t = -1, 0

for x_index,x_elem in enumerate(x):
   y_indexes = y_matchlist.get(x_elem, [])
   for y_index in reversed(y_indexes):               
       k_start = bisect.bisect_left(THRESH, y_index, 1, t)
       k_end   = bisect.bisect_right(THRESH, y_index, k_start, t)
       for k in xrange(k_start, k_end+2):
           if THRESH[k-1] < y_index and y_index < THRESH[k]:
               THRESH[k] = y_index
               LINK_x[k] = (x_index, LINK_x[k-1])
           if k > t:
               t = k

После этого нам остается только собрать из LINK_x саму последовательность.
Время работы этого алгоритма равно O((r + n) log n), где n — длина большей последовательности, а r — количество совпадений, в худшем случае при r = n*n, мы получаем время работы хуже чем в предыдущем варианте решения. Требования по памяти остаются порядка O(r+n).

Итого

Алгоритмов решающих данную проблему довольно много, асимптотически, самый эффективный алгоритм (Masek and Paterson) имеет оценку по времени O(n*n/log n). Учитывая общую небольшую эффективность при вычислениях LCS, на практике перед работой алгоритма выполняются простейшие подготовки, вроде отбрасывания одинаковых элементов в начале и в конце последовательностей и поиск тривиальных отличий между последовательностями. Также, существуют некоторые оптимизации с использованием битовых операций, не влияющие на асимптотику времени работы.
скачать весь код с примерами