# Chapter 8 trial. Part 1: the simulated resource game, played by two simple rule-following
# bots (a MODEL of the hypothesis, not data about people). Part 2: sample size.
# Hidden stock dynamics per round (participants are not told the model or its fold):
#   x <- x + 0.5 [x (1 - x/10) - E x^2/(x^2 + 1)] + 0.1 N(0,1); fold at E = 2.604.
# Profit per round = catch - 0.2 E. Collapse (primary outcome): x < 2.0 in any round; by the
# game's rules the fishery is then closed for good. 80 rounds; effort E reviewed every 4 rounds.
import math, random
from statistics import NormalDist
rng = random.Random(16)

def play(dashboard, rounds=80, heed=0.3):
    """A bot that nudges effort up or down by its own step, reversing when the block's profit
    falls, and stops exploring once profit per round reaches its own aspiration. The dashboard
    bot also looks at lag-1 autocorrelation of the last 20 stock values (with probability heed
    at each review) and cuts effort by 0.3 if it exceeds 0.6."""
    x, E, step = 8.9, rng.uniform(0.5, 2.0), rng.uniform(0.05, 0.3)
    asp, stopped, last, hold = rng.uniform(0.6, 0.9), False, None, 0
    xs, block, total = [], 0.0, 0.0
    for t in range(1, rounds + 1):
        catch = 0.5 * E * x * x / (x * x + 1)
        x += 0.5 * (x * (1 - x / 10)) - catch + 0.1 * rng.gauss(0, 1)
        xs.append(x); block += catch - 0.2 * E; total += catch - 0.2 * E
        if x < 2.0: return True, total
        if t % 4: continue
        if not stopped and block / 4 >= asp: stopped = True
        if last is not None and block < last: step = -step       # profit fell: reverse
        last, block = block, 0.0
        if dashboard and len(xs) >= 20 and rng.random() < heed:
            w = xs[-20:]; m = sum(w) / 20; v = sum((u - m) ** 2 for u in w)
            a = sum((w[i] - m) * (w[i + 1] - m) for i in range(19)) / v
            if a > 0.6: E, step, hold = max(0.0, E - 0.3), -abs(step), 3
        if hold: hold -= 1; continue
        if not stopped: E = max(0.0, E + step)
    return False, total

for dash, lab in ((False, "stock and profit only"), (True, "plus recovery dashboard")):
    res = [play(dash) for _ in range(2000)]
    print(f"bots, {lab:24s}: collapse in {sum(c for c, _ in res) / 2000:.0%} of games;"
          f" mean profit {sum(p for _, p in res) / 2000:.1f}")

def n_per_arm(p1, p2, alpha=0.05, power=0.90):
    z = NormalDist().inv_cdf; pb = (p1 + p2) / 2
    num = (z(1 - alpha / 2) * math.sqrt(2 * pb * (1 - pb))
           + z(power) * math.sqrt(p1 * (1 - p1) + p2 * (1 - p2)))
    return math.ceil(num ** 2 / (p1 - p2) ** 2)

print("assumed collapse share: control -> dashboard   analysable per arm (90 % power)")
for p1, p2 in ((0.40, 0.25), (0.40, 0.28), (0.50, 0.35), (0.30, 0.18)):
    print(f"                        {p1:.2f} -> {p2:.2f}            {n_per_arm(p1, p2)}")
n = n_per_arm(0.40, 0.25); rec = math.ceil(n / 0.9)
print(f"chosen 0.40 -> 0.25: {n} per arm; allowing 10 % incomplete sessions,")
print(f"  recruit {rec} per arm, {2 * rec} in all")
hits = 0
for _ in range(4000):
    a = sum(rng.random() < 0.40 for _ in range(n)); b = sum(rng.random() < 0.25 for _ in range(n))
    pa, pb = a / n, b / n; pp = (a + b) / (2 * n)
    hits += abs(pa - pb) / math.sqrt(pp * (1 - pp) * 2 / n) > 1.96
print(f"simulation check of power with {n} per arm: {hits / 4000:.2f}")
# bots, stock and profit only   : collapse in 40% of games; mean profit 42.3
# bots, plus recovery dashboard : collapse in 7% of games; mean profit 35.3
# assumed collapse share: control -> dashboard   analysable per arm (90 % power)
#                         0.40 -> 0.25            203
#                         0.40 -> 0.28            326
#                         0.50 -> 0.35            227
#                         0.30 -> 0.18            265
# chosen 0.40 -> 0.25: 203 per arm; allowing 10 % incomplete sessions,
#   recruit 226 per arm, 452 in all
# simulation check of power with 203 per arm: 0.90

# Decision rule for the Chapter 8 trial, fixed before any data are collected.
# diff = collapse share in the control arm (stock and profit only) minus the dashboard arm
# (plus rolling recovery time and variance); POSITIVE diff = fewer collapses with the dashboard.
# 95 % CI from the pre-specified model (logistic regression, arm + site stratum, reported as a
# risk difference). Smallest difference of interest: 0.10 (10 percentage points).
SESOI = 0.10

def verdict(diff, lo, hi, shown, profit_change):
    if shown < 0.90:
        return (f"diff {diff:+.2f}: NOT INTERPRETABLE - dashboard shown as specified"
                f" in only {shown:.0%} of sessions")
    if lo >= SESOI:    v = "dashboard cuts collapses by >= 10 points"
    elif lo > 0:       v = "dashboard cuts collapses; >= 10 points not shown"
    elif hi < 0:       v = "dashboard INCREASES collapses"
    elif lo > -SESOI and hi < SESOI: v = "no difference worth acting on"
    else:              v = "unclear: the interval is too wide to decide"
    return f"diff {diff:+.2f} ({lo:+.2f} to {hi:+.2f}): {v}; profit {profit_change:+.0%}"

for case in ((0.19, 0.11, 0.27, 0.96, -0.08), (0.08, 0.01, 0.15, 0.95, -0.03),
             (0.02, -0.05, 0.09, 0.94, 0.00), (-0.09, -0.17, -0.01, 0.95, 0.05),
             (0.21, 0.13, 0.29, 0.80, -0.10), (0.05, -0.08, 0.18, 0.93, -0.02)):
    print(verdict(*case))
print("profit and time to collapse are secondary: reported with intervals, never used to")
print("  override the primary verdict")
# diff +0.19 (+0.11 to +0.27): dashboard cuts collapses by >= 10 points; profit -8%
# diff +0.08 (+0.01 to +0.15): dashboard cuts collapses; >= 10 points not shown; profit -3%
# diff +0.02 (-0.05 to +0.09): no difference worth acting on; profit +0%
# diff -0.09 (-0.17 to -0.01): dashboard INCREASES collapses; profit +5%
# diff +0.21: NOT INTERPRETABLE - dashboard shown as specified in only 80% of sessions
# diff +0.05 (-0.08 to +0.18): unclear: the interval is too wide to decide; profit -2%
# profit and time to collapse are secondary: reported with intervals, never used to
#   override the primary verdict