Hands-on lab — Day 113: Probability: Events, Rules, and Intuition

"""Exercise 1 -- the sample space, events as sets, and probability as counting.

Two fair dice. Enumerate every outcome, define an event as the subset of
outcomes where something is true, and read a probability off the ratio of
two counts -- exactly, with a Fraction, never a float that might be
0.16666666666666663 instead of 1/6.
"""

from fractions import Fraction

import dataset as D
import probability as P

checks_held = []


def check(label: str, condition: bool) -> None:
    checks_held.append((label, condition))
    print(f"  {'ok' if condition else 'FAIL'}: {label}")


print("The sample space of two dice")
print("-" * 60)

space = P.sample_space_two_dice()
print(f"  built with itertools.product: {len(space)} outcomes")
print(f"  first five: {space[:5]}")
check("the space has 36 outcomes", len(space) == 36)
check("it matches the outcomes dataset.py enumerates", set(space) == set(D.TWO_DICE_SPACE))

print()
print("An event is just a subset of the space")
print("-" * 60)

sum_seven = P.event(space, D.is_sum(7))
print(f"  'sum == 7': {sorted(sum_seven)}")
check("6 outcomes sum to 7", len(sum_seven) == 6)

p_sum_seven = P.probability(sum_seven, space)
print(f"  P(sum == 7) = {p_sum_seven} = {float(p_sum_seven):.6f}")
check("P(sum == 7) is exactly Fraction(1, 6)", p_sum_seven == Fraction(1, 6))
check("probability() returns a Fraction, not a float", isinstance(p_sum_seven, Fraction))

double = P.event(space, D.is_double)
p_double = P.probability(double, space)
print(f"  P(both dice match) = {p_double} = {float(p_double):.6f}")
check("P(double) is exactly Fraction(1, 6)", p_double == Fraction(1, 6))

p_whole = P.probability(space, space)
p_empty = P.probability(frozenset(), space)
print(f"  P(whole space) = {p_whole},  P(empty set) = {p_empty}")
check("the three axioms hold here: P(space) = 1", p_whole == 1)
check("and P(empty) = 0", p_empty == 0)
check("and every probability is non-negative", p_sum_seven >= 0 and p_double >= 0)

print()
if all(ok for _, ok in checks_held):
    print(f"01_sample_space_and_events.py: every assertion held. ({len(checks_held)} checks)")
else:
    failed = [label for label, ok in checks_held if not ok]
    raise SystemExit(f"FAILED: {failed}")

"""Exercise 2 -- the addition rule, and exactly how the naive sum lies.

A = "the dice sum to 7" (6 outcomes). B = "the first die shows 6" (6
outcomes). Someone in a hurry adds P(A) + P(B) and gets 1/3. The true answer
is 11/36, and the gap is exactly P(A and B) -- the one outcome, (6, 1), that
belongs to both events and got counted twice.
"""

from fractions import Fraction

import dataset as D
import probability as P

checks_held = []


def check(label: str, condition: bool) -> None:
    checks_held.append((label, condition))
    print(f"  {'ok' if condition else 'FAIL'}: {label}")


space = P.sample_space_two_dice()
a = P.event(space, D.ADDITION_EVENT_A)
b = P.event(space, D.ADDITION_EVENT_B)
a_and_b = a & b

print("A = 'sum is 7', B = 'first die is 6'")
print("-" * 60)
print(f"  A = {sorted(a)}")
print(f"  B = {sorted(b)}")
print(f"  A and B = {sorted(a_and_b)}   <- exactly one outcome, double-counted")

p_a = P.probability(a, space)
p_b = P.probability(b, space)
p_a_and_b = P.probability(a_and_b, space)
print(f"  P(A) = {p_a},  P(B) = {p_b},  P(A and B) = {p_a_and_b}")

naive = P.naive_sum(p_a, p_b)
true_union = P.addition_rule(p_a, p_b, p_a_and_b)
print()
print(f"  naive P(A) + P(B)              = {naive}  = {float(naive):.6f}   WRONG")
print(f"  true  P(A) + P(B) - P(A and B) = {true_union}  = {float(true_union):.6f}   correct")

# Verify the true union directly by counting the union set, independent of
# the formula -- a second, unrelated route to the same number.
union_by_counting = P.probability(a | b, space)
print(f"  counting |A union B| directly  = {union_by_counting}   <- agrees")

check("A has 6 outcomes", len(a) == 6)
check("B has 6 outcomes", len(b) == 6)
check("A and B overlap in exactly one outcome", len(a_and_b) == 1)
check("the naive sum is 1/3", naive == Fraction(1, 3))
check("the true union is 11/36", true_union == Fraction(11, 36))
check("counting the union directly agrees with the formula", union_by_counting == true_union)
check(
    "the naive sum overstates the truth by exactly P(A and B)",
    naive - true_union == p_a_and_b,
)
check("the naive sum is measurably wrong", naive != true_union)

print()
if all(ok for _, ok in checks_held):
    print(f"02_addition_rule.py: every assertion held. ({len(checks_held)} checks)")
else:
    failed = [label for label, ok in checks_held if not ok]
    raise SystemExit(f"FAILED: {failed}")

"""Exercise 3 -- the Chevalier de Méré's two bets, exact and simulated.

In the 1650s de Méré believed these two bets were equally good:

  bet 1: at least one 6 in 4 rolls of one die
  bet 2: at least one double-six in 24 rolls of two dice

The reasoning was seductive: a double six is 1/6 as likely as a six, so
rolling 6x as many times should even things out. It does not. This script
derives both exactly with the complement rule, then confirms both by
simulation -- the exact arithmetic and the measurement have to agree, or one
of them is wrong.
"""

from fractions import Fraction

import numpy as np

import dataset as D
import probability as P
import simulate as S

checks_held = []


def check(label: str, condition: bool) -> None:
    checks_held.append((label, condition))
    print(f"  {'ok' if condition else 'FAIL'}: {label}")


print("Bet 1: at least one 6 in 4 rolls of one die")
print("-" * 60)

# The complement rule collapses "at least one" to one line: 1 minus the
# probability that every single roll fails.
p_no_six_one_roll = P.complement(Fraction(1, 6))
p_no_six_all_four = p_no_six_one_roll**D.DE_MERE_SINGLE_ROLLS
p_bet_one = P.complement(p_no_six_all_four)
print(f"  P(no 6 in one roll)       = {p_no_six_one_roll}")
print(f"  P(no 6 in all 4 rolls)    = {p_no_six_one_roll}^4 = {p_no_six_all_four}")
print(f"  P(at least one 6)         = 1 - {p_no_six_all_four} = {p_bet_one}")
print(f"                            = {float(p_bet_one):.10f}  ~ {round(float(p_bet_one), 4)}")

via_at_least_one = P.at_least_one(Fraction(1, 6), D.DE_MERE_SINGLE_ROLLS)
check("at_least_one() agrees with the step-by-step derivation", via_at_least_one == p_bet_one)
check("bet 1 rounds to 0.5177", round(float(p_bet_one), 4) == 0.5177)
check("bet 1 matches dataset.py's exact value", p_bet_one == D.DE_MERE_SINGLE_EXACT)

print()
print("Bet 2: at least one double-six in 24 rolls of two dice")
print("-" * 60)

p_no_double_one_roll = P.complement(Fraction(1, 36))
p_no_double_all_24 = p_no_double_one_roll**D.DE_MERE_DOUBLE_ROLLS
p_bet_two = P.complement(p_no_double_all_24)
print(f"  P(no double-six in one roll)     = {p_no_double_one_roll}")
print(f"  P(no double-six in all 24 rolls) = {p_no_double_one_roll}^24")
print(f"  P(at least one double-six)       = {float(p_bet_two):.10f}  ~ {round(float(p_bet_two), 4)}")

via_at_least_one_2 = P.at_least_one(Fraction(1, 36), D.DE_MERE_DOUBLE_ROLLS)
check("at_least_one() agrees with the step-by-step derivation", via_at_least_one_2 == p_bet_two)
check("bet 2 rounds to 0.4914", round(float(p_bet_two), 4) == 0.4914)
check("bet 2 matches dataset.py's exact value", p_bet_two == D.DE_MERE_DOUBLE_EXACT)

print()
print("The two bets are NOT equally good")
print("-" * 60)
print(f"  bet 1 (one die,  4 rolls):  {float(p_bet_one):.6f}  -- above 0.5, favours the player")
print(f"  bet 2 (two dice, 24 rolls): {float(p_bet_two):.6f}  -- below 0.5, favours the house")
check("bet 1 is favourable to the player", p_bet_one > Fraction(1, 2))
check("bet 2 is NOT favourable to the player", p_bet_two < Fraction(1, 2))
check("the two bets differ, contrary to the naive '6x the rolls' reasoning", p_bet_one != p_bet_two)

print()
print(f"Simulated at n = {D.DE_MERE_SIM_TRIALS:,} trials per bet, seed 42")
print("-" * 60)

rng = np.random.default_rng(D.REPRODUCIBILITY_SEED_A)
sim_one = S.simulate_at_least_one_six(rng, D.DE_MERE_SIM_TRIALS)
sim_two = S.simulate_at_least_one_double_six(rng, D.DE_MERE_SIM_TRIALS)

print(f"  bet 1: simulated {sim_one:.6f}  vs exact {float(p_bet_one):.6f}"
      f"  (gap {abs(sim_one - float(p_bet_one)):.6f}, tolerance {D.DE_MERE_SINGLE_TOL:.6f})")
print(f"  bet 2: simulated {sim_two:.6f}  vs exact {float(p_bet_two):.6f}"
      f"  (gap {abs(sim_two - float(p_bet_two)):.6f}, tolerance {D.DE_MERE_DOUBLE_TOL:.6f})")

check(
    "bet 1's simulation lands within 3 standard errors of the exact value",
    abs(sim_one - float(p_bet_one)) < D.DE_MERE_SINGLE_TOL,
)
check(
    "bet 2's simulation lands within 3 standard errors of the exact value",
    abs(sim_two - float(p_bet_two)) < D.DE_MERE_DOUBLE_TOL,
)

print()
if all(ok for _, ok in checks_held):
    print(f"03_de_mere.py: every assertion held. ({len(checks_held)} checks)")
else:
    failed = [label for label, ok in checks_held if not ok]
    raise SystemExit(f"FAILED: {failed}")

"""Exercise 4 -- independence versus dependence, the most common conflation
in the subject.

Two events are independent when P(A and B) == P(A) x P(B) -- knowing one
happened tells you nothing about the other. This script finds one genuinely
independent pair and one genuinely dependent pair inside the same sample
space, so the difference is not abstract.
"""

from fractions import Fraction

import dataset as D
import probability as P

checks_held = []


def check(label: str, condition: bool) -> None:
    checks_held.append((label, condition))
    print(f"  {'ok' if condition else 'FAIL'}: {label}")


space = P.sample_space_two_dice()


def summarize(name_a, pred_a, name_b, pred_b):
    a = P.event(space, pred_a)
    b = P.event(space, pred_b)
    p_a = P.probability(a, space)
    p_b = P.probability(b, space)
    p_ab = P.probability(a & b, space)
    independent = P.is_independent(p_a, p_b, p_ab)
    print(f"  {name_a}: P = {p_a}   {name_b}: P = {p_b}")
    print(f"  P(both) = {p_ab}   P(A) x P(B) = {p_a * p_b}   independent? {independent}")
    return p_a, p_b, p_ab, independent


print("A genuinely independent pair")
print("-" * 60)
print("  A = 'sum is 7', B = 'first die is 3'")
print("  Reasoning: whatever the first die shows, exactly one value of the")
print("  second die makes the sum 7 -- so P(sum=7 | first die=v) = 1/6 for")
print("  EVERY v. Conditioning on the first die changes nothing.")
p_a1, p_b1, p_ab1, indep1 = summarize("sum=7", D.is_sum(7), "first=3", D.is_first_die(3))
check("P(sum=7) is 1/6", p_a1 == Fraction(1, 6))
check("P(first=3) is 1/6", p_b1 == Fraction(1, 6))
check("P(A and B) equals P(A) x P(B) exactly", p_ab1 == p_a1 * p_b1)
check("is_independent() reports True", indep1 is True)

print()
print("A genuinely dependent pair")
print("-" * 60)
print("  A = 'sum is 2', B = 'first die is 1'")
print("  Reasoning: sum=2 is ONLY possible when both dice show 1 -- so")
print("  P(sum=2 | first die=1) = 1/6, but P(sum=2 | first die != 1) = 0.")
print("  Conditioning on the first die changes everything.")
p_a2, p_b2, p_ab2, indep2 = summarize("sum=2", D.is_sum(2), "first=1", D.is_first_die(1))
check("P(A and B) does NOT equal P(A) x P(B)", p_ab2 != p_a2 * p_b2)
check("is_independent() reports False", indep2 is False)

print()
print("The conflation this exercise exists to prevent")
print("-" * 60)
print("  | Property             | Independent          | Mutually exclusive        |")
print("  | --------------------- | --------------------- | -------------------------- |")
print("  | P(A and B)             | P(A) x P(B), generally > 0 | exactly 0                    |")
print("  | knowing A happened      | tells you nothing about B | tells you B did NOT happen |")
print("  | can both have P > 0?    | yes                    | yes, but never together    |")
print("  independent events with non-zero probability CAN both happen at once.")
print("  mutually exclusive events with non-zero probability NEVER can --")
print("  which, as exercise 5 shows next, makes them dependent.")

print()
if all(ok for _, ok in checks_held):
    print(f"04_independence_vs_dependence.py: every assertion held. ({len(checks_held)} checks)")
else:
    failed = [label for label, ok in checks_held if not ok]
    raise SystemExit(f"FAILED: {failed}")

"""Exercise 5 -- mutually exclusive events with non-zero probability are
NECESSARILY dependent.

This is the sharpest possible case of dependence, and it is the one most
people get backwards: "mutually exclusive" sounds like it should mean
"unrelated", and it means the opposite. If A and B cannot both happen, then
knowing B happened tells you everything about A -- specifically, that A did
not happen, even though A was possible before you knew anything.
"""

from fractions import Fraction

import dataset as D
import probability as P

checks_held = []


def check(label: str, condition: bool) -> None:
    checks_held.append((label, condition))
    print(f"  {'ok' if condition else 'FAIL'}: {label}")


space = P.sample_space_two_dice()
pred_a, pred_b = D.MUTUALLY_EXCLUSIVE_PAIR
a = P.event(space, pred_a)
b = P.event(space, pred_b)

print("A = 'sum is 2' (only (1,1)), B = 'sum is 12' (only (6,6))")
print("-" * 60)
print(f"  A = {sorted(a)}")
print(f"  B = {sorted(b)}")
print(f"  A and B = {sorted(a & b)}   <- empty. The dice cannot sum to both 2 and 12.")

p_a = P.probability(a, space)
p_b = P.probability(b, space)
p_a_and_b = P.probability(a & b, space)
print(f"  P(A) = {p_a},  P(B) = {p_b},  P(A and B) = {p_a_and_b}")

p_a_given_b = P.conditional(p_a_and_b, p_b)
print()
print(f"  P(A | B) = P(A and B) / P(B) = {p_a_and_b} / {p_b} = {p_a_given_b}")
print(f"  but P(A) on its own is {p_a}")
print(f"  P(A | B) != P(A)  ->  knowing B happened changed what you believe about A")
print(f"  ->  A and B are DEPENDENT, despite being unable to co-occur")

check("A and B are mutually exclusive: their intersection is empty", len(a & b) == 0)
check("P(A) is non-zero", p_a != 0)
check("P(B) is non-zero", p_b != 0)
check("P(A | B) is exactly zero", p_a_given_b == 0)
check("P(A | B) does not equal P(A) -- the events are dependent", p_a_given_b != p_a)

print()
print("Contrast this with the independent pair from exercise 4")
print("-" * 60)
indep_a, indep_b = D.INDEPENDENT_PAIR
ia = P.event(space, indep_a)
ib = P.event(space, indep_b)
p_ia = P.probability(ia, space)
p_ib = P.probability(ib, space)
p_ia_given_ib = P.conditional(P.probability(ia & ib, space), p_ib)
print(f"  independent pair: P(sum=7 | first=3) = {p_ia_given_ib} = P(sum=7) = {p_ia}")
check("for a genuinely independent pair, conditioning changes nothing", p_ia_given_ib == p_ia)
check(
    "mutual exclusivity and independence pull P(A | B) in opposite "
    "directions here: one collapses it to 0, the other leaves it unchanged",
    p_a_given_b == 0 and p_ia_given_ib == p_ia and p_a_given_b != p_ia_given_ib,
)

print()
if all(ok for _, ok in checks_held):
    print(f"05_mutual_exclusivity_implies_dependence.py: every assertion held. ({len(checks_held)} checks)")
else:
    failed = [label for label, ok in checks_held if not ok]
    raise SystemExit(f"FAILED: {failed}")

"""Exercise 6 -- conditioning as restriction: throwing away rows.

P(sum = 8 | first die is even), computed two ways that have to agree: the
formula P(A and B) / P(B), and literally filtering the sample space down to
the rows where B is true and asking what fraction of THOSE rows satisfy A.
That second method is what conditioning IS -- not a formula to apply, but a
smaller table to look at.
"""

from fractions import Fraction

import dataset as D
import probability as P

checks_held = []


def check(label: str, condition: bool) -> None:
    checks_held.append((label, condition))
    print(f"  {'ok' if condition else 'FAIL'}: {label}")


space = P.sample_space_two_dice()
a = P.event(space, D.CONDITIONING_EVENT_A)  # sum == 8
b = P.event(space, D.CONDITIONING_EVENT_B)  # first die even

print("Method 1: the formula, P(A | B) = P(A and B) / P(B)")
print("-" * 60)
p_a_and_b = P.probability(a & b, space)
p_b = P.probability(b, space)
by_formula = P.conditional(p_a_and_b, p_b)
print(f"  A and B = {sorted(a & b)}   P(A and B) = {p_a_and_b}")
print(f"  B = {sorted(b)}   P(B) = {p_b}")
print(f"  P(A | B) = {p_a_and_b} / {p_b} = {by_formula}")

print()
print("Method 2: restriction -- throw away every row where B is false")
print("-" * 60)
restricted_space = b  # the 18 outcomes where the first die is even
restricted_event = P.event(restricted_space, D.CONDITIONING_EVENT_A)
by_filtering = P.probability(restricted_event, restricted_space)
print(f"  restricted sample space (first die even): {len(restricted_space)} outcomes")
print(f"  of those, sum == 8: {sorted(restricted_event)}")
print(f"  P(A | B) by filtering = {len(restricted_event)}/{len(restricted_space)} = {by_filtering}")

check("both methods give exactly 1/6", by_formula == Fraction(1, 6) and by_filtering == Fraction(1, 6))
check("the two methods agree exactly, not approximately", by_formula == by_filtering)

print()
print("Compare against the UNCONDITIONED probability")
print("-" * 60)
p_a_unconditioned = P.probability(a, space)
print(f"  P(sum == 8), no conditioning: {p_a_unconditioned}")
print(f"  P(sum == 8 | first die even): {by_formula}")
check(
    "conditioning on 'first die even' changed the probability of sum=8",
    p_a_unconditioned != by_formula,
)

print()
if all(ok for _, ok in checks_held):
    print(f"06_conditioning_by_restriction.py: every assertion held. ({len(checks_held)} checks)")
else:
    failed = [label for label, ok in checks_held if not ok]
    raise SystemExit(f"FAILED: {failed}")

"""Exercise 7 -- the law of total probability, two urns.

Urn 1 has 3 red and 7 blue balls. Urn 2 has 6 red and 4 blue. A fair coin
picks the urn, then a ball is drawn. What is P(red), overall? Weight each
urn's conditional probability of red by the probability of picking that urn,
and sum. Then check the answer a completely different way: enumerate the
combined 20-outcome experiment directly and count.

This rule is the one Day 115's Bayes' theorem runs backwards.
"""

from fractions import Fraction

import dataset as D
import probability as P

checks_held = []


def check(label: str, condition: bool) -> None:
    checks_held.append((label, condition))
    print(f"  {'ok' if condition else 'FAIL'}: {label}")


print("The setup")
print("-" * 60)
print(f"  urn 1: {D.URN_1_RED} red, {D.URN_1_BLUE} blue")
print(f"  urn 2: {D.URN_2_RED} red, {D.URN_2_BLUE} blue")
print(f"  P(urn 1) = P(urn 2) = 1/2, chosen by a fair coin")

print()
print("Method 1: the law of total probability")
print("-" * 60)
total = P.total_probability(D.URN_PRIOR, D.URN_CONDITIONAL_RED)
for i, (prior, cond) in enumerate(zip(D.URN_PRIOR, D.URN_CONDITIONAL_RED), start=1):
    print(f"  P(urn {i}) x P(red | urn {i}) = {prior} x {cond} = {prior * cond}")
print(f"  P(red) = sum of those = {total} = {float(total)}")

print()
print("Method 2: enumerate the combined 20-outcome experiment directly")
print("-" * 60)
urn1 = ["red"] * D.URN_1_RED + ["blue"] * D.URN_1_BLUE
urn2 = ["red"] * D.URN_2_RED + ["blue"] * D.URN_2_BLUE
combined = [("urn1", ball) for ball in urn1] + [("urn2", ball) for ball in urn2]
print(f"  combined space: {len(combined)} equally likely (urn, ball) pairs")
reds = [outcome for outcome in combined if outcome[1] == "red"]
enumerated = Fraction(len(reds), len(combined))
print(f"  {len(reds)} of them are red: P(red) = {len(reds)}/{len(combined)} = {enumerated}")

check("the weighted total is exactly 9/20", total == Fraction(9, 20))
check("the enumeration agrees exactly with the weighted total", enumerated == total)
check("both equal 0.45", float(total) == 0.45 and float(enumerated) == 0.45)

print()
print("A structural note for Day 115")
print("-" * 60)
print("  P(red) was built by summing P(urn) x P(red | urn) over every urn.")
print("  Bayes' theorem asks the reverse question -- given that the ball WAS")
print("  red, how likely is it that it came from urn 2? -- and its denominator")
print("  is exactly this rule, computed for the event that was observed.")

print()
if all(ok for _, ok in checks_held):
    print(f"07_law_of_total_probability.py: every assertion held. ({len(checks_held)} checks)")
else:
    failed = [label for label, ok in checks_held if not ok]
    raise SystemExit(f"FAILED: {failed}")

"""Exercise 8 -- Monte Carlo error scaling: a hundred times the samples buys
a tenth of the error, not a hundredth.

Estimate P(two dice sum to 7) -- exactly 1/6 -- by simulation, at four sample
sizes four decades apart. The error shrinks like 1/sqrt(n): multiplying n by
10 should shrink the average error by about sqrt(10), roughly 3.16x, not by
10x. Day 117 explains why with the central limit theorem; this script only
measures that it is true.

The assertion is about the SHAPE of the trend, averaged over many seeds --
never a single sampled value, which would be flaky on someone else's machine.
"""

import numpy as np

import dataset as D
import simulate as S

checks_held = []


def check(label: str, condition: bool) -> None:
    checks_held.append((label, condition))
    print(f"  {'ok' if condition else 'FAIL'}: {label}")


target = float(D.MONTE_CARLO_TARGET)
print(f"True probability: P(sum == 7) = {D.MONTE_CARLO_TARGET} = {target:.6f}")
print(f"Averaging over {len(D.MONTE_CARLO_SEEDS)} seeds at each sample size")
print("-" * 60)

mean_errors = []
for n in D.MONTE_CARLO_SAMPLE_SIZES:
    errors = []
    for seed in D.MONTE_CARLO_SEEDS:
        rng = np.random.default_rng(seed)
        estimate = S.simulate_sum_seven(rng, n)
        errors.append(abs(estimate - target))
    mean_error = sum(errors) / len(errors)
    mean_errors.append(mean_error)
    predicted_se = D.standard_error(target, n)
    print(
        f"  n = {n:>7,}   mean |error| = {mean_error:.6f}   "
        f"predicted standard error = {predicted_se:.6f}"
    )

print()
print("Does the error shrink, and does it shrink like 1/sqrt(n)?")
print("-" * 60)
first, last = mean_errors[0], mean_errors[-1]
n_first, n_last = D.MONTE_CARLO_SAMPLE_SIZES[0], D.MONTE_CARLO_SAMPLE_SIZES[-1]
n_ratio = n_last / n_first
error_ratio = first / last
sqrt_prediction = n_ratio**0.5
linear_prediction = n_ratio
print(f"  n grew by a factor of {n_ratio:.0f} (from {n_first:,} to {n_last:,})")
print(f"  the mean error shrank by a factor of {error_ratio:.2f}")
print(f"  a 1/sqrt(n) law predicts a shrink of  {sqrt_prediction:.2f}x")
print(f"  a 1/n law predicts a shrink of        {linear_prediction:.0f}x")
print(f"  observed shrink ({error_ratio:.2f}x) sits far closer to the sqrt(n) "
      "prediction than the 1/n one")

check("the error is monotonically smaller at every larger n", all(
    mean_errors[i + 1] < mean_errors[i] for i in range(len(mean_errors) - 1)
))
check("the error at the largest n is well below the error at the smallest n", last < first / 5.0)
check(
    "the observed shrink lands within a factor of 3 of the sqrt(n) prediction",
    0.33 * sqrt_prediction < error_ratio < 3.0 * sqrt_prediction,
)
check(
    "the observed shrink is nowhere near the 1/n prediction",
    error_ratio < linear_prediction / 10.0,
)

print()
if all(ok for _, ok in checks_held):
    print(f"08_monte_carlo_error_scaling.py: every assertion held. ({len(checks_held)} checks)")
else:
    failed = [label for label, ok in checks_held if not ok]
    raise SystemExit(f"FAILED: {failed}")

"""Exercise 9 -- reproducibility: same seed, byte-identical results.

`numpy.random.default_rng(seed)` builds an independent, stateful Generator.
Two Generators built from the same seed produce identical sequences, so a
whole simulation reproduces exactly. Two Generators built from different
seeds diverge -- but both still land close to the true probability, because
they are both honest estimates of the same thing.
"""

import numpy as np

import dataset as D
import simulate as S

checks_held = []


def check(label: str, condition: bool) -> None:
    checks_held.append((label, condition))
    print(f"  {'ok' if condition else 'FAIL'}: {label}")


target = float(D.MONTE_CARLO_TARGET)
n = D.REPRODUCIBILITY_TRIALS

print(f"Two independent Generators, same seed ({D.REPRODUCIBILITY_SEED_A})")
print("-" * 60)
rng_a1 = np.random.default_rng(D.REPRODUCIBILITY_SEED_A)
rng_a2 = np.random.default_rng(D.REPRODUCIBILITY_SEED_A)
result_a1 = S.simulate_sum_seven(rng_a1, n)
result_a2 = S.simulate_sum_seven(rng_a2, n)
print(f"  run 1: {result_a1}")
print(f"  run 2: {result_a2}")
check("the same seed gives byte-identical results", result_a1 == result_a2)

print()
print(f"A different seed ({D.REPRODUCIBILITY_SEED_B})")
print("-" * 60)
rng_b = np.random.default_rng(D.REPRODUCIBILITY_SEED_B)
result_b = S.simulate_sum_seven(rng_b, n)
print(f"  seed {D.REPRODUCIBILITY_SEED_A}: {result_a1}")
print(f"  seed {D.REPRODUCIBILITY_SEED_B}: {result_b}")
check("a different seed gives a different result", result_a1 != result_b)

tol = 4.0 * D.standard_error(target, n)
print()
print(f"Both still estimate the true probability, {target:.6f}, within tolerance")
print(f"(4 standard errors at n = {n:,}: {tol:.6f})")
print(f"  seed {D.REPRODUCIBILITY_SEED_A}: gap {abs(result_a1 - target):.6f}")
print(f"  seed {D.REPRODUCIBILITY_SEED_B}: gap {abs(result_b - target):.6f}")
check("seed A's estimate is within tolerance of the true value", abs(result_a1 - target) < tol)
check("seed B's estimate is within tolerance of the true value", abs(result_b - target) < tol)

print()
print("Why this matters, and why numpy.random.seed is the wrong tool")
print("-" * 60)
print("  numpy.random.seed(n) mutates ONE GLOBAL state shared by every piece")
print("  of code that calls numpy.random.* -- importing a library that seeds")
print("  it, or calling a function twice, silently changes your results.")
print("  default_rng(seed) hands back an independent object: two Generators")
print("  never interfere with each other, and reproducibility does not")
print("  depend on what else your program happened to do first.")

print()
if all(ok for _, ok in checks_held):
    print(f"09_reproducibility.py: every assertion held. ({len(checks_held)} checks)")
else:
    failed = [label for label, ok in checks_held if not ok]
    raise SystemExit(f"FAILED: {failed}")

"""Make this directory's own modules the ones its tests import.

Both `examples/` and `starter/` contain modules called `probability`,
`simulate`, `dataset` and `answers`, and pytest imports test files by putting
their directory on `sys.path`. Without this file, running `pytest` across both
directories at once would import whichever `probability` was seen first and
then reuse it for the other suite -- so these starter tests would silently
pass against the reference solution instead of skipping. That is a wrong
answer with a green tick on it, which is the worst kind.

So: put this directory first on the import path, and drop any already-imported
module of those names that came from somewhere else.
"""

import sys
from pathlib import Path

HERE = str(Path(__file__).parent.resolve())

if HERE in sys.path:
    sys.path.remove(HERE)
sys.path.insert(0, HERE)

for name in ("probability", "simulate", "dataset", "answers"):
    module = sys.modules.get(name)
    origin = getattr(module, "__file__", "") or ""
    if module is not None and not origin.startswith(HERE):
        del sys.modules[name]

"""The sample space, the events, the urns, and every tolerance this lab
compares against.

Read this file. Nothing here is tuned: every probability below is either
computed by enumeration (via `itertools.product` and `fractions.Fraction`) or
derived algebraically and then checked against the enumeration in the tests.
The simulation tolerance is derived from the standard error of a proportion,
`sqrt(p(1-p)/n)`, with the arithmetic written out beside it -- not chosen by
running a test and loosening the number until it passed.
"""

import itertools
import math
from fractions import Fraction
from typing import Callable

# --------------------------------------------------------------------------
# The sample space: two fair six-sided dice, 36 equally likely outcomes
# --------------------------------------------------------------------------

#: Every ordered pair (first die, second die), each equally likely. Built by
#: enumeration, not written down as a literal list.
TWO_DICE_SPACE: tuple[tuple[int, int], ...] = tuple(
    itertools.product(range(1, 7), range(1, 7))
)

assert len(TWO_DICE_SPACE) == 36

# --------------------------------------------------------------------------
# Events, as predicates over one outcome -- the enumerated set is built by
# filtering the space with these, in probability.py and in the scripts below.
# --------------------------------------------------------------------------


def is_sum(target: int) -> Callable[[tuple[int, int]], bool]:
    """An event predicate: the two dice sum to `target`."""

    def predicate(outcome: tuple[int, int]) -> bool:
        return outcome[0] + outcome[1] == target

    return predicate


def is_first_die(value: int) -> Callable[[tuple[int, int]], bool]:
    """An event predicate: the first die shows `value`."""

    def predicate(outcome: tuple[int, int]) -> bool:
        return outcome[0] == value

    return predicate


def first_die_even(outcome: tuple[int, int]) -> bool:
    """An event predicate: the first die shows an even number."""
    return outcome[0] % 2 == 0


def is_double(outcome: tuple[int, int]) -> bool:
    """An event predicate: both dice show the same value."""
    return outcome[0] == outcome[1]


# --------------------------------------------------------------------------
# The addition rule, on a concrete pair of events
# --------------------------------------------------------------------------

#: A = "the dice sum to 7" (6 outcomes). B = "the first die shows 6" (6
#: outcomes). Their overlap is exactly one outcome, (6, 1), because a first
#: die of 6 needs a second die of 1 to sum to 7. The naive sum P(A) + P(B)
#: double-counts that one outcome, so it overstates P(A union B) by exactly
#: P(A intersect B) = 1/36.
ADDITION_EVENT_A = is_sum(7)
ADDITION_EVENT_B = is_first_die(6)

# --------------------------------------------------------------------------
# The Chevalier de Méré's two bets
# --------------------------------------------------------------------------

#: Bet 1: at least one six in four rolls of one die.
DE_MERE_SINGLE_ROLLS: int = 4
DE_MERE_SINGLE_FACES: int = 6

#: Bet 2: at least one double-six in twenty-four rolls of two dice.
DE_MERE_DOUBLE_ROLLS: int = 24
DE_MERE_DOUBLE_FACES: int = 36  # 6 x 6 possible pairs, one of which is (6, 6)

#: Exact answers, by the complement rule: 1 - P(none of the trials succeed).
#: Fraction keeps every digit; the lesson rounds these to four places.
DE_MERE_SINGLE_EXACT: Fraction = 1 - Fraction(5, 6) ** DE_MERE_SINGLE_ROLLS
DE_MERE_DOUBLE_EXACT: Fraction = 1 - Fraction(35, 36) ** DE_MERE_DOUBLE_ROLLS

assert round(float(DE_MERE_SINGLE_EXACT), 4) == 0.5177
assert round(float(DE_MERE_DOUBLE_EXACT), 4) == 0.4914

# --------------------------------------------------------------------------
# Independence and dependence: two named pairs of events in TWO_DICE_SPACE
# --------------------------------------------------------------------------

#: A genuinely independent pair. "Sum is 7" is independent of "first die is
#: 3" -- for ANY value the first die shows, exactly one value of the second
#: die makes the sum 7, so P(sum = 7 | first die = v) = 1/6 for every v, which
#: is P(sum = 7) itself. That is what independence means.
INDEPENDENT_PAIR = (is_sum(7), is_first_die(3))

#: A genuinely dependent pair. "Sum is 2" requires BOTH dice to show 1, so it
#: is only possible when the first die shows 1 -- P(sum = 2 | first die = 1)
#: = 1/6, but P(sum = 2 | first die != 1) = 0. Those are not equal, so the
#: events are dependent.
DEPENDENT_PAIR = (is_sum(2), is_first_die(1))

# --------------------------------------------------------------------------
# Mutual exclusivity implies dependence
# --------------------------------------------------------------------------

#: "Sum is 2" and "sum is 12" cannot both happen (2 needs (1,1), 12 needs
#: (6,6)), so they are mutually exclusive. Knowing one occurred makes the
#: other impossible -- the sharpest form of dependence there is.
MUTUALLY_EXCLUSIVE_PAIR = (is_sum(2), is_sum(12))

# --------------------------------------------------------------------------
# Conditioning by restriction
# --------------------------------------------------------------------------

#: P(sum = 8 | first die is even), computed two ways in the lab: by the
#: formula P(A and B) / P(B), and by filtering TWO_DICE_SPACE down to the
#: rows where the first die is even and asking what fraction of THOSE rows
#: sum to 8.
CONDITIONING_EVENT_A = is_sum(8)
CONDITIONING_EVENT_B = first_die_even

# --------------------------------------------------------------------------
# The law of total probability: two urns, drawn from with a fair coin
# --------------------------------------------------------------------------

#: Urn 1: 3 red, 7 blue, out of 10 balls.
URN_1_RED: int = 3
URN_1_BLUE: int = 7

#: Urn 2: 6 red, 4 blue, out of 10 balls.
URN_2_RED: int = 6
URN_2_BLUE: int = 4

assert URN_1_RED + URN_1_BLUE == 10
assert URN_2_RED + URN_2_BLUE == 10

#: A fair coin decides which urn to draw from -- P(urn 1) = P(urn 2) = 1/2.
URN_PRIOR: tuple[Fraction, Fraction] = (Fraction(1, 2), Fraction(1, 2))

#: P(red | urn 1) and P(red | urn 2), read straight off each urn's contents.
URN_CONDITIONAL_RED: tuple[Fraction, Fraction] = (
    Fraction(URN_1_RED, URN_1_RED + URN_1_BLUE),
    Fraction(URN_2_RED, URN_2_RED + URN_2_BLUE),
)

# --------------------------------------------------------------------------
# Monte Carlo error scaling
# --------------------------------------------------------------------------

#: The event whose probability the Monte Carlo experiment estimates: two
#: fair dice sum to 7. The exact answer is 1/6, established in exercise 1.
MONTE_CARLO_TARGET: Fraction = Fraction(1, 6)

#: Sample sizes to sweep, four decades apart at the ends.
MONTE_CARLO_SAMPLE_SIZES: tuple[int, ...] = (100, 1_000, 10_000, 100_000)

#: Independent seeds to average over at each sample size, so the assertion is
#: about the SHAPE of the error trend across many runs rather than about one
#: sampled value, which would be flaky on someone else's machine.
MONTE_CARLO_SEEDS: tuple[int, ...] = tuple(range(20))


def standard_error(p: float, n: int) -> float:
    """The standard error of a proportion estimated from n trials.

    This is the quantity that governs how far a Monte Carlo estimate can be
    expected to land from the true probability it estimates: about 68% of
    estimates land within one standard error, and about 99.7% within three.
    """
    return math.sqrt(p * (1.0 - p) / n)


#: Three standard errors, at the sample sizes de Méré's simulations use in
#: exercise 3. This is not a guessed tolerance: it is the width inside which
#: 99.7% of repeated simulations should land, derived from the formula above.
DE_MERE_SIM_TRIALS: int = 200_000
DE_MERE_SINGLE_SE: float = standard_error(float(DE_MERE_SINGLE_EXACT), DE_MERE_SIM_TRIALS)
DE_MERE_DOUBLE_SE: float = standard_error(float(DE_MERE_DOUBLE_EXACT), DE_MERE_SIM_TRIALS)
DE_MERE_SINGLE_TOL: float = 3.0 * DE_MERE_SINGLE_SE
DE_MERE_DOUBLE_TOL: float = 3.0 * DE_MERE_DOUBLE_SE

# --------------------------------------------------------------------------
# Reproducibility
# --------------------------------------------------------------------------

REPRODUCIBILITY_SEED_A: int = 42
REPRODUCIBILITY_SEED_B: int = 43
REPRODUCIBILITY_TRIALS: int = 10_000

"""Exercise 1, 2, 4, 5, 6 and 7: exact probability, computed two ways.

Every function here returns a `fractions.Fraction` wherever the answer is
rational, so an assertion against it is exact rather than "close enough".
"""

import itertools
from fractions import Fraction
from typing import Callable, Iterable


# ---------------------------------------------------------------------------
# Exercise 1: the sample space, and probability as counting
# ---------------------------------------------------------------------------


def sample_space_two_dice() -> tuple[tuple[int, int], ...]:
    """Every ordered pair (first die, second die), 36 outcomes."""
    return tuple(itertools.product(range(1, 7), range(1, 7)))


def event(
    space: Iterable[tuple[int, int]], predicate: Callable[[tuple[int, int]], bool]
) -> frozenset[tuple[int, int]]:
    """The subset of `space` for which `predicate` is true."""
    return frozenset(outcome for outcome in space if predicate(outcome))


def probability(
    outcome_set: Iterable[tuple[int, int]], space: Iterable[tuple[int, int]]
) -> Fraction:
    """P(event), for equally likely outcomes: |event| / |space|, exactly."""
    outcome_set = frozenset(outcome_set)
    space = tuple(space)
    return Fraction(len(outcome_set), len(space))


# ---------------------------------------------------------------------------
# Exercise 2: the addition rule
# ---------------------------------------------------------------------------


def addition_rule(p_a: Fraction, p_b: Fraction, p_a_and_b: Fraction) -> Fraction:
    """P(A or B) = P(A) + P(B) - P(A and B)."""
    return p_a + p_b - p_a_and_b


def naive_sum(p_a: Fraction, p_b: Fraction) -> Fraction:
    """The wrong shortcut: P(A) + P(B), which double-counts the overlap."""
    return p_a + p_b


# ---------------------------------------------------------------------------
# Exercise 3 helper: the complement rule
# ---------------------------------------------------------------------------


def complement(p: Fraction) -> Fraction:
    """P(not A) = 1 - P(A)."""
    return 1 - p


def at_least_one(p_single_success: Fraction, trials: int) -> Fraction:
    """P(at least one success in `trials` independent tries).

    Collapsed to one line with the complement rule: 1 minus the probability
    that every single trial fails.
    """
    p_single_failure = complement(p_single_success)
    return complement(p_single_failure**trials)


# ---------------------------------------------------------------------------
# Exercise 4: independence
# ---------------------------------------------------------------------------


def is_independent(p_a: Fraction, p_b: Fraction, p_a_and_b: Fraction) -> bool:
    """True exactly when P(A and B) == P(A) * P(B)."""
    return p_a_and_b == p_a * p_b


# ---------------------------------------------------------------------------
# Exercise 5 and 6: conditional probability
# ---------------------------------------------------------------------------


def conditional(p_a_and_b: Fraction, p_b: Fraction) -> Fraction:
    """P(A | B) = P(A and B) / P(B)."""
    if p_b == 0:
        raise ValueError("conditional probability given a zero-probability event")
    return p_a_and_b / p_b


# ---------------------------------------------------------------------------
# Exercise 7: the law of total probability
# ---------------------------------------------------------------------------


def total_probability(
    priors: Iterable[Fraction], conditionals: Iterable[Fraction]
) -> Fraction:
    """P(A) = sum over i of P(A | condition_i) * P(condition_i)."""
    total = Fraction(0)
    for prior, cond in zip(priors, conditionals):
        total += prior * cond
    return total

"""Exercise 3, 8 and 9: simulation, and the error that comes with it.

Every simulation here takes an explicit `numpy.random.Generator` -- built by
`numpy.random.default_rng(seed)` -- rather than touching global random state.
Two calls with the same seed return byte-identical arrays; that is
`default_rng`'s whole advantage over the legacy `numpy.random.seed` global,
and exercise 9 asserts it directly.
"""

import numpy as np


# ---------------------------------------------------------------------------
# Exercise 3: de Méré's two bets, simulated
# ---------------------------------------------------------------------------


def simulate_at_least_one_six(rng: np.random.Generator, trials: int) -> float:
    """Simulate de Méré's first bet: at least one 6 in 4 rolls of one die."""
    rolls = rng.integers(1, 7, size=(trials, 4))
    at_least_one = (rolls == 6).any(axis=1)
    return float(np.mean(at_least_one))


def simulate_at_least_one_double_six(rng: np.random.Generator, trials: int) -> float:
    """Simulate de Méré's second bet: at least one double-six in 24 rolls of
    two dice."""
    first = rng.integers(1, 7, size=(trials, 24))
    second = rng.integers(1, 7, size=(trials, 24))
    double_six = (first == 6) & (second == 6)
    at_least_one = double_six.any(axis=1)
    return float(np.mean(at_least_one))


# ---------------------------------------------------------------------------
# Exercise 8: Monte Carlo error scaling
# ---------------------------------------------------------------------------


def simulate_sum_seven(rng: np.random.Generator, trials: int) -> float:
    """Simulate P(two fair dice sum to 7) by rolling `trials` pairs of dice."""
    first = rng.integers(1, 7, size=trials)
    second = rng.integers(1, 7, size=trials)
    return float(np.mean(first + second == 7))

"""The reference suite: real values, real exceptions, real simulations.

Run from the lab directory:

    .venv/bin/pytest examples -q -p no:cacheprovider
"""

import itertools
from fractions import Fraction

import numpy as np
import pytest

import dataset as D
import probability as P
import simulate as S

SPACE = P.sample_space_two_dice()

# ---------------------------------------------------------------------------
# Exercise 1 -- the sample space and probability as counting
# ---------------------------------------------------------------------------


def test_sample_space_has_36_outcomes():
    assert len(SPACE) == 36


def test_sample_space_has_no_duplicates():
    assert len(set(SPACE)) == len(SPACE)


def test_sample_space_matches_direct_itertools_product():
    assert set(SPACE) == set(itertools.product(range(1, 7), range(1, 7)))


def test_event_of_sum_seven_has_six_outcomes():
    ev = P.event(SPACE, D.is_sum(7))
    assert ev == frozenset({(1, 6), (2, 5), (3, 4), (4, 3), (5, 2), (6, 1)})


def test_event_of_sum_two_has_one_outcome():
    ev = P.event(SPACE, D.is_sum(2))
    assert ev == frozenset({(1, 1)})


def test_event_of_sum_twelve_has_one_outcome():
    ev = P.event(SPACE, D.is_sum(12))
    assert ev == frozenset({(6, 6)})


@pytest.mark.parametrize(
    "target,count",
    [(2, 1), (3, 2), (4, 3), (5, 4), (6, 5), (7, 6), (8, 5), (9, 4), (10, 3), (11, 2), (12, 1)],
)
def test_event_counts_for_every_possible_sum(target, count):
    assert len(P.event(SPACE, D.is_sum(target))) == count


def test_probability_of_sum_seven_is_exactly_one_sixth():
    ev = P.event(SPACE, D.is_sum(7))
    assert P.probability(ev, SPACE) == Fraction(1, 6)


def test_probability_returns_a_fraction_type():
    ev = P.event(SPACE, D.is_sum(7))
    assert isinstance(P.probability(ev, SPACE), Fraction)


def test_probability_of_the_whole_space_is_one():
    assert P.probability(SPACE, SPACE) == 1


def test_probability_of_the_empty_event_is_zero():
    assert P.probability(frozenset(), SPACE) == 0


def test_probability_of_a_double_is_one_sixth():
    ev = P.event(SPACE, D.is_double)
    assert P.probability(ev, SPACE) == Fraction(1, 6)


# ---------------------------------------------------------------------------
# Exercise 2 -- the addition rule
# ---------------------------------------------------------------------------


def test_addition_rule_matches_enumerated_union():
    a = P.event(SPACE, D.ADDITION_EVENT_A)
    b = P.event(SPACE, D.ADDITION_EVENT_B)
    p_a, p_b = P.probability(a, SPACE), P.probability(b, SPACE)
    p_ab = P.probability(a & b, SPACE)
    formula = P.addition_rule(p_a, p_b, p_ab)
    enumerated = P.probability(a | b, SPACE)
    assert formula == enumerated == Fraction(11, 36)


def test_naive_sum_overstates_by_exactly_the_intersection():
    a = P.event(SPACE, D.ADDITION_EVENT_A)
    b = P.event(SPACE, D.ADDITION_EVENT_B)
    p_a, p_b = P.probability(a, SPACE), P.probability(b, SPACE)
    p_ab = P.probability(a & b, SPACE)
    naive = P.naive_sum(p_a, p_b)
    true_union = P.addition_rule(p_a, p_b, p_ab)
    assert naive - true_union == p_ab


def test_addition_rule_reduces_to_naive_sum_for_disjoint_events():
    a = P.event(SPACE, D.is_sum(2))
    b = P.event(SPACE, D.is_sum(12))
    assert len(a & b) == 0
    p_a, p_b = P.probability(a, SPACE), P.probability(b, SPACE)
    p_ab = P.probability(a & b, SPACE)
    assert P.addition_rule(p_a, p_b, p_ab) == P.naive_sum(p_a, p_b)


@pytest.mark.parametrize(
    "target_a,target_b",
    [(4, 6), (5, 9), (7, 11)],
)
def test_addition_rule_on_further_pairs_matches_enumeration(target_a, target_b):
    a = P.event(SPACE, D.is_sum(target_a))
    b = P.event(SPACE, D.is_sum(target_b))
    p_a, p_b = P.probability(a, SPACE), P.probability(b, SPACE)
    p_ab = P.probability(a & b, SPACE)
    assert P.addition_rule(p_a, p_b, p_ab) == P.probability(a | b, SPACE)


# ---------------------------------------------------------------------------
# Exercise 3 -- the complement rule and de Méré, exact
# ---------------------------------------------------------------------------


def test_complement_of_a_sixth_is_five_sixths():
    assert P.complement(Fraction(1, 6)) == Fraction(5, 6)


def test_complement_is_its_own_inverse():
    p = Fraction(7, 20)
    assert P.complement(P.complement(p)) == p


def test_complement_of_zero_is_one():
    assert P.complement(Fraction(0)) == 1


def test_complement_of_one_is_zero():
    assert P.complement(Fraction(1)) == 0


def test_at_least_one_of_zero_trials_is_zero():
    assert P.at_least_one(Fraction(1, 6), 0) == 0


def test_at_least_one_of_a_certain_event_is_certain():
    assert P.at_least_one(Fraction(1, 1), 5) == 1


def test_de_mere_bet_one_matches_dataset_exact_value():
    assert P.at_least_one(Fraction(1, 6), D.DE_MERE_SINGLE_ROLLS) == D.DE_MERE_SINGLE_EXACT


def test_de_mere_bet_two_matches_dataset_exact_value():
    assert P.at_least_one(Fraction(1, 36), D.DE_MERE_DOUBLE_ROLLS) == D.DE_MERE_DOUBLE_EXACT


def test_de_mere_bet_one_rounds_to_the_historical_figure():
    assert round(float(D.DE_MERE_SINGLE_EXACT), 4) == 0.5177


def test_de_mere_bet_two_rounds_to_the_historical_figure():
    assert round(float(D.DE_MERE_DOUBLE_EXACT), 4) == 0.4914


def test_de_mere_bet_one_favours_the_player():
    assert D.DE_MERE_SINGLE_EXACT > Fraction(1, 2)


def test_de_mere_bet_two_does_not_favour_the_player():
    assert D.DE_MERE_DOUBLE_EXACT < Fraction(1, 2)


def test_de_mere_bets_are_not_equal_despite_the_six_to_one_scaling():
    # 24 = 6 x 4, matching the 6x smaller probability of a double six --
    # and the two bets are still not equal. That gap IS the lesson.
    assert D.DE_MERE_SINGLE_EXACT != D.DE_MERE_DOUBLE_EXACT


def test_de_mere_simulation_bet_one_within_three_standard_errors():
    rng = np.random.default_rng(D.REPRODUCIBILITY_SEED_A)
    got = S.simulate_at_least_one_six(rng, D.DE_MERE_SIM_TRIALS)
    assert abs(got - float(D.DE_MERE_SINGLE_EXACT)) < D.DE_MERE_SINGLE_TOL


def test_de_mere_simulation_bet_two_within_three_standard_errors():
    rng = np.random.default_rng(D.REPRODUCIBILITY_SEED_A)
    got = S.simulate_at_least_one_double_six(rng, D.DE_MERE_SIM_TRIALS)
    assert abs(got - float(D.DE_MERE_DOUBLE_EXACT)) < D.DE_MERE_DOUBLE_TOL


@pytest.mark.parametrize("seed", [1, 2, 3, 4, 5])
def test_de_mere_simulation_bet_one_is_stable_across_seeds(seed):
    rng = np.random.default_rng(seed)
    got = S.simulate_at_least_one_six(rng, D.DE_MERE_SIM_TRIALS)
    assert abs(got - float(D.DE_MERE_SINGLE_EXACT)) < D.DE_MERE_SINGLE_TOL


@pytest.mark.parametrize("seed", [1, 2, 3, 4, 5])
def test_de_mere_simulation_bet_two_is_stable_across_seeds(seed):
    rng = np.random.default_rng(seed)
    got = S.simulate_at_least_one_double_six(rng, D.DE_MERE_SIM_TRIALS)
    assert abs(got - float(D.DE_MERE_DOUBLE_EXACT)) < D.DE_MERE_DOUBLE_TOL


def test_simulate_at_least_one_six_returns_a_plain_float():
    rng = np.random.default_rng(0)
    got = S.simulate_at_least_one_six(rng, 1000)
    assert isinstance(got, float)
    assert 0.0 <= got <= 1.0


def test_simulate_at_least_one_double_six_returns_a_plain_float():
    rng = np.random.default_rng(0)
    got = S.simulate_at_least_one_double_six(rng, 1000)
    assert isinstance(got, float)
    assert 0.0 <= got <= 1.0


# ---------------------------------------------------------------------------
# Exercise 4 -- independence
# ---------------------------------------------------------------------------


def _pair_stats(pred_a, pred_b):
    a = P.event(SPACE, pred_a)
    b = P.event(SPACE, pred_b)
    p_a, p_b = P.probability(a, SPACE), P.probability(b, SPACE)
    p_ab = P.probability(a & b, SPACE)
    return p_a, p_b, p_ab


def test_the_independent_pair_satisfies_the_product_rule_exactly():
    p_a, p_b, p_ab = _pair_stats(*D.INDEPENDENT_PAIR)
    assert p_ab == p_a * p_b
    assert P.is_independent(p_a, p_b, p_ab) is True


def test_the_dependent_pair_fails_the_product_rule():
    p_a, p_b, p_ab = _pair_stats(*D.DEPENDENT_PAIR)
    assert p_ab != p_a * p_b
    assert P.is_independent(p_a, p_b, p_ab) is False


def test_is_independent_returns_a_python_bool():
    p_a, p_b, p_ab = _pair_stats(*D.INDEPENDENT_PAIR)
    assert isinstance(P.is_independent(p_a, p_b, p_ab), bool)


@pytest.mark.parametrize("first_die_value", [1, 2, 4, 5, 6])
def test_sum_seven_is_independent_of_every_value_of_the_first_die(first_die_value):
    # The independence of "sum = 7" from the first die is not special to the
    # value 3 -- it holds for every value, because exactly one second-die
    # value completes the sum to 7 regardless.
    p_a, p_b, p_ab = _pair_stats(D.is_sum(7), D.is_first_die(first_die_value))
    assert p_ab == p_a * p_b


# ---------------------------------------------------------------------------
# Exercise 5 -- mutual exclusivity implies dependence
# ---------------------------------------------------------------------------


def test_mutually_exclusive_pair_has_empty_intersection():
    pred_a, pred_b = D.MUTUALLY_EXCLUSIVE_PAIR
    a, b = P.event(SPACE, pred_a), P.event(SPACE, pred_b)
    assert (a & b) == frozenset()


def test_mutually_exclusive_pair_has_zero_conditional_probability():
    pred_a, pred_b = D.MUTUALLY_EXCLUSIVE_PAIR
    a, b = P.event(SPACE, pred_a), P.event(SPACE, pred_b)
    p_ab = P.probability(a & b, SPACE)
    p_b = P.probability(b, SPACE)
    assert P.conditional(p_ab, p_b) == 0


def test_mutually_exclusive_pair_is_therefore_dependent():
    pred_a, pred_b = D.MUTUALLY_EXCLUSIVE_PAIR
    a, b = P.event(SPACE, pred_a), P.event(SPACE, pred_b)
    p_a = P.probability(a, SPACE)
    p_ab = P.probability(a & b, SPACE)
    p_b = P.probability(b, SPACE)
    p_a_given_b = P.conditional(p_ab, p_b)
    assert p_a != 0
    assert p_a_given_b != p_a


def test_conditional_refuses_to_divide_by_a_zero_probability_event():
    with pytest.raises(ValueError):
        P.conditional(Fraction(0), Fraction(0))


# ---------------------------------------------------------------------------
# Exercise 6 -- conditioning by restriction
# ---------------------------------------------------------------------------


def test_conditional_by_formula_matches_filtering_the_space():
    a = P.event(SPACE, D.CONDITIONING_EVENT_A)
    b = P.event(SPACE, D.CONDITIONING_EVENT_B)
    p_ab = P.probability(a & b, SPACE)
    p_b = P.probability(b, SPACE)
    by_formula = P.conditional(p_ab, p_b)

    restricted_event = P.event(b, D.CONDITIONING_EVENT_A)
    by_filtering = P.probability(restricted_event, b)

    assert by_formula == by_filtering == Fraction(1, 6)


@pytest.mark.parametrize("sum_target", [4, 6, 8, 10])
def test_conditioning_on_first_die_even_matches_filtering_for_several_sums(sum_target):
    b = P.event(SPACE, D.first_die_even)
    a = P.event(SPACE, D.is_sum(sum_target))
    p_ab = P.probability(a & b, SPACE)
    p_b = P.probability(b, SPACE)
    by_formula = P.conditional(p_ab, p_b)
    by_filtering = P.probability(P.event(b, D.is_sum(sum_target)), b)
    assert by_formula == by_filtering


def test_conditioning_on_the_whole_space_changes_nothing():
    a = P.event(SPACE, D.is_sum(7))
    p_a = P.probability(a, SPACE)
    p_a_and_space = P.probability(a & frozenset(SPACE), SPACE)
    p_space = P.probability(SPACE, SPACE)
    assert P.conditional(p_a_and_space, p_space) == p_a


# ---------------------------------------------------------------------------
# Exercise 7 -- the law of total probability
# ---------------------------------------------------------------------------


def test_total_probability_matches_the_documented_answer():
    assert P.total_probability(D.URN_PRIOR, D.URN_CONDITIONAL_RED) == Fraction(9, 20)


def test_total_probability_matches_the_combined_enumeration():
    urn1 = ["red"] * D.URN_1_RED + ["blue"] * D.URN_1_BLUE
    urn2 = ["red"] * D.URN_2_RED + ["blue"] * D.URN_2_BLUE
    combined = [("urn1", b) for b in urn1] + [("urn2", b) for b in urn2]
    reds = [o for o in combined if o[1] == "red"]
    enumerated = Fraction(len(reds), len(combined))
    assert P.total_probability(D.URN_PRIOR, D.URN_CONDITIONAL_RED) == enumerated


def test_total_probability_with_a_certain_prior_reduces_to_one_conditional():
    priors = (Fraction(1), Fraction(0))
    conditionals = (Fraction(3, 10), Fraction(9, 10))
    assert P.total_probability(priors, conditionals) == Fraction(3, 10)


def test_total_probability_of_blue_plus_red_is_one():
    urn_conditional_blue = (
        Fraction(D.URN_1_BLUE, 10),
        Fraction(D.URN_2_BLUE, 10),
    )
    p_red = P.total_probability(D.URN_PRIOR, D.URN_CONDITIONAL_RED)
    p_blue = P.total_probability(D.URN_PRIOR, urn_conditional_blue)
    assert p_red + p_blue == 1


# ---------------------------------------------------------------------------
# Exercise 8 -- Monte Carlo error scaling
# ---------------------------------------------------------------------------


def test_standard_error_shrinks_as_n_grows():
    small = D.standard_error(0.5, 100)
    large = D.standard_error(0.5, 100_000)
    assert large < small


def test_standard_error_of_a_certain_event_is_zero():
    assert D.standard_error(1.0, 1000) == 0.0
    assert D.standard_error(0.0, 1000) == 0.0


def test_standard_error_is_maximised_at_p_one_half():
    assert D.standard_error(0.5, 1000) > D.standard_error(0.1, 1000)
    assert D.standard_error(0.5, 1000) > D.standard_error(0.9, 1000)


def test_monte_carlo_mean_error_decreases_across_the_four_sample_sizes():
    target = float(D.MONTE_CARLO_TARGET)
    means = []
    for n in D.MONTE_CARLO_SAMPLE_SIZES:
        errors = [
            abs(S.simulate_sum_seven(np.random.default_rng(seed), n) - target)
            for seed in D.MONTE_CARLO_SEEDS
        ]
        means.append(sum(errors) / len(errors))
    assert all(means[i + 1] < means[i] for i in range(len(means) - 1))


def test_monte_carlo_error_shrink_is_closer_to_sqrt_n_than_to_n():
    target = float(D.MONTE_CARLO_TARGET)
    n_small, n_large = D.MONTE_CARLO_SAMPLE_SIZES[0], D.MONTE_CARLO_SAMPLE_SIZES[-1]

    def mean_error(n):
        errors = [
            abs(S.simulate_sum_seven(np.random.default_rng(seed), n) - target)
            for seed in D.MONTE_CARLO_SEEDS
        ]
        return sum(errors) / len(errors)

    ratio = mean_error(n_small) / mean_error(n_large)
    n_ratio = n_large / n_small
    sqrt_prediction = n_ratio**0.5
    assert 0.33 * sqrt_prediction < ratio < 3.0 * sqrt_prediction
    assert ratio < n_ratio / 10.0


def test_simulate_sum_seven_is_close_to_one_sixth_at_a_large_sample():
    rng = np.random.default_rng(123)
    got = S.simulate_sum_seven(rng, 500_000)
    assert abs(got - 1.0 / 6.0) < 4.0 * D.standard_error(1.0 / 6.0, 500_000)


# ---------------------------------------------------------------------------
# Exercise 9 -- reproducibility
# ---------------------------------------------------------------------------


def test_same_seed_gives_byte_identical_results():
    a1 = S.simulate_sum_seven(np.random.default_rng(D.REPRODUCIBILITY_SEED_A), 10_000)
    a2 = S.simulate_sum_seven(np.random.default_rng(D.REPRODUCIBILITY_SEED_A), 10_000)
    assert a1 == a2


def test_different_seeds_give_different_results():
    a = S.simulate_sum_seven(np.random.default_rng(D.REPRODUCIBILITY_SEED_A), 10_000)
    b = S.simulate_sum_seven(np.random.default_rng(D.REPRODUCIBILITY_SEED_B), 10_000)
    assert a != b


def test_different_seeds_both_still_estimate_the_true_value():
    target = 1.0 / 6.0
    tol = 4.0 * D.standard_error(target, D.REPRODUCIBILITY_TRIALS)
    for seed in (D.REPRODUCIBILITY_SEED_A, D.REPRODUCIBILITY_SEED_B):
        got = S.simulate_sum_seven(np.random.default_rng(seed), D.REPRODUCIBILITY_TRIALS)
        assert abs(got - target) < tol


@pytest.mark.parametrize("seed", [0, 1, 2, 3, 4])
def test_reproducibility_holds_across_several_seeds(seed):
    a1 = S.simulate_at_least_one_six(np.random.default_rng(seed), 5_000)
    a2 = S.simulate_at_least_one_six(np.random.default_rng(seed), 5_000)
    assert a1 == a2


def test_two_generators_from_the_same_seed_produce_identical_raw_draws():
    rng1 = np.random.default_rng(7)
    rng2 = np.random.default_rng(7)
    draws1 = rng1.integers(1, 7, size=100)
    draws2 = rng2.integers(1, 7, size=100)
    assert np.array_equal(draws1, draws2)


def test_a_single_generator_does_not_repeat_its_own_draws():
    rng = np.random.default_rng(7)
    first = rng.integers(1, 7, size=100)
    second = rng.integers(1, 7, size=100)
    assert not np.array_equal(first, second)

lesson_id: D113
day: 113
kind: guided-build
languages: [python, bash]
setup_commands:
  - cd labs/sections/math-statistics-and-data/day-113-probability-events-rules-and-intuition
  - python3 -m venv .venv
  - .venv/bin/pip install -r requirements/requirements.txt
  - .venv/bin/python3 -c "import numpy; print(numpy.__version__)"
run_commands:
  - 'cd examples && ../.venv/bin/python3 01_sample_space_and_events.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 02_addition_rule.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 03_de_mere.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 04_independence_vs_dependence.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 05_mutual_exclusivity_implies_dependence.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 06_conditioning_by_restriction.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 07_law_of_total_probability.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 08_monte_carlo_error_scaling.py && cd ..'
  - 'cd examples && ../.venv/bin/python3 09_reproducibility.py && cd ..'
  - .venv/bin/pytest examples -q -p no:cacheprovider
  - .venv/bin/pytest starter -q -p no:cacheprovider
test_commands:
  - bash tests/run_tests.sh
cleanup_commands:
  - "find . -path ./.venv -prune -o -type d -name '__pycache__' -print -exec rm -rf -- {} +"
  - rm -rf .pytest_cache
  - 'rm -rf .venv  # optional: removes the lab virtual environment'
  - 'git checkout -- starter/  # optional: reset your work'
requires_network: true
requires_api_key: false
estimated_minutes: 30
last_executed: '2026-08-17'
executed_on: 'macOS 26.5.2 (Apple Silicon, arm64), Python 3.14.0, numpy 2.5.2, pytest 9.1.1, bash 3.2.57 -- bash tests/run_tests.sh -> 57 checks, 0 failure(s), exit 0; pytest examples -> 93 passed; pytest starter -> 3 passed, 43 skipped on an untouched checkout, and 46 passed against a fully solved copy of starter/ (verified by temporarily copying the reference probability.py and simulate.py into starter/ and filling in every answers.py prediction, then restoring the blank skeletons -- the skip counts before and after that restore were confirmed identical). All nine reference scripts exit 0 with every internal assertion holding. Everything was run through a real lab-local .venv created by the documented setup commands, not through an authoring environment; pip install used the network exactly once, as documented. Section 6 of the harness re-runs the reference pytest suite with one assertion (the addition-rule error-amount check) temporarily replaced by a deliberately wrong one, confirms the run exits non-zero with exactly one failure named in the output, and restores the original file -- so the suite is demonstrated to be capable of failing rather than merely claimed to be. Two honesty notes from this run. FIRST: scipy and pandas are not installed in this environment. scipy.stats is described from its public documentation in the lesson''s Tools section and explicitly marked as not run here; no output attributed to it anywhere in this lab or its lesson was actually produced by it. SECOND: the Monte Carlo error-scaling exercise (08) reports real, freshly measured numbers rather than fixed literals -- on this run the mean absolute error across 20 seeds fell from 0.023000 at n=100 to 0.000977 at n=100,000, a 23.55x shrink against a sqrt(1000)=31.62x prediction, comfortably within the suite''s stated band of 0.33x to 3.0x the sqrt(n) prediction and nowhere near the false 1/n prediction of 1000x; these exact figures will differ slightly on another machine or NumPy version, and the tests assert the trend shape and the tolerance band rather than these specific digits, as documented in expected-output/FIELDS.md. De Méré''s two exact probabilities, 671/1296 and the double-six figure, are exact rational arithmetic via fractions.Fraction and are therefore identical on any correct Python implementation; only the simulated confirmations of them carry machine-dependent noise, and both landed within their three-standard-error tolerance (0.00335 at 200,000 trials) on this run: bet 1 simulated at 0.515865 against an exact 0.517747, bet 2 at 0.492350 against an exact 0.491404.'

# What is installed, why, and what it costs

Two packages, both free and open source, both installed into a lab-local
virtual environment that `rm -rf .venv` completely undoes.

| Package | Version pinned | Licence | What this lab uses it for |
| --- | --- | --- | --- |
| `numpy` | 2.5.2 | BSD 3-Clause | The `numpy.random.Generator` built by `default_rng(seed)`, used for every simulation in exercises 3, 8 and 9 — de Méré's two bets, the Monte Carlo error-scaling sweep, and the reproducibility checks. |
| `pytest` | 9.1.1 | MIT | The reference suite (93 tests) and your running score in `starter/`. |

There is no paid tier of anything in this lab, no account, no key and no
signup, personally or commercially.

## The one time the network is needed

```bash
.venv/bin/pip install -r requirements/requirements.txt
```

That is the only command in the lab that opens a connection. Section 5 of
`tests/run_tests.sh` greps every source file in `examples/` and `starter/` to
prove that nothing else does.

## If you cannot install anything at all

Every exact-probability exercise — 1, 2, 4, 5, 6 and 7 — needs only
`itertools`, `fractions` and `statistics` from the standard library, and
none of them touch NumPy. Only the three simulation exercises (3, 8, 9) need
`numpy.random.Generator`. If NumPy is unavailable, Python's own `random`
module can stand in for the vectorised roll:

```python
import random

def simulate_sum_seven_stdlib(seed: int, trials: int) -> float:
    rng = random.Random(seed)
    hits = sum(
        1 for _ in range(trials)
        if rng.randint(1, 6) + rng.randint(1, 6) == 7
    )
    return hits / trials
```

It is slower — a Python-level loop instead of a vectorised NumPy call — but
it is the same statistics. `troubleshooting.md` shows the substitution for
the other two simulation functions. What you lose is `pytest`, so no running
score and no skip-versus-fail distinction; you would read the numbers back
yourself instead.

## What is deliberately *not* installed

`scipy.stats` does this job too — and considerably more, once random
variables and distributions arrive on Day 114. It is **not installed in this
environment, and no output from it is reproduced anywhere** in this lab or
its lesson. The lesson's Tools section describes it from its documentation
and marks it as not run here.

`pandas` is also not installed and is not needed for anything in this lab;
every table here is small enough to enumerate and print directly.

That is not a limitation to apologise for. Every probability this lab
computes is either counted exactly from an enumerated sample space or
estimated by a simulation you wrote yourself with the standard library and
NumPy's random module — nothing here depends on a statistics package to be
correct.

numpy==2.5.2
pytest==9.1.1

# The nine exercises

Work through these in order. Predict the answer to each `answers.py` question
*before* running anything — the two that catch almost everyone (the
addition-rule error amount and de Méré's favourable bet) only catch you if
you commit to a guess first.

Check yourself as you go:

```bash
.venv/bin/pytest starter -q
```

Unattempted work reports as **skipped**, never failed. Wrong work **fails**
with your answer printed beside the correct one.

## 1. The sample space (`probability.py`)

Write `sample_space_two_dice()` with `itertools.product`, `event()` to
filter a space by a predicate, and `probability()` to return the exact
`Fraction` of an event. Assert `P(sum == 7) == Fraction(1, 6)` — not
`0.16666...`, the exact fraction.

## 2. The addition rule (`probability.py`)

`addition_rule(p_a, p_b, p_a_and_b)` and `naive_sum(p_a, p_b)`. Verify on
`A = "sum is 7"`, `B = "first die is 6"` that the naive sum is wrong by
exactly `P(A and B)`.

## 3. De Méré, exact and simulated (`probability.py`, `simulate.py`)

`complement()` and `at_least_one()` collapse both of de Méré's bets to one
line each. Then `simulate_at_least_one_six()` and
`simulate_at_least_one_double_six()` confirm both by simulation, within a
tolerance derived from the standard error of a proportion.

## 4. Independence (`probability.py`)

`is_independent(p_a, p_b, p_a_and_b)`. Confirmed against one genuinely
independent pair of dice events and one genuinely dependent pair.

## 5. Mutual exclusivity implies dependence

Uses `conditional()` from exercise 6 on a pair of mutually exclusive events.
No new function — the point is what the existing tools reveal.

## 6. Conditioning by restriction (`probability.py`)

`conditional(p_a_and_b, p_b)`. Computed by formula and by literally filtering
the sample space down to the rows where the condition holds; both must agree
exactly.

## 7. The law of total probability (`probability.py`)

`total_probability(priors, conditionals)`, checked against a direct
enumeration of the combined two-urn experiment.

## 8. Monte Carlo error scaling (`simulate.py`)

`simulate_sum_seven()`, called at four sample sizes four decades apart,
averaged over twenty seeds at each. The error should shrink like
`1/sqrt(n)`, not `1/n`.

## 9. Reproducibility (`simulate.py`)

The same `numpy.random.default_rng(seed)` gives byte-identical results
across two calls; a different seed gives a different but still-close result.

"""Exercises 1 through 9 -- eighteen predictions.

Replace each `None` with the value you think is correct. A `None` is a skip,
not a failure: `pytest starter -q` counts only what you have attempted. When
you are wrong it prints both your answer and the real one, so a wrong guess
is worth more than a blank.

Predict BEFORE you run anything. The two that catch almost everyone are the
addition-rule error amount (exercise 2) and de Méré's favourable bet
(exercise 3) -- and they only catch you if you commit first.

Every answer is a number, a Python bool, or (for one question) an int naming
a bet.
"""

ANSWERS: dict[str, object] = {
    # ----------------------------------------------------------------------
    # Exercise 1 -- the sample space
    # ----------------------------------------------------------------------
    # 1.1 How many outcomes are in the sample space of two dice?
    "sample_space_size": None,
    # 1.2 P(the two dice sum to 7), as a decimal.
    "p_sum_seven": None,
    # ----------------------------------------------------------------------
    # Exercise 2 -- the addition rule
    # ----------------------------------------------------------------------
    # 2.1 A = "sum is 7" (6 outcomes), B = "first die is 6" (6 outcomes).
    #     The WRONG naive sum P(A) + P(B), as a decimal.
    "addition_naive_sum": None,
    # 2.2 The TRUE P(A or B), as a decimal.
    "addition_true_union": None,
    # 2.3 By exactly how much does the naive sum overstate the truth?
    "addition_error_amount": None,
    # ----------------------------------------------------------------------
    # Exercise 3 -- de Méré's two bets
    # ----------------------------------------------------------------------
    # 3.1 P(at least one 6 in 4 rolls of one die), as a decimal.
    "de_mere_single_bet_probability": None,
    # 3.2 P(at least one double-six in 24 rolls of two dice), as a decimal.
    "de_mere_double_bet_probability": None,
    # 3.3 Which bet is favourable to the player (probability above 0.5)?
    #     Answer 1 or 2.
    "de_mere_favorable_bet": None,
    # ----------------------------------------------------------------------
    # Exercise 4 -- independence
    # ----------------------------------------------------------------------
    # 4.1 Does P(A and B) == P(A) * P(B) hold for the INDEPENDENT_PAIR?
    "independent_pair_holds": None,
    # 4.2 Does it hold for the DEPENDENT_PAIR?
    "dependent_pair_holds": None,
    # ----------------------------------------------------------------------
    # Exercise 5 -- mutual exclusivity implies dependence
    # ----------------------------------------------------------------------
    # 5.1 For the mutually exclusive pair (sum = 2, sum = 12), is
    #     P(sum=2 | sum=12) == 0 while P(sum=2) != 0 -- i.e. are they
    #     dependent?
    "mutually_exclusive_implies_dependent": None,
    # ----------------------------------------------------------------------
    # Exercise 6 -- conditioning by restriction
    # ----------------------------------------------------------------------
    # 6.1 P(sum = 8 | first die is even), as a decimal.
    "conditional_p_sum8_given_first_even": None,
    # 6.2 Does the formula method and the filter-the-space method agree
    #     exactly?
    "conditional_formula_matches_filter": None,
    # ----------------------------------------------------------------------
    # Exercise 7 -- the law of total probability
    # ----------------------------------------------------------------------
    # 7.1 P(red), across both urns, weighted by the fair coin. As a decimal.
    "urn_total_probability_red": None,
    # 7.2 Does the weighted-total answer match the answer from enumerating
    #     the combined 20-outcome experiment directly?
    "urn_enumeration_matches_formula": None,
    # ----------------------------------------------------------------------
    # Exercise 8 -- Monte Carlo error scaling
    # ----------------------------------------------------------------------
    # 8.1 Does the average simulation error shrink as the sample size grows
    #     from 100 to 100,000?
    "monte_carlo_error_shrinks_with_n": None,
    # 8.2 Multiplying n by 10 should shrink the error by roughly sqrt(10),
    #     not by 10. Is the OBSERVED shrink much closer to sqrt(10) than
    #     to 10?
    "monte_carlo_error_ratio_near_sqrt10": None,
    # ----------------------------------------------------------------------
    # Exercise 9 -- reproducibility
    # ----------------------------------------------------------------------
    # 9.1 Two simulations with the SAME seed: byte-identical results?
    "reproducibility_same_seed_identical": None,
    # 9.2 Two simulations with DIFFERENT seeds: different results that both
    #     still land within tolerance of the true probability?
    "reproducibility_different_seed_differs": None,
}

"""Make this directory's own modules the ones its tests import.

Both `examples/` and `starter/` contain modules called `probability`,
`simulate`, `dataset` and `answers`, and pytest imports test files by putting
their directory on `sys.path`. Without this file, running `pytest` across both
directories at once would import whichever `probability` was seen first and
then reuse it for the other suite -- so these starter tests would silently
pass against the reference solution instead of skipping. That is a wrong
answer with a green tick on it, which is the worst kind.

So: put this directory first on the import path, and drop any already-imported
module of those names that came from somewhere else.
"""

import sys
from pathlib import Path

HERE = str(Path(__file__).parent.resolve())

if HERE in sys.path:
    sys.path.remove(HERE)
sys.path.insert(0, HERE)

for name in ("probability", "simulate", "dataset", "answers"):
    module = sys.modules.get(name)
    origin = getattr(module, "__file__", "") or ""
    if module is not None and not origin.startswith(HERE):
        del sys.modules[name]

"""The sample space, the events, the urns, and every tolerance this lab
compares against.

Read this file. Nothing here is tuned: every probability below is either
computed by enumeration (via `itertools.product` and `fractions.Fraction`) or
derived algebraically and then checked against the enumeration in the tests.
The simulation tolerance is derived from the standard error of a proportion,
`sqrt(p(1-p)/n)`, with the arithmetic written out beside it -- not chosen by
running a test and loosening the number until it passed.
"""

import itertools
import math
from fractions import Fraction
from typing import Callable

# --------------------------------------------------------------------------
# The sample space: two fair six-sided dice, 36 equally likely outcomes
# --------------------------------------------------------------------------

#: Every ordered pair (first die, second die), each equally likely. Built by
#: enumeration, not written down as a literal list.
TWO_DICE_SPACE: tuple[tuple[int, int], ...] = tuple(
    itertools.product(range(1, 7), range(1, 7))
)

assert len(TWO_DICE_SPACE) == 36

# --------------------------------------------------------------------------
# Events, as predicates over one outcome -- the enumerated set is built by
# filtering the space with these, in probability.py and in the scripts below.
# --------------------------------------------------------------------------


def is_sum(target: int) -> Callable[[tuple[int, int]], bool]:
    """An event predicate: the two dice sum to `target`."""

    def predicate(outcome: tuple[int, int]) -> bool:
        return outcome[0] + outcome[1] == target

    return predicate


def is_first_die(value: int) -> Callable[[tuple[int, int]], bool]:
    """An event predicate: the first die shows `value`."""

    def predicate(outcome: tuple[int, int]) -> bool:
        return outcome[0] == value

    return predicate


def first_die_even(outcome: tuple[int, int]) -> bool:
    """An event predicate: the first die shows an even number."""
    return outcome[0] % 2 == 0


def is_double(outcome: tuple[int, int]) -> bool:
    """An event predicate: both dice show the same value."""
    return outcome[0] == outcome[1]


# --------------------------------------------------------------------------
# The addition rule, on a concrete pair of events
# --------------------------------------------------------------------------

#: A = "the dice sum to 7" (6 outcomes). B = "the first die shows 6" (6
#: outcomes). Their overlap is exactly one outcome, (6, 1), because a first
#: die of 6 needs a second die of 1 to sum to 7. The naive sum P(A) + P(B)
#: double-counts that one outcome, so it overstates P(A union B) by exactly
#: P(A intersect B) = 1/36.
ADDITION_EVENT_A = is_sum(7)
ADDITION_EVENT_B = is_first_die(6)

# --------------------------------------------------------------------------
# The Chevalier de Méré's two bets
# --------------------------------------------------------------------------

#: Bet 1: at least one six in four rolls of one die.
DE_MERE_SINGLE_ROLLS: int = 4
DE_MERE_SINGLE_FACES: int = 6

#: Bet 2: at least one double-six in twenty-four rolls of two dice.
DE_MERE_DOUBLE_ROLLS: int = 24
DE_MERE_DOUBLE_FACES: int = 36  # 6 x 6 possible pairs, one of which is (6, 6)

#: Exact answers, by the complement rule: 1 - P(none of the trials succeed).
#: Fraction keeps every digit; the lesson rounds these to four places.
DE_MERE_SINGLE_EXACT: Fraction = 1 - Fraction(5, 6) ** DE_MERE_SINGLE_ROLLS
DE_MERE_DOUBLE_EXACT: Fraction = 1 - Fraction(35, 36) ** DE_MERE_DOUBLE_ROLLS

assert round(float(DE_MERE_SINGLE_EXACT), 4) == 0.5177
assert round(float(DE_MERE_DOUBLE_EXACT), 4) == 0.4914

# --------------------------------------------------------------------------
# Independence and dependence: two named pairs of events in TWO_DICE_SPACE
# --------------------------------------------------------------------------

#: A genuinely independent pair. "Sum is 7" is independent of "first die is
#: 3" -- for ANY value the first die shows, exactly one value of the second
#: die makes the sum 7, so P(sum = 7 | first die = v) = 1/6 for every v, which
#: is P(sum = 7) itself. That is what independence means.
INDEPENDENT_PAIR = (is_sum(7), is_first_die(3))

#: A genuinely dependent pair. "Sum is 2" requires BOTH dice to show 1, so it
#: is only possible when the first die shows 1 -- P(sum = 2 | first die = 1)
#: = 1/6, but P(sum = 2 | first die != 1) = 0. Those are not equal, so the
#: events are dependent.
DEPENDENT_PAIR = (is_sum(2), is_first_die(1))

# --------------------------------------------------------------------------
# Mutual exclusivity implies dependence
# --------------------------------------------------------------------------

#: "Sum is 2" and "sum is 12" cannot both happen (2 needs (1,1), 12 needs
#: (6,6)), so they are mutually exclusive. Knowing one occurred makes the
#: other impossible -- the sharpest form of dependence there is.
MUTUALLY_EXCLUSIVE_PAIR = (is_sum(2), is_sum(12))

# --------------------------------------------------------------------------
# Conditioning by restriction
# --------------------------------------------------------------------------

#: P(sum = 8 | first die is even), computed two ways in the lab: by the
#: formula P(A and B) / P(B), and by filtering TWO_DICE_SPACE down to the
#: rows where the first die is even and asking what fraction of THOSE rows
#: sum to 8.
CONDITIONING_EVENT_A = is_sum(8)
CONDITIONING_EVENT_B = first_die_even

# --------------------------------------------------------------------------
# The law of total probability: two urns, drawn from with a fair coin
# --------------------------------------------------------------------------

#: Urn 1: 3 red, 7 blue, out of 10 balls.
URN_1_RED: int = 3
URN_1_BLUE: int = 7

#: Urn 2: 6 red, 4 blue, out of 10 balls.
URN_2_RED: int = 6
URN_2_BLUE: int = 4

assert URN_1_RED + URN_1_BLUE == 10
assert URN_2_RED + URN_2_BLUE == 10

#: A fair coin decides which urn to draw from -- P(urn 1) = P(urn 2) = 1/2.
URN_PRIOR: tuple[Fraction, Fraction] = (Fraction(1, 2), Fraction(1, 2))

#: P(red | urn 1) and P(red | urn 2), read straight off each urn's contents.
URN_CONDITIONAL_RED: tuple[Fraction, Fraction] = (
    Fraction(URN_1_RED, URN_1_RED + URN_1_BLUE),
    Fraction(URN_2_RED, URN_2_RED + URN_2_BLUE),
)

# --------------------------------------------------------------------------
# Monte Carlo error scaling
# --------------------------------------------------------------------------

#: The event whose probability the Monte Carlo experiment estimates: two
#: fair dice sum to 7. The exact answer is 1/6, established in exercise 1.
MONTE_CARLO_TARGET: Fraction = Fraction(1, 6)

#: Sample sizes to sweep, four decades apart at the ends.
MONTE_CARLO_SAMPLE_SIZES: tuple[int, ...] = (100, 1_000, 10_000, 100_000)

#: Independent seeds to average over at each sample size, so the assertion is
#: about the SHAPE of the error trend across many runs rather than about one
#: sampled value, which would be flaky on someone else's machine.
MONTE_CARLO_SEEDS: tuple[int, ...] = tuple(range(20))


def standard_error(p: float, n: int) -> float:
    """The standard error of a proportion estimated from n trials.

    This is the quantity that governs how far a Monte Carlo estimate can be
    expected to land from the true probability it estimates: about 68% of
    estimates land within one standard error, and about 99.7% within three.
    """
    return math.sqrt(p * (1.0 - p) / n)


#: Three standard errors, at the sample sizes de Méré's simulations use in
#: exercise 3. This is not a guessed tolerance: it is the width inside which
#: 99.7% of repeated simulations should land, derived from the formula above.
DE_MERE_SIM_TRIALS: int = 200_000
DE_MERE_SINGLE_SE: float = standard_error(float(DE_MERE_SINGLE_EXACT), DE_MERE_SIM_TRIALS)
DE_MERE_DOUBLE_SE: float = standard_error(float(DE_MERE_DOUBLE_EXACT), DE_MERE_SIM_TRIALS)
DE_MERE_SINGLE_TOL: float = 3.0 * DE_MERE_SINGLE_SE
DE_MERE_DOUBLE_TOL: float = 3.0 * DE_MERE_DOUBLE_SE

# --------------------------------------------------------------------------
# Reproducibility
# --------------------------------------------------------------------------

REPRODUCIBILITY_SEED_A: int = 42
REPRODUCIBILITY_SEED_B: int = 43
REPRODUCIBILITY_TRIALS: int = 10_000

"""Exercise 1, 2, 4, 5, 6 and 7: exact probability, computed two ways.

Every function here returns a `fractions.Fraction` wherever the answer is
rational, so an assertion against it is exact rather than "close enough".
Fill in the bodies marked `# YOUR CODE HERE`. `dataset.py` has everything you
need: the sample space, the event predicates, and the urn compositions.
"""

import itertools
from fractions import Fraction
from typing import Callable, Iterable


# ---------------------------------------------------------------------------
# Exercise 1: the sample space, and probability as counting
# ---------------------------------------------------------------------------


def sample_space_two_dice() -> tuple[tuple[int, int], ...]:
    """Every ordered pair (first die, second die), 36 outcomes.

    Build it with `itertools.product(range(1, 7), range(1, 7))` -- do not
    write the 36 pairs out by hand.
    """
    # YOUR CODE HERE
    raise NotImplementedError


def event(
    space: Iterable[tuple[int, int]], predicate: Callable[[tuple[int, int]], bool]
) -> frozenset[tuple[int, int]]:
    """The subset of `space` for which `predicate` is true."""
    # YOUR CODE HERE
    raise NotImplementedError


def probability(
    outcome_set: Iterable[tuple[int, int]], space: Iterable[tuple[int, int]]
) -> Fraction:
    """P(event), for equally likely outcomes: |event| / |space|, exactly."""
    # YOUR CODE HERE
    raise NotImplementedError


# ---------------------------------------------------------------------------
# Exercise 2: the addition rule
# ---------------------------------------------------------------------------


def addition_rule(p_a: Fraction, p_b: Fraction, p_a_and_b: Fraction) -> Fraction:
    """P(A or B) = P(A) + P(B) - P(A and B)."""
    # YOUR CODE HERE
    raise NotImplementedError


def naive_sum(p_a: Fraction, p_b: Fraction) -> Fraction:
    """The wrong shortcut: P(A) + P(B), which double-counts the overlap."""
    # YOUR CODE HERE
    raise NotImplementedError


# ---------------------------------------------------------------------------
# Exercise 3 helper: the complement rule
# ---------------------------------------------------------------------------


def complement(p: Fraction) -> Fraction:
    """P(not A) = 1 - P(A)."""
    # YOUR CODE HERE
    raise NotImplementedError


def at_least_one(p_single_success: Fraction, trials: int) -> Fraction:
    """P(at least one success in `trials` independent tries).

    Collapse it to one line with the complement rule: 1 minus the
    probability that every single trial fails.
    """
    # YOUR CODE HERE
    raise NotImplementedError


# ---------------------------------------------------------------------------
# Exercise 4: independence
# ---------------------------------------------------------------------------


def is_independent(p_a: Fraction, p_b: Fraction, p_a_and_b: Fraction) -> bool:
    """True exactly when P(A and B) == P(A) * P(B)."""
    # YOUR CODE HERE
    raise NotImplementedError


# ---------------------------------------------------------------------------
# Exercise 5 and 6: conditional probability
# ---------------------------------------------------------------------------


def conditional(p_a_and_b: Fraction, p_b: Fraction) -> Fraction:
    """P(A | B) = P(A and B) / P(B)."""
    # YOUR CODE HERE
    raise NotImplementedError


# ---------------------------------------------------------------------------
# Exercise 7: the law of total probability
# ---------------------------------------------------------------------------


def total_probability(
    priors: Iterable[Fraction], conditionals: Iterable[Fraction]
) -> Fraction:
    """P(A) = sum over i of P(A | condition_i) * P(condition_i).

    `priors` and `conditionals` are matched by position: priors[i] is
    P(condition_i), conditionals[i] is P(A | condition_i).
    """
    # YOUR CODE HERE
    raise NotImplementedError

"""Exercise 3, 8 and 9: simulation, and the error that comes with it.

Every simulation here takes an explicit `numpy.random.Generator` -- built by
`numpy.random.default_rng(seed)` -- rather than touching global random state.
Two calls with the same seed must return byte-identical arrays; that is
`default_rng`'s whole advantage over the legacy `numpy.random.seed` global,
and exercise 9 asserts it directly.
"""

import numpy as np


# ---------------------------------------------------------------------------
# Exercise 3: de Méré's two bets, simulated
# ---------------------------------------------------------------------------


def simulate_at_least_one_six(rng: np.random.Generator, trials: int) -> float:
    """Simulate de Méré's first bet: at least one 6 in 4 rolls of one die.

    Roll a (trials, 4) array of dice, and return the fraction of rows that
    contain at least one 6. Vectorised: no Python-level loop over trials.
    """
    # YOUR CODE HERE
    raise NotImplementedError


def simulate_at_least_one_double_six(rng: np.random.Generator, trials: int) -> float:
    """Simulate de Méré's second bet: at least one double-six in 24 rolls of
    two dice.

    Roll two (trials, 24) arrays -- one per die -- and return the fraction of
    rows in which some roll shows 6 on both dice at once.
    """
    # YOUR CODE HERE
    raise NotImplementedError


# ---------------------------------------------------------------------------
# Exercise 8: Monte Carlo error scaling
# ---------------------------------------------------------------------------


def simulate_sum_seven(rng: np.random.Generator, trials: int) -> float:
    """Simulate P(two fair dice sum to 7) by rolling `trials` pairs of dice.

    Returns the fraction of pairs that summed to 7. The exact answer,
    established in exercise 1, is exactly 1/6.
    """
    # YOUR CODE HERE
    raise NotImplementedError

"""Your running score. Unattempted work SKIPS; wrong work FAILS with both values.

Run from the lab directory:

    .venv/bin/pytest starter -q

On an untouched checkout this reports one pass and everything else skipped.
A skip means "not attempted". A failure means "attempted and wrong", and the
message shows your answer next to the real one so you can see the gap rather
than guess at it.

Nothing in here checks that a function exists or that a file is present.
Every test runs your code and compares a value.
"""

import itertools
from fractions import Fraction

import numpy as np
import pytest

import answers
import dataset as D
import probability as P
import simulate as S

# --------------------------------------------------------------------------
# The skip machinery
# --------------------------------------------------------------------------


def need(value, what):
    """Skip if the exercise has not been attempted, otherwise hand it back."""
    if value is None:
        pytest.skip(f"not attempted yet: {what}")
    return value


def attempt(fn, what):
    """Call something that may not be written yet, and skip if it is not."""
    try:
        result = fn()
    except (TypeError, AttributeError, NotImplementedError):
        pytest.skip(f"not attempted yet: {what}")
    if result is None:
        pytest.skip(f"not attempted yet: {what}")
    return result


def close(got, want, tol, what):
    assert abs(float(got) - float(want)) < tol, (
        f"{what}: your answer {got!r}, expected {want!r} "
        f"(difference {abs(float(got) - float(want)):.3e}, tolerance {tol:g})"
    )


def test_the_suite_itself_runs():
    """One test that always passes, so a green run is distinguishable from
    a collection error that quietly ran nothing at all."""
    assert len(D.TWO_DICE_SPACE) == 36


# --------------------------------------------------------------------------
# Exercise 1 -- the sample space and probability as counting
# --------------------------------------------------------------------------


def test_1_sample_space_has_36_outcomes():
    space = attempt(P.sample_space_two_dice, "sample_space_two_dice")
    assert len(space) == 36


def test_1_sample_space_matches_the_dataset_space():
    space = attempt(P.sample_space_two_dice, "sample_space_two_dice")
    assert set(space) == set(D.TWO_DICE_SPACE)


def test_1_event_filters_by_predicate():
    space = attempt(P.sample_space_two_dice, "sample_space_two_dice")
    ev = attempt(lambda: P.event(space, D.is_sum(7)), "event")
    assert len(ev) == 6
    assert all(a + b == 7 for a, b in ev)


def test_1_probability_of_sum_seven_is_exactly_one_sixth():
    space = attempt(P.sample_space_two_dice, "sample_space_two_dice")
    ev = attempt(lambda: P.event(space, D.is_sum(7)), "event")
    p = attempt(lambda: P.probability(ev, space), "probability")
    assert p == Fraction(1, 6), "P(sum == 7) must be exactly Fraction(1, 6)"


def test_1_probability_returns_a_fraction():
    space = attempt(P.sample_space_two_dice, "sample_space_two_dice")
    ev = attempt(lambda: P.event(space, D.is_sum(7)), "event")
    p = attempt(lambda: P.probability(ev, space), "probability")
    assert isinstance(p, Fraction), (
        "probability() must return a Fraction, not a float -- that is the "
        "whole point of using it here"
    )


def test_1_probability_of_the_whole_space_is_one():
    space = attempt(P.sample_space_two_dice, "sample_space_two_dice")
    p = attempt(lambda: P.probability(space, space), "probability")
    assert p == 1


def test_1_probability_of_the_empty_set_is_zero():
    space = attempt(P.sample_space_two_dice, "sample_space_two_dice")
    p = attempt(lambda: P.probability(frozenset(), space), "probability")
    assert p == 0


# --------------------------------------------------------------------------
# Exercise 2 -- the addition rule
# --------------------------------------------------------------------------


def _addition_setup():
    space = D.TWO_DICE_SPACE
    a = frozenset(o for o in space if D.ADDITION_EVENT_A(o))
    b = frozenset(o for o in space if D.ADDITION_EVENT_B(o))
    p_a = Fraction(len(a), len(space))
    p_b = Fraction(len(b), len(space))
    p_ab = Fraction(len(a & b), len(space))
    return p_a, p_b, p_ab


def test_2_addition_rule_gives_the_true_union():
    p_a, p_b, p_ab = _addition_setup()
    got = attempt(lambda: P.addition_rule(p_a, p_b, p_ab), "addition_rule")
    assert got == Fraction(11, 36)


def test_2_naive_sum_double_counts_the_overlap():
    p_a, p_b, p_ab = _addition_setup()
    naive = attempt(lambda: P.naive_sum(p_a, p_b), "naive_sum")
    assert naive == Fraction(1, 3), "P(A) + P(B) should be 6/36 + 6/36 = 1/3"


def test_2_the_naive_sum_overstates_by_exactly_the_intersection():
    p_a, p_b, p_ab = _addition_setup()
    naive = attempt(lambda: P.naive_sum(p_a, p_b), "naive_sum")
    true_union = attempt(lambda: P.addition_rule(p_a, p_b, p_ab), "addition_rule")
    assert naive - true_union == p_ab, (
        "the naive sum's error should be exactly P(A and B) = 1/36"
    )


# --------------------------------------------------------------------------
# Exercise 3 -- de Méré, exact and simulated
# --------------------------------------------------------------------------


def test_3a_complement_of_one_sixth_is_five_sixths():
    got = attempt(lambda: P.complement(Fraction(1, 6)), "complement")
    assert got == Fraction(5, 6)


def test_3b_at_least_one_six_matches_the_exact_de_mere_answer():
    got = attempt(
        lambda: P.at_least_one(Fraction(1, 6), D.DE_MERE_SINGLE_ROLLS), "at_least_one"
    )
    assert got == D.DE_MERE_SINGLE_EXACT


def test_3b_at_least_one_double_six_matches_the_exact_de_mere_answer():
    got = attempt(
        lambda: P.at_least_one(Fraction(1, 36), D.DE_MERE_DOUBLE_ROLLS),
        "at_least_one",
    )
    assert got == D.DE_MERE_DOUBLE_EXACT


def test_3c_bet_one_is_favourable_and_bet_two_is_not():
    assert D.DE_MERE_SINGLE_EXACT > Fraction(1, 2)
    assert D.DE_MERE_DOUBLE_EXACT < Fraction(1, 2)


def test_3d_simulation_of_bet_one_lands_within_three_standard_errors():
    rng = np.random.default_rng(D.REPRODUCIBILITY_SEED_A)
    got = attempt(
        lambda: S.simulate_at_least_one_six(rng, D.DE_MERE_SIM_TRIALS),
        "simulate_at_least_one_six",
    )
    close(got, float(D.DE_MERE_SINGLE_EXACT), D.DE_MERE_SINGLE_TOL, "de Méré bet 1")


def test_3d_simulation_of_bet_two_lands_within_three_standard_errors():
    rng = np.random.default_rng(D.REPRODUCIBILITY_SEED_A)
    got = attempt(
        lambda: S.simulate_at_least_one_double_six(rng, D.DE_MERE_SIM_TRIALS),
        "simulate_at_least_one_double_six",
    )
    close(got, float(D.DE_MERE_DOUBLE_EXACT), D.DE_MERE_DOUBLE_TOL, "de Méré bet 2")


# --------------------------------------------------------------------------
# Exercise 4 -- independence
# --------------------------------------------------------------------------


def test_4_the_independent_pair_is_reported_independent():
    a_pred, b_pred = D.INDEPENDENT_PAIR
    space = D.TWO_DICE_SPACE
    a = frozenset(o for o in space if a_pred(o))
    b = frozenset(o for o in space if b_pred(o))
    p_a, p_b = Fraction(len(a), 36), Fraction(len(b), 36)
    p_ab = Fraction(len(a & b), 36)
    got = attempt(lambda: P.is_independent(p_a, p_b, p_ab), "is_independent")
    assert got is True


def test_4_the_dependent_pair_is_reported_dependent():
    a_pred, b_pred = D.DEPENDENT_PAIR
    space = D.TWO_DICE_SPACE
    a = frozenset(o for o in space if a_pred(o))
    b = frozenset(o for o in space if b_pred(o))
    p_a, p_b = Fraction(len(a), 36), Fraction(len(b), 36)
    p_ab = Fraction(len(a & b), 36)
    got = attempt(lambda: P.is_independent(p_a, p_b, p_ab), "is_independent")
    assert got is False


# --------------------------------------------------------------------------
# Exercise 5 -- mutual exclusivity implies dependence
# --------------------------------------------------------------------------


def test_5_mutually_exclusive_events_have_zero_conditional():
    a_pred, b_pred = D.MUTUALLY_EXCLUSIVE_PAIR
    space = D.TWO_DICE_SPACE
    a = frozenset(o for o in space if a_pred(o))
    b = frozenset(o for o in space if b_pred(o))
    p_a = Fraction(len(a), 36)
    p_ab = Fraction(len(a & b), 36)
    p_b = Fraction(len(b), 36)
    got = attempt(lambda: P.conditional(p_ab, p_b), "conditional")
    assert got == 0
    assert p_a != 0, "P(A) must be non-zero for this to demonstrate dependence"
    assert got != p_a, "P(A | B) == 0 while P(A) != 0 -- the events are dependent"


# --------------------------------------------------------------------------
# Exercise 6 -- conditioning by restriction
# --------------------------------------------------------------------------


def test_6_conditional_by_formula_matches_conditional_by_filtering():
    space = D.TWO_DICE_SPACE
    a = frozenset(o for o in space if D.CONDITIONING_EVENT_A(o))
    b = frozenset(o for o in space if D.CONDITIONING_EVENT_B(o))
    p_ab = Fraction(len(a & b), 36)
    p_b = Fraction(len(b), 36)
    by_formula = attempt(lambda: P.conditional(p_ab, p_b), "conditional")

    restricted = frozenset(o for o in b if D.CONDITIONING_EVENT_A(o))
    by_filtering = Fraction(len(restricted), len(b))

    assert by_formula == by_filtering == Fraction(1, 6)


# --------------------------------------------------------------------------
# Exercise 7 -- the law of total probability
# --------------------------------------------------------------------------


def test_7_total_probability_matches_the_combined_enumeration():
    got = attempt(
        lambda: P.total_probability(D.URN_PRIOR, D.URN_CONDITIONAL_RED),
        "total_probability",
    )
    assert got == Fraction(9, 20)

    # Enumerate the combined 20-outcome experiment directly: 10 balls in
    # each of 2 urns, chosen with equal prior, so all 20 (urn, ball) pairs
    # are equally likely.
    urn1 = ["red"] * D.URN_1_RED + ["blue"] * D.URN_1_BLUE
    urn2 = ["red"] * D.URN_2_RED + ["blue"] * D.URN_2_BLUE
    combined = [("urn1", b) for b in urn1] + [("urn2", b) for b in urn2]
    reds = [o for o in combined if o[1] == "red"]
    enumerated = Fraction(len(reds), len(combined))
    assert got == enumerated


# --------------------------------------------------------------------------
# Exercise 8 -- Monte Carlo error scaling
# --------------------------------------------------------------------------


def test_8_error_shrinks_across_four_decades_of_sample_size():
    target = float(D.MONTE_CARLO_TARGET)
    mean_errors = []
    for n in D.MONTE_CARLO_SAMPLE_SIZES:
        errors = []
        for seed in D.MONTE_CARLO_SEEDS:
            rng = np.random.default_rng(seed)
            got = attempt(lambda: S.simulate_sum_seven(rng, n), "simulate_sum_seven")
            errors.append(abs(got - target))
        mean_errors.append(sum(errors) / len(errors))
    assert mean_errors[-1] < mean_errors[0], (
        "the average error at the largest n must be smaller than at the "
        "smallest n"
    )
    assert mean_errors[-1] < mean_errors[0] / 5.0, (
        "a thousandfold increase in n should shrink the error by roughly "
        "sqrt(1000) = ~31x, well below a factor of 5"
    )


def test_8_the_shrink_looks_like_one_over_sqrt_n_not_one_over_n():
    target = float(D.MONTE_CARLO_TARGET)
    small_n, large_n = D.MONTE_CARLO_SAMPLE_SIZES[0], D.MONTE_CARLO_SAMPLE_SIZES[-1]
    ratio_n = large_n / small_n  # 1000

    def mean_error(n):
        errors = []
        for seed in D.MONTE_CARLO_SEEDS:
            rng = np.random.default_rng(seed)
            got = attempt(lambda: S.simulate_sum_seven(rng, n), "simulate_sum_seven")
            errors.append(abs(got - target))
        return sum(errors) / len(errors)

    error_ratio = mean_error(small_n) / mean_error(large_n)
    sqrt_prediction = ratio_n**0.5  # ~31.6
    linear_prediction = ratio_n  # 1000
    # The observed shrink must land far closer to the sqrt(n) prediction
    # than to the 1/n prediction -- checked as a wide band, not a point, so
    # the test is not flaky on a different machine.
    assert 0.3 * sqrt_prediction < error_ratio < 3.0 * sqrt_prediction
    assert error_ratio < linear_prediction / 10.0


# --------------------------------------------------------------------------
# Exercise 9 -- reproducibility
# --------------------------------------------------------------------------


def test_9_the_same_seed_gives_byte_identical_results():
    rng_a = np.random.default_rng(D.REPRODUCIBILITY_SEED_A)
    rng_b = np.random.default_rng(D.REPRODUCIBILITY_SEED_A)
    result_a = attempt(
        lambda: S.simulate_sum_seven(rng_a, D.REPRODUCIBILITY_TRIALS),
        "simulate_sum_seven",
    )
    result_b = attempt(
        lambda: S.simulate_sum_seven(rng_b, D.REPRODUCIBILITY_TRIALS),
        "simulate_sum_seven",
    )
    assert result_a == result_b, (
        "two Generators built from the same seed must produce identical "
        "results -- that is the entire point of default_rng(seed)"
    )


def test_9_a_different_seed_gives_a_different_but_still_close_result():
    rng_a = np.random.default_rng(D.REPRODUCIBILITY_SEED_A)
    rng_b = np.random.default_rng(D.REPRODUCIBILITY_SEED_B)
    result_a = attempt(
        lambda: S.simulate_sum_seven(rng_a, D.REPRODUCIBILITY_TRIALS),
        "simulate_sum_seven",
    )
    result_b = attempt(
        lambda: S.simulate_sum_seven(rng_b, D.REPRODUCIBILITY_TRIALS),
        "simulate_sum_seven",
    )
    assert result_a != result_b, "different seeds should not coincide exactly"
    target = float(D.MONTE_CARLO_TARGET)
    tol = 4.0 * D.standard_error(target, D.REPRODUCIBILITY_TRIALS)
    close(result_a, target, tol, "seed A")
    close(result_b, target, tol, "seed B")


# --------------------------------------------------------------------------
# The eighteen predictions
# --------------------------------------------------------------------------

EXPECTED: dict[str, object] = {
    "sample_space_size": 36,
    "p_sum_seven": float(Fraction(1, 6)),
    "addition_naive_sum": float(Fraction(1, 3)),
    "addition_true_union": float(Fraction(11, 36)),
    "addition_error_amount": float(Fraction(1, 36)),
    "de_mere_single_bet_probability": float(D.DE_MERE_SINGLE_EXACT),
    "de_mere_double_bet_probability": float(D.DE_MERE_DOUBLE_EXACT),
    "de_mere_favorable_bet": 1,
    "independent_pair_holds": True,
    "dependent_pair_holds": False,
    "mutually_exclusive_implies_dependent": True,
    "conditional_p_sum8_given_first_even": float(Fraction(1, 6)),
    "conditional_formula_matches_filter": True,
    "urn_total_probability_red": 0.45,
    "urn_enumeration_matches_formula": True,
    "monte_carlo_error_shrinks_with_n": True,
    "monte_carlo_error_ratio_near_sqrt10": True,
    "reproducibility_same_seed_identical": True,
    "reproducibility_different_seed_differs": True,
}

HINTS: dict[str, str] = {
    "addition_error_amount": (
        "This is exactly P(A and B), the region the naive sum counted twice."
    ),
    "de_mere_favorable_bet": (
        "A bet is favourable to the player when its probability is above "
        "0.5. Only one of the two is."
    ),
    "mutually_exclusive_implies_dependent": (
        "P(A | B) collapses to 0 for a mutually exclusive pair, but P(A) is "
        "not 0 -- so knowing B happened changed what you believe about A."
    ),
    "monte_carlo_error_ratio_near_sqrt10": (
        "Multiplying n by 10 should shrink the error by about sqrt(10), "
        "roughly 3.16x -- not by 10x."
    ),
}


@pytest.mark.parametrize("key", sorted(EXPECTED))
def test_predictions(key):
    got = need(answers.ANSWERS.get(key), f"answers.ANSWERS[{key!r}]")
    want = EXPECTED[key]
    hint = HINTS.get(key, "")
    if isinstance(want, bool) or isinstance(want, int):
        assert got == want, f"{key}: your answer {got!r}, expected {want!r}. {hint}"
    else:
        assert abs(float(got) - want) < 1e-9, (
            f"{key}: your answer {got!r}, expected {want!r}. {hint}"
        )


def test_every_answer_key_is_still_present():
    missing = sorted(set(EXPECTED) - set(answers.ANSWERS))
    assert not missing, f"answers.py is missing these keys: {missing}"

#!/usr/bin/env bash
# Tests for the Day 113 lab. Run from the lab directory:
#   bash tests/run_tests.sh
#
# The harness proves the lesson's claims by running code and reading real
# values, never by reading source:
#
#   * P(two dice sum to 7) is exactly Fraction(1, 6), by enumeration;
#   * the naive addition-rule sum overstates the true union by exactly
#     P(A and B), and the exact size of that error is checked;
#   * de Méré's two bets -- 1 - (5/6)^4 and 1 - (35/36)^24 -- are computed
#     exactly and both round to the historical figures, 0.5177 and 0.4914,
#     and both are confirmed by simulation within three standard errors;
#   * a genuinely independent pair of dice events satisfies P(A and B) ==
#     P(A) x P(B) exactly, and a genuinely dependent pair does not;
#   * a mutually exclusive pair with non-zero probabilities has P(A | B) ==
#     0 while P(A) != 0, proving mutual exclusivity implies dependence;
#   * conditioning by formula and conditioning by filtering the sample space
#     agree exactly;
#   * the law of total probability over two urns matches a direct
#     enumeration of the combined experiment;
#   * Monte Carlo error shrinks like 1/sqrt(n), asserted as a trend averaged
#     over twenty seeds, never a single sampled value;
#   * the same seed gives byte-identical simulated results, and a different
#     seed gives a different but still-close one;
#   * nothing is left behind on disk.
#
# Everything after the one-time install runs offline. Nothing binds a port,
# nothing writes outside the lab, nothing needs a key. Deterministic,
# non-interactive, exits 0 only if every check passes.
set -u

export PYTHONDONTWRITEBYTECODE=1

lab_dir="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)"

# Bytecode left by an EARLIER command is not this run's litter. The README
# documents `pytest starter -q`, and running it writes .pyc files that would
# then fail the cleanliness check at the end of this script -- failing the
# reader for following the instructions. Clearing them here makes that final
# check measure what it claims to: what THIS run left behind. `.venv` is
# untouched, because the packages' own bytecode is theirs, not ours.
find "${lab_dir}" -name '.venv' -prune -o -type d -name '__pycache__' -exec rm -rf {} + 2>/dev/null || true
find "${lab_dir}" -name '.venv' -prune -o -type d -name '.pytest_cache' -exec rm -rf {} + 2>/dev/null || true

failures=0
checks=0

check() {
  local label="$1" ok="$2"
  checks=$((checks + 1))
  if [ "${ok}" = "yes" ]; then
    echo "  ok: ${label}"
  else
    echo "  FAIL: ${label}"
    failures=$((failures + 1))
  fi
}

check_eq() {
  # check_eq <label> <expected> <actual>
  if [ "$2" = "$3" ]; then
    check "$1" "yes"
  else
    check "$1 (expected [$2], got [$3])" "no"
  fi
}

# Resolve pytest: an explicit override, then this lab's .venv, then PATH.
# Fails loudly with instructions rather than silently skipping checks.
resolve_tool() {
  local tool="$1" override="$2"
  if [ -n "${override}" ] && [ -x "${override}" ]; then echo "${override}"; return 0; fi
  if [ -x "${lab_dir}/.venv/bin/${tool}" ]; then echo "${lab_dir}/.venv/bin/${tool}"; return 0; fi
  if command -v "${tool}" >/dev/null 2>&1; then command -v "${tool}"; return 0; fi
  return 1
}

pytest_bin="$(resolve_tool pytest "${PYTEST:-}")" || {
  echo "FAIL: pytest not found." >&2
  echo "  Install the lab's dependencies with:" >&2
  echo "    python3 -m venv .venv" >&2
  echo "    .venv/bin/pip install -r requirements/requirements.txt" >&2
  echo "  Or point this suite at an existing pytest:" >&2
  echo "    PYTEST=/path/to/pytest bash tests/run_tests.sh" >&2
  exit 1
}

python_bin="$(dirname "${pytest_bin}")/python3"
if [ ! -x "${python_bin}" ]; then
  python_bin="$(command -v python3 || true)"
fi
if [ -z "${python_bin}" ]; then
  echo "FAIL: python3 not found on PATH." >&2
  exit 1
fi

if ! "${python_bin}" -c "import numpy" >/dev/null 2>&1; then
  echo "FAIL: numpy is not importable from ${python_bin}." >&2
  echo "  Install the lab's dependencies with:" >&2
  echo "    python3 -m venv .venv" >&2
  echo "    .venv/bin/pip install -r requirements/requirements.txt" >&2
  exit 1
fi

echo "Day 113 — Probability You Can Count"
echo

# --------------------------------------------------------------------------
echo "1. The tools and the versions this lab was written against"
# --------------------------------------------------------------------------

versions="$("${python_bin}" - <<'PY'
import platform
import sys
from importlib.metadata import version

print(f"python   {platform.python_version()}")
for name in ("numpy", "pytest"):
    print(f"{name:<8} {version(name)}")
print(f"platform {platform.platform()}")
print(f"exe      {sys.executable.rsplit('/', 3)[-1]}")
PY
)"
echo "${versions}" | sed 's/^/  /'

pinned_numpy="$(grep -E '^numpy==' "${lab_dir}/requirements/requirements.txt" | cut -d= -f3)"
installed_numpy="$("${python_bin}" -c "from importlib.metadata import version; print(version('numpy'))")"
check_eq "installed numpy matches requirements.txt" "${pinned_numpy}" "${installed_numpy}"

major="$("${python_bin}" -c "import numpy; print(numpy.__version__.split('.')[0])")"
check_eq "numpy is version 2 or later" "2" "${major}"

# --------------------------------------------------------------------------
echo
echo "2. Every reference script runs and every assertion inside it holds"
# --------------------------------------------------------------------------

for script in 01_sample_space_and_events 02_addition_rule 03_de_mere \
              04_independence_vs_dependence 05_mutual_exclusivity_implies_dependence \
              06_conditioning_by_restriction 07_law_of_total_probability \
              08_monte_carlo_error_scaling 09_reproducibility; do
  out="$(cd "${lab_dir}/examples" && "${python_bin}" "${script}.py" 2>&1)"
  status=$?
  if [ "${status}" -ne 0 ]; then
    check "${script}.py exits 0" "no"
    echo "${out}" | tail -5 | sed 's/^/      /'
  else
    check "${script}.py exits 0" "yes"
  fi
  case "${out}" in
    *"${script}.py: every assertion held."*)
      check "${script}.py reports every assertion held" "yes" ;;
    *) check "${script}.py reports every assertion held" "no" ;;
  esac
done

# --------------------------------------------------------------------------
echo
echo "3. The reference pytest suite: real values, real exceptions"
# --------------------------------------------------------------------------

ref_out="$(cd "${lab_dir}" && "${pytest_bin}" examples -q -p no:cacheprovider 2>&1)"
ref_status=$?
echo "${ref_out}" | tail -3 | sed 's/^/  /'
if [ "${ref_status}" -eq 0 ]; then
  check "pytest examples exits 0" "yes"
else
  check "pytest examples exits 0" "no"
fi
case "${ref_out}" in
  *" failed"*) check "no test in the reference suite failed" "no" ;;
  *)           check "no test in the reference suite failed" "yes" ;;
esac
ref_passed="$(printf '%s\n' "${ref_out}" | grep -o '[0-9][0-9]* passed' | head -1 | cut -d' ' -f1)"
if [ "${ref_passed:-0}" -ge 80 ]; then
  check "the reference suite ran at least 80 tests (ran ${ref_passed})" "yes"
else
  check "the reference suite ran at least 80 tests (ran ${ref_passed:-0})" "no"
fi

# --------------------------------------------------------------------------
echo
echo "4. The starter suite skips unattempted work instead of failing it"
# --------------------------------------------------------------------------

start_out="$(cd "${lab_dir}" && "${pytest_bin}" starter -q -p no:cacheprovider 2>&1)"
start_status=$?
echo "${start_out}" | tail -3 | sed 's/^/  /'
if [ "${start_status}" -eq 0 ]; then
  check "pytest starter exits 0 on an untouched checkout" "yes"
else
  check "pytest starter exits 0 on an untouched checkout" "no"
fi
case "${start_out}" in
  *" failed"*) check "the starter suite reports no failures" "no" ;;
  *)           check "the starter suite reports no failures" "yes" ;;
esac
case "${start_out}" in
  *skipped*) check "unwritten exercises are reported as skipped, not passed" "yes" ;;
  *) check "unwritten exercises are reported as skipped, not passed" "no" ;;
esac

# The import guard. Both directories contain modules called `probability`,
# `simulate`, `dataset` and `answers`, and pytest imports test files by
# putting their directory on sys.path -- so collecting both suites at once
# would otherwise let the starter tests import the REFERENCE solution and
# report unwritten exercises as passing. Each directory's conftest.py
# prevents that. This check proves it still does: across both suites, the
# skip count must be unchanged.
both_out="$(cd "${lab_dir}" && "${pytest_bin}" -q -p no:cacheprovider 2>&1)"
start_skipped="$(printf '%s\n' "${start_out}" | grep -o '[0-9][0-9]* skipped' | head -1 | cut -d' ' -f1)"
both_skipped="$(printf '%s\n' "${both_out}" | grep -o '[0-9][0-9]* skipped' | head -1 | cut -d' ' -f1)"
check_eq "collecting both suites at once does not turn skips into passes" \
  "${start_skipped:-none}" "${both_skipped:-none}"

# --------------------------------------------------------------------------
echo
echo "5. The lesson's claims, checked one value at a time"
# --------------------------------------------------------------------------

facts="$(cd "${lab_dir}/examples" && "${python_bin}" - <<'PY'
from fractions import Fraction

import numpy as np

import dataset as D
import probability as P
import simulate as S

space = P.sample_space_two_dice()
print("space_size", len(space))

a = P.event(space, D.ADDITION_EVENT_A)
b = P.event(space, D.ADDITION_EVENT_B)
p_a, p_b = P.probability(a, space), P.probability(b, space)
p_ab = P.probability(a & b, space)
naive = P.naive_sum(p_a, p_b)
true_union = P.addition_rule(p_a, p_b, p_ab)
print("addition_naive", naive)
print("addition_true", true_union)
print("addition_error", naive - true_union)
print("addition_error_equals_intersection", (naive - true_union) == p_ab)

single = P.at_least_one(Fraction(1, 6), D.DE_MERE_SINGLE_ROLLS)
double = P.at_least_one(Fraction(1, 36), D.DE_MERE_DOUBLE_ROLLS)
print("de_mere_single_rounded", round(float(single), 4))
print("de_mere_double_rounded", round(float(double), 4))
print("de_mere_single_favourable", single > Fraction(1, 2))
print("de_mere_double_favourable", double > Fraction(1, 2))

rng = np.random.default_rng(D.REPRODUCIBILITY_SEED_A)
sim_single = S.simulate_at_least_one_six(rng, D.DE_MERE_SIM_TRIALS)
sim_double = S.simulate_at_least_one_double_six(rng, D.DE_MERE_SIM_TRIALS)
print("de_mere_single_sim_within_tol", abs(sim_single - float(single)) < D.DE_MERE_SINGLE_TOL)
print("de_mere_double_sim_within_tol", abs(sim_double - float(double)) < D.DE_MERE_DOUBLE_TOL)

ip_a, ip_b = D.INDEPENDENT_PAIR
ia, ib = P.event(space, ip_a), P.event(space, ip_b)
p_ia, p_ib = P.probability(ia, space), P.probability(ib, space)
p_iab = P.probability(ia & ib, space)
print("independent_holds", P.is_independent(p_ia, p_ib, p_iab))

dp_a, dp_b = D.DEPENDENT_PAIR
da, db = P.event(space, dp_a), P.event(space, dp_b)
p_da, p_db = P.probability(da, space), P.probability(db, space)
p_dab = P.probability(da & db, space)
print("dependent_holds", P.is_independent(p_da, p_db, p_dab))

me_a, me_b = D.MUTUALLY_EXCLUSIVE_PAIR
ma, mb = P.event(space, me_a), P.event(space, me_b)
p_ma = P.probability(ma, space)
p_mab = P.probability(ma & mb, space)
p_mb = P.probability(mb, space)
p_ma_given_mb = P.conditional(p_mab, p_mb)
print("mutually_exclusive_conditional_is_zero", p_ma_given_mb == 0)
print("mutually_exclusive_implies_dependent", p_ma_given_mb != p_ma)

ca = P.event(space, D.CONDITIONING_EVENT_A)
cb = P.event(space, D.CONDITIONING_EVENT_B)
p_ca_cb = P.probability(ca & cb, space)
p_cb = P.probability(cb, space)
by_formula = P.conditional(p_ca_cb, p_cb)
by_filtering = P.probability(P.event(cb, D.CONDITIONING_EVENT_A), cb)
print("conditioning_agrees", by_formula == by_filtering == Fraction(1, 6))

total = P.total_probability(D.URN_PRIOR, D.URN_CONDITIONAL_RED)
urn1 = ["red"] * D.URN_1_RED + ["blue"] * D.URN_1_BLUE
urn2 = ["red"] * D.URN_2_RED + ["blue"] * D.URN_2_BLUE
combined = [("u1", x) for x in urn1] + [("u2", x) for x in urn2]
reds = [o for o in combined if o[1] == "red"]
enumerated = Fraction(len(reds), len(combined))
print("urn_total", total)
print("urn_matches_enumeration", total == enumerated)

target = float(D.MONTE_CARLO_TARGET)
means = []
for n in D.MONTE_CARLO_SAMPLE_SIZES:
    errors = [
        abs(S.simulate_sum_seven(np.random.default_rng(seed), n) - target)
        for seed in D.MONTE_CARLO_SEEDS
    ]
    means.append(sum(errors) / len(errors))
print("mc_error_monotone_decreasing", all(means[i + 1] < means[i] for i in range(len(means) - 1)))
n_ratio = D.MONTE_CARLO_SAMPLE_SIZES[-1] / D.MONTE_CARLO_SAMPLE_SIZES[0]
error_ratio = means[0] / means[-1]
sqrt_pred = n_ratio ** 0.5
print("mc_error_ratio_near_sqrt_n", 0.3 * sqrt_pred < error_ratio < 3.0 * sqrt_pred)
print("mc_error_ratio_far_from_linear", error_ratio < n_ratio / 10.0)

a1 = S.simulate_sum_seven(np.random.default_rng(D.REPRODUCIBILITY_SEED_A), D.REPRODUCIBILITY_TRIALS)
a2 = S.simulate_sum_seven(np.random.default_rng(D.REPRODUCIBILITY_SEED_A), D.REPRODUCIBILITY_TRIALS)
b1 = S.simulate_sum_seven(np.random.default_rng(D.REPRODUCIBILITY_SEED_B), D.REPRODUCIBILITY_TRIALS)
print("same_seed_identical", a1 == a2)
print("different_seed_differs", a1 != b1)
tol = 4.0 * D.standard_error(target, D.REPRODUCIBILITY_TRIALS)
print("both_seeds_within_tolerance", abs(a1 - target) < tol and abs(b1 - target) < tol)
PY
)"

get() { printf '%s\n' "${facts}" | grep "^$1 " | cut -d' ' -f2-; }

check_eq "the sample space has 36 outcomes" "36" "$(get space_size)"
check_eq "the naive addition-rule sum is 1/3" "1/3" "$(get addition_naive)"
check_eq "the true union is 11/36" "11/36" "$(get addition_true)"
check_eq "the naive sum's error is exactly 1/36" "1/36" "$(get addition_error)"
check_eq "and that error equals P(A and B) exactly" "True" "$(get addition_error_equals_intersection)"
check_eq "de Méré's single-die bet rounds to 0.5177" "0.5177" "$(get de_mere_single_rounded)"
check_eq "de Méré's double-dice bet rounds to 0.4914" "0.4914" "$(get de_mere_double_rounded)"
check_eq "the single-die bet favours the player" "True" "$(get de_mere_single_favourable)"
check_eq "the double-dice bet does NOT favour the player" "False" "$(get de_mere_double_favourable)"
check_eq "the single-die simulation lands within 3 standard errors" "True" "$(get de_mere_single_sim_within_tol)"
check_eq "the double-dice simulation lands within 3 standard errors" "True" "$(get de_mere_double_sim_within_tol)"
check_eq "the independent pair satisfies P(A and B) == P(A) x P(B)" "True" "$(get independent_holds)"
check_eq "the dependent pair does NOT satisfy it" "False" "$(get dependent_holds)"
check_eq "a mutually exclusive pair has P(A | B) exactly 0" "True" "$(get mutually_exclusive_conditional_is_zero)"
check_eq "so mutual exclusivity implies dependence" "True" "$(get mutually_exclusive_implies_dependent)"
check_eq "conditioning by formula and by filtering agree exactly at 1/6" "True" "$(get conditioning_agrees)"
check_eq "the urns' weighted total probability of red is 9/20" "9/20" "$(get urn_total)"
check_eq "and it matches the direct enumeration of the combined experiment" "True" "$(get urn_matches_enumeration)"
check_eq "Monte Carlo error falls monotonically across four decades of n" "True" "$(get mc_error_monotone_decreasing)"
check_eq "the observed shrink is close to the sqrt(n) prediction" "True" "$(get mc_error_ratio_near_sqrt_n)"
check_eq "and far from the false 1/n prediction" "True" "$(get mc_error_ratio_far_from_linear)"
check_eq "the same seed gives byte-identical simulated results" "True" "$(get same_seed_identical)"
check_eq "a different seed gives a different result" "True" "$(get different_seed_differs)"
check_eq "and both seeds still land within tolerance of the truth" "True" "$(get both_seeds_within_tolerance)"

# --------------------------------------------------------------------------
echo
echo "6. The harness can actually fail"
# --------------------------------------------------------------------------

# A green test suite proves nothing until you have watched it go red. This
# section re-runs the reference pytest suite with one assertion deliberately
# broken -- the addition-rule error amount replaced with the wrong value --
# and asserts that the run reports the failure and exits non-zero.
sentinel_file="${lab_dir}/examples/test_reference.py"
if grep -q "def test_naive_sum_overstates_by_exactly_the_intersection" "${sentinel_file}"; then
  backup="$(mktemp)"
  cp "${sentinel_file}" "${backup}"
  python3 - "${sentinel_file}" <<'PY'
import sys
path = sys.argv[1]
text = open(path).read()
needle = "assert naive - true_union == p_ab"
replacement = "assert naive - true_union == p_ab * 2  # DELIBERATELY WRONG: self-test only"
assert needle in text, "sentinel assertion not found"
open(path, "w").write(text.replace(needle, replacement, 1))
PY
  self_out="$(cd "${lab_dir}" && "${pytest_bin}" examples -q -p no:cacheprovider 2>&1)"
  self_status=$?
  cp "${backup}" "${sentinel_file}"
  rm -f "${backup}"
  find "${lab_dir}" -name '.venv' -prune -o -type d -name '__pycache__' -exec rm -rf {} + 2>/dev/null || true

  if [ "${self_status}" -ne 0 ]; then
    check "a deliberately broken assertion makes the suite exit non-zero (${self_status})" "yes"
  else
    check "a deliberately broken assertion makes the suite exit non-zero" "no"
  fi
  case "${self_out}" in
    *"test_naive_sum_overstates_by_exactly_the_intersection"*"failed"*|*"FAILED"*"test_naive_sum_overstates_by_exactly_the_intersection"*)
      check "the failing test is named in the output" "yes" ;;
    *) check "the failing test is named in the output" "no" ;;
  esac
  case "${self_out}" in
    *"1 failed"*) check "exactly one test failed" "yes" ;;
    *) check "exactly one test failed" "no" ;;
  esac
else
  check "sentinel assertion found in examples/test_reference.py" "no"
fi

# --------------------------------------------------------------------------
echo
echo "7. Nothing was left behind"
# --------------------------------------------------------------------------

# `.venv` is pruned from both searches below. The virtual environment ships
# NumPy's and pytest's own precompiled bytecode -- hundreds of __pycache__
# directories that came with the packages and have nothing to do with
# whether THIS lab tidied up after itself.

if find "${lab_dir}" -name '.venv' -prune -o -type d -name '__pycache__' -print -quit 2>/dev/null | grep -q .; then
  check "no __pycache__ directory left by the lab's own code" "no"
else
  check "no __pycache__ directory left by the lab's own code" "yes"
fi

if find "${lab_dir}" -name '.venv' -prune -o -type d -name '.pytest_cache' -print -quit 2>/dev/null | grep -q .; then
  check "no .pytest_cache directory left under the lab" "no"
else
  check "no .pytest_cache directory left under the lab" "yes"
fi

if grep -rqE 'urlopen|requests\.|socket\.|http://|https://' \
     "${lab_dir}/examples" "${lab_dir}/starter" 2>/dev/null; then
  check "no lab source opens a network connection" "no"
else
  check "no lab source opens a network connection" "yes"
fi

echo
echo "${checks} checks, ${failures} failure(s)."
[ "${failures}" -eq 0 ]