"""H7 — Digital-assortment complexity indicator with validation (V1 bootstrap stability, V2 store-reshuffle null, V4 community stability).
V3 (incremental information vs external covariates) is computed in src/models/h2_scale.py once external data exist.
Matrix M_rc = 1[RAA_rc >= 1], scheme M on S_R. Indicators: diversity (row sums), ECI (method of reflections / eigenvector of
D_r^-1 M D_c^-1 M^T, sign fixed so that corr(ECI, diversity) >= 0), Fitness (Tacchella et al. 2012, 200 iterations).
Also: category relatedness (Hidalgo proximity across STORES), category network, Louvain communities.
"""
from __future__ import annotations
import json
import numpy as np, pandas as pd, networkx as nx
from scipy import stats
from src.common import ROOT, SEED, get_logger, ledger, register_output, write_json, PROC_DIR, EXP_DIR, P
from src.models.h1_raa import weights, table, raa

log = get_logger("models.h7")
PROC = PROC_DIR


def eci(M):
    kr = M.sum(1); kc = M.sum(0)
    ok_r = kr > 0; ok_c = kc > 0
    Mm = M[np.ix_(ok_r, ok_c)]; kr = Mm.sum(1); kc = Mm.sum(0)
    Mt = (Mm / kr[:, None]) @ (Mm.T / kc[:, None]).T if False else (Mm / kr[:, None]) @ (Mm / kc).T
    w, v = np.linalg.eig(Mt); order = np.argsort(-w.real); vec = v[:, order[1]].real
    vec = (vec - vec.mean()) / vec.std()
    if np.corrcoef(vec, kr)[0, 1] < 0: vec = -vec
    out = np.full(M.shape[0], np.nan); out[ok_r] = vec; return out


def lambda2(M):
    kr = M.sum(1); kc = M.sum(0); ok_r = kr > 0; ok_c = kc > 0; Mm = M[np.ix_(ok_r, ok_c)]
    Mt = (Mm / Mm.sum(1)[:, None]) @ (Mm / Mm.sum(0)).T; w = np.sort(np.linalg.eigvals(Mt).real)[::-1]; return float(w[1]) if len(w) > 1 else np.nan


def spectral_radius(M):
    R, C = M.shape; A = np.zeros((R + C, R + C)); A[:R, R:] = M; A[R:, :R] = M.T; return float(np.max(np.linalg.eigvalsh(A)))


def fitness(M, iters=200):
    F = np.ones(M.shape[0]); Q = np.ones(M.shape[1])
    for _ in range(iters):
        F1 = M @ Q; F1 = F1 / F1.mean()
        Q1 = 1.0 / (M.T @ (1.0 / np.maximum(F, 1e-12))); Q1 = Q1 / Q1.mean(); F, Q = F1, Q1
    return F, Q


def matrix(d, regions, cats):
    rcode = d.region_final.map({r: i for i, r in enumerate(regions)}).values; ccode = d.category.map({c: i for i, c in enumerate(cats)}).values
    t = table(rcode, ccode, weights(d, "M"), len(regions), len(cats)); A = raa(np.where(t == 0, 1e-12, t))
    return (A >= 1).astype(float), t, A


def proximity_stores(sc):
    pres = sc.assign(v=1).pivot_table(index="store_id", columns="category", values="v", fill_value=0)
    co = pres.T @ pres; k = np.diag(co).astype(float)
    phi = co.values / np.maximum(np.minimum.outer(k, k), 1); np.fill_diagonal(phi, 0)
    return pd.DataFrame(phi, index=pres.columns, columns=pres.columns), pres


def main():
    rng = np.random.default_rng(SEED)
    sc = pd.read_parquet(PROC / "store_category.parquet"); d = sc[sc.region_final != "Inconnue"].reset_index(drop=True)
    regions = sorted(d.region_final.unique()); cats = sorted(d.category.unique())
    M, t, A = matrix(d, regions, cats)
    div = M.sum(1); ubiq = M.sum(0); e = eci(M); F, Q = fitness(M)
    base = pd.DataFrame({"region": regions, "diversity": div, "eci": e, "fitness": F, "n_stores": d.groupby("region_final").store_id.nunique().reindex(regions).values})
    # V1: store bootstrap
    stores = d.store_id.unique(); sidx = d.store_id.map({s: i for i, s in enumerate(stores)}).values
    rows_by_store = [np.where(sidx == i)[0] for i in range(len(stores))]
    rhos, rhos_f, divs = [], [], []
    for b in range(500):
        draw = rng.integers(0, len(stores), len(stores)); idx = np.concatenate([rows_by_store[i] for i in draw])
        Mb, _, _ = matrix(d.iloc[idx], regions, cats); eb = eci(Mb); Fb, _ = fitness(Mb)
        m = ~np.isnan(eb) & ~np.isnan(e)
        rhos.append(stats.spearmanr(e[m], eb[m]).correlation); rhos_f.append(stats.spearmanr(F, Fb).correlation); divs.append(stats.spearmanr(div, Mb.sum(1)).correlation)
    # V2: store-level region reshuffle null
    sreg = d.drop_duplicates("store_id").set_index("store_id").region_final.reindex(stores).values
    null_l2, null_sr, null_corr_div = [], [], []
    for k in range(2000):
        perm = rng.permutation(sreg); dd = d.assign(region_final=perm[sidx])
        Mp, _, _ = matrix(dd, regions, cats); ep = eci(Mp)
        null_l2.append(lambda2(Mp)); null_sr.append(spectral_radius(Mp)); m = ~np.isnan(ep); null_corr_div.append(stats.spearmanr(ep[m], Mp.sum(1)[m]).correlation)
    obs_corr_div = float(stats.spearmanr(e, div).correlation); obs_l2 = lambda2(M); obs_sr = spectral_radius(M)
    p_l2 = float((np.sum(np.array(null_l2) >= obs_l2) + 1) / 2001); p_sr = float((np.sum(np.array(null_sr) >= obs_sr) + 1) / 2001)
    # relatedness network and communities (V4)
    phi, pres = proximity_stores(sc)
    G = nx.Graph()
    for i, a in enumerate(cats):
        for j, b in enumerate(cats):
            if j > i and phi.loc[a, b] > 0: G.add_edge(a, b, weight=float(phi.loc[a, b]))
    try:
        import community as louvain
        parts = [louvain.best_partition(G, weight="weight", random_state=s) for s in range(200)]
        from sklearn.metrics import normalized_mutual_info_score as nmi
        labs = [[p[c] for c in cats] for p in parts]
        nmis = [nmi(labs[0], labs[k]) for k in range(1, 200)]
        nmods = [len(set(p.values())) for p in parts]
        v4 = {"nmi_median_vs_run0": round(float(np.median(nmis)), 3), "nmi_p10": round(float(np.percentile(nmis, 10)), 3), "n_communities_mode": int(stats.mode(nmods, keepdims=False).mode),
              "communities_run0": {c: int(parts[0][c]) for c in cats}}
    except Exception as ex:  # noqa: BLE001
        v4 = {"error": str(ex)[:200]}
    # maximum spanning tree for the figure
    mst = nx.maximum_spanning_tree(G, weight="weight")
    nx.write_weighted_edgelist(G, PROC / "category_proximity_edges.txt"); phi.to_csv(PROC / "category_proximity.csv")
    base.to_csv(PROC / "h7_region_indicators.csv", index=False); register_output(PROC / "h7_region_indicators.csv", "h7")
    pd.DataFrame(A, index=regions, columns=cats).to_csv(PROC / "h7_raa_M.csv"); pd.DataFrame(M, index=regions, columns=cats).to_csv(PROC / "h7_M.csv")
    res = {"regions": len(regions), "categories": len(cats), "M_density": round(float(M.mean()), 3), "diversity": dict(zip(regions, div.astype(int).tolist())),
           "eci": {r: round(float(x), 3) for r, x in zip(regions, e)}, "fitness": {r: round(float(x), 3) for r, x in zip(regions, F)},
           "V1_bootstrap_rank_stability": {"eci_rho_median": round(float(np.nanmedian(rhos)), 3), "eci_rho_p10": round(float(np.nanpercentile(rhos, 10)), 3), "fitness_rho_median": round(float(np.nanmedian(rhos_f)), 3), "diversity_rho_median": round(float(np.nanmedian(divs)), 3), "pass_eci": bool(np.nanmedian(rhos) >= 0.8)},
           "V2_reshuffle_null": {"obs_lambda2": round(obs_l2, 4), "null_lambda2_p95": round(float(np.percentile(null_l2, 95)), 4), "p_lambda2": p_l2,
                                 "obs_spectral_radius": round(obs_sr, 4), "null_spectral_radius_p95": round(float(np.percentile(null_sr, 95)), 4), "p_spectral_radius": p_sr,
                                 "obs_corr_eci_diversity": round(obs_corr_div, 3), "null_corr_eci_diversity_p5_p95": [round(float(np.nanpercentile(null_corr_div, 5)), 3), round(float(np.nanpercentile(null_corr_div, 95)), 3)],
                                 "pass_lambda2": bool(p_l2 < 0.05), "pass_spectral_radius": bool(p_sr < 0.05), "pass": bool(p_l2 < 0.05 and p_sr < 0.05)},
           "V4_community_stability": v4, "mst_edges": [(a, b, round(w, 3)) for a, b, w in mst.edges(data="weight")]}
    write_json(res, EXP_DIR / "h7_complexity.json", "h7"); ledger("models.h7_complexity", **{k: v for k, v in res.items() if k not in ("mst_edges",)})
    log.info(json.dumps({k: v for k, v in res.items() if k in ("M_density", "V1_bootstrap_rank_stability", "V2_reshuffle_null", "V4_community_stability")}, ensure_ascii=False)[:900])


if __name__ == "__main__":
    main()