# -*- coding: utf-8 -*-
# ------------------------------------------------------------------------------
# Küresel bir deniz mikroplastik arşivinde birim hataları, yöntem sınıfları ve eğilim tahminleri
# Unit errors, method classes and trend estimates in a global archive of marine microplastic concentrations
#
# GOMD Yayınlar · Makale 27 — https://gomd.pirireis.edu.tr/publications.php?id=rp27
#
# 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 rp27_microplastic_archive_units_methods_trends.py
#
# Microplastic data: NOAA National Centers for Environmental Information (NCEI), Marine Microplastics database (Nyadjro et al. 2023), https://www.ncei.noaa.gov/products/microplastics ; all records retrieved from the public map service on 14 November 2026. The original publications of the contributing studies should be cited when individual records are used.
# 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 = ['mp/mp_analysis.py', 'mp/mp_deep.py', 'paper/fig28en.py', 'mp/mp_flag.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 = ['mp', 'mp/out', 'paper', 'paper/out', 'paper/fig']

# ==============================================================================
# mp/mp_analysis.py
# ==============================================================================
MODULES['mp/mp_analysis.py'] = r'''
# P28: NCEI deniz mikroplastik arsivi - olcumler ne kadar karsilastirilabilir?
import json, itertools, warnings, numpy as np, pandas as pd
warnings.filterwarnings('ignore')
rng = np.random.default_rng(28)
d = pd.DataFrame([json.loads(l) for l in open('mp/ncei.jsonl')])
R = {'n_all': int(len(d)), 'medium': d.Medium.value_counts().to_dict(), 'unit': d.Unit.value_counts().to_dict(),
     'n_ref': int(d.Short_Reference.nunique()), 'n_org': int(d.ORGANIZATION.nunique())}
w = d[(d.Medium == 'Ocean water') & (d.Unit == 'pieces/m3')].copy()
w['date'] = pd.to_datetime(w.Date_m_d_yyyy, unit='ms'); w['yr'] = w.date.dt.year; w['mon'] = w.date.dt.month
w['c'] = w.Microplastics_measurement.astype(float)
w = w.dropna(subset=['c', 'Latitude__degree_', 'Longitude_degree_', 'yr']).rename(columns={'Latitude__degree_': 'lat', 'Longitude_degree_': 'lon'})
m = w.Sampling_Method.fillna('unknown').str.lower().str.strip()
cls = np.where(m.str.contains('neuston|manta|avani|plankton'), 'net', np.where(m.str.contains('grab|jar'), 'grab',
      np.where(m.str.contains('pump|intake'), 'pump', 'bulk')))   # bulk: kova, Niskin, Van Dorn, CTD rozet, ornekleyici
w['cls'] = cls; w['meth'] = m.replace({'stainless steel bucket': 'steel bucket', 'stainless-steel sampler': 'steel bucket'})
w['mesh'] = w.Mesh_size__mm_.round(3)
w['proto'] = w.meth + '|' + w.mesh.fillna(-1).astype(str)
w['study'] = w.Short_Reference.fillna(w.ORGANIZATION).str.strip()
w['cell'] = (np.floor(w.lat / 10) * 10).astype(int).astype(str) + '_' + (np.floor(w.lon / 10) * 10).astype(int).astype(str)
w['cell5'] = (np.floor(w.lat / 5) * 5).astype(int).astype(str) + '_' + (np.floor(w.lon / 5) * 5).astype(int).astype(str)
R['water'] = {'n': int(len(w)), 'studies': int(w.study.nunique()), 'orgs': int(w.ORGANIZATION.nunique()), 'protocols': int(w.proto.nunique()), 'methods': int(w.meth.nunique()),
              'years': [int(w.yr.min()), int(w.yr.max())], 'zero_share': float((w.c == 0).mean()), 'cells10': int(w.cell.nunique()),
              'ocean': w.Location_Oceans.value_counts().to_dict(), 'mesh_missing': float(w.mesh.isna().mean()),
              'q': {str(q): float(w.c.quantile(q)) for q in (0.05, 0.25, 0.5, 0.75, 0.95, 0.99)}, 'max': float(w.c.max()),
              'top_org_share': float(w.ORGANIZATION.value_counts(normalize=True).iloc[0]), 'top_org': w.ORGANIZATION.value_counts().index[0]}
pos = w[w.c > 0].copy(); pos['y'] = np.log10(pos.c)
R['water']['pos_n'] = int(len(pos)); R['water']['orders_of_magnitude_p1_p99'] = float(pos.y.quantile(0.99) - pos.y.quantile(0.01))
# 1) yontem sinifi / yontem ozeti
def summ(g):
    p = g[g.c > 0]
    return pd.Series({'n': len(g), 'studies': g.study.nunique(), 'zero': (g.c == 0).mean(), 'median_all': g.c.median(),
                      'gm_pos': 10 ** np.log10(p.c).mean() if len(p) else np.nan, 'p10': p.c.quantile(.1) if len(p) else np.nan, 'p90': p.c.quantile(.9) if len(p) else np.nan,
                      'mesh_med': g.mesh.median()})
bycls = w.groupby('cls').apply(summ); bymeth = w.groupby('meth').apply(summ).sort_values('n', ascending=False)
bycls.to_csv('mp/out/by_class.csv'); bymeth.to_csv('mp/out/by_method.csv')
R['by_class'] = json.loads(bycls.to_json(orient='index')); R['by_method'] = json.loads(bymeth.to_json(orient='index'))
# ag gozu egimi (yalniz ag ornekleri, calisma duzeyinde)
net = pos[(pos.cls == 'net') & pos.mesh.notna()]
sm = net.groupby('study').agg(y=('y', 'mean'), mesh=('mesh', 'median'), n=('y', 'size'))
sm = sm[sm.n >= 5]; x = np.log10(sm.mesh); b = np.polyfit(x, sm.y, 1)
bs = []
for _ in range(4000):
    i = rng.integers(0, len(sm), len(sm)); 
    if np.ptp(x.values[i]) > 0: bs.append(np.polyfit(x.values[i], sm.y.values[i], 1)[0])
R['mesh_slope'] = {'n_studies': int(len(sm)), 'slope_log_log': float(b[0]), 'ci': [float(np.percentile(bs, 2.5)), float(np.percentile(bs, 97.5))],
                   'r': float(np.corrcoef(x, sm.y)[0, 1]), 'mesh_range': [float(sm.mesh.min()), float(sm.mesh.max())]}
# 2) varyans ayristirma (LMG / Shapley): protokol, calisma, konum (10 derece), yil, ay
def r2(df, facs):
    if not facs: return 0.0
    X = [np.ones((len(df), 1))]
    for f in facs:
        X.append(pd.get_dummies(df[f], drop_first=True).values.astype(float))
    X = np.hstack(X); y = df.y.values
    beta, *_ = np.linalg.lstsq(X, y, rcond=None); res = y - X @ beta
    return 1 - res.var() / y.var()
def lmg(df, facs):
    cache = {}
    def g(s):
        s = tuple(sorted(s))
        if s not in cache: cache[s] = r2(df, list(s))
        return cache[s]
    sh = {f: 0.0 for f in facs}; perms = list(itertools.permutations(facs))
    for p in perms:
        for i, f in enumerate(p): sh[f] += g(p[:i + 1]) - g(p[:i])
    return {f: sh[f] / len(perms) for f in facs}, g(facs), {f: [g([f]), g(facs) - g([x for x in facs if x != f])] for f in facs}
pos['yrs'] = pos.yr.astype(str); pos['mons'] = pos.mon.astype(str)
F = ['proto', 'study', 'cell', 'yrs', 'mons']
sh, tot, fl = lmg(pos, F)
R['variance_all'] = {'n': int(len(pos)), 'total_r2': float(tot), 'lmg': {k: float(v) for k, v in sh.items()}, 'first_last': {k: [float(a), float(b_)] for k, (a, b_) in fl.items()},
                     'sd_log10': float(pos.y.std())}
pn = pos[pos.cls == 'net']
sh, tot, fl = lmg(pn, F)
R['variance_net'] = {'n': int(len(pn)), 'total_r2': float(tot), 'lmg': {k: float(v) for k, v in sh.items()}, 'first_last': {k: [float(a), float(b_)] for k, (a, b_) in fl.items()}, 'sd_log10': float(pn.y.std())}
# tek program (SEA) icinde konum+yil+ay
sea = pos[pos.ORGANIZATION == 'Sea Education Association']
sh, tot, fl = lmg(sea, ['cell', 'yrs', 'mons'])
R['variance_sea'] = {'n': int(len(sea)), 'total_r2': float(tot), 'lmg': {k: float(v) for k, v in sh.items()}, 'sd_log10': float(sea.y.std())}
# 3) es-konumlu ciftler: ayni 1 derece hucresi, ayni mevsim (+-1 ay), +-3 yil; calisma ici ve calismalar arasi
w['c1'] = np.floor(w.lat).astype(int).astype(str) + '_' + np.floor(w.lon).astype(int).astype(str)
pp = w[w.c > 0].copy(); pp['y'] = np.log10(pp.c)
def pairs(df, maxper=40):
    out = []
    for c1, g in df.groupby('c1'):
        if len(g) < 2: continue
        if len(g) > maxper: g = g.sample(maxper, random_state=1)
        a = g[['y', 'yr', 'mon', 'study', 'cls', 'proto']].values
        for i in range(len(a)):
            for j in range(i + 1, len(a)):
                if abs(a[i][1] - a[j][1]) <= 3 and min(abs(a[i][2] - a[j][2]), 12 - abs(a[i][2] - a[j][2])) <= 1:
                    out.append((abs(a[i][0] - a[j][0]), a[i][3] == a[j][3], a[i][4] == a[j][4], a[i][5] == a[j][5], a[i][4], a[j][4], c1))
    return pd.DataFrame(out, columns=['dy', 'same_study', 'same_cls', 'same_proto', 'c_a', 'c_b', 'c1'])
PR = pairs(pp)
def ps(x):
    return {'n': int(len(x)), 'cells': int(x.c1.nunique()), 'median_factor': float(10 ** x.dy.median()), 'p90_factor': float(10 ** x.dy.quantile(.9)), 'share_gt10x': float((x.dy > 1).mean()), 'share_lt2x': float((x.dy < np.log10(2)).mean())} if len(x) else {'n': 0}
R['pairs'] = {'within_study': ps(PR[PR.same_study]), 'between_study_same_protocol': ps(PR[~PR.same_study & PR.same_proto]),
              'between_study_same_class': ps(PR[~PR.same_study & PR.same_cls]), 'between_study_net_net': ps(PR[~PR.same_study & (PR.c_a == 'net') & (PR.c_b == 'net')]),
              'between_study_diff_class': ps(PR[~PR.same_study & ~PR.same_cls]), 'between_study_all': ps(PR[~PR.same_study])}
# hucre-kume bootstrap: calismalar arasi / calisma ici medyan orani
def med_ratio(x):
    a = x[~x.same_study & (x.c_a == 'net') & (x.c_b == 'net')].dy; b_ = x[x.same_study & (x.c_a == 'net')].dy
    return 10 ** (a.median() - b_.median()) if len(a) and len(b_) else np.nan
cells = PR.c1.unique(); grp = {c: g for c, g in PR.groupby('c1')}; bs = []
for _ in range(500):
    s = rng.choice(cells, len(cells)); bs.append(med_ratio(pd.concat([grp[c] for c in s])))
R['pairs']['net_between_over_within'] = {'ratio': float(med_ratio(PR)), 'ci': [float(np.nanpercentile(bs, 2.5)), float(np.nanpercentile(bs, 97.5))]}
# 4) karsilastirilabilirlik agi: calismalar ayni 5 derece hucresini paylasiyorsa bagli (yalniz ayni sinif)
def components(df, key):
    st = sorted(df.study.unique()); par = {s: s for s in st}
    def f(a):
        while par[a] != a: par[a] = par[par[a]]; a = par[a]
        return a
    for c, g in df.groupby(key):
        u = g.study.unique()
        for s in u[1:]: par[f(s)] = f(u[0])
    comp = pd.Series({s: f(s) for s in st}); sz = df.study.map(comp).value_counts()
    alone = sum(1 for s in st if (comp == comp[s]).sum() == 1)
    return {'studies': len(st), 'components': int(comp.nunique()), 'largest_component_studies': int(comp.value_counts().iloc[0]), 'largest_component_record_share': float(sz.iloc[0] / len(df)), 'isolated_studies': int(alone)}
w['clscell5'] = w.cls + '#' + w.cell5; w['clscell10'] = w.cls + '#' + w.cell
R['network'] = {'any_method_5deg': components(w, 'cell5'), 'same_class_5deg': components(w, 'clscell5'), 'same_class_10deg': components(w, 'clscell10'), 'net_only_5deg': components(w[w.cls == 'net'], 'cell5')}
ov = w.groupby('cell5').study.nunique(); R['network']['cells5'] = int(len(ov)); R['network']['cells5_multi_study'] = int((ov >= 2).sum()); R['network']['cells5_one_study_share'] = float((ov == 1).mean())
# 5) egilim: havuzlanmis ve program-ici (Kuzey Atlantik, ag ornekleri)
def trend(df, fe=None, label=''):
    df = df[df.c > 0].copy(); df['y'] = np.log10(df.c); X = [np.ones(len(df)), (df.yr - 2000).values / 10.0]
    cols = []
    for f in (fe or []): D = pd.get_dummies(df[f], drop_first=True).values.astype(float); cols.append(D)
    X = np.column_stack(X + cols) if cols else np.column_stack(X)
    beta, *_ = np.linalg.lstsq(X, df.y.values, rcond=None); res = df.y.values - X @ beta
    # yil-kumeli bootstrap
    ys = df.yr.unique(); idx = {y_: np.where(df.yr.values == y_)[0] for y_ in ys}; bs = []
    for _ in range(300):
        ii = np.concatenate([idx[y_] for y_ in rng.choice(ys, len(ys))])
        try: bs.append(np.linalg.lstsq(X[ii], df.y.values[ii], rcond=None)[0][1])
        except Exception: pass
    return {'n': int(len(df)), 'years': [int(df.yr.min()), int(df.yr.max())], 'slope_log10_per_decade': float(beta[1]), 'factor_per_decade': float(10 ** beta[1]),
            'ci_factor': [float(10 ** np.percentile(bs, 2.5)), float(10 ** np.percentile(bs, 97.5))]}
na = w[(w.Location_Oceans == 'Atlantic Ocean') & (w.lat > 0) & (w.lon < -10) & (w.cls == 'net')]
allnet = w[w.cls == 'net']
R['trend'] = {'global_all_methods_pooled': trend(w), 'global_net_pooled': trend(allnet), 'global_net_study_fe': trend(allnet, ['study']), 'global_net_cell_fe': trend(allnet, ['cell']),
              'global_net_study_cell_fe': trend(allnet, ['study', 'cell']), 'natl_net_pooled': trend(na), 'natl_net_cell_fe': trend(na, ['cell']), 'natl_net_study_cell_fe': trend(na, ['study', 'cell']),
              'sea_only_natl_cell_month_fe': trend(na[na.ORGANIZATION == 'Sea Education Association'], ['cell', 'mon']),
              'global_all_study_cell_fe': trend(w, ['study', 'cell'])}
# sifir paylari zaman icinde (SEA)
s2 = w[w.ORGANIZATION == 'Sea Education Association']; s2 = s2.assign(dec=(s2.yr // 10) * 10)
R['sea_decade'] = json.loads(s2.groupby('dec').agg(n=('c', 'size'), zero=('c', lambda x: (x == 0).mean()), median=('c', 'median'), p90=('c', lambda x: x.quantile(.9))).to_json(orient='index'))
# bilesim: on yila gore yontem sinifi ve calisma sayisi
w['dec'] = (w.yr // 10) * 10
comp = w.groupby('dec').agg(n=('c', 'size'), studies=('study', 'nunique'), net=('cls', lambda x: (x == 'net').mean()), median=('c', 'median'), gm_pos=('c', lambda x: 10 ** np.log10(x[x > 0]).mean()), protocols=('proto', 'nunique'))
comp.to_csv('mp/out/by_decade.csv'); R['by_decade'] = json.loads(comp.to_json(orient='index'))
# 6) uyumlastirilmis indeks: log10(c) = calisma etkisi + hucre(10) etkisi + yil ; hucre etkisi = protokolden arindirilmis seviye (yalniz ag, en buyuk bagli bilesen)
def harmon(df):
    df = df[df.c > 0].copy(); df['y'] = np.log10(df.c)
    S = pd.get_dummies(df.study, drop_first=True).astype(float); C = pd.get_dummies(df.cell).astype(float)
    X = np.hstack([C.values, S.values]); beta, *_ = np.linalg.lstsq(X, df.y.values, rcond=None)
    ce = pd.Series(beta[:C.shape[1]], index=C.columns); raw = df.groupby('cell').y.mean(); n = df.groupby('cell').size(); ns = df.groupby('cell').study.nunique()
    ref = df.study.value_counts().index[0]
    return pd.DataFrame({'adj': ce, 'raw': raw, 'n': n, 'studies': ns}), ref, pd.Series(beta[C.shape[1]:], index=S.columns)
H, ref, se = harmon(allnet)
H = H[H.n >= 10]; H.to_csv('mp/out/harmonised_cells.csv')
R['harmonised'] = {'cells': int(len(H)), 'reference_study': ref, 'corr_raw_adj': float(np.corrcoef(H.raw, H.adj)[0, 1]), 'spearman': float(H.raw.rank().corr(H.adj.rank())),
                   'median_abs_shift_factor': float(10 ** (H.adj - H.raw).abs().median()), 'share_shift_gt3x': float(((H.adj - H.raw).abs() > np.log10(3)).mean()),
                   'top10_raw_in_top10_adj': int(len(set(H.raw.nlargest(10).index) & set(H.adj.nlargest(10).index))),
                   'study_effect_sd_log10': float(se.std()), 'study_effect_range_factor': float(10 ** (se.max() - se.min())), 'n_study_effects': int(len(se))}
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('mp/out/results.json', 'w'), indent=1, ensure_ascii=False)
for k in ('water', 'mesh_slope', 'variance_all', 'variance_net', 'variance_sea', 'pairs', 'network', 'trend', 'sea_decade', 'by_decade', 'harmonised', 'by_class'): print(k, json.dumps(clean(R[k]), ensure_ascii=False))
'''

# ==============================================================================
# mp/mp_deep.py
# ==============================================================================
MODULES['mp/mp_deep.py'] = r'''
# P28 derinlestirme: birim hatasi sinamasi, engel (hurdle) + PPML egilimleri, birini-disarida-birak, cift tasarimi duyarliligi
import json, warnings, numpy as np, pandas as pd, statsmodels.api as sm
warnings.filterwarnings('ignore'); rng = np.random.default_rng(280)
d = pd.DataFrame([json.loads(l) for l in open('mp/ncei.jsonl')])
w = d[(d.Medium == 'Ocean water') & (d.Unit == 'pieces/m3')].copy()
w['date'] = pd.to_datetime(w.Date_m_d_yyyy, unit='ms'); w['yr'] = w.date.dt.year; w['mon'] = w.date.dt.month; w['c'] = w.Microplastics_measurement.astype(float)
w = w.rename(columns={'Latitude__degree_': 'lat', 'Longitude_degree_': 'lon'}).dropna(subset=['c', 'lat', 'lon', 'yr'])
m = w.Sampling_Method.fillna('unknown').str.lower().str.strip()
w['cls'] = np.where(m.str.contains('neuston|manta|avani|plankton'), 'net', np.where(m.str.contains('grab|jar'), 'grab', np.where(m.str.contains('pump|intake'), 'pump', 'bulk')))
w['study'] = w.Short_Reference.fillna(w.ORGANIZATION).str.strip(); w['cell'] = (np.floor(w.lat / 10) * 10).astype(int).astype(str) + '_' + (np.floor(w.lon / 10) * 10).astype(int).astype(str)
R = {}
# 1) birim sinamasi: Barrows vd. 2018 (1 L elle ornek; yayimlanan ortalama 11.8 +- 24.0 adet/L, n = 1393)
g = w[w.study.str.contains('Barrows')]
p = g[g.c > 0].c
R['unit_check'] = {'n_archive': int(len(g)), 'n_published': 1393, 'archive_mean': float(g.c.mean()), 'archive_sd': float(g.c.std()), 'archive_mean_x1000': float(1000 * g.c.mean()), 'archive_sd_x1000': float(1000 * g.c.std()),
                   'published_mean_per_L': 11.8, 'published_sd_per_L': 24.0, 'min_pos': float(p.min()), 'share_integer_thousandths': float(np.isclose((p * 1000) % 1, 0, atol=1e-6).mean() + np.isclose((p * 1000) % 1, 1, atol=1e-6).mean()),
                   'implied_volume_m3_if_unit_true': float(1 / p.min()), 'true_volume_m3': 0.001, 'correct_value_factor': 1e6, 'archive_median': float(g.c.median()), 'corrected_median_per_m3': float(g.c.median() * 1e6),
                   'class_text': g.Concentration_class_text.value_counts().to_dict(), 'oceans': g.Location_Oceans.value_counts().to_dict(), 'years': [int(g.yr.min()), int(g.yr.max())]}
# genel tutarlilik taramasi: calisma bazinda en kucuk pozitif deger -> ima edilen ornek hacmi
rows = []
for s, q in w.groupby('study'):
    pp = q[q.c > 0].c
    if len(q) < 10 or not len(pp): continue
    rows.append({'study': s, 'n': len(q), 'cls': q.cls.mode().iloc[0], 'method': q.Sampling_Method.mode().iloc[0], 'mesh': float(q.Mesh_size__mm_.median()), 'min_pos': float(pp.min()), 'implied_volume_m3': float(1 / pp.min()), 'median_pos': float(pp.median())})
V = pd.DataFrame(rows).sort_values('implied_volume_m3'); V.to_csv('mp/out/implied_volume.csv', index=False)
R['implied_volume'] = {c: {'studies': int(len(q)), 'median_m3': float(q.implied_volume_m3.median()), 'min_m3': float(q.implied_volume_m3.min()), 'max_m3': float(q.implied_volume_m3.max())} for c, q in V.groupby('cls')}
R['implied_volume_rows'] = json.loads(V.to_json(orient='records'))
# duzeltmenin etkisi: arsiv "sinif" etiketleri ve yontemler arasi fark
w2 = w.copy(); w2.loc[w2.study.str.contains('Barrows'), 'c'] *= 1e6
def gm(x): x = x[x > 0]; return float(10 ** np.log10(x).mean())
R['correction_effect'] = {'grab_gm_before': gm(w[w.cls == 'grab'].c), 'grab_gm_after': gm(w2[w2.cls == 'grab'].c), 'net_gm': gm(w[w.cls == 'net'].c), 'bulk_gm': gm(w[w.cls == 'bulk'].c), 'pump_gm': gm(w[w.cls == 'pump'].c),
                          'global_median_before': float(w.c.median()), 'global_median_after': float(w2.c.median()), 'global_mean_before': float(w.c.mean()), 'global_mean_after': float(w2.c.mean()),
                          'share_records_affected': float(w.study.str.contains('Barrows').mean())}
# 2) egilimler: ag ornekleri; (a) tespit olasiligi (dogrusal olasilik), (b) pozitif kisim log10, (c) PPML (sifirlar dahil duzey)
net = w[w.cls == 'net'].copy(); net['t'] = (net.yr - 2000) / 10.0; net['pos'] = (net.c > 0).astype(float)
def design(df, fe):
    X = [np.ones(len(df)), df.t.values]
    for f in fe: X.append(pd.get_dummies(df[f], drop_first=True).values.astype(float))
    return np.column_stack(X)
def boot_year(df, fn, nb=200):
    ys = df.yr.unique(); idx = {y: np.where(df.yr.values == y)[0] for y in ys}; out = []
    for _ in range(nb):
        ii = np.concatenate([idx[y] for y in rng.choice(ys, len(ys))])
        try: out.append(fn(df.iloc[ii]))
        except Exception: pass
    return [float(np.percentile(out, 2.5)), float(np.percentile(out, 97.5))]
def lpm(df, fe):
    X = design(df, fe); return float(np.linalg.lstsq(X, df.pos.values, rcond=None)[0][1])
def logpos(df, fe):
    q = df[df.c > 0]; X = design(q, fe); return float(np.linalg.lstsq(X, np.log10(q.c.values), rcond=None)[0][1])
def ppml(df, fe):
    X = design(df, fe).astype(np.float32).astype(float); keep = X.std(0) > 0; keep[0] = True; X = X[:, keep]; y = df.c.values
    beta = np.zeros(X.shape[1]); beta[0] = np.log(max(y.mean(), 1e-9))
    for _ in range(60):
        eta = np.clip(X @ beta, -30, 30); mu = np.exp(eta); z = eta + (y - mu) / mu; WX = X * mu[:, None]
        nb_ = np.linalg.solve(X.T @ WX + 1e-8 * np.eye(X.shape[1]), WX.T @ z)
        if np.max(np.abs(nb_ - beta)) < 1e-7: beta = nb_; break
        beta = nb_
    return float(beta[1])
def allthree(df, fe, nb=120):
    o = {'n': int(len(df)), 'zero_share': float(1 - df.pos.mean()), 'years': [int(df.yr.min()), int(df.yr.max())], 'studies': int(df.study.nunique())}
    b = lpm(df, fe); o['detect_pp_per_decade'] = 100 * b; o['detect_ci'] = [100 * v for v in boot_year(df, lambda q: lpm(q, fe), nb)]
    b = logpos(df, fe); o['pos_factor_per_decade'] = 10 ** b; o['pos_ci'] = [10 ** v for v in boot_year(df, lambda q: logpos(q, fe), nb)]
    try:
        b = ppml(df, fe); o['ppml_factor_per_decade'] = float(np.exp(b)); o['ppml_ci'] = [float(np.exp(v)) for v in boot_year(df, lambda q: ppml(q, fe), 60)] if (len(fe) < 2 or df.study.nunique() <= 2) else None
    except Exception as e: o['ppml_error'] = str(e)[:80]
    return o
sea = net[net.ORGANIZATION == 'Sea Education Association']; natl = sea[(sea.Location_Oceans == 'Atlantic Ocean') & (sea.lat > 0)]; pac = sea[sea.Location_Oceans == 'Pacific Ocean']
R['trend_models'] = {'global_pooled': allthree(net, []), 'global_study_fe': allthree(net, ['study']), 'global_study_cell_fe': allthree(net, ['study', 'cell']),
                     'sea_natl_cell_month': allthree(natl, ['cell', 'mon']), 'sea_pacific_cell_month': allthree(pac, ['cell', 'mon']) if len(pac) > 200 else None, 'sea_all_cell_month': allthree(sea, ['cell', 'mon'])}
R['sea_spans'] = {'natl': [int(natl.yr.min()), int(natl.yr.max()), int(len(natl))], 'pacific': [int(pac.yr.min()), int(pac.yr.max()), int(len(pac))] if len(pac) else None}
# 3) bilesim payi ve birini-disarida-birak
def slope(df, fe): q = df[df.c > 0]; X = design(q, fe); return float(np.linalg.lstsq(X, np.log10(q.c.values), rcond=None)[0][1])
sp, sw = slope(net, []), slope(net, ['study'])
R['composition'] = {'pooled_log10': sp, 'within_log10': sw, 'between_share': float(1 - sw / sp)}
loo = []
for s in net.study.unique():
    q = net[net.study != s]; loo.append({'study': s, 'n': int((net.study == s).sum()), 'pooled': 10 ** slope(q, []), 'within': 10 ** slope(q, ['study'])})
L = pd.DataFrame(loo); L.to_csv('mp/out/leave_one_out.csv', index=False)
R['leave_one_out'] = {'pooled_min': float(L.pooled.min()), 'pooled_max': float(L.pooled.max()), 'within_min': float(L.within.min()), 'within_max': float(L.within.max()), 'within_most_influential': L.loc[(L.within - 10 ** sw).abs().idxmax()].to_dict()}
# calisma surelerinin dagilimi: egilimi kim tanimliyor?
sp_ = net.groupby('study').agg(n=('c', 'size'), y0=('yr', 'min'), y1=('yr', 'max')); sp_['span'] = sp_.y1 - sp_.y0 + 1
R['study_spans'] = {'studies': int(len(sp_)), 'span_le2': int((sp_.span <= 2).sum()), 'span_ge5': int((sp_.span >= 5).sum()), 'span_ge10': int((sp_.span >= 10).sum()), 'records_in_span_ge10': float(sp_[sp_.span >= 10].n.sum() / sp_.n.sum()),
                    'post2010_studies': int((sp_.y0 >= 2010).sum()), 'post2010_span_ge5': int(((sp_.y0 >= 2010) & (sp_.span >= 5)).sum()), 'long': sp_[sp_.span >= 5].reset_index().values.tolist()}
# 4) cift tasarimi duyarliligi
pp_ = w[(w.c > 0) & (w.cls == 'net')].copy(); pp_['y'] = np.log10(pp_.c)
def pairs(df, deg, dyr, maxper=40):
    df = df.assign(k=(np.floor(df.lat / deg)).astype(int).astype(str) + '_' + (np.floor(df.lon / deg)).astype(int).astype(str)); W = []; B = []; cells = 0
    for k, q in df.groupby('k'):
        if len(q) < 2: continue
        if len(q) > maxper: q = q.sample(maxper, random_state=1)
        a = q[['y', 'yr', 'mon', 'study']].values; hit = False
        for i in range(len(a)):
            for j in range(i + 1, len(a)):
                if abs(a[i][1] - a[j][1]) <= dyr and min(abs(a[i][2] - a[j][2]), 12 - abs(a[i][2] - a[j][2])) <= 1:
                    (W if a[i][3] == a[j][3] else B).append(abs(a[i][0] - a[j][0])); hit = hit or a[i][3] != a[j][3]
        cells += hit
    return {'within_n': len(W), 'between_n': len(B), 'between_cells': int(cells), 'within_factor': float(10 ** np.median(W)) if W else None, 'between_factor': float(10 ** np.median(B)) if B else None, 'between_p90': float(10 ** np.percentile(B, 90)) if B else None}
R['pairs_sensitivity'] = {'%g_deg_%d_yr' % (dg, dy): pairs(pp_, dg, dy) for dg in (0.5, 1, 2) for dy in (1, 3, 5)}
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('mp/out/deep.json', 'w'), indent=1, ensure_ascii=False)
for k, v in clean(R).items():
    if k != 'implied_volume_rows': print(k, json.dumps(v, ensure_ascii=False)[:1600])
'''

# ==============================================================================
# paper/fig28en.py
# ==============================================================================
MODULES['paper/fig28en.py'] = r'''
import sys; sys.path.insert(0, '.')
from figstyle import *
import json, numpy as np, pandas as pd
R = json.load(open('mp/out/results.json')); DP = json.load(open('mp/out/deep.json')); EX = json.load(open('mp/out/extra.json'))
R['deep'] = DP; R['extra'] = {k: EX[k] for k in ('sigma_log10_within', 'n_per_group')}; json.dump(R, open('paper/out/results.json', 'w'))
d = pd.DataFrame([json.loads(l) for l in open('mp/ncei.jsonl')]); w = d[(d.Medium == 'Ocean water')].copy(); w['c'] = w.Microplastics_measurement.astype(float)
m = w.Sampling_Method.fillna('').str.lower(); w['cls'] = np.where(m.str.contains('neuston|manta|avani|plankton'), 'net', np.where(m.str.contains('grab|jar'), 'grab', np.where(m.str.contains('pump|intake'), 'pump', 'bulk')))
pos = w[w.c > 0]; names = {'net': 'Net tow', 'pump': 'Pump / intake', 'bulk': 'Bottle / bucket', 'grab': 'Grab sample'}; cols = {'net': C['b'], 'pump': C['g'], 'bulk': C['r'], 'grab': C['o']}
fig, ax = plt.subplots(1, 2, figsize=(7.4, 3.3)); bins = np.linspace(-4, 5, 46)
for k in ('net', 'grab', 'pump', 'bulk'):
    v = np.log10(pos[pos.cls == k].c); ax[0].hist(v, bins=bins, color=cols[k], alpha=.6, label='%s (n = %s)' % (names[k], format(len(v), ',')), density=True)
v = np.log10(pos[pos.cls == 'grab'].c * 1e6); ax[0].hist(v, bins=bins, histtype='step', color=C['o'], lw=1.3, ls='--', density=True, label='Grab sample, values × 10$^6$')
ax[0].set_xlabel('log$_{10}$ concentration as archived (pieces m$^{-3}$)'); ax[0].set_ylabel('Density'); ax[0].set_title('(a) Concentration by method class'); ax[0].legend(fontsize=6.5)
V = pd.DataFrame(DP['implied_volume_rows']); V = V.sort_values('implied_volume_m3'); jit = np.random.default_rng(3).uniform(-.25, .25, len(V)); order = ['bulk', 'pump', 'net', 'grab']
for i, k in enumerate(order):
    q = V[V.cls == k]; ax[1].scatter(q.implied_volume_m3, i + jit[:len(q)], s=10 + 3 * np.sqrt(q.n), color=cols[k], alpha=.8, lw=0)
ax[1].set_xscale('log'); ax[1].set_yticks(range(4)); ax[1].set_yticklabels([names[k] for k in order]); ax[1].set_xlabel('Implied sample volume, 1 / smallest non-zero value (m$^3$)'); ax[1].set_title('(b) Detection-limit consistency check')
ax[1].axvline(0.001, color='k', lw=.6, ls=':'); ax[1].text(0.0013, -0.42, '1 L', fontsize=7)
fig.tight_layout(); fig.savefig('paper/fig/Figure_1.png'); plt.close()
fig, ax = plt.subplots(1, 2, figsize=(7.4, 3.1))
lab = {'proto': 'Protocol', 'study': 'Study', 'cell': 'Location (10°)', 'yrs': 'Year', 'mons': 'Month'}; ks = ['proto', 'study', 'cell', 'yrs', 'mons']; x = np.arange(len(ks))
ax[0].bar(x - .2, [100 * R['variance_all']['lmg'][k] for k in ks], .4, color=C['b'], label='All methods'); ax[0].bar(x + .2, [100 * R['variance_net']['lmg'][k] for k in ks], .4, color=C['lb'], label='Net tows only')
ax[0].set_xticks(x); ax[0].set_xticklabels([lab[k] for k in ks], fontsize=7.5); ax[0].set_ylabel('Share of variance explained (%)'); ax[0].set_title('(a) Variance partitioning'); ax[0].legend()
P = R['pairs']; ks = [('within_study', 'Same study'), ('between_study_same_protocol', 'Different study,\nsame protocol'), ('between_study_net_net', 'Different study,\nnet vs net'), ('between_study_diff_class', 'Different\nmethod class')]
ax[1].bar(range(4), [P[k]['median_factor'] for k, _ in ks], color=[C['gr'], C['lb'], C['b'], C['r']]); ax[1].set_yscale('log'); ax[1].set_xticks(range(4)); ax[1].set_xticklabels([l for _, l in ks], fontsize=7)
for i, (k, _) in enumerate(ks): ax[1].text(i, P[k]['median_factor'] * 1.12, '×%.1f\nn = %s' % (P[k]['median_factor'], format(P[k]['n'], ',')), ha='center', fontsize=7)
ax[1].set_ylim(1, 900); ax[1].set_ylabel('Median fold difference within pairs'); ax[1].set_title('(b) Co-located sample pairs')
fig.tight_layout(); fig.savefig('paper/fig/Figure_2.png'); plt.close()
T = DP['trend_models']; ks = [('global_pooled', 'All net tows, pooled'), ('global_study_fe', '+ study effects'), ('global_study_cell_fe', '+ study and location effects'), ('sea_natl_cell_month', 'Single programme, N Atlantic\n(1986–2008)'), ('sea_all_cell_month', 'Single programme, all cruises\n(1986–2012)')]
fig, ax = plt.subplots(1, 2, figsize=(7.4, 3.1), sharey=True); y = np.arange(len(ks))[::-1]
for a, key, cik, ttl in ((ax[0], 'pos_factor_per_decade', 'pos_ci', '(a) Geometric mean, non-zero tows'), (ax[1], 'ppml_factor_per_decade', 'ppml_ci', '(b) Arithmetic mean, all tows')):
    for i, (k, l) in zip(y, ks):
        t = T[k]
        if t.get(cik): a.plot(t[cik], [i, i], color=C['b'], lw=1.5)
        a.plot(t[key], i, 'o', color=C['b'], mfc=(C['b'] if t.get(cik) else 'white'))
    a.axvline(1, color='k', lw=.7); a.set_xscale('log'); a.set_xticks([0.8, 1, 1.5, 2, 3, 4, 6]); a.set_xticklabels(['0.8', '1', '1.5', '2', '3', '4', '6']); a.set_xlim(0.7, 7); a.set_title(ttl); a.set_xlabel('Fold change per decade')
ax[0].set_yticks(y); ax[0].set_yticklabels([l for _, l in ks], fontsize=7.5)
fig.tight_layout(); fig.savefig('paper/fig/Figure_3.png'); plt.close()
D = R['by_decade']; dec = sorted(D); fig, ax = plt.subplots(1, 2, figsize=(7.4, 3))
ax[0].bar([int(x) for x in dec], [D[x]['studies'] for x in dec], width=7, color=C['lb'], edgecolor=C['b']); ax[0].set_ylabel('Number of studies', color=C['b']); a2 = ax[0].twinx(); a2.plot([int(x) for x in dec], [D[x]['protocols'] for x in dec], color=C['r'], marker='o'); a2.set_ylabel('Number of protocols', color=C['r']); a2.spines['right'].set_visible(True); a2.grid(False)
ax[0].set_xlabel('Decade'); ax[0].set_title('(a) Composition of the archive')
L = pd.read_csv('mp/out/leave_one_out.csv').sort_values('within'); ax[1].scatter(L.pooled, L.within, s=8 + 2 * np.sqrt(L.n), color=C['b'], alpha=.7, lw=0); ax[1].axhline(1, color='k', lw=.6); ax[1].set_xlabel('Pooled trend (fold per decade)'); ax[1].set_ylabel('Within-study trend (fold per decade)'); ax[1].set_title('(b) Leave-one-study-out')
r = L.loc[L.within.idxmax()]; ax[1].annotate('without the\nlongest programme', (r.pooled, r.within), fontsize=7, xytext=(-75, -28), textcoords='offset points', arrowprops=dict(arrowstyle='-', lw=.5))
fig.tight_layout(); fig.savefig('paper/fig/Figure_4.png'); plt.close()
'''

# ==============================================================================
# mp/mp_flag.py
# ==============================================================================
MODULES['mp/mp_flag.py'] = r'''
# Birim tutarsizligi icin isaretlenen iki calismanin (Barrows vd. 2018; Dehm vd. 2020) ozet sayilari -> paper/out/results.json ('flag')
import json, numpy as np, pandas as pd
d = pd.DataFrame([json.loads(l) for l in open('mp/ncei.jsonl')]); w = d[d.Medium == 'Ocean water'].copy(); w['c'] = w.Microplastics_measurement.astype(float); sr = w.Short_Reference.fillna('')
bar = w[sr.str.contains('Barrows')]; de = w[sr.str.contains('Dehm')]; c = de.c; pp = c[c > 0]; yr = pd.to_datetime(de.Date_m_d_yyyy, unit='ms').dt.year
dehm = {'n': int(len(de)), 'mean': float(c.mean()), 'median': float(c.median()), 'min_pos': float(pp.min()), 'max': float(c.max()), 'mean_x1000': float(1000 * c.mean()), 'published_per_L': [1.6, 2.0], 'zero': float((c == 0).mean()),
        'unit': de.Unit.value_counts().to_dict(), 'class': de.Concentration_class_text.value_counts().to_dict(), 'years': [int(yr.min()), int(yr.max())]}
json.dump(dehm, open('mp/out/dehm.json', 'w'))
x = {'barrows_nonzero': int((bar.c > 0).sum()), 'dehm_nonzero': int((de.c > 0).sum()), 'very_high_now': int((w.Concentration_class_text == 'Very High').sum()), 'flagged_records': int(len(bar) + len(de)), 'flagged_share_water': float((len(bar) + len(de)) / len(w)),
     'barrows_class': bar.Concentration_class_text.value_counts().to_dict(), 'water_n': int(len(w)), 'barrows_min_corrected': float(bar.c[bar.c > 0].min() * 1e6), 'dehm_min_corrected': float(de.c[de.c > 0].min() * 1e6),
     'p99_before': float(w.c.quantile(.99)), 'share_ge10_before': float((w.c >= 10).mean())}
w2 = w.copy(); w2.loc[sr.str.contains('Barrows|Dehm'), 'c'] *= 1e6
x.update({'share_ge10_after': float((w2.c >= 10).mean()), 'median_after': float(w2.c.median()), 'mean_after': float(w2.c.mean()), 'mean_before': float(w.c.mean()), 'median_before': float(w.c.median()),
          'gm_pos_before': float(10 ** np.log10(w.c[w.c > 0]).mean()), 'gm_pos_after': float(10 ** np.log10(w2.c[w2.c > 0]).mean())})
x['dehm'] = dehm
R = json.load(open('paper/out/results.json')); R['flag'] = x; json.dump(R, open('paper/out/results.json', 'w'))
print('flag', x['flagged_records'], round(x['flagged_share_water'], 4))
'''

# ==============================================================================
# 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')