"""Independent numerical answer key using Fraction arithmetic and source XML.

No imports from ligninqc, no reading the generated CSVs, and no fitted package
outputs. Display decimals from the source are treated as exact input values;
this is computational verification, not an uncertainty model for experiments.
"""
import argparse
from decimal import Decimal, localcontext
from fractions import Fraction as F
from pathlib import Path
import hashlib
import json
import re
import xml.etree.ElementTree as ET

PARSER = argparse.ArgumentParser(description="Recompute the real-case answer key independently with exact rational arithmetic. Reads only source XML; does not import ligninqc.")
PARSER.add_argument("--source-dir", type=Path, default=Path(__file__).resolve().parents[1] / "data/sources", help="Directory containing PMC6539611.xml and PMC9331396.xml")
PARSER.add_argument("--out", type=Path, default=Path("independent-case-analysis.json"), help="Output JSON path")
ARGS = PARSER.parse_args()
EVIDENCE = ARGS.source_dir

def text(element):
    return "".join(element.itertext()).strip()

def rows(tree, table):
    return [[text(c) for c in r] for r in tree.findall(f'.//table-wrap[@id="{table}"]//tbody/tr')]

def sqrt(value):
    with localcontext() as ctx:
        ctx.prec = 50
        return float((Decimal(value.numerator) / Decimal(value.denominator)).sqrt())

def avg(xs):
    return sum(xs, F(0)) / len(xs)

def ranks(xs):
    return [F(sum(y < x for y in xs)) + F(sum(y == x for y in xs) + 1, 2) for x in xs]

def fit(xs, ys):
    mx, my = avg(xs), avg(ys)
    xx = sum((x - mx) ** 2 for x in xs)
    yy = sum((y - my) ** 2 for y in ys)
    xy = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
    return {"slope": xy / xx, "intercept": my - xy / xx * mx,
            "r": (1 if xy >= 0 else -1) * sqrt(xy * xy / (xx * yy)),
            "r2": xy * xy / (xx * yy)}

def analyze(points):
    ids, xs, ys = map(list, zip(*points))
    model = fit(xs, ys)
    rankfit = fit(ranks(xs), ranks(ys))
    slopes = []
    folds = []
    model_errors, mean_errors = [], []
    for i in range(len(xs)):
        trainx, trainy = xs[:i] + xs[i+1:], ys[:i] + ys[i+1:]
        training = fit(trainx, trainy)
        prediction = training["intercept"] + training["slope"] * xs[i]
        baseline = avg(trainy)
        model_errors.append(prediction - ys[i])
        mean_errors.append(baseline - ys[i])
        slopes.append({"excluded_compound_id": ids[i], "remaining_r2": float(training["r2"])})
        folds.append({"held_out_compound_id": ids[i], "observed_rdi": float(ys[i]),
                      "linear_prediction": float(prediction), "training_mean_prediction": float(baseline)})
    positive, negative, tie_x, tie_y, tie_both = 0, 0, 0, 0, 0
    reversal_pairs = []
    for i in range(len(xs)):
        for j in range(i + 1, len(xs)):
            dx, dy = xs[j] - xs[i], ys[j] - ys[i]
            if dx == 0 and dy == 0:
                tie_both += 1
            elif dx == 0:
                tie_x += 1
            elif dy == 0:
                tie_y += 1
            elif dx * dy > 0:
                positive += 1
                reversal_pairs.append([ids[i], ids[j]])
            else:
                negative += 1
    model_mse, mean_mse = avg([e*e for e in model_errors]), avg([e*e for e in mean_errors])
    return {
        "n": len(xs), "compound_ids": ids,
        "pearson_r": model["r"], "r2": float(model["r2"]),
        "slope": float(model["slope"]), "intercept": float(model["intercept"]),
        "spearman_rho": rankfit["r"],
        "rank_pairs": {"total": len(xs) * (len(xs)-1)//2,
            "consistent_with_lower_bde_higher_rdi": negative,
            "discordant_with_lower_bde_higher_rdi": positive,
            "tie_bde_only": tie_x, "tie_rdi_only": tie_y, "tie_both": tie_both,
            "untied_pairs": negative + positive,
            "discordance_fraction_among_untied": positive / (negative + positive),
            "discordant_compound_pairs": reversal_pairs},
        "loo": {"linear_mae": float(avg([abs(e) for e in model_errors])), "linear_rmse": sqrt(model_mse),
            "training_mean_mae": float(avg([abs(e) for e in mean_errors])), "training_mean_rmse": sqrt(mean_mse),
            "mse_skill_vs_training_mean": float(1 - model_mse / mean_mse), "folds": folds},
        "leave_one_compound_out_association_sensitivity": {
            "min_r2": min(x["remaining_r2"] for x in slopes),
            "max_r2": max(x["remaining_r2"] for x in slopes), "each_exclusion": slopes},
    }

source = ET.parse(EVIDENCE / "PMC6539611.xml")
records = []
for r in rows(source, "molecules-24-01794-t001"):
    def outcome(cell):
        return F(re.match(r"[0-9.]+", cell).group()) if re.match(r"[0-9.]+", cell) else None
    records.append({"id": int(r[0]), "neutral": F(r[2].split('/')[0]),
                    "carboxylate": F(r[2].split('/')[-1]), "dpph": outcome(r[3]), "abts": outcome(r[5]),
                    "phenolic_oh": 2 if int(r[0]) in (16, 17) else 1})

definitions = [
    ("dpph_neutral_all", "dpph", "neutral", "all", False),
    ("dpph_author_exclusions", "dpph", "neutral", "author_exclusions", False),
    ("dpph_neutral_monomers", "dpph", "neutral", "monomers", False),
    ("dpph_per_phenolic_oh", "dpph", "neutral", "all", True),
    ("dpph_carboxylate_diagnostic", "dpph", "carboxylate", "all", False),
    ("dpph_exclude_imported", "dpph", "neutral", "exclude_imported", False),
    ("abts_neutral_diagnostic", "abts", "neutral", "all", False),
    ("abts_carboxylate_all", "abts", "carboxylate", "all", False),
    ("abts_carboxylate_author_exclusions", "abts", "carboxylate", "author_exclusions", False),
    ("abts_carboxylate_monomers", "abts", "carboxylate", "monomers", False),
    ("abts_carboxylate_monomers_exclude5", "abts", "carboxylate", "monomers_exclude5", False),
    ("abts_per_phenolic_oh", "abts", "carboxylate", "all", True),
]
results = {}
for key, assay, state, subset, per_oh in definitions:
    included = [r for r in records if r[assay] is not None]
    if subset.startswith("monomers"):
        included = [r for r in included if r["id"] <= 14]
    if subset == "monomers_exclude5":
        included = [r for r in included if r["id"] != 5]
    if subset == "author_exclusions":
        included = [r for r in included if r["id"] not in (5, 16)]
    if subset == "exclude_imported":
        included = [r for r in included if r["id"] not in (15, 16)]
    points = [(r["id"], r[state], r[assay] / (r["phenolic_oh"] if per_oh else 1)) for r in included]
    results[key] = {"assay": assay, "bde_state": state, "subset": subset,
                    "outcome_normalization": "per_phenolic_oh_reanalysis" if per_oh else "per_molecule_as_published",
                    **analyze(points)}

abts_prose = next(text(p) for p in source.findall('.//sec[@id="sec2dot2-molecules-24-01794"]/p') if '1.45 and 1.44' in text(p))
abts_prose_match = re.search(r'RDI values ([0-9.]+) and ([0-9.]+) for 5 and 6', abts_prose)
assert abts_prose_match is not None
prose_replacements = {5: F(abts_prose_match.group(1)), 6: F(abts_prose_match.group(2))}
post_hoc_abts = {
    "status": "post_hoc_source_discrepancy_sensitivity; not a prespecified main analysis or a correction",
    "changed_rows": [{"compound_id": i, "table_abts_rdi": float(next(r['abts'] for r in records if r['id'] == i)), "prose_abts_rdi": float(v)} for i, v in prose_replacements.items()],
    "held_fixed": "All 14 ABTS observations and carboxylate-mapped BDE values; change only response means for compounds 5 and 6; do not transfer table confidence intervals to prose values",
    "source": {"section": "2.2", "xml_id": "sec2dot2-molecules-24-01794", "publisher_pdf_page": 6, "table_pdf_page": 4},
    "result": analyze([(r['id'], r['carboxylate'], prose_replacements.get(r['id'], r['abts'])) for r in records if r['abts'] is not None]),
    "interpretation": "The association is sensitive to which conflicting source values are used. Neither representation is experimentally adjudicated. This does not establish how the authors calculated their published R2."
}

thermal_source = ET.parse(EVIDENCE / "PMC9331396.xml")
t = thermal_source.find('.//table-wrap[@id="polymers-14-03000-t002"]')
temperatures = [int(text(c)) for c in t.findall('.//thead/tr')[1]]
values = rows(thermal_source, 'polymers-14-03000-t002')
thermal = [{"temperature_c": temp, "initial_mass_loss_pct": float(F(initial)),
            "sulfated_mass_loss_pct": float(F(sulfated)),
            "sulfated_minus_initial_percentage_points": float(F(sulfated)-F(initial)),
            "initial_residue_pct": float(100-F(initial)), "sulfated_residue_pct": float(100-F(sulfated))}
           for temp, initial, sulfated in zip(temperatures, values[0][1:], values[1][1:])]
gpc_source = rows(thermal_source, 'polymers-14-03000-t001')
gpc = [{"row_id": int(r[0]), "catalyst": r[1],
        "mn_ratio_to_initial": float(F(r[4]) / F(gpc_source[0][4])),
        "mw_ratio_to_initial": float(F(r[5]) / F(gpc_source[0][5])),
        "pdi_ratio_to_initial": float(F(r[6]) / F(gpc_source[0][6]))} for r in gpc_source[1:]]

answer = {
    "schema_version": 1, "date": "2026-09-05",
    "method": "Independent source XML parser. Fraction arithmetic for centered OLS, R², average-tie ranks and leave-one-out predictions; Decimal square root at 50-digit precision. No release-code imports or generated-data reads.",
    "source_hashes": {p: hashlib.sha256((EVIDENCE / p).read_bytes()).hexdigest() for p in ['PMC6539611.xml', 'PMC9331396.xml']},
    "interpretation_limits": ["Rounded published aggregate values; no replicate-level uncertainty or confidence intervals for these reanalyses.", "Rows are purposefully selected related compounds, not independent random population draws.", "Leave-one-compound-out is internal exploratory predictive validation, not an external test or prospective chemical-series generalization.", "Post-hoc author exclusions and monomer restrictions change the evaluation set; compare each model to the baseline on its same folds, not raw errors across different subsets.", "A smaller BDE thermodynamically favors H-atom donation; RDI is a steady-state radical-consumption endpoint, not a kinetic rate constant.", "All DFT input values use methanol continuum; selecting carboxylates for aqueous ABTS fixes a protonation-state mismatch only, not solvent or full-speciation differences."],
    "lauberte2019": {"variants": results, "post_hoc_abts_prose_substitution": post_hoc_abts},
    "kazachenko2022": {"thermal": thermal,
        "full_table": {"temperature_count": 13, "supporting_count": 11, "counterexample_temperatures_c": [200, 250]},
        "prose_from_250": {"temperature_count": 12, "supporting_count": 11, "counterexample_temperatures_c": [250]},
        "from_300": {"temperature_count": 11, "supporting_count": 11, "counterexample_temperatures_c": []},
        "prose_800_value_conflict": {"prose_initial": 76.6, "table_initial": 61.9, "prose_sulfated": 61.9, "table_sulfated": 76.6, "classification": "exact reversal of sample assignments between prose and Table 2"},
        "gpc_observed_ratios": gpc,
        "gpc_scope": "Apparent measured-size comparisons, not uniquely identified chain-scission fractions; eluent/calibration and sample transformations need interpretation."},
}
out = ARGS.out
out.parent.mkdir(parents=True, exist_ok=True)
out.write_text(json.dumps(answer, indent=2, ensure_ascii=False) + '\n')
for key, value in results.items():
    print(f"{key}: n={value['n']} r2={value['r2']:.9f} rho={value['spearman_rho']:.9f} LOO_RMSE={value['loo']['linear_rmse']:.9f} mean_RMSE={value['loo']['training_mean_rmse']:.9f} skill={value['loo']['mse_skill_vs_training_mean']:.9f}")
print(out)
