#!/usr/bin/env python3
"""Mechanistic load-sharing statistics at the first failure event, computed from the exported per-member data
(utilisation per unit displacement of the intact structure x displacement at first failure)."""
import os, json, sys
import numpy as np, pandas as pd
from hfx_analysis import load_log, ROOT, PROC

def stats_for_run(dir_rel, dFirst, **kwargs):
    m = pd.read_csv(os.path.join(ROOT, dir_rel, 'members.csv'))
    m = m[m.absent == 0]
    w = m['utilInitial_perUnitDisp'].values * dFirst   # utilisation of every member at the first failure event
    w = np.clip(w, 0, None)
    N = len(w); wmax = w.max() if N else np.nan
    out = {}
    out['ls_mean_over_max'] = w.mean() / wmax if wmax > 0 else np.nan
    out['ls_participation'] = (w.sum() ** 2) / (N * (w ** 2).sum()) if (w ** 2).sum() > 0 else np.nan
    out['ls_frac_ge0.5'] = float((w >= 0.5 * wmax).mean())
    out['ls_frac_ge0.8'] = float((w >= 0.8 * wmax).mean())
    out['ls_cv'] = w.std() / w.mean() if w.mean() > 0 else np.nan
    out['ls_gini'] = gini(w)
    # dominant load path: largest single-member share of the reaction force (axial force per unit displacement / K0)
    Nax = m['NInitial_perUnitDisp'].values
    K0 = kwargs.get('K0', np.nan)
    out['lp_max_share'] = float(np.abs(Nax).max() / K0) if K0 > 0 else np.nan
    # share of the reaction carried by level-0 members crossing the mid-height line y = H/2
    y1, y2 = m['y1'].values, m['y2'].values
    yc = 15.0 / 32.0  # cut line between grid lines of the 16x16-equivalent ordered designs
    cross = ((y1 - yc) * (y2 - yc) < 0)
    lvl = m['level'].values
    tot = np.abs(Nax[cross]).sum()
    out['lp_level0_share'] = float(np.abs(Nax[cross & (lvl == 0)]).sum() / tot) if tot > 0 else np.nan
    out['lp_top3_share'] = float(np.sort(np.abs(Nax[cross]))[-3:].sum() / tot) if tot > 0 else np.nan
    out['ls_frac_tension'] = float((Nax > 1e-3 * np.abs(Nax).max()).mean()) if np.abs(Nax).max() > 0 else np.nan
    # bending share of the fibre strain at the most utilised member: util - axial strain / eps_c
    return out

def gini(x):
    x = np.sort(np.asarray(x, dtype=float)); n = len(x)
    if n == 0 or x.sum() == 0: return np.nan
    return (2 * np.sum(np.arange(1, n + 1) * x) / (n * x.sum())) - (n + 1) / n

if __name__ == '__main__':
    df = load_log()
    rows = []
    for _, r in df[df.status == 'ok'].iterrows():
        try:
            s = stats_for_run(r['dir'], r['dFirst'], K0=r['K0']); s['runId'] = r['runId']; rows.append(s)
        except Exception as e:
            print('skip', r['runId'], e)
    ls = pd.DataFrame(rows).set_index('runId')
    ls.to_csv(os.path.join(PROC, 'load_sharing.csv'))
    d = df.set_index('runId').join(ls)
    from scipy.stats import spearmanr, pearsonr
    print('n =', len(d))
    for k in ['ls_mean_over_max', 'ls_participation', 'ls_frac_ge0.5', 'ls_frac_ge0.8', 'ls_cv', 'ls_gini', 'ls_frac_tension', 'lp_max_share', 'lp_level0_share', 'lp_top3_share']:
        ok = d[[k, 'abruptness']].dropna()
        print('%-20s Spearman rho vs abruptness = %+.3f (p=%.1e), vs prePeakFailures = %+.3f, vs Wf = %+.3f' % (k, spearmanr(ok[k], ok['abruptness'])[0], spearmanr(ok[k], ok['abruptness'])[1], spearmanr(d[[k,'prePeakFailures']].dropna()[k], d[[k,'prePeakFailures']].dropna()['prePeakFailures'])[0], spearmanr(d[[k,'Wf']].dropna()[k], d[[k,'Wf']].dropna()['Wf'])[0]))
    g = d.groupby('designId')[['ls_frac_ge0.8', 'lp_max_share', 'lp_top3_share', 'abruptness', 'prePeakFailures', 'Wf']].mean().sort_values('abruptness')
    pd.set_option('display.width', 200); pd.set_option('display.max_rows', 200)
    print(g.round(3).to_string())