STR Строки и тексты Алгоритм
Префикс-функция и КМП
Также: Knuth–Morris–Pratt algorithm · KMP · prefix function · алгоритм Кнута — Морриса — Пратта · функция отказов · failure function
Поиск образца в тексте, при котором каждая буква текста читается один раз и указатель никогда не идёт назад. Всё, что нужно знать при несовпадении, заранее записано в префикс-функции образца: какой длины его начало повторяется в конце. Отсюда линейное время в худшем случае, поиск в потоке и простой способ найти период строки.
В каждой лаборатории молекулярной биологии есть пробирка с надписью «λ/HindIII». В ней ДНК бактериального вируса — фага лямбда, разрезанная ферментом HindIII. На геле эта смесь расходится на восемь полос известной длины: 23 130, 9416, 6557, 4361, 2322, 2027, 564 и 125 пар оснований, и по ним, как по линейке, меряют другие молекулы. Предскажем эту лестницу программой. Фермент режет ДНК везде, где встречается мотив AAGCTT, поэтому достаточно найти все вхождения шести букв в 48 502. Ниже их ищет алгоритм Кнута — Морриса — Пратта, по одной букве, слева направо. Выберите другой геном, другой фермент или впишите свой мотив.
Программа нашла шесть сайтов, кусков получилось семь, а у маркера полос восемь: вместо 6557 и 125 у нас один кусок 6682. Ошибки в поиске нет. В базе NCBI лежит геном дикого фага, а в пробирках продают ДНК лабораторного штамма — мутанта cI857ind1 Sam7. Его мутация ind1 — замена C на T в позиции 37 589, и последовательность AAGCTC в этом месте становится AAGCTT. Поставьте галочку «λ из пробирки»: кусок 6682 распадётся на 6557 и 125, и лестница совпадёт с маркером до последней полосы. Одна буква из 48 502 видна на геле.
Теперь посмотрите на счётчики. На лямбде с HindIII поиск сделал 60 083 сравнения букв, по 1,24 на букву, а наивный перебор всех позиций — 65 502, по 1,35. Разница небольшая, и это честный результат: в настоящей ДНК первая же буква образца обычно не совпадает, и наивному поиску почти нечего терять. Алгоритм КМП ценят за другое. Он ни разу не возвращается назад по тексту, поэтому годится для потока, который нельзя перемотать, и даже в самом неудачном тексте делает меньше 2n сравнений. Как выглядит неудачный текст, покажет вторая сцена.
Идея
Пусть образец совпал с текстом на первых k буквах, а на следующей разошёлся. Эти k букв текста мы уже знаем, не перечитывая: это начало образца. Значит, куда сдвинуть образец, можно решить, вообще не глядя в текст. Надо найти самое длинное начало образца, которым заканчиваются его первые k букв, и продолжить сравнение с его конца. Такие длины вычисляются заранее, один раз для каждого k, и указатель в тексте больше никогда не идёт назад.
Начало строки, которое одновременно её конец, но не вся строка, называют гранью. У слова «колокол» самая длинная грань — «кол»: им слово начинается и заканчивается. У «абракадабры» — «абра», у «GAATTC» граней нет. Префикс-функция строки p — это массив π, где π[i] — длина самой длинной грани начала p[0..i]:
к о л о к о л 0 0 0 0 1 2 3
Под последней «л» стоит 3: у всего слова грань «кол». Под второй «к» — 1: у «колок» грань «к». Таблица зависит только от образца, а поиск, как сейчас увидим, сводится к одному числу k и этой таблице. Префикс-функцию с нуля и шаг за шагом строит глава курса «Иголка в стоге»; здесь — доказательства, поведение на худших текстах, код для потока и задачи, где поиска нет вовсе.
Поиск без шагов назад
Идём по тексту слева направо и держим одно число k — сколько первых букв образца совпадает с концом прочитанного. Читаем букву c:
- если c равна p[k], совпадение продлилось: k растёт на единицу;
- если не равна и k > 0, образец сдвигается: k становится π[k − 1], и та же буква c сравнивается снова, уже с другой буквой образца;
- если не равна и k = 0, совпавших букв нет — переходим к следующей букве текста.
Когда k дорастает до длины образца m, вхождение найдено, и k становится π[m − 1]: следующее вхождение может начаться внутри этого. Вот оба поиска на одних и тех же строках, по одному сравнению за шаг. Слева у графика — номер буквы текста, на которую смотрит поиск, по горизонтали — сколько сравнений уже сделано.
На теломере — так называют концы хромосом, где у человека повторяется звено TTAGGG, — образец из пяти звеньев почти совпадает с текстом в каждой шестой позиции, и только последняя буква C его отвергает. Наивный поиск каждый раз проходит 29 совпадающих букв, спотыкается на тридцатой, сдвигается на одну позицию и перечитывает их снова: на графике это пила. КМП после несовпадения сдвигает образец ровно на одно звено: у совпавших 29 букв самая длинная грань — 23 буквы, на звено короче. И сравнение продолжается с той же буквы текста. Его линия — лестница, которая никогда не идёт вниз. На 300 буквах теломеры наивный поиск делает 1650 сравнений, КМП — 346.
Переключите на «AAAA…» — текст из одних A и образец из девятнадцати A и одной C. Здесь КМП приходится тяжелее всего: на каждой букве он сравнивает дважды — с C и, после отката, с A. Это 581 сравнение против 5620 у наивного. А на куске настоящей ДНК лямбды две линии почти сливаются: в обычном тексте наивный поиск хорош, и разница видна только на худших случаях.
Как считать префикс-функцию
Перебирать для каждого начала все возможные длины грани — порядка m³ действий. Но префикс-функция строится из самой себя, тем же движением, что и поиск. Пусть уже известно k = π[i − 1]. Если p[i] = p[k], грань продлевается: π[i] = k + 1. Если нет, нужна грань покороче, которая продлится. Какие ещё грани есть у строки, кроме самой длинной, подсказывает лемма из следующего раздела: это грань грани, потом грань грани грани и так далее. Их длины уже записаны в таблице: π[k − 1], π[π[k − 1] − 1], …
Нажмите на любую букву: под строкой появится цепочка граней начала, которое этой буквой заканчивается, а оранжевые черты отметят его наименьший период.
Посмотрите на «abaababaabaababaaba». Это строка Фибоначчи без двух последних букв: f₁ = b, f₂ = a, а дальше каждая следующая строка — склейка предыдущей и той, что перед ней, fk = fk−1fk−2. У неё необычно длинная для такой короткой строки цепочка граней: 11, 6, 3, 1. А у «aaaaaaaa» цепочка самая длинная из возможных — каждая следующая грань короче предыдущей ровно на одну букву.
Есть и короткий путь, без отдельной функции поиска. Склейте образец, разделитель, которого нет в строках, и текст: p + «#» + t. Префикс-функция склейки равна m ровно там, где в тексте кончается вхождение образца: длиннее m грань быть не может, потому что через «#» не перешагнёт. Поиск КМП — это то же самое вычисление, только без хранения склейки.
Почему это верно
Пусть β — самая длинная грань строки s. Тогда грани s, которые короче β, — это в точности грани самой β. Поэтому все грани s по убыванию длины — это π[n − 1], π[π[n − 1] − 1], … , пока не получится ноль.
Пусть γ — грань s и |γ| < |β|. И γ, и β — начала s, причём γ короче, значит, γ — начало β. И γ, и β — концы s, γ короче, значит, γ — конец β. Это и есть грань β. Обратно: грань β — начало β, а потому и начало s, и конец β, а потому и конец s. По той же причине самая длинная грань β — это следующая по длине грань s, и так далее вниз по цепочке.
Алгоритм построения вычисляет π правильно. После чтения буквы t[i] поиск хранит в k длину самого длинного начала образца, которое заканчивается в позиции i текста, и поэтому сообщает о каждом вхождении ровно один раз.
Посмотрим, откуда берутся грани начала s = p[0..i]. Если у s есть грань длины b > 0, то, отрезав у неё последнюю букву, получим грань длины b − 1 у строки p[0..i − 1], причём p[b − 1] = p[i]. Значит, каждая грань s — это продлённая на букву p[i] грань предыдущего начала (возможно, пустая). Обратное тоже верно: если у p[0..i − 1] есть грань длины c и p[c] = p[i], то у s есть грань длины c + 1.
Цикл перебирает грани p[0..i − 1] от самой длинной, π[i − 1], вниз по цепочке, и по лемме о гранях не пропускает ни одной. Первая же, которая продлевается буквой p[i], и даёт самую длинную грань s. Если не продлевается ни одна, грани у s нет, и π[i] = 0.
Поиск — то же вычисление для склейки p + «#» + t: значение k после буквы t[i] — это π склейки в этой позиции. Значение m означает вхождение, а следующее k = π[m − 1] — следующая по длине грань, как и требует лемма.
Корректность держится на одном наблюдении: когда сравнение сорвалось, единственное, что стоит помнить о прочитанном, — его самый длинный конец, совпадающий с началом образца. Всё остальное следует из таблицы. Поэтому буквы текста можно не хранить вовсе.
Код
Поиск оформлен как объект, которому текст скармливают кусками: между кусками он помнит только k и сколько букв прочитано. Так удобно читать большой файл по строкам, сеть по пакетам, ленту по блокам. Для сравнения рядом наивный поиск, и оба считают сравнения букв — те же, что показывает сцена выше.
Посмотрите на последние строки вывода. На теломере из 5000 звеньев образец «30 звеньев и TTAGGC» наивный поиск проверяет почти миллион раз, КМП — 34 970. А когда образец ровно 30 звеньев, КМП находит 4971 вхождение и делает ровно по одному сравнению на букву: после каждого найденного вхождения откат на одно звено оставляет совпавшими 174 буквы, и следующая буква сразу их продлевает.
Теперь на настоящем файле. Геном лямбды лежит в формате FASTA, по 70 букв в строке. Программа не склеивает строки в одну, а подаёт их по очереди девяти поискам сразу, по одному на фермент:
31 sites, split between two lines of the file: NdeI at 38357-38362
Сайт NdeI CATATG начинается в одной строке файла и заканчивается в следующей. Ни одна строка его целиком не содержит, а поиск всё равно нашёл: состояние k пережило перевод строки. Поиск, который на несовпадении откатывается назад по тексту, пришлось бы снабдить буфером. Именно эту трудность Джеймс Моррис и решал, когда придумал алгоритм.
Сколько это стоит
Построение префикс-функции делает меньше 2m сравнений букв, поиск в тексте длины n — меньше 2n. Время всего — O(n + m), память — O(m) на таблицу.
Каждое сравнение при поиске кончается одним из трёх исходов. Совпадение — переход к следующей букве текста. Несовпадение при k = 0 — тоже переход к следующей букве. Таких сравнений вместе не больше n. Несовпадение при k > 0 уменьшает k хотя бы на единицу, потому что грань короче строки. Но k растёт только при совпадении и только на единицу, то есть в сумме не больше чем на n, а ниже нуля не опускается. Уменьшений поэтому не больше, чем увеличений, и всех сравнений меньше 2n. Для построения π тот же счёт даёт меньше 2m.
Это амортизированная оценка: на одной букве текста откатов может быть много. Образец из m − 1 букв A и буквы C на тексте A…AB откатывается на букве B m − 1 раз подряд, зато эти откаты оплачены m − 1 буквами, которые им предшествовали. В среднем на букву выходит не больше двух сравнений, и дорогая буква не может повториться, пока k снова не накопится.
В статье 1977 года авторы улучшают таблицу: если после отката образец сравнивал бы ту же букву, что только что не совпала, откат делается сразу дальше. С такой таблицей откатов на одной букве текста не больше 1 + logφ m, где φ ≈ 1,618 — золотое сечение. Худший случай для неё дают как раз строки Фибоначчи. А Цви Галиль позже показал, как сделать поиск, который тратит на каждую букву ограниченное число действий, — это важно, когда букву надо обработать до прихода следующей.
Замер: текст — тысяча звеньев теломеры, образец — несколько звеньев и TTAGGC в конце, так что почти каждое шестое прикладывание доходит до последней буквы. Длина образца растёт, длина текста та же.
Наивный поиск дорожает вместе с образцом: от 96 букв до 384 образец удлинился вчетверо, и время выросло примерно вчетверо — это слагаемое n · m/6 от прикладываний к началу каждого звена. На коротких образцах его ещё заслоняет обычный проход по тексту, поэтому общий наклон меньше единицы. Линия КМП горизонтальная: длина образца на время поиска почти не влияет. Встроенный поиск Python написан на C и в этом замере быстрее всех; для длинных образцов он с версии 3.10 использует двусторонний алгоритм Крошмора и Перрена, который, как и КМП, линеен в худшем случае. Встроенный поиск Go для коротких образцов перебирает позиции машинными командами, сравнивающими сразу много байтов, а для длинных, если перебор часто ошибается, переходит на алгоритм Рабина — Карпа. В обычной работе пишите find, in и strings.Index; свой КМП нужен там, где встроенный поиск не подходит: в потоке, на последовательностях не из символов, с подсчётом откатов, как в задачах ниже.
Префикс-функция без поиска: период
Число q называют периодом строки s, если каждая буква повторяет ту, что стоит на q левее: s[i] = s[i + q] для всех i. У «колокол» период 4: «коло|кол». Длина строки не обязана делиться на период.
Наименьший период строки длины n равен n − π[n − 1].
Равенство s[i] = s[i + q] для всех i от 0 до n − q − 1 означает ровно то, что начало длины n − q совпадает с концом той же длины. Значит, q — период тогда и только тогда, когда у строки есть грань длины n − q. Самый маленький период соответствует самой длинной грани.
По лемме о гранях заодно получаются все периоды: n минус длины всех граней цепочки. Если n делится на наименьший период q, строка — это n/q копий своего начала длины q. Если не делится, составить её из одинаковых кусков короче всей строки нельзя; это следует из леммы о двух периодах Файна и Уилфа, которую доказывают и в статье Кнута, Морриса и Пратта.
Так префикс-функция отвечает на вопросы, где поиска нет вовсе: из какого звена собран повтор в ДНК, сжимается ли строка в «звено × число», какой самый короткий кусок надо повторять, чтобы получить данную строку. Попробуйте в сцене «TTAGGGTTAGGGTTA»: период 6, но 15 на 6 не делится — рид оборвался посреди звена. Это задача о теломере ниже.
Где это работает
Сам КМП в стандартных библиотеках встречается редко. В glibc, стандартной библиотеке C в Linux, длинные образцы ищет двусторонний алгоритм Крошмора и Перрена: он, как и КМП, линеен в худшем случае, но вместо таблицы на m чисел обходится несколькими переменными. GNU grep ищет фиксированные строки вариантами алгоритма Бойера — Мура, который на обычном тексте перепрыгивает через буквы, не читая их. Зато идея КМП — откат по заранее построенной таблице, не трогая текст, — работает везде, где образцов много или текст нельзя перемотать.
- Много образцов сразу. Автомат Ахо — Корасик (1975) — это префикс-функция, перенесённая с одной строки на бор из тысяч слов: у каждого узла ссылка отката на самый длинный суффикс, который тоже есть в боре. На нём была сделана первая версия fgrep в Unix, на нём же работают фильтры, которые ищут в трафике и файлах тысячи сигнатур.
- Автомат. Если алфавит маленький, как у ДНК, можно заранее посчитать для каждого k и каждой буквы, куда попадёт поиск с учётом всех откатов. Получится конечный автомат, который тратит на букву ровно один переход, без цикла. Так устроены движки регулярных выражений, которые строят автоматы; подробнее — в главе «Автоматы и регулярки». Собрать автомат КМП для своего слова можно в главе курса.
- Поток. КМП читает каждую букву один раз и помнит одно число, поэтому подходит для данных, которые проходят мимо и не сохраняются: пакеты в сети, показания датчика, лента.
- Структура строки. Периоды, повторы, сжатие, число различных подстрок, сравнение строки с её поворотами — префикс-функция и родственная ей Z-функция решают это за линейное время, и на олимпиадах по программированию это рабочие лошадки.
Редактор, автомат и задача о палиндромах
Летом 1969 года Джеймс Моррис писал текстовый редактор для машины CDC 6400. Буферизация файла была сложной, и Моррис хотел найти образец, ни разу не возвращаясь назад по тексту. Он придумал способ, опираясь на теорию конечных автоматов, но в его записи не было видно, что время линейно. В статье 1977 года соавторы рассказывают, что процедура Морриса оказалась слишком сложной для других разработчиков системы, и через несколько месяцев он обнаружил, что непрошеные «исправления» превратили её в развалины.
Независимо от него Дональд Кнут в начале 1970 года разбирал теорему Стивена Кука о двусторонних автоматах с магазинной памятью. Из неё следовало, что задача, которую Кнут не умел решить быстрее чем примерно за n² шагов, — распознать строку, начинающуюся с палиндрома чётной длины, — решается за линейное время. Кнут вручную проследил конструкцию Кука, выделил из неё механизм и обобщил его до программы, которая находит самое длинное начало одной строки, входящее в другую. По его словам, это был первый случай, когда теория автоматов научила его решать практическую задачу лучше, чем он умел раньше. Вон Пратт сделал время работы не зависящим от размера алфавита, а когда рассказал алгоритм Моррису, тот узнал в нём свой. Моррис и Пратт описали его в техническом отчёте Беркли в 1970 году, общая статья трёх авторов «Fast pattern matching in strings» вышла в журнале SIAM Journal on Computing в 1977-м.
Ещё раньше, в 1969 году, линейный поиск для текстов из нулей и единиц построил Юрий Матиясевич — в виде машины Тьюринга с двумерной памятью; об этом рассказывает глава курса.
Задачи
В первой задаче геном замкнут в кольцо, и вхождения могут перекрываться и проходить через место, где кольцо разрезали при записи в файл. Во второй ничего искать не надо: по куску повтора, который начинается и обрывается посреди звена, нужно найти само звено, и это задача о периоде.
Сайты на кольце Звено повтораДанные: Геномы фагов λ и φX174 и митохондрии человека: NCBI RefSeq NC_001416.1, NC_001422.1, NC_012920.1