# -*- coding: utf-8 -*-
# ------------------------------------------------------------------------------
# Büyük gemiler limanlarda yük yığılmasını ne zaman artırır? 1.264 limanda bileşik Poisson öngörüsünün sınanması ve uğrak başına yükteki artışın ayrıştırılması, 2019–2026
# When do large ships increase cargo bunching at ports? A test of the compound Poisson prediction at 1,264 ports and a decomposition of the rise in cargo per call, 2019–2026
#
# GOMD Yayınlar · Makale 25 — https://gomd.pirireis.edu.tr/publications.php?id=rp25
#
# Bu dosya makalenin bütün modüllerini içerir. Replikasyon paketi klasöründe çalıştırıldığında modülleri `s/` klasörüne yazar ve sırayla, her birini ayrı süreçte çalıştırır;
# veriyi `data/` ve `s/out/` klasörlerinden okur, sonuçları `s/out/` ve şekilleri `s/fig/` klasörüne yazar.
# This file contains all modules of the paper. Run in the replication-package folder, it writes the modules to `s/` and runs them in order, each in its own process;
# it reads data from `data/` and `s/out/`, and writes results to `s/out/` and figures to `s/fig/`.
#
# Kurulum / Setup:  pip install pandas numpy matplotlib scipy statsmodels pyarrow
# Çalıştırma / Run:  python rp25_large_ships_cargo_bunching_compound_poisson.py
#
# Port data: IMF PortWatch / UN Global Platform / University of Oxford (https://portwatch.imf.org), daily port activity and trade estimates derived from AIS; retrieved from the official open ArcGIS API on 13 November 2026.
# 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 = ['s/s_analysis.py', 's/s_extra.py', 's/s_extra2.py', 's/s_disp.py', 's/s_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

# ==============================================================================
# s/s_analysis.py
# ==============================================================================
MODULES['s/s_analysis.py'] = r'''
# Daha buyuk gemi, daha az ugrak mi? Liman yuk hacmi = ugrak sayisi x ugrak basina yuk (PortWatch gunluk)
import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore')
import statsmodels.formula.api as smf
d=pd.read_parquet('data/pwd/daily.parquet'); P=pd.read_csv('data/pw/ports_global.csv').set_index('portid'); wb=pd.read_csv('data/wb_country.csv').set_index('iso3')
d['y']=d.date.dt.year; d['thr']=d['import']+d['export']; d['thr_c']=d.import_container+d.export_container; d['thr_n']=d.thr-d.thr_c; d['calls_n']=d.portcalls-d.portcalls_container
R={}
Y=list(range(2019,2026)); full=d[d.y<=2025]
A=full.groupby(['portid','y'])[['portcalls','thr','portcalls_container','thr_c','calls_n','thr_n','import','export']].sum().reset_index()
# 1. dunya toplami
W=A.groupby('y').sum(numeric_only=True); nd=full.groupby('y').date.nunique()
def tab(c,t):
    o={}
    for y in Y: o[y]=dict(calls=float(W.loc[y,c]/nd[y]),thr=float(W.loc[y,t]/nd[y]),q=float(W.loc[y,t]/W.loc[y,c]))
    o['dln']={k:float(np.log(o[2025][k]/o[2019][k])) for k in('calls','thr','q')}
    o['dln2022']={k:float(np.log(o[2025][k]/o[2022][k])) for k in('calls','thr','q')}
    return o
R['world']={'all':tab('portcalls','thr'),'cont':tab('portcalls_container','thr_c'),'non':tab('calls_n','thr_n')}
# ilk 10 ay karsilastirmasi 2019 - 2026
m10=d[d.date.dt.month<=10].groupby('y')[['portcalls','thr','portcalls_container','thr_c']].sum()
R['jan_oct']={str(y):dict(q=float(m10.loc[y,'thr']/m10.loc[y,'portcalls']),qc=float(m10.loc[y,'thr_c']/m10.loc[y,'portcalls_container']),calls=float(m10.loc[y,'portcalls']),thr=float(m10.loc[y,'thr'])) for y in m10.index}
# 2. ugrak basina yukteki degisimin kaydirma-pay ayristirmasi (liman ici / limanlar arasi), 2019->2025
def shiftshare(c,t,y0=2019,y1=2025):
    a=A[A.y==y0].set_index('portid'); b=A[A.y==y1].set_index('portid'); ids=a.index.union(b.index); a=a.reindex(ids).fillna(0); b=b.reindex(ids).fillna(0)
    s0=a[c]/a[c].sum(); s1=b[c]/b[c].sum(); q0=(a[t]/a[c]).replace([np.inf,-np.inf],np.nan); q1=(b[t]/b[c]).replace([np.inf,-np.inf],np.nan)
    Q0=a[t].sum()/a[c].sum(); Q1=b[t].sum()/b[c].sum(); both=q0.notna()&q1.notna()
    within=((s0+s1)/2*(q1-q0))[both].sum(); between=(((q0+q1)/2-(Q0+Q1)/2)*(s1-s0))[both].sum(); entry=(s1*(q1-(Q0+Q1)/2))[q0.isna()&q1.notna()].sum()-(s0*(q0-(Q0+Q1)/2))[q0.notna()&q1.isna()].sum()
    return dict(Q0=float(Q0),Q1=float(Q1),d=float(Q1-Q0),within=float(within),between=float(between),entry=float(entry),within_share=float(within/(Q1-Q0)),between_share=float(between/(Q1-Q0)))
R['ss']={'all':shiftshare('portcalls','thr'),'cont':shiftshare('portcalls_container','thr_c'),'non':shiftshare('calls_n','thr_n')}
# gemi tipi bilesimi: q = sc*qc + (1-sc)*qn
w=R['world']; sc0=W.loc[2019,'portcalls_container']/W.loc[2019,'portcalls']; sc1=W.loc[2025,'portcalls_container']/W.loc[2025,'portcalls']
qc0,qc1,qn0,qn1=w['cont'][2019]['q'],w['cont'][2025]['q'],w['non'][2019]['q'],w['non'][2025]['q']
R['mix']=dict(sc0=float(sc0),sc1=float(sc1),within_c=float((sc0+sc1)/2*(qc1-qc0)),within_n=float((1-(sc0+sc1)/2)*(qn1-qn0)),mix=float((sc1-sc0)*((qc0+qc1)/2-(qn0+qn1)/2)),d=float(w['all'][2025]['q']-w['all'][2019]['q']))
# 3. liman duzeyi: 2019-20 ortalamasi -> 2024-25 ortalamasi
def portchg(c,t,minc):
    a=A[A.y.isin([2019,2020])].groupby('portid')[[c,t]].sum(); b=A[A.y.isin([2024,2025])].groupby('portid')[[c,t]].sum(); ok=(a[c]>=minc)&(b[c]>=minc)&(a[t]>0)&(b[t]>0); a=a[ok]; b=b[ok]
    D=pd.DataFrame({'dT':np.log(b[t]/a[t]),'dC':np.log(b[c]/a[c])}); D['dq']=D.dT-D.dC; D['lnT0']=np.log(a[t]); D['lnq0']=np.log(a[t]/a[c]); D['w']=(a[t]+b[t])/2; D['calls0']=a[c]/2
    D['cont']=P.continent.reindex(D.index); D['iso']=P.ISO3.reindex(D.index); D['inc']=wb.income.reindex(D.iso.values).values
    return D
out={}
for nm,(c,t,mc) in {'all':('portcalls','thr',200),'cont':('portcalls_container','thr_c',100),'non':('calls_n','thr_n',200)}.items():
    D=portchg(c,t,mc); D.to_csv(f's/out/port_{nm}.csv'); v=D.dT.var()
    o=dict(n=int(len(D)),med={k:float(D[k].median()) for k in('dT','dC','dq')},mean={k:float(D[k].mean()) for k in('dT','dC','dq')},wmean={k:float(np.average(D[k],weights=D.w)) for k in('dT','dC','dq')},
        share_q_up=float((D.dq>0).mean()),share_c_down=float((D.dC<0).mean()),share_fewer_bigger=float(((D.dC<0)&(D.dq>0)).mean()),share_more_smaller=float(((D.dC>0)&(D.dq<0)).mean()),share_both_up=float(((D.dC>0)&(D.dq>0)).mean()),share_both_down=float(((D.dC<0)&(D.dq<0)).mean()),
        var_share_calls=float(np.cov(D.dT,D.dC)[0,1]/v),var_share_q=float(np.cov(D.dT,D.dq)[0,1]/v),corr_dC_dq=float(np.corrcoef(D.dC,D.dq)[0,1]),sd={k:float(D[k].std()) for k in('dT','dC','dq')})
    # kosullu yakinsama: baslangictaki ugrak basina yuk
    m=smf.ols('dq~lnq0+lnT0+C(cont)',D).fit(cov_type='HC1'); o['conv']=dict(b_lnq0=float(m.params['lnq0']),se=float(m.bse['lnq0']),b_lnT0=float(m.params['lnT0']),se_T=float(m.bse['lnT0']),r2=float(m.rsquared))
    o['by_cont']={str(k):dict(n=int(len(g)),dT=float(np.average(g.dT,weights=g.w)),dC=float(np.average(g.dC,weights=g.w)),dq=float(np.average(g.dq,weights=g.w)),med_dq=float(g.dq.median()),up=float((g.dq>0).mean())) for k,g in D.groupby('cont')}
    o['by_inc']={str(k):dict(n=int(len(g)),dT=float(np.average(g.dT,weights=g.w)),dC=float(np.average(g.dC,weights=g.w)),dq=float(np.average(g.dq,weights=g.w)),med_dq=float(g.dq.median())) for k,g in D.groupby('inc')}
    D['qs']=pd.qcut(D.lnT0,5,labels=False); o['by_size']=[dict(q=int(k)+1,n=int(len(g)),dT=float(g.dT.median()),dC=float(g.dC.median()),dq=float(g.dq.median()),lnq0=float(np.exp(g.lnq0.median()))) for k,g in D.groupby('qs')]
    out[nm]=o
R['port']=out
# 4. yillik panel: liman ici esneklikler (ln C ve ln q'nun ln T'ye gore), liman ve yil sabit etkili
Pn=A[(A.portcalls>=100)&(A.thr>0)].copy(); n=Pn.groupby('portid').y.count(); Pn=Pn[Pn.portid.isin(n[n==7].index)]; Pn['lT']=np.log(Pn.thr); Pn['lC']=np.log(Pn.portcalls); Pn['lq']=Pn.lT-Pn.lC
def dm(df,cols):
    o=df[cols]-df.groupby('portid')[cols].transform('mean'); o=o-o.groupby(df.y).transform('mean'); return o
Z=dm(Pn,['lT','lC','lq']); bC=float((Z.lT*Z.lC).sum()/(Z.lT**2).sum()); bq=float((Z.lT*Z.lq).sum()/(Z.lT**2).sum())
cs=Pn.groupby('portid')[['lT','lC','lq']].mean(); xC=float(np.polyfit(cs.lT,cs.lC,1)[0]); xq=float(np.polyfit(cs.lT,cs.lq,1)[0])
R['elas']=dict(n_ports=int(Pn.portid.nunique()),within_calls=bC,within_q=bq,cross_calls=xC,cross_q=xq,cross_r2_q=float(np.corrcoef(cs.lT,cs.lq)[0,1]**2))
# 5. gunluk yuk hacminin yigilmasi: liman-yil bazinda CV ve tepe orani
g=full.groupby(['portid','y']).thr
S=pd.DataFrame({'mu':g.mean(),'sd':g.std(),'p95':g.quantile(.95),'p99':g.quantile(.99),'mx':g.max(),'zero':full.assign(z=(full.thr==0)).groupby(['portid','y']).z.mean()}).reset_index()
r7=full.sort_values(['portid','date']).groupby('portid').thr.rolling(7).sum().reset_index(level=0,drop=True); full=full.assign(r7=r7); S['mx7']=full.groupby(['portid','y']).r7.max().values
S=S.merge(A[['portid','y','portcalls','thr']],on=['portid','y']); S=S[(S.portcalls>=100)&(S.mu>0)]; n=S.groupby('portid').y.count(); S=S[S.portid.isin(n[n==7].index)].copy()
S['cv']=S.sd/S.mu; S['pk']=S.p95/S.mu; S['pk7']=S.mx7/(7*S.mu); S['lcv']=np.log(S.cv); S['lpk']=np.log(S.pk); S['lpk7']=np.log(S.pk7); S['lC']=np.log(S.portcalls); S['lT']=np.log(S.thr); S['lq']=S.lT-S.lC
S.to_csv('s/out/port_year_peak.csv',index=False)
def fe(dep,rhs):
    Zd=dm(S,[dep]+rhs); X=Zd[rhs].values; yv=Zd[dep].values; b=np.linalg.lstsq(X,yv,rcond=None)[0]; e=yv-X@b
    # liman kumeli standart hata
    XtXi=np.linalg.inv(X.T@X); meat=np.zeros((len(rhs),len(rhs)))
    for pid,idx in S.groupby('portid').indices.items():
        s=X[idx].T@e[idx]; meat+=np.outer(s,s)
    G=S.portid.nunique(); V=XtXi@meat@XtXi*G/(G-1); return {r:dict(b=float(b[i]),se=float(np.sqrt(V[i,i]))) for i,r in enumerate(rhs)}|{'r2w':float(1-(e**2).sum()/(yv**2).sum())}
R['peak']=dict(n_ports=int(S.portid.nunique()),n_obs=int(len(S)),med_cv=float(S.cv.median()),med_pk=float(S.pk.median()),med_pk7=float(S.pk7.median()),
    cv_TQ=fe('lcv',['lT','lq']),pk7_TQ=fe('lpk7',['lT','lq']),cv_C=fe('lcv',['lC']),cv_Cq=fe('lcv',['lC','lq']))
# kesit: limanlar arasi
cs=S.groupby('portid')[['lcv','lpk','lC','lT','lq']].mean(); m=smf.ols('lcv~lT+lq',cs).fit(cov_type='HC1'); R['peak']['cross_cv']=dict(bT=float(m.params['lT']),seT=float(m.bse['lT']),bq=float(m.params['lq']),seq=float(m.bse['lq']),r2=float(m.rsquared))
m=smf.ols('lcv~lC',cs).fit(cov_type='HC1'); R['peak']['cross_cv_C']=dict(b=float(m.params['lC']),se=float(m.bse['lC']),r2=float(m.rsquared))
# bilesik Poisson karsilastirmasi: CV^2 * (gunluk ugrak) = 1 + cv_boyut^2 ; medyan
S['k']=S.cv**2*S.portcalls/365.25; R['peak']['k_med']=float(S.k.median()); R['peak']['k_q']=[float(S.k.quantile(x)) for x in(.1,.25,.75,.9)]
# zaman icinde: ayni limanlarin medyan CV ve tepe orani
R['peak']['by_year']={int(y):dict(cv=float(g.cv.median()),pk=float(g.pk.median()),pk7=float(g.pk7.median()),wcv=float(np.average(g.cv,weights=g.thr)),wpk=float(np.average(g.pk,weights=g.thr))) for y,g in S.groupby('y')}
# karsi-olgusal: 2025 hacmi 2019 ugrak basina yukuyle tasinsaydi CV ne olurdu (liman ici esneklikle)
bq_=R['peak']['cv_TQ']['lq']['b']; a=S[S.y==2019].set_index('portid'); b=S[S.y==2025].set_index('portid'); dq=(b.lq-a.lq); R['peak']['cf']=dict(med_dq=float(dq.median()),med_effect=float(np.median(bq_*dq)),w_effect=float(np.average(bq_*dq,weights=b.thr)),p90_effect=float(np.quantile(bq_*dq,.9)),share_gt5=float((bq_*dq>np.log(1.05)).mean()))
json.dump(R,open('s/out/results.json','w'),ensure_ascii=False,indent=1)
def rd(o,n=3):
    if isinstance(o,dict): return {k:rd(v,n) for k,v in o.items()}
    if isinstance(o,list): return [rd(v,n) for v in o]
    return round(o,n) if isinstance(o,float) else o
for k in R: print('##',k,rd(R[k]))
'''

# ==============================================================================
# s/s_extra.py
# ==============================================================================
MODULES['s/s_extra.py'] = r'''
import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore')
d=pd.read_parquet('data/pwd/daily.parquet'); P=pd.read_csv('data/pw/ports_global.csv').set_index('portid')
d['y']=d.date.dt.year; d['thr']=d['import']+d['export']; d['thr_c']=d.import_container+d.export_container; R=json.load(open('s/out/results.json')); X={}
full=d[d.y<=2025]; A=full.groupby(['portid','y'])[['portcalls','thr','portcalls_container','thr_c']].sum().reset_index(); A['cont']=P.continent.reindex(A.portid).values
# kita x yil
g=A.groupby(['cont','y']).sum(numeric_only=True); X['cont_year']={c:{int(y):dict(q=float(r.thr/r.portcalls),qc=float(r.thr_c/r.portcalls_container),calls=float(r.portcalls),thr=float(r.thr)) for y,r in g.loc[c].iterrows()} for c in g.index.levels[0]}
# liman ozyukleme: dunya dln (2019->2025) guven araligi
a=A[A.y==2019].set_index('portid')[['portcalls','thr','portcalls_container','thr_c']]; b=A[A.y==2025].set_index('portid')[['portcalls','thr','portcalls_container','thr_c']]; ids=a.index.union(b.index); a=a.reindex(ids).fillna(0).values; b=b.reindex(ids).fillna(0).values
rng=np.random.default_rng(3); bs=[]
for k in range(2000):
    i=rng.integers(0,len(ids),len(ids)); s0=a[i].sum(0); s1=b[i].sum(0); bs.append([np.log(s1[1]/s0[1]),np.log(s1[0]/s0[0]),np.log((s1[1]/s1[0])/(s0[1]/s0[0])),np.log((s1[3]/s1[2])/(s0[3]/s0[2]))])
bs=np.array(bs); X['boot']={n:[float(np.percentile(bs[:,j],2.5)),float(np.percentile(bs[:,j],97.5))] for j,n in enumerate(('thr','calls','q','qc'))}; X['boot_share_q']=[float(np.percentile(bs[:,2]/bs[:,0],2.5)),float(np.percentile(bs[:,2]/bs[:,0],97.5))]
# aylik dunya serisi: ugrak basina yuk (12 aylik hareketli)
d['m']=d.date.dt.to_period('M'); M=d[d.date<='2026-10-31'].groupby('m')[['portcalls','thr','portcalls_container','thr_c']].sum(); r=M.rolling(12).sum().dropna(); pd.DataFrame({'m':r.index.astype(str),'q':r.thr/r.portcalls,'qc':r.thr_c/r.portcalls_container,'calls':r.portcalls/365.25,'thr':r.thr/365.25}).to_csv('s/out/world_roll12.csv',index=False)
lt=np.log(r.thr/r.portcalls); t=np.arange(len(lt))/12; X['trend_q']=float(np.polyfit(t,lt,1)[0]); lt=np.log(r.thr_c/r.portcalls_container); X['trend_qc']=float(np.polyfit(t,lt,1)[0])
# bolme yanliligi denetimi: tek ve cift gunlerden bagimsiz olcum
full=full.assign(odd=(full.date.dt.dayofyear%2))
def chg(sub):
    a=sub[sub.y.isin([2019,2020])].groupby('portid')[['portcalls','thr']].sum(); b=sub[sub.y.isin([2024,2025])].groupby('portid')[['portcalls','thr']].sum(); return pd.DataFrame({'dC':np.log(b.portcalls/a.portcalls),'dT':np.log(b.thr/a.thr)})
D=pd.read_csv('s/out/port_all.csv',index_col=0); o=chg(full[full.odd==1]).reindex(D.index); e=chg(full[full.odd==0]).reindex(D.index); ok=np.isfinite(o).all(1)&np.isfinite(e).all(1)
X['split']=dict(n=int(ok.sum()),corr_same=float(np.corrcoef(o.dC[ok],(o.dT-o.dC)[ok])[0,1]),corr_cross=float(np.corrcoef(o.dC[ok],(e.dT-e.dC)[ok])[0,1]),rel_dC=float(np.corrcoef(o.dC[ok],e.dC[ok])[0,1]),rel_dq=float(np.corrcoef((o.dT-o.dC)[ok],(e.dT-e.dC)[ok])[0,1]))
# tepe orani regresyonu (p99/ortalama ve en yogun gun)
S=pd.read_csv('s/out/port_year_peak.csv'); S['lp99']=np.log(S.p99/S.mu); S['lmx']=np.log(S.mx/S.mu)
def dm(df,cols):
    o=df[cols]-df.groupby('portid')[cols].transform('mean'); return o-o.groupby(df.y).transform('mean')
def fe(dep,rhs,S=S):
    Zd=dm(S,[dep]+rhs); Xm=Zd[rhs].values; yv=Zd[dep].values; b=np.linalg.lstsq(Xm,yv,rcond=None)[0]; e_=yv-Xm@b; XtXi=np.linalg.inv(Xm.T@Xm); meat=np.zeros((len(rhs),)*2)
    for pid,idx in S.groupby('portid').indices.items():
        s=Xm[idx].T@e_[idx]; meat+=np.outer(s,s)
    G=S.portid.nunique(); V=XtXi@meat@XtXi*G/(G-1); return {r_:dict(b=float(b[i]),se=float(np.sqrt(V[i,i]))) for i,r_ in enumerate(rhs)}|{'r2w':float(1-(e_**2).sum()/(yv**2).sum())}
X['mx_TQ']=fe('lmx',['lT','lq'])
S2=S[S.p95>0]; n=S2.groupby('portid').y.count(); S2=S2[S2.portid.isin(n[n==7].index)].reset_index(drop=True); X['pk_TQ']=fe('lpk',['lT','lq'],S2)|{'n_ports':int(S2.portid.nunique())}
# buyuk limanlar (ust beste bir) icin CV regresyonu
big=S.groupby('portid').lT.mean(); S3=S[S.portid.isin(big[big>=big.quantile(.8)].index)].reset_index(drop=True); X['cv_TQ_big']=fe('lcv',['lT','lq'],S3)|{'n_ports':int(S3.portid.nunique()),'med_cv':float(S3.cv.median()),'med_pk':float(S3.pk.median())}
# benzetim: bilesik Poisson altinda beklenen katsayilar (lT=-0.5, lq=+0.5)
# en buyuk degisimler: ugrak basina yuk (buyuk limanlar)
D['name']=P.fullname.reindex(D.index); Db=D[D.w>=D.w.quantile(.8)]
X['top_up']=[[r['name'],float(r.dT),float(r.dC),float(r.dq)] for _,r in Db.sort_values('dq',ascending=False).head(8).iterrows()]; X['top_dn']=[[r['name'],float(r.dT),float(r.dC),float(r.dq)] for _,r in Db.sort_values('dq').head(8).iterrows()]
top=D.sort_values('w',ascending=False).head(15); X['top15']=[[r['name'],float(r.w/730.5),float(np.exp(r.lnq0)),float(r.dT),float(r.dC),float(r.dq)] for _,r in top.iterrows()]
R['x']=X; json.dump(R,open('s/out/results.json','w'),ensure_ascii=False,indent=1)
for k in('boot','boot_share_q','trend_q','trend_qc','split','mx_TQ','pk_TQ','cv_TQ_big','top_up','top_dn','top15'): print(k,X[k])
print({c:{y:round(v['q']) for y,v in X['cont_year'][c].items() if y in(2019,2022,2025)} for c in X['cont_year']}); print({c:{y:round(v['qc']) for y,v in X['cont_year'][c].items() if y in(2019,2022,2025)} for c in X['cont_year']})
'''

# ==============================================================================
# s/s_extra2.py
# ==============================================================================
MODULES['s/s_extra2.py'] = r'''
import pandas as pd, numpy as np, json
d=pd.read_parquet('data/pwd/daily.parquet'); P=pd.read_csv('data/pw/ports_global.csv').set_index('portid'); d['y']=d.date.dt.year; d['thr']=d['import']+d['export']
A=d[d.y<=2025].groupby(['portid','y'])[['portcalls','thr','import','export','portcalls_dry_bulk','portcalls_tanker','portcalls_general_cargo','portcalls_roro','portcalls_container']].sum().reset_index(); A['iso']=P.ISO3.reindex(A.portid).values
R=json.load(open('s/out/results.json')); Z={}
def dec(S,y0=2019,y1=2025,t='thr'):
    a=S[S.y==y0][['portcalls',t]].sum(); b=S[S.y==y1][['portcalls',t]].sum(); dT=np.log(b[t]/a[t]); dC=np.log(b.portcalls/a.portcalls); return dict(dT=float(dT),dC=float(dC),dq=float(dT-dC),share=float((dT-dC)/dT))
Z['ex_china']=dec(A[A.iso!='CHN']); Z['china']=dec(A[A.iso=='CHN'])
top=A.groupby('portid').thr.sum().sort_values(ascending=False).head(20).index; Z['ex_top20']=dec(A[~A.portid.isin(top)]); Z['top20']=dec(A[A.portid.isin(top)])
Z['y2019_2024']=dec(A,2019,2024); Z['y2020_2025']=dec(A,2020,2025); Z['y2021_2025']=dec(A,2021,2025); Z['import']=dec(A,t='import'); Z['export']=dec(A,t='export')
# dengeli panel: her yil faaliyeti olan limanlar
n=A[A.portcalls>0].groupby('portid').y.count(); Z['balanced']=dec(A[A.portid.isin(n[n==7].index)])|{'n':int((n==7).sum())}
# gemi tipi ugrak paylari
W=A.groupby('y').sum(numeric_only=True); Z['type_share']={str(y):{k:float(W.loc[y,'portcalls_'+k]/W.loc[y,'portcalls']) for k in('container','dry_bulk','tanker','general_cargo','roro')} for y in(2019,2025)}
Z['type_dln']={k:float(np.log(W.loc[2025,'portcalls_'+k]/W.loc[2019,'portcalls_'+k])) for k in('container','dry_bulk','tanker','general_cargo','roro')}
R['z']=Z; json.dump(R,open('s/out/results.json','w'),ensure_ascii=False,indent=1,allow_nan=False)
for k,v in Z.items(): print(k,{a:(round(b,3) if isinstance(b,float) else b) for a,b in v.items()} if not isinstance(list(v.values())[0],dict) else v)
print(R['port']['all']['by_size'])
'''

# ==============================================================================
# s/s_disp.py
# ==============================================================================
MODULES['s/s_disp.py'] = r'''
# Gunluk ugrak sayisinin varyans/ortalama orani (Poisson gelislerinde 1) ve ugrak buyuklugu esitsizligi
import pandas as pd, numpy as np, json
d=pd.read_parquet('data/pwd/daily.parquet'); d['y']=d.date.dt.year; d=d[d.y<=2025]
S=pd.read_csv('s/out/port_year_peak.csv'); g=d.groupby(['portid','y']).portcalls; D=pd.DataFrame({'m':g.mean(),'v':g.var()}).reset_index(); D['disp']=D.v/D.m; S=S.merge(D[['portid','y','disp']],on=['portid','y'])
R=json.load(open('s/out/results.json')); S['k']=S.cv**2*S.portcalls/365.25; S['v2']=S.k-S.disp
big=S.groupby('portid').lT.transform('mean'); top=big>=big.quantile(.8)
R['disp']=dict(med=float(S.disp.median()),q10=float(S.disp.quantile(.1)),q90=float(S.disp.quantile(.9)),share_gt1=float((S.disp>1).mean()),med_big=float(S[top].disp.median()),k_med=float(S.k.median()),v2_med=float(S.v2.median()),cv_size=float(np.sqrt(max(S.v2.median(),0))),k_big=float(S[top].k.median()))
json.dump(R,open('s/out/results.json','w'),ensure_ascii=False,indent=1,allow_nan=False); print(R['disp'])
'''

# ==============================================================================
# s/s_fig.py
# ==============================================================================
MODULES['s/s_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('s/out/results.json')); X=R['x']; B='#1b4f72'; O='#b9770e'; Rd='#b03a2e'; Gn='#1e8449'; Gy='#909497'
TR={'Africa':'Afrika','Asia & Pacific':'Asya ve\nPasifik','Europe':'Avrupa','North America':'Kuzey\nAmerika','South America':'Güney\nAmerika'}
# S1
w=pd.read_csv('s/out/world_roll12.csv'); w.index=pd.PeriodIndex(w.m,freq='M').to_timestamp(); b=w.iloc[0]
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1))
a=ax[0]; a.plot(w.index,100*w.thr/b.thr,color=B,lw=1.4,label='Yük hacmi'); a.plot(w.index,100*w.calls/b.calls,color=O,lw=1.4,label='Uğrak sayısı'); a.plot(w.index,100*w.q/b.q,color=Rd,lw=1.4,label='Uğrak başına yük'); a.axhline(100,color='#999',lw=.6)
a.set_ylabel('Endeks (2019 = 100)'); a.legend(frameon=False,fontsize=7.5); a.set_title('(a) Dünya toplamı, 12 aylık hareketli',fontsize=9,loc='left')
a=ax[1]; a.plot(w.index,w.q/1000,color=B,lw=1.4,label='Tüm gemiler'); a.plot(w.index,w.qc/1000,color=O,lw=1.4,label='Konteyner gemileri'); a.set_ylabel('Uğrak başına yük (bin ton)'); a.legend(frameon=False,fontsize=7.5); a.set_title('(b) Uğrak başına yük',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('s/fig/Figure_1.png'); plt.close()
# S2
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1)); ks=['all','cont','non']; lab=['Tüm gemiler','Konteyner','Konteyner dışı']
a=ax[0]; x=np.arange(3); c=[100*R['world'][k]['dln']['calls'] for k in ks]; q=[100*R['world'][k]['dln']['q'] for k in ks]
a.bar(x,c,.55,color=O,label='Uğrak sayısı'); a.bar(x,q,.55,bottom=c,color=Rd,label='Uğrak başına yük')
for i in range(3): a.text(i,c[i]+q[i]+.25,f"{c[i]+q[i]:.1f}".replace('.',','),ha='center',fontsize=8)
a.set_xticks(x); a.set_xticklabels(lab); a.set_ylabel('Logaritmik değişim × 100, 2019–2025'); a.legend(frameon=False,fontsize=7.5,loc='upper right'); a.set_ylim(0,17); a.set_title('(a) Yük hacmi büyümesinin iki kaynağı',fontsize=9,loc='left')
a=ax[1]; wi=[100*R['ss'][k]['within_share'] for k in ks]; be=[100*R['ss'][k]['between_share'] for k in ks]; en=[100-wi[i]-be[i] for i in range(3)]
a.bar(x,wi,.55,color=B,label='Liman içi'); a.bar(x,be,.55,bottom=wi,color='#85c1e9',label='Limanlar arası'); a.bar(x,en,.55,bottom=np.array(wi)+np.array(be),color='#d5d8dc',label='Giriş–çıkış')
a.set_xticks(x); a.set_xticklabels(lab); a.set_ylabel('Uğrak başına yük artışındaki pay (%)'); a.set_ylim(0,122); a.legend(frameon=False,fontsize=7,loc='upper center',ncol=3,columnspacing=.8,handlelength=1.2); a.set_title('(b) Kaydırma–pay ayrıştırması',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('s/fig/Figure_2.png'); plt.close()
# S3
D=pd.read_csv('s/out/port_all.csv',index_col=0); P=R['port']['all']
fig,ax=plt.subplots(1,2,figsize=(7.4,3.3))
a=ax[0]; a.scatter(100*D.dC,100*D.dq,s=3+np.sqrt(D.w)/900,color=B,alpha=.35,lw=0); a.axhline(0,color='#999',lw=.6); a.axvline(0,color='#999',lw=.6); a.set_xlim(-120,120); a.set_ylim(-90,90)
for (xx,yy,k,ha) in((-115,80,'share_fewer_bigger','left'),(115,80,'share_both_up','right'),(-115,-85,'share_both_down','left'),(115,-85,'share_more_smaller','right')): a.text(xx,yy,'%'+f"{100*P[k]:.0f}",fontsize=9,ha=ha,color=Rd,weight='bold')
a.set_xlabel('Uğrak sayısındaki değişim (log × 100)'); a.set_ylabel('Uğrak başına yükteki değişim (log × 100)'); a.set_title('(a) Limanlar, 2019–20 → 2024–25',fontsize=9,loc='left')
a=ax[1]; cs=list(TR); x=np.arange(len(cs)); a.bar(x-.2,[100*P['by_cont'][c]['dC'] for c in cs],.4,color=O,label='Uğrak sayısı'); a.bar(x+.2,[100*P['by_cont'][c]['dq'] for c in cs],.4,color=Rd,label='Uğrak başına yük'); a.axhline(0,color='#999',lw=.6)
a.set_xticks(x); a.set_xticklabels([TR[c] for c in cs],fontsize=7.5); a.set_ylabel('Ağırlıklı ortalama değişim (log × 100)'); a.legend(frameon=False,fontsize=7.5); a.set_title('(b) Bölgelere göre',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('s/fig/Figure_3.png'); plt.close()
# S4
S=pd.read_csv('s/out/port_year_peak.csv'); cs=S.groupby('portid')[['lcv','lC','lT','lq']].mean()
fig,ax=plt.subplots(1,2,figsize=(7.4,3.2))
a=ax[0]; a.scatter(np.exp(cs.lC)/365.25,np.exp(cs.lcv),s=4,color=B,alpha=.35,lw=0); xx=np.logspace(-0.6,2.2,50); k=R['peak']['k_med']; a.plot(xx,np.sqrt(k/xx),color=Rd,lw=1.2,label='Eğim −0,5 (bileşik Poisson)'); a.plot(xx,np.sqrt(1/xx),color=Gy,lw=1,ls='--',label='Eşit yüklü Poisson')
a.set_xscale('log'); a.set_yscale('log'); a.set_xlabel('Günlük ortalama uğrak sayısı'); a.set_ylabel('Günlük yük hacminin değişim katsayısı'); a.legend(frameon=False,fontsize=7.5); a.set_title('(a) Limanlar arası',fontsize=9,loc='left')
a=ax[1]; names=['Değişim\nkatsayısı','%95 günü /\nortalama','En yoğun gün /\nortalama','En yoğun hafta /\nortalama']; ks=[R['peak']['cv_TQ'],X['pk_TQ'],X['mx_TQ'],R['peak']['pk7_TQ']]; x=np.arange(4)
a.bar(x-.2,[k['lT']['b'] for k in ks],.4,yerr=[1.96*k['lT']['se'] for k in ks],color=B,label='Yük hacmi esnekliği',error_kw=dict(lw=.8)); a.bar(x+.2,[k['lq']['b'] for k in ks],.4,yerr=[1.96*k['lq']['se'] for k in ks],color=Rd,label='Uğrak başına yük esnekliği',error_kw=dict(lw=.8))
a.axhline(0,color='#999',lw=.6); a.axhline(.5,color=Gy,lw=.7,ls=':'); a.axhline(-.5,color=Gy,lw=.7,ls=':'); a.set_xticks(x); a.set_xticklabels(names,fontsize=6.2); a.set_ylim(-1.0,0.62); a.set_ylabel('Liman içi esneklik'); a.legend(frameon=False,fontsize=7.5,loc='lower left'); a.set_title('(b) Liman içi, yıldan yıla',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('s/fig/Figure_4.png'); plt.close()
# S5
a0=S[S.y==2019].set_index('portid'); b0=S[S.y==2025].set_index('portid'); bq=R['peak']['cv_TQ']['lq']['b']; bT=R['peak']['cv_TQ']['lT']['b']; eq=100*bq*(b0.lq-a0.lq); eT=100*bT*(b0.lT-a0.lT)
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1))
a=ax[0]; bins=np.arange(-40,41,2.5); a.hist(eq.clip(-40,40),bins=bins,color=Rd,alpha=.7,label='Uğrak başına yük etkisi'); a.hist(eT.clip(-40,40),bins=bins,histtype='step',color=B,lw=1.2,label='Yük hacmi etkisi'); a.axvline(0,color='#999',lw=.6)
a.set_xlabel('Değişim katsayısına etki (log × 100), 2019–2025'); a.set_ylabel('Liman sayısı'); a.set_ylim(0,185); a.legend(frameon=False,fontsize=7.5,loc='upper left'); a.set_title('(a) Yığılma üzerindeki iki karşıt etki',fontsize=9,loc='left')
a=ax[1]; by=R['peak']['by_year']; ys=sorted(by,key=int); a.plot([int(y) for y in ys],[by[y]['cv'] for y in ys],'o-',color=B,label='Ortanca liman'); a.plot([int(y) for y in ys],[by[y]['wcv'] for y in ys],'s-',color=O,label='Yük ağırlıklı ortalama'); a.set_ylabel('Günlük yük hacminin değişim katsayısı'); a.set_ylim(0.5,1.5); a.legend(frameon=False,fontsize=7.5); a.set_title('(b) Zaman içinde',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('s/fig/Figure_5.png'); plt.close()
'''

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