# -*- coding: utf-8 -*-
# ------------------------------------------------------------------------------
# Kitlesel ağarmalardan sonra resif balığı avı düşüyor mu? Bildirilen ve yeniden yapılandırılmış av verileriyle 51 ülke ve 95 münhasır ekonomik bölgede sınama, 1986–2024
# Do reef fish catches fall after mass bleaching? A test with reported and reconstructed catch data in 51 countries and 95 exclusive economic zones, 1986–2024
#
# GOMD Yayınlar · Makale 24 — https://gomd.pirireis.edu.tr/publications.php?id=rp24
#
# Bu dosya makalenin bütün modüllerini içerir. Replikasyon paketi klasöründe çalıştırıldığında modülleri `r/` klasörüne yazar ve sırayla, her birini ayrı süreçte çalıştırır;
# veriyi `data/` ve `r/out/` klasörlerinden okur, sonuçları `r/out/` ve şekilleri `r/fig/` klasörüne yazar.
# This file contains all modules of the paper. Run in the replication-package folder, it writes the modules to `r/` and runs them in order, each in its own process;
# it reads data from `data/` and `r/out/`, and writes results to `r/out/` and figures to `r/fig/`.
#
# Kurulum / Setup:  pip install pandas numpy matplotlib scipy statsmodels pyarrow shapely
# Çalıştırma / Run:  python rp24_reef_fish_catch_after_mass_bleaching.py
#
# Heat-stress data: NOAA Coral Reef Watch, Regional Virtual Stations, daily 5 km product suite v3.1 (https://coralreefwatch.noaa.gov). Catch data: FAO Global Capture Production (FishStat, release 2026.1.0), © FAO, https://www.fao.org/fishery/en/collection/capture . Reconstructed catch: Sea Around Us (https://www.seaaroundus.org; Pauly, Zeller and Palomares, eds.), retrieved through the public API by the script itself; these data are not redistributed in this package.
# 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 = ['r/r_prep.py', 'r/r_analysis.py', 'r/r_extra.py', 'r/r_extra2.py', 'r/r_extra3.py', 'r/r_sau_prep.py', 'r/r_sau.py', 'r/r_sau2.py', 'r/r_fig.py']
SKIP_IF_EXISTS = {'r/r_sau_prep.py': 'r/out/sau_raw.parquet'}   # hazırlık modülü: çıktısı pakette varsa atlanır / preparation module: skipped if its output is in the package

# ==============================================================================
# r/r_prep.py
# ==============================================================================
MODULES['r/r_prep.py'] = r'''
# Resif isi stresi -> ulusal resif baligi avi: veri hazirligi
import pandas as pd, numpy as np, json, re
s=pd.read_csv('c/out/stations.csv'); P=pd.read_csv('data/pw/ports_global.csv')
la=np.radians(s.lat.values)[:,None]; lo=np.radians(s.lon.values)[:,None]; pa=np.radians(P.lat.values)[None]; po=np.radians(P.lon.values)[None]
d=6371*2*np.arcsin(np.sqrt(np.sin((pa-la)/2)**2+np.cos(la)*np.cos(pa)*np.sin((po-lo)/2)**2)); s['iso']=P.ISO3.values[d.argmin(1)]
fix={'bermuda':'BMU','dongsha':'TWN','east_gulf_of_thailand':'THA','egypt':'EGY','el_salvador':'SLV','great_nicobar':'IND','nicaragua':'NIC','northern_cook_islands':'COK','northern_mauritius':'MUS','northern_vietnam':'VNM','phoenix_islands':'KIR','singapore':'SGP','socotra':'YEM','southern_line_islands':'KIR','southern_myanmar':'MMR','torres_strait':'AUS','wallis_futuna':'WLF',
 'leeward_caribbean':'ATG,KNA,MSR,AIA,GLP,SXM,MAF,BLM','windward_caribbean':'LCA,VCT,GRD,BRB,DMA,MTQ','abc_islands':'ABW,CUW,BES','usvi':'VIR,VGB','samoas':'WSM,ASM','hong_kong':'HKG','tromelin_reunion':'REU','mayotte':'MYT','gulf_of_aqaba':'JOR,EGY','eastern_persian_gulf':'IRN,ARE,QAT','western_persian_gulf':'SAU,KWT,QAT','gulf_of_oman':'OMN,ARE','southern_hispaniola':'HTI,DOM','northern_hispaniola':'DOM,HTI','honduras':'HND,GTM'}
drop=['bassas_da_india','chagos_archipelago','christmas_island','clipperton_island','easter_island','galapagos','glorieuse_geyser','howland_baker','isla_del_coco','johnston_atoll','norfolk_island','nwhi','paracel_islands','pitcairn','pulu_keeling','spratly_islands','lord_howe_island','ogasawara','trinadade_martim_vaz','wake_atoll']
s=s[~s.sid.isin(drop)].copy(); s['iso']=[fix.get(r.sid,r.iso) for r in s.itertuples()]
M=s.assign(iso=s.iso.str.split(',')).explode('iso')[['sid','iso']]; M.to_csv('r/out/station_country.csv',index=False)
c=pd.read_parquet('c/out/crw.parquet',columns=['sid','date','dhw','ssta']); c['yr']=pd.to_datetime(c.date).dt.year; c=c[(c.yr>=1985)&(c.yr<=2025)]
sy=c.groupby(['sid','yr']).agg(mx=('dhw','max'),ssta=('ssta','mean')).reset_index()
cy=M.merge(sy,on='sid').groupby(['iso','yr']).agg(dhw=('mx','mean'),dhw_max=('mx','max'),ssta=('ssta','mean'),nst=('sid','nunique')).reset_index(); cy.to_csv('r/out/country_dhw.csv',index=False)
m=pd.read_parquet('f/out/marine.parquet'); m=m[m.Major_Group.astype(str).str.upper().str.contains('PISCES')|True]
kw=r'snapper|jobfish|grouper|hind|coral ?trout|emperor|fusilier|spinefe|rabbitfish|parrotfish|surgeonfish|unicornfish|wrasse|hogfish|triggerfish|goatfish|sweetlips|grunt|squirrelfish|soldierfish|angelfish|butterflyfish|coralgroup|seabasses nei|coney|graysby|tang|humphead|bumphead|rockcod'
nm=m.sp_name.astype(str); isreef=nm.str.contains(kw,case=False,regex=True)&m.isscaap.isin(['Miscellaneous coastal fishes','Miscellaneous demersal fishes'])&~nm.str.contains('Marbled rockcod|Striped-eyed|Humped rockcod|Antarctic|Patagonian|icefish',case=False)
print(m[isreef].groupby('sp_name').val.sum().sort_values(ascending=False).head(45).round(0).to_string())
m['reef']=isreef; m['grp']=np.where(isreef,'reef',np.where(m.isscaap=='Tunas, bonitos, billfishes','tuna',np.where(m.isscaap.isin(['Herrings, sardines, anchovies','Miscellaneous pelagic fishes']),'pel',np.where(m.isscaap=='Miscellaneous coastal fishes','coast','oth'))))
iso=cy.iso.unique(); F=m[m.iso.isin(iso)&(m.yr>=1975)]
C=F.groupby(['iso','country','yr','grp']).val.sum().unstack('grp').fillna(0).reset_index(); C['tot']=C[['reef','tuna','pel','coast','oth']].sum(axis=1)
off=F[F.reef].groupby(['iso','yr']).apply(lambda g:(g.val*g.official).sum()/max(g.val.sum(),1e-9)).rename('off').reset_index(); C=C.merge(off,on=['iso','yr'],how='left')
C.to_csv('r/out/country_catch.csv',index=False); print(C.iso.nunique(),len(C)); print(sorted(set(iso)-set(C.iso)))
'''

# ==============================================================================
# r/r_analysis.py
# ==============================================================================
MODULES['r/r_analysis.py'] = r'''
# Resif isi stresi ulusal resif baligi avinda gorunuyor mu? Yerel izdusum (local projection) paneli
import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore')
C=pd.read_csv('r/out/country_catch.csv'); D=pd.read_csv('r/out/country_dhw.csv'); R={}
C=C[(C.yr>=1980)&(C.yr<=2024)]
ok=C.groupby('iso').apply(lambda g:(g.reef>0).sum()>=30 and g.reef.mean()>=20); keep=ok[ok].index; C=C[C.iso.isin(keep)].copy()
C=C.merge(D,on=['iso','yr'],how='left'); C=C.sort_values(['iso','yr']).reset_index(drop=True)
for g in('reef','tuna','pel','coast','tot'): C['l'+g]=np.log(C[g].where(C[g]>0))
C['ev4']=(C.dhw>=4).astype(float); C['ev8']=(C.dhw>=8).astype(float); C['evm8']=(C.dhw_max>=8).astype(float)
R['sample']=dict(n_cty=int(C.iso.nunique()),yrs=[int(C.yr.min()),int(C.yr.max())],n_obs=int(len(C)),reef_share=float(C[C.yr>=2015].reef.sum()/C[C.yr>=2015].tot.sum()),reef_2015_24=float(C[C.yr>=2015].groupby('yr').reef.sum().mean()),
    n_ev4=int(C[C.yr>=1986].ev4.sum()),n_ev8=int(C[C.yr>=1986].ev8.sum()),n_cty_ev4=int(C[C.ev4==1].iso.nunique()),n_cty_ev8=int(C[C.ev8==1].iso.nunique()),med_dhw=float(C[C.yr>=1986].dhw.median()),p90_dhw=float(C[C.yr>=1986].dhw.quantile(.9)),
    off=float(C.off.mean()),ev_by_dec={str(k):[int(g.ev4.sum()),int(g.ev8.sum())] for k,g in C[C.yr>=1986].groupby((C.yr//10)*10)})
def lp(dep,x,H=range(-4,11),df=C,ctrl=True,y0=1986,wt=None,trend=False):
    out={}
    for h in H:
        d=df.copy(); g=d.groupby('iso')[dep]
        d['y']=g.shift(-h)-g.shift(1) if h>=0 else g.shift(-h)-g.shift(1)   # h<0: ornek oncesi (plasebo), t-1'e gore
        if h==-1: continue
        d['g1']=g.shift(1)-g.shift(2); d['g2']=g.shift(2)-g.shift(3)
        cols=[x]+(['g1','g2'] if ctrl else [])
        d=d[(d.yr>=y0)].dropna(subset=['y']+cols)
        d=d[np.isfinite(d.y)]
        # iki yonlu sabit etki: ulke ve yil kuklalari (dengesiz panel -> acik kuklalarla)
        X=pd.concat([d[cols],pd.get_dummies(d.iso,drop_first=True,dtype=float),pd.get_dummies(d.yr,drop_first=True,dtype=float)],axis=1); X.insert(0,'c',1.0)
        Xv=X.values; yv=d.y.values; w=np.ones(len(d)) if wt is None else d[wt].values
        sw=np.sqrt(w); b=np.linalg.lstsq(Xv*sw[:,None],yv*sw,rcond=None)[0]; e=yv-Xv@b
        XtXi=np.linalg.pinv((Xv*w[:,None]).T@Xv); meat=np.zeros((Xv.shape[1],)*2)
        for iso,idx in d.reset_index(drop=True).groupby('iso').indices.items():
            s=(Xv[idx]*w[idx,None]).T@e[idx]; meat+=np.outer(s,s)
        G=d.iso.nunique(); V=XtXi@meat@XtXi*G/(G-1); out[h]=dict(b=float(b[1]),se=float(np.sqrt(V[1,1])),n=int(len(d)),G=int(G))
    return out
C['dhw4']=C.dhw/4
R['lp_reef_dhw']=lp('lreef','dhw4'); R['lp_reef_ev4']=lp('lreef','ev4'); R['lp_reef_ev8']=lp('lreef','ev8'); R['lp_reef_evm8']=lp('lreef','evm8')
R['lp_tuna_ev4']=lp('ltuna','ev4'); R['lp_pel_ev4']=lp('lpel','ev4'); R['lp_tot_ev4']=lp('ltot','ev4'); R['lp_coast_ev4']=lp('lcoast','ev4')
R['lp_reef_ev4_noctrl']=lp('lreef','ev4',ctrl=False)
# resmi veri agirlikli ulkeler; resif bagimliligi yuksek ulkeler; ada ulkeleri
offc=C.groupby('iso').off.mean(); hi=offc[offc>=0.8].index; R['lp_reef_ev4_official']=lp('lreef','ev4',df=C[C.iso.isin(hi)]); R['n_official']=int(len(hi))
dep=C[C.yr>=2000].groupby('iso').apply(lambda g:g.reef.sum()/g.tot.sum()); hid=dep[dep>=dep.median()].index; R['lp_reef_ev4_dep']=lp('lreef','ev4',df=C[C.iso.isin(hid)]); R['dep_median']=float(dep.median())
big=C.groupby('iso').reef.mean(); hb=big[big>=big.median()].index; R['lp_reef_ev4_big']=lp('lreef','ev4',df=C[C.iso.isin(hb)])
# 1998 ve 2016 olaylari ayri
for y in(1998,2010,2016):
    d=C.copy(); d['e']=((d.yr==y)&(d.dhw>=4)).astype(float); R[f'lp_reef_{y}']=lp('lreef','e',H=range(0,min(9,2024-y+1)),df=d); R[f'n_{y}']=int(d.e.sum())
# doz-tepki: 5 yil sonraki degisim, DHW siniflari
d=C.copy(); g=d.groupby('iso').lreef; d['y5']=(g.shift(-3)+g.shift(-4)+g.shift(-5))/3-(g.shift(1)+g.shift(2)+g.shift(3))/3; d=d[d.yr>=1986].dropna(subset=['y5','dhw']); d['y5r']=d.y5-d.groupby('yr').y5.transform('mean')
bins=[-1,1,2,4,8,100]; d['bin']=pd.cut(d.dhw,bins,labels=['0-1','1-2','2-4','4-8','8+']); R['dose']={str(k):dict(n=int(len(g)),mean=float(g.y5r.mean()),se=float(g.y5r.std()/np.sqrt(len(g))),med=float(g.y5r.median())) for k,g in d.groupby('bin')}
d.to_csv('r/out/dose.csv',index=False)
# gurultu ve guc: ulke ici yillik log degisim sd; tespit edilebilir en kucuk etki (h=3..5 ortalamasi icin)
gg=C.groupby('iso').lreef.diff(); R['noise']=dict(sd_dlog=float(gg.std()),med_sd=float(C.assign(g=gg).groupby('iso').g.std().median()),ac1=float(C.assign(g=gg).groupby('iso').g.apply(lambda s:s.autocorr()).median()))
for k in('lp_reef_ev4','lp_reef_ev8','lp_reef_dhw'): R['mde_'+k]={str(h):float(2.8*R[k][h]['se']) for h in(0,3,5,8)}
# yil-ulke permutasyonu: ulke ici isi stresi yillarinin rastgele kaydirilmasi (h=3..5 ortalama katsayi)
def avg35(df,x): 
    o=lp('lreef',x,H=[3,4,5],df=df); return np.mean([o[h]['b'] for h in(3,4,5)])
obs=avg35(C,'ev4'); rng=np.random.default_rng(7); null=[]
for k in range(300):
    d=C.copy(); d['ev4']=d.groupby('iso').ev4.transform(lambda s:np.roll(s.values,rng.integers(1,len(s)))); null.append(avg35(d,'ev4'))
R['perm']=dict(obs=float(obs),p=float((np.abs(null)>=abs(obs)).mean()),null_sd=float(np.std(null)),lo=float(np.percentile(null,2.5)),hi=float(np.percentile(null,97.5)))
C.to_csv('r/out/panel.csv',index=False); json.dump(R,open('r/out/results.json','w'),ensure_ascii=False,indent=1)
print(R['sample']); print(R['noise'],R['perm'],R['dose'])
for k in R:
    if k.startswith('lp_'): print(k,' '.join(f"{h}:{100*v['b']:.1f}({100*v['se']:.1f})" for h,v in R[k].items()),'n=',list(R[k].values())[0]['n'],list(R[k].values())[0]['G'])
print({k:v for k,v in R.items() if k.startswith('mde') or k.startswith('n_')})
'''

# ==============================================================================
# r/r_extra.py
# ==============================================================================
MODULES['r/r_extra.py'] = r'''
import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore')
import statsmodels.formula.api as smf
Cc=pd.read_csv('r/out/panel.csv'); R=json.load(open('r/out/results.json')); X={}
# 1. doz-tepki, ulke kumeli (yil ortalamasindan arindirilmis 5 yillik degisim)
d=pd.read_csv('r/out/dose.csv'); d['d4']=d.dhw/4
m=smf.ols('y5~d4+C(yr)',d).fit(cov_type='cluster',cov_kwds={'groups':d.iso}); X['dose_lin']=dict(b=float(m.params['d4']),se=float(m.bse['d4']),p=float(m.pvalues['d4']),n=int(m.nobs))
d['b8']=(d.dhw>=8).astype(float); d['b48']=((d.dhw>=4)&(d.dhw<8)).astype(float)
m=smf.ols('y5~b48+b8+C(yr)',d).fit(cov_type='cluster',cov_kwds={'groups':d.iso}); X['dose_bins']={k:dict(b=float(m.params[k]),se=float(m.bse[k]),p=float(m.pvalues[k])) for k in('b48','b8')}
m=smf.ols('y5~b48+b8+C(yr)+C(iso)',d).fit(cov_type='cluster',cov_kwds={'groups':d.iso}); X['dose_bins_fe']={k:dict(b=float(m.params[k]),se=float(m.bse[k]),p=float(m.pvalues[k])) for k in('b48','b8')}
# 2. uzun fark: 1990-99 -> 2015-24 resif avi degisimi ve 2000-2019 birikimli isi stresi
a=Cc[(Cc.yr>=1990)&(Cc.yr<=1999)].groupby('iso').agg(r0=('reef','mean'),t0=('tot','mean')); b=Cc[(Cc.yr>=2015)&(Cc.yr<=2024)].groupby('iso').agg(r1=('reef','mean'),t1=('tot','mean'))
e=Cc[(Cc.yr>=1998)&(Cc.yr<=2019)].groupby('iso').agg(n8=('ev8','sum'),n4=('ev4','sum'),mdhw=('dhw','mean')); L=a.join(b).join(e).dropna(); L=L[(L.r0>0)&(L.r1>0)]; L['dl']=np.log(L.r1/L.r0); L['dsh']=L.r1/L.t1-L.r0/L.t0; L['dlo']=np.log((L.t1-L.r1)/(L.t0-L.r0))
L.to_csv('r/out/longdiff.csv')
for dep in('dl','dsh'):
    for x in('n8','mdhw'):
        m=smf.ols(f'{dep}~{x}',L).fit(cov_type='HC1'); X[f'ld_{dep}_{x}']=dict(b=float(m.params[x]),se=float(m.bse[x]),p=float(m.pvalues[x]),r2=float(m.rsquared),n=int(m.nobs))
m=smf.ols('dl~n8+dlo',L).fit(cov_type='HC1'); X['ld_dl_n8_ctrl']=dict(b=float(m.params['n8']),se=float(m.bse['n8']),p=float(m.pvalues['n8']),b_o=float(m.params['dlo']))
X['ld_desc']=dict(n=int(len(L)),med_dl=float(L.dl.median()),share_up=float((L.dl>0).mean()),med_n8=float(L.n8.median()),max_n8=float(L.n8.max()),sp=float(L[['dl','n8']].corr(method='spearman').iloc[0,1]))
hi=L[L.n8>=L.n8.median()]; lo=L[L.n8<L.n8.median()]; X['ld_groups']=dict(hi_n=int(len(hi)),lo_n=int(len(lo)),hi_med=float(hi.dl.median()),lo_med=float(lo.dl.median()))
# 3. toplam resif avi ve kuresel isi stresi
W=Cc.groupby('yr').agg(reef=('reef','sum'),tot=('tot','sum'),dhw=('dhw','mean'),ev8=('ev8','mean')).reset_index(); W.to_csv('r/out/world.csv',index=False)
X['world']={str(int(r.yr)):dict(reef=float(r.reef),share=float(r.reef/r.tot),dhw=float(r.dhw) if r.dhw==r.dhw else None,ev8=float(r.ev8)) for _,r in W.iterrows() if r.yr in(1980,1990,2000,2010,2016,2020,2024)}
X['reef_growth']=float(np.log(W[W.yr>=2015].reef.mean()/W[(W.yr>=1990)&(W.yr<=1999)].reef.mean()))
# 4. en buyuk ureticiler
top=Cc[Cc.yr>=2015].groupby(['iso','country']).reef.mean().sort_values(ascending=False).head(12).reset_index(); tot=Cc[Cc.yr>=2015].groupby('yr').reef.sum().mean()
X['top']=[[r.country,float(r.reef),float(r.reef/tot),int(Cc[(Cc.iso==r.iso)&(Cc.yr>=1986)].ev8.sum()),float(L.dl.get(r.iso,np.nan)) if r.iso in L.index else None] for r in top.itertuples()]
X['top5_share']=float(top.reef.head(5).sum()/tot)
# 5. bilinen olaylar: Seyseller 1998, Maldivler 1998/2016, vb. olay oncesi 3 yil / sonrasi 3-5 yil
def case(iso,y):
    g=Cc[Cc.iso==iso].set_index('yr').reef; pre=g.loc[y-3:y-1].mean(); post=g.loc[y+3:y+5].mean(); return dict(pre=float(pre),post=float(post),d=float(np.log(post/pre)),dhw=float(Cc[(Cc.iso==iso)&(Cc.yr==y)].dhw.iloc[0]))
cs={}
for iso,y in(('SYC',1998),('MDV',1998),('MDV',2016),('KEN',1998),('LKA',1998),('PLW',1998),('AUS',2016),('FJI',2000),('KIR',2015),('BLZ',1998)):
    if iso in set(Cc.iso):
        try: cs[f'{iso}_{y}']=case(iso,y)
        except Exception as ex: pass
X['cases']=cs
# 6. guc: n olay ve h=3-5 ortalama MDE (ev8)
X['isos']=sorted(Cc.iso.unique())
R['x']=X; json.dump(R,open('r/out/results.json','w'),ensure_ascii=False,indent=1,allow_nan=False)
for k,v in X.items(): print(k,v)
'''

# ==============================================================================
# r/r_extra2.py
# ==============================================================================
MODULES['r/r_extra2.py'] = r'''
import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore')
exec(open('r/r_analysis.py').read().split("R['lp_reef_dhw']")[0].replace("R={}","R={}").replace("C=pd.read_csv('r/out/country_catch.csv')","C=pd.read_csv('r/out/country_catch.csv')"))
R=json.load(open('r/out/results.json')); M=pd.read_csv('r/out/station_country.csv'); S=pd.read_csv('c/out/station_year.csv')[['sid','region']].drop_duplicates(); M=M.merge(S,on='sid',how='left')
reg=M.groupby('iso').region.agg(lambda s:s.mode().iloc[0] if s.notna().any() else None); print(reg.value_counts())
C['region']=C.iso.map(reg); atl=C.region.astype(str).str.contains('Atlantik'); H=[0,1,3,5,8,10]
R['lp_reef_ev4_atl']=lp('lreef','ev4',H=H,df=C[atl]); R['lp_reef_ev4_indopac']=lp('lreef','ev4',H=H,df=C[~atl]); R['n_atl']=int(C[atl].iso.nunique()); R['n_indopac']=int(C[~atl].iso.nunique())
R['lp_reef_ev8_indopac']=lp('lreef','ev8',H=H,df=C[~atl])
C['w']=C.groupby('iso').reef.transform('mean'); R['lp_reef_ev4_w']=lp('lreef','ev4',H=H,wt='w'); R['lp_reef_ev8_w']=lp('lreef','ev8',H=H,wt='w')
R['lp_reef_ev4_to2016']=lp('lreef','ev4',H=H,df=C[C.yr<=2016]); R['lp_reef_ev8_to2016']=lp('lreef','ev8',H=H,df=C[C.yr<=2016])
R['lp_reef_ev4_exIDN']=lp('lreef','ev4',H=H,df=C[C.iso!='IDN'])
# avin resif payi (pay olarak sonuc)
C['lsh']=np.log((C.reef/C.tot).where(C.reef>0)); R['lp_share_ev4']=lp('lsh','ev4',H=H); R['lp_share_ev8']=lp('lsh','ev8',H=H)
# ardisik olay: son 5 yilda baska siddetli olay yasanmis mi
g=C.groupby('iso').ev8; C['prev8']=sum(g.shift(k).fillna(0) for k in range(1,6)); C['first8']=((C.ev8==1)&(C.prev8==0)).astype(float); C['rep8']=((C.ev8==1)&(C.prev8>0)).astype(float)
R['lp_reef_first8']=lp('lreef','first8',H=H); R['lp_reef_rep8']=lp('lreef','rep8',H=H); R['n_first8']=int(C[C.yr>=1986].first8.sum()); R['n_rep8']=int(C[C.yr>=1986].rep8.sum())
R['sample']['n_st']=int(M[M.iso.isin(C.iso.unique())].sid.nunique()); R['sample']['n_st_all']=int(M.sid.nunique())
json.dump(R,open('r/out/results.json','w'),ensure_ascii=False,indent=1,allow_nan=False)
for k in('lp_reef_ev4_atl','lp_reef_ev4_indopac','lp_reef_ev8_indopac','lp_reef_ev4_w','lp_reef_ev8_w','lp_reef_ev4_to2016','lp_reef_ev8_to2016','lp_reef_ev4_exIDN','lp_share_ev4','lp_share_ev8','lp_reef_first8','lp_reef_rep8'): print(k,' '.join(f"{h}:{100*v['b']:.1f}({100*v['se']:.1f})" for h,v in R[k].items()))
print(R['n_atl'],R['n_indopac'],R['n_first8'],R['n_rep8'],R['sample']['n_st'])
'''

# ==============================================================================
# r/r_extra3.py
# ==============================================================================
MODULES['r/r_extra3.py'] = r'''
import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore')
exec(open('r/r_analysis.py').read().split("R['lp_reef_dhw']")[0])
R=json.load(open('r/out/results.json')); N={}
# --- (a) ilk / yinelenen olaylar
def mk(C,win,thr='ev8'):
    g=C.groupby('iso')[thr]; prev=sum(g.shift(k).fillna(0) for k in range(1,win+1)); return ((C[thr]==1)&(prev==0)).astype(float),((C[thr]==1)&(prev>0)).astype(float),prev
def lp2(dep,xs,H,df):   # iki olay degiskeni birlikte; katsayilar ve farkin kumeli se'si
    out={}
    for h in H:
        d=df.copy(); g=d.groupby('iso')[dep]; d['y']=g.shift(-h)-g.shift(1); d['g1']=g.shift(1)-g.shift(2); d['g2']=g.shift(2)-g.shift(3); cols=xs+['g1','g2']
        d=d[d.yr>=1986].dropna(subset=['y']+cols); d=d[np.isfinite(d.y)]
        X=pd.concat([d[cols],pd.get_dummies(d.iso,drop_first=True,dtype=float),pd.get_dummies(d.yr,drop_first=True,dtype=float)],axis=1); X.insert(0,'c',1.0); Xv=X.values; yv=d.y.values
        b=np.linalg.lstsq(Xv,yv,rcond=None)[0]; e=yv-Xv@b; XtXi=np.linalg.pinv(Xv.T@Xv); meat=np.zeros((Xv.shape[1],)*2)
        for iso,idx in d.reset_index(drop=True).groupby('iso').indices.items():
            s=Xv[idx].T@e[idx]; meat+=np.outer(s,s)
        G=d.iso.nunique(); V=XtXi@meat@XtXi*G/(G-1); k=len(xs); o={x:dict(b=float(b[1+i]),se=float(np.sqrt(V[1+i,1+i]))) for i,x in enumerate(xs)}
        if k==2: o['diff']=dict(b=float(b[2]-b[1]),se=float(np.sqrt(V[1,1]+V[2,2]-2*V[1,2])))
        out[h]=o
    return out
H=[0,1,2,3,4,5,6,8,10]
for win in(3,5,8):
    C['first'],C['rep'],C['prev']=mk(C,win); o=lp2('lreef',['first','rep'],H,C); N[f'fr_w{win}']={str(h):v for h,v in o.items()}; N[f'n_w{win}']=[int(C[C.yr>=1986]['first'].sum()),int(C[C.yr>=1986].rep.sum())]
C['first'],C['rep'],C['prev']=mk(C,5)
def avg(o,k,hs=(3,4,5)): return float(np.mean([o[h][k]['b'] for h in hs]))
o=lp2('lreef',['first','rep'],[3,4,5],C); obs={k:avg(o,k) for k in('first','rep','diff')}
rng=np.random.default_rng(21); nul=[]
for it in range(500):
    d=C.copy(); d['ev8']=d.groupby('iso').ev8.transform(lambda s:np.roll(s.values,rng.integers(1,len(s)))); d['first'],d['rep'],_=mk(d,5); oo=lp2('lreef',['first','rep'],[3,4,5],d); nul.append([avg(oo,k) for k in('first','rep','diff')])
nul=np.array(nul); N['fr_perm']={k:dict(obs=obs[k],p_two=float((np.abs(nul[:,i])>=abs(obs[k])).mean()),p_low=float((nul[:,i]<=obs[k]).mean()),sd=float(nul[:,i].std())) for i,k in enumerate(('first','rep','diff'))}
# olaylar arasi sure (surekli): onceki siddetli olaydan bu yana gecen yil
def gap(s):
    out=[];last=None
    for i,v in enumerate(s.values):
        out.append(np.nan if last is None else i-last)
        if v==1: last=i
    return pd.Series(out,index=s.index)
C['gapyr']=C.groupby('iso').ev8.transform(gap); E=C[(C.ev8==1)&(C.yr>=1986)].copy(); N['gap_desc']=dict(n=int(len(E)),n_with_prev=int(E.gapyr.notna().sum()),med=float(E.gapyr.median()),share_le5=float((E.gapyr<=5).mean()),by_dec={str(k):[int(len(g)),float((g.gapyr<=5).mean())] for k,g in E.groupby((E.yr//10)*10)})
# plasebo: orkinos ve toplam av icin ayni ayrim
for dep in('ltuna','ltot','lcoast'): N['fr_'+dep]={str(h):v for h,v in lp2(dep,['first','rep'],[1,3,5,8],C).items()}
# Hint-Pasifik
M=pd.read_csv('r/out/station_country.csv').merge(pd.read_csv('c/out/station_year.csv')[['sid','region']].drop_duplicates(),on='sid',how='left'); reg=M.groupby('iso').region.agg(lambda s:s.mode().iloc[0] if s.notna().any() else None); atl=C.iso.map(reg).astype(str).str.contains('Atlantik')
N['fr_indopac']={str(h):v for h,v in lp2('lreef',['first','rep'],[1,3,5,8,10],C[~atl]).items()}; N['fr_atl']={str(h):v for h,v in lp2('lreef',['first','rep'],[1,3,5,8,10],C[atl]).items()}
# 4 C-hafta esigiyle
C['first4'],C['rep4'],_=mk(C,5,'ev4'); N['fr4']={str(h):v for h,v in lp2('lreef',['first4','rep4'],[1,3,5,8,10],C).items()}
# --- (c) degiskenlik: |dlog| sonuc
C['absd']=C.groupby('iso').lreef.diff().abs()
def lpv(x,H=[1,2,3,4,5]):
    out={}
    for h in H:
        d=C.copy(); g=d.groupby('iso').absd; d['y']=sum(g.shift(-k) for k in range(1,h+1))/h-(g.shift(1)+g.shift(2)+g.shift(3))/3; d=d[d.yr>=1986].dropna(subset=['y',x])
        X=pd.concat([d[[x]],pd.get_dummies(d.iso,drop_first=True,dtype=float),pd.get_dummies(d.yr,drop_first=True,dtype=float)],axis=1); X.insert(0,'c',1.0); Xv=X.values; yv=d.y.values; b=np.linalg.lstsq(Xv,yv,rcond=None)[0]; e=yv-Xv@b; XtXi=np.linalg.pinv(Xv.T@Xv); meat=np.zeros((Xv.shape[1],)*2)
        for iso,idx in d.reset_index(drop=True).groupby('iso').indices.items():
            s=Xv[idx].T@e[idx]; meat+=np.outer(s,s)
        G=d.iso.nunique(); V=XtXi@meat@XtXi*G/(G-1); out[str(h)]=dict(b=float(b[1]),se=float(np.sqrt(V[1,1])))
    return out
N['vol_ev8']=lpv('ev8'); N['vol_ev4']=lpv('ev4'); N['vol_base']=float(C.absd.median())
R['n']=N; json.dump(R,open('r/out/results.json','w'),ensure_ascii=False,indent=1,allow_nan=False); C.to_csv('r/out/panel.csv',index=False)
for w in(3,5,8): print('win',w,N[f'n_w{w}']); [print('  ',h,{k:(round(100*v['b'],1),round(100*v['se'],1)) for k,v in o.items()}) for h,o in N[f'fr_w{w}'].items()]
print(N['fr_perm']); print(N['gap_desc'])
for k in('fr_ltuna','fr_ltot','fr_indopac','fr_atl','fr4'): print(k,{h:{a:(round(100*v['b'],1),round(100*v['se'],1)) for a,v in o.items()} for h,o in N[k].items()})
print('vol',N['vol_ev8'],N['vol_ev4'],N['vol_base'])
'''

# ==============================================================================
# r/r_sau_prep.py
# ==============================================================================
MODULES['r/r_sau_prep.py'] = r'''
import os, urllib.request
os.makedirs('data/sau', exist_ok=True)
if not os.path.exists('data/sau/eez.json'):
    with urllib.request.urlopen(urllib.request.Request('https://api.seaaroundus.org/api/v1/eez/?geojson=true', headers={'User-Agent': 'gomd-replication/1.0'}), timeout=120) as _r, open('data/sau/eez.json', 'wb') as _f: _f.write(_r.read())
# Sea Around Us yeniden yapilandirilmis av verisi: CRW istasyonlarinin MEB'lere atanmasi ve resifle iliskili balik avi
import json, sys, time, urllib.request, os, pandas as pd, numpy as np
from shapely.geometry import shape, Point
E=json.load(open('data/sau/eez.json'))['features']; G=[(f['properties']['region_id'],f['properties']['title'],shape(f['geometry']).buffer(0)) for f in E]
s=pd.read_csv('c/out/stations.csv'); rows=[]
for r in s.itertuples():
    p=Point(r.lon,r.lat); hit=[(i,t) for i,t,g in G if g.contains(p)]
    if hit: i,t=hit[0]; dist=0.0
    else:
        dd=[(g.distance(p),i,t) for i,t,g in G]; dist,i,t=min(dd)
    rows.append((r.sid,r.name,i,t,round(dist,2)))
M=pd.DataFrame(rows,columns=['sid','name','eez','title','deg']); M=M[M.deg<=1.5]; M.to_csv('r/out/station_eez.csv',index=False); print(len(M),M.eez.nunique()); print(M[M.deg>0].to_string())
os.makedirs('data/sau/fg',exist_ok=True)
def get(dim,i):
    fn=f'data/sau/fg/{dim}_{i}.json'
    if os.path.exists(fn): return json.load(open(fn))
    for k in range(4):
        try:
            with urllib.request.urlopen(urllib.request.Request(f'https://api.seaaroundus.org/api/v1/eez/tonnage/{dim}/?format=json&limit=50&sciname=false&region_id={i}',headers={'User-Agent':'research/1.0'}),timeout=60) as r: d=json.load(r); json.dump(d,open(fn,'w')); time.sleep(.4); return d
        except Exception as e: time.sleep(3*(k+1))
    return None
out=[]
for i,t in M[['eez','title']].drop_duplicates().itertuples(index=False):
    for dim in('functionalgroup','sector','reporting-status'):
        d=get(dim,i)
        if not d or not d.get('data'): print('FAIL',dim,i,t); continue
        for x in d['data']:
            for y,v in x['values']: out.append((i,t,dim,x['key'],y,v))
D=pd.DataFrame(out,columns=['eez','title','dim','key','yr','val']); D.to_parquet('r/out/sau_raw.parquet'); print(D.groupby('dim').key.unique().to_dict()); print(D.eez.nunique(),D.yr.min(),D.yr.max())
'''

# ==============================================================================
# r/r_sau.py
# ==============================================================================
MODULES['r/r_sau.py'] = r'''
# Sea Around Us yeniden yapilandirilmis resif baligi avi ile ayni sinama (MEB duzeyi, 1986-2019)
import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore')
D=pd.read_parquet('r/out/sau_raw.parquet'); M=pd.read_csv('r/out/station_eez.csv')
c=pd.read_parquet('c/out/crw.parquet',columns=['sid','date','dhw']); c['yr']=pd.to_datetime(c.date).dt.year; sy=c[(c.yr>=1985)&(c.yr<=2019)].groupby(['sid','yr']).dhw.max().rename('mx').reset_index()
H=M.merge(sy,on='sid').groupby(['eez','yr']).agg(dhw=('mx','mean'),dhw_max=('mx','max')).reset_index()
fg=D[D.dim=='functionalgroup']; reef=fg[fg.key.str.contains('reef assoc')].groupby(['eez','title','yr']).val.sum().rename('reef'); tot=fg.groupby(['eez','yr']).val.sum().rename('tot'); lp_=fg[fg.key.str.contains('Large pelagics')].groupby(['eez','yr']).val.sum().rename('tuna')
sec=D[D.dim=='sector'].pivot_table(index=['eez','yr'],columns='key',values='val',aggfunc='sum').fillna(0); sec['ssf']=sec.get('Artisanal',0)+sec.get('Subsistence',0)
rs=D[D.dim=='reporting-status'].pivot_table(index=['eez','yr'],columns='key',values='val',aggfunc='sum').fillna(0); print(rs.columns.tolist())
C=reef.reset_index().merge(tot.reset_index(),on=['eez','yr']).merge(lp_.reset_index(),on=['eez','yr'],how='left').merge(sec[['ssf']].reset_index(),on=['eez','yr'],how='left').merge(rs.reset_index(),on=['eez','yr'],how='left')
C=C[(C.yr>=1980)&(C.yr<=2019)]; ok=C.groupby('eez').apply(lambda g:(g.reef>0).sum()>=35 and g.reef.mean()>=20); C=C[C.eez.isin(ok[ok].index)].merge(H,on=['eez','yr'],how='left').sort_values(['eez','yr']).reset_index(drop=True); C['iso']=C.eez.astype(str)
for g in('reef','tot','tuna','ssf'): C['l'+g]=np.log(C[g].where(C[g]>0))
C['ev4']=(C.dhw>=4).astype(float); C['ev8']=(C.dhw>=8).astype(float)
src=open('r/r_analysis.py').read(); a=src.index('def lp('); b=src.index("C['dhw4']=C.dhw/4"); exec(src[a:b])
C['dhw4']=C.dhw/4; R={}; Hh=[0,1,2,3,4,5,6,8,10]
R['sample']=dict(n=int(C.eez.nunique()),n_obs=int(len(C)),n_ev4=int(C[C.yr>=1986].ev4.sum()),n_ev8=int(C[C.yr>=1986].ev8.sum()),reef_2010_19=float(C[C.yr>=2010].groupby('yr').reef.sum().mean()),reef_share=float(C[C.yr>=2010].reef.sum()/C[C.yr>=2010].tot.sum()),sd_dlog=float(C.groupby('eez').lreef.diff().groupby(C.eez).std().median()),ac1=float(C.assign(g=C.groupby('eez').lreef.diff()).groupby('eez').g.apply(lambda s:s.autocorr()).median()),
   unrep=float(C[C.yr>=2010].get('Unreported',pd.Series([0])).sum()/C[C.yr>=2010].tot.sum()) if 'Unreported' in C else None, ssf_share=float(C[C.yr>=2010].ssf.sum()/C[C.yr>=2010].tot.sum()))
for k,(dep,x) in {'reef_ev4':('lreef','ev4'),'reef_ev8':('lreef','ev8'),'reef_dhw':('lreef','dhw4'),'tuna_ev4':('ltuna','ev4'),'tot_ev4':('ltot','ev4'),'ssf_ev4':('lssf','ev4'),'ssf_ev8':('lssf','ev8')}.items(): R[k]={str(h):v for h,v in lp(dep,x,H=Hh,df=C).items()}
R['reef_ev4_pre']={str(h):v for h,v in lp('lreef','ev4',H=[-4,-3,-2],df=C,ctrl=False).items()}
g=C.groupby('iso').ev8; prev=sum(g.shift(k).fillna(0) for k in range(1,6)); C['first8']=((C.ev8==1)&(prev==0)).astype(float); C['rep8']=((C.ev8==1)&(prev>0)).astype(float)
R['reef_first8']={str(h):v for h,v in lp('lreef','first8',H=Hh,df=C).items()}; R['reef_rep8']={str(h):v for h,v in lp('lreef','rep8',H=Hh,df=C).items()}; R['n_first8']=int(C[C.yr>=1986].first8.sum()); R['n_rep8']=int(C[C.yr>=1986].rep8.sum())
big=C.groupby('eez').apply(lambda g:g.reef.sum()/g.tot.sum()); hid=big[big>=big.median()].index; R['reef_ev8_dep']={str(h):v for h,v in lp('lreef','ev8',H=Hh,df=C[C.eez.isin(hid)]).items()}
R['mde_ev8']={str(h):float(2.8*R['reef_ev8'][str(h)]['se']) for h in(0,3,5,8)}; R['mde_ev4']={str(h):float(2.8*R['reef_ev4'][str(h)]['se']) for h in(0,3,5,8)}
C.to_csv('r/out/sau_panel.csv',index=False); J=json.load(open('r/out/results.json')); J['sau']=R; json.dump(J,open('r/out/results.json','w'),ensure_ascii=False,indent=1,allow_nan=False)
print(R['sample'])
for k in R:
    if isinstance(R[k],dict) and isinstance(list(R[k].values())[0],dict): print(k,' '.join(f"{h}:{100*v['b']:.1f}({100*v['se']:.1f})" for h,v in R[k].items()))
print(R['mde_ev8'],R['mde_ev4'],R['n_first8'],R['n_rep8'])
'''

# ==============================================================================
# r/r_sau2.py
# ==============================================================================
MODULES['r/r_sau2.py'] = r'''
import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore')
C=pd.read_csv('r/out/sau_panel.csv'); C['iso']=C.eez.astype(str); M=pd.read_csv('r/out/station_eez.csv').merge(pd.read_csv('c/out/station_year.csv')[['sid','region']].drop_duplicates(),on='sid',how='left'); reg=M.groupby('eez').region.agg(lambda s:s.mode().iloc[0]); C['region']=C.eez.map(reg)
C['nonreef']=C.tot-C.reef; C['lnon']=np.log(C.nonreef.where(C.nonreef>0)); C['lsh']=C.lreef-C.ltot
def lp(dep,x,H,df,fe='yr',ctrl=True,extra=[]):
    out={}
    for h in H:
        d=df.copy(); g=d.groupby('iso')[dep]; d['y']=g.shift(-h)-g.shift(1); d['g1']=g.shift(1)-g.shift(2); d['g2']=g.shift(2)-g.shift(3); cols=[x]+extra+(['g1','g2'] if ctrl else [])
        d=d[d.yr>=1986].dropna(subset=['y']+cols); d=d[np.isfinite(d.y)]
        tfe=pd.get_dummies(d.yr,drop_first=True,dtype=float) if fe=='yr' else pd.get_dummies(d.region+'_'+d.yr.astype(str),drop_first=True,dtype=float)
        X=pd.concat([d[cols],pd.get_dummies(d.iso,drop_first=True,dtype=float),tfe],axis=1); X.insert(0,'c',1.0); Xv=X.values; yv=d.y.values; b=np.linalg.lstsq(Xv,yv,rcond=None)[0]; e=yv-Xv@b; XtXi=np.linalg.pinv(Xv.T@Xv); meat=np.zeros((Xv.shape[1],)*2)
        for iso,idx in d.reset_index(drop=True).groupby('iso').indices.items():
            s=Xv[idx].T@e[idx]; meat+=np.outer(s,s)
        G=d.iso.nunique(); V=XtXi@meat@XtXi*G/(G-1); out[str(h)]=dict(b=float(b[1]),se=float(np.sqrt(V[1,1])),n=int(len(d)),G=int(G))
    return out
H=[0,1,2,3,4,5,6,8,10]; S={}
S['non_ev4']=lp('lnon','ev4',H,C); S['share_ev4']=lp('lsh','ev4',H,C); S['share_ev8']=lp('lsh','ev8',H,C); S['reef_ev4_ry']=lp('lreef','ev4',H,C,fe='ry'); S['reef_ev8_ry']=lp('lreef','ev8',H,C,fe='ry'); S['share_ev4_ry']=lp('lsh','ev4',H,C,fe='ry'); S['non_ev4_ry']=lp('lnon','ev4',H,C,fe='ry')
C['dnon']=C.groupby('iso').lnon.diff(); 
exC=C[~C.title.str.contains('China|Viet')]; S['reef_ev4_exCN']=lp('lreef','ev4',H,exC); S['n_exCN']=int(exC.eez.nunique())
dep=C.groupby('eez').apply(lambda g:g.reef.sum()/g.tot.sum()); hi=dep[dep>=dep.median()].index; S['reef_ev4_dep']=lp('lreef','ev4',H,C[C.eez.isin(hi)]); S['share_ev4_dep']=lp('lsh','ev4',H,C[C.eez.isin(hi)])
atl=C.region.astype(str).str.contains('Atlantik'); S['reef_ev4_indopac']=lp('lreef','ev4',H,C[~atl]); S['reef_ev4_atl']=lp('lreef','ev4',H,C[atl]); S['n_atl']=int(C[atl].eez.nunique())
S['reef_ev4_pre_ry']=lp('lreef','ev4',[-4,-3,-2],C,fe='ry',ctrl=False); S['share_ev4_pre']=lp('lsh','ev4',[-4,-3,-2],C,ctrl=False)
# permutasyon: h=3..5 ortalamasi, resif ve pay
def avg(df,dep,fe='yr'): o=lp(dep,'ev4',[3,4,5],df,fe=fe); return np.mean([o[str(h)]['b'] for h in(3,4,5)])
rng=np.random.default_rng(9); obs={'reef':avg(C,'lreef'),'share':avg(C,'lsh'),'reef_ry':avg(C,'lreef','ry')}; nul={k:[] for k in obs}
for it in range(300):
    d=C.copy(); d['ev4']=d.groupby('iso').ev4.transform(lambda s:np.roll(s.values,rng.integers(1,len(s)))); nul['reef'].append(avg(d,'lreef')); nul['share'].append(avg(d,'lsh')); nul['reef_ry'].append(avg(d,'lreef','ry'))
S['perm']={k:dict(obs=float(obs[k]),p=float((np.abs(nul[k])>=abs(obs[k])).mean()),sd=float(np.std(nul[k])),null_mean=float(np.mean(nul[k]))) for k in obs}
J=json.load(open('r/out/results.json')); J['sau2']=S; json.dump(J,open('r/out/results.json','w'),ensure_ascii=False,indent=1,allow_nan=False)
for k,v in S.items():
    if isinstance(v,dict) and 'obs' not in str(list(v.values())[0]): print(k,' '.join(f"{h}:{100*o['b']:.1f}({100*o['se']:.1f})" for h,o in v.items()))
print(S['perm'],S['n_exCN'],S['n_atl'])
'''

# ==============================================================================
# r/r_fig.py
# ==============================================================================
MODULES['r/r_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
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('r/out/results.json')); X=R['x']; C=pd.read_csv('r/out/panel.csv'); B='#1b4f72'; O='#b9770e'; Rd='#b03a2e'; Gn='#1e8449'; Gy='#909497'
def irf(a,key,c,lab,H=range(0,11),off=0,fill=True,ls='-'):
    o=R[key]; h=[k for k in H if str(k) in o]; b=np.array([100*o[str(k)]['b'] for k in h]); s=np.array([100*o[str(k)]['se'] for k in h]); x=np.array(h)+off
    a.plot(x,b,ls,color=c,lw=1.4,marker='o',ms=3,label=lab)
    if fill: a.fill_between(x,b-1.96*s,b+1.96*s,color=c,alpha=.13,lw=0)
# S1
W=C.groupby('yr').agg(reef=('reef','sum'),tot=('tot','sum')).reset_index(); E=C[C.yr>=1986].groupby('yr').agg(e4=('ev4','mean'),e8=('ev8','mean')).reset_index()
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1))
a=ax[0]; a.plot(W.yr,W.reef/1e6,color=B,lw=1.5); a.set_ylabel('Resif balığı avı (milyon ton)',color=B); a2=a.twinx(); a2.plot(W.yr,100*W.reef/W.tot,color=O,lw=1.2,ls='--'); a2.set_ylabel('Toplam deniz avındaki pay (%)',color=O); a2.spines['top'].set_visible(False); a.set_title('(a) Bildirilen resif balığı avı',fontsize=9,loc='left')
a=ax[1]; a.bar(E.yr,100*E.e4,color='#f5b7b1',width=.85,label='≥ 4 °C-hafta'); a.bar(E.yr,100*E.e8,color=Rd,width=.85,label='≥ 8 °C-hafta'); a.set_ylabel('Isı stresi yaşayan ülkelerin payı (%)'); a.legend(frameon=False,fontsize=7.5,loc='upper left'); a.set_title('(b) Isı stresinin yaygınlığı',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('r/fig/Figure_1.png'); plt.close()
# S2
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1),sharey=True)
a=ax[0]; irf(a,'lp_reef_ev4',B,'≥ 4 °C-hafta'); o=R['lp_reef_ev4_noctrl']; hh=[-4,-3,-2]; a.errorbar(hh,[100*o[str(k)]['b'] for k in hh],yerr=[196*o[str(k)]['se'] for k in hh],fmt='s',color=Gy,ms=3.5,lw=.8,label='Olay öncesi (plasebo)'); a.axhline(0,color='#999',lw=.6); a.axvline(-.5,color='#999',lw=.5,ls=':')
a.set_xlabel('Isı stresi yılından sonraki yıl'); a.set_ylabel('Resif balığı avındaki değişim (log × 100)'); a.legend(frameon=False,fontsize=7.5,loc='lower left'); a.set_title('(a) Orta düzey ısı stresi',fontsize=9,loc='left')
a=ax[1]; irf(a,'lp_reef_ev8',Rd,'≥ 8 °C-hafta (ülke ortalaması)'); irf(a,'lp_reef_evm8',O,'≥ 8 °C-hafta (en az bir istasyon)',fill=False,ls='--',off=.12); a.axhline(0,color='#999',lw=.6); a.set_xlabel('Isı stresi yılından sonraki yıl'); a.legend(frameon=False,fontsize=7.5,loc='lower left'); a.set_title('(b) Şiddetli ısı stresi',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('r/fig/Figure_2.png'); plt.close()
# S3
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1),sharey=True)
a=ax[0]; irf(a,'lp_reef_ev4',B,'Resif balıkları'); irf(a,'lp_tuna_ev4',Rd,'Orkinos ve benzerleri',fill=False,off=.1); irf(a,'lp_pel_ev4',Gn,'Küçük pelajikler',fill=False,off=-.1,ls='--'); irf(a,'lp_tot_ev4',Gy,'Toplam deniz avı',fill=False,ls=':'); a.axhline(0,color='#999',lw=.6)
a.set_xlabel('Isı stresi yılından sonraki yıl'); a.set_ylabel('Avdaki değişim (log × 100)'); a.legend(frameon=False,fontsize=7,loc='lower left'); a.set_title('(a) Tür gruplarına göre',fontsize=9,loc='left')
a=ax[1]; irf(a,'lp_reef_ev4',B,'Tüm ülkeler',fill=False); irf(a,'lp_reef_ev4_official',O,'Resmî veri ağırlıklı',fill=False,off=.1); irf(a,'lp_reef_ev4_dep',Rd,'Resif bağımlılığı yüksek',fill=True,off=-.1); irf(a,'lp_reef_ev4_big',Gn,'Büyük üreticiler',fill=False,ls='--'); a.axhline(0,color='#999',lw=.6)
a.set_xlabel('Isı stresi yılından sonraki yıl'); a.legend(frameon=False,fontsize=7,loc='lower left'); a.set_title('(b) Ülke gruplarına göre (resif balıkları)',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('r/fig/Figure_3.png'); plt.close()
# S4
L=pd.read_csv('r/out/longdiff.csv',index_col=0)
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1))
a=ax[0]; ks=['0-1','1-2','2-4','4-8','8+']; a.bar(range(5),[100*R['dose'][k]['mean'] for k in ks],yerr=[196*R['dose'][k]['se'] for k in ks],color=[B,B,B,O,Rd],width=.6,error_kw=dict(lw=.8)); a.axhline(0,color='#999',lw=.6)
a.set_xticks(range(5)); a.set_xticklabels(['0–1','1–2','2–4','4–8','8+']); a.set_xlabel('Yılın en yüksek ısı stresi (°C-hafta)'); a.set_ylabel('Beş yıllık değişim, yıl ortalamasından sapma (log × 100)'); a.yaxis.label.set_size(7.5); a.set_title('(a) Doz–tepki',fontsize=9,loc='left')
a=ax[1]; a.scatter(L.n8,100*L.dl,s=8+np.sqrt(L.r1)/6,color=B,alpha=.6,lw=0); xx=np.array([0,L.n8.max()]); m=X['ld_dl_n8']; ic=100*(L.dl.mean()-m['b']*L.n8.mean()); a.plot(xx,ic+100*m['b']*xx,color=Rd,lw=1.1); a.axhline(0,color='#999',lw=.6)
for iso in('IDN','IND','PHL','FJI','SYC','EGY'):
    if iso in L.index: a.annotate(iso,(L.n8[iso],100*L.dl[iso]),fontsize=6.5,xytext=(3,3),textcoords='offset points')
a.set_xlabel('Şiddetli ısı stresi yılı sayısı, 1998–2019'); a.set_ylabel('Resif avı değişimi, 1990–99 → 2015–24 (log × 100)'); a.yaxis.label.set_size(7.5); a.set_title('(b) Uzun dönem farkı',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('r/fig/Figure_4.png'); plt.close()
# S5
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1))
a=ax[0]
for y,c in((1998,B),(2010,O),(2016,Rd)): irf(a,f'lp_reef_{y}',c,str(y),H=range(0,9),fill=False,off=(y-2010)/90)
a.axhline(0,color='#999',lw=.6); a.set_xlabel('Olaydan sonraki yıl'); a.set_ylabel('Resif balığı avındaki değişim (log × 100)'); a.legend(frameon=False,fontsize=7.5); a.set_title('(a) Üç küresel ağarma olayı',fontsize=9,loc='left')
a=ax[1]; ev=pd.Series({int(h):v for h,v in R['mde_lp_reef_ev4'].items()}); e8=pd.Series({int(h):v for h,v in R['mde_lp_reef_ev8'].items()}); a.plot(ev.index,100*ev.values,'o-',color=B,label='≥ 4 °C-hafta'); a.plot(e8.index,100*e8.values,'s-',color=Rd,label='≥ 8 °C-hafta'); a.set_ylim(0,25)
a.set_xlabel('Isı stresi yılından sonraki yıl'); a.set_ylabel('En küçük saptanabilir etki (log × 100)'); a.legend(frameon=False,fontsize=7.5); a.set_title('(b) İstatistiksel güç',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('r/fig/Figure_5.png'); plt.close()
# S6 Sea Around Us
R=json.load(open('r/out/results.json')); SA=R['sau']; S2=R['sau2']
def irf2(a,o,c,lab,off=0,fill=True,ls='-'):
    h=sorted(int(k) for k in o if int(k)>=0); b=np.array([100*o[str(k)]['b'] for k in h]); s=np.array([100*o[str(k)]['se'] for k in h]); x=np.array(h)+off; a.plot(x,b,ls,color=c,lw=1.4,marker='o',ms=3,label=lab)
    if fill: a.fill_between(x,b-1.96*s,b+1.96*s,color=c,alpha=.13,lw=0)
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1),sharey=True)
a=ax[0]; irf2(a,SA['reef_ev4'],B,'≥ 4 °C-hafta'); irf2(a,SA['reef_ev8'],Rd,'≥ 8 °C-hafta',off=.12); o=SA['reef_ev4_pre']; hh=[-4,-3,-2]; a.errorbar(hh,[100*o[str(k)]['b'] for k in hh],yerr=[196*o[str(k)]['se'] for k in hh],fmt='s',color=Gy,ms=3.5,lw=.8,label='Olay öncesi (plasebo)'); a.axhline(0,color='#999',lw=.6); a.axvline(-.5,color='#999',lw=.5,ls=':')
a.set_xlabel('Isı stresi yılından sonraki yıl'); a.set_ylabel('Avdaki değişim (log × 100)'); a.legend(frameon=False,fontsize=7.5,loc='upper left'); a.set_title('(a) Yeniden yapılandırılmış resif balığı avı',fontsize=9,loc='left')
a=ax[1]; irf2(a,SA['reef_ev4'],B,'Resifle ilişkili balıklar',fill=False); irf2(a,S2['non_ev4'],Gy,'Öteki türler',fill=False,ls='--',off=.1); irf2(a,S2['share_ev4'],O,'Resif balıklarının payı',off=-.1); a.axhline(0,color='#999',lw=.6); a.set_xlabel('Isı stresi yılından sonraki yıl'); a.legend(frameon=False,fontsize=7.5,loc='upper left'); a.set_title('(b) Resife özgü mü? (≥ 4 °C-hafta)',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('r/fig/Figure_6.png'); plt.close()
'''

if __name__ == "__main__":
    for _d in ('r', 'r/out', 'r/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: r/out/results.json')