用GST解决Longest Common Substring Problem
最长公共子字符串问题,是指给定N个字符串,在每个字符串中找最长的公共子字符串。该问题最常出现在DNA测序的计算中。给定一些列DNA片段,找到重复的地方将其对接起来。解决这个问题需要我们前面提到的后缀树Suffix Tree(ST)和Generalised Suffix Tree(GST)。

创建GST的方法和创建ST的方法类似,首先我们为每一个字符串Sn加上一个终结符号,根据这个终结符号,可以在GST中区别出该后缀来自哪个字符串。然后我们先将每个添加终结符后的字符串Sn’创建一个ST,并且将Sn’添加到GST中去。
如图得到的GST中,可以直观的看到我们要的结果。但程序应该如何判断呢?
首先我们从每个ST(Sn)的每个后缀开始迭代。在GST的底部开始向根节点查找。每遇到一个新的支节点,则将其标记做上属于字符串Sn的标记。如果遇到已经属于本字符串的标记,开始迭代下一个字符串。
然后,我们将获得了所有Sn个字串标记的支节点。取的其代表的公共字符串L(Sn)。如图,得到了A, CTG, TG和G。很显然,其中最长的是CTG。
这样就得到了两个字符串的公共子串CTG。
本文参考了[Inoke Lee 2005] Linear Time Algorithm for the Longest Common Repeat Problem。事实上论文中论述的是最长公共重复子串问题,比本问题还要复杂。由于篇幅限制不在这里做详细介绍,有兴趣的读者可以Google下这篇论文。本文可以看作是GST应用的简单例子。

构建Suffix Tree

先了解什么是字符串的后缀:字符串任意一个字符到尾字符形成自字符串。比如banana,所有a,na,ana,nana,anana,banana都算作该字符串的后缀。同理可以知道前缀的定义。 有了前
后缀的定义,再来看后缀树的定义.
一个长度为n的字符串S的后缀树是一种这样的树:
从根节点到叶节点的路径同字符串的后缀是一一对应关系;
叶节点所拼的是非空字符串;
所有的内节点(除root外)最少含有两个孩子。
后缀树的构建:
一个字符串S的后缀树的构建的方法是:
1)从字符串S最短的前缀P1(即{S1})按照由短到长的顺序开始,
2)从根节点开始查找能够表示该Pi前缀的每个后缀Psi的路径;
3)如果存在完整的路径,到5);
4)如果不存在完整路径,则在树中够在出一条能够表示该后缀的路径;
5)处理下一个后缀Psi+1
举例来说,我们用banana$构建一个后缀树。注意字符串是banana+一个字符$。至于为什么多了一个$字符下一篇文章再做介绍。这里先将$当作和a,b,n一样的字符。
1)首先处理b
从根节点开始,因为找不到能够表示b的路径。我们插入一个子节点,并将其路径表示为b;
|(1:b)|leaf
由于b只有一个后缀,我们进行下一步。
2)处理ba
|(1:ba)|leaf
tree:|
|(2:a)|leaf
ba的后缀有a,和ba,目前树中不存在表示a的路径,又根目录有子节点。所以创建一个新的子节点并将路径表示为a。对于ba,因为存在一个路径表示为b的节点。因此在其后添加a即可。
3)直到banana都采用相同的算法:
|(1:banana)|leaf
tree:|
|(2:anana)|leaf
|
|(3:nana)|leaf
4)处理banana$
|(1:banana$)|leaf
tree:|
| | |(5:na$)|leaf
| |(3:na)|
| | |(7:$)|leaf
|(2:a)|
| |(7:$)|leaf
|
| |(5:na$)|leaf
|(3:na)|
| |(7:$)|leaf
|
|(7:$)|leaf
先处理$,直接在跟节点后插入。
再处理a$,由于树中已经有了anana,我们需要将这条路径拆分,他们享有公共的路径a,和两条分支$,和nana。
同理,对于na$,我们找到nana,将其拆分成na->$和na->na。
再处理ana$,我们查到有a->nana这条路径,因此,我们将其拆分成a->na->$和a->na->na
最后处理banana$,anana$, nana$由于已经有了banana,anana, nana在其后添加$即可。
这样,一个字符串banana$的后缀树就构造完成了。我们观察到。通过这个后缀树我们可以轻松的获得字符串中重复的子字符串:na和ana。其中ana最长,我们称之为最长重复子字符串。事实上,我们还可以用把多个字符串放到一个Suffix Tree中,构造出通用的Suffix Tree。利用这种树可以得出多个字符串中公共的子串。具体方法将在下一篇文章中介绍。
参考链接

Y Combinator
我们常常遇到函数的自变量取值和函数值相等的情况,比如y = x^2当,x取0时,y=0, x = 1是y = 1。我们称这样的x的取值为fixed point。我们把这个概念推广到一般的情况,函数变成higher order function,自变量变成first-oder value。我们可以称一个高阶函数f的不动点p为满足等式f(p)=p的函数p。而Fixed point combinator,则是生成不动点p的函数g:
p = g(f), f(p)=p 或者 f(g(f)) = g(f) ;; 参见维基百科
Y combinator是一种最常见的不动点组合算子。我们先看它的应用:
(define combined (lambda (n)
((lambda (f) (f f n))
(lambda (g n)
(if (= n 0)
1
(* n (g g (- n 1))))))))
这是典型的匿名递归函数。实现的是n的阶乘n!。看起来让人发蒙吧?我们可以把求阶乘的部分用F(n)代替:
(((lambda (f) F(f f)) ((lambda (f) F(f f)))
[code lang="list"]
我们把这样的段代码一般化,就得到了Y combinator的一般形式:
Y = λf·(λx·f (x x)) (λx·f (x x))
最后,我们可以得到一个通用的combinator:
[code lang="lisp"]
(define combinator (lambda (g) (
(lambda (f) (g (lambda (arg) ((f f) arg))))
(lambda (f) (g (lambda (arg) ((f f) arg)))) )))
(define factor (lambda (f) (lambda (x)
(if (= x 0)
1
(* x (f (- x 1))) ))))
>((combinator factor) 5)
120
>
Python的版本:
Y = lambda g: (lambda f: g(lambda arg: f(f)(arg))) (lambda f: g(lambda arg: f(f)(arg))) fac = lambda f: lambda n: (1 if n<2 else n*f(n-1)) print Y(fac)(7)

在Hadoop上进行朴素贝叶斯分类器的训练 Part 1
基础知识可以看:http://www.cnblogs.com/phinecos/archive/2008/10/21/1315948.html
只需要注意关于贝叶斯公式的推到有处笔误,应该改成:P(Ai|B) = P(B|Ai)P(Ai)/sigma_j(B|Aj)P(Aj)
问题的分解
不考虑各类平滑的技术。我们只需要计算:1)likelihood(P(cj) = M(C = cj)/ N 即训练文本中属于cj类别的文本数量/训练文本的总数量;2) likelihood(P(xj|cj) = N(Xi = xi, C = cj) / (N(C = cj) 即类别cj中包含属性xi训练文本数量/类别中的训练文本数量。
需要从训练集中提取的信息有:
1)属于每个类别的文本的数量
2)训练文本总数
3)类别cj中包含属性xi的训练文本的数量
4) 类别的个数
其中只有3)需要进行大量的查找和匹配,适合放在Hadoop上进行计算。
Map
先将训练文本拆分成tokenlist,然后按照为参数传递给Mapper, Mapper使用分类和每个token组成一个新的key,并且赋值为1。即该token和该分类同时出现的次数。
public static class Map extends MapReduceBase implements Mapper {
private final static IntWritable one = new IntWritable( 1 );
public void map( Text key ,ArrayWritable value,
OutputCollector output, Reporter reporter )
throws IOException {
Writable[] lines = value.get();
for ( int i = 0; i < lines.length; i++ ) {
output.collect( new Text( key.toString() + " " + ( ( Text )lines[ i ] ).toString() ), one );
}
}
}
Reduce
Reducer按照class+token为key进行计数,统计出该训练文本中token属于该类别的次数。
public static class Reduce extends MapReduceBase implements Reducer {
public void reduce( Text key, Iterator values, OutputCollector output, Reporter reporter ) throws IOException {
int sum = 0;
while ( values.hasNext() ) {
sum += values.next().get();
}
output.collect( key, new IntWritable( sum ) );
}
}
实测采用了lingpipe的1800份英文文本,在公司的单台机器本地上本地测试跑了3分多钟。感觉数据量不大的话训练的工作完全可以由单机完成。真正需要hadoop出马的是后续的分类工作。准备放在Part II进行介绍。
-------------------------------------------------------------------------
本文所介绍的贝叶斯和最大似然公式均来自网上,实际应用可能需要进行一定的变换和添加平滑的参数,才能在程序中进行准确的计算。参考请自行承担风险。

再谈尾递归
C/PASCAL出身的程序员刚接触函数式编程时都会对递归和尾递归的大量运用有些不适应。 那么能否找到一种通用的方法,将迭代算法转换为尾递归算法呢?个人觉得是可以的。以迭代为例,迭代可以认为是影射函数Fn(x),源集合和目标集合组成的三元组:
destination []
for x in source:
y = Fn(x)
destination.append(y)
我们观察每次迭代,都是在源集合中取出一个元素x,进行Fn(x)运算后,将Fn(x)的值插入到目标集合中去。于是for迭代可以转化成一下的形式:
def helper( Fn, source, destination): x = source.pop_back() destination.append(Fn(x)) return destination def iterate(Fn, source): destination = [] while true: if (source.emtpy()) return destination helper(Fn, source, destination)
我们观察到iterate函数里,引入了一个destination容器来存储中间运算结果。所以在尾递归的helper函数中,我们要使用三个参数:低阶函数fn(), source集合,destination集合。只要参照第二段代码,在循环处使用递归,就可以将其转化成尾递归的写法:
(define (foreach fn x) (define (helper fn pending done) (if (empty? pending) done (helper fn (butlast pending) (sentence (fn (last pending)) done)) )) (helper fn x '()) )
See?这就是所谓的线性递归.对于更复杂的递归方式,我们一样也可以参照这种思路生成进行转换.有的童鞋可能会想N个参数的迭代函数可以转化为N+1个helper参数的尾递归函数.其实不见得.helper的参数个数是由高介函数的定义所决定的。
First Class Data
SICP中由First Order Function引出了First Class Data.所谓First Class Data是指可以在程序中“没有使用限制”的数据。如何理解没有限制的数据,我们可以通过看几个例子。First Class Data可以被用来作为:
变量的值;
一个过程的参数;
一个过程的返回值;
集合的成员等等。
个人理解是一个语言的First Class Data可以比喻一个集合,把First Order Function比喻成集合上的操作。如果函数的值域也是First Class Data,那么就可以形成一个数据和函数的field。
一门语言允许哪些数据作为First Class Data可以衡量这们语言的动态性。

关于hashCode算法
Java中选择了这样的实现:
s[0]*31^(n-1) + s[1]*31^(n-2) + ... + s[n-1]
这个算法我在之前提到过。今天在地铁上看书想起来,当时对31的理解是错误的。31很容易让人想起32-1,32位计算机算术中的常见数值。其实是一种以X为底,字符串中的ascii值作为系数的多项式。为了防止化简造成的系数的冲突,所以X应该取质数。这里用的是31,还有建议用100003的。
Sn的再用2^64取模,就得到了int型的hashCode()。

尾递归vs迭代
在很多支持函数式特性语言中,经常鼓励使用尾递归。网上看到很多文章引用SICP中的例子作为支撑:
(define (factorial n)
(fact-iter 1 1 n))
(define (fact-iter product counter max-count)
(if (> counter max-count)
product
(fact-iter (* counter product)
(+ counter 1)
max-count)))
尾递归是否像描述的那样,可以用在任何场合来替代迭代语法呢?不是这样! 对于基于stack的语言,函数调用必然会带来额外的空间消耗,即使尾递归,空间增长也是O(n)的。解决办法是通过”trampine“这样的技术来为 尾递归做优化,尾递归函数调用自身时,先返回给trampine要调用的参数,由tarmpine再次调用自身。这样递归可以无限进行下去。理论上讲和迭 代的效率是一样的。
所以对于适合放心使用尾递归的语言,应当是在编译/解释时进行优化的语言。只有在了解相应语言的实现方法的前提下,才可以考虑是否放心的使用尾递归。
关于内核中的红黑树
内核文件linux/include/linux/rbtree.h中有这样的定义:
#define rb_parent(r) ((struct rb_node *)((r)->rb_parent_color & ~3))
static inline void rb_set_parent(struct rb_node *rb, struct rb_node *p)
{
rb->rb_parent_color = (rb->rb_parent_color & 3) | (unsigned long)p;
}
想了半天才发现作者用了一个十分tricky的技巧。我们看rb_node的定义:
struct rb_node
{
unsigned long rb_parent_color;
#define RB_RED 0
#define RB_BLACK 1
struct rb_node *rb_right;
struct rb_node *rb_left;
} __attribute__((aligned(sizeof(long))));
/* The alignment might seem pointless, but allegedly CRIS needs it */
我们知道rb_node占3个计算机字长的空间。如果分配了N个rb_node的连续空间。于是我们计算第k个rb_node地址的时候。其地址应当为: 数组的起始地址 + k*3。
那么,当我们能够保证这段内存起始地址的最后两个位是0(也就是~3的最后两位)时,就可以保证每一个rb_node的地址的最后两位都是0。这相当与对于所有用户管理这段内存的指针,指针的最后两位是空闲的。在内核的rbtree代码里面。这两位就被用来保存红黑树节点的颜色了(事实上还浪费了一位)。
因此里的rb_parent_color不是父亲的颜色,而是父节点的指针+颜色。只要在使用时将后两位置零,就获得了父节点的指针。
可惜这个技巧在现实中应用的有限。我们很少遇到申请一段较小的连续内存,又要十分节省的使用的情况。
关于外部FIFO的实现方法
最近有个项目的数据量较大,需要将原本在内存中的队列改写为外部存储。粗略的分析了一下,这个队列的要求如下。
1、First in first out
不管队列有多长,必须保证先进先出的特性;
2、Blocking
由于生产消费模型的需要,队列必须有阻塞能力;
3、外部存储
队列必须利用磁盘作为存储媒介;
4、缓存
必须拥有缓存机制来提高入队和出队操作的速度;
5、封装性;
实现必须保证界面同原来内存中的队列一致;
其中外部存储我想到两种方案,一种是类似larbin项目中的分段文件的方式。另一种就是使用sqlite。
前者的特点是可控性强,由于代码均有自己编写。另外一个特点是针对性强,对于队列的拆分和通过文件名寻址等等都由自己控制,采用什么样的参数也可以自己设置。此外利用Java的串行化技术,可以省去存储格式转换的过程。后者的特点是实现简单,单一文件存储,以及队列中的数据可以通过标准的工具查看和导出。缺点是存储的实现被掩盖,表中数据的访问也可能不是顺序的,性能不便优化。同时,数据的读出和写入都需要进行转换。
那么到底选择哪种呢?有没有更好的方法呢?今天晚上回家要好好想一下。