import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans
from sklearn.compose import make_column_transformer
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score, roc_curve, silhouette_score
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import OneHotEncoder, StandardScaler

# - Load data/songs_2023.csv as songs (parse the release_date column) and data/catalog.csv as catalog.
# - Make tracks: one row per track_id, leaving out the tracks whose danceability is blank (the data notebook
#   set it to blank where Spotify's rhythm analysis returned nothing).
# - Show, as one small table, the number of songs, catalog rows, unique tracks in the catalog, and
#   tracks whose sound was measured.
import sys; sys.path.append(".claude/skills/fresh-figures/scripts")
import fresh_style as ms  # one style for every table and figure

ms.use_theme()  # the theme theme.json names

SEED = 6120  # every random step in the notebook uses this seed

songs = pd.read_csv("data/songs_2023.csv", parse_dates=["release_date"], low_memory=False)
catalog = pd.read_csv("data/catalog.csv", low_memory=False)
# One row per track, keeping its first genre, minus the tracks Spotify could not measure.
tracks = catalog.drop_duplicates("track_id").dropna(subset=["danceability_pct"]).reset_index(drop=True)
assert catalog["danceability_pct"].isna().sum() == 157  # the data notebook's missing rhythm analyses

display(ms.table(pd.Series({
    "songs on the 2023 list": len(songs),
    "catalog rows": len(catalog),
    "unique tracks in the catalog": catalog["track_id"].nunique(),
    "tracks whose sound was measured": len(tracks),
}, name="count").to_frame(), index_label="what we counted"))

# Meet the measure before testing it: hear the scale, then see where genres, artists and the list's music sit.
from mixtape import danceability_primer

display(danceability_primer(songs, tracks, catalog))

# The median danceability of the 2023 list and of every track, as one small table.
medians = pd.Series({
    "2023 list": songs["danceability_pct"].median(),
    "every track in the catalog": tracks["danceability_pct"].median(),
}, name="median danceability (0-100)")
assert medians.round(1).tolist() == [69.0, 57.6]  # the numbers the introduction quotes
display(ms.table(medians.round(1).to_frame(), index_label="group"))

# Figure: danceability on the 2023 list compared with every track in the
#   catalog, from songs and tracks.
import sys; sys.path.append(".claude/skills/fresh-figures/scripts")
import fresh_style as ms

ms.use_theme()
fig, ax = plt.subplots(figsize=(8, 4.5))
bins = np.arange(0, 101, 5)
# Shares of each group, not counts: the catalog has about 94 times as many tracks as the list.
share = lambda s: np.ones(len(s)) / len(s)
ax.hist(tracks["danceability_pct"], bins=bins, weights=share(tracks), color=ms.CONTEXT)
ax.hist(songs["danceability_pct"], bins=bins, weights=share(songs), histtype="step", linewidth=2.2,
        color=ms.HIGHLIGHT)
ax.yaxis.set_major_formatter(lambda v, _: f"{v:.0%}")
ax.grid(axis="y")
ax.set_axisbelow(True)
top = ax.get_ylim()[1] * 1.12
ax.set_ylim(0, top)

# Each group's median as a labeled line over the bars, solid for the list and dashed for the catalog.
for group, data, color, style, side in [
    ("Catalog", tracks, ms.INK, "--", "right"),
    ("2023 list", songs, ms.HIGHLIGHT, "-", "left"),
]:
    m = data["danceability_pct"].median()
    ax.axvline(m, color=color, linestyle=style, linewidth=1.4, zorder=3)
    ax.annotate(f"{group} median: {m:.1f}", xy=(m, top), xytext=(-6 if side == "right" else 6, -4),
                textcoords="offset points", ha=side, va="top", fontsize=10, color=color,
                fontweight="bold")

# Direct labels in place of a legend.
ax.annotate(f"2023 list\n{len(songs):,} songs", xy=(87, 0.10), ha="left", fontsize=10,
            color=ms.HIGHLIGHT)
ax.annotate(f"Every track in\nthe catalog\n{len(tracks):,} tracks", xy=(3, 0.06), ha="left",
            fontsize=10, color=ms.INK)

ax.set_xlim(0, 100)
ax.set_xlabel("Danceability (0–100)")
ax.set_ylabel("Share of the group (%)")
gap = songs["danceability_pct"].median() - tracks["danceability_pct"].median()
ms.finish(ax, f"The 2023 hits' median is {gap:.1f} points more danceable than the catalog's",
          f"Danceability in 5-point bins · share of each group · {len(songs):,} songs, "
          f"{len(tracks):,} catalog tracks",
          "Source: data/songs_2023.csv, data/catalog.csv (one row per track; the 157 tracks with no "
          "rhythm analysis are left out)", tag=ms.tag(1, 17))
ms.save(fig, n=1, slug="danceability-list-vs-catalog")
plt.show()

# Play the list's five most-streamed songs.
from mixtape import mixtape

display(mixtape(songs.nlargest(5, "streams_millions"), "The 2023 list's 5 most-streamed songs", "danceability_pct"))

# Make your calls before any test runs; each verdict below shows your call beside the data's.
from IPython.display import HTML

CALLS = "hit-song"  # where this notebook's calls are saved in the browser
display(HTML(ms.predictions(CALLS, [
    ("Q1·H1", "More danceable songs have more streams."),
    ("Q1·H2", "More upbeat songs, those with higher valence, have more streams."),
    ("Q1·H3", "Something other than how a song sounds explains its streams."),
    ("Q2·H1", "Danceability alone would have singled out the hits."),
    ("Q2·H2", "All the audio features together would have singled out the hits."),
    ("Q2·H3", "Adding genre, known before release, to the audio features would have singled out the hits."),
])))

# - The median streams in each quarter of the list by danceability: group songs by
#   pd.qcut(songs["danceability_pct"], 4) (with observed=True) and take the median of streams.
#   Call it dance, and show it as a table.
# - The slope (np.polyfit) and r (.corr) of streams on danceability across every song with a stream
#   count, as a small table.
def quarter_medians(df, x, y):
    """Median of y in each of four equal-sized groups of df by x, with each group's size."""
    groups = df.groupby(pd.qcut(df[x], 4), observed=True)[y]
    return groups.median(), groups.count()


def slope_and_r(df, x, y):
    """Slope (np.polyfit) and r (.corr) of y on x, over the rows where both are known."""
    known = df.dropna(subset=[x, y])
    slope = np.polyfit(known[x], known[y], 1)[0]
    return pd.Series({"slope": slope, "r": known[x].corr(known[y]), "songs": len(known)})


dance, dance_n = quarter_medians(songs, "danceability_pct", "streams_millions")
display(ms.table(pd.DataFrame({"median streams (millions)": dance.round(1), "songs": dance_n})))
display(ms.table(slope_and_r(songs, "danceability_pct", "streams_millions").to_frame("streams on danceability").round(3), index_label="measure"))

# Figure: the test of hypothesis 1, from dance and songs.
fig, (left, right) = plt.subplots(1, 2, figsize=(11, 4.5))

# Left: the median streams of each danceability quarter, with each bar's median and number of songs.
known = songs.dropna(subset=["streams_millions"])
per_quarter = known.groupby(pd.qcut(songs["danceability_pct"], 4), observed=True).size()
ranges = [f"{int(np.floor(q.left)) + 1}–{int(q.right)}" for q in dance.index]
bars = left.bar(ranges, dance.values, color=ms.CONTEXT, width=0.65)
for bar, median, n in zip(bars, dance.values, per_quarter.values):
    middle = bar.get_x() + bar.get_width() / 2
    left.annotate(f"{median:.1f}", xy=(middle, median), xytext=(0, 4), textcoords="offset points",
                  ha="center", va="bottom", fontsize=10, color=ms.INK, fontweight="bold")
    left.annotate(f"n = {n}", xy=(middle, 0), xytext=(0, 6), textcoords="offset points",
                  ha="center", va="bottom", fontsize=9, color=ms.on(ms.CONTEXT))
left.set_ylim(0, dance.max() * 1.18)
left.set_xlabel("Danceability quarter (0–100)")
left.set_ylabel("Median streams (millions)")
left.set_title("Median streams by quarter", color=ms.MUTED, fontsize=10)

# Right: every song with a stream count, and the regression line of streams on danceability.
slope, intercept, r = ms.regression(known["danceability_pct"], known["streams_millions"])
right.scatter(known["danceability_pct"], known["streams_millions"], s=12, color=ms.CONTEXT,
              alpha=0.6, linewidths=0)
x = np.array([known["danceability_pct"].min(), known["danceability_pct"].max()])
right.plot(x, intercept + slope * x, color=ms.HIGHLIGHT, linewidth=2.2)
right.annotate(f"Regression line\nslope = {slope:.2f} million streams per point\nr = {r:.2f}",
               xy=(0.03, 0.90), xycoords="axes fraction", ha="left", va="top", fontsize=10,
               color=ms.HIGHLIGHT, fontweight="bold")
right.set_ylim(0, known["streams_millions"].max() * 1.05)
right.yaxis.set_major_formatter(lambda v, _: f"{v:,.0f}")
right.set_xlabel("Danceability (0–100)")
right.set_ylabel("Streams (millions)")
right.set_title(f"Every song ({len(known):,})", color=ms.MUTED, fontsize=10)

for ax in (left, right):
    ax.grid(axis="y")
    ax.set_axisbelow(True)

ms.finish([left, right],
          f"More danceable songs do not get more streams: the least danceable quarter leads, r = {r:.2f}",
          f"Streams since release · millions · {len(known):,} songs with a stream count, in four "
          "equal-sized quarters by danceability",
          f"Source: data/songs_2023.csv ({len(songs) - len(known)} of the {len(songs):,} songs have no "
          "stream count and are left out)", tag=ms.tag(2, 17))
ms.save(fig, n=2, slug="streams-by-danceability")
plt.show()

# Reveal: your call on Q1·H1 beside the data.
display(HTML(ms.reveal(CALLS, "Q1·H1", holds=False, verdict="The least danceable quarter leads (366.8 million vs 253.9), and r = −0.11.")))

# Play the songs this test evaluated.
from mixtape import mixtape, span

quarter = pd.qcut(songs["danceability_pct"], 4)
for q, which in [(quarter.cat.categories[0], "Least"), (quarter.cat.categories[-1], "Most")]:
    picks = songs[quarter == q].nlargest(5, "streams_millions")
    display(mixtape(picks, f"{which} danceable quarter ({span(q)}): its 5 most-streamed songs", "danceability_pct"))

# The same computations with valence in place of danceability. Call the medians upbeat.
upbeat, upbeat_n = quarter_medians(songs, "valence_pct", "streams_millions")
display(ms.table(pd.DataFrame({"median streams (millions)": upbeat.round(1), "songs": upbeat_n})))
display(ms.table(slope_and_r(songs, "valence_pct", "streams_millions").to_frame("streams on valence").round(3), index_label="measure"))

# Figure: the test of hypothesis 2, from upbeat and songs.
fig, (left, right) = plt.subplots(1, 2, figsize=(11, 4.5))

# Left: the median streams of each valence quarter, with each bar's median and number of songs.
known = songs.dropna(subset=["streams_millions"])
ranges = [f"{int(np.floor(q.left)) + 1}–{int(q.right)}" for q in upbeat.index]
bars = left.bar(ranges, upbeat.values, color=ms.CONTEXT, width=0.65)
for bar, median, n in zip(bars, upbeat.values, upbeat_n.values):
    middle = bar.get_x() + bar.get_width() / 2
    left.annotate(f"{median:.1f}", xy=(middle, median), xytext=(0, 4), textcoords="offset points",
                  ha="center", va="bottom", fontsize=10, color=ms.INK, fontweight="bold")
    left.annotate(f"n = {n}", xy=(middle, 0), xytext=(0, 6), textcoords="offset points",
                  ha="center", va="bottom", fontsize=9, color=ms.on(ms.CONTEXT))
left.set_ylim(0, upbeat.max() * 1.18)
left.set_xlabel("Valence quarter (0–100)")
left.set_ylabel("Median streams (millions)")
left.set_title("Median streams by quarter", color=ms.MUTED, fontsize=10)

# Right: every song with a stream count, and the regression line of streams on valence.
slope, intercept, r = ms.regression(known["valence_pct"], known["streams_millions"])
right.scatter(known["valence_pct"], known["streams_millions"], s=12, color=ms.CONTEXT,
              alpha=0.6, linewidths=0)
x = np.array([known["valence_pct"].min(), known["valence_pct"].max()])
right.plot(x, intercept + slope * x, color=ms.HIGHLIGHT, linewidth=2.2)
right.annotate(f"Regression line\nslope = {slope:.2f} million streams per point\nr = {r:.2f}",
               xy=(0.03, 0.90), xycoords="axes fraction", ha="left", va="top", fontsize=10,
               color=ms.HIGHLIGHT, fontweight="bold")
right.set_ylim(0, known["streams_millions"].max() * 1.05)
right.yaxis.set_major_formatter(lambda v, _: f"{v:,.0f}")
right.set_xlabel("Valence (0–100, sad to happy)")
right.set_ylabel("Streams (millions)")
right.set_title(f"Every song ({len(known):,})", color=ms.MUTED, fontsize=10)

for ax in (left, right):
    ax.grid(axis="y")
    ax.set_axisbelow(True)

ms.finish([left, right],
          f"More upbeat songs do not get more streams: the most upbeat quarter is lowest, r = {r:.2f}",
          f"Streams since release · millions · {len(known):,} songs with a stream count, in four "
          "equal-sized quarters by valence",
          f"Source: data/songs_2023.csv ({len(songs) - len(known)} of the {len(songs):,} songs have no "
          "stream count and are left out)", tag=ms.tag(3, 17))
ms.save(fig, n=3, slug="streams-by-valence")
plt.show()

# Reveal: your call on Q1·H2 beside the data.
display(HTML(ms.reveal(CALLS, "Q1·H2", holds=False, verdict="The happiest quarter is lowest (261.4 million), and r = −0.04.")))

# Play the songs this test evaluated.
from mixtape import mixtape, span

quarter = pd.qcut(songs["valence_pct"], 4)
for q, which in [(quarter.cat.categories[0], "Least"), (quarter.cat.categories[-1], "Most")]:
    picks = songs[quarter == q].nlargest(5, "streams_millions")
    display(mixtape(picks, f"{which} upbeat quarter ({span(q)}): its 5 most-streamed songs", "valence_pct"))

# Group the songs by release year with pd.cut into a column called released: before 2000, 2000-09,
# 2010-14, 2015-18, 2019-20, 2021, 2022 and 2023. The next section uses the same groups.
# Each edge is the last year of a group: 1930-1999, 2000-2009, ..., then one group per year.
RELEASE_EDGES = [1929, 1999, 2009, 2014, 2018, 2020, 2021, 2022, 2023]
RELEASE_GROUPS = ["before 2000", "2000-09", "2010-14", "2015-18", "2019-20", "2021", "2022", "2023"]
songs["released"] = pd.cut(songs["release_year"], bins=RELEASE_EDGES, labels=RELEASE_GROUPS)
assert songs["released"].notna().sum() == songs["release_year"].notna().sum()  # every known year has a group
display(ms.table(songs["released"].value_counts(sort=False).rename("songs").to_frame()))

# - The median streams of each release group. Call it by_release.
# - A table of each group's number of songs and median streams.
release = songs.groupby("released", observed=True)
by_release = release["streams_millions"].median()
display(ms.table(pd.DataFrame({"songs": release.size(), "median streams (millions)": by_release.round(1)})))

# Figure: median streams by release group, from by_release.
import sys; sys.path.append(".claude/skills/fresh-figures/scripts")
import fresh_style as ms

ms.use_theme()
# The bar panel is wider so its eight year ranges fit side by side.
fig, (left, right) = plt.subplots(1, 2, figsize=(11, 4.5), width_ratios=[1.85, 1])

# Left: the median streams of each release group, with the 2021-2023 releases highlighted and each
# bar's median and number of songs.
known = songs.dropna(subset=["streams_millions"])
per_group = known.groupby("released", observed=True).size()
recent = ["2021", "2022", "2023"]
colors = [ms.HIGHLIGHT if group in recent else ms.CONTEXT for group in by_release.index]
# Each tick names the group's real years (the oldest song is from 1930) and its number of songs.
first_year = int(songs["release_year"].min())
years = [f"{first_year}–99" if g == "before 2000" else g.replace("-", "–") for g in by_release.index]
ticks = [f"{y}\nn = {n}" for y, n in zip(years, per_group.values)]
bars = left.bar(ticks, by_release.values, color=colors, width=0.7)
for bar, median in zip(bars, by_release.values):
    left.annotate(f"{median:,.1f}", xy=(bar.get_x() + bar.get_width() / 2, median), xytext=(0, 4),
                  textcoords="offset points", ha="center", va="bottom", fontsize=9, color=ms.INK,
                  fontweight="bold")
# A bracket over the three newest groups, so the highlight is named in words, not only in color.
first, last = bars[-3], bars[-1]
x0, x1 = first.get_x(), last.get_x() + last.get_width()
top = by_release[recent].max() * 1.3
left.plot([x0, x0, x1, x1], [top * 0.95, top, top, top * 0.95], color=ms.HIGHLIGHT, linewidth=1.2)
left.annotate("released 2021–23", xy=((x0 + x1) / 2, top), xytext=(0, 4), textcoords="offset points",
              ha="center", va="bottom", fontsize=10, color=ms.HIGHLIGHT, fontweight="bold")
left.set_ylim(0, by_release.max() * 1.15)
left.yaxis.set_major_formatter(lambda v, _: f"{v:,.0f}")
left.tick_params(axis="x", labelsize=9)
left.set_xlabel("Release year (n = songs in the group)")
left.set_ylabel("Median streams (millions)")
left.set_title("Median streams by release group", color=ms.MUTED, fontsize=10)

# Right: every song with a stream count, and the regression line of streams on release year.
slope, intercept, r = ms.regression(known["release_year"], known["streams_millions"])
right.scatter(known["release_year"], known["streams_millions"], s=12, color=ms.CONTEXT, alpha=0.6,
              linewidths=0)
x = np.array([known["release_year"].min(), known["release_year"].max()])
right.plot(x, intercept + slope * x, color=ms.HIGHLIGHT, linewidth=2.2)
right.annotate(f"Regression line\nslope = {slope:.2f} million\nstreams per year\nr = {r:.2f}",
               xy=(0.03, 0.90), xycoords="axes fraction", ha="left", va="top", fontsize=10,
               color=ms.HIGHLIGHT, fontweight="bold")
right.set_ylim(0, known["streams_millions"].max() * 1.05)
right.yaxis.set_major_formatter(lambda v, _: f"{v:,.0f}")
right.set_xlabel("Release year")
right.set_ylabel("Streams (millions)")
right.set_title(f"Every song ({len(known):,})", color=ms.MUTED, fontsize=10)

for ax in (left, right):
    ax.grid(axis="y")
    ax.set_axisbelow(True)

ms.finish([left, right],
          f"Newer songs have far fewer streams: the typical 2023 release has "
          f"{by_release['2023']:,.1f} million",
          f"Streams since release · millions · {len(known):,} songs with a stream count, grouped by "
          "release year",
          f"Source: data/songs_2023.csv ({len(songs) - len(known)} of the {len(songs):,} songs have no "
          "stream count and are left out; streams were read once, in mid-2023)", tag=ms.tag(4, 17))
# The two-line ticks push the axis label down, so the caption moves below it.
left.texts[-1].xyann = (0, -58)
ms.save(fig, n=4, slug="streams-by-release-year")
plt.show()

# Play the most-streamed song from each release group, oldest first.
from mixtape import mixtape

picks = songs.sort_values("streams_millions", ascending=False).drop_duplicates("released").sort_values("released")
display(mixtape(picks, "The biggest song of each era on the 2023 list", "release_year", limit=8))

# - The median danceability of each release group. Call it dance_by_release.
# - A table of each group's median danceability and its share of the list.
dance_by_release = release["danceability_pct"].median()
display(ms.table(pd.DataFrame({
    "median danceability (0-100)": dance_by_release,
    "share of the list": (release.size() / len(songs)).round(3),
})))

# Figure: median danceability by release group, from dance_by_release.
import sys; sys.path.append(".claude/skills/fresh-figures/scripts")
import fresh_style as ms

ms.use_theme()
# The bar panel is wider so its eight year ranges fit side by side.
fig, (left, right) = plt.subplots(1, 2, figsize=(11, 4.5), width_ratios=[1.75, 1])

# Left: the median danceability of each release group, with the 2021-2023 releases highlighted and
# each bar's median and number of songs.
known = songs.dropna(subset=["danceability_pct"])
per_group = known.groupby("released", observed=True).size()
recent = ["2021", "2022", "2023"]
colors = [ms.HIGHLIGHT if group in recent else ms.CONTEXT for group in dance_by_release.index]
# Each tick names the group's real years (the oldest song is from 1930) and its number of songs.
first_year = int(songs["release_year"].min())
years = [f"{first_year}–99" if g == "before 2000" else g.replace("-", "–")
         for g in dance_by_release.index]
ticks = [f"{y}\nn = {n}" for y, n in zip(years, per_group.values)]
bars = left.bar(ticks, dance_by_release.values, color=colors, width=0.7)
for bar, median in zip(bars, dance_by_release.values):
    left.annotate(f"{median:.1f}", xy=(bar.get_x() + bar.get_width() / 2, median), xytext=(0, 4),
                  textcoords="offset points", ha="center", va="bottom", fontsize=9, color=ms.INK,
                  fontweight="bold")
# A bracket over the three newest groups, so the highlight is named in words, not only in color.
first, last = bars[-3], bars[-1]
x0, x1 = first.get_x(), last.get_x() + last.get_width()
# Raised clear of the 80 gridline, so the bracket does not merge with it.
top = dance_by_release[recent].max() * 1.17
left.plot([x0, x0, x1, x1], [top * 0.97, top, top, top * 0.97], color=ms.HIGHLIGHT, linewidth=1.2)
left.annotate("released 2021–23", xy=((x0 + x1) / 2, top), xytext=(0, 4), textcoords="offset points",
              ha="center", va="bottom", fontsize=10, color=ms.HIGHLIGHT, fontweight="bold")
left.set_ylim(0, 100)
left.tick_params(axis="x", labelsize=9)
left.set_xlabel("Release year (n = songs in the group)")
left.set_ylabel("Median danceability (0–100)")
left.set_title("Median danceability by release group", color=ms.MUTED, fontsize=10)

# Right: every song, and the regression line of danceability on release year.
slope, intercept, r = ms.regression(known["release_year"], known["danceability_pct"])
right.scatter(known["release_year"], known["danceability_pct"], s=12, color=ms.CONTEXT, alpha=0.6,
              linewidths=0)
x = np.array([known["release_year"].min(), known["release_year"].max()])
right.plot(x, intercept + slope * x, color=ms.HIGHLIGHT, linewidth=2.2)
right.annotate(f"Regression line\nslope = {slope:.2f} points per year\nr = {r:.2f}",
               xy=(0.03, 0.06), xycoords="axes fraction", ha="left", va="bottom", fontsize=10,
               color=ms.HIGHLIGHT, fontweight="bold")
right.set_ylim(0, 100)
right.set_xlabel("Release year")
right.set_ylabel("Danceability (0–100)")
right.set_title(f"Every song ({len(known):,})", color=ms.MUTED, fontsize=10)

for ax in (left, right):
    ax.grid(axis="y")
    ax.set_axisbelow(True)

older = dance_by_release.drop(recent)
ms.finish([left, right],
          f"The newest songs are the most danceable: {dance_by_release[recent].min():.0f}+ for "
          f"2021–23, {older.min():.1f}–{older.max():.1f} before",
          f"Danceability · Spotify's 0–100 score · {len(known):,} songs, grouped by release year",
          "Source: data/songs_2023.csv", tag=ms.tag(5, 17))
# The two-line ticks push the axis label down, so the caption moves below it.
left.texts[-1].xyann = (0, -58)
ms.save(fig, n=5, slug="danceability-by-release-year")
plt.show()

# Play the most-streamed songs released 2021-23, the newest and most danceable group.
from mixtape import mixtape

recent = songs[songs["release_year"] >= 2021].nlargest(5, "streams_millions")
display(mixtape(recent, "Released 2021–23: the 5 most-streamed", "danceability_pct"))

# Hypothesis 1's computations again, on the songs released in 2023 only. Call those songs
# released_2023, and their medians dance_2023.
released_2023 = songs[songs["release_year"] == 2023]
dance_2023, dance_2023_n = quarter_medians(released_2023, "danceability_pct", "streams_millions")
display(ms.table(pd.DataFrame({"median streams (millions)": dance_2023.round(1), "songs": dance_2023_n})))
display(ms.table(slope_and_r(released_2023, "danceability_pct", "streams_millions").to_frame("2023 releases").round(3), index_label="measure"))

# Figure: the test of hypothesis 1 on the 2023 releases, from dance_2023 and
#   released_2023.
import sys; sys.path.append(".claude/skills/fresh-figures/scripts")
import fresh_style as ms

ms.use_theme()
fig, (left, right) = plt.subplots(1, 2, figsize=(11, 4.5))

# Left: the median streams of each danceability quarter, with each bar's median and number of songs.
known = released_2023.dropna(subset=["streams_millions"])
per_quarter = known.groupby(pd.qcut(released_2023["danceability_pct"], 4), observed=True).size()
ranges = [f"{int(np.floor(q.left)) + 1}–{int(q.right)}" for q in dance_2023.index]
bars = left.bar(ranges, dance_2023.values, color=ms.CONTEXT, width=0.65)
for bar, median, n in zip(bars, dance_2023.values, per_quarter.values):
    middle = bar.get_x() + bar.get_width() / 2
    left.annotate(f"{median:.1f}", xy=(middle, median), xytext=(0, 4), textcoords="offset points",
                  ha="center", va="bottom", fontsize=10, color=ms.INK, fontweight="bold")
    left.annotate(f"n = {n}", xy=(middle, 0), xytext=(0, 6), textcoords="offset points",
                  ha="center", va="bottom", fontsize=9, color=ms.on(ms.CONTEXT))
left.set_ylim(0, dance_2023.max() * 1.18)
left.set_xlabel("Danceability quarter (0–100)")
left.set_ylabel("Median streams (millions)")
left.set_title("Median streams by quarter", color=ms.MUTED, fontsize=10)

# Right: every 2023 release with a stream count, and the regression line of streams on danceability.
slope, intercept, r = ms.regression(known["danceability_pct"], known["streams_millions"])
right.scatter(known["danceability_pct"], known["streams_millions"], s=14, color=ms.CONTEXT,
              alpha=0.7, linewidths=0)
x = np.array([known["danceability_pct"].min(), known["danceability_pct"].max()])
right.plot(x, intercept + slope * x, color=ms.HIGHLIGHT, linewidth=2.2)
right.annotate(f"Regression line\nslope = {slope:.2f} million streams per point\nr = {r:.2f}",
               xy=(0.03, 0.90), xycoords="axes fraction", ha="left", va="top", fontsize=10,
               color=ms.HIGHLIGHT, fontweight="bold")
right.set_ylim(0, known["streams_millions"].max() * 1.05)
right.yaxis.set_major_formatter(lambda v, _: f"{v:,.0f}")
right.set_xlabel("Danceability (0–100)")
right.set_ylabel("Streams (millions)")
right.set_title(f"Every 2023 release ({len(known):,})", color=ms.MUTED, fontsize=10)

for ax in (left, right):
    ax.grid(axis="y")
    ax.set_axisbelow(True)

ms.finish([left, right],
          f"Among songs of the same age, danceability still does not bring streams: r = {r:.2f}",
          f"Streams since release · millions · {len(known):,} songs released in 2023, in four "
          "equal-sized quarters by danceability",
          f"Source: data/songs_2023.csv (the {len(known):,} of the {len(songs):,} songs on the list "
          "that were released in 2023)", tag=ms.tag(6, 17))
ms.save(fig, n=6, slug="streams-by-danceability-2023")
plt.show()

# Play the songs this test evaluated: 2023's most-streamed new releases.
from mixtape import mixtape

display(mixtape(released_2023.nlargest(5, "streams_millions"), "Released in 2023: the 5 most-streamed", "danceability_pct"))

# - A matching name for every song and every track: the track name plus the first artist (split
#   the artists on "," or ";"), in lower case, with spaces trimmed.
# - hit_line: the popularity of the catalog's top 10%, the 0.9 quantile of the tracks' popularity.
def match_name(track, artists):
    """Track name plus first artist, lower case, trimmed: the key shared by both files."""
    first_artist = artists.str.split(r"[,;]").str[0]
    return (track.str.strip() + " | " + first_artist.str.strip()).str.lower()


songs["match_name"] = match_name(songs["track"], songs["artists"])
tracks["match_name"] = match_name(tracks["track"], tracks["artists"])
hit_line = tracks["popularity"].quantile(0.9)
display(ms.table(pd.Series({"hit_line: popularity of the catalog's top 10%": hit_line}, name="popularity"),
         index_label="measure"))

# - Give each list song the highest popularity among the tracks with its matching name; the songs
#   that get one are the ones found in the catalog. Call their popularity found.
# - A table: songs on the list, how many have a damaged name (the character "\ufffd"), how many were
#   found, the median popularity of the found songs and of every track, hit_line, and the share of
#   found songs at or above it.
best_popularity = tracks.groupby("match_name")["popularity"].max()
songs["popularity"] = songs["match_name"].map(best_popularity)
found = songs.dropna(subset=["popularity"])
damaged = songs["track"].str.contains("\ufffd") | songs["artists"].str.contains("\ufffd")
damaged_in_file = songs["track_in_file"].str.contains("\ufffd") | songs["artists_in_file"].str.contains("\ufffd")
display(ms.table(pd.Series({
    "songs on the list": len(songs),
    "songs with a damaged name in the file": damaged_in_file.sum(),
    "songs with a damaged name after the approved repairs": damaged.sum(),
    "songs found in the catalog": len(found),
    "not found: same artists in the catalog, in a different order (flagged in the data notebook)":
        (songs["artist_order_differs_from_catalog"] & songs["popularity"].isna()).sum(),
    "median popularity of the found songs": found["popularity"].median(),
    "median popularity of every track": tracks["popularity"].median(),
    "hit_line": hit_line,
    "share of found songs at or above hit_line": round((found["popularity"] >= hit_line).mean(), 3),
}, name="value"), index_label="measure"))

# Figure: the popularity of the list songs found in the catalog compared with
#   every track, from found, tracks and hit_line.
import sys; sys.path.append(".claude/skills/fresh-figures/scripts")
import fresh_style as ms

ms.use_theme()
fig, ax = plt.subplots(figsize=(8, 4.5))
bins = np.arange(0, 101, 5)
# Shares of each group, not counts: the catalog has hundreds of times as many tracks as were found.
share = lambda s: np.ones(len(s)) / len(s)
ax.hist(tracks["popularity"], bins=bins, weights=share(tracks), color=ms.CONTEXT)
ax.hist(found["popularity"], bins=bins, weights=share(found), histtype="step", linewidth=2.2,
        color=ms.HIGHLIGHT)
ax.yaxis.set_major_formatter(lambda v, _: f"{v:.0%}")
ax.grid(axis="y")
ax.set_axisbelow(True)
top = ax.get_ylim()[1] * 1.22
ax.set_ylim(0, top)

# The hit line and each group's median as labeled lines, each label just left of its own line.
above = (found["popularity"] >= hit_line).mean()
ax.axvline(hit_line, color=ms.MUTED, linestyle=":", linewidth=1.6, zorder=3)
ax.annotate(f"Top 10% line: {hit_line:.0f}\n{above:.0%} of found\nsongs at or above",
            xy=(hit_line, top * 0.80), xytext=(-6, 0), textcoords="offset points", ha="right",
            va="top", fontsize=10, color=ms.MUTED, fontweight="bold")
for group, data, color, style in [
    ("Catalog median", tracks, ms.INK, "--"),
    ("Found median", found, ms.HIGHLIGHT, "-"),
]:
    m = data["popularity"].median()
    ax.axvline(m, color=color, linestyle=style, linewidth=1.4, zorder=3)
    ax.annotate(f"{group}: {m:.0f}", xy=(m, top), xytext=(-6, -4), textcoords="offset points",
                ha="right", va="top", fontsize=10, color=color, fontweight="bold")

# Direct labels in place of a legend.
ax.annotate(f"List songs\nfound in\nthe catalog\n{len(found):,} songs", xy=(63, 0.12), ha="left",
            va="bottom", fontsize=10, color=ms.HIGHLIGHT)
ax.annotate(f"Every track in\nthe catalog\n{len(tracks):,} tracks", xy=(2, top * 0.55), ha="left",
            va="top", fontsize=10, color=ms.INK)

ax.set_xlim(0, 100)
ax.set_xlabel("Spotify popularity (0–100)")
ax.set_ylabel("Share of the group (%)")
ms.finish(ax, f"The hits sit far above the catalog: {above:.0%} clear its top-10% line",
          f"Popularity in 5-point bins · share of each group · {len(found):,} found list songs, "
          f"{len(tracks):,} catalog tracks",
          f"Source: data/songs_2023.csv, data/catalog.csv (only the {len(found):,} of {len(songs):,} "
          "list songs matched by track name and first artist)", tag=ms.tag(7, 17))
ms.save(fig, n=7, slug="popularity-found-vs-catalog")
plt.show()

# Play the 5 most popular list songs found in the catalog.
from mixtape import mixtape

display(mixtape(found.nlargest(5, "popularity"), "The list's 5 most popular songs in the catalog", "popularity"))

# Hypothesis 1's computations on every track, with popularity in place of streams. Call the medians
# dance_catalog.
dance_catalog, dance_catalog_n = quarter_medians(tracks, "danceability_pct", "popularity")
display(ms.table(pd.DataFrame({"median popularity (0-100)": dance_catalog, "tracks": dance_catalog_n})))
display(ms.table(slope_and_r(tracks, "danceability_pct", "popularity").rename({"songs": "tracks"})
         .to_frame("popularity on danceability").round(3), index_label="measure"))

# Figure: the test of hypothesis 1 on every track in the catalog, from
#   dance_catalog and tracks.
import sys; sys.path.append(".claude/skills/fresh-figures/scripts")
import fresh_style as ms

ms.use_theme()
fig, (left, right) = plt.subplots(1, 2, figsize=(11, 4.5))

# Left: the median popularity of each danceability quarter, with each bar's median and number of tracks.
known = tracks.dropna(subset=["popularity"])
per_quarter = known.groupby(pd.qcut(tracks["danceability_pct"], 4), observed=True).size()
ranges = [f"{int(np.floor(q.left)) + 1}–{int(q.right)}" for q in dance_catalog.index]
# The lowest quarter starts at the lowest score in the catalog.
ranges[0] = f"{int(tracks['danceability_pct'].min())}–{int(dance_catalog.index[0].right)}"
bars = left.bar(ranges, dance_catalog.values, color=ms.CONTEXT, width=0.65)
for bar, median, n in zip(bars, dance_catalog.values, per_quarter.values):
    middle = bar.get_x() + bar.get_width() / 2
    left.annotate(f"{median:.0f}", xy=(middle, median), xytext=(0, 4), textcoords="offset points",
                  ha="center", va="bottom", fontsize=10, color=ms.INK, fontweight="bold")
    # n on two lines, so five-digit counts leave room between neighboring bars.
    left.annotate(f"{n:,}\ntracks", xy=(middle, 0), xytext=(0, 6), textcoords="offset points",
                  ha="center", va="bottom", fontsize=9, color=ms.on(ms.CONTEXT))
left.set_ylim(0, 100)
left.set_xlabel("Danceability quarter (0–100)")
left.set_ylabel("Median popularity (0–100)")
left.set_title("Median popularity by quarter", color=ms.MUTED, fontsize=10)

# Right: every track, and the regression line of popularity on danceability. The dots are small and
# faint so the dense middle still shows where most tracks sit.
slope, intercept, r = ms.regression(known["danceability_pct"], known["popularity"])
right.scatter(known["danceability_pct"], known["popularity"], s=3, color=ms.CONTEXT, alpha=0.15,
              linewidths=0, rasterized=True)
x = np.array([known["danceability_pct"].min(), known["danceability_pct"].max()])
right.plot(x, intercept + slope * x, color=ms.HIGHLIGHT, linewidth=2.4)
right.annotate(f"Regression line\nslope = {slope:.2f} points\nper danceability point\nr = {r:.2f}",
               xy=(0.03, 0.97), xycoords="axes fraction", ha="left", va="top", fontsize=10,
               color=ms.HIGHLIGHT, fontweight="bold",
               bbox=dict(facecolor=ms.BG, edgecolor="none", alpha=0.9, pad=3))
right.set_ylim(0, 100)
right.set_xlim(0, 100)
right.set_xlabel("Danceability (0–100)")
right.set_ylabel("Popularity (0–100)")
right.set_title(f"Every track ({len(known):,})", color=ms.MUTED, fontsize=10)

for ax in (left, right):
    ax.grid(axis="y")
    ax.set_axisbelow(True)

ms.finish([left, right],
          f"Across the whole catalog, danceability barely moves popularity: r = {r:.2f}",
          f"Spotify popularity · 0–100 · {len(known):,} catalog tracks whose sound was measured, in four "
          "equal-sized quarters by danceability",
          f"Source: data/catalog.csv (one row per track; the {catalog['danceability_pct'].isna().sum()} "
          "rows with no rhythm analysis are left out)", tag=ms.tag(8, 17))
ms.save(fig, n=8, slug="popularity-by-danceability-catalog")
plt.show()

# Play the most popular tracks in the least and the most danceable quarters of the catalog.
from mixtape import mixtape, span

quarter = pd.qcut(tracks["danceability_pct"], 4)
for q, which in [(quarter.cat.categories[0], "Least"), (quarter.cat.categories[-1], "Most")]:
    picks = tracks[quarter == q].nlargest(5, "popularity")
    display(mixtape(picks, f"{which} danceable quarter ({span(q)}): its 5 most popular tracks", "danceability_pct"))

# - The mean popularity of each genre, from every catalog row, from least to most popular. Call it
#   by_genre.
# - Show the five least and five most popular genres.
by_genre = catalog.groupby("genre")["popularity"].mean().sort_values()
ends = pd.concat([by_genre.head(5), by_genre.tail(5)]).round(1)
display(ms.table(ends.rename("mean popularity (0-100)").to_frame()
         .assign(end=["least popular"] * 5 + ["most popular"] * 5), index_label="genre"))

# To compare like with like, the danceability quarters' popularity as means too (the test above used
# medians, as planned).
dance_catalog_means = tracks.groupby(pd.qcut(tracks["danceability_pct"], 4), observed=True)["popularity"].mean()
display(ms.table(dance_catalog_means.round(1).rename("mean popularity (0-100)").to_frame(),
                 index_label="danceability quarter"))

# Figure: the mean popularity of every genre, from by_genre.
import sys; sys.path.append(".claude/skills/fresh-figures/scripts")
import fresh_style as ms

ms.use_theme()
# One named bar per genre: horizontal, and tall enough (0.16 inch a bar) that every name fits.
fig, ax = plt.subplots(figsize=(8, 0.16 * len(by_genre) + 1.6))

# Every genre from least popular (bottom) to most popular (top), with the five at each end (the ten
# in the table above) highlighted, and each bar's mean printed at its tip.
ends = list(by_genre.index[:5]) + list(by_genre.index[-5:])
colors = [ms.HIGHLIGHT if genre in ends else ms.CONTEXT for genre in by_genre.index]
rows = np.arange(len(by_genre))
bars = ax.barh(rows, by_genre.values, color=colors, height=0.72)
for bar, genre, mean in zip(bars, by_genre.index, by_genre.values):
    ax.annotate(f"{mean:.1f}", xy=(mean, bar.get_y() + bar.get_height() / 2), xytext=(3, 0),
                textcoords="offset points", ha="left", va="center", fontsize=9,
                color=ms.HIGHLIGHT if genre in ends else ms.MUTED,
                fontweight="bold" if genre in ends else "normal", zorder=4,
                # A background patch, so the dashed mean line never runs through a number.
                bbox=dict(boxstyle="square,pad=0.05", facecolor=ms.BG, edgecolor="none"))
ax.set_yticks(rows, by_genre.index, fontsize=9)
ax.set_ylim(-0.8, len(by_genre) - 0.2)

# The range of the four danceability quarters' means, shaded behind the bars, so the comparison the
# section asks about (does genre spread popularity more than danceability?) is on the chart.
low, high = dance_catalog_means.min(), dance_catalog_means.max()
band = ms.palette(2)[1]  # the theme's second color, so the band is not read as a highlight
ax.axvspan(low, high, color=band, alpha=0.2, zorder=0.5, linewidth=0)
ax.annotate(f"Danceability quarters' means:\n{low:.1f} to {high:.1f}", xy=(high, 10), xytext=(6, 0),
            textcoords="offset points", ha="left", va="center", fontsize=10, color=band,
            fontweight="bold")

# The mean over every catalog row as a labeled reference line.
overall = catalog["popularity"].mean()
ax.axvline(overall, color=ms.INK, linestyle="--", linewidth=1.2, zorder=3)
ax.annotate(f"Mean over every row: {overall:.1f}", xy=(overall, 20), xytext=(6, 0),
            textcoords="offset points", ha="left", va="center", fontsize=10, color=ms.INK,
            fontweight="bold")

# The highlight named in words beside each end.
ax.annotate("Five most\npopular", xy=(by_genre.iloc[-1], len(by_genre) - 3), xytext=(32, 0),
            textcoords="offset points", ha="left", va="center", fontsize=10, color=ms.HIGHLIGHT,
            fontweight="bold")
# On two lines, so it ends well short of the dashed mean line.
ax.annotate("Five least\npopular", xy=(by_genre.iloc[4], 2), xytext=(28, 0),
            textcoords="offset points", ha="left", va="center", fontsize=10, color=ms.HIGHLIGHT,
            fontweight="bold")

ax.set_xlim(0, 75)
ax.set_xlabel("Mean popularity (0–100)")
ax.grid(axis="x")
ax.grid(axis="y", visible=False)
ax.set_axisbelow(True)

per_genre = catalog.groupby("genre").size()
ms.finish(ax,
          f"Genre splits popularity wide: means run from {by_genre.min():.0f} to {by_genre.max():.0f}",
          f"Mean Spotify popularity of each of {len(by_genre)} genres · 0–100 · {len(catalog):,} "
          f"catalog rows, {per_genre.min():,}–{per_genre.max():,} per genre",
          "Source: data/catalog.csv (one row per track and genre search; a track found by two "
          "searches counts in both genres)", tag=ms.tag(9, 17))
ms.save(fig, n=9, slug="popularity-by-genre")
plt.show()

# Play the most popular track of each of the five most and five least popular genres.
from mixtape import mixtape

for genres, which in [(by_genre.index[::-1][:5], "most"), (by_genre.index[:5], "least")]:
    picks = pd.concat([catalog[catalog["genre"] == g].nlargest(1, "popularity") for g in genres])
    display(mixtape(picks, f"The top track of each of the five {which} popular genres", "genre"))

# - The six sound features: danceability, energy, acousticness, valence, speechiness and
#   instrumentalness. Standardize them with StandardScaler.
# - Fit k-means for k = 2 to 8 (n_init=10, random_state=6120). For each k, keep the within-cluster
#   spread (the fit's inertia_) and the silhouette score, computed on a random sample of 10,000
#   tracks (np.random.default_rng(6120)).
# - Show both as a table indexed by k. Call it choosing_k.
SOUND = ["danceability_pct", "energy_pct", "acousticness_pct", "valence_pct", "speechiness_pct",
         "instrumentalness_pct"]
sound = StandardScaler().fit_transform(tracks[SOUND])  # each feature in standard deviations from its mean
sample = np.random.default_rng(SEED).choice(len(sound), 10_000, replace=False)

fits, rows = {}, {}
for k in range(2, 9):
    fits[k] = KMeans(n_clusters=k, n_init=10, random_state=SEED).fit(sound)
    rows[k] = {"within-cluster spread": fits[k].inertia_,
               "silhouette": silhouette_score(sound[sample], fits[k].labels_[sample])}
choosing_k = pd.DataFrame(rows).T.rename_axis("k")
display(ms.table(choosing_k.round({"within-cluster spread": 0, "silhouette": 3}), index_label="k"))

# Figure: the within-cluster spread and the silhouette score for each k, from
#   choosing_k.
fig, (left, right) = plt.subplots(1, 2, figsize=(11, 4.5), sharex=True)
ks = choosing_k.index.to_numpy()
spread, silhouette = choosing_k["within-cluster spread"], choosing_k["silhouette"]
picked = 4  # the fit the next section uses, highlighted in both panels

# Left: the within-cluster spread at each k, in thousands, with each point's value printed above
# and to the right, clear of the falling line (the last one centered below, so it stays in the panel).
left.plot(ks, spread / 1_000, color=ms.CONTEXT, linewidth=2, marker="o", markersize=7, zorder=3)
left.plot(picked, spread[picked] / 1_000, "o", color=ms.HIGHLIGHT, markersize=11, zorder=4)
for k, value in spread.items():
    last = k == ks.max()
    left.annotate(f"{value:,.0f}", xy=(k, value / 1_000), xytext=(0, -12) if last else (7, 6),
                  textcoords="offset points", ha="center" if last else "left",
                  va="top" if last else "bottom", fontsize=9,
                  color=ms.HIGHLIGHT if k == picked else ms.MUTED,
                  fontweight="bold" if k == picked else "normal")
left.set_ylim(0, spread.max() / 1_000 * 1.18)
left.set_ylabel("Within-cluster spread (thousands)")
left.set_title("Spread falls steadily, with no sharp elbow", color=ms.MUTED, fontsize=10)
# The highlight named in words, in the open space under the line.
left.annotate(f"k = {picked}, the fit used next", xy=(picked, spread[picked] / 1_000),
              xytext=(0, -58), textcoords="offset points", ha="center", va="top", fontsize=10,
              color=ms.HIGHLIGHT, fontweight="bold",
              arrowprops=dict(arrowstyle="-", color=ms.HIGHLIGHT, linewidth=1))

# Right: the silhouette score at each k, on an axis tall enough that its low level shows, with a
# labeled reference line at 0.5, a common bar for clearly separate groups.
right.plot(ks, silhouette, color=ms.CONTEXT, linewidth=2, marker="o", markersize=7, zorder=3)
right.plot(picked, silhouette[picked], "o", color=ms.HIGHLIGHT, markersize=11, zorder=4)
for k, value in silhouette.items():
    right.annotate(f"{value:.3f}", xy=(k, value), xytext=(0, 9), textcoords="offset points",
                   ha="center", va="bottom", fontsize=9,
                   color=ms.HIGHLIGHT if k == picked else ms.MUTED,
                   fontweight="bold" if k == picked else "normal")
right.axhline(0.5, color=ms.INK, linestyle="--", linewidth=1.2)
right.annotate("0.5: a common bar for clearly separate groups", xy=(ks.max(), 0.5), xytext=(0, 5),
               textcoords="offset points", ha="right", va="bottom", fontsize=10, color=ms.INK,
               fontweight="bold")
right.set_ylim(0, 0.6)
right.set_ylabel("Silhouette score (-1 to 1)")
right.set_title("Silhouette stays near 0.25 at every k", color=ms.MUTED, fontsize=10)

for ax in (left, right):
    ax.set_xticks(ks)
    ax.set_xlim(ks.min() - 0.5, ks.max() + 0.7)
    ax.set_xlabel("Number of groups, k")
    ax.grid(axis="y")
    ax.grid(axis="x", visible=False)
    ax.set_axisbelow(True)

ms.finish([left, right],
          "Tracks do not fall into clear sound groups at any k",
          f"K-means on six standardized sound features · k = 2 to 8 · {len(tracks):,} tracks "
          f"(silhouette on a sample of {len(sample):,})",
          "Source: data/catalog.csv, one row per track (duplicates across genre searches dropped)",
          tag=ms.tag(10, 17))
ms.save(fig, n=10, slug="choosing-k")
plt.show()

# - Label each track with its group from the k = 4 fit, numbered 1 to 4.
# - For each group: the median of each sound feature and of popularity, and its number of tracks.
#   Call it profile, and show it as a table.
tracks["sound_group"] = fits[4].labels_ + 1
by_group = tracks.groupby("sound_group")
profile = by_group[SOUND + ["popularity"]].median().round(1).assign(tracks=by_group.size())
display(ms.table(profile, index_label="sound group"))

# Figure: the median popularity of each sound group, from profile.
import sys; sys.path.append(".claude/skills/fresh-figures/scripts")
import fresh_style as ms

ms.use_theme()
fig, ax = plt.subplots(figsize=(8, 4.5))

# Each group named by its number and by what sets it apart: how far its median on each feature sits
# from the median of the four groups, in standard deviations over every track, so features with
# different spreads compare fairly. Up to two features more than half a standard deviation above,
# and the one more than half below.
WORDS = {"danceability_pct": "danceable", "energy_pct": "energetic", "acousticness_pct": "acoustic",
         "valence_pct": "upbeat", "speechiness_pct": "wordy", "instrumentalness_pct": "instrumental"}
apart = (profile[SOUND] - profile[SOUND].median()) / tracks[SOUND].std()
sounds = []
for g in profile.index:
    highs = [WORDS[f] for f in apart.loc[g].nlargest(2).index if apart.loc[g, f] > 0.5]
    lows = [f"not {WORDS[f]}" for f in apart.loc[g].nsmallest(1).index if apart.loc[g, f] < -0.5]
    sounds.append(",\n".join(highs + lows))
ticks = [f"Group {g}\n{s}" for g, s in zip(profile.index, sounds)]  # one word a line, so neighbors keep apart

# One series, so one color: no group is singled out before the data speak.
overall = tracks["popularity"].median()
bars = ax.bar(ticks, profile["popularity"], color=ms.CONTEXT, width=0.62)
for bar, median, n in zip(bars, profile["popularity"], profile["tracks"]):
    middle = bar.get_x() + bar.get_width() / 2
    # A bar whose top sits just under the dashed median line gets its number inside the bar, so the
    # line never runs into it; every other number sits on top of its bar.
    inside = overall - 6 < median < overall
    ax.annotate(f"{median:.0f}", xy=(middle, median), xytext=(0, -4 if inside else 4),
                textcoords="offset points", ha="center", va="top" if inside else "bottom",
                fontsize=11, color=ms.on(ms.CONTEXT) if inside else ms.INK, fontweight="bold",
                zorder=4)
    ax.annotate(f"{n:,}\ntracks", xy=(middle, 0), xytext=(0, 6), textcoords="offset points",
                ha="center", va="bottom", fontsize=9, color=ms.on(ms.CONTEXT))

# The median over every track as a labeled reference line, its label in open space past the last bar.
ax.axhline(overall, color=ms.INK, linestyle="--", linewidth=1.2, zorder=3)
ax.annotate(f"Median of\nevery track: {overall:.0f}", xy=(len(profile) - 0.55, overall), xytext=(0, 0),
            textcoords="offset points", ha="left", va="center", fontsize=10, color=ms.INK,
            fontweight="bold", bbox=dict(boxstyle="square,pad=0.2", facecolor=ms.BG, edgecolor="none"))

ax.set_xlim(-0.55, len(profile) + 0.45)
ax.set_ylim(0, 100)
ax.set_xlabel("Sound group from k-means, k = 4 (what sets its sound apart)")
ax.set_ylabel("Median popularity (0–100)")
ax.grid(axis="y")
ax.grid(axis="x", visible=False)
ax.set_axisbelow(True)

low, high = profile["popularity"].min(), profile["popularity"].max()
ms.finish(ax, f"How a track sounds moves its popularity only a little: medians run {low:.0f} to {high:.0f}",
          f"Median Spotify popularity of each sound group · 0–100 · {profile['tracks'].sum():,} tracks",
          "Source: data/catalog.csv, one row per track (duplicates across genre searches dropped); "
          "groups from k-means on six standardized sound features", tag=ms.tag(11, 17))
# The three-line ticks push the axis label down, so the caption moves below it.
ax.texts[-1].xyann = (0, -70)
ms.save(fig, n=11, slug="popularity-by-sound-group")
plt.show()

# For each genre, the share of its tracks in its most common sound group (pd.crosstab with
# normalize="index", then the largest share in each row). Show the three lowest, the three highest,
# and the median genre.
genre_rows = catalog.merge(tracks[["track_id", "sound_group"]], on="track_id")  # every row: one genre search
top_share = pd.crosstab(genre_rows["genre"], genre_rows["sound_group"], normalize="index").max(axis=1).sort_values()
middle = top_share.index[len(top_share) // 2]
display(ms.table(pd.concat([top_share.head(3), top_share[[middle]], top_share.tail(3)]).round(3)
                 .rename("share in its most common sound group").to_frame()
                 .assign(where=["least concentrated"] * 3 + ["median genre"] + ["most concentrated"] * 3),
                 index_label="genre"))

# Hear the four sound groups: the most popular track in each.
from mixtape import mixtape

for group, picks in tracks.sort_values("popularity", ascending=False).groupby("sound_group").head(2).groupby("sound_group"):
    display(mixtape(picks, f"Sound group {group}: its 2 most popular tracks", "popularity", limit=2))

# Reveal: your call on Q1·H3 beside the data.
display(HTML(ms.reveal(CALLS, "Q1·H3", holds=True, verdict="Age stands in for sound: 2000–09 songs have 1,229M streams, 2023 releases 108M; with age held fixed, r = 0.01.")))

# - A hit column: 1 when a track's popularity is at or above hit_line, 0 otherwise.
# - Split tracks into train (70%) and test (30%) with train_test_split (random_state=6120,
#   stratified on hit). Show the hit line, the size of each set, and the share of test tracks that
#   are hits.
N_PICKS = 15  # a label can promote a handful of songs: each model is judged by its 15 top picks

tracks["hit"] = (tracks["popularity"] >= hit_line).astype(int)
train, test = train_test_split(tracks, test_size=0.3, random_state=SEED, stratify=tracks["hit"])
display(ms.table(pd.Series({
    "hit_line (popularity)": hit_line,
    "training tracks (70%)": len(train),
    "test tracks (30%)": len(test),
    "share of test tracks that are hits": round(test["hit"].mean(), 4),
}, name="value"), index_label="measure"))

# - 10,000 times, pick 15 test tracks at random without replacement and count the hits among them
#   (use np.random.default_rng(6120)). Call the counts random_picks.
# - Show the average count, the chance of 5 or more, and the chance of 8 or more.
rng = np.random.default_rng(SEED)
test_hits = test["hit"].to_numpy()
random_picks = np.array([test_hits[rng.choice(len(test_hits), N_PICKS, replace=False)].sum()
                         for _ in range(10_000)])
assert round(random_picks.mean(), 2) == 1.65  # the luck yardstick
display(ms.table(pd.Series({
    "average hits among 15 random picks": round(random_picks.mean(), 2),
    "chance of 5 or more": round((random_picks >= 5).mean(), 4),
    "chance of 8 or more": round((random_picks >= 8).mean(), 4),
}, name="value"), index_label="measure"))

# Figure: how often each number of hits came up in the 10,000 random picks, from
#   random_picks.
import sys; sys.path.append(".claude/skills/fresh-figures/scripts")
import fresh_style as ms

ms.use_theme()
fig, ax = plt.subplots(figsize=(8, 4.5))

# One bar for every possible count, 0 to 15 hits, so the empty right side shows as empty.
counts = np.bincount(random_picks, minlength=N_PICKS + 1)
hits = np.arange(N_PICKS + 1)
ax.bar(hits, counts, color=[ms.HIGHLIGHT if k >= 8 else ms.CONTEXT for k in hits], width=0.6)
top = counts.max() * 1.3
ax.set_ylim(0, top)
mean = random_picks.mean()
for k, c in zip(hits, counts):
    if c:
        # A label next to the dashed average line steps away from it, so the two never touch.
        nudge = 0 if abs(k - mean) > 0.6 else (10 if k > mean else -10)
        ax.annotate(f"{c:,}", xy=(k, c), xytext=(nudge, 4), textcoords="offset points", ha="center",
                    va="bottom", fontsize=9, color=ms.INK, fontweight="bold", zorder=4,
                    # A background patch, so the dashed average line breaks around a number.
                    bbox=dict(boxstyle="square,pad=0.3", facecolor=ms.BG if k < 8 else "none",
                              edgecolor="none"))

# The zone a model must reach (8 or more hits), shaded and named with its chance by luck.
eight = (random_picks >= 8).mean()
lucky = int((random_picks >= 8).sum())
ax.axvspan(7.5, N_PICKS + 0.5, color=ms.GRID, alpha=0.35, zorder=0, linewidth=0)
ax.annotate(f"8 or more hits: a model that\nsingles out the hits\nBy luck: {lucky} of 10,000 picks ({eight:.2%})",
            xy=(11.5, top * 0.55), ha="center", va="center", fontsize=10, color=ms.HIGHLIGHT,
            fontweight="bold")

# The average count as a labeled reference line.
ax.axvline(mean, color=ms.MUTED, linestyle="--", linewidth=1.4, zorder=3)
# The label sits between gridlines, so none runs through it.
ax.annotate(f"Average: {mean:.2f} hits", xy=(mean, counts.max() * 1.17), xytext=(6, 0),
            textcoords="offset points", ha="left", va="center", fontsize=10, color=ms.MUTED, fontweight="bold")

ax.set_xlim(-0.6, N_PICKS + 0.5)
ax.set_xticks(hits)
ax.set_xlabel(f"Hits among {N_PICKS} randomly picked test tracks")
ax.set_ylabel("Number of picks (of 10,000)")
ax.yaxis.set_major_formatter(lambda v, _: f"{v:,.0f}")
ax.grid(axis="y")
ax.grid(axis="x", visible=False)
ax.set_axisbelow(True)

ms.finish(ax, f"By luck, 15 picks reach 8 hits only {lucky} times in 10,000",
          f"Hits among {N_PICKS} test tracks picked at random, 10,000 times · count of picks · "
          f"{len(test):,} test tracks, {test['hit'].mean():.0%} hits",
          "Source: data/catalog.csv, test set (30% of tracks); a hit is a track in the top 10% by "
          "popularity; picks without replacement, seed 6120", tag=ms.tag(12, 17))
ms.save(fig, n=12, slug="random-picks")
plt.show()

# - The list songs with a stream count. Call them listed.
# - The slope (np.polyfit) and r of streams on Spotify playlists, as a small table.
listed = songs.dropna(subset=["streams_millions"])
display(ms.table(slope_and_r(listed, "spotify_playlists", "streams_millions")
                 .to_frame("streams on Spotify playlists").round(4), index_label="measure"))

# Figure: streams against Spotify playlists, from listed.
import sys; sys.path.append(".claude/skills/fresh-figures/scripts")
import fresh_style as ms

ms.use_theme()
fig, ax = plt.subplots(figsize=(8, 4.5))

# Every song with a stream count, and the regression line of streams on playlists.
slope, intercept, r = ms.regression(listed["spotify_playlists"], listed["streams_millions"])
ax.scatter(listed["spotify_playlists"] / 1_000, listed["streams_millions"], s=12, color=ms.CONTEXT,
           alpha=0.6, linewidths=0)
x = np.array([listed["spotify_playlists"].min(), listed["spotify_playlists"].max()])
ax.plot(x / 1_000, intercept + slope * x, color=ms.HIGHLIGHT, linewidth=2.2)
# The slope and r in the empty top-left corner, clear of the dots.
ax.annotate(f"Regression line\nslope = {slope * 1_000:.1f} million streams\nper 1,000 playlists\n"
            f"r = {r:.2f}", xy=(0.03, 0.97), xycoords="axes fraction", ha="left", va="top",
            fontsize=10, color=ms.HIGHLIGHT, fontweight="bold")

# Two real songs named: the most streamed, and the one on the most playlists.
# The second label hangs below its dot, clear of the dots above it.
for song, offset, va in [(listed.loc[listed["streams_millions"].idxmax()], (-8, 8), "center"),
                         (listed.loc[listed["spotify_playlists"].idxmax()], (-4, -4), "top")]:
    ax.annotate(f"{song['track'].split(' - ')[0]}\n{song['spotify_playlists']:,} playlists, "
                f"{song['streams_millions']:,.0f}M streams",
                xy=(song["spotify_playlists"] / 1_000, song["streams_millions"]), xytext=offset,
                textcoords="offset points", ha="right", va=va, fontsize=9, color=ms.INK)

ax.set_xlim(0, listed["spotify_playlists"].max() / 1_000 * 1.04)
ax.set_ylim(0, listed["streams_millions"].max() * 1.08)
ax.xaxis.set_major_formatter(lambda v, _: f"{v:,.0f}")
ax.yaxis.set_major_formatter(lambda v, _: f"{v:,.0f}")
ax.set_xlabel("Spotify playlists the song is on (thousands)")
ax.set_ylabel("Streams (millions)")
ax.grid(axis="y")
ax.grid(axis="x", visible=False)
ax.set_axisbelow(True)

ms.finish(ax, f"Songs on more playlists have far more streams, r = {r:.2f}",
          f"Streams since release · millions · against Spotify playlists · {len(listed):,} songs "
          "with a stream count",
          f"Source: data/songs_2023.csv ({len(songs) - len(listed)} of the {len(songs):,} songs have "
          "no stream count and are left out)", tag=ms.tag(13, 17))
ms.save(fig, n=13, slug="streams-by-playlists")
plt.show()

# Play the most-playlisted songs on the list.
from mixtape import mixtape

display(mixtape(listed.nlargest(5, "spotify_playlists"), "The list's 5 most-playlisted songs", "spotify_playlists"))

# Write a function fit_and_score(title, numeric, categorical) that:
# - fits a logistic regression of hit on train, in a pipeline that standardizes the numeric columns
#   and one-hot encodes the categorical ones (make_column_transformer, make_pipeline);
# - scores every test track (predict_proba(test)[:, 1]);
# - returns (1) a Series named with title, holding the AUC on the test set (roc_auc_score) and the
#   number of hits among the 15 test tracks with the highest scores, and (2) the ROC curve as a
#   DataFrame of the share of non-hits flagged and the share of hits flagged (roc_curve).
# Also make audio, the ten audio features: danceability, energy, valence, acousticness,
# instrumentalness, liveness and speechiness (the _pct columns), loudness_db, tempo_bpm and duration_s.
def fit_and_score(title, numeric, categorical):
    """Fit a logistic regression of hit on train; score every test track; return its scores and ROC curve."""
    parts = [(StandardScaler(), numeric)] + ([(OneHotEncoder(handle_unknown="ignore"), categorical)] if categorical else [])
    model = make_pipeline(make_column_transformer(*parts), LogisticRegression(max_iter=1000))
    model.fit(train, train["hit"])
    score = model.predict_proba(test)[:, 1]  # each test track's estimated chance of being a hit
    top_picks = test.assign(score=score).nlargest(N_PICKS, "score")
    false_share, true_share, _ = roc_curve(test["hit"], score)
    scores = pd.Series({"AUC": roc_auc_score(test["hit"], score),
                        f"hits among its {N_PICKS} top picks": top_picks["hit"].sum()}, name=title)
    roc = pd.DataFrame({"share of non-hits flagged": false_share, "share of hits flagged": true_share})
    return scores, roc


audio = ["danceability_pct", "energy_pct", "valence_pct", "acousticness_pct", "instrumentalness_pct",
         "liveness_pct", "speechiness_pct", "loudness_db", "tempo_bpm", "duration_s"]

# Score a model with danceability as its only feature. Call its results dance_scores and
# dance_roc.
dance_scores, dance_roc = fit_and_score("danceability alone", ["danceability_pct"], [])
display(ms.table(dance_scores.round(3).to_frame(), index_label="measure"))

# Figure: the ROC curve of the danceability model, from dance_roc and dance_scores.
import sys; sys.path.append(".claude/skills/fresh-figures/scripts")
import fresh_style as ms

ms.use_theme()
fig, ax = plt.subplots(figsize=(8, 4.5))
auc = dance_scores["AUC"]
top_hits = int(dance_scores[f"hits among its {N_PICKS} top picks"])
are_hits = "is a hit" if top_hits == 1 else "are hits"

# The coin-flip diagonal as context, the model's curve as the highlight.
ax.plot([0, 1], [0, 1], color=ms.CONTEXT, linestyle="--", linewidth=1.4)
ax.annotate("Coin flip\nAUC = 0.50", xy=(0.80, 0.80), xytext=(10, -6), textcoords="offset points",
            ha="left", va="top", fontsize=10, color=ms.MUTED)
ax.plot(dance_roc["share of non-hits flagged"], dance_roc["share of hits flagged"],
        color=ms.HIGHLIGHT, linewidth=2.2)
# The deciding numbers in the empty top-left corner, clear of both lines.
ax.annotate(f"Danceability model\nAUC = {auc:.3f}\n{top_hits} of its {N_PICKS} top picks {are_hits}\n"
            "(8 or more would single out the hits)", xy=(0.03, 0.97), xycoords="axes fraction",
            ha="left", va="top", fontsize=10, color=ms.HIGHLIGHT, fontweight="bold")

ax.set_xlim(0, 1)
ax.set_ylim(0, 1.02)
ax.xaxis.set_major_formatter(lambda v, _: f"{v:.0%}")
ax.yaxis.set_major_formatter(lambda v, _: f"{v:.0%}")
ax.set_xlabel("Share of all non-hits flagged so far (%)")
ax.set_ylabel("Share of all hits flagged so far (%)")
ax.grid(axis="y")
ax.grid(axis="x", visible=False)
ax.set_axisbelow(True)

ms.finish(ax, f"Danceability alone barely beats a coin flip at spotting hits, AUC = {auc:.2f}",
          f"ROC curve · test tracks flagged from the highest score down · {len(test):,} tracks, "
          f"{int(test['hit'].sum()):,} of them hits",
          "Source: data/catalog.csv (test set: the 30% of tracks the model never saw in training)",
          tag=ms.tag(14, 17))
ms.save(fig, n=14, slug="danceability-roc")
plt.show()

# Reveal: your call on Q2·H1 beside the data.
display(HTML(ms.reveal(CALLS, "Q2·H1", holds=False, verdict="AUC 0.552, and 1 of its 15 top picks is a hit (luck: 1.65; the bar: 8).")))

# Score a model with all ten audio features. Call its results audio_scores and audio_roc.
audio_scores, audio_roc = fit_and_score("all ten audio features", audio, [])
display(ms.table(audio_scores.round(3).to_frame(), index_label="measure"))

# Figure: the ROC curve of the audio-features model, from audio_roc and audio_scores.
import sys; sys.path.append(".claude/skills/fresh-figures/scripts")
import fresh_style as ms

ms.use_theme()
fig, ax = plt.subplots(figsize=(8, 4.5))
auc = audio_scores["AUC"]
top_hits = int(audio_scores[f"hits among its {N_PICKS} top picks"])
are_hits = "is a hit" if top_hits == 1 else "are hits"

# The coin-flip diagonal as context, the model's curve as the highlight.
ax.plot([0, 1], [0, 1], color=ms.CONTEXT, linestyle="--", linewidth=1.4)
ax.annotate("Coin flip\nAUC = 0.50", xy=(0.80, 0.80), xytext=(10, -6), textcoords="offset points",
            ha="left", va="top", fontsize=10, color=ms.MUTED)
ax.plot(audio_roc["share of non-hits flagged"], audio_roc["share of hits flagged"],
        color=ms.HIGHLIGHT, linewidth=2.2)
# The deciding numbers in the empty top-left corner, clear of both lines.
ax.annotate(f"Audio-features model\nAUC = {auc:.3f}\n{top_hits} of its {N_PICKS} top picks {are_hits}\n"
            "(8 or more would single out the hits)", xy=(0.03, 0.97), xycoords="axes fraction",
            ha="left", va="top", fontsize=10, color=ms.HIGHLIGHT, fontweight="bold")

ax.set_xlim(0, 1)
ax.set_ylim(0, 1.02)
ax.xaxis.set_major_formatter(lambda v, _: f"{v:.0%}")
ax.yaxis.set_major_formatter(lambda v, _: f"{v:.0%}")
ax.set_xlabel("Share of all non-hits flagged so far (%)")
ax.set_ylabel("Share of all hits flagged so far (%)")
ax.grid(axis="y")
ax.grid(axis="x", visible=False)
ax.set_axisbelow(True)

ms.finish(ax, f"All ten audio features rank hits better than a coin flip, AUC = {auc:.2f}",
          f"ROC curve · test tracks flagged from the highest score down · {len(test):,} tracks, "
          f"{int(test['hit'].sum()):,} of them hits",
          "Source: data/catalog.csv (test set: the 30% of tracks the model never saw in training)",
          tag=ms.tag(15, 17))
ms.save(fig, n=15, slug="audio-features-roc")
plt.show()

# Reveal: your call on Q2·H2 beside the data.
display(HTML(ms.reveal(CALLS, "Q2·H2", holds=False, verdict="AUC 0.644, and 3 of its 15 top picks are hits (the bar: 8).")))

# Score a model with the ten audio features plus genre, as a category. Call its results
# genre_scores and genre_roc.
genre_scores, genre_roc = fit_and_score("audio features and genre", audio, ["genre"])
display(ms.table(pd.concat([dance_scores, audio_scores, genre_scores], axis=1).round(3), index_label="measure"))

# Figure: the ROC curve of the audio-features-and-genre model, from genre_roc and genre_scores.
import sys; sys.path.append(".claude/skills/fresh-figures/scripts")
import fresh_style as ms

ms.use_theme()
fig, ax = plt.subplots(figsize=(8, 4.5))
auc = genre_scores["AUC"]
top_hits = int(genre_scores[f"hits among its {N_PICKS} top picks"])
are_hits = "is a hit" if top_hits == 1 else "are hits"

# The coin-flip diagonal as context, the model's curve as the highlight.
ax.plot([0, 1], [0, 1], color=ms.CONTEXT, linestyle="--", linewidth=1.4)
ax.annotate("Coin flip\nAUC = 0.50", xy=(0.80, 0.80), xytext=(10, -6), textcoords="offset points",
            ha="left", va="top", fontsize=10, color=ms.MUTED)
ax.plot(genre_roc["share of non-hits flagged"], genre_roc["share of hits flagged"],
        color=ms.HIGHLIGHT, linewidth=2.2)
# The deciding numbers in the empty space below the diagonal, between gridlines and clear of both lines.
ax.annotate(f"Audio-features-and-genre model\nAUC = {auc:.3f}\n{top_hits} of its {N_PICKS} top picks {are_hits}\n"
            "(8 or more would single out the hits)", xy=(0.97, 0.22), xycoords="axes fraction",
            ha="right", va="bottom", fontsize=10, color=ms.HIGHLIGHT, fontweight="bold")

ax.set_xlim(0, 1)
ax.set_ylim(0, 1.02)
ax.xaxis.set_major_formatter(lambda v, _: f"{v:.0%}")
ax.yaxis.set_major_formatter(lambda v, _: f"{v:.0%}")
ax.set_xlabel("Share of all non-hits flagged so far (%)")
ax.set_ylabel("Share of all hits flagged so far (%)")
ax.grid(axis="y")
ax.grid(axis="x", visible=False)
ax.set_axisbelow(True)

ms.finish(ax, f"Adding genre ranks hits far better than a coin flip, AUC = {auc:.2f}",
          f"ROC curve · test tracks flagged from the highest score down · {len(test):,} tracks, "
          f"{int(test['hit'].sum()):,} of them hits",
          "Source: data/catalog.csv (test set: the 30% of tracks the model never saw in training)",
          tag=ms.tag(16, 17))
ms.save(fig, n=16, slug="audio-genre-roc")
plt.show()

# Reveal: your call on Q2·H3 beside the data.
display(HTML(ms.reveal(CALLS, "Q2·H3", holds=True, verdict="AUC 0.855, and 10 of its 15 top picks are hits; but the next section shows how the catalog's collection inflates genre.")))

# - For each genre, the share of its catalog rows that are hits (popularity at or above hit_line).
#   Call it hit_share.
# - Show the rows per genre, how many genres have no hits, and which have hits in over half their
#   rows.
rows_per_genre = catalog["genre"].value_counts()
hit_share = (catalog["popularity"] >= hit_line).groupby(catalog["genre"]).mean().sort_values()
display(ms.table(pd.Series({
    "rows per genre (fewest to most)": f"{rows_per_genre.min():,} to {rows_per_genre.max():,}",
    "genres with no hits": int((hit_share == 0).sum()),
    "genres with hits in over half their rows": ", ".join(hit_share[hit_share > 0.5].index),
    "share of hits over all tracks": f"{tracks['hit'].mean():.1%}",
}, name="value"), index_label="measure"))

# Figure: the share of each genre's rows that are hits, compared with 10%, the
#   share of hits over all tracks, from hit_share.
import sys; sys.path.append(".claude/skills/fresh-figures/scripts")
import fresh_style as ms

ms.use_theme()
# One named bar per genre: horizontal, and tall enough (0.16 inch a bar) that every name fits.
fig, ax = plt.subplots(figsize=(8, 0.16 * len(hit_share) + 1.6))

# Every genre from fewest hits (bottom) to most (top). The genres with no hits and the genres with
# hits in over half their rows (the table above) are highlighted, and each bar's share is printed
# at its tip, so the genres with no bar at all still read "0%" (shares under 1% get a decimal, so
# they don't read as 0% too).
no_hits = hit_share.index[hit_share == 0]
over_half = hit_share.index[hit_share > 0.5]
ends = set(no_hits) | set(over_half)
colors = [ms.HIGHLIGHT if genre in ends else ms.CONTEXT for genre in hit_share.index]
rows = np.arange(len(hit_share))
bars = ax.barh(rows, hit_share.values * 100, color=colors, height=0.72)
for bar, genre, share in zip(bars, hit_share.index, hit_share.values):
    label = f"{share:.1%}" if 0 < share < 0.01 else f"{share:.0%}"
    ax.annotate(label, xy=(share * 100, bar.get_y() + bar.get_height() / 2), xytext=(3, 0),
                textcoords="offset points", ha="left", va="center", fontsize=9,
                color=ms.HIGHLIGHT if genre in ends else ms.MUTED,
                fontweight="bold" if genre in ends else "normal",
                # Above the dashed line, on a background patch, so the line never runs
                # through a number (R7).
                zorder=5,
                bbox=dict(boxstyle="square,pad=0.05", facecolor=ms.BG, edgecolor="none"))
ax.set_yticks(rows, hit_share.index, fontsize=9)
ax.set_ylim(-0.8, len(hit_share) - 0.2)

# The share of hits over all tracks as a labeled reference line.
overall = tracks["hit"].mean()
ax.axvline(overall * 100, color=ms.INK, linestyle="--", linewidth=1.2, zorder=3)
ax.annotate(f"Hits over all tracks: {overall:.1%}", xy=(overall * 100, len(hit_share) * 0.27),
            xytext=(30, 0), textcoords="offset points", ha="left", va="center", fontsize=10,
            color=ms.INK, fontweight="bold",
            bbox=dict(boxstyle="square,pad=0.15", facecolor=ms.BG, edgecolor="none"),
            arrowprops=dict(arrowstyle="-", color=ms.INK, linewidth=1))

# The highlight named in words beside each end.
ax.annotate(f"{len(over_half)} genres with hits\nin over half their rows",
            xy=(hit_share.max() * 100, len(hit_share) - len(over_half) / 2 - 0.5), xytext=(40, 0),
            textcoords="offset points", ha="left", va="center", fontsize=10, color=ms.HIGHLIGHT,
            fontweight="bold")
ax.annotate(f"{len(no_hits)} genres with no hits at all", xy=(0, len(no_hits) / 2 - 0.5),
            xytext=(60, 0), textcoords="offset points", ha="left", va="center", fontsize=10,
            color=ms.HIGHLIGHT, fontweight="bold")

ax.set_xlim(0, 100)
ax.xaxis.set_major_formatter(lambda v, _: f"{v:.0f}%")
ax.set_xlabel(f"Share of the genre's rows that are hits, popularity {hit_line:g} or more (%)")
ax.grid(axis="x")
ax.grid(axis="y", visible=False)
ax.set_axisbelow(True)

ms.finish(ax,
          f"Hit shares run from {hit_share.min():.0%} to {hit_share.max():.0%} by genre: "
          f"{len(no_hits)} have none, {len(over_half)} are mostly hits",
          f"Share of hits in each of {len(hit_share)} genres · % · {len(catalog):,} catalog rows, "
          f"{rows_per_genre.min():,}–{rows_per_genre.max():,} per genre",
          "Source: data/catalog.csv (one row per track and genre search; the overall share counts "
          "each track once)", tag=ms.tag(17, 17))
ms.save(fig, n=17, slug="hit-share-by-genre")
plt.show()

danceability_pct	median streams (millions)	songs
(22.999, 57.0]	366.8	244
(57.0, 69.0]	286.4	235
(69.0, 78.0]	302.1	238
(78.0, 96.0]	253.9	230

valence_pct	median streams (millions)	songs
(3.999, 32.0]	304.1	239
(32.0, 51.0]	307.4	237
(51.0, 70.0]	290.8	239
(70.0, 97.0]	261.4	232

released	median danceability (0-100)	share of the list
before 2000	60	0.051
2000-09	62.5	0.021
2010-14	61	0.056
2015-18	59.5	0.065
2019-20	67	0.077
2021	71.5	0.124
2022	71	0.42
2023	71	0.184

danceability_pct	median streams (millions)	songs
(42.999, 62.0]	117.7	47
(62.0, 71.0]	118	46
(71.0, 79.0]	84.8	42
(79.0, 93.0]	112.7	40

danceability_pct	median popularity (0-100)	tracks
(5.129, 45.0]	28	22,401
(45.0, 57.6]	35	22,401
(57.6, 69.2]	36	22,406
(69.2, 98.5]	35	22,375

What Makes a Hit Song?¶

SIDE A · THE ANALYSIS · SHINING AND GLEANING
Shining and gleaning
The data arrives defined, cleaned and graded in the data notebook. Now we test the hunches and see what the data has to say.
¶

SIDE A · INTRO
The data
Two files: the 949 songs everyone streamed in 2023, and a catalog of nearly 90,000 tracks to hold them up against.
¶

SIDE A · TRACK 01
The observation
Everyone says the hits are danceable. Are they really more danceable than music in general?
¶

SIDE A · INTERLUDE
Hypotheses
The label's hunch, turned into six claims that could each turn out wrong.
¶

SIDE A · TRACK 02
Hypothesis 1 · Streams by danceability
If danceable songs are what people play, the most danceable songs should have the most streams. Do they?
¶

SIDE A · TRACK 03
Hypothesis 2 · Streams by valence
Danceability didn't do it. Maybe songs that sound happier get played more?
¶

SIDE A · INTERLUDE
Hypothesis 3 · Something else?
Neither sound explains streams. So what does a stream count actually measure?
¶

SIDE A · TRACK 04
Digging deeper · Streams by release year
A stream count is a running total, and older songs have had years longer to run it up.
¶

SIDE A · TRACK 05
Digging deeper · Danceability by release year
If the newest songs are also the most danceable, danceability was set up to look bad.
¶

SIDE A · TRACK 06
Digging deeper · Songs released in 2023 only
Line up songs of the same age, and danceability finally gets a fair test.
¶

SIDE A · TRACK 07
Digging deeper · What the list leaves out
A list of hits can't show what makes a hit. The catalog holds hits and misses alike.
¶

SIDE A · TRACK 08
Digging deeper · Danceability across the whole catalog
Danceability's fairest test yet: every track, hits and misses, on one popularity scale.
¶

SIDE A · TRACK 09
Digging deeper · Popularity by genre
If how a track sounds doesn't separate hits from misses, does the kind of music it is?
¶

SIDE B · TRACK 10
Digging deeper · Do tracks group by sound?
Genre is a label people give. Set it aside: does sound alone sort tracks into natural groups?
¶

SIDE B · TRACK 11
Digging deeper · The sound groups, popularity, and genre
Four groups by sound: what each sounds like, how popular each is, and whether genres live in one.
¶

SIDE B · INTERLUDE
Synthesis · What stood in for danceability
Every finding so far, side by side: what really goes with streams and success?
¶

SIDE B · INTERLUDE
Question 2 · Could we have predicted a hit before release?
The label can't wait for streams to pile up. Could anything known on release day have picked the hits?
¶

SIDE B · TRACK 12
Question 2 · The yardstick: 15 tracks picked at random
Before judging any model, measure what luck alone achieves.
¶

SIDE B · TRACK 13
Question 2 · The strongest predictor we cannot use
One column tracks streams far better than any sound. It just arrives too late.
¶

SIDE B · TRACK 14
Question 2 · Hypothesis 1 · Danceability alone
If the label's belief were right, danceability alone should pick the hits.
¶

SIDE B · TRACK 15
Question 2 · Hypothesis 2 · All the audio features
Maybe one measure isn't enough. Does everything Spotify measures about the sound do better?
¶

SIDE B · TRACK 16
Question 2 · Hypothesis 3 · The audio features and genre
Genre is known before release, and it split popularity 57 points. Does it pick the hits?
¶

SIDE B · TRACK 17
Question 2 · How the catalog was collected
Before trusting genre, check how the genres got into the data.
¶

SIDE B · OUTRO
Limits · What this analysis cannot tell us
Every finding has edges. Here they are, so nobody leans on a number harder than it can hold.
¶

SIDE B · OUTRO
The memo · To the Head of A&R
Both answers on one page: what we found, what to do, and what to get next.
¶

SIDE B · OUTRO
Recap · What each hypothesis showed
Six claims, six verdicts, each with the number that decides it.
¶

what we counted	count
songs on the 2023 list	949
catalog rows	113,999
unique tracks in the catalog	89,740
tracks whose sound was measured	89,583

released	songs	median streams (millions)
before 2000	48	663.8
2000-09	20	1,229.3
2010-14	53	1,131.1
2015-18	62	948.4
2019-20	73	807
2021	118	505.3
2022	399	222.1
2023	175	107.8

measure	value
songs on the list	949
songs with a damaged name in the file	109
songs with a damaged name after the approved repairs	40
songs found in the catalog	308
not found: same artists in the catalog, in a different order (flagged in the data notebook)	36
median popularity of the found songs	84
median popularity of every track	33
hit_line	60
share of found songs at or above hit_line	0.925

genre	mean popularity (0-100)	end
iranian	2.2	least popular
romance	3.2	least popular
latin	8.3	least popular
detroit-techno	11.2	least popular
chicago-house	12.3	least popular
grunge	49.6	most popular
sad	52.4	most popular
chill	53.7	most popular
k-pop	57	most popular
pop-film	59.3	most popular

danceability quarter	mean popularity (0-100)
(5.129, 45.0]	30.4
(45.0, 57.6]	34.4
(57.6, 69.2]	34.4
(69.2, 98.5]	33.6

k	within-cluster spread	silhouette
2	405,247	0.285
3	336,671	0.233
4	285,604	0.247
5	247,364	0.25
6	211,683	0.273
7	192,076	0.279
8	176,983	0.254

sound group	danceability_pct	energy_pct	acousticness_pct	valence_pct	speechiness_pct	instrumentalness_pct	popularity	tracks
1	63.7	66.8	38	50.7	39.9	0	26	3,884
2	68.2	72.3	17.3	70.8	5.1	0	35	36,196
3	46.5	29.4	83.2	27.4	3.8	0.1	30	21,198
4	49.1	83.9	1	29	5.4	0.3	33	28,305

genre	share in its most common sound group	where
german	0.369	least concentrated
brazil	0.393	least concentrated
anime	0.402	least concentrated
goth	0.606	median genre
black-metal	0.959	most concentrated
death-metal	0.978	most concentrated
grindcore	0.98	most concentrated

measure	value
hit_line (popularity)	60
training tracks (70%)	62,708
test tracks (30%)	26,875
share of test tracks that are hits	0.1094

measure	value
average hits among 15 random picks	1.65
chance of 5 or more	0.0177
chance of 8 or more	0.0002

measure	value
rows per genre (fewest to most)	999 to 1,000
genres with no hits	8
genres with hits in over half their rows	k-pop, pop-film, pop
share of hits over all tracks	10.9%

Question	Hypothesis	What the data showed
1	More danceable songs have more streams.	Does not hold. The least danceable quarter leads (366.8 million vs 253.9 million), r = −0.11; among same-age songs, r = 0.01.
1	More upbeat songs, those with higher valence, have more streams.	Does not hold. The happiest quarter is lowest (261.4 million), r = −0.04.
1	Something other than how a song sounds explains its streams.	Holds. Age: 2000–09 songs have 1,229 million streams, 2023 releases 108 million; the newest songs are the most danceable.
2	Danceability alone would have singled out the hits.	Does not hold. AUC 0.552; 1 of its 15 top picks is a hit (luck: 1.65; the bar: 8).
2	All the audio features together would have singled out the hits.	Does not hold. AUC 0.644; 3 of 15 top picks are hits.
2	Adding genre, which is known before release, to the audio features would have singled out the hits.	Holds on this test, not proven. AUC 0.855; 10 of 15 top picks are hits, but genres' hit shares (0% to 64%) come partly from how the catalog was collected.

What Makes a Hit Song?¶

SIDE A · THE ANALYSIS · SHINING AND GLEANINGShining and gleaningThe data arrives defined, cleaned and graded in the data notebook. Now we test the hunches and see what the data has to say.¶

SIDE A · INTROThe dataTwo files: the 949 songs everyone streamed in 2023, and a catalog of nearly 90,000 tracks to hold them up against.¶

SIDE A · TRACK 01The observationEveryone says the hits are danceable. Are they really more danceable than music in general?¶

SIDE A · INTERLUDEHypothesesThe label's hunch, turned into six claims that could each turn out wrong.¶

SIDE A · TRACK 02Hypothesis 1 · Streams by danceabilityIf danceable songs are what people play, the most danceable songs should have the most streams. Do they?¶

SIDE A · TRACK 03Hypothesis 2 · Streams by valenceDanceability didn't do it. Maybe songs that sound happier get played more?¶

SIDE A · INTERLUDEHypothesis 3 · Something else?Neither sound explains streams. So what does a stream count actually measure?¶

SIDE A · TRACK 04Digging deeper · Streams by release yearA stream count is a running total, and older songs have had years longer to run it up.¶

SIDE A · TRACK 05Digging deeper · Danceability by release yearIf the newest songs are also the most danceable, danceability was set up to look bad.¶

SIDE A · TRACK 06Digging deeper · Songs released in 2023 onlyLine up songs of the same age, and danceability finally gets a fair test.¶

SIDE A · TRACK 07Digging deeper · What the list leaves outA list of hits can't show what makes a hit. The catalog holds hits and misses alike.¶

SIDE A · TRACK 08Digging deeper · Danceability across the whole catalogDanceability's fairest test yet: every track, hits and misses, on one popularity scale.¶

SIDE A · TRACK 09Digging deeper · Popularity by genreIf how a track sounds doesn't separate hits from misses, does the kind of music it is?¶

SIDE B · TRACK 10Digging deeper · Do tracks group by sound?Genre is a label people give. Set it aside: does sound alone sort tracks into natural groups?¶

SIDE B · TRACK 11Digging deeper · The sound groups, popularity, and genreFour groups by sound: what each sounds like, how popular each is, and whether genres live in one.¶

SIDE B · INTERLUDESynthesis · What stood in for danceabilityEvery finding so far, side by side: what really goes with streams and success?¶

SIDE B · INTERLUDEQuestion 2 · Could we have predicted a hit before release?The label can't wait for streams to pile up. Could anything known on release day have picked the hits?¶

SIDE B · TRACK 12Question 2 · The yardstick: 15 tracks picked at randomBefore judging any model, measure what luck alone achieves.¶

SIDE B · TRACK 13Question 2 · The strongest predictor we cannot useOne column tracks streams far better than any sound. It just arrives too late.¶

SIDE B · TRACK 14Question 2 · Hypothesis 1 · Danceability aloneIf the label's belief were right, danceability alone should pick the hits.¶

SIDE B · TRACK 15Question 2 · Hypothesis 2 · All the audio featuresMaybe one measure isn't enough. Does everything Spotify measures about the sound do better?¶

SIDE B · TRACK 16Question 2 · Hypothesis 3 · The audio features and genreGenre is known before release, and it split popularity 57 points. Does it pick the hits?¶

SIDE B · TRACK 17Question 2 · How the catalog was collectedBefore trusting genre, check how the genres got into the data.¶

SIDE B · OUTROLimits · What this analysis cannot tell usEvery finding has edges. Here they are, so nobody leans on a number harder than it can hold.¶

SIDE B · OUTROThe memo · To the Head of A&RBoth answers on one page: what we found, what to do, and what to get next.¶

SIDE B · OUTRORecap · What each hypothesis showedSix claims, six verdicts, each with the number that decides it.¶

SIDE A · THE ANALYSIS · SHINING AND GLEANING
Shining and gleaning
The data arrives defined, cleaned and graded in the data notebook. Now we test the hunches and see what the data has to say.
¶

SIDE A · INTRO
The data
Two files: the 949 songs everyone streamed in 2023, and a catalog of nearly 90,000 tracks to hold them up against.
¶

SIDE A · TRACK 01
The observation
Everyone says the hits are danceable. Are they really more danceable than music in general?
¶

SIDE A · INTERLUDE
Hypotheses
The label's hunch, turned into six claims that could each turn out wrong.
¶

SIDE A · TRACK 02
Hypothesis 1 · Streams by danceability
If danceable songs are what people play, the most danceable songs should have the most streams. Do they?
¶

SIDE A · TRACK 03
Hypothesis 2 · Streams by valence
Danceability didn't do it. Maybe songs that sound happier get played more?
¶

SIDE A · INTERLUDE
Hypothesis 3 · Something else?
Neither sound explains streams. So what does a stream count actually measure?
¶

SIDE A · TRACK 04
Digging deeper · Streams by release year
A stream count is a running total, and older songs have had years longer to run it up.
¶

SIDE A · TRACK 05
Digging deeper · Danceability by release year
If the newest songs are also the most danceable, danceability was set up to look bad.
¶

SIDE A · TRACK 06
Digging deeper · Songs released in 2023 only
Line up songs of the same age, and danceability finally gets a fair test.
¶

SIDE A · TRACK 07
Digging deeper · What the list leaves out
A list of hits can't show what makes a hit. The catalog holds hits and misses alike.
¶

SIDE A · TRACK 08
Digging deeper · Danceability across the whole catalog
Danceability's fairest test yet: every track, hits and misses, on one popularity scale.
¶

SIDE A · TRACK 09
Digging deeper · Popularity by genre
If how a track sounds doesn't separate hits from misses, does the kind of music it is?
¶

SIDE B · TRACK 10
Digging deeper · Do tracks group by sound?
Genre is a label people give. Set it aside: does sound alone sort tracks into natural groups?
¶

SIDE B · TRACK 11
Digging deeper · The sound groups, popularity, and genre
Four groups by sound: what each sounds like, how popular each is, and whether genres live in one.
¶

SIDE B · INTERLUDE
Synthesis · What stood in for danceability
Every finding so far, side by side: what really goes with streams and success?
¶

SIDE B · INTERLUDE
Question 2 · Could we have predicted a hit before release?
The label can't wait for streams to pile up. Could anything known on release day have picked the hits?
¶

SIDE B · TRACK 12
Question 2 · The yardstick: 15 tracks picked at random
Before judging any model, measure what luck alone achieves.
¶

SIDE B · TRACK 13
Question 2 · The strongest predictor we cannot use
One column tracks streams far better than any sound. It just arrives too late.
¶

SIDE B · TRACK 14
Question 2 · Hypothesis 1 · Danceability alone
If the label's belief were right, danceability alone should pick the hits.
¶

SIDE B · TRACK 15
Question 2 · Hypothesis 2 · All the audio features
Maybe one measure isn't enough. Does everything Spotify measures about the sound do better?
¶

SIDE B · TRACK 16
Question 2 · Hypothesis 3 · The audio features and genre
Genre is known before release, and it split popularity 57 points. Does it pick the hits?
¶

SIDE B · TRACK 17
Question 2 · How the catalog was collected
Before trusting genre, check how the genres got into the data.
¶

SIDE B · OUTRO
Limits · What this analysis cannot tell us
Every finding has edges. Here they are, so nobody leans on a number harder than it can hold.
¶

SIDE B · OUTRO
The memo · To the Head of A&R
Both answers on one page: what we found, what to do, and what to get next.
¶

SIDE B · OUTRO
Recap · What each hypothesis showed
Six claims, six verdicts, each with the number that decides it.
¶