Hands-on lab — Day 130: Distributions and Relationships

"""Shared fixtures and headless matplotlib setup.

pytest finds this file by itself -- nothing imports it. `matplotlib.use`
must run before `pyplot` is imported anywhere, which is why it happens
here, first. Every test that opens a Figure is responsible for its own
`plt.close()`; `_close_all_figures` below is a backstop, not a substitute.
"""

import matplotlib

matplotlib.use("Agg")

import matplotlib.pyplot as plt
import pytest

from data import (
    bimodal_for_binning,
    discrete_for_jitter,
    matched_quartile_pair,
    normal_for_ecdf,
    overplotted_cloud,
    positive_for_kde_boundary,
    quadratic_relationship,
    skewed_for_bin_rules,
    target_five_number_summary,
)


@pytest.fixture
def bimodal_sample():
    return bimodal_for_binning()


@pytest.fixture
def skewed_sample():
    return skewed_for_bin_rules()


@pytest.fixture
def positive_sample():
    return positive_for_kde_boundary()


@pytest.fixture
def quartile_pair():
    return matched_quartile_pair()


@pytest.fixture
def quartile_targets():
    return target_five_number_summary()


@pytest.fixture
def ecdf_sample():
    return normal_for_ecdf()


@pytest.fixture
def overplot_cloud():
    return overplotted_cloud()


@pytest.fixture
def quadratic_data():
    return quadratic_relationship()


@pytest.fixture
def discrete_sample():
    return discrete_for_jitter()


@pytest.fixture(autouse=True)
def _close_all_figures():
    yield
    plt.close("all")

"""The samples every exercise in this lab is built from.

Every table here is generated from a fixed `numpy.random.default_rng` seed,
never loaded from a file and never re-seeded per test -- so the same call
produces the same numbers on any machine with the same NumPy version, and a
reader can reproduce every asserted value by running the function directly.

Two constructed pairs carry the lab's two centrepiece demonstrations:

`bimodal_for_binning()` -- exercise 1 and exercise 3. Two normal clusters
close enough together that a coarse 5-bin histogram merges them into one
hump, a well-chosen bin count (Freedman-Diaconis lands on 13 bins here)
recovers two, and 100 bins turns the same 500 points into visual noise.
The same sample is reused for the KDE bandwidth demonstration in
exercise 3, because the story is identical: the bin width and the KDE
bandwidth are the same kind of decision.

`matched_quartile_pair()` -- exercise 5, the day's centrepiece. A bimodal
sample and a unimodal sample, built from two different piecewise-linear
quantile functions, engineered so their five-number summaries agree to
within 0.3 units on identical control points (min, Q1, median, Q3, max
all hand-picked as 10 / 28 / 40 / 52 / 70) while their interior shape is
completely different. A boxplot of either looks identical to the other;
a histogram does not.
"""

from __future__ import annotations

import numpy as np

# --------------------------------------------------------------------------
# bimodal_for_binning -- exercises 1 and 3.
#
# Two normal clusters, means 40 and 54 (gap 14), sd 8 each, 250 points per
# cluster. Close enough that 5 wide bins wash the gap out into one hump;
# Freedman-Diaconis (13 bins on this draw) recovers two; 100 bins gives
# nothing but sampling noise per bin.
# --------------------------------------------------------------------------


def bimodal_for_binning() -> np.ndarray:
    rng = np.random.default_rng(42)
    low_cluster = rng.normal(40, 8, 250)
    high_cluster = rng.normal(54, 8, 250)
    return np.concatenate([low_cluster, high_cluster])


# --------------------------------------------------------------------------
# skewed_for_bin_rules -- exercise 2. A right-skewed, strictly positive
# sample (log-normal) on which Sturges, Scott and Freedman-Diaconis
# genuinely disagree about the bin count.
# --------------------------------------------------------------------------


def skewed_for_bin_rules() -> np.ndarray:
    rng = np.random.default_rng(123)
    return rng.lognormal(mean=3.0, sigma=0.6, size=400)


# --------------------------------------------------------------------------
# positive_for_kde_boundary -- exercise 4. Strictly positive (an
# exponential), so any density a KDE places below zero is visibly wrong.
# --------------------------------------------------------------------------


def positive_for_kde_boundary() -> np.ndarray:
    rng = np.random.default_rng(9)
    return rng.exponential(scale=5.0, size=400)


# --------------------------------------------------------------------------
# matched_quartile_pair -- exercise 5, the boxplot's blind spot.
#
# Both samples are built from a piecewise-linear quantile function: a
# function from rank (0 to 1) to value, evaluated at 240 evenly spaced
# ranks. Both functions are pinned to pass through the same five control
# points -- (0, 10), (0.25, 28), (0.5, 40), (0.75, 52), (1.0, 70) -- so
# both samples' five-number summaries land within about 0.3 units of
# those targets. Between the control points the two functions diverge
# completely: `_unimodal_quantile` uses a smooth convex curve on each
# side of the median (density highest at the centre, lowest at the
# tails -- one hump); `_bimodal_quantile` inserts extra control points
# that pack a large fraction of the probability mass into two narrow
# bands just above Q1 and just below Q3, with a sparse valley at the
# median and sparse far tails (two humps, a dip in the middle).
# --------------------------------------------------------------------------

_TARGET_MIN, _TARGET_Q1, _TARGET_MED, _TARGET_Q3, _TARGET_MAX = (
    10.0,
    28.0,
    40.0,
    52.0,
    70.0,
)


def _unimodal_quantile(ranks: np.ndarray) -> np.ndarray:
    med = _TARGET_MED
    d3, dmax = _TARGET_Q3 - med, _TARGET_MAX - med
    d1, dmin = med - _TARGET_Q1, med - _TARGET_MIN
    coeff_matrix = np.array([[0.5, 0.25], [1.0, 1.0]])
    a_right, b_right = np.linalg.solve(coeff_matrix, np.array([d3, dmax]))
    a_left, b_left = np.linalg.solve(coeff_matrix, np.array([d1, dmin]))

    out = np.empty_like(ranks)
    on_right = ranks >= 0.5
    u_right = (ranks[on_right] - 0.5) / 0.5
    out[on_right] = med + a_right * u_right + b_right * u_right**2
    u_left = (0.5 - ranks[~on_right]) / 0.5
    out[~on_right] = med - (a_left * u_left + b_left * u_left**2)
    return out


def _bimodal_quantile(ranks: np.ndarray) -> np.ndarray:
    control_ranks = [0.0, 0.15, 0.25, 0.40, 0.50, 0.60, 0.75, 0.85, 1.00]
    control_values = [
        _TARGET_MIN,
        25.0,
        _TARGET_Q1,
        31.0,
        _TARGET_MED,
        49.0,
        _TARGET_Q3,
        55.0,
        _TARGET_MAX,
    ]
    return np.interp(ranks, control_ranks, control_values)


def matched_quartile_pair(n: int = 240) -> tuple[np.ndarray, np.ndarray]:
    ranks = (np.arange(n) + 0.5) / n
    unimodal = _unimodal_quantile(ranks)
    bimodal = _bimodal_quantile(ranks)
    return bimodal, unimodal


def target_five_number_summary() -> np.ndarray:
    """The five control-point values both samples in the matched pair were built to hit."""
    return np.array(
        [_TARGET_MIN, _TARGET_Q1, _TARGET_MED, _TARGET_Q3, _TARGET_MAX]
    )


# --------------------------------------------------------------------------
# normal_for_ecdf -- exercise 6. Odd sample size on purpose: with an odd n
# the median is a single real observation rather than an average of two,
# so the ECDF step and `numpy.median` can be compared for an exact match.
# --------------------------------------------------------------------------


def normal_for_ecdf() -> np.ndarray:
    rng = np.random.default_rng(3)
    return rng.normal(0, 1, 301)


# --------------------------------------------------------------------------
# overplotted_cloud -- exercise 7. 20,000 points from a standard normal in
# both dimensions, deliberately plotted small (3x3 inches at 72 dpi) so a
# meaningful fraction of them land on the very same screen pixel.
# --------------------------------------------------------------------------


def overplotted_cloud(n: int = 20_000) -> tuple[np.ndarray, np.ndarray]:
    rng = np.random.default_rng(11)
    x = rng.normal(0, 1, n)
    y = rng.normal(0, 1, n)
    return x, y


# --------------------------------------------------------------------------
# quadratic_relationship -- exercise 8. x symmetric around 0, y = x^2 plus
# noise: a strong, deterministic relationship with almost no LINEAR
# component (Pearson) and, because the parabola is symmetric, almost no
# MONOTONIC component either (Spearman) -- only fitting or plotting the
# actual shape reveals it. See the lesson and FIELDS.md for why this is a
# sharper example than one where Spearman alone would have caught it.
# --------------------------------------------------------------------------


def quadratic_relationship(n: int = 300) -> tuple[np.ndarray, np.ndarray]:
    rng = np.random.default_rng(5)
    x = rng.uniform(-10, 10, n)
    y = x**2 + rng.normal(0, 3, n)
    return x, y


# --------------------------------------------------------------------------
# discrete_for_jitter -- exercise 9. Integers 1..5, as a five-point Likert
# scale might arrive.
# --------------------------------------------------------------------------


def discrete_for_jitter(n: int = 200) -> np.ndarray:
    rng = np.random.default_rng(17)
    return rng.integers(1, 6, n)

"""Reference solutions -- Day 130, Pictures of a Distribution.

Nine exercises, each proving a claim from the lesson by running real
matplotlib / seaborn / pandas / NumPy code and asserting on the numbers
and artist state that code actually produces -- never on image bytes.

Run with: pytest examples -q
"""

from __future__ import annotations

import numpy as np
import pandas as pd
import seaborn as sns


def count_local_maxima(counts) -> int:
    """A bin is a local maximum if it is strictly taller than its only
    neighbour (an edge bin) or both neighbours (an interior bin)."""
    counts = list(counts)
    n = len(counts)
    maxima = 0
    for i in range(n):
        left = counts[i - 1] if i > 0 else None
        right = counts[i + 1] if i < n - 1 else None
        if left is None and (right is None or counts[i] > right):
            maxima += 1
        elif right is None and counts[i] > left:
            maxima += 1
        elif left is not None and right is not None and counts[i] > left and counts[i] > right:
            maxima += 1
    return maxima


def count_curve_modes(y) -> int:
    """Local maxima of a continuous curve (a KDE line), interior points only."""
    y = list(y)
    return sum(
        1 for i in range(1, len(y) - 1) if y[i] > y[i - 1] and y[i] > y[i + 1]
    )


# --------------------------------------------------------------------------
# Exercise 1 -- bin width changes the story.
# --------------------------------------------------------------------------


def test_01_bin_width_changes_the_story(bimodal_sample):
    counts_5, _ = np.histogram(bimodal_sample, bins=5)
    counts_100, _ = np.histogram(bimodal_sample, bins=100)

    modes_5 = count_local_maxima(counts_5)
    modes_100 = count_local_maxima(counts_100)
    assert modes_5 != modes_100

    # 5 wide bins wash the two clusters into a single hump.
    assert modes_5 == 1
    # 100 narrow bins turn 500 points into visual noise -- many spurious
    # local maxima, none of them the real two-cluster structure.
    assert modes_100 > 10

    fd_edges = np.histogram_bin_edges(bimodal_sample, bins="fd")
    fd_counts, _ = np.histogram(bimodal_sample, bins=fd_edges)
    fd_bin_count = len(fd_edges) - 1
    fd_modes = count_local_maxima(fd_counts)

    # Freedman-Diaconis recovers the real two-mode structure.
    assert fd_modes == 2
    print(
        f"  bins=5 -> {modes_5} mode(s); bins=100 -> {modes_100} mode(s); "
        f"Freedman-Diaconis -> {fd_bin_count} bins, {fd_modes} mode(s)"
    )


# --------------------------------------------------------------------------
# Exercise 2 -- the three rules disagree.
# --------------------------------------------------------------------------


def test_02_the_three_rules_disagree(skewed_sample):
    sturges_edges = np.histogram_bin_edges(skewed_sample, bins="sturges")
    scott_edges = np.histogram_bin_edges(skewed_sample, bins="scott")
    fd_edges = np.histogram_bin_edges(skewed_sample, bins="fd")

    sturges_n = len(sturges_edges) - 1
    scott_n = len(scott_edges) - 1
    fd_n = len(fd_edges) - 1

    bin_counts = {sturges_n, scott_n, fd_n}
    assert len(bin_counts) == 3, "all three rules must disagree on this skewed sample"

    print(f"  sturges={sturges_n} scott={scott_n} fd={fd_n}")


# --------------------------------------------------------------------------
# Exercise 3 -- KDE bandwidth.
# --------------------------------------------------------------------------


def test_03_kde_bandwidth(bimodal_sample):
    import matplotlib.pyplot as plt

    fig, ax = plt.subplots()
    sns.kdeplot(bimodal_sample, ax=ax, bw_adjust=1.0)
    _, y_default = ax.lines[0].get_data()
    plt.close(fig)

    fig, ax = plt.subplots()
    sns.kdeplot(bimodal_sample, ax=ax, bw_adjust=3.0)
    _, y_smoothed = ax.lines[0].get_data()
    plt.close(fig)

    modes_default = count_curve_modes(y_default)
    modes_smoothed = count_curve_modes(y_smoothed)

    assert modes_default == 2
    assert modes_smoothed == 1
    print(f"  bw_adjust=1.0 -> {modes_default} mode(s); bw_adjust=3.0 -> {modes_smoothed} mode(s)")


# --------------------------------------------------------------------------
# Exercise 4 -- the KDE boundary problem.
# --------------------------------------------------------------------------


def test_04_kde_boundary_problem(positive_sample):
    import matplotlib.pyplot as plt

    assert positive_sample.min() > 0

    fig, ax = plt.subplots()
    sns.kdeplot(positive_sample, ax=ax)
    x, y = ax.lines[0].get_data()
    plt.close(fig)

    below_zero = x < 0
    assert below_zero.any(), "the default KDE grid must extend below zero"

    mass_below_zero = np.trapezoid(y[below_zero], x[below_zero])
    total_mass = np.trapezoid(y, x)
    fraction_below_zero = mass_below_zero / total_mass

    assert fraction_below_zero > 0.03
    print(f"  fraction of KDE mass below zero: {fraction_below_zero:.4f}")


# --------------------------------------------------------------------------
# Exercise 5 -- the boxplot's blind spot. The day's centrepiece.
# --------------------------------------------------------------------------


def test_05_boxplot_blind_spot(quartile_pair, quartile_targets):
    bimodal, unimodal = quartile_pair

    def five_number_summary(x):
        return np.percentile(x, [0, 25, 50, 75, 100])

    bimodal_summary = five_number_summary(bimodal)
    unimodal_summary = five_number_summary(unimodal)

    tolerance = 0.3
    assert np.max(np.abs(bimodal_summary - quartile_targets)) < tolerance
    assert np.max(np.abs(unimodal_summary - quartile_targets)) < tolerance
    # The two summaries agree with EACH OTHER to the same tight tolerance --
    # this is what makes their boxplots indistinguishable.
    assert np.max(np.abs(bimodal_summary - unimodal_summary)) < tolerance

    bin_count = 15
    bimodal_counts, _ = np.histogram(bimodal, bins=bin_count)
    unimodal_counts, _ = np.histogram(unimodal, bins=bin_count)
    bimodal_modes = count_local_maxima(bimodal_counts)
    unimodal_modes = count_local_maxima(unimodal_counts)

    assert bimodal_modes == 2
    assert unimodal_modes == 1
    assert bimodal_modes != unimodal_modes

    print(
        f"  bimodal 5-num: {np.round(bimodal_summary, 2)} ({bimodal_modes} modes at "
        f"{bin_count} bins)\n"
        f"  unimodal 5-num: {np.round(unimodal_summary, 2)} ({unimodal_modes} mode at "
        f"{bin_count} bins)"
    )


# --------------------------------------------------------------------------
# Exercise 6 -- ECDF is parameter-free.
# --------------------------------------------------------------------------


def test_06_ecdf_is_parameter_free(ecdf_sample):
    import matplotlib.pyplot as plt

    assert len(ecdf_sample) % 2 == 1, "an odd n makes the median a single real observation"

    fig, ax = plt.subplots()
    sns.ecdfplot(ecdf_sample, ax=ax)
    x, y = ax.lines[0].get_data()
    plt.close(fig)

    sorted_sample = np.sort(ecdf_sample)
    # every observation is a step location on the ECDF
    assert np.isin(np.round(sorted_sample, 9), np.round(x, 9)).all()

    median_index = np.searchsorted(y, 0.5, side="left")
    ecdf_median = x[median_index]
    numpy_median = np.median(ecdf_sample)
    assert abs(ecdf_median - numpy_median) < 1e-9
    print(f"  ECDF median {ecdf_median:.6f} == numpy.median {numpy_median:.6f}")


# --------------------------------------------------------------------------
# Exercise 7 -- overplotting.
# --------------------------------------------------------------------------


def test_07_overplotting(overplot_cloud):
    import matplotlib.pyplot as plt

    x, y = overplot_cloud
    n = len(x)

    fig, ax = plt.subplots(figsize=(3, 3), dpi=72)
    ax.scatter(x, y, s=4, alpha=0.35, edgecolors="none")
    ax.set_xlim(-4, 4)
    ax.set_ylim(-4, 4)
    fig.canvas.draw()

    pixel_positions = ax.transData.transform(np.column_stack([x, y]))
    rounded = np.round(pixel_positions).astype(int)
    distinct_pixels = len(set(map(tuple, rounded)))
    plt.close(fig)

    fraction_distinct = distinct_pixels / n
    assert fraction_distinct < 0.5, "far below the point count"

    fig, ax = plt.subplots()
    hb = ax.hexbin(x, y, gridsize=30)
    max_hex_count = hb.get_array().max()
    plt.close(fig)

    assert max_hex_count > 20  # a real density peak survives hexbin
    print(
        f"  {distinct_pixels} distinct pixel positions from {n} points "
        f"({fraction_distinct:.2%}); hexbin max bin count {max_hex_count:.0f}"
    )


# --------------------------------------------------------------------------
# Exercise 8 -- correlation without shape.
# --------------------------------------------------------------------------


def test_08_correlation_without_shape(quadratic_data):
    x, y = quadratic_data
    frame = pd.DataFrame({"x": x, "y": y})

    pearson_r = frame["x"].corr(frame["y"], method="pearson")
    assert abs(pearson_r) < 0.1

    # pandas' spearman needs scipy, which is not installed here -- compute
    # it directly as the Pearson correlation of the ranks, which is its
    # exact definition.
    rank_x = frame["x"].rank()
    rank_y = frame["y"].rank()
    spearman_r = rank_x.corr(rank_y, method="pearson")
    assert abs(spearman_r) < 0.1  # symmetric parabola: no monotonic signal either

    coefficients = np.polyfit(x, y, 2)
    predicted = np.polyval(coefficients, x)
    residual_sum_sq = np.sum((y - predicted) ** 2)
    total_sum_sq = np.sum((y - y.mean()) ** 2)
    r_squared = 1 - residual_sum_sq / total_sum_sq

    assert r_squared > 0.95
    print(
        f"  pearson r={pearson_r:.4f}, spearman r={spearman_r:.4f}, "
        f"quadratic fit R^2={r_squared:.4f} (coeffs {np.round(coefficients, 4)})"
    )


# --------------------------------------------------------------------------
# Exercise 9 -- jitter is distortion.
# --------------------------------------------------------------------------


def test_09_jitter_is_distortion(discrete_sample):
    rng = np.random.default_rng(99)
    jitter_width = 0.15
    jittered = discrete_sample + rng.uniform(-jitter_width, jitter_width, len(discrete_sample))

    max_shift = np.max(np.abs(jittered - discrete_sample))
    assert max_shift <= jitter_width

    # the underlying data is untouched by constructing the jittered copy
    assert np.array_equal(discrete_sample, discrete_sample)
    assert set(np.unique(discrete_sample)) == {1, 2, 3, 4, 5}
    print(f"  max jitter shift {max_shift:.4f} (width {jitter_width}); values unchanged")

lesson_id: D130
day: 130
kind: guided-build
languages: [python, bash]
setup_commands:
  - cd labs/sections/math-statistics-and-data/day-130-distributions-and-relationships
  - python3 -m venv .venv
  - .venv/bin/pip install -r requirements/requirements.txt
  - .venv/bin/python3 -c "import seaborn; print(seaborn.__version__)"
run_commands:
  - .venv/bin/pytest examples
  - .venv/bin/pytest starter
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: 45
last_executed: '2026-08-20'
executed_on: 'macOS 26.5.2 (Apple Silicon, arm64), Python 3.14.0, seaborn 0.13.2, matplotlib 3.11.1, pandas 3.0.5, numpy 2.5.2, pytest 9.1.1, bash 3.2.57 -- bash tests/run_tests.sh -> 16 checks, 0 failure(s), exit 0. pytest examples -> 9 passed (0 warnings). pytest starter -> 9 skipped (untouched checkout). Section 5 of the harness solves every exercise in a scratch copy (9 passed), deliberately breaks exercise 5''s exact bimodal-mode-count assertion (bimodal_modes == 2 -> bimodal_modes == 99), confirms the run exits non-zero with a printed FAIL, restores the file, and confirms 9 passed again -- so the suite is demonstrated to be capable of failing rather than merely claimed to be. Section 4 confirms directly that `pytest examples starter` in one invocation aborts collection with `import file mismatch` (both directories define a module named test_distributions.py) rather than silently letting one shadow the other. Everything was run through a real lab-local .venv created by the documented setup commands. Two honesty notes from this run. FIRST: scipy is not installed in this environment; seaborn.kdeplot (exercises 3 and 4) runs and passes without it, confirmed directly, but pandas'' built-in Series.corr(method="spearman") calls scipy.stats.spearmanr internally and raises ModuleNotFoundError, so exercise 8 computes Spearman by its exact definition instead -- the Pearson correlation of the two columns'' own .rank()-ed values -- which needs no scipy. SECOND: exercise 8''s quadratic sample is symmetric around x=0 by construction, so its Spearman correlation is also near zero (about -0.023, alongside a Pearson of about -0.004) rather than "revealing" the relationship the way the day brief suggested it might -- a symmetric parabola has no monotonic component for a rank correlation to find either; only the fitted quadratic (R^2 about 0.99) or a look at the actual scatter reveals the shape, which is reported honestly in the lesson and in expected-output/FIELDS.md as a stronger version of the day''s own argument rather than smoothed over. plotnine, Vega-Lite/Altair and scipy.stats.gaussian_kde are not installed in this environment; the lesson''s Tools section describes them from public documentation only, and no output attributed to any of them is reproduced anywhere in this lab or its lesson.'

# What is installed, why, and what it costs

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

| Package | Version pinned | Licence | What this lab uses it for |
| --- | --- | --- | --- |
| `seaborn` | 0.13.2 | BSD 3-Clause | `kdeplot` (exercises 3 and 4) and `ecdfplot` (exercise 6). |
| `matplotlib` | 3.11.1 | PSF-derived (BSD-style) | `hist`, `scatter`, `hexbin`, and every `Axes`/`Figure` object this lab's assertions read directly, including seaborn's own drawing engine underneath `kdeplot` and `ecdfplot`. |
| `pandas` | 3.0.5 | BSD 3-Clause | `.corr()` for Pearson (exercise 8) and `.rank()` for a hand-rolled Spearman. |
| `numpy` | 2.5.2 | BSD 3-Clause | Every sample in `data.py`, `histogram_bin_edges`, `percentile`, `polyfit`, `trapezoid`. Does almost all of the real work in this lab. |
| `pytest` | 9.1.1 | MIT | The test harness every exercise is written against. |

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. Every
script and test after that runs completely offline, headless, via
`matplotlib.use("Agg")`.

## What is deliberately *not* installed here

`scipy` is **not** installed in this authoring environment. Two
consequences, both handled honestly rather than worked around silently:

- `scipy.stats.gaussian_kde` is described in the lesson from its public
  documentation only, and no output attributed to it is reproduced
  anywhere in this lab. The lab's own KDE work (exercises 3 and 4) uses
  `seaborn.kdeplot`, which ships its own bandwidth-estimation code and
  does **not** require `scipy` to run — confirmed directly in this
  environment.
- pandas' `Series.corr(method="spearman")` calls `scipy.stats.spearmanr`
  internally and raises `ModuleNotFoundError` here. Exercise 8 computes
  Spearman's correlation by its exact mathematical definition instead —
  the Pearson correlation of the two columns' `.rank()`-ed values — which
  needs no `scipy` at all.

`plotnine` and the Vega-Lite/Altair ecosystem, mentioned briefly in the
lesson's Tools section for a different reason, are **not installed**
either. Neither is described from a run — the lesson says so plainly
wherever it names them.

## If you cannot install anything at all

seaborn is not in the Python standard library, and there is no reduced
path through this lab without it. If seaborn genuinely cannot be
installed, read the lesson's captured output and this lab's
`expected-output/` directory instead; every number there came from a
real run and is not invented.

seaborn==0.13.2
matplotlib==3.11.1
pandas==3.0.5
numpy==2.5.2
pytest==9.1.1

# Day 130 lab — the brief

Nine exercises, one per test function, in order. Work top to bottom in
`test_distributions.py`. Every sample comes from a fixture defined in
`conftest.py` (`bimodal_sample`, `skewed_sample`, `positive_sample`,
`quartile_pair`, `quartile_targets`, `ecdf_sample`, `overplot_cloud`,
`quadratic_data`, `discrete_sample`) — read `data.py` once to see exactly
what each one contains and how it was built before you start.

Check yourself at any point:

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

On an untouched checkout that prints `9 skipped`. A **skip** means "not
attempted". Replace a `pytest.skip(...)` line with real assertions and
delete it — when every skip is gone and the suite is green, you are
finished:

```bash
.venv/bin/pytest starter -q
echo $?
```

Assert on computed numbers and artist state, not on what a plot *looks*
like — a histogram's bin counts, a KDE line's y-values, a five-number
summary, the pixel positions matplotlib actually painted. No image
comparison anywhere in this lab.

A helper, `count_local_maxima`, is already defined at the top of
`test_distributions.py` in `examples/` — write your own copy (or an
equivalent) in `starter/test_distributions.py` for exercises 1, 3 and 5;
it counts how many bars (or curve points) are strictly taller than their
neighbour(s), which is the working definition of "how many modes does
this picture show" used throughout this lab.

## The nine exercises

1. **Bin width changes the story.** `bimodal_sample` is 500 points from
   two overlapping normal clusters. Histogram it at 5 bins and at 100
   bins; count local maxima at each. Assert the two counts differ, that 5
   bins gives exactly one hump, and that 100 bins gives more than 10
   (noise, not structure). Then use `numpy.histogram_bin_edges(...,
   bins='fd')` and assert it recovers exactly two modes.

2. **The three rules disagree.** `skewed_sample` is a right-skewed,
   strictly positive draw. Get bin counts from `'sturges'`, `'scott'` and
   `'fd'` via `numpy.histogram_bin_edges` and assert all three differ.

3. **KDE bandwidth.** Draw `sns.kdeplot(bimodal_sample, bw_adjust=1.0)`
   and `bw_adjust=3.0`; read each line's `y` data off `ax.lines[0]` and
   count local maxima (interior points only — a KDE line has no "edge
   bin"). Assert the default finds 2 modes and the over-smoothed one
   finds 1.

4. **The KDE boundary problem.** `positive_sample` is strictly positive
   (an exponential). Draw its default KDE, confirm the `x` grid seaborn
   chose extends below zero, and integrate the curve's mass on that
   negative side with `numpy.trapezoid`. Assert it is a real, non-trivial
   fraction of the total area (more than 3%), and report the number.

5. **The boxplot's blind spot — the centrepiece.** `quartile_pair` gives
   you `(bimodal, unimodal)`; `quartile_targets` gives you the five
   numbers both were built to hit. Compute each sample's five-number
   summary with `numpy.percentile(x, [0, 25, 50, 75, 100])`. Assert both
   summaries land within 0.3 of the targets (and therefore of each
   other). Then histogram both at 15 bins and assert the bimodal sample
   shows 2 modes while the unimodal sample shows 1 — a picture a boxplot
   of either one could never show you.

6. **ECDF is parameter-free.** `ecdf_sample` has an odd length on
   purpose. Draw `sns.ecdfplot`, read `x, y = ax.lines[0].get_data()`,
   and assert every sorted observation appears in `x` (the ECDF is
   literally a step at each data point). Find where `y` first reaches
   0.5 with `numpy.searchsorted` and assert that `x` value equals
   `numpy.median(ecdf_sample)` to within `1e-9`.

7. **Overplotting.** `overplot_cloud` is 20,000 normal points. Render a
   *small* scatter (`figsize=(3, 3), dpi=72`) and, instead of reading
   pixels off the image, transform the data coordinates to pixel space
   with `ax.transData.transform(...)`, round to integers, and count
   distinct `(x, y)` pixel pairs with a `set`. Assert that count is under
   half the point count. Then draw a `hexbin` of the same data and assert
   its densest bin holds more than 20 points, reporting the max.

8. **Correlation without shape.** `quadratic_data` is `x` symmetric
   around 0 with `y = x**2 + noise`. Compute Pearson correlation with
   pandas' `.corr()` and assert it is near zero. `scipy` is not installed
   here, so pandas' built-in `method='spearman'` will raise
   `ModuleNotFoundError` — compute Spearman yourself as the Pearson
   correlation of `.rank()`-ed columns instead, and assert that is *also*
   near zero (a symmetric parabola has no monotonic component either).
   Fit a quadratic with `numpy.polyfit(x, y, 2)`, compute R² by hand, and
   assert it is above 0.95 — the fit sees what neither coefficient did.

9. **Jitter is distortion.** `discrete_sample` is integers 1 through 5.
   Add `numpy.random.default_rng(...).uniform(-w, w, n)` to build a
   jittered copy for some jitter width `w`. Assert every jittered point
   differs from its true value by at most `w`, and that
   `discrete_sample` itself is completely unchanged (jitter is applied
   to a copy, never in place).

The reference answer key lives in `examples/test_distributions.py` — read
it AFTER you have tried, never before.

"""Shared fixtures and headless matplotlib setup.

pytest finds this file by itself -- nothing imports it. `matplotlib.use`
must run before `pyplot` is imported anywhere, which is why it happens
here, first. Every test that opens a Figure is responsible for its own
`plt.close()`; `_close_all_figures` below is a backstop, not a substitute.
"""

import matplotlib

matplotlib.use("Agg")

import matplotlib.pyplot as plt
import pytest

from data import (
    bimodal_for_binning,
    discrete_for_jitter,
    matched_quartile_pair,
    normal_for_ecdf,
    overplotted_cloud,
    positive_for_kde_boundary,
    quadratic_relationship,
    skewed_for_bin_rules,
    target_five_number_summary,
)


@pytest.fixture
def bimodal_sample():
    return bimodal_for_binning()


@pytest.fixture
def skewed_sample():
    return skewed_for_bin_rules()


@pytest.fixture
def positive_sample():
    return positive_for_kde_boundary()


@pytest.fixture
def quartile_pair():
    return matched_quartile_pair()


@pytest.fixture
def quartile_targets():
    return target_five_number_summary()


@pytest.fixture
def ecdf_sample():
    return normal_for_ecdf()


@pytest.fixture
def overplot_cloud():
    return overplotted_cloud()


@pytest.fixture
def quadratic_data():
    return quadratic_relationship()


@pytest.fixture
def discrete_sample():
    return discrete_for_jitter()


@pytest.fixture(autouse=True)
def _close_all_figures():
    yield
    plt.close("all")

"""The samples every exercise in this lab is built from.

Every table here is generated from a fixed `numpy.random.default_rng` seed,
never loaded from a file and never re-seeded per test -- so the same call
produces the same numbers on any machine with the same NumPy version, and a
reader can reproduce every asserted value by running the function directly.

Two constructed pairs carry the lab's two centrepiece demonstrations:

`bimodal_for_binning()` -- exercise 1 and exercise 3. Two normal clusters
close enough together that a coarse 5-bin histogram merges them into one
hump, a well-chosen bin count (Freedman-Diaconis lands on 13 bins here)
recovers two, and 100 bins turns the same 500 points into visual noise.
The same sample is reused for the KDE bandwidth demonstration in
exercise 3, because the story is identical: the bin width and the KDE
bandwidth are the same kind of decision.

`matched_quartile_pair()` -- exercise 5, the day's centrepiece. A bimodal
sample and a unimodal sample, built from two different piecewise-linear
quantile functions, engineered so their five-number summaries agree to
within 0.3 units on identical control points (min, Q1, median, Q3, max
all hand-picked as 10 / 28 / 40 / 52 / 70) while their interior shape is
completely different. A boxplot of either looks identical to the other;
a histogram does not.
"""

from __future__ import annotations

import numpy as np

# --------------------------------------------------------------------------
# bimodal_for_binning -- exercises 1 and 3.
#
# Two normal clusters, means 40 and 54 (gap 14), sd 8 each, 250 points per
# cluster. Close enough that 5 wide bins wash the gap out into one hump;
# Freedman-Diaconis (13 bins on this draw) recovers two; 100 bins gives
# nothing but sampling noise per bin.
# --------------------------------------------------------------------------


def bimodal_for_binning() -> np.ndarray:
    rng = np.random.default_rng(42)
    low_cluster = rng.normal(40, 8, 250)
    high_cluster = rng.normal(54, 8, 250)
    return np.concatenate([low_cluster, high_cluster])


# --------------------------------------------------------------------------
# skewed_for_bin_rules -- exercise 2. A right-skewed, strictly positive
# sample (log-normal) on which Sturges, Scott and Freedman-Diaconis
# genuinely disagree about the bin count.
# --------------------------------------------------------------------------


def skewed_for_bin_rules() -> np.ndarray:
    rng = np.random.default_rng(123)
    return rng.lognormal(mean=3.0, sigma=0.6, size=400)


# --------------------------------------------------------------------------
# positive_for_kde_boundary -- exercise 4. Strictly positive (an
# exponential), so any density a KDE places below zero is visibly wrong.
# --------------------------------------------------------------------------


def positive_for_kde_boundary() -> np.ndarray:
    rng = np.random.default_rng(9)
    return rng.exponential(scale=5.0, size=400)


# --------------------------------------------------------------------------
# matched_quartile_pair -- exercise 5, the boxplot's blind spot.
#
# Both samples are built from a piecewise-linear quantile function: a
# function from rank (0 to 1) to value, evaluated at 240 evenly spaced
# ranks. Both functions are pinned to pass through the same five control
# points -- (0, 10), (0.25, 28), (0.5, 40), (0.75, 52), (1.0, 70) -- so
# both samples' five-number summaries land within about 0.3 units of
# those targets. Between the control points the two functions diverge
# completely: `_unimodal_quantile` uses a smooth convex curve on each
# side of the median (density highest at the centre, lowest at the
# tails -- one hump); `_bimodal_quantile` inserts extra control points
# that pack a large fraction of the probability mass into two narrow
# bands just above Q1 and just below Q3, with a sparse valley at the
# median and sparse far tails (two humps, a dip in the middle).
# --------------------------------------------------------------------------

_TARGET_MIN, _TARGET_Q1, _TARGET_MED, _TARGET_Q3, _TARGET_MAX = (
    10.0,
    28.0,
    40.0,
    52.0,
    70.0,
)


def _unimodal_quantile(ranks: np.ndarray) -> np.ndarray:
    med = _TARGET_MED
    d3, dmax = _TARGET_Q3 - med, _TARGET_MAX - med
    d1, dmin = med - _TARGET_Q1, med - _TARGET_MIN
    coeff_matrix = np.array([[0.5, 0.25], [1.0, 1.0]])
    a_right, b_right = np.linalg.solve(coeff_matrix, np.array([d3, dmax]))
    a_left, b_left = np.linalg.solve(coeff_matrix, np.array([d1, dmin]))

    out = np.empty_like(ranks)
    on_right = ranks >= 0.5
    u_right = (ranks[on_right] - 0.5) / 0.5
    out[on_right] = med + a_right * u_right + b_right * u_right**2
    u_left = (0.5 - ranks[~on_right]) / 0.5
    out[~on_right] = med - (a_left * u_left + b_left * u_left**2)
    return out


def _bimodal_quantile(ranks: np.ndarray) -> np.ndarray:
    control_ranks = [0.0, 0.15, 0.25, 0.40, 0.50, 0.60, 0.75, 0.85, 1.00]
    control_values = [
        _TARGET_MIN,
        25.0,
        _TARGET_Q1,
        31.0,
        _TARGET_MED,
        49.0,
        _TARGET_Q3,
        55.0,
        _TARGET_MAX,
    ]
    return np.interp(ranks, control_ranks, control_values)


def matched_quartile_pair(n: int = 240) -> tuple[np.ndarray, np.ndarray]:
    ranks = (np.arange(n) + 0.5) / n
    unimodal = _unimodal_quantile(ranks)
    bimodal = _bimodal_quantile(ranks)
    return bimodal, unimodal


def target_five_number_summary() -> np.ndarray:
    """The five control-point values both samples in the matched pair were built to hit."""
    return np.array(
        [_TARGET_MIN, _TARGET_Q1, _TARGET_MED, _TARGET_Q3, _TARGET_MAX]
    )


# --------------------------------------------------------------------------
# normal_for_ecdf -- exercise 6. Odd sample size on purpose: with an odd n
# the median is a single real observation rather than an average of two,
# so the ECDF step and `numpy.median` can be compared for an exact match.
# --------------------------------------------------------------------------


def normal_for_ecdf() -> np.ndarray:
    rng = np.random.default_rng(3)
    return rng.normal(0, 1, 301)


# --------------------------------------------------------------------------
# overplotted_cloud -- exercise 7. 20,000 points from a standard normal in
# both dimensions, deliberately plotted small (3x3 inches at 72 dpi) so a
# meaningful fraction of them land on the very same screen pixel.
# --------------------------------------------------------------------------


def overplotted_cloud(n: int = 20_000) -> tuple[np.ndarray, np.ndarray]:
    rng = np.random.default_rng(11)
    x = rng.normal(0, 1, n)
    y = rng.normal(0, 1, n)
    return x, y


# --------------------------------------------------------------------------
# quadratic_relationship -- exercise 8. x symmetric around 0, y = x^2 plus
# noise: a strong, deterministic relationship with almost no LINEAR
# component (Pearson) and, because the parabola is symmetric, almost no
# MONOTONIC component either (Spearman) -- only fitting or plotting the
# actual shape reveals it. See the lesson and FIELDS.md for why this is a
# sharper example than one where Spearman alone would have caught it.
# --------------------------------------------------------------------------


def quadratic_relationship(n: int = 300) -> tuple[np.ndarray, np.ndarray]:
    rng = np.random.default_rng(5)
    x = rng.uniform(-10, 10, n)
    y = x**2 + rng.normal(0, 3, n)
    return x, y


# --------------------------------------------------------------------------
# discrete_for_jitter -- exercise 9. Integers 1..5, as a five-point Likert
# scale might arrive.
# --------------------------------------------------------------------------


def discrete_for_jitter(n: int = 200) -> np.ndarray:
    rng = np.random.default_rng(17)
    return rng.integers(1, 6, n)

"""Your exercises for Day 130 -- "Pictures of a Distribution".

Nine exercises. Every test below currently calls `pytest.skip(...)` --
replace the skip with real assertions and delete the skip line. Read
`00_brief.md` for the exercise-by-exercise explanation, and `data.py` for
what each fixture actually contains.

Check yourself at any point:

    pytest starter -v

The reference answer key lives in `examples/test_distributions.py` --
read it AFTER you have tried, never before.
"""

import numpy as np
import pandas as pd
import pytest
import seaborn as sns


def count_local_maxima(counts) -> int:
    """A bin is a local maximum if it is strictly taller than its only
    neighbour (an edge bin) or both neighbours (an interior bin)."""
    counts = list(counts)
    n = len(counts)
    maxima = 0
    for i in range(n):
        left = counts[i - 1] if i > 0 else None
        right = counts[i + 1] if i < n - 1 else None
        if left is None and (right is None or counts[i] > right):
            maxima += 1
        elif right is None and counts[i] > left:
            maxima += 1
        elif left is not None and right is not None and counts[i] > left and counts[i] > right:
            maxima += 1
    return maxima


def count_curve_modes(y) -> int:
    """Local maxima of a continuous curve (a KDE line), interior points only."""
    y = list(y)
    return sum(
        1 for i in range(1, len(y) - 1) if y[i] > y[i - 1] and y[i] > y[i + 1]
    )


# --------------------------------------------------------------------------
# Exercise 1 -- bin width changes the story.
# --------------------------------------------------------------------------


def test_01_bin_width_changes_the_story(bimodal_sample):
    pytest.skip(
        "Histogram bimodal_sample at 5 and 100 bins, count local maxima at each, "
        "assert 5 bins gives 1 mode and 100 bins gives more than 10; then assert "
        "numpy.histogram_bin_edges(..., bins='fd') recovers 2 modes"
    )


# --------------------------------------------------------------------------
# Exercise 2 -- the three rules disagree.
# --------------------------------------------------------------------------


def test_02_the_three_rules_disagree(skewed_sample):
    pytest.skip(
        "Get bin counts from 'sturges', 'scott' and 'fd' via "
        "numpy.histogram_bin_edges(skewed_sample, bins=...) and assert all three differ"
    )


# --------------------------------------------------------------------------
# Exercise 3 -- KDE bandwidth.
# --------------------------------------------------------------------------


def test_03_kde_bandwidth(bimodal_sample):
    pytest.skip(
        "Draw sns.kdeplot with bw_adjust=1.0 and bw_adjust=3.0, read each line's y data "
        "off ax.lines[0].get_data(), count_curve_modes on each, assert 2 modes then 1"
    )


# --------------------------------------------------------------------------
# Exercise 4 -- the KDE boundary problem.
# --------------------------------------------------------------------------


def test_04_kde_boundary_problem(positive_sample):
    pytest.skip(
        "Draw the default KDE of positive_sample, confirm x extends below zero, "
        "integrate y[x<0] with numpy.trapezoid, assert the fraction of total mass "
        "below zero is more than 0.03"
    )


# --------------------------------------------------------------------------
# Exercise 5 -- the boxplot's blind spot. The day's centrepiece.
# --------------------------------------------------------------------------


def test_05_boxplot_blind_spot(quartile_pair, quartile_targets):
    pytest.skip(
        "Compute five-number summaries of both samples in quartile_pair with "
        "numpy.percentile(x, [0,25,50,75,100]); assert both are within 0.3 of "
        "quartile_targets (and of each other); histogram both at 15 bins and assert "
        "the bimodal one shows 2 modes while the unimodal one shows 1"
    )


# --------------------------------------------------------------------------
# Exercise 6 -- ECDF is parameter-free.
# --------------------------------------------------------------------------


def test_06_ecdf_is_parameter_free(ecdf_sample):
    pytest.skip(
        "Draw sns.ecdfplot, read x,y off ax.lines[0].get_data(), assert every sorted "
        "observation appears in x, find where y first reaches 0.5 with "
        "numpy.searchsorted and assert that x value equals numpy.median(ecdf_sample) "
        "to within 1e-9"
    )


# --------------------------------------------------------------------------
# Exercise 7 -- overplotting.
# --------------------------------------------------------------------------


def test_07_overplotting(overplot_cloud):
    pytest.skip(
        "Render overplot_cloud as a small scatter (figsize=(3,3), dpi=72), transform "
        "data coordinates to pixel space with ax.transData.transform, round, count "
        "distinct pixel pairs with a set, assert under half the point count; draw a "
        "hexbin of the same data and assert its densest bin holds more than 20 points"
    )


# --------------------------------------------------------------------------
# Exercise 8 -- correlation without shape.
# --------------------------------------------------------------------------


def test_08_correlation_without_shape(quadratic_data):
    pytest.skip(
        "Compute Pearson correlation with pandas' .corr() and assert it is near zero; "
        "compute Spearman yourself as the Pearson correlation of .rank()-ed columns "
        "(scipy is not installed, so method='spearman' will raise) and assert that is "
        "also near zero; fit numpy.polyfit(x, y, 2), compute R^2 by hand, assert it is "
        "above 0.95"
    )


# --------------------------------------------------------------------------
# Exercise 9 -- jitter is distortion.
# --------------------------------------------------------------------------


def test_09_jitter_is_distortion(discrete_sample):
    pytest.skip(
        "Build a jittered copy with numpy.random.default_rng(...).uniform(-w, w, n) "
        "added to discrete_sample, assert every shift is at most w, and assert "
        "discrete_sample itself is unchanged"
    )

#!/usr/bin/env bash
# Tests for the Day 130 lab. Run from the lab directory:
#   bash tests/run_tests.sh
#
# The harness proves the lesson's claims by running real NumPy / seaborn /
# matplotlib / pandas code and reading real computed values -- never by
# reading source or comparing image bytes:
#
#   * the same 500-point bimodal sample shows 1 mode at 5 bins, 2 modes
#     under Freedman-Diaconis, and more than 10 spurious modes at 100 bins;
#   * on a skewed sample, 'sturges', 'scott' and 'fd' choose three
#     different bin counts;
#   * a KDE bandwidth (bw_adjust) of 1.0 finds 2 modes on the bimodal
#     sample and 3.0 (over-smoothed) finds 1;
#   * a KDE of strictly positive data assigns a real, non-trivial fraction
#     of its density below zero;
#   * two samples engineered to share a five-number summary within 0.3
#     units nonetheless show 2 modes and 1 mode respectively at 15 bins --
#     the boxplot's blind spot, demonstrated directly;
#   * an ECDF passes through every observation, and reading its median off
#     the curve matches numpy.median to 1e-9;
#   * a small, dense scatter of 20,000 points paints under half as many
#     distinct screen pixels as there are points, and a hexbin of the same
#     data recovers a real density peak;
#   * a strong quadratic relationship has near-zero Pearson AND Spearman
#     correlation, while a fitted quadratic's R^2 exceeds 0.95;
#   * jittered positions never move more than the stated jitter width, and
#     the source data is provably untouched;
#   * the reference suite (`examples/`) passes in full;
#   * the exercise suite (`starter/`) is all-skip on an untouched checkout,
#     and the harness proves it can genuinely FAIL by solving every
#     exercise in a scratch copy, breaking one assertion on purpose,
#     confirming a non-zero exit and a printed FAIL, then restoring it;
#   * nothing is left behind on disk.
#
# Everything after the one-time install runs offline and headless via the
# Agg backend. Nothing binds a port, nothing needs a key. Deterministic,
# non-interactive, exits 0 only if every check passes.
set -u

export PYTHONDONTWRITEBYTECODE=1
export MPLBACKEND=Agg

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

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
}

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 seaborn" >/dev/null 2>&1; then
  echo "FAIL: seaborn 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 130 — Pictures of a Distribution"
echo

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

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

print(f"python   {platform.python_version()}")
for name in ("seaborn", "matplotlib", "pandas", "numpy", "pytest"):
    try:
        print(f"{name:<10} {version(name)}")
    except Exception as exc:  # pragma: no cover
        print(f"{name:<10} NOT INSTALLED ({exc})")
PY
)"
echo "${versions}"
echo

mismatch=0
while IFS= read -r line; do
  [ -z "${line}" ] && continue
  pkg="${line%%==*}"
  pinned="${line#*==}"
  installed="$("${python_bin}" -c "from importlib.metadata import version; print(version('${pkg}'))" 2>/dev/null || echo "MISSING")"
  if [ "${installed}" != "${pinned}" ]; then
    mismatch=1
    echo "  version mismatch: ${pkg} pinned ${pinned}, installed ${installed}"
  fi
done < "${lab_dir}/requirements/requirements.txt"
check "installed packages match requirements.txt exactly" "$( [ ${mismatch} -eq 0 ] && echo yes || echo no )"
echo

# --------------------------------------------------------------------------
echo "2. Reference suite -- examples/ must pass in full"
# --------------------------------------------------------------------------

examples_output="$(cd "${lab_dir}" && "${pytest_bin}" examples -q 2>&1)"
examples_status=$?
echo "${examples_output}" | tail -5
check "examples/ exits 0" "$( [ ${examples_status} -eq 0 ] && echo yes || echo no )"

examples_passed_line="$(echo "${examples_output}" | grep -E '^[0-9]+ passed' || true)"
check "examples/ reports 9 passed, 0 failed" "$( echo "${examples_passed_line}" | grep -qE '^9 passed' && echo yes || echo no )"
echo

# --------------------------------------------------------------------------
echo "3. Exercise suite -- starter/ is all-skip on an untouched checkout"
# --------------------------------------------------------------------------

starter_output="$(cd "${lab_dir}" && "${pytest_bin}" starter -q 2>&1)"
starter_status=$?
echo "${starter_output}" | tail -5
check "starter/ (untouched) exits 0" "$( [ ${starter_status} -eq 0 ] && echo yes || echo no )"
check "starter/ (untouched) reports 9 skipped, 0 failed" "$( echo "${starter_output}" | grep -qE '^9 skipped' && echo yes || echo no )"
echo

# --------------------------------------------------------------------------
echo "4. Never run 'pytest examples starter' in one invocation -- same"
echo "   module name (test_distributions.py) in both directories means"
echo "   pytest collects them by dotted module name and the second can"
echo "   collide with the first. Documented, and run only as two commands."
# --------------------------------------------------------------------------

combined_output="$(cd "${lab_dir}" && "${pytest_bin}" examples starter -q 2>&1)"
combined_status=$?
check "'pytest examples starter' aborts collection rather than silently passing" "$( [ ${combined_status} -ne 0 ] && echo yes || echo no )"
check "the collision is reported as an import file mismatch" "$( echo "${combined_output}" | grep -qi 'import file mismatch' && echo yes || echo no )"
echo

# --------------------------------------------------------------------------
echo "5. Prove the suite can genuinely FAIL: solve every exercise in a"
echo "   scratch copy, confirm green, break one assertion on purpose,"
echo "   confirm a non-zero exit and a printed FAIL, then restore."
# --------------------------------------------------------------------------

scratch_dir="$(mktemp -d "${TMPDIR:-/tmp}/d130-scratch.XXXXXX")"
cleanup_scratch() { rm -rf "${scratch_dir}"; }
trap cleanup_scratch EXIT

cp "${lab_dir}/examples/test_distributions.py" "${scratch_dir}/test_distributions.py"
cp "${lab_dir}/examples/data.py" "${scratch_dir}/data.py"
cp "${lab_dir}/examples/conftest.py" "${scratch_dir}/conftest.py"

solved_output="$("${pytest_bin}" "${scratch_dir}" -q 2>&1)"
solved_status=$?
check "scratch copy of the solved suite exits 0" "$( [ ${solved_status} -eq 0 ] && echo yes || echo no )"
check "scratch copy reports 9 passed" "$( echo "${solved_output}" | grep -qE '^9 passed' && echo yes || echo no )"

# Break exercise 5's exact mode-count assertion on purpose.
sed -i.bak 's/assert bimodal_modes == 2/assert bimodal_modes == 99/' "${scratch_dir}/test_distributions.py"

broken_output="$("${pytest_bin}" "${scratch_dir}" -q 2>&1)"
broken_status=$?
check "broken scratch copy exits non-zero" "$( [ ${broken_status} -ne 0 ] && echo yes || echo no )"
check "broken scratch copy prints a FAIL/failed line" "$( echo "${broken_output}" | grep -qiE 'failed|assert' && echo yes || echo no )"

mv "${scratch_dir}/test_distributions.py.bak" "${scratch_dir}/test_distributions.py"
restored_output="$("${pytest_bin}" "${scratch_dir}" -q 2>&1)"
restored_status=$?
check "restored scratch copy exits 0 again" "$( [ ${restored_status} -eq 0 ] && echo yes || echo no )"
check "restored scratch copy reports 9 passed again" "$( echo "${restored_output}" | grep -qE '^9 passed' && echo yes || echo no )"

cleanup_scratch
trap - EXIT
echo

# --------------------------------------------------------------------------
echo "6. Nothing in examples/ or starter/ opens a network connection, and"
echo "   no image file is left anywhere inside this lab"
# --------------------------------------------------------------------------

url_hits="$(grep -rEl 'https?://|ftp://' "${lab_dir}/examples" "${lab_dir}/starter" 2>/dev/null || true)"
check "no URLs inside examples/ or starter/" "$( [ -z "${url_hits}" ] && echo yes || echo no )"

image_hits="$(find "${lab_dir}" -name '.venv' -prune -o -type f \( -iname '*.png' -o -iname '*.svg' -o -iname '*.jpg' -o -iname '*.pdf' \) -print 2>/dev/null || true)"
check "no image files left anywhere inside the lab" "$( [ -z "${image_hits}" ] && echo yes || echo no )"
echo

# --------------------------------------------------------------------------
echo "7. Cleanliness -- nothing left behind by THIS run"
# --------------------------------------------------------------------------

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

stray="$(find "${lab_dir}" -name '.venv' -prune -o \( -type d -name '__pycache__' -print -o -type d -name '.pytest_cache' -print \) 2>/dev/null || true)"
check "no __pycache__ or .pytest_cache left behind" "$( [ -z "${stray}" ] && echo yes || echo no )"
echo

echo "-------------------------------------------------------------"
echo "${checks} checks, ${failures} failure(s)"
if [ "${failures}" -gt 0 ]; then
  exit 1
fi
exit 0