如何将Python实现的动态时间规整(DTW)代码移植为C语言版本用于序列比对

现有基础 #

• 已收集到大量可正常运行的Python版DTW实现，但本人已有一年多未编写C代码，需要将对应功能转换为C语言版本。
• C代码中已定义两个二维整型数组，DTW计算仅需对比数组右列元素codeLock(n)(1)与code(m)(1)，左列的1~5编号值无需参与计算，现有数组定义如下：

static int codeLock[6][2] = {
    {1, 0},
    {2, 671},
    {3, 1400},
    {4, 2000},
    {5, 2800},
};

static int code[6][2] = {
    {1, 0},
    {2, 600},
    {3, 1360},
    {4, 1990},
    {5, 2800},
};

已验证的Python端实现 #

欧氏距离矩阵计算 #

实现函数为compute_euclidean_distance_matrix，功能为计算两个序列的逐点欧氏距离平方，支持长度不同的序列输入，代码如下：

def compute_euclidean_distance_matrix(x, y) -> np.array:
    """Calculate distance matrix
    This method calcualtes the pairwise Euclidean distance between two sequences.
    The sequences can have different lengths.
    """
    dist = np.zeros((len(y), len(x)))
    for i in range(len(y)):
        for j in range(len(x)):
            dist[i,j] = (x[j]-y[i])**2
    return dist

累积成本矩阵计算 #

实现函数为compute_accumulated_cost_matrix，功能为基于欧氏距离矩阵，计算DTW对齐路径所需的累积成本，代码如下：

def compute_accumulated_cost_matrix(x, y) -> np.array:
    """Compute accumulated cost matrix for warp path using Euclidean distance
    """
    distances = compute_euclidean_distance_matrix(x, y)
    # Initialization
    cost = np.zeros((len(y), len(x)))
    cost[0,0] = distances[0,0]

    for i in range(1, len(y)):
        cost[i, 0] = distances[i, 0] + cost[i-1, 0]  
        
    for j in range(1, len(x)):
        cost[0, j] = distances[0, j] + cost[0, j-1]  

    # Accumulated warp path cost
    for i in range(1, len(y)):
        for j in range(1, len(x)):
            cost[i, j] = min(
                cost[i-1, j],    # insertion
                cost[i, j-1],    # deletion
                cost[i-1, j-1]   # match
            ) + distances[i, j] 
            
    return cost

测试验证结果 #

基础实现测试 #

测试输入序列为：

x = [0, 671, 1400, 2000, 2800]
y = [0, 600, 1360, 1990, 2800]

运行测试代码后，会依次输出欧氏距离矩阵、累积成本矩阵、上下翻转后的累积成本矩阵，对应控制台输出如下：

fastdtw库交叉验证 #

使用fastdtw库对同组输入进行DTW距离计算，代码如下：

x = [0, 671, 1400, 2000, 2800]
y = [0, 600, 1360, 1990, 2800]

dtw_distance, warp_path = fastdtw(x, y, dist=euclidean)

print("\ndtw_distance")
print(dtw_distance)

运行后输出的DTW距离结果如下：

实现诉求 #

• 需要可参考的DTW C语言实现逻辑
• 希望将现有两个Python函数的逻辑准确翻译为C语言版本，保证计算结果与已验证的Python实现一致

内容的提问来源于stack exchange，提问作者WRD BAC