---
name: street-orientation-entropy
description: >
  Compares street orientation entropy between 10 planned-grid cities and 10 organic cities
  using OpenStreetMap Overpass API data. Uses exact combinatorial testing (all 184,756
  group assignments enumerated), bootstrap confidence intervals, and sensitivity analysis
  across bin sizes and leave-one-out city removal.
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget"
tags: ["claw4s-2026", "urban-science", "entropy", "openstreetmap", "permutation-test", "exact-enumeration"]
python_version: ">=3.8"
dependencies: []
---

# Street Orientation Entropy: Planned vs Organic Cities

## Overview

Downloads road network data for 20 cities (10 grid-planned, 10 organic) from the
OpenStreetMap Overpass API. Computes Shannon entropy of street orientation distributions.
Tests whether planned cities have systematically lower orientation entropy using exact
combinatorial enumeration of all C(20,10) = 184,756 possible group assignments — no
Monte Carlo approximation needed.

**Methodological hook:** Exact enumeration eliminates sampling error in the p-value.
Combined with stdlib-only implementation (no OSMnx/NetworkX), this is fully reproducible
on any Python 3.8+ system with zero dependency installation.

**Data source:** OpenStreetMap via https://overpass-api.de/api/interpreter

## Step 1: Create workspace

```bash
mkdir -p /tmp/claw4s_auto_openstreetmap-city-grid-entropy
```

**Expected output:** Directory created (exit code 0).

## Step 2: Write analysis script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_openstreetmap-city-grid-entropy/script.py
#!/usr/bin/env python3
"""
Street Orientation Entropy: Do Planned Cities Have Lower Entropy Than Organic Cities?

Downloads road networks for 20 cities from OpenStreetMap Overpass API,
computes Shannon entropy of street orientations, and performs exact
combinatorial testing to compare planned vs organic city groups.

Author: Claw, David Austin
License: MIT
"""

import urllib.request
import urllib.error
import urllib.parse
import json
import hashlib
import os
import math
import time
import random
import itertools
import sys

# ─── Configuration ───────────────────────────────────────────────────────────
SEED = 42
N_BINS = 36              # 10-degree bins over 0-180 range
N_BOOTSTRAP = 2000
N_PERM = 3000
OVERPASS_URL = "https://overpass-api.de/api/interpreter"
REQUEST_DELAY = 8        # seconds between API calls (be polite)
REQUEST_TIMEOUT = 180    # seconds
MAX_RETRIES = 5
CACHE_DIR = "cache"
MANIFEST_FILE = os.path.join("cache", "manifest.json")

# ─── City Definitions ────────────────────────────────────────────────────────
# Each bbox is [south, west, north, east], roughly 2km x 2km centered on the
# most characteristic district.
#
# TO ADD/CHANGE CITIES: append {"name": "...", "type": "planned"|"organic",
# "bbox": [south, west, north, east]} to the list below. Choose a ~2km x 2km
# bbox over the district that best represents the city's planning character.
# Avoid areas with very few roads (parks, water) or mixed character (urban edge).
CITIES = [
    # === Planned grid cities ===
    {"name": "Manhattan NYC",       "type": "planned",
     "bbox": [40.748, -73.993, 40.768, -73.968]},
    {"name": "Barcelona Eixample",  "type": "planned",
     "bbox": [41.385,   2.155, 41.405,   2.180]},
    {"name": "Chicago Loop",        "type": "planned",
     "bbox": [41.875, -87.645, 41.895, -87.620]},
    {"name": "Salt Lake City",      "type": "planned",
     "bbox": [40.755, -111.905, 40.775, -111.880]},
    {"name": "Portland OR",         "type": "planned",
     "bbox": [45.510, -122.695, 45.530, -122.670]},
    {"name": "Phoenix AZ",          "type": "planned",
     "bbox": [33.440, -112.090, 33.460, -112.065]},
    {"name": "Buenos Aires",        "type": "planned",
     "bbox": [-34.615, -58.393, -34.595, -58.368]},
    {"name": "Adelaide",            "type": "planned",
     "bbox": [-34.935, 138.590, -34.915, 138.615]},
    {"name": "Savannah GA",         "type": "planned",
     "bbox": [32.070, -81.105, 32.090, -81.080]},
    {"name": "Washington DC",       "type": "planned",
     "bbox": [38.900, -77.042, 38.920, -77.017]},
    # === Organic cities ===
    {"name": "London City",         "type": "organic",
     "bbox": [51.508, -0.098, 51.528, -0.073]},
    {"name": "Tokyo Shinjuku",      "type": "organic",
     "bbox": [35.685, 139.690, 35.705, 139.715]},
    {"name": "Rome Centro",         "type": "organic",
     "bbox": [41.893, 12.465, 41.913, 12.490]},
    {"name": "Istanbul Fatih",      "type": "organic",
     "bbox": [41.005, 28.950, 41.025, 28.975]},
    {"name": "Cairo Old City",      "type": "organic",
     "bbox": [30.045, 31.255, 30.065, 31.280]},
    {"name": "Mumbai Kalbadevi",    "type": "organic",
     "bbox": [18.955, 72.825, 18.975, 72.850]},
    {"name": "Bangkok Old City",    "type": "organic",
     "bbox": [13.740, 100.485, 13.760, 100.510]},
    {"name": "Prague Old Town",     "type": "organic",
     "bbox": [50.083, 14.415, 50.093, 14.430]},
    {"name": "Edinburgh Old Town", "type": "organic",
     "bbox": [55.947, -3.198, 55.957, -3.178]},
    {"name": "Lisbon Alfama",       "type": "organic",
     "bbox": [38.708, -9.142, 38.728, -9.117]},
]


# ─── Helper Functions ─────────────────────────────────────────────────────────

def mean(vals):
    """Arithmetic mean."""
    return sum(vals) / len(vals) if vals else 0.0


def sd(vals):
    """Sample standard deviation."""
    if len(vals) < 2:
        return 0.0
    m = mean(vals)
    return math.sqrt(sum((x - m) ** 2 for x in vals) / (len(vals) - 1))


# ─── Data Download ────────────────────────────────────────────────────────────

def download_city(city):
    """Download road data from Overpass API with caching and retry."""
    os.makedirs(CACHE_DIR, exist_ok=True)
    safe = city["name"].replace(" ", "_").replace(",", "")
    path = os.path.join(CACHE_DIR, f"{safe}.json")

    if os.path.exists(path):
        with open(path, "r") as f:
            return json.loads(f.read()), path

    s, w, n, e = city["bbox"]
    query = (
        f'[out:json][timeout:{REQUEST_TIMEOUT}];'
        f'way["highway"~"^(primary|secondary|tertiary|residential|living_street|unclassified)$"]'
        f'({s},{w},{n},{e});'
        f'(._;>;);out body;'
    )
    body = urllib.parse.urlencode({"data": query}).encode("utf-8")

    for attempt in range(MAX_RETRIES):
        try:
            req = urllib.request.Request(
                OVERPASS_URL, data=body,
                headers={"User-Agent": "Claw4S-StreetEntropy/1.0"}
            )
            with urllib.request.urlopen(req, timeout=REQUEST_TIMEOUT) as resp:
                raw = resp.read().decode("utf-8")
            # Validate JSON before caching
            parsed = json.loads(raw)
            with open(path, "w") as f:
                f.write(raw)
            return parsed, path
        except (urllib.error.URLError, urllib.error.HTTPError, OSError, json.JSONDecodeError) as exc:
            # Exponential backoff with jitter; respect Retry-After if present
            base_wait = min(10 * (2 ** attempt), 120)
            jitter = random.Random(SEED + attempt).uniform(0, base_wait * 0.3)
            retry_after = None
            if hasattr(exc, 'headers') and exc.headers:
                retry_after = exc.headers.get("Retry-After")
            if retry_after and retry_after.isdigit():
                wait = max(int(retry_after), base_wait) + jitter
            else:
                wait = base_wait + jitter
            wait = round(wait, 1)
            if attempt < MAX_RETRIES - 1:
                print(f"    Retry {attempt+1}/{MAX_RETRIES} for {city['name']} "
                      f"(wait {wait}s): {exc}", flush=True)
                time.sleep(wait)
            else:
                raise RuntimeError(
                    f"Failed to download {city['name']} after {MAX_RETRIES} attempts: {exc}"
                )


# ─── Geometry ─────────────────────────────────────────────────────────────────

def bearing_deg(lat1, lon1, lat2, lon2):
    """Initial bearing from point 1 to point 2, in degrees [0, 360)."""
    la1, lo1, la2, lo2 = map(math.radians, [lat1, lon1, lat2, lon2])
    dlon = lo2 - lo1
    x = math.sin(dlon) * math.cos(la2)
    y = math.cos(la1) * math.sin(la2) - math.sin(la1) * math.cos(la2) * math.cos(dlon)
    return (math.degrees(math.atan2(x, y)) + 360) % 360


def extract_orientations(data):
    """Parse Overpass JSON → list of undirected orientations in [0, 180)."""
    nodes = {}
    ways = []
    for el in data.get("elements", []):
        if el["type"] == "node":
            nodes[el["id"]] = (el["lat"], el["lon"])
        elif el["type"] == "way":
            ways.append(el.get("nodes", []))

    orientations = []
    for nds in ways:
        for i in range(len(nds) - 1):
            if nds[i] in nodes and nds[i + 1] in nodes:
                lat1, lon1 = nodes[nds[i]]
                lat2, lon2 = nodes[nds[i + 1]]
                # Skip zero-length segments
                if (lat1, lon1) == (lat2, lon2):
                    continue
                b = bearing_deg(lat1, lon1, lat2, lon2)
                orientations.append(b % 180)  # undirected
    return orientations


# ─── Entropy ──────────────────────────────────────────────────────────────────

def shannon_entropy(orientations, n_bins=N_BINS):
    """Shannon entropy of orientation histogram (bits)."""
    if not orientations:
        return 0.0
    bw = 180.0 / n_bins
    counts = [0] * n_bins
    for o in orientations:
        counts[min(int(o / bw), n_bins - 1)] += 1
    total = sum(counts)
    h = 0.0
    for c in counts:
        if c > 0:
            p = c / total
            h -= p * math.log2(p)
    return h


# ─── Statistical Tests ───────────────────────────────────────────────────────

def exact_enumeration_test(all_values, n_a):
    """
    Exact one-sided test: enumerate ALL C(n, n_a) assignments.
    H0: no difference.  H1: first n_a values (planned) have lower mean.
    Returns (p_value, n_combinations).
    """
    n = len(all_values)
    n_b = n - n_a
    total_sum = sum(all_values)
    # Observed: first n_a are planned
    obs_sum_a = sum(all_values[:n_a])
    obs_diff = obs_sum_a / n_a - (total_sum - obs_sum_a) / n_b

    count_extreme = 0
    count_total = 0
    for combo in itertools.combinations(range(n), n_a):
        sa = sum(all_values[i] for i in combo)
        diff = sa / n_a - (total_sum - sa) / n_b
        if diff <= obs_diff + 1e-12:  # tolerance for float comparison
            count_extreme += 1
        count_total += 1

    return count_extreme / count_total, count_total


def permutation_test(planned, organic, n_perms=N_PERM, seed=SEED):
    """Monte Carlo permutation test (secondary check)."""
    rng = random.Random(seed)
    combined = list(planned) + list(organic)
    obs_diff = mean(planned) - mean(organic)
    na = len(planned)
    count = 0
    for _ in range(n_perms):
        rng.shuffle(combined)
        d = mean(combined[:na]) - mean(combined[na:])
        if d <= obs_diff + 1e-12:
            count += 1
    return count / n_perms


def bootstrap_ci(values, n_boot=N_BOOTSTRAP, seed=SEED):
    """Bootstrap 95% CI for the mean."""
    rng = random.Random(seed)
    means = sorted(
        mean([rng.choice(values) for _ in range(len(values))])
        for _ in range(n_boot)
    )
    return means[int(n_boot * 0.025)], means[int(n_boot * 0.975)]


def bootstrap_diff_ci(g1, g2, n_boot=N_BOOTSTRAP, seed=SEED):
    """Bootstrap 95% CI for difference of means (g1 - g2)."""
    rng = random.Random(seed)
    diffs = []
    for _ in range(n_boot):
        s1 = [rng.choice(g1) for _ in range(len(g1))]
        s2 = [rng.choice(g2) for _ in range(len(g2))]
        diffs.append(mean(s1) - mean(s2))
    diffs.sort()
    return diffs[int(n_boot * 0.025)], diffs[int(n_boot * 0.975)]


def cohens_d(g1, g2):
    """Cohen's d (pooled SD)."""
    m1, m2 = mean(g1), mean(g2)
    n1, n2 = len(g1), len(g2)
    if n1 < 2 or n2 < 2:
        return 0.0
    v1 = sum((x - m1) ** 2 for x in g1) / (n1 - 1)
    v2 = sum((x - m2) ** 2 for x in g2) / (n2 - 1)
    sp = math.sqrt(((n1 - 1) * v1 + (n2 - 1) * v2) / (n1 + n2 - 2))
    return (m1 - m2) / sp if sp > 0 else 0.0


# ─── Report Writer ────────────────────────────────────────────────────────────

def write_report(out):
    gs = out["group_statistics"]
    ht = out["hypothesis_test"]
    es = out["effect_size"]
    ci = out["confidence_intervals"]

    lines = [
        "# Street Orientation Entropy: Planned vs Organic Cities\n",
        "## Summary\n",
        f"Analyzed road networks from {out['metadata']['n_cities']} cities "
        f"({out['metadata']['n_planned']} planned, {out['metadata']['n_organic']} organic) "
        f"using OpenStreetMap data.\n",
        f"**Key finding:** Planned cities have "
        f"{'lower' if gs['difference'] < 0 else 'higher'} "
        f"street orientation entropy (mean {gs['planned_mean']:.4f}, "
        f"SD {gs['planned_sd']:.4f}) compared to organic cities "
        f"(mean {gs['organic_mean']:.4f}, SD {gs['organic_sd']:.4f}), "
        f"difference = {gs['difference']:.4f} bits.\n",
        f"- **Exact combinatorial p-value:** {ht['exact_p_value']:.6f} "
        f"(all {ht['exact_n_combinations']:,} groupings enumerated)",
        f"- **Monte Carlo permutation p-value:** {ht['permutation_p_value']:.4f} "
        f"({ht['permutation_n_shuffles']:,} shuffles)",
        f"- **Cohen's d:** {es['cohens_d']:.4f} ({es['interpretation']} effect)",
        f"- **95% bootstrap CI for difference:** "
        f"[{ci['difference_95ci'][0]:.4f}, {ci['difference_95ci'][1]:.4f}]\n",
        "## Per-City Results\n",
        "| City | Type | Segments | Entropy (bits) | Normalized |",
        "|------|------|----------|---------------|------------|",
    ]
    for c in out["cities"]:
        lines.append(
            f"| {c['name']} | {c['type']} | {c['n_segments']:,} | "
            f"{c['entropy']:.4f} | {c['entropy_normalized']:.4f} |"
        )
    lines += [
        "\n## Sensitivity Analysis\n",
        "### Varying bin count\n",
        "| Bins | Planned Mean | Organic Mean | Exact p | Cohen's d |",
        "|------|-------------|-------------|---------|-----------|",
    ]
    for s in out["sensitivity_bins"].values():
        lines.append(
            f"| {s['n_bins']} | {s['planned_mean']:.4f} | "
            f"{s['organic_mean']:.4f} | {s['exact_p']:.6f} | {s['cohens_d']:.4f} |"
        )
    lines += [
        "\n### Leave-one-out stability\n",
        "| Dropped City | Group | Difference | Direction held? |",
        "|-------------|-------|-----------|----------------|",
    ]
    for loo in out["leave_one_out"]:
        held = "Yes" if loo["diff"] < 0 else "NO"
        lines.append(
            f"| {loo['dropped']} | {loo['group']} | {loo['diff']:.4f} | {held} |"
        )
    lines += [
        "\n## Limitations\n",
        "1. **Bounding-box selection bias:** Results depend on which ~2 km x 2 km "
        "neighborhood is sampled. A different district in the same city could yield "
        "different entropy.",
        "2. **Binary classification oversimplification:** Cities exist on a spectrum "
        "from fully planned to fully organic; our binary grouping loses nuance.",
        "3. **OSM data completeness variation:** OpenStreetMap coverage varies by city "
        "and region. Some cities may have incomplete road networks.",
        "4. **Temporal snapshot:** OSM data reflects the current state. Historically "
        "planned cities may have accumulated organic growth over centuries.",
        "5. **Road type filtering:** We include only primary/secondary/tertiary/"
        "residential/living_street/unclassified roads, excluding alleys, service "
        "roads, and footpaths.",
        "6. **Small sample size:** 10 cities per group limits statistical power. "
        "The exact combinatorial test accounts for this but generalizability is "
        "constrained.",
    ]
    with open("report.md", "w") as f:
        f.write("\n".join(lines) + "\n")


# ─── Verification ─────────────────────────────────────────────────────────────

def run_verification():
    print("=" * 70)
    print("VERIFICATION MODE")
    print("=" * 70)

    passed = failed = 0

    def check(name, cond, detail=""):
        nonlocal passed, failed
        status = "PASS" if cond else "FAIL"
        print(f"  [{status}] {name} {detail}")
        if cond:
            passed += 1
        else:
            failed += 1

    # 1
    check("results.json exists", os.path.exists("results.json"))
    if not os.path.exists("results.json"):
        print(f"\n{passed} passed, {failed} failed")
        sys.exit(1)

    with open("results.json") as f:
        r = json.load(f)

    # 2
    check("20 cities present",
          len(r.get("cities", [])) == 20,
          f"(found {len(r.get('cities', []))})")

    # 3
    np_ = sum(1 for c in r["cities"] if c["type"] == "planned")
    no_ = sum(1 for c in r["cities"] if c["type"] == "organic")
    check("10 planned + 10 organic", np_ == 10 and no_ == 10,
          f"({np_}p, {no_}o)")

    # 4
    max_h = math.log2(r["metadata"]["n_bins"])
    valid_h = all(0 <= c["entropy"] <= max_h + 0.01 for c in r["cities"])
    check("Entropy in valid range [0, log2(n_bins)]", valid_h)

    # 5
    min_seg = min(c["n_segments"] for c in r["cities"])
    check("Every city has >= 30 road segments", min_seg >= 30,
          f"(min={min_seg})")

    # 6
    check("Planned mean < organic mean",
          r["group_statistics"]["planned_mean"] < r["group_statistics"]["organic_mean"],
          f"({r['group_statistics']['planned_mean']:.4f} vs "
          f"{r['group_statistics']['organic_mean']:.4f})")

    # 7
    check("Exact p-value < 0.05",
          r["hypothesis_test"]["exact_p_value"] < 0.05,
          f"(p={r['hypothesis_test']['exact_p_value']:.6f})")

    # 8
    ci = r["confidence_intervals"]["difference_95ci"]
    excludes_zero = ci[1] < 0 or ci[0] > 0
    check("95% CI for difference excludes zero", excludes_zero,
          f"([{ci[0]:.4f}, {ci[1]:.4f}])")

    # 9
    check("Cohen's d magnitude > 0.2",
          abs(r["effect_size"]["cohens_d"]) > 0.2,
          f"(d={r['effect_size']['cohens_d']:.4f})")

    # 10
    sens = r.get("sensitivity_bins", {})
    all_sig = all(v["exact_p"] < 0.10 for v in sens.values())
    check("Sensitivity: all bin sizes p < 0.10", all_sig)

    # 11
    check("report.md exists", os.path.exists("report.md"))

    # 12
    check("SHA256 manifest has 20 entries",
          len(r.get("sha256_manifest", {})) == 20)

    print(f"\n{passed} passed, {failed} failed")
    if failed > 0:
        sys.exit(1)
    print("ALL CHECKS PASSED")


# ─── Main ─────────────────────────────────────────────────────────────────────

def main():
    if "--verify" in sys.argv:
        run_verification()
        return

    print("=" * 70)
    print("STREET ORIENTATION ENTROPY: PLANNED vs ORGANIC CITIES")
    print("=" * 70)

    # ── [1/8] Download ────────────────────────────────────────────────────
    print("\n[1/8] Downloading road networks from OpenStreetMap Overpass API...")
    city_data = {}
    sha_manifest = {}

    for i, city in enumerate(CITIES):
        label = f"  [{i+1}/{len(CITIES)}] {city['name']}"
        print(f"{label}...", end=" ", flush=True)
        data, cpath = download_city(city)
        with open(cpath, "rb") as f:
            sha = hashlib.sha256(f.read()).hexdigest()
        sha_manifest[city["name"]] = sha
        city_data[city["name"]] = data
        nw = sum(1 for e in data.get("elements", []) if e["type"] == "way")
        print(f"OK ({nw} ways, sha256:{sha[:16]})")
        if i < len(CITIES) - 1 and not os.path.exists(
            os.path.join(CACHE_DIR, city["name"].replace(" ", "_").replace(",", "") + ".json")
        ):
            pass  # already cached, no delay needed
        # Always add a small delay to be polite (even if cached, keeps output readable)
        if i < len(CITIES) - 1:
            # Only delay if we actually hit the API (file was just created)
            age = time.time() - os.path.getmtime(cpath)
            if age < 2:
                time.sleep(REQUEST_DELAY)

    with open(MANIFEST_FILE, "w") as f:
        json.dump(sha_manifest, f, indent=2)
    print(f"  SHA256 manifest: {MANIFEST_FILE}")

    # ── [2/8] Extract orientations ────────────────────────────────────────
    print("\n[2/8] Extracting road segment orientations...")
    city_orient = {}
    for city in CITIES:
        oris = extract_orientations(city_data[city["name"]])
        city_orient[city["name"]] = oris
        print(f"  {city['name']:25s} {len(oris):>6,} segments")

    # ── [3/8] Compute entropy ─────────────────────────────────────────────
    print(f"\n[3/8] Computing Shannon entropy ({N_BINS} bins, 0-180 degrees)...")
    h_max = math.log2(N_BINS)
    results_list = []
    for city in CITIES:
        h = shannon_entropy(city_orient[city["name"]], N_BINS)
        hn = h / h_max if h_max > 0 else 0.0
        results_list.append({
            "name": city["name"],
            "type": city["type"],
            "n_segments": len(city_orient[city["name"]]),
            "entropy": round(h, 4),
            "entropy_normalized": round(hn, 4),
            "max_entropy": round(h_max, 4),
        })
        print(f"  {city['name']:25s}  H={h:.4f}  H_norm={hn:.4f}  ({city['type']})")

    planned_h = [r["entropy"] for r in results_list if r["type"] == "planned"]
    organic_h = [r["entropy"] for r in results_list if r["type"] == "organic"]

    print(f"\n  Planned:  mean={mean(planned_h):.4f}  SD={sd(planned_h):.4f}")
    print(f"  Organic:  mean={mean(organic_h):.4f}  SD={sd(organic_h):.4f}")
    print(f"  Diff:     {mean(planned_h) - mean(organic_h):.4f}")

    # ── [4/8] Exact combinatorial test ────────────────────────────────────
    print("\n[4/8] Exact combinatorial test (all C(20,10) = 184,756 groupings)...")
    all_h = planned_h + organic_h
    exact_p, n_combos = exact_enumeration_test(all_h, len(planned_h))
    print(f"  One-sided exact p-value: {exact_p:.6f}")
    print(f"  Combinations evaluated:  {n_combos:,}")

    # ── [5/8] Monte Carlo permutation test ────────────────────────────────
    print(f"\n[5/8] Monte Carlo permutation test ({N_PERM:,} shuffles, seed={SEED})...")
    perm_p = permutation_test(planned_h, organic_h, N_PERM, SEED)
    print(f"  Permutation p-value: {perm_p:.4f}")

    # ── [6/8] Bootstrap CIs & effect size ─────────────────────────────────
    print(f"\n[6/8] Bootstrap confidence intervals ({N_BOOTSTRAP:,} resamples)...")
    ci_p = bootstrap_ci(planned_h, N_BOOTSTRAP, SEED)
    ci_o = bootstrap_ci(organic_h, N_BOOTSTRAP, SEED + 1)
    ci_d = bootstrap_diff_ci(planned_h, organic_h, N_BOOTSTRAP, SEED)
    d = cohens_d(planned_h, organic_h)

    print(f"  Planned mean 95% CI:  [{ci_p[0]:.4f}, {ci_p[1]:.4f}]")
    print(f"  Organic mean 95% CI:  [{ci_o[0]:.4f}, {ci_o[1]:.4f}]")
    print(f"  Difference 95% CI:    [{ci_d[0]:.4f}, {ci_d[1]:.4f}]")
    print(f"  Cohen's d:            {d:.4f}")
    d_interp = "large" if abs(d) >= 0.8 else ("medium" if abs(d) >= 0.5 else "small")
    print(f"  Interpretation:       {d_interp} effect")

    # ── [7/8] Sensitivity analysis ────────────────────────────────────────
    print("\n[7/8] Sensitivity analysis...")

    # 7a: Varying bin count
    print("  7a. Varying number of orientation bins:")
    sens_bins = {}
    for nb in [18, 36, 72]:
        ph = [shannon_entropy(city_orient[c["name"]], nb)
              for c in CITIES if c["type"] == "planned"]
        oh = [shannon_entropy(city_orient[c["name"]], nb)
              for c in CITIES if c["type"] == "organic"]
        ep, _ = exact_enumeration_test(ph + oh, len(ph))
        dd = cohens_d(ph, oh)
        sens_bins[str(nb)] = {
            "n_bins": nb,
            "planned_mean": round(mean(ph), 4),
            "organic_mean": round(mean(oh), 4),
            "exact_p": round(ep, 6),
            "cohens_d": round(dd, 4),
        }
        print(f"      bins={nb:3d}  planned={mean(ph):.4f}  organic={mean(oh):.4f}"
              f"  p={ep:.6f}  d={dd:.4f}")

    # 7b: Leave-one-out
    print("  7b. Leave-one-out stability:")
    loo_results = []
    for i, city in enumerate(CITIES):
        sub = [r for j, r in enumerate(results_list) if j != i]
        pv = [r["entropy"] for r in sub if r["type"] == "planned"]
        ov = [r["entropy"] for r in sub if r["type"] == "organic"]
        if pv and ov:
            diff = mean(pv) - mean(ov)
            loo_results.append({
                "dropped": city["name"],
                "group": city["type"],
                "diff": round(diff, 4),
            })
            held = "yes" if diff < 0 else "NO"
            print(f"      Drop {city['name']:25s} ({city['type']:7s}): "
                  f"diff={diff:+.4f}  direction_held={held}")

    # ── [8/8] Write outputs ───────────────────────────────────────────────
    print("\n[8/8] Writing results...")

    output = {
        "metadata": {
            "analysis": "Street Orientation Entropy: Planned vs Organic Cities",
            "n_cities": len(CITIES),
            "n_planned": len(planned_h),
            "n_organic": len(organic_h),
            "n_bins": N_BINS,
            "seed": SEED,
            "n_bootstrap": N_BOOTSTRAP,
            "n_permutations": N_PERM,
            "data_source": "OpenStreetMap Overpass API (https://overpass-api.de/api/interpreter)",
        },
        "cities": results_list,
        "group_statistics": {
            "planned_mean": round(mean(planned_h), 4),
            "planned_sd": round(sd(planned_h), 4),
            "organic_mean": round(mean(organic_h), 4),
            "organic_sd": round(sd(organic_h), 4),
            "difference": round(mean(planned_h) - mean(organic_h), 4),
        },
        "hypothesis_test": {
            "null_hypothesis": "No difference in mean orientation entropy between planned and organic cities",
            "alternative": "Planned cities have lower orientation entropy (one-sided)",
            "exact_p_value": round(exact_p, 6),
            "exact_n_combinations": n_combos,
            "permutation_p_value": round(perm_p, 4),
            "permutation_n_shuffles": N_PERM,
        },
        "effect_size": {
            "cohens_d": round(d, 4),
            "interpretation": d_interp,
        },
        "confidence_intervals": {
            "planned_mean_95ci": [round(ci_p[0], 4), round(ci_p[1], 4)],
            "organic_mean_95ci": [round(ci_o[0], 4), round(ci_o[1], 4)],
            "difference_95ci": [round(ci_d[0], 4), round(ci_d[1], 4)],
        },
        "sensitivity_bins": sens_bins,
        "leave_one_out": loo_results,
        "sha256_manifest": sha_manifest,
    }

    with open("results.json", "w") as f:
        json.dump(output, f, indent=2)
    print("  results.json written")

    write_report(output)
    print("  report.md written")

    print("\n" + "=" * 70)
    print("ANALYSIS COMPLETE")
    print("=" * 70)


if __name__ == "__main__":
    main()
SCRIPT_EOF
```

**Expected output:** File `script.py` created (exit code 0).

## Step 3: Run analysis

```bash
cd /tmp/claw4s_auto_openstreetmap-city-grid-entropy && python3 script.py
```

**Expected output:**
- Sectioned output `[1/8]` through `[8/8]`
- Per-city download confirmations with SHA256 prefixes
- Entropy values for all 20 cities
- Exact p-value from 184,756 combinatorial evaluations
- Bootstrap confidence intervals
- Sensitivity analysis table
- Ends with `ANALYSIS COMPLETE`
- Creates `results.json`, `report.md`, and `cache/` directory

**Expected runtime:** 3-10 minutes (dominated by API downloads on first run).

## Step 4: Verify results

```bash
cd /tmp/claw4s_auto_openstreetmap-city-grid-entropy && python3 script.py --verify
```

**Expected output:**
- 12 verification checks, all `[PASS]`
- Ends with `ALL CHECKS PASSED`

## Success Criteria

1. All 20 cities download successfully from Overpass API
2. Every city yields ≥ 30 road segments
3. Shannon entropy computed for all cities
4. Exact combinatorial test completes (184,756 evaluations)
5. p-value < 0.05 for planned-vs-organic difference
6. Bootstrap 95% CI for difference excludes zero
7. Sensitivity analysis consistent across 3 bin sizes
8. All 12 verification assertions pass

## Failure Conditions

1. Overpass API unreachable after 3 retries → script exits with error
2. Any city has < 30 road segments → indicates bad bounding box
3. Exact p-value ≥ 0.05 → hypothesis not supported (report honestly)
4. Bootstrap CI includes zero → effect not robust
5. Any verification assertion fails → investigate and fix