序列比對（25）編輯距離
- 2019 年 10 月 4 日
- 筆記
本文介紹兩個字元串的編輯距離並給出程式碼。
編輯距離
編輯距離的求解過程和全局比對是十分相似的（關於全局比對，可以參見前文《序列比對（一）全局比對Needleman-Wunsch演算法》），都需要全部符號參與比對，都允許插入、缺失和錯配。所以，編輯距離可以用動態規劃演算法求解，其迭代公式是：
效果如下：
編輯距離與最長公共子序列
在只允許插入和缺失而不允許錯配的情況下，兩個字元串的編輯距離可以通過最長公共子序列的長度（關於最長公共子序列，可以參看前文《序列比對（24）最長公共子序列》）間接算出來。
解編輯距離的程式碼
#include <stdio.h> #include <stdlib.h> #include <string.h> #define MAXSEQ 1000 #define GAP_CHAR '-' #define min(a, b) ((a) < (b) ? (a) : (b)) #define diff(a, b) ((a) == (b) ? 0 : 1) struct Unit { int W1; // 是否往上回溯一格 int W2; // 是否往左上回溯一格 int W3; // 是否往左回溯一格 int M; // 得分矩陣第(i, j)這個單元的分值，即序列s(1,...,i)與序列r(1,...,j)比對的最低得分 }; typedef struct Unit *pUnit; void strUpper(char *s); void printAlign(pUnit** a, const int i, const int j, char* s, char* r, char* saln, char* raln, int n); void align(char *s, char *r); int main() { char s[MAXSEQ]; char r[MAXSEQ]; printf("The 1st seq: "); scanf("%s", s); printf("The 2nd seq: "); scanf("%s", r); align(s, r); return 0; } void strUpper(char *s) { while (*s != ' ') { if (*s >= 'a' && *s <= 'z') { *s -= 32; } s++; } } void printAlign(pUnit** a, const int i, const int j, char* s, char* r, char* saln, char* raln, int n) { int k; pUnit p = a[i][j]; if (! i && ! j) { // 到達(0, 0) for (k = n - 1; k >= 0; k--) printf("%c", saln[k]); printf("n"); for (k = n - 1; k >= 0; k--) printf("%c", raln[k]); printf("nn"); return; } if (p->W1) { // 向上回溯一格 saln[n] = s[i - 1]; raln[n] = GAP_CHAR; printAlign(a, i - 1, j, s, r, saln, raln, n + 1); } if (p->W2) { // 向左上回溯一格 saln[n] = s[i - 1]; raln[n] = r[j - 1]; printAlign(a, i - 1, j - 1, s, r, saln, raln, n + 1); } if (p->W3) { // 向左回溯一格 saln[n] = GAP_CHAR; raln[n] = r[j - 1]; printAlign(a, i, j - 1, s, r, saln, raln, n + 1); } } void align(char *s, char *r) { int i, j; int m = strlen(s); int n = strlen(r); int m1, m2, m3, minm; pUnit **aUnit; char* salign; char* ralign; // 初始化 if ((aUnit = (pUnit **) malloc(sizeof(pUnit*) * (m + 1))) == NULL) { fputs("Error: Out of space!n", stderr); exit(1); } for (i = 0; i <= m; i++) { if ((aUnit[i] = (pUnit *) malloc(sizeof(pUnit) * (n + 1))) == NULL) { fputs("Error: Out of space!n", stderr); exit(1); } for (j = 0; j <= n; j++) { if ((aUnit[i][j] = (pUnit) malloc(sizeof(struct Unit))) == NULL) { fputs("Error: Out of space!n", stderr); exit(1); } aUnit[i][j]->W1 = 0; aUnit[i][j]->W2 = 0; aUnit[i][j]->W3 = 0; } } for (i = 0; i <= m; i++) { aUnit[i][0]->M = i; aUnit[i][0]->W1 = 1; } for (j = 1; j <= n; j++) { aUnit[0][j]->M = j; aUnit[0][j]->W3 = 1; } // 將字元串都變成大寫 strUpper(s); strUpper(r); // 動態規劃演算法計算得分矩陣每個單元的分值 for (i = 1; i <= m; i++) { for (j = 1; j <= n; j++) { m1 = aUnit[i - 1][j]->M + 1; m2 = aUnit[i - 1][j - 1]->M + diff(s[i - 1], r[j - 1]); m3 = aUnit[i][j - 1]->M + 1; minm = min(min(m1, m2), m3); aUnit[i][j]->M = minm; if (m1 == minm) aUnit[i][j]->W1 = 1; if (m2 == minm) aUnit[i][j]->W2 = 1; if (m3 == minm) aUnit[i][j]->W3 = 1; } } /* // 列印得分矩陣 for (i = 0; i <= m; i++) { for (j = 0; j <= n; j++) printf("%d ", aUnit[i][j]->M); printf("n"); } */ printf("min score: %dn", aUnit[m][n]->M); // 列印最優比對結果，如果有多個，全部列印 // 遞歸法 if (aUnit[m][n]->M == 0) { fputs("Two seqs are totally the same.n", stdout); } else { if ((salign = (char*) malloc(sizeof(char) * (m + n))) == NULL || (ralign = (char*) malloc(sizeof(char) * (m + n))) == NULL) { fputs("Error: Out of space!n", stderr); exit(1); } printAlign(aUnit, m, n, s, r, salign, ralign, 0); // 釋放記憶體 free(salign); free(ralign); } for (i = 0; i <= m; i++) { for (j = 0; j <= n; j++) free(aUnit[i][j]); free(aUnit[i]); } free(aUnit); }
