#!/usr/bin/env python3 """Run the precommitted fixed-capacity queue experiment using Python stdlib only. Usage: python3 run_experiment.py --plan plan.json --output-dir rerun-results The output directory may exist, but computed result files must not already exist. This code makes no network requests and does not read external scientific data. """ import argparse import csv from datetime import datetime, timezone from fractions import Fraction import hashlib import json import math from pathlib import Path import platform import random import statistics import sys import time def sha(data): return hashlib.sha256(data).hexdigest() def save_json(path, data): path.write_text(json.dumps(data, indent=2) + '\n', encoding='utf-8') def save_csv(path, rows): with path.open('w', newline='', encoding='utf-8') as stream: writer = csv.DictWriter(stream, fieldnames=list(rows[0])) writer.writeheader() writer.writerows(rows) def deterministic(plan): rows = [] summaries = [] cycles = plan['deterministic']['cycles'] for name, pattern in [('constant_fluid', [Fraction(4, 5)]), ('periodic_burst', list(map(Fraction, [4, 0, 0, 0, 0])))]: q = Fraction(0) arrivals, queues = [], [] for t in range(1, cycles + 1): arrival = pattern[(t - 1) % len(pattern)] q = max(Fraction(0), q + arrival - 1) arrivals.append(arrival) queues.append(q) rows.append({'pattern': name, 'cycle': t, 'arrival_exact': str(arrival), 'arrival': float(arrival), 'after_service_queue_exact': str(q), 'after_service_queue': float(q)}) mean_a = sum(arrivals) / cycles variance_a = sum((a - mean_a) ** 2 for a in arrivals) / cycles summaries.append({'pattern': name, 'cycles': cycles, 'mean_arrival_exact': str(mean_a), 'mean_arrival': float(mean_a), 'arrival_variance_exact': str(variance_a), 'arrival_variance': float(variance_a), 'mean_queue_exact': str(sum(queues) / cycles), 'mean_queue': float(sum(queues) / cycles), 'max_queue': float(max(queues)), 'end_queue': float(q)}) maps = [] for initial in plan['deterministic']['period_map_initial_queues']: q = Fraction(initial) for arrival in [4, 0, 0, 0, 0]: q = max(Fraction(0), q + arrival - 1) expected = max(Fraction(0), Fraction(initial) - 1) if q != expected: raise AssertionError('Exact periodic map does not match its derivation') maps.append({'initial_queue': initial, 'queue_after_five_cycles': float(q), 'expected_max_0_q_minus_1': float(expected)}) if summaries[0]['mean_queue'] != 0 or summaries[1]['mean_queue_exact'] != '6/5': raise AssertionError('Exact deterministic control failed') return rows, summaries, maps def simulate(rho_string, batch, replicate, plan): cfg = plan['stochastic'] rho = float(rho_string) probability = rho / batch if not 0 <= probability <= 1: raise ValueError('Infeasible Bernoulli probability') seed_text = f'congestion-v1|rho={rho_string}|batch={batch}|replicate={replicate}' seed = int.from_bytes(hashlib.sha256(seed_text.encode()).digest()[:8], 'big') rng = random.Random(seed) draw = rng.random horizon = cfg['cycles_per_replicate'] burn = cfg['burn_in_cycles'] measured = horizon - burn split = burn + measured // 2 every = cfg['checkpoint_every_cycles'] q = total_arrivals = sum_q = sum_a = early_sum = late_sum = 0 above10 = above100 = empty = max_q = half_queue = 0 checkpoints = [{'rho': rho_string, 'batch': batch, 'replicate': replicate, 'cycle': 0, 'queue': 0, 'cumulative_arrivals': 0}] for t in range(1, horizon + 1): arrival = batch if draw() < probability else 0 total_arrivals += arrival q += arrival if q > 0: q -= 1 if t == horizon // 2: half_queue = q if t > burn: sum_q += q sum_a += arrival above10 += q >= 10 above100 += q >= 100 empty += q == 0 if q > max_q: max_q = q if t <= split: early_sum += q else: late_sum += q if t % every == 0: checkpoints.append({'rho': rho_string, 'batch': batch, 'replicate': replicate, 'cycle': t, 'queue': q, 'cumulative_arrivals': total_arrivals}) early = early_sum / (measured // 2) late = late_sum / (measured // 2) result = {'rho': rho_string, 'batch': batch, 'replicate': replicate, 'seed': seed, 'seed_text': seed_text, 'arrival_probability': probability, 'regime': 'below_capacity' if rho < 1 else 'overload', 'mean_queue': sum_q / measured, 'tail_ge_10': above10 / measured, 'tail_ge_100': above100 / measured, 'empty_fraction': empty / measured, 'measured_mean_arrivals': sum_a / measured, 'max_queue': max_q, 'queue_at_half_horizon': half_queue, 'queue_at_end': q, 'late_endpoint_growth': (q - half_queue) / (horizon / 2), 'early_half_mean_queue': early, 'late_half_mean_queue': late, 'half_window_growth': (late - early) / (measured / 2)} return result, checkpoints def aggregate(rows): metrics = ['mean_queue', 'tail_ge_10', 'tail_ge_100', 'empty_fraction', 'measured_mean_arrivals', 'max_queue', 'queue_at_end', 'late_endpoint_growth', 'early_half_mean_queue', 'late_half_mean_queue', 'half_window_growth'] first = rows[0] rho = Fraction(first['rho']) batch = first['batch'] expected = rho * (batch - 1) / (2 * (1 - rho)) if rho < 1 else None result = {'rho': first['rho'], 'batch': batch, 'replicates': len(rows), 'regime': first['regime'], 'arrival_probability': first['arrival_probability'], 'analytic_stationary_mean_queue': float(expected) if expected is not None else None, 'analytic_stationary_mean_queue_exact': str(expected) if expected is not None else None, 'theoretical_arrival_variance': float(rho * batch - rho * rho)} for metric in metrics: values = [r[metric] for r in rows] mean = statistics.mean(values) standard_error = statistics.stdev(values) / math.sqrt(len(values)) low, high = mean - 2 * standard_error, mean + 2 * standard_error if 'growth' not in metric: low = max(0, low) if metric in ['tail_ge_10', 'tail_ge_100', 'empty_fraction']: high = min(1, high) result[metric + '_mean'] = mean result[metric + '_se'] = standard_error result[metric + '_ci_low'] = low result[metric + '_ci_high'] = high return result def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument('--plan', type=Path, default=Path(__file__).with_name('plan.json')) parser.add_argument('--output-dir', type=Path, default=Path(__file__).with_name('rerun-results')) args = parser.parse_args() out = args.output_dir.resolve() out.mkdir(parents=True, exist_ok=True) if (out / 'results.json').exists() or (out / 'replicate-results.csv').exists(): raise SystemExit('Computed output already exists; choose a new directory.') plan_bytes = args.plan.read_bytes() plan = json.loads(plan_bytes) start = datetime.now(timezone.utc).isoformat() clock = time.monotonic() controls, control_summary, maps = deterministic(plan) save_csv(out / 'deterministic-trace.csv', controls) save_csv(out / 'period-map-check.csv', maps) save_json(out / 'deterministic-results.json', {'patterns': control_summary, 'period_map': maps}) all_replicates, all_checkpoints, summaries, infeasible = [], [], [], [] cfg = plan['stochastic'] for rho in cfg['rho_values']: for batch in cfg['batch_sizes']: if float(rho) / batch > 1: infeasible.append({'rho': rho, 'batch': batch, 'status': 'not_simulated_infeasible_probability', 'arrival_probability': float(rho) / batch}) continue runs = [] for replicate in range(cfg['replicates_per_condition']): result, checkpoints = simulate(rho, batch, replicate, plan) runs.append(result) all_replicates.append(result) all_checkpoints.extend(checkpoints) summary = aggregate(runs) summaries.append(summary) print(f"rho={rho} batch={batch}: mean_Q={summary['mean_queue_mean']:.6f}, " f"late_growth={summary['late_endpoint_growth_mean']:.6f}", flush=True) if len(summaries) != cfg['valid_conditions']: raise AssertionError('Unexpected number of valid conditions') save_csv(out / 'replicate-results.csv', all_replicates) save_csv(out / 'condition-summary.csv', summaries) save_csv(out / 'queue-checkpoints.csv', all_checkpoints) save_json(out / 'replicate-results.json', all_replicates) save_json(out / 'queue-checkpoints.json', all_checkpoints) save_json(out / 'condition-summary.json', summaries) metadata = {'plan_sha256': sha(plan_bytes), 'code_sha256': sha(Path(__file__).read_bytes()), 'started_utc': start, 'finished_utc': datetime.now(timezone.utc).isoformat(), 'elapsed_seconds': time.monotonic() - clock, 'python': sys.version, 'platform': platform.platform(), 'condition_count': len(summaries), 'replicate_count': len(all_replicates), 'simulated_cycles': len(all_replicates) * cfg['cycles_per_replicate'], 'infeasible_conditions': infeasible, 'plan_modified_during_run': sha(args.plan.read_bytes()) != sha(plan_bytes)} save_json(out / 'results.json', {'metadata': metadata, 'deterministic': control_summary, 'period_map': maps, 'conditions': summaries}) save_json(out / 'run-metadata.json', metadata) print(json.dumps(metadata, indent=2), flush=True) if __name__ == '__main__': main()