# 序列比对（24）最长公共子序列

### 动态规划求解最长公共子序列的代码

```#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MAXSEQ 1000
#define max(a, b) ((a) > (b) ? (a) : (b))
#define min(a, b) ((a) < (b) ? (a) : (b))

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, int* is, int* ir, 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 != '\0') {
if (*s >= 'a' && *s <= 'z') {
*s -= 32;
}
s++;
}
}

void printAlign(pUnit** a, const int i, const int j, char* s, char* r, int* is, int* ir, int n) {
int k;
pUnit p = a[i][j];
if (p->M == 0) {   // 到值为0的矩阵单元就结束
printf("index of common seq in seq1: ");
for (k = n - 1; k >= 0; k--)
printf("%d ", is[k]);
printf("\n");
printf("index of common seq in seq2: ");
for (k = n - 1; k >= 0; k--)
printf("%d ", ir[k]);
printf("\n");
for (k = n - 1; k >= 0; k--)
printf("%c", s[is[k] - 1]);
printf("\n");
for (k = n - 1; k >= 0; k--)
printf("%c", r[ir[k] - 1]);
printf("\n\n");
return;
}
if (p->W1)     // 向上回溯一格
printAlign(a, i - 1, j, s, r, is, ir, n);
if (p->W2) {    // 向左上回溯一格
is[n] = i;
ir[n] = j;
printAlign(a, i - 1, j - 1, s, r, is, ir, n + 1);
}
if (p->W3)     // 向左回溯一格
printAlign(a, i, j - 1, s, r, is, ir, n);
}

void align(char *s, char *r) {
int i, j;
int m = strlen(s);
int n = strlen(r);
int ml = min(m, n);
int m1, m2, m3, maxm;
pUnit **aUnit;
int* sidx;
int* ridx;
// 初始化
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 = 0;
for (j = 1; j <= n; j++)
aUnit[0][j]->M = 0;
// 将字符串都变成大写
strUpper(s);
strUpper(r);
// 动态规划算法计算得分矩阵每个单元的分值
for (i = 1; i <= m; i++) {
for (j = 1; j <= n; j++) {
m1 = aUnit[i - 1][j]->M;
m2 = s[i - 1] == r[j - 1] ? aUnit[i - 1][j - 1]->M + 1 : -1;
m3 = aUnit[i][j - 1]->M;
maxm = max(max(m1, m2), m3);
aUnit[i][j]->M = maxm;
if (m1 == maxm) aUnit[i][j]->W1 = 1;
if (m2 == maxm) aUnit[i][j]->W2 = 1;
if (m3 == maxm) 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("max score: %d\n", aUnit[m][n]->M);
// 打印最优比对结果，如果有多个，全部打印
// 递归法
if (aUnit[m][n]->M == 0) {
fputs("No common seq found.\n", stdout);
} else {
if ((sidx = (int*) malloc(sizeof(int) * ml)) == NULL || \
(ridx = (int*) malloc(sizeof(int) * ml)) == NULL) {
fputs("Error: Out of space!\n", stderr);
exit(1);
}
printAlign(aUnit, m, n, s, r, sidx, ridx, 0);
// 释放内存
free(sidx);
free(ridx);
}
for (i = 0; i <= m; i++) {
for (j = 0; j <= n; j++)
free(aUnit[i][j]);
free(aUnit[i]);
}
free(aUnit);
}```

0 条评论

• ### 序列比对（九）从掷骰子说起HMM

前面八篇文章介绍了动态规划在序列比对中的基础应用。从本文开始，开始介绍HMM（隐马尔可夫模型）。为什么要介绍它呢？因为该模型在发现Motif、预测CpG岛等多方...

• ### 序列比对（14）viterbi算法和后验解码的比较

前文《序列比对（十）viterbi算法求解最可能路径》介绍了用viterbi算法求解最可能路径：在符号序列已知而状态序列未知时，最可能路径是：

• ### 序列比对（22）中间字符串分支定界方法中更紧的界

前文《序列比对（21）中间字符串问题的算法及实现代码》介绍了中间字符串的算法和代码，但是使用分支定界策略时所使用的界限是很宽松的。分支定界法的伪代码如下：

• ### Spring 复盘 | AOP

Spring AOP 基础 Java 动态代理实现，阅读文章之前，你最好有以下基础：

• ### 【Basic algorithm学习笔记】 排序

原理： 其实原理很简单 大白话来讲 就是 先定义很多 容器 即这里形象的表达 为桶，之后每次输入的数字则为选择哪个桶，将对应的 桶中扔一个标志物表示有这个 桶...

• ### Fence Repair （PKU 3253）

思路：用朴素的方法实现的话，时间复杂度为O（n^n）。因为只需要从板的集合中取出最短的两块，并且把长度为两块长度之和的板加入集合中即可，所有使用优先队列就可以高...

• ### 聊聊artemis的NetworkHealthCheck

activemq-artemis-2.11.0/artemis-commons/src/main/java/org/apache/activemq/artemi...

• ### LeetCode 1390. 四因数

来源：力扣（LeetCode） 链接：https://leetcode-cn.com/problems/four-divisors 著作权归领扣网络所有。商...

• ### 聊聊artemis的NetworkHealthCheck

activemq-artemis-2.11.0/artemis-commons/src/main/java/org/apache/activemq/artemi...