#!/usr/bin/env python3
"""
Reproducible structured characterization of malaria diagnostic
Target Product Profile (TPP) requirements across two use-cases.

Source dataset: "Target product profiles for malaria diagnostics"
(supplementary Table 1 from the malERA consultative group on diagnoses
and diagnostics), DOI 10.1371/journal.pmed.1000396.t001, licence CC BY 4.0.

Re-running this script reproduces every figure, table, and number in
results.json. All randomness is seeded.

Pipeline:
  1. Acquire Table_1.xls: attempt canonical figshare download + md5 verify,
     fall back to a local md5-gated copy (hard-fail on mismatch).
  2. Parse: strip HTML, split section headers from substantive criteria,
     extract per-cell priority-tier tags (E/D/O).
  3. Pre-registered Test 1: priority-tier composition x use-case
     (Freeman-Halton exact; Cramer's V; chi-square for reference) + sensitivity.
  4. Pre-registered Test 2: cross-use-case concordance per criterion
     (exact binomial vs 0.5; Wilson CI) + sensitivity.
  5. Descriptive requirement-map (most-stringent tier per criterion x use-case).
  6. Emit figures/, tables/ (+ Frictionless datapackage.json) and results.json.
"""
from __future__ import annotations
import os, re, sys, json, hashlib, collections, datetime
from math import exp as mexp
import numpy as np
import pandas as pd
from scipy import stats
from scipy.special import gammaln
import matplotlib as mpl
mpl.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Patch

# ----------------------------------------------------------------------------
# Reproducibility
# ----------------------------------------------------------------------------
SEED = 20240517
np.random.seed(SEED)

HERE = os.path.dirname(os.path.abspath(__file__))
FIG_DIR = os.path.join(HERE, "figures")
TBL_DIR = os.path.join(HERE, "tables")
os.makedirs(FIG_DIR, exist_ok=True)
os.makedirs(TBL_DIR, exist_ok=True)

# ----------------------------------------------------------------------------
# Data source constants
# ----------------------------------------------------------------------------
CANONICAL_URL = "https://ndownloader.figshare.com/files/803775"
EXPECTED_MD5 = "2bbf6174e03e39ce858e4fcdc0ad287f"
LOCAL_FALLBACK = "/Users/markhahnel/datasetpapers-handfeed/malaria/Table_1.xls"
SOURCE_DOI = "10.1371/journal.pmed.1000396.t001"
SOURCE_LICENCE = "CC BY 4.0"


def _md5(path: str) -> str:
    h = hashlib.md5()
    with open(path, "rb") as fh:
        for chunk in iter(lambda: fh.read(65536), b""):
            h.update(chunk)
    return h.hexdigest()


def acquire_data() -> tuple[str, str]:
    """Return (path, provenance) for an md5-verified Table_1.xls.

    Attempts the canonical figshare download first; on any failure or md5
    mismatch, falls back to the local pinned copy. The local copy is md5-gated:
    a mismatch is a hard failure.
    """
    dest = os.path.join(HERE, "Table_1.xls")
    # 1. Attempt canonical download.
    try:
        import urllib.request
        with urllib.request.urlopen(CANONICAL_URL, timeout=30) as resp:
            data = resp.read()
        with open(dest, "wb") as fh:
            fh.write(data)
        got = _md5(dest)
        if got == EXPECTED_MD5:
            return dest, f"downloaded from canonical URL {CANONICAL_URL} (md5 verified)"
        else:
            sys.stderr.write(
                f"[warn] canonical download md5 {got} != expected {EXPECTED_MD5}; "
                "falling back to local copy\n")
    except Exception as exc:  # network blocked / unreachable
        sys.stderr.write(f"[warn] canonical download failed ({exc!r}); "
                         "falling back to local copy\n")

    # 2. Local md5-gated fallback.
    if not os.path.exists(LOCAL_FALLBACK):
        raise FileNotFoundError(
            f"Canonical download unavailable and local fallback not found at "
            f"{LOCAL_FALLBACK}")
    got = _md5(LOCAL_FALLBACK)
    if got != EXPECTED_MD5:
        raise ValueError(
            f"Local fallback md5 {got} != expected {EXPECTED_MD5}; refusing to "
            "proceed with an unverified source file.")
    return LOCAL_FALLBACK, f"local md5-verified copy ({LOCAL_FALLBACK})"


# ----------------------------------------------------------------------------
# Parsing
# ----------------------------------------------------------------------------
def strip_html(s):
    if pd.isna(s):
        return None
    s = re.sub(r"<sup>.*?</sup>", "", str(s))   # drop footnote superscripts
    s = re.sub(r"<[^>]+>", "", s)                # drop remaining tags
    return re.sub(r"\s+", " ", s).strip()


def is_bold(s):
    return isinstance(s, str) and "<b>" in s


# Priority tier tokens: standalone E / D / O (not embedded in words).
TIER_RE = re.compile(r"(?<![A-Za-z0-9])([EDO])(?![A-Za-z0-9])")
TIERS = ["E", "D", "O"]
TIER_NAMES = {"E": "Essential", "D": "Desirable", "O": "Optional"}


def tier_tokens(s):
    return TIER_RE.findall(s) if s else []


def parse_table(path: str) -> tuple[pd.DataFrame, str, str]:
    """Parse Table_1.xls into a long-format criterion table.

    Returns (dataframe, use_case_1_name, use_case_2_name).
    """
    raw = pd.read_excel(path, sheet_name=0, header=None)
    uc1 = strip_html(raw.iloc[0, 1])
    uc2 = strip_html(raw.iloc[0, 2])
    rows = []
    section = None
    for i in range(1, raw.shape[0]):
        c0_raw = raw.iloc[i, 0]
        c0 = strip_html(c0_raw)
        c1 = strip_html(raw.iloc[i, 1])
        c2 = strip_html(raw.iloc[i, 2])
        is_header = (c1 is None and c2 is None)
        if is_header:
            if is_bold(c0_raw):          # bold => top-level section header
                section = c0
            # non-bold header (e.g. "Species detection/differentiation:") is a
            # sub-label with no data of its own; skip but keep the section.
            continue
        rows.append(dict(
            row_index=i, section=section, criterion=c0,
            cm_text=c1, ss_text=c2,
            cm_tiers=tier_tokens(c1), ss_tiers=tier_tokens(c2),
        ))
    return pd.DataFrame(rows), uc1, uc2


# ----------------------------------------------------------------------------
# Statistics helpers
# ----------------------------------------------------------------------------
def freeman_halton_exact(table) -> float:
    """Exact p-value for a 2x3 contingency table with fixed margins
    (Freeman-Halton extension of Fisher's exact test).

    Enumerates all tables consistent with the observed margins and sums the
    probability of every table no more likely than the observed one.
    """
    table = np.asarray(table, dtype=int)
    assert table.shape == (2, 3), "expects a 2x3 table"
    rs = table.sum(1)
    cs = table.sum(0)
    N = table.sum()

    def logprob(t):
        t = np.asarray(t)
        return (sum(gammaln(r + 1) for r in rs)
                + sum(gammaln(c + 1) for c in cs)
                - gammaln(N + 1)
                - sum(gammaln(x + 1) for x in t.ravel()))

    obs_lp = logprob(table)
    tol = 1e-7
    total = 0.0
    R0 = rs[0]
    for a in range(0, min(R0, cs[0]) + 1):
        for b in range(0, min(R0 - a, cs[1]) + 1):
            c = R0 - a - b
            if c < 0 or c > cs[2]:
                continue
            t = [[a, b, c], [cs[0] - a, cs[1] - b, cs[2] - c]]
            if any(x < 0 for row in t for x in row):
                continue
            lp = logprob(t)
            if lp <= obs_lp + tol:
                total += mexp(lp)
    return float(min(total, 1.0))


def cramers_v(table) -> float:
    table = np.asarray(table, dtype=float)
    chi2 = stats.chi2_contingency(table, correction=False)[0]
    n = table.sum()
    k = min(table.shape)
    return float(np.sqrt((chi2 / n) / (k - 1)))


def holm_correction(pvals):
    """Holm-Bonferroni corrected p-values, returned in the input order."""
    p = np.asarray(pvals, dtype=float)
    m = len(p)
    order = np.argsort(p)
    corrected = np.empty(m)
    running = 0.0
    for rank, idx in enumerate(order):
        val = (m - rank) * p[idx]
        running = max(running, val)
        corrected[idx] = min(running, 1.0)
    return corrected.tolist()


def norm_text(s):
    return re.sub(r"[^a-z0-9]+", " ", (s or "").lower()).strip()


# ----------------------------------------------------------------------------
# Test 1: priority-tier composition across use-cases
# ----------------------------------------------------------------------------
def tier_count_table(df, tiers=TIERS, presence_absence=False):
    """2x3 counts (rows = [use-case 1, use-case 2], cols = tiers).

    When presence_absence=True each criterion contributes at most one count per
    tier (collapses within-cell repeats), removing pseudoreplication.
    """
    def counts(col):
        c = collections.Counter()
        for lst in col:
            items = set(lst) if presence_absence else lst
            for t in items:
                c[t] += 1
        return c
    cm = counts(df.cm_tiers)
    ss = counts(df.ss_tiers)
    return np.array([[cm.get(t, 0) for t in tiers],
                     [ss.get(t, 0) for t in tiers]], dtype=int)


def run_test1(df):
    ct = tier_count_table(df, presence_absence=False)
    n = int(ct.sum())
    fh_p = freeman_halton_exact(ct)
    chi2, chi_p, dof, expected = stats.chi2_contingency(ct, correction=False)
    v = cramers_v(ct)
    # sensitivity: presence/absence (one count per criterion per tier)
    ct_pa = tier_count_table(df, presence_absence=True)
    fh_p_pa = freeman_halton_exact(ct_pa)
    v_pa = cramers_v(ct_pa)
    return {
        "test": "Freeman-Halton exact test (2x3), priority tier x use-case",
        "contingency_table": ct.tolist(),
        "contingency_rows": ["use_case_1", "use_case_2"],
        "contingency_cols": TIERS,
        "n_items": n,
        "statistic_chi2": float(chi2),
        "chi2_p": float(chi_p),
        "chi2_dof": int(dof),
        "min_expected_count": float(expected.min()),
        "exact_p": fh_p,
        "effect_size_cramers_v": v,
        "sensitivity_presence_absence": {
            "contingency_table": ct_pa.tolist(),
            "n_items": int(ct_pa.sum()),
            "exact_p": fh_p_pa,
            "effect_size_cramers_v": v_pa,
        },
    }


# ----------------------------------------------------------------------------
# Test 2: cross-use-case concordance per criterion
# ----------------------------------------------------------------------------
def run_test2(df):
    df = df.copy()
    df["identical_raw"] = (df.cm_text.astype(str).str.strip()
                           == df.ss_text.astype(str).str.strip())
    df["identical_norm"] = df.apply(
        lambda r: norm_text(r.cm_text) == norm_text(r.ss_text), axis=1)
    n = len(df)
    k_diff = int((~df.identical_raw).sum())
    bt = stats.binomtest(k_diff, n, 0.5)
    ci = bt.proportion_ci(method="wilson")
    # sensitivity: normalized comparison
    k_diff_norm = int((~df.identical_norm).sum())
    bt2 = stats.binomtest(k_diff_norm, n, 0.5)
    ci2 = bt2.proportion_ci(method="wilson")
    return {
        "test": "Exact binomial test, proportion of criteria differing between "
                "use-cases vs 0.5",
        "n_criteria": n,
        "n_differ": k_diff,
        "n_identical": int(df.identical_raw.sum()),
        "proportion_differ": k_diff / n,
        "wilson_95ci": [float(ci.low), float(ci.high)],
        "exact_p": float(bt.pvalue),
        "sensitivity_normalized": {
            "n_differ": k_diff_norm,
            "n_identical": int(df.identical_norm.sum()),
            "proportion_differ": k_diff_norm / n,
            "wilson_95ci": [float(ci2.low), float(ci2.high)],
            "exact_p": float(bt2.pvalue),
        },
    }, df


# ----------------------------------------------------------------------------
# Descriptive requirement map
# ----------------------------------------------------------------------------
# Ordinal stringency: Essential > Desirable > Optional > None-required.
STRINGENCY = {"E": 3, "D": 2, "O": 1, "None": 0}
STRINGENCY_LABEL = {3: "Essential", 2: "Desirable", 1: "Optional", 0: "Not required"}


def most_stringent(tiers):
    """Most stringent tier present in a cell, or 'None' when no E/D/O tag."""
    if not tiers:
        return "None"
    return max(tiers, key=lambda t: STRINGENCY[t])


def build_requirement_map(df, uc1, uc2):
    rows = []
    for _, r in df.iterrows():
        rows.append(dict(
            section=r.section, criterion=r.criterion,
            use_case_1_tier=most_stringent(r.cm_tiers),
            use_case_2_tier=most_stringent(r.ss_tiers),
        ))
    m = pd.DataFrame(rows)
    m["use_case_1_stringency"] = m.use_case_1_tier.map(STRINGENCY)
    m["use_case_2_stringency"] = m.use_case_2_tier.map(STRINGENCY)
    return m


# ----------------------------------------------------------------------------
# Figures  (colourblind-safe Okabe-Ito palette; no tool/platform names)
# ----------------------------------------------------------------------------
UC_SHORT = ("Case management\n(elimination)", "Screening /\nsurveillance")
# Okabe-Ito colourblind-safe hues.
TIER_COLOR = {"E": "#0072B2", "D": "#E69F00", "O": "#009E73", "None": "#D0D0D0"}


def _style():
    mpl.rcParams.update({
        "figure.dpi": 150, "savefig.dpi": 150,
        "font.size": 11, "axes.titlesize": 12, "axes.labelsize": 11,
        "xtick.labelsize": 10, "ytick.labelsize": 10, "legend.fontsize": 10,
        "axes.spines.top": False, "axes.spines.right": False,
        "font.family": "DejaVu Sans", "figure.autolayout": False,
    })


def fig1_tier_composition(t1, path):
    ct = np.array(t1["contingency_table"])
    _style()
    fig, ax = plt.subplots(figsize=(7.2, 4.6))
    x = np.arange(len(TIERS))
    w = 0.38
    b1 = ax.bar(x - w / 2, ct[0], w, label=UC_SHORT[0].replace("\n", " "),
                color="#0072B2")
    b2 = ax.bar(x + w / 2, ct[1], w, label=UC_SHORT[1].replace("\n", " "),
                color="#E69F00")
    for bars in (b1, b2):
        for rect in bars:
            h = rect.get_height()
            ax.annotate(f"{int(h)}", (rect.get_x() + rect.get_width() / 2, h),
                        ha="center", va="bottom", fontsize=9,
                        xytext=(0, 1), textcoords="offset points")
    ax.set_xticks(x)
    ax.set_xticklabels([f"{TIER_NAMES[t]}\n({t})" for t in TIERS])
    ax.set_ylabel("Number of tagged requirement items")
    ax.set_xlabel("Priority tier")
    ax.set_title("Priority-tier composition is similar across use-cases\n"
                 f"(Freeman-Halton exact p = {t1['exact_p']:.2f}, "
                 f"Cramer's V = {t1['effect_size_cramers_v']:.2f}, "
                 f"n = {t1['n_items']} items)", fontsize=11)
    ax.set_ylim(0, ct.max() * 1.28)
    ax.legend(frameon=False, loc="upper center", ncol=2,
              bbox_to_anchor=(0.5, 1.0))
    fig.tight_layout()
    fig.savefig(path)
    plt.close(fig)


def fig2_concordance(t2, path):
    _style()
    fig, ax = plt.subplots(figsize=(6.2, 4.6))
    cats = ["Differ", "Identical"]
    vals = [t2["n_differ"], t2["n_identical"]]
    colors = ["#0072B2", "#CC79A7"]
    bars = ax.bar(cats, vals, color=colors, width=0.6)
    for rect, v in zip(bars, vals):
        ax.annotate(f"{v}", (rect.get_x() + rect.get_width() / 2, v),
                    ha="center", va="bottom", fontsize=10,
                    xytext=(0, 1), textcoords="offset points")
    lo, hi = t2["wilson_95ci"]
    ax.set_ylabel("Number of criteria")
    ax.set_title(
        "Criteria differ between use-cases about as often as not\n"
        f"proportion differing = {t2['proportion_differ']:.2f} "
        f"(Wilson 95% CI {lo:.2f}\u2013{hi:.2f})\n"
        f"exact binomial vs 0.5: p = {t2['exact_p']:.2f}, "
        f"n = {t2['n_criteria']}", fontsize=10)
    ax.set_ylim(0, max(t2["n_differ"], t2["n_identical"]) * 1.18)
    fig.tight_layout()
    fig.savefig(path)
    plt.close(fig)


def fig3_requirement_map(m, uc1, uc2, path):
    _style()
    order = {"Technical specifications": 0,
             "Health systems and technical specifications": 1}
    m = m.copy()
    m["_o"] = m.section.map(lambda s: order.get(s, 9))
    m = m.sort_values(["_o", "row_index"] if "row_index" in m else ["_o"]).reset_index(drop=True)
    n = len(m)
    grid = np.array([m.use_case_1_stringency.values,
                     m.use_case_2_stringency.values]).T  # n x 2
    fig, ax = plt.subplots(figsize=(7.6, 0.42 * n + 1.8))
    from matplotlib.colors import ListedColormap, BoundaryNorm
    cmap = ListedColormap(["#D0D0D0", "#009E73", "#E69F00", "#0072B2"])
    norm = BoundaryNorm([-0.5, 0.5, 1.5, 2.5, 3.5], cmap.N)
    ax.imshow(grid, cmap=cmap, norm=norm, aspect="auto")
    ax.set_xticks([0, 1])
    ax.set_xticklabels([UC_SHORT[0], UC_SHORT[1]])
    ax.set_yticks(np.arange(n))
    ax.set_yticklabels(m.criterion.values, fontsize=8)
    # annotate each cell with tier letter
    letter = {3: "E", 2: "D", 1: "O", 0: "-"}
    for i in range(n):
        for j, col in enumerate(["use_case_1_stringency", "use_case_2_stringency"]):
            val = grid[i, j]
            ax.text(j, i, letter[val], ha="center", va="center",
                    color="white" if val in (3,) else "black", fontsize=8)
    # section separators
    sec_boundaries = m.groupby("_o").size().cumsum().values[:-1]
    for b in sec_boundaries:
        ax.axhline(b - 0.5, color="black", lw=1.2)
    ax.set_title("Most-stringent required tier per criterion, by use-case",
                 fontsize=12)
    handles = [Patch(facecolor="#0072B2", label="Essential (E)"),
               Patch(facecolor="#E69F00", label="Desirable (D)"),
               Patch(facecolor="#009E73", label="Optional (O)"),
               Patch(facecolor="#D0D0D0", label="Not required (-)")]
    ax.legend(handles=handles, frameon=False, ncol=4,
              loc="upper center", bbox_to_anchor=(0.5, -0.06), fontsize=9)
    ax.set_xticks(np.arange(-0.5, 2, 1), minor=True)
    ax.set_yticks(np.arange(-0.5, n, 1), minor=True)
    ax.grid(which="minor", color="white", lw=1.0)
    ax.tick_params(which="minor", length=0)
    fig.tight_layout()
    fig.savefig(path, bbox_inches="tight")
    plt.close(fig)


# ----------------------------------------------------------------------------
# Table / datapackage emission
# ----------------------------------------------------------------------------
def write_datapackage(tbl_specs, path):
    pkg = {
        "name": "malaria-diagnostic-tpp-characterization",
        "title": "Malaria diagnostic Target Product Profile characterization tables",
        "licenses": [{"name": "CC-BY-4.0",
                      "title": "Creative Commons Attribution 4.0",
                      "path": "https://creativecommons.org/licenses/by/4.0/"}],
        "sources": [{"title": "Target product profiles for malaria diagnostics "
                     "(malERA consultative group), supplementary Table 1",
                     "path": f"https://doi.org/{SOURCE_DOI}"}],
        "resources": tbl_specs,
    }
    with open(path, "w") as fh:
        json.dump(pkg, fh, indent=2)


def _fields(df):
    typemap = {"int64": "integer", "float64": "number", "bool": "boolean",
               "object": "string"}
    return [{"name": c, "type": typemap.get(str(df[c].dtype), "string")}
            for c in df.columns]


# ----------------------------------------------------------------------------
# Main
# ----------------------------------------------------------------------------
def main():
    path, provenance = acquire_data()
    src_md5 = _md5(path)
    df, uc1, uc2 = parse_table(path)

    t1 = run_test1(df)
    t2, df2 = run_test2(df)
    req_map = build_requirement_map(df, uc1, uc2)
    req_map["row_index"] = df["row_index"].values

    # Multiple-comparison correction across the 2 primary tests.
    holm = holm_correction([t1["exact_p"], t2["exact_p"]])
    t1["holm_corrected_p"] = holm[0]
    t2["holm_corrected_p"] = holm[1]

    # ---- figures ----
    f1 = os.path.join(FIG_DIR, "fig-1-tier-composition.png")
    f2 = os.path.join(FIG_DIR, "fig-2-concordance.png")
    f3 = os.path.join(FIG_DIR, "fig-3-requirement-map.png")
    fig1_tier_composition(t1, f1)
    fig2_concordance(t2, f2)
    fig3_requirement_map(req_map, uc1, uc2, f3)

    # ---- tables ----
    tbl1 = pd.DataFrame(t1["contingency_table"],
                        columns=[f"tier_{t}" for t in TIERS])
    tbl1.insert(0, "use_case", ["use_case_1", "use_case_2"])
    tbl1_path = os.path.join(TBL_DIR, "tbl-1-tier-counts.csv")
    tbl1.to_csv(tbl1_path, index=False)

    tbl2 = df2[["section", "criterion", "cm_text", "ss_text",
                "identical_raw", "identical_norm"]].copy()
    tbl2 = tbl2.rename(columns={"cm_text": "use_case_1_text",
                                "ss_text": "use_case_2_text",
                                "identical_raw": "identical",
                                "identical_norm": "identical_normalized"})
    tbl2_path = os.path.join(TBL_DIR, "tbl-2-criterion-concordance.csv")
    tbl2.to_csv(tbl2_path, index=False)

    tbl3 = req_map[["section", "criterion", "use_case_1_tier", "use_case_2_tier",
                    "use_case_1_stringency", "use_case_2_stringency"]].copy()
    tbl3_path = os.path.join(TBL_DIR, "tbl-3-requirement-map.csv")
    tbl3.to_csv(tbl3_path, index=False)

    write_datapackage([
        {"name": "tbl-1-tier-counts", "path": "tbl-1-tier-counts.csv",
         "profile": "tabular-data-resource", "format": "csv",
         "title": "Priority-tier counts (E/D/O) by use-case",
         "schema": {"fields": _fields(tbl1)}},
        {"name": "tbl-2-criterion-concordance",
         "path": "tbl-2-criterion-concordance.csv",
         "profile": "tabular-data-resource", "format": "csv",
         "title": "Per-criterion cross-use-case concordance",
         "schema": {"fields": _fields(tbl2)}},
        {"name": "tbl-3-requirement-map", "path": "tbl-3-requirement-map.csv",
         "profile": "tabular-data-resource", "format": "csv",
         "title": "Most-stringent required tier per criterion by use-case",
         "schema": {"fields": _fields(tbl3)}},
    ], os.path.join(TBL_DIR, "datapackage.json"))

    # ---- descriptive summary counts for results.json ----
    map_summary = {
        "use_case_1": req_map.use_case_1_tier.value_counts().to_dict(),
        "use_case_2": req_map.use_case_2_tier.value_counts().to_dict(),
    }
    n_more_stringent_cm = int((req_map.use_case_1_stringency
                               > req_map.use_case_2_stringency).sum())
    n_more_stringent_ss = int((req_map.use_case_2_stringency
                               > req_map.use_case_1_stringency).sum())
    n_equal = int((req_map.use_case_1_stringency
                   == req_map.use_case_2_stringency).sum())

    results = {
        "meta": {
            "seed": SEED,
            "source_doi": SOURCE_DOI,
            "source_licence": SOURCE_LICENCE,
            "source_provenance": provenance,
            "source_md5": src_md5,
            "expected_md5": EXPECTED_MD5,
            "generated_utc": datetime.datetime.utcnow().isoformat() + "Z",
            "use_case_1": uc1,
            "use_case_2": uc2,
            "n_substantive_criteria": int(len(df)),
            "n_sections": int(df.section.nunique()),
            "multiple_comparison_correction": "Holm-Bonferroni across 2 primary tests",
        },
        "test1_tier_composition": t1,
        "test2_concordance": t2,
        "descriptive_requirement_map": {
            "tier_distribution": map_summary,
            "n_criteria_more_stringent_use_case_1": n_more_stringent_cm,
            "n_criteria_more_stringent_use_case_2": n_more_stringent_ss,
            "n_criteria_equal_stringency": n_equal,
        },
    }
    with open(os.path.join(HERE, "results.json"), "w") as fh:
        json.dump(results, fh, indent=2)

    print(json.dumps(results, indent=2))
    return results


if __name__ == "__main__":
    main()