Общее продолжение суффиксов: поиск с ошибками, джокеры и тандемные повторы

Если ты знаешь, как суффиксный массив с LCP и разреженной таблицей отвечают на LCE за O(1) и как на этом строят поиск с k ошибками и поиск повторов, листай до «Теперь сам».

Как думаешь, как искать с опечатками быстро

Ищем в геноме фрагмент из 300 букв, допуская две замены: так выглядят мутации или ошибки чтения. Наивно на каждом сдвиге сравниваем по букве, пока не наберётся третья ошибка. На случайном тексте это дёшево - ошибки набираются за пару букв. Но геном полон повторов: на сдвигах, где фрагмент почти совпадает, сравнение идёт до конца - m шагов.

Прежде чем читать дальше, ответь себе: на совпадающем куске нам ведь не важно, ЧТО там написано, - нам важно, ГДЕ кончается совпадение. Есть ли структура, которая сразу говорит, сколько символов совпадает начиная с двух заданных мест?

На что это похоже

Сверка двух длинных договоров. Юрист не перечитывает каждую строку: у него есть программа, которая мгновенно показывает, до какого места тексты совпадают, начиная с выбранных строк. Он прыгает от расхождения к расхождению: ткнул - «совпадает 412 слов», пропустил, отметил правку, ткнул снова. Число прыжков - число правок, а не длина договора.

Механизм: LCE через суффиксный массив, и три применения

LCE(i, j) = длина общего начала суффиксов s[i:] и s[j:]
    ранги r_i < r_j в суффиксном массиве →
    LCE = min LCP[r_i + 1 .. r_j]       (соседи по порядку
                                          совпадают не меньше)
    минимум на отрезке - разреженная таблица: O(1) на запрос
    предподсчёт: суффиксный массив, Касаи, таблица - O(n log n)

текст и образец склеивают через разделитель: s = образец # текст
и LCE между позицией в образце и позицией в тексте - тот же запрос

k ОШИБОК: на сдвиге t
    p = 0; повторять: d = LCE(p, t + p);  p += d
           если p < m: ошибка, p += 1;  ошибок > k → сдвиг не подходит
    не больше k + 1 запроса на сдвиг: O(n · k) вместо O(n · m)

ДЖОКЕРЫ «?» в образце: тот же прыжок, но на джокере ошибка
    не считается; k джокеров - k + 1 запрос

ТАНДЕМНЫЙ ПОВТОР ww с |w| = p, начало в i:  LCE(i, i + p) ≥ p

Почему минимум LCP: суффиксы в массиве отсортированы, и общее начало двух далёких суффиксов не может быть длиннее общего начала любой пары соседей между ними - иначе порядок нарушился бы. А оно и не короче: если все соседи между ними совпадают в первых h символах, совпадают и крайние.

Правило

LCE двух позиций - минимум массива LCP между их рангами в суффиксном массиве, и разреженная таблица даёт его за O(1). Поиск с не более чем k несовпадениями прыгает через совпадения: k + 1 запрос на сдвиг, O(n · k) вместо O(n · m); джокеры обрабатываются так же. Позиция i начинает повтор ww с |w| = p, когда LCE(i, i + p) ≥ p.

Разбор: LCE в слове banana

суффиксы по порядку и LCP (с соседом сверху):
    ранг 0: a       (5)   LCP 0
    ранг 1: ana     (3)   LCP 1
    ранг 2: anana   (1)   LCP 3
    ранг 3: banana  (0)   LCP 0
    ранг 4: na      (4)   LCP 0
    ранг 5: nana    (2)   LCP 2

LCE(1, 3): позиции 1 (anana, ранг 2) и 3 (ana, ранг 1)
    min LCP[2..2] = 3 → «ana» общее
LCE(1, 5): anana (ранг 2) и a (ранг 0): min LCP[1..2] = min(1, 3) = 1
повтор с p = 2 в позиции 1: LCE(1, 3) = 3 ≥ 2 → «anan» = (an)(an)

Повтор «anan» найден одним запросом. Так же находится «nana» в позиции 2: LCE(2, 4) = 2 ≥ 2.

Разбор: геном с повторами

Программа строит текст из повторяющегося блока длиной 50 с мутациями, берёт образец из 300 букв с двумя заменами и ищет его с не более чем двумя ошибками - прыжками LCE и наивно. Затем считает, в скольких местах начинаются повторы ww разной длины.

import random


def postroit(s):                      # суффиксный массив, LCP, таблица
    n = len(s)
    sa = sorted(range(n), key=lambda i: s[i:])
    rang = [0] * n
    for r, i in enumerate(sa):
        rang[i] = r
    lcp, h = [0] * n, 0               # Касаи: lcp[r] - с соседом сверху
    for i in range(n):
        if rang[i] > 0:
            j = sa[rang[i] - 1]
            while i + h < n and j + h < n and s[i + h] == s[j + h]:
                h += 1
            lcp[rang[i]] = h
            h = max(h - 1, 0)
        else:
            h = 0
    tabl = [lcp]                      # разреженная таблица минимумов
    k = 1
    while 2 * k <= n:
        p = tabl[-1]
        tabl.append([min(p[i], p[i + k]) for i in range(n - 2 * k + 1)])
        k *= 2
    return rang, tabl


def lce(rang, tabl, i, j, n):         # общее начало суффиксов i и j
    if i == j:
        return n - i
    a, b = sorted((rang[i], rang[j]))
    a += 1                            # минимум lcp на (a, b]
    k = (b - a + 1).bit_length() - 1
    return min(tabl[k][a], tabl[k][b - 2 ** k + 1])


rng = random.Random(31)
blok = "".join(rng.choice("acgt") for _ in range(50))
tekst = list(blok * 60)               # повторы, как в ДНК, с мутациями
for _ in range(60):
    tekst[rng.randrange(len(tekst))] = rng.choice("acgt")
tekst = "".join(tekst)
obrazec = list(tekst[1000:1300])      # образец из 300 букв, две ошибки
obrazec[70], obrazec[200] = "g", "t"
obrazec = "".join(obrazec)
s = obrazec + "#" + tekst             # склейка: образец и текст
rang, tabl = postroit(s)
m, n = len(obrazec), len(s)

for _ in range(2000):                 # LCE против наивного сравнения
    i, j = rng.randrange(n), rng.randrange(n)
    h = 0
    while i + h < n and j + h < n and s[i + h] == s[j + h]:
        h += 1
    assert lce(rang, tabl, i, j, n) == h

k, nashli, zaprosov = 2, [], 0       # не больше k несовпадений
for start in range(len(tekst) - m + 1):
    p, t, oshibok = 0, m + 1 + start, 0
    while p < m and oshibok <= k:
        zaprosov += 1
        d = lce(rang, tabl, p, t, n)  # прыжок через все совпадения
        p, t = p + d, t + d
        if p < m:
            oshibok, p, t = oshibok + 1, p + 1, t + 1
    if oshibok <= k:
        nashli.append(start)

naivno, sravnenij = [], 0            # наивно: до k + 1 ошибки
for st in range(len(tekst) - m + 1):
    oshibok = 0
    for a, b in zip(obrazec, tekst[st:st + m]):
        sravnenij += 1
        oshibok += a != b
        if oshibok > k:
            break
    if oshibok <= k:
        naivno.append(st)
assert nashli == naivno
print("вхождений с ≤ 2 ошибками:", len(nashli))
print("запросов LCE", zaprosov, " сравнений у наивного", sravnenij)

kvadraty = {}                         # повторы ww: lce(i, i + p) ≥ p
for p in (1, 2, 3, 25, 50, 100):
    kvadraty[p] = sum(1 for i in range(m + 1, n - 2 * p + 1)
                      if lce(rang, tabl, i, i + p, n) >= p)
for p, c in kvadraty.items():
    print(f"повторов ww с |w| = {p}: начинается в {c} местах")

вхождений с ≤ 2 ошибками: 1
запросов LCE 8103  сравнений у наивного 15246
повторов ww с |w| = 1: начинается в 893 местах
повторов ww с |w| = 2: начинается в 247 местах
повторов ww с |w| = 3: начинается в 0 местах
повторов ww с |w| = 25: начинается в 0 местах
повторов ww с |w| = 50: начинается в 630 местах
повторов ww с |w| = 100: начинается в 105 местах

Ответы совпали, а запросов LCE почти вдвое меньше, чем сравнений у наивного поиска, - и это на тексте, где почти все сдвиги отбрасываются за пару букв. Разрыв набирается на «почти совпадающих» сдвигах: там наивный проходит все 300 букв, а LCE - три прыжка. Повторы ww с |w| = 50 начинаются в 630 местах - это и есть структура текста: блок длиной 50, повторённый с мутациями.

Что ломается без этого

Первая ошибка - забыть разделитель при склейке образца и текста: общее начало «перетечёт» из образца в текст, и LCE окажется длиннее реального. Разделитель должен быть символом, которого нет нигде. Вторая - искать минимум LCP на [r_i, r_j], включив LCP[r_i]: это общее начало с соседом СВЕРХУ от меньшего ранга, а не с нужным суффиксом. Третья - строить суффиксный массив сортировкой суффиксов как строк на длинном тексте: это O(n² log n) по времени и копирование всех суффиксов; на миллионах символов нужны линейные алгоритмы построения. Четвёртая - искать с ошибками, где важны вставки и удаления: LCE считает только замены; для вставок нужна динамика расстояния редактирования.

Зачем это в работе

Биоинформатика: поиск прочтений с ошибками в геноме, поиск повторов (тандемные повторы связаны с болезнями и используются в анализе ДНК). Поиск по шаблонам с подстановками в логах и коде. Сравнение версий и поиск дубликатов в больших текстах. Общий урок: после построения индекса над текстом многие задачи «сравнить куски» решаются запросом за константу, и цена переходит с длины совпадения на число мест, где оно кончается.

Теперь сам

Строка abcabcabd. Есть ли тандемный повтор с p = 3, начинающийся в позиции 0? В позиции 3? Посчитай LCE напрямую.