#!/usr/bin/env python3
"""
Orrery Lab · H4 · Índice cíclico de Barbault × guerra (análise confirmatória)

Uso:
  python3 h4_analise.py H4_indice.csv UCDP_ACD.csv COW_InterState.csv

Só usa a biblioteca padrão do Python (sem numpy), para ser reprodutível em qualquer máquina.
Tudo o que este script faz está descrito no protocolo H4.md e não pode ser alterado depois do selo.
"""
import csv, sys, math

K_MIN, K_MAX = 30, 600          # desfasamentos do nulo, em anos (± , excluindo |k|<30)
ALFA = 0.025                     # por primário (Bonferroni, 2 primários)

# ---------- utilitários estatísticos ----------
def ranks(v):
    o = sorted(range(len(v)), key=lambda i: v[i])
    r = [0.0] * len(v); i = 0
    while i < len(o):
        j = i
        while j + 1 < len(o) and v[o[j + 1]] == v[o[i]]:
            j += 1
        m = (i + j) / 2 + 1
        for k in range(i, j + 1):
            r[o[k]] = m
        i = j + 1
    return r

def pearson(a, b):
    n = len(a); ma = sum(a) / n; mb = sum(b) / n
    sab = sum((x - ma) * (y - mb) for x, y in zip(a, b))
    saa = sum((x - ma) ** 2 for x in a); sbb = sum((y - mb) ** 2 for y in b)
    return sab / math.sqrt(saa * sbb) if saa > 0 and sbb > 0 else float('nan')

def detrend(years, v):
    """Resíduos da regressão linear (OLS) de v no ano."""
    n = len(v); mx = sum(years) / n; my = sum(v) / n
    sxx = sum((x - mx) ** 2 for x in years)
    b = sum((x - mx) * (y - my) for x, y in zip(years, v)) / sxx
    a = my - b * mx
    return [y - (a + b * x) for x, y in zip(years, v)]

def stat(years, idx, out):
    """Spearman entre os resíduos destendenciados do índice e do resultado."""
    return pearson(ranks(detrend(years, idx)), ranks(detrend(years, out)))

# ---------- leitura ----------
def ler_indice(path):
    I = {}
    with open(path, newline='') as f:
        for r in csv.DictReader(f):
            I[int(r['ano'])] = float(r['indice_media_anual'])
    return I

def get(r, *nomes):
    for n in nomes:
        if n in r and r[n] not in (None, ''):
            return r[n]
    raise KeyError(nomes)

def resultado_ucdp(path):
    """P1: n.º de conflitos distintos com intensity_level == 2 (guerra, ≥1000 mortes) por ano."""
    guerras = {}; anos = set()
    with open(path, newline='', encoding='utf-8-sig') as f:
        for r in csv.DictReader(f):
            y = int(get(r, 'year', 'Year')); anos.add(y)
            if int(get(r, 'intensity_level', 'IntensityLevel')) == 2:
                guerras.setdefault(y, set()).add(get(r, 'conflict_id', 'ConflictID'))
    a0, a1 = min(anos), max(anos)
    return {y: float(len(guerras.get(y, ()))) for y in range(a0, a1 + 1)}

def resultado_cow(path):
    """P2: log10(1 + mortes em combate) por ano, guerras interestatais COW v4.0.
    As mortes de cada participante são repartidas por igual pelos anos civis em que
    esteve em combate (período 1 e, se existir, período 2). -7 (em curso) = 2007.
    Participantes com mortes desconhecidas (-9) contam 0."""
    D = {y: 0.0 for y in range(1816, 2008)}
    with open(path, newline='', encoding='latin-1') as f:
        for r in csv.DictReader(f):
            d = float(get(r, 'BatDeath', 'BatDeaths'))
            if d < 0:
                continue
            anos = []
            for s, e in (('StartYear1', 'EndYear1'), ('StartYear2', 'EndYear2')):
                a = int(get(r, s)); b = int(get(r, e))
                if a < 0:
                    continue
                if b == -7:
                    b = 2007
                if b < 0:
                    b = a
                anos.extend(range(a, b + 1))
            anos = sorted(set(y for y in anos if 1816 <= y <= 2007))
            for y in anos:
                D[y] += d / len(anos)
    return {y: math.log10(1 + v) for y, v in D.items()}

# ---------- teste ----------
def teste(nome, I, R):
    anos = sorted(R)
    out = [R[y] for y in anos]
    obs = stat(anos, [I[y] for y in anos], out)
    nulos = []
    for k in list(range(-K_MAX, -K_MIN + 1)) + list(range(K_MIN, K_MAX + 1)):
        nulos.append(stat(anos, [I[y + k] for y in anos], out))
    p = (1 + sum(1 for x in nulos if x <= obs)) / (1 + len(nulos))
    nulos.sort()
    print(f"{nome}: anos {anos[0]}–{anos[-1]} (n={len(anos)})")
    print(f"  rho observado = {obs:+.4f}  (esperado negativo: índice baixo ↔ mais guerra)")
    print(f"  nulo: {len(nulos)} desfasamentos; mediana {nulos[len(nulos)//2]:+.4f}; "
          f"percentil 2,5% {nulos[int(0.025*len(nulos))]:+.4f}")
    print(f"  p unilateral = {p:.4f}  →  {'SUCESSO' if p < ALFA else 'sem efeito'} (α = {ALFA})")
    return obs, p

if __name__ == '__main__':
    I = ler_indice(sys.argv[1])
    r1 = teste('P1 · UCDP/PRIO guerras por ano', I, resultado_ucdp(sys.argv[2]))
    r2 = teste('P2 · COW mortes interestatais (log10)', I, resultado_cow(sys.argv[3]))
    print('H4 CONFIRMADA' if min(r1[1], r2[1]) < ALFA else 'H4 NÃO CONFIRMADA')