DTW算法解决时序问题的解析实践

d(i,j) = |xᵢ - yⱼ|  (一维数据)

d(i,j) = √[(xᵢ₁-yⱼ₁)² + (xᵢ₂-yⱼ₂)² + ... + (xᵢₚ-yⱼₚ)²]  (p维数据)

D(i,j) = d(i,j) + min{
    D(i-1,j),     // 垂直移动（X序列前进，Y序列保持不变）
    D(i,j-1),     // 水平移动（Y序列前进，X序列保持不变）
    D(i-1,j-1)    // 对角移动（X和Y序列同时前进）
}

D(1,1) = d(1,1)
D(i,1) = d(i,1) + D(i-1,1), for i > 1
D(1,j) = d(1,j) + D(1,j-1), for j > 1

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import euclidean
import seaborn as sns
from typing import Tuple, List, Optional
import time
import warnings
warnings.filterwarnings('ignore')

class DTWAnalyzer:
    """
    DTW算法分析器
    
    该类实现了基础DTW算法和带约束的DTW算法，提供了计算距离矩阵、
    构建累积成本矩阵和回溯最优路径等核心功能。
    """
    
    def __init__(self):
        """
        初始化DTW分析器
        
        初始化三个主要属性：
        - distance_matrix: 存储两个序列之间的局部距离矩阵
        - cost_matrix: 存储累积成本矩阵
        - path: 存储最优对齐路径
        """
        self.distance_matrix = None
        self.cost_matrix = None
        self.path = None
        
    def compute_distance_matrix(self, x: np.ndarray, y: np.ndarray, 
                               distance_func=euclidean) -> np.ndarray:
        """
        计算两个序列之间的距离矩阵
        
        参数:
            x: 第一个时间序列，形状为(m,)或(m,d)
            y: 第二个时间序列，形状为(n,)或(n,d)
            distance_func: 距离函数，默认为欧几里得距离
            
        返回:
            distance_matrix: 形状为(m,n)的距离矩阵
        """
        m, n = len(x), len(y)
        distance_matrix = np.zeros((m, n))
        
        # 计算每对元素之间的距离
        for i in range(m):
            for j in range(n):
                # 对一维序列使用简单的绝对差值
                if x.ndim == 1 and y.ndim == 1:
                    distance_matrix[i, j] = abs(x[i] - y[j])
                # 对多维序列使用提供的距离函数
                else:
                    distance_matrix[i, j] = distance_func(x[i], y[j])
                    
        self.distance_matrix = distance_matrix
        return distance_matrix
    
    def dtw_basic(self, x: np.ndarray, y: np.ndarray) -> Tuple[float, np.ndarray]:
        """
        基础DTW算法实现
        
        使用动态规划方法计算两个序列之间的DTW距离和累积成本矩阵
        
        参数:
            x: 第一个时间序列
            y: 第二个时间序列
            
        返回:
            distance: DTW距离值
            cost_matrix: 累积成本矩阵
        """
        m, n = len(x), len(y)
        
        # 步骤1: 计算局部距离矩阵
        distance_matrix = self.compute_distance_matrix(x, y)
        
        # 步骤2: 初始化累积成本矩阵
        cost_matrix = np.full((m, n), np.inf)  # 初始化为无穷大
        cost_matrix[0, 0] = distance_matrix[0, 0]  # 设置起点
        
        # 步骤3: 填充第一行和第一列（边界条件）
        for i in range(1, m):
            cost_matrix[i, 0] = distance_matrix[i, 0] + cost_matrix[i-1, 0]
        for j in range(1, n):
            cost_matrix[0, j] = distance_matrix[0, j] + cost_matrix[0, j-1]
            
        # 步骤4: 使用动态规划填充整个矩阵
        for i in range(1, m):
            for j in range(1, n):
                # 当前点的累积成本 = 当前点的局部距离 + 三个可能前驱点中的最小累积成本
                cost_matrix[i, j] = distance_matrix[i, j] + min(
                    cost_matrix[i-1, j],      # 垂直移动（X序列前进，Y序列保持不变）
                    cost_matrix[i, j-1],      # 水平移动（Y序列前进，X序列保持不变）
                    cost_matrix[i-1, j-1]     # 对角移动（X和Y序列同时前进）
                )
        
        # 保存累积成本矩阵供后续使用
        self.cost_matrix = cost_matrix
        
        # 返回DTW距离（右下角元素）和累积成本矩阵
        return cost_matrix[m-1, n-1], cost_matrix
    
    def find_optimal_path(self, cost_matrix: np.ndarray) -> List[Tuple[int, int]]:
        """
        从累积成本矩阵回溯找到最优对齐路径
        
        参数:
            cost_matrix: DTW累积成本矩阵
            
        返回:
            path: 最优对齐路径，格式为[(i1,j1), (i2,j2), ...]
        """
        m, n = cost_matrix.shape
        path = []
        i, j = m-1, n-1  # 从右下角开始回溯
        
        # 回溯直到到达左上角
        while i > 0 or j > 0:
            path.append((i, j))
            
            # 处理边界情况
            if i == 0:
                j -= 1  # 只能水平移动
            elif j == 0:
                i -= 1  # 只能垂直移动
            else:
                # 选择成本最小的前驱点
                costs = [
                    cost_matrix[i-1, j-1],  # 对角移动
                    cost_matrix[i-1, j],    # 垂直移动
                    cost_matrix[i, j-1]     # 水平移动
                ]
                min_idx = np.argmin(costs)
                
                # 根据最小成本方向移动
                if min_idx == 0:
                    i, j = i-1, j-1  # 对角移动
                elif min_idx == 1:
                    i = i-1  # 垂直移动
                else:
                    j = j-1  # 水平移动
                    
        # 添加起点并反转路径（从起点到终点）
        path.append((0, 0))
        path.reverse()
        
        # 保存路径供后续使用
        self.path = path
        return path
    
    def dtw_with_constraints(self, x: np.ndarray, y: np.ndarray, 
                           window: Optional[int] = None) -> Tuple[float, np.ndarray]:
        """
        带约束的DTW算法（Sakoe-Chiba带约束）
        
        通过限制搜索窗口提高计算效率
        
        参数:
            x: 第一个时间序列
            y: 第二个时间序列
            window: 窗口大小，表示对角线两侧的最大偏移量
            
        返回:
            distance: DTW距离值
            cost_matrix: 累积成本矩阵
        """
        m, n = len(x), len(y)
        
        # 如果未指定窗口大小，使用默认值
        if window is None:
            window = max(m, n)  # 默认为无约束
            
        # 计算距离矩阵
        distance_matrix = self.compute_distance_matrix(x, y)
        
        # 初始化累积成本矩阵
        cost_matrix = np.full((m, n), np.inf)
        cost_matrix[0, 0] = distance_matrix[0, 0]
        
        # 动态规划填充矩阵（带窗口约束）
        for i in range(m):
            # 只计算窗口内的元素，窗口定义为对角线两侧的带状区域
            # max(0, i-window)确保下界不小于0
            # min(n, i+window+1)确保上界不超过n
            for j in range(max(0, i-window), min(n, i+window+1)):
                # 跳过已初始化的起点
                if i == 0 and j == 0:
                    continue
                    
                # 收集可能的前驱点的累积成本
                candidates = []
                if i > 0 and abs(i-1-j) <= window:
                    candidates.append(cost_matrix[i-1, j])  # 垂直移动
                if j > 0 and abs(i-j+1) <= window:
                    candidates.append(cost_matrix[i, j-1])  # 水平移动
                if i > 0 and j > 0:
                    candidates.append(cost_matrix[i-1, j-1])  # 对角移动
                    
                # 如果有有效的前驱点，更新当前点的累积成本
                if candidates:
                    cost_matrix[i, j] = distance_matrix[i, j] + min(candidates)
        
        # 保存累积成本矩阵供后续使用
        self.cost_matrix = cost_matrix
        
        # 返回DTW距离和累积成本矩阵
        return cost_matrix[m-1, n-1], cost_matrix

def run_comprehensive_experiments():
    """
    运行综合实验
    
    该函数执行一系列实验来评估不同DTW算法变体的性能和效果。
    实验包括:
    1. 生成多种类型的测试数据
    2. 对每种数据应用不同的DTW算法
    3. 收集性能指标和结果
    4. 进行详细的算法分析
    
    返回:
        experiment_results: 详细的实验结果字典
        performance_results: 性能评估结果
        visualizer: 可视化工具对象
    """
    print("=== DTW算法综合实验 ===\n")
    
    # 初始化实验所需的各个组件
    data_generator = ExperimentDataGenerator()  # 数据生成器
    dtw_analyzer = DTWAnalyzer()                # 基础DTW分析器
    fast_dtw = FastDTW(radius=2)                # 快速DTW算法（搜索半径为2）
    visualizer = DTWVisualizer()                # 可视化工具
    evaluator = DTWPerformanceEvaluator()       # 性能评估工具
    
    # 第一阶段: 生成多种类型的测试数据
    print("1. 生成测试数据...")
    test_data = []    # 存储测试数据对
    test_names = []   # 存储测试案例名称
    
    # 1.1 生成正弦波数据 - 模拟简单的周期性信号
    x_sine, y_sine = data_generator.generate_sine_waves(100, 0.1)
    test_data.append((x_sine, y_sine))
    test_names.append("Sine Waves")
    
    # 1.2 生成股票数据 - 模拟金融时间序列
    x_stock, y_stock = data_generator.generate_stock_like_data(150)
    test_data.append((x_stock, y_stock))
    test_names.append("Stock-like Data")
    
    # 1.3 生成语音数据 - 模拟语音信号
    x_speech, y_speech = data_generator.generate_speech_like_data(120)
    test_data.append((x_speech, y_speech))
    test_names.append("Speech-like Data")
    
    # 1.4 生成生物信号数据 - 模拟心电图等生物医学信号
    x_bio, y_bio = data_generator.generate_biological_signal(200)
    test_data.append((x_bio, y_bio))
    test_names.append("Biological Signal")
    
    print(f"生成了 {len(test_data)} 组测试数据\n")
    
    # 第二阶段: 算法性能评估
    print("2. 算法性能评估...")
    # 定义要评估的算法
    algorithms = {
        'Basic DTW': dtw_analyzer,    # 基础DTW算法
        'FastDTW': fast_dtw           # 快速DTW算法
    }
    
    # 执行性能评估，收集执行时间、内存使用等指标
    performance_results = evaluator.evaluate_algorithm_performance(
        algorithms, test_data, test_names
    )
    
    # 第三阶段: 详细实验分析
    print("\n3. 详细实验分析...")
    experiment_results = {}  # 存储详细实验结果
    
    # 对每个测试案例进行详细分析
    for i, (test_name, (x, y)) in enumerate(zip(test_names, test_data)):
        print(f"\n--- {test_name} 详细分析 ---")
        
        # 3.1 基础DTW分析
        distance, cost_matrix = dtw_analyzer.dtw_basic(x, y)
        path = dtw_analyzer.find_optimal_path(cost_matrix)
        
        # 3.2 约束DTW分析（使用窗口约束）
        # 窗口大小设为序列长度的1/4，平衡效率和准确性
        window_size = min(len(x), len(y)) // 4
        constrained_distance, constrained_matrix = dtw_analyzer.dtw_with_constraints(
            x, y, window=window_size
        )
        constrained_path = dtw_analyzer.find_optimal_path(constrained_matrix)
        
        # 3.3 FastDTW分析
        fast_distance, fast_path = fast_dtw.compute(x, y)
        
        # 3.4 存储分析结果
        experiment_results[test_name] = {
            'sequences': (x, y),                # 原始序列对
            'basic_dtw': {                      # 基础DTW结果
                'distance': distance,
                'cost_matrix': cost_matrix,
                'path': path
            },
            'constrained_dtw': {                # 约束DTW结果
                'distance': constrained_distance,
                'cost_matrix': constrained_matrix,
                'path': constrained_path
            },
            'fast_dtw': {                       # 快速DTW结果
                'distance': fast_distance,
                'path': fast_path
            }
        }
        
        # 3.5 输出关键结果指标
        print(f"基础DTW距离: {distance:.4f}")
        print(f"约束DTW距离: {constrained_distance:.4f}")
        print(f"快速DTW距离: {fast_distance:.4f}")
        print(f"路径长度比较: 基础={len(path)}, 约束={len(constrained_path)}, 快速={len(fast_path)}")
    
    # 返回实验结果、性能结果和可视化工具
    return experiment_results, performance_results, visualizer

# 执行综合实验
experiment_results, performance_results, visualizer = run_comprehensive_experiments()

def create_comprehensive_visualizations(experiment_results, visualizer):
    """
    创建综合可视化结果
    
    为每个测试案例生成一系列可视化图表，展示DTW算法的工作原理和结果。
    包括序列对比图、成本矩阵热图、对齐结果图和路径对比图等。
    
    参数:
        experiment_results: 实验结果字典
        visualizer: 可视化工具对象
    """
    
    print("\n4. 生成可视化结果...")
    
    # 为每个测试案例创建可视化
    for test_name, results in experiment_results.items():
        print(f"生成 {test_name} 的可视化图表...")
        
        # 提取实验数据
        x, y = results['sequences']                # 原始序列对
        basic_results = results['basic_dtw']       # 基础DTW结果
        constrained_results = results['constrained_dtw']  # 约束DTW结果
        fast_results = results['fast_dtw']         # 快速DTW结果
        
        # 1. 序列对比图 - 展示原始时间序列的形态和差异
        fig1 = visualizer.plot_sequences(
            x, y, 
            f"{test_name} - 时间序列对比"
        )
        
        # 2. 基础DTW成本矩阵和路径图 - 展示DTW的核心计算过程和最优路径
        fig2 = visualizer.plot_cost_matrix(
            basic_results['cost_matrix'], 
            basic_results['path'],
            f"{test_name} - 基础DTW成本矩阵与路径"
        )
        
        # 3. 约束DTW成本矩阵 - 展示窗口约束如何影响计算过程
        fig3 = visualizer.plot_cost_matrix(
            constrained_results['cost_matrix'],
            constrained_results['path'],
            f"{test_name} - 约束DTW成本矩阵与路径"
        )
        
        # 4. 序列对齐结果图 - 直观展示DTW如何对齐两个序列
        fig4 = visualizer.plot_alignment(
            x, y, 
            basic_results['path'],
            f"{test_name} - DTW序列对齐结果"
        )
        
        # 5. 时间规整路径对比图 - 展示时间变形的模式和程度
        fig5 = visualizer.plot_warping_path(
            basic_results['path'], 
            len(x), len(y),
            f"{test_name} - DTW时间规整路径"
        )
        
        # 在实际应用中，可以将图表保存为图片文件
        # 这里注释掉保存代码，实际使用时可以取消注释
        """
        # 创建输出目录
        output_dir = f"visualization_results/{test_name}"
        os.makedirs(output_dir, exist_ok=True)
        
        # 保存图片，使用高分辨率和紧凑布局
        fig1.savefig(f'{output_dir}/sequences.png', dpi=300, bbox_inches='tight')
        fig2.savefig(f'{output_dir}/basic_dtw_matrix.png', dpi=300, bbox_inches='tight')
        fig3.savefig(f'{output_dir}/constrained_dtw_matrix.png', dpi=300, bbox_inches='tight')
        fig4.savefig(f'{output_dir}/alignment.png', dpi=300, bbox_inches='tight')
        fig5.savefig(f'{output_dir}/warping_path.png', dpi=300, bbox_inches='tight')
        """
        
        # 关闭所有图形，释放内存
        plt.close('all')

# 生成可视化结果
create_comprehensive_visualizations(experiment_results, visualizer)

def quantitative_analysis(experiment_results, performance_results):
    """
    定量分析实验结果
    
    对实验结果进行系统的定量分析，包括距离度量分析、性能分析和路径特征分析。
    通过这些分析，可以客观评估不同DTW算法变体的优缺点。
    
    参数:
        experiment_results: 详细的实验结果字典
        performance_results: 性能评估结果
        
    返回:
        distance_analysis: 距离度量分析结果
    """
    
    print("\n=== 定量分析结果 ===\n")
    
    # 1. 距离度量分析 - 比较不同算法计算的DTW距离
    print("1. 距离度量分析:")
    print("-" * 50)
    
    distance_analysis = {}  # 存储距离分析结果
    
    # 对每个测试案例进行分析
    for test_name, results in experiment_results.items():
        # 提取各算法计算的距离
        basic_dist = results['basic_dtw']['distance']
        constrained_dist = results['constrained_dtw']['distance']
        fast_dist = results['fast_dtw']['distance']
        
        # 计算约束DTW和快速DTW相对于基础DTW的距离比率
        # 比率接近1表示近似算法的精度较高
        distance_analysis[test_name] = {
            'basic': basic_dist,                                 # 基础DTW距离
            'constrained': constrained_dist,                     # 约束DTW距离
            'fast': fast_dist,                                   # 快速DTW距离
            'constrained_ratio': constrained_dist / basic_dist,  # 约束DTW/基础DTW比率
            'fast_ratio': fast_dist / basic_dist                 # 快速DTW/基础DTW比率
        }
        
        # 输出分析结果
        print(f"{test_name}:")
        print(f"  基础DTW距离: {basic_dist:.4f} (基准值)")
        print(f"  约束DTW距离: {constrained_dist:.4f} (相对基准比率: {constrained_dist/basic_dist:.3f})")
        print(f"  快速DTW距离: {fast_dist:.4f} (相对基准比率: {fast_dist/basic_dist:.3f})")
        
        # 解释距离差异
        if constrained_dist/basic_dist > 1.05:
            print("  注: 约束DTW距离明显大于基础DTW，表明窗口约束可能导致次优路径")
        if fast_dist/basic_dist > 1.10:
            print("  注: 快速DTW距离显著大于基础DTW，表明多分辨率近似带来了一定精度损失")
        print()
    
    # 2. 性能分析 - 评估算法的执行效率
    print("2. 性能分析:")
    print("-" * 50)
    
    # 按算法分组计算平均性能指标
    algorithms = list(set(performance_results['algorithm']))
    for alg in algorithms:
        # 找出该算法的所有测试结果
        alg_indices = [i for i, a in enumerate(performance_results['algorithm']) if a == alg]
        
        # 计算平均执行时间和内存使用
        avg_time = np.mean([performance_results['time_cost'][i] for i in alg_indices])
        avg_memory = np.mean([performance_results['memory_usage'][i] for i in alg_indices])
        
        # 计算性能指标的标准差，评估算法在不同数据上的稳定性
        std_time = np.std([performance_results['time_cost'][i] for i in alg_indices])
        
        # 输出性能分析结果
        print(f"{alg}:")
        print(f"  平均执行时间: {avg_time:.6f}秒 (标准差: {std_time:.6f}秒)")
        print(f"  平均内存使用: {avg_memory/1024:.2f}KB")
        
        # 如果是快速DTW，计算相对于基础DTW的加速比
        if alg == 'FastDTW' and 'Basic DTW' in algorithms:
            basic_indices = [i for i, a in enumerate(performance_results['algorithm']) if a == 'Basic DTW']
            basic_avg_time = np.mean([performance_results['time_cost'][i] for i in basic_indices])
            speedup = basic_avg_time / avg_time if avg_time > 0 else float('inf')
            print(f"  相对基础DTW加速比: {speedup:.2f}倍")
        print()
    
    # 3. 路径特征分析 - 评估对齐路径的质量和特点
    print("3. 路径特征分析:")
    print("-" * 50)
    
    # 对每个测试案例分析路径特征
    for test_name, results in experiment_results.items():
        # 提取各算法生成的路径
        basic_path = results['basic_dtw']['path']
        constrained_path = results['constrained_dtw']['path']
        fast_path = results['fast_dtw']['path']
        
        # 定义路径分析函数
        def analyze_path(path):
            """分析路径的特征，包括长度、对角移动比例和扭曲程度"""
            if not path:
                return {'length': 0, 'diagonal_ratio': 0, 'warping_degree': 0}
            
            # 计算对角移动的比例（对角移动表示两个序列同步前进）
            diagonal_moves = sum(1 for i in range(1, len(path)) 
                               if path[i][0] - path[i-1][0] == 1 and 
                                  path[i][1] - path[i-1][1] == 1)
            diagonal_ratio = diagonal_moves / (len(path) - 1) if len(path) > 1 else 0
            
            # 计算扭曲程度（路径偏离对角线的平均距离）
            x_coords = [p[0] for p in path]
            y_coords = [p[1] for p in path]
            # 计算理想对角线上的预期y坐标
            expected_y = [x * len(y_coords) / len(x_coords) for x in x_coords]
            # 计算实际y坐标与预期y坐标的平均偏差
            warping_degree = np.mean([abs(y - ey) for y, ey in zip(y_coords, expected_y)])
            
            # 返回路径特征指标
            return {
                'length': len(path),              # 路径长度
                'diagonal_ratio': diagonal_ratio,  # 对角移动比例
                'warping_degree': warping_degree   # 扭曲程度
            }
        
        # 分析各算法的路径特征
        basic_analysis = analyze_path(basic_path)
        constrained_analysis = analyze_path(constrained_path)
        fast_analysis = analyze_path(fast_path)
        
        # 输出路径分析结果
        print(f"{test_name}:")
        print(f"  基础DTW路径: 长度={basic_analysis['length']}, "
              f"对角比例={basic_analysis['diagonal_ratio']:.3f}, "
              f"扭曲度={basic_analysis['warping_degree']:.3f}")
        print(f"  约束DTW路径: 长度={constrained_analysis['length']}, "
              f"对角比例={constrained_analysis['diagonal_ratio']:.3f}, "
              f"扭曲度={constrained_analysis['warping_degree']:.3f}")
        print(f"  快速DTW路径: 长度={fast_analysis['length']}, "
              f"对角比例={fast_analysis['diagonal_ratio']:.3f}, "
              f"扭曲度={fast_analysis['warping_degree']:.3f}")
        
        # 解释路径特征差异
        if constrained_analysis['diagonal_ratio'] > basic_analysis['diagonal_ratio'] * 1.1:
            print("  注: 约束DTW的对角比例明显高于基础DTW，表明窗口约束倾向于更直接的对齐")
        if abs(fast_analysis['warping_degree'] - basic_analysis['warping_degree']) > basic_analysis['warping_degree'] * 0.2:
            print("  注: 快速DTW的扭曲度与基础DTW有显著差异，表明多分辨率方法影响了对齐模式")
        print()
    
    # 返回距离分析结果，供后续使用
    return distance_analysis

# 执行定量分析
distance_analysis = quantitative_analysis(experiment_results, performance_results)