"""Phase C — effect-size harmonisation (protocol §3–4, CLAUDE.md §7). Every transformation is explicit and recorded in the
`formula` column; nothing is converted without the information the formula requires.

Families: pct_premium (binary attribute, log price), semi_elasticity (% change in price per unit of X, log price, X in level),
elasticity (log-log), monetary (level price, level X: currency per unit), partial_r (unit-free, from t and df).
"""
import math
import numpy as np, pandas as pd

AREA_TO_M2 = {"m2": 1.0, "sqm": 1.0, "sqft": 0.09290304, "ft2": 0.09290304, "100sqft": 9.290304, "1000sqft": 92.90304, "ha": 10000.0, "acre": 4046.8564224, "are": 100.0, "100m2": 100.0, "pyeong": 3.3058, "tsubo": 3.3058}
DIST_TO_KM = {"km": 1.0, "m": 0.001, "100m": 0.1, "mile": 1.609344, "miles": 1.609344, "ft": 0.0003048, "1000ft": 0.3048, "minutes": None}


def pct_premium(beta, se, kennedy=True):
    """Binary X in a log-price model: 100[exp(b)-1]; Kennedy (1981) correction 100[exp(b - V/2)-1]; delta-method SE 100 exp(b) se."""
    if kennedy and se is not None and not np.isnan(se):
        est = 100 * (math.exp(beta - se ** 2 / 2) - 1); formula = "100*(exp(b - se^2/2) - 1)"
    else:
        est = 100 * (math.exp(beta) - 1); formula = "100*(exp(b) - 1)"
    se_out = 100 * math.exp(beta) * se if se is not None and not np.isnan(se) else np.nan
    return est, se_out, formula


def semi_elasticity(beta, se, unit_scale=1.0):
    """log price, level X: 100*b per unit of X (after rescaling X to the harmonised unit: b_h = b / scale)."""
    return 100 * beta / unit_scale, 100 * se / unit_scale if se is not None else np.nan, f"100*b/{unit_scale}"


def elasticity(beta, se):
    return beta, se, "b (log-log)"


def partial_r(t, df):
    if t is None or df is None or np.isnan(t) or np.isnan(df) or df <= 2:
        return np.nan, np.nan, ""
    r = t / math.sqrt(t * t + df); se = math.sqrt((1 - r * r) / df)
    return r, se, "t/sqrt(t^2+df); se=sqrt((1-r^2)/df)"


def recover_se(beta, se, t, p, ci_lo, ci_hi, df=None):
    """Recover a standard error transparently; returns (se, formula) or (nan, '')."""
    if se is not None and not np.isnan(se) and se > 0:
        return se, "reported"
    if t is not None and not np.isnan(t) and t != 0:
        return abs(beta / t), "|b/t|"
    if ci_lo is not None and ci_hi is not None and not np.isnan(ci_lo) and not np.isnan(ci_hi):
        return (ci_hi - ci_lo) / (2 * 1.959964), "(ci_hi-ci_lo)/(2*1.96)"
    if p is not None and not np.isnan(p) and 0 < p < 1 and beta != 0:
        from scipy.stats import norm
        z = norm.isf(p / 2); return abs(beta / z), "|b/z(p)| (two-sided normal)"
    return np.nan, ""


UNIT_ALIASES = {"m2": "m2", "sqm": "m2", "sq m": "m2", "sq. m": "m2", "m²": "m2", "square meters": "m2", "square metres": "m2", "square meter": "m2", "sqft": "sqft", "sq ft": "sqft", "sq. ft": "sqft",
                "square feet": "sqft", "ft2": "sqft", "ft²": "sqft", "feet2": "sqft", "100sqft": "100sqft", "1000sqft": "1000sqft", "1,000 sqft": "1000sqft", "acres": "acre", "acre": "acre", "hectare": "ha", "ha": "ha",
                "km": "km", "kilometers": "km", "kilometres": "km", "kilometer": "km", "m": "m", "meters": "m", "metres": "m", "meter": "m", "metre": "m", "100m": "100m", "100 m": "100m", "miles": "miles", "mile": "mile",
                "ft": "ft", "feet": "ft", "foot": "ft", "1000ft": "1000ft", "1,000 ft": "1000ft", "minutes": "minutes", "min": "minutes", "years": "years", "year": "years", "yrs": "years", "count": "count", "number": "count",
                "units": "count", "index": "index", "share": "share", "percent": "percent", "%": "percent", "db": "db", "dba": "db", "db(a)": "db", "ug/m3": "ug/m3", "µg/m3": "ug/m3", "μg/m3": "ug/m3", "currency": "currency"}


AREA_LEAVES = {"living_area", "lot_size", "garden"}
DIST_LEAVES = {"cbd_distance", "transit_distance", "school_distance", "university_distance", "hospital_distance", "highway_distance", "commercial_distance", "airport_distance", "water_proximity", "green_space", "hazard_site"}
COUNT_LEAVES = {"bedrooms", "bathrooms", "storeys", "floor_level", "garage", "parking", "fireplace", "terrace_balcony"}
AREA_UNITS = {"m2": 1.0, "sqft": 0.09290304, "1000sqft": 92.90304, "100sqft": 9.290304, "100m2": 100.0, "acre": 4046.8564224, "ha": 10000.0, "pyeong": 3.3058, "are": 100.0, "1000m2": 1000.0}
DIST_UNITS = {"km": 1.0, "m": 0.001, "100m": 0.1, "miles": 1.609344, "ft": 0.0003048, "1000ft": 0.3048, "100km": 100.0, "10km": 10.0, "100ft": 0.03048}
PLAUSIBLE = {"area_per_m2": 5.0, "dist_per_km": 60.0, "years": 10.0, "count": 60.0, "db": 5.0, "ug/m3": 10.0}


def norm_unit(u):
    """Parse a free-text unit into a canonical token: strips parentheticals, spaces, plurals; recognises scaled units."""
    import re as _re
    u = (u or "").strip().lower(); u = _re.sub(r"\(.*?\)", "", u).strip(); u = u.replace("square", "sq").replace(",", "")
    u = _re.sub(r"\s+", " ", u)
    # scaled units like "1000 sqft", "100 m", "100 m2", "hundred years"
    m = _re.match(r"^(1000|100|10)\s*(sq ?ft|sqft|ft2|m2|m²|sqm|m|ft|km|feet)$", u)
    if m:
        n, base = m.group(1), m.group(2)
        base = {"sq ft": "sqft", "sqft": "sqft", "ft2": "sqft", "m2": "m2", "m²": "m2", "sqm": "m2", "feet": "ft", "ft": "ft", "m": "m", "km": "km"}[base]
        return f"{n}{base}"
    if u in ("hundred_years", "hundred years", "100 years", "centuries", "century"): return "hundred_years"
    if u.startswith("ping") or u.startswith("pyeong") or u.startswith("tsubo"): return "pyeong"
    if u in ("sq ft", "sqft", "sq. ft", "ft2", "ft²", "feet2", "sq feet", "sqf"): return "sqft"
    if u in ("m2", "m²", "sqm", "sq m", "sq. m", "sq meters", "sq metres", "sq meter", "square meters", "square metres", "sq.m"): return "m2"
    if u in ("acre", "acres"): return "acre"
    if u in ("hectare", "hectares", "ha"): return "ha"
    if u in ("km", "kilometer", "kilometers", "kilometre", "kilometres"): return "km"
    if u in ("m", "meter", "meters", "metre", "metres"): return "m"
    if u in ("mile", "miles"): return "miles"
    if u in ("ft", "feet", "foot"): return "ft"
    if u in ("year", "years", "yrs", "yr"): return "years"
    if u in ("count", "number", "units", "n", "integer", "rooms"): return "count"
    if u in ("db", "dba", "db(a)", "decibel", "decibels", "dB"): return "db"
    if u in ("ug/m3", "µg/m3", "μg/m3", "ug/m³", "micrograms per cubic meter"): return "ug/m3"
    if u in ("minutes", "min", "minute"): return "minutes"
    return UNIT_ALIASES.get(u, u)


def unit_scale(tax, unit):
    """Return (scale to the harmonised unit, unit family, note) or (None, None, reason) if not harmonisable for this leaf."""
    if tax in AREA_LEAVES:
        return (AREA_UNITS[unit], "area_per_m2", f"X per {unit}->m2") if unit in AREA_UNITS else (None, None, "area unit not harmonisable")
    if tax in DIST_LEAVES:
        return (DIST_UNITS[unit], "dist_per_km", f"X per {unit}->km") if unit in DIST_UNITS else (None, None, "distance unit not harmonisable (counts, shares, indices excluded)")
    if tax == "age":
        return (1.0, "years", "per year") if unit == "years" else ((100.0, "years", "hundred years->year") if unit == "hundred_years" else (None, None, "age unit not harmonisable"))
    if tax in COUNT_LEAVES:
        return (1.0, "count", "per unit count") if unit == "count" else (None, None, "count expected")
    if tax == "noise":
        return (1.0, "db", "per dB") if unit == "db" else (None, None, "noise unit not dB")
    if tax == "air_pollution":
        return (1.0, "ug/m3", "per ug/m3") if unit == "ug/m3" else (None, None, "pollution unit not ug/m3")
    return (None, None, "no harmonised unit defined for this characteristic")


def harmonise_row(r):
    """r: dict with keys beta, se, t, p, ci_lo, ci_hi, dep_transform (log|level|boxcox|other), x_transform (level|log|binary|
    dummy|other), x_unit, taxonomy_id, n_obs, k_params, var_family (from taxonomy). Returns dict with family, es, es_se, formula."""
    beta = r["beta"]; se, se_formula = recover_se(beta, r.get("se"), r.get("t"), r.get("p"), r.get("ci_lo"), r.get("ci_hi"))
    dep = (r.get("dep_transform") or "").lower(); xt = (r.get("x_transform") or "").lower(); unit = norm_unit(r.get("x_unit"))
    if dep == "per_m2":
        dep = "level"  # price per m2 in levels: monetary family (currency per m2 per unit of X)
    name = (r.get("var_name_printed") or "").lower()
    year_built = r.get("taxonomy_id") == "age" and any(k in name for k in ("built", "construct", "vintage", "year of", "yr built", "baujahr", "année"))
    fam = es = es_se = None; formula = ""
    if dep == "log" and xt in ("binary", "dummy"):
        es, es_se, formula = pct_premium(beta, se); fam = "pct_premium"
    elif dep == "log" and xt == "log":
        es, es_se, formula = elasticity(beta, se); fam = "elasticity"
    elif dep == "log" and xt == "level":
        sc, ufam, note = unit_scale(r.get("taxonomy_id"), unit)
        es, es_se, formula = semi_elasticity(beta, se, 1.0); fam = "semi_elasticity"
        if sc is None:
            formula = f"100*b (unit as reported: {unit or 'not stated'}; {note})"; unit_ok = False; ufam_ = None
        else:
            # b is per reported unit; per harmonised unit = b * (reported units per harmonised unit) = b / sc  ... sc = harmonised units per reported unit
            es, es_se = es / sc, (es_se / sc if es_se is not None and not np.isnan(es_se) else np.nan); formula = f"100*b/{sc:.6g} ({note})"; unit_ok = True; ufam_ = ufam
    elif dep == "level" and xt in ("level", "binary", "dummy"):
        es, es_se, formula, fam = beta, se, "b (currency per unit, as reported)", "monetary"
    else:
        fam = "other"; es, es_se, formula = beta, se, "not harmonised"
    df = (r.get("n_obs") or np.nan) - (r.get("k_params") or 10)
    t = r.get("t") if r.get("t") is not None and not np.isnan(r.get("t")) else (beta / se if se and se > 0 else np.nan)
    pr, pr_se, pr_formula = partial_r(t, df)
    flags = []
    if year_built: flags.append("year_built_not_age")
    if fam == "semi_elasticity" and not locals().get("unit_ok", False): flags.append("unit_not_harmonisable")
    if fam == "semi_elasticity" and es is not None and locals().get("ufam_") and abs(es) > PLAUSIBLE.get(locals().get("ufam_"), 100): flags.append("implausible_magnitude")
    if fam == "pct_premium" and es is not None and abs(es) > 200: flags.append("implausible_magnitude")
    if fam == "elasticity" and es is not None and abs(es) > 5: flags.append("implausible_magnitude")
    if xt in ("squared", "other"): flags.append("nonstandard_x_transform")
    return {"family": fam, "es": es, "es_se": es_se, "es_formula": formula, "se_source": se_formula, "partial_r": pr, "partial_r_se": pr_se, "partial_r_formula": pr_formula, "t_used": t, "df_used": df,
            "x_unit_norm": unit, "harm_flags": ";".join(flags)}