benford pig breeding sim

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()

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()

I integrated this version in. yours now powers the pinned version and I've narrowed mine down to just a naive idealized version.