Needleman-Wunsch 的维基百科伪代码中的路径是否为 1 索引?
Are the paths 1-indexed in the wikipedia pseudocode for Needleman-Wunsch?
https://en.wikipedia.org/wiki/Needleman%E2%80%93Wunsch_algorithm中的伪代码显示
for i = 1 to length(A)
for j = 1 to length(B)
{
Match ← F(i−1, j−1) + S(Ai, Bj)
Delete ← F(i−1, j) + d
Insert ← F(i, j−1) + d
F(i,j) ← max(Match, Insert, Delete)
}
对于路径 A
和 B
,评分矩阵 F
将有 A.length+1
行和 B.length+1
列。 i,j
指的是分数矩阵中的位置。但是,在相似度函数中,它使用路径的 i
和 j
索引。对于这个循环,路径的第一个元素永远不会被访问。路径不是 0 索引的吗?
是的,序列 A
和 B
是 1 索引的,to length(A)
应该读作 <= length(A)
。
这在动态规划伪代码中很常见,例如维基百科上给出的Longest Common Subsequence伪代码是:
function LCSLength(X[1..m], Y[1..n])
C = array(0..m, 0..n)
for i := 0..m
C[i,0] = 0
for j := 0..n
C[0,j] = 0
for i := 1..m
for j := 1..n
if X[i] = Y[j]
C[i,j] := C[i-1,j-1] + 1
else
C[i,j] := max(C[i,j-1], C[i-1,j])
return C[m,n]
(其中 ..
包含在内,如 Ruby)
这是为了方便起见——如果序列是从 0 开始索引的,则有两个典型的更改会使代码更加冗长:
- 所有序列访问索引从
Sequence[i][j]
变为Sequence[i-1][j-1]
- 如果端点不包含在您的语言中,运行
length(Sequence)
应该 运行 length(Sequence) + 1
的所有循环。这是必要的,因为查找 table 的大小为 length(Sequence) + 1
.
...或者,用您的 0 索引语言填充序列,使它们成为 1 索引,然后在不进行上述调整的情况下应用算法(这可能令人困惑且效率低下,但可能)。
然后,维基百科提供了按照上述 2 条规则翻译成 C#(0 索引语言)的上述代码:
static int[,] LcsLength(string a, string b)
{
int[,] C = new int[a.Length + 1, b.Length + 1]; // (a, b).Length + 1
for (int i = 0; i < a.Length; i++)
C[i, 0] = 0;
for (int j = 0; j < b.Length; j++)
C[0, j] = 0;
for (int i = 1; i <= a.Length; i++)
for (int j = 1; j <= b.Length; j++)
{
if (a[i - 1] == b[j - 1])//i-1,j-1
C[i, j] = C[i - 1, j - 1] + 1;
else
C[i, j] = Math.Max(C[i, j - 1], C[i - 1, j]);
}
return C;
}
我将使用相同的规则将 Needleman-Wunsch 的维基百科伪代码翻译成 Python。
我从 here 中抓取了一些测试用例,但认为除此之外还没有经过测试。
此外,正如您在链接要点中看到的那样,将新元素附加到列表并反转它们比在前面添加字符串更有效,如下所示,但我选择了维基百科的逐字翻译。
def needleman_wunsch_make_F(A, B, S, d):
F = [[0] * (len(B) + 1) for _ in range(len(A) + 1)]
for i in range(0, len(A) + 1):
F[i][0] = d * i
for j in range(0, len(B) + 1):
F[0][j] = d * j
for i in range(1, len(A) + 1):
for j in range(1, len(B) + 1):
match = F[i-1][j-1] + S(A[i-1], B[j-1])
delete = F[i-1][j] + d
insert = F[i][j-1] + d
F[i][j] = max(match, insert, delete)
return F
def needleman_wunsch_align(A, B, F, S, d):
AlignmentA = "" # TODO use lists for efficiency
AlignmentB = ""
i = len(F) - 1
j = len(F[0]) - 1
while i > 0 or j > 0:
if i > 0 and j > 0 and F[i][j] == F[i-1][j-1] + S(A[i-1], B[j-1]):
AlignmentA = A[i-1] + AlignmentA
AlignmentB = B[j-1] + AlignmentB
i -= 1
j -= 1
elif i > 0 and F[i][j] == F[i-1][j] + d:
AlignmentA = A[i-1] + AlignmentA
AlignmentB = "-" + AlignmentB
i -= 1
else:
AlignmentA = "-" + AlignmentA
AlignmentB = B[j-1] + AlignmentB
j -= 1
return AlignmentA, AlignmentB
def needleman_wunsch(A, B, S=lambda a, b: 1 if a == b else -1, d=-1):
F = needleman_wunsch_make_F(A, B, S, d)
return needleman_wunsch_align(A, B, F, S, d)
if __name__ == "__main__":
import collections
Test = collections.namedtuple("Test", "A B exp_A exp_B d")
tests = (
Test(A="GATTACA", B="GCATGCU", exp_A="G-ATTACA", exp_B="GCA-TGCU", d=-1),
Test(
A="GCAGGCAAGTGGGGCACCCGTATCCTTTCCAACTTACAAGGGTCCCCGTT",
B="GTGCGCCAGAGGAAGTCACTTTATATCCGCGCACGGTACTCCTTTTTCTA",
exp_A="----G-C--AGGCAAGTGGGGCACCCGTATCCT-T-T-C-C-AACTTACAAGGGT-C-CC-----CGT-T",
exp_B="GTGCGCCAGAGG-AAGT----CA--C-T-T--TATATCCGCG--C--AC---GGTACTCCTTTTTC-TA-",
d=0
),
Test(
A="GCAGGCAAGTGGGGCACCCGTATCCTTTCCAACTTACAAGGGTCCCCGTT",
B="GTGCGCCAGAGGAAGTCACTTTATATCCGCGCACGGTACTCCTTTTTCTA",
exp_A="GCAG-GCAAGTGG--GGCAC-CCGTATCCTTTC-CAAC-TTACAAGGGTCC-CCGT-T-",
exp_B="G-TGCGCCAGAGGAAGTCACTTTATATCC--GCGC-ACGGTAC-----TCCTTTTTCTA",
d=-1
),
Test(
A="GCAGGCAAGTGGGGCACCCGTATCCTTTCCAACTTACAAGGGTCCCCGTT",
B="GTGCGCCAGAGGAAGTCACTTTATATCCGCGCACGGTACTCCTTTTTCTA",
exp_A="GCAGGCAAGTGG--GGCAC-CCGTATCCTTTCCAACTTACAAGGGTCCCCGTT",
exp_B="GTGCGCCAGAGGAAGTCACTTTATATCC-GCGCACGGTAC-TCCTTTTTC-TA",
d=-2
),
)
for A, B, exp_A, exp_B, d in tests:
act_A, act_B = needleman_wunsch(A, B, d=d)
assert act_A == exp_A
assert act_B == exp_B
顺便说一句,我写了一个小应用程序,将用于编写 DP 算法的 1 索引伪代码语言转换为 0 索引 Python。如果你好奇的话,它在 GitHub 上。
https://en.wikipedia.org/wiki/Needleman%E2%80%93Wunsch_algorithm中的伪代码显示
for i = 1 to length(A)
for j = 1 to length(B)
{
Match ← F(i−1, j−1) + S(Ai, Bj)
Delete ← F(i−1, j) + d
Insert ← F(i, j−1) + d
F(i,j) ← max(Match, Insert, Delete)
}
对于路径 A
和 B
,评分矩阵 F
将有 A.length+1
行和 B.length+1
列。 i,j
指的是分数矩阵中的位置。但是,在相似度函数中,它使用路径的 i
和 j
索引。对于这个循环,路径的第一个元素永远不会被访问。路径不是 0 索引的吗?
是的,序列 A
和 B
是 1 索引的,to length(A)
应该读作 <= length(A)
。
这在动态规划伪代码中很常见,例如维基百科上给出的Longest Common Subsequence伪代码是:
function LCSLength(X[1..m], Y[1..n])
C = array(0..m, 0..n)
for i := 0..m
C[i,0] = 0
for j := 0..n
C[0,j] = 0
for i := 1..m
for j := 1..n
if X[i] = Y[j]
C[i,j] := C[i-1,j-1] + 1
else
C[i,j] := max(C[i,j-1], C[i-1,j])
return C[m,n]
(其中 ..
包含在内,如 Ruby)
这是为了方便起见——如果序列是从 0 开始索引的,则有两个典型的更改会使代码更加冗长:
- 所有序列访问索引从
Sequence[i][j]
变为Sequence[i-1][j-1]
- 如果端点不包含在您的语言中,运行
length(Sequence)
应该 运行length(Sequence) + 1
的所有循环。这是必要的,因为查找 table 的大小为length(Sequence) + 1
.
...或者,用您的 0 索引语言填充序列,使它们成为 1 索引,然后在不进行上述调整的情况下应用算法(这可能令人困惑且效率低下,但可能)。
然后,维基百科提供了按照上述 2 条规则翻译成 C#(0 索引语言)的上述代码:
static int[,] LcsLength(string a, string b)
{
int[,] C = new int[a.Length + 1, b.Length + 1]; // (a, b).Length + 1
for (int i = 0; i < a.Length; i++)
C[i, 0] = 0;
for (int j = 0; j < b.Length; j++)
C[0, j] = 0;
for (int i = 1; i <= a.Length; i++)
for (int j = 1; j <= b.Length; j++)
{
if (a[i - 1] == b[j - 1])//i-1,j-1
C[i, j] = C[i - 1, j - 1] + 1;
else
C[i, j] = Math.Max(C[i, j - 1], C[i - 1, j]);
}
return C;
}
我将使用相同的规则将 Needleman-Wunsch 的维基百科伪代码翻译成 Python。
我从 here 中抓取了一些测试用例,但认为除此之外还没有经过测试。
此外,正如您在链接要点中看到的那样,将新元素附加到列表并反转它们比在前面添加字符串更有效,如下所示,但我选择了维基百科的逐字翻译。
def needleman_wunsch_make_F(A, B, S, d):
F = [[0] * (len(B) + 1) for _ in range(len(A) + 1)]
for i in range(0, len(A) + 1):
F[i][0] = d * i
for j in range(0, len(B) + 1):
F[0][j] = d * j
for i in range(1, len(A) + 1):
for j in range(1, len(B) + 1):
match = F[i-1][j-1] + S(A[i-1], B[j-1])
delete = F[i-1][j] + d
insert = F[i][j-1] + d
F[i][j] = max(match, insert, delete)
return F
def needleman_wunsch_align(A, B, F, S, d):
AlignmentA = "" # TODO use lists for efficiency
AlignmentB = ""
i = len(F) - 1
j = len(F[0]) - 1
while i > 0 or j > 0:
if i > 0 and j > 0 and F[i][j] == F[i-1][j-1] + S(A[i-1], B[j-1]):
AlignmentA = A[i-1] + AlignmentA
AlignmentB = B[j-1] + AlignmentB
i -= 1
j -= 1
elif i > 0 and F[i][j] == F[i-1][j] + d:
AlignmentA = A[i-1] + AlignmentA
AlignmentB = "-" + AlignmentB
i -= 1
else:
AlignmentA = "-" + AlignmentA
AlignmentB = B[j-1] + AlignmentB
j -= 1
return AlignmentA, AlignmentB
def needleman_wunsch(A, B, S=lambda a, b: 1 if a == b else -1, d=-1):
F = needleman_wunsch_make_F(A, B, S, d)
return needleman_wunsch_align(A, B, F, S, d)
if __name__ == "__main__":
import collections
Test = collections.namedtuple("Test", "A B exp_A exp_B d")
tests = (
Test(A="GATTACA", B="GCATGCU", exp_A="G-ATTACA", exp_B="GCA-TGCU", d=-1),
Test(
A="GCAGGCAAGTGGGGCACCCGTATCCTTTCCAACTTACAAGGGTCCCCGTT",
B="GTGCGCCAGAGGAAGTCACTTTATATCCGCGCACGGTACTCCTTTTTCTA",
exp_A="----G-C--AGGCAAGTGGGGCACCCGTATCCT-T-T-C-C-AACTTACAAGGGT-C-CC-----CGT-T",
exp_B="GTGCGCCAGAGG-AAGT----CA--C-T-T--TATATCCGCG--C--AC---GGTACTCCTTTTTC-TA-",
d=0
),
Test(
A="GCAGGCAAGTGGGGCACCCGTATCCTTTCCAACTTACAAGGGTCCCCGTT",
B="GTGCGCCAGAGGAAGTCACTTTATATCCGCGCACGGTACTCCTTTTTCTA",
exp_A="GCAG-GCAAGTGG--GGCAC-CCGTATCCTTTC-CAAC-TTACAAGGGTCC-CCGT-T-",
exp_B="G-TGCGCCAGAGGAAGTCACTTTATATCC--GCGC-ACGGTAC-----TCCTTTTTCTA",
d=-1
),
Test(
A="GCAGGCAAGTGGGGCACCCGTATCCTTTCCAACTTACAAGGGTCCCCGTT",
B="GTGCGCCAGAGGAAGTCACTTTATATCCGCGCACGGTACTCCTTTTTCTA",
exp_A="GCAGGCAAGTGG--GGCAC-CCGTATCCTTTCCAACTTACAAGGGTCCCCGTT",
exp_B="GTGCGCCAGAGGAAGTCACTTTATATCC-GCGCACGGTAC-TCCTTTTTC-TA",
d=-2
),
)
for A, B, exp_A, exp_B, d in tests:
act_A, act_B = needleman_wunsch(A, B, d=d)
assert act_A == exp_A
assert act_B == exp_B
顺便说一句,我写了一个小应用程序,将用于编写 DP 算法的 1 索引伪代码语言转换为 0 索引 Python。如果你好奇的话,它在 GitHub 上。