广义后缀自动机

1).自动机的介绍
      首先我们先来介绍一下什么是自动机，有限状态自动机的功能是识别字符串，令一个自动机A，若他能识别字符串S，就记为A(S)=Ture,否则A(S)=False。
自动机由五个部分组成，alpha:字符集，state:状态集合，init:初始状态，end：结束状态集合，trans：状态转移函数。
      令trans(s,ch)表示当前状态是s，在读入字符ch之后，所到达的状态。如果trans(s,ch)这个转移不存在，为了方便，设其为null，同时null只能转移到null。null表示不存在的状态。同时另trans(s,str)表示当前状态是s，在读入str之后，所到达的状态。

trans(s,str):
Cur=s;
For i=0 to length(str)-1
Cur=trans(Cur,str[i]);

Trans(s,str)就是Cur。
      那么自动机A能识别的字符串就是所有使得trans(init,x)∈end的字符串x。令其为Reg(A)。从状态s开始能识别的字符串，就是所有使得的字符串x。令其为Reg(s)。
2).后缀自动机模型的建立过程
      后缀自动机的定义：给定字符串S，S的后缀自动机suffix automaton(以后简记为SAM)是一个能够识别S的所有后缀的自动机。即SAM(x)=Ture,当且仅当x是S的后缀，后面会说明，后缀自动机也能用来识别所有子串。
      最简单的实现:考虑字符串”aabbabd”,我们可以将该字符串的所有后缀插入一个Trie(字典树)中，就像下图那样。那么初始状态就是根，状态转移函数就是这颗树的边，结束状态集合就是所有的叶子。注意到这个结构对于长度为N的串，会有O(N2)的节点。节点数量太多了，怎么办？

      最简状态自后缀自动机：顾名思义，就是状态数最少的后缀自动机，在后面可以证明它的大小是线性的，我们先介绍一些性质。假如我们得到了这个最简状态后缀自动机SAM，我们令ST(str)表示trans(init,str)。就是初始状态开始读入字符串str之后，所能到达的状态。
      分析：令母串为S，他的后缀的集合为Suf，他的连续子串的集合为Fac，从位置a开始的后缀为Suffix(a)，S[l,r)表示S中[l,r)这个区间构成的子串。下标从0开始，对于一个字符串s，如果它不属于Fac，那么ST(s)=null。因为s后面加上任何字符串都不可能是s的后缀了，没有理由浪费空间。同时如果字符串s属于Fac，那么ST(s)≠null。因为s既然是S的子串，就可以在后面加上一些字符使得其变成s的后缀，我们既然要识别所有的后缀，就不能放过这种可能性。我们不能对每个s∈Fac都建立一个状态，因为Fac的大小是O(N2)的。我们考虑ST(a)能识别哪些字符串，即Reg(ST(a))。字符串x能被自动机识别，当且仅当x∈Suf。ST(a)能够识别字符串x，当且仅当ax∈Suf。因为我们已经读入了字符串a了。也就是说ax是s的后缀，那么x也是s的后缀，Reg(ST(a))是一些后缀集合。对于一个状态s，我们唯一关心的是Reg(s)。
      如果a在s中的[l,r)位置出现，那么他就能识别s从r开始的后缀。例子：那么如果a在s中的出现位置集合是{[l1,r1),[l2,r2),...,[ln,rn)}那么Reg(ST(a))就是{Suffix(r1),Suffix(r2),...,Suffix(r3)}。令Right(a)={r1,r2,...,rn}那么Reg(ST(a))就完全由Right(a)决定。
      对于两个子串a,b∈Fac如果Right(a)=Right(b),那么ST(a)=ST(b)。所以一个状态s，由所有Right集合是Right(s)的字符串组成。不妨令r∈Right(s),那么只要给定一个长度就可以确定子串了。考虑对于一个Right集合，容易证明如果长度l,r合适，那么长度l≤m≤r的m也一定合适。所以合适长度必然是一个区间，不妨令s的区间是[min(s),max(s)]。
      关于子串的性质：由于每个子串都必然包含在SAM的某个状态里，那么一个字符串a是S的子串，当且仅当，ST(s)≠null,那么我们就可以用SAM来解决子串判定问题，同时也可以求出这个子串出现个数，就是所在状态Right集合的大小。在一个状态中直接保存Right集合会消耗过多的空间，我们可以发现状态Right就是它Parent树中所有孩子Right集合的并集，进一步的话，就是Parent树中它所有后代中叶子节点的Right集合的并集。那么如果按dfs序排序，一个状态的Right集合就是一段连续的区间中的叶子节点的Right集合的并集，那么我们也就可以快速求出一个子串的所有出现的位置了。树的dfs序列：所有子树中节点组成一个区间。
      线性构造算法：我们的构造算法是Online(在线)的，也就是从左到右逐个添加字符串中的字符，依次构造SAM，这个算法实现相比后缀树来说要简单很多，尽管可能不是非常好理解。让我们先回顾一下性质：
      定义和性质的回顾：状态s,转移trans,初始状态init，结束状态集合end，母串S,S的后缀自动机SAM(Suffix Automaton的缩写)，Right(str)表示str在母串S中所有出现的结束位置集合。

      一个状态s表示的所有子串Right集合相同，为Right(s)。
      Parent(s)表示使得Right(s)是Right(x)的真子集，并且是Right(x)的大小最小的状态x。
      Parent函数可以表示一个树形结构，叫他Parent树。
      一个Right集合和一个长度，定义了一个子串。
      对于状态s，使得Right(s)合法的子串长度是一个区间，为[Min(s),Max(s)]。
      Max(Parent(s)) = Min(s)-1。
      SMA的状态数量和边的数量，都是O(N)的。
      trans(s,ch)==null表示从s出发没有标号为ch的边。
      考虑一个状态s,他的Right(s)={r1,r2,...,rn},假如有一条s→t标号为c的边，考虑r的Right集合，由于多了一个字符，s的Right集合中，只有S[ri]=c的符合要求。那么t的Right集合就是{ri+1|S[ri]=c }。
      那么如果s出发有标号为x的边，那么Parent(s)出发必然也有。
      同时令f=Parent(s)，Right(trans(s,c))∈Right(trans(f,c))。
      有一个很显然的推论是Max(t)>Max(s)。
      我们每次添加一个字符，并且更新当前的SAM使得他成为包含这个新字符的SAM令当前的字符串为T,新字符为x,令T的长度为L，SAM(T)→SAM(Tx)
      那么我们新增加了一些子串，他们都是串Tx的后缀，Tx的后缀，就是T的后面添加一个x，那么我们考虑所有表示T的后缀(也就是Right集合中包含L)的节点v1,v2,v3,...
由于必然存在一个Right(p)={L}的节点p(ST(T))。那么v1,v2,...,vk,由于Right集合都含有L，那么他们在Parent树中必然全是p的祖先。可以使用Parent函数得到他们，同时我们添加一个字符x后，另np表示ST(Tx)，则Right(np)={L+1}。不妨让他们的后代到祖先排位v1=p,v2,...,vk=root，考虑其中一个v的Right集合={r1,r2,...,rn=L}，那么在他的后面添加一个新字符x的话，形成新的状态nv的话，只有S[ri]=x的ri那些是符合要求的。
      同时在之前我们知道，如果从v出发没有标号为x的边(我么你先不看rn),那v的Right集合内就没有满足这个要求的ri，
那么由于v1,v2,v3,...的Right集合逐渐扩大，如果vi出发有标号为x的边，那么vi+1出发也肯定有。
对于出发没有标号为x的边的v，他的Right集合内只有rn是满足要求的，所以根据之前提到的转移规则，让它连一条到np标号为x的边。
令vp为v1,v2,...,vk中第一有标号为x的变得状态。
考虑vp的Right集合={r1,r2,...,rn},令trans(vp,x)=q,那么q的Right集合就是{ri+1|S[ri]=x}的集合(注意到这是更新之前的情况，所以rn是不算的)。
注意：我们不一定能直接在q的Right集合中插入L+1。
接下来考虑节点nq，在转移的过程中，结束位置L+1是不起作用的，所以trans(nq)就跟原来的trans(q)是一样的，拷贝即可。
接下来，如果已经存在节点nq，我们还得处理。
回忆：v1,v2,...,vk 是多有Right集合包含{L}的节点按后代到祖先排序，其中vp 是第一个由标号为x的边的祖先。x是这轮新加入的字符。
由于vp,...,vk 都有标号为x的边，并且到达的点的Right集合，随着出发点Right集合的变大，也会变大，那么只有一段vp,...,ve ,通过标号为x的边，原来是到节点q的。(q=Trans(vp,x))。
那么由于在这里q节点已经被替换成了nq，我们只要把vp,...,ve 的Trans(*,x)设为nq即可。
每个阶段：回顾：
令当前串为T,加入字符为x.
令p=ST(T)，Right（p）={Length(T)}的字点
新建np=ST(Tx),Right(np)={Length(T)+1}的节点。
对p的所有没有标号x的边的祖先v,trans(x,v)=np
找到p的第一个祖先vp ,它有标号x的边，如果没有这样的vp，那么
Parent(p)=root,结束该阶段。
令q=trans(vp,x)，
若Max(q)=Max(vp)+1令Parent=(np)=q,结束该阶段。
否则新建节点nq，trans(nq,*)=trans(q,*)
Parent(nq) = Parent(q) (先前的)
Parent(q) = nq
Parent(np)=nq
对所有trans(v,x) == q的p的祖先v，trans(v,x) 改成nq。
图解：（来自：http://hi.baidu.com/myidea/item/142c5cd45901a51820e25039）
我们要将字符串”aabbab”建立成后缀自动机模型。
1.起始只有一个空节点root。

2.现在添加字符’a’,即字符串”a”的后缀自动机为：

3.再添加字符’a’，即字符串”aa”的后缀自动机为：

4.再添加字符’b’,即字符串”aab”的后缀自动机为：




5.再添加字符’b’，即字符串”aabb”的后缀自动机为：

(错误！)
所有结束点的parent树中点都是结束态，所以结点3是结束态，所以此图的”ab”也是后缀，但”ab”并不是”aabb”的后缀，我们还需要一个点辅助，所以我们要这么建立：

(正确！)
找到0(S)时,发现转移"b"已经有节点占了,所以新建节点5,将3号所有信息Copy(包括Parent),然后更新len值,就是node[5]->len=node[5]->parent->len+1,所以5号节点可以代表后缀空串(0号代表的串)+字符"b"=后缀"b",节点3成了中间态,所以将节点为原接受态的节点指向3的转移改为指向5,这时,我们发现指向3的原接受态节点一定是当前节点0(S)及当前未访问的原接受态节点,所以可以直接沿着Parent往上更新，然后节点5的Parent及祖先加入了现在的接受态，再次重申一点:一个节点及其父辈的代表的串有相同的后缀,且代表串长度递减,由于5号节点是接受态,所以他的父辈也是接受态,同时反过来也一样,与任意接受态拥有相同后缀的长度小于当前节点的未访问节点一定是当前节点的父辈,如与5号节点有相同后缀的长度小于5号节点的未访问的节点一定是5号的父辈,一定可以作为接受态。因此为了维护这个性质,我们应该将3号节点的父亲重定义为5，到这里基本上应该明白了。
再分别插入’a’,’b’的后缀自动机为:

广义后缀自动机：

      传统后缀自动机是解决单个主串的匹配问题，广义后缀自动机可以用来解决多个主串的匹配问题。

      如何将多个主串构建成广义后缀自动机？先将一个主串建立成后缀自动机，让后将重置last，令last=root，下一个字符串再从头节点开始建立，下一状态如果不存在，则以后缀自动机的规则进行建立新节点。

      如果下一状态已经建立，我们直接转移到该状态即可，既然到达该状态，说明已经匹配成功的字符串的所有后缀都在该状态及该状态的父节点，父节点的父节点...直到root，所以我们需要对该状态以及他的父节点，他的父节点的父节点。。。直到root进行内容更新，更新的内容当然因题目而异，如果求某个字符串出现的个数，就cnt++，如果求某个字符串出现的位置，就将位置存在结点的一个数组里。

举例如下：(给n个长度为100的字符串，有m个长度为K的目标字符串，问每个目标字符串在主串中哪些位置出现过。)

      通用性很强，如果目标字符串是任意长度，就把此代码所有用到K的地方删去；如果求目标字符串的个数，就在结构体中加一个计数cnt，然后将下方xx.push_back()改成cnt++；如果主串长度任意，可以用if(last==root&&weizhi！=0)来记录行数。其他用处于后缀自动机相同，网上也有好多例子与题解，在此不多叙述~。~

      优点：再也不用将多个主串连起来建立后缀自动机啦~！空间复杂度大大降低~代码难度不高，就是理解有点复杂。

#include <cstdio>
#include <cstring>
#include <iostream>
#include<vector>
using namespace std;
int weizhi,n,K,m;/*全局变量wezhi方便传递此时坐标，K代表K值，m代表有m个k-mer index*/
int f[30];/*分别将对A,C,G,T映射到0，1，2，3，可以大量减少单个结点的开销*/
struct Complex{/*定义节点*/
    Complex* tranc[4],*father;/*通过father和tranc来链接Parent树,tranc进行对下一个状态询问*/
    int len;/*代表次状态Right的MAX*/
    vector<int> ans;/*存储答案*/
}*last,*root,none,*nil=&none;

Complex *NewComplex(int _)/*申请新节点*/
{
    Complex *C=new Complex;
    C->len = _;
    fill(C->tranc,C->tranc + 4,nil);
    C->father = nil;
    return C;
}

void init()/*初始化*/
{
    f['A'-'A']=0,f['C'-'A']=1,f['G'-'A']=2,f['T'-'A']=3;/*对转移函数f的初始化*/
    root = last = NewComplex(0);/*建立节点root,初始化last*/
}

inline void Add(int c)/*对广义SAM中此时last后面添加状态c*/
{
    Complex *np,*p = last;
    if(last->tranc[f[c]]==nil){
        np= NewComplex(last->len + 1);
        for(; p != nil && p->tranc[f[c]] == nil; p = p->father)  p->tranc[f[c]] = np;
        if(p == nil)    np->father = root;
        else {
            Complex *q = p->tranc[f[c]];
            if(q->len == p->len + 1)  np->father = q;
            else {
                Complex *nq = NewComplex(p->len + 1);
                memcpy(nq->tranc,q->tranc,sizeof(nq->tranc));
                nq->father = q->father;
                np->father = q->father = nq;
                if(nq->len>=K&&nq->father->len<K){
                    nq->ans=q->ans;
                }
                for(;p!=nil&&p->tranc[f[c]]==q;p=p->father)    p->tranc[f[c]] = nq;
            }
        }
    }
    else{
        np=last->tranc[f[c]];
    }
    last = np;
    for(; np != root; np = np->father){
        if(np->len<K) break;
        if(np->father->len<K){
            np->ans.push_back(weizhi);
        }
    }
}

inline int Ask(char s[])/*插入新字符*/
{
    Complex *now = root;
    int length = strlen(s);
    for(int i = 0; i < length; ++i)
        if(now->tranc[f[s[i] - 'A']] != nil)
            now = now->tranc[f[s[i] - 'A']];
        else    return 0;
    int p=0;
    for(int i=0;i<now->ans.size();i++){
        if(now->ans[i]%100>=length-1){
            printf("NO.%d:\tX:%d\t,Y:%d\n",p+1,now->ans[i]/100+1,now->ans[i]%100+2-K);/*此处输出anser，请在此处修改输出格式。*/
        }
        else
            p--;
        p++;
    }
    return now->ans.size();
}
void scan_DNA(){
    for(int j=0;j<n;j++){
        for(int i=0;i<100;i++){
            char c;
            if(c=getchar(),c==EOF) break;
            if(c>='A'&&c<'Z'){
                Add(c - 'A');weizhi++;
            }
            else{
                i--;
                if(last!=root&&weizhi%100){
                    weizhi+=100-weizhi%100;
                }
                last=root;
            }
        }
    }

}
void scan_k_mer_print(){

    //printf("%lld",(long long)sizeof(Complex)); /*结点消耗内存测试*/

}
int k(){
    int k;
    printf("please input n and K:\n");
    scanf("%d%d",&n,&k);
    return k;
}
void scan_m_s(){
    fclose(stdin);
    freopen("CON","r",stdin);
    printf("Please input m:\n");
    scanf("%d",&m);
    char ss[105];
    for(int i = 0; i < m; ++i) {
        printf("Please input the %dth lenght is %d k-mer index:\n",i+1,K);
        scanf("%s",ss);
	printf("%s",ss);
        int last_ans = Ask(ss);
        printf("ALL:%d\n",last_ans);
    }
}

void autoscan_m_s(){
    scanf("%d",&m);
    char ss[105];
    for(int i = 0; i < m; ++i) {
        scanf("%s",ss);
	printf("%s",ss);
        int last_ans = Ask(ss);
        printf("ALL:%d\n",last_ans);
    }
}
int main(){
    K=k();/*读入K值*/
    freopen("f1out.txt","r",stdin);/*测试数据文件*/
	/*数据输出文件*/
    init();/*初始化*/
    scan_DNA();/*读入DNA序列并建立广义后缀自动机模型*/
    printf("The establishment of success!\n");
    //freopen("out_2.txt","w",stdout);/*想要将答案输出文件请取消此行注释*/
    /*想要控制台输入k-mer index请选择scan_s();并且注释掉autoscan_m_s();*/
    /*想要文件输入k-mer index，请选择autoscan_m_s();并且注释掉scan_m_s();*/
    scan_m_s();/*手动输入*/
    //autoscan_m_s();/*文件读入，注意:请保证文件格式正确!*/
    return 0;
}

压缩后缀自动机：

      答题思路就是将没有支链的一小段节点压缩成一个节点，这个压缩节点的父节点的next应记录此压缩节点中的第一个结点的状态。

      图解：

      因为一直没有写，实现就交给你们啦~！

谢谢CLJ大神的PPT指导！谢谢小彭玉、YYN、叉姐的指导！谢谢建模让我深入思考！