#!/usr/bin/env python3
"""
Correlation metrics for human DA ↔ COMET-QE (and ΔQE ↔ ΔDA).

No third-party deps (stdlib only). Compatible with precision_sweep.py labels.

Examples:
  from correlation_metrics import pearson, spearman, kendall_tau, correlate_report

  report = correlate_report(da_scores, qe_scores)
  print(report)

  # CLI:
  python correlation_metrics.py --x da.txt --y qe.txt
  python correlation_metrics.py --jsonl labels.jsonl --x-key human_da --y-key comet_qe
"""
from __future__ import annotations

import argparse
import json
import math
from typing import Iterable, List, Optional, Sequence, Tuple


def _clean_pairs(
    x: Sequence[float], y: Sequence[float]
) -> Tuple[List[float], List[float]]:
    xs, ys = [], []
    for a, b in zip(x, y):
        if a is None or b is None:
            continue
        try:
            fa, fb = float(a), float(b)
        except (TypeError, ValueError):
            continue
        if math.isnan(fa) or math.isnan(fb) or math.isinf(fa) or math.isinf(fb):
            continue
        xs.append(fa)
        ys.append(fb)
    return xs, ys


def pearson(x: Sequence[float], y: Sequence[float]) -> float:
    """Pearson product-moment correlation coefficient r ∈ [-1, 1]."""
    xs, ys = _clean_pairs(x, y)
    n = len(xs)
    if n < 3:
        return float("nan")
    mx = sum(xs) / n
    my = sum(ys) / n
    num = sum((a - mx) * (b - my) for a, b in zip(xs, ys))
    denx = math.sqrt(sum((a - mx) ** 2 for a in xs))
    deny = math.sqrt(sum((b - my) ** 2 for b in ys))
    if denx == 0.0 or deny == 0.0:
        return float("nan")
    return num / (denx * deny)


def _rankdata(vals: Sequence[float]) -> List[float]:
    """Average ranks for ties; ranks start at 1."""
    ordered = sorted(enumerate(vals), key=lambda t: t[1])
    ranks = [0.0] * len(vals)
    i = 0
    while i < len(ordered):
        j = i
        while j + 1 < len(ordered) and ordered[j + 1][1] == ordered[i][1]:
            j += 1
        avg = (i + j) / 2.0 + 1.0
        for k in range(i, j + 1):
            ranks[ordered[k][0]] = avg
        i = j + 1
    return ranks


def spearman(x: Sequence[float], y: Sequence[float]) -> float:
    """Spearman rank correlation ρ (Pearson on ranks)."""
    xs, ys = _clean_pairs(x, y)
    if len(xs) < 3:
        return float("nan")
    return pearson(_rankdata(xs), _rankdata(ys))


def kendall_tau_detail(x: Sequence[float], y: Sequence[float]) -> dict:
    """
    Kendall's τ-b with full tie breakdown.

    Pair counts (i < j):
      concordant  — same order in x and y
      discordant  — opposite order
      tie_x       — tied in x only
      tie_y       — tied in y only
      tie_both    — tied in both (excluded from τ-b denominator)

    τ-b = (C - D) / sqrt((C + D + T_x) * (C + D + T_y))
    """
    xs, ys = _clean_pairs(x, y)
    n = len(xs)
    out = {
        "n": n,
        "concordant": 0,
        "discordant": 0,
        "tie_x": 0,
        "tie_y": 0,
        "tie_both": 0,
        "tau_b": float("nan"),
    }
    if n < 3:
        return out
    c = d = tx = ty = tb = 0
    for i in range(n):
        for j in range(i + 1, n):
            dx = xs[i] - xs[j]
            dy = ys[i] - ys[j]
            if dx == 0.0 and dy == 0.0:
                tb += 1
            elif dx == 0.0:
                tx += 1
            elif dy == 0.0:
                ty += 1
            elif (dx > 0 and dy > 0) or (dx < 0 and dy < 0):
                c += 1
            else:
                d += 1
    out.update(
        concordant=c, discordant=d, tie_x=tx, tie_y=ty, tie_both=tb
    )
    den = math.sqrt((c + d + tx) * (c + d + ty))
    if den == 0.0:
        return out
    out["tau_b"] = (c - d) / den
    return out


def kendall_tau(x: Sequence[float], y: Sequence[float]) -> float:
    """Kendall's τ-b (tie-adjusted). See kendall_tau_detail() for counts."""
    return kendall_tau_detail(x, y)["tau_b"]


def r_squared(x: Sequence[float], y: Sequence[float]) -> float:
    """Coefficient of determination from Pearson r."""
    r = pearson(x, y)
    if math.isnan(r):
        return float("nan")
    return r * r


def bootstrap_ci(
    x: Sequence[float],
    y: Sequence[float],
    stat_fn=spearman,
    n_boot: int = 1000,
    alpha: float = 0.05,
    seed: int = 42,
) -> Tuple[float, float, float]:
    """
    Bootstrap CI for a bivariate statistic.
    Returns (point_estimate, ci_low, ci_high).
    """
    xs, ys = _clean_pairs(x, y)
    n = len(xs)
    point = stat_fn(xs, ys)
    if n < 5 or n_boot < 1:
        return point, float("nan"), float("nan")

    # simple LCG for reproducibility without numpy
    state = seed & 0xFFFFFFFF

    def rnd() -> float:
        nonlocal state
        state = (1664525 * state + 1013904223) & 0xFFFFFFFF
        return state / 0x100000000

    samples = []
    for _ in range(n_boot):
        bx, by = [], []
        for _ in range(n):
            k = int(rnd() * n)
            if k >= n:
                k = n - 1
            bx.append(xs[k])
            by.append(ys[k])
        v = stat_fn(bx, by)
        if not math.isnan(v):
            samples.append(v)
    if not samples:
        return point, float("nan"), float("nan")
    samples.sort()
    lo_i = int(math.floor(alpha / 2 * len(samples)))
    hi_i = int(math.ceil((1 - alpha / 2) * len(samples))) - 1
    hi_i = max(0, min(hi_i, len(samples) - 1))
    return point, samples[lo_i], samples[hi_i]


def correlate_report(
    x: Sequence[float],
    y: Sequence[float],
    *,
    x_name: str = "x",
    y_name: str = "y",
    bootstrap: bool = True,
    n_boot: int = 1000,
) -> dict:
    """Full correlation report for DA ↔ QE style pairs."""
    xs, ys = _clean_pairs(x, y)
    n = len(xs)
    report = {
        "n": n,
        "x_name": x_name,
        "y_name": y_name,
        "pearson_r": pearson(xs, ys),
        "spearman_rho": spearman(xs, ys),
        "kendall_tau_b": kendall_tau(xs, ys),
        "r_squared": r_squared(xs, ys),
    }
    if bootstrap and n >= 5:
        for name, fn in (
            ("pearson_r", pearson),
            ("spearman_rho", spearman),
            ("kendall_tau_b", kendall_tau),
        ):
            pt, lo, hi = bootstrap_ci(xs, ys, stat_fn=fn, n_boot=n_boot)
            report[f"{name}_ci95"] = {"point": pt, "low": lo, "high": hi}
    return report


def delta_correlation(
    s_trunc: Sequence[float],
    s_complete: Sequence[float],
    da_trunc: Sequence[float],
    da_complete: Sequence[float],
) -> dict:
    """
    Correlation of truncation impact:
      ΔQE = s_complete - s_trunc
      ΔDA = da_complete - da_trunc
    """
    d_qe = [c - t for c, t in zip(s_complete, s_trunc)]
    d_da = [c - t for c, t in zip(da_complete, da_trunc)]
    rep = correlate_report(d_qe, d_da, x_name="delta_qe", y_name="delta_da")
    rep["mean_delta_qe"] = sum(d_qe) / len(d_qe) if d_qe else float("nan")
    rep["mean_delta_da"] = sum(d_da) / len(d_da) if d_da else float("nan")
    return rep


def tdi(
    s_trunc: Sequence[float],
    s_complete: Sequence[float],
    t_star: float,
) -> dict:
    """Truncation Decision Impact: P(s_trunc < t_star <= s_complete)."""
    st, sc = _clean_pairs(s_trunc, s_complete)
    if not st:
        return {"tdi": float("nan"), "n": 0, "flips": 0, "t_star": t_star}
    flips = sum(1 for a, b in zip(st, sc) if a < t_star <= b)
    return {
        "tdi": flips / len(st),
        "n": len(st),
        "flips": flips,
        "t_star": t_star,
        "definition": "P(s_trunc < t_star <= s_complete)",
    }



def goodman_kruskal_gamma(x: Sequence[float], y: Sequence[float]) -> float:
    """Goodman–Kruskal γ = (C - D) / (C + D); ties ignored in denominator."""
    d = kendall_tau_detail(x, y)
    c, disc = d["concordant"], d["discordant"]
    if c + disc == 0:
        return float("nan")
    return (c - disc) / (c + disc)


def somers_d(
    x: Sequence[float], y: Sequence[float], dependent: str = "y"
) -> float:
    """
    Somers' D.
    dependent='y' → D_yx = (C-D)/(C+D+T_y)  (y ordered by x)
    dependent='x' → D_xy = (C-D)/(C+D+T_x)
    """
    d = kendall_tau_detail(x, y)
    c, disc = d["concordant"], d["discordant"]
    if dependent == "y":
        den = c + disc + d["tie_y"]
    else:
        den = c + disc + d["tie_x"]
    if den == 0:
        return float("nan")
    return (c - disc) / den


def _inv_norm_approx(p: float) -> float:
    """Approximate Φ^{-1} (Acklam/rational approximation)."""
    if p <= 0.0:
        return -8.0
    if p >= 1.0:
        return 8.0
    # coefficients
    a = [
        -3.969683028665376e01,
        2.209460984245205e02,
        -2.759285104469687e02,
        1.383577518672690e02,
        -3.066479806614736e01,
        2.506628277459239e00,
    ]
    b = [
        -5.447609879822406e01,
        1.615858368580409e02,
        -1.556989798598866e02,
        6.680131188771972e01,
        -1.328068155288572e01,
    ]
    c = [
        -7.784894002430293e-03,
        -3.223964580411365e-01,
        -2.400758277161838e00,
        -2.549732539343734e00,
        4.374664141464968e00,
        2.938163982698783e00,
    ]
    d = [
        7.784695709041462e-03,
        3.224671290700398e-01,
        2.445134137142996e00,
        3.754408661907416e00,
    ]
    plow = 0.02425
    phigh = 1 - plow
    if p < plow:
        q = math.sqrt(-2 * math.log(p))
        return (
            (((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5])
            / ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1)
        )
    if p > phigh:
        q = math.sqrt(-2 * math.log(1 - p))
        return -(
            (((((c[0] * q + c[1]) * q + c[2]) * q + c[3]) * q + c[4]) * q + c[5])
            / ((((d[0] * q + d[1]) * q + d[2]) * q + d[3]) * q + 1)
        )
    q = p - 0.5
    r = q * q
    return (
        (((((a[0] * r + a[1]) * r + a[2]) * r + a[3]) * r + a[4]) * r + a[5]) * q
        / (((((b[0] * r + b[1]) * r + b[2]) * r + b[3]) * r + b[4]) * r + 1)
    )


def _norm_cdf_approx(z: float) -> float:
    """Approximate Φ(z) via erf-like rational."""
    # Abramowitz-Stegun style
    t = 1.0 / (1.0 + 0.2316419 * abs(z))
    d = 0.3989422804 * math.exp(-0.5 * z * z)
    p = d * t * (
        0.3193815
        + t * (-0.3565638 + t * (1.781478 + t * (-1.8212560 + t * 1.3302740)))
    )
    return 1.0 - p if z > 0 else p


_STAT_FN = {
    "pearson": pearson,
    "spearman": spearman,
    "kendall": kendall_tau,
    "gamma": goodman_kruskal_gamma,
}


def bootstrap_report(
    x: Sequence[float],
    y: Sequence[float],
    stat: str = "spearman",
    n_boot: int = 2000,
    alpha: float = 0.05,
    seed: int = 42,
) -> dict:
    """
    Bootstrap bias estimate + percentile / basic / BCa confidence intervals.

    Returns:
      point, mean_boot, bias_boot, point_bc,
      ci_percentile, ci_basic, ci_bca, n, n_boot, stat
    """
    xs, ys = _clean_pairs(x, y)
    n = len(xs)
    fn = _STAT_FN.get(stat, spearman)
    point = fn(xs, ys)
    empty = {
        "stat": stat,
        "n": n,
        "n_boot": n_boot,
        "point": point,
        "mean_boot": float("nan"),
        "bias_boot": float("nan"),
        "point_bc": float("nan"),
        "ci_percentile": {"low": float("nan"), "high": float("nan")},
        "ci_basic": {"low": float("nan"), "high": float("nan")},
        "ci_bca": {"low": float("nan"), "high": float("nan")},
        "z0": float("nan"),
        "acceleration": float("nan"),
    }
    if n < 5 or math.isnan(point):
        return empty

    state = seed & 0xFFFFFFFF

    def rnd() -> float:
        nonlocal state
        state = (1664525 * state + 1013904223) & 0xFFFFFFFF
        return state / 0x100000000

    boots = []
    for _ in range(n_boot):
        bx, by = [], []
        for _ in range(n):
            k = int(rnd() * n)
            if k >= n:
                k = n - 1
            bx.append(xs[k])
            by.append(ys[k])
        v = fn(bx, by)
        if not math.isnan(v):
            boots.append(v)
    if len(boots) < 50:
        return empty

    boots.sort()
    mean_boot = sum(boots) / len(boots)
    bias = mean_boot - point
    point_bc = point - bias  # = 2*point - mean_boot

    def percentile(p: float) -> float:
        if not boots:
            return float("nan")
        i = p * (len(boots) - 1)
        lo = int(math.floor(i))
        hi = int(math.ceil(i))
        if lo == hi:
            return boots[lo]
        w = i - lo
        return boots[lo] * (1 - w) + boots[hi] * w

    lo_p, hi_p = percentile(alpha / 2), percentile(1 - alpha / 2)
    # basic (reverse percentile)
    lo_b, hi_b = 2 * point - hi_p, 2 * point - lo_p

    # z0 bias correction for BCa
    prop = sum(1 for v in boots if v < point) / len(boots)
    # avoid 0/1
    prop = min(max(prop, 1.0 / (len(boots) + 1)), 1.0 - 1.0 / (len(boots) + 1))
    z0 = _inv_norm_approx(prop)

    # jackknife acceleration
    jack = []
    for i in range(n):
        xj = xs[:i] + xs[i + 1 :]
        yj = ys[:i] + ys[i + 1 :]
        jv = fn(xj, yj)
        if not math.isnan(jv):
            jack.append(jv)
    a = 0.0
    if len(jack) >= 3:
        jmean = sum(jack) / len(jack)
        diffs = [jmean - j for j in jack]
        s2 = sum(d * d for d in diffs)
        s3 = sum(d * d * d for d in diffs)
        if s2 > 0:
            a = s3 / (6.0 * (s2 ** 1.5))

    z_alpha = _inv_norm_approx(alpha / 2)
    z_one = _inv_norm_approx(1 - alpha / 2)

    def bca_alpha(z_a: float) -> float:
        num = z0 + z_a
        den = 1.0 - a * num
        if den == 0:
            return 0.5
        return _norm_cdf_approx(z0 + num / den)

    a1, a2 = bca_alpha(z_alpha), bca_alpha(z_one)
    a1 = min(max(a1, 0.0), 1.0)
    a2 = min(max(a2, 0.0), 1.0)
    if a1 > a2:
        a1, a2 = a2, a1
    lo_bca, hi_bca = percentile(a1), percentile(a2)

    return {
        "stat": stat,
        "n": n,
        "n_boot": len(boots),
        "point": point,
        "mean_boot": mean_boot,
        "bias_boot": bias,
        "point_bc": point_bc,
        "ci_percentile": {"low": lo_p, "high": hi_p},
        "ci_basic": {"low": lo_b, "high": hi_b},
        "ci_bca": {"low": lo_bca, "high": hi_bca},
        "z0": z0,
        "acceleration": a,
    }




def jackknife_acceleration(
    x: Sequence[float],
    y: Sequence[float],
    stat_fn=spearman,
) -> dict:
    """
    Jackknife estimate of acceleration `a` for BCa intervals.

    a = sum(d_i^3) / (6 * (sum(d_i^2))^(3/2))
    where d_i = mean(theta_(·)) - theta_(i), theta_(i) = stat without observation i.

    Also returns leave-one-out values and jackknife bias estimate.
    """
    xs, ys = _clean_pairs(x, y)
    n = len(xs)
    out = {
        "n": n,
        "acceleration": float("nan"),
        "jackknife_values": [],
        "jackknife_mean": float("nan"),
        "jackknife_bias": float("nan"),
        "point": float("nan"),
    }
    if n < 4:
        return out
    point = stat_fn(xs, ys)
    out["point"] = point
    jack = []
    for i in range(n):
        xj = xs[:i] + xs[i + 1 :]
        yj = ys[:i] + ys[i + 1 :]
        jv = stat_fn(xj, yj)
        if not math.isnan(jv):
            jack.append(jv)
    if len(jack) < 3:
        return out
    jmean = sum(jack) / len(jack)
    diffs = [jmean - j for j in jack]
    s2 = sum(d * d for d in diffs)
    s3 = sum(d * d * d for d in diffs)
    a = 0.0
    if s2 > 0:
        a = s3 / (6.0 * (s2 ** 1.5))
    # jackknife bias ≈ (n-1) * (jmean - point)
    jbias = (n - 1) * (jmean - point) if not math.isnan(point) else float("nan")
    out.update(
        acceleration=a,
        jackknife_values=jack,
        jackknife_mean=jmean,
        jackknife_bias=jbias,
    )
    return out


def permutation_test(
    x: Sequence[float],
    y: Sequence[float],
    stat: str = "spearman",
    n_perm: int = 5000,
    alternative: str = "two-sided",
    seed: int = 42,
) -> dict:
    """
    Permutation test of association: shuffle y, recompute statistic.

    Under H0: x and y independent (no ordinal/linear association).
    p-value = fraction of |perm| as extreme as |observed| (two-sided),
    or one-sided per `alternative` in {'two-sided','greater','less'}.

    Returns observed, p_value, n_perm, null_mean, null_std, etc.
    """
    xs, ys = _clean_pairs(x, y)
    n = len(xs)
    fn = _STAT_FN.get(stat, spearman)
    observed = fn(xs, ys)
    out = {
        "stat": stat,
        "n": n,
        "n_perm": n_perm,
        "alternative": alternative,
        "observed": observed,
        "p_value": float("nan"),
        "null_mean": float("nan"),
        "null_std": float("nan"),
        "null_extreme_count": 0,
    }
    if n < 4 or math.isnan(observed):
        return out

    state = seed & 0xFFFFFFFF

    def rnd() -> float:
        nonlocal state
        state = (1664525 * state + 1013904223) & 0xFFFFFFFF
        return state / 0x100000000

    def shuffle(arr: List[float]) -> List[float]:
        a = list(arr)
        # Fisher–Yates
        for i in range(len(a) - 1, 0, -1):
            j = int(rnd() * (i + 1))
            if j > i:
                j = i
            a[i], a[j] = a[j], a[i]
        return a

    nulls = []
    extreme = 0
    for _ in range(n_perm):
        y_p = shuffle(ys)
        v = fn(xs, y_p)
        if math.isnan(v):
            continue
        nulls.append(v)
        if alternative == "greater":
            if v >= observed:
                extreme += 1
        elif alternative == "less":
            if v <= observed:
                extreme += 1
        else:  # two-sided
            if abs(v) >= abs(observed):
                extreme += 1

    if not nulls:
        return out
    # add-one smoothing for p-value
    p = (extreme + 1) / (len(nulls) + 1)
    mu = sum(nulls) / len(nulls)
    var = sum((v - mu) ** 2 for v in nulls) / max(len(nulls) - 1, 1)
    out.update(
        p_value=p,
        null_mean=mu,
        null_std=math.sqrt(var),
        null_extreme_count=extreme,
        n_perm_effective=len(nulls),
    )
    return out


def permutation_ci(
    x: Sequence[float],
    y: Sequence[float],
    stat: str = "spearman",
    n_perm: int = 5000,
    alpha: float = 0.05,
    seed: int = 42,
) -> dict:
    """
    Permutation-based confidence interval via studentized-ish percentile of
    the null-centered distribution mapped back to the observed scale.

    Practical approach used here (correlation-friendly):
      1. Build null distribution of stat under shuffled y (H0).
      2. Also bootstrap (paired resample) the sampling distribution of the stat.
      3. Report:
         - permutation p-value (from permutation_test)
         - percentile CI from paired bootstrap (same as bootstrap_report percentile)
         - 'permutation_test_ci' as the range of null stats at alpha/2 (diagnostic
           under H0 — not a CI for the parameter, but a null reference band)

    For a true parameter CI, prefer BCa from bootstrap_report.
    This function emphasizes *significance* (permutation p) + null band clarity.
    """
    pt = permutation_test(x, y, stat=stat, n_perm=n_perm, seed=seed)
    br = bootstrap_report(x, y, stat=stat, n_boot=max(1000, n_perm // 2), alpha=alpha, seed=seed)

    # Null reference band from permutation distribution
    xs, ys = _clean_pairs(x, y)
    fn = _STAT_FN.get(stat, spearman)
    state = (seed ^ 0xA5A5A5A5) & 0xFFFFFFFF

    def rnd() -> float:
        nonlocal state
        state = (1664525 * state + 1013904223) & 0xFFFFFFFF
        return state / 0x100000000

    nulls = []
    ys_list = list(ys)
    n = len(xs)
    for _ in range(n_perm):
        # shuffle
        a = list(ys_list)
        for i in range(n - 1, 0, -1):
            j = int(rnd() * (i + 1))
            if j > i:
                j = i
            a[i], a[j] = a[j], a[i]
        v = fn(xs, a)
        if not math.isnan(v):
            nulls.append(v)
    nulls.sort()

    def perc(sorted_vals, p):
        if not sorted_vals:
            return float("nan")
        i = p * (len(sorted_vals) - 1)
        lo = int(math.floor(i))
        hi = int(math.ceil(i))
        if lo == hi:
            return sorted_vals[lo]
        w = i - lo
        return sorted_vals[lo] * (1 - w) + sorted_vals[hi] * w

    null_band = {
        "low": perc(nulls, alpha / 2),
        "high": perc(nulls, 1 - alpha / 2),
        "note": "H0 reference band for the statistic under shuffled y — not a CI for rho",
    }

    return {
        "stat": stat,
        "n": pt["n"],
        "observed": pt["observed"],
        "permutation_p_value": pt["p_value"],
        "permutation_null_mean": pt["null_mean"],
        "permutation_null_std": pt["null_std"],
        "null_reference_band": null_band,
        "bootstrap_point_bc": br.get("point_bc"),
        "bootstrap_ci_bca": br.get("ci_bca"),
        "bootstrap_ci_percentile": br.get("ci_percentile"),
        "guidance": (
            "Use permutation_p_value for significance under independence. "
            "Use bootstrap_ci_bca for confidence interval of the correlation parameter. "
            "null_reference_band is the H0 sampling band, not a parameter CI."
        ),
    }



# ── CLI ───────────────────────────────────────────────────────────────

def _read_column(path: str) -> List[float]:
    vals = []
    with open(path, encoding="utf-8") as f:
        for line in f:
            line = line.strip()
            if not line:
                continue
            vals.append(float(line.split()[0]))
    return vals


def _read_jsonl(path: str, x_key: str, y_key: str):
    xs, ys = [], []
    with open(path, encoding="utf-8") as f:
        for line in f:
            line = line.strip()
            if not line:
                continue
            row = json.loads(line)
            xs.append(row[x_key])
            ys.append(row[y_key])
    return xs, ys


def main():
    ap = argparse.ArgumentParser(description="Correlation metrics (Pearson / Spearman / Kendall)")
    ap.add_argument("--x", type=str, help="Text file: one number per line")
    ap.add_argument("--y", type=str, help="Text file: one number per line")
    ap.add_argument("--jsonl", type=str, help="JSONL with numeric fields")
    ap.add_argument("--x-key", default="human_da")
    ap.add_argument("--y-key", default="comet_qe")
    ap.add_argument("--no-bootstrap", action="store_true")
    ap.add_argument("--n-boot", type=int, default=1000)
    args = ap.parse_args()

    if args.jsonl:
        xs, ys = _read_jsonl(args.jsonl, args.x_key, args.y_key)
        x_name, y_name = args.x_key, args.y_key
    elif args.x and args.y:
        xs, ys = _read_column(args.x), _read_column(args.y)
        x_name, y_name = "x", "y"
    else:
        # self-demo
        xs = [10, 20, 30, 40, 50, 60, 70, 80, 90]
        ys = [0.2, 0.35, 0.4, 0.55, 0.5, 0.7, 0.75, 0.85, 0.9]
        x_name, y_name = "demo_da", "demo_qe"
        print("No input — running demo series\n")

    report = correlate_report(
        xs,
        ys,
        x_name=x_name,
        y_name=y_name,
        bootstrap=not args.no_bootstrap,
        n_boot=args.n_boot,
    )
    print(json.dumps(report, indent=2))


if __name__ == "__main__":
    main()