# =============================================================
# DISEÑO 1: FACTORIAL SIMPLE EN ARREGLO COMPLETAMENTE AL AZAR
# Efecto del CaCl2 sobre la firmeza poscosecha del tomate de árbol
# (Solanum betaceum Cav.) — día 20
# =============================================================

import numpy as np
import pandas as pd
from scipy import stats
from scipy.stats import studentized_range
from itertools import combinations
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

# ── 1. DATOS ──────────────────────────────────────────────────
# Firmeza (% del peso inicial) a los 20 días poscosecha
# 15 frutos por tratamiento — 3 tratamientos — N = 45

datos = {
    'Control (0 mM)': [
        32.2352, 29.3778, 32.9146, 36.8536, 28.9463,
        28.9464, 37.1065, 33.4535, 27.8874, 32.4415,
        27.9146, 27.9042, 31.0888, 21.3902, 22.2379
    ],
    '570 mM CaCl2': [
        47.3010, 45.1384, 51.5084, 45.6415, 43.2209,
        57.0351, 48.9163, 50.3241, 43.1612, 47.3870,
        50.5324, 44.4752, 51.8034, 47.1169, 48.5999
    ],
    '862 mM CaCl2': [
        49.2321, 60.5205, 51.9379, 47.1345, 55.7837,
        46.3841, 52.9608, 42.9855, 45.8903, 52.9056,
        55.3969, 52.7883, 51.4680, 50.6149, 45.1988
    ]
}

trts = list(datos.keys())

# ── 2. DATAFRAME LARGO ────────────────────────────────────────
registros = []
for trt, valores in datos.items():
    for i, v in enumerate(valores, 1):
        registros.append({'tratamiento': trt, 'fruto': i, 'firmeza': v})
df = pd.DataFrame(registros)
print('=' * 60)
print('DISEÑO COMPLETAMENTE AL AZAR (DCA) — DATOS')
print('=' * 60)
print(df.to_string(index=False))

# ── 3. ESTADÍSTICAS DESCRIPTIVAS ──────────────────────────────
print('\n' + '=' * 60)
print('ESTADÍSTICAS DESCRIPTIVAS POR TRATAMIENTO')
print('=' * 60)
desc = df.groupby('tratamiento')['firmeza'].agg(
    n='count', media='mean', mediana='median',
    DE=lambda x: x.std(ddof=1),
    minimo='min', maximo='max',
    CV=lambda x: x.std(ddof=1) / x.mean() * 100
).round(4)
print(desc.to_string())

# ── 4. ANOVA ──────────────────────────────────────────────────
grupos     = [np.array(datos[t]) for t in trts]
N          = sum(len(g) for g in grupos)
k          = len(grupos)
grand_mean = np.mean(np.concatenate(grupos))
means      = [g.mean() for g in grupos]
ni         = [len(g) for g in grupos]

SS_trt = sum(n * (m - grand_mean)**2 for n, m in zip(ni, means))
SS_err = sum(np.sum((g - m)**2) for g, m in zip(grupos, means))
SS_tot = SS_trt + SS_err

df_trt = k - 1
df_err = N - k
df_tot = N - 1

MS_trt = SS_trt / df_trt
MS_err = SS_err / df_err
F_val  = MS_trt / MS_err
p_val  = stats.f.sf(F_val, df_trt, df_err)
F_crit = stats.f.ppf(0.95, df_trt, df_err)

print('\n' + '=' * 60)
print('TABLA DE ANÁLISIS DE VARIANZA (ANOVA)')
print('=' * 60)
print(f"{'Fuente':<28} {'SC':>10} {'gl':>4} {'CM':>10} {'F':>8} {'p-valor':>12}")
print('-' * 76)
print(f"{'Tratamientos (CaCl2)':<28} {SS_trt:>10.4f} {df_trt:>4} "
      f"{MS_trt:>10.4f} {F_val:>8.4f} {p_val:>12.6e}")
print(f"{'Error experimental':<28} {SS_err:>10.4f} {df_err:>4} "
      f"{MS_err:>10.4f} {'—':>8} {'—':>12}")
print(f"{'Total':<28} {SS_tot:>10.4f} {df_tot:>4}")
print(f'\nF crítico (α=0.05): {F_crit:.4f}')

# ── 5. COMPARACIONES POR PARES DE TUKEY ───────────────────────
print('\n' + '=' * 60)
print('COMPARACIONES POR PARES DE TUKEY')
print('(p-valor ajustado por distribución del rango estudentizado)')
print('=' * 60)
print(f"{'Comparación':<42} {'Diferencia':>10} {'q obs':>8} {'p-adj':>14}")
print('-' * 78)

for i, j in combinations(range(k), 2):
    diff  = means[i] - means[j]
    se    = np.sqrt(MS_err / ni[i])
    q_obs = abs(diff) / se
    p_adj = float(studentized_range.sf(q_obs, k, df_err))
    comp  = f'{trts[i]}  vs  {trts[j]}'
    print(f"{comp:<42} {diff:>10.4f} {q_obs:>8.4f} {p_adj:>14.4e}")

# ── 6. SUPUESTOS ──────────────────────────────────────────────
residuals = np.concatenate([g - m for g, m in zip(grupos, means)])
sw_stat, sw_p   = stats.shapiro(residuals)
bart_stat, bart_p = stats.bartlett(*grupos)

print('\n' + '=' * 60)
print('VERIFICACIÓN DE SUPUESTOS')
print('=' * 60)
print(f'Shapiro-Wilk (residuales)  W = {sw_stat:.4f},  p = {sw_p:.4f}')
print(f'Bartlett (homocedasticidad) K² = {bart_stat:.4f},  p = {bart_p:.4f}')

# ── 7. GRÁFICAS ───────────────────────────────────────────────
fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
fig.suptitle('DCA — Firmeza poscosecha del tomate de árbol (día 20)', fontsize=11)

# Boxplot
colores = ['#B5D4F4', '#9FE1CB', '#F5C4B3']
bp = axes[0].boxplot([datos[t] for t in trts], patch_artist=True, widths=0.45,
                     medianprops=dict(color='#2C2C2A', linewidth=2))
for patch, col in zip(bp['boxes'], colores):
    patch.set_facecolor(col); patch.set_edgecolor('#444441')
axes[0].set_xticks([1, 2, 3])
axes[0].set_xticklabels(trts, fontsize=9)
axes[0].set_ylabel('Firmeza (% peso inicial)')
axes[0].set_title('Boxplots por tratamiento')
axes[0].yaxis.grid(True, linestyle='--', alpha=0.4)

# Q-Q residuales
(osm, osr), (slope, intercept, _) = stats.probplot(residuals, dist='norm')
axes[1].scatter(osm, osr, color='#378ADD', s=18, alpha=0.7)
axes[1].plot([osm[0], osm[-1]],
             [slope*osm[0]+intercept, slope*osm[-1]+intercept],
             color='#D85A30', linewidth=1.5)
axes[1].set_xlabel('Cuantiles teóricos')
axes[1].set_ylabel('Residuales ordenados')
axes[1].set_title(f'Q-Q Normal  W={sw_stat:.4f}  p={sw_p:.4f}')
axes[1].yaxis.grid(True, linestyle='--', alpha=0.3)

plt.tight_layout()
plt.savefig('dca_graficas.png', dpi=150, bbox_inches='tight')
plt.show()
print('\nGráfica guardada: dca_graficas.png')