#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
# Copyright (C) 2026 World in Transition
# This program is free software: you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version. It is distributed WITHOUT ANY WARRANTY. Full text:
# https://worldintransition.net/downloads/models/LICENSE-GPL-3.0.txt
"""
Is Money Necessary? A model of coordination with and without currency.
=====================================================================

Companion script for
  https://worldintransition.net/explore/research/money-and-coordination/

What it does
------------
1. FEASIBILITY DEMO. Poses the core coordination problem -- allocate finite
   production capacity to real needs -- as a linear program and solves it with
   no prices, no money and no purchasing power anywhere in the input. It then
   prints the dual variables ("shadow values"), which are the internal
   valuations Kantorovich (1939) called "objectively determined valuations":
   numbers the solver computes, not money anyone holds.

2. COMPARISON MODEL. Compares the same population, needs and production
   capacity under two allocation regimes:

     Currency-based:  W_cur  = (1-L1)(1-L2)(1-L3)(1-L4)(1-L5)(1-L6)(1-L7)
     Currency-free:   W_free = (1-C_coord)(1-C_fed)(1-C_local)(1-C_agg)(1-C_dyn)(1-C_inc)(1-C_err)

   W is welfare-equivalent output as a share of what the capacity K could
   deliver if allocated to needs with no losses at all (W = 1).
   L1..L7 are the loss channels of allocation by purchasing power; the C terms
   are the costs of coordinating without currency (including the Hayek
   knowledge problem and incentive risk). All are fractions of K.

   Currency-free allocation is more efficient iff
        W_free > W_cur
   <=>  sum_k -ln(1-L_k)  >  sum_j -ln(1-C_j)
   (approximately: L1+...+L7 > C_coord+C_fed+C_local+C_agg+C_dyn+C_inc+C_err)

   The currency-free alternative is modelled as POLYCENTRIC coordination, not
   central planning: local nodes match local needs to local capacity using
   feedback signals (stock levels, waiting lists), and exchange surpluses and
   node-level aggregates with neighbouring nodes / federations. Its costs are
   split into local terms (C_coord, C_local) and cross-node terms (C_fed, C_agg).

   Two comparisons are reported:
     * OUTPUT-ONLY: L1, L2, L4, L5, L6, L7  vs  C_coord, C_fed, C_local, C_agg, C_dyn, C_inc
       (counts only goods and labour; ignores who receives them)
     * WELLBEING:   adds L3 (allocation by income rather than need) and C_err
       (need-assessment error in the currency-free system). Both are measured
       with the same concave wellbeing function, so the comparison is symmetric.

3. MONTE CARLO + SENSITIVITY. Samples every parameter over a documented
   plausible range (fixed seed), reports the share of parameter space in which
   currency-free allocation comes out ahead, a one-at-a-time tornado analysis,
   rank correlations, and break-even values.

4. TWO REGIMES (printed as section 5). Section 1a poses the hard case,
   capacity physically below need. This part runs the same construction from
   shortage (0.8x need) to sufficiency (1.0x-3.0x need) and compares, for the
   SAME needs and capacity, the needs-based allocation with a budget-
   constrained market at three levels of income inequality (low, US-like,
   high; equal incomes as a reference). It also re-runs the polycentric
   network of 1b with spare capacity. It uses its own random generator, so
   it changes none of the numbers in the other sections.

5. CHARTS. Writes four SVG charts: to public/images/models/ when run inside
   the website's repository, otherwise to a new folder money-model-charts/ in
   the directory you run it from.

Every parameter below carries its source or is explicitly marked ASSUMPTION.
Change them and re-run: the conclusions are meant to be checked, not trusted.

Requirements: Python 3.9+, numpy, matplotlib; scipy (optional, for the
feasibility and two-regimes demos -- skipped with a message if absent).

Run:  python3 money-vs-moneyless-model.py            (from anywhere)
      python3 money-vs-moneyless-model.py --no-charts
"""

import math
import os
import sys
import time

import numpy as np

SEED = 20261011           # fixed seed: results are exactly reproducible
N_DRAWS = 200_000         # Monte Carlo draws

# ---------------------------------------------------------------------------
# Helper: closed-form Atkinson index for a log-normal distribution
# ---------------------------------------------------------------------------
# For identical CRRA wellbeing u(c) = c^(1-e)/(1-e)  (u = ln c when e = 1),
# Jensen's inequality implies: for a fixed total of goods, an equal (need-
# matched) allocation gives strictly higher total wellbeing than any unequal
# allocation. The size of the gap is the Atkinson (1970) index A_e: the share
# of total goods that could be given up under the need-matched allocation
# while still reaching the same total wellbeing. For a log-normal with log-sd
# sigma, A_e = 1 - exp(-e * sigma^2 / 2)  (exact, all e > 0).
# Gini of a log-normal: G = 2*Phi(sigma/sqrt(2)) - 1.


def _norm_ppf(p):
    """Inverse standard normal CDF (Acklam's rational approximation, |err|<1e-9)."""
    a = [-3.969683028665376e01, 2.209460984245205e02, -2.759285104469687e02,
         1.383577518672690e02, -3.066479806614716e01, 2.506628277459239e00]
    b = [-5.447609879822406e01, 1.615858368580409e02, -1.556989798598866e02,
         6.680131188771972e01, -1.328068155288572e01]
    c = [-7.784894002430293e-03, -3.223964580411365e-01, -2.400758277161838e00,
         -2.549732539343734e00, 4.374664141464968e00, 2.938163982698783e00]
    d = [7.784695709041462e-03, 3.224671290700398e-01, 2.445134137142996e00,
         3.754408661907416e00]
    p = np.asarray(p, dtype=float)
    out = np.empty_like(p)
    lo, hi = 0.02425, 1 - 0.02425
    m = p < lo
    q = np.sqrt(-2 * np.log(p[m]))
    out[m] = (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) / ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1)
    m = (p >= lo) & (p <= hi)
    q = p[m] - 0.5
    r = q*q
    out[m] = (((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r+a[5])*q / (((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r+1)
    m = p > hi
    q = np.sqrt(-2 * np.log(1 - p[m]))
    out[m] = -(((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) / ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1)
    return out


def sigma_from_gini(g):
    return math.sqrt(2) * _norm_ppf((np.asarray(g) + 1) / 2)


def atkinson_lognormal(sigma, eps):
    return 1 - np.exp(-np.asarray(eps) * np.asarray(sigma) ** 2 / 2)


def atkinson_sample(x, eps):
    """Atkinson index of an empirical sample x (numerical check)."""
    x = np.asarray(x, dtype=float)
    mu = x.mean()
    if abs(eps - 1) < 1e-12:
        ede = np.exp(np.mean(np.log(x)))
    else:
        ede = np.mean(x ** (1 - eps)) ** (1 / (1 - eps))
    return 1 - ede / mu


def gini_sample(x):
    x = np.sort(np.asarray(x, dtype=float))
    n = len(x)
    i = np.arange(1, n + 1)
    return (2 * np.sum(i * x) / (n * np.sum(x))) - (n + 1) / n


# ---------------------------------------------------------------------------
# PARAMETERS
# Each entry: key -> (low, mode, high, label, source)
# Sampled as triangular(low, mode, high); uniform where low..high and mode=None.
# All values are fractions (0.039 = 3.9 %).
# ---------------------------------------------------------------------------

P = {
    # ---- L1 Overhead: labour devoted purely to handling money -------------
    "fin":   (0.035, 0.039, 0.045, "Finance & insurance employment share",
              "BLS 2024: 6.70 M of 170.0 M jobs = 3.9 % (Rieley & Colato, MLR Jan 2026)"),
    "hlth":  (0.015, 0.023, 0.033, "Excess health-care administration",
              "Himmelstein et al. 2020: US admin 34.2 % of health spending vs 17.0 % Canada; "
              "excess ~3.3 % of GDP (2017). Lower values remove overlap with insurer staff "
              "already counted in finance & insurance (ASSUMPTION on overlap)"),
    "tax":   (0.018, 0.022, 0.026, "Tax-compliance labour",
              "Tax Foundation 2024: 7.9 bn hours = 3.8 M full-time workers; / 170 M jobs = 2.2 %"),
    "adv":   (0.013, 0.015, 0.020, "Advertising",
              "St. Louis Fed (Marto & Le 2024): US advertising 1.3-2 % of GDP over 20 years"),
    "other": (0.005, 0.015, 0.030, "Other money handling (billing, collections, payments, "
              "pricing, transactional legal/real-estate work)",
              "ASSUMPTION: broad range; the site's projections page estimates more"),
    "rho":   (0.25, None, 0.50, "Residual share still needed without money",
              "ASSUMPTION: record-keeping, risk pooling, product information and logistics "
              "accounting persist in any system; 25-50 % of gross overhead kept"),

    # ---- L2 Access-gating: capacity idle because need lacks purchasing power
    "food_gate": (0.10, 0.25, 0.40, "Share of surplus food lost because of price/purchasing-power gating",
                  "ASSUMPTION applied to ReFED 2023: surplus food 31 % of supply, $382 bn = 1.4 % of GDP"),
    "house_use": (0.10, 0.25, 0.50, "Usable share of 'other vacant' homes",
                  "ASSUMPTION applied to Census 2020: 3.36 M 'other vacant' of 140.5 M units; "
                  "HUD 2024: 771,480 people homeless on one night"),
    "idle":  (0.005, 0.010, 0.020, "Idle labour and plant (output below potential)",
              "FRED GDPC1/GDPPOT 1980-2026: mean shortfall below CBO potential = 1.0-1.1 % "
              "(computed for this page)"),

    # ---- L3 Allocation by income rather than need ---------------------------
    "gini":  (0.30, 0.35, 0.40, "Gini of household resources (US-like)",
              "OECD: US disposable-income Gini 0.39 (2022); consumption inequality is lower "
              "(Krueger & Perri 2006; Attanasio & Pistaferri 2016) -> range 0.30-0.40"),
    "eps":   (0.70, 1.30, 2.00, "Curvature of wellbeing (elasticity of marginal utility)",
              "Layard, Mayraz & Nickell 2008: 1.26 from wellbeing data; "
              "Acland & Greenberg 2023 meta-analysis: 1.61 (95 % CI 1.18-2.05)"),

    # ---- L4 Market power ---------------------------------------------------
    "L4":    (0.015, 0.035, 0.080, "Market-power (markup) losses",
              "Edmond, Midrigan & Xu 2023 (JPE): total cost of markups 7.5 % of consumption, "
              "about 1/3 from misallocation; Baqaee & Farhi 2020: up to ~15-20 % TFP (heavy imputation); "
              "De Loecker et al. 2020: avg markup 1.21 (1980) -> 1.61 (2016)"),

    # ---- L5 Competitive duplication & marketing-driven demand ---------------
    "mkt_ind": (0.005, 0.02, 0.05, "Production for advertising-induced demand that serves no need",
                "ASSUMPTION informed by Molinari & Turino 2018 (Economic Journal): advertising raises "
                "aggregate consumption AND hours worked, leaving households worse off on balance; "
                "Benhabib & Bisin 2001. No verified point estimate is used - wide range"),
    "dup":     (0.000, 0.01, 0.03, "Competitive duplication (parallel lines/infrastructure) net of "
                "the value of variety and redundancy",
                "ASSUMPTION: Mankiw & Whinston 1986 show free entry can be excessive; but Edmond, "
                "Midrigan & Xu 2023 find inefficient-entry costs small -> low central value"),

    # ---- L6 Premature obsolescence / shortened lifespans --------------------
    "life_ext": (1.05, 1.20, 1.50, "Achievable average lifetime extension factor for durables & clothing",
                 "UBA/Oeko-Institut 2016: washing machines replaced within 5 years due to a defect "
                 "rose from 3.5 % (2004) to 8.3 % (2012/13); Wieser & Troger 2018: most phones replaced "
                 "while still (partly) working; EEB 2019 Coolproducts. Range is an ASSUMPTION"),
    "obs_resid": (0.20, 0.35, 0.50, "Share of the saving used up by repair labour and sturdier design",
                  "ASSUMPTION"),

    # ---- L7 Market coordination errors --------------------------------------
    "L7":    (0.003, 0.008, 0.020, "Unsold / destroyed / returned non-food goods, failed-venture losses",
              "ASSUMPTION informed by: US returns $685 bn (13.2 % of retail sales, Appriss/Deloitte 2024), "
              "~9.5 bn lb of returns landfilled (Optoro); France: EUR 630 M of unsold non-food goods "
              "destroyed per year before the 2020 AGEC ban. Excludes food (in L2) and the output gap (in L2)"),

    # ---- Costs of currency-free coordination ---------------------------------
    # The alternative modelled here is NOT central planning. It is polycentric
    # coordination: many local nodes match local needs to local capacity using
    # feedback signals (stock levels, waiting lists, usage data), and exchange
    # information and surpluses with neighbouring nodes and federations.
    # Costs are therefore split into a LOCAL part and a CROSS-NODE part.
    "C_coord": (0.02, 0.04, 0.08, "Local coordination labour (matching, deliberation, records)",
                "ASSUMPTION: computation itself is cheap (see demo); the cost is people's time "
                "for gathering needs data, deliberating and keeping open records"),
    "C_fed":   (0.005, 0.015, 0.04, "Federation labour (cross-node information exchange, clearing)",
                "ASSUMPTION: regional/interregional coordination of node-level aggregates only"),
    "C_local": (0.01, 0.03, 0.08, "Local mismatch (wrong goods / wrong time within a node)",
                "ASSUMPTION: GROSS (the market's own mismatch is now charged separately as L7); "
                "lower than under central planning because decisions use local and tacit knowledge"),
    "C_agg":   (0.005, 0.02, 0.06, "Cross-node aggregation error (surpluses and shortfalls between nodes)",
                "ASSUMPTION: the part of Hayek's knowledge problem that survives decentralisation; "
                "upper end reflects documented shortage/surplus problems of planned economies (Kornai 1980)"),
    "C_dyn":   (0.00, 0.02, 0.08, "Slower discovery and innovation (level-equivalent)",
                "ASSUMPTION: Hayek/Schumpeter critique; partly offset by publicly funded and "
                "open-source innovation"),
    "C_inc":   (0.00, 0.04, 0.15, "Participation / effort shortfall",
                "Vivalt et al. 2024 (NBER w32719): $1,000/month for 3 years -> 1.3-1.4 h/week "
                "less work (5-6 % of earnings); Jones & Marinescu 2022: no employment drop from "
                "Alaska dividend. Upper end allows for much larger effects"),
    "s_err":   (0.08, None, 0.35, "Need-assessment error (log-sd)",
                "ASSUMPTION: errors in recognising needs re-create dispersion; s = 0.35 alone "
                "corresponds to a Gini of ~0.20"),
}

# Fixed structural constants (sources as above)
FOOD_SURPLUS_GDP = 0.014           # ReFED 2023: $382 bn = 1.4 % of GDP
OTHER_VACANT_SHARE = 3.363 / 140.5  # Census 2020
DURABLE_CLOTHING_SHARE = 0.093    # FRED/BEA: durable-goods PCE $2.44 tn + clothing & footwear $0.56 tn ~ 9.3 % of GDP
HOUSING_SHARE = 0.15               # ASSUMPTION: housing services ~15 % of consumption
LOW_INEQ_GINI = 0.26               # lowest-inequality OECD monetary economies (Nordic / Central Europe)
HIGH_INEQ_GINI = (0.50, 0.55)      # World Bank SI.POV.GINI 2024: Brazil 0.503, Colombia 0.544
                                   # (South Africa ~0.60-0.65 in 2005-2014 surveys; 0.541 in 2022)


def central(key):
    lo, mode, hi = P[key][:3]
    return (lo + hi) / 2 if mode is None else mode


def sample(key, rng, n):
    lo, mode, hi = P[key][:3]
    if mode is None:
        return rng.uniform(lo, hi, n)
    return rng.triangular(lo, mode, hi, n)


# ---------------------------------------------------------------------------
# Channel formulas
# ---------------------------------------------------------------------------

def channels(v):
    """v: dict of parameter arrays (or scalars). Returns dict of channel values."""
    gross = v["fin"] + v["hlth"] + v["tax"] + v["adv"] + v["other"]
    L1 = gross * (1 - v["rho"])
    L2 = (FOOD_SURPLUS_GDP * v["food_gate"]
          + OTHER_VACANT_SHARE * HOUSING_SHARE * v["house_use"]
          + v["idle"])
    L3 = atkinson_lognormal(sigma_from_gini(v["gini"]), v["eps"])
    L4 = v["L4"]
    L5 = v["mkt_ind"] + v["dup"]
    L6 = DURABLE_CLOTHING_SHARE * (1 - 1 / v["life_ext"]) * (1 - v["obs_resid"])
    L7 = v["L7"]
    C_err = atkinson_lognormal(v["s_err"], v["eps"])
    return dict(gross_overhead=gross, L1=L1, L2=L2, L3=L3, L4=L4, L5=L5, L6=L6, L7=L7,
                C_coord=v["C_coord"], C_fed=v["C_fed"], C_local=v["C_local"], C_agg=v["C_agg"],
                C_dyn=v["C_dyn"], C_inc=v["C_inc"], C_err=C_err)


COST_KEYS = ["C_coord", "C_fed", "C_local", "C_agg", "C_dyn", "C_inc"]
OUT_LOSS_KEYS = ["L1", "L2", "L4", "L5", "L6", "L7"]


def welfare(ch):
    """Return (W_cur, W_free) for output-only and wellbeing comparisons."""
    W_cur_out = 1.0
    for k in OUT_LOSS_KEYS:
        W_cur_out = W_cur_out * (1 - ch[k])
    W_free_out = 1.0
    for k in COST_KEYS:
        W_free_out = W_free_out * (1 - ch[k])
    W_cur = W_cur_out * (1 - ch["L3"])
    W_free = W_free_out * (1 - ch["C_err"])
    return W_cur_out, W_free_out, W_cur, W_free


def net_advantage(ch):
    """Percent by which currency-free welfare exceeds currency-based (output-only, wellbeing)."""
    a, b, c, d = welfare(ch)
    return (b / a - 1), (d / c - 1)


def make_needs_problem(rng, n_people, n_goods, n_res):
    """Random needs n[i,g] >= 0 and resource requirements a[r,g] (used by 1a and 5)."""
    needs = rng.gamma(2.0, 1.0, (n_people, n_goods)) * (rng.random((n_people, n_goods)) < 0.6)
    # each good draws on 1-2 of the resources (labour types, materials, energy...)
    a = rng.uniform(0.1, 1.0, (n_res, n_goods)) * (rng.random((n_res, n_goods)) < 0.35)
    a[rng.integers(0, n_res, n_goods), np.arange(n_goods)] = rng.uniform(0.1, 1.0, n_goods)
    return needs, a


def solve_needs_lp(needs, a, K):
    """Needs-based allocation as a linear program (no prices, no incomes).

    Variables: x[i,g] (only where need > 0), then t (minimum coverage share).
    Objective (lexicographic via weights): first raise the minimum coverage t that
    every person gets for every need, then cover as much total need as possible.
    """
    from scipy.optimize import linprog
    from scipy.sparse import hstack, vstack, identity, csr_matrix
    n_res = a.shape[0]
    I, G = np.nonzero(needs)
    nv = len(I)
    nij = needs[I, G]
    A_res = csr_matrix(a[:, G])                       # sum_i a[r,g] x[i,g] <= K[r]
    A_fair = hstack([-identity(nv, format="csr"),      # t*n[i,g] - x[i,g] <= 0
                     csr_matrix(nij.reshape(-1, 1))])
    A_ub = vstack([hstack([A_res, csr_matrix((n_res, 1))]), A_fair]).tocsr()
    b_ub = np.concatenate([K, np.zeros(nv)])
    bounds = [(0, u) for u in nij] + [(0, 1)]
    c = np.concatenate([-np.ones(nv) / nij.sum(), [-10.0]])
    t0 = time.perf_counter()
    r = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method="highs-ipm")
    dt = time.perf_counter() - t0
    x = r.x[:-1]
    used = A_res @ x                                   # capacity used, per resource
    return dict(t=r.x[-1], covered=x.sum() / nij.sum(), seconds=dt, variables=nv,
                used=used, x=x, I=I, G=G)


# ---------------------------------------------------------------------------
# 1. FEASIBILITY DEMO: allocation as a linear program, no prices anywhere
# ---------------------------------------------------------------------------

def feasibility_demo(rng):
    print("=" * 76)
    print("1a. FEASIBILITY: allocating capacity to needs as a linear program (no prices)")
    print("    (computational scale of one node or one federation tier)")
    print("=" * 76)
    try:
        from scipy.optimize import linprog  # noqa: F401  (used in solve_needs_lp)
    except ImportError:
        print("scipy not installed - feasibility demo skipped (pip install scipy).")
        return None

    results = []
    for n_people, n_goods, n_res in [(1_000, 50, 4), (10_000, 60, 6)]:
        needs, a = make_needs_problem(rng, n_people, n_goods, n_res)
        total_need = needs.sum(axis=0)
        # capacities: each resource covers 70-100 % of what meeting every need would require
        K = rng.uniform(0.70, 1.00, n_res) * (a @ total_need)
        r_lp = solve_needs_lp(needs, a, K)
        t_star, covered, dt, nv = r_lp["t"], r_lp["covered"], r_lp["seconds"], r_lp["variables"]
        print(f"  {n_people:>6,} people x {n_goods} goods x {n_res} resources "
              f"({nv:,} variables): solved in {dt:.2f} s")
        print(f"     every need met to at least {t_star:.1%}; total need covered {covered:.1%}")
        results.append(dict(people=n_people, goods=n_goods, resources=n_res, variables=nv,
                            seconds=dt, min_cov=t_star, covered=covered))
    print("  -> The allocation uses only needs, capacities and technical coefficients.")
    print("     (Its dual variables are internal valuations - Kantorovich's 'objectively")
    print("     determined valuations' - computed by the solver, not money anyone holds.)")
    return results


def polycentric_demo(rng, n_nodes=12, n_goods=20, periods=300, burn_in=100,
                     loss=0.10, gain=0.5, cap=0.95, verbose=True):
    """Decentralised, feedback-controlled provisioning on a ring of local nodes.

    Each node knows only its own capacity, its own stocks and its own waiting
    list, plus the waiting lists its two neighbours post. Every period it:
      1. produces with its current labour shares,
      2. meets local needs from stock,
      3. offers leftover stock to neighbours that report unmet needs
         (a share `loss` is lost in transfer),
      4. shifts labour toward goods with waiting lists and away from goods
         piling up in stock (proportional feedback, gain `gain`).
    No prices, no money, no central allocator.

    Benchmark: a full-information oracle that, every period, knows every
    node's actual needs and solves one global LP (any node may ship to any
    other with the same loss). The ratio of needs met is a rough, toy-scale
    indicator of how much is lost to keeping decisions local.

    `cap` is each node's labour capacity relative to what meeting its mean
    needs locally would take (0.95 in section 1b; section 5 re-runs the same
    network with spare capacity). `verbose=False` suppresses the printout.
    """
    if verbose:
        print()
        print("=" * 76)
        print("1b. POLYCENTRIC FEEDBACK COORDINATION vs FULL-INFORMATION OPTIMUM")
        print("=" * 76)
    try:
        from scipy.optimize import linprog
    except ImportError:
        print("scipy not installed - polycentric demo skipped.")
        return None
    M, G = n_nodes, n_goods
    # Labour per unit of good g at node m (nodes differ: local strengths)
    a = rng.uniform(0.5, 1.5, (M, G))
    lam = rng.uniform(2, 10, (M, G)) * rng.uniform(0.6, 1.4, (M, 1))  # mean needs
    # Capacity: 95 % of the labour needed to meet mean needs locally (cap = 0.95)
    Lab = cap * (a * lam).sum(axis=1)
    w = (a * lam) / (a * lam).sum(axis=1, keepdims=True)   # initial labour shares
    stock = np.zeros((M, G))
    met_dec = met_orc = need_tot = 0.0
    # Oracle LP (built once): variables y[m,g] production, z[m,n,g] shipments m->n,
    # r[n,g] needs met at n (bounded by actual needs each period). Maximise sum r.
    from scipy.sparse import lil_matrix
    ny, nz = M * G, M * M * G
    nvar = ny + nz + ny
    A_orc = lil_matrix((M + 2 * M * G, nvar))
    b_orc = np.zeros(M + 2 * M * G)
    row = 0
    for m in range(M):                                   # labour: sum_g a y <= Lab
        A_orc[row, m*G:(m+1)*G] = a[m]; b_orc[row] = Lab[m]; row += 1
    for m in range(M):                                   # shipped out <= produced
        for g in range(G):
            A_orc[row, ny + (m*M + np.arange(M))*G + g] = 1; A_orc[row, m*G + g] = -1; row += 1
    for n in range(M):                                   # met <= received (with transfer loss)
        for g in range(G):
            eff = np.full(M, 1 - loss); eff[n] = 1.0
            A_orc[row, ny + (np.arange(M)*M + n)*G + g] = -eff; A_orc[row, ny + nz + n*G + g] = 1; row += 1
    A_orc = A_orc.tocsr()
    c_orc = np.concatenate([np.zeros(ny + nz), -np.ones(ny)])
    bounds_fixed = [(0, None)] * (ny + nz)
    nbrs = [((m - 1) % M, (m + 1) % M) for m in range(M)]
    for t in range(periods):
        D = rng.poisson(lam).astype(float)
        prod = w * Lab[:, None] / a
        avail = stock + prod
        served = np.minimum(avail, D)
        left = avail - served
        unmet = D - served
        # neighbours share leftover stock with posted waiting lists
        for m in range(M):
            for n in nbrs[m]:
                give = np.minimum(left[m], unmet[n] / (1 - loss))
                left[m] -= give
                got = give * (1 - loss)
                served[n] += got
                unmet[n] -= got
        stock = left * 0.9                         # 10 % of stock spoils per period
        # feedback: move labour toward waiting lists, away from excess stock
        signal = (unmet - 0.5 * stock) / np.maximum(lam, 1e-9)
        w = w * np.exp(gain * signal * 0.2)
        w = w / w.sum(axis=1, keepdims=True)
        if t >= burn_in:
            met_dec += served.sum()
            need_tot += D.sum()
            r = linprog(c_orc, A_ub=A_orc, b_ub=b_orc,
                        bounds=bounds_fixed + [(0, d) for d in D.ravel()], method="highs")
            met_orc += -r.fun
    ratio = met_dec / met_orc
    if not verbose:
        return dict(nodes=M, goods=G, dec=met_dec/need_tot, orc=met_orc/need_tot, ratio=ratio)
    print(f"  {M} nodes x {G} goods, {periods - burn_in} measured periods, random needs each period")
    print(f"  Needs met, decentralised feedback (local info + neighbours): {met_dec/need_tot:.1%}")
    print(f"  Needs met, full-information global optimum (oracle):          {met_orc/need_tot:.1%}")
    print(f"  Decentralised result as share of the oracle:                  {ratio:.1%}")
    print("  -> Toy model, not a calibration: it shows that simple local feedback with")
    print("     neighbour sharing gets close to the full-information optimum, and that the")
    print("     remaining gap is a real cost (represented by C_local and C_agg).")
    return dict(nodes=M, goods=G, dec=met_dec/need_tot, orc=met_orc/need_tot, ratio=ratio)


# ---------------------------------------------------------------------------
# 5. TWO REGIMES: SHORTAGE AND SUFFICIENCY
# ---------------------------------------------------------------------------
# Section 1a deliberately poses the hard case: capacity physically below total
# need. Here the same construction is run across capacity levels from shortage
# to sufficiency, and the needs-based allocation (no prices) is compared with an
# allocation gated by purchasing power, for the SAME needs and the SAME capacity.
# It uses its own random generator (seed SEED + 5), so adding it changes none of
# the numbers in sections 1-4.

REGIME_LEVELS = (0.8, 1.0, 1.2, 1.5, 2.0, 3.0)
# 0.8: middle of the 0.70-1.00 shortage range used in section 1a.
# 1.0: capacity exactly equal to what meeting every need requires.
# 1.2, 1.5: modest and ample slack. 2.0, 3.0: large slack, to show how far
# purchasing-power gating persists. ASSUMPTION: illustrative levels, not
# calibrated to any real economy's ratio of capacity to need.
REGIME_GINIS = (("Equal incomes (reference)", 0.0),
                ("Low inequality (Gini 0.26)", LOW_INEQ_GINI),
                ("US-like (Gini 0.35)", 0.35),
                ("High inequality (Gini 0.525)", sum(HIGH_INEQ_GINI) / 2))


def market_allocation(needs_res, K, income, max_iter=20_000, tol=1e-10):
    """Allocation gated by purchasing power: a budget-constrained market.

    needs_res[i,r]: capacity of resource r embodied in person i's FULL need bundle.
    K[r]: capacity of resource r.  income[i]: income shares (mean 1).

    Mechanism (all assumptions chosen to be FAVOURABLE to the market):
      * Every unit of capacity is paid for and all income is spent: total income
        equals the value of all capacity (full employment, no unsold stock,
        no market power, no money-handling overhead).
      * Prices: one market-clearing price per resource (lam[r]); a good's price
        is the value of the capacity it embodies.
      * Each person spends on needs FIRST, buying the same fraction
        f_i = min(1, income_i / cost of i's full need bundle) of every need.
      * Income left after all needs are met is spent on uses beyond need: a
        composite that spends a fixed share of that income on each resource,
        equal to the resource's share of total capacity (Cobb-Douglas).
      * Incomes are independent of needs.
    Prices are found by multiplicative tatonnement until demand = capacity for
    every resource (to `tol`).
    """
    n, R = needs_res.shape
    wshare = K / K.sum()
    lam = np.ones(R)
    lam *= n / (lam @ K)
    for it in range(max_iter):
        m = income * (lam @ K) / n
        cost = needs_res @ lam
        f = np.minimum(1.0, m / cost)
        beyond = np.maximum(0.0, m - cost)
        beyond_res = beyond.sum() * wshare / lam
        D = f @ needs_res + beyond_res
        ex = D / K
        if np.max(np.abs(ex - 1)) < tol:
            break
        lam = lam * ex ** 0.5
        lam *= n / (lam @ K)
    gap_res = (1 - f) @ needs_res
    return dict(f=f, beyond_share=float(np.mean(beyond_res / K)),
                gap_share=float(np.mean(gap_res / K)), residual=float(np.max(np.abs(ex - 1))),
                iters=it)


def regimes_demo(n_people=2_000, n_goods=50, n_res=4):
    print()
    print("=" * 76)
    print("5. TWO REGIMES: SHORTAGE AND SUFFICIENCY (same needs, same capacity)")
    print("=" * 76)
    try:
        from scipy.optimize import linprog  # noqa: F401
    except ImportError:
        print("scipy not installed - regimes demo skipped.")
        return None
    rng = np.random.default_rng(SEED + 5)
    needs, a = make_needs_problem(rng, n_people, n_goods, n_res)
    need_tot = needs.sum()
    person_need = needs.sum(axis=1)
    R_req = a @ needs.sum(axis=0)          # capacity needed to meet every need, per resource
    needs_res = needs @ a.T                # (people x resources)
    z = rng.standard_normal(n_people)      # ONE draw, rescaled for each Gini (common random numbers)
    print(f"  {n_people:,} people x {n_goods} goods x {n_res} resources; capacity level s = every")
    print("  resource at s times what meeting every need requires. Needs as in section 1a.")
    print()
    print("  (a) NEEDS-BASED allocation (linear program, no prices, no incomes)")
    print(f"      {'level':>5}  {'lowest share met':>16}  {'total need met':>14}  {'capacity unused':>15}")
    nb = []
    for s in REGIME_LEVELS:
        r = solve_needs_lp(needs, a, s * R_req)
        idle = float(np.mean(1 - r["used"] / (s * R_req)))
        nb.append(dict(level=s, min=float(r["t"]), total=float(r["covered"]), idle=idle))
        # with every capacity scaled by the same factor the optimum is exactly min(s, 1)
        assert abs(r["t"] - min(s, 1.0)) < 1e-6 and abs(r["covered"] - min(s, 1.0)) < 1e-6
        print(f"      {s:>4.1f}x  {r['t']:>16.1%}  {r['covered']:>14.1%}  {idle:>15.1%}")
    print("      -> Below 1.0x every need is met to the same share; from 1.0x every need is")
    print("         met in full and the rest of the capacity is spare (slack).")
    print()
    print("  (b) PURCHASING-POWER-GATED allocation (budget-constrained market; see")
    print("      market_allocation() for the mechanism - every assumption favours the market)")
    mk = {}
    for lab, g in REGIME_GINIS:
        if g == 0:
            inc = np.ones(n_people)
        else:
            inc = np.exp(float(sigma_from_gini(g)) * z)
            inc = inc / inc.mean()
        g_s = gini_sample(inc) if g > 0 else 0.0
        print(f"      {lab} - sample Gini {g_s:.3f}")
        print(f"      {'level':>5}  {'people short':>12}  {'lowest met':>10}  {'poorest 10%':>11}"
              f"  {'total met':>9}  {'to beyond-need':>14}  {'would close gaps':>16}")
        rows = []
        for s in REGIME_LEVELS:
            r = market_allocation(needs_res, s * R_req, inc)
            f = r["f"]
            short = float(np.mean(f < 1 - 1e-9))
            low = float(f.min())
            poor = float(np.mean(f[inc <= np.quantile(inc, 0.10)])) if g > 0 else float("nan")
            total = float((f * person_need).sum() / need_tot)
            rows.append(dict(level=s, short=short, min=low, poor10=poor, total=total,
                             beyond=r["beyond_share"], gap=r["gap_share"], residual=r["residual"]))
            poor_s = f"{poor:>11.1%}" if g > 0 else f"{'-':>11}"
            print(f"      {s:>4.1f}x  {short:>12.1%}  {low:>10.1%}  {poor_s}  {total:>9.1%}"
                  f"  {r['beyond_share']:>14.1%}  {r['gap_share']:>16.1%}")
        # finer grid of capacity levels for the chart (Figure 4)
        grid = np.round(np.arange(0.80, 3.0001, 0.05), 2)
        curve = []
        for s in grid:
            r = market_allocation(needs_res, s * R_req, inc)
            f = r["f"]
            rows_ok = r["residual"] < 1e-6
            assert rows_ok, "market prices did not converge"
            poor = float(np.mean(f[inc <= np.quantile(inc, 0.10)])) if g > 0 else float("nan")
            curve.append((float(s), float((f * person_need).sum() / need_tot), poor,
                          float(np.mean(f < 1 - 1e-9))))
        mk[lab] = dict(gini=g, sample_gini=g_s, rows=rows, curve=curve)
        assert max(x["residual"] for x in rows) < 1e-6, "market prices did not converge"
    print("      Columns: share of people with at least one need not fully met; lowest share")
    print("      of need met by anyone; average share met by the 10 % with the lowest incomes;")
    print("      total need met; share of capacity spent on uses beyond anyone's need; share")
    print("      of capacity that would close every remaining gap. Unused capacity is 0 by")
    print("      construction (all income is spent, so all capacity is bought).")
    print()
    print("  (c) POLYCENTRIC NETWORK (section 1b) with spare capacity - same code, own seed")
    poly = []
    for cap in (0.95, 1.10, 1.25):
        r = polycentric_demo(np.random.default_rng(SEED + 6), cap=cap, verbose=False)
        poly.append(dict(cap=cap, dec=r["dec"], orc=r["orc"], ratio=r["ratio"]))
        print(f"      capacity {cap:.2f}x mean need: decentralised {r['dec']:.1%}, "
              f"oracle {r['orc']:.1%}, ratio {r['ratio']:.1%}")
    print("      (0.95x repeats section 1b's setting with a different random seed.)")
    print("  -> Under sufficiency the needs-based allocation meets every need with capacity")
    print("     to spare, while allocation by purchasing power leaves needs unmet even though")
    print("     the capacity to meet them exists and is being used for other things.")
    return dict(levels=list(REGIME_LEVELS), needs_based=nb, market=mk, poly=poly)


# ---------------------------------------------------------------------------
# 2. Numerical check of L3 with a Pareto top tail
# ---------------------------------------------------------------------------

def l3_numerical_check(rng):
    print()
    print("=" * 76)
    print("2. L3 CHECK: log-normal closed form vs log-normal body + Pareto top 1 % (same Gini)")
    print("=" * 76)
    rows = []
    z = rng.standard_normal(400_000)
    u = rng.random(400_000)

    def tailed(sig):
        # log-normal body; top 1 % replaced by a Pareto tail (alpha = 1.6, roughly the
        # US top-income tail; Atkinson, Piketty & Saez 2011)
        x = np.exp(sig * z)
        cut = np.quantile(x, 0.99)
        top = x >= cut
        x[top] = cut * (1 - u[top]) ** (-1 / 1.6)
        return x

    for g in (0.30, 0.35, 0.39):
        sig = float(sigma_from_gini(g))
        # re-fit the body so the tailed distribution has the SAME overall Gini
        lo_s, hi_s = 0.05, sig
        for _ in range(40):
            mid = (lo_s + hi_s) / 2
            if gini_sample(tailed(mid)) > g:
                hi_s = mid
            else:
                lo_s = mid
        x = tailed((lo_s + hi_s) / 2)
        for e in (0.7, 1.3, 2.0):
            cf = float(atkinson_lognormal(sig, e))
            num = atkinson_sample(x, e)
            rows.append((g, e, cf, num, gini_sample(x)))
            print(f"  Gini {g:.2f}  e={e:.1f}:  log-normal {cf:6.1%}   log-normal + Pareto tail, same Gini {num:6.1%}")
    print("  -> At the same overall Gini, adding a heavy top tail changes L3 by about two")
    print("     percentage points or less (up for low curvature, down for high curvature).")
    print("     The log-normal closed form used in the model is a good approximation.")
    return rows


# ---------------------------------------------------------------------------
# 3. Central estimate, Monte Carlo, sensitivity
# ---------------------------------------------------------------------------

def main(make_charts=True):
    rng = np.random.default_rng(SEED)
    feas = feasibility_demo(rng)
    poly = polycentric_demo(rng)
    l3rows = l3_numerical_check(rng)

    keys = list(P.keys())
    cen = {k: central(k) for k in keys}
    ch_c = channels(cen)
    out_c, wb_c = net_advantage(ch_c)

    print()
    print("=" * 76)
    print("3. CENTRAL ESTIMATES (share of total productive capacity K)")
    print("=" * 76)
    print(f"  Gross money-handling overhead (before residual)  {ch_c['gross_overhead']:6.1%}")
    for k, lab in [("L1", "L1 Overhead (net of residual coordination)"),
                   ("L2", "L2 Access-gating (idle capacity)"),
                   ("L3", "L3 Allocation by income, not need (wellbeing-equiv.)"),
                   ("L4", "L4 Market power"),
                   ("L5", "L5 Duplication & advertising-induced demand"),
                   ("L6", "L6 Premature obsolescence"),
                   ("L7", "L7 Market mismatch (unsold, returns, failures)")]:
        print(f"  {lab:<50}{ch_c[k]:6.1%}")
    for k, lab in [("C_coord", "C_coord Local coordination labour"),
                   ("C_fed", "C_fed Federation (cross-node) labour"),
                   ("C_local", "C_local Local mismatch"),
                   ("C_agg", "C_agg Cross-node aggregation error"),
                   ("C_dyn", "C_dyn Slower discovery"),
                   ("C_inc", "C_inc Participation shortfall"),
                   ("C_err", "C_err Need-assessment error (wellbeing-equiv.)")]:
        print(f"  {lab:<50}{ch_c[k]:6.1%}")
    print(f"  Net advantage of currency-free, OUTPUT ONLY      {out_c:+6.1%}")
    print(f"  Net advantage of currency-free, WELLBEING        {wb_c:+6.1%}")

    # ---- Monte Carlo ----
    draws = {k: sample(k, rng, N_DRAWS) for k in keys}
    ch = channels(draws)
    adv_out, adv_wb = net_advantage(ch)
    # Variant: currency economy with low inequality (Nordic-like Gini 0.26)
    d_low = dict(draws); d_low["gini"] = np.full(N_DRAWS, LOW_INEQ_GINI)
    adv_out_low, adv_wb_low = net_advantage(channels(d_low))
    # Variant: high-inequality currency economy (Brazil-/Colombia-like), drawn AFTER all
    # other draws so the main results are unchanged by its addition
    d_high = dict(draws); d_high["gini"] = rng.uniform(*HIGH_INEQ_GINI, N_DRAWS)
    adv_out_high, adv_wb_high = net_advantage(channels(d_high))
    L3_high = channels(d_high)["L3"]
    print(f"  L3 in the high-inequality case: median {np.median(L3_high):.1%}"
          f" (90 %: {np.percentile(L3_high,5):.1%}-{np.percentile(L3_high,95):.1%})")

    def pct(a, q):
        return np.percentile(a, q)

    print()
    print("=" * 76)
    print(f"4. MONTE CARLO ({N_DRAWS:,} draws, seed {SEED})")
    print("=" * 76)
    for name in ["L1", "L2", "L3"] + OUT_LOSS_KEYS[2:] + COST_KEYS + ["C_err"]:
        a = ch[name]
        print(f"  {name:<8} median {np.median(a):6.1%}   90 % range {pct(a,5):6.1%} - {pct(a,95):6.1%}")
    tot_L_out = sum(ch[k] for k in OUT_LOSS_KEYS)
    tot_C_out = sum(ch[k] for k in COST_KEYS)
    print(f"  Sum L1+L2+L4+L5+L6+L7 median {np.median(tot_L_out):.1%}; sum of output costs median {np.median(tot_C_out):.1%}")
    summary = {}
    for lab, arr in [("Output only (US-like)", adv_out),
                     ("Wellbeing vs high-inequality currency economy (Gini 0.50-0.55)", adv_wb_high),
                     ("Wellbeing (US-like inequality)", adv_wb),
                     ("Wellbeing vs low-inequality currency economy (Gini 0.26)", adv_wb_low)]:
        share = float(np.mean(arr > 0))
        summary[lab] = (share, float(np.median(arr)), float(pct(arr, 5)), float(pct(arr, 95)))
        print(f"  {lab:<58} currency-free ahead in {share:6.1%} of draws;"
              f" median {np.median(arr):+.1%} (90 %: {pct(arr,5):+.1%} to {pct(arr,95):+.1%})")

    # ---- Rank correlations (Spearman) with wellbeing advantage ----
    def ranks(a):
        r = np.empty(len(a)); r[np.argsort(a)] = np.arange(len(a)); return r
    ry = ranks(adv_wb)
    spear = {k: float(np.corrcoef(ranks(draws[k]), ry)[0, 1]) for k in keys}
    print()
    print("  Spearman rank correlation of each parameter with the WELLBEING advantage:")
    for k, r in sorted(spear.items(), key=lambda kv: -abs(kv[1])):
        print(f"    {k:<10} {r:+.2f}   {P[k][3]}")

    # ---- Tornado: one-at-a-time, 5th / 95th percentile, others central ----
    tornado = []
    for k in keys:
        lo_v, hi_v = np.percentile(draws[k], [5, 95])
        vals = []
        for val in (lo_v, hi_v):
            v = dict(cen); v[k] = val
            vals.append(net_advantage(channels(v))[1])
        tornado.append((k, vals[0], vals[1]))
    tornado.sort(key=lambda t: -abs(t[2] - t[1]))
    tornado_out = []
    for k in keys:
        lo_v, hi_v = np.percentile(draws[k], [5, 95])
        vals = []
        for val in (lo_v, hi_v):
            v = dict(cen); v[k] = val
            vals.append(net_advantage(channels(v))[0])
        tornado_out.append((k, vals[0], vals[1]))
    tornado_out.sort(key=lambda t: -abs(t[2] - t[1]))
    ry_o = ranks(adv_out)
    print()
    print("  OUTPUT-ONLY: Spearman rank correlations (top 8):")
    sp_o = {k: float(np.corrcoef(ranks(draws[k]), ry_o)[0, 1]) for k in keys}
    for k, r in sorted(sp_o.items(), key=lambda kv: -abs(kv[1]))[:8]:
        print(f"    {k:<10} {r:+.2f}   {P[k][3]}")
    print("  OUTPUT-ONLY tornado (top 8):")
    for k, a, b in tornado_out[:8]:
        print(f"    {k:<10} {a:+6.1%}  {b:+6.1%}   swing {abs(b-a):5.1%}")
    print()
    print("  Tornado (wellbeing advantage when one parameter moves to its 5th / 95th pct):")
    for k, a, b in tornado:
        print(f"    {k:<10} {a:+6.1%}  {b:+6.1%}   swing {abs(b-a):5.1%}")

    # ---- Break-even values (others at central values) ----
    print()
    print("  Break-even values (all other parameters central):")
    W_cur_out, W_free_out, W_cur, W_free = welfare(ch_c)
    # participation shortfall at which wellbeing comparison ties
    rest = W_free / (1 - ch_c["C_inc"])
    cinc_star = 1 - W_cur / rest
    print(f"    Wellbeing: currency-free stays ahead while participation shortfall C_inc < {cinc_star:.1%}")
    rest_o = W_free_out / (1 - ch_c["C_inc"])
    cinc_star_o = 1 - W_cur_out / rest_o
    print(f"    Output only: break-even participation shortfall C_inc = {cinc_star_o:+.1%}"
          f" (negative = would need effort to RISE)")
    # Gini at which wellbeing ties
    g_grid = np.linspace(0.05, 0.60, 2201)
    adv_g = []
    for g in g_grid:
        v = dict(cen); v["gini"] = g
        adv_g.append(net_advantage(channels(v))[1])
    adv_g = np.array(adv_g)
    g_star = float(g_grid[np.argmax(adv_g > 0)]) if np.any(adv_g > 0) else float("nan")
    print(f"    Wellbeing: currency-free ahead when resource Gini > {g_star:.3f}"
          f" (US disposable income 0.39; lowest OECD ~0.23-0.26)")
    # total output-side cost at which wellbeing ties
    total_cost_star = 1 - W_cur / (1 - ch_c["C_err"])
    print(f"    Wellbeing: currency-free ahead while combined output costs (all six C terms,"
          f" multiplicative) < {total_cost_star:.1%}; central = {1 - W_free_out:.1%}")

    # ---- Consistency with the site's projections page ----
    print()
    print("  Consistency with /explore/data/projections/ :")
    g = ch["gross_overhead"]
    print(f"    Gross money-handling overhead here: median {np.median(g):.1%}"
          f" (90 %: {pct(g,5):.1%}-{pct(g,95):.1%}); projections page: 13-18 % of the workforce.")
    # net capacity freed = output-only advantage (all output losses avoided minus all output costs)
    W_o_cur = np.ones(N_DRAWS); W_o_free = np.ones(N_DRAWS)
    for k in OUT_LOSS_KEYS:
        W_o_cur = W_o_cur * (1 - ch[k])
    for k in COST_KEYS:
        W_o_free = W_o_free * (1 - ch[k])
    # labour needed to deliver today's useful output: today's hours * W_o_cur / W_o_free
    hours_ratio = W_o_cur / W_o_free
    net_lab = 1 - hours_ratio
    print(f"    Net share of working hours freed (all output channels, all coordination costs): median"
          f" {np.median(net_lab):+.1%} (90 %: {pct(net_lab,5):+.1%} to {pct(net_lab,95):+.1%}).")
    base_h = 38.6   # ASSUMPTION: average usual weekly hours of US workers (BLS CPS, ~38-39)
    print(f"    Implied work week for today's workers: median {base_h*np.median(hours_ratio):.1f} h"
          f" (90 %: {base_h*pct(hours_ratio,5):.1f}-{base_h*pct(hours_ratio,95):.1f} h) from {base_h} h.")
    print("    -> A 15-17 hour week would need a ~56-61 % cut in hours per worker; channels modelled")
    print("       here do not reach that. Any remaining gap would have to come from sharing work")
    print("       among more adults and from channels not modelled here.")

    # ---- Two regimes (own random generator: changes none of the numbers above) ----
    regimes = regimes_demo()

    if make_charts:
        write_charts(ch_c, ch, adv_out, adv_wb, adv_wb_low, tornado, tornado_out, summary, adv_wb_high,
                     regimes)
    return dict(central=ch_c, summary=summary, tornado=tornado, spear=spear,
                cinc_star=cinc_star, g_star=g_star, feas=feas, poly=poly, l3=l3rows,
                regimes=regimes)


# ---------------------------------------------------------------------------
# 4. Charts (SVG, colour-blind-safe Okabe-Ito blue/vermillion, light surface)
# ---------------------------------------------------------------------------

LABELS = {
    "fin": "Finance employment share", "hlth": "Excess health admin", "tax": "Tax-compliance labour",
    "adv": "Advertising", "other": "Other money handling", "rho": "Residual share still needed",
    "food_gate": "Food lost to price gating", "house_use": "Usable vacant homes",
    "idle": "Idle capacity (output gap)", "gini": "Inequality (Gini)",
    "eps": "Wellbeing curvature (ε)", "L4": "Market-power loss",
    "mkt_ind": "Advertising-induced demand", "dup": "Competitive duplication",
    "life_ext": "Achievable lifetime extension", "obs_resid": "Repair/durability offset",
    "L7": "Unsold, returned, failed output",
    "C_coord": "Local coordination labour", "C_fed": "Federation labour",
    "C_local": "Local mismatch", "C_agg": "Cross-node aggregation error", "C_dyn": "Slower discovery",
    "C_inc": "Participation shortfall", "s_err": "Need-assessment error",
}


def write_charts(ch_c, ch, adv_out, adv_wb, adv_wb_low, tornado, tornado_out, summary, adv_wb_high,
                 regimes=None):
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    # Inside the website's repository (public/downloads/models/) the charts go to
    # public/images/models/; anywhere else, to ./money-model-charts/ in the
    # directory you run the script from.
    here = os.path.dirname(os.path.abspath(__file__))
    root = os.path.abspath(os.path.join(here, "..", "..", ".."))
    if os.path.basename(here) == "models" and os.path.isdir(os.path.join(root, "public", "images")):
        outdir = os.path.join(root, "public", "images", "models")
    else:
        outdir = os.path.join(os.getcwd(), "money-model-charts")
    os.makedirs(outdir, exist_ok=True)

    BLUE, VERM, INK, MUTED, GRID, BG = "#0072B2", "#D55E00", "#1f1f1f", "#555555", "#dddddd", "#fcfcfb"
    plt.rcParams.update({
        "font.family": "DejaVu Sans", "font.size": 11, "svg.fonttype": "none",
        "axes.edgecolor": MUTED, "axes.labelcolor": INK, "xtick.color": MUTED, "ytick.color": MUTED,
        "axes.spines.top": False, "axes.spines.right": False,
        "figure.facecolor": BG, "axes.facecolor": BG, "savefig.facecolor": BG,
        # fixed element ids and no timestamp: re-running writes byte-identical SVG files
        "svg.hashsalt": "money-vs-moneyless-model",
    })
    META = {"Date": None}

    # ---- Chart 1: channel breakdown (medians with 90 % ranges) ----
    items = [("L1 Money-handling overhead", "L1", BLUE), ("L2 Idle capacity (access-gating)", "L2", BLUE),
             ("L3 Allocation by income, not need*", "L3", BLUE), ("L4 Market power", "L4", BLUE),
             ("L5 Duplication & ad-induced demand", "L5", BLUE), ("L6 Premature obsolescence", "L6", BLUE),
             ("L7 Unsold, returned, failed output", "L7", BLUE),
             ("Local coordination labour", "C_coord", VERM), ("Federation (cross-node) labour", "C_fed", VERM),
             ("Local mismatch", "C_local", VERM), ("Cross-node aggregation error", "C_agg", VERM),
             ("Slower discovery", "C_dyn", VERM), ("Participation shortfall", "C_inc", VERM),
             ("Need-assessment error*", "C_err", VERM)]
    fig, ax = plt.subplots(figsize=(9, 7.8))
    y = np.arange(len(items))[::-1].astype(float)
    y[:7] += 0.6  # gap between groups
    for yi, (lab, k, col) in zip(y, items):
        med = np.median(ch[k]) * 100
        lo, hi = np.percentile(ch[k], [5, 95]) * 100
        ax.barh(yi, med, color=col, height=0.62)
        ax.plot([lo, hi], [yi, yi], color=INK, lw=1.5)
        ax.plot([lo, lo], [yi - 0.15, yi + 0.15], color=INK, lw=1.5)
        ax.plot([hi, hi], [yi - 0.15, yi + 0.15], color=INK, lw=1.5)
        ax.text(hi + 0.6, yi, f"{med:.1f}%", va="center", color=INK, fontsize=10)
    ax.set_yticks(y)
    ax.set_yticklabels([i[0] for i in items])
    ax.set_xlabel("Share of total productive capacity (%), median and 90% range")
    ax.xaxis.grid(True, color=GRID, lw=0.8); ax.set_axisbelow(True)
    ax.text(0, y[0] + 0.75, "Losses of currency-based allocation", color=BLUE, fontweight="bold")
    ax.text(0, y[7] + 0.75, "Costs of currency-free coordination", color=VERM, fontweight="bold")
    ax.set_title("Where each system loses capacity", loc="left", color=INK, fontsize=13, pad=22)
    fig.text(0.01, 0.01, "* wellbeing-equivalent (Atkinson index, same wellbeing curvature on both sides). "
             "Monte Carlo, 200,000 draws.", fontsize=8.5, color=MUTED)
    fig.tight_layout(rect=(0, 0.03, 1, 1))
    fig.savefig(os.path.join(outdir, "money-channels.svg"), metadata=META)
    plt.close(fig)

    # ---- Chart 2: tornado, wellbeing (top) and output-only (bottom) ----
    from matplotlib.patches import Patch
    fig, axes = plt.subplots(2, 1, figsize=(9, 8.6))
    panels = [(axes[0], tornado, net_advantage(ch_c)[1] * 100, "Wellbeing comparison (vs US-like economy)"),
              (axes[1], tornado_out, net_advantage(ch_c)[0] * 100, "Output-only comparison")]
    for ax, tor, base, title in panels:
        top = tor[:8][::-1]
        for i, (k, a_, b_) in enumerate(top):
            ax.barh(i, a_ * 100 - base, left=base, color=VERM, height=0.6)
            ax.barh(i, b_ * 100 - base, left=base, color=BLUE, height=0.6)
        ax.axvline(base, color=INK, lw=1)
        ax.axvline(0, color=MUTED, lw=1, ls="--")
        ax.text(base, len(top) - 0.35, f" central {base:+.1f}%", color=INK, fontsize=9)
        ax.text(0, -0.85, " break-even", color=MUTED, fontsize=9)
        ax.set_yticks(range(len(top)))
        ax.set_yticklabels([LABELS[t[0]] for t in top])
        ax.set_ylim(-1.1, len(top))
        ax.xaxis.grid(True, color=GRID, lw=0.8); ax.set_axisbelow(True)
        ax.set_title(title, loc="left", color=INK, fontsize=11)
    axes[1].set_xlabel("Advantage of currency-free allocation (%), one parameter varied at a time")
    axes[0].legend(handles=[Patch(color=BLUE, label="parameter at its 95th percentile"),
                            Patch(color=VERM, label="parameter at its 5th percentile")],
                   loc="upper left", frameon=False, fontsize=9)
    fig.suptitle("Which assumptions matter most", x=0.01, ha="left", fontsize=13, color=INK)
    fig.tight_layout()
    fig.savefig(os.path.join(outdir, "money-tornado.svg"), metadata=META)
    plt.close(fig)

    # ---- Chart 3: Monte Carlo distributions ----
    fig, axes = plt.subplots(4, 1, figsize=(9, 9.2), sharex=True)
    bins = np.linspace(-30, 150, 361)
    sets = [("Output only (goods and labour, ignoring who receives them)", adv_out, "Output only (US-like)"),
            ("Wellbeing vs a HIGH-inequality currency economy (Gini 0.50-0.55)", adv_wb_high,
             "Wellbeing vs high-inequality currency economy (Gini 0.50-0.55)"),
            ("Wellbeing vs a US-LIKE currency economy (Gini 0.30-0.40)", adv_wb, "Wellbeing (US-like inequality)"),
            ("Wellbeing vs a LOW-inequality currency economy (Gini 0.26)", adv_wb_low,
             "Wellbeing vs low-inequality currency economy (Gini 0.26)")]
    for ax, (title, arr, key) in zip(axes, sets):
        a = arr * 100
        beyond = float(np.mean(a > bins[-1]))
        h, e = np.histogram(a, bins=bins)   # values beyond the axis are not drawn
        centers = (e[:-1] + e[1:]) / 2
        w = e[1] - e[0]
        ax.bar(centers[centers < 0], h[centers < 0], width=w * 0.92, color=BLUE)
        ax.bar(centers[centers >= 0], h[centers >= 0], width=w * 0.92, color=VERM)
        ax.axvline(0, color=INK, lw=1)
        share = summary[key][0]
        ax.set_title(title, loc="left", fontsize=10.5, color=INK)
        fmt = lambda x: f"{x:.1%}" if (0 < x < 0.01 or 0.99 < x < 1) else f"{x:.0%}"
        ax.text(0.99, 0.85, f"currency-free ahead: {fmt(share)}", transform=ax.transAxes,
                ha="right", color=VERM, fontsize=10, fontweight="bold")
        ax.text(0.01, 0.85, f"currency ahead: {fmt(1-share)}", transform=ax.transAxes,
                ha="left", color=BLUE, fontsize=10, fontweight="bold")
        if beyond > 0:
            ax.text(0.01, 0.62, f"{beyond:.1%} of scenarios lie beyond +{bins[-1]:.0f}% (not drawn)",
                    transform=ax.transAxes, ha="left", color=MUTED, fontsize=8.5)
            print(f"  Chart note: {beyond:.2%} of '{key}' draws beyond +{bins[-1]:.0f}%")
        ax.set_yticks([])
        ax.spines["left"].set_visible(False)
    axes[-1].set_xlabel("Advantage of currency-free allocation (% of welfare-equivalent output)")
    fig.suptitle("How often each system comes out ahead across 200,000 plausible scenarios",
                 x=0.01, ha="left", fontsize=13, color=INK)
    fig.tight_layout()
    fig.savefig(os.path.join(outdir, "money-monte-carlo.svg"), metadata=META)
    plt.close(fig)

    # ---- Chart 4: two regimes - share of need met vs capacity level ----
    names = ["money-channels.svg", "money-tornado.svg", "money-monte-carlo.svg"]
    if regimes is not None:
        SKY, NAVY = "#56B4E9", "#003F63"
        fig, axes = plt.subplots(2, 1, figsize=(9, 7.6), sharex=True)
        xs = np.round(np.arange(0.80, 3.0001, 0.05), 2)
        nb_line = np.minimum(xs, 1.0) * 100     # equals the LP optimum (asserted in regimes_demo)
        styles = [("Low inequality (Gini 0.26)", SKY, "Purchasing power, low inequality (Gini 0.26)"),
                  ("US-like (Gini 0.35)", BLUE, "Purchasing power, US-like inequality (Gini 0.35)"),
                  ("High inequality (Gini 0.525)", NAVY, "Purchasing power, high inequality (Gini 0.525)")]
        panels = [(axes[0], 1, "Total need met (%)", (50, 103)),
                  (axes[1], 2, "Share of need met for the 10% with the lowest incomes (%)", (0, 105))]
        for ax, col_i, title, ylim in panels:
            ax.axvspan(0.8, 1.0, color="#eeeeee", lw=0)
            ax.axvline(1.0, color=MUTED, lw=1, ls="--")
            ax.plot(xs, nb_line, color=VERM, lw=2.6, label="Needs-based coordination (no prices)")
            for lab, colr, legend_lab in styles:
                cur = regimes["market"][lab]["curve"]
                ax.plot([c[0] for c in cur], [c[col_i] * 100 for c in cur], color=colr, lw=2,
                        label=legend_lab)
            ax.set_ylim(*ylim)
            ax.set_xlim(0.8, 3.0)
            ax.yaxis.grid(True, color=GRID, lw=0.8); ax.set_axisbelow(True)
            ax.set_title(title, loc="left", color=INK, fontsize=11)
        axes[0].legend(loc="lower right", frameon=False, fontsize=9)
        axes[1].text(0.81, 2, "shortage", color=MUTED, fontsize=9)
        axes[1].text(1.03, 2, "sufficiency: capacity \u2265 need", color=MUTED, fontsize=9)
        axes[1].set_xlabel("Capacity, as a multiple of what meeting every need requires")
        axes[1].set_xticks([0.8, 1.0, 1.2, 1.5, 2.0, 2.5, 3.0])
        axes[1].set_xticklabels(["0.8\u00d7", "1\u00d7", "1.2\u00d7", "1.5\u00d7", "2\u00d7",
                                 "2.5\u00d7", "3\u00d7"])
        fig.suptitle("Same needs, same capacity: allocation by need vs by purchasing power",
                     x=0.01, ha="left", fontsize=13, color=INK)
        fig.text(0.01, 0.01, "Blue lines: budget-constrained market (full employment, no market power, no public "
                 "provision). 2,000 people, 50 goods, 4 resources.",
                 fontsize=8.5, color=MUTED)
        fig.tight_layout(rect=(0, 0.03, 1, 1))
        fig.savefig(os.path.join(outdir, "money-regimes.svg"), metadata=META)
        plt.close(fig)
        names.append("money-regimes.svg")
    # let browsers fall back to a sans-serif font if DejaVu Sans is not installed
    for name in names:
        fp = os.path.join(outdir, name)
        with open(fp, encoding="utf-8") as fh:
            svg = fh.read()
        svg = svg.replace("font-family: 'DejaVu Sans'",
                          "font-family: 'DejaVu Sans', 'Helvetica Neue', Arial, sans-serif")
        with open(fp, "w", encoding="utf-8") as fh:
            fh.write(svg)
    print()
    print(f"  Charts written to {outdir}/" + ", ".join(names))


if __name__ == "__main__":
    main(make_charts="--no-charts" not in sys.argv)
