"""Small exact-rational inquiry, not an experimental reproduction. Park-bench height updates follow Paulsen & Keim (2019), section 2(a,b). The iid driving protocol is an explicit variation; alpha is restricted to [0,1]. No simulation sampling, third-party dependencies, network or background work. Run: python grass.py """ from fractions import Fraction as F from math import comb import json def kernel(height, alpha, visit): """Rows: grow by one with alpha*(H-h)/H, then possibly wear by one.""" if height < 1 or not 0 <= alpha <= 1 or not 0 <= visit <= 1: raise ValueError("Require H>=1 and alpha, visit in [0,1]") rows = [] for h in range(height + 1): row = [F(0)] * (height + 1) growth = alpha * F(height - h, height) for grow, pg in ((0, 1 - growth), (1, growth)): for wear, pw in ((0, 1 - visit), (1, visit)): if pg * pw: after = max(0, h + grow - wear) row[after] += pg * pw assert sum(row) == 1 and min(row) >= 0 rows.append(row) return rows def advance(dist, matrix): return [sum(p * row[j] for p, row in zip(dist, matrix)) for j in range(len(dist))] def compose(first, second): return [advance(row, second) for row in first] def stationary(matrix): """Solve pi P = pi and sum(pi)=1; check every residual exactly.""" n = len(matrix) equations = [[matrix[j][i] - int(i == j) for j in range(n)] + [F(0)] for i in range(n - 1)] equations.append([F(1)] * n + [F(1)]) for col in range(n): pivot = next(i for i in range(col, n) if equations[i][col]) equations[col], equations[pivot] = equations[pivot], equations[col] scale = equations[col][col] equations[col] = [v / scale for v in equations[col]] for i in range(n): if i != col: scale = equations[i][col] equations[i] = [a - scale * b for a, b in zip(equations[i], equations[col])] answer = [row[-1] for row in equations] assert sum(answer) == 1 and min(answer) >= 0 assert advance(answer, matrix) == answer return answer def summary(dist): mean = sum(h * p for h, p in enumerate(dist)) var = sum((h - mean) ** 2 * p for h, p in enumerate(dist)) return {"probabilities": list(map(str, dist)), "mean": str(mean), "variance": str(var), "interior_probability": str(sum(dist[1:-1]))} def main(): # H=2, two benches, three patches. Patch 1 is always visited and in # stationarity has height 0 at alpha<=1; patch 3 stays pristine at height 2. # All distributions below concern the middle patch after a complete cycle. H = 2 wear = kernel(H, F(1), F(1)) rest = kernel(H, F(1), F(0)) iid = kernel(H, F(1), F(1, 2)) # Hand-derived H=2 matrices independently fix the update order and edges. assert wear == [[F(1), F(0), F(0)], [F(1, 2), F(1, 2), F(0)], [F(0), F(1), F(0)]] assert rest == [[F(0), F(1), F(0)], [F(0), F(1, 2), F(1, 2)], [F(0), F(0), F(1)]] iid_pi = stationary(iid) before_visit = stationary(compose(wear, rest)) after_visit = advance(before_visit, wear) phase_average = [(a + b) / 2 for a, b in zip(before_visit, after_visit)] assert advance(after_visit, rest) == before_visit # Independent closed-form oracle for iid alpha=1, for several sizes/rates. cases = 0 for hmax in (1, 2, 3, 5, 10): for q in (F(1, 3), F(1, 2), F(2, 3)): expected = [F(comb(hmax, h)) * (1 - q) ** h * q ** (hmax - h) for h in range(hmax + 1)] assert stationary(kernel(hmax, F(1), q)) == expected cases += 1 no_noise = stationary(kernel(H, F(0), F(1, 2))) weaker_noise = stationary(kernel(H, F(1, 2), F(1, 2))) assert no_noise == [F(1), F(0), F(0)] assert iid_pi == [F(1, 4), F(1, 2), F(1, 4)] assert phase_average == [F(1, 6), F(2, 3), F(1, 6)] # Stop all visits but retain regrowth. Contrast saturated grass with a # previously trained distribution. Their mean difference must decay. paused = iid_pi[:] pause_rows = [] for t in range(21): mean = sum(h * p for h, p in enumerate(paused)) difference = H - mean assert difference == F(1, 2) ** t if t in (0, 1, 2, 5, 10, 20): pause_rows.append({"growth_cycles": t, "mean": str(mean), "difference_from_pristine": str(difference)}) paused = advance(paused, rest) # A different continuation: visit only bench 2, without changing noise. new_protocol = iid_pi[:] for _ in range(20): new_protocol = advance(new_protocol, wear) # A separate, deterministic variant: add the expected regrowth as a real # height instead of randomly adding an integer. This tests whether the # restoring part alone suffices to keep an intermediate grass height. def deterministic(h, visited): return max(F(0), h + (H - h) / H - int(visited)) deterministic_before = F(4, 3) deterministic_after = deterministic(deterministic_before, True) assert deterministic_after == F(2, 3) assert deterministic(deterministic_after, False) == deterministic_before deterministic_start = F(2) for _ in range(20): deterministic_start = deterministic(deterministic_start, True) deterministic_start = deterministic(deterministic_start, False) assert deterministic_start - deterministic_before == F(2, 3) * F(1, 4) ** 20 output = { "model": "middle grass patch, H=2, growth before driving", "deterministic_regrowth_control": { "after_visit": str(deterministic_after), "after_no_visit": str(deterministic_before), "twenty_period_residual": str(deterministic_start - deterministic_before), "note": "Real-valued deterministic variant, not the stochastic source model." }, "iid_no_regrowth": summary(no_noise), "iid_alpha_half": summary(weaker_noise), "iid_alpha_one": summary(iid_pi), "alternating_alpha_one_after_visit": summary(after_visit), "alternating_alpha_one_after_no_visit": summary(before_visit), "alternating_alpha_one_phase_average": summary(phase_average), "stop_visits_keep_growth": pause_rows, "switch_to_only_bench_two_after_20": summary(new_protocol), "binomial_oracle_cases": cases, "checks": "All transition rows, stationary residuals, periodic phase closure, binomial oracles and pause mean identities exact." } print(json.dumps(output, indent=2)) if __name__ == "__main__": main()