Алгоритм k различий Ландау-Вишкина [Landau, Vishkin, 1986b, 1989] основан на подходе, близком методу динамического программирования для вычисления расстояния между строками, который предложил Укконен [Ukkonen, 1983, 1985a]. Перед тем, как перейти к этому алгоритму, опишем метод динамического программирования и его адаптацию в стиле Укконена.

Вспомним, что элемент d_i,j, определяемый (9), сообщает расстояние между префиксами с длинами i и j строк, соответственно, x и y. Чтобы решить задачу k различий, матрицу расстояний надо преобразовать таким образом, чтобы d_i,j представлял минимальное расстояние между x(1, i) и любой подстрокой y, заканчивающейся символом y_j (Селлерс, 1980). Для этого достаточно изменить граничные условия (задаваемые (12)), на

так как минимальное расстояние между и любой подстрокой y равно 0.

Остальная часть матрицы вычисляется как раньше, с использованием цен редактирования расстояния Левенштейна, (10), и рекуррентного соотношения для d_i,j (11). По завершении, каждое значение, не превосходящее k, в конечной строке указывает позицию в тексте, в которой заканчивается строка, имеющая не больше k отличий от образца. Это иллюстрируется приведенным ниже примером, в котором показана матрица расстояния для случая x = ABCDE и y = ACEABPCQDEABCR. Из строки 5 этой матрицы можно видеть, что вхождения образца с точностью до 2 отличий, заканчиваются в позициях 3, 10, 13 и 1 Соответствующими подстроками являются ACE, ABPCQDE, ABC и ABCR.

Вспомним, что, как описано в алгоритме Укконена, при вычисления расстояний между строками диагонали матрицы можно пронумеровать целыми числами p[-m, n], таким образом, чтобы диагональ p состояла из элементов (i, j), у которых j - i = p. Пусть r_p,q представляет наибольшую строку i, у которой d_i,j = q и (i, j) лежит на диагонали p. Таким образом, q – это минимальное число различий между x(1, r_p,q) и любой подстрокой текста, заканчивающейся, и, кроме того, так как в противном случае требуемым условиям удовлетворяла бы строка r_p,q+1, которая больше строки r_p,q.

Вернемся к нашему примеру. Ниже представлена диагональ 9 матрицы расстояний. Мы видим, что значениями r_p,q в этом случае являются r_9,0 = 0, r_9,1 = 4 и r_9,2 = 5.

Значение m в строке r_p,q, для q < k, указывает, что в тексте имеется вхождение образца с точностью до k отличий, заканчивающееся в y_m+p. Таким образом, чтобы решить задачу k различий, достаточно вычислить значения r_p,q для q < k.

Чтобы объяснить, как можно вычислить значения r_p,q, рассмотрим вычисление входов d_i,j в матрице динамического программирования. Если d_i,j = q, этот вход находится на диагонали p, и i является наибольшей строкой, удовлетворяющей этим требованиям, тогда по определению r_p,q = i. Значение d_i,j будет выведено по одному из предшествующих ему соседних значений в матрице в соответствии с (11). При этом существуют следующие возможности:

1. d_i-1,j-1 = q и x_i = y_j
2. d_i-1,j-1 = q-1 и x_i =/=y_j
3. d_i,j-1 = q-1
4. d_i-1,j = q-1

Если последнее появление значения q-1 на диагонали уже найдено, последнее появление значения q на этой диагонали можно найти, проходя через последующие элементы диагонали, пока x_i = y_j.

Обратите внимание, что первое возможное вхождение образца с точностью до k отличий встречается в позиции текста j = m-k. Эта позиция в строке m лежит на диагонали -k, поэтому диагонали p < -k при вычислениях мы можем не рассматривать. Заметим также, что смежные элементы в матрице могут отличаться только на 0 или 1. Алгоритм, приведенный на рисунке 29, можно использовать для вычисления значений r_p,q. Начальное значение строки r получают из предыдущих значений r_p,q-1, r_p-1,q-1 и r_p+1,q-1, и последовательно увеличивают на единицу за один раз, пока не будет достигнуто правильное значение r_p,q.

То, что алгоритм правильно вычисляет значения r_p,q, можно показать следующим образом. Рассмотрим случай, когда q = 0 и вычисляются значения r_p,0 (p > 0). Из граничных условий, начальное значение r равно 0. Цикл while находит, что x(1, r_p,0) = y(p+1, p+r_p,0) и, устанавливая, таким образом, r_p,0 равным правильному значению. Для диагоналей p < 0 правильные значения для r_p,|p|-1 и r_p,|p|-2 будут установлены во время инициализации. Пусть q = t, и предположим, что все значения r_p,t-1 правильны. Рассмотрим вычисление r_p,t (p > - t). Переменная r устанавливается равной наибольшей строке на диагонали p, такой, что d_r,r+p может достичь величины t посредством возможностей 2, 3 или 4, описанных ранее. Цикл while после этого станет перебирать последующие позиции на диагонали, в которых символы шаблона и текста совпадают, пока не придет к правильном значению r_p,t.

- инициализация

rp,-1 = -1 for 0 < p < n

rp,|p|-1 = |p| - 1, rp,|p|-2 = |p| - 2 for -(k+1) < p < -1

rn+1,q = -1 for -1 * q * k

- вычисление значений rp,q

for q = 0 to k

   for p = -q to n

      r = max{rp,q-1 + 1, rp-1,q-1, rp+1,q-1 + 1}

      r = min{r, m}

      while (r < m) and (r + p < n) and (xr+1 = yr+1+p)

         r = r + 1

      rp,q = r

      if rp,q = m then

         имеется вхождение с k отличиями, заканчивающееся в yp+m

Алгоритм вычисляет значения r_p,t на n+k+1 диагоналях. Для каждой диагонали переменной строки r можно присвоить не больше m различных значений, что приводит ко времени вычислений O(mn). Таким образом, временные затраты те же, что и у прямого метода динамического программирования, описанного ранее. Однако, теперь мы покажем, как структуру этого метода можно использовать для разработки более эффективного алгоритма.

Улучшение эффективности достигается за счет изменения метода вычисления r_p,q. В описанном выше методе переменная r последовательно увеличивается на единичных шагах, пока не будет достигнуто правильное значение r_p,q. Однако, благодаря подходящим предварительным вычислениям, значение r_p,q можно находить за фиксированное время. Как это сделать, мы расскажем позже, однако сначала рассмотрим требуемые для этого предварительные вычисления.

Стадия предварительных вычислений включает построение так называемого суффиксного дерева, которое можно описать следующим образом. Рассмотрим строку s = s(1, q), где s_q = $ – маркер конца, не являющийся частью алфавита, над которым построена строка s(1, q-1). Каждый суффикс s(i, q), 1 < i < q, строки s определяет один из листьев суффиксного дерева для s. Все ребра в дереве направлены от корня, и количество выходящих из каждого узла либо равно нулю (у листьев), либо > 2. Если s(i, i+p) – самый длинный общий префикс двух суффиксов s(i, q) и s(j, q), то есть s(i, i+p) = s(j, j+p) и s_i+p+1 =/=s_j+p+1, то s(i, i+p) определяет внутренний узел дерева. Два узла, представляющих подстроку b и соответствующий префикс из нее a, соединены ребром только в том случае, если не существует узла, представляющего подстроку c, такого, что c является префиксом b и a является префиксом c. Прерывающий символ $ используется для того, чтобы разделять в дереве суффиксы s(i, q-1) и s(j, q-1), когда один из них является префиксом другого. Обратите внимание, что суффиксное дерево для данной строки уникально с точностью до изоморфизма графов. Структура данных суффиксных деревьев более подробно рассматривается в последующих разделах этой главы.

В качестве примера ниже приведено суффиксное дерево для строки EWEW$, в котором каждый узел помечен определяющей его подстрокой и номером в скобках.

корень --------------------------- (1) $

|  |

|  |

|   ------------- (2) W ---------- (4) W$

|                                |

|                                |

|                                  ----------- (5) WEW$

|

|

  --------------- (3)  EW -------- (6) EW$

                                 |

                                 |

                                   ------------ (7) EWEW$

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

start1 = 5 lenght1 = 1

start2 = 4 lenght2 = 1

start3 = 3 lenght3 = 2

start5 = 2 lenght5 = 4

start6 = 3 lenght6 = 3

start7 = 1 length7 = 5

На этапе предварительной обработки, суффиксное дерево строки y#x$, где # и $ – символы, не принадлежащие алфавиту, над которыми построены строки x и y, строится с помощью алгоритма Вейнера [Weiner, 1973] (который подробно описан Ченом и Сейферасом [Chen, Seiferas, 1985). Этот алгоритм требует линейных затрат памяти, и, для алфавита фиксированного размера, линейного времени. Для неограниченных алфавитов этот алфавит можно преобразовать, так что он будет выполняться за время O(n log), где – число различающихся символов образца. Стадия предварительной обработки требует время O(n) и O(n log m) для постоянного и неограниченного алфавитов, соответственно. Альтернативные способы предварительной обработки с использованием суффиксных деревьев рассматривались Галилом и Джанкарло (1988).

Теперь можно модифицировать алгоритм, приведенный на рисунке 29, чтобы воспользоваться преимуществами вышеуказанной предварительной обработки. Непосредственно перед циклом while для диагонали p, r было присвоено значение, такое, чтобы x(1, r) сопоставлялось с точностью до k различий с некоторой подстрокой текста, заканчивающейся y_r+p. Тогда функция цикла while находит максимальное значение, назовем его h, для которого x(r+1, r+h) = y(r+p+1, r+p+h). Это эквивалентно нахождению длины самого длинного общего префикса суффиксов x(r+1, m)$ и y(r+p+1,n)#x$ предварительно вычисленной конкатенированной строки. Символ # используется для предотвращения ситуаций, в которых может ошибочно рассматриваться префикс, состоящий из символов как y, так и x. Если lca(r,p) определяется как самый низкий общий предок (LCA) в суффиксном дереве с листьями, определенными вышеуказанными суффиксами, нужное значение h задается length_lca(r,p). Таким образом, алгоритм самого низкого общего предка Харела и Тарьяна [Harel, Tarjan, 1984] или Шибера и Вишкина [Schieber, Vishkin, 1988] может использоваться в алгоритме k различий в качестве средства, альтернативного выводу длин нужных сопоставляемых подстрок.

Суффиксное дерево имеет O(n) узлов. Для поддержки определения самого низкого общего предка за линейное время, алгоритмам LCA (самого низкого общего предка) требуется преобразование дерева, проводимое за линейное время. Как упоминалось ранее, значения r_p,q вычисляются на n+k+1 диагоналях. Более того, для каждой диагонали надо вычислить k+1 таких значений, что в общей сложности дает O(kn) запросов. Эти вычисления с временем O(kn) дают наибольший вклад в вычисление значений r_p,q. Таким образом, общее время прогона для этого алгоритма k различий составляет O(kn) для алфавитов фиксированного размера, и O(n * log(m) +kn) для неограниченных алфавитов.

Кроме того, Ландау и Вишкин [Landau,Vishkin, 1986b, 1989] разработали параллельную версию вышеприведенного алгоритма. Суффиксное дерево может быть вычислено за время O(log n) с использованием n процессоров (Ландау, Шибер и Вишкин [Landau, Schieber, Vishkin, 1987]; Апостолико et al. [Apostolico et al., 1988]). Параллельный алгоритм Шибера и Вишкина может использоваться для предварительной обработки дерева за время O(log n) с использованием n/log n процессоров, при условии, что отдельный процессор выполняет последовательный запрос самого низкого общего предка за время O(1). Таким образом, использование n+k+1 процессоров позволяет выполнять все требуемые запросы самого низкого общего предка за время O(k), – граница, которая также может быть достигнута моделированием алгоритма на n процессорах. Итак, общее требуемое время равно O(log n+k) для n процессоров. Кроме того, Ландау и Вишкин обсуждают, как геометрическая декомпозиция задачи может привести к адаптации алгоритма, по прежнему требующей n процессоров, но выполняющейся за время O(log m).