#!/usr/bin/env python3
"""Independent stdlib logistic MLE. Never emits or persists participant-level data."""
import argparse, csv, hashlib, json, math, pathlib, platform

NAMES = ['Intercept', 'Knowledge_Score_Total', 'Q43', 'Q44', 'Q15']

def pdf_score(row):
    items = [row['Q'+str(i)].strip() for i in range(23, 42)]
    if any(v and float(v) not in (0, 1) for v in items):
        raise ValueError('Expected binary correctness items or blanks')
    return '' if any(not v for v in items) else str(sum(
        1-int(float(v)) if i in (31, 37) else int(float(v)) for i, v in enumerate(items, 23)))

def solve(a, b):
    n = len(b)
    m = [list(row) + [v] for row, v in zip(a, b)]
    for k in range(n):
        p = max(range(k, n), key=lambda i: abs(m[i][k]))
        if abs(m[p][k]) < 1e-14:
            raise ValueError('Singular information matrix')
        m[p], m[k] = m[k], m[p]
        d = m[k][k]
        m[k] = [v/d for v in m[k]]
        for i in range(n):
            if i != k:
                f = m[i][k]
                m[i] = [v-f*w for v, w in zip(m[i], m[k])]
    return [r[-1] for r in m]

def sigmoid(z):
    return 1/(1+math.exp(-z)) if z >= 0 else math.exp(z)/(1+math.exp(z))

def evaluate(x, y, beta):
    n = len(beta)
    gradient = [0.0]*n
    info = [[0.0]*n for _ in range(n)]
    terms = []
    for row, outcome in zip(x, y):
        z = math.fsum(a*b for a, b in zip(row, beta))
        p = sigmoid(z)
        terms.append(outcome*z-max(z, 0)-math.log1p(math.exp(-abs(z))))
        for j in range(n):
            gradient[j] += row[j]*(outcome-p)
            for k in range(n):
                info[j][k] += row[j]*row[k]*p*(1-p)
    return math.fsum(terms), gradient, info

def fit(x, y, initial=None):
    beta = list(initial or [0.0]*len(x[0]))
    history = []
    for iteration in range(1, 101):
        ll, g, information = evaluate(x, y, beta)
        step = solve(information, g)
        scale = 1.0
        while True:
            candidate = [b+scale*d for b, d in zip(beta, step)]
            new_ll, _, _ = evaluate(x, y, candidate)
            if new_ll >= ll-1e-12:
                break
            scale /= 2
            if scale < 2**-30:
                raise ValueError('Line search failed')
        beta = candidate
        history.append({'iteration': iteration, 'log_likelihood': new_ll,
                        'step_max_abs': max(abs(scale*d) for d in step), 'step_scale': scale})
        if max(abs(scale*d) for d in step) < 1e-10:
            break
    else:
        raise ValueError('No convergence within 100 iterations')
    ll, g, info = evaluate(x, y, beta)
    inverse_columns = [solve(info, [float(j == k) for j in range(len(beta))]) for k in range(len(beta))]
    covariance = [list(row) for row in zip(*inverse_columns)]
    return beta, covariance, {'converged': True, 'iterations': iteration, 'log_likelihood': ll,
                              'gradient_max_abs': max(map(abs, g)), 'history': history}

def main():
    p = argparse.ArgumentParser()
    p.add_argument('--csv')
    p.add_argument('--output')
    p.add_argument('--pdf-key', action='store_true', help='Correct historical Q31/Q37 correctness bits and sum all 19 items before fitting')
    p.add_argument('--self-test', action='store_true', help='Run a synthetic closed-form check without study data')
    args = p.parse_args()
    if args.self_test:
        x = [[1., 0.]]*100 + [[1., 1.]]*100
        y = [1.]*30 + [0.]*70 + [1.]*60 + [0.]*40
        beta, cov, _ = fit(x, y)
        expected = [math.log(30/70), math.log(60/40)-math.log(30/70)]
        assert max(abs(a-b) for a,b in zip(beta, expected)) < 1e-10
        assert abs(cov[0][0]-(1/30+1/70)) < 1e-10
        assert abs(cov[1][1]-(1/30+1/70+1/60+1/40)) < 1e-10
        assert abs(cov[0][1]+(1/30+1/70)) < 1e-10
        for bit, expected_score in [('0', '2'), ('1', '17')]:
            row = {'Q'+str(i): bit for i in range(23, 42)}
            assert pdf_score(row) == expected_score
            for i in range(23, 42):
                assert pdf_score({**row, 'Q'+str(i): ''}) == ''
        row = {'Q'+str(i): '0' for i in range(23, 42)}
        assert pdf_score({**row, 'Q31': '1', 'Q37': '1'}) == '0'
        try:
            pdf_score({**row, 'Q31': '999'})
        except ValueError:
            pass
        else:
            raise AssertionError('Invalid score accepted')
        print('PASS: synthetic logistic coefficients/covariance and PDF-key score correction, missingness and invalid-item checks.')
        return
    if not args.csv or not args.output:
        p.error('--csv and --output are required unless --self-test is used')
    source = pathlib.Path(args.csv)
    if args.pdf_key and hashlib.sha256(source.read_bytes()).hexdigest() != '501042c4be0c3c7f40eab1cc1907d0e23607b609212af8fbcd846209c88c7e37':
        raise ValueError('PDF-key correction requires the identified historical input, not an already corrected derivative')
    x, y = [], []
    total = excluded = invalid = missing = 0
    with source.open(newline='', encoding='utf-8-sig') as f:
        reader = csv.DictReader(f)
        required = ['Q21'] + NAMES[1:]
        assert all(k in reader.fieldnames for k in required), 'Missing required column'
        for row in reader:
            total += 1
            if args.pdf_key:
                row['Knowledge_Score_Total'] = pdf_score(row)
            if any(row[k] is None or not row[k].strip() for k in required):
                excluded += 1
                missing += 1
                continue
            try:
                vals = [float(row[k]) for k in required]
            except ValueError:
                excluded += 1
                invalid += 1
                continue
            if not all(math.isfinite(v) for v in vals) or vals[0] not in (0, 1):
                excluded += 1
                invalid += 1
                continue
            y.append(vals[0])
            x.append([1.0]+vals[1:])
    beta, cov, diagnostics = fit(x, y)
    # A different initial point checks optimizer stability without using any Sutrix code.
    beta2, _, _ = fit(x, y, [0.25, -0.1, 0.1, -0.2, 0.05])
    diagnostics['alternate_start_max_abs_coefficient_difference'] = max(abs(a-b) for a,b in zip(beta,beta2))
    _, g, info = evaluate(x, y, beta)
    eps = 1e-5
    gradient_error = []
    hessian_error = []
    for j in range(len(beta)):
        plus, minus = list(beta), list(beta)
        plus[j] += eps
        minus[j] -= eps
        lp, gp, _ = evaluate(x, y, plus)
        lm, gm, _ = evaluate(x, y, minus)
        gradient_error.append(abs((lp-lm)/(2*eps)-g[j]))
        hessian_error.extend(abs((gp[k]-gm[k])/(2*eps)+info[k][j]) for k in range(len(beta)))
    diagnostics['finite_difference_gradient_max_abs_error'] = max(gradient_error)
    diagnostics['finite_difference_information_max_abs_error'] = max(hessian_error)
    rows = []
    for j, name in enumerate(NAMES):
        se = math.sqrt(cov[j][j])
        z = beta[j]/se
        lo, hi = beta[j]-1.959963984540054*se, beta[j]+1.959963984540054*se
        rows.append({'term': name, 'coefficient': beta[j], 'standard_error': se, 'odds_ratio': math.exp(beta[j]),
                     'coefficient_wald95ci': [lo, hi], 'odds_ratio_wald95ci': [math.exp(lo), math.exp(hi)],
                     'wald_z': z, 'p_two_sided': math.erfc(abs(z)/math.sqrt(2))})
    report = {'scope': 'Same-cleaned-input implementation cross-check; not independent sample replication, clinical validation, or raw-cleaning proof.',
              'scoring_version': 'pdf-key-v2' if args.pdf_key else 'historical-v1',
              'input_role': 'Historical cleaned CSV; Q31/Q37 correction and 19-item sum are applied in memory before fitting. No derivative CSV is emitted.' if args.pdf_key else 'Historical cleaned CSV, using its saved total without correction.',
              'method': 'Unweighted, unpenalized Bernoulli logistic maximum likelihood; intercept; predictors treated as supplied numeric values; Newton steps with backtracking; model-based inverse-information standard errors; two-sided normal Wald tests and 95% confidence intervals.',
              'outcome': 'Q21=1 versus Q21=0', 'predictors': NAMES[1:],
              'csv_sha256': hashlib.sha256(source.read_bytes()).hexdigest(), 'python': platform.python_version(),
              'counts': {'total': total, 'complete_cases': len(y), 'excluded': excluded, 'missing_required_field': missing, 'invalid_numeric_or_outcome': invalid},
              'coefficients': rows, 'diagnostics': diagnostics,
              'limitations': ['Uses the existing cleaned input; does not verify raw recoding, consent, sampling, construct validity, or cleaning provenance.',
                              'No survey weighting, clustering, robust variance, interactions, nonlinearity assessment, out-of-sample evaluation, or causal interpretation.',
                              'Agreement would establish this model calculation only; it would not validate the full engine or every reported analysis.']}
    pathlib.Path(args.output).write_text(json.dumps(report, indent=2)+'\n')
    print(json.dumps(report, indent=2))

if __name__ == '__main__':
    main()