import tempfile
import time
import warnings
from pathlib import Path

import numpy as np
import pandas as pd

import wingallom as w

warnings.filterwarnings("ignore", module="rdata")
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
import logging

started = time.time()
for name, digest in w.check_staged().items():
    print(f"{name}: SHA-256 {digest} (matches Dryad's file list, version {w.DRYAD_VERSION} of {w.DRYAD_ID})")
t0 = time.time()
print("CmdStan", w.CMDSTAN_VERSION, "at", w.install_cmdstan(Path(tempfile.gettempdir()) / "cmdstan"),
      f"({time.time() - t0:.0f} s to install or find)")
logging.getLogger("cmdstanpy").setLevel(logging.ERROR)   # sampler messages; divergences are counted below
code_zip = w.code_zip()
original = code_zip.read("Code/lm_mod.stan").decode()
model_code = w.port_stan(original)
workdir = Path(tempfile.mkdtemp())
t0 = time.time()
model = w.compile_stan(model_code, workdir)
cpu = next((line.split(":", 1)[1].strip() for line in open("/proc/cpuinfo") if line.startswith("model name")), "unknown")
print("CPU:", cpu, "(the model is compiled with -march=native, as Script01 does, so the draws depend on the CPU)")
print(f"lm_mod.stan compiled in {time.time() - t0:.0f} s")

READING = "c"

Data.zip: SHA-256 2185ac50e9a59c3a1dff59473216de2ab001d0d0044b1fa15442ae18a520d554 (matches Dryad's file list, version 333377 of bvq83bkk9)
README.md: SHA-256 6c19de67b1c6c52802be42a6593c175a02bee0dab0ca2f6f958f08a0a3fa0482 (matches Dryad's file list, version 333377 of bvq83bkk9)
CmdStan 2.40.0 at /tmp/cmdstan/cmdstan-2.40.0 (0 s to install or find)

/home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

03:39:44 - cmdstanpy - INFO - compiling stan file /tmp/tmprecf6m42/lm_mod.stan to exe file /tmp/tmprecf6m42/lm_mod

03:40:08 - cmdstanpy - INFO - compiled model executable: /tmp/tmprecf6m42/lm_mod

CPU: AMD EPYC 7763 64-Core Processor (the model is compiled with -march=native, as Script01 does, so the draws depend on the CPU)
lm_mod.stan compiled in 24 s

morph = w.morph_prepared(READING)
removed = w.outlier_mask(w.read_csv("morphdata.csv"), READING)
p = w.prepare(morph, "autumn AHY", "main")
data = w.stan_data(p, w.MAIN_PREDICTORS["autumn AHY"])
md_ = p["md"]
print(f"reading ({READING}): {int(removed.sum())} specimens removed by the filter; {len(morph)} kept")
print(f"autumn AHY: {data['N']} birds in {data['J']} species-years of {data['K']} species (paper: 3250 and 171)")
print(f"species-years dropped for a missing climate value: rows {p['missing_rows']} of {p['n_merged']}")
print("\nBirds per species and sex:")
print(pd.crosstab(md_.Taxon, md_.Sex).to_string())
print(f"\nYears {p['G'].year.min()} to {p['G'].year.max()}; year column {data['G'][:, -1].min():.2f} to {data['G'][:, -1].max():.2f} decades")
print(f"scaled wing {data['WL'].min():.2f} to {data['WL'].max():.2f}, scaled tarsus {data['X'][:, 1].min():.2f} to {data['X'][:, 1].max():.2f}")
print("climate columns of G (scaled within species), min and max:")
print(pd.DataFrame(data["G"][:, 1:-1], columns=w.MAIN_PREDICTORS["autumn AHY"][:-1]).agg(["min", "max"]).round(2).to_string())

reading (c): 17 specimens removed by the filter; 27469 kept
autumn AHY: 3249 birds in 171 species-years of 5 species (paper: 3250 and 171)
species-years dropped for a missing climate value: rows [np.int64(1), np.int64(36), np.int64(106), np.int64(111), np.int64(112), np.int64(142)] of 177

Birds per species and sex:
Sex                       f    m
Taxon                           
Junco_hyemalis          250  279
Melospiza_georgiana     353  440
Melospiza_melodia        78  112
Seiurus_aurocapilla     260  314
Zonotrichia_albicollis  627  536

Years 1982 to 2016; year column -1.60 to 1.80 decades
scaled wing -4.02 to 3.04, scaled tarsus -4.28 to 3.79
climate columns of G (scaled within species), min and max:
     PrevPrevPrecipWinMo  PrevTempBreedBase  PrevPrecipBreedBase  PrevMinTempFallMig  MinTempFallMig
min                -2.74              -1.92                -2.24               -2.95           -2.95
max                 2.65               2.00                 1.77                2.23            2.47

t0 = time.time()
fit = w.fit_stan(model, data, w.STAN_SEEDS["autumn AHY"], workdir)
fit_minutes = (time.time() - t0) / 60
print(f"fitted in {fit_minutes:.1f} minutes")
print(f"divergent transitions after warm-up: {fit['divergent']} of 8000 draws (chains 1 to 4: {fit['divergent_per_chain']}); "
      f"largest R-hat over mu_i, mu_s, mu_stl and the standard deviations: {fit['max_rhat']:.3f}; smallest bulk ESS: {fit['min_ess_bulk']:.0f}")
names = ["Intercept", *w.MAIN_PREDICTORS["autumn AHY"]]
diag = fit["mu_summary"].copy()
diag.index = [f"{'mu_i' if i.startswith('mu_i') else 'mu_s'} {names[int(i.split('[')[1][:-1]) - 1]}" for i in diag.index]
print("\nCmdStan's diagnostics of mu_i (intercept) and mu_s (slope):")
print(diag.round(4).to_string())
table3 = w.effect_table(fit, w.MAIN_PREDICTORS["autumn AHY"])
print("\nTable 3 as fitted here (mean, s.d. and pd of mu_i, effects on the allometric intercept, and mu_s, on the slope):")
print(table3[["int_mean", "int_sd", "int_pd", "slope_mean", "slope_sd", "slope_pd"]].round(3).to_string())
print("\nMonte Carlo error of each pd (from the effective sample size of the sign indicator) and the pd of each chain:")
mc = pd.DataFrame({"int_pd": table3.int_pd, "int_pd_mcse": table3.int_pd_mcse, "int_pd_ess": table3.int_pd_ess,
                   "int_pd_chains": [[round(x, 3) for x in v] for v in table3.int_pd_chains],
                   "slope_pd": table3.slope_pd, "slope_pd_mcse": table3.slope_pd_mcse, "slope_pd_ess": table3.slope_pd_ess,
                   "slope_pd_chains": [[round(x, 3) for x in v] for v in table3.slope_pd_chains]})
print(mc.round(4).to_string())
stuck = w.stuck_chains(fit)
if stuck:
    moving = w.effect_table(w.drop_chains(fit, stuck), w.MAIN_PREDICTORS["autumn AHY"])
    print(f"\nChain(s) {', '.join(str(c + 1) for c in stuck)} stuck (every kept draw divergent). The pd over the "
          f"{4 - len(stuck)} chains that move, pooled, with its Monte Carlo error:")
    print(moving[["int_mean", "int_pd", "int_pd_mcse", "slope_mean", "slope_pd", "slope_pd_mcse"]].round(4).to_string())
else:
    moving = None
    print("\nNo chain is stuck (none has every kept draw divergent).")


def show(term, paper):
    """The claim's numbers for one term: mean, s.d. and pd with its Monte Carlo error and chains,
    against the paper's mean, tolerance and pd."""
    row = table3.loc[term]
    print(f"{term}, reading ({READING}):")
    for part in ("int", "slope"):
        m, sd, pd_, se = row[part + "_mean"], row[part + "_sd"], row[part + "_pd"], row[part + "_pd_mcse"]
        pm, ptol, ppd = paper[part]
        print(f"  {'intercept' if part == 'int' else 'slope':9s} mean {m:+.3f} (s.d. {sd:.3f}), pd {pd_:.2%};"
              f" paper {pm:+.2f} (tolerance {ptol}), pd {ppd}%; within tolerance: {abs(m - pm) <= ptol + 1e-9};"
              f" sign as paper: {np.sign(m) == np.sign(pm) or pm == 0}; pd >= 90%: {pd_ >= 0.9}")
        print(f"            Monte Carlo error of pd {se:.2%} (ESS {row[part + '_pd_ess']:.0f}), so pd - 90% is "
              f"{(pd_ - 0.9) / se:+.1f} errors; pd of chains 1 to 4: {', '.join(f'{x:.2%}' for x in row[part + '_pd_chains'])}")
        if moving is not None:
            mrow = moving.loc[term]
            mpd, mse = mrow[part + "_pd"], mrow[part + "_pd_mcse"]
            print(f"            chains that move, pooled: mean {mrow[part + '_mean']:+.3f}, pd {mpd:.2%}, Monte Carlo error {mse:.2%},"
                  f" so pd - 90% is {(mpd - 0.9) / mse:+.1f} errors; pd >= 90%: {mpd >= 0.9}")

03:40:08 - cmdstanpy - INFO - CmdStan start processing

03:40:08 - cmdstanpy - INFO - Chain [1] start processing

03:40:08 - cmdstanpy - INFO - Chain [2] start processing

03:40:08 - cmdstanpy - INFO - Chain [3] start processing

03:40:08 - cmdstanpy - INFO - Chain [4] start processing

03:49:10 - cmdstanpy - INFO - Chain [1] done processing

03:50:19 - cmdstanpy - INFO - Chain [4] done processing

03:51:18 - cmdstanpy - INFO - Chain [2] done processing

03:52:22 - cmdstanpy - INFO - Chain [3] done processing

03:52:22 - cmdstanpy - WARNING - Non-fatal error during sampling:
Exception: normal_lpdf: Scale parameter[1] is 0, but must be positive! (in 'lm_mod.stan', line 43, column 6 to column 35)
Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in 'lm_mod.stan', line 38, column 4 to column 29)
	Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in 'lm_mod.stan', line 38, column 4 to column 29)
Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in 'lm_mod.stan', line 38, column 4 to column 29)
	Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in 'lm_mod.stan', line 38, column 4 to column 29)
	Exception: normal_lpdf: Scale parameter[2] is 0, but must be positive! (in 'lm_mod.stan', line 43, column 6 to column 35)
Consider re-running with show_console=True if the above output is unclear!

paper = {"int": (-0.09, 0.03, 99.5), "slope": (-0.01, 0.03, 75.1)}
show("PrevTempBreedBase", paper)

PrevTempBreedBase, reading (c):
  intercept mean -0.091 (s.d. 0.029), pd 99.39%; paper -0.09 (tolerance 0.03), pd 99.5%; within tolerance: True; sign as paper: True; pd >= 90%: True
            Monte Carlo error of pd 0.17% (ESS 2041), so pd - 90% is +54.4 errors; pd of chains 1 to 4: 99.85%, 99.45%, 99.45%, 98.80%
  slope     mean -0.026 (s.d. 0.045), pd 76.51%; paper -0.01 (tolerance 0.03), pd 75.1%; within tolerance: True; sign as paper: True; pd >= 90%: False
            Monte Carlo error of pd 1.16% (ESS 1338), so pd - 90% is -11.6 errors; pd of chains 1 to 4: 82.80%, 75.70%, 74.55%, 73.00%

paper = {"int": (-0.02, 0.07, 78.0), "slope": (0.01, 0.05, 70.9)}
show("PrevMinTempFallMig", paper)

PrevMinTempFallMig, reading (c):
  intercept mean -0.019 (s.d. 0.045), pd 74.17%; paper -0.02 (tolerance 0.07), pd 78.0%; within tolerance: True; sign as paper: True; pd >= 90%: False
            Monte Carlo error of pd 1.14% (ESS 1473), so pd - 90% is -13.9 errors; pd of chains 1 to 4: 74.60%, 74.00%, 73.80%, 74.30%
  slope     mean +0.008 (s.d. 0.031), pd 63.35%; paper +0.01 (tolerance 0.05), pd 70.9%; within tolerance: True; sign as paper: True; pd >= 90%: False
            Monte Carlo error of pd 1.96% (ESS 603), so pd - 90% is -13.6 errors; pd of chains 1 to 4: 55.30%, 67.45%, 62.75%, 67.90%

paper = {"int": (0.03, 0.04, 88.6), "slope": (-0.03, 0.05, 90.3)}
show("MinTempFallMig", paper)

MinTempFallMig, reading (c):
  intercept mean +0.027 (s.d. 0.029), pd 88.92%; paper +0.03 (tolerance 0.04), pd 88.6%; within tolerance: True; sign as paper: True; pd >= 90%: False
            Monte Carlo error of pd 0.65% (ESS 2333), so pd - 90% is -1.7 errors; pd of chains 1 to 4: 90.60%, 88.10%, 88.75%, 88.25%
  slope     mean -0.032 (s.d. 0.035), pd 91.24%; paper -0.03 (tolerance 0.05), pd 90.3%; within tolerance: True; sign as paper: True; pd >= 90%: True
            Monte Carlo error of pd 0.54% (ESS 2712), so pd - 90% is +2.3 errors; pd of chains 1 to 4: 92.85%, 90.70%, 90.95%, 90.45%

paper = {"int": (0.03, 0.03, 90.0), "slope": (0.02, 0.03, 88.7)}
show("PrevPrevPrecipWinMo", paper)

PrevPrevPrecipWinMo, reading (c):
  intercept mean +0.029 (s.d. 0.027), pd 89.15%; paper +0.03 (tolerance 0.03), pd 90.0%; within tolerance: True; sign as paper: True; pd >= 90%: False
            Monte Carlo error of pd 0.66% (ESS 2198), so pd - 90% is -1.3 errors; pd of chains 1 to 4: 92.15%, 89.05%, 88.90%, 86.50%
  slope     mean +0.020 (s.d. 0.035), pd 81.46%; paper +0.02 (tolerance 0.03), pd 88.7%; within tolerance: True; sign as paper: True; pd >= 90%: False
            Monte Carlo error of pd 10.47% (ESS 14), so pd - 90% is -0.8 errors; pd of chains 1 to 4: 59.80%, 87.15%, 90.10%, 88.80%

print(table3.loc[["y_id", "PrevPrecipBreedBase"], ["int_mean", "int_sd", "int_pd", "int_pd_mcse", "slope_mean", "slope_sd", "slope_pd", "slope_pd_mcse"]].round(3).to_string())
print("paper, Table 3: year intercept 0.19 (pd 99.8%), slope 0.02 (pd 75.9%); PrevPrecipBreedBase intercept -0.05 (pd 93.4%), slope -0.01 (pd 65.3%)")
climate = [t for t in w.MAIN_PREDICTORS["autumn AHY"] if t != "y_id"]
year = abs(table3.loc["y_id", "int_mean"])
print(f"\nclimate-explains-little, reading (a) of its criterion, autumn AHY part: |year effect on the intercept| {year:.3f}; "
      f"|climate effects on the intercept|: " + ", ".join(f"{t} {abs(table3.loc[t, 'int_mean']):.3f}" for t in climate))
print(f"every climate effect smaller than the year effect: {all(abs(table3.loc[t, 'int_mean']) < year for t in climate)}")
print(f"\nNotebook time: {(time.time() - started) / 60:.1f} minutes")

                     int_mean  int_sd  int_pd  int_pd_mcse  slope_mean  slope_sd  slope_pd  slope_pd_mcse
term                                                                                                     
y_id                    0.197   0.049   0.995        0.002       0.014     0.028     0.731          0.016
PrevPrecipBreedBase    -0.056   0.039   0.938        0.005      -0.007     0.043     0.562          0.067
paper, Table 3: year intercept 0.19 (pd 99.8%), slope 0.02 (pd 75.9%); PrevPrecipBreedBase intercept -0.05 (pd 93.4%), slope -0.01 (pd 65.3%)

climate-explains-little, reading (a) of its criterion, autumn AHY part: |year effect on the intercept| 0.197; |climate effects on the intercept|: PrevPrevPrecipWinMo 0.029, PrevTempBreedBase 0.091, PrevPrecipBreedBase 0.056, PrevMinTempFallMig 0.019, MinTempFallMig 0.027
every climate effect smaller than the year effect: True

Notebook time: 12.7 minutes

Dias et al. (2025): the autumn AHY main model (Table 3), filtering reading (c)¶

Data check¶

Fit¶

table3-autumn-ahy-prior-breeding-temperature¶

table3-autumn-ahy-prior-autumn-temperature¶

table3-autumn-ahy-recent-autumn-temperature¶

table3-autumn-ahy-prior-wintering-precipitation¶

Autumn AHY parts of claims that also need the HY and spring AHY models¶