Carcion Octo Fest - Raising an Octopus to level 9/8/7/6 : Expected time

import numpy as np
import matplotlib.pyplot as plt

# Simulation parameters
N_TRIALS = 10000       # increase if you want tighter percentile precision
RNG = np.random.default_rng(42)

# Percentiles to report (in feeds)
PCTS = [1, 5, 10, 20, 30, 50, 75, 85, 95, 99]

# Success/Failure/Fled probabilities by level
# fled rate = 0% for Levels 1–5 and 3% for Levels 6–8.
rates = {
    1: {"success": 1.00, "failure": 0.00, "fled": 0.00},
    2: {"success": 0.70, "failure": 0.30, "fled": 0.00},
    3: {"success": 0.55, "failure": 0.45, "fled": 0.00},
    4: {"success": 0.40, "failure": 0.60, "fled": 0.00},
    5: {"success": 0.307, "failure": 0.693, "fled": 0.00},
    6: {"success": 0.205, "failure": 0.765, "fled": 0.03},
    7: {"success": 0.153, "failure": 0.817, "fled": 0.03},
    8: {"success": 0.10,  "failure": 0.87,  "fled": 0.03},
}

def one_trial_feeds_to_goal(goal_level: int, max_feeds: int = None, rng=RNG):
    """
    Simulate one run starting at level 1.
    Count feeds only (no minutes). Startup has zero feeds.
    Each feed:
      - Success: level += 1
      - Failure: level -= 1, but not below Lv.2 (if level is 2, stays at 2)
      - Fled: reset to Lv.2
    Returns:
      - feeds_used (int) if goal reached
      - If max_feeds is provided and reached without hitting goal, returns None and final_level_after_max_feeds
    """
    level = 1
    feeds = 0

    def step_from_level(lvl):
        p_s = rates[lvl]["success"]
        p_f = rates[lvl]["failure"]
        r = rng.random()
        if r < p_s:
            return "success"
        elif r < p_s + p_f:
            return "failure"
        else:
            return "fled"

    if goal_level <= 1:
        return 0, None  # already at goal with zero feeds

    while True:
        if max_feeds is not None and feeds >= max_feeds:
            return None, level

        feeds += 1
        outcome = step_from_level(level)

        if outcome == "success":
            level += 1
        elif outcome == "failure":
            if level > 2:
                level -= 1
            else:
                level = 2
        else:  # fled
            level = 2

        if level >= goal_level:
            return feeds, None

def simulate(goal_level: int, n_trials: int = N_TRIALS, rng=RNG):
    feeds_arr = np.empty(n_trials, dtype=np.int32)
    for i in range(n_trials):
        f, _ = one_trial_feeds_to_goal(goal_level, max_feeds=None, rng=rng)
        feeds_arr[i] = f
    return feeds_arr

def simulate_censored(goal_level: int, feeds_limit: int, n_trials: int = N_TRIALS, rng=RNG):
    """
    Run Monte Carlo with a cap on number of feeds.
    Returns:
      - feeds_to_goal for runs that hit goal within feeds_limit
      - final levels for runs that do NOT hit goal within feeds_limit
    """
    hit_feeds = []
    not_hit_levels = []
    for _ in range(n_trials):
        f, lvl = one_trial_feeds_to_goal(goal_level, max_feeds=feeds_limit, rng=rng)
        if f is not None:
            hit_feeds.append(f)
        else:
            not_hit_levels.append(lvl)
    return np.array(hit_feeds, dtype=np.int32), np.array(not_hit_levels, dtype=np.int32)

def percentiles(feeds: np.ndarray, ps=PCTS):
    return {p: np.percentile(feeds, p) for p in ps}

def report(goal_level: int, feeds_limit: int = 100, n_trials: int = N_TRIALS, rng=RNG, show_plot=True):
    print(f"=== Goal Level {goal_level} ===")
    # Uncensored distribution
    feeds_to_goal = simulate(goal_level, n_trials=n_trials, rng=rng)
    # Percentiles in feeds
    pct = percentiles(feeds_to_goal)
    for p in PCTS:
        print(f"{p}% chance within {pct[p]:.0f} feeds")

    # Censored at feeds_limit for median level not reaching goal
    hit_feeds, not_hit_levels = simulate_censored(goal_level, feeds_limit, n_trials=n_trials, rng=rng)
    if len(not_hit_levels) > 0:
        median_level_not_hit = int(np.median(not_hit_levels))
        print(f"Median level reached if you don't reach level {goal_level} after {feeds_limit} feeds: {median_level_not_hit}")
    else:
        print(f"All runs reached level {goal_level} within {feeds_limit} feeds (no censored cases).")

    # Plot PDF, mark 100 feeds
    if show_plot:
        plt.figure(figsize=(9,5))
        bins = 80
        plt.hist(feeds_to_goal, bins=bins, density=True, alpha=0.6, color='steelblue', edgecolor='black')
        plt.axvline(feeds_limit, color='red', linestyle='--', linewidth=2, label=f"{feeds_limit} feeds")
        plt.title(f"Probability Distribution Function of Feeds to Reach Level {goal_level}")
        plt.xlabel("Feeds")
        plt.ylabel("Probability Density")
        plt.legend()
        plt.tight_layout()
        plt.show()

    return {
        "feeds_to_goal": feeds_to_goal,
        "percentiles": pct,
        "hit_feeds_within_limit": hit_feeds,
        "not_hit_levels_within_limit": not_hit_levels
    }

if __name__ == "__main__":
    # Goal = level 9
    res9 = report(goal_level=9, feeds_limit=100, n_trials=N_TRIALS, rng=RNG, show_plot=True)
    # Goal = level 8
    res8 = report(goal_level=8, feeds_limit=100, n_trials=N_TRIALS, rng=RNG, show_plot=True)
    # Goal = level 7
    res7 = report(goal_level=7, feeds_limit=100, n_trials=N_TRIALS, rng=RNG, show_plot=True)
    # Goal = level 6
    res6 = report(goal_level=6, feeds_limit=100, n_trials=N_TRIALS, rng=RNG, show_plot=True)