最近在学习后缀树的构造,在网上找了好久发觉国内详解它的构造的文章胜少,在苦苦寻觅了许久,终于发现了一个网友翻译的一篇文章,很好,于是我转帖出来,希望能有更多的人受益,也希望国内多一些英文高手多翻译一些国外的技术文章,好让我们这些英文很烂的人受益,呵呵!
后缀树
Fast String Searching With Suffix Trees
原著
Mark Nelson. Fast string searching with suffix trees. 1996.
构造法
E. Ukkonen. On-line construction of suffix trees. 1995.
翻译
3xian / 三鲜 in GDUT
三鲜 序
原来是打算翻译Sartaj Sahni的Suffix tree, 并专注地进行了一周, 连复习备考的时间也不惜占去. 我希望给国产的同好者提供更通俗易懂的资料, 在翻译的同时对原文进行了删改, 并加入了许多自己的心得. 然而后来发现了Mark Nelson的这篇文章, 相比之下更有亲和力, 于是老实地尽弃前功来翻译这篇. 更重要一个原因是, Mark Nelson介绍的是Ukkonen的构造法O(n), 它比Sartaj Sahni的构造法O(nr), r为字母表大小 在时间上更有优势. 但我们不能说Sartaj Sahni的算法慢, 因为r往往会很小, 因此实际效率也接近线性, 两种构造法在思想上均有可取之处.
本文偏重于阐述后缀树的构造过程, 而并没有直接介绍后缀树除了匹配以外还能做什么. 其实后缀树是一种功能非常强大的数据结构, 你可以去搜索引擎了解一下它还有多少功能, 当然我最希望的是你在阅读本文之后已经足以体会后缀树的妙处, 日后遇到诸多问题的时候都能随心随意地用上.
最后唠叨一句. 我所见过的各种介绍后缀树的论文都难免使初学者陷入混乱, 本文估计也好不到哪里去. 这在一定程度上说明了后缀树的原理是不太浅显的, 理解它需要在整体上把握, 建议希望读者先不要纠结于细节, 思路不清则反复阅读.
问题的来源
字符串匹配问题是程序员经常要面对的问题. 字符串匹配算法的改进可以使许多工程受益良多, 比如数据压缩和DNA排列. 这篇文章讨论的是一种相对鲜为人知的数据结构 --- 后缀树, 并介绍它是如何通过自身的特性去解决一些复杂的匹配问题.
你可以把自己想象成一名工作于DNA排列工程的程序员. 那些基因研究者们天天忙着分切病毒的基因材料, 制造出一段一段的核苷酸序列. 他们把这些序列发到你的服务器里, 指望你在基因数据库中定位. 要知道, 你的数据库里有数百种病毒的数据, 而一个特定的病毒可以有成千上万的碱基. 你的程序必须像C/S工程那样实时向博士们反馈信息, 这需要一个很好的方案.
很明显, 在这个问题上采取暴力算法是极其低效的. 这种方法需要你在基因数据库里对比每一个核苷酸, 测试一个较长的基因段基本会把你的C/S系统变成一台古老的批处理机.
直觉上的解决方法
由于基因数据库一般是不变的, 通过预处理来把搜索简化或许是个好主意. 一种预处理的方法是建立一棵Trie. 我们通过Trie引申出一种东西叫作后缀Trie. (后缀Trie离后缀树仅一步之遥.) 首先, Trie是一种n叉树, n为字母表大小, 每个节点表示从根节点到此节点所经过的所有字符组成的字符串. 而后缀Trie的 “后缀” 说明这棵Trie包含了所给字段的所有后缀 (也许正是一个病毒基因).
图1
BANANAS的后缀Trie
图1展示了文本BANANAS的后缀Trie. 关于这棵Trie有两个地方需要注意. 第一, 从根节点开始, BANANAS的每一个后缀都插入到Trie中, 包括BANANAS, ANANAS, NANAS, ANAS, NAS, AS, S. 第二, 鉴于这种结构, 你可以通过从根节点往下匹配的方式搜索到单词的任何一个子串.
这里所说的第二点正是我们认为后缀Trie优秀的原因. 如果你输入一个长度为N的文本并想在其中搜索一个长度为M的串, 传统的暴力匹配需要进行N*M次字符对比, 而一些改进过的匹配技术, 比如像Boyer-Moore算法, 可以在O(N+M)的时间开销内解决问题, 平均效率更是令人满意. 然而, 后缀Trie亮出了O(M)的牌子, 彻底鄙视了其他算法的成绩, 后缀Trie对比的次数仅仅相当于被搜索串的长度!
这确实是可圈可点的威力, 这意味着你能通过仅仅7次对比便在莎士比亚所有作品中找出BANANAS. 但有一点我们可不能忘了, 构造后缀Trie也是需要时间的.
后缀Trie之所以没有家喻户晓正是因为构造它需要O(n2)的时间和空间. 平方级的开销使它在最需要它的领域 --- 长串搜索 中被拒之门外.
横空出世
直到1976年, Edward McCreigh发表了一篇论文, 咱们的后缀树问世了. 后缀Trie的困境被彻底打破.
后缀树跟后缀Trie有着一样的布局, 但它把只有一个儿子的节点给剔除了. 这个过程被称为路径压缩, 这意味着树上的某些边将表示一个序列而不是单独的字符.
图2
BANANAS的后缀树
图2是由图1的后缀Trie转化而来的后缀树. 你会发现这树基本还是那个形状, 只是节点变少了. 在剔除了只有一个儿子的节点之后, 总节点数由23降为11. 经过证明, 在最坏情况下, 后缀树的节点数也不会超过2N (N为文本的长度). 这使构造后缀树的线性时空开销成为可能.
然而, McCreight最初的构造法是有些缺陷的, 原则上它要按逆序构造, 也就是说字符要从末端开始插入. 如此一来, 便不能作为在线算法, 它变得更加难以应用于实际问题, 如数据压缩.
20年后, 来自赫尔辛基理工大学的Esko Ukkonen把原算法作了一些改动, 把它变成了从左往右. 本文接下来的所有描述和代码都是基于Esko Ukkonen的成果.
对于所给的文本T, Esko Ukkonen的算法是由一棵空树开始, 逐步构造T的每个前缀的后缀树. 比如我们构造BANANAS的后缀树, 先由B开始, 接着是BA, 然后BAN, … . 不断更新直到构造出BANANAS的后缀树.
图3
逐步构造后缀树
初窥门径
加入一个新的前缀需要访问树中已有的后缀. 我们从最长的一个后缀开始(图3中的BAN), 一直访问到最短的后缀(空后缀). 每个后缀会在以下三种节点的其中一种结束.
l 一个叶节点. 这个是常识了, 图4中标号为1, 2, 4, 5的就是叶节点.
l 一个显式节点. 图4中标号为0, 3的是显式节点, 它表示该节点之后至少有两条边.
l 一个隐式节点. 图4中, 前缀BO, BOO, 或者非前缀OO, 它们都在某条表示序列的边上结束, 这些位置就叫作隐式节点. 它表示后缀Trie中存在的由于路径压缩而剔除的节点. 在后缀树的构造过程中, 有时要把一些隐式节点转化为显式节点.
图4
加入BOOK之后的BOOKKEEPER
(也就是BOOK的后缀树)
如图4, 在加入BOOK之后, 树中有5个后缀(包括空后缀). 那么要构造下一个前缀BOOKK的后缀树的话, 只需要访问树中已存在的每一个后缀, 然后在它们的末尾加上K.
前4个后缀BOOK, OOK, OK和K都在叶节点上结束. 由于我们要路径压缩, 只需要在通往叶节点的边上直接加一个字符, 而不需要创建一个新节点.
在所有叶节点更新之后, 我们还需要在空后缀后面加上K. 这时候我们发现已经存在一条从0节点出发的边的首字符为K, 没必要画蛇添足了. 换句话说, 新加入的后缀K可以在0节点和2节点之间的隐式节点中找到. 最终形态见图5.
图5加入BOOKK之后的BOOKKEEPER
相比图4, 树的结构没有发生变化
如果你是一位敏感的读者, 可能要发问了, 如果加入K我们什么都不做的话, 在查找的时候如何知道它到底是一个后缀呢还是某个后缀的一截? 如果你同时又是一位熟悉字符串算法的朋友, 心里可能马上就有答案了 --- 我们只需要在文本后面加个字母表以外的字符, 比如$或者#. 那我们查找到K$或K#的话就说明这是一个后缀了.
稍微麻烦一点的事情
从图4到图5这个更新过程是相对简单的, 其中我们执行了两种更新: 一种是将某条边延长, 另一种是啥都不做. 但接下来往图5继续加入BOOKKE, 我们则会遇到另外两种更新:
1. 创建一个新节点来割开某一隐式节点所处的边, 并在其后加一条新边.
2. 在显式节点后加一条新边.
图6
先分割, 再添加
当我们往图5的树中加入BOOKKE的时候, 我们是从已存在的最长后缀BOOKK开始, 一直操作到最短的后缀空后缀. 更新最长的后缀必然是更新叶节点, 之前提到了, 非常简单. 除此之外, 图5中结束在叶节点上的后缀还有OOKK, OKK, KK. 图6的第一棵树展示了这一类节点的更新.
图5中首个不是结束在叶节点上的后缀是K. 这里我们先引入一个定义:
在每次更新后缀树的过程中, 第一个非叶节点称为激活节点. 它有以下性质:
1. 所有比激活节点长的后缀都在叶节点上结束.
2. 所有在激活节点之后加入的后缀都不在叶节点上结束.
后缀K在边KKE上的隐式节点结束. 在后缀树中我们要判断一个节点是不是非叶节点需要看它是否有跟待加入字符相同的儿子, 即本例中的E.
一眼可以看出, KKE中的第一个K只有一个儿子: K. 所以它是非叶节点(这里同时也是激活节点), 我们要给他加一个儿子来表示E. 这个过程有两个步骤:
1. 在第一个K和第二个K之间把边分割开, 于是第一个K(隐式节点)成了一个显式节点, 如图6第二棵树.
2. 在刚刚变身而来的显式节点后加一个新节点表示E, 如图6第三棵树. 由此我们又多了一个叶节点.
后缀K更新之后, 别忘了还有空后缀. 空后缀在根节点(节点0)结束, 显然此时根节点是一个显式节点. 我们看一下它后面有没有以E开头的边---没有, 那么加入一个新的叶节点(如果存在以E开头的边, 则不用任何操作). 最终如图7.
图7
归纳, 反思, 优化
借助后缀树的特性, 我们可以做出一个相当有效的算法. 首先一个重要的特性是: 一朝为叶, 终生为叶. 一个叶节点自诞生以后绝不会有子孙. 更重要的是, 每当我们往树上加入一个新的前缀, 每一条通往叶节点的边都会延长一个字符(新前缀的最后一个字符). 这使得处理通往叶节点的边变得异常简单, 我们完全可以在创建叶节点的时候就把当前字符到文本末的所有字符一股脑塞进去. 是的, 我们不需要知道后面的字符是啥, 但我们知道它们最终都要被加进去. 因此, 一个叶节点诞生的时候, 也正是它可以被我们遗忘的时候. 你可能会担心通往叶节点的边被分割了怎么办, 那也不要紧, 分割之后只是起点变了, 尾部该怎么着还是怎么着.
如此一来, 我们只需要关心显式节点和隐式节点上的更新.
还要提到一个节约时间的方法. 当我们遍历所有后缀时, 如果某个后缀的某个儿子跟待加字符(新前缀最后一个字符)相同, 那么我们当前前缀的所有更新就可以停止了. 如果你理解了后缀树的本质, 你会知道一旦待加字符跟某个后缀的某个儿子相同, 那么更短的后缀必然也有这个儿子. 我们不妨把首个这样的节点定义为结束节点. 比结束节点长的后缀必然是叶节点, 这一点很好解释, 要么本来就是叶节点, 要么就是新创建的节点(新创建的必然是叶节点). 这意味着, 每一个前缀更新完之后, 当前的结束节点将成为下一轮更新的激活节点.
好了, 现在我们可以把后缀树的更新限制在激活节点和结束节点之间, 效率有了很大的改善. 整理成伪代码如下:
C:
1. Update( 新前缀 )
2. {
3. 当前后缀 = 激活节点
4. 待加字符 = 新前缀最后一个字符
5. done = false;
6. while ( !done ) {
7. if ( 当前后缀在显式节点结束 ) {
8. if ( 当前节点后没有以待加字符开始的边 )
9. 在当前节点后创建一个新的叶节点
10. else
11. done = true;
12. } else {
13. if ( 当前隐式节点的下一个字符不是待加字符 ) {
14. 从隐式节点后分割此边
15. 在分割处创建一个新的叶节点
16. } else
17. done = true;
18. if ( 当前后缀是空后缀 )
19. done = true;
20. else
21. 当前后缀 = 下一个更短的后缀
22. }
23. 激活节点 = 当前后缀
24. }
后缀指针
上面的伪代码看上去很完美, 但它掩盖了一个问题. 注意到第21行, “下一个更短的后缀”, 如果呆板地沿着树枝去搜索我们想要的后缀, 那这种算法就不是线性的了. 要解决此问题, 我们得附加一种指针: 后缀指针. 后缀指针存在于每个结束在非叶节点的后缀上, 它指向“下一个更短的后缀”. 即, 如果一个后缀表示文本的第0到第N个字符, 那么它的后缀指针指向的节点表示文本的第1到第N个字符.
图8是文本ABABABC的后缀树. 第一个后缀指针在表示ABAB的节点上. ABAB的后缀指针指向表示BAB的节点. 同样地, BAB也有它的后缀指针, 指向AB. 如此这般.
图8
加上后缀指针(虚线)的ABABABC的后缀树
介绍一下如何创建后缀指针. 后缀指针的创建是跟后缀树的更新同步的. 随着我们从激活节点移动到结束节点, 我把每个新的叶节点的父亲的路径保存下来. 每当创建一条新边, 我同时也在上一个叶节点的父亲那儿创建一个后缀指针来指向当前新边开始的节点. (显然, 我们不能在第一条新边上做这样的操作, 但除此之外都可以这么做.)
有了后缀指针, 就可以方便地一个后缀跳到另一个后缀. 这个关键性的附加品使得算法的时间上限成功降为O(N).
参考文献
E.M. McCreight. A space-economical suffix tree construction algorithm. Journal of the ACM, 23:262-272, 1976.
E. Ukkonen. On-line construction of suffix trees. Algorithmica, 14(3):249-260, September 1995.
1 #include <iostream> 2 #include <stdio.h> 3 #include <string> 4 5 const int oo = 1<<25; 6 const int ALPHABET_SIZE = 256; 7 const int MAXN = 5000; 8 9 using namespace std; 10 11 int root, last_added, pos, needSL, remainder, 12 active_node, active_e, active_len; 13 14 struct node { 15 /* 16 There is no need to create an "Edge" struct. 17 Information about the edge is stored right in the node. 18 [start; end) interval specifies the edge, 19 by which the node is connected to its parent node. 20 */ 21 22 int start, end, slink; 23 int next[ALPHABET_SIZE]; 24 25 int edge_length() { 26 return min(end, pos + 1) - start; 27 } 28 }; 29 30 node tree[2*MAXN]; 31 char text[MAXN]; 32 33 int new_node(int start, int end = oo) { 34 node nd; 35 nd.start = start; 36 nd.end = end; 37 nd.slink = 0; 38 for (int i = 0; i < ALPHABET_SIZE; i++) 39 nd.next[i] = 0; 40 tree[++last_added] = nd; 41 return last_added; 42 } 43 44 char active_edge() { 45 return text[active_e]; 46 } 47 48 void add_SL(int node) { 49 if (needSL > 0) tree[needSL].slink = node; 50 needSL = node; 51 } 52 53 bool walk_down(int node) { 54 if (active_len >= tree[node].edge_length()) { 55 active_e += tree[node].edge_length(); 56 active_len -= tree[node].edge_length(); 57 active_node = node; 58 return true; 59 } 60 return false; 61 } 62 63 void st_init() { 64 needSL = 0, last_added = 0, pos = -1, 65 remainder = 0, active_node = 0, active_e = 0, active_len = 0; 66 root = active_node = new_node(-1, -1); 67 } 68 69 void st_extend(char c) { 70 text[++pos] = c; 71 needSL = 0; 72 remainder++; 73 while(remainder > 0) { 74 if (active_len == 0) active_e = pos; 75 if (tree[active_node].next[active_edge()] == 0) { 76 int leaf = new_node(pos); 77 tree[active_node].next[active_edge()] = leaf; 78 add_SL(active_node); //rule 2 79 } else { 80 int nxt = tree[active_node].next[active_edge()]; 81 if (walk_down(nxt)) continue; //observation 2 82 if (text[tree[nxt].start + active_len] == c) { //observation 1 83 active_len++; 84 add_SL(active_node); //observation 3 85 break; 86 } 87 int split = new_node(tree[nxt].start, tree[nxt].start + active_len); 88 tree[active_node].next[active_edge()] = split; 89 int leaf = new_node(pos); 90 tree[split].next[c] = leaf; 91 tree[nxt].start += active_len; 92 tree[split].next[text[tree[nxt].start]] = nxt; 93 add_SL(split); //rule 2 94 } 95 remainder--; 96 if (active_node == root && active_len > 0) { //rule 1 97 active_len--; 98 active_e = pos - remainder + 1; 99 } else 100 active_node = tree[active_node].slink > 0 ? tree[active_node].slink : root; //rule 3 101 } 102 } 103 104 int main() { 105 // 106 return 0; 107 }