#!/usr/bin/env python3
"""Do Estonia's cereal yields follow the thermal growing season, year to year?

Our check for https://factchecked.science/estonia-growing-season/ (Fact-Checked Science, October 26, 2026).
Inputs, next to this script: gsl-rerun-numbers.json (written by gsl_rerun.py: the 11-station network series of
season length, start and end) and estonia-crops.csv (Statistics Estonia, PM0281, CC BY-SA 4.0: the national
cereal yield). Yields are detrended over 1992-2025, the years after the post-Soviet collapse (2026 is left out,
being a preliminary estimate), with a straight line and, separately, by first differences; the script prints
Pearson and Spearman correlations with the season's length, start and end. A screening, not an attribution:
it cannot rule out a slow climate contribution to the long-term rise in yields.

Usage: python3 yield_vs_gsl.py      Standard library only (Python 3.8 or later).
"""
import csv
import json
import math
import os
import statistics
import sys

HERE = os.path.dirname(os.path.abspath(__file__))


def pearson(x, y):
    mx, my = statistics.fmean(x), statistics.fmean(y)
    sxy = sum((a - mx) * (b - my) for a, b in zip(x, y))
    return sxy / math.sqrt(sum((a - mx) ** 2 for a in x) * sum((b - my) ** 2 for b in y))


def ranks(v):
    order = sorted(range(len(v)), key=lambda i: v[i])
    r = [0.0] * len(v)
    i = 0
    while i < len(v):
        j = i
        while j + 1 < len(v) and v[order[j + 1]] == v[order[i]]:
            j += 1
        for k in range(i, j + 1):
            r[order[k]] = (i + j) / 2 + 1
        i = j + 1
    return r


def p_approx(r, n):
    """Two-sided p from Fisher's z, an approximation adequate for a screening statement."""
    return math.erfc(abs(0.5 * math.log((1 + r) / (1 - r)) * math.sqrt(n - 3)) / math.sqrt(2))


def main():
    if hasattr(sys.stdout, "reconfigure"):
        sys.stdout.reconfigure(encoding="utf-8")
    net = json.load(open(os.path.join(HERE, "gsl-rerun-numbers.json"), encoding="utf-8"))["network"]
    series = {k: {int(y): v for y, v in net[k + "_series"].items()} for k in ("GSL", "SOS", "EOS")}
    yi = {}
    with open(os.path.join(HERE, "estonia-crops.csv"), encoding="utf-8") as f:
        for row in csv.DictReader(f):
            if row["crop_code"] == "CRC0000" and row["indicator_code"] == "YI":
                yi[int(row["year"])] = float(row["value"])
    years = [y for y in range(1992, 2026) if y in yi and y in series["GSL"]]
    ys = [yi[y] for y in years]
    mx, my = statistics.fmean(years), statistics.fmean(ys)
    b = sum((x - mx) * (v - my) for x, v in zip(years, ys)) / sum((x - mx) ** 2 for x in years)
    resid = [v - (my + b * (x - mx)) for x, v in zip(years, ys)]
    print(f"Cereal yield {years[0]}-{years[-1]}: trend {b:+.1f} kg/ha a year; {len(years)} years with a season value")
    names = {"GSL": "season length", "SOS": "season start", "EOS": "season end"}
    for k in ("GSL", "SOS", "EOS"):
        x = [series[k][y] for y in years]
        r = pearson(x, resid)
        rs = pearson(ranks(x), ranks(resid))
        pairs = [(series[k][years[i]] - series[k][years[i - 1]], yi[years[i]] - yi[years[i - 1]])
                 for i in range(1, len(years)) if years[i] - years[i - 1] == 1]
        r2 = pearson([p[0] for p in pairs], [p[1] for p in pairs])
        print(f"  {names[k]:13s}: detrended r = {r:+.2f} (Spearman {rs:+.2f}, p about {p_approx(r, len(years)):.2f});"
              f" first differences r = {r2:+.2f} ({len(pairs)} pairs, p about {p_approx(r2, len(pairs)):.2f})")


if __name__ == "__main__":
    main()