Поиск наибольшей общей последовательности
Вместо вступления
Задача нахождения общей подпоследовательности двух последовательностей широко ис-пользуется во многих областях, например: проверки правописания, сравнения файлов, изучение мутации генов. И если для первых двух областей достаточно сказать одинаковые последова-тельности или нет (или количество общих элементов) – для этого подойдет алгоритм 2.а, то в задачах, связанных с генетикой, необходимо найти общую подпоследовательность максималь-ной длинны. Для этой цели мог бы подойти любой алгоритм, O(N2) времени и использующий двумерный массив, но для его хранения памяти может быть недостаточно, т.к. длина цепочки ДНК может достигать нескольких сотен тысяч, а иногда и миллионов символов (обозначающих кислоты).
Для подобных задач, как нельзя лучше, подойдет алгоритм 2.б.
В статье кроме алгоритмов 2.а и 2.б, будут рассмотрены алгоритмы, им предшествовав-шие, работающие над той же проблемой, но с меньшей производительностью.
Итак, начнем!
Алгоритм 1. Малого объема памяти, но обратно пропорционального ему времени
Данный алгоритм достаточно прост и может быть применен на небольших похожих по-следовательностях, в приложениях не часто использующих данную операцию. Суть метода за-ключается в том, что программа использует символы двух исходных последовательностей для построения графа (возможно циклического), одновременно ищет путь до конечной точки (на рис.1 – большая красная точка) с минимальным количеством несовпадений (на рис.1 несовпаде-ния – не вертикальные ребра).
Временная оценка метода: Временная сложность данного алгоритма –((2*(N+M))!/((N+M)!(N+M)!)) времени. Данный алгоритм нахождения максимальной общей подпоследовательности похож на переборный вариант решения задачи о минимальном пути че-репашки, поэтому за подробностями оценки можно обращаться к [2].
Суть метода: Метод основан на обходе графа в глубину с возвратом. На каждом шаге про-граммы прибавляется только одно на совпадающее ребро и следующее за ним совпадающие ребра (если таковые существуют). По достижении конца графа, возвращаемся на начальное зна-чение предыдущего шага, и если первый раз на данном шаге поиска в графе путь «уходил» вправо, то при возврате путь «уходит» влево, создавая дальнейший путь. (Построение графа по-казано на рис.1).
Усовершенствование: Для некоторого ускорения в работе программы можно учитывать еще количество шагов по первому и дальнейшему достижению конечной точки. И если в теку-щем пути количество невертикальных ребер превосходит количество невертикальных ребер лучшего пути, то построение данной ветки в графе закончено.
Часть алгоритма с некоторым усовершенствованием представлена ниже:
Procedure Solve(x,y:integer{length(a),length(b)}); // A,B – последовательности, представлен-ные в виде массива
begin
delete(Res,((x+y-step)div 2)+1,length(Res));
// проходим по вертикальным ребрам
while (A[x+1]=B[y+1])and(x
inc(x);
inc(y);
Res:=Res+A[x];
end;
//вызов рекурсии с новыми значениями индекса сравниваемых символов
if ((x
if x
if y
end;
//улучшение и хранение оптимального пути
end else
begin
if (x>=length(A))and(y>=Length(B)) then
begin
if length(Res)>Length(BestRes) then
begin
BestStep:=Step;
BestRes:=Res;
end;
end;
end;
dec(step);
end;
Алгоритм 2. O(ND) временной алгоритм
Временных O(ND) алгоритмов существует несколько, мы не будем затрагивать все из них, только укажем, что данный алгоритм 2.а (алгоритм Myers’а [1]) несколько напоминает алгоритм Нудельмана-Вунша [2] и по времени работы и по объему занимаемой памяти (если требуется вывод всех путей), но в некоторых случаях он может быть гораздо производительней – при не-больших различиях в последовательности (это существенной заметно, если количество разли-чий не превосходит 25% символов). В случае если требуется только посчитать количество эле-ментов в подпоследовательности, то данный алгоритм может быть очень хорош, т.к. в этом слу-чае он занимает линейный объем памяти (особенно эффективен при больших длинах примерно одинаковых последовательностей).
Алгоритм 2.а O(ND) временной алгоритм
Предварительные разъяснения. Редактируемый граф.
Пусть A=a1a2…an и B=b1b2…by — последовательность длины N и M соответственно. Ре-дактируемый граф A и B обладает вершинами в каждой точке сети (x,y), где x[0,N], y[0,M]. Редактируемый граф обладает вертикальными, горизонтальными и диагональными ребрами. Вертикальные ребра соединены с соседними вершинами справа, т.е. (x-1,y)(x,y), где x[1,N], y[0,M]. Горизонтальные ребра соединены с соседними вершинами сверху, т.е. (x,y-1)(x,y), где x[0,N], y[1,M]. Если ax=by, тогда образуется диагональное ребро, соединяющее диаго-нальных соседей так, что (x-1,y-1)(x,y). Точки (x,y), для которых ax=by назовем равными точ-ками. На рис.2 представлен граф для двух последовательностей (последовательности отличны от последовательностей алгоритма 1).
Путь длинной L это последовательность из L равных точек, (x1,y1)(x2,y2)…(xL,yL), таких, что xi
Суть: Пусть D-путь начинается в (0,0), который содержит D недиагональных ребер. 0 путь доложен содержать только диагональные ребра. На следующем шаге мы добавляем к имеюще-муся пути еще одно либо вертикальное, либо диагональное ребро (вызываем либо (x+1,y), либо (x,y+1)) и следующие за ним диагональные ребра (если таковые имеются). Таким образом, D путь должен состоять из D-1 следующего недиагонального ребра и возможно пустой последова-тельностью диагональных ребер, названных змейкой. Прохождение по графу будет продолжать-ся до тех пор, пока (x,y) не достигнет правого нижнего угла (если смотреть на рисунок так, как задумывал автор).
Необходимо упомянуть два важных момента для создания алгоритма:
1. за один шаг невозможно пройти более 1 недиагонального ребра
2. при добавлении змейки к D-пути, необходимо брать наибольшее число диагональных ребер, иначе полученный путь может оказаться не наибольшей общей подпоследовательностью
3. из 1 следует, что любой D-путь содержит ровно (D-1) не диагональное ребро (т.к. всего D шагов на каждом из которых, кроме первого, содержалось 1 недиагональное ребро)
Представленный ниже процедура предлагают нахождение длины максимальной общей по-следовательности.
Procedure Solve;
var x,y,k,D,Dout:integer;
begin
Res[1]:=0; //вспомогательное значение для 0 пути
//внешний цикл For – расширяет область диагоналей на каждом шаге
For D:=0 to length(A)+Length(B)do
begin
//внутренний цикл For – начинает проход с самой левой диагонали и
//заканчивает на самой правой
For k:=-d to d do
//проход по диагоналям, находящимся через 1 диагональ с самой левой
//об этом позже в леммах 1 и 2
If (abs(k) mod 2)=(d mod 2)then begin
//выбор диагонали (D-1)-пути (пути-предка)
If (k=-d)or(k<>D)and(Res[k-1]
else x:=Res[k-1]+1;
y:=x-k;
//создание змейки максимальной длины
While x
inc(y); end;
Res[k]:=x;
//проверка условия выхода из процедуры
If (x>=length(A))and(y>=length(B))then
begin
Res[k] – длина максимальной последовательности};
exit
end
end;
end;
На практике алгоритм ищет D-пути, где D<= максимуму и если такой путь достигает (N,M) тогда любая редактируемая последовательность для A и B должна по строке exit быть больше чем максимальная длина последовательностей (далее просто максимум). Полагая, что максимум эта сумма двух длин последовательностей, алгоритм гарантированно найдет макси-мальную длину общей последовательности. Рис.3 ниже показывает найденный D-путь, когда алгоритм применен к примеру на рис.1. Заметим, что несуществующая точка (0,-1) была введена для нахождения змейки 0-пути. Также заметим, что D-пути не увеличиваются, достигая правой и нижней границ в редактируемом графе при работе алгоритма. Эта пограничная ситуация об-работана тем условием, что после граничных ситуаций в одной из последовательностей нет символов.
Что бы не быть голословным докажем 1 и 3 утверждения (они также упоминаются в опи-сании процедуры Solve):
Лемма 1: D путь должен заканчиваться на диагонали k{-D,-D+2,…D-2,D}.
Доказательство: 0 путь состоит исключительно из диагональных ребер и начинается на диагонали 0. Отсюда он должен заканчиваться на диагонали 0. Предположим, что D путь дол-жен заканчиваться на диагонали k{-D,-D+2,…D-2,D}. Каждый (D+1) путь состоит из преды-дущего D пути, заканчивающегося, как сказано, на диагонали k, недиагонального ребра закан-чивающегося на диагонали k-1 (k+1) и змейки, которая также должна заканчиваться на диагона-ли k-1 (k+1). Из этого следует, что каждый (D+1) путь должен заканчиваться на {(-D)1,(-D+2) 1,… (D-2) 1,D1}={-D-1,-D+1,-D+3,…D-1,D+1}. Таким образом, можно продолжить рассуж-дения для следующих D+1+l-путей и для предыдущих D-l-путей, где l-произвольное целое чис-ло. Лемма подразумевает, что D путь заканчивается исключительно на нечетной диагонали, ко-гда D нечетное и на четной, когда D четная.
D путь является самым длинным на диагонали k, тогда и только тогда, когда один из D пу-тей заканчивается на диагонали k, чья конечная точка содержит наибольшее возможное число рядов (столбцов) всех таких путей. Другими словами из всех D путей, заканчивающихся на диа-гонали k, он заканчивается на самом длинном пути из всех начинавшихся в (0,0). Следующая лемма дает характеристику длиннейшего D пути, и реализует жадный принцип: длиннейший D путь получается путем увеличения длиннейшего (D-1) пути.
Лемма 2: Самый длинный 0 путь заканчивается в (х,х), где х – минимум (z-1II axbz или z>M или z>N). Самый длинный D путь на диагонали k может быть заменен без потерь на самый длинный D-1 путь на диагонали k-1, следующее горизонтальное ребро, и следующая наиболее длинная змейка или он может быть заменен на самый длинный D-1 путь на диагонали k+1, сле-дующее вертикальное ребро, и следующая наиболее длинная змейка.
Доказательство: Основа 0 путей проста. Как замечалось ранее, D путь содержит D-1 путь, недиагонального ребра и змейки. Если D путь заканчивается на диагонали k, из этого следует, что D-1 путь должен заканчиваться на диагонали k1 в зависимости от того, горизонтальное или вертикальное ребро было выбрано перед финальной змейкой. Финальная змейка должна быть максимальной, т.к. D путь не будет максимальным, если змейку можно будет увеличить. Пред-положим D-1 путь не самый дальний на этой диагонали. Но тогда самый дальний D-1 путь мо-жет быть связан с последней змейкой D пути данным недиагональным шагом. Таким образом, D путь по желанию может быть всегда разложен на D-1-путь, недиагональное ребро и змейку (ес-ли хранить предыдущее значение).
Время выполнения процедуры: Данная процедура обладает временем O((M+N)max(N,M)). Складывающимся из:
• внешний цикл повторяется в худшем случае M+N раз
• внутренний цикл For k:=-D to D повторяется 2D+1 раз
• обращение ко всем внутренним операциям во внутреннем цикле For k:=-D to D происхо-дит за счет условия (abs(k) mod 2)=(d mod 2), которое определяет шаг ((D div 2)+(D mod 2)) раз
• таким образом внутренние операции, кроме цикла While повторяются D/2 раз (округле-ние)
• если цикл While повториться хотя бы один раз, то время выполнения процедуры будет улучшено, поэтому его повторение мы будем учитывать только для худшего случая, т.о. он будет пропорционален O(1)
• исходя из всех вышеперечисленных пунктов общее время складывается только из двух циклов For и времени затрачиваемом на внутренние операции, и в худшем случае на данную процедуру уйдет O((M+N)max(N,M)) времени.
Замечание: Как было сказано выше, для вывода обратного пути итогового массива D не-достаточно, необходимо знать еще и предыдущее значение. Поэтому для вывода обратного D-пути следует хранить все исходные массивы, что иногда может и не удовлетворять условиям задачи (на больших последовательностях), но в некоторых случаях он может быть полезен.
Алгоритм 2.б Улучшенный алгоритм 2.а с линейным объемом памяти
Усовершенствованный алгоритм Mayers’а [1].
Как было сказано выше, для вывода пути необходимо запоминать массив D на каждом по-вторении внешнего цикла For. Таким образом, для вывода общей последовательности необхо-дим квадратичный объем памяти. Поэтому попробуем уменьшить объем до линейного. Для это-го немного изменим структуру алгоритма. Алгоритмом, представленным выше, определяет ко-личество не диагональных ребер в максимальной общей последовательности.
Для дальнейшей сборки алгоритма используем тип алгоритма “разделяй и властвуй”. Раз-делим максимальную последовательность пополам, пусть мы остановились на точке (x,y), затем найдем максимальную последовательность от (0,0) до (x,y) и от (x’,y’) до конечной точки. Сле-дующим шагом повторим вышеописанные действия. Так повторяем до тех пор, пока в общей последовательности не останется недиагональных ребер 1. Для доказательства возможности использования данного алгоритма воспользуемся леммой 3 и ее доказательством.
Лемма 3: D-путь из (0,0) до (N,M) существует тогда и только тогда, когда существует [D/2]-путь из (0,0) в (x,y) и [D/2]-путь из некоторой точки (u,v) до (N,M), такой, что:
{выполнимость} (u + v) [D/2] и (x + y) (N+M-[D/2]) и
{частичное совпадение} x – y = u – v и x u
Более того, оба [D/2]-пути составляют D-путь из (0,0) до (N,M)
Доказательство: Предположим, что существует D-путь из (0,0) до (N,M). Он может быть разделен от начала (0,0), до точки (x,y),что является серединой змейки [D/2]-пути из (0,0) до (x,y) и [D/2]-пути от (u,v) до (N,M), где (u,v)=(x,y). Путь из (0,0) до (u,v) может содержать мак-симум u + v не диагональных ребер и [D/2]-путь до (u,v) предполагает, что u + v [D/2]. Путь их (x,y) до (N,M) может содержать (N+M)-(x+y) не диагональных ребер и [D/2]-путь подразумева-ет, что x+yN+M-[D/2]. И, наконец, u-v=x-y и ux, т.к. (x, y)=(u,v).
Обратно, предположим [D/2] и [D/2]-пути существуют. Но ux подразумевает, что сущест-вует k-путь из (0,0) в (u,v), где k[D/2]. По лемме 1, =[D/2]-k кратно 2, т.к. оба пути и k и [D/2] заканчиваются на одной диагонали. Более того, k-путь содержит (u + v - k)/2/2 диагоналей, т.к. u + v [D/2]. Заменив каждую /2 диагонали в k-пути с парой вертикальных и горизонталь-ных ребер, [D/2] путь из (0,0) до (u,v) найдется. Но тогда будет существовать D-путь из (0,0) до (N,M) состоящий из этого [D/2]-путь до (u,v) и и данного [D/2]-путь из (u,v) до (N,M). Замете, что [D/2]-путь это часть D-пути. И симметричный [D/2]-путь это также часть D-пути из (0,0) до (N,M).
Ниже приведены процедуры для алгоритма, использующего вышеописанный метод.
//процедура нахождения серединного пути для двух
//подпоследовательностей
function MiddlePoint(var x,y,u,v:integer):Integer;
var D,k,l,delta:integer;
way,back:MyArray;
begin
delta:=length(AA)-Length(BB);
FillChar(Way,SizeOf(Way),0);
FillChar(Back,SizeOf(Back),10);
Way[1]:=0;
Back[delta+1]:=Maxim(Length(AA),Length(BB))+1;
For D:=0 to ((length(AA)+Length(BB))div 2)+1 do
begin
//прямой путь
For k:=-d to d do
If (abs(k) mod 2)=(d mod 2)then
begin
If (k=-d)or(k<>D)and(Way[k-1]
else x:=Way[k-1]+1;
y:=x-k;
While (x
inc(y); end;
Way[k]:=x;
//выход из функции нахождения срединной точки
If (x>=length(AA))and(y>=length(BB))or (Way[k]>=Back[k])
then begin
u:=Back[k];
v:=Back[k]-k;
MiddlePoint:=2*D-1;
exit;end;
end;
//обратный путь
For k:=-d to d do
If (abs(k) mod 2)=(d mod 2)then
begin
// переопределение нумерации диагоналей в обратном пути
// для боле легкого понимания
l:=k+delta;
If (k=-d)or(k<>D)and(Back[l-1]>Back[l+1])
then u:=Back[l+1]-1
else u:=Back[l-1];
v:=u-l;
While (u>0)and(v>0)and(AA[u]=BB[v])do
begin dec(u);
dec(v); end;
Back[l]:=u;
//выход из функции нахождения срединной точки
If (u<=0)and(v>=0)or (Way[l]>=Back[l])
then begin
x:=Back[l];
y:=Back[l]-l;
MiddlePoint:=2*D;
exit;end;
end;
end;
end;
// функция нахождения подпоследовательности из имеющейся
// последовательности, координат начала и конца подпоследовательности
function out(Script:string;x,u:integer):string;
var i:integer;
st:string;
begin
st:='';
for i:=x to u do
st:=st+Script[i];
out:=st;
end;
// основное тело процедуры решения
Procedure Solve(AA,BB:string);
var x,y,u,v:integer;
begin
if (length(AA)>0)and(Length(BB)>0) then
begin
// если в пути до срединной точки более одного недиагонального ребра,
// то делим путь дальше
if MiddlePoint(x,y,u,v)>1 then
begin
Solve(out(AA,1,u),out(BB,1,v));
if (x>u)and(y>v)then write(Out(AA,u+1,x));
Solve(out(AA,x+1,length(AA)),out(BB,y+1,length(BB)));
end else
// если путь содержит менее двух различий, то выводим
// последовательность минимальной длинны
if length(AA)<=length(BB) then
write(AA) else write(BB)
end;
end;
По примерным подсчетам времени на работу данной процедуры требуется в два раза больше, чем для описанной выше процедуры с циклами For и While, также время тратиться на проверку всех условий внутри рекурсии. Время работы программы, как было сказано выше не намного хуже лучших вариантов реализации. Теперь подсчитаем объем памяти. На хранение всех последовательностей в стеке рекурсии требуется объем равный 2 объема, необходимых для хранения исходных последовательностей (т.к. алгоритм делит последовательность в среднем случае пополам, в худшем сначала сохраняется большая часть, затем меньшая). Так же для хра-нения координат точек в стеке требуется два массива на D элементов. В итоге необходимый объем памяти составляет: O((M+N)lgS+(M+N)lg(M+N)), где
(M+N)lgS – объем, необходимый для хранения последовательностей, (M+N) – длинны по-следовательностей, lgS – байты необходимые для хранения символа
(M+N)lg(M+N) – объем, необходимый для хранения координат точек, (M+N) – общее ко-личество точек в худшем случае, lg(M+N) – байты необходимые для хранения координат симво-лов.
На рис.3 представлены две последовательности, разбитые алгоритмом.
Замечания:
Приведенные выше алгоритмы достаточно широко используются для решения разнооб-разных задач. Но ни один из приведенных алгоритмов не является универсальным для различ-ных последовательностей (кстати, элементами последовательности могут быть не только симво-лы, но и слова, а также в роли последовательности может выступать массив чисел или состоя-ний). Поэтому при решении задач следует комбинировать алгоритмы в зависимости от резуль-тата и входных данных задачи.
В случае если набор входных данных неизвестен или вы уверены, что около 30 – 50 % символов совпадут, то лучше всего воспользоваться квадратичным алгоритмом. В этом случае он будет оптимальным алгоритмом по времени.
Если вам известно, что последовательности слишком большие, для хранения их в двумер-ном массиве, то для решения задачи придется пожертвовать временем (существуют измененный алгоритм Нудельмана-Вунша, работающий за большее время, но с линейным объемом – обра-щаться к автору статьи).
Алгоритм 2.б хорошо подходит для последовательностей любой длинны, но при условии, что две последовательности достаточно сильно совпадают (он использовался для нахождения мутации генов ДНК длинной более 200 000 элементов, но различия между двумя последова-тельностями не превышали 15 %), в противоположном случае (совпадений менее 70 % ), на по-следовательностях большой длины, алгоритм существенно замедляется.
Для подобного случая хорошо подходит другой алгоритм Hunt’а и Szimansk’ого[3], кото-рый хорошо работает на последовательностях различающихся друг от друга или последователь-ностях, алфавит которых достаточно широк.
Список литературы:
1. Eugene M., An O(ND) Difference Algorithm and its Variations, www.algolist.manual.ru
2. Окулов С.М. Программирование в алгоритмах. М.: БИНОМ. Лаборатория знаний, 2002.
3. Алгоритм Ханта-Шиманского. www.algolist.manual.ru/Поиск/Общие подпоследова-тельности. Дистанция/ Алгоритм Ханта-Шиманского
Рассмотрим граф на рисунке 1. Матрица