序列比对（26）精准匹配之KMP算法、Trie树以及AC自动机

一只羊

发布于 2019-09-27 18:41:53

1K0

发布于 2019-09-27 18:41:53

文章被收录于专栏：生信了

前文已经介绍过KMP算法和Trie树，本文将在此基础上介绍AC自动机。

之前的序列比对文章大都在利用动态规划算法解决字符串的非精准匹配（允许错配、插入和缺失），比如全局比对和局部比对问题。当然，后来我们还介绍了模序发现和中间字符串问题，并初次学习了如何用分支定界法解决这一类问题。

既然有非精准匹配问题，就有精准匹配问题。所谓精准匹配，就是两个字符串在比对时，不允许错配、插入和缺失。其实，笔者之前已经有过介绍：

KMP算法

《算法（四）（转载）KMP算法》一文介绍了KMP算法，该算法常用来解决如何在一个长字符串（模板串）中寻找一个短字符串（模式串）的问题。其特点就是要预先对模式串进行处理，对其构建一个next数组，该数组中存放着模式串每一个位置所对应的跳转指针，这样做的效果是模板串的指针不用回溯，一直在递增，从而增加效率。

KMP算法的完整代码：

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MAXLEN 1000

void getNext(const char *pat, int *next, const int n) {
/* assume n >= 2 */
    int i, j;
    for (next[0] = -1, next[1] = 0, i = 1, j = 0; i < n - 1;) {
        while (j != -1 && pat[i] != pat[j])
            j = next[j];
        next[++i] = ++j;
    }
}

int KMP(const char *temp, const int m, const char *pat, const int n) {
    int i, j;
    int *next = (int*) malloc(sizeof(int) * n);
    getNext(pat, next, n);
    for (i = 0, j = 0; i < m && j < n;) {
        if (j == -1 || temp[i] == pat[j]) {
            i++;
            j++;
        } else {
            j = next[j];
        }
    }
    free(next);
    if (j < n)
        return -1;
    return i - j;
}

int main(void) {
    char t[MAXLEN], p[MAXLEN];
    int index;
    printf("the temp str: ");
    scanf("%s", t);
    printf("the pattern str: ");
    scanf("%s", p);
    if ((index = KMP(t, strlen(t), p, strlen(p))) == -1)
        printf("Not found!\n");
    else
        printf("Index is %d.\n", index);
    return 0;
}

其效果如下：

Trie树（字典树）

《算法（五）字典树算法快速查找单词前缀》一文介绍了Trie树，该算法常用来解决如何在多个模板串中寻找一个模式串的问题。其特点就是将多个模板串装填进一个树结构中，虽然使用了较多的空间，但是查找效率得到了提升。

Trie树的完整代码：

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define N 26
#define MAXLEN 1000

struct Node {
    struct Node *a[N];
    int is_end;
};
typedef struct Node *tnode;

tnode init() {
    tnode p;
    int i;
    if ((p = (tnode) malloc(sizeof(struct Node))) == NULL) {
        fputs("Error: out of space!\n", stderr);
        exit(1);
    }
    for (i = 0; i < N; i++)
        p->a[i] = NULL;
    p->is_end = 0;
    return p;
}

void delete(tnode p) {
    int i;
    for (i = 0; i < N; i++)
        if (p->a[i] != NULL)
            delete(p->a[i]);
    free(p);
}

void insert(const char *s, tnode p) {
    int i;
    while (*s != '\n' && *s != '\0') {
        i = *s - 'a';
        if (p->a[i] == NULL) {
            p->a[i] = init();
        }
        p = p->a[i];
        s++;
    }
    p->is_end = 1;
}

int find(const char *s, tnode p) {
    int i;
    while (*s != '\n' && *s != '\0') {
        i = *s - 'a';
        if (p->a[i] == NULL)
            return 0;
        p = p->a[i];
        s++;
    }
    if (p->is_end)
        return 1;
    return 0;
}

void strLower(char *s) {
    while (*s != '\0') {
        if (*s >= 'A' && *s <= 'Z') {
            *s += 32;
        }
        s++;
    }
}

int main(int argc, char *argv[]) {
    FILE *fpp, *fpw;
    char buf[MAXLEN];
    tnode root = init();
    if (argc < 3) {
        fputs("Usage: trie.exe <pattern.txt> <words.txt>\n", stderr);
        exit(2);
    }
    if ((fpw = fopen(argv[2], "r")) == NULL) {
        fprintf(stderr, "Error: cannot open %s\n", argv[2]);
        exit(4);
    }
    while (fgets(buf, MAXLEN, fpw) != NULL) {
        strLower(buf);
        insert(buf, root);
    }
    if (! feof(fpw)) {
        fprintf(stderr, "Error: error while reading %s\n", argv[2]);
        exit(8);
    }
    fclose(fpw);
    if ((fpp = fopen(argv[1], "r")) == NULL) {
        fprintf(stderr, "Error: cannot open %s\n", argv[1]);
        exit(4);
    }
    while (fgets(buf, MAXLEN, fpp) != NULL) {
        strLower(buf);
        if (find(buf, root)) 
            printf("%s", buf);
    }
    if (! feof(fpp)) {
        fprintf(stderr, "Error: error while reading %s\n", argv[1]);
        exit(8);
    }    
    fclose(fpp);
    delete(root);
}

其效果如下：

其中，pat.txt文件内容（里面是多个模式串）如下：

how are you

words.txt文件内容（里面是多个模板串）如下：

lily we are friends

AC自动机

本文要介绍的AC自动机算法，是用来解决如何有效地在一个模板串中查找多个模式串的问题。该算法是建立在KMP算法和Trie树基础上的。其特点是：

预先将多个模式串装填进Trie树中。
在进行查找（匹配）之前，对Trie树中的每个有效节点构建fail指针，该指针的作用类似KMP算法中next数组的作用：如果到该节点时匹配失败，就可以利用fail指针跳转到下一个可利用节点继续进行匹配，从而避免了模板串的回溯而提升查找效率。

如同KMP算法的next数组充分利用了模式串的内部信息，AC自动机中的fail指针也是充分利用了多个模式串的内部信息，每一次跳转都是跳到“最大可利用后缀子字符串”的节点。

笔者在学习该算法时，主要参考了这一篇文章（https://www.jianshu.com/p/641142c6d477），看里面的图非常有助于理解AC自动机原理。

关于AC自动机的实现有几点说明：

在对Trie树里所有有效节点构建fail指针的过程中，采用了BFS（广度优先搜索）的策略，具体的实现是利用了队列这种数据结构。采用BFS的原因是Trie树中某一层节点的fail指针指向的节点一定是来自于上面的层（可能是上面第一层，也有可能是上面第二层等等，甚至是根节点）。
无论是Trie树还是队列，都是利用一种链表实现的，不得不感叹，对链表进行操作真是很容易出错，经常要填坑。

AC自动机的完整代码如下：

#include <stdio.h>
#include <stdlib.h>
#define ALPHA_NUM 26
#define QUEUE_EMPTY NULL
#define BUFSIZE 64
#define MAXTEMP 10000

/* Data structure for trie */
struct TNode {
    struct TNode *a[ALPHA_NUM];
    struct TNode *fail;
    struct TNode *parent;
    int c;
    int is_end;
};
typedef struct TNode *tnode;

tnode initTrie(int c, tnode parent) {
    tnode p;
    int i;
    if ((p = (tnode) malloc(sizeof(struct TNode))) == NULL) {
        fputs("Error: out of space!\n", stderr);
        exit(1);
    }
    for (i = 0; i < ALPHA_NUM; i++)
        p->a[i] = NULL;
    p->fail = NULL;
    p->parent = parent;
    p->c = c;
    p->is_end = 0;
    return p;
}

void deleteTrie(tnode p) {
    int i;
    for (i = 0; i < ALPHA_NUM; i++)
        if (p->a[i] != NULL)
            deleteTrie(p->a[i]);
    free(p);
}

void insertTrie(const char *s, tnode p) {
    int i;
    while (*s != '\n' && *s != '\0') {
        i = *s - 'a';
        if (p->a[i] == NULL) {
            p->a[i] = initTrie(*s, p);
        }
        p = p->a[i];
        s++;
    }
    p->is_end = 1;
}

/*
void printTrie(tnode p) {
    int i;
    printf("%c ", p->c);
    for (i = 0; i < ALPHA_NUM; i++) {
        if (p->a[i] != NULL)
            printTrie(p->a[i]);
    }
}
*/

void backtrack(char *s, const int n, tnode p) {
    int i, j, c;
    for (i = 0; i < n - 1 && p != NULL && p->parent != NULL; i++) {
        s[i] = p->c;
        p = p->parent;
    }
    if (i >= n - 1) {
        fputs("Error: pattern word is too long!\n", stderr);
        exit(2);
    }
    for (j = 0; j < i / 2; j++) {
        c = s[j];
        s[j] = s[i - 1 - j];
        s[i - 1 - j] = c;
    }
    s[i] = '\0';
}

/* Data structure for queue */
typedef tnode QNodeType;

struct QNode {
    QNodeType value;
    struct QNode *next;
};
typedef struct QNode *qnode;

qnode initQueue(QNodeType v) {
    qnode p;
    if ((p = (qnode) malloc(sizeof(struct QNode))) == NULL) {
        fputs("Error: out of space!\n", stderr);
        exit(1);
    }
    p->value = v;
    p->next = NULL;
    return p;
}

void addQueue(QNodeType v, qnode head) {
    qnode p;
    for (p = head; p->next != NULL; p = p->next)
        ;
    p->next = initQueue(v);
}

int isEmptyQueue(qnode head) {
    return head->next == NULL ? 1 : 0;
}

QNodeType popQueue(qnode head) {
    qnode p;
    QNodeType q;
    if (isEmptyQueue(head))
        return QUEUE_EMPTY;
    else {
        p = head->next;
        head->next = p->next;
        q = p->value;
        free(p);
        return q;
    }
}

void destroyQueue(qnode head) {
    while (popQueue(head) != QUEUE_EMPTY)
        ;
    free(head);
}

void strLower(char *s) {
    while (*s != '\0') {
        if (*s >= 'A' && *s <= 'Z') {
            *s += 32;
        }
        s++;
    }
}

void getFailPointer(tnode root) {
    int i;
    qnode head = initQueue(NULL);
    tnode this, tmp;
    addQueue(root, head);
    while (! isEmptyQueue(head)) {
        this = popQueue(head);
        for (i = 0; i < ALPHA_NUM; i++) {
            if (this->a[i] != NULL) {
                addQueue(this->a[i], head);
                tmp = this;
                while (tmp->fail != NULL) {
                    if (tmp->fail->a[i] != NULL) {
                        this->a[i]->fail = tmp->fail->a[i];
                        break;
                    } else {
                        tmp = tmp->fail;
                    }
                }
                if (tmp->fail == NULL)
                    this->a[i]->fail = root;
            }
        }       
    }
    destroyQueue(head);
}

void ACauto(char *buf, const int n, char *temp, tnode root) {
    char *s = temp;
    tnode p, tmp;
    int i;
    for (p = root; *s != '\0'; s++) {
        i = *s - 'a';
        while (p->a[i] == NULL && p->fail != NULL)  /* 这里必须先判断p->a[i]，否则有可能错过root->a[i] != NULL的情况 */
            p = p->fail;
        if ((p = p->a[i]) == NULL)
            p = root;
        for (tmp = p; tmp->fail != NULL; tmp = tmp->fail)
            if (tmp->is_end) {
                backtrack(buf, BUFSIZE, tmp);
                printf("%s\n", buf);
            }
    }
}

int main(int argc, char *argv[]) {
    FILE *fpp, *fpt;
    char buf[BUFSIZE];
    char temp[MAXTEMP];
    tnode root = initTrie('\0', NULL);
    int n;
    if (argc < 3) {
        fputs("Usage: trie.exe <pattern.txt> <template.txt>\n", stderr);
        exit(2);
    }
    if ((fpp = fopen(argv[1], "r")) == NULL) {
        fprintf(stderr, "Error: cannot open %s\n", argv[1]);
        exit(4);
    }
    while (fgets(buf, BUFSIZE, fpp) != NULL) {
        strLower(buf);
        insertTrie(buf, root);
    }
    if (! feof(fpp)) {
        fprintf(stderr, "Error: error while reading %s\n", argv[1]);
        exit(8);
    }
    fclose(fpp);
    if ((fpt = fopen(argv[2], "r")) == NULL) {
        fprintf(stderr, "Error: cannot open %s\n", argv[2]);
        exit(4);
    }
    n = fread(temp, 1, MAXTEMP - 1, fpt);
    temp[n] = '\0';
    strLower(temp);
    if (! feof(fpt)) {
        fprintf(stderr, "Error: error while reading %s\n", argv[2]);
        exit(8);
    }    
    fclose(fpt);
    getFailPointer(root);
    ACauto(buf, BUFSIZE, temp, root);
    deleteTrie(root);
}

其效果如下：