Boyer-Moore
本章节内容需要以 《前缀函数与 KMP 算法》 作为前置章节。
之前的 KMP 算法将前缀匹配的信息用到了极致,
而 BM 算法背后的基本思想是通过后缀匹配获得比前缀匹配更多的信息来实现更快的字符跳转。
基础介绍¶
想象一下,如果我们的的模式字符串
在这里做定义,往后不赘述:
假如我们知道了
观察 1:
如果我们知道
观察 2:
更一般地,如果出现在
那么就可以不用匹配,直接将
因此除非
注意:显然这个表只需计算到
现在假设
如果是,就继续回退直到整个模式串
或者,我们也可能会在匹配完
观察 3(a):
在 观察 2 中提到,当匹配完
需要把
而
所以我们的注意力应该沿着
然而,我们有机会跳过更多的字符,请继续看下去。
观察 3(b):
如果我们知道
我们还知道在
使得失配字符
假设此时
假定
所以有:
于是我们在失配时,可以把把
实例说明:¶
箭头指向失配字符
根据 观察 2,我们需要将
现在char:
根据 观察 3(a),
这时
显然直观上看,此时根据 观察 3(b),将
而从形式化逻辑看,此时,
现在我们发现了
算法设计¶
最初的匹配算法¶
现在看这样一个利用
如果上面的算法
然后让我们更精细地描述下计算
根据前文定义,
也就是说需要找到一个最好的
- 当
k<0 pat delta_2 - 当
k>0 pat[k-1]=pat[j] pat[k\dots k+patlastpos-j-1] subpat pat[j] pat k pat[k-1]
还要注意两个限制条件:
k < j k=j pat[k]=pat[j] pat[j] pat[k] - 考虑到
delta_2(patlastpos)= 0 rpr(patlastpos) = patlastpos
由于理解
对于
对于
对于
对于
对于
对于
对于
对于
对于
现在再看一下另一个例子:
对于
对于
对于
对于
对于
对于
对于
对于
对于
对匹配算法的一个改进¶
最后,实践过程中考虑到搜索过程中估计有 80% 的时间用在了基于 观察 1 的跳转上,也就是
于是,可以为此进行特别的优化:
我们定义一个
用
其中
经过改进,比起原算法,在做 观察 1 跳转时不必每次进行
delta_2 构建细节¶
历史细节¶
说起
而构造
-
1969 年夏天 Morris 为某个大型机编写文本编辑器时利用有限自动机的理论发明了等价于 KMP 算法的字符串匹配算法,而他的算法由于过于复杂,被不理解他算法的同事当做 bug 修改得一团糟,哈哈。
-
1970 年 KMP 中的“带头人”Knuth 在研究 Cook 的关于两路确定下推自动机(two-way deterministic pushdown automaton)的理论时受到启发,也独立发明了 KMP 算法的雏形,并把它展示给他的同事 Pratt,Pratt 改进了算法的数据结构。
-
1974 年 Boyer、Moor 发现通过更快地跳过不可能匹配的文本能实现比 KMP 更快的字符串匹配,(Gosper 也独立地发现了这一点),而一个只有原始
delta_1 -
1975 年 Boyer、Moor 提出了原始的
delta_2 delta_2 delta_2 delta_1 delta_2 -
1976 年 1 月 Knuth 证明了关于
delta_2 delta_2 -
Standish 又为 Boyer、Moor 提供了现在的匹配算法的改进。
-
1977 年 6 月 Knuth、Morris、Pratt 正式联合发表了 KMP 算法的论文,其中在提及比 KMP 表现更好的算法中提出了
delta_2
这个 BM 算法的发展的故事,切实地向我们展示了团结、友谊、协作,以及谦虚好学不折不挠“在平凡中实现伟大”!😂😂😂
时间复杂度为 O(n^3) 的构建 delta_2 的朴素算法¶
在介绍 Knuth 的
- 对于
[0, patlen)
区间的每一个位置i
,根据subpat
的长度确定其重现位置的区间,也就是[-subpatlen, i]
; - 可能的重现位置按照从右到左进行逐字符比较,寻找符合
delta_2 subpat - 最后别忘了令
delta_2(lastpos)= 0
use std::cmp::PartialEq;
pub fn build_delta_2_table_naive(p: &[impl PartialEq]) -> Vec<usize> {
let patlen = p.len();
let lastpos = patlen - 1;
let mut delta_2 = vec![];
for i in 0..patlen {
let subpatlen = (lastpos - i) as isize;
if subpatlen == 0 {
delta_2.push(0);
break;
}
for j in (-subpatlen..(i + 1) as isize).rev() {
// subpat 匹配
if (j..j + subpatlen)
.zip(i + 1..patlen)
.all(|(rpr_index, subpat_index)| {
if rpr_index < 0 {
return true;
}
if p[rpr_index as usize] == p[subpat_index] {
return true;
}
false
})
&& (j <= 0 || p[(j - 1) as usize] != p[i])
{
delta_2.push((lastpos as isize - j) as usize);
break;
}
}
}
delta_2
}
特别地,对 Rust 语言特性进行必要地解释,下不赘述:
usize
和isize
是和内存指针同字节数的无符号整数和有符号整数,在 32 位机上相当于u32
和i32
,64 位机上相当于u64
和i64
。- 索引数组、向量、分片时使用
usize
类型的数字(因为在做内存上的随机访问并且下标不能为负值),所以如果需要处理负值要用isize
,而进行索引时又要用usize
,这就看到使用as
关键字进行二者之间的显式转换。 impl PartialEq
只是用作泛型,可以同时支持Unicode
编码的char
和二进制的u8
。
显见这是个时间复杂度为
时间复杂度为 O(n) 的构建 delta_2 的高效算法¶
下面我们要介绍的是时间复杂度为
需要指出的是,虽然 1977 年 Knuth 提出了这个构建方法,然而他的原始版本的构建算法存在一个缺陷,实际上对于某些
Rytter 在 1980 年SIAM Journal on Computing上发表的文章3对此提出了修正,但是 Rytter 的这篇文章在细节上有些令人疑惑的地方,包括不限于:
- 示例中奇怪的
delta_2 delta_2 delta_2 - 明显的在复述 Knuth 算法时的笔误、算法上错误的缩进(可能是文章录入时的问题?)
- 奇妙的变量命名(考虑到那个时代的标签:
goto
语句、汇编语言、大型机,随性的变量命名也很合理)
总之就是你绝对不想看他的那个修正算法的具体实现,不过好在他在用文字描述的时候比用伪代码清晰多了呢,现在我们用更清晰的思路和代码结构整理这么一个
首先考虑到
按照重现位置由远到近,也就是偏移量由大到小,分成如下几类:
-
整个
subpat pat \texttt{[(EYX)]ABYXCDEYX} delta_2(j) = patlastpos\times 2 - j -
subpat pat pat \texttt{[(XX)ABC]XXXABC} patlastpos < delta_2(j) < patlastpos\times 2 - j subpat pat \texttt{[ABC]XXXABC} patlastpos = delta_2(j) -
subpat pat \texttt{AB[YX]CDEYX} delta_2(j) < patlastpos
现在来讨论如何高效地计算这三种情况:
第一种情况¶
这是最简单的情况,只需一次遍历并且可以顺便将
第二种情况¶
我们观察什么时候会出现
比如之前的例子:
说到这个拗口的前缀后缀相等,此时看过之前《前缀函数与 KMP 算法》的小伙伴们可能已经有所悟了,
没错,实际上对第二种和第三种情况的计算的关键都离不开前缀函数的计算和和应用
那么只要
而当
可以计算此时的
设此时这对相等的前后缀长度为
而
那么问题到这儿是不是结束了呢,并不是,因为可能会有多对相等的前缀和后缀,比如:
在
之前提到的 Knuth 算法的缺陷就是只考虑了最长的那一对的情况。
所以实际上我们要考虑所有
不同的
从
第三种情况¶
也就是按照从右到左的顺序,在
开启脑洞:既然是个字符串搜索的问题,那么当然可以用著名的 BM 算法本身解决,于是我们就得到了一个 BM 的递归实现的第三种情况,结束条件是
而且根据
这就很好地启发了我们(起码很好地启发了 Knuth)使用类似于计算前缀函数的过程计算第三种情况,只不过是左右反过来的前缀函数:
- 两个指针分别指向子串的左端点和子串最长公共前后缀的“前缀”位置,从右向左移动,在发现指向的两个字符相等时继续移动,此时相当于“前缀”变大;
- 当两个字符不相等时,之前相等的部分就满足了
delta_2
同前缀函数一样,需要一个辅助数组,用于回退,可以使用之前计算第二种情况所生成的前缀数组的空间。
上述实现¶
use std::cmp::PartialEq;
use std::cmp::min;
pub fn build_delta_2_table_improved_minghu6(p: &[impl PartialEq]) -> Vec<usize> {
let patlen = p.len();
let lastpos = patlen - 1;
let mut delta_2 = Vec::with_capacity(patlen);
// 第一种情况
// delta_2[j] = lastpos * 2 - j
for i in 0..patlen {
delta_2.push(lastpos * 2 - i);
}
// 第二种情况
// lastpos <= delata2[j] = lastpos * 2 - j
let pi = compute_pi(p); // 计算前缀函数
let mut i = lastpos;
let mut last_i = lastpos; // 只是为了初始化
while pi[i] > 0 {
let start;
let end;
if i == lastpos {
start = 0;
} else {
start = patlen - pi[last_i];
}
end = patlen - pi[i];
for j in start..end {
delta_2[j] = lastpos * 2 - j - pi[i];
}
last_i = i;
i = pi[i] - 1;
}
// 第三种情况
// delata2[j] < lastpos
let mut j = lastpos;
let mut t = patlen;
let mut f = pi;
loop {
f[j] = t;
while t < patlen && p[j] != p[t] {
delta_2[t] = min(delta_2[t], lastpos - 1 - j); // 使用min函数保证后面可能的回退不会覆盖前面的数据
t = f[t];
}
t -= 1;
if j == 0 {
break;
}
j -= 1;
}
// 没有实际意义,只是为了完整定义
delta_2[lastpos] = 0;
delta_2
}
Galil 规则对多次匹配时最坏情况的改善¶
关于后缀匹配算法的多次匹配问题¶
之前的搜索算法只涉及到在
如何利用之前匹配成功的字符的信息,将最坏情况下的时间复杂度降为线性。
在原始的成功匹配后,简单的
比如一个极端的例子:
对此 Knuth 提出来的一个方法是用一个“数量有限”的状态的集合来记录
而 Knuth 提出的另一个方法,嗯这里就不介绍了,同在上面的 Knuth 那篇“宝藏”论文里被介绍,缺点是除了过于复杂以外主要是构建辅助的数据结构需要的预处理时间太大:
于是在这个背景下就有了下面介绍的思路简单,不需要额外预处理开销的 Galil 算法4。
Galil 规则¶
原理很简单,假定一个
比如,
当然其实
事实上,广义地讲,
我们规定这个最短的周期为
在搜索过程中,假如我们的
而如何计算这个最短周期的长度呢,假如我们知道
而从数学的角度看这个公式,显然我们已经有了长度为
而这个最长相等的前后缀长度,
结合上述优化的 BM 的搜索算法最终实现¶
#[cfg(target_pointer_width = "64")]
const LARGE: usize = 10_000_000_000_000_000_000;
#[cfg(not(target_pointer_width = "64"))]
const LARGE: usize = 2_000_000_000;
pub struct BMPattern<'a> {
pat_bytes: &'a [u8],
delta_1: [usize; 256],
delta_2: Vec<usize>,
k: usize // pat的最短周期长度
}
impl<'a> BMPattern<'a> {
// ...
pub fn find_all(&self, string: &str) -> Vec<usize> {
let mut result = vec![];
let string_bytes = string.as_bytes();
let stringlen = string_bytes.len();
let patlen = self.pat_bytes.len();
let pat_last_pos = patlen - 1;
let mut string_index = pat_last_pos;
let mut pat_index;
let l0 = patlen - self.k;
let mut l = 0;
while string_index < stringlen {
let old_string_index = string_index;
while string_index < stringlen {
string_index += self.delta0(string_bytes[string_index]);
}
if string_index < LARGE {
break;
}
string_index -= LARGE;
// 如果string_index发生移动,意味着自从上次成功匹配后发生了至少一次的失败匹配。
// 此时需要将Galil规则的二次匹配的偏移量归零。
if old_string_index < string_index {
l = 0;
}
pat_index = pat_last_pos;
while pat_index > l && string_bytes[string_index] == self.pat_bytes[pat_index] {
string_index -= 1;
pat_index -= 1;
}
if pat_index == l && string_bytes[string_index] == self.pat_bytes[pat_index] {
result.push(string_index - l);
string_index += pat_last_pos - l + self.k;
l = l0;
} else {
l = 0;
string_index += max(
self.delta_1[string_bytes[string_index] as usize],
self.delta_2[pat_index],
);
}
}
result
}
}
最坏情况在实践中性能影响¶
从实践的角度上说,理论上的最坏情况并不容易影响性能表现,哪怕是很小的只有 4 的字符集的随机文本测试下这种最坏情况的影响也小到难以观察。
也因此如果没有很好地设计,使用 Galil 法则会拖累一点平均的性能表现,但对于一些极端特殊的
实践及后续¶
这个部分要讨论实践中的具体问题。
尽管前面给出了一些算法的实现代码,但并没有真正讨论过完整实现可能面临的一些“小问题”。
字符类型的考虑¶
在英语环境下,特别是上世纪 70 年代那个时候,人们考虑字符,默认的前提是它是 ASCII 码,通用字符表是容易通过一个固定大小的数组来确定的。
而在尝试用 Rust 实现上述算法的时候,第一个遇到的问题是字符的问题,用一门很新的 2010 + 发展起来的语言来实现 1970 + 时代的算法,是一件很有意思的事情。
会观察到一些因时代发展而产生的一些变化,现代的编程语言,内生的 char
类型就是 Unicode,首先不可能用一个全字符集大小的数组来计算
但更严重的问题是 Unicode 使用的都是变长的字节编码方案,所以没办法直接按照字符个数计算跳转的字节数,当然,如果限定文本是简单的 ASCII 字符集,我们仍然可以按照 1 字符宽 1 字节来进行快速跳转,但这样的实现根本就没啥卵用!😠
在思考的过程中,首先的一个想法是直接将字符串转为按字符索引的向量数组,但这意味着啥都不用做就先有了一个遍历字符串的时间开销,和额外的大于等于字符串字节数的额外空间开销(因为 char
类型是 Unicode 字面值,采用固定 4 字节大小保存)。
于是我改进了思路,对于变长编码字符串,至少要完全遍历一遍,才能完成字符串匹配,那么在遍历过程中,我使用一个基于可增长的环结构实现的双头队列作为滑动窗口,保存过去
于是挫折使我困惑,困惑使我思考,终于一束阳光照进了石头缝里:
- 字符串匹配算法高效的关键在于字符索引的快速跳转
- 字符索引一定要建立在等宽字符的基础上,
基于这两条原则思考,我就发现二进制字节本身:1 字节等宽、字符全集大小是 256,就是符合条件的完美字符!在这个基础上完成了一系列后缀匹配算法的高效实现。
Simplified Boyer-Moore 算法¶
BM 算法最复杂的地方就在于
Boyer-Moore-Horspol 算法¶
Horspol 算法同样是基于坏字符的规则,不过是在与
pub struct HorspoolPattern<'a> {
pat_bytes: &'a [u8],
bm_bc: [usize; 256],
}
impl<'a> HorspoolPattern<'a> {
// ...
pub fn find_all(&self, string: &str) -> Vec<usize> {
let mut result = vec![];
let string_bytes = string.as_bytes();
let stringlen = string_bytes.len();
let pat_last_pos = self.pat_bytes.len() - 1;
let mut string_index = pat_last_pos;
while string_index < stringlen {
if &string_bytes[string_index-pat_last_pos..string_index+1] == self.pat_bytes {
result.push(string_index-pat_last_pos);
}
string_index += self.bm_bc[string_bytes[string_index] as usize];
}
result
}
}
Boyer-Moore-Sunday 算法¶
Sunday 算法同样是利用坏字符规则,只不过相比 Horspool 它更进一步,直接关注
实现它只需要稍微修改一下
Sunday 算法通常用作一般情况下实现最简单而且平均表现最好之一的实用算法,通常表现比 Horspool、BM 都要快一点。
pub struct SundayPattern<'a> {
pat_bytes: &'a [u8],
sunday_bc: [usize; 256],
}
impl<'a> SundayPattern<'a> {
// ...
fn build_sunday_bc(p: &'a [u8]) -> [usize; 256] {
let mut sunday_bc_table = [p.len() + 1; 256];
for i in 0..p.len() {
sunday_bc_table[p[i] as usize] = p.len() - i;
}
sunday_bc_table
}
pub fn find_all(&self, string: &str) -> Vec<usize> {
let mut result = vec![];
let string_bytes = string.as_bytes();
let pat_last_pos = self.pat_bytes.len() - 1;
let stringlen = string_bytes.len();
let mut string_index = pat_last_pos;
while string_index < stringlen {
if &string_bytes[string_index - pat_last_pos..string_index+1] == self.pat_bytes {
result.push(string_index - pat_last_pos);
}
if string_index + 1 == stringlen {
break;
}
string_index += self.sunday_bc[string_bytes[string_index + 1] as usize];
}
result
}
}
BMHBNFS 算法¶
该算法结合了 Horspool 和 Sunday,是 CPython 实现 stringlib
模块时用到的 find
的算法5,似乎国内更有名气,不清楚为何叫这个名字,怎么就“AKA”了?
以下简称 B5S。
B5S 基本想法是:
-
按照后缀匹配的思路,首先比较
patlastpos 0\dots patlastpos-1 -
如果任何一个阶段发生不匹配,就进入跳转阶段;
-
在跳转阶段,首先观察
patlastpos pat patlen+1 如果这个字符在
pat patlastpos delta_1
而这个算法根据时间节省还是空间节省为第一目标,会有差别巨大的不同实现。
时间节省版本¶
pub struct B5STimePattern<'a> {
pat_bytes: &'a [u8],
alphabet: [bool;256],
bm_bc: [usize;256],
k: usize
}
impl<'a> B5STimePattern<'a> {
pub fn new(pat: &'a str) -> Self {
assert_ne!(pat.len(), 0);
let pat_bytes = pat.as_bytes();
let (alphabet, bm_bc, k) = B5STimePattern::build(pat_bytes);
B5STimePattern { pat_bytes, alphabet, bm_bc, k }
}
fn build(p: &'a [u8]) -> ([bool;256], [usize;256], usize) {
let mut alphabet = [false;256];
let mut bm_bc = [p.len(); 256];
let lastpos = p.len() - 1;
for i in 0..lastpos {
alphabet[p[i] as usize] = true;
bm_bc[p[i] as usize] = lastpos - i;
}
alphabet[p[lastpos] as usize] = true;
(alphabet, bm_bc, compute_k(p))
}
pub fn find_all(&self, string: &str) -> Vec<usize> {
let mut result = vec![];
let string_bytes = string.as_bytes();
let pat_last_pos = self.pat_bytes.len() - 1;
let patlen = self.pat_bytes.len();
let stringlen = string_bytes.len();
let mut string_index = pat_last_pos;
let mut offset = pat_last_pos;
let offset0 = self.k - 1;
while string_index < stringlen {
if string_bytes[string_index] == self.pat_bytes[pat_last_pos] {
if &string_bytes[string_index-offset..string_index] == &self.pat_bytes[pat_last_pos-offset..pat_last_pos] {
result.push(string_index-pat_last_pos);
offset = offset0;
// Galil rule
string_index += self.k;
continue;
}
}
if string_index + 1 == stringlen {
break;
}
offset = pat_last_pos;
if !self.alphabet[string_bytes[string_index+1] as usize] {
string_index += patlen + 1; // sunday
} else {
string_index += self.bm_bc[string_bytes[string_index] as usize]; // horspool
}
}
result
}
}
这个版本的 B5S 性能表现非常理想,是通常情况下,目前介绍的后缀匹配系列算法中最快的。
空间节省版本¶
这也是 CPython stringlib
中实现的版本,使用了两个整数近似取代了字符表和
用一个简单的 Bloom 过滤器取代字符表(alphabet)¶
pub struct BytesBloomFilter {
mask: u64,
}
impl BytesBloomFilter {
pub fn new() -> Self {
SimpleBloomFilter {
mask: 0,
}
}
fn insert(&mut self, byte: &u8) {
(self.mask) |= 1u64 << (byte & 63);
}
fn contains(&self, char: &u8) -> bool {
(self.mask & (1u64 << (byte & 63))) != 0
}
}
Bloom 过滤器设设计通过牺牲准确率(实际还有运行时间)来极大地节省存储空间的 Set
类型的数据结构,它的特点是会将集合中不存在的项误判为存在(False Positives,简称 FP),但不会把集合中存在的项判断为不存在(False Negatives,简称 FN),因此使用它可能会因为 FP 而没有得到最大的字符跳转,但不会因为 FN 而跳过本应匹配的字符。
理论上分析,上述“Bloom 过滤器”的实现在
当然这不是一个标准的 Bloom 过滤器,首先它实际上没有使用一个真正的哈希函数,实际上它只是一个字符映射,就是将 0-255 的字节映射为它的前六位构成的数,考虑到我们在做内存上的字符搜索,这样的简化就非常重要,因为即使用目前已知最快的非加密哈希算法 xxHash,计算所需要的时间都要比它高一个数量级。
另外,按照计算,当 pat 在 30 字节以下时,为了达到最佳的 FP 概率,需要超过一个哈希函数,但这么做意义不大,因为用装有两个 u128
数字的数组就已经可以构建字符表的全字符集。
使用 delta_1(pat[patlastpos]) 代替整个 delta_1 ¶
观察 skip = delta1(pat[patlastpos])
,
在第一阶段不匹配时,直接向下滑动 skip
个字符;但当第二阶段不配时,因为缺乏整个
pub struct B5SSpacePattern<'a> {
pat_bytes: &'a [u8],
alphabet: BytesBloomFilter,
skip: usize,
}
impl<'a> B5SSpacePattern<'a> {
pub fn new(pat: &'a str) -> Self {
assert_ne!(pat.len(), 0);
let pat_bytes = pat.as_bytes();
let (alphabet, skip) = B5SSpacePattern::build(pat_bytes);
B5SSpacePattern { pat_bytes, alphabet, skip}
}
fn build(p: &'a [u8]) -> (BytesBloomFilter, usize) {
let mut alphabet = BytesBloomFilter::new();
let lastpos = p.len() - 1;
let mut skip = p.len();
for i in 0..p.len()-1 {
alphabet.insert(&p[i]);
if p[i] == p[lastpos] {
skip = lastpos - i;
}
}
alphabet.insert(&p[lastpos]);
(alphabet, skip)
}
pub fn find_all(&self, string: &'a str) -> Vec<usize> {
let mut result = vec![];
let string_bytes = string.as_bytes();
let pat_last_pos = self.pat_bytes.len() - 1;
let patlen = self.pat_bytes.len();
let stringlen = string_bytes.len();
let mut string_index = pat_last_pos;
while string_index < stringlen {
if string_bytes[string_index] == self.pat_bytes[pat_last_pos] {
if &string_bytes[string_index-pat_last_pos..string_index] == &self.pat_bytes[..patlen-1] {
result.push(string_index-pat_last_pos);
}
if string_index + 1 == stringlen {
break;
}
if !self.alphabet.contains(&string_bytes[string_index+1]) {
string_index += patlen + 1; // sunday
} else {
string_index += self.skip; // horspool
}
} else {
if string_index + 1 == stringlen {
break;
}
if !self.alphabet.contains(&string_bytes[string_index+1]) {
string_index += patlen + 1; // sunday
} else {
string_index += 1;
}
}
}
result
}
}
这个版本的算法相对于前面的后缀匹配算法不够快,但差距并不大,仍然比 KMP 这种快得多,特别是考虑到它极为优秀的空间复杂度:至多两个 u64
的整数,这确实是极为实用的适合作为标准库实现的一种算法!
理论分析¶
现在我们通过一个简单的概率模型来做一些绝不枯燥的理论上的分析,借此可以发现一些有趣而更深入的事实。
建立模型¶
想象一下,我们滑动字符串
假如这个代价是用对
而如果这个代价是用机器指令衡量,那我们可以知道平均每个字符需要多少条机器指令;
当然也可以有其他的衡量方式,这并不影响什么,这里我们采用对
同时为我们的概率模型提出一个假设:
显然,假如全字母表的大小为
现在可以更准确地刻画这个比率,
其中,
计算 BoyerMoore 算法的 skip(m,k) ¶
delta_1 ¶
首先考虑
而对于
-
当
k = 1 -
失配字符对应位置的下一个字符恰好等于失配字符;
-
失配字符已经是
pat
-
-
当
1 < k < patlen-m pat -
当
k = patlen - m pat pat -
当
k > patlen - m delta_1
于是有计算
delta_2 ¶
对于
于是
汇总¶
前面已经独立讨论了
当
因此针对
于是通过组合
分析比较¶
为了结构清晰、书写简单、演示方便,我们使用 Python 平台的 Lisp 方言 Hy 来进行实际计算:
myprob.hy
(require [hy.contrib.sequences [defseq seq]])
(import [hy.contrib.sequences [Sequence end-sequence]])
(import [hy.models [HySymbol]])
(defmacro simplify-probfn [patlen p probfn-list]
"(prob-xxx patlen p m k) -> (prob-xxx-s m k)"
(lfor probfn probfn-list
[(setv simplified-probfn-symbol (HySymbol (.format "{}-s" (name probfn))))
`(defn ~simplified-probfn-symbol [&rest args] (~probfn patlen p #*args))]))
(defn map-sum [range-args func]
(setv [start end] range-args)
(-> func (map (range start (inc end))) sum))
(defn cost [m] (+ m 1))
(defn prob-m [patlen p m]
(/
(* (** p m) (- 1 p))
(- 1 (** p patlen))))
(defn prob-delta1 [patlen p m k]
(cond [(= 1 k) (* (** (- 1 p) m)
(if (= (inc m) patlen) 1 p))]
[(< k (- patlen m)) (* p (** (- 1 p) (dec (+ k m))))]
[(= k (- patlen m)) (** (- 1 p) (dec patlen))]
[(> k (- patlen m)) 0]))
(defn prob-delta1-worthless [p m] (- 1 (** (- 1 p) m)))
(defn prob-pr [patlen p m k] (if (< patlen (+ m k))
(* (- 1 p) (** p m))
(** p (- patlen k))))
(defn prob-delta2 [patlen p m]
"prob-delta2(_, k) = prob-delta2-seq[k]"
(defseq prob-delta2-seq [n]
(cond [(< n 1) 0]
[(= n 1) (prob-pr patlen p m 1)]
[(> n 1) (* (prob-pr patlen p m n) (- 1 (sum (take n prob-delta2-seq))))]))
prob-delta2-seq)
(defn prob-delta2-1 [patlen p m]
(defseq prob-delta2-1-seq [n]
(cond [(< n 2) 0]
[(>= n 2) (* (prob-pr patlen p m n) (- 1 (sum (take n prob-delta2-1-seq))))]))
prob-delta2-1-seq)
(defn skip [patlen p m k prob-delta2-seq prob-delta2-1-seq]
(simplify-probfn patlen p [prob-delta1])
(if (= k 1) (* (prob-delta1-s m 1) (get prob-delta2-seq 1))
(sum [(* (prob-delta1-worthless p m) (get prob-delta2-1-seq k))
(* (prob-delta1-s m k) (sum (take k prob-delta2-seq)))
(* (get prob-delta2-seq k) (map-sum [1 (- k 1)]
(fn [n] (prob-delta1-s m n))))
(* (prob-delta1-s m k) (get prob-delta2-seq k))])))
(defn bm-rate [patlen p]
(simplify-probfn patlen p [prob-m prob-delta2 prob-delta2-1 skip])
(/
(map-sum [0 (dec patlen)]
(fn [m] (* (cost m) (prob-m-s m))))
(map-sum [0 (dec patlen)]
(fn [m]
(setv prob-delta2-seq (prob-delta2-s m)
prob-delta2-1-seq (prob-delta2-1-s m))
(* (prob-m-s m) (map-sum [1 patlen]
(fn [k] (* k (skip-s m k prob-delta2-seq prob-delta2-1-seq)))))))))
并且为了进行比较,还额外计算了简化 BM 算法:
myprob.hy
(defn simplified-bm-skip [patlen p m k]
(simplify-probfn patlen p [prob-delta1])
(if (= k 1) (+ (prob-delta1-s m 1) (prob-delta1-worthless p m))
(prob-delta1-s m k)))
(defn sbm-rate [patlen p]
(simplify-probfn patlen p [prob-m simplified-bm-skip])
(/
(map-sum [0 (dec patlen)]
(fn [m] (* (cost m) (prob-m-s m))))
(map-sum [0 (dec patlen)]
(fn [m] (* (prob-m-s m) (map-sum [1 patlen]
(fn [k] (* k (simplified-bm-skip-s m k)))))))))
和 KMP 算法:
myprob.hy
(defn getone [&rest body &kwonly [default None]]
(try
(get #*body)
(except [_ [IndexError KeyError TypeError]]
default)))
(defn prob-pi [patlen p]
"prob-pi-s(m, l) = prob-pi-seq[m][l]"
(defseq prob-pi-seq [n]
(cond [(= n 0) []]
[(= n 1) [1]]
[(= n 2) [(- 1 p) p]]
[(> n 2) (lfor i (range n)
(+
(* (getone prob-pi-seq (dec n) i :default 0) (- 1 p))
(* (get prob-pi-seq (dec n) (dec i)))))]))
prob-pi-seq)
(defn skip [m k prob-pi-seq]
(if (= m 0) (if (= k 1) 1 0)
(get prob-pi-seq m (- m k))))
(defn at-least-1 [n]
(if (< n 1) 1 n))
(defn kmp-rate [patlen p]
(simplify-probfn patlen p [prob-m prob-pi])
(setv prob-pi-seq (prob-pi-s))
(/
(map-sum [0 (dec patlen)]
(fn [m] (* (cost m) (prob-m-s m))))
(map-sum [0 (dec patlen)]
(fn [m] (* (prob-m-s m) (map-sum [1 (at-least-1 (dec m))]
(fn [k] (* k (skip m k prob-pi-seq)))))))))
然后我们就可以通过 Python 上的 plotnine
图形包看一下计算的数据(并用高斯过程回归拟合曲线):
import pandas as pd
from plotnine import *
import hy
from my_prob import bm_rate, sbm_rate, kmp_rate
theme_update(text=element_text(family="SimHei"))
def plot(p, title, N=30):
model_range = range(1, N+1)
data = {'rate':[], 'alg':[], 'patlen':[]}
categories_list = [(bm_rate, 'BoyerMoore'),
(sbm_rate, 'S BoyerMoore'),
(kmp_rate, 'KMP'),
(lambda patlen, p: 1/patlen, '$\\frac{1}{patlen}$')]
for alg_fun, label in categories_list:
data['rate'].extend([alg_fun(patlen, p) for patlen in model_range])
data['alg'].extend([label for _ in model_range])
data['patlen'].extend(model_range)
df = pd.DataFrame(data)
return (ggplot(df, aes(x='patlen', y='rate', color='alg'))
+ geom_point()
+ geom_smooth(method='gpr')
+ labs(color='Algs', title=title, x='$patlen$', y='$\\frac{cost}{skip}$'))
plot(1/256, '$p= \\frac{1}{256}$')
:
观察这个图像,令人印象深刻的首先就是抬头的一条大兰线,几乎笔直地画出了算法性能的下限,不愧是 KMP 算法,
接着会发现 BoyerMoore 算法与简化版 BoyerMoore 算法高度重叠的这条红绿紫曲线,同时也是
这就是在一般字符集下随机文本搜索能达到的
另外此时可以绝大多数的字符跳转依靠
接着我们可以看一下在经典的小字符集,比如在 DNA {A, C, T, G} 碱基对序列中算法的性能表现(plot(1/4, '$p= \\frac{1}{4}$')
):
曲线出现了明显的分化,当然 KMP 还是一如既往地稳定(゚▽゚)/,如果此时在测试中监控一下一下
总结一下,通过概率模型的计算,一方面看到了在较大的字符集,比如日常搜索的过程中 BoyerMoore 系列算法的优越表现,其中主要依赖
引用¶
buildLast update and/or translate time of this article,Check the history
editFound smelly bugs? Translation outdated? Wanna contribute with us? Edit this Page on Github
peopleContributor of this article minghu6
translateTranslator of this article Visit the original article!
copyrightThe article is available under CC BY-SA 4.0 & SATA ; additional terms may apply.