使用后缀数组和 LCP(-LR) 实现字符串模式匹配
Implementation of string pattern matching using Suffix Array and LCP(-LR)
在过去的几周里,我试图找出如何有效地在另一个字符串中找到一个字符串模式。
我发现很长一段时间以来,最有效的方法一直是使用后缀树。然而,由于这种数据结构在空间上非常昂贵,我进一步研究了后缀数组的使用(它使用的空间要少得多)。 "Suffix Arrays: A new method for on-line string searches" (Manber & Myers, 1993) 等不同的论文指出,搜索子串可以在 O(P+log(N)) 中实现(其中 P 是模式的长度,N 是模式的长度字符串)通过使用二进制搜索和后缀数组以及 LCP 数组。
为了理解搜索算法,特意研究了后面的论文。 This answer did a great job in helping me understand the algorithm (and incidentally made it into the LCP Wikipedia Page).
但我仍在寻找实现此算法的方法。尤其是上面提到的LCP-LR阵列的构建,显得非常复杂。
参考文献:
曼伯和迈尔斯,1993 年:曼伯、乌迪; Myers, Gene, SIAM Journal on Computing, 1993, Vol.22(5), pp.935-948, http://epubs.siam.org/doi/pdf/10.1137/0222058
更新 1
只是强调一下我感兴趣的内容:我了解 LCP 阵列并找到了实现它们的方法。但是,"plain" LCP 数组不适用于高效的模式匹配(如参考资料中所述)。因此,我对实现 LCP-LR 阵列很感兴趣,这似乎比仅实现 LCP 阵列复杂得多
更新 2
已将 link 添加到参考论文
Here 是一个很好的 post,包括一些代码以帮助您更好地理解 LCP 阵列和比较实现。
我了解您的需求是代码,而不是实现您自己的代码。
虽然写在Javathis is an implementation of Suffix Array with LCP by Sedgewick and Wayne from their Algorithms booksite。它应该可以为您节省一些时间,并且应该不会很难移植到 C/C++。
LCP array construction 对于那些可能需要有关该算法的更多信息的人来说是伪造的。
这是一个相当简单的 C++ 实现,尽管 build()
过程在 O(N lg^2 N)
时间内构建后缀数组。 lcp_compute()
过程具有线性复杂度。我在许多编程比赛中使用过这段代码,它从未让我失望:)
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;
const int MAX = 200005;
char str[MAX];
int N, h, sa[MAX], pos[MAX], tmp[MAX], lcp[MAX];
bool compare(int i, int j) {
if(pos[i] != pos[j]) return pos[i] < pos[j]; // compare by the first h chars
i += h, j += h; // if prefvious comparing failed, use 2*h chars
return (i < N && j < N) ? pos[i] < pos[j] : i > j; // return results
}
void build() {
N = strlen(str);
for(int i=0; i<N; ++i) sa[i] = i, pos[i] = str[i]; // initialize variables
for(h=1;;h<<=1) {
sort(sa, sa+N, compare); // sort suffixes
for(int i=0; i<N-1; ++i) tmp[i+1] = tmp[i] + compare(sa[i], sa[i+1]); // bucket suffixes
for(int i=0; i<N; ++i) pos[sa[i]] = tmp[i]; // update pos (reverse mapping of suffix array)
if(tmp[N-1] == N-1) break; // check if done
}
}
void lcp_compute() {
for(int i=0, k=0; i<N; ++i)
if(pos[i] != N-1) {
for(int j=sa[pos[i]+1]; str[i+k] == str[j+k];) k++;
lcp[pos[i]] = k;
if(k) k--;
}
}
int main() {
scanf("%s", str);
build();
for(int i=0; i<N; ++i) printf("%d\n", sa[i]);
return 0;
}
注:如果想让build()
过程的复杂度变成O(N lg N)
,可以用基数排序代替STL排序,但是这样会使代码复杂化。
编辑: 抱歉,我误解了你的问题。虽然我还没有用后缀数组实现字符串匹配,但我想我可以为您描述一个简单的非标准但相当有效的字符串匹配算法。给定两个字符串,text
和 pattern
。给定这些字符串你创建一个新的,我们称之为 concat
,它是两个给定字符串的连接(首先是 text
,然后是 pattern
)。你 运行 concat
上的后缀数组构建算法,你构建了普通的 lcp 数组。然后,你在刚刚构建的后缀数组中搜索长度为 pattern.size()
的后缀。让我们调用它在后缀数组中的位置 pos
。然后你需要两个指针 lo
和 hi
。开始 lo = hi = pos
。在 lcp(lo, pos) = pattern.size()
时减少 lo
,在 lcp(hi, pos) = pattern.size()
时增加 hi
。然后,您在 [lo, hi]
范围内搜索长度至少为 2*pattern.size()
的后缀。如果你找到它,你就找到了一个匹配项。否则,不存在匹配项。
编辑[2]:一旦我有一个实现我会回来...
编辑[3]:
这里是:
// It works assuming you have builded the concatenated string and
// computed the suffix and the lcp arrays
// text.length() ---> tlen
// pattern.length() ---> plen
// concatenated string: str
bool match(int tlen, int plen) {
int total = tlen + plen;
int pos = -1;
for(int i=0; i<total; ++i)
if(total-sa[i] == plen)
{ pos = i; break; }
if(pos == -1) return false;
int lo, hi;
lo = hi = pos;
while(lo-1 >= 0 && lcp[lo-1] >= plen) lo--;
while(hi+1 < N && lcp[hi] >= plen) hi++;
for(int i=lo; i<=hi; ++i)
if(total-sa[i] >= 2*plen)
return true;
return false;
}
对你有帮助的术语:enchanced suffix array
,用来描述后缀数组和其他各种数组,以取代后缀树(lcp,child)。
这些可以是一些示例:
https://code.google.com/p/esaxx/ ESAXX
http://bibiserv.techfak.uni-bielefeld.de/mkesa/ MKESA
esaxx 似乎可以满足您的要求,此外,它还有示例 enumSubstring.cpp 如何使用它。
如果您查看引用的 paper,它会提到一个有用的 属性 (4.2)
。由于SO不支持数学,这里没有意义复制它。
我已经快速实现了,它使用线段树:
// note that arrSize is O(n)
// int arrSize = 2 * 2 ^ (log(N) + 1) + 1; // start from 1
// LCP = new int[N];
// fill the LCP...
// LCP_LR = new int[arrSize];
// memset(LCP_LR, maxValueOfInteger, arrSize);
//
// init: buildLCP_LR(1, 1, N);
// LCP_LR[1] == [1..N]
// LCP_LR[2] == [1..N/2]
// LCP_LR[3] == [N/2+1 .. N]
// rangeI = LCP_LR[i]
// rangeILeft = LCP_LR[2 * i]
// rangeIRight = LCP_LR[2 * i + 1]
// ..etc
void buildLCP_LR(int index, int low, int high)
{
if(low == high)
{
LCP_LR[index] = LCP[low];
return;
}
int mid = (low + high) / 2;
buildLCP_LR(2*index, low, mid);
buildLCP_LR(2*index+1, mid + 1, high);
LCP_LR[index] = min(LCP_LR[2*index], LCP_LR[2*index + 1]);
}
我认为@Erti-Chris Eelmaa 的算法是错误的。
L ... 'M ... M ... M' ... R
|-----|-----|
左子域和右子域应该都包含M。因此我们不能对LCP-LR阵列进行正常的线段树划分。
代码应该看起来像
def lcp_from_i_j(i, j): # means [i, j] not [i, j)
if (j-i<1) return lcp_2_elem(i, j)
return lcp_merge(lcp_from_i_j(i, (i+j)/2), lcp_from_i_j((i+j)/2, j)
左右子范围重叠。线段树支持范围最小查询。但是,[a,b] 之间的范围 min 不等于 [a,b] 之间的 lcp。 LCP 阵列是连续的,简单的 range-min 是行不通的!
在过去的几周里,我试图找出如何有效地在另一个字符串中找到一个字符串模式。
我发现很长一段时间以来,最有效的方法一直是使用后缀树。然而,由于这种数据结构在空间上非常昂贵,我进一步研究了后缀数组的使用(它使用的空间要少得多)。 "Suffix Arrays: A new method for on-line string searches" (Manber & Myers, 1993) 等不同的论文指出,搜索子串可以在 O(P+log(N)) 中实现(其中 P 是模式的长度,N 是模式的长度字符串)通过使用二进制搜索和后缀数组以及 LCP 数组。
为了理解搜索算法,特意研究了后面的论文。 This answer did a great job in helping me understand the algorithm (and incidentally made it into the LCP Wikipedia Page).
但我仍在寻找实现此算法的方法。尤其是上面提到的LCP-LR阵列的构建,显得非常复杂。
参考文献:
曼伯和迈尔斯,1993 年:曼伯、乌迪; Myers, Gene, SIAM Journal on Computing, 1993, Vol.22(5), pp.935-948, http://epubs.siam.org/doi/pdf/10.1137/0222058
更新 1
只是强调一下我感兴趣的内容:我了解 LCP 阵列并找到了实现它们的方法。但是,"plain" LCP 数组不适用于高效的模式匹配(如参考资料中所述)。因此,我对实现 LCP-LR 阵列很感兴趣,这似乎比仅实现 LCP 阵列复杂得多
更新 2
已将 link 添加到参考论文
Here 是一个很好的 post,包括一些代码以帮助您更好地理解 LCP 阵列和比较实现。
我了解您的需求是代码,而不是实现您自己的代码。 虽然写在Javathis is an implementation of Suffix Array with LCP by Sedgewick and Wayne from their Algorithms booksite。它应该可以为您节省一些时间,并且应该不会很难移植到 C/C++。
LCP array construction 对于那些可能需要有关该算法的更多信息的人来说是伪造的。
这是一个相当简单的 C++ 实现,尽管 build()
过程在 O(N lg^2 N)
时间内构建后缀数组。 lcp_compute()
过程具有线性复杂度。我在许多编程比赛中使用过这段代码,它从未让我失望:)
#include <stdio.h>
#include <string.h>
#include <algorithm>
using namespace std;
const int MAX = 200005;
char str[MAX];
int N, h, sa[MAX], pos[MAX], tmp[MAX], lcp[MAX];
bool compare(int i, int j) {
if(pos[i] != pos[j]) return pos[i] < pos[j]; // compare by the first h chars
i += h, j += h; // if prefvious comparing failed, use 2*h chars
return (i < N && j < N) ? pos[i] < pos[j] : i > j; // return results
}
void build() {
N = strlen(str);
for(int i=0; i<N; ++i) sa[i] = i, pos[i] = str[i]; // initialize variables
for(h=1;;h<<=1) {
sort(sa, sa+N, compare); // sort suffixes
for(int i=0; i<N-1; ++i) tmp[i+1] = tmp[i] + compare(sa[i], sa[i+1]); // bucket suffixes
for(int i=0; i<N; ++i) pos[sa[i]] = tmp[i]; // update pos (reverse mapping of suffix array)
if(tmp[N-1] == N-1) break; // check if done
}
}
void lcp_compute() {
for(int i=0, k=0; i<N; ++i)
if(pos[i] != N-1) {
for(int j=sa[pos[i]+1]; str[i+k] == str[j+k];) k++;
lcp[pos[i]] = k;
if(k) k--;
}
}
int main() {
scanf("%s", str);
build();
for(int i=0; i<N; ++i) printf("%d\n", sa[i]);
return 0;
}
注:如果想让build()
过程的复杂度变成O(N lg N)
,可以用基数排序代替STL排序,但是这样会使代码复杂化。
编辑: 抱歉,我误解了你的问题。虽然我还没有用后缀数组实现字符串匹配,但我想我可以为您描述一个简单的非标准但相当有效的字符串匹配算法。给定两个字符串,text
和 pattern
。给定这些字符串你创建一个新的,我们称之为 concat
,它是两个给定字符串的连接(首先是 text
,然后是 pattern
)。你 运行 concat
上的后缀数组构建算法,你构建了普通的 lcp 数组。然后,你在刚刚构建的后缀数组中搜索长度为 pattern.size()
的后缀。让我们调用它在后缀数组中的位置 pos
。然后你需要两个指针 lo
和 hi
。开始 lo = hi = pos
。在 lcp(lo, pos) = pattern.size()
时减少 lo
,在 lcp(hi, pos) = pattern.size()
时增加 hi
。然后,您在 [lo, hi]
范围内搜索长度至少为 2*pattern.size()
的后缀。如果你找到它,你就找到了一个匹配项。否则,不存在匹配项。
编辑[2]:一旦我有一个实现我会回来...
编辑[3]:
这里是:
// It works assuming you have builded the concatenated string and
// computed the suffix and the lcp arrays
// text.length() ---> tlen
// pattern.length() ---> plen
// concatenated string: str
bool match(int tlen, int plen) {
int total = tlen + plen;
int pos = -1;
for(int i=0; i<total; ++i)
if(total-sa[i] == plen)
{ pos = i; break; }
if(pos == -1) return false;
int lo, hi;
lo = hi = pos;
while(lo-1 >= 0 && lcp[lo-1] >= plen) lo--;
while(hi+1 < N && lcp[hi] >= plen) hi++;
for(int i=lo; i<=hi; ++i)
if(total-sa[i] >= 2*plen)
return true;
return false;
}
对你有帮助的术语:enchanced suffix array
,用来描述后缀数组和其他各种数组,以取代后缀树(lcp,child)。
这些可以是一些示例:
https://code.google.com/p/esaxx/ ESAXX
http://bibiserv.techfak.uni-bielefeld.de/mkesa/ MKESA
esaxx 似乎可以满足您的要求,此外,它还有示例 enumSubstring.cpp 如何使用它。
如果您查看引用的 paper,它会提到一个有用的 属性 (4.2)
。由于SO不支持数学,这里没有意义复制它。
我已经快速实现了,它使用线段树:
// note that arrSize is O(n)
// int arrSize = 2 * 2 ^ (log(N) + 1) + 1; // start from 1
// LCP = new int[N];
// fill the LCP...
// LCP_LR = new int[arrSize];
// memset(LCP_LR, maxValueOfInteger, arrSize);
//
// init: buildLCP_LR(1, 1, N);
// LCP_LR[1] == [1..N]
// LCP_LR[2] == [1..N/2]
// LCP_LR[3] == [N/2+1 .. N]
// rangeI = LCP_LR[i]
// rangeILeft = LCP_LR[2 * i]
// rangeIRight = LCP_LR[2 * i + 1]
// ..etc
void buildLCP_LR(int index, int low, int high)
{
if(low == high)
{
LCP_LR[index] = LCP[low];
return;
}
int mid = (low + high) / 2;
buildLCP_LR(2*index, low, mid);
buildLCP_LR(2*index+1, mid + 1, high);
LCP_LR[index] = min(LCP_LR[2*index], LCP_LR[2*index + 1]);
}
我认为@Erti-Chris Eelmaa 的算法是错误的。
L ... 'M ... M ... M' ... R
|-----|-----|
左子域和右子域应该都包含M。因此我们不能对LCP-LR阵列进行正常的线段树划分。 代码应该看起来像
def lcp_from_i_j(i, j): # means [i, j] not [i, j)
if (j-i<1) return lcp_2_elem(i, j)
return lcp_merge(lcp_from_i_j(i, (i+j)/2), lcp_from_i_j((i+j)/2, j)
左右子范围重叠。线段树支持范围最小查询。但是,[a,b] 之间的范围 min 不等于 [a,b] 之间的 lcp。 LCP 阵列是连续的,简单的 range-min 是行不通的!