沙子原地不动能零运费，为什么平滑算法还报了一个正数？

"""Tiny entropic-transport checks. KL convention: KL(pi || a tensor b)."""
import json
import math
from pathlib import Path
import platform


def cost_matrix(x, y):
    # Squared distance, with no factor 1/2.
    return [[(u - v) ** 2 for v in y] for u in x]


def objective(plan, a, b, costs, epsilon):
    transport = sum(plan[i][j] * costs[i][j] for i in range(len(a)) for j in range(len(b)))
    kl = sum(plan[i][j] * math.log(plan[i][j] / (a[i] * b[j]))
             for i in range(len(a)) for j in range(len(b)) if plan[i][j] > 0)
    return transport, kl, transport + epsilon * kl


def negative_entropy_objective(plan, costs, epsilon):
    transport = sum(plan[i][j] * costs[i][j] for i in range(len(plan)) for j in range(len(plan[i])))
    negative_entropy = sum(p * math.log(p) for row in plan for p in row if p > 0)
    return transport + epsilon * negative_entropy


def sinkhorn(a, b, costs, epsilon, tolerance=1e-13, max_iterations=20000):
    if not math.isfinite(epsilon) or epsilon <= 0:
        raise ValueError('positive finite epsilon required')
    if not a or not b or any(not math.isfinite(v) or v <= 0 for v in a + b):
        raise ValueError('strictly positive finite probability weights required')
    if abs(sum(a) - 1) > tolerance or abs(sum(b) - 1) > tolerance:
        raise ValueError('unit mass required')
    if len(costs) != len(a) or any(len(row) != len(b) for row in costs):
        raise ValueError('cost shape mismatch')
    if any(not math.isfinite(c) or c < 0 for row in costs for c in row):
        raise ValueError('finite nonnegative costs required')
    kernel = [[math.exp(-c / epsilon) for c in row] for row in costs]
    if any(v == 0 for row in kernel for v in row):
        raise ValueError('kernel underflow: use a log-domain solver outside this prototype')
    v = [1.0] * len(b)
    for iteration in range(1, max_iterations + 1):
        u = [a[i] / sum(kernel[i][j] * v[j] for j in range(len(b))) for i in range(len(a))]
        v = [b[j] / sum(kernel[i][j] * u[i] for i in range(len(a))) for j in range(len(b))]
        plan = [[u[i] * kernel[i][j] * v[j] for j in range(len(b))] for i in range(len(a))]
        residual = max([abs(sum(plan[i]) - a[i]) for i in range(len(a))]
                       + [abs(sum(plan[i][j] for i in range(len(a))) - b[j]) for j in range(len(b))])
        if residual < tolerance:
            return plan, residual, iteration
    raise RuntimeError('iteration cap reached')


def reference_2x2(a, b, costs, epsilon):
    """Independent 1D convex optimization via the monotone analytic derivative."""
    lo, hi = max(0.0, a[0] + b[0] - 1), min(a[0], b[0])
    delta = costs[0][0] - costs[0][1] - costs[1][0] + costs[1][1]
    for _ in range(100):
        t = (lo + hi) / 2
        cells = (t, a[0] - t, b[0] - t, 1 - a[0] - b[0] + t)
        derivative = delta + epsilon * (math.log(cells[0]) + math.log(cells[3])
                                          - math.log(cells[1]) - math.log(cells[2]))
        if derivative > 0:
            hi = t
        else:
            lo = t
        if hi - lo <= 2e-15:
            break
    t = (lo + hi) / 2
    return [[t, a[0] - t], [b[0] - t, 1 - a[0] - b[0] + t]]


def main():
    checks = 0
    def check(condition, label):
        nonlocal checks
        checks += 1
        if not condition:
            raise AssertionError(label)
    a = [0.5, 0.5]
    epsilon = 1 / math.log(3)
    costs = cost_matrix([0, 1], [0, 1])
    plan, residual, iterations = sinkhorn(a, a, costs, epsilon)
    expected = [[3 / 8, 1 / 8], [1 / 8, 3 / 8]]
    check(max(abs(plan[i][j] - expected[i][j]) for i in range(2) for j in range(2)) < 1e-14, 'closed-form plan')
    transport, kl, raw_self = objective(plan, a, a, costs, epsilon)
    closed_raw = math.log(1.5) / math.log(3)
    check(abs(raw_self - closed_raw) < 1e-13, 'closed-form regularized cost')
    check(abs(transport - 0.25) < 1e-13 and raw_self > 0, 'positive self transport and objective')
    identity = [[0.5, 0], [0, 0.5]]
    check(objective(identity, a, a, costs, epsilon)[0] == 0, 'unregularized self cost')
    check(raw_self < objective(identity, a, a, costs, epsilon)[2], 'splitting reduces regularized objective')
    collapsed_plan, _, _ = sinkhorn(a, [1.0], cost_matrix([0, 1], [0.5]), epsilon)
    raw_collapsed = objective(collapsed_plan, a, [1.0], cost_matrix([0, 1], [0.5]), epsilon)[2]
    debiased_collapsed = raw_collapsed - 0.5 * raw_self  # delta_0.5 self cost is exactly zero.
    check(abs(raw_collapsed - 0.25) < 1e-13 and raw_collapsed < raw_self, 'raw cost favors collapsed target in this case')
    check(debiased_collapsed > 0 and raw_self - 0.5 * raw_self - 0.5 * raw_self == 0, 'debiased comparison')

    convention_table = []
    for eps in (0.1, 0.3, epsilon, 1.5, 3.0):
        for label, bb, yy in [('same_two_piles', a, [0, 1]), ('one_pile_at_mean', [1.0], [0.5])]:
            cc = cost_matrix([0, 1], yy)
            pp, res, its = sinkhorn(a, bb, cc, eps)
            pure, relative_kl, full_kl = objective(pp, a, bb, cc, eps)
            negative_entropy_value = negative_entropy_objective(pp, cc, eps)
            marginal_entropy = -sum(w * math.log(w) for w in a + bb)
            check(abs((full_kl - negative_entropy_value) - eps * marginal_entropy) < 1e-12,
                  'KL versus negative entropy marginal constant')
            if label == 'same_two_piles':
                exact_pure = 1 / (1 + math.exp(1 / eps))
                exact_full = eps * (math.log(2) - math.log1p(math.exp(-1 / eps)))
                check(abs(pure - exact_pure) < 1e-12, 'epsilon sweep pure transport reference')
                check(abs(full_kl - exact_full) < 1e-12, 'epsilon sweep KL objective reference')
            convention_table.append({'epsilon': eps, 'target': label, 'pure_transport': pure,
                                     'KL_pi_vs_a_tensor_b': relative_kl, 'full_KL_objective': full_kl,
                                     'negative_entropy_objective': negative_entropy_value,
                                     'KL_minus_negative_entropy_constant': eps * marginal_entropy})
    special = [row for row in convention_table if row['epsilon'] == epsilon]
    check(special[0]['full_KL_objective'] > special[1]['full_KL_objective'], 'KL objective comparison across different target marginals')
    check(special[0]['negative_entropy_objective'] < special[1]['negative_entropy_objective'], 'negative entropy comparison can reverse across target marginals')

    comparisons = 0
    worst_plan_error = 0.0
    worst_residual = 0.0
    max_used_iterations = 0
    for mass_a in (0.2, 0.35, 0.5, 0.8):
        for mass_b in (0.2, 0.4, 0.5, 0.7):
            aa, bb = [mass_a, 1 - mass_a], [mass_b, 1 - mass_b]
            for target in ([0.0, 1.0], [0.2, 1.4], [-0.3, 0.7]):
                cc = cost_matrix([0, 1], target)
                for eps in (0.3, 0.7, 1.5):
                    pp, res, its = sinkhorn(aa, bb, cc, eps)
                    ref = reference_2x2(aa, bb, cc, eps)
                    err = max(abs(pp[i][j] - ref[i][j]) for i in range(2) for j in range(2))
                    check(err < 5e-12, 'independent 1D reference')
                    check(abs(objective(pp, aa, bb, cc, eps)[2] - objective(ref, aa, bb, cc, eps)[2]) < 1e-11, 'reference objective')
                    check(res < 1e-13, 'marginal residual')
                    worst_plan_error = max(worst_plan_error, err)
                    worst_residual = max(worst_residual, res)
                    max_used_iterations = max(max_used_iterations, its)
                    comparisons += 1
    rejected = 0
    for aa, bb, cc, eps in [(a, a, costs, 0), ([0, 1], a, costs, 1),
                             ([0.6, 0.6], a, costs, 1), (a, a, [[0]], 1),
                             (a, a, costs, 1e-6), (a, a, costs, float('nan'))]:
        try:
            sinkhorn(aa, bb, cc, eps)
        except ValueError:
            rejected += 1
    check(rejected == 6, 'invalid inputs and underflow rejected')
    report = {'status': 'PASS', 'runtime': platform.python_version(), 'checks': checks,
              'independent_2x2_comparisons': comparisons, 'worst_plan_error': worst_plan_error,
              'worst_marginal_residual': worst_residual, 'maximum_iterations_used': max_used_iterations,
              'epsilon': epsilon, 'self_plan': plan, 'self_transport': transport, 'self_KL': kl,
              'self_regularized_objective': raw_self, 'collapsed_raw_objective': raw_collapsed,
              'self_sinkhorn_divergence': 0.0, 'collapsed_sinkhorn_divergence': debiased_collapsed,
              'epsilon_and_convention_table': convention_table,
              'convention_identity': 'J_KL = J_negative_entropy + epsilon*(H(a)+H(b)); the added constant is fixed only when both marginals are fixed.',
              'conventions': 'C(x,y)=(x-y)^2; no factor 1/2. OT_eps=min <C,pi>+eps*KL(pi||a tensor b).',
              'boundary': 'Only tiny strictly-positive finite discrete distributions. No log-domain solver, GPU, GeomLoss run, training, or model-quality claim.'}
    Path('report.json').write_text(json.dumps(report, indent=2) + '\n', encoding='utf-8')
    print('epsilon = 1/log(3); self plan = [[3/8,1/8],[1/8,3/8]]')
    print(f'self: transport={transport:.9f}; KL={kl:.9f}; full objective={raw_self:.9f}')
    print(f'collapsed mean: raw={raw_collapsed:.9f}; debiased={debiased_collapsed:.9f}; debiased self=0')
    print('epsilon | target | pure transport | full KL objective | negative entropy objective')
    for row in convention_table:
        print(f"{row['epsilon']:.9f} | {row['target']} | {row['pure_transport']:.9f} | {row['full_KL_objective']:.9f} | {row['negative_entropy_objective']:.9f}")
    print(f'PASS: {checks} checks; {comparisons} independent 2x2 comparisons; {rejected} rejected inputs')


if __name__ == '__main__':
    main()