run_subset_sum_sweep.py

#!/usr/bin/env python3
"""Subset-sum ambiguity sweep: partition entropy vs. output type and count.

Sweeps a synthetic CoinJoin over three dimensions:
  - n: number of participants (2..8)
  - output_type: equal / standard-denomination / interval (channel-open)
  - fee_homogeneity: distinct vs. homogenized maker fees

For each cell, computes:
  - n_valid_partitions
  - partition_entropy_bits   (Shannon entropy over valid partitions)
  - partition_min_entropy_bits  (worst-case / min-entropy)
  - partition_bayes_vulnerability (max partition probability)

Also empirically recovers the Lagarias-Odlyzko density threshold by sweeping
the subset-sum density d = n / log2(max_value) and finding where
partition_min_entropy crosses 0 (unique solution).

Outputs subset_sum_sweep_results.json.
"""

from __future__ import annotations

import json
import math
import time
from pathlib import Path

from coinjoin_simulator.anonymity import (
    compute_change_anonymity,
    compute_change_anonymity_with_intervals,
    interval_valid_partition_count,
)
from coinjoin_simulator.models import (
    UTXO,
    AnonymityMetrics,
    CoinJoinTransaction,
    Participant,
    Role,
)

OUTPUT_PATH = Path("subset_sum_sweep_results.json")

# Standard Wasabi/WabiSabi denomination values (sats) for the denomination test.
_STD_DENOMS = (
    100_000, 131_072, 200_000, 262_144, 500_000, 524_288,
    1_000_000, 1_048_576, 2_000_000, 5_000_000,
)


def _make_tx(
    n_makers: int,
    *,
    cj_amount: int = 1_000_000,
    use_std_denom: bool = False,
    distinct_fees: bool = True,
    interval_changes: bool = False,
) -> CoinJoinTransaction:
    """Build a synthetic CoinJoin for the sweep."""
    participants: list[Participant] = []
    for i in range(n_makers):
        change = _nearest_std_denom(i * 80_000 + 50_000) if use_std_denom else (i + 1) * 80_000
        fee = 500 if distinct_fees else 250
        if distinct_fees:
            fee = 100 * (i + 1)
        p = Participant(
            role=Role.MAKER,
            entity_id=f"maker_{i}",
            utxos_in=[UTXO(value_sats=cj_amount + change + fee, owner_id=f"maker_{i}")],
            equal_output=UTXO(value_sats=cj_amount, owner_id=f"maker_{i}", is_equal_output=True),
            change_output=UTXO(value_sats=change, owner_id=f"maker_{i}", is_change=True),
            cj_fee_sats=fee,
        )
        participants.append(p)

    taker_fee = sum(p.cj_fee_sats for p in participants) + 2_000
    taker_change = 600_000
    taker = Participant(
        role=Role.TAKER,
        entity_id="taker",
        utxos_in=[
            UTXO(value_sats=cj_amount + taker_change + taker_fee, owner_id="taker")
        ],
        equal_output=UTXO(value_sats=cj_amount, owner_id="taker", is_equal_output=True),
        change_output=UTXO(value_sats=taker_change, owner_id="taker", is_change=True),
        cj_fee_sats=-taker_fee,
    )
    participants.append(taker)

    return CoinJoinTransaction(
        cj_amount=cj_amount,
        participants=participants,
        total_mining_fee=2_000,
    )


def _nearest_std_denom(amount: int) -> int:
    """Round amount to the nearest standard denomination (simulates WabiSabi)."""
    return min(_STD_DENOMS, key=lambda d: abs(d - amount))


def _summary(metrics: dict[str, AnonymityMetrics]) -> dict[str, float]:
    """Aggregate metrics across change outputs."""
    if not metrics:
        return {
            "mean_partition_entropy": 0.0,
            "min_partition_entropy": 0.0,
            "mean_partition_min_entropy": 0.0,
            "min_partition_min_entropy": 0.0,
            "mean_bayes_vulnerability": 1.0,
            "mean_n_valid_mappings": 1.0,
        }
    entropies = [m.partition_entropy_bits for m in metrics.values()]
    min_entropies = [m.partition_min_entropy_bits for m in metrics.values()]
    bvulns = [m.partition_bayes_vulnerability for m in metrics.values()]
    n_maps = [float(m.n_valid_mappings) for m in metrics.values()]
    return {
        "mean_partition_entropy": sum(entropies) / len(entropies),
        "min_partition_entropy": min(entropies),
        "mean_partition_min_entropy": sum(min_entropies) / len(min_entropies),
        "min_partition_min_entropy": min(min_entropies),
        "mean_bayes_vulnerability": sum(bvulns) / len(bvulns),
        "mean_n_valid_mappings": sum(n_maps) / len(n_maps),
    }


def run_output_type_sweep() -> list[dict[str, object]]:
    """Sweep partition entropy over (n_makers, output_type, fee_homogeneity)."""
    rows: list[dict[str, object]] = []
    for n in range(1, 8):
        for use_std_denom in (False, True):
            for distinct_fees in (True, False):
                output_type = "standard_denom" if use_std_denom else "exact_change"
                tx = _make_tx(n, use_std_denom=use_std_denom, distinct_fees=distinct_fees)
                metrics_exact = compute_change_anonymity(tx)
                s = _summary(metrics_exact)

                # Interval variant: treat all change outputs as interval-constrained
                # (simulating all-channel-open mini-CoinJoinXT round).
                interval_idx = frozenset(range(n + 1))  # n makers + taker
                metrics_interval = compute_change_anonymity_with_intervals(
                    tx, interval_change_indices=interval_idx
                )
                s_interval = _summary(metrics_interval)

                rows.append({
                    "n_participants": n + 1,  # makers + 1 taker
                    "output_type": output_type,
                    "distinct_fees": distinct_fees,
                    "exact": s,
                    "interval": s_interval,
                    "interval_entropy_gain": (
                        s_interval["mean_partition_entropy"] - s["mean_partition_entropy"]
                    ),
                    "interval_min_entropy_gain": (
                        s_interval["mean_partition_min_entropy"] - s["mean_partition_min_entropy"]
                    ),
                })
    return rows


def run_density_sweep() -> list[dict[str, object]]:
    """Sweep partition min-entropy vs Lagarias-Odlyzko subset-sum density.

    Density d = n / log2(V) where V is the max output value. For d < 0.94 the
    subset sum is solvable in polynomial time; for d >= 1.0 it becomes
    combinatorially hard. Here we measure where the transition appears in
    partition min-entropy.
    """
    rows: list[dict[str, object]] = []
    # Fix n=4 participants; vary V by scaling the CoinJoin amount.
    n_fixed = 3  # 3 makers + 1 taker = 4 participants
    for log2_v in range(15, 27):
        cj_amount = 2**log2_v
        density = (n_fixed + 1) / log2_v
        tx = _make_tx(n_fixed, cj_amount=cj_amount, distinct_fees=True)
        metrics = compute_change_anonymity(tx)
        s = _summary(metrics)
        rows.append({
            "n_participants": n_fixed + 1,
            "log2_max_value": log2_v,
            "density_d": round(density, 3),
            "cj_amount_sats": cj_amount,
            "mean_partition_min_entropy": s["mean_partition_min_entropy"],
            "mean_n_valid_mappings": s["mean_n_valid_mappings"],
        })
    return rows


def run_interval_ratio_sweep() -> list[dict[str, object]]:
    """Sweep exact-vs-interval partition count ratio over n_makers.

    For each n, computes:
      - exact partition count (normal CoinJoin, all exact amounts)
      - interval partition count (all channel-open outputs)
      - ratio = interval / exact  (the mini-CoinJoinXT partition-entropy gain)
    """
    rows: list[dict[str, object]] = []
    for n in range(1, 7):
        tx = _make_tx(n, distinct_fees=True)
        participants = tx.participants
        input_totals = [sum(u.value_sats for u in p.utxos_in) for p in participants]
        change_values = [
            p.change_output.value_sats if p.change_output else None
            for p in participants
        ]

        exact_count = interval_valid_partition_count(
            input_totals, change_values, tx.cj_amount, participants, tx.total_mining_fee,
            is_interval=[False] * len(change_values),
        )
        interval_count = interval_valid_partition_count(
            input_totals, change_values, tx.cj_amount, participants, tx.total_mining_fee,
            is_interval=[True] * len(change_values),
        )
        ratio = interval_count / max(1, exact_count)
        rows.append({
            "n_participants": n + 1,
            "exact_partition_count": exact_count,
            "interval_partition_count": interval_count,
            "interval_to_exact_ratio": round(ratio, 2),
            "exact_min_entropy": round(math.log2(max(1, exact_count)), 3),
            "interval_min_entropy": round(math.log2(max(1, interval_count)), 3),
        })
    return rows


def main() -> None:
    start = time.time()
    print("Running subset-sum ambiguity sweep...")

    output_type = run_output_type_sweep()
    density = run_density_sweep()
    interval_ratio = run_interval_ratio_sweep()

    payload: dict[str, object] = {
        "description": (
            "Partition entropy vs. output type, count, and interval constraint. "
            "Tests: exact/standard-denom/interval change outputs across n=2..8 "
            "participants; Lagarias-Odlyzko density threshold; "
            "interval-to-exact partition-count ratio (the mini-CoinJoinXT gain)."
        ),
        "output_type_sweep": output_type,
        "density_sweep": density,
        "interval_ratio_sweep": interval_ratio,
        "generated_at": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime()),
    }
    OUTPUT_PATH.write_text(json.dumps(payload, indent=2) + "\n")
    elapsed = time.time() - start
    print(f"\nWrote {OUTPUT_PATH} in {elapsed:.2f}s\n")

    print("Interval-to-exact partition ratio (mini-CoinJoinXT gain):")
    print(f"  {'n':>3} {'exact':>8} {'interval':>9} {'ratio':>7} {'H_inf gain':>11}")
    for r in interval_ratio:
        h_gain = r["interval_min_entropy"] - r["exact_min_entropy"]
        print(
            f"  {r['n_participants']:>3} {r['exact_partition_count']:>8} "
            f"{r['interval_partition_count']:>9} {r['interval_to_exact_ratio']:>7.1f}x "
            f"{h_gain:>10.2f} bits"
        )

    print("\nDensity sweep (n=4 fixed, varying cj_amount):")
    print(f"  {'log2(V)':>7} {'density d':>10} {'mean H_inf partition':>21}")
    for r in density:
        print(
            f"  {r['log2_max_value']:>7} {r['density_d']:>10.3f} "
            f"{r['mean_partition_min_entropy']:>21.3f}"
        )


if __name__ == "__main__":
    main()