1"""Phase 2 orchestrator: build the three panels, write processed files, coverage tables, missingness diagnostics,2temporal-distribution figures and the data dictionary. Usage: python -m wp35.harmonization.build_panels3No predictor-outcome association is computed here.4"""5from __future__ import annotations6import json7from datetime import datetime, timezone8from pathlib import Path9import matplotlib10matplotlib.use("Agg")11import matplotlib.pyplot as plt12import numpy as np13import pandas as pd14from wp35.harmonization import seshat as S, dplace as D, modern as M1516ROOT = Path(__file__).resolve().parents[3]17PROC = ROOT / "data" / "processed"; META = ROOT / "data" / "metadata"; DIAG = ROOT / "results" / "diagnostics"18FIG = ROOT / "paper" / "figures"; TAB = ROOT / "paper" / "tables"19for p in (PROC, META, DIAG, FIG, TAB):20 p.mkdir(parents=True, exist_ok=True)2122K_BLOCKS = {23 "admin": ["full_time_bureaucrat", "specialized_government_building", "merit_promotion", "examination_system", "administrative_level"],24 "rules_money": ["formal_legal_code", "judge", "court", "professional_lawyer", "indigenous_coin", "foreign_coin",25 "paper_currency", "token", "precious_metal", "article"],26 "comm_infra": ["courier", "postal_station", "general_postal_service", "road", "bridge", "canal", "port"],27 "records": ["written_record", "script", "phonetic_alphabetic_writing", "non_phonetic_writing",28 "lists_tables_and_classification", "calendar", "practical_literature", "nonwritten_record", "mnemonic_device"],29}30STANDARDS = ["length_measurement_system", "area_measurement_system", "volume_measurement_system", "weight_measurement_system", "time_measurement_system"]31D_VARS = ["polity_population", "polity_territory", "population_of_the_largest_settlement", "settlement_hierarchy"]32PROD_INFRA = ["irrigation_system", "drinking_water_supply_system", "market", "food_storage_site", "mines_or_quarry"]33MILITARY = ["professional_soldier", "professional_military_officer", "military_level"]343536def seshat_panels():37 out = S.build_all()38 facts, agg, meta, crisis, pre = out["facts"], out["polity_agg"], out["polity_meta"], out["crisis"], out["precrisis_agg"]39 facts.to_parquet(PROC / "seshat_facts.parquet", index=False)40 agg.to_parquet(PROC / "seshat_polity_agg.parquet", index=False)41 meta.to_parquet(PROC / "seshat_polity_meta.parquet", index=False)42 crisis.to_parquet(PROC / "seshat_crisis.parquet", index=False)43 pre.to_parquet(PROC / "seshat_precrisis_agg.parquet", index=False)44 wide = S.wide(agg, "value_mean"); wide_log = S.wide(agg, "logv_mean")45 wide.to_parquet(PROC / "seshat_polity_wide.parquet")46 # coverage table: per variable, polities coded, facts, dated facts, inferred share47 cov = agg.groupby("variable").agg(polities_coded=("polity_id", "nunique"), facts=("n_facts", "sum"),48 dated_facts=("n_dated", "sum"), inferred_facts=("n_inferred", "sum"),49 disputed_facts=("n_disputed", "sum")).reset_index()50 sect = facts.drop_duplicates("variable").set_index("variable").section51 cov["section"] = cov.variable.map(sect)52 cov["block"] = cov.variable.map({v: b for b, vs in K_BLOCKS.items() for v in vs}).fillna(53 cov.variable.map({v: "D_scale" for v in D_VARS}).fillna(cov.variable.map({v: "standards" for v in STANDARDS})))54 cov = cov.sort_values(["block", "polities_coded"], ascending=[True, False])55 cov.to_csv(DIAG / "seshat_coverage_by_variable.csv", index=False)56 # Warfare and other sheets: coverage57 warfare_vars = facts[facts.sheet == "Warfare"].variable.unique().tolist()58 # missingness among the pre-registered indicator set, per polity59 ind = sum(K_BLOCKS.values(), []) + D_VARS + STANDARDS + PROD_INFRA + MILITARY60 pres = wide.reindex(columns=ind)61 miss = pd.DataFrame({"indicator": ind, "share_polities_missing": pres.isna().mean().values})62 miss.to_csv(DIAG / "seshat_missingness_by_indicator.csv", index=False)63 per_pol = pd.DataFrame({"polity_id": pres.index, "n_K_indicators_coded": pres[sum(K_BLOCKS.values(), [])].notna().sum(axis=1).values,64 "n_D_coded": pres[D_VARS].notna().sum(axis=1).values})65 per_pol = per_pol.merge(meta[["polity_id", "macro_region", "era", "start_year", "end_year"]], on="polity_id")66 per_pol.to_csv(DIAG / "seshat_missingness_by_polity.csv", index=False)67 # temporal distribution figure (polity time spans by macro-region) and crisis start years68 fig, axes = plt.subplots(1, 2, figsize=(11, 4.2), gridspec_kw={"width_ratios": [2, 1]})69 regs = sorted(meta.macro_region.unique())70 for i, r in enumerate(regs):71 sub = meta[meta.macro_region == r]72 axes[0].hlines([i + np.random.RandomState(1).uniform(-0.35, 0.35) for _ in range(len(sub))], sub.start_year, sub.end_year, lw=1.2, alpha=0.6)73 axes[0].set_yticks(range(len(regs))); axes[0].set_yticklabels([r.replace("_", " ") for r in regs], fontsize=8)74 axes[0].set_xlabel("Year (negative = BCE)"); axes[0].set_title(f"Polity time spans, Polaris-2026 (n = {len(meta)})", fontsize=10)75 axes[1].hist(crisis.crisis_start.dropna(), bins=30, color="#555"); axes[1].set_xlabel("Crisis start year"); axes[1].set_title(f"CrisisDB records (n = {len(crisis)})", fontsize=10)76 fig.tight_layout(); fig.savefig(FIG / "fig_seshat_temporal_coverage.pdf"); plt.close(fig)77 # missingness heatmap for pre-registered indicators78 fig, ax = plt.subplots(figsize=(10, 5))79 order = pres[sum(K_BLOCKS.values(), []) + D_VARS].notna()80 order = order.loc[order.sum(axis=1).sort_values(ascending=False).index]81 ax.imshow(order.values.T, aspect="auto", cmap="Greys", interpolation="nearest")82 ax.set_yticks(range(order.shape[1])); ax.set_yticklabels(order.columns, fontsize=6); ax.set_xlabel("Polities (sorted by number of coded indicators)")83 ax.set_title("Coded (dark) vs missing (white) indicators, Seshat Polaris-2026", fontsize=10)84 fig.tight_layout(); fig.savefig(FIG / "fig_seshat_missingness.pdf"); plt.close(fig)85 # coverage summary numbers86 summary = {87 "polities": int(len(meta)), "regions": int(meta.seshat_region.nunique()), "macro_regions": int(meta.macro_region.nunique()),88 "facts_total": int(len(facts)), "facts_dated": int(facts.dated.sum()), "facts_unparsed": int(facts.unparsed.sum()),89 "facts_inferred": int(facts.inferred.sum()), "facts_disputed_code": int(facts.disputed.sum()),90 "polities_with_population": int(wide.polity_population.notna().sum()), "polities_with_territory": int(wide.polity_territory.notna().sum()),91 "polities_with_any_D": int(pres[D_VARS].notna().any(axis=1).sum()),92 "polities_with_ge50pct_K": int((pres[sum(K_BLOCKS.values(), [])].notna().mean(axis=1) >= 0.5).sum()),93 "crisis_records": int(len(crisis)), "crisis_polities": int(crisis.polity_id.nunique()),94 "crisis_y1_coded": int(crisis.y1_severity.notna().sum()), "crisis_y1b_coded": int(crisis.y1b_severe.notna().sum()),95 "crisis_y1b_positive": int((crisis.y1b_severe == 1).sum()),96 "y3_coded": int(meta.y3_discontinuous.notna().sum()), "y3_positive": int((meta.y3_discontinuous == 1).sum()),97 "duration_censored_ge1900": int(meta.censored.sum()), "warfare_indicators": len(warfare_vars),98 "era_counts": meta.era.value_counts().to_dict(), "macro_region_counts": meta.macro_region.value_counts().to_dict(),99 }100 return summary, cov101102103def dplace_panel():104 out = D.build()105 soc = out["societies"]106 soc.to_parquet(PROC / "dplace_ea_societies.parquet", index=False)107 out["codes"].to_csv(META / "dplace_ea_codes_used.csv", index=False)108 cov = pd.DataFrame({"variable": [c for c in soc.columns if c.startswith("EA")], })109 cov["societies_coded"] = [int(soc[c].notna().sum()) for c in cov.variable]110 cov.to_csv(DIAG / "dplace_coverage_by_variable.csv", index=False)111 fig, ax = plt.subplots(figsize=(6, 3.5))112 ax.hist(soc.main_focal_year.dropna(), bins=40, color="#555"); ax.set_xlabel("Main focal year of the ethnographic description"); ax.set_title(f"D-PLACE EA v3.3 societies (n = {len(soc)})", fontsize=10)113 fig.tight_layout(); fig.savefig(FIG / "fig_dplace_focal_years.pdf"); plt.close(fig)114 return {"societies": int(len(soc)), "with_EA031": int(soc.ea031_code.notna().sum()), "with_EA202": int(soc.ea202_code.notna().sum()),115 "with_EA033": int(soc.EA033.notna().sum()), "with_glottolog_family": int(soc.Family_ID.notna().sum()),116 "families": int(soc.Family_ID.nunique()), "regions": int(soc.region.nunique()),117 "focal_year_range": [int(soc.main_focal_year.min()), int(soc.main_focal_year.max())]}118119120def modern_panel():121 out = M.build(); p = out["panel"]122 p.to_parquet(PROC / "modern_panel.parquet", index=False)123 cov = pd.DataFrame({"variable": [c for c in p.columns if c not in ("country_name", "country_text_id", "country_id", "year")]})124 cov["country_years_coded"] = [int(p[c].notna().sum()) for c in cov.variable]125 cov["first_year"] = [int(p.year[p[c].notna()].min()) if p[c].notna().any() else None for c in cov.variable]126 cov["last_year"] = [int(p.year[p[c].notna()].max()) if p[c].notna().any() else None for c in cov.variable]127 cov["countries"] = [int(p.country_text_id[p[c].notna()].nunique()) for c in cov.variable]128 cov.to_csv(DIAG / "modern_coverage_by_variable.csv", index=False)129 fig, ax = plt.subplots(figsize=(7, 3.5))130 for c, lab in [("maddison_pop", "Maddison population"), ("v2clrspct", "V-Dem K items"), ("e_civil_war", "civil war (COW)"), ("e_coups", "coups"), ("v2regendtype", "regime end type")]:131 ax.plot(p.groupby("year")[c].apply(lambda s: s.notna().sum()), label=lab, lw=1.2)132 ax.set_ylabel("countries coded"); ax.set_xlabel("year"); ax.legend(fontsize=7); ax.set_title("Modern panel coverage by year", fontsize=10)133 fig.tight_layout(); fig.savefig(FIG / "fig_modern_coverage.pdf"); plt.close(fig)134 cw = out["crosswalk"]135 return {"country_years": int(len(p)), "countries": int(p.country_text_id.nunique()), "years": [int(p.year.min()), int(p.year.max())],136 "with_maddison_pop": int(p.maddison_pop.notna().sum()), "countries_matched_maddison": int(cw[cw.in_vdem & cw.in_maddison].shape[0]),137 "vdem_countries_unmatched": int(cw[cw.in_vdem & ~cw.in_maddison].shape[0]),138 "y4_onsets": int((p.y4_civilwar_onset == 1).sum()), "y4_coded": int(p.y4_civilwar_onset.notna().sum()),139 "y5_coups": int((p.y5_coup == 1).sum()), "y5_coded": int(p.y5_coup.notna().sum()),140 "y6_irregular_ends": int((p.y6_irregular_regime_end == 1).sum()), "y6_coded": int(p.y6_irregular_regime_end.notna().sum())}141142143def main():144 s, cov = seshat_panels(); d = dplace_panel(); m = modern_panel()145 summary = {"built_at_utc": datetime.now(timezone.utc).isoformat(timespec="seconds"), "seshat": s, "dplace": d, "modern": m}146 (DIAG / "coverage_summary.json").write_text(json.dumps(summary, indent=2))147 print(json.dumps(summary, indent=2))148149150if __name__ == "__main__":151 main()152