DP 在工业界：资源调度、广告投放与路径规划

class PacingPIDController:
    """PID 控制器实现广告预算 pacing。

    通过动态调整 lambda（花费惩罚系数），
    使得广告主的花费曲线尽量贴合理想的 pacing 曲线。
    """

    def __init__(self, budget: float, total_periods: int,
                 kp: float = 0.5, ki: float = 0.1, kd: float = 0.05):
        self.budget = budget
        self.total_periods = total_periods
        self.kp = kp  # 比例系数
        self.ki = ki  # 积分系数
        self.kd = kd  # 微分系数

        self.lambda_val = 1.0       # 初始 lambda
        self.integral_error = 0.0   # 累积误差
        self.prev_error = 0.0       # 上一次误差
        self.spent = 0.0            # 已花费

    def ideal_spend(self, period: int) -> float:
        """理想花费曲线：均匀分配。"""
        return self.budget * (period / self.total_periods)

    def update(self, period: int, period_spend: float) -> float:
        """根据当前花费情况更新 lambda。

        返回值:
            更新后的 lambda 值。
        """
        self.spent += period_spend

        # 误差 = 实际花费 - 理想花费（正值意味着花多了）
        error = self.spent - self.ideal_spend(period)

        # PID 更新
        self.integral_error += error
        derivative = error - self.prev_error
        self.prev_error = error

        delta = (self.kp * error +
                 self.ki * self.integral_error +
                 self.kd * derivative)

        # lambda 越大，越倾向于不花钱
        self.lambda_val = max(0.01, self.lambda_val + delta / self.budget)

        return self.lambda_val

    def should_bid(self, value: float, cost: float) -> bool:
        """根据当前 lambda 决定是否参与竞价。"""
        return value - self.lambda_val * cost > 0

def offline_dp_lambda_table(
    budget: float,
    num_periods: int,
    num_budget_levels: int,
    historical_value_distributions: list
) -> list:
    """离线预计算 lambda 查找表。

    dp[t][b] = 从时段 t 开始、剩余预算为 b 时的最优期望价值。
    lambda_table[t][b] = 对应的最优 lambda。
    """
    budget_step = budget / num_budget_levels

    # dp[t][b]: 从时段 t 开始，剩余预算 b 的最优期望价值
    dp = [[0.0] * (num_budget_levels + 1) for _ in range(num_periods + 1)]
    lambda_table = [[0.0] * (num_budget_levels + 1)
                    for _ in range(num_periods)]

    for t in range(num_periods - 1, -1, -1):
        dist = historical_value_distributions[t]
        for b in range(num_budget_levels + 1):
            remaining_budget = b * budget_step
            best_val = dp[t + 1][b]
            best_lambda = 0.0

            # 尝试不同的 lambda 值
            for lam in [i * 0.1 for i in range(1, 51)]:
                expected_val = 0.0
                expected_cost = 0.0

                for v, c, prob in dist:
                    if v - lam * c > 0:
                        expected_val += v * prob
                        expected_cost += c * prob

                # 预计花费后的剩余预算
                new_b = max(0, int((remaining_budget - expected_cost)
                                   / budget_step))
                total = expected_val + dp[t + 1][new_b]

                if total > best_val:
                    best_val = total
                    best_lambda = lam

            dp[t][b] = best_val
            lambda_table[t][b] = best_lambda

    return lambda_table

class OnlineLambdaLearner:
    """在线梯度下降更新 lambda。"""

    def __init__(self, lambda_table: list, learning_rate: float = 0.01):
        self.lambda_table = lambda_table
        self.lr = learning_rate
        self.lambda_correction = 0.0

    def get_lambda(self, period: int, remaining_budget_level: int) -> float:
        """获取修正后的 lambda。"""
        base_lambda = self.lambda_table[period][remaining_budget_level]
        return max(0.01, base_lambda + self.lambda_correction)

    def update(self, gradient: float):
        """根据梯度更新修正量。

        gradient > 0 表示花费偏多，需要增大 lambda。
        """
        self.lambda_correction += self.lr * gradient

import numpy as np
from typing import Tuple


def dynamic_pricing_value_iteration(
    num_periods: int,
    max_inventory: int,
    prices: list,
    demand_pmf: callable,
    discount: float = 1.0
) -> Tuple[np.ndarray, np.ndarray]:
    """有限 horizon 动态定价的价值迭代。

    参数:
        num_periods: 总时段数。
        max_inventory: 最大库存量。
        prices: 可选价格列表。
        demand_pmf: demand_pmf(price, d) 返回价格为 price 时需求为 d 的概率。
        discount: 折扣因子。

    返回:
        (V, policy): V[t][s] 是最优值函数，policy[t][s] 是最优价格。
    """
    V = np.zeros((num_periods + 1, max_inventory + 1))
    policy = np.zeros((num_periods, max_inventory + 1))

    # 从最后一个时段往前迭代
    for t in range(num_periods - 1, -1, -1):
        for s in range(max_inventory + 1):
            if s == 0:
                V[t][s] = 0.0
                continue

            best_val = -np.inf
            best_price = prices[0]

            for p in prices:
                expected = 0.0
                for d in range(s + 1):
                    prob = demand_pmf(p, d)
                    sold = min(d, s)
                    revenue = p * sold
                    future = V[t + 1][s - sold]
                    expected += prob * (revenue + discount * future)

                # 需求可能超过库存
                for d in range(s + 1, s + 50):
                    prob = demand_pmf(p, d)
                    if prob < 1e-10:
                        break
                    revenue = p * s  # 全部卖完
                    future = V[t + 1][0]
                    expected += prob * (revenue + discount * future)

                if expected > best_val:
                    best_val = expected
                    best_price = p

            V[t][s] = best_val
            policy[t][s] = best_price

    return V, policy

def cdn_cache_placement_dp(
    items: list,
    capacity: int
) -> Tuple[float, list]:
    """CDN 缓存放置的 0-1 背包 DP。

    参数:
        items: [(content_id, size, value), ...] 内容列表。
        capacity: 缓存容量（归一化为整数单位）。

    返回:
        (max_value, selected_ids): 最大价值和选中的内容 ID。
    """
    n = len(items)
    dp = [0.0] * (capacity + 1)
    choice = [[False] * (capacity + 1) for _ in range(n)]

    for i in range(n):
        cid, size, value = items[i]
        for c in range(capacity, size - 1, -1):
            if dp[c - size] + value > dp[c]:
                dp[c] = dp[c - size] + value
                choice[i][c] = True

    # 回溯
    selected = []
    c = capacity
    for i in range(n - 1, -1, -1):
        if choice[i][c]:
            selected.append(items[i][0])
            c -= items[i][1]

    return dp[capacity], selected

class OnlineCDNCache:
    """基于阈值的在线缓存策略。

    使用类似广告 pacing 的 lambda 阈值机制。
    """

    def __init__(self, capacity: int, threshold: float = 1.0):
        self.capacity = capacity
        self.threshold = threshold  # 类似 lambda 的阈值
        self.cache = {}  # content_id -> (size, estimated_value)
        self.used = 0

    def _value_density(self, value: float, size: int) -> float:
        """单位存储的价值密度。"""
        return value / max(size, 1)

    def request(self, content_id: str, size: int,
                estimated_value: float) -> bool:
        """处理一个内容请求。

        如果内容不在缓存中，决定是否缓存。
        返回是否命中缓存。
        """
        if content_id in self.cache:
            # 命中，更新价值估计
            self.cache[content_id] = (size, estimated_value)
            return True

        density = self._value_density(estimated_value, size)
        if density < self.threshold:
            return False  # 价值密度不够，不缓存

        # 如果空间不够，淘汰低密度内容
        while self.used + size > self.capacity and self.cache:
            min_id = min(self.cache,
                         key=lambda k: self._value_density(
                             self.cache[k][1], self.cache[k][0]))
            min_density = self._value_density(
                self.cache[min_id][1], self.cache[min_id][0])
            if min_density >= density:
                break  # 淘汰代价太高
            self.used -= self.cache[min_id][0]
            del self.cache[min_id]

        if self.used + size <= self.capacity:
            self.cache[content_id] = (size, estimated_value)
            self.used += size
            return False  # 未命中但已缓存

        return False  # 未命中且未缓存

import numpy as np
from typing import List, Tuple


class HiddenMarkovModel:
    """隐马尔可夫模型的完整实现。"""

    def __init__(self, num_states: int, num_observations: int):
        self.N = num_states
        self.M = num_observations
        self.pi = np.ones(num_states) / num_states
        self.A = np.ones((num_states, num_states)) / num_states
        self.B = np.ones((num_states, num_observations)) / num_observations

    def set_params(self, pi: np.ndarray, A: np.ndarray, B: np.ndarray):
        """设置模型参数。"""
        self.pi = pi.copy()
        self.A = A.copy()
        self.B = B.copy()

    def viterbi(self, observations: List[int]) -> Tuple[List[int], float]:
        """Viterbi 算法：求最优隐状态序列。

        参数:
            observations: 观测序列（整数索引）。

        返回:
            (best_path, log_prob): 最优路径和对应的对数概率。
        """
        T = len(observations)
        N = self.N

        # 使用对数概率避免下溢
        log_pi = np.log(self.pi + 1e-300)
        log_A = np.log(self.A + 1e-300)
        log_B = np.log(self.B + 1e-300)

        # delta[t][j]: 在时刻 t 处于状态 j 的最大对数概率
        delta = np.full((T, N), -np.inf)
        # psi[t][j]: 回溯指针
        psi = np.zeros((T, N), dtype=int)

        # 初始化
        delta[0] = log_pi + log_B[:, observations[0]]

        # 递推
        for t in range(1, T):
            for j in range(N):
                # delta[t][j] = max_i(delta[t-1][i] + log_A[i][j]) + log_B[j][o_t]
                scores = delta[t - 1] + log_A[:, j]
                psi[t][j] = np.argmax(scores)
                delta[t][j] = scores[psi[t][j]] + log_B[j, observations[t]]

        # 终止
        best_last = np.argmax(delta[T - 1])
        log_prob = delta[T - 1][best_last]

        # 回溯
        path = [0] * T
        path[T - 1] = best_last
        for t in range(T - 2, -1, -1):
            path[t] = psi[t + 1][path[t + 1]]

        return path, log_prob

    def forward(self, observations: List[int]) -> float:
        """前向算法：计算观测序列的对数似然。"""
        T = len(observations)
        log_alpha = np.log(self.pi + 1e-300) + np.log(
            self.B[:, observations[0]] + 1e-300)

        for t in range(1, T):
            new_log_alpha = np.full(self.N, -np.inf)
            for j in range(self.N):
                log_terms = log_alpha + np.log(self.A[:, j] + 1e-300)
                max_term = np.max(log_terms)
                new_log_alpha[j] = (max_term
                                    + np.log(np.sum(np.exp(log_terms - max_term)))
                                    + np.log(self.B[j, observations[t]] + 1e-300))
            log_alpha = new_log_alpha

        max_alpha = np.max(log_alpha)
        return max_alpha + np.log(np.sum(np.exp(log_alpha - max_alpha)))

import math
from typing import List, Tuple, Dict


class GPSMapMatcher:
    """基于 HMM + Viterbi 的 GPS 轨迹匹配。"""

    def __init__(self, sigma_z: float = 4.07, beta: float = 3.0):
        """
        参数:
            sigma_z: GPS 测量噪声的标准差（米）。
            beta: 转移概率的参数。
        """
        self.sigma_z = sigma_z
        self.beta = beta

    def emission_prob(self, gps_point: Tuple[float, float],
                      candidate_point: Tuple[float, float]) -> float:
        """发射概率：GPS 点与候选路段投影点之间的距离。"""
        dist = self._haversine(gps_point, candidate_point)
        return (1.0 / (math.sqrt(2 * math.pi) * self.sigma_z)
                * math.exp(-0.5 * (dist / self.sigma_z) ** 2))

    def transition_prob(self, route_dist: float,
                        great_circle_dist: float) -> float:
        """转移概率：路网距离与直线距离的差异。"""
        dt = abs(route_dist - great_circle_dist)
        return (1.0 / self.beta) * math.exp(-dt / self.beta)

    def match(
        self,
        gps_points: List[Tuple[float, float]],
        candidates: List[List[Dict]],
        route_distances: callable
    ) -> List[int]:
        """Viterbi 解码匹配 GPS 轨迹到道路网络。

        参数:
            gps_points: GPS 采样点序列。
            candidates: candidates[t] 是第 t 个 GPS 点的候选路段列表。
                每个候选是 {'segment_id': int, 'proj_point': (lat, lon)}。
            route_distances: 函数，输入两个路段 ID，返回路网最短距离。

        返回:
            匹配到的路段 ID 序列。
        """
        T = len(gps_points)

        # 初始化
        prev_delta = {}
        prev_psi = {}

        for i, cand in enumerate(candidates[0]):
            prob = self.emission_prob(gps_points[0], cand['proj_point'])
            prev_delta[i] = math.log(max(prob, 1e-300))
            prev_psi[i] = -1

        # 记录回溯路径
        all_psi = [{}]

        # 递推
        for t in range(1, T):
            curr_delta = {}
            curr_psi = {}
            gc_dist = self._haversine(gps_points[t - 1], gps_points[t])

            for j, cand_j in enumerate(candidates[t]):
                emit = self.emission_prob(gps_points[t], cand_j['proj_point'])
                log_emit = math.log(max(emit, 1e-300))

                best_score = -math.inf
                best_prev = 0

                for i, cand_i in enumerate(candidates[t - 1]):
                    if i not in prev_delta:
                        continue
                    rd = route_distances(cand_i['segment_id'],
                                         cand_j['segment_id'])
                    trans = self.transition_prob(rd, gc_dist)
                    score = prev_delta[i] + math.log(max(trans, 1e-300))

                    if score > best_score:
                        best_score = score
                        best_prev = i

                curr_delta[j] = best_score + log_emit
                curr_psi[j] = best_prev

            prev_delta = curr_delta
            all_psi.append(curr_psi)

        # 终止
        best_j = max(prev_delta, key=prev_delta.get)

        # 回溯
        path_indices = [0] * T
        path_indices[T - 1] = best_j
        for t in range(T - 2, -1, -1):
            path_indices[t] = all_psi[t + 1][path_indices[t + 1]]

        return [candidates[t][path_indices[t]]['segment_id']
                for t in range(T)]

    @staticmethod
    def _haversine(p1: Tuple[float, float],
                   p2: Tuple[float, float]) -> float:
        """Haversine 公式计算两点之间的距离（米）。"""
        R = 6371000.0
        lat1, lon1 = math.radians(p1[0]), math.radians(p1[1])
        lat2, lon2 = math.radians(p2[0]), math.radians(p2[1])
        dlat = lat2 - lat1
        dlon = lon2 - lon1
        a = (math.sin(dlat / 2) ** 2
             + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2) ** 2)
        return R * 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))


## 六、序列对齐：生物信息学的 Needleman-Wunsch / Smith-Waterman

### 6.1 全局对齐：Needleman-Wunsch

生物信息学中的序列对齐是 DP 最早的工业应用之一。给定两条蛋白质或 DNA 序列，我们需要找到它们之间的最佳对齐方式。

Needleman-Wunsch 算法求解**全局对齐**——将两条序列从头到尾完整对齐：

**状态**：$F(i, j)$ 表示将序列 $A$ 的前 $i$ 个字符与序列 $B$ 的前 $j$ 个字符对齐的最优得分。

**转移方程**：

$$
F(i, j) = \max \begin{cases}
F(i-1, j-1) + S(A_i, B_j) & \text{（匹配/替换）} \\
F(i-1, j) + d & \text{（A 中插入空位）} \\
F(i, j-1) + d & \text{（B 中插入空位）}
\end{cases}
$$

其中 $S(A_i, B_j)$ 是替换矩阵（如 BLOSUM62），$d$ 是空位罚分。

```python
def needleman_wunsch(
    seq_a: str,
    seq_b: str,
    match_score: int = 1,
    mismatch_score: int = -1,
    gap_penalty: int = -2
) -> Tuple[str, str, int]:
    """Needleman-Wunsch 全局序列对齐。

    返回:
        (aligned_a, aligned_b, score): 对齐后的序列和得分。
    """
    m, n = len(seq_a), len(seq_b)

    # 初始化得分矩阵
    F = [[0] * (n + 1) for _ in range(m + 1)]
    for i in range(m + 1):
        F[i][0] = i * gap_penalty
    for j in range(n + 1):
        F[0][j] = j * gap_penalty

    # 填表
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            if seq_a[i - 1] == seq_b[j - 1]:
                diag = F[i - 1][j - 1] + match_score
            else:
                diag = F[i - 1][j - 1] + mismatch_score
            up = F[i - 1][j] + gap_penalty
            left = F[i][j - 1] + gap_penalty
            F[i][j] = max(diag, up, left)

    # 回溯
    aligned_a, aligned_b = [], []
    i, j = m, n

    while i > 0 or j > 0:
        if i > 0 and j > 0:
            if seq_a[i - 1] == seq_b[j - 1]:
                score_diag = F[i - 1][j - 1] + match_score
            else:
                score_diag = F[i - 1][j - 1] + mismatch_score

        if i > 0 and j > 0 and F[i][j] == score_diag:
            aligned_a.append(seq_a[i - 1])
            aligned_b.append(seq_b[j - 1])
            i -= 1
            j -= 1
        elif i > 0 and F[i][j] == F[i - 1][j] + gap_penalty:
            aligned_a.append(seq_a[i - 1])
            aligned_b.append('-')
            i -= 1
        else:
            aligned_a.append('-')
            aligned_b.append(seq_b[j - 1])
            j -= 1

    return ''.join(reversed(aligned_a)), ''.join(reversed(aligned_b)), F[m][n]

def smith_waterman(
    seq_a: str,
    seq_b: str,
    match_score: int = 2,
    mismatch_score: int = -1,
    gap_penalty: int = -1
) -> Tuple[str, str, int]:
    """Smith-Waterman 局部序列对齐。"""
    m, n = len(seq_a), len(seq_b)

    F = [[0] * (n + 1) for _ in range(m + 1)]
    max_score = 0
    max_pos = (0, 0)

    for i in range(1, m + 1):
        for j in range(1, n + 1):
            if seq_a[i - 1] == seq_b[j - 1]:
                diag = F[i - 1][j - 1] + match_score
            else:
                diag = F[i - 1][j - 1] + mismatch_score
            up = F[i - 1][j] + gap_penalty
            left = F[i][j - 1] + gap_penalty
            F[i][j] = max(0, diag, up, left)

            if F[i][j] > max_score:
                max_score = F[i][j]
                max_pos = (i, j)

    # 回溯（直到得分为 0）
    aligned_a, aligned_b = [], []
    i, j = max_pos

    while i > 0 and j > 0 and F[i][j] > 0:
        if seq_a[i - 1] == seq_b[j - 1]:
            score_diag = F[i - 1][j - 1] + match_score
        else:
            score_diag = F[i - 1][j - 1] + mismatch_score

        if F[i][j] == score_diag:
            aligned_a.append(seq_a[i - 1])
            aligned_b.append(seq_b[j - 1])
            i -= 1
            j -= 1
        elif F[i][j] == F[i - 1][j] + gap_penalty:
            aligned_a.append(seq_a[i - 1])
            aligned_b.append('-')
            i -= 1
        else:
            aligned_a.append('-')
            aligned_b.append(seq_b[j - 1])
            j -= 1

    return ''.join(reversed(aligned_a)), ''.join(reversed(aligned_b)), max_score

class SimpleValueNetwork:
    """简单的值函数近似网络（线性函数近似）。

    实际系统中通常用深度神经网络替代。
    """

    def __init__(self, feature_dim: int, learning_rate: float = 0.01):
        self.weights = np.zeros(feature_dim)
        self.lr = learning_rate

    def predict(self, features: np.ndarray) -> float:
        """预测状态值。"""
        return np.dot(self.weights, features)

    def update(self, features: np.ndarray, target: float):
        """TD 更新。"""
        prediction = self.predict(features)
        td_error = target - prediction
        self.weights += self.lr * td_error * features


class QLearningAgent:
    """基于表格的 Q-learning 实现。"""

    def __init__(self, num_states: int, num_actions: int,
                 alpha: float = 0.1, gamma: float = 0.99,
                 epsilon: float = 0.1):
        self.Q = np.zeros((num_states, num_actions))
        self.alpha = alpha
        self.gamma = gamma
        self.epsilon = epsilon

    def choose_action(self, state: int) -> int:
        """epsilon-greedy 策略选择动作。"""
        if np.random.random() < self.epsilon:
            return np.random.randint(self.Q.shape[1])
        return int(np.argmax(self.Q[state]))

    def update(self, state: int, action: int, reward: float,
               next_state: int, done: bool):
        """Q-learning 更新。"""
        if done:
            target = reward
        else:
            target = reward + self.gamma * np.max(self.Q[next_state])

        self.Q[state, action] += self.alpha * (target - self.Q[state, action])

import random
from dataclasses import dataclass
from typing import List


@dataclass
class Impression:
    """一次展示机会。"""
    value: float   # 预估价值（CTR * bid）
    cost: float    # 花费（CPC）


def generate_traffic(num_periods: int, imps_per_period: int,
                     seed: int = 42) -> List[List[Impression]]:
    """生成日内流量，上午/晚间高峰，下午低谷。"""
    rng = random.Random(seed)
    traffic = []
    for t in range(num_periods):
        hour = (t * 24) // num_periods
        if 8 <= hour <= 11 or 19 <= hour <= 22:
            mult = 1.5
        elif 12 <= hour <= 17:
            mult = 0.8
        else:
            mult = 0.5
        n = int(imps_per_period * mult)
        period = []
        for _ in range(n):
            v = rng.expovariate(1.0) * 0.5
            c = v * rng.uniform(0.3, 0.8)
            period.append(Impression(value=v, cost=c))
        traffic.append(period)
    return traffic


def simulate_pacing(budget: float, traffic: List[List[Impression]],
                    ctrl: PacingPIDController) -> dict:
    """模拟完整 pacing 过程并返回统计信息。"""
    total_val, total_cost, total_imp, won = 0.0, 0.0, 0, 0
    for t, imps in enumerate(traffic):
        p_cost = 0.0
        for imp in imps:
            total_imp += 1
            if total_cost >= budget:
                break
            if ctrl.should_bid(imp.value, imp.cost):
                spend = min(imp.cost, budget - total_cost)
                total_val += imp.value
                total_cost += spend
                p_cost += spend
                won += 1
        ctrl.update(t + 1, p_cost)

    return {'value': total_val, 'cost': total_cost,
            'impressions': total_imp, 'won': won,
            'utilization': total_cost / budget}


# 运行演示
budget, periods = 100.0, 24
traffic = generate_traffic(periods, 200)
ctrl = PacingPIDController(budget, periods, kp=0.3, ki=0.05, kd=0.02)
result = simulate_pacing(budget, traffic, ctrl)
print(f"价值={result['value']:.2f}  花费={result['cost']:.2f}/{budget}"
      f"  利用率={result['utilization']:.1%}  胜率={result['won']/result['impressions']:.1%}")

def run_viterbi_demo():
    """天气-活动 HMM 解码演示。"""
    # 隐状态：0=晴天, 1=阴天, 2=雨天；观测：0=散步, 1=购物, 2=清洁
    hmm = HiddenMarkovModel(num_states=3, num_observations=3)
    hmm.set_params(
        pi=np.array([0.6, 0.3, 0.1]),
        A=np.array([[0.7, 0.2, 0.1],
                     [0.3, 0.4, 0.3],
                     [0.2, 0.3, 0.5]]),
        B=np.array([[0.6, 0.3, 0.1],   # 晴天
                     [0.3, 0.4, 0.3],   # 阴天
                     [0.1, 0.4, 0.5]])  # 雨天
    )
    observations = [0, 1, 2, 2, 1, 0, 0, 1, 2]
    names = ['晴天', '阴天', '雨天']
    path, lp = hmm.viterbi(observations)
    print("最优路径:", ' -> '.join(names[s] for s in path))
    print(f"对数概率: {lp:.4f}  概率: {np.exp(lp):.2e}")

run_viterbi_demo()