# -*- coding: utf-8 -*-

# ------------------------------------------------------------------------------
# # Küresel Argo dizisi ne kadar sürdürülebilir? 20.619 şamandıranın sağkalım analizi, tasarım hedefinin kapsanması ve filo yenileme modeli, 2000–2026
#
# *How sustainable is the global Argo array? Survival analysis of 20,619 floats, coverage of the design target and a fleet renewal model, 2000–2026*
#
# **GOMD Yayınlar · Makale 17** — https://gomd.pirireis.edu.tr/publications.php?id=rp17
#
# Bu not defteri makaledeki **tüm sayıları, tabloları ve şekilleri** üreten kodun Colab sürümüdür. Veri, replikasyon paketinde (12 Kasım 2026 indirmesi) varsa oradan; yoksa özgün sağlayıcılardan (Argo GDAC (Coriolis/Ifremer)) indirilir. Bütün sayılar `out/results.json` dosyasına yazılır.
#
# This notebook reproduces **every number, table and figure** of the paper. Data are taken from the replication package (retrieved 12 November 2026) if present; otherwise they are downloaded from the original providers. All numbers are written to `out/results.json`.
#
# *Makale hazırlık aşamasındadır; dergi adı ve tam metin kabulden sonra eklenecektir. / The paper is in preparation; journal name and full text will be added after acceptance.*
# ------------------------------------------------------------------------------

# ------------------------------------------------------------------------------
# ## 0. Kurulum / Setup
# ------------------------------------------------------------------------------
# [Colab Magic] !pip install -q pandas numpy matplotlib lifelines pyarrow

# ------------------------------------------------------------------------------
# ## 1. Veri / Data
# Colab'da birebir yeniden üretim için: GOMD Yayınlar sayfasındaki replikasyon paketini (.zip) yükleyin. / For exact reproduction in Colab: upload the replication package (.zip) from the GOMD Publications page.
# ------------------------------------------------------------------------------
VERI_KAYNAGI = 'kaynak'   # Colab'da 'zip' secenegi icin not defterine bakin / see the notebook for the 'zip' option in Colab
# Veri / Data. Replikasyon paketindeki veri varsa onu kullanır (12 Kasım 2026 indirmesi, birebir yeniden üretim);
# yoksa özgün kaynaktan indirir (güncel sürüm; sağlayıcı güncellemeleri nedeniyle sonuçlar farklılaşabilir).
# Uses the data shipped in the replication package if present (retrieved 12 Nov 2026, exact reproduction);
# otherwise downloads from the original providers (current release; results may differ after provider updates).
import os, json, zipfile
import pandas as pd, numpy as np
def _get(url, path, timeout=600):
    if os.path.exists(path): return
    print('indiriliyor / downloading:', url)
    import urllib.request
    req = urllib.request.Request(url, headers={'User-Agent': 'gomd-replication/1.0'})
    with urllib.request.urlopen(req, timeout=timeout) as r, open(path + '.tmp', 'wb') as f:
        while True:
            b = r.read(1 << 20)
            if not b: break
            f.write(b)
    os.replace(path + '.tmp', path)
for _d in ('data', 'out', 'fig'): os.makedirs(_d, exist_ok=True)
_get('https://data-argo.ifremer.fr/ar_index_global_prof.txt.gz', 'data/argo_index.txt.gz', 1800)
_get('https://coastwatch.pfeg.noaa.gov/erddap/griddap/etopo180.csv?altitude[(-89.75):30:(89.75)][(-179.75):30:(179.75)]', 'data/etopo_30min.csv')
if not os.path.exists('data/argo.parquet'):   # profil dizini -> sıkıştırılmış tablo / profile index -> compact table
    _d = pd.read_csv('data/argo_index.txt.gz', comment='#', dtype={'date': str, 'date_update': str, 'profiler_type': str})
    _p = _d['file'].str.split('/', expand=True)
    _d['dac'] = _p[0]; _d['wmo'] = _p[1]; _d['fn'] = _p[3]; _d['mode'] = _d['fn'].str[0]
    _d['t'] = pd.to_datetime(_d['date'], format='%Y%m%d%H%M%S', errors='coerce')
    _d.drop(columns=['file', 'date', 'date_update', 'fn']).to_parquet('data/argo.parquet')
    del _d, _p

# ------------------------------------------------------------------------------
# ## 2. Yardımcı işlevler / Helper functions
# ------------------------------------------------------------------------------
# -*- coding: utf-8 -*-
import pandas as pd, numpy as np
IDX=pd.Timestamp('2026-11-12'); DEAD=120
FAM={'APEX':['845','846','847','877','0846','0845','0847'],'SOLO':['850','851','852','0851','0852'],'SOLO-II/S2A':['853','854','879','880','0853','0854'],
     'ARVOR':['844','878','837','0844'],'PROVOR':['840','841','842','836','834','835','839','0841'],'NAVIS':['863','869','0863'],'Derin':['838','862','864','874','881'],}
def fam(p):
    p=str(p).strip()
    for k,v in FAM.items():
        if p in v: return k
    return 'Diğer'
def floats():
    f=pd.read_csv('out/floats.csv',dtype={'wmo':str,'ptype':str},parse_dates=['first','last'])
    f['fam']=f.ptype.map(fam); f['event']=(f.gap>DEAD).astype(int)
    f['dur']=f.life.clip(lower=1/365.25); f['y0']=f['first'].dt.year
    f['coh']=pd.cut(f.y0,[1996,2004,2009,2014,2019,2026],labels=['2000–2004','2005–2009','2010–2014','2015–2019','2020–2026'])
    f['band']=pd.cut(f.lat0.abs(),[-1,30,60,91],labels=['0–30°','30–60°','>60°'])
    return f
def mask():
    """3x3 derece hucreler: 0.5 derecelik ETOPO noktalarinin en az yarisi 2000 m'den derin ise tasarim hucresi."""
    e=pd.read_csv('data/etopo_30min.csv',skiprows=[1])
    e['ci']=np.floor((e.latitude+90)/3).astype(int); e['cj']=np.floor((e.longitude+180)/3).astype(int)
    e['deep']=(e.altitude<-2000).astype(float); e['sea']=(e.altitude<0).astype(float)
    m=e.groupby(['ci','cj']).agg(deep=('deep','mean'),sea=('sea','mean')).reset_index()
    m['lat']=-90+3*m.ci+1.5; m['lon']=-180+3*m.cj+1.5
    m['area']=np.cos(np.deg2rad(m.lat))*m.deep   # derin kisimla agirliklandirilmis goreli alan
    m['design']=m.deep>=0.5
    return m

# ------------------------------------------------------------------------------
# ## 3. Şamandıra tablosu ve hücre sayımları / Float table and cell counts
# ------------------------------------------------------------------------------
# -*- coding: utf-8 -*-
import pandas as pd, numpy as np
d=pd.read_parquet('data/argo.parquet')
IDX=pd.Timestamp('2026-11-12')
d=d[d.t.notna()&d.latitude.notna()&d.longitude.notna()&(d.latitude.abs()<=90)&(d.longitude.abs()<=180)]
d=d[(d.t>='1997-01-01')&(d.t<=IDX)]
d=d.sort_values(['wmo','t'])
g=d.groupby('wmo')
f=g.agg(dac=('dac','first'),ptype=('profiler_type','first'),first=('t','min'),last=('t','max'),nprof=('t','size'),lat0=('latitude','first'),lon0=('longitude','first'),ocean0=('ocean','first'),latm=('latitude','median'))
f['life']=(f['last']-f['first']).dt.days/365.25
f['gap']=(IDX-f['last']).dt.days
f.to_csv('out/floats.csv'); print(len(f)); print(f.describe().round(2).to_string())
print(f.ptype.value_counts().head(25).to_dict())
# 3x3 aylik hucre sayimlari
d['ci']=(np.floor((d.latitude+90)/3)).clip(0,59).astype(int); d['cj']=(np.floor((d.longitude+180)/3)).clip(0,119).astype(int)
d['ym']=d.t.dt.year*12+d.t.dt.month-1
c=d.groupby(['ym','ci','cj']).agg(n=('t','size'),nf=('wmo','nunique')).reset_index()
c.to_parquet('out/cell_month.parquet'); print(len(c))
# cycle interval
dt=g['t'].diff().dt.total_seconds()/86400
print(dt.describe(percentiles=[.05,.25,.5,.75,.95]).round(2).to_string())

# ------------------------------------------------------------------------------
# ## 4. Sağkalım, filo, kapsama ve yenileme modeli / Survival, fleet, coverage and renewal model
# ------------------------------------------------------------------------------
# -*- coding: utf-8 -*-
import warnings; warnings.filterwarnings('ignore')
import json
from lifelines import KaplanMeierFitter, WeibullFitter, WeibullAFTFitter, CoxPHFitter
from lifelines.statistics import multivariate_logrank_test
from lifelines.utils import restricted_mean_survival_time
R={}
f=floats(); f=f[(f.y0>=2000)]
R['n_floats']=len(f); R['n_active']=int((f.event==0).sum()); R['n_profiles']=int(f.nprof.sum())
def km(sub,label=''):
    k=KaplanMeierFitter().fit(sub.dur,sub.event)
    s=lambda t: float(k.predict(t))
    return dict(n=len(sub),ev=int(sub.event.sum()),med=float(k.median_survival_time_),s2=s(2),s4=s(4),s5=s(5),s6=s(6),s8=s(8),rmst10=float(restricted_mean_survival_time(k,t=10))),k
R['km_all'],KALL=km(f)
R['km_coh']={c:km(f[f.coh==c])[0] for c in f.coh.cat.categories if (f.coh==c).sum()>0}
R['km_fam']={c:km(f[f.fam==c])[0] for c in f.fam.unique()}
R['km_dac']={c:km(f[f.dac==c])[0] for c in f.dac.unique() if (f.dac==c).sum()>=200}
R['km_band']={c:km(f[f.band==c])[0] for c in f.band.cat.categories}
lr=multivariate_logrank_test(f.dur,f.coh.astype(str),f.event); R['logrank_coh']=dict(stat=float(lr.test_statistic),p=float(lr.p_value))
# erken ariza: ilk 1 yil
R['early_fail']={c:float(1-km(f[f.coh==c])[1].predict(1)) for c in f.coh.cat.categories}
# yillik kohort medyan/5-yil sagkalim
yr={}
for y in range(2001,2022):
    s=f[f.y0==y]; r,_=km(s); yr[y]=dict(n=r['n'],med=r['med'],s5=r['s5'],s1=float(_.predict(1)))
R['km_year']=yr
# Weibull (toplu ve kohort)
wb={}
for c in list(f.coh.cat.categories)+['all']:
    s=f if c=='all' else f[f.coh==c]
    w=WeibullFitter().fit(s.dur,s.event); wb[c]=dict(lam=float(w.lambda_),rho=float(w.rho_),rho_se=float(w.summary.loc['rho_','se(coef)']))
R['weibull']=wb
# Cox
X=f[['dur','event']].copy()
for c in ['2005–2009','2010–2014','2015–2019','2020–2026']: X['coh_'+c]=(f.coh==c).astype(int)
for c in ['SOLO','SOLO-II/S2A','ARVOR','PROVOR','NAVIS','Derin','Diğer']: X['fam_'+c]=(f.fam==c).astype(int)
X['lat30_60']=(f.band=='30–60°').astype(int); X['lat60']=(f.band=='>60°').astype(int)
cx=CoxPHFitter().fit(X,'dur','event')
R['cox']={k:dict(hr=float(np.exp(v['coef'])),lo=float(np.exp(v['coef lower 95%'])),hi=float(np.exp(v['coef upper 95%'])),p=float(v['p'])) for k,v in cx.summary.iterrows()}
R['cox_conc']=float(cx.concordance_index_)
# ---- filo: aylik etkin samandira, birakma, yas
d=pd.read_parquet('data/argo.parquet'); d=d[d.t.notna()&d.latitude.notna()&d.longitude.notna()&(d.latitude.abs()<=90)&(d.longitude.abs()<=180)&(d.t<=IDX)]
d['ym']=d.t.dt.year*12+d.t.dt.month-1
d=d[d.wmo.isin(f.wmo)]
d=d.merge(f[['wmo','first','dac','fam']],on='wmo',suffixes=('','_f'))
d['age']=(d.t-d['first']).dt.days/365.25
act=d.groupby('ym').agg(nf=('wmo','nunique'),nprof=('t','size')).reset_index()
fa=d.drop_duplicates(['ym','wmo']).groupby('ym').age.agg(['mean','median',lambda x:(x>5).mean(),lambda x:(x>7).mean()]).reset_index(); fa.columns=['ym','age_mean','age_med','sh5','sh7']
act=act.merge(fa,on='ym'); act['date']=pd.to_datetime((act.ym//12).astype(str)+'-'+((act.ym%12)+1).astype(str)+'-15')
act=act[(act.date>='2000-01-01')&(act.date<'2026-11-01')]
act.to_csv('out/fleet_month.csv',index=False)
ann=act.groupby(act.date.dt.year).agg(nf=('nf','mean'),age=('age_mean','mean'),sh5=('sh5','mean'),sh7=('sh7','mean'),nprof=('nprof','sum'))
dep=f.groupby('y0').size(); ann['dep']=dep; ann.to_csv('out/fleet_year.csv')
R['fleet_year']={int(y):{k:float(v) for k,v in r.items()} for y,r in ann.iterrows()}
R['peak_fleet']=dict(val=int(act.nf.max()),date=str(act.loc[act.nf.idxmax(),'date'].date()))
# birakmalarda ulke (DAC) payi
dd=f.groupby(['y0','dac']).size().unstack(fill_value=0)
sh=dd.div(dd.sum(1),axis=0); R['dac_share']={int(y):{k:float(v) for k,v in r.items()} for y,r in sh.loc[2005:].iterrows()}
R['hhi_dep']={int(y):float((r**2).sum()) for y,r in sh.loc[2005:].iterrows()}
lastm=act.ym.max(); cur=d[d.ym==lastm].drop_duplicates('wmo')
R['active_dac_share']=cur.dac.value_counts(normalize=True).round(4).to_dict(); R['active_fam_share']=cur.fam.value_counts(normalize=True).round(4).to_dict(); R['active_n_lastmonth']=len(cur)
# ---- kapsama
m=mask(); D=m[m.design].copy(); D60=D[D.lat.abs()<60]
c=pd.read_parquet('out/cell_month.parquet')
def cover(Dm,thr=1):
    A=Dm.area.sum(); x=c.merge(Dm[['ci','cj','area']],on=['ci','cj'])
    return x[x.n>=thr].groupby('ym').area.sum()/A, x.groupby('ym').nf.sum()
cov1,_=cover(D60,1); cov3,_=cover(D60,3); covg1,_=cover(D,1); covg3,_=cover(D,3)
cv=pd.DataFrame(dict(c1=cov1,c3=cov3,g1=covg1,g3=covg3)).reset_index(); cv['date']=pd.to_datetime((cv.ym//12).astype(str)+'-'+((cv.ym%12)+1).astype(str)+'-15')
cv=cv[(cv.date>='2000-01-01')&(cv.date<'2026-11-01')]
# tasarim alanindaki samandira sayisi (ay icinde tasarim hucresinde profili olan)
d['ci']=np.floor((d.latitude+90)/3).clip(0,59).astype(int); d['cj']=np.floor((d.longitude+180)/3).clip(0,119).astype(int)
din=d.merge(D60[['ci','cj']],on=['ci','cj']); nin=din.groupby('ym').wmo.nunique(); cv['n_in']=cv.ym.map(nin)
# rastgele (alanla orantili bagimsiz yerlesim) ve ideal (her hucreye bir) kiyaslari
a=(D60.area/D60.area.sum()).values
def crand(N): return float((a*(1-(1-a)**N)).sum())
def cideal(N):
    o=np.sort(a)[::-1]; return float(o[:int(min(N,len(o)))].sum())
cv['rand']=cv.n_in.map(crand); cv['ideal']=cv.n_in.map(cideal); cv['eff']=(cv.c1-cv['rand'])/(cv['ideal']-cv['rand'])
cv.to_csv('out/coverage_month.csv',index=False)
ca=cv.groupby(cv.date.dt.year)[['c1','c3','g1','g3','n_in','rand','ideal','eff']].mean(); ca.to_csv('out/coverage_year.csv')
R['cover_year']={int(y):{k:float(v) for k,v in r.items()} for y,r in ca.iterrows()}
R['n_design60']=len(D60); R['n_design']=len(D)
# enlem kusagi kapsama (yillik ort., c1) 2010-14, 2019, 2023-26
bands=[(-90,-60),(-60,-30),(-30,0),(0,30),(30,60),(60,90)]
bc={}
for lo,hi in bands:
    Db=D[(D.lat>lo)&(D.lat<hi)]; x1,_=cover(Db,1); x1=x1.reindex(cv.ym).fillna(0); x1.index=cv.date.values
    bc[f'{lo}..{hi}']=dict(n=len(Db),**{str(y):float(x1[x1.index.year==y].mean()) for y in (2005,2010,2015,2019,2021,2023,2025,2026)})
R['cover_band']=bc
# kalici bosluklar: 2019-2026 aylarinin >=%50'sinde profilsiz tasarim hucreleri
per=c[(c.ym>=2019*12)&(c.ym<2026*12+10)].groupby(['ci','cj']).ym.nunique().rename('mo').reset_index(); nm=(2026*12+10)-2019*12
P=D.merge(per,on=['ci','cj'],how='left').fillna({'mo':0}); P['frac']=P.mo/nm; P.to_csv('out/cell_persistence.csv',index=False)
P60=P[P.lat.abs()<60]
R['persist']=dict(nm=nm,gap50_share60=float(P60[P60.frac<0.5].area.sum()/P60.area.sum()),gap50_n60=int((P60.frac<0.5).sum()),never60=int((P60.mo==0).sum()),
   full90_share60=float(P60[P60.frac>=0.9].area.sum()/P60.area.sum()),gap50_share_all=float(P[P.frac<0.5].area.sum()/P.area.sum()),
   polar_gap50=float(P[(P.lat.abs()>=60)&(P.frac<0.5)].area.sum()/P[P.lat.abs()>=60].area.sum()))
# Gini (hucre basina yillik profil, tasarim alani <60)
def gini(x,w):
    o=np.argsort(x); x=x[o]; w=w[o]; cw=np.cumsum(w)/w.sum(); cx=np.cumsum(x*w)/(x*w).sum()
    return float(1-np.sum((cx[1:]+cx[:-1])*np.diff(cw))-cx[0]*cw[0])
gy={}
for y in range(2004,2027):
    s=c[(c.ym>=y*12)&(c.ym<y*12+12)].groupby(['ci','cj']).n.sum().reset_index()
    x=D60.merge(s,on=['ci','cj'],how='left').fillna({'n':0}); gy[y]=gini((x.n/x.area).values,x.area.values)
R['gini']=gy
# ---- yenileme modeli
# guncel kohort sagkalimi: 2015-2021 birakmalari
rec=f[(f.y0>=2014)&(f.y0<=2021)]; kr=KaplanMeierFitter().fit(rec.dur,rec.event); wr=WeibullFitter().fit(rec.dur,rec.event)
lam,rho=float(wr.lambda_),float(wr.rho_)
from math import gamma
meanL=lam*gamma(1+1/rho); R['renew']=dict(n=len(rec),lam=lam,rho=rho,meanL=meanL,rmst12=float(restricted_mean_survival_time(kr,t=12)),med=float(kr.median_survival_time_))
S=lambda t: np.exp(-(np.maximum(t,0)/lam)**rho)
# tum samandiralar (her yerde) icin etkin filo = D * ortL ; hedefler
recdep=float(dep.loc[2023:2025].mean()); R['renew']['dep_recent']=recdep; R['renew']['dep_2010_19']=float(dep.loc[2010:2019].mean())
for tgt in (3000,3900,4700): R['renew'][f'need_{tgt}']=tgt/R['renew']['rmst12']
R['renew']['steady_recent']=recdep*R['renew']['rmst12']
# projeksiyon: mevcut etkin samandiralar (yas a) kosullu sagkalim + yeni birakmalar
A=f[f.event==0].copy(); ages=A.life.values; R['renew']['n_active_now']=len(A); R['renew']['age_mean_now']=float(ages.mean()); R['renew']['age_sh5_now']=float((ages>5).mean())
def project(Dyr,years=8,us_cut=0.0,us_share=None):
    out=[]
    for h in np.arange(0,years+0.001,0.25):
        old=(S(ages+h)/S(ages)).sum()
        # yeni birakmalar: yil icinde duzgun
        taus=np.arange(0.125,h,0.25); new=sum(Dyr*0.25*S(h-t) for t in taus)
        out.append((2026.87+h,old+new))
    return out
sc={'son3yil':recdep,'eksi20':0.8*recdep,'eksi40':0.6*recdep,'hedef4700':4700/meanL}
us=float(np.mean([R['dac_share'][y]['aoml'] for y in (2023,2024,2025)])); R['renew']['aoml_share_recent']=us
sc['aoml_yarim']=recdep*(1-0.5*us)
R['proj']={k:[(float(a),float(b)) for a,b in project(v)] for k,v in sc.items()}; R['proj_dep']=sc
# geriye donuk dogrulama: 2019 sonundaki filo icin 2015-oncesi sagkalimla 2025 tahmini? -> basit: model etkin sayiyi ne kadar tutturuyor
def backcast(year):
    # year sonunda: tum birakmalar x Weibull(kohorta ozgu degil, guncel) sagkalim
    t=pd.Timestamp(f'{year}-12-31'); g=f[f['first']<=t]; a=(t-g['first']).dt.days/365.25
    return float(S(a.values).sum()), int(((g['last']>=t-pd.Timedelta(days=DEAD))&(g['first']<=t)).sum())
R['backcast']={y:backcast(y) for y in (2016,2018,2020,2022,2024,2025)}
json.dump(R,open('out/results.json','w'),ensure_ascii=False,indent=1,default=float)
for k in ['n_floats','n_active','n_profiles','km_all','km_coh','km_fam','km_band','early_fail','weibull','cox_conc','peak_fleet','active_dac_share','persist','renew','proj_dep','backcast','n_design60']:
    print(k, json.dumps(R[k],ensure_ascii=False,default=float)[:900])
print('cox',{k:(round(v['hr'],2),round(v['p'],4)) for k,v in R['cox'].items()})
print(ann.round(2).to_string()); print(ca.round(3).to_string()); print({k:round(v,3) for k,v in gy.items()})
print('band',json.dumps(bc)); print('hhi',{k:round(v,3) for k,v in R['hhi_dep'].items()})
print('proj',{k:[round(x[1]) for x in v[::4]] for k,v in R['proj'].items()})

# ------------------------------------------------------------------------------
# ## 5. Ek ölçüler / Additional measures
# ------------------------------------------------------------------------------
import warnings; warnings.filterwarnings('ignore')
import json
R=json.load(open('out/results.json')); f=floats(); f=f[f.y0>=2000]
m=mask(); D=m[m.design&(m.lat.abs()<60)]; a=(D.area/D.area.sum()).values
crand=lambda N: float((a*(1-(1-a)**N)).sum()); cideal=lambda N: float(np.sort(a)[::-1][:int(min(N,len(a)))].sum())
cv=pd.read_csv('out/coverage_month.csv',parse_dates=['date']); rec=cv[cv.date>='2019-01-01']; eff=float(rec.eff.mean())
def cobs(N): return crand(N)+eff*(cideal(N)-crand(N))
def need(fn,tgt):
    for N in range(500,40000,25):
        if fn(N)>=tgt: return N
    return None
R['elastic']=dict(eff=eff,n_in_recent=float(rec.n_in.mean()),c1_recent=float(rec.c1.mean()),c3_recent=float(rec.c3.mean()),rand_recent=float(rec['rand'].mean()),
  need80_rand=need(crand,0.8),need90_rand=need(crand,0.9),need80_obs=need(cobs,0.8),need90_obs=need(cobs,0.9),need95_obs=need(cobs,0.95),
  marg_100=float(cobs(rec.n_in.mean()+100)-cobs(rec.n_in.mean())),share_in=float((rec.n_in/ pd.read_csv('out/fleet_month.csv',parse_dates=['date']).set_index('ym').reindex(rec.ym).nf.values).mean()))
# hucre basina samandira dagilimi: varyans/ortalama (son ay ornekleri)
c=pd.read_parquet('out/cell_month.parquet'); vm=[]
for ym in range(2019*12,2026*12+10):
    x=D.merge(c[c.ym==ym][['ci','cj','nf']],on=['ci','cj'],how='left').fillna({'nf':0}); vm.append(x.nf.var()/x.nf.mean())
R['elastic']['var_mean']=float(np.mean(vm))
x=D.merge(c[c.ym==2026*12+8][['ci','cj','nf']],on=['ci','cj'],how='left').fillna({'nf':0}); R['elastic']['dist_sep2026']={str(k):float((x.nf==k).mean()) for k in range(0,4)}; R['elastic']['dist_ge4']=float((x.nf>=4).mean())
# kutup mevsimselligi
Dp=m[m.design&(m.lat<=-60)]; xs=c.merge(Dp[['ci','cj','area']],on=['ci','cj']); xs=xs[(xs.ym>=2019*12)&(xs.ym<2026*12)]; s=xs.groupby('ym').area.sum()/Dp.area.sum(); mo=(s.index%12)+1
R['polar_season']=dict(south_summer=float(s[np.isin(mo,[1,2,3])].mean()),south_winter=float(s[np.isin(mo,[7,8,9])].mean()))
Dn=m[m.design&(m.lat>=60)]; xs=c.merge(Dn[['ci','cj','area']],on=['ci','cj']); xs=xs[(xs.ym>=2019*12)&(xs.ym<2026*12)]; s=xs.groupby('ym').area.sum()/Dn.area.sum(); mo=(s.index%12)+1
R['polar_season'].update(north_summer=float(s[np.isin(mo,[7,8,9])].mean()),north_winter=float(s[np.isin(mo,[1,2,3])].mean()))
# olu samandiralarda omur boyu profil
dead=f[f.event==1]; R['prof_life']={c:dict(med=float(s.nprof.median()),mean=float(s.nprof.mean())) for c,s in dead.groupby('coh')}
R['prof_all']=dict(med_dead=float(dead.nprof.median()),top10_share=float(f.sort_values('nprof',ascending=False).nprof.head(len(f)//10).sum()/f.nprof.sum()))
# erken ariza: <90 gun, <1 yil olenlerin payi (2010-2021 birakmalari)
g=f[(f.y0>=2010)&(f.y0<=2021)]; R['infant']=dict(n=len(g),lt90=float(((g.event==1)&(g.life<0.25)).mean()),lt1y=float(((g.event==1)&(g.life<1)).mean()),lost_floatyears=float(((g.event==1)&(g.life<1)).mean()))
# DAC bazinda birakma ve medyan omur tablosu
names={'aoml':'AOML (ABD)','coriolis':'Coriolis (Fransa/Avrupa)','jma':'JMA (Japonya)','csiro':'CSIRO (Avustralya)','bodc':'BODC (Birleşik Krallık)','meds':'MEDS (Kanada)','incois':'INCOIS (Hindistan)','csio':'CSIO (Çin)','kma':'KMA (Güney Kore)','kiost':'KIOST (Güney Kore)','nmdis':'NMDIS (Çin)'}
act=f[f.event==0]
R['dac_tab']=[dict(dac=k,name=names.get(k,k),n=int((f.dac==k).sum()),share=float((f.dac==k).mean()),dep2325=float(((f.dac==k)&(f.y0>=2023)&(f.y0<=2025)).sum()/3),active=int((act.dac==k).sum()),active_share=float((act.dac==k).mean()),med=R['km_dac'].get(k,{}).get('med'),s5=R['km_dac'].get(k,{}).get('s5')) for k in f.dac.value_counts().index]
R['hhi_active']=float((act.dac.value_counts(normalize=True)**2).sum())
# projeksiyon: sabit durum degerleri (Weibull ort.)
rn=R['renew']; R['renew']['need_w']={str(t):t/rn['meanL'] for t in (3000,3900,4700)}; R['renew']['steady_w']=rn['dep_recent']*rn['meanL']
json.dump(R,open('out/results.json','w'),ensure_ascii=False,indent=1,default=float)
for k in ('elastic','polar_season','prof_life','prof_all','infant','dac_tab','hhi_active'): print(k,json.dumps(R[k],ensure_ascii=False,default=lambda x:round(float(x),3))[:1800])
print(R['renew'])
import math
def _c(o):
    if isinstance(o,dict): return {k:_c(v) for k,v in o.items()}
    if isinstance(o,(list,tuple)): return [_c(v) for v in o]
    if isinstance(o,float) and (math.isnan(o) or math.isinf(o)): return None
    return o
json.dump(_c(json.load(open('out/results.json'))),open('out/results.json','w'),ensure_ascii=False,indent=1)

# ------------------------------------------------------------------------------
# ## 6. Kalan ömür, yenileme dalgası, havzalar / Residual life, renewal wave, basins
# ------------------------------------------------------------------------------
import warnings; warnings.filterwarnings('ignore')
import json, math
R=json.load(open('out/results.json')); f=floats(); f=f[f.y0>=2000]
lam,rho=R['renew']['lam'],R['renew']['rho']; S=lambda t: np.exp(-(np.maximum(t,0)/lam)**rho)
tt=np.arange(0,40,0.01)
def mrl(a): return float((S(a+tt)/S(a)).sum()*0.01)
R['mrl']={str(a):mrl(a) for a in (0,1,3,5,7,9)}
A=f[f.event==0]; ages=A.life.values
R['wave']={str(h):float((1-S(ages+h)/S(ages)).sum()) for h in (1,2,3,5)}; R['wave']['n']=len(A)
R['wave']['age_q']=[float(np.percentile(ages,q)) for q in (25,50,75,90)]; R['wave']['old7']=int((ages>7).sum()); R['wave']['old10']=int((ages>10).sum())
# havza
d=pd.read_parquet('data/argo.parquet'); d=d[d.t.notna()&d.latitude.notna()&d.longitude.notna()&(d.latitude.abs()<=90)&(d.longitude.abs()<=180)]
d['ci']=np.floor((d.latitude+90)/3).clip(0,59).astype(int); d['cj']=np.floor((d.longitude+180)/3).clip(0,119).astype(int)
b=d.groupby(['ci','cj']).ocean.agg(lambda x:x.value_counts().index[0]).rename('basin').reset_index()
m=mask(); D=m[m.design&(m.lat.abs()<60)].merge(b,on=['ci','cj'],how='left'); D['basin']=D.basin.fillna('?')
c=pd.read_parquet('out/cell_month.parquet'); out={}
for k,nm in (('A','Atlantik'),('P','Pasifik'),('I','Hint')):
    Db=D[D.basin==k]; x=c.merge(Db[['ci','cj','area']],on=['ci','cj']); cov=x.groupby('ym').area.sum()/Db.area.sum(); nf=x.groupby('ym').nf.sum()
    yr=lambda y:(cov[(cov.index>=y*12)&(cov.index<y*12+12)].mean(), nf[(nf.index>=y*12)&(nf.index<y*12+12)].mean())
    out[nm]=dict(n=len(Db),area=float(Db.area.sum()/D.area.sum()),**{f'c{y}':float(yr(y)[0]) for y in (2010,2019,2025)},dens2025=float(yr(2025)[1]/len(Db)))
R['basin']=out
dt=d.sort_values(['wmo','t']).groupby('wmo').t.diff().dt.total_seconds()/86400; dt=dt.dropna()
R['cycle']=dict(med=float(dt.median()),sh9_11=float(((dt>=9)&(dt<=11)).mean()),sh_lt5=float((dt<5).mean()))
def _c(o):
    if isinstance(o,dict): return {k:_c(v) for k,v in o.items()}
    if isinstance(o,(list,tuple)): return [_c(v) for v in o]
    if isinstance(o,float) and (math.isnan(o) or math.isinf(o)): return None
    return o
json.dump(_c(R),open('out/results.json','w'),ensure_ascii=False,indent=1)
for k in ('mrl','wave','basin','cycle'): print(k,json.dumps(R[k],ensure_ascii=False,default=float))

# ------------------------------------------------------------------------------
# ## 7. Hücre boyutu ve konum tanımına duyarlılık / Sensitivity to cell size and position definition
# ------------------------------------------------------------------------------
import warnings; warnings.filterwarnings('ignore')
import json, math
R=json.load(open('out/results.json'))
d=pd.read_parquet('data/argo.parquet'); d=d[d.t.notna()&d.latitude.notna()&d.longitude.notna()&(d.latitude.abs()<=90)&(d.longitude.abs()<=180)&(d.t>='2019-01-01')&(d.t<'2026-11-01')]
d['ym']=d.t.dt.year*12+d.t.dt.month-1
e=pd.read_csv('data/etopo_30min.csv',skiprows=[1]); e['deep']=(e.altitude<-2000).astype(float)
def run(size,one):
    nlat=int(180/size); nlon=int(360/size)
    e['ci']=np.floor((e.latitude+90)/size).astype(int); e['cj']=np.floor((e.longitude+180)/size).astype(int)
    m=e.groupby(['ci','cj']).deep.mean().reset_index(); m['lat']=-90+size*m.ci+size/2; m['area']=np.cos(np.deg2rad(m.lat))*m.deep
    D=m[(m.deep>=0.5)&(m.lat.abs()<60)]; a=(D.area/D.area.sum()).values
    x=d.sort_values('t').drop_duplicates(['ym','wmo']) if one else d
    x=x.assign(ci=np.floor((x.latitude+90)/size).clip(0,nlat-1).astype(int),cj=np.floor((x.longitude+180)/size).clip(0,nlon-1).astype(int)).merge(D[['ci','cj','area']],on=['ci','cj'])
    N=x.groupby('ym').wmo.nunique(); cov=x.drop_duplicates(['ym','ci','cj']).groupby('ym').area.sum()/D.area.sum()
    rand=N.map(lambda n: float((a*(1-(1-a)**n)).sum())); ideal=N.map(lambda n: float(np.sort(a)[::-1][:int(min(n,len(a)))].sum()))
    eff=((cov-rand)/(ideal-rand)).replace([np.inf,-np.inf],np.nan)
    return dict(cells=len(D),N=float(N.mean()),cov=float(cov.mean()),rand=float(rand.mean()),ideal=float(ideal.mean()),eff=float(eff.mean()))
R['rob_cov']={'3_all':run(3,False),'3_one':run(3,True),'6_all':run(6,False),'6_one':run(6,True),'2_one':run(2,True)}
json.dump(R,open('out/results.json','w'),ensure_ascii=False,indent=1)
for k,v in R['rob_cov'].items(): print(k,{a:round(b,3) for a,b in v.items()})

# ------------------------------------------------------------------------------
# ## 8. Anlık kapsama ve dağılım etkinliği / Snapshot coverage and placement efficiency
# ------------------------------------------------------------------------------
# -*- coding: utf-8 -*-
"""Anlik (samandira basina ayda tek konum) kapsama ve dagilim etkinligi. Rastgele/kusursuz kiyaslar her samandirayi
tek hucreye yerlestirdiginden, gozlenen kapsama da samandiranin aydaki ilk profilinin konumuyla hesaplanir."""
import warnings; warnings.filterwarnings('ignore')
import json
R=json.load(open('out/results.json'))
d=pd.read_parquet('data/argo.parquet'); d=d[d.t.notna()&d.latitude.notna()&d.longitude.notna()&(d.latitude.abs()<=90)&(d.longitude.abs()<=180)&(d.t>='2000-01-01')&(d.t<'2026-11-01')]
d['ym']=d.t.dt.year*12+d.t.dt.month-1
x=d.sort_values('t').drop_duplicates(['ym','wmo'])
x['ci']=np.floor((x.latitude+90)/3).clip(0,59).astype(int); x['cj']=np.floor((x.longitude+180)/3).clip(0,119).astype(int)
m=mask(); D=m[m.design&(m.lat.abs()<60)]; a=(D.area/D.area.sum()).values
crand=lambda N: float((a*(1-(1-a)**N)).sum()); cideal=lambda N: float(np.sort(a)[::-1][:int(min(N,len(a)))].sum())
xi=x.merge(D[['ci','cj','area']],on=['ci','cj'])
N=xi.groupby('ym').wmo.nunique(); s1=xi.drop_duplicates(['ym','ci','cj']).groupby('ym').area.sum()/D.area.sum()
cnt=xi.groupby(['ym','ci','cj']).size().rename('k').reset_index()
cv=pd.read_csv('out/coverage_month.csv',parse_dates=['date'])
cv['n_in']=cv.ym.map(N); cv['s1']=cv.ym.map(s1); cv['rand']=cv.n_in.map(crand); cv['ideal']=cv.n_in.map(cideal); cv['eff']=(cv.s1-cv['rand'])/(cv['ideal']-cv['rand'])
cv.to_csv('out/coverage_month.csv',index=False)
ca=cv.groupby(cv.date.dt.year)[['c1','c3','g1','g3','s1','n_in','rand','ideal','eff']].mean(); ca.to_csv('out/coverage_year.csv')
R['cover_year']={str(int(y)):{k:float(v) for k,v in r.items()} for y,r in ca.iterrows()}
rec=cv[cv.date>='2019-01-01']; eff=float(rec.eff.mean())
cobs=lambda n: crand(n)+eff*(cideal(n)-crand(n))
def need(fn,tgt):
    for n in range(500,60000,25):
        if fn(n)>=tgt: return n
fm=pd.read_csv('out/fleet_month.csv').set_index('ym')
E=R['elastic']; E.update(eff=eff,n_in_recent=float(rec.n_in.mean()),s1_recent=float(rec.s1.mean()),c1_recent=float(rec.c1.mean()),c3_recent=float(rec.c3.mean()),rand_recent=float(rec['rand'].mean()),
  need80_rand=need(crand,0.8),need90_rand=need(crand,0.9),need80_obs=need(cobs,0.8),need90_obs=need(cobs,0.9),marg_100=float(cobs(rec.n_in.mean()+100)-cobs(rec.n_in.mean())),
  share_in=float((rec.n_in/fm.reindex(rec.ym).nf.values).mean()),eff_min1625=float(ca.loc[2016:2025,'eff'].min()),eff_max1625=float(ca.loc[2016:2025,'eff'].max()),
  s1_min1625=float(ca.loc[2016:2025,'s1'].min()),s1_max1625=float(ca.loc[2016:2025,'s1'].max()))
# hucre basina samandira dagilimi, Eylul 2026 (anlik) ve ayni N icin rastgele beklenti (benzetim)
ym=2026*12+8; k=D.merge(cnt[cnt.ym==ym][['ci','cj','k']],on=['ci','cj'],how='left').fillna({'k':0}).k.values; n=int(N[ym])
rng=np.random.default_rng(1); sim=np.zeros(5)
for _ in range(300):
    z=np.bincount(rng.choice(len(a),size=n,p=a),minlength=len(a)); sim+=np.array([(z==0).mean(),(z==1).mean(),(z==2).mean(),(z==3).mean(),(z>=4).mean()])
sim/=300
E['dist_sep2026']={str(i):float((k==i).mean()) for i in range(4)}; E['dist_ge4']=float((k>=4).mean()); E['dist_rand']=[float(v) for v in sim]; E['n_sep2026']=n
# varyans/ortalama: gozlenen ve rastgele (2019-2026 ort.)
vm=[]; 
for y in range(2019*12,2026*12+10):
    kk=D.merge(cnt[cnt.ym==y][['ci','cj','k']],on=['ci','cj'],how='left').fillna({'k':0}).k.values; vm.append(kk.var()/kk.mean())
z=np.array([np.bincount(rng.choice(len(a),size=n,p=a),minlength=len(a)) for _ in range(200)]); E['var_mean']=float(np.mean(vm)); E['var_mean_rand']=float((z.var(axis=1)/z.mean(axis=1)).mean())
json.dump(R,open('out/results.json','w'),ensure_ascii=False,indent=1)
print({k:(round(v,3) if isinstance(v,float) else v) for k,v in E.items()})
print(ca[['c1','c3','s1','n_in','rand','eff']].round(3).to_string())

# ------------------------------------------------------------------------------
# ## 9. Şekiller / Figures
# ------------------------------------------------------------------------------
# -*- coding: utf-8 -*-
import warnings; warnings.filterwarnings('ignore')
import json, matplotlib; matplotlib.use('Agg'); import matplotlib.pyplot as plt
from lifelines import KaplanMeierFitter
plt.rcParams.update({'font.family':'DejaVu Sans','font.size':9,'axes.spines.top':False,'axes.spines.right':False,'savefig.dpi':300})
C=['#1b4f72','#b03a2e','#1e8449','#b9770e','#6c3483','#566573']
R=json.load(open('out/results.json')); f=floats(); f=f[f.y0>=2000]
# Sekil 1: KM kohort + aile
fig,ax=plt.subplots(1,2,figsize=(7.4,3.2))
for i,c in enumerate(f.coh.cat.categories):
    s=f[f.coh==c]; k=KaplanMeierFitter().fit(s.dur,s.event); t=k.survival_function_
    ax[0].step(t.index,t.iloc[:,0],where='post',color=C[i],lw=1.4,label=c)
ax[0].set_xlim(0,12); ax[0].set_ylim(0,1); ax[0].set_xlabel('Bırakmadan sonra geçen süre (yıl)'); ax[0].set_ylabel('Sağkalım olasılığı'); ax[0].legend(title='Bırakma dönemi',frameon=False,fontsize=8); ax[0].set_title('(a) Bırakma dönemine göre',fontsize=9,loc='left')
for i,c in enumerate(['APEX','SOLO-II/S2A','ARVOR','NAVIS','PROVOR','Derin']):
    s=f[(f.fam==c)&(f.y0>=2010)]; k=KaplanMeierFitter().fit(s.dur,s.event); t=k.survival_function_
    ax[1].step(t.index,t.iloc[:,0],where='post',color=C[i],lw=1.4,label=f'{c} (n={len(s)})'.replace(',','.'))
ax[1].set_xlim(0,12); ax[1].set_ylim(0,1); ax[1].set_xlabel('Bırakmadan sonra geçen süre (yıl)'); ax[1].legend(title='Şamandıra ailesi, 2010 ve sonrası',frameon=False,fontsize=7.5); ax[1].set_title('(b) Şamandıra ailesine göre',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('fig/Figure_1.png'); plt.close()
# Sekil 2: filo, birakma, yas
a=pd.read_csv('out/fleet_month.csv',parse_dates=['date']); y=pd.read_csv('out/fleet_year.csv',index_col=0)
fig,ax=plt.subplots(1,2,figsize=(7.4,3.0))
ax[0].plot(a.date,a.nf,color=C[0],lw=1.2); ax[0].set_ylabel('Ay içinde profil gönderen şamandıra',color=C[0]); ax0=ax[0].twinx(); ax0.spines['right'].set_visible(True)
yy=y.loc[2000:2025]; ax0.bar(pd.to_datetime(yy.index.astype(str)+'-07-01'),yy.dep,width=250,color=C[3],alpha=.45); ax0.set_ylabel('Yıllık bırakma (adet)',color=C[3]); ax0.set_ylim(0,2600)
ax[0].set_title('(a) Filo büyüklüğü ve bırakmalar',fontsize=9,loc='left')
ax[1].plot(a.date,a.age_mean,color=C[1],lw=1.2,label='Ortalama yaş (yıl)'); ax[1].set_ylabel('Etkin şamandıraların ortalama yaşı (yıl)',color=C[1]); ax1=ax[1].twinx(); ax1.spines['right'].set_visible(True)
ax1.plot(a.date,100*a.sh5,color=C[4],lw=1.2); ax1.set_ylabel('5 yaşından büyük şamandıra payı (%)',color=C[4]); ax[1].set_title('(b) Filonun yaşlanması',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('fig/Figure_2.png'); plt.close()
# Sekil 3: kapsama
cv=pd.read_csv('out/coverage_month.csv',parse_dates=['date']); cv=cv[cv.date>='2002-01-01']
fig,ax=plt.subplots(1,2,figsize=(7.4,3.0))
ax[0].plot(cv.date,100*cv.ideal,color='#999',lw=1,ls='--',label='Kusursuz yerleşim'); ax[0].plot(cv.date,100*cv.c1,color=C[2],lw=1.0,label='Ay içinde ≥1 profil alınan alan'); ax[0].plot(cv.date,100*cv['rand'],color=C[1],lw=1.2,label='Rastgele yerleşim beklentisi'); ax[0].plot(cv.date,100*cv.s1,color=C[0],lw=1.4,label='Anlık kapsama (şamandıra başına tek konum)')
ax[0].set_ylabel('Kapsanan tasarım alanı (%)'); ax[0].set_ylim(0,105); ax[0].legend(frameon=False,fontsize=7,loc='lower right'); ax[0].set_title('(a) 60°G–60°K tasarım alanının aylık kapsaması',fontsize=9,loc='left')
ax[1].scatter(cv.n_in,100*cv.s1,s=5,color=C[0],alpha=.6,label='Gözlenen aylar (anlık kapsama)'); n=np.arange(0,4300,50)
D=mask(); D=D[D.design&(D.lat.abs()<60)]; aa=(D.area/D.area.sum()).values
ax[1].plot(n,[100*(aa*(1-(1-aa)**k)).sum() for k in n],color=C[1],lw=1.2,label='Rastgele yerleşim'); ax[1].plot(n,[100*np.sort(aa)[::-1][:int(min(k,len(aa)))].sum() for k in n],color='#999',ls='--',lw=1,label='Kusursuz yerleşim')
ax[1].set_xlabel('Tasarım alanındaki şamandıra sayısı'); ax[1].set_ylabel('Kapsanan tasarım alanı (%)'); ax[1].legend(frameon=False,fontsize=7.5,loc='lower right'); ax[1].set_title('(b) Şamandıra sayısı–kapsama ilişkisi',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('fig/Figure_3.png'); plt.close()
# Sekil 4: harita
P=pd.read_csv('out/cell_persistence.csv'); m=mask()
G=np.full((60,120),np.nan); 
for r in P.itertuples(): G[r.ci,r.cj]=100*r.frac
L=np.full((60,120),np.nan); 
for r in m.itertuples():
    if r.sea<0.5: L[r.ci,r.cj]=1
fig,ax=plt.subplots(figsize=(7.4,3.7))
ax.pcolormesh(np.arange(-180,181,3),np.arange(-90,91,3),L,cmap='Greys',vmin=0,vmax=2.2)
pc=ax.pcolormesh(np.arange(-180,181,3),np.arange(-90,91,3),G,cmap='viridis',vmin=0,vmax=100)
cb=fig.colorbar(pc,ax=ax,fraction=0.025,pad=0.02); cb.set_label('Profil alınan ayların payı (%), 2019–2026')
ax.axhline(60,color='k',lw=.5,ls=':'); ax.axhline(-60,color='k',lw=.5,ls=':'); ax.set_xlabel('Boylam'); ax.set_ylabel('Enlem'); ax.set_xlim(-180,180); ax.set_ylim(-81,90)
for s in ax.spines.values(): s.set_visible(True)
fig.tight_layout(); fig.savefig('fig/Figure_4.png'); plt.close()
# Sekil 5: projeksiyon
fig,ax=plt.subplots(figsize=(5.6,3.2))
lab={'hedef4700':'Hedefe göre bırakma','son3yil':'2023–2025 ortalaması','eksi20':'%20 azalma','aoml_yarim':'En büyük veri merkezinin bırakmaları yarıya iner','eksi40':'%40 azalma'}
for i,(k,l) in enumerate(lab.items()):
    p=np.array(R['proj'][k]); ax.plot(p[:,0],p[:,1],color=C[i],lw=1.4,label=f"{l} ({R['proj_dep'][k]:.0f}/yıl)")
ax.axhline(4700,color='#999',ls='--',lw=.8); ax.text(2027,4730,'OneArgo hedefi: 4.700',fontsize=7.5,color='#666'); ax.axhline(3000,color='#999',ls=':',lw=.8); ax.text(2027,3030,'İlk tasarım: 3.000',fontsize=7.5,color='#666')
ax.set_ylabel('Etkin şamandıra sayısı'); ax.set_ylim(2500,5300); ax.legend(frameon=False,fontsize=7.5,loc='lower left',bbox_to_anchor=(0,0.08)); fig.tight_layout(); fig.savefig('fig/Figure_5.png'); plt.close()
# duyarlilik: olum esigi
out={}
for dead in (60,120,180,365):
    g=f.copy(); g['event']=(g.gap>dead).astype(int); r=g[(g.y0>=2014)&(g.y0<=2021)]; k=KaplanMeierFitter().fit(r.dur,r.event)
    from lifelines.utils import restricted_mean_survival_time as rm
    out[dead]=dict(active=int((g.event==0).sum()),med=float(k.median_survival_time_),rmst12=float(rm(k,t=12)))
# ayni samandiranin 120 gunden uzun sessizlikten sonra geri donmesi
d=pd.read_parquet('data/argo.parquet'); d=d[d.t.notna()].sort_values(['wmo','t']); gp=d.groupby('wmo').t.diff().dt.days
out['ret120']=float((gp>120).groupby(d.wmo).any().mean()); out['ret365']=float((gp>365).groupby(d.wmo).any().mean())
R['sens_dead']=out; json.dump(R,open('out/results.json','w'),ensure_ascii=False,indent=1,default=float); print(out)

# ------------------------------------------------------------------------------
# ## Şekiller / Figures
# ------------------------------------------------------------------------------
import glob
print('Sekiller / figures:', sorted(glob.glob('fig/*.png')))