"""
Statistical analysis for:
"Strength reduction and strength-scaling relations of M40 concrete with
untreated bamboo chips as partial coarse-aggregate replacement"

Reproduces Tables 2-5 and 7 and the regressions reported in the paper:
  - mean, SD, COV and strength retained for each mix (Table 3)
  - one-way ANOVA with Levene and Shapiro-Wilk checks (Table 4)
  - Tukey HSD pairwise comparisons (Table 5)
  - linear regression of strength on bamboo content
  - power-law relations f = a (f'c)^b and k = f / sqrt(f'c) (Table 7)
  - mix proportions per m3 by the absolute-volume method (Table 2)

Requires: Python 3.8+, numpy, scipy
Run:      python statistical_analysis.py
"""
import numpy as np
from scipy import stats

from bamboo_data import LEVELS, COMPRESSIVE, SPLITTING, FLEXURAL


def descriptive(name, data):
    means = np.array([np.mean(data[k]) for k in LEVELS])
    sds = np.array([np.std(data[k], ddof=1) for k in LEVELS])
    print(f"\n=== {name} ===")
    print("Bamboo %   Mean     SD    COV %   Retained %")
    for k, m, s in zip(LEVELS, means, sds):
        print(f"{k:8d} {m:7.2f} {s:6.3f} {100 * s / m:7.1f} {100 * m / means[0]:11.1f}")
    return means, sds


def anova(data):
    groups = [np.asarray(data[k], float) for k in LEVELS]
    values = np.concatenate(groups)
    grand_mean, n, k = values.mean(), len(values), len(groups)
    ss_between = sum(len(g) * (g.mean() - grand_mean) ** 2 for g in groups)
    ss_within = sum(((g - g.mean()) ** 2).sum() for g in groups)
    ss_total = ((values - grand_mean) ** 2).sum()
    df_between, df_within = k - 1, n - k
    ms_between, ms_within = ss_between / df_between, ss_within / df_within
    f_value = ms_between / ms_within
    p_value = stats.f.sf(f_value, df_between, df_within)

    print("One-way ANOVA")
    print(f"  {'Source':<22}{'SS':>11}{'df':>5}{'MS':>11}{'F':>9}{'p':>11}")
    print(f"  {'Between mixes':<22}{ss_between:11.4f}{df_between:5d}{ms_between:11.4f}"
          f"{f_value:9.1f}{p_value:11.2e}")
    print(f"  {'Within mixes (error)':<22}{ss_within:11.4f}{df_within:5d}{ms_within:11.5f}")
    print(f"  {'Total':<22}{ss_total:11.4f}{n - 1:5d}")
    print(f"  eta^2 = {ss_between / ss_total:.3f}")

    if min(len(g) for g in groups) >= 3:
        levene = stats.levene(*groups, center="median")
        shapiro = stats.shapiro(np.concatenate([g - g.mean() for g in groups]))
        print(f"  Levene: W = {levene.statistic:.2f}, p = {levene.pvalue:.3f}; "
              f"Shapiro-Wilk: W = {shapiro.statistic:.3f}, p = {shapiro.pvalue:.3f}")
    else:
        print("  Assumption tests not performed (fewer than 3 specimens per mix).")


def tukey(data):
    result = stats.tukey_hsd(*[data[k] for k in LEVELS])
    ci = result.confidence_interval(0.95)
    print("Tukey HSD (mean difference = lower minus higher bamboo content)")
    for i in range(len(LEVELS)):
        for j in range(i + 1, len(LEVELS)):
            flag = "" if result.pvalue[i, j] < 0.05 else "  (n.s.)"
            print(f"  {LEVELS[i]:2d} vs {LEVELS[j]:2d}%: {result.statistic[i, j]:7.3f} "
                  f"(95% CI {ci.low[i, j]:.3f} to {ci.high[i, j]:.3f}), "
                  f"p = {result.pvalue[i, j]:.4f}{flag}")


def linear_fit(data):
    x = np.concatenate([[k] * len(data[k]) for k in LEVELS])
    y = np.concatenate([data[k] for k in LEVELS])
    r = stats.linregress(x, y)
    print(f"Linear fit: f = {r.intercept:.3f} {r.slope:+.4f} p "
          f"(SE of slope {r.stderr:.4f}), R^2 = {r.rvalue ** 2:.4f}, n = {len(x)}")


def power_law(fc, f, name):
    x, y = np.log(fc), np.log(f)
    r = stats.linregress(x, y)
    t = stats.t.ppf(0.975, len(x) - 2)
    k = f / np.sqrt(fc)
    print(f"{name}: f = {np.exp(r.intercept):.3f} (f'c)^{r.slope:.3f}, "
          f"R^2 = {r.rvalue ** 2:.3f}, 95% CI of exponent "
          f"{r.slope - t * r.stderr:.2f} to {r.slope + t * r.stderr:.2f}")
    print(f"  f / f'c       : {np.round(f / fc, 3)}")
    print(f"  f / sqrt(f'c) : {np.round(k, 2)}  (mean {k.mean():.2f}, "
          f"SD {k.std(ddof=1):.2f}, COV {100 * k.std(ddof=1) / k.mean():.1f}%)")


def mix_proportions(air=0.02):
    """Absolute-volume method; mass ratio cement : sand : stone = 1 : 0.99 : 1.80, w/c = 0.40."""
    specific_gravity = {"cement": 3.15, "sand": 2.50, "stone": 2.65, "water": 1.00}
    ratio = {"cement": 1.00, "sand": 0.99, "stone": 1.80, "water": 0.40}
    litres_per_kg_cement = sum(ratio[m] / specific_gravity[m] for m in ratio)
    cement = 1000 * (1 - air) / litres_per_kg_cement
    mix = {m: cement * ratio[m] for m in ratio}
    stone_volume = mix["stone"] / 2650                  # m3 per m3 of concrete
    bamboo_density = 650                                # kg/m3 (assumed)
    print(f"\n=== Mix proportions per m3 ({100 * air:.0f}% air) ===")
    print("Bamboo %   Cement   Sand   Stone   Bamboo   Water  (kg)")
    for p in LEVELS:
        v = stone_volume * p / 100
        print(f"{p:8d} {mix['cement']:8.0f} {mix['sand']:6.0f} {mix['stone'] - v * 2650:7.0f} "
              f"{v * bamboo_density:8.1f} {mix['water']:7.0f}")


if __name__ == "__main__":
    results = {}
    for name, data in [("Compressive strength (MPa)", COMPRESSIVE),
                       ("Splitting tensile strength (MPa)", SPLITTING),
                       ("Flexural strength (MPa)", FLEXURAL)]:
        means, _ = descriptive(name, data)
        anova(data)
        tukey(data)
        linear_fit(data)
        results[name] = means

    fc = results["Compressive strength (MPa)"]
    fsp = results["Splitting tensile strength (MPa)"]
    fr = results["Flexural strength (MPa)"]
    print("\n=== Strength scaling ===")
    power_law(fc, fr, "Flexural")
    power_law(fc, fsp, "Splitting tensile")
    print(f"Measured flexural / ACI 318 (0.62 sqrt(f'c)): {np.round(fr / (0.62 * np.sqrt(fc)), 2)}")

    mix_proportions()