# (c) 2026 William J. Dwyer. Free and open-source software under the MIT License (see LICENSE).
# If you use or build upon this software, please retain this notice and cite the accompanying
# methodology paper together with the archived software (Zenodo DOI: 10.5281/zenodo.20709963).
"""
cvscreen_core — pure statistical functions for the CV-ratio screen.
No GUI dependencies; importable from the Streamlit app or any script.

Companion to: Dwyer WJ, "The Coefficient of Variation as a Cutoff for
Detecting Non-Normal Distributions from Summary Statistics."
"""
import numpy as np
from scipy import stats


def lognormal_params(mu_n, cv):
    """Natural-domain mean and CV -> log-domain (mu_e, sigma_e). Paper eqs 10, 13."""
    s2 = np.log(cv ** 2 + 1.0)
    return np.log(mu_n) - s2 / 2.0, np.sqrt(s2)


def draw_lognormal(mu_n, cv, n, reps, rng, logshift=0.0):
    """Draw a reps x n array of log-normal samples with given natural mean & CV."""
    mu_e, s_e = lognormal_params(mu_n, cv)
    return rng.lognormal(mu_e + logshift, s_e, size=(reps, n))


def permutation_t_p(a, b, B=2000, use_log=True, rng=None):
    """Two-sided permutation test on the difference in (log) means. Size is exact by construction."""
    rng = rng or np.random.default_rng()
    a = np.asarray(a, float); b = np.asarray(b, float)
    if use_log:
        a, b = np.log(a), np.log(b)
    obs = abs(a.mean() - b.mean())
    pool = np.concatenate([a, b]); n1 = len(a); count = 0
    for _ in range(B):
        p = rng.permutation(pool)
        if abs(p[:n1].mean() - p[n1:].mean()) >= obs - 1e-12:
            count += 1
    return (count + 1) / (B + 1)


def cv_ratio_screen(mu1, sd1, mu2, sd2):
    """Apply the eq-17 / CV-ratio screen to two groups' summary statistics."""
    cv1, cv2 = sd1 / mu1, sd2 / mu2
    cv_hi, cv_lo = max(cv1, cv2), min(cv1, cv2)
    ratio = cv_hi / cv_lo if cv_lo > 0 else float("inf")
    larger_mean_has_larger_var = (sd1 >= sd2) == (mu1 >= mu2)
    # CV ~ 0.5 is where the log transform begins to matter (paper Result 4);
    # a between-group CV ratio > 1.5 with the larger mean also flags (eq 17).
    suspect = (cv_hi >= 0.5) or (ratio > 1.5 and larger_mean_has_larger_var)
    if cv_hi < 0.3:
        verdict = "Effectively normal - raw-scale t-test is appropriate."
    elif cv_hi < 0.5:
        verdict = "Borderline - log transform begins to matter; check the raw data if possible."
    elif cv_hi < 1.0:
        verdict = "Suspect non-normality - prefer a log-scale analysis."
    else:
        verdict = "Strongly suspect log-normality - a raw-scale t-test is likely deficient; use the log scale."
    return dict(cv1=cv1, cv2=cv2, cv_ratio=ratio, suspect=bool(suspect), verdict=verdict)


def analyze_two_groups(x, y, n_perm=2000, seed=0):
    """Full analysis of two raw samples: screen, scale choice, and competing tests."""
    rng = np.random.default_rng(seed)
    x = np.asarray(x, float); y = np.asarray(y, float)
    R = {}
    R["n1"], R["n2"] = len(x), len(y)
    R["mean1"], R["sd1"] = x.mean(), x.std(ddof=1)
    R["mean2"], R["sd2"] = y.mean(), y.std(ddof=1)
    positive = bool((x > 0).all() and (y > 0).all())
    R["positive"] = positive
    R["cv1"] = R["sd1"] / R["mean1"] if R["mean1"] else np.nan
    R["cv2"] = R["sd2"] / R["mean2"] if R["mean2"] else np.nan
    pooled = np.concatenate([x, y])
    R["pooled_cv"] = pooled.std(ddof=1) / pooled.mean() if pooled.mean() else np.nan
    R["sw_p_raw_1"] = stats.shapiro(x).pvalue if len(x) >= 3 else np.nan
    R["sw_p_raw_2"] = stats.shapiro(y).pvalue if len(y) >= 3 else np.nan
    R["p_raw_t"] = stats.ttest_ind(x, y, equal_var=False).pvalue
    R["p_mw"] = stats.mannwhitneyu(x, y, alternative="two-sided").pvalue
    R["p_log_t"] = R["p_perm_log"] = R["sw_p_log_1"] = R["sw_p_log_2"] = None
    if positive:
        R["sw_p_log_1"] = stats.shapiro(np.log(x)).pvalue if len(x) >= 3 else np.nan
        R["sw_p_log_2"] = stats.shapiro(np.log(y)).pvalue if len(y) >= 3 else np.nan
        R["p_log_t"] = stats.ttest_ind(np.log(x), np.log(y), equal_var=False).pvalue
        R["p_perm_log"] = permutation_t_p(x, y, B=n_perm, use_log=True, rng=rng)
    scr = cv_ratio_screen(R["mean1"], R["sd1"], R["mean2"], R["sd2"])
    R.update({"screen_" + k: v for k, v in scr.items()})
    if not positive:
        R["recommended"] = ("Data contain non-positive values: log transform not applicable. "
                            "Use the raw-scale t-test if approximately normal, else Mann-Whitney or a permutation test.")
    elif scr["suspect"]:
        R["recommended"] = ("Log-scale t-test (most powerful under log-normality); report the permutation "
                            "t-test on the log scale as an assumption-light confirmation.")
    else:
        R["recommended"] = "Raw-scale t-test is appropriate."
    return R


if __name__ == "__main__":
    # self-test against the paper's worked example
    print("lognormal_params(0.5, 1.0)  ->", tuple(round(v, 3) for v in lognormal_params(0.5, 1.0)), "expect (-1.040, 0.833)")
    print("lognormal_params(1.5, 3.4/1.5) ->", tuple(round(v, 3) for v in lognormal_params(1.5, 3.4 / 1.5)), "expect (-0.502, 1.347)")
    s = cv_ratio_screen(0.5, 0.5, 1.5, 3.4)
    print("screen(0.5±0.5 vs 1.5±3.4) -> CVs", round(s["cv1"], 3), round(s["cv2"], 3), "ratio", round(s["cv_ratio"], 3), "suspect", s["suspect"])
    rng = np.random.default_rng(0)
    x = rng.lognormal(*lognormal_params(1.0, 0.8), size=20)
    y = rng.lognormal(*lognormal_params(1.6, 0.8), size=20)
    R = analyze_two_groups(x, y, n_perm=1000)
    print("analyze demo -> raw_t={:.4f} log_t={:.4f} perm={:.4f} MW={:.4f}".format(R["p_raw_t"], R["p_log_t"], R["p_perm_log"], R["p_mw"]))
    print("recommended:", R["recommended"])
