# -*- coding: utf-8 -*-
# ------------------------------------------------------------------------------
# Zenginlik kestiricileri sonradan bulunanı öngörüyor mu? 160 deniz alanında deniz canlısı kayıtlarıyla bir geri-test
# Do richness estimators anticipate what is found later? A back-test with marine occurrence records from 160 sea areas
#
# GOMD Yayınlar · Makale 28 — https://gomd.pirireis.edu.tr/publications.php?id=rp28
#
# Bu dosya makalenin bütün modüllerini içerir. Replikasyon paketi klasöründe çalıştırıldığında modülleri yerlerine yazar ve sırayla, her birini ayrı süreçte çalıştırır;
# sonuçları `paper/out/` ve şekilleri `paper/fig/` klasörüne yazar.
# This file contains all modules of the paper. Run in the replication-package folder, it writes the modules to their places and runs them in order, each in its own process;
# it writes results to `paper/out/` and figures to `paper/fig/`.
#
# Kurulum / Setup:  pip install pandas numpy matplotlib
# Çalıştırma / Run:  python rp28_richness_estimators_backtest_marine_inventories.py
#
# Occurrence data: OBIS (2026) Ocean Biodiversity Information System. Intergovernmental Oceanographic Commission of UNESCO. obis.org ; area checklists retrieved from the public service on 13–14 November 2026. OBIS datasets are licensed CC0, CC BY or CC BY-NC; the checklists in this package are aggregated counts of records per species and area and are provided for non-commercial research use with attribution to OBIS and its data providers.
# 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 = ['obis/obis_analysis.py', 'obis/obis_deep.py', 'paper/fig30en.py', 'obis/obis_thr.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
DIRS = ['obis', 'obis/out', 'obis/c', 'paper', 'paper/out', 'paper/fig', 'raw']

# ==============================================================================
# obis/obis_analysis.py
# ==============================================================================
MODULES['obis/obis_analysis.py'] = r'''
# P30: OBIS tur envanterinin tamligi (66 LME + 101 IHO denizi), Chao1 ve 1999 geri-testi
import json, glob, os, warnings, numpy as np, pandas as pd
warnings.filterwarnings('ignore')
FISH = {'Teleostei', 'Elasmobranchii', 'Holocephali', 'Myxini', 'Petromyzonti', 'Coelacanthi', 'Chondrostei', 'Holostei', 'Cladistii', 'Dipneusti', 'Actinopteri', 'Actinopterygii'}
TETRA = {'Mammalia', 'Aves', 'Reptilia', 'Testudines', 'Squamata', 'Crocodylia'}
def group(ph, cl):
    if cl in FISH: return 'fish'
    if cl in TETRA: return 'tetrapod'
    if ph == 'Mollusca': return 'mollusc'
    if ph == 'Arthropoda': return 'arthropod'
    if ph == 'Annelida': return 'annelid'
    if ph == 'Cnidaria': return 'cnidarian'
    if ph == 'Echinodermata': return 'echinoderm'
    if ph == 'Porifera': return 'sponge'
    return 'other'
def chao(n):
    n = np.asarray(n); n = n[n > 0]; S = len(n); N = int(n.sum()); f1 = int((n == 1).sum()); f2 = int((n == 2).sum())
    if S == 0: return dict(S=0, N=0, f1=0, f2=0, chao=np.nan, C=np.nan, cov=np.nan, lo=np.nan, hi=np.nan)
    est = S + f1 * (f1 - 1) / (2.0 * (f2 + 1))
    # Chao 1987 log-normal aralik
    if f2 > 0 and f1 > 0:
        r = f1 / f2; var = f2 * (0.5 * r ** 2 + r ** 3 + 0.25 * r ** 4); T = est - S
        K = np.exp(1.96 * np.sqrt(np.log(1 + var / max(T, 1e-9) ** 2))) if T > 0 else 1.0
        lo, hi = S + T / K, S + T * K
    else: lo = hi = est
    cov = 1 - (f1 / N) * (((N - 1) * f1) / ((N - 1) * f1 + 2 * f2)) if (f1 > 0 and N > 1) else 1.0
    return dict(S=S, N=N, f1=f1, f2=f2, chao=est, C=S / est, cov=cov, lo=lo, hi=hi)
rows = []; SP = {}
for f in sorted(glob.glob('obis/c/*.json')):
    d = json.load(open(f)); a = d['area']; tag = d['tag']
    df = pd.DataFrame(d['rows'], columns=['id', 'rank', 'n', 'phylum', 'cls', 'kingdom'])
    df = df[(df['rank'] == 'Species') & (df.kingdom == 'Animalia')].copy(); dup = int(df.id.duplicated().sum()); df = df.groupby('id', as_index=False).agg(n=('n', 'sum'), phylum=('phylum', 'first'), cls=('cls', 'first')); df['g'] = [group(p, c) for p, c in zip(df.phylum, df.cls)]
    SP[(a['id'], tag)] = df.set_index('id')[['n', 'g']]
    for g, sub in [('all', df)] + list(df.groupby('g')):
        r = chao(sub.n.values); r.update(area=a['id'], name=a['name'], type=a['type'], tag=tag, group=g, records_total=d['stats'].get('records'), datasets=d['stats'].get('datasets')); rows.append(r)
T = pd.DataFrame(rows); T.to_csv('obis/out/completeness.csv', index=False)
R = {}
cur = T[(T.tag == 'all') & (T.group == 'all')].copy()
R['n_areas'] = {'lme': int((cur.type == 'lme').sum()), 'iho': int((cur.type == 'iho').sum())}
def dist(x): x = x.dropna(); return {'n': int(len(x)), 'median': float(x.median()), 'p10': float(x.quantile(.1)), 'p25': float(x.quantile(.25)), 'p75': float(x.quantile(.75)), 'p90': float(x.quantile(.9)), 'min': float(x.min()), 'max': float(x.max())}
ok = cur[cur.S >= 100]
R['now'] = {t: {'completeness': dist(ok[ok.type == t].C), 'coverage': dist(ok[ok.type == t]['cov']), 'singleton_share': dist((ok[ok.type == t].f1 / ok[ok.type == t].S)), 'S': dist(ok[ok.type == t].S), 'N': dist(ok[ok.type == t].N),
                    'missing_species_sum': float((ok[ok.type == t].chao - ok[ok.type == t].S).sum())} for t in ('lme', 'iho')}
R['now']['excluded_lt100'] = cur[cur.S < 100][['name', 'type', 'S']].values.tolist()
R['now']['share_below_70'] = {t: float((ok[ok.type == t].C < 0.7).mean()) for t in ('lme', 'iho')}
R['now']['share_above_90'] = {t: float((ok[ok.type == t].C > 0.9).mean()) for t in ('lme', 'iho')}
lst = ok.sort_values('C'); cols = ['name', 'type', 'S', 'N', 'f1', 'f2', 'chao', 'C', 'cov']
R['lowest15'] = json.loads(lst.head(15)[cols].to_json(orient='records')); R['highest15'] = json.loads(lst.tail(15)[cols].to_json(orient='records'))
R['corr'] = {'spearman_C_logN': float(ok.C.rank().corr(np.log(ok.N).rank())), 'spearman_C_logS': float(ok.C.rank().corr(np.log(ok.S).rank())), 'spearman_C_rec_per_sp': float(ok.C.rank().corr((ok.N / ok.S).rank()))}
# gruplar
gg = T[(T.tag == 'all') & (T.group != 'all') & (T.S >= 30)]
R['groups'] = {g: {'areas': int(len(s)), 'C_median': float(s.C.median()), 'C_p10': float(s.C.quantile(.1)), 'C_p90': float(s.C.quantile(.9)), 'singleton_share_median': float((s.f1 / s.S).median()), 'S_median': float(s.S.median()),
                   'rec_per_species_median': float((s.N / s.S).median())} for g, s in gg.groupby('group')}
# zaman: kumulatif kesitler
for tag in ('to1999', 'to2009', 'to2019', 'dated', 'all'):
    s = T[(T.tag == tag) & (T.group == 'all') & T.area.isin(ok.area)]
    R.setdefault('by_cut', {})[tag] = {'C_median': float(s.C.median()), 'S_median': float(s.S.median()), 'N_median': float(s.N.median()), 'S_sum': float(s.S.sum()), 'singleton_share_median': float((s.f1 / s.S).median()), 'cov_median': float(s['cov'].median())}
# geri-test: 1999'a kadar toplanan kayitlarla Chao1 -> 2026'da gozlenen
A = T[(T.group == 'all')].pivot(index='area', columns='tag', values=['S', 'chao', 'lo', 'hi', 'C', 'N', 'f1']).loc[ok.area]
bt = pd.DataFrame({'S99': A['S']['to1999'], 'chao99': A['chao']['to1999'], 'hi99': A['hi']['to1999'], 'S26': A['S']['dated'], 'N99': A['N']['to1999'], 'N26': A['N']['dated'], 'C99': A['C']['to1999'], 'C26': A['C']['dated'], 'Sall': A['S']['all'], 'Nall': A['N']['all'], 'S09': A['S']['to2009'], 'chao09': A['chao']['to2009'], 'hi09': A['hi']['to2009'], 'S19': A['S']['to2019'], 'chao19': A['chao']['to2019'], 'hi19': A['hi']['to2019']})
bt = bt[bt.S99 >= 100]; bt['name'] = ok.set_index('area').name.reindex(bt.index); bt['type'] = ok.set_index('area').type.reindex(bt.index)
bt['growth'] = bt.S26 / bt.S99; bt['pred_growth'] = bt.chao99 / bt.S99; bt['exceed'] = bt.S26 > bt.chao99; bt['exceed_hi'] = bt.S26 > bt.hi99
bt.to_csv('obis/out/backtest.csv')
R['undated'] = {'species_share_only_undated_median': float((1 - bt.S26 / bt.Sall).median()), 'p90': float((1 - bt.S26 / bt.Sall).quantile(.9)), 'record_share_undated_median': float((1 - bt.N26 / bt.Nall).median())}
R['backtest99'] = {'n': int(len(bt)), 'share_observed_exceeds_chao': float(bt.exceed.mean()), 'share_exceeds_upper95': float(bt.exceed_hi.mean()), 'median_S26_over_chao99': float((bt.S26 / bt.chao99).median()),
                   'p90_S26_over_chao99': float((bt.S26 / bt.chao99).quantile(.9)), 'median_growth': float(bt.growth.median()), 'median_pred_growth': float(bt.pred_growth.median()),
                   'spearman_pred_vs_real_growth': float(bt.pred_growth.rank().corr(bt.growth.rank())), 'pearson_log': float(np.corrcoef(np.log(bt.pred_growth), np.log(bt.growth))[0, 1]),
                   'median_C99': float(bt.C99.median()), 'median_true_share99': float((bt.S99 / bt.S26).median()), 'median_record_growth': float((bt.N26 / bt.N99).median())}
b = np.polyfit(np.log(bt.pred_growth), np.log(bt.growth), 1); R['backtest99']['loglog_slope_intercept'] = [float(b[0]), float(b[1])]
# kayit artisi kontrol edilerek: log buyume ~ log ongorulen + log kayit artisi
X = np.column_stack([np.ones(len(bt)), np.log(bt.pred_growth), np.log(bt.N26 / bt.N99)]); beta, *_ = np.linalg.lstsq(X, np.log(bt.growth), rcond=None)
res = np.log(bt.growth) - X @ beta; R['backtest99']['multi'] = {'b_pred': float(beta[1]), 'b_records': float(beta[2]), 'r2': float(1 - res.var() / np.log(bt.growth).var())}
for a_, b_, tg in (('S09', 'chao09', 'hi09'), ('S19', 'chao19', 'hi19')):
    q = bt[bt[a_] >= 100]; R['backtest_' + a_] = {'n': int(len(q)), 'share_observed_exceeds_chao': float((q.S26 > q[b_]).mean()), 'share_exceeds_upper95': float((q.S26 > q[tg]).mean()), 'median_S26_over_chao': float((q.S26 / q[b_]).median())}
# 1999 teklerinin akibeti ve yeni turlerin bollugu
fate = []
for aid in bt.index:
    s0 = SP[(aid, 'to1999')]; s1 = SP[(aid, 'dated')]
    one = s0.index[s0.n == 1]; n1 = s1.n.reindex(one)
    new = s1.index.difference(s0.index); nn = s1.n.reindex(new)
    fate.append({'area': aid, 'single99': len(one), 'still_single': float((n1 == 1).mean()) if len(one) else np.nan, 'new_species': len(new), 'new_single': float((nn == 1).mean()) if len(new) else np.nan,
                 'new_ge10': float((nn >= 10).mean()) if len(new) else np.nan, 'new_per_1000rec': 1000 * len(new) / max(s1.n.sum() - s0.n.sum(), 1),
                 'new_fish_share': float((s1.g.reindex(new) == 'fish').mean()) if len(new) else np.nan})
FA = pd.DataFrame(fate).set_index('area'); FA.to_csv('obis/out/fate.csv')
R['fate'] = {'still_single_median': float(FA.still_single.median()), 'new_single_median': float(FA.new_single.median()), 'new_ge10_median': float(FA.new_ge10.median()), 'new_per_1000rec_median': float(FA.new_per_1000rec.median()),
             'new_species_sum': int(FA.new_species.sum()), 'single99_sum': int(FA.single99.sum())}
# kesif orani: donemler arasi 1000 yeni tur-duzeyi kayit basina yeni tur
dr = {}
for t0, t1 in (('to1999', 'to2009'), ('to2009', 'to2019'), ('to2019', 'dated')):
    v = []
    for aid in ok.area:
        s0 = SP[(aid, t0)]; s1 = SP[(aid, t1)]; dn = s1.n.sum() - s0.n.sum()
        if dn >= 1000: v.append(1000 * (len(s1) - len(s0)) / dn)
    dr[t0 + '->' + t1] = {'n': len(v), 'median': float(np.median(v)), 'p25': float(np.percentile(v, 25)), 'p75': float(np.percentile(v, 75))}
R['discovery_rate'] = dr
# gruplarda geri-test
gb = {}
for g in ('fish', 'tetrapod', 'mollusc', 'arthropod', 'annelid', 'cnidarian', 'echinoderm', 'sponge'):
    q0 = T[(T.group == g) & (T.tag == 'to1999')].set_index('area'); q1 = T[(T.group == g) & (T.tag == 'dated')].set_index('area'); ix = q0.index[(q0.S >= 30)].intersection(q1.index)
    gb[g] = {'n': int(len(ix)), 'share_exceeds_chao': float((q1.S[ix] > q0.chao[ix]).mean()), 'median_S26_over_chao99': float((q1.S[ix] / q0.chao[ix]).median()), 'median_growth': float((q1.S[ix] / q0.S[ix]).median()), 'C99_median': float(q0.C[ix].median()), 'C26_median': float(q1.C[ix].median())}
R['backtest_groups'] = gb
# Turkiye'yi cevreleyen denizler
tr = cur[cur.name.str.contains('Black Sea|Marmara|Aegean|Mediterranean|Azov|Levant', case=False, regex=True)]
R['turkiye_seas'] = json.loads(tr[['name', 'type', 'S', 'N', 'f1', 'f2', 'chao', 'lo', 'hi', 'C', 'cov']].to_json(orient='records'))
R['turkiye_backtest'] = json.loads(bt[bt.name.str.contains('Black Sea|Marmara|Aegean|Mediterranean|Azov', case=False, regex=True)][['name', 'type', 'S99', 'chao99', 'hi99', 'S26', 'growth', 'pred_growth']].to_json(orient='records'))
def clean(o):
    if isinstance(o, dict): return {str(k): clean(v) for k, v in o.items()}
    if isinstance(o, (list, tuple)): return [clean(v) for v in o]
    if isinstance(o, (float, np.floating)): return None if (np.isnan(o) or np.isinf(o)) else float(o)
    if isinstance(o, np.integer): return int(o)
    return o
json.dump(clean(R), open('obis/out/results.json', 'w'), indent=1, ensure_ascii=False)
for k, v in clean(R).items(): print(k, json.dumps(v, ensure_ascii=False)[:1800])
'''

# ==============================================================================
# obis/obis_deep.py
# ==============================================================================
MODULES['obis/obis_deep.py'] = r'''
# P30 derinlestirme: dort kestirici, rastgele inceltme (null) geri-testi, ornek-disi kalibrasyon, yeni turlerin bilesimi
import json, glob, warnings, numpy as np, pandas as pd
warnings.filterwarnings('ignore'); rng = np.random.default_rng(300)
FISH = {'Teleostei', 'Elasmobranchii', 'Holocephali', 'Myxini', 'Petromyzonti', 'Coelacanthi', 'Chondrostei', 'Holostei', 'Cladistii', 'Dipneusti', 'Actinopteri', 'Actinopterygii'}; TETRA = {'Mammalia', 'Aves', 'Reptilia', 'Testudines', 'Squamata', 'Crocodylia'}
def group(ph, cl):
    if cl in FISH: return 'fish'
    if cl in TETRA: return 'tetrapod'
    return {'Mollusca': 'mollusc', 'Arthropoda': 'arthropod', 'Annelida': 'annelid', 'Cnidaria': 'cnidarian', 'Echinodermata': 'echinoderm', 'Porifera': 'sponge'}.get(ph, 'other')
def est(n):
    n = np.asarray(n); n = n[n > 0]; S = len(n); N = n.sum(); f = [int((n == k).sum()) for k in range(1, 11)]; f1, f2, f3, f4 = f[:4]
    chao1 = S + f1 * (f1 - 1) / (2.0 * (f2 + 1))
    ich = chao1 + (f3 / (4.0 * max(f4, 1))) * max(f1 - f2 * f3 / (2.0 * max(f4, 1)), 0)                      # Chiu vd. 2014 iChao1
    # ACE (nadir esigi 10)
    rare = n[n <= 10]; Sr = len(rare); Nr = rare.sum(); Sa = S - Sr
    if Nr > 0 and Nr != f1:
        Ca = 1 - f1 / Nr; g2 = max(Sr / Ca * sum(k * (k - 1) * f[k - 1] for k in range(1, 11)) / (Nr * (Nr - 1)) - 1, 0); ace = Sa + Sr / Ca + f1 / Ca * g2
    else: ace = np.nan
    jack1 = S + f1; jack2 = S + 2 * f1 - f2
    return {'S': S, 'N': int(N), 'f1': f1, 'f2': f2, 'chao1': chao1, 'ichao1': ich, 'ace': ace, 'jack1': jack1, 'jack2': jack2}
SP = {}; AR = {}
for f in sorted(glob.glob('obis/c/*.json')):
    d = json.load(open(f)); a = d['area']
    if d['tag'] not in ('to1999', 'to2009', 'to2019', 'dated'): continue
    df = pd.DataFrame(d['rows'], columns=['id', 'rank', 'n', 'phylum', 'cls', 'kingdom']); df = df[(df['rank'] == 'Species') & (df.kingdom == 'Animalia')]
    df = df.groupby('id', as_index=False).agg(n=('n', 'sum'), phylum=('phylum', 'first'), cls=('cls', 'first')); df['g'] = [group(p, c) for p, c in zip(df.phylum, df.cls)]
    SP[(a['id'], d['tag'])] = df.set_index('id'); AR[a['id']] = a
rows = []
for (aid, tag), df in SP.items():
    r = est(df.n.values); r.update(area=aid, tag=tag, name=AR[aid]['name'], type=AR[aid]['type']); rows.append(r)
T = pd.DataFrame(rows); T.to_csv('obis/out/estimators.csv', index=False)
P = T.pivot(index='area', columns='tag'); ok = P.index[P['S']['to1999'] >= 100]
E_ = ['chao1', 'ichao1', 'ace', 'jack1', 'jack2']; R = {'n_areas': int(len(ok))}
def bt(t0, t1='dated'):
    o = {}
    for e in E_:
        x = P[e][t0][ok]; y = P['S'][t1][ok]; s0 = P['S'][t0][ok]
        o[e] = {'share_exceeded': float((y > x).mean()), 'median_obs_over_est': float((y / x).median()), 'median_pred_growth': float((x / s0).median()), 'spearman': float((x / s0).rank().corr((y / s0).rank())),
                'median_claimed_completeness': float((s0 / x).median())}
    o['median_real_growth'] = float((P['S'][t1][ok] / P['S'][t0][ok]).median()); return o
R['backtest_1999_2026'] = bt('to1999'); R['backtest_2009_2026'] = bt('to2009'); R['backtest_1999_2009'] = bt('to1999', 'to2009'); R['backtest_2009_2019'] = bt('to2009', 'to2019'); R['backtest_1999_2019'] = bt('to1999', 'to2019')
# rastgele inceltme: 2026 tarihli envanteri 1999 kayit sayisina rastgele indir -> ayni kestirici 2026 turunu bulur mu?
thin = []
for aid in ok:
    full = SP[(aid, 'dated')].n.values.astype(np.int64); N0 = int(SP[(aid, 'to1999')].n.sum()); S26 = len(full); S99 = len(SP[(aid, 'to1999')])
    if N0 >= full.sum(): continue
    res = {e: [] for e in E_}; Sth = []
    for _ in range(30):
        x = rng.multivariate_hypergeometric(full, N0) if full.sum() < 2e9 and len(full) < 60000 else rng.binomial(full, N0 / full.sum())
        r = est(x); Sth.append(r['S'])
        for e in E_: res[e].append(r[e])
    thin.append({'area': aid, 'S26': S26, 'S99': S99, 'S_thin': float(np.mean(Sth)), **{e + '_thin': float(np.nanmean(res[e])) for e in E_}, 'chao1_real99': float(P['chao1']['to1999'][aid])})
TH = pd.DataFrame(thin).set_index('area'); TH.to_csv('obis/out/thinning.csv')
R['thinning'] = {'n': int(len(TH)), 'median_S_thin_over_S99': float((TH.S_thin / TH.S99).median()), 'share_S_thin_gt_S99': float((TH.S_thin > TH.S99).mean()),
                 **{e: {'share_exceeded': float((TH.S26 > TH[e + '_thin']).mean()), 'median_S26_over_est': float((TH.S26 / TH[e + '_thin']).median())} for e in E_},
                 'real_chao1_share_exceeded_same_areas': float((TH.S26 > TH.chao1_real99).mean()), 'real_chao1_median_ratio_same_areas': float((TH.S26 / TH.chao1_real99).median()),
                 'median_real_over_thin_chao': float((TH.chao1_real99 / TH.chao1_thin).median())}
# tur devri: 1999'da olup 2026 tarihli kume icinde kalanlar hep kalir (kumulatif). Yeni turlerin grup bilesimi
comp = []; newab = []
for aid in ok:
    a = SP[(aid, 'to1999')]; b = SP[(aid, 'dated')]; new = b.index.difference(a.index)
    comp.append(pd.concat([b.g.reindex(new).value_counts(normalize=True).rename('new'), a.g.value_counts(normalize=True).rename('old')], axis=1).fillna(0).assign(area=aid))
    nn = b.n.reindex(new); newab.append({'area': aid, 'new': len(new), 'new_single': float((nn == 1).mean()) if len(new) else np.nan, 'new_ge10': float((nn >= 10).mean()) if len(new) else np.nan, 'new_ge100': float((nn >= 100).mean()) if len(new) else np.nan,
                                         'old_single_share': float((a.n == 1).mean())})
CM = pd.concat(comp).reset_index().rename(columns={'index': 'g'}); cm = CM.groupby('g')[['new', 'old']].median(); R['new_species_composition'] = json.loads(cm.to_json(orient='index'))
NA_ = pd.DataFrame(newab); R['new_species_abundance'] = {'new_single_median': float(NA_.new_single.median()), 'new_ge10_median': float(NA_.new_ge10.median()), 'new_ge100_median': float(NA_.new_ge100.median()), 'total_new': int(NA_.new.sum())}
# ornek-disi kalibrasyon (10 yillik ufuk): 1999->2009 ile kur, 2009->2019 ile sina
def xy(t0, t1):
    s0 = P['S'][t0][ok]; return pd.DataFrame({'g': np.log(P['S'][t1][ok] / s0), 'm': np.log(P['chao1'][t0][ok] / s0), 'r': np.log(P['N'][t1][ok] / P['N'][t0][ok]), 'f1': P['f1'][t0][ok] / s0}).replace([np.inf, -np.inf], np.nan).dropna()
tr_, te = xy('to1999', 'to2009'), xy('to2009', 'to2019')
def fit(cols): X = np.column_stack([np.ones(len(tr_))] + [tr_[c] for c in cols]); return np.linalg.lstsq(X, tr_.g, rcond=None)[0]
def score(beta, cols):
    X = np.column_stack([np.ones(len(te))] + [te[c] for c in cols]); pr = X @ beta; return {'r2_oos': float(1 - ((te.g - pr) ** 2).sum() / ((te.g - te.g.mean()) ** 2).sum()), 'mae_log': float((te.g - pr).abs().mean()), 'bias_log': float((pr - te.g).mean())}
R['calibration'] = {'train_n': int(len(tr_)), 'test_n': int(len(te)), 'naive_chao_as_forecast': {'r2_oos': float(1 - ((te.g - te.m) ** 2).sum() / ((te.g - te.g.mean()) ** 2).sum()), 'mae_log': float((te.g - te.m).abs().mean()), 'bias_log': float((te.m - te.g).mean())},
                    'm_only': {'beta': fit(['m']).tolist(), **score(fit(['m']), ['m'])}, 'm_plus_records': {'beta': fit(['m', 'r']).tolist(), **score(fit(['m', 'r']), ['m', 'r'])}, 'records_only': {'beta': fit(['r']).tolist(), **score(fit(['r']), ['r'])},
                    'mean_growth_only': {'mae_log': float((te.g - tr_.g.mean()).abs().mean())}}
# efor dagilimi: yeni kayitlarin ne kadari yeni veri kumelerinden? (istatistiklerden: veri kumesi sayisi)
ds = []
for f in glob.glob('obis/c/*_to1999.json'):
    a = json.load(open(f)); b = json.load(open(f.replace('to1999', 'dated'))); ds.append((a['area']['id'], a['stats'].get('datasets'), b['stats'].get('datasets')))
DS = pd.DataFrame(ds, columns=['area', 'd99', 'd26']).set_index('area').loc[ok]; R['datasets'] = {'median_1999': float(DS.d99.median()), 'median_2026': float(DS.d26.median()), 'median_ratio': float((DS.d26 / DS.d99).median())}
g = np.log(P['S']['dated'][ok] / P['S']['to1999'][ok]); R['datasets']['spearman_growth_vs_dataset_growth'] = float(g.rank().corr(np.log(DS.d26 / DS.d99).rank())); R['datasets']['spearman_growth_vs_record_growth'] = float(g.rank().corr(np.log(P['N']['dated'][ok] / P['N']['to1999'][ok]).rank()))
def clean(o):
    if isinstance(o, dict): return {str(k): clean(v) for k, v in o.items()}
    if isinstance(o, (list, tuple)): return [clean(v) for v in o]
    if isinstance(o, (float, np.floating)): return None if (np.isnan(o) or np.isinf(o)) else float(o)
    if isinstance(o, np.integer): return int(o)
    return o
json.dump(clean(R), open('obis/out/deep.json', 'w'), indent=1, ensure_ascii=False)
for k, v in clean(R).items(): print(k, json.dumps(v, ensure_ascii=False)[:1500])
'''

# ==============================================================================
# paper/fig30en.py
# ==============================================================================
MODULES['paper/fig30en.py'] = r'''
import sys; sys.path.insert(0, '.')
from figstyle import *
import json, numpy as np, pandas as pd
R = json.load(open('obis/out/results.json')); DP = json.load(open('obis/out/deep.json')); R['deep'] = DP; json.dump(R, open('paper/out/results.json', 'w'))
T = pd.read_csv('obis/out/completeness.csv'); B = pd.read_csv('obis/out/backtest.csv'); TH = pd.read_csv('obis/out/thinning.csv')
cur = T[(T.tag == 'all') & (T.group == 'all') & (T.S >= 100)]
fig, ax = plt.subplots(1, 2, figsize=(7.4, 3.1)); bins = np.linspace(45, 95, 21)
ax[0].hist([100 * cur[cur.type == 'lme'].C, 100 * cur[cur.type == 'iho'].C], bins=bins, stacked=True, color=[C['b'], C['lb']], label=['Large marine ecosystems (66)', 'IHO sea areas (100)'])
ax[0].set_xlabel('Estimated inventory completeness (%)'); ax[0].set_ylabel('Number of areas'); ax[0].set_title('(a) Chao1 completeness, all records'); ax[0].legend(fontsize=7)
ax[1].scatter(cur.N, 100 * cur.C, s=10, c=np.where(cur.type == 'lme', C['b'], C['lb']), lw=0); ax[1].set_xscale('log'); ax[1].set_xlabel('Species-level records'); ax[1].set_ylabel('Estimated completeness (%)'); ax[1].set_title('(b) Completeness and sampling effort')
fig.tight_layout(); fig.savefig('paper/fig/Figure_1.png'); plt.close()
fig, ax = plt.subplots(1, 2, figsize=(7.4, 3.3))
ax[0].scatter(B.chao99, B.S26, s=9, c=np.where(B.S26 > B.hi99, C['r'], np.where(B.S26 > B.chao99, C['o'], C['b'])), lw=0, zorder=3); l = [80, 40000]; ax[0].plot(l, l, color='k', lw=.7); ax[0].set_xscale('log'); ax[0].set_yscale('log')
ax[0].set_xlabel('Chao1 estimate from records collected to 1999'); ax[0].set_ylabel('Species observed by 2026'); ax[0].set_title('(a) Estimate versus later observation')
for c, lab in ((C['r'], 'above upper 95% bound'), (C['o'], 'above point estimate'), (C['b'], 'within estimate')): ax[0].scatter([], [], s=12, c=c, label=lab)
ax[0].legend(fontsize=7, loc='upper left')
es = ['chao1', 'ichao1', 'ace', 'jack1', 'jack2']; nm = ['Chao1', 'iChao1', 'ACE', 'Jack-1', 'Jack-2']; x = np.arange(5)
ax[1].bar(x - .2, [100 * DP['backtest_1999_2026'][e]['share_exceeded'] for e in es], .4, color=C['r'], label='Records as collected to 1999'); ax[1].bar(x + .2, [100 * DP['thinning'][e]['share_exceeded'] for e in es], .4, color=C['gr'], label='Random subsample of same size')
ax[1].set_xticks(x); ax[1].set_xticklabels(nm, fontsize=7.5); ax[1].set_ylabel('Areas where 2026 richness exceeds estimate (%)'); ax[1].set_ylim(0, 100); ax[1].set_title('(b) Exceedance by estimator'); ax[1].legend(fontsize=7)
fig.tight_layout(); fig.savefig('paper/fig/Figure_2.png'); plt.close()
fig, ax = plt.subplots(1, 2, figsize=(7.4, 3.2))
ax[0].scatter(100 * (B.pred_growth - 1), 100 * (B.growth - 1), s=9, c=C['b'], lw=0); l = [0, 160]; ax[0].plot(l, l, color='k', lw=.7); ax[0].set_xlim(0, 160); ax[0].set_ylim(0, 300)
ax[0].set_xlabel('Predicted undetected species (% of 1999 inventory)'); ax[0].set_ylabel('Realised increase 1999–2026 (%)'); ax[0].set_title('(a) Predicted and realised increase')
cal = DP['calibration']; ks = [('naive_chao_as_forecast', 'Chao1,\nuncalibrated'), ('m_only', 'Chao1,\nrecalibrated'), ('records_only', 'Record growth\nonly'), ('m_plus_records', 'Chao1 +\nrecord growth')]
ax[1].bar(range(4), [cal[k]['r2_oos'] for k, _ in ks], color=[C['gr'], C['lb'], C['lb'], C['b']]); ax[1].axhline(0, color='k', lw=.6); ax[1].set_xticks(range(4)); ax[1].set_xticklabels([l for _, l in ks], fontsize=7); ax[1].set_ylabel('Out-of-sample R$^2$ (2009–2019)'); ax[1].set_title('(b) Forecasting ten-year inventory growth')
fig.tight_layout(); fig.savefig('paper/fig/Figure_3.png'); plt.close()
G = R['groups']; BG = R['backtest_groups']; nm = {'fish': 'Fishes', 'tetrapod': 'Mammals, birds, reptiles', 'mollusc': 'Molluscs', 'arthropod': 'Arthropods', 'annelid': 'Annelids', 'cnidarian': 'Cnidarians', 'echinoderm': 'Echinoderms', 'sponge': 'Sponges'}
ks = sorted(nm, key=lambda k: G[k]['C_median']); fig, ax = plt.subplots(1, 2, figsize=(7.4, 3), sharey=True); y = np.arange(len(ks))
ax[0].barh(y, [100 * G[k]['C_median'] for k in ks], color=C['lb'], edgecolor=C['b']); ax[0].errorbar([100 * G[k]['C_median'] for k in ks], y, xerr=[[100 * (G[k]['C_median'] - G[k]['C_p10']) for k in ks], [100 * (G[k]['C_p90'] - G[k]['C_median']) for k in ks]], fmt='none', ecolor='k', lw=.7, capsize=2)
ax[0].set_yticks(y); ax[0].set_yticklabels([nm[k] for k in ks]); ax[0].set_xlabel('Completeness (%), median and 10–90%'); ax[0].set_title('(a) Completeness by group')
ax[1].barh(y, [100 * BG[k]['share_exceeds_chao'] for k in ks], color=C['o']); ax[1].set_xlabel('Areas exceeding the 1999 estimate (%)'); ax[1].set_title('(b) Back-test by group')
fig.tight_layout(); fig.savefig('paper/fig/Figure_4.png'); plt.close()
'''

# ==============================================================================
# obis/obis_thr.py
# ==============================================================================
MODULES['obis/obis_thr.py'] = r'''
# 1999'da Chao1 tamligi %70 ve uzeri olan envanterlerin sonraki durumu -> paper/out/results.json ('thr70')
import json, pandas as pd
B = pd.read_csv('obis/out/backtest.csv'); q = B[B.C99 >= 0.7]
x = {'n_ge70': int(len(q)), 'share_ge70': float(len(q) / len(B)), 'median_growth_ge70': float(q.growth.median()), 'share_exceeded_ge70': float((q.S26 > q.chao99).mean()), 'n_lt70': int((B.C99 < 0.7).sum()), 'median_growth_lt70': float(B[B.C99 < 0.7].growth.median())}
R = json.load(open('paper/out/results.json')); R['thr70'] = x; json.dump(R, open('paper/out/results.json', 'w')); print('thr70', x)
'''

# ==============================================================================
# figstyle.py
# ==============================================================================
MODULES['figstyle.py'] = r'''
import matplotlib; matplotlib.use('Agg')
import matplotlib.pyplot as plt
plt.rcParams.update({'font.family': 'DejaVu Sans', 'font.size': 9, 'axes.spines.top': False, 'axes.spines.right': False, 'axes.grid': True, 'grid.alpha': .25, 'grid.linewidth': .5,
                     'axes.titlesize': 9.5, 'axes.titleweight': 'bold', 'axes.titlelocation': 'left', 'legend.frameon': False, 'legend.fontsize': 8, 'savefig.dpi': 300, 'savefig.bbox': 'tight', 'figure.dpi': 110})
C = {'b': '#1f5a94', 'r': '#c0392b', 'g': '#2e7d5b', 'o': '#d98200', 'k': '#333333', 'p': '#6a4c93', 'lb': '#8fb8de', 'gr': '#9aa0a6'}
'''

# ==============================================================================
# dl_all.py
# ==============================================================================
MODULES['dl_all.py'] = r'''
# n5: bes makale icin ham veri indirme. Kullanim: dl_all.py tc|ice|mp|ices|obis
import sys, os, json, time, urllib.request, urllib.parse, concurrent.futures as cf
UA = {'User-Agent': 'gomd-research/1.0 (gomd.pirireis.edu.tr)'}
def get(url, timeout=120, tries=4):
    last = None
    for k in range(tries):
        try:
            with urllib.request.urlopen(urllib.request.Request(url, headers=UA), timeout=timeout) as r: return r.read()
        except Exception as e:
            last = e; time.sleep(2 + 3 * k)
    raise last
what = sys.argv[1]
if what == 'tc':
    u = 'https://www.ncei.noaa.gov/data/international-best-track-archive-for-climate-stewardship-ibtracs/v04r01/access/csv/ibtracs.since1980.list.v04r01.csv'
    b = get(u, 900); open('tc/ibtracs.since1980.csv.part', 'wb').write(b); os.replace('tc/ibtracs.since1980.csv.part', 'tc/ibtracs.since1980.csv'); print('tc', len(b))
elif what == 'ice':
    for h in 'NS':
        f = '%s_Sea_Ice_Index_Regional_Daily_Data_G02135_v4.0.xlsx' % h
        b = get('https://noaadata.apps.nsidc.org/NOAA/G02135/seaice_analysis/' + f, 600); open('ice/' + f, 'wb').write(b); print(f, len(b))
    for h, d in (('N', 'north'), ('S', 'south')):
        f = '%s_seaice_extent_daily_v4.0.csv' % h
        b = get('https://noaadata.apps.nsidc.org/NOAA/G02135/%s/daily/data/%s' % (d, f), 600); open('ice/' + f, 'wb').write(b); print(f, len(b))
elif what == 'mp':
    base = 'https://services2.arcgis.com/C8EMgrsFcRFL6LrL/arcgis/rest/services/Marine_Microplastics_WGS84/FeatureServer/0/query'
    out = open('mp/ncei.jsonl.part', 'w', encoding='utf-8'); n = 0; off = 0
    while True:
        q = urllib.parse.urlencode(dict(where='1=1', outFields='*', returnGeometry='false', orderByFields='OBJECTID', resultOffset=off, resultRecordCount=2000, f='json'))
        d = json.loads(get(base + '?' + q)); fs = d.get('features', [])
        for x in fs: out.write(json.dumps(x['attributes'], ensure_ascii=False) + '\n')
        n += len(fs); off += len(fs)
        if len(fs) < 2000 and not d.get('exceededTransferLimit'): break
        if not fs: break
    out.close(); os.replace('mp/ncei.jsonl.part', 'mp/ncei.jsonl'); print('mp', n)
elif what == 'ices':
    A = 'https://sag.ices.dk/SAG_API/api/'
    lists = {}
    for y in range(2014, 2027):
        try: lists[y] = json.loads(get(A + 'StockList?year=%d' % y))
        except Exception as e: print('liste', y, e); lists[y] = []
        print(y, len(lists[y]), flush=True)
    json.dump(lists, open('ices/stocklists.json', 'w'))
    keys = sorted({s['AssessmentKey'] for y in lists for s in lists[y]})
    def one(k):
        f = 'ices/a/%d.json' % k
        if os.path.isfile(f): return 0
        try:
            s = json.loads(get(A + 'SummaryTable?assessmentKey=%d' % k, 90)); r = json.loads(get(A + 'FishStockReferencePoints?assessmentKey=%d' % k, 90))
            json.dump({'s': s, 'r': r}, open(f + '.part', 'w')); os.replace(f + '.part', f); return 1
        except Exception as e: return -1
    os.makedirs('ices/a', exist_ok=True)
    with cf.ThreadPoolExecutor(6) as ex: res = list(ex.map(one, keys))
    print('ices anahtar', len(keys), 'indirilen', res.count(1), 'hata', res.count(-1))
elif what == 'obis':
    A = 'https://api.obis.org/v3/'
    areas = [a for a in json.load(open('raw/obis_areas.json'))['results'] if a['type'] in ('lme', 'iho')]
    os.makedirs('obis/c', exist_ok=True)
    def one(arg):
        a, tag, extra = arg; f = 'obis/c/%s_%s.json' % (a['id'], tag)
        if os.path.isfile(f): return 0
        try:
            rows = []; skip = 0
            while True:
                d = json.loads(get(A + 'checklist?areaid=%s&size=10000&skip=%d%s' % (a['id'], skip, extra), 300))
                rs = d.get('results', [])
                rows += [(r.get('taxonID'), r.get('taxonRank'), r.get('records'), r.get('phylum'), r.get('class'), r.get('kingdom')) for r in rs]
                skip += len(rs)
                if len(rs) < 10000 or skip >= d.get('total', 0): break
            st = json.loads(get(A + 'statistics?areaid=%s%s' % (a['id'], extra), 120))
            json.dump({'area': a, 'tag': tag, 'stats': st, 'rows': rows}, open(f + '.part', 'w')); os.replace(f + '.part', f); return 1
        except Exception as e: return -1
    jobs = [(a, tag, extra) for a in areas for tag, extra in (('dated', '&enddate=2026-12-31'), ('all', ''), ('to1999', '&enddate=1999-12-31'), ('to2009', '&enddate=2009-12-31'), ('to2019', '&enddate=2019-12-31'))]
    with cf.ThreadPoolExecutor(4) as ex: res = list(ex.map(one, jobs))
    print('obis is', len(jobs), 'indirilen', res.count(1), 'hata', res.count(-1))
'''

if __name__ == "__main__":
    for _d in DIRS: os.makedirs(_d, exist_ok=True)
    for _m in MODULES:
        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, env=dict(os.environ, MPLBACKEND='Agg'))   # her modul ayri surecte calisir / each module runs in its own process
    print('Bitti / done. Sonuclar / results: paper/out/results.json')