# -*- coding: utf-8 -*-
# ------------------------------------------------------------------------------
# Hangi iklim kipi ulusal deniz avında iz bırakır? On bir iklim kipi ve 146 ülke için alan anlamlılığı, El Niño'nun çeşidi ve naif taramanın yanılgısı, 1952–2024
# Which climate mode leaves a footprint in national marine catches? Field significance for eleven climate modes and 146 countries, the flavour of El Niño and the fallacy of naive screening, 1952–2024
#
# GOMD Yayınlar · Makale 23 — https://gomd.pirireis.edu.tr/publications.php?id=rp23
#
# Bu dosya makalenin bütün modüllerini içerir. Replikasyon paketi klasöründe çalıştırıldığında modülleri `e/` klasörüne yazar ve sırayla, her birini ayrı süreçte çalıştırır;
# veriyi `data/` ve `e/out/` klasörlerinden okur, sonuçları `e/out/` ve şekilleri `e/fig/` klasörüne yazar.
# This file contains all modules of the paper. Run in the replication-package folder, it writes the modules to `e/` and runs them in order, each in its own process;
# it reads data from `data/` and `e/out/`, and writes results to `e/out/` and figures to `e/fig/`.
#
# Kurulum / Setup:  pip install pandas numpy matplotlib scipy statsmodels pyarrow
# Çalıştırma / Run:  python rp23_climate_modes_national_catch_field_significance.py
#
# Catch data: FAO Global Capture Production (FishStat, release 2026.1.0), © FAO, https://www.fao.org/fishery/en/collection/capture . Climate indices: NOAA Physical Sciences Laboratory (https://psl.noaa.gov/data/climateindices/) and NOAA Climate Prediction Center (Oceanic Niño Index).
# GOMD, Piri Reis Üniversitesi ile T.C. Çevre, Şehircilik ve İklim Değişikliği Bakanlığı Çevre Yönetimi Genel Müdürlüğü tarafından geliştirilen açık veri platformudur.
# GOMD is an open data platform developed by Piri Reis University together with the Directorate General of Environmental Management of the Ministry of Environment, Urbanisation and Climate Change of Türkiye.
# ------------------------------------------------------------------------------
import os, subprocess, sys

MODULES = {}
ORDER = ['e/e_analysis.py', 'e/e_extra.py', 'e/e_extra_b.py', 'e/e_modes.py', 'e/e_fig.py']
SKIP_IF_EXISTS = {}   # hazırlık modülü: çıktısı pakette varsa atlanır / preparation module: skipped if its output is in the package

# ==============================================================================
# e/e_analysis.py
# ==============================================================================
MODULES['e/e_analysis.py'] = r'''
# ENSO'nun ulusal deniz avindaki izi: FAO av (1950-2024) + ONI
import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore')
import statsmodels.api as sm
from statsmodels.nonparametric.smoothers_lowess import lowess
from scipy import stats
m=pd.read_parquet('f/out/marine.parquet')
oni=pd.read_csv('data/oni.data',sep=r'\s+',skiprows=1,header=None,nrows=77).set_index(0); oni=oni.where(oni>-90)
# ENSO yili: Tem(t-1)-Haz(t) ortalamasi (kis tepe merkezli); ONI sutunlari 1=DJF..12=NDJ (uc aylik merkezler Oca..Ara)
mo=oni.stack(); mo.index=pd.PeriodIndex([f'{int(y)}-{int(k):02d}' for y,k in mo.index],freq='M'); mo=mo.sort_index()
E={y:float(mo[f'{y-1}-07':f'{y}-06'].mean()) for y in range(1951,2025)}; Ecal={y:float(mo[f'{y}-01':f'{y}-12'].mean()) for y in range(1950,2025)}
ES=pd.Series(E); R={'E':E,'years':[1952,2024],'n_el':int((ES>=0.5).sum()),'n_la':int((ES<=-0.5).sum()),'sd_E':float(ES.std()),'top_el':[[int(y),float(v)] for y,v in ES.sort_values(ascending=False).head(6).items()]}
YRS=list(range(1952,2025))
def detr(s):
    """log av - LOWESS egilimi (frac 0.3)"""
    s=s.reindex(range(1950,2025)); 
    if s.isna().any() or (s<=0).any(): return None
    l=np.log(s.values); tr=lowess(l,np.arange(len(l)),frac=0.3,return_sorted=False); return pd.Series(l-tr,index=s.index)
def fit(y,lag=True,nw=3):
    y=y.loc[1952:2024]; X=pd.DataFrame({'e0':[E[t] for t in y.index],'e1':[E[t-1] for t in y.index]},index=y.index)
    mm=sm.OLS(y.values,sm.add_constant(X.values)).fit(cov_type='HAC',cov_kwds={'maxlags':nw})
    r=dict(b0=float(mm.params[1]),se0=float(mm.bse[1]),p0=float(mm.pvalues[1]),b1=float(mm.params[2]),se1=float(mm.bse[2]),p1=float(mm.pvalues[2]),r2=float(mm.rsquared),sd=float(y.std()))
    w=mm.wald_test(np.array([[0,1,0],[0,0,1]]),scalar=True); r['pj']=float(w.pvalue); r['bsum']=float(mm.params[1]+mm.params[2])
    c=mm.cov_params(); r['se_sum']=float(np.sqrt(c[1,1]+c[2,2]+2*c[1,2])); return r
def asym(y):
    y=y.loc[1952:2024]; e=np.array([E[t] for t in y.index]); X=np.column_stack([np.maximum(e,0),np.minimum(e,0)]); mm=sm.OLS(y.values,sm.add_constant(X)).fit(cov_type='HAC',cov_kwds={'maxlags':3})
    w=mm.wald_test(np.array([[0,1,-1]]),scalar=True); return dict(bp=float(mm.params[1]),bn=float(mm.params[2]),p_diff=float(w.pvalue),pp=float(mm.pvalues[1]),pn=float(mm.pvalues[2]))
# ---- dunya toplami
tot=m.groupby('yr').val.sum(); anc=m[m.sp=='VET'].groupby('yr').val.sum().reindex(tot.index).fillna(0)
R['world']={k:dict(**fit(detr(s)),**asym(detr(s))) for k,s in(('total',tot),('ex_anchoveta',tot-anc),('anchoveta',anc.replace(0,np.nan).fillna(anc[anc>0].min())))}
R['world']['mean_2000_24']=float(tot.loc[2000:2024].mean()); R['world']['anch_mean_2000_24']=float(anc.loc[2000:2024].mean())
# ---- FAO alanlari
AR={}
for (a,nm),g in m.groupby(['area','area_name']):
    y=detr(g.groupby('yr').val.sum())
    if y is None: continue
    AR[nm]=dict(**fit(y),mean=float(g[g.yr>=2000].groupby('yr').val.sum().mean()))
R['area']=AR
# ---- ISSCAAP gruplari (dunya)
GR={}
for k,g in m.groupby('isscaap'):
    s=g.groupby('yr').val.sum(); y=detr(s)
    if y is None or s.loc[2000:2024].mean()<50000: continue
    GR[k]=dict(**fit(y),mean=float(s.loc[2000:2024].mean()))
R['group']=GR
# ---- ulkeler: tam seri 1950-2024, ortalama >= 1000 t
CO={}
for un,g in m.groupby('un'):
    s=g.groupby('yr').val.sum(); y=detr(s)
    if y is None or s.mean()<1000: continue
    i=g.iloc[0]; r=fit(y); r.update(asym(y)); r.update(country=i.country,iso=i.iso,cont=i.cont,eco=i.eco,mean=float(s.loc[2000:2024].mean()),smallpel=float(g[(g.isscaap=='Herrings, sardines, anchovies')&(g.yr>=2000)].val.sum()/g[g.yr>=2000].val.sum()),tuna=float(g[(g.isscaap=='Tunas, bonitos, billfishes')&(g.yr>=2000)].val.sum()/g[g.yr>=2000].val.sum()))
    # Pasifik payi
    r['pacific']=float(g[g.area.isin(['61','67','71','77','81','87'])&(g.yr>=2000)].val.sum()/g[g.yr>=2000].val.sum())
    # donem ayrimi
    for nm,(a,b) in(('early',(1952,1987)),('late',(1988,2024))):
        yy=y.loc[a:b]; X=np.column_stack([[E[t] for t in yy.index],[E[t-1] for t in yy.index]]); mm=sm.OLS(yy.values,sm.add_constant(X)).fit(cov_type='HAC',cov_kwds={'maxlags':3}); r['b0_'+nm]=float(mm.params[1]); r['p0_'+nm]=float(mm.pvalues[1]); r['r2_'+nm]=float(mm.rsquared)
    CO[un]=r
C=pd.DataFrame(CO).T; 
for c in C.columns:
    try: C[c]=pd.to_numeric(C[c])
    except Exception: pass
# FDR (Benjamini-Hochberg) ortak test p'leri
def bh(p,q=0.05):
    p=np.asarray(p); o=np.argsort(p); n=len(p); thr=q*np.arange(1,n+1)/n; ok=p[o]<=thr; k=np.max(np.where(ok)[0])+1 if ok.any() else 0; sig=np.zeros(n,bool); sig[o[:k]]=True; return sig
C['sig_fdr']=bh(C.pj.values,0.05); C['sig_fdr10']=bh(C.pj.values,0.10); C['sig05']=C.pj<0.05; C['sig0']=C.p0<0.05
C.to_csv('e/out/country_enso.csv')
R['country']=dict(n=int(len(C)),sig05=int(C.sig05.sum()),sig_fdr=int(C.sig_fdr.sum()),sig_fdr10=int(C.sig_fdr10.sum()),sig0=int(C.sig0.sum()),neg0=int((C.b0<0).sum()),neg_sig0=int(((C.b0<0)&C.sig0).sum()),pos_sig0=int(((C.b0>0)&C.sig0).sum()),
   r2_med=float(C.r2.median()),r2_p90=float(C.r2.quantile(.9)),r2_w=float(np.average(C.r2,weights=C['mean'])),absb_med=float(C.b0.abs().median()),sd_med=float(C.sd.median()),share_catch_sig=float(C[C.sig05]['mean'].sum()/C['mean'].sum()),share_catch_fdr=float(C[C.sig_fdr]['mean'].sum()/C['mean'].sum()),
   asym_sig=int((C.p_diff<0.05).sum()))
f=lambda r:[r.country,float(r['mean']),float(r.b0),float(r.se0),float(r.p0),float(r.b1),float(r.p1),float(r.r2),float(r.pj),bool(r.sig_fdr),float(r.smallpel),float(r.pacific)]
big=C[C['mean']>=50000]
R['country']['n_big']=int(len(big)); R['country']['most_neg']=[f(r) for _,r in big.sort_values('b0').head(12).iterrows()]; R['country']['most_pos']=[f(r) for _,r in big.sort_values('b0',ascending=False).head(12).iterrows()]
R['country']['top_r2']=[f(r) for _,r in big.sort_values('r2',ascending=False).head(12).iterrows()]; R['country']['largest']=[f(r) for _,r in C.sort_values('mean',ascending=False).head(20).iterrows()]
R['country']['fdr_list']=[f(r) for _,r in C[C.sig_fdr].sort_values('pj').iterrows()]
# ---- permutasyon: ONI serisini dairesel kaydirarak tesadufi anlamlilik sayisi
Y={un:detr(g.groupby('yr').val.sum()) for un,g in m.groupby('un') if un in CO}
Ymat=np.column_stack([Y[u].loc[1952:2024].values for u in C.index]); ev=np.array([E[t] for t in range(1951,2025)])
def count_sig(e):
    X=sm.add_constant(np.column_stack([e[1:],e[:-1]])); n=0; r2=[]
    XtXi=np.linalg.inv(X.T@X); B=XtXi@X.T@Ymat; res=Ymat-X@B; 
    # hizli: OLS F-testi (HAC yerine) ile ortak anlamlilik
    ssr=(res**2).sum(0); sst=((Ymat-Ymat.mean(0))**2).sum(0); r2=1-ssr/sst; Fst=(r2/2)/((1-r2)/(len(e)-1-3)); p=stats.f.sf(Fst,2,len(e)-1-3); return int((p<0.05).sum()),float(np.median(r2)),float(np.average(r2,weights=C['mean'].values))
obs=count_sig(ev); rng=np.random.default_rng(7); perm=[]
for k in range(1000):
    # faz rastgelelestirme (vekil seri): ayni guc spektrumu
    f_=np.fft.rfft(ev-ev.mean()); ph=np.exp(1j*rng.uniform(0,2*np.pi,len(f_))); ph[0]=1; s_=np.fft.irfft(np.abs(f_)*ph,n=len(ev)); perm.append(count_sig(s_))
perm=np.array(perm); R['perm']=dict(obs_sig=obs[0],obs_r2_med=obs[1],obs_r2_w=obs[2],null_sig_mean=float(perm[:,0].mean()),null_sig_p95=float(np.percentile(perm[:,0],95)),p_sig=float((perm[:,0]>=obs[0]).mean()),null_r2w_mean=float(perm[:,2].mean()),null_r2w_p95=float(np.percentile(perm[:,2],95)),p_r2w=float((perm[:,2]>=obs[2]).mean()),null_r2med_mean=float(perm[:,1].mean()),p_r2med=float((perm[:,1]>=obs[1]).mean()))
# ---- belirleyiciler: |b0| ve r2 ~ Pasifik payi, kucuk pelajik, ton baligi, log buyukluk
import statsmodels.formula.api as smf
C['absb']=C.b0.abs(); C['lm']=np.log(C['mean'].clip(lower=1)); C['lr2']=np.log(C.r2.clip(lower=1e-4))
for dep in('absb','r2'):
    mm=smf.ols(dep+' ~ pacific + smallpel + tuna + lm',C).fit(cov_type='HC1'); R['det_'+dep]={k:dict(b=float(mm.params[k]),se=float(mm.bse[k]),p=float(mm.pvalues[k])) for k in('pacific','smallpel','tuna','lm')}; R['det_'+dep]['r2']=float(mm.rsquared)
R['by_cont']={str(k):dict(n=int(len(g)),sig=int(g.sig05.sum()),r2_med=float(g.r2.median()),b0_med=float(g.b0.median()),b0_w=float(np.average(g.b0,weights=g['mean']))) for k,g in C.groupby('cont')}
pac=C.pacific>=0.5; R['pacific_split']={k:dict(n=int(len(g)),sig=float(g.sig05.mean()),r2_med=float(g.r2.median()),absb=float(g.absb.median())) for k,g in(('pacific',C[pac]),('other',C[~pac]))}
# ---- donem karsilastirmasi
R['period']=dict(corr=float(np.corrcoef(C.b0_early,C.b0_late)[0,1]),spearman=float(stats.spearmanr(C.b0_early,C.b0_late)[0]),r2_early=float(C.r2_early.median()),r2_late=float(C.r2_late.median()),sig_early=int((C.p0_early<0.05).sum()),sig_late=int((C.p0_late<0.05).sum()),same_sign=float((np.sign(C.b0_early)==np.sign(C.b0_late)).mean()),
   wilc=float(stats.wilcoxon(C.r2_late-C.r2_early).pvalue))
# ---- guclu El Nino (E=+1.5) icin beklenen degisim (ton): kazananlar/kaybedenler
C['dton']=C['mean']*(np.exp(1.5*C.b0)-1); R['tonnage']=dict(E=1.5,loss=float(C[C.dton<0].dton.sum()),gain=float(C[C.dton>0].dton.sum()),net=float(C.dton.sum()),loss_sig=float(C[(C.dton<0)&C.sig0].dton.sum()),gain_sig=float(C[(C.dton>0)&C.sig0].dton.sum()),top_loss=[[r.country,float(r.dton),float(r.b0)] for _,r in C.sort_values('dton').head(8).iterrows()],top_gain=[[r.country,float(r.dton),float(r.b0)] for _,r in C.sort_values('dton',ascending=False).head(8).iterrows()],world_mean=float(C['mean'].sum()),peru_share_loss=float(C[C.country=='Peru'].dton.sum()/C[C.dton<0].dton.sum()))
# ---- stok duzeyi: ulke x takson (ortalama >= 50 bin ton, tam seri) 
ST=[]
for (un,sp),g in m.groupby(['un','sp']):
    s=g.groupby('yr').val.sum()
    if len(s)<75 or s.loc[2000:2024].mean()<50000: continue
    y=detr(s)
    if y is None: continue
    r=fit(y); i=g.iloc[0]; r.update(country=i.country,sp=i.sp_name,isscaap=i.isscaap,mean=float(s.loc[2000:2024].mean())); ST.append(r)
S=pd.DataFrame(ST); S['sig_fdr']=bh(S.pj.values,0.05); S.to_csv('e/out/stock_enso.csv',index=False)
R['stock']=dict(n=int(len(S)),sig05=int((S.pj<0.05).sum()),sig_fdr=int(S.sig_fdr.sum()),r2_med=float(S.r2.median()),top=[[r.country,r.sp,float(r['mean']),float(r.b0),float(r.b1),float(r.r2),float(r.pj)] for _,r in S.sort_values('r2',ascending=False).head(15).iterrows()],
   by_group={k:dict(n=int(len(g)),sig=float((g.pj<0.05).mean()),r2=float(g.r2.median())) for k,g in S.groupby('isscaap') if len(g)>=8})
def clean(o):
    if isinstance(o,dict): return {str(k):clean(v) for k,v in o.items()}
    if isinstance(o,(list,tuple)): return [clean(v) for v in o]
    if isinstance(o,(float,np.floating)): return None if not np.isfinite(o) else float(o)
    if isinstance(o,(np.integer,)): return int(o)
    if isinstance(o,(np.bool_,)): return bool(o)
    return o
json.dump(clean(R),open('e/out/results.json','w'),ensure_ascii=False,indent=1)
pr=lambda d:{k:(round(v,3) if isinstance(v,float) else v) for k,v in d.items()}
print('E',R['n_el'],R['n_la'],round(R['sd_E'],2),R['top_el'])
for k in('total','ex_anchoveta','anchoveta'): print('WORLD',k,pr(R['world'][k]))
for k,v in sorted(R['area'].items(),key=lambda x:x[1]['pj']): print('AREA',k[:28],pr({a:v[a] for a in('b0','p0','b1','p1','r2','pj')}),round(v['mean']/1e6,2))
for k,v in sorted(R['group'].items(),key=lambda x:x[1]['pj'])[:14]: print('GROUP',k[:30],pr({a:v[a] for a in('b0','p0','b1','p1','r2','pj')}),round(v['mean']/1e6,2))
print('COUNTRY',{k:v for k,v in R['country'].items() if not isinstance(v,list)})
g=lambda L:[[x[0][:22],round(x[1]/1e3),round(x[2],3),round(x[4],3),round(x[5],3),round(x[7],2),round(x[8],4),x[9]] for x in L]
print('NEG',g(R['country']['most_neg'])); print('POS',g(R['country']['most_pos'])); print('R2',g(R['country']['top_r2'])); print('FDR',g(R['country']['fdr_list']))
print('PERM',pr(R['perm'])); print('DET',R['det_absb'],R['det_r2']); print('CONT',R['by_cont']); print('PAC',R['pacific_split']); print('PERIOD',pr(R['period'])); print('TON',R['tonnage'])
print('STOCK',{k:v for k,v in R['stock'].items() if k not in('top','by_group')},[[x[0][:14],x[1][:26],round(x[2]/1e3),round(x[3],2),round(x[4],2),round(x[5],2)] for x in R['stock']['top']],R['stock']['by_group'])
'''

# ==============================================================================
# e/e_extra.py
# ==============================================================================
MODULES['e/e_extra.py'] = r'''
import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore')
import statsmodels.api as sm
from statsmodels.nonparametric.smoothers_lowess import lowess
from scipy import stats
m=pd.read_parquet('f/out/marine.parquet'); R=json.load(open('e/out/results.json')); C=pd.read_csv('e/out/country_enso.csv',dtype={'Unnamed: 0':str}).rename(columns={'Unnamed: 0':'un'}).set_index('un')
E={int(k):v for k,v in R['E'].items()}; X={}
oni=pd.read_csv('data/oni.data',sep=r'\s+',skiprows=1,header=None,nrows=77).set_index(0); oni=oni.where(oni>-90); mo=oni.stack(); mo.index=pd.PeriodIndex([f'{int(y)}-{int(k):02d}' for y,k in mo.index],freq='M'); mo=mo.sort_index()
Ecal={y:float(mo[f'{y}-01':f'{y}-12'].mean()) for y in range(1950,2025)}
def detr(s,frac=0.3):
    s=s.reindex(range(1950,2025))
    if s.isna().any() or (s<=0).any(): return None
    l=np.log(s.values); return pd.Series(l-lowess(l,np.arange(len(l)),frac=frac,return_sorted=False),index=s.index)
tot={un:g.groupby('yr').val.sum() for un,g in m.groupby('un') if un in C.index}
# 1. rastgele etkiler meta-analizi (DerSimonian-Laird) b0 icin
def meta(b,se):
    b=np.asarray(b,float); se=np.asarray(se,float); w=1/se**2; mu=(w*b).sum()/w.sum(); Q=(w*(b-mu)**2).sum(); k=len(b); c=w.sum()-(w**2).sum()/w.sum(); tau2=max(0,(Q-(k-1))/c); wr=1/(se**2+tau2); mur=(wr*b).sum()/wr.sum(); ser=np.sqrt(1/wr.sum())
    return dict(k=int(k),fixed=float(mu),random=float(mur),se=float(ser),p=float(2*stats.norm.sf(abs(mur/ser))),Q=float(Q),pQ=float(stats.chi2.sf(Q,k-1)),I2=float(max(0,(Q-(k-1))/Q)),tau=float(np.sqrt(tau2)))
X['meta_all']=meta(C.b0,C.se0); X['meta_b1']=meta(C.b1,C.se1)
ex=C[C.country!='Peru']; X['meta_ex_peru']=meta(ex.b0,ex.se0)
X['meta_cont']={k:meta(g.b0,g.se0) for k,g in C.groupby('cont')}
EP=['Peru','Chile','Ecuador','Colombia','Panama','Costa Rica','Mexico','El Salvador','Guatemala','Nicaragua']; g=C[C.country.isin(EP)]; X['meta_epac']=dict(meta(g.b0,g.se0),names=list(g.country))
WP=C[(C.pacific>=0.5)&~C.country.isin(EP)]; X['meta_wpac']=meta(WP.b0,WP.se0)
X['meta_big']=meta(C[C['mean']>=500000].b0,C[C['mean']>=500000].se0)
# 2. guc: en kucuk saptanabilir etki (%80 guc, %5 duzey) = 2.8*se
C['mde']=2.8*C.se0; X['power']=dict(mde_med=float(C.mde.median()),mde_p25=float(C.mde.quantile(.25)),mde_p75=float(C.mde.quantile(.75)),mde_big=float(C[C['mean']>=500000].mde.median()),share_mde_lt5=float((C.mde<0.05).mean()),share_mde_lt10=float((C.mde<0.10).mean()),se_med=float(C.se0.median()),
   ci_excl_10=float(((C.b0-1.96*C.se0>-0.10)&(C.b0+1.96*C.se0<0.10)).mean()),ci_excl_5=float(((C.b0-1.96*C.se0>-0.05)&(C.b0+1.96*C.se0<0.05)).mean()),catch_share_excl5=float(C[((C.b0-1.96*C.se0>-0.05)&(C.b0+1.96*C.se0<0.05))]['mean'].sum()/C['mean'].sum()),catch_share_excl10=float(C[((C.b0-1.96*C.se0>-0.10)&(C.b0+1.96*C.se0<0.10))]['mean'].sum()/C['mean'].sum()))
# 3. saglamlik: farkli belirlemeler icin anlamli ulke sayisi ve dunya katsayisi
def run(spec):
    out=[]; 
    for un,s in tot.items():
        if spec=='growth':
            y=np.log(s.reindex(range(1950,2025))).diff().loc[1952:2024]; e0=np.array([E[t]-E[t-1] for t in y.index]); e1=np.array([E[t-1]-E[t-2] if t-2 in E else np.nan for t in y.index]); ok=~np.isnan(e1); Xm=sm.add_constant(np.column_stack([e0[ok],e1[ok]])); yy=y.values[ok]
        else:
            fr={'lowess02':0.2,'lowess05':0.5}.get(spec,0.3); y=detr(s,fr).loc[1952:2024]; ee=Ecal if spec=='calendar' else E
            if spec=='lag2': Xm=sm.add_constant(np.column_stack([[E[t] for t in y.index],[E[t-1] for t in y.index],[E.get(t-2,np.nan) for t in y.index]])); ok=~np.isnan(Xm).any(axis=1); Xm=Xm[ok]; yy=y.values[ok]
            elif spec=='calendar': Xm=sm.add_constant(np.column_stack([[ee[t] for t in y.index],[ee[t-1] for t in y.index]])); yy=y.values
            else: Xm=sm.add_constant(np.column_stack([[E[t] for t in y.index],[E[t-1] for t in y.index]])); yy=y.values
        mm=sm.OLS(yy,Xm).fit(cov_type='HAC',cov_kwds={'maxlags':3}); k=Xm.shape[1]-1; w=mm.wald_test(np.eye(k+1)[1:],scalar=True); out.append((float(w.pvalue),float(mm.pvalues[1]),float(mm.params[1]),float(mm.rsquared)))
    o=np.array(out); return dict(n=int(len(o)),sig_joint=int((o[:,0]<0.05).sum()),sig_b0=int((o[:,1]<0.05).sum()),r2_med=float(np.median(o[:,3])),corr_b0=float(np.corrcoef(o[:,2],C.loc[list(tot.keys())].b0)[0,1]))
X['specs']={s:run(s) for s in('base','growth','lowess02','lowess05','calendar','lag2')}
# 4. buyuk El Nino'lar icin ust uste bindirilmis donem analizi (olay yili 0: ENSO yili)
EV=[1958,1966,1973,1983,1987,1992,1998,2010,2016,2024]; EVs=[1973,1983,1998,2016]
def sea(y,ev,lags=range(-2,4)):
    return {str(k):[float(np.mean([y.get(t+k,np.nan) for t in ev if (t+k) in y.index])),float(np.std([y.get(t+k,np.nan) for t in ev if (t+k) in y.index],ddof=1)/np.sqrt(sum((t+k) in y.index for t in ev)))] for k in lags}
W=m.groupby('yr').val.sum(); A=m[m.sp=='VET'].groupby('yr').val.sum().reindex(W.index).fillna(0)
X['sea']={'events':EV,'strong':EVs,'world':sea(detr(W),EV),'world_ex':sea(detr(W-A),EV),'world_strong':sea(detr(W),EVs),'world_ex_strong':sea(detr(W-A),EVs)}
for c in('Peru','Chile','Ecuador','Japan','Indonesia','Norway','United States of America','China','Morocco'):
    un=C[C.country==c].index[0]; X['sea'][c]=sea(detr(tot[un]),EV); X['sea'][c+'_strong']=sea(detr(tot[un]),EVs)
# ulkeler arasi: guclu El Nino yilinda (0) negatif sapma gosteren ulke payi ve ortanca
Y=pd.DataFrame({un:detr(s) for un,s in tot.items()})
X['sea']['share_neg']={str(t):float((Y.loc[t]<0).mean()) for t in EVs}; X['sea']['med_dev']={str(t):float(Y.loc[t].median()) for t in EVs}
Yz=Y/Y.std(); X['sea']['share_neg_all']=float((Y.loc[EV]<0).values.mean()); X['sea']['z_mean_ev']=float(Yz.loc[EV].values.mean()); 
# binom testi: olay yillarinda negatif sapma payi 0.5'ten farkli mi (yil kumeleri ile bootstrap)
rng=np.random.default_rng(3); yrs=np.array(Y.index); base=[(Y.loc[rng.choice(yrs[(yrs>=1952)],len(EV),replace=False)]<0).values.mean() for _ in range(2000)]
X['sea']['null_share_neg']=[float(np.percentile(base,2.5)),float(np.percentile(base,97.5))]
# 5. pozitif kontroller: bilinen ENSO duyarli stoklar
S=pd.read_csv('e/out/stock_enso.csv'); 
def one(country,spname,minmean=0):
    g=m[(m.country==country)&(m.sp_name.str.contains(spname,regex=False))]
    if g.empty: return None
    s=g.groupby('yr').val.sum(); y=detr(s)
    if y is None: return None
    y=y.loc[1952:2024]; Xm=sm.add_constant(np.column_stack([[E[t] for t in y.index],[E[t-1] for t in y.index]])); mm=sm.OLS(y.values,Xm).fit(cov_type='HAC',cov_kwds={'maxlags':3}); w=mm.wald_test(np.eye(3)[1:],scalar=True)
    return dict(b0=float(mm.params[1]),p0=float(mm.pvalues[1]),b1=float(mm.params[2]),p1=float(mm.pvalues[2]),r2=float(mm.rsquared),pj=float(w.pvalue),mean=float(s.loc[2000:2024].mean()))
PC=[('Peru','Anchoveta'),('Chile','Anchoveta'),('Chile','Araucanian herring'),('United States of America','Opalescent inshore squid'),('Ecuador','Skipjack tuna'),('Peru','Jumbo flying squid'),('Peru','Chilean jack mackerel'),('Chile','Chilean jack mackerel'),('Mexico','Pacific sardine'),('United States of America','Pacific sardine'),('Japan','Skipjack tuna'),('Papua New Guinea','Skipjack tuna'),('Kiribati','Skipjack tuna')]
X['controls']={c+' | '+s:one(c,s) for c,s in PC}
# 6. Dogu Pasifik bolgeleri icin gecikme yapisi: FAO 87 alan toplam, anchoveta haric
a87=m[m.area=='87']; s=a87.groupby('yr').val.sum(); s2=a87[a87.sp!='VET'].groupby('yr').val.sum()
def lagfit(s):
    y=detr(s).loc[1953:2024]; Xm=sm.add_constant(np.column_stack([[E[t] for t in y.index],[E[t-1] for t in y.index],[E[t-2] for t in y.index]])); mm=sm.OLS(y.values,Xm).fit(cov_type='HAC',cov_kwds={'maxlags':3}); return dict(b=[float(v) for v in mm.params[1:]],p=[float(v) for v in mm.pvalues[1:]],r2=float(mm.rsquared))
X['area87']=dict(all=lagfit(s),ex_anch=lagfit(s2),share_anch=float(1-s2.loc[2000:2024].mean()/s.loc[2000:2024].mean()))
R['x']=X; json.dump(R,open('e/out/results.json','w'),ensure_ascii=False,indent=1)
rd=lambda d:{k:(round(v,4) if isinstance(v,float) else v) for k,v in d.items()}
for k in('meta_all','meta_b1','meta_ex_peru','meta_epac','meta_wpac','meta_big','power'): print(k,rd(X[k]))
print({k:rd(v) for k,v in X['meta_cont'].items()}); print(X['specs'])
for k in('world','world_ex','world_strong','world_ex_strong','Peru','Peru_strong','Chile_strong','Ecuador_strong'): print('SEA',k,{a:[round(100*b[0],1),round(100*b[1],1)] for a,b in X['sea'][k].items()})
print(X['sea']['share_neg'],X['sea']['med_dev'],X['sea']['share_neg_all'],X['sea']['null_share_neg'],X['sea']['z_mean_ev'])
for k,v in X['controls'].items(): print('PC',k,rd(v) if v else None)
print(X['area87'])
'''

# ==============================================================================
# e/e_extra_b.py
# ==============================================================================
MODULES['e/e_extra_b.py'] = r'''
# Ek ozetler: evre asimetrisi (El Nino / La Nina) ve coklu sinama denetiminden gecen stoklar
import pandas as pd, json
R=json.load(open('e/out/results.json')); C=pd.read_csv('e/out/country_enso.csv'); S=pd.read_csv('e/out/stock_enso.csv')
R['x']['asym']=dict(n_diff=int((C.p_diff<.05).sum()),n_pos=int((C.pp<.05).sum()),n_neg=int((C.pn<.05).sum()),bp_med=float(C.bp.median()),bn_med=float(C.bn.median()))
F=S[S.sig_fdr].sort_values('pj'); R['x']['stock_fdr']=[[r.country,r.sp,float(r['mean']),float(r.b0),float(r.p0),float(r.b1),float(r.p1),float(r.r2),float(r.pj)] for _,r in F.iterrows()]
json.dump(R,open('e/out/results.json','w'),ensure_ascii=False,indent=1); print(R['x']['asym'],len(R['x']['stock_fdr']))
'''

# ==============================================================================
# e/e_modes.py
# ==============================================================================
MODULES['e/e_modes.py'] = r'''
# Iklim kipleri atlasi: on bir iklim kipinin ulusal deniz avindaki izi (1952-2022), ortak yontem, alan anlamliligi
import pandas as pd, numpy as np, json, warnings, io, re; warnings.filterwarnings('ignore')
from statsmodels.nonparametric.smoothers_lowess import lowess
from scipy import stats
m=pd.read_parquet('f/out/marine.parquet'); Y0,Y1=1952,2022; YRS=np.arange(Y0,Y1+1); T=len(YRS)
def rd(fn,skip=1):
    rows=[]
    for l in open(fn):
        p=l.split()
        if len(p)==13 and re.fullmatch(r'\d{4}',p[0]):
            try: rows.append([int(p[0])]+[float(x) for x in p[1:]])
            except ValueError: pass
    d=pd.DataFrame(rows).set_index(0); d=d.where((d>-90)&(d.abs()<9000)); s=d.stack(); s.index=pd.PeriodIndex([f'{y}-{k:02d}' for y,k in s.index],freq='M'); return s.sort_index()
def win(s,m0,m1,prev):   # prev=True: pencere onceki yilin m0 ayindan bu yilin m1 ayina
    out={}
    for y in range(1950,Y1+1):
        a=f'{y-1 if prev else y}-{m0:02d}'; b=f'{y}-{m1:02d}'; v=s[a:b]
        n=(12-m0+1+m1) if prev else (m1-m0+1)
        out[y]=float(v.mean()) if len(v)==n else np.nan
    return pd.Series(out)
oni=pd.read_csv('data/oni.data',sep=r'\s+',skiprows=1,header=None,nrows=77).set_index(0); oni=oni.where(oni>-90); mo=oni.stack(); mo.index=pd.PeriodIndex([f'{int(y)}-{int(k):02d}' for y,k in mo.index],freq='M'); mo=mo.sort_index()
I={'ENSO (ONI)':win(mo,7,6,True),'Doğu Pasifik El Niño (Niño 1+2)':win(rd('data/idx/nina1.anom.data'),7,6,True),'Orta Pasifik El Niño (Niño 4)':win(rd('data/idx/nina4.anom.data'),7,6,True),
   'PDO':win(rd('data/idx/pdo.data'),1,12,False),'Kuzey Pasifik (NP)':win(rd('data/idx/np.data'),11,3,True),'PNA':win(rd('data/idx/pna.data'),12,2,True),'NAO':win(rd('data/idx/nao.data'),12,3,True),'AMO':win(rd('data/idx/amon.us.data'),1,12,False),
   'Tropikal Kuzey Atlantik (TNA)':win(rd('data/idx/tna.data'),1,12,False),'Tropikal Güney Atlantik (TSA)':win(rd('data/idx/tsa.data'),1,12,False),'Hint Okyanusu Dipolü (DMI)':win(rd('data/idx/dmi.data'),6,11,False)}
def hp(v): v=np.asarray(v,float); return v-lowess(v,np.arange(len(v)),frac=0.3,return_sorted=False)
IX={}
for k,s in I.items():
    s=s.loc[1951:Y1]; assert s.notna().all(),(k,s[s.isna()].index.tolist()); z=hp(s.values); z=z/z.std(); IX[k]=pd.Series(z,index=s.index)
# ulke serileri (egilimden arindirilmis log av), 1950-2024 tam, ortalama >= 1000 t
def detr(s):
    s=s.reindex(range(1950,2025))
    if s.isna().any() or (s<=0).any(): return None
    l=np.log(s.values); return pd.Series(l-lowess(l,np.arange(len(l)),frac=0.3,return_sorted=False),index=s.index)
cols={}; meta={}
for un,g in m.groupby('un'):
    s=g.groupby('yr').val.sum(); y=detr(s)
    if y is None or s.mean()<1000: continue
    cols[un]=y.loc[Y0:Y1].values; i=g.iloc[0]; a=g[g.yr>=2000].groupby('area_name').val.sum(); meta[un]=dict(country=i.country,cont=i.cont,mean=float(s.loc[2000:2024].mean()),area=a.idxmax())
Y=np.column_stack(list(cols.values())); uns=list(cols); N=Y.shape[1]; W=np.array([meta[u]['mean'] for u in uns]); W=W/W.sum(); areas=np.array([meta[u]['area'] for u in uns])
def basin(a):
    if 'Pacific' in a: return 'Pasifik'
    if 'Indian' in a: return 'Hint'
    return 'Atlantik ve Akdeniz'
bas=np.array([basin(a) for a in areas])
def fitall(e0,e1,Y,L=3):
    X=np.column_stack([np.ones(T),e0,e1]); XtXi=np.linalg.inv(X.T@X); B=XtXi@X.T@Y; U=Y-X@B
    # Newey-West, sutun bazinda vektorize
    se=np.zeros((3,Y.shape[1])); S=np.einsum('ti,tj,tn->ijn',X,X,U*U)
    for l in range(1,L+1):
        w=1-l/(L+1); G=np.einsum('ti,tj,tn->ijn',X[l:],X[:-l],U[l:]*U[:-l]); S+=w*(G+G.transpose(1,0,2))
    V=np.einsum('ab,bcn,cd->adn',XtXi,S,XtXi)*T/(T-3)
    t0=B[1]/np.sqrt(V[1,1]); t1=B[2]/np.sqrt(V[2,2])
    # ortak Wald (b0=b1=0)
    Vi=V[1:,1:]; det=Vi[0,0]*Vi[1,1]-Vi[0,1]*Vi[1,0]; wald=(B[1]**2*Vi[1,1]-2*B[1]*B[2]*Vi[0,1]+B[2]**2*Vi[0,0])/det
    return B[1],B[2],t0,t1,wald,np.sqrt(V[1,1])
def summ(e,Y):
    e0=e.loc[Y0:Y1].values; e1=e.loc[Y0-1:Y1-1].values; b0,b1,t0,t1,w,se0=fitall(e0,e1,Y); p0=2*stats.t.sf(np.abs(t0),T-3); pj=stats.chi2.sf(w,2); return b0,b1,p0,pj,se0
def bh(p,q=0.05):
    o=np.sort(p); k=np.arange(1,len(p)+1); ok=o<=q*k/len(p); return int(k[ok].max()) if ok.any() else 0
rng=np.random.default_rng(42); R={}; rows=[]
for k,e in IX.items():
    b0,b1,p0,pj,se0=summ(e,Y); n0=int((p0<.05).sum()); nj=int((pj<.05).sum())
    nul0=[];nulj=[];nulw=[]
    for it in range(1000):
        sh=rng.integers(5,T-5); ev=pd.Series(np.roll(e.values,sh),index=e.index); _,_,q0,qj,_=summ(ev,Y); nul0.append((q0<.05).sum()); nulj.append((qj<.05).sum()); nulw.append(W[qj<.05].sum())
    # rastgele etkiler meta (DL)
    w=1/se0**2; mu=(w*b0).sum()/w.sum(); Q=(w*(b0-mu)**2).sum(); tau2=max(0,(Q-(N-1))/(w.sum()-(w**2).sum()/w.sum())); wr=1/(se0**2+tau2); mur=(wr*b0).sum()/wr.sum(); ser=np.sqrt(1/wr.sum())
    o=dict(n0=n0,nj=nj,exp=0.05*N,fdr0=bh(p0),fdrj=bh(pj),p_field0=float((np.array(nul0)>=n0).mean()),p_fieldj=float((np.array(nulj)>=nj).mean()),null_mean_j=float(np.mean(nulj)),null_p95_j=float(np.percentile(nulj,95)),wshare_j=float(W[pj<.05].sum()),p_wshare=float((np.array(nulw)>=W[pj<.05].sum()).mean()),
        meta_mu=float(mur),meta_se=float(ser),tau=float(np.sqrt(tau2)),med_abs_b0=float(np.median(np.abs(b0))),ac1=float(pd.Series(e.values).autocorr()),
        by_basin={b:dict(n=int((bas==b).sum()),nj=int(((pj<.05)&(bas==b)).sum())) for b in np.unique(bas)})
    top=np.argsort(pj)[:6]; o['top']=[[meta[uns[i]]['country'],float(b0[i]),float(b1[i]),float(pj[i]),float(meta[uns[i]]['mean'])] for i in top]
    R[k]=o
    for i,u in enumerate(uns): rows.append((k,u,meta[u]['country'],bas[i],b0[i],b1[i],p0[i],pj[i],meta[u]['mean']))
D=pd.DataFrame(rows,columns=['mode','un','country','basin','b0','b1','p0','pj','mean']); D.to_csv('e/out/modes_country.csv',index=False)
# kipler arasi korelasyon
CM=pd.DataFrame(IX).loc[Y0:Y1].corr(); CM.to_csv('e/out/modes_corr.csv')
# cok kipli ortak model: her ulke icin 12 kipin (es zamanli) birlikte aciklama gucu ve en iyi kip
E0=pd.DataFrame(IX).loc[Y0:Y1].values; Xf=np.column_stack([np.ones(T),E0]); Bf=np.linalg.lstsq(Xf,Y,rcond=None)[0]; r2=1-((Y-Xf@Bf)**2).sum(0)/((Y-Y.mean(0))**2).sum(0); adj=1-(1-r2)*(T-1)/(T-Xf.shape[1])
nulr=[]
for it in range(300):
    Es=np.column_stack([np.roll(E0[:,j],rng.integers(5,T-5)) for j in range(E0.shape[1])]); Xs=np.column_stack([np.ones(T),Es]); Bs=np.linalg.lstsq(Xs,Y,rcond=None)[0]; nulr.append(np.median(1-((Y-Xs@Bs)**2).sum(0)/((Y-Y.mean(0))**2).sum(0)))
best=D.loc[D.groupby('un').pj.idxmin()]; 
ALL=dict(N=int(N),T=int(T),years=[Y0,Y1],r2_med=float(np.median(r2)),adj_med=float(np.median(adj)),r2_null_med=float(np.mean(nulr)),r2_null_p95=float(np.percentile(nulr,95)),p_r2=float((np.array(nulr)>=np.median(r2)).mean()),share_adj_gt10=float((adj>0.10).mean()),
    best_mode_counts={k:int(v) for k,v in best['mode'].value_counts().items()},n_any_sig=int((D.groupby('un').pj.min()<0.05).sum()),n_any_bonf=int((D.groupby('un').pj.min()<0.05/len(IX)).sum()),total_tests=int(len(D)),total_sig=int((D.pj<.05).sum()),total_fdr=bh(D.pj.values),
    corr_ENSO={k:float(CM.loc['ENSO (ONI)',k]) for k in IX})
json.dump({'modes':R,'all':ALL},open('e/out/modes.json','w'),ensure_ascii=False,indent=1)
for k,o in R.items(): print(f"{k[:28]:28s} n0={o['n0']:3d} nj={o['nj']:3d} fdrj={o['fdrj']:2d} pfield={o['p_fieldj']:.3f} null={o['null_mean_j']:.1f}/{o['null_p95_j']:.0f} wshare={o['wshare_j']:.2f}(p={o['p_wshare']:.2f}) mu={100*o['meta_mu']:.2f}({100*o['meta_se']:.2f}) tau={100*o['tau']:.1f} ac1={o['ac1']:.2f} {o['by_basin']}  top={[(t[0][:14],round(100*t[1],1),round(t[3],4)) for t in o['top'][:3]]}")
print(ALL)
# --- ortak sifir dagilimi: butun kipler ayni miktarda kaydirilir (kipler arasi iliski korunur)
def allp(shift):
    out=[]
    for k,e in IX.items():
        ev=e if shift==0 else pd.Series(np.roll(e.values,shift),index=e.index); out.append(summ(ev,Y)[3])
    return np.array(out)   # kip x ulke
P0=allp(0); obs_tot=int((P0<.05).sum()); obs_any=int((P0.min(0)<.05).sum()); obs_min=np.sort(P0.min(0))
nt=[];na=[]
for sh in range(5,T-5):
    Ps=allp(sh); nt.append((Ps<.05).sum()); na.append((Ps.min(0)<.05).sum())
J=json.load(open('e/out/modes.json')); J['all'].update(null_total_mean=float(np.mean(nt)),null_total_p95=float(np.percentile(nt,95)),p_total=float((np.array(nt)>=obs_tot).mean()),null_any_mean=float(np.mean(na)),null_any_p95=float(np.percentile(na,95)),p_any=float((np.array(na)>=obs_any).mean()),n_shifts=len(nt),naive_expected_total=0.05*P0.size,naive_any=float(N*(1-0.95**len(IX))))
# kip basina yanlis pozitif orani (kaydirma sifir dagiliminda %5 duzeyinde anlamli bulunan ulke payi)
for k in J['modes']: J['modes'][k]['fpr']=J['modes'][k]['null_mean_j']/N
# Peru: ENSO cesidi
pe=[i for i,u in enumerate(uns) if meta[u]['country']=='Peru'][0]; e0=pd.DataFrame(IX).loc[Y0:Y1]; X=np.column_stack([np.ones(T),e0['Doğu Pasifik El Niño (Niño 1+2)'],e0['Orta Pasifik El Niño (Niño 4)']]); import statsmodels.api as sm; mm=sm.OLS(Y[:,pe],X).fit(cov_type='HAC',cov_kwds={'maxlags':3}); J['peru']=dict(b_ep=float(mm.params[1]),t_ep=float(mm.tvalues[1]),b_cp=float(mm.params[2]),t_cp=float(mm.tvalues[2]),r2=float(mm.rsquared))
tot=np.log(m.groupby('yr').val.sum().reindex(range(1950,2025)).values); wy=(tot-lowess(tot,np.arange(75),frac=0.3,return_sorted=False))[Y0-1950:Y1-1950+1]; mm=sm.OLS(wy,X).fit(cov_type='HAC',cov_kwds={'maxlags':3}); J['world_flavor']=dict(b_ep=float(mm.params[1]),t_ep=float(mm.tvalues[1]),b_cp=float(mm.params[2]),t_cp=float(mm.tvalues[2]),r2=float(mm.rsquared))
J['world_modes']={}
for k,e in IX.items():
    Xk=np.column_stack([np.ones(T),e.loc[Y0:Y1].values,e.loc[Y0-1:Y1-1].values]); mm=sm.OLS(wy,Xk).fit(cov_type='HAC',cov_kwds={'maxlags':3}); J['world_modes'][k]=dict(b0=float(mm.params[1]),p0=float(mm.pvalues[1]),b1=float(mm.params[2]),p1=float(mm.pvalues[2]),r2=float(mm.rsquared))
json.dump(J,open('e/out/modes.json','w'),ensure_ascii=False,indent=1)
print({k:v for k,v in J['all'].items() if 'null' in k or k.startswith('p_') or 'naive' in k},obs_tot,obs_any); print(J['peru'],J['world_flavor']); print({k[:12]:(round(100*v['b0'],2),round(v['p0'],3),round(v['r2'],3)) for k,v in J['world_modes'].items()})
'''

# ==============================================================================
# e/e_fig.py
# ==============================================================================
MODULES['e/e_fig.py'] = r'''
# -*- coding: utf-8 -*-
import warnings; warnings.filterwarnings('ignore')
import json, matplotlib; matplotlib.use('Agg'); import matplotlib.pyplot as plt
import numpy as np, pandas as pd
from scipy import stats
from statsmodels.nonparametric.smoothers_lowess import lowess
plt.rcParams.update({'font.family':'DejaVu Sans','font.size':9,'axes.spines.top':False,'axes.spines.right':False,'savefig.dpi':300})
R=json.load(open('e/out/results.json')); X=R['x']; C=pd.read_csv('e/out/country_enso.csv'); m=pd.read_parquet('f/out/marine.parquet')
B='#1b4f72'; O='#b9770e'; Rd='#b03a2e'; Gn='#1e8449'
E=pd.Series({int(k):v for k,v in R['E'].items()})
def detr(s):
    l=np.log(s.values); return pd.Series(l-lowess(l,np.arange(len(l)),frac=0.3,return_sorted=False),index=s.index)
W=m.groupby('yr').val.sum(); A=m[m.sp=='VET'].groupby('yr').val.sum().reindex(W.index).fillna(0)
# S1 ENSO ve dunya avi sapmasi
fig,ax=plt.subplots(2,1,figsize=(7.4,4.2),sharex=True)
a=ax[0]; a.bar(E.index,E.values,color=[Rd if v>=0.5 else (B if v<=-0.5 else '#bbbbbb') for v in E.values],width=.85); a.axhline(0,color='#999',lw=.6); a.set_ylabel('ENSO indeksi (°C)'); a.set_title('(a) ENSO yılı indeksi (Temmuz–Haziran ortalaması)',fontsize=9,loc='left')
a=ax[1]; d1=100*detr(W); d2=100*detr(W-A); a.plot(d1.index,d1,color=B,lw=1.3,label='Dünya toplam avı'); a.plot(d2.index,d2,color=O,lw=1.1,label='Perulu hamsi dışında'); a.axhline(0,color='#999',lw=.6)
for y in X['sea']['strong']: a.axvline(y,color=Rd,lw=.6,ls=':')
a.set_ylabel('Eğilimden sapma (%)'); a.legend(frameon=False,fontsize=7.5,ncol=2,loc='lower left'); a.set_title('(b) Dünya deniz avının eğilimden sapması',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('e/fig/Figure_1.png'); plt.close()
# S2 ulke t-istatistikleri ve huni grafigi
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1))
a=ax[0]; t=C.b0/C.se0; a.hist(t,bins=np.arange(-5,5.01,.5),density=True,color='#aed6f1',edgecolor='white',lw=.3,label='Ülkeler'); g=np.linspace(-5,5,200); a.plot(g,stats.norm.pdf(g),color=Rd,lw=1.1,label='Standart normal'); a.axvline(-1.96,color='#999',lw=.6,ls='--'); a.axvline(1.96,color='#999',lw=.6,ls='--')
a.set_xlabel('t istatistiği (eş zamanlı ENSO katsayısı)'); a.set_ylabel('Yoğunluk'); a.legend(frameon=False,fontsize=7.5); a.set_title('(a) Ülke katsayılarının dağılımı',fontsize=9,loc='left')
a=ax[1]; a.scatter(100*C.b0,C.se0*100,s=4+np.sqrt(C['mean'])/60,color=B,alpha=.55,lw=0); yy=np.linspace(0.3,20,50); a.plot(1.96*yy,yy,color='#999',lw=.7,ls='--'); a.plot(-1.96*yy,yy,color='#999',lw=.7,ls='--'); a.axvline(0,color='#999',lw=.6)
for nm in('Peru','Ecuador','Chile'):
    r=C[C.country==nm].iloc[0]; a.annotate({'Peru':'Peru','Ecuador':'Ekvador','Chile':'Şili'}[nm],(100*r.b0,100*r.se0),fontsize=7,xytext=(3,3),textcoords='offset points')
a.set_ylim(20,0); a.set_xlim(-30,30); a.set_xlabel('ENSO katsayısı (ONI birimi başına %)'); a.set_ylabel('Standart hata (%)'); a.set_title('(b) Huni grafiği',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('e/fig/Figure_2.png'); plt.close()
# S3 ust uste bindirilmis donem analizi
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1)); ks=[-2,-1,0,1,2,3]
a=ax[0]
for key,c,l in(('world_strong',B,'Dünya toplamı'),('world_ex_strong',O,'Perulu hamsi dışında')):
    v=np.array([X['sea'][key][str(k)] for k in ks])*100; a.errorbar(ks,v[:,0],yerr=v[:,1],fmt='o-',color=c,ms=3.5,lw=1.1,capsize=2,label=l)
a.axhline(0,color='#999',lw=.6); a.axvline(0,color='#999',lw=.6,ls=':'); a.set_xlabel('Güçlü El Niño yılına göre yıl'); a.set_ylabel('Eğilimden sapma (%)'); a.legend(frameon=False,fontsize=7.5); a.set_title('(a) Dünya',fontsize=9,loc='left')
a=ax[1]
for key,c,l in(('Peru_strong',Rd,'Peru'),('Chile_strong',O,'Şili'),('Ecuador_strong',Gn,'Ekvador'),('Japan_strong',B,'Japonya'),('Norway_strong','#6c3483','Norveç')):
    v=np.array([X['sea'][key][str(k)] for k in ks]); a.plot(ks,100*(np.exp(v[:,0])-1),'o-',color=c,ms=3,lw=1.1,label=l)
a.axhline(0,color='#999',lw=.6); a.axvline(0,color='#999',lw=.6,ls=':'); a.set_xlabel('Güçlü El Niño yılına göre yıl'); a.legend(frameon=False,fontsize=7,ncol=2); a.set_title('(b) Seçilmiş ülkeler',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('e/fig/Figure_3.png'); plt.close()
# S4 FAO alanlari orman grafigi (b0 ve b1)
AR=R['area']; TRA={'Pacific, Southeast':'Güneydoğu Pasifik','Pacific, Eastern Central':'Doğu Orta Pasifik','Pacific, Northeast':'Kuzeydoğu Pasifik','Pacific, Northwest':'Kuzeybatı Pasifik','Pacific, Western Central':'Batı Orta Pasifik','Pacific, Southwest':'Güneybatı Pasifik','Indian Ocean, Eastern':'Doğu Hint Okyanusu','Indian Ocean, Western':'Batı Hint Okyanusu','Atlantic, Northwest':'Kuzeybatı Atlantik','Atlantic, Northeast':'Kuzeydoğu Atlantik','Atlantic, Western Central':'Batı Orta Atlantik','Atlantic, Eastern Central':'Doğu Orta Atlantik','Atlantic, Southwest':'Güneybatı Atlantik','Atlantic, Southeast':'Güneydoğu Atlantik','Mediterranean and Black Sea':'Akdeniz ve Karadeniz'}
ks_=[k for k in TRA if k in AR]; y=np.arange(len(ks_))[::-1]
fig,ax=plt.subplots(figsize=(6.6,4.0))
for i,(b,s,c,l) in enumerate((('b0','se0',B,'Aynı yıl'),('b1','se1',O,'Bir yıl gecikmeli'))):
    ax.errorbar([100*AR[k][b] for k in ks_],y+0.17-i*0.34,xerr=[196*AR[k][s] for k in ks_],fmt='o',color=c,ms=3.5,capsize=2,lw=.9,label=l)
ax.axvline(0,color='#999',lw=.6); ax.set_yticks(y); ax.set_yticklabels([TRA[k] for k in ks_],fontsize=8); ax.set_xlabel('ENSO katsayısı (ONI birimi başına %, %95 güven aralığı)'); ax.legend(frameon=False,fontsize=8)
fig.tight_layout(); fig.savefig('e/fig/Figure_4.png'); plt.close()
# S5 en kucuk saptanabilir etki dagilimi ve pozitif kontroller
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1))
a=ax[0]; a.hist(100*2.8*C.se0.clip(upper=0.25),bins=30,color='#aed6f1',edgecolor='white',lw=.3); a.axvline(100*X['power']['mde_med'],color=Rd,lw=1,ls='--'); a.set_xlabel('En küçük saptanabilir etki (ONI birimi başına %)'); a.set_ylabel('Ülke sayısı'); a.set_title('(a) İstatistiksel güç',fontsize=9,loc='left')
a=ax[1]; PC=[('Chile | Anchoveta','Şili: Perulu hamsi'),('United States of America | Opalescent inshore squid','ABD: Kaliforniya kalamarı'),('Chile | Araucanian herring','Şili: Arauko ringası'),('Peru | Anchoveta','Peru: Perulu hamsi'),('Mexico | Pacific sardine','Meksika: Pasifik sardalyası'),('Ecuador | Skipjack tuna','Ekvador: çizgili orkinos'),('Japan | Skipjack tuna','Japonya: çizgili orkinos'),('Chile | Chilean jack mackerel','Şili: istavrit')]
S=pd.read_csv('e/out/stock_enso.csv').set_index(['country','sp'])
yy=np.arange(len(PC))[::-1]; v=[100*X['controls'][k]['b0'] for k,_ in PC]; pv=[X['controls'][k]['p0'] for k,_ in PC]
a.barh(yy,v,color=[Rd if p<0.05 else ('#e59866' if p<0.1 else '#bbbbbb') for p in pv],height=.6); a.set_yticks(yy); a.set_yticklabels([l for _,l in PC],fontsize=7.5); a.axvline(0,color='#999',lw=.6); a.set_xlabel('ENSO katsayısı (ONI birimi başına %)'); a.set_title('(b) Bilinen ENSO duyarlı balıkçılıklar',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('e/fig/Figure_5.png'); plt.close()
# S6 iklim kipleri atlasi
J=json.load(open('e/out/modes.json')); MO=J['modes']; ks=list(MO); short={'ENSO (ONI)':'ENSO (ONI)','Doğu Pasifik El Niño (Niño 1+2)':'Niño 1+2 (doğu)','Orta Pasifik El Niño (Niño 4)':'Niño 4 (orta)','PDO':'PDO','Kuzey Pasifik (NP)':'NP','PNA':'PNA','NAO':'NAO','AMO':'AMO','Tropikal Kuzey Atlantik (TNA)':'TNA','Tropikal Güney Atlantik (TSA)':'TSA','Hint Okyanusu Dipolü (DMI)':'DMI'}
fig,ax=plt.subplots(1,2,figsize=(7.4,3.4))
a=ax[0]; y=np.arange(len(ks))[::-1]; a.barh(y,[MO[k]['nj'] for k in ks],color=B,height=.55,label='Gözlenen'); a.scatter([MO[k]['null_mean_j'] for k in ks],y,marker='|',color='#222',s=90,zorder=3,label='Sıfır dağılımı ortalaması'); a.scatter([MO[k]['null_p95_j'] for k in ks],y,marker='|',color=Rd,s=90,zorder=3,label='Sıfır dağılımı, %95'); a.axvline(0.05*J['all']['N'],color='#999',lw=.7,ls=':')
a.set_yticks(y); a.set_yticklabels([short[k] for k in ks],fontsize=7.5); a.set_xlabel('Anlamlı ilişki bulunan ülke sayısı (p < 0,05)'); a.legend(frameon=False,fontsize=6.8,loc='lower right'); a.set_title('(a) Ülke sayısı ve alan anlamlılığı',fontsize=9,loc='left')
a=ax[1]; wm=J['world_modes']; a.barh(y,[100*wm[k]['r2'] for k in ks],color=[Rd if wm[k]['p0']<0.05 else '#aab7b8' for k in ks],height=.55); a.set_yticks(y); a.set_yticklabels([]); a.set_xlabel('Dünya toplam avında açıklanan değişkenlik (R², %)'); a.set_title('(b) Dünya toplam avı',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('e/fig/Figure_6.png'); plt.close()
'''

if __name__ == "__main__":
    for _d in ('e', 'e/out', 'e/fig'): os.makedirs(_d, exist_ok=True)
    for _m in ORDER:
        with open(_m, 'w', encoding='utf-8') as _f: _f.write(MODULES[_m].lstrip('\n'))
    for _m in ORDER:
        if _m in SKIP_IF_EXISTS and os.path.exists(SKIP_IF_EXISTS[_m]):
            print('== atlandi / skipped (cikti mevcut / output present):', _m); continue
        print('==', _m, flush=True); subprocess.run([sys.executable, _m], check=True)   # her modul ayri surecte calisir / each module runs in its own process
    print('Bitti / done. Sonuclar / results: e/out/results.json')