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)
    }

对于路径 AB,评分矩阵 F 将有 A.length+1 行和 B.length+1 列。 i,j 指的是分数矩阵中的位置。但是,在相似度函数中,它使用路径的 ij 索引。对于这个循环,路径的第一个元素永远不会被访问。路径不是 0 索引的吗?

是的,序列 AB 是 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 上。