Thuật toán Optimal Mismatch

Đặc điểm chính

Ý tưởng

Trong một phép thử, mục tiêu là phát hiện sai lệch càng sớm càng tốt để loại bỏ vị trí hiện tại và trượt sang vị trí mới. Nếu ta so các vị trí mẫu theo đúng thứ tự trái-phải, có thể phải so nhiều ký tự "phổ biến" (dễ khớp) trước khi chạm tới ký tự thực sự phân biệt.

Optimal Mismatch của Sunday đảo ngược trực giác đó: hãy so ký tự hiếm trước. Nếu một ký tự của mẫu ít xuất hiện trong văn bản, thì tại một vị trí ngẫu nhiên xác suất nó khớp là thấp — nghĩa là khả năng cao ta sẽ tìm được sai lệch ngay lập tức và bỏ qua vị trí đó với chi phí chỉ một phép so sánh. Do đó thuật toán xếp các vị trí của mẫu theo tần suất ký tự tăng dần và so theo thứ tự đó.

Mô tả chi tiết

Thuật toán cần ba cấu trúc tiền xử lý.

1. Thứ tự so sánh (orderPattern). Với mỗi vị trí i của mẫu ta lập một cặp (loc = i, c = x[i]). Các cặp được sắp xếp theo tần suất của ký tự c trong văn bản tăng dần; khi hai ký tự có cùng tần suất, ta phá hoà bằng vị trí lớn hơn xếp trước. Kết quả là mảng pat[]: pat[0] là vị trí "đáng so nhất" (ký tự hiếm nhất), pat[m-1] là vị trí ít ưu tiên nhất.

2. Bảng dịch Quick Search (qsBc). Giống thuật toán Quick Search, bảng này cho biết dịch bao xa dựa vào ký tự văn bản ngay sau cửa sổ, tức y[j+m]. Với mỗi ký tự a: qsBc[a] = m - i với i là vị trí xuất hiện phải nhất của a trong mẫu, hoặc m + 1 nếu a không có trong mẫu.

3. Bảng adaptedGs. Đây là biến thể của quy tắc "good suffix" nhưng thích ứng với thứ tự so sánh đã đảo. adaptedGs[k] cho biết cần dịch cửa sổ bao nhiêu khi ta đã khớp thành công k vị trí đầu tiên (theo thứ tự pat) rồi mới gặp sai lệch (hoặc khớp toàn bộ khi k = m). Bảng được dựng bằng hàm matchShift, tính lượng dịch nhỏ nhất sao cho các vị trí đã kiểm tra vẫn nhất quán khi cửa sổ dời chỗ.

Trong pha tìm kiếm, tại vị trí cửa sổ j ta so lần lượt pat[0], pat[1], …: kiểm tra pat[i].c với y[j + pat[i].loc]. Nếu tất cả m vị trí khớp, ghi nhận xuất hiện tại j. Lượng dịch là giá trị lớn hơn giữa adaptedGs[i] (i là số vị trí đã khớp) và qsBc[y[j+m]].

Cài đặt bằng C

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define ASIZE 256
#define XSIZE 1024
#define MAX(a, b) ((a) > (b) ? (a) : (b))

typedef struct {
    int loc;
    unsigned char c;
} Pair;

static int freq[ASIZE];      /* tần suất ký tự trong văn bản */

/* Đếm tần suất ký tự của văn bản để định thứ tự so sánh */
void computeFreq(char *y, int n) {
    int i;
    memset(freq, 0, sizeof(freq));
    for (i = 0; i < n; ++i)
        ++freq[(unsigned char)y[i]];
}

/* So sánh: ký tự tần suất thấp xếp trước (so trước) */
int optimalPcmp(const void *a, const void *b) {
    const Pair *p1 = (const Pair *)a, *p2 = (const Pair *)b;
    int d = freq[p1->c] - freq[p2->c];
    return d != 0 ? d : (p2->loc - p1->loc);
}

void orderPattern(char *x, int m,
                  int (*pcmp)(const void *, const void *), Pair *pat) {
    int i;
    for (i = 0; i < m; ++i) {
        pat[i].loc = i;
        pat[i].c = (unsigned char)x[i];
    }
    qsort(pat, m, sizeof(Pair), pcmp);
}

/* Bảng dịch kiểu Quick Search dựa vào ký tự y[j+m] */
void preQsBc(char *x, int m, int qsBc[]) {
    int i;
    for (i = 0; i < ASIZE; ++i)
        qsBc[i] = m + 1;
    for (i = 0; i < m; ++i)
        qsBc[(unsigned char)x[i]] = m - i;
}

int matchShift(char *x, int m, int ploc, int lshift, Pair *pat) {
    int i, j;
    for (; lshift < m; ++lshift) {
        i = ploc;
        while (--i >= 0) {
            j = pat[i].loc - lshift;
            if (j < 0)
                continue;
            if (pat[i].c != (unsigned char)x[j])
                break;
        }
        if (i < 0)
            break;
    }
    return lshift;
}

void preAdaptedGs(char *x, int m, int adaptedGs[], Pair *pat) {
    int lshift, i, ploc;

    adaptedGs[0] = lshift = 1;
    for (ploc = 1; ploc <= m; ++ploc) {
        lshift = matchShift(x, m, ploc, lshift, pat);
        adaptedGs[ploc] = lshift;
    }
    for (ploc = 0; ploc < m; ++ploc) {
        lshift = adaptedGs[ploc];
        while (lshift < m) {
            i = pat[ploc].loc - lshift;
            if (i < 0 || pat[ploc].c != (unsigned char)x[i])
                break;
            ++lshift;
            lshift = matchShift(x, m, ploc, lshift, pat);
        }
        adaptedGs[ploc] = lshift;
    }
}

void search(char *x, int m, char *y, int n) {
    int i, j;
    int adaptedGs[XSIZE + 1], qsBc[ASIZE];
    Pair pat[XSIZE];

    /* Tiền xử lý */
    computeFreq(y, n);
    orderPattern(x, m, optimalPcmp, pat);
    preQsBc(x, m, qsBc);
    preAdaptedGs(x, m, adaptedGs, pat);

    /* Tìm kiếm */
    j = 0;
    while (j <= n - m) {
        i = 0;
        while (i < m && pat[i].c == (unsigned char)y[j + pat[i].loc])
            ++i;
        if (i >= m)
            printf("Tim thay tai vi tri %d\n", j);
        j += MAX(adaptedGs[i], qsBc[(unsigned char)y[j + m]]);
    }
}

Ví dụ minh họa

Xét mẫu x = "GCAGAGAG" (m = 8) trên văn bản y = "GCATCGCAGAGAGTATACAGTACG".

Trước hết đếm tần suất trong văn bản: A và G xuất hiện nhiều, còn C, T ít hơn. Khi sắp thứ tự so sánh, các vị trí mang ký tự hiếm (ví dụ vị trí 1 mang C) được đẩy lên đầu mảng pat, còn các vị trí mang A/G xuống cuối.

Vị trí cửa sổ j = 0 (y[0..7] = "GCATCGCA"): thay vì so x[0]=G trước, thuật toán so vị trí ưu tiên nhất — chẳng hạn vị trí 1 (C). Ở đây y[0+1] = C khớp, nên nó so tiếp các vị trí sau; đến vị trí mang G tại chỉ số 3, y[0+3] = T sai lệch. Cửa sổ dịch theo max(adaptedGs, qsBc[y[8]]).

Vị trí cửa sổ j = 5 (y[5..12] = "GCAGAGAG"): mọi vị trí trong pat đều khớp, thuật toán duyệt hết m vị trí và ghi nhận xuất hiện tại vị trí 5. Sau đó dịch tiếp bằng adaptedGs[m].

Nhờ so ký tự hiếm trước, ở đa số vị trí không khớp, thuật toán chỉ tốn một hoặc hai phép so sánh là đã loại được vị trí đó.

Tài liệu tham khảo