#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ voyagermon degrau-0 — três medidas, declaradas antes de rodar: A. CODEC — tube -pack (v4s+índice GOP) vs gzip -9 na mesma série, por canal. ε por canal (a lição da 1ª volta ESTÁ no código): dens/T/fluxo caem décadas com a distância radial — bound ABSOLUTO calibrado a 1 UA (convenção BYO-ε) vira ε=1,62 n/cc num sinal de 0,002 n/cc a 119 UA e o canal degenera em saída quase constante (0,13 bit/amostra — o bug catalogado de "tamanho constante"). Canal multiplicativo pede bound RELATIVO: empacotamos log10 e o ε_log certifica erro relativo (ε_rel ≈ ln(10)·ε_log). B : 0,02 e 0,05 nT — acurácia DOCUMENTADA no vy2mgd.txt (validação dual-magnetômetro pós-1989: "accuracy of 0.02 nT - 0.05 nT") V : 0,05 km/s (meio-quantum do formato F7.1 — fidelidade ao arquivo) e o BYO-ε de 28,3 km/s mantido SÓ como teto refutado (40% do sinal a 119 UA) dens : log10, ε_rel 1% e 10% (sem acurácia no arquivo; PLS no heliosheath tem incerteza de dezenas de %) T : log10, ε_rel 1% e 10% lecp : log10, ε_rel 10% (contagem Poisson; fluxo baixo → incerteza alta) B. CRUZAMENTO DA HELIOPAUSA (5/11/2018 = dia 309, Burlaga et al. 2019) — veredito em 3 valores sobre a média diária de |B| (≥4 h válidas): limiar L = ponto médio entre mediana heliosheath (dias 250–305/2018) e mediana VLISM (dias 315–365/2018) — o evento é CONHECIDO; a pergunta não é descobri-lo, é: quantos dias de INDECISÃO cada ε compra? Varredura: ε ∈ {0,02, 0,05, 0,10, 0,13, 0,20} nT (instrumento⊕codec). "Certo acima" sustentado = 3 dias válidos consecutivos com média−ε > L. C. PROVA DE AUSÊNCIA — para 2013–2017 (heliosheath), empacotar o |B| horário do ano e medir o artefato que prova "nenhum cruzamento": índice GOP (min/max) + manifesto, com max(reconstruído)+ε_codec+ε_instr < L. Saída: resultados.json + impressão. stdlib + binário ../tube. """ import csv, gzip, json, math, statistics as st, subprocess from collections import defaultdict from pathlib import Path AQUI = Path(__file__).parent DADOS = AQUI / "dados" PACK = DADOS / "pack" TUBE = AQUI.parent / "tube" GOP = 600 EPS_SWEEP = [0.02, 0.05, 0.10, 0.13, 0.20] JAN_SHEATH = (250, 305) # baseline heliosheath (2018) JAN_VLISM = (315, 365) # baseline VLISM (2018) MIN_H_DIA = 4 SUSTENTA = 3 # dias válidos consecutivos p/ "certo" ANOS_AUSENCIA = range(2013, 2018) BPS_LINK = 160 # downlink de ciência (ordem de grandeza) def roda(cmd): p = subprocess.run([str(c) for c in cmd], capture_output=True, text=True, timeout=600) if p.returncode != 0: raise RuntimeError(f"{cmd}\n{p.stdout}{p.stderr}") return p.stdout def le_canal(nome): with gzip.open(DADOS / f"canal_{nome}.csv.gz", "rt") as f: return [float(l) for l in f if l.strip()] def byo_eps(x, horas_cal=17520): """P98 de |Δ| na janela de calibração (2 primeiros anos válidos). MANTIDO só para exibir a armadilha: a 1 UA ele superestima o bound do resto da missão (decaimento radial).""" cal = x[:horas_cal] d = sorted(abs(cal[i+1] - cal[i]) for i in range(len(cal) - 1)) e = d[int(len(d) * 0.98)] return e if e > 0 else 1e-6 LOG10_1PCT = math.log10(1.01) # ε_log p/ erro relativo certificado de 1% LOG10_10PCT = math.log10(1.10) # idem 10% # nome -> (transforma_log10?, [(eps, rótulo), ...]) — 1º eps = o "honesto" # usado na régua do link CANAIS_SPEC = { "B": (False, [(0.05, "doc dual-mag teto"), (0.02, "doc dual-mag piso")]), "V": (False, [(0.05, "meio-quantum F7.1"), (28.3, "BYO-ε 1 UA (TETO REFUTADO)")]), "dens": (True, [(LOG10_10PCT, "10% rel"), (LOG10_1PCT, "1% rel")]), "T": (True, [(LOG10_10PCT, "10% rel"), (LOG10_1PCT, "1% rel")]), "lecp": (False, [(LOG10_10PCT, "10% rel")]), # CSV já vem em log10 do fetch } def pack_e_verifica(nome, eps, seq=None, verifica=True): PACK.mkdir(parents=True, exist_ok=True) csv_p = DADOS / f"canal_{nome}.csv.gz" if seq is not None: csv_p = PACK / f"tmp_{nome}.csv" csv_p.write_text("".join(f"{v:.6g}\n" for v in seq)) pref = PACK / f"{nome}_{eps:g}" roda([TUBE, "-pack", csv_p, eps, pref, GOP, 0]) if verifica: v = roda([TUBE, "-verify", str(pref) + ".manifest.json", csv_p, 0]) assert "FALHA" not in v.upper(), v stream = (PACK / (pref.name + ".v4s")).stat().st_size idx = (PACK / (pref.name + ".idx")).stat().st_size man = (PACK / (pref.name + ".manifest.json")).stat().st_size return stream, idx, man def gz9_bytes(raw): return len(raw), len(gzip.compress(raw, 9)) # ---------- B. série diária de 2018 ------------------------------------- def dias_2018(): soma, n = defaultdict(float), defaultdict(int) with gzip.open(DADOS / "serie.csv.gz", "rt") as f: for r in csv.DictReader(f): if r["ano"] == "2018" and r["B"]: d = int(r["dia"]) soma[d] += float(r["B"]) n[d] += 1 return {d: soma[d] / n[d] for d in sorted(soma) if n[d] >= MIN_H_DIA} def veredito_cruzamento(dias, L, eps): """último dia CERTO abaixo, primeiro dia CERTO acima (sustentado), e nº de dias-calendário indecidíveis entre eles.""" ds = sorted(dias) def certo_acima(i): seq = [dias[d] for d in ds[i:i + SUSTENTA]] return len(seq) == SUSTENTA and all(v - eps > L for v in seq) def certo_abaixo(i): seq = [dias[d] for d in ds[i:i + SUSTENTA]] return len(seq) == SUSTENTA and all(v + eps < L for v in seq) ult_abaixo = max((ds[i] for i in range(len(ds)) if certo_abaixo(i) and ds[i] < 309 + 20), default=None) prim_acima = min((ds[i] for i in range(len(ds)) if certo_acima(i) and ds[i] > 250), default=None) if ult_abaixo is None or prim_acima is None: return None return ult_abaixo, prim_acima, prim_acima - ult_abaixo - 1 def main(): res = {} # ---- A. codec por canal ------------------------------------------- print("== A. codec (tube v4s+idx vs gzip -9, mesma série) ==") res["codec"] = {} for nome, (translog, eps_lista) in CANAIS_SPEC.items(): x = le_canal(nome) if translog: x = [math.log10(v) for v in x if v > 0] seq_txt = "".join(f"{v:.6g}\n" for v in x).encode() raw_csv, gz = gz9_bytes(seq_txt) bruto = 8 * len(x) linha = {"n": len(x), "log10": translog, "bruto_f64_B": bruto, "csv_B": raw_csv, "gzip9_B": gz, "eps": {}} for e, rotulo in eps_lista: s, i, m = pack_e_verifica(nome, e, seq=x) bits = 8 * (s + i) / len(x) linha["eps"][f"{e:g}"] = { "rotulo": rotulo, "stream_B": s, "idx_B": i, "manifest_B": m, "bits_por_amostra": round(bits, 2), "vs_gzip": round(gz / (s + i), 2), "vs_f64": round(bruto / (s + i), 2)} print(f" {nome:5s} ε={e:<8.4g} [{rotulo}] n={len(x):>7,} " f"tube {s+i:>8,} B ({bits:.2f} b/am) gzip9 {gz:>8,} B " f"→ {gz/(s+i):5.2f}× gzip, {bruto/(s+i):5.1f}× f64") res["codec"][nome] = linha # ---- B. cruzamento ------------------------------------------------ print("\n== B. heliopausa 2018: dias de indecisão por ε ==") dias = dias_2018() sheath = [v for d, v in dias.items() if JAN_SHEATH[0] <= d <= JAN_SHEATH[1] and d < 306] vlism = [v for d, v in dias.items() if JAN_VLISM[0] <= d <= JAN_VLISM[1]] m_sh, m_vl = st.median(sheath), st.median(vlism) L = (m_sh + m_vl) / 2 print(f" mediana heliosheath={m_sh:.3f} nT, VLISM={m_vl:.3f} nT, " f"L={L:.3f} nT, degrau={m_vl-m_sh:.3f} nT") res["cruzamento"] = {"mediana_sheath_nT": round(m_sh, 3), "mediana_vlism_nT": round(m_vl, 3), "L_nT": round(L, 3), "por_eps": {}} for e in EPS_SWEEP: v = veredito_cruzamento(dias, L, e) if v is None: print(f" ε={e:g}: INDECIDÍVEL — nunca há {SUSTENTA} dias certos") res["cruzamento"]["por_eps"][str(e)] = None else: ua, pa, ind = v print(f" ε={e:g}: último certo-abaixo dia {ua}, primeiro " f"certo-acima dia {pa} → {ind} dias indecidíveis") res["cruzamento"]["por_eps"][str(e)] = { "ultimo_abaixo": ua, "primeiro_acima": pa, "dias_indecidiveis": ind} # ---- C. prova de ausência ----------------------------------------- print("\n== C. prova de ausência de cruzamento (heliosheath 2013–2017) ==") res["ausencia"] = {} with gzip.open(DADOS / "serie.csv.gz", "rt") as f: por_ano = defaultdict(list) for r in csv.DictReader(f): if r["B"] and int(r["ano"]) in ANOS_AUSENCIA: por_ano[int(r["ano"])].append(float(r["B"])) eps_c, eps_i = 0.02, 0.05 # codec ⊕ instrumento (soma sound) for ano in sorted(por_ano): seq = por_ano[ano] s, i, m = pack_e_verifica(f"B{ano}", eps_c, seq=seq) # reconstrói e checa o predicado da ausência pref = PACK / f"B{ano}_{eps_c:g}" rec_p = PACK / f"B{ano}.rec.csv" roda([TUBE, "-recon", str(pref) + ".manifest.json", rec_p]) mx = max(float(l) for l in rec_p.read_text().splitlines() if l.strip()) prova = mx + eps_c + eps_i < L cob = len(seq) / 8760 # a prova só cobre as horas MEDIDAS print(f" {ano}: {len(seq):>5,} h ({100*cob:4.1f}% do ano) " f"máx_rec={mx:.3f} nT máx+ε={mx+eps_c+eps_i:.3f} " f"{'<' if prova else '≥'} L={L:.3f} → idx {i} B + manifest {m} B " f"({'AUSÊNCIA PROVADA nas horas medidas' if prova else 'NÃO PROVA'})") res["ausencia"][str(ano)] = {"horas": len(seq), "cobertura_ano": round(cob, 3), "max_rec_nT": round(mx, 3), "prova": prova, "idx_B": i, "manifest_B": m} # ---- D. régua do link --------------------------------------------- # 1º ε de cada canal (o "honesto": doc p/ B, quantum p/ V, 10% rel p/ resto) tot_tube = 0 for nome, (_, eps_lista) in CANAIS_SPEC.items(): v = res["codec"][nome]["eps"][f"{eps_lista[0][0]:g}"] tot_tube += v["stream_B"] + v["idx_B"] tot_gz = sum(c["gzip9_B"] for c in res["codec"].values()) masc = sum((DADOS / f"presenca_{n}.bin.gz").stat().st_size for n in CANAIS_SPEC) seg_tube = (tot_tube + masc) * 8 / BPS_LINK seg_gz = (tot_gz + masc) * 8 / BPS_LINK print(f"\n== D. régua do link ({BPS_LINK} bps) ==") print(f" missão inteira, 5 canais: tube {tot_tube:,} B (+máscaras {masc:,} B) " f"= {seg_tube/3600:.1f} h de link; gzip9 {tot_gz:,} B = {seg_gz/3600:.1f} h") res["link"] = {"tube_B": tot_tube, "gzip9_B": tot_gz, "mascaras_B": masc, "bps": BPS_LINK, "tube_h_link": round(seg_tube / 3600, 1), "gzip9_h_link": round(seg_gz / 3600, 1)} (AQUI / "resultados.json").write_text( json.dumps(res, ensure_ascii=False, indent=1)) print("\nresultados.json gravado") if __name__ == "__main__": main()