ferro-ta/benchmarks/bench_derivatives_compare.py

"""
Compare ferro_ta derivatives analytics against analytical references and
optional third-party libraries available in the current environment.

The suite focuses on selected core workflows:

- Black-Scholes-Merton call pricing
- implied-volatility recovery
- first-order Greeks
- Black-76 call pricing

Outputs include:

- speed timings with per-run samples and variance stats
- analytical accuracy metrics
- Python-tracked peak allocation snapshots
- machine, runtime, build, and package metadata
"""

from __future__ import annotations

import argparse
import importlib.util
import json
import math
import sys
import time
import tracemalloc
from collections.abc import Callable
from dataclasses import dataclass
from pathlib import Path
from typing import Any

import numpy as np

from ferro_ta.analysis.options import (
    black_76_price as ft_black_76_price,
)
from ferro_ta.analysis.options import (
    greeks as ft_greeks,
)
from ferro_ta.analysis.options import (
    implied_volatility as ft_implied_volatility,
)
from ferro_ta.analysis.options import (
    option_price as ft_option_price,
)

try:
    from benchmarks.metadata import benchmark_metadata, package_versions
except ModuleNotFoundError:  # pragma: no cover - script execution fallback
    from metadata import benchmark_metadata, package_versions


N_WARMUP = 1
N_RUNS = 7
DEFAULT_SIZES = [1_000, 10_000]
DEFAULT_ACCURACY_SIZE = 512
DEFAULT_SEED = 42
INV_SQRT_2PI = 1.0 / math.sqrt(2.0 * math.pi)


@dataclass(frozen=True)
class Case:
    name: str
    label: str
    expected_key: str
    component_names: tuple[str, ...] | None = None
    accuracy_target: str = "expected_output"


@dataclass(frozen=True)
class Provider:
    name: str
    kind: str
    note: str
    functions: dict[str, Callable[[dict[str, np.ndarray]], np.ndarray]]
    max_speed_size: int | None = None

    def supports(self, case_name: str) -> bool:
        return case_name in self.functions


CASES = [
    Case("bsm_call_price", "BSM Call Price", "call_price"),
    Case(
        "bsm_call_iv",
        "BSM Call IV Recovery",
        "volatility",
        accuracy_target="reconstructed_price",
    ),
    Case(
        "bsm_call_greeks",
        "BSM Call Greeks",
        "call_greeks",
        component_names=("delta", "gamma", "vega", "theta", "rho"),
    ),
    Case("black76_call_price", "Black-76 Call Price", "black76_call_price"),
]


def _median(values: list[float]) -> float:
    ordered = sorted(values)
    mid = len(ordered) // 2
    if len(ordered) % 2:
        return ordered[mid]
    return (ordered[mid - 1] + ordered[mid]) / 2.0


def _summary_stats(samples_ms: list[float]) -> dict[str, float]:
    if not samples_ms:
        return {
            "median_ms": 0.0,
            "mean_ms": 0.0,
            "min_ms": 0.0,
            "max_ms": 0.0,
            "stddev_ms": 0.0,
            "cv_pct": 0.0,
        }

    mean_ms = sum(samples_ms) / len(samples_ms)
    variance = (
        sum((sample - mean_ms) ** 2 for sample in samples_ms) / (len(samples_ms) - 1)
        if len(samples_ms) > 1
        else 0.0
    )
    stddev_ms = math.sqrt(variance)
    cv_pct = (stddev_ms / mean_ms * 100.0) if mean_ms else 0.0
    return {
        "median_ms": round(_median(samples_ms), 4),
        "mean_ms": round(mean_ms, 4),
        "min_ms": round(min(samples_ms), 4),
        "max_ms": round(max(samples_ms), 4),
        "stddev_ms": round(stddev_ms, 4),
        "cv_pct": round(cv_pct, 3),
    }


def _timed_runs_ms(
    fn: Callable[[dict[str, np.ndarray]], np.ndarray],
    chain: dict[str, np.ndarray],
) -> list[float]:
    for _ in range(N_WARMUP):
        fn(chain)

    samples_ms: list[float] = []
    for _ in range(N_RUNS):
        t0 = time.perf_counter()
        fn(chain)
        samples_ms.append((time.perf_counter() - t0) * 1000.0)
    return samples_ms


def _python_peak_bytes(
    fn: Callable[[dict[str, np.ndarray]], np.ndarray],
    chain: dict[str, np.ndarray],
) -> int | None:
    try:
        tracemalloc.start()
        tracemalloc.reset_peak()
        fn(chain)
        _, peak = tracemalloc.get_traced_memory()
        return int(peak)
    except Exception:
        return None
    finally:
        tracemalloc.stop()


def _throughput_contracts_s(size: int, median_ms: float) -> float:
    if median_ms <= 0:
        return 0.0
    return size / (median_ms / 1000.0)


def _normal_pdf_numpy(x: np.ndarray) -> np.ndarray:
    return INV_SQRT_2PI * np.exp(-0.5 * x * x)


def _normal_cdf_numpy(x: np.ndarray) -> np.ndarray:
    abs_x = np.abs(x)
    t = 1.0 / (1.0 + 0.2316419 * abs_x)
    poly = (
        (((((1.330274429 * t) - 1.821255978) * t) + 1.781477937) * t - 0.356563782) * t
        + 0.319381530
    ) * t
    cdf = 1.0 - _normal_pdf_numpy(abs_x) * poly
    return np.where(x >= 0.0, cdf, 1.0 - cdf)


def _normal_cdf_scalar(x: float) -> float:
    return 0.5 * (1.0 + math.erf(x / math.sqrt(2.0)))


def _normal_pdf_scalar(x: float) -> float:
    return INV_SQRT_2PI * math.exp(-0.5 * x * x)


def _bsm_price_numpy(
    spot: np.ndarray,
    strike: np.ndarray,
    rate: np.ndarray,
    time_to_expiry: np.ndarray,
    volatility: np.ndarray,
    *,
    option_type: str,
    carry: np.ndarray,
) -> np.ndarray:
    sqrt_t = np.sqrt(time_to_expiry)
    sigma_sqrt_t = volatility * sqrt_t
    d1 = (
        np.log(spot / strike)
        + (rate - carry + 0.5 * volatility * volatility) * time_to_expiry
    ) / sigma_sqrt_t
    d2 = d1 - sigma_sqrt_t
    spot_df = np.exp(-carry * time_to_expiry)
    strike_df = np.exp(-rate * time_to_expiry)
    if option_type == "call":
        out = spot * spot_df * _normal_cdf_numpy(
            d1
        ) - strike * strike_df * _normal_cdf_numpy(d2)
    else:
        out = strike * strike_df * _normal_cdf_numpy(
            -d2
        ) - spot * spot_df * _normal_cdf_numpy(-d1)
    return np.ascontiguousarray(out, dtype=np.float64)


def _black76_price_numpy(
    forward: np.ndarray,
    strike: np.ndarray,
    rate: np.ndarray,
    time_to_expiry: np.ndarray,
    volatility: np.ndarray,
    *,
    option_type: str,
) -> np.ndarray:
    sqrt_t = np.sqrt(time_to_expiry)
    sigma_sqrt_t = volatility * sqrt_t
    d1 = (
        np.log(forward / strike) + 0.5 * volatility * volatility * time_to_expiry
    ) / sigma_sqrt_t
    d2 = d1 - sigma_sqrt_t
    discount = np.exp(-rate * time_to_expiry)
    if option_type == "call":
        out = discount * (
            forward * _normal_cdf_numpy(d1) - strike * _normal_cdf_numpy(d2)
        )
    else:
        out = discount * (
            strike * _normal_cdf_numpy(-d2) - forward * _normal_cdf_numpy(-d1)
        )
    return np.ascontiguousarray(out, dtype=np.float64)


def _bsm_greeks_numpy(
    spot: np.ndarray,
    strike: np.ndarray,
    rate: np.ndarray,
    time_to_expiry: np.ndarray,
    volatility: np.ndarray,
    *,
    option_type: str,
    carry: np.ndarray,
) -> np.ndarray:
    sqrt_t = np.sqrt(time_to_expiry)
    sigma_sqrt_t = volatility * sqrt_t
    d1 = (
        np.log(spot / strike)
        + (rate - carry + 0.5 * volatility * volatility) * time_to_expiry
    ) / sigma_sqrt_t
    d2 = d1 - sigma_sqrt_t
    pdf = _normal_pdf_numpy(d1)
    carry_df = np.exp(-carry * time_to_expiry)
    strike_df = np.exp(-rate * time_to_expiry)

    if option_type == "call":
        delta = carry_df * _normal_cdf_numpy(d1)
        theta = (
            -(spot * carry_df * pdf * volatility) / (2.0 * sqrt_t)
            - rate * strike * strike_df * _normal_cdf_numpy(d2)
            + carry * spot * carry_df * _normal_cdf_numpy(d1)
        )
        rho = strike * time_to_expiry * strike_df * _normal_cdf_numpy(d2)
    else:
        delta = carry_df * (_normal_cdf_numpy(d1) - 1.0)
        theta = (
            -(spot * carry_df * pdf * volatility) / (2.0 * sqrt_t)
            + rate * strike * strike_df * _normal_cdf_numpy(-d2)
            - carry * spot * carry_df * _normal_cdf_numpy(-d1)
        )
        rho = -strike * time_to_expiry * strike_df * _normal_cdf_numpy(-d2)

    gamma = carry_df * pdf / (spot * sigma_sqrt_t)
    vega = spot * carry_df * pdf * sqrt_t
    return np.ascontiguousarray(
        np.column_stack([delta, gamma, vega, theta, rho]),
        dtype=np.float64,
    )


def _bsm_price_scalar(
    spot: float,
    strike: float,
    rate: float,
    time_to_expiry: float,
    volatility: float,
    *,
    option_type: str,
    carry: float,
) -> float:
    sqrt_t = math.sqrt(time_to_expiry)
    sigma_sqrt_t = volatility * sqrt_t
    d1 = (
        math.log(spot / strike)
        + (rate - carry + 0.5 * volatility * volatility) * time_to_expiry
    ) / sigma_sqrt_t
    d2 = d1 - sigma_sqrt_t
    spot_df = math.exp(-carry * time_to_expiry)
    strike_df = math.exp(-rate * time_to_expiry)
    if option_type == "call":
        return spot * spot_df * _normal_cdf_scalar(
            d1
        ) - strike * strike_df * _normal_cdf_scalar(d2)
    return strike * strike_df * _normal_cdf_scalar(
        -d2
    ) - spot * spot_df * _normal_cdf_scalar(-d1)


def _black76_price_scalar(
    forward: float,
    strike: float,
    rate: float,
    time_to_expiry: float,
    volatility: float,
    *,
    option_type: str,
) -> float:
    sqrt_t = math.sqrt(time_to_expiry)
    sigma_sqrt_t = volatility * sqrt_t
    d1 = (
        math.log(forward / strike) + 0.5 * volatility * volatility * time_to_expiry
    ) / sigma_sqrt_t
    d2 = d1 - sigma_sqrt_t
    discount = math.exp(-rate * time_to_expiry)
    if option_type == "call":
        return discount * (
            forward * _normal_cdf_scalar(d1) - strike * _normal_cdf_scalar(d2)
        )
    return discount * (
        strike * _normal_cdf_scalar(-d2) - forward * _normal_cdf_scalar(-d1)
    )


def _bsm_greeks_scalar(
    spot: float,
    strike: float,
    rate: float,
    time_to_expiry: float,
    volatility: float,
    *,
    option_type: str,
    carry: float,
) -> tuple[float, float, float, float, float]:
    sqrt_t = math.sqrt(time_to_expiry)
    sigma_sqrt_t = volatility * sqrt_t
    d1 = (
        math.log(spot / strike)
        + (rate - carry + 0.5 * volatility * volatility) * time_to_expiry
    ) / sigma_sqrt_t
    d2 = d1 - sigma_sqrt_t
    pdf = _normal_pdf_scalar(d1)
    carry_df = math.exp(-carry * time_to_expiry)
    strike_df = math.exp(-rate * time_to_expiry)

    if option_type == "call":
        delta = carry_df * _normal_cdf_scalar(d1)
        theta = (
            -(spot * carry_df * pdf * volatility) / (2.0 * sqrt_t)
            - rate * strike * strike_df * _normal_cdf_scalar(d2)
            + carry * spot * carry_df * _normal_cdf_scalar(d1)
        )
        rho = strike * time_to_expiry * strike_df * _normal_cdf_scalar(d2)
    else:
        delta = carry_df * (_normal_cdf_scalar(d1) - 1.0)
        theta = (
            -(spot * carry_df * pdf * volatility) / (2.0 * sqrt_t)
            + rate * strike * strike_df * _normal_cdf_scalar(-d2)
            - carry * spot * carry_df * _normal_cdf_scalar(-d1)
        )
        rho = -strike * time_to_expiry * strike_df * _normal_cdf_scalar(-d2)

    gamma = carry_df * pdf / (spot * sigma_sqrt_t)
    vega = spot * carry_df * pdf * sqrt_t
    return delta, gamma, vega, theta, rho


def _implied_vol_bisection_numpy(
    price: np.ndarray,
    spot: np.ndarray,
    strike: np.ndarray,
    rate: np.ndarray,
    time_to_expiry: np.ndarray,
    *,
    option_type: str,
    carry: np.ndarray,
    lower: float = 1e-6,
    upper: float = 5.0,
    tolerance: float = 1e-8,
    max_iterations: int = 100,
) -> np.ndarray:
    lo = np.full_like(price, lower, dtype=np.float64)
    hi = np.full_like(price, upper, dtype=np.float64)
    mid = np.full_like(price, 0.2, dtype=np.float64)
    for _ in range(max_iterations):
        mid = (lo + hi) / 2.0
        estimate = _bsm_price_numpy(
            spot,
            strike,
            rate,
            time_to_expiry,
            mid,
            option_type=option_type,
            carry=carry,
        )
        too_low = estimate < price
        lo = np.where(too_low, mid, lo)
        hi = np.where(too_low, hi, mid)
        if float(np.max(np.abs(estimate - price))) < tolerance:
            break
    return np.ascontiguousarray(mid, dtype=np.float64)


def _implied_vol_bisection_scalar(
    price: float,
    spot: float,
    strike: float,
    rate: float,
    time_to_expiry: float,
    *,
    option_type: str,
    carry: float,
    lower: float = 1e-6,
    upper: float = 5.0,
    tolerance: float = 1e-10,
    max_iterations: int = 100,
) -> float:
    lo = lower
    hi = upper
    mid = 0.2
    for _ in range(max_iterations):
        mid = (lo + hi) / 2.0
        estimate = _bsm_price_scalar(
            spot,
            strike,
            rate,
            time_to_expiry,
            mid,
            option_type=option_type,
            carry=carry,
        )
        if abs(estimate - price) < tolerance:
            return mid
        if estimate < price:
            lo = mid
        else:
            hi = mid
    return mid


def _reference_python_loop(
    chain: dict[str, np.ndarray],
    fn: Callable[..., float | tuple[float, ...]],
    *,
    include_forward: bool = False,
    include_price: bool = False,
) -> np.ndarray:
    rows: list[Any] = []
    for idx in range(len(chain["strike"])):
        kwargs: dict[str, float | str] = {
            "strike": float(chain["strike"][idx]),
            "rate": float(chain["rate"][idx]),
            "time_to_expiry": float(chain["time_to_expiry"][idx]),
            "volatility": float(chain["volatility"][idx]),
            "carry": float(chain["carry"][idx]),
        }
        if include_forward:
            kwargs["forward"] = float(chain["forward"][idx])
        else:
            kwargs["spot"] = float(chain["spot"][idx])
        if include_price:
            kwargs["price"] = float(chain["call_price"][idx])
        rows.append(fn(**kwargs))
    return np.asarray(rows, dtype=np.float64)


def _build_chain(n: int, *, seed: int) -> dict[str, np.ndarray]:
    rng = np.random.default_rng(seed)
    spot = np.ascontiguousarray(rng.uniform(80.0, 120.0, size=n), dtype=np.float64)
    strike = np.ascontiguousarray(rng.uniform(70.0, 130.0, size=n), dtype=np.float64)
    rate = np.ascontiguousarray(rng.uniform(0.0, 0.07, size=n), dtype=np.float64)
    carry = np.ascontiguousarray(rng.uniform(0.0, 0.03, size=n), dtype=np.float64)
    time_to_expiry = np.ascontiguousarray(
        rng.uniform(7.0 / 365.0, 2.0, size=n), dtype=np.float64
    )
    volatility = np.ascontiguousarray(rng.uniform(0.08, 0.65, size=n), dtype=np.float64)
    forward = np.ascontiguousarray(
        spot * np.exp((rate - carry) * time_to_expiry),
        dtype=np.float64,
    )

    chain = {
        "spot": spot,
        "strike": strike,
        "rate": rate,
        "carry": carry,
        "time_to_expiry": time_to_expiry,
        "volatility": volatility,
        "forward": forward,
    }

    call_price = _reference_python_loop(
        chain,
        lambda *, spot, strike, rate, time_to_expiry, volatility, carry: (
            _bsm_price_scalar(
                spot,
                strike,
                rate,
                time_to_expiry,
                volatility,
                option_type="call",
                carry=carry,
            )
        ),
    )
    call_greeks = _reference_python_loop(
        chain,
        lambda *, spot, strike, rate, time_to_expiry, volatility, carry: (
            _bsm_greeks_scalar(
                spot,
                strike,
                rate,
                time_to_expiry,
                volatility,
                option_type="call",
                carry=carry,
            )
        ),
    )
    black76_call_price = _reference_python_loop(
        chain,
        lambda *, forward, strike, rate, time_to_expiry, volatility, carry: (
            _black76_price_scalar(
                forward,
                strike,
                rate,
                time_to_expiry,
                volatility,
                option_type="call",
            )
        ),
        include_forward=True,
    )

    chain["call_price"] = np.ascontiguousarray(call_price, dtype=np.float64)
    chain["call_greeks"] = np.ascontiguousarray(call_greeks, dtype=np.float64)
    chain["black76_call_price"] = np.ascontiguousarray(
        black76_call_price, dtype=np.float64
    )
    return chain


def _ferro_ta_call_price(chain: dict[str, np.ndarray]) -> np.ndarray:
    return np.asarray(
        ft_option_price(
            chain["spot"],
            chain["strike"],
            chain["rate"],
            chain["time_to_expiry"],
            chain["volatility"],
            option_type="call",
            model="bsm",
            carry=chain["carry"],
        ),
        dtype=np.float64,
    )


def _ferro_ta_call_iv(chain: dict[str, np.ndarray]) -> np.ndarray:
    return np.asarray(
        ft_implied_volatility(
            chain["call_price"],
            chain["spot"],
            chain["strike"],
            chain["rate"],
            chain["time_to_expiry"],
            option_type="call",
            model="bsm",
            carry=chain["carry"],
        ),
        dtype=np.float64,
    )


def _ferro_ta_call_greeks(chain: dict[str, np.ndarray]) -> np.ndarray:
    result = ft_greeks(
        chain["spot"],
        chain["strike"],
        chain["rate"],
        chain["time_to_expiry"],
        chain["volatility"],
        option_type="call",
        model="bsm",
        carry=chain["carry"],
    )
    return np.ascontiguousarray(
        np.column_stack(
            [
                np.asarray(result.delta, dtype=np.float64),
                np.asarray(result.gamma, dtype=np.float64),
                np.asarray(result.vega, dtype=np.float64),
                np.asarray(result.theta, dtype=np.float64),
                np.asarray(result.rho, dtype=np.float64),
            ]
        ),
        dtype=np.float64,
    )


def _ferro_ta_black76_call_price(chain: dict[str, np.ndarray]) -> np.ndarray:
    return np.asarray(
        ft_black_76_price(
            chain["forward"],
            chain["strike"],
            chain["rate"],
            chain["time_to_expiry"],
            chain["volatility"],
            option_type="call",
        ),
        dtype=np.float64,
    )


def _reference_numpy_call_price(chain: dict[str, np.ndarray]) -> np.ndarray:
    return _bsm_price_numpy(
        chain["spot"],
        chain["strike"],
        chain["rate"],
        chain["time_to_expiry"],
        chain["volatility"],
        option_type="call",
        carry=chain["carry"],
    )


def _reference_numpy_call_iv(chain: dict[str, np.ndarray]) -> np.ndarray:
    return _implied_vol_bisection_numpy(
        chain["call_price"],
        chain["spot"],
        chain["strike"],
        chain["rate"],
        chain["time_to_expiry"],
        option_type="call",
        carry=chain["carry"],
    )


def _reference_numpy_call_greeks(chain: dict[str, np.ndarray]) -> np.ndarray:
    return _bsm_greeks_numpy(
        chain["spot"],
        chain["strike"],
        chain["rate"],
        chain["time_to_expiry"],
        chain["volatility"],
        option_type="call",
        carry=chain["carry"],
    )


def _reference_numpy_black76_call_price(chain: dict[str, np.ndarray]) -> np.ndarray:
    return _black76_price_numpy(
        chain["forward"],
        chain["strike"],
        chain["rate"],
        chain["time_to_expiry"],
        chain["volatility"],
        option_type="call",
    )


def _reference_python_call_price(chain: dict[str, np.ndarray]) -> np.ndarray:
    return _reference_python_loop(
        chain,
        lambda *, spot, strike, rate, time_to_expiry, volatility, carry: (
            _bsm_price_scalar(
                spot,
                strike,
                rate,
                time_to_expiry,
                volatility,
                option_type="call",
                carry=carry,
            )
        ),
    )


def _reference_python_call_iv(chain: dict[str, np.ndarray]) -> np.ndarray:
    return _reference_python_loop(
        chain,
        lambda *, price, spot, strike, rate, time_to_expiry, volatility, carry: (
            _implied_vol_bisection_scalar(
                price,
                spot,
                strike,
                rate,
                time_to_expiry,
                option_type="call",
                carry=carry,
            )
        ),
        include_price=True,
    )


def _reference_python_call_greeks(chain: dict[str, np.ndarray]) -> np.ndarray:
    return _reference_python_loop(
        chain,
        lambda *, spot, strike, rate, time_to_expiry, volatility, carry: (
            _bsm_greeks_scalar(
                spot,
                strike,
                rate,
                time_to_expiry,
                volatility,
                option_type="call",
                carry=carry,
            )
        ),
    )


def _reference_python_black76_call_price(chain: dict[str, np.ndarray]) -> np.ndarray:
    return _reference_python_loop(
        chain,
        lambda *, forward, strike, rate, time_to_expiry, volatility, carry: (
            _black76_price_scalar(
                forward,
                strike,
                rate,
                time_to_expiry,
                volatility,
                option_type="call",
            )
        ),
        include_forward=True,
    )


def _reprice_bsm_call_from_iv(
    chain: dict[str, np.ndarray],
    implied_vols: np.ndarray,
) -> np.ndarray:
    rows = [
        _bsm_price_scalar(
            float(spot),
            float(strike),
            float(rate),
            float(time_to_expiry),
            max(float(iv), 1e-12),
            option_type="call",
            carry=float(carry),
        )
        for spot, strike, rate, time_to_expiry, iv, carry in zip(
            chain["spot"],
            chain["strike"],
            chain["rate"],
            chain["time_to_expiry"],
            np.asarray(implied_vols, dtype=np.float64),
            chain["carry"],
        )
    ]
    return np.asarray(rows, dtype=np.float64)


def _py_vollib_provider() -> Provider | None:
    if importlib.util.find_spec("py_vollib") is None:
        return None

    from py_vollib.black_scholes_merton import black_scholes_merton as py_vollib_bsm
    from py_vollib.black_scholes_merton.implied_volatility import (
        implied_volatility as py_vollib_iv,
    )

    def _price(chain: dict[str, np.ndarray]) -> np.ndarray:
        return np.asarray(
            [
                py_vollib_bsm(
                    "c",
                    float(s),
                    float(k),
                    float(t),
                    float(r),
                    float(vol),
                    float(q),
                )
                for s, k, r, t, vol, q in zip(
                    chain["spot"],
                    chain["strike"],
                    chain["rate"],
                    chain["time_to_expiry"],
                    chain["volatility"],
                    chain["carry"],
                )
            ],
            dtype=np.float64,
        )

    def _iv(chain: dict[str, np.ndarray]) -> np.ndarray:
        return np.asarray(
            [
                py_vollib_iv(
                    float(price),
                    "c",
                    float(s),
                    float(k),
                    float(t),
                    float(r),
                    float(q),
                )
                for price, s, k, r, t, q in zip(
                    chain["call_price"],
                    chain["spot"],
                    chain["strike"],
                    chain["rate"],
                    chain["time_to_expiry"],
                    chain["carry"],
                )
            ],
            dtype=np.float64,
        )

    return Provider(
        name="py_vollib",
        kind="third_party",
        note="Scalar Black-Scholes-Merton baseline from py_vollib.",
        functions={
            "bsm_call_price": _price,
            "bsm_call_iv": _iv,
        },
        max_speed_size=1_000,
    )


def available_providers() -> list[Provider]:
    providers = [
        Provider(
            name="ferro_ta",
            kind="project",
            note="Rust-backed vectorized implementation.",
            functions={
                "bsm_call_price": _ferro_ta_call_price,
                "bsm_call_iv": _ferro_ta_call_iv,
                "bsm_call_greeks": _ferro_ta_call_greeks,
                "black76_call_price": _ferro_ta_black76_call_price,
            },
        ),
        Provider(
            name="reference_numpy",
            kind="reference",
            note="Pure NumPy analytical formulas with vectorized IV bisection.",
            functions={
                "bsm_call_price": _reference_numpy_call_price,
                "bsm_call_iv": _reference_numpy_call_iv,
                "bsm_call_greeks": _reference_numpy_call_greeks,
                "black76_call_price": _reference_numpy_black76_call_price,
            },
        ),
        Provider(
            name="reference_python_loop",
            kind="reference",
            note="Scalar math-loop analytical baseline; useful for accuracy sanity checks.",
            functions={
                "bsm_call_price": _reference_python_call_price,
                "bsm_call_iv": _reference_python_call_iv,
                "bsm_call_greeks": _reference_python_call_greeks,
                "black76_call_price": _reference_python_black76_call_price,
            },
            max_speed_size=1_000,
        ),
    ]
    optional = _py_vollib_provider()
    if optional is not None:
        providers.append(optional)
    return providers


def _accuracy_metrics(actual: np.ndarray, expected: np.ndarray) -> dict[str, Any]:
    actual_arr = np.asarray(actual, dtype=np.float64)
    expected_arr = np.asarray(expected, dtype=np.float64)
    abs_error = np.abs(actual_arr - expected_arr)
    rel_error = abs_error / np.maximum(np.abs(expected_arr), 1e-12)
    return {
        "output_shape": list(actual_arr.shape),
        "max_abs_error": round(float(np.max(abs_error)), 12),
        "mean_abs_error": round(float(np.mean(abs_error)), 12),
        "rmse": round(
            float(np.sqrt(np.mean(np.square(actual_arr - expected_arr)))), 12
        ),
        "max_rel_error": round(float(np.max(rel_error)), 12),
    }


def _accuracy_summary(rows: list[dict[str, Any]]) -> list[dict[str, Any]]:
    summary: list[dict[str, Any]] = []
    for case in CASES:
        case_rows = [
            row for row in rows if row["case"] == case.name and "max_abs_error" in row
        ]
        if not case_rows:
            continue
        best = min(case_rows, key=lambda row: float(row["max_abs_error"]))
        worst = max(case_rows, key=lambda row: float(row["max_abs_error"]))
        summary.append(
            {
                "case": case.name,
                "best_provider": best["provider"],
                "best_max_abs_error": best["max_abs_error"],
                "worst_provider": worst["provider"],
                "worst_max_abs_error": worst["max_abs_error"],
            }
        )
    return summary


def _speed_summary(rows: list[dict[str, Any]]) -> list[dict[str, Any]]:
    summary: list[dict[str, Any]] = []
    for case in CASES:
        for size in sorted(
            {
                int(row["size"])
                for row in rows
                if row["case"] == case.name and "median_ms" in row
            }
        ):
            case_rows = [
                row
                for row in rows
                if row["case"] == case.name
                and row.get("size") == size
                and "median_ms" in row
            ]
            if not case_rows:
                continue
            fastest = min(case_rows, key=lambda row: float(row["median_ms"]))
            ranking = [
                {
                    "provider": row["provider"],
                    "median_ms": row["median_ms"],
                    "contracts_per_s": row["contracts_per_s"],
                }
                for row in sorted(case_rows, key=lambda row: float(row["median_ms"]))
            ]
            summary.append(
                {
                    "case": case.name,
                    "size": size,
                    "fastest_provider": fastest["provider"],
                    "fastest_median_ms": fastest["median_ms"],
                    "ranking": ranking,
                }
            )
    return summary


def _print_provider_inventory(providers: list[Provider]) -> None:
    print("Providers:")
    for provider in providers:
        supported = ", ".join(
            case.name for case in CASES if provider.supports(case.name)
        )
        cap = (
            f" (speed cap {provider.max_speed_size})" if provider.max_speed_size else ""
        )
        print(f" - {provider.name} [{provider.kind}] {provider.note}{cap}")
        print(f"   supported: {supported}")
    print()


def _print_accuracy_table(rows: list[dict[str, Any]], accuracy_size: int) -> None:
    print(f"Accuracy ({accuracy_size} contracts)")
    print(
        "IV accuracy is measured as price reconstruction error from the recovered IV."
    )
    header = f"{'Case':<22} {'Provider':<24} {'Max abs err':<14} {'RMSE':<14} {'Max rel err':<14}"
    print(header)
    print("-" * len(header))
    for row in rows:
        if "max_abs_error" not in row:
            continue
        print(
            f"{row['label']:<22} {row['provider']:<24} "
            f"{row['max_abs_error']:<14.6g} {row['rmse']:<14.6g} {row['max_rel_error']:<14.6g}"
        )
    print()


def _print_speed_table(rows: list[dict[str, Any]]) -> None:
    print(f"Speed (median of {N_RUNS} measured runs after {N_WARMUP} warmup)")
    header = f"{'Case':<22} {'Size':<8} {'Provider':<24} {'Median ms':<12} {'Contracts/s':<14} {'Peak alloc':<12}"
    print(header)
    print("-" * len(header))
    for row in rows:
        if "median_ms" not in row:
            continue
        peak = row.get("python_peak_allocation_bytes")
        peak_label = "n/a" if peak is None else str(peak)
        print(
            f"{row['label']:<22} {row['size']:<8} {row['provider']:<24} "
            f"{row['median_ms']:<12.4f} {row['contracts_per_s']:<14.2f} {peak_label:<12}"
        )
    print()


def run_benchmark(
    *,
    sizes: list[int],
    accuracy_size: int,
    json_path: str | None,
) -> dict[str, Any]:
    providers = available_providers()
    speed_chains = {
        size: _build_chain(size, seed=DEFAULT_SEED + size) for size in sizes
    }
    accuracy_chain = _build_chain(accuracy_size, seed=DEFAULT_SEED)

    accuracy_rows: list[dict[str, Any]] = []
    speed_rows: list[dict[str, Any]] = []

    _print_provider_inventory(providers)

    for case in CASES:
        for provider in providers:
            if not provider.supports(case.name):
                continue
            fn = provider.functions[case.name]
            actual = fn(accuracy_chain)
            if case.name == "bsm_call_iv":
                compared_actual = _reprice_bsm_call_from_iv(accuracy_chain, actual)
                expected = np.asarray(accuracy_chain["call_price"], dtype=np.float64)
            else:
                compared_actual = actual
                expected = np.asarray(
                    accuracy_chain[case.expected_key], dtype=np.float64
                )
            row = {
                "case": case.name,
                "label": case.label,
                "provider": provider.name,
                "provider_kind": provider.kind,
                "sample_size": accuracy_size,
                "accuracy_target": case.accuracy_target,
            }
            row.update(_accuracy_metrics(compared_actual, expected))
            if case.component_names is not None:
                row["component_names"] = list(case.component_names)
            accuracy_rows.append(row)

    _print_accuracy_table(accuracy_rows, accuracy_size)

    for case in CASES:
        for size in sizes:
            chain = speed_chains[size]
            for provider in providers:
                if not provider.supports(case.name):
                    continue
                if (
                    provider.max_speed_size is not None
                    and size > provider.max_speed_size
                ):
                    continue

                fn = provider.functions[case.name]
                samples_ms = _timed_runs_ms(fn, chain)
                stats = _summary_stats(samples_ms)
                median_ms = float(stats["median_ms"])
                contracts_per_s = _throughput_contracts_s(size, median_ms)
                peak_bytes = _python_peak_bytes(fn, chain)

                speed_rows.append(
                    {
                        "case": case.name,
                        "label": case.label,
                        "size": size,
                        "provider": provider.name,
                        "provider_kind": provider.kind,
                        "median_ms": round(median_ms, 4),
                        "contracts_per_s": round(contracts_per_s, 2),
                        "runs_ms": [round(sample, 4) for sample in samples_ms],
                        "stats": stats,
                        "python_peak_allocation_bytes": peak_bytes,
                        "input_layout": {
                            "dtype": "float64",
                            "contiguous": True,
                        },
                    }
                )

    _print_speed_table(speed_rows)

    package_names = ["numpy", "ferro-ta", "py_vollib"]
    metadata = benchmark_metadata(
        "benchmark_derivatives_compare",
        extra={
            "dataset": {
                "generator": "synthetic_option_chain",
                "speed_sizes": sizes,
                "accuracy_size": accuracy_size,
                "dtype": "float64",
                "array_layout": "C-contiguous",
                "seed": DEFAULT_SEED,
                "ranges": {
                    "spot": [80.0, 120.0],
                    "strike": [70.0, 130.0],
                    "rate": [0.0, 0.07],
                    "carry": [0.0, 0.03],
                    "time_to_expiry_years": [7.0 / 365.0, 2.0],
                    "volatility": [0.08, 0.65],
                },
            },
            "methodology": {
                "warmup_runs": N_WARMUP,
                "measured_runs": N_RUNS,
                "reported_metric": "median_ms",
                "speed_metric": "contracts_per_second",
                "accuracy_reference": (
                    "Scalar analytical Black-Scholes-Merton and Black-76 formulas "
                    "using math.erf; IV accuracy is measured as repriced error "
                    "from the recovered volatility because direct volatility "
                    "differences can be unstable on low-vega contracts."
                ),
                "input_layout_notes": (
                    "Benchmarks use contiguous float64 arrays. If your workload "
                    "passes non-contiguous arrays or mixed dtypes, benchmark that "
                    "path separately."
                ),
                "allocation_notes": (
                    "python_peak_allocation_bytes is a tracemalloc snapshot of "
                    "Python-tracked allocations only; it does not measure native RSS."
                ),
                "provider_notes": (
                    "reference_python_loop and py_vollib are scalar baselines and "
                    "are size-capped in the speed table to keep runtime reasonable."
                ),
            },
            "providers": [
                {
                    "name": provider.name,
                    "kind": provider.kind,
                    "note": provider.note,
                    "max_speed_size": provider.max_speed_size,
                    "supported_cases": [
                        case.name for case in CASES if provider.supports(case.name)
                    ],
                }
                for provider in providers
            ],
            "packages": package_versions(*package_names),
        },
    )

    result = {
        "schema_version": 1,
        "command": " ".join(["python", *sys.argv]),
        "n_warmup": N_WARMUP,
        "n_runs": N_RUNS,
        "accuracy_size": accuracy_size,
        "sizes": sizes,
        "metadata": metadata,
        "accuracy": {
            "summary": _accuracy_summary(accuracy_rows),
            "results": accuracy_rows,
        },
        "speed": {
            "summary": _speed_summary(speed_rows),
            "results": speed_rows,
        },
    }

    if json_path:
        output_path = Path(json_path)
        output_path.parent.mkdir(parents=True, exist_ok=True)
        output_path.write_text(json.dumps(result, indent=2), encoding="utf-8")
        print(f"Results written to {output_path}")

    return result


def main() -> int:
    parser = argparse.ArgumentParser(
        description="Compare ferro_ta derivatives analytics against reference implementations"
    )
    parser.add_argument(
        "--json",
        default=None,
        help="Write the benchmark artifact to JSON",
    )
    parser.add_argument(
        "--sizes",
        type=int,
        nargs="+",
        default=DEFAULT_SIZES,
        help="Contract counts to benchmark (default: 1000 10000)",
    )
    parser.add_argument(
        "--accuracy-size",
        type=int,
        default=DEFAULT_ACCURACY_SIZE,
        help="Contract count used for the accuracy pass (default: 512)",
    )
    args = parser.parse_args()
    run_benchmark(
        sizes=args.sizes,
        accuracy_size=args.accuracy_size,
        json_path=args.json,
    )
    return 0


if __name__ == "__main__":
    raise SystemExit(main())