Protonk · May 6, 2026 04:20 · srtdog64 · May 6, 2026 · srtdog64 · May 6, 2026
diff --git a/pig_breeding_naive.py b/pig_breeding_naive.py
 """
 Naive mean-field Minecraft pig-breeding simulation.

 Deterministic two-compartment continuous-time model. Idealized regime only:
 infinite carrots, no pen cap. For carrot-limited, pen-capped, or
 player-scheduled regimes, see the stochastic per-individual file.

 Mechanics approximated:
  * 5-min cooldown per pig: max breeding rate = A/10 events per minute
    (A adult pigs form A/2 pairs; each pair breeds once per cooldown)
  * 20-min maturation: piglets born now become breeders 20 min later
  * Carrots are not modeled in the idealized regime.

 State:
  A(t) = adult pigs  (age >= 20 min)
  B(t) = baby pigs   (age <  20 min)
  N(t) = A + B

 Dynamics:
  dB/dt = births(t) - maturations(t)
  dA/dt = maturations(t)
  births(t)      = A(t) / 10
  maturations(t) = births(t - 20)

 In steady-state exponential growth, r solves 10r = exp(-20r), giving
 r ≈ 0.04263 /min, doubling time ≈ 16.26 min, adult fraction ≈ 0.4263.

 Run as `python pig_breeding_naive.py`.
 """
 from __future__ import annotations

 import math

 import numpy as np


 # ---- Vanilla mechanics ----
 COOLDOWN_MIN = 5.0
 MATURATION_MIN = 20.0

 # ---- Sim defaults ----
 T_TOTAL_MIN = 150.0
 DT_MIN = 0.5
 WARMUP_MIN = 30.0
 A0_DEFAULT = 2.0

 # ---- Numerical tolerances ----
 TIME_TOL = 1e-12
 SLOPE_TOL = 1e-10


 def _exact_step_count(duration: float, dt: float, name: str) -> int:
    """Return duration / dt as an integer, or raise if it is not integral."""
    if duration < 0:
        raise ValueError(f"{name} must be nonnegative")
    steps = int(round(duration / dt))
    if not math.isclose(steps * dt, duration, rel_tol=0.0, abs_tol=TIME_TOL):
        raise ValueError(f"{name}={duration!r} must be an integer multiple of dt={dt!r}")
    return steps


 def solve_steady_state_r() -> float:
    """Solve 10r = exp(-20r) numerically for the idealized growth rate."""
    lo, hi = 0.0, 0.1
    for _ in range(100):
        mid = 0.5 * (lo + hi)
        if 10 * mid - math.exp(-20 * mid) < 0:
            lo = mid
        else:
            hi = mid
    return 0.5 * (lo + hi)


 R_THEORY = solve_steady_state_r()
 T_DOUBLE_THEORY = math.log(2) / R_THEORY
 ADULT_FRAC_THEORY = math.exp(-MATURATION_MIN * R_THEORY)


 def simulate(
    T_total: float = T_TOTAL_MIN,
    dt: float = DT_MIN,
    A0: float = A0_DEFAULT,
    B0: float = 0.0,
 ) -> dict[str, np.ndarray]:
    """
    Run the two-compartment mean-field idealized simulation.

    Returns a dict with keys 'times', 'N', 'A', 'B' (numpy arrays of length
    n_steps + 1; index 0 holds the initial state).
    """
    if dt <= 0:
        raise ValueError("dt must be positive")
    if A0 < 0 or B0 < 0:
        raise ValueError("A0 and B0 must be nonnegative")

    n_steps = _exact_step_count(T_total, dt, "T_total")
    delay_steps = _exact_step_count(MATURATION_MIN, dt, "MATURATION_MIN")

    maturation_due = np.zeros(n_steps + delay_steps + 1)
    A = float(A0)
    B = float(B0)

    times = np.arange(n_steps + 1, dtype=float) * dt
    N_rec = np.zeros(n_steps + 1)
    A_rec = np.zeros(n_steps + 1)
    B_rec = np.zeros(n_steps + 1)

    A_rec[0] = A
    B_rec[0] = B
    N_rec[0] = A + B

    for step in range(n_steps):
        maturations = maturation_due[step]
        if maturations > B and maturations - B < 1e-9:
            maturations = B
        if maturations > B + 1e-9:
            raise RuntimeError("scheduled maturations exceeded baby population")

        A += maturations
        B -= maturations

        rate = A / (2.0 * COOLDOWN_MIN)
        births = rate * dt

        B += births
        maturation_due[step + delay_steps] += births

        A_rec[step + 1] = A
        B_rec[step + 1] = B
        N_rec[step + 1] = A + B

    return {"times": times, "N": N_rec, "A": A_rec, "B": B_rec}


 # ---- Shared analysis helpers (definitions match the stochastic file) ----
 BENFORD = [math.log10(1.0 + 1.0 / d) for d in range(1, 10)]


 def leading_digit(x: float) -> int:
    """Return the base-10 leading digit of x, or 0 for nonpositive/non-finite x."""
    if x <= 0 or not math.isfinite(x):
        return 0
    exponent = math.floor(math.log10(x))
    digit = int(x / (10 ** exponent))
    return max(1, min(9, digit))


 def histogram(values) -> list[float]:
    """Leading-digit histogram from a flat sequence of population values."""
    digits = [leading_digit(float(v)) for v in values]
    n = len(digits)
    if n == 0:
        raise ValueError("no values to histogram")
    return [digits.count(d) / n for d in range(1, 10)]


 def l1_to_benford(hist: list[float]) -> float:
    """L1 distance between a leading-digit histogram and Benford's law."""
    return sum(abs(h - b) for h, b in zip(hist, BENFORD))


 # ---- File-specific helper ----
 def empirical_logfit_doubling_time(
    result: dict[str, np.ndarray], warmup_min: float = WARMUP_MIN
 ) -> float:
    """Fit slope of log N vs t after warmup; return implied doubling time."""
    times = result["times"]
    N = result["N"]
    warmup_idx = int(np.searchsorted(times, warmup_min, side="left"))
    if warmup_idx >= len(N) - 1:
        raise ValueError("warmup_min leaves too few samples")

    fit_t = times[warmup_idx:]
    fit_N = N[warmup_idx:]
    positive = fit_N > 0
    if positive.sum() < 2:
        return float("inf")

    slope, _ = np.polyfit(fit_t[positive], np.log(fit_N[positive]), 1)
    if slope <= SLOPE_TOL:
        return float("inf")
    return float(math.log(2) / slope)


 def main() -> None:
    print(f"Theoretical r:                {R_THEORY:.5f} /min")
    print(f"Theoretical doubling time:    {T_DOUBLE_THEORY:.3f} min")
    print(f"Theoretical adult fraction:   {ADULT_FRAC_THEORY:.4f}")
    print()

    result = simulate()
    times = result["times"]
    N = result["N"]
    warmup_idx = int(np.searchsorted(times, WARMUP_MIN, side="left"))
    values = N[warmup_idx:]

    hist = histogram(values)
    l1 = l1_to_benford(hist)
    final_N = float(N[-1])
    decades = math.log10(max(final_N, 1.0) / 2.0)
    Td_emp = empirical_logfit_doubling_time(result, WARMUP_MIN)

    print("=" * 58)
    print(" NAIVE MEAN-FIELD IDEALIZED")
    print("=" * 58)
    print(f" Run time:                    {T_TOTAL_MIN:.0f} min")
    print(f" Warmup discarded:            {WARMUP_MIN:.0f} min")
    print(f" Final N:                     {final_N:.3e}")
    print(f" Decades of growth:           {decades:.2f}")
    print(f" Empirical doubling time:     {Td_emp:.3f} min")
    print(f" L1 to Benford:               {l1:.4f}")
    print()
    print(" Leading-digit histogram (uniform-in-time samples after warmup)")
    print(f"   {'D':>2} | {'Benford':>8} | {'Empirical':>10}")
    print("   " + "-" * 28)
    for d in range(1, 10):
        print(f"   {d:>2} | {BENFORD[d - 1]:>8.4f} | {hist[d - 1]:>10.4f}")


 if __name__ == "__main__":
    main()
diff --git a/pig_breeding_stochastic_improved.py b/pig_breeding_stochastic_improved.py
 """
 Stochastic/discrete Minecraft pig-breeding simulator.

 This is a per-pig, integer-population simulator intended to complement the
 naive mean-field idealized model. It is not a full Minecraft implementation.
 It models the breeding-relevant mechanics only, with explicit time-step and
 resource semantics.

 Mechanics represented:
  * 5-minute breeding cooldown per adult pig.
  * 20-minute maturation delay from baby to adult.
  * 2 carrots consumed per breeding event.
  * Optional per-minute activation probability for each breedable adult.
  * Random pairing among activated adults.
  * Optional experimental population cap, counting adults plus babies.

 Time semantics:
  * snapshots[0] is the initial state at t = 0, before any breeding.
  * Each step advances one interval [t, t + step_min).
  * Carrot supply is credited at the start of the interval.
  * Breeding, if possible, occurs at the start of the interval.
  * snapshots[k] records the state at the end of the kth interval.
  * The final snapshot is exactly at duration_min; no events are processed at
    the endpoint after the run has already ended.

 Activation semantics:
  activation_prob_per_min is the probability that a breedable adult activates
  during one minute. For a step of length step_min, the interval activation
  probability is 1 - (1 - activation_prob_per_min) ** step_min.

 Regimes compared:
  Idealized-max          unlimited carrots, no cap, all breedable adults activate
  Stochastic-activation  unlimited carrots, no cap, activation probability < 1
  Carrot-limited         steady carrot supply caps the event rate
  Pen-capped             experimental hard population ceiling pins late behavior
  Player-scheduled       sparse carrot bursts; nothing in between

 Run as `python pig_breeding_stochastic_improved.py`.
 """

 from __future__ import annotations

 import math
 import random
 import statistics
 from dataclasses import dataclass
 from typing import Callable


 # ---- Breeding mechanics ----
 COOLDOWN_MIN = 5
 MATURATION_MIN = 20
 CARROTS_PER_BREEDING = 2

 # A stochastic behavior parameter, not a vanilla cooldown/maturation constant.
 DEFAULT_ACTIVATION_PROB_PER_MIN = 0.5


 # ---- Domain types ----

 @dataclass(frozen=True)
 class Pig:
    id: int
    mature_at: int
    cooldown_until: int

    def is_adult(self, time_min: int) -> bool:
        return time_min >= self.mature_at

    def can_breed(self, time_min: int) -> bool:
        return self.is_adult(time_min) and time_min >= self.cooldown_until


 @dataclass(frozen=True)
 class Snapshot:
    """
    State at a snapshot time.

    breeding_events is the number of events during the immediately preceding
    interval. For the initial snapshot at t = 0, it is 0.
    """

    time_min: int
    total_pigs: int
    adult_pigs: int
    baby_pigs: int
    carrots: int
    breeding_events: int
    cumulative_breeding_events: int


 CarrotSupply = Callable[[int, int], int]


 @dataclass(frozen=True)
 class SimulationConfig:
    duration_min: int
    step_min: int
    initial_adults: int
    initial_carrots: int
    carrot_supply: CarrotSupply
    pen_cap: int | None
    random_seed: int | None
    activation_prob_per_min: float = 1.0


 # ---- Carrot supply factories ----

 def constant_supply(rate_per_min: int) -> CarrotSupply:
    """
    Return a supply function that adds rate_per_min * step_min carrots per step.

    The result must be integral because the simulator keeps carrots discrete.
    """
    if rate_per_min < 0:
        raise ValueError("rate_per_min must be nonnegative")

    def supply(time_min: int, step_min: int) -> int:
        del time_min
        return rate_per_min * step_min

    return supply


 def burst_supply(interval_min: int, amount: int) -> CarrotSupply:
    """
    Return a supply function that adds amount carrots at t = interval_min,
    2 * interval_min, 3 * interval_min, ... . No burst occurs at t = 0.
    """
    if interval_min <= 0:
        raise ValueError("interval_min must be positive")
    if amount < 0:
        raise ValueError("amount must be nonnegative")

    def supply(time_min: int, step_min: int) -> int:
        if interval_min % step_min != 0:
            raise ValueError("step_min must divide burst interval_min")
        if time_min <= 0:
            return 0
        if time_min % interval_min != 0:
            return 0
        return amount

    return supply


 # ---- Simulation ----

 class PigBreedingSimulation:
    def __init__(self, config: SimulationConfig) -> None:
        self.validate_config(config)

        self.config = config
        self.rng = random.Random(config.random_seed)
        self.next_pig_id = 1
        self.pigs: list[Pig] = []
        self.carrots = config.initial_carrots
        self.cumulative_breeding_events = 0

        for _ in range(config.initial_adults):
            self.pigs.append(self.create_adult_pig())

    @staticmethod
    def validate_config(config: SimulationConfig) -> None:
        if config.duration_min < 0:
            raise ValueError("duration_min must be nonnegative")
        if config.step_min <= 0:
            raise ValueError("step_min must be positive")
        if config.initial_adults < 0 or config.initial_carrots < 0:
            raise ValueError("initial_adults and initial_carrots must be nonnegative")
        if config.duration_min % config.step_min != 0:
            raise ValueError("duration_min must be a multiple of step_min")
        if COOLDOWN_MIN % config.step_min != 0:
            raise ValueError("step_min must divide COOLDOWN_MIN")
        if MATURATION_MIN % config.step_min != 0:
            raise ValueError("step_min must divide MATURATION_MIN")
        if not 0.0 <= config.activation_prob_per_min <= 1.0:
            raise ValueError("activation_prob_per_min must lie in [0, 1]")
        if config.pen_cap is not None and config.pen_cap < config.initial_adults:
            raise ValueError("pen_cap is below initial_adults")

    def create_adult_pig(self) -> Pig:
        pig = Pig(id=self.next_pig_id, mature_at=0, cooldown_until=0)
        self.next_pig_id += 1
        return pig

    def create_baby_pig(self, birth_time_min: int) -> Pig:
        mature_at = birth_time_min + MATURATION_MIN
        pig = Pig(
            id=self.next_pig_id,
            mature_at=mature_at,
            cooldown_until=mature_at,
        )
        self.next_pig_id += 1
        return pig

    def run(self) -> list[Snapshot]:
        snapshots: list[Snapshot] = [self.create_snapshot(time_min=0, breeding_events=0)]

        for start_min in range(0, self.config.duration_min, self.config.step_min):
            breeding_events = self.step_interval(start_min)
            end_min = start_min + self.config.step_min
            snapshots.append(self.create_snapshot(end_min, breeding_events))

        return snapshots

    def step_interval(self, start_min: int) -> int:
        """Advance the interval [start_min, start_min + step_min)."""
        self.add_carrots(start_min)

        breeders = self.get_activated_breeders(start_min)
        pair_count = len(breeders) // 2
        event_count = self.get_event_limit(pair_count)

        if event_count <= 0:
            return 0

        self.apply_breeding(start_min, breeders, event_count)
        self.cumulative_breeding_events += event_count
        return event_count

    def add_carrots(self, start_min: int) -> None:
        added = self.config.carrot_supply(start_min, self.config.step_min)
        if added < 0:
            raise ValueError("carrot_supply returned a negative amount")
        self.carrots += added

    def interval_activation_probability(self) -> float:
        p = self.config.activation_prob_per_min
        return 1.0 - (1.0 - p) ** self.config.step_min

    def get_activated_breeders(self, time_min: int) -> list[Pig]:
        candidates = [pig for pig in self.pigs if pig.can_breed(time_min)]
        p_activate = self.interval_activation_probability()
        activated = [pig for pig in candidates if self.rng.random() < p_activate]
        self.rng.shuffle(activated)
        return activated

    def get_event_limit(self, pair_count: int) -> int:
        carrot_limit = self.carrots // CARROTS_PER_BREEDING
        capacity_limit = self.get_capacity_limit()
        return min(pair_count, carrot_limit, capacity_limit)

    def get_capacity_limit(self) -> int:
        # This is an experimental cap used to study constrained behavior. It
        # counts adults plus babies and is not meant to claim exact equivalence
        # to Minecraft's passive-mob spawning rules.
        if self.config.pen_cap is None:
            return 10**18
        return max(0, self.config.pen_cap - len(self.pigs))

    def apply_breeding(
        self, time_min: int, breeders: list[Pig], event_count: int
    ) -> None:
        breeder_ids = {pig.id for pig in breeders[: event_count * 2]}

        self.pigs = [
            self.with_cooldown(pig, time_min) if pig.id in breeder_ids else pig
            for pig in self.pigs
        ]

        for _ in range(event_count):
            self.pigs.append(self.create_baby_pig(time_min))

        self.carrots -= event_count * CARROTS_PER_BREEDING
        if self.carrots < 0:
            raise RuntimeError("carrot count became negative")

    @staticmethod
    def with_cooldown(pig: Pig, time_min: int) -> Pig:
        return Pig(
            id=pig.id,
            mature_at=pig.mature_at,
            cooldown_until=time_min + COOLDOWN_MIN,
        )

    def create_snapshot(self, time_min: int, breeding_events: int) -> Snapshot:
        adult_count = sum(1 for pig in self.pigs if pig.is_adult(time_min))
        baby_count = len(self.pigs) - adult_count
        return Snapshot(
            time_min=time_min,
            total_pigs=len(self.pigs),
            adult_pigs=adult_count,
            baby_pigs=baby_count,
            carrots=self.carrots,
            breeding_events=breeding_events,
            cumulative_breeding_events=self.cumulative_breeding_events,
        )


 # ---- Shared analysis helpers ----

 BENFORD = [math.log10(1.0 + 1.0 / d) for d in range(1, 10)]


 def leading_digit(x: float) -> int:
    """Return the base-10 leading digit of x, or 0 for nonpositive/non-finite x."""
    if x <= 0 or not math.isfinite(x):
        return 0
    exponent = math.floor(math.log10(x))
    digit = int(x / (10 ** exponent))
    return max(1, min(9, digit))


 def histogram(values: list[int] | list[float]) -> list[float]:
    """Leading-digit histogram from a flat sequence of population values."""
    digits = [leading_digit(float(v)) for v in values]
    if not digits:
        raise ValueError("no values to histogram")
    return [digits.count(d) / len(digits) for d in range(1, 10)]


 def l1_to_benford(hist: list[float]) -> float:
    """L1 distance between a leading-digit histogram and Benford's law."""
    return sum(abs(h - b) for h, b in zip(hist, BENFORD))


 def values_from_snapshots(
    snapshots: list[Snapshot], warmup_min: int, *, include_warmup_time: bool = True
 ) -> list[int]:
    """Pull total_pigs from snapshots after the warmup cutoff."""
    if include_warmup_time:
        return [s.total_pigs for s in snapshots if s.time_min >= warmup_min]
    return [s.total_pigs for s in snapshots if s.time_min > warmup_min]


 # ---- Multi-seed runner ----

 @dataclass(frozen=True)
 class RegimeResult:
    name: str
    final_n_mean: float
    final_n_std: float
    l1_mean: float
    l1_std: float
    l1_of_average_hist: float
    average_hist: list[float]


 def run_regime(
    name: str,
    config_builder: Callable[[int], SimulationConfig],
    seeds: list[int],
    warmup_min: int,
 ) -> RegimeResult:
    if not seeds:
        raise ValueError("seeds must be nonempty")

    finals: list[int] = []
    l1_values: list[float] = []
    hist_accum = [0.0] * 9

    for seed in seeds:
        config = config_builder(seed)
        snapshots = PigBreedingSimulation(config).run()
        finals.append(snapshots[-1].total_pigs)

        values = values_from_snapshots(snapshots, warmup_min)
        hist = histogram(values)
        l1_values.append(l1_to_benford(hist))
        for digit_index, value in enumerate(hist):
            hist_accum[digit_index] += value

    seed_count = len(seeds)
    average_hist = [h / seed_count for h in hist_accum]

    return RegimeResult(
        name=name,
        final_n_mean=statistics.mean(finals),
        final_n_std=statistics.pstdev(finals) if seed_count > 1 else 0.0,
        l1_mean=statistics.mean(l1_values),
        l1_std=statistics.pstdev(l1_values) if seed_count > 1 else 0.0,
        l1_of_average_hist=l1_to_benford(average_hist),
        average_hist=average_hist,
    )


 # ---- Regime definitions ----

 DURATION_MIN = 150
 STEP_MIN = 1
 INITIAL_ADULTS = 2
 WARMUP_MIN = 30
 SEEDS = list(range(20))


 def build_idealized_max(seed: int) -> SimulationConfig:
    return SimulationConfig(
        duration_min=DURATION_MIN,
        step_min=STEP_MIN,
        initial_adults=INITIAL_ADULTS,
        initial_carrots=10**9,
        carrot_supply=constant_supply(10**6),
        pen_cap=None,
        random_seed=seed,
        activation_prob_per_min=1.0,
    )


 def build_stochastic_activation(seed: int) -> SimulationConfig:
    return SimulationConfig(
        duration_min=DURATION_MIN,
        step_min=STEP_MIN,
        initial_adults=INITIAL_ADULTS,
        initial_carrots=10**9,
        carrot_supply=constant_supply(10**6),
        pen_cap=None,
        random_seed=seed,
        activation_prob_per_min=DEFAULT_ACTIVATION_PROB_PER_MIN,
    )


 def build_carrot_limited(seed: int) -> SimulationConfig:
    return SimulationConfig(
        duration_min=DURATION_MIN,
        step_min=STEP_MIN,
        initial_adults=INITIAL_ADULTS,
        initial_carrots=20,
        carrot_supply=constant_supply(1),
        pen_cap=None,
        random_seed=seed,
        activation_prob_per_min=1.0,
    )


 def build_pen_capped(seed: int) -> SimulationConfig:
    return SimulationConfig(
        duration_min=DURATION_MIN,
        step_min=STEP_MIN,
        initial_adults=INITIAL_ADULTS,
        initial_carrots=10**9,
        carrot_supply=constant_supply(10**6),
        pen_cap=200,
        random_seed=seed,
        activation_prob_per_min=1.0,
    )


 def build_player_scheduled(seed: int) -> SimulationConfig:
    return SimulationConfig(
        duration_min=DURATION_MIN,
        step_min=STEP_MIN,
        initial_adults=INITIAL_ADULTS,
        initial_carrots=0,
        carrot_supply=burst_supply(interval_min=30, amount=30),
        pen_cap=None,
        random_seed=seed,
        activation_prob_per_min=1.0,
    )


 REGIMES: list[tuple[str, Callable[[int], SimulationConfig]]] = [
    ("Idealized-max", build_idealized_max),
    ("Stochastic-activation", build_stochastic_activation),
    ("Carrot-limited", build_carrot_limited),
    ("Pen-capped", build_pen_capped),
    ("Player-scheduled", build_player_scheduled),
 ]


 # ---- Reporting ----

 def print_summary(results: list[RegimeResult]) -> None:
    print("=" * 106)
    print(" REGIME COMPARISON (mean +/- population std over seeds)")
    print("=" * 106)
    print(
        f"{'Regime':<23} | {'Final N (mean +/- std)':>30} | "
        f"{'mean L1 to B':>22} | {'L1(avg hist)':>12}"
    )
    print("-" * 106)
    for r in results:
        final_col = f"{r.final_n_mean:.3e} +/- {r.final_n_std:.3e}"
        l1_col = f"{r.l1_mean:.4f} +/- {r.l1_std:.4f}"
        print(
            f"{r.name:<23} | {final_col:>30} | "
            f"{l1_col:>22} | {r.l1_of_average_hist:>12.4f}"
        )


 def print_histograms(results: list[RegimeResult]) -> None:
    short = {
        "Idealized-max": "Max",
        "Stochastic-activation": "Stoch",
        "Carrot-limited": "Carrot",
        "Pen-capped": "Pen",
        "Player-scheduled": "Player",
    }
    print()
    print("=" * 106)
    print(" AVERAGED LEADING-DIGIT HISTOGRAMS (uniform-in-time snapshots after warmup)")
    print("=" * 106)
    head = f"{'D':>2} | {'Benford':>7}"
    for r in results:
        head += f" | {short[r.name]:>8}"
    print(head)
    print("-" * len(head))
    for d in range(1, 10):
        row = f"{d:>2} | {BENFORD[d - 1]:>7.4f}"
        for r in results:
            row += f" | {r.average_hist[d - 1]:>8.4f}"
        print(row)


 def main() -> None:
    print(f"Seeds per regime: {len(SEEDS)}")
    print(f"Duration: {DURATION_MIN} min, warmup: {WARMUP_MIN} min, step: {STEP_MIN} min")
    print(
        "Snapshot convention: t=0 is initial state; each later snapshot follows "
        "one completed interval."
    )
    print()

    results = [
        run_regime(name, builder, SEEDS, WARMUP_MIN)
        for name, builder in REGIMES
    ]

    print_summary(results)
    print_histograms(results)


 if __name__ == "__main__":
    main()
	"""
	Naive mean-field Minecraft pig-breeding simulation.

	Deterministic two-compartment continuous-time model. Idealized regime only:
	infinite carrots, no pen cap. For carrot-limited, pen-capped, or
	player-scheduled regimes, see the stochastic per-individual file.

	Mechanics approximated:
	* 5-min cooldown per pig: max breeding rate = A/10 events per minute
	(A adult pigs form A/2 pairs; each pair breeds once per cooldown)
	* 20-min maturation: piglets born now become breeders 20 min later
	* Carrots are not modeled in the idealized regime.

	State:
	A(t) = adult pigs (age >= 20 min)
	B(t) = baby pigs (age < 20 min)
	N(t) = A + B

	Dynamics:
	dB/dt = births(t) - maturations(t)
	dA/dt = maturations(t)
	births(t) = A(t) / 10
	maturations(t) = births(t - 20)

	In steady-state exponential growth, r solves 10r = exp(-20r), giving
	r ≈ 0.04263 /min, doubling time ≈ 16.26 min, adult fraction ≈ 0.4263.

	Run as `python pig_breeding_naive.py`.
	"""
	from __future__ import annotations

	import math

	import numpy as np


	# ---- Vanilla mechanics ----
	COOLDOWN_MIN = 5.0
	MATURATION_MIN = 20.0

	# ---- Sim defaults ----
	T_TOTAL_MIN = 150.0
	DT_MIN = 0.5
	WARMUP_MIN = 30.0
	A0_DEFAULT = 2.0

	# ---- Numerical tolerances ----
	TIME_TOL = 1e-12
	SLOPE_TOL = 1e-10


	def _exact_step_count(duration: float, dt: float, name: str) -> int:
	"""Return duration / dt as an integer, or raise if it is not integral."""
	if duration < 0:
	raise ValueError(f"{name} must be nonnegative")
	steps = int(round(duration / dt))
	if not math.isclose(steps * dt, duration, rel_tol=0.0, abs_tol=TIME_TOL):
	raise ValueError(f"{name}={duration!r} must be an integer multiple of dt={dt!r}")
	return steps


	def solve_steady_state_r() -> float:
	"""Solve 10r = exp(-20r) numerically for the idealized growth rate."""
	lo, hi = 0.0, 0.1
	for _ in range(100):
	mid = 0.5 * (lo + hi)
	if 10 * mid - math.exp(-20 * mid) < 0:
	lo = mid
	else:
	hi = mid
	return 0.5 * (lo + hi)


	R_THEORY = solve_steady_state_r()
	T_DOUBLE_THEORY = math.log(2) / R_THEORY
	ADULT_FRAC_THEORY = math.exp(-MATURATION_MIN * R_THEORY)


	def simulate(
	T_total: float = T_TOTAL_MIN,
	dt: float = DT_MIN,
	A0: float = A0_DEFAULT,
	B0: float = 0.0,
	) -> dict[str, np.ndarray]:
	"""
	Run the two-compartment mean-field idealized simulation.

	Returns a dict with keys 'times', 'N', 'A', 'B' (numpy arrays of length
	n_steps + 1; index 0 holds the initial state).
	"""
	if dt <= 0:
	raise ValueError("dt must be positive")
	if A0 < 0 or B0 < 0:
	raise ValueError("A0 and B0 must be nonnegative")

	n_steps = _exact_step_count(T_total, dt, "T_total")
	delay_steps = _exact_step_count(MATURATION_MIN, dt, "MATURATION_MIN")

	maturation_due = np.zeros(n_steps + delay_steps + 1)
	A = float(A0)
	B = float(B0)

	times = np.arange(n_steps + 1, dtype=float) * dt
	N_rec = np.zeros(n_steps + 1)
	A_rec = np.zeros(n_steps + 1)
	B_rec = np.zeros(n_steps + 1)

	A_rec[0] = A
	B_rec[0] = B
	N_rec[0] = A + B

	for step in range(n_steps):
	maturations = maturation_due[step]
	if maturations > B and maturations - B < 1e-9:
	maturations = B
	if maturations > B + 1e-9:
	raise RuntimeError("scheduled maturations exceeded baby population")

	A += maturations
	B -= maturations

	rate = A / (2.0 * COOLDOWN_MIN)
	births = rate * dt

	B += births
	maturation_due[step + delay_steps] += births

	A_rec[step + 1] = A
	B_rec[step + 1] = B
	N_rec[step + 1] = A + B

	return {"times": times, "N": N_rec, "A": A_rec, "B": B_rec}


	# ---- Shared analysis helpers (definitions match the stochastic file) ----
	BENFORD = [math.log10(1.0 + 1.0 / d) for d in range(1, 10)]


	def leading_digit(x: float) -> int:
	"""Return the base-10 leading digit of x, or 0 for nonpositive/non-finite x."""
	if x <= 0 or not math.isfinite(x):
	return 0
	exponent = math.floor(math.log10(x))
	digit = int(x / (10 ** exponent))
	return max(1, min(9, digit))


	def histogram(values) -> list[float]:
	"""Leading-digit histogram from a flat sequence of population values."""
	digits = [leading_digit(float(v)) for v in values]
	n = len(digits)
	if n == 0:
	raise ValueError("no values to histogram")
	return [digits.count(d) / n for d in range(1, 10)]


	def l1_to_benford(hist: list[float]) -> float:
	"""L1 distance between a leading-digit histogram and Benford's law."""
	return sum(abs(h - b) for h, b in zip(hist, BENFORD))


	# ---- File-specific helper ----
	def empirical_logfit_doubling_time(
	result: dict[str, np.ndarray], warmup_min: float = WARMUP_MIN
	) -> float:
	"""Fit slope of log N vs t after warmup; return implied doubling time."""
	times = result["times"]
	N = result["N"]
	warmup_idx = int(np.searchsorted(times, warmup_min, side="left"))
	if warmup_idx >= len(N) - 1:
	raise ValueError("warmup_min leaves too few samples")

	fit_t = times[warmup_idx:]
	fit_N = N[warmup_idx:]
	positive = fit_N > 0
	if positive.sum() < 2:
	return float("inf")

	slope, _ = np.polyfit(fit_t[positive], np.log(fit_N[positive]), 1)
	if slope <= SLOPE_TOL:
	return float("inf")
	return float(math.log(2) / slope)


	def main() -> None:
	print(f"Theoretical r: {R_THEORY:.5f} /min")
	print(f"Theoretical doubling time: {T_DOUBLE_THEORY:.3f} min")
	print(f"Theoretical adult fraction: {ADULT_FRAC_THEORY:.4f}")
	print()

	result = simulate()
	times = result["times"]
	N = result["N"]
	warmup_idx = int(np.searchsorted(times, WARMUP_MIN, side="left"))
	values = N[warmup_idx:]

	hist = histogram(values)
	l1 = l1_to_benford(hist)
	final_N = float(N[-1])
	decades = math.log10(max(final_N, 1.0) / 2.0)
	Td_emp = empirical_logfit_doubling_time(result, WARMUP_MIN)

	print("=" * 58)
	print(" NAIVE MEAN-FIELD IDEALIZED")
	print("=" * 58)
	print(f" Run time: {T_TOTAL_MIN:.0f} min")
	print(f" Warmup discarded: {WARMUP_MIN:.0f} min")
	print(f" Final N: {final_N:.3e}")
	print(f" Decades of growth: {decades:.2f}")
	print(f" Empirical doubling time: {Td_emp:.3f} min")
	print(f" L1 to Benford: {l1:.4f}")
	print()
	print(" Leading-digit histogram (uniform-in-time samples after warmup)")
	print(f" {'D':>2} \| {'Benford':>8} \| {'Empirical':>10}")
	print(" " + "-" * 28)
	for d in range(1, 10):
	print(f" {d:>2} \| {BENFORD[d - 1]:>8.4f} \| {hist[d - 1]:>10.4f}")


	if __name__ == "__main__":
	main()
	"""
	Stochastic/discrete Minecraft pig-breeding simulator.

	This is a per-pig, integer-population simulator intended to complement the
	naive mean-field idealized model. It is not a full Minecraft implementation.
	It models the breeding-relevant mechanics only, with explicit time-step and
	resource semantics.

	Mechanics represented:
	* 5-minute breeding cooldown per adult pig.
	* 20-minute maturation delay from baby to adult.
	* 2 carrots consumed per breeding event.
	* Optional per-minute activation probability for each breedable adult.
	* Random pairing among activated adults.
	* Optional experimental population cap, counting adults plus babies.

	Time semantics:
	* snapshots[0] is the initial state at t = 0, before any breeding.
	* Each step advances one interval [t, t + step_min).
	* Carrot supply is credited at the start of the interval.
	* Breeding, if possible, occurs at the start of the interval.
	* snapshots[k] records the state at the end of the kth interval.
	* The final snapshot is exactly at duration_min; no events are processed at
	the endpoint after the run has already ended.

	Activation semantics:
	activation_prob_per_min is the probability that a breedable adult activates
	during one minute. For a step of length step_min, the interval activation
	probability is 1 - (1 - activation_prob_per_min) ** step_min.

	Regimes compared:
	Idealized-max unlimited carrots, no cap, all breedable adults activate
	Stochastic-activation unlimited carrots, no cap, activation probability < 1
	Carrot-limited steady carrot supply caps the event rate
	Pen-capped experimental hard population ceiling pins late behavior
	Player-scheduled sparse carrot bursts; nothing in between

	Run as `python pig_breeding_stochastic_improved.py`.
	"""

	from __future__ import annotations

	import math
	import random
	import statistics
	from dataclasses import dataclass
	from typing import Callable


	# ---- Breeding mechanics ----
	COOLDOWN_MIN = 5
	MATURATION_MIN = 20
	CARROTS_PER_BREEDING = 2

	# A stochastic behavior parameter, not a vanilla cooldown/maturation constant.
	DEFAULT_ACTIVATION_PROB_PER_MIN = 0.5


	# ---- Domain types ----

	@dataclass(frozen=True)
	class Pig:
	id: int
	mature_at: int
	cooldown_until: int

	def is_adult(self, time_min: int) -> bool:
	return time_min >= self.mature_at

	def can_breed(self, time_min: int) -> bool:
	return self.is_adult(time_min) and time_min >= self.cooldown_until


	@dataclass(frozen=True)
	class Snapshot:
	"""
	State at a snapshot time.

	breeding_events is the number of events during the immediately preceding
	interval. For the initial snapshot at t = 0, it is 0.
	"""

	time_min: int
	total_pigs: int
	adult_pigs: int
	baby_pigs: int
	carrots: int
	breeding_events: int
	cumulative_breeding_events: int


	CarrotSupply = Callable[[int, int], int]


	@dataclass(frozen=True)
	class SimulationConfig:
	duration_min: int
	step_min: int
	initial_adults: int
	initial_carrots: int
	carrot_supply: CarrotSupply
	pen_cap: int \| None
	random_seed: int \| None
	activation_prob_per_min: float = 1.0


	# ---- Carrot supply factories ----

	def constant_supply(rate_per_min: int) -> CarrotSupply:
	"""
	Return a supply function that adds rate_per_min * step_min carrots per step.

	The result must be integral because the simulator keeps carrots discrete.
	"""
	if rate_per_min < 0:
	raise ValueError("rate_per_min must be nonnegative")

	def supply(time_min: int, step_min: int) -> int:
	del time_min
	return rate_per_min * step_min

	return supply


	def burst_supply(interval_min: int, amount: int) -> CarrotSupply:
	"""
	Return a supply function that adds amount carrots at t = interval_min,
	2 * interval_min, 3 * interval_min, ... . No burst occurs at t = 0.
	"""
	if interval_min <= 0:
	raise ValueError("interval_min must be positive")
	if amount < 0:
	raise ValueError("amount must be nonnegative")

	def supply(time_min: int, step_min: int) -> int:
	if interval_min % step_min != 0:
	raise ValueError("step_min must divide burst interval_min")
	if time_min <= 0:
	return 0
	if time_min % interval_min != 0:
	return 0
	return amount

	return supply


	# ---- Simulation ----

	class PigBreedingSimulation:
	def __init__(self, config: SimulationConfig) -> None:
	self.validate_config(config)

	self.config = config
	self.rng = random.Random(config.random_seed)
	self.next_pig_id = 1
	self.pigs: list[Pig] = []
	self.carrots = config.initial_carrots
	self.cumulative_breeding_events = 0

	for _ in range(config.initial_adults):
	self.pigs.append(self.create_adult_pig())

	@staticmethod
	def validate_config(config: SimulationConfig) -> None:
	if config.duration_min < 0:
	raise ValueError("duration_min must be nonnegative")
	if config.step_min <= 0:
	raise ValueError("step_min must be positive")
	if config.initial_adults < 0 or config.initial_carrots < 0:
	raise ValueError("initial_adults and initial_carrots must be nonnegative")
	if config.duration_min % config.step_min != 0:
	raise ValueError("duration_min must be a multiple of step_min")
	if COOLDOWN_MIN % config.step_min != 0:
	raise ValueError("step_min must divide COOLDOWN_MIN")
	if MATURATION_MIN % config.step_min != 0:
	raise ValueError("step_min must divide MATURATION_MIN")
	if not 0.0 <= config.activation_prob_per_min <= 1.0:
	raise ValueError("activation_prob_per_min must lie in [0, 1]")
	if config.pen_cap is not None and config.pen_cap < config.initial_adults:
	raise ValueError("pen_cap is below initial_adults")

	def create_adult_pig(self) -> Pig:
	pig = Pig(id=self.next_pig_id, mature_at=0, cooldown_until=0)
	self.next_pig_id += 1
	return pig

	def create_baby_pig(self, birth_time_min: int) -> Pig:
	mature_at = birth_time_min + MATURATION_MIN
	pig = Pig(
	id=self.next_pig_id,
	mature_at=mature_at,
	cooldown_until=mature_at,
	)
	self.next_pig_id += 1
	return pig

	def run(self) -> list[Snapshot]:
	snapshots: list[Snapshot] = [self.create_snapshot(time_min=0, breeding_events=0)]

	for start_min in range(0, self.config.duration_min, self.config.step_min):
	breeding_events = self.step_interval(start_min)
	end_min = start_min + self.config.step_min
	snapshots.append(self.create_snapshot(end_min, breeding_events))

	return snapshots

	def step_interval(self, start_min: int) -> int:
	"""Advance the interval [start_min, start_min + step_min)."""
	self.add_carrots(start_min)

	breeders = self.get_activated_breeders(start_min)
	pair_count = len(breeders) // 2
	event_count = self.get_event_limit(pair_count)

	if event_count <= 0:
	return 0

	self.apply_breeding(start_min, breeders, event_count)
	self.cumulative_breeding_events += event_count
	return event_count

	def add_carrots(self, start_min: int) -> None:
	added = self.config.carrot_supply(start_min, self.config.step_min)
	if added < 0:
	raise ValueError("carrot_supply returned a negative amount")
	self.carrots += added

	def interval_activation_probability(self) -> float:
	p = self.config.activation_prob_per_min
	return 1.0 - (1.0 - p) ** self.config.step_min

	def get_activated_breeders(self, time_min: int) -> list[Pig]:
	candidates = [pig for pig in self.pigs if pig.can_breed(time_min)]
	p_activate = self.interval_activation_probability()
	activated = [pig for pig in candidates if self.rng.random() < p_activate]
	self.rng.shuffle(activated)
	return activated

	def get_event_limit(self, pair_count: int) -> int:
	carrot_limit = self.carrots // CARROTS_PER_BREEDING
	capacity_limit = self.get_capacity_limit()
	return min(pair_count, carrot_limit, capacity_limit)

	def get_capacity_limit(self) -> int:
	# This is an experimental cap used to study constrained behavior. It
	# counts adults plus babies and is not meant to claim exact equivalence
	# to Minecraft's passive-mob spawning rules.
	if self.config.pen_cap is None:
	return 10**18
	return max(0, self.config.pen_cap - len(self.pigs))

	def apply_breeding(
	self, time_min: int, breeders: list[Pig], event_count: int
	) -> None:
	breeder_ids = {pig.id for pig in breeders[: event_count * 2]}

	self.pigs = [
	self.with_cooldown(pig, time_min) if pig.id in breeder_ids else pig
	for pig in self.pigs
	]

	for _ in range(event_count):
	self.pigs.append(self.create_baby_pig(time_min))

	self.carrots -= event_count * CARROTS_PER_BREEDING
	if self.carrots < 0:
	raise RuntimeError("carrot count became negative")

	@staticmethod
	def with_cooldown(pig: Pig, time_min: int) -> Pig:
	return Pig(
	id=pig.id,
	mature_at=pig.mature_at,
	cooldown_until=time_min + COOLDOWN_MIN,
	)

	def create_snapshot(self, time_min: int, breeding_events: int) -> Snapshot:
	adult_count = sum(1 for pig in self.pigs if pig.is_adult(time_min))
	baby_count = len(self.pigs) - adult_count
	return Snapshot(
	time_min=time_min,
	total_pigs=len(self.pigs),
	adult_pigs=adult_count,
	baby_pigs=baby_count,
	carrots=self.carrots,
	breeding_events=breeding_events,
	cumulative_breeding_events=self.cumulative_breeding_events,
	)


	# ---- Shared analysis helpers ----

	BENFORD = [math.log10(1.0 + 1.0 / d) for d in range(1, 10)]


	def leading_digit(x: float) -> int:
	"""Return the base-10 leading digit of x, or 0 for nonpositive/non-finite x."""
	if x <= 0 or not math.isfinite(x):
	return 0
	exponent = math.floor(math.log10(x))
	digit = int(x / (10 ** exponent))
	return max(1, min(9, digit))


	def histogram(values: list[int] \| list[float]) -> list[float]:
	"""Leading-digit histogram from a flat sequence of population values."""
	digits = [leading_digit(float(v)) for v in values]
	if not digits:
	raise ValueError("no values to histogram")
	return [digits.count(d) / len(digits) for d in range(1, 10)]


	def l1_to_benford(hist: list[float]) -> float:
	"""L1 distance between a leading-digit histogram and Benford's law."""
	return sum(abs(h - b) for h, b in zip(hist, BENFORD))


	def values_from_snapshots(
	snapshots: list[Snapshot], warmup_min: int, *, include_warmup_time: bool = True
	) -> list[int]:
	"""Pull total_pigs from snapshots after the warmup cutoff."""
	if include_warmup_time:
	return [s.total_pigs for s in snapshots if s.time_min >= warmup_min]
	return [s.total_pigs for s in snapshots if s.time_min > warmup_min]


	# ---- Multi-seed runner ----

	@dataclass(frozen=True)
	class RegimeResult:
	name: str
	final_n_mean: float
	final_n_std: float
	l1_mean: float
	l1_std: float
	l1_of_average_hist: float
	average_hist: list[float]


	def run_regime(
	name: str,
	config_builder: Callable[[int], SimulationConfig],
	seeds: list[int],
	warmup_min: int,
	) -> RegimeResult:
	if not seeds:
	raise ValueError("seeds must be nonempty")

	finals: list[int] = []
	l1_values: list[float] = []
	hist_accum = [0.0] * 9

	for seed in seeds:
	config = config_builder(seed)
	snapshots = PigBreedingSimulation(config).run()
	finals.append(snapshots[-1].total_pigs)

	values = values_from_snapshots(snapshots, warmup_min)
	hist = histogram(values)
	l1_values.append(l1_to_benford(hist))
	for digit_index, value in enumerate(hist):
	hist_accum[digit_index] += value

	seed_count = len(seeds)
	average_hist = [h / seed_count for h in hist_accum]

	return RegimeResult(
	name=name,
	final_n_mean=statistics.mean(finals),
	final_n_std=statistics.pstdev(finals) if seed_count > 1 else 0.0,
	l1_mean=statistics.mean(l1_values),
	l1_std=statistics.pstdev(l1_values) if seed_count > 1 else 0.0,
	l1_of_average_hist=l1_to_benford(average_hist),
	average_hist=average_hist,
	)


	# ---- Regime definitions ----

	DURATION_MIN = 150
	STEP_MIN = 1
	INITIAL_ADULTS = 2
	WARMUP_MIN = 30
	SEEDS = list(range(20))


	def build_idealized_max(seed: int) -> SimulationConfig:
	return SimulationConfig(
	duration_min=DURATION_MIN,
	step_min=STEP_MIN,
	initial_adults=INITIAL_ADULTS,
	initial_carrots=10**9,
	carrot_supply=constant_supply(10**6),
	pen_cap=None,
	random_seed=seed,
	activation_prob_per_min=1.0,
	)


	def build_stochastic_activation(seed: int) -> SimulationConfig:
	return SimulationConfig(
	duration_min=DURATION_MIN,
	step_min=STEP_MIN,
	initial_adults=INITIAL_ADULTS,
	initial_carrots=10**9,
	carrot_supply=constant_supply(10**6),
	pen_cap=None,
	random_seed=seed,
	activation_prob_per_min=DEFAULT_ACTIVATION_PROB_PER_MIN,
	)


	def build_carrot_limited(seed: int) -> SimulationConfig:
	return SimulationConfig(
	duration_min=DURATION_MIN,
	step_min=STEP_MIN,
	initial_adults=INITIAL_ADULTS,
	initial_carrots=20,
	carrot_supply=constant_supply(1),
	pen_cap=None,
	random_seed=seed,
	activation_prob_per_min=1.0,
	)


	def build_pen_capped(seed: int) -> SimulationConfig:
	return SimulationConfig(
	duration_min=DURATION_MIN,
	step_min=STEP_MIN,
	initial_adults=INITIAL_ADULTS,
	initial_carrots=10**9,
	carrot_supply=constant_supply(10**6),
	pen_cap=200,
	random_seed=seed,
	activation_prob_per_min=1.0,
	)


	def build_player_scheduled(seed: int) -> SimulationConfig:
	return SimulationConfig(
	duration_min=DURATION_MIN,
	step_min=STEP_MIN,
	initial_adults=INITIAL_ADULTS,
	initial_carrots=0,
	carrot_supply=burst_supply(interval_min=30, amount=30),
	pen_cap=None,
	random_seed=seed,
	activation_prob_per_min=1.0,
	)


	REGIMES: list[tuple[str, Callable[[int], SimulationConfig]]] = [
	("Idealized-max", build_idealized_max),
	("Stochastic-activation", build_stochastic_activation),
	("Carrot-limited", build_carrot_limited),
	("Pen-capped", build_pen_capped),
	("Player-scheduled", build_player_scheduled),
	]


	# ---- Reporting ----

	def print_summary(results: list[RegimeResult]) -> None:
	print("=" * 106)
	print(" REGIME COMPARISON (mean +/- population std over seeds)")
	print("=" * 106)
	print(
	f"{'Regime':<23} \| {'Final N (mean +/- std)':>30} \| "
	f"{'mean L1 to B':>22} \| {'L1(avg hist)':>12}"
	)
	print("-" * 106)
	for r in results:
	final_col = f"{r.final_n_mean:.3e} +/- {r.final_n_std:.3e}"
	l1_col = f"{r.l1_mean:.4f} +/- {r.l1_std:.4f}"
	print(
	f"{r.name:<23} \| {final_col:>30} \| "
	f"{l1_col:>22} \| {r.l1_of_average_hist:>12.4f}"
	)


	def print_histograms(results: list[RegimeResult]) -> None:
	short = {
	"Idealized-max": "Max",
	"Stochastic-activation": "Stoch",
	"Carrot-limited": "Carrot",
	"Pen-capped": "Pen",
	"Player-scheduled": "Player",
	}
	print()
	print("=" * 106)
	print(" AVERAGED LEADING-DIGIT HISTOGRAMS (uniform-in-time snapshots after warmup)")
	print("=" * 106)
	head = f"{'D':>2} \| {'Benford':>7}"
	for r in results:
	head += f" \| {short[r.name]:>8}"
	print(head)
	print("-" * len(head))
	for d in range(1, 10):
	row = f"{d:>2} \| {BENFORD[d - 1]:>7.4f}"
	for r in results:
	row += f" \| {r.average_hist[d - 1]:>8.4f}"
	print(row)


	def main() -> None:
	print(f"Seeds per regime: {len(SEEDS)}")
	print(f"Duration: {DURATION_MIN} min, warmup: {WARMUP_MIN} min, step: {STEP_MIN} min")
	print(
	"Snapshot convention: t=0 is initial state; each later snapshot follows "
	"one completed interval."
	)
	print()

	results = [
	run_regime(name, builder, SEEDS, WARMUP_MIN)
	for name, builder in REGIMES
	]

	print_summary(results)
	print_histograms(results)


	if __name__ == "__main__":
	main()