Shared Emotional Atlas — Participant × Emotion Explorer

This notebook uses the already fitted shared 3D UMAP and audited metadata to build:

1. an interactive participant × emotion explorer;
2. participant hull and centroid overlays;
3. a shared-space occupancy map showing how many participants visit each local region;
4. a local participant-entropy map showing where the shared embedding is highly individual versus broadly shared;
5. summary tables for participant spread and peripheral occupation.

The notebook does not refit UMAP.

from pathlib import Path
import json
import numpy as np
import pandas as pd
import plotly.graph_objects as go

try:
    from scipy.spatial import ConvexHull
    SCIPY_AVAILABLE = True
except Exception:
    SCIPY_AVAILABLE = False

print("SciPy available:", SCIPY_AVAILABLE)

SciPy available: True

1. Configuration

BASE_DIR = Path("shared_umap_76D_SOT_exports")

Z_PATH = BASE_DIR / "Zumap3_76D_SOT_AUDITED.npy"
SUBJECT_PATH = BASE_DIR / "frame_subject_76D_SOT_AUDITED.npy"
EMOTION_PATH = BASE_DIR / "frame_emotion_76D_SOT_AUDITED.npy"
SEGMENTS_PATH = BASE_DIR / "segments_76D_SOT_AUDITED.csv"

OUT_DIR = BASE_DIR / "shared_atlas_outputs"
OUT_DIR.mkdir(parents=True, exist_ok=True)

RANDOM_SEED = 42
MAX_POINTS_PER_GROUP = 5000
POINT_SIZE = 2.0
POINT_OPACITY = 0.45
CONTEXT_OPACITY = 0.025

VOXEL_BINS = 22
VOXEL_MIN_COUNT = 20

PARTICIPANT_COLOURS = {
    "P1": "#1f77b4",
    "P2": "#ff7f0e",
    "P3": "#2ca02c",
    "P4": "#d62728",
    "P5": "#9467bd",
    "P6": "#17becf",
    "P7": "#bcbd22",
}

EMOTION_COLOURS = {
    "anger": "#1f77b4",
    "contempt": "#d62728",
    "disgust": "#2ca02c",
    "fear": "#9467bd",
    "flow": "#ff7f0e",
    "happiness": "#17becf",
    "neutral": "#8c564b",
    "sadness": "#e377c2",
    "surprise": "#7f7f7f",
}

2. Load and validate

Z = np.load(Z_PATH)
subjects_raw = np.load(SUBJECT_PATH, allow_pickle=True)
emotions_raw = np.load(EMOTION_PATH, allow_pickle=True)
segments = pd.read_csv(SEGMENTS_PATH)

def normalise_subject(value):
    text = str(value).strip().upper()
    if text.startswith("P"):
        return text
    try:
        return f"P{int(float(text))}"
    except Exception:
        return text

def normalise_emotion(value):
    return str(value).strip().lower()

subjects = np.array([normalise_subject(v) for v in subjects_raw], dtype=object)
emotions = np.array([normalise_emotion(v) for v in emotions_raw], dtype=object)

assert Z.ndim == 2 and Z.shape[1] == 3, Z.shape
assert len(subjects) == len(Z)
assert len(emotions) == len(Z)

participant_order = [p for p in PARTICIPANT_COLOURS if p in set(subjects)]
emotion_order = [e for e in EMOTION_COLOURS if e in set(emotions)]

print("Shared UMAP:", Z.shape)
print("Participants:", participant_order)
print("Emotions:", emotion_order)
print("Segments:", segments.shape)

Shared UMAP: (228574, 3)
Participants: ['P1', 'P2', 'P3', 'P4', 'P5', 'P6', 'P7']
Emotions: ['anger', 'contempt', 'disgust', 'fear', 'flow', 'happiness', 'neutral', 'sadness', 'surprise']
Segments: (630, 6)

def clean_scene(background="black"):
    grid = "#222222" if background == "black" else "#e8e8e8"
    return dict(
        xaxis=dict(title="", showticklabels=False, ticks="", showgrid=True, gridcolor=grid,
                   zeroline=False, showbackground=False),
        yaxis=dict(title="", showticklabels=False, ticks="", showgrid=True, gridcolor=grid,
                   zeroline=False, showbackground=False),
        zaxis=dict(title="", showticklabels=False, ticks="", showgrid=True, gridcolor=grid,
                   zeroline=False, showbackground=False),
        bgcolor=background,
        aspectmode="data",
        camera=dict(eye=dict(x=1.45, y=1.45, z=0.9)),
    )

rng = np.random.default_rng(RANDOM_SEED)

3. Participant × emotion explorer

This creates one interactive HTML file with dropdowns for:

• participant;
• emotion;
• colouring mode;
• context visibility.

The selected subset remains in colour while all other frames can be faded into grey context.

# Stratified sample by participant × emotion
sample_idx = []
for participant in participant_order:
    for emotion in emotion_order:
        idx = np.flatnonzero((subjects == participant) & (emotions == emotion))
        if len(idx) == 0:
            continue
        n = min(MAX_POINTS_PER_GROUP, len(idx))
        chosen = rng.choice(idx, size=n, replace=False)
        sample_idx.append(np.sort(chosen))

sample_idx = np.concatenate(sample_idx)
sample_idx.sort()

Zs = Z[sample_idx]
Ss = subjects[sample_idx]
Es = emotions[sample_idx]

print("Displayed frames:", len(sample_idx), "of", len(Z))

Displayed frames: 223142 of 228574

# One trace per participant × emotion group
fig_explorer = go.Figure()
trace_meta = []

for participant in participant_order:
    for emotion in emotion_order:
        mask = (Ss == participant) & (Es == emotion)
        if not np.any(mask):
            continue

        trace_meta.append((participant, emotion))
        fig_explorer.add_trace(go.Scatter3d(
            x=Zs[mask, 0],
            y=Zs[mask, 1],
            z=Zs[mask, 2],
            mode="markers",
            name=f"{participant} · {emotion.title()}",
            legendgroup=participant,
            showlegend=False,
            marker=dict(
                size=POINT_SIZE,
                color=PARTICIPANT_COLOURS[participant],
                opacity=POINT_OPACITY,
            ),
            customdata=np.column_stack([
                np.repeat(participant, mask.sum()),
                np.repeat(emotion, mask.sum()),
            ]),
            hovertemplate=(
                "<b>%{customdata[0]}</b><br>"
                "Emotion: %{customdata[1]}<br>"
                "UMAP: (%{x:.3f}, %{y:.3f}, %{z:.3f})"
                "<extra></extra>"
            ),
        ))

def explorer_state(selected_participant="All", selected_emotion="All",
                   colour_mode="Participant", show_context=True):
    opacities, colours, sizes = [], [], []

    for participant, emotion in trace_meta:
        selected = (
            (selected_participant == "All" or participant == selected_participant)
            and
            (selected_emotion == "All" or emotion == selected_emotion)
        )

        if selected:
            opacity = POINT_OPACITY
            size = POINT_SIZE + 0.5
            colour = (
                PARTICIPANT_COLOURS[participant]
                if colour_mode == "Participant"
                else EMOTION_COLOURS[emotion]
            )
        else:
            opacity = CONTEXT_OPACITY if show_context else 0.0
            size = POINT_SIZE
            colour = "#9a9a9a"

        opacities.append(opacity)
        colours.append(colour)
        sizes.append(size)

    return {
        "marker.opacity": opacities,
        "marker.color": colours,
        "marker.size": sizes,
    }

# Independent dropdowns are generated as preset combinations.
# This keeps the exported HTML self-contained and avoids custom JavaScript.
buttons = []

buttons.append(dict(
    label="All participants · all emotions",
    method="restyle",
    args=[explorer_state()]
))

for participant in participant_order:
    buttons.append(dict(
        label=f"{participant} · all emotions",
        method="restyle",
        args=[explorer_state(selected_participant=participant)]
    ))

for emotion in emotion_order:
    buttons.append(dict(
        label=f"All participants · {emotion.title()}",
        method="restyle",
        args=[explorer_state(selected_emotion=emotion, colour_mode="Emotion")]
    ))

for participant in participant_order:
    for emotion in ["flow", "neutral", "fear"]:
        if (participant, emotion) not in trace_meta:
            continue
        buttons.append(dict(
            label=f"{participant} · {emotion.title()}",
            method="restyle",
            args=[explorer_state(
                selected_participant=participant,
                selected_emotion=emotion,
                colour_mode="Emotion",
            )]
        ))

fig_explorer.update_layout(
    title=dict(
        text="Shared emotional atlas — participant × emotion explorer",
        x=0.5,
        xanchor="center",
    ),
    scene=clean_scene("black"),
    paper_bgcolor="black",
    font=dict(color="white", family="Arial, sans-serif"),
    showlegend=False,
    updatemenus=[dict(
        type="dropdown",
        x=0.01,
        y=0.99,
        xanchor="left",
        yanchor="top",
        buttons=buttons,
        bgcolor="#111111",
        font=dict(color="white"),
    )],
    annotations=[dict(
        text=f"{len(Z):,} total frames · {len(Zs):,} displayed · one shared 76D → 3D UMAP",
        x=0.5, y=0.01, xref="paper", yref="paper",
        showarrow=False, font=dict(size=11, color="#bbbbbb"),
    )],
    width=1150,
    height=850,
    margin=dict(l=0, r=0, t=70, b=30),
)

explorer_out = OUT_DIR / "shared_atlas_participant_emotion_explorer.html"
fig_explorer.write_html(explorer_out, include_plotlyjs="cdn")
print("Wrote:", explorer_out.resolve())

fig_explorer.show()

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/shared_atlas_participant_emotion_explorer.html

4. Participant centroids, spread, and peripheral occupation

global_centroid = Z.mean(axis=0)
global_radius = np.linalg.norm(Z - global_centroid, axis=1)
global_r95 = np.quantile(global_radius, 0.95)

summary_rows = []

for participant in participant_order:
    pts = Z[subjects == participant]
    centroid = pts.mean(axis=0)
    radial = np.linalg.norm(pts - centroid, axis=1)
    global_radial = np.linalg.norm(pts - global_centroid, axis=1)

    row = {
        "participant": participant,
        "frames": len(pts),
        "mean_distance_to_own_centroid": radial.mean(),
        "median_distance_to_own_centroid": np.median(radial),
        "r95_own_centroid": np.quantile(radial, 0.95),
        "max_distance_to_own_centroid": radial.max(),
        "mean_distance_to_global_centroid": global_radial.mean(),
        "share_beyond_global_95pct_radius": np.mean(global_radial > global_r95),
    }

    if SCIPY_AVAILABLE:
        hull_sample_n = min(12000, len(pts))
        hull_idx = rng.choice(len(pts), size=hull_sample_n, replace=False)
        hull_pts = pts[hull_idx]
        try:
            row["sampled_convex_hull_volume"] = ConvexHull(hull_pts).volume
        except Exception:
            row["sampled_convex_hull_volume"] = np.nan
    else:
        row["sampled_convex_hull_volume"] = np.nan

    summary_rows.append(row)

participant_summary = pd.DataFrame(summary_rows).sort_values(
    "mean_distance_to_global_centroid", ascending=False
)

participant_summary

	participant	frames	mean_distance_to_own_centroid	median_distance_to_own_centroid	r95_own_centroid	max_distance_to_own_centroid	mean_distance_to_global_centroid	share_beyond_global_95pct_radius	sampled_convex_hull_volume
1	P2	27505	10.050427	9.724325	14.306516	16.264431	10.807572	0.205344	7352.704609
5	P6	29910	7.945413	8.114986	11.820291	13.902251	7.962641	0.071983	5477.428029
2	P3	23979	7.705568	8.368054	12.207925	13.629522	7.823942	0.036032	5736.489296
6	P7	41688	7.346733	7.439608	11.506913	14.477082	7.588328	0.035886	5852.489985
0	P1	37682	7.003020	7.088642	11.680943	14.167866	7.474209	0.004538	5763.610374
3	P4	26262	6.520965	6.031725	11.665685	15.162266	7.391629	0.026883	5484.790872
4	P5	41548	6.898783	7.196618	11.196620	13.459812	7.042881	0.009411	5935.642833

summary_csv = OUT_DIR / "participant_shared_umap_spread_summary.csv"
participant_summary.to_csv(summary_csv, index=False)
print("Wrote:", summary_csv.resolve())

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/participant_shared_umap_spread_summary.csv

5. Participant hull and centroid overlays

MAX_HULL_POINTS = 9000

fig_hulls = go.Figure()

# Grey context
context_n = min(50000, len(Z))
context_idx = rng.choice(len(Z), size=context_n, replace=False)

fig_hulls.add_trace(go.Scatter3d(
    x=Z[context_idx, 0],
    y=Z[context_idx, 1],
    z=Z[context_idx, 2],
    mode="markers",
    name="Shared context",
    marker=dict(size=1.5, color="#8a8a8a", opacity=0.035),
    hoverinfo="skip",
    showlegend=False,
))

hull_trace_indices = {}

for participant in participant_order:
    pts = Z[subjects == participant]
    centroid = pts.mean(axis=0)

    # Centroid
    fig_hulls.add_trace(go.Scatter3d(
        x=[centroid[0]], y=[centroid[1]], z=[centroid[2]],
        mode="markers+text",
        text=[participant],
        textposition="top center",
        marker=dict(
            size=8,
            color=PARTICIPANT_COLOURS[participant],
            line=dict(color="white", width=1),
        ),
        name=f"{participant} centroid",
        visible=False,
        showlegend=False,
    ))
    centroid_trace_idx = len(fig_hulls.data) - 1

    mesh_trace_idx = None
    if SCIPY_AVAILABLE:
        n = min(MAX_HULL_POINTS, len(pts))
        idx = rng.choice(len(pts), size=n, replace=False)
        hp = pts[idx]
        try:
            hull = ConvexHull(hp)
            fig_hulls.add_trace(go.Mesh3d(
                x=hp[:, 0], y=hp[:, 1], z=hp[:, 2],
                i=hull.simplices[:, 0],
                j=hull.simplices[:, 1],
                k=hull.simplices[:, 2],
                color=PARTICIPANT_COLOURS[participant],
                opacity=0.10,
                name=f"{participant} sampled hull",
                visible=False,
                showlegend=False,
                hoverinfo="skip",
            ))
            mesh_trace_idx = len(fig_hulls.data) - 1
        except Exception as exc:
            print(f"Hull failed for {participant}:", exc)

    hull_trace_indices[participant] = (centroid_trace_idx, mesh_trace_idx)

buttons = []

all_centroids_visible = [True] + [False] * (len(fig_hulls.data) - 1)
for participant, (centroid_idx, mesh_idx) in hull_trace_indices.items():
    all_centroids_visible[centroid_idx] = True

buttons.append(dict(
    label="All centroids",
    method="update",
    args=[{"visible": all_centroids_visible}, {"title.text": "Participant centroids in the shared UMAP"}]
))

for participant, (centroid_idx, mesh_idx) in hull_trace_indices.items():
    visible = [False] * len(fig_hulls.data)
    visible[0] = True
    visible[centroid_idx] = True
    if mesh_idx is not None:
        visible[mesh_idx] = True

    buttons.append(dict(
        label=f"{participant} hull",
        method="update",
        args=[
            {"visible": visible},
            {"title.text": f"{participant} occupied volume in the shared UMAP"}
        ]
    ))

fig_hulls.update_layout(
    title=dict(text="Participant hulls in the shared UMAP", x=0.5, xanchor="center"),
    scene=clean_scene("black"),
    paper_bgcolor="black",
    font=dict(color="white", family="Arial, sans-serif"),
    showlegend=False,
    updatemenus=[dict(
        type="dropdown",
        x=0.01, y=0.99,
        xanchor="left", yanchor="top",
        buttons=buttons,
        bgcolor="#111111",
        font=dict(color="white"),
    )],
    width=1150,
    height=850,
    margin=dict(l=0, r=0, t=70, b=20),
)

hulls_out = OUT_DIR / "shared_atlas_participant_hulls.html"
fig_hulls.write_html(hulls_out, include_plotlyjs="cdn")
print("Wrote:", hulls_out.resolve())

fig_hulls.show()

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/shared_atlas_participant_hulls.html

6. Shared-space occupancy map

The shared UMAP is partitioned into a 3D voxel grid.

For every occupied voxel, the notebook counts:

• total frames;
• number of distinct participants;
• dominant participant;
• participant entropy.

Interpretation:

• 1 participant = highly individual territory;
• 7 participants = broadly shared territory;
• low entropy = dominated by one participant;
• high entropy = participants contribute more evenly.

def build_voxel_table(Z, subjects, bins=22, min_count=20):
    mins = Z.min(axis=0)
    maxs = Z.max(axis=0)
    edges = [np.linspace(mins[d], maxs[d], bins + 1) for d in range(3)]

    ijk = np.column_stack([
        np.clip(np.digitize(Z[:, d], edges[d]) - 1, 0, bins - 1)
        for d in range(3)
    ])

    records = {}
    for idx, key_arr in enumerate(ijk):
        key = tuple(int(v) for v in key_arr)
        if key not in records:
            records[key] = {"count": 0, "participant_counts": {}}
        records[key]["count"] += 1
        participant = subjects[idx]
        pc = records[key]["participant_counts"]
        pc[participant] = pc.get(participant, 0) + 1

    rows = []
    max_entropy = np.log2(len(participant_order))

    for (i, j, k), rec in records.items():
        if rec["count"] < min_count:
            continue

        center = np.array([
            0.5 * (edges[0][i] + edges[0][i + 1]),
            0.5 * (edges[1][j] + edges[1][j + 1]),
            0.5 * (edges[2][k] + edges[2][k + 1]),
        ])

        counts = rec["participant_counts"]
        values = np.array(list(counts.values()), dtype=float)
        probs = values / values.sum()
        entropy = -(probs * np.log2(probs)).sum()
        norm_entropy = entropy / max_entropy if max_entropy > 0 else 0.0

        dominant = max(counts, key=counts.get)
        dominant_share = counts[dominant] / rec["count"]

        rows.append({
            "x": center[0],
            "y": center[1],
            "z": center[2],
            "count": rec["count"],
            "participant_count": len(counts),
            "dominant_participant": dominant,
            "dominant_share": dominant_share,
            "participant_entropy": entropy,
            "participant_entropy_normalised": norm_entropy,
        })

    return pd.DataFrame(rows)

voxels = build_voxel_table(
    Z, subjects,
    bins=VOXEL_BINS,
    min_count=VOXEL_MIN_COUNT,
)

print("Occupied voxels retained:", len(voxels))
voxels.head()

Occupied voxels retained: 2936

	x	y	z	count	participant_count	dominant_participant	dominant_share	participant_entropy	participant_entropy_normalised
0	1.270103	6.826300	-6.807026	45	1	P1	1.000000	-0.00000	-0.000000
1	7.473257	0.234696	11.694004	56	1	P1	1.000000	-0.00000	-0.000000
2	6.232626	0.234696	10.537690	68	3	P2	0.573529	1.07663	0.383503
3	3.751365	1.333297	10.537690	50	1	P1	1.000000	-0.00000	-0.000000
4	2.510734	1.333297	10.537690	120	1	P1	1.000000	-0.00000	-0.000000

voxel_csv = OUT_DIR / "shared_umap_participant_voxels.csv"
voxels.to_csv(voxel_csv, index=False)
print("Wrote:", voxel_csv.resolve())

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/shared_umap_participant_voxels.csv

6a. Number of participants per occupied voxel

fig_occupancy = go.Figure(go.Scatter3d(
    x=voxels["x"],
    y=voxels["y"],
    z=voxels["z"],
    mode="markers",
    marker=dict(
        size=np.clip(np.sqrt(voxels["count"]) / 2.5, 2.5, 9),
        color=voxels["participant_count"],
        colorscale="Viridis",
        cmin=1,
        cmax=len(participant_order),
        opacity=0.72,
        colorbar=dict(
            title="Participants<br>per voxel",
            tickmode="array",
            tickvals=list(range(1, len(participant_order) + 1)),
        ),
    ),
    customdata=np.column_stack([
        voxels["count"],
        voxels["participant_count"],
        voxels["dominant_participant"],
        voxels["dominant_share"],
    ]),
    hovertemplate=(
        "Frames: %{customdata[0]}<br>"
        "Participants: %{customdata[1]}<br>"
        "Dominant: %{customdata[2]}<br>"
        "Dominant share: %{customdata[3]:.1%}"
        "<extra></extra>"
    )
))

fig_occupancy.update_layout(
    title=dict(
        text="How many participants occupy each local region?",
        x=0.5, xanchor="center",
    ),
    scene=clean_scene("black"),
    paper_bgcolor="black",
    font=dict(color="white", family="Arial, sans-serif"),
    width=1150,
    height=850,
    margin=dict(l=0, r=0, t=70, b=20),
)

occupancy_out = OUT_DIR / "shared_atlas_participant_occupancy.html"
fig_occupancy.write_html(occupancy_out, include_plotlyjs="cdn")
print("Wrote:", occupancy_out.resolve())

fig_occupancy.show()

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/shared_atlas_participant_occupancy.html

6b. Local participant entropy

fig_entropy = go.Figure(go.Scatter3d(
    x=voxels["x"],
    y=voxels["y"],
    z=voxels["z"],
    mode="markers",
    marker=dict(
        size=np.clip(np.sqrt(voxels["count"]) / 2.5, 2.5, 9),
        color=voxels["participant_entropy_normalised"],
        colorscale="Plasma",
        cmin=0,
        cmax=1,
        opacity=0.72,
        colorbar=dict(title="Normalised<br>participant entropy"),
    ),
    customdata=np.column_stack([
        voxels["count"],
        voxels["participant_count"],
        voxels["dominant_participant"],
        voxels["dominant_share"],
        voxels["participant_entropy_normalised"],
    ]),
    hovertemplate=(
        "Frames: %{customdata[0]}<br>"
        "Participants: %{customdata[1]}<br>"
        "Dominant: %{customdata[2]}<br>"
        "Dominant share: %{customdata[3]:.1%}<br>"
        "Normalised entropy: %{customdata[4]:.3f}"
        "<extra></extra>"
    )
))

fig_entropy.update_layout(
    title=dict(
        text="Local participant entropy — individual versus shared territory",
        x=0.5, xanchor="center",
    ),
    scene=clean_scene("black"),
    paper_bgcolor="black",
    font=dict(color="white", family="Arial, sans-serif"),
    width=1150,
    height=850,
    margin=dict(l=0, r=0, t=70, b=20),
)

entropy_out = OUT_DIR / "shared_atlas_participant_entropy.html"
fig_entropy.write_html(entropy_out, include_plotlyjs="cdn")
print("Wrote:", entropy_out.resolve())

fig_entropy.show()

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/shared_atlas_participant_entropy.html

Interpretation guardrails

These plots describe the geometry of the shared UMAP embedding, not exact distances in the original 76-dimensional feature space.

Recommended language:

• “participant-specific density and peripheral occupation”;
• “local regions dominated by one participant”;
• “shared versus individually occupied regions of the common embedding”;
• “descriptive evidence motivating quantitative follow-up.”

Avoid claiming that:

• outer position has inherent psychological meaning;
• convex hulls are exact high-dimensional volumes;
• visible overlap proves equivalence;
• visible separation alone proves categorical difference.

8. Quantitative atlas refinements

The exploratory atlas is now supplemented with quantitative summaries designed to reduce over-interpretation of the 3D display.

This section adds:

1. voxel-occupancy percentages for regions visited by 1–7 participants;
2. frame-weighted occupancy percentages;
3. entropy summary statistics;
4. exclusive-volume share by participant;
5. dominance-strength distributions;
6. peripheral ownership;
7. sensitivity checks across voxel resolutions and minimum-count thresholds;
8. emotion-specific occupancy summaries, especially Flow and Neutral.

These metrics describe the shared UMAP embedding, not exact geometry in the original 76-dimensional feature space.

8.1 Occupancy percentages

# Unweighted: each retained occupied voxel contributes equally
occupancy_summary = (
    voxels["participant_count"]
    .value_counts()
    .reindex(range(1, len(participant_order) + 1), fill_value=0)
    .rename_axis("participants_in_voxel")
    .reset_index(name="voxel_count")
)

occupancy_summary["voxel_percent"] = (
    100.0
    * occupancy_summary["voxel_count"]
    / occupancy_summary["voxel_count"].sum()
)

# Frame-weighted: dense voxels contribute in proportion to their retained frame count
frame_weighted_occupancy = (
    voxels.groupby("participant_count", as_index=False)["count"]
    .sum()
    .rename(columns={
        "participant_count": "participants_in_voxel",
        "count": "frames_in_voxels",
    })
    .set_index("participants_in_voxel")
    .reindex(range(1, len(participant_order) + 1), fill_value=0)
    .reset_index()
)

frame_weighted_occupancy["frame_percent"] = (
    100.0
    * frame_weighted_occupancy["frames_in_voxels"]
    / frame_weighted_occupancy["frames_in_voxels"].sum()
)

occupancy_combined = occupancy_summary.merge(
    frame_weighted_occupancy,
    on="participants_in_voxel",
    how="left",
)

occupancy_combined

	participants_in_voxel	voxel_count	voxel_percent	frames_in_voxels	frame_percent
0	1	1492	50.817439	85042	38.544377
1	2	943	32.118529	76538	34.690030
2	3	383	13.044959	40819	18.500775
3	4	96	3.269755	14087	6.384782
4	5	21	0.715259	4034	1.828367
5	6	1	0.034060	114	0.051669
6	7	0	0.000000	0	0.000000

occupancy_csv = OUT_DIR / "voxel_occupancy_percentages.csv"
occupancy_combined.to_csv(occupancy_csv, index=False)
print("Wrote:", occupancy_csv.resolve())

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/voxel_occupancy_percentages.csv

fig_occupancy_percent = go.Figure()

fig_occupancy_percent.add_trace(go.Bar(
    x=occupancy_combined["participants_in_voxel"],
    y=occupancy_combined["voxel_percent"],
    name="Occupied voxels",
))

fig_occupancy_percent.add_trace(go.Bar(
    x=occupancy_combined["participants_in_voxel"],
    y=occupancy_combined["frame_percent"],
    name="Frames within retained voxels",
))

fig_occupancy_percent.update_layout(
    title="Shared-space occupancy by number of participants",
    xaxis_title="Participants occupying voxel",
    yaxis_title="Percentage",
    barmode="group",
    template="plotly_white",
    width=900,
    height=520,
)

occupancy_percent_out = OUT_DIR / "voxel_occupancy_percentages.html"
fig_occupancy_percent.write_html(occupancy_percent_out, include_plotlyjs="cdn")
print("Wrote:", occupancy_percent_out.resolve())

fig_occupancy_percent.show()

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/voxel_occupancy_percentages.html

Interpret both summaries together:

• voxel percentage describes the proportion of retained spatial regions;
• frame-weighted percentage describes where the retained observations actually occur.

A large singleton-voxel percentage accompanied by a smaller singleton-frame percentage would indicate many sparse participant-specific peripheral regions surrounding a more densely shared core.

8.2 Entropy and dominance summaries

entropy_values = voxels["participant_entropy_normalised"].to_numpy()
dominance_values = voxels["dominant_share"].to_numpy()

entropy_summary = pd.DataFrame([{
    "retained_voxels": len(voxels),
    "mean_normalised_entropy": float(np.mean(entropy_values)),
    "median_normalised_entropy": float(np.median(entropy_values)),
    "entropy_p05": float(np.quantile(entropy_values, 0.05)),
    "entropy_p25": float(np.quantile(entropy_values, 0.25)),
    "entropy_p75": float(np.quantile(entropy_values, 0.75)),
    "entropy_p95": float(np.quantile(entropy_values, 0.95)),
    "voxel_percent_entropy_lt_0_2": float(100 * np.mean(entropy_values < 0.2)),
    "voxel_percent_entropy_gt_0_8": float(100 * np.mean(entropy_values > 0.8)),
    "mean_dominant_share": float(np.mean(dominance_values)),
    "median_dominant_share": float(np.median(dominance_values)),
    "voxel_percent_dominant_share_ge_0_75": float(
        100 * np.mean(dominance_values >= 0.75)
    ),
    "voxel_percent_dominant_share_ge_0_90": float(
        100 * np.mean(dominance_values >= 0.90)
    ),
}])

entropy_summary

	retained_voxels	mean_normalised_entropy	median_normalised_entropy	entropy_p05	entropy_p25	entropy_p75	entropy_p95	voxel_percent_entropy_lt_0_2	voxel_percent_entropy_gt_0_8	mean_dominant_share	median_dominant_share	voxel_percent_dominant_share_ge_0_75	voxel_percent_dominant_share_ge_0_90
0	2936	0.160336	0.0	-0.0	-0.0	0.331744	0.519296	59.264305	0.0	0.84796	1.0	68.69891	57.561308

entropy_csv = OUT_DIR / "participant_entropy_and_dominance_summary.csv"
entropy_summary.to_csv(entropy_csv, index=False)
print("Wrote:", entropy_csv.resolve())

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/participant_entropy_and_dominance_summary.csv

fig_dominance = go.Figure()

fig_dominance.add_trace(go.Histogram(
    x=voxels["dominant_share"],
    nbinsx=25,
    name="Dominant participant share",
))

fig_dominance.update_layout(
    title="Strength of local participant dominance",
    xaxis_title="Dominant participant share",
    yaxis_title="Retained occupied voxels",
    template="plotly_white",
    width=900,
    height=500,
)

dominance_out = OUT_DIR / "local_dominance_distribution.html"
fig_dominance.write_html(dominance_out, include_plotlyjs="cdn")
print("Wrote:", dominance_out.resolve())

fig_dominance.show()

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/local_dominance_distribution.html

8.3 Exclusive local territory by participant

exclusive_voxels = voxels[voxels["participant_count"] == 1].copy()

exclusive_counts = (
    exclusive_voxels["dominant_participant"]
    .value_counts()
    .reindex(participant_order, fill_value=0)
    .rename_axis("participant")
    .reset_index(name="exclusive_voxel_count")
)

exclusive_counts["share_of_all_exclusive_voxels_percent"] = (
    100.0
    * exclusive_counts["exclusive_voxel_count"]
    / max(1, exclusive_counts["exclusive_voxel_count"].sum())
)

# Share of each participant's dominant voxels that are exclusive
dominant_counts = (
    voxels["dominant_participant"]
    .value_counts()
    .reindex(participant_order, fill_value=0)
    .rename_axis("participant")
    .reset_index(name="all_dominant_voxel_count")
)

exclusive_by_participant = exclusive_counts.merge(
    dominant_counts,
    on="participant",
    how="left",
)

exclusive_by_participant["exclusive_share_of_own_dominant_voxels_percent"] = (
    100.0
    * exclusive_by_participant["exclusive_voxel_count"]
    / exclusive_by_participant["all_dominant_voxel_count"].replace(0, np.nan)
)

exclusive_by_participant

	participant	exclusive_voxel_count	share_of_all_exclusive_voxels_percent	all_dominant_voxel_count	exclusive_share_of_own_dominant_voxels_percent
0	P1	180	12.064343	448	40.178571
1	P2	387	25.938338	484	79.958678
2	P3	163	10.924933	333	48.948949
3	P4	141	9.450402	303	46.534653
4	P5	221	14.812332	505	43.762376
5	P6	159	10.656836	331	48.036254
6	P7	241	16.152815	532	45.300752

exclusive_csv = OUT_DIR / "exclusive_voxel_share_by_participant.csv"
exclusive_by_participant.to_csv(exclusive_csv, index=False)
print("Wrote:", exclusive_csv.resolve())

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/exclusive_voxel_share_by_participant.csv

8.4 Peripheral ownership

# Peripheral regions are defined here using voxel-centre distance from the
# global shared-UMAP centroid. The outer 10% can be adjusted below.
OUTER_QUANTILE = 0.90

voxel_xyz = voxels[["x", "y", "z"]].to_numpy()
voxel_radius = np.linalg.norm(voxel_xyz - global_centroid, axis=1)
outer_threshold = np.quantile(voxel_radius, OUTER_QUANTILE)

peripheral_voxels = voxels.loc[voxel_radius >= outer_threshold].copy()

peripheral_ownership = (
    peripheral_voxels["dominant_participant"]
    .value_counts()
    .reindex(participant_order, fill_value=0)
    .rename_axis("participant")
    .reset_index(name="outer_voxel_count")
)

peripheral_ownership["outer_voxel_percent"] = (
    100.0
    * peripheral_ownership["outer_voxel_count"]
    / max(1, peripheral_ownership["outer_voxel_count"].sum())
)

peripheral_ownership["outer_quantile"] = OUTER_QUANTILE
peripheral_ownership["outer_radius_threshold"] = outer_threshold

peripheral_ownership

	participant	outer_voxel_count	outer_voxel_percent	outer_quantile	outer_radius_threshold
0	P1	10	3.401361	0.9	12.06182
1	P2	160	54.421769	0.9	12.06182
2	P3	27	9.183673	0.9	12.06182
3	P4	17	5.782313	0.9	12.06182
4	P5	14	4.761905	0.9	12.06182
5	P6	34	11.564626	0.9	12.06182
6	P7	32	10.884354	0.9	12.06182

peripheral_csv = OUT_DIR / "peripheral_ownership_by_participant.csv"
peripheral_ownership.to_csv(peripheral_csv, index=False)
print("Wrote:", peripheral_csv.resolve())

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/peripheral_ownership_by_participant.csv

fig_periphery = go.Figure(go.Bar(
    x=peripheral_ownership["participant"],
    y=peripheral_ownership["outer_voxel_percent"],
))

fig_periphery.update_layout(
    title=f"Dominant ownership of the outer {int((1-OUTER_QUANTILE)*100)}% of retained voxels",
    xaxis_title="Participant",
    yaxis_title="Percentage of peripheral voxels",
    template="plotly_white",
    width=850,
    height=500,
)

periphery_out = OUT_DIR / "peripheral_ownership.html"
fig_periphery.write_html(periphery_out, include_plotlyjs="cdn")
print("Wrote:", periphery_out.resolve())

fig_periphery.show()

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/peripheral_ownership.html

8.5 Sensitivity to voxel resolution and minimum-count threshold

SENSITIVITY_BINS = [16, 20, 22, 26, 30]
SENSITIVITY_MIN_COUNTS = [10, 20, 40]

sensitivity_rows = []

for bins in SENSITIVITY_BINS:
    for min_count in SENSITIVITY_MIN_COUNTS:
        v = build_voxel_table(Z, subjects, bins=bins, min_count=min_count)
        if len(v) == 0:
            continue

        sensitivity_rows.append({
            "bins": bins,
            "min_count": min_count,
            "retained_voxels": len(v),
            "singleton_voxel_percent": float(
                100 * np.mean(v["participant_count"] == 1)
            ),
            "mean_participants_per_voxel": float(
                v["participant_count"].mean()
            ),
            "mean_normalised_entropy": float(
                v["participant_entropy_normalised"].mean()
            ),
            "median_normalised_entropy": float(
                v["participant_entropy_normalised"].median()
            ),
            "mean_dominant_share": float(
                v["dominant_share"].mean()
            ),
            "dominant_share_ge_0_75_percent": float(
                100 * np.mean(v["dominant_share"] >= 0.75)
            ),
        })

sensitivity = pd.DataFrame(sensitivity_rows)
sensitivity

	bins	min_count	retained_voxels	singleton_voxel_percent	mean_participants_per_voxel	mean_normalised_entropy	median_normalised_entropy	mean_dominant_share	dominant_share_ge_0_75_percent
0	16	10	1909	35.882661	2.209534	0.258468	0.293749	0.762089	51.440545
1	16	20	1748	30.320366	2.317506	0.281340	0.315573	0.741219	47.196796
2	16	40	1491	23.004695	2.494299	0.315547	0.338900	0.711781	41.582830
3	20	10	2944	52.513587	1.747283	0.164930	0.000000	0.844615	67.730978
4	20	20	2569	46.321526	1.848190	0.187159	0.104934	0.823588	63.254185
5	20	40	1944	36.574074	2.038580	0.225130	0.246947	0.790034	56.790123
6	22	10	3404	56.345476	1.626322	0.141280	0.000000	0.866389	72.679201
7	22	20	2936	50.817439	1.710490	0.160336	0.000000	0.847960	68.698910
8	22	40	2076	39.643545	1.899807	0.200329	0.206904	0.812049	61.608863
9	26	10	4357	67.339913	1.428965	0.095643	0.000000	0.909964	82.097774
10	26	20	3520	62.073864	1.504261	0.112124	0.000000	0.894157	78.750000
11	26	40	2142	52.054155	1.663399	0.144467	0.000000	0.864811	72.735761
12	30	10	5035	73.406157	1.327110	0.072792	0.000000	0.932470	86.454816
13	30	20	3834	68.596766	1.391497	0.086501	0.000000	0.919516	83.672405
14	30	40	2119	60.122699	1.518169	0.109822	0.000000	0.898921	79.424257

sensitivity_csv = OUT_DIR / "voxel_parameter_sensitivity.csv"
sensitivity.to_csv(sensitivity_csv, index=False)
print("Wrote:", sensitivity_csv.resolve())

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/voxel_parameter_sensitivity.csv

A pattern is more trustworthy when its direction remains stable across reasonable voxel resolutions and minimum-count thresholds. The precise percentages may change; the key question is whether conclusions such as “most retained regions are locally participant-dominated” persist.

8.6 Emotion-specific occupancy atlas

emotion_atlas_rows = []

for emotion in emotion_order:
    mask = emotions == emotion
    Ze = Z[mask]
    Se = subjects[mask]

    if len(Ze) == 0:
        continue

    ev = build_voxel_table(
        Ze,
        Se,
        bins=VOXEL_BINS,
        min_count=max(8, VOXEL_MIN_COUNT // 2),
    )

    if len(ev) == 0:
        continue

    emotion_atlas_rows.append({
        "emotion": emotion,
        "frames": int(mask.sum()),
        "retained_voxels": len(ev),
        "singleton_voxel_percent": float(
            100 * np.mean(ev["participant_count"] == 1)
        ),
        "mean_participants_per_voxel": float(
            ev["participant_count"].mean()
        ),
        "mean_normalised_participant_entropy": float(
            ev["participant_entropy_normalised"].mean()
        ),
        "median_normalised_participant_entropy": float(
            ev["participant_entropy_normalised"].median()
        ),
        "mean_dominant_share": float(
            ev["dominant_share"].mean()
        ),
        "dominant_share_ge_0_75_percent": float(
            100 * np.mean(ev["dominant_share"] >= 0.75)
        ),
    })

emotion_atlas_summary = (
    pd.DataFrame(emotion_atlas_rows)
    .sort_values("singleton_voxel_percent", ascending=False)
)

emotion_atlas_summary

	emotion	frames	retained_voxels	singleton_voxel_percent	mean_participants_per_voxel	mean_normalised_participant_entropy	median_normalised_participant_entropy	mean_dominant_share	dominant_share_ge_0_75_percent
4	flow	23390	443	97.065463	1.029345	0.007081	0.0	0.993172	98.419865
5	happiness	23327	510	93.529412	1.070588	0.014934	0.0	0.986444	97.450980
7	sadness	28301	661	92.586989	1.077156	0.014647	0.0	0.987161	97.881997
3	fear	31442	752	92.021277	1.082447	0.022214	0.0	0.978301	95.345745
0	anger	27466	633	91.469194	1.094787	0.021967	0.0	0.979758	96.366509
1	contempt	25408	638	90.909091	1.109718	0.024055	0.0	0.978525	96.081505
6	neutral	27744	687	90.829694	1.094614	0.024547	0.0	0.975786	94.905386
2	disgust	22535	555	89.009009	1.124324	0.030249	0.0	0.971510	94.414414
8	surprise	18961	473	89.006342	1.124736	0.024376	0.0	0.978708	96.405920

emotion_atlas_csv = OUT_DIR / "emotion_specific_shared_atlas_summary.csv"
emotion_atlas_summary.to_csv(emotion_atlas_csv, index=False)
print("Wrote:", emotion_atlas_csv.resolve())

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/emotion_specific_shared_atlas_summary.csv

fig_emotion_exclusivity = go.Figure(go.Bar(
    x=emotion_atlas_summary["emotion"].str.title(),
    y=emotion_atlas_summary["singleton_voxel_percent"],
    customdata=np.column_stack([
        emotion_atlas_summary["mean_normalised_participant_entropy"],
        emotion_atlas_summary["mean_dominant_share"],
    ]),
    hovertemplate=(
        "<b>%{x}</b><br>"
        "Singleton voxels: %{y:.1f}%<br>"
        "Mean participant entropy: %{customdata[0]:.3f}<br>"
        "Mean dominant share: %{customdata[1]:.3f}"
        "<extra></extra>"
    ),
))

fig_emotion_exclusivity.update_layout(
    title="Participant-specific local territory by emotion",
    xaxis_title="Elicitation state",
    yaxis_title="Singleton occupied voxels (%)",
    template="plotly_white",
    width=950,
    height=520,
)

emotion_exclusivity_out = OUT_DIR / "emotion_specific_exclusive_territory.html"
fig_emotion_exclusivity.write_html(
    emotion_exclusivity_out,
    include_plotlyjs="cdn",
)
print("Wrote:", emotion_exclusivity_out.resolve())

fig_emotion_exclusivity.show()

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/emotion_specific_exclusive_territory.html

8.7 Dominant-participant ownership map

Each retained voxel is coloured by the participant contributing the largest number of frames. Plotly's 3D scatter markers accept only one scalar opacity per trace, so local dominance strength is encoded through marker size while the exact dominant share remains available on hover.

ownership_colours = [
    PARTICIPANT_COLOURS[p]
    for p in voxels["dominant_participant"]
]

dominance_scaled_size = np.clip(
    (np.sqrt(voxels["count"].to_numpy()) / 2.5)
    * (0.55 + voxels["dominant_share"].to_numpy()),
    2.5,
    11,
)

fig_ownership_v2 = go.Figure()

fig_ownership_v2.add_trace(go.Scatter3d(
    x=voxels["x"],
    y=voxels["y"],
    z=voxels["z"],
    mode="markers",
    name="Dominant local territory",
    showlegend=False,
    marker=dict(
        size=dominance_scaled_size,
        color=ownership_colours,
        opacity=0.72,
        line=dict(width=0),
    ),
    customdata=np.column_stack([
        voxels["count"],
        voxels["participant_count"],
        voxels["dominant_participant"],
        voxels["dominant_share"],
        voxels["participant_entropy_normalised"],
    ]),
    hovertemplate=(
        "Frames: %{customdata[0]}<br>"
        "Participants: %{customdata[1]}<br>"
        "Dominant participant: %{customdata[2]}<br>"
        "Dominant share: %{customdata[3]:.1%}<br>"
        "Normalised entropy: %{customdata[4]:.3f}"
        "<extra></extra>"
    )
))

for participant in participant_order:
    fig_ownership_v2.add_trace(go.Scatter3d(
        x=[None],
        y=[None],
        z=[None],
        mode="markers",
        marker=dict(
            size=7,
            color=PARTICIPANT_COLOURS[participant],
        ),
        name=participant,
        showlegend=True,
    ))

fig_ownership_v2.update_layout(
    title=dict(
        text="Dominant participant in each occupied local region",
        x=0.5,
        xanchor="center",
    ),
    scene=clean_scene("black"),
    paper_bgcolor="black",
    font=dict(color="white", family="Arial, sans-serif"),
    legend=dict(
        title="Dominant participant",
        bgcolor="rgba(0,0,0,0)",
    ),
    width=1150,
    height=850,
    margin=dict(l=0, r=0, t=70, b=20),
)

ownership_v2_out = OUT_DIR / "shared_atlas_dominant_participant_v2.html"
fig_ownership_v2.write_html(
    ownership_v2_out,
    include_plotlyjs="cdn",
)
print("Wrote:", ownership_v2_out.resolve())

fig_ownership_v2.show()

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/shared_atlas_dominant_participant_v2.html

10. Local nearest-neighbour purity

Voxel occupancy depends on the selected grid resolution. A complementary, boundary-free measure asks:

For each frame, what proportion of its k nearest neighbours in the shared UMAP belong to the same participant?

High purity indicates locally participant-specific organisation even when trajectories appear globally intermixed.

This section reports:

• overall participant-neighbour purity for several values of k;
• per-participant purity;
• emotion-neighbour purity for comparison;
• a shuffled-label baseline;
• optional bootstrap confidence intervals.

The calculation is performed in the shared 3D UMAP, so it measures local organisation in that embedding rather than in the original 76-dimensional feature space.

from sklearn.neighbors import NearestNeighbors

# Use all frames when memory permits. Reduce this on machines with limited RAM.
NN_MAX_FRAMES = 120_000
NN_K_VALUES = [5, 10, 25, 50]
NN_METRIC = "euclidean"
NN_RANDOM_SEED = 42
NN_BOOTSTRAP_REPS = 500

nn_rng = np.random.default_rng(NN_RANDOM_SEED)

if len(Z) > NN_MAX_FRAMES:
    # Stratify by participant × emotion so all groups remain represented.
    nn_indices_parts = []
    groups = [
        (participant, emotion)
        for participant in participant_order
        for emotion in emotion_order
        if np.any((subjects == participant) & (emotions == emotion))
    ]
    target_per_group = max(1, NN_MAX_FRAMES // len(groups))

    for participant, emotion in groups:
        idx = np.flatnonzero(
            (subjects == participant) & (emotions == emotion)
        )
        n = min(target_per_group, len(idx))
        nn_indices_parts.append(
            nn_rng.choice(idx, size=n, replace=False)
        )

    nn_indices = np.concatenate(nn_indices_parts)
    nn_indices.sort()
else:
    nn_indices = np.arange(len(Z))

Z_nn = Z[nn_indices]
S_nn = subjects[nn_indices]
E_nn = emotions[nn_indices]

print(f"kNN sample: {len(Z_nn):,} of {len(Z):,} frames")
print(pd.Series(S_nn).value_counts().sort_index())

kNN sample: 118,985 of 228,574 frames
P1    17136
P2    16626
P3    17136
P4    17136
P5    17136
P6    16679
P7    17136
Name: count, dtype: int64

max_k = max(NN_K_VALUES)

nn_model = NearestNeighbors(
    n_neighbors=max_k + 1,  # +1 because the first neighbour is the point itself
    metric=NN_METRIC,
    n_jobs=-1,
)

nn_model.fit(Z_nn)
nn_distances, nn_indices_local = nn_model.kneighbors(Z_nn)

# Remove self-neighbour in column 0
nn_indices_local = nn_indices_local[:, 1:]
nn_distances = nn_distances[:, 1:]

print("Neighbour matrix:", nn_indices_local.shape)

Neighbour matrix: (118985, 50)

def bootstrap_mean_ci(values, reps=500, seed=42, ci=0.95):
    values = np.asarray(values, dtype=float)
    rng = np.random.default_rng(seed)

    means = np.empty(reps, dtype=float)
    n = len(values)

    for r in range(reps):
        sample = rng.choice(values, size=n, replace=True)
        means[r] = sample.mean()

    alpha = (1.0 - ci) / 2.0
    return (
        float(values.mean()),
        float(np.quantile(means, alpha)),
        float(np.quantile(means, 1.0 - alpha)),
    )


purity_rows = []
frame_level_purity = {}

for k in NN_K_VALUES:
    neighbours = nn_indices_local[:, :k]

    participant_matches = (
        S_nn[neighbours] == S_nn[:, None]
    )
    emotion_matches = (
        E_nn[neighbours] == E_nn[:, None]
    )

    participant_purity = participant_matches.mean(axis=1)
    emotion_purity = emotion_matches.mean(axis=1)

    frame_level_purity[k] = {
        "participant": participant_purity,
        "emotion": emotion_purity,
    }

    p_mean, p_low, p_high = bootstrap_mean_ci(
        participant_purity,
        reps=NN_BOOTSTRAP_REPS,
        seed=NN_RANDOM_SEED + k,
    )
    e_mean, e_low, e_high = bootstrap_mean_ci(
        emotion_purity,
        reps=NN_BOOTSTRAP_REPS,
        seed=NN_RANDOM_SEED + 1000 + k,
    )

    purity_rows.append({
        "k": k,
        "participant_purity_mean": p_mean,
        "participant_purity_ci_low": p_low,
        "participant_purity_ci_high": p_high,
        "emotion_purity_mean": e_mean,
        "emotion_purity_ci_low": e_low,
        "emotion_purity_ci_high": e_high,
    })

knn_purity_summary = pd.DataFrame(purity_rows)
knn_purity_summary

	k	participant_purity_mean	participant_purity_ci_low	participant_purity_ci_high	emotion_purity_mean	emotion_purity_ci_low	emotion_purity_ci_high
0	5	0.987550	0.987074	0.988042	0.976200	0.975540	0.976810
1	10	0.966931	0.966218	0.967591	0.949944	0.949057	0.950834
2	25	0.882586	0.881421	0.883859	0.850717	0.849322	0.851970
3	50	0.758862	0.757186	0.760453	0.711859	0.710301	0.713297

10.1 Shuffled-label baseline

The baseline preserves the geometry and class frequencies but randomly permutes participant and emotion labels. This estimates the purity expected if labels were unrelated to local neighbourhood structure.

shuffled_rows = []

S_shuffled = nn_rng.permutation(S_nn)
E_shuffled = nn_rng.permutation(E_nn)

for k in NN_K_VALUES:
    neighbours = nn_indices_local[:, :k]

    participant_baseline = (
        S_shuffled[neighbours] == S_shuffled[:, None]
    ).mean(axis=1)

    emotion_baseline = (
        E_shuffled[neighbours] == E_shuffled[:, None]
    ).mean(axis=1)

    shuffled_rows.append({
        "k": k,
        "participant_purity_shuffled": participant_baseline.mean(),
        "emotion_purity_shuffled": emotion_baseline.mean(),
    })

knn_shuffled_baseline = pd.DataFrame(shuffled_rows)

knn_purity_comparison = knn_purity_summary.merge(
    knn_shuffled_baseline,
    on="k",
    how="left",
)

knn_purity_comparison["participant_purity_above_baseline"] = (
    knn_purity_comparison["participant_purity_mean"]
    - knn_purity_comparison["participant_purity_shuffled"]
)

knn_purity_comparison["emotion_purity_above_baseline"] = (
    knn_purity_comparison["emotion_purity_mean"]
    - knn_purity_comparison["emotion_purity_shuffled"]
)

knn_purity_comparison

	k	participant_purity_mean	participant_purity_ci_low	participant_purity_ci_high	emotion_purity_mean	emotion_purity_ci_low	emotion_purity_ci_high	participant_purity_shuffled	emotion_purity_shuffled	participant_purity_above_baseline	emotion_purity_above_baseline
0	5	0.987550	0.987074	0.988042	0.976200	0.975540	0.976810	0.142956	0.111204	0.844594	0.864996
1	10	0.966931	0.966218	0.967591	0.949944	0.949057	0.950834	0.142840	0.111484	0.824091	0.838460
2	25	0.882586	0.881421	0.883859	0.850717	0.849322	0.851970	0.143150	0.111298	0.739436	0.739419
3	50	0.758862	0.757186	0.760453	0.711859	0.710301	0.713297	0.143124	0.111412	0.615738	0.600447

knn_summary_csv = OUT_DIR / "knn_local_purity_summary.csv"
knn_purity_comparison.to_csv(knn_summary_csv, index=False)
print("Wrote:", knn_summary_csv.resolve())

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/knn_local_purity_summary.csv

fig_knn = go.Figure()

fig_knn.add_trace(go.Scatter(
    x=knn_purity_comparison["k"],
    y=100 * knn_purity_comparison["participant_purity_mean"],
    mode="lines+markers",
    name="Participant purity",
    error_y=dict(
        type="data",
        symmetric=False,
        array=100 * (
            knn_purity_comparison["participant_purity_ci_high"]
            - knn_purity_comparison["participant_purity_mean"]
        ),
        arrayminus=100 * (
            knn_purity_comparison["participant_purity_mean"]
            - knn_purity_comparison["participant_purity_ci_low"]
        ),
    ),
))

fig_knn.add_trace(go.Scatter(
    x=knn_purity_comparison["k"],
    y=100 * knn_purity_comparison["participant_purity_shuffled"],
    mode="lines+markers",
    name="Participant shuffled baseline",
    line=dict(dash="dash"),
))

fig_knn.add_trace(go.Scatter(
    x=knn_purity_comparison["k"],
    y=100 * knn_purity_comparison["emotion_purity_mean"],
    mode="lines+markers",
    name="Emotion purity",
))

fig_knn.add_trace(go.Scatter(
    x=knn_purity_comparison["k"],
    y=100 * knn_purity_comparison["emotion_purity_shuffled"],
    mode="lines+markers",
    name="Emotion shuffled baseline",
    line=dict(dash="dash"),
))

fig_knn.update_layout(
    title="Local nearest-neighbour purity in the shared UMAP",
    xaxis_title="Number of nearest neighbours (k)",
    yaxis_title="Same-label neighbours (%)",
    template="plotly_white",
    width=950,
    height=560,
)

knn_plot_out = OUT_DIR / "knn_local_purity.html"
fig_knn.write_html(knn_plot_out, include_plotlyjs="cdn")
print("Wrote:", knn_plot_out.resolve())

fig_knn.show()

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/knn_local_purity.html

10.2 Per-participant purity

PER_PARTICIPANT_K = 25

participant_purity_k = frame_level_purity[PER_PARTICIPANT_K]["participant"]

per_participant_rows = []

for participant in participant_order:
    mask = S_nn == participant
    values = participant_purity_k[mask]

    mean, low, high = bootstrap_mean_ci(
        values,
        reps=NN_BOOTSTRAP_REPS,
        seed=NN_RANDOM_SEED + int(participant[1:]),
    )

    per_participant_rows.append({
        "participant": participant,
        "frames": int(mask.sum()),
        "k": PER_PARTICIPANT_K,
        "mean_same_participant_neighbour_purity": mean,
        "ci_low": low,
        "ci_high": high,
        "median_purity": float(np.median(values)),
        "frames_with_purity_ge_0_75_percent": float(
            100 * np.mean(values >= 0.75)
        ),
        "frames_with_purity_eq_1_percent": float(
            100 * np.mean(values == 1.0)
        ),
    })

per_participant_purity = (
    pd.DataFrame(per_participant_rows)
    .sort_values(
        "mean_same_participant_neighbour_purity",
        ascending=False,
    )
)

per_participant_purity

	participant	frames	k	mean_same_participant_neighbour_purity	ci_low	ci_high	median_purity	frames_with_purity_ge_0_75_percent	frames_with_purity_eq_1_percent
3	P4	17136	25	0.932085	0.929616	0.934731	1.0	87.914332	80.287115
5	P6	16679	25	0.918530	0.915784	0.921076	1.0	87.277415	73.997242
2	P3	17136	25	0.909914	0.907145	0.912636	1.0	84.535481	72.619048
1	P2	16626	25	0.886296	0.882937	0.889693	1.0	80.753037	72.705401
0	P1	17136	25	0.855852	0.852448	0.859102	1.0	76.184641	57.913165
6	P7	17136	25	0.840075	0.836530	0.843613	1.0	73.739496	58.327498
4	P5	17136	25	0.836422	0.832967	0.840103	1.0	73.150093	55.345472

per_participant_csv = OUT_DIR / "knn_purity_by_participant.csv"
per_participant_purity.to_csv(per_participant_csv, index=False)
print("Wrote:", per_participant_csv.resolve())

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/knn_purity_by_participant.csv

fig_participant_purity = go.Figure(go.Bar(
    x=per_participant_purity["participant"],
    y=100 * per_participant_purity["mean_same_participant_neighbour_purity"],
    error_y=dict(
        type="data",
        symmetric=False,
        array=100 * (
            per_participant_purity["ci_high"]
            - per_participant_purity["mean_same_participant_neighbour_purity"]
        ),
        arrayminus=100 * (
            per_participant_purity["mean_same_participant_neighbour_purity"]
            - per_participant_purity["ci_low"]
        ),
    ),
    customdata=np.column_stack([
        per_participant_purity["median_purity"],
        per_participant_purity["frames_with_purity_ge_0_75_percent"],
        per_participant_purity["frames_with_purity_eq_1_percent"],
    ]),
    hovertemplate=(
        "<b>%{x}</b><br>"
        "Mean purity: %{y:.1f}%<br>"
        "Median purity: %{customdata[0]:.3f}<br>"
        "Frames ≥75% pure: %{customdata[1]:.1f}%<br>"
        "Frames 100% pure: %{customdata[2]:.1f}%"
        "<extra></extra>"
    ),
))

fig_participant_purity.update_layout(
    title=f"Participant-local neighbourhood purity (k={PER_PARTICIPANT_K})",
    xaxis_title="Participant",
    yaxis_title="Same-participant neighbours (%)",
    template="plotly_white",
    width=900,
    height=520,
)

participant_purity_out = OUT_DIR / "knn_purity_by_participant.html"
fig_participant_purity.write_html(
    participant_purity_out,
    include_plotlyjs="cdn",
)
print("Wrote:", participant_purity_out.resolve())

fig_participant_purity.show()

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/knn_purity_by_participant.html

10.3 Per-emotion participant purity

# This asks whether particular elicitation states occupy especially
# participant-specific neighbourhoods in the shared embedding.
emotion_participant_rows = []

for emotion in emotion_order:
    mask = E_nn == emotion
    values = participant_purity_k[mask]

    mean, low, high = bootstrap_mean_ci(
        values,
        reps=NN_BOOTSTRAP_REPS,
        seed=NN_RANDOM_SEED + 2000 + emotion_order.index(emotion),
    )

    emotion_participant_rows.append({
        "emotion": emotion,
        "frames": int(mask.sum()),
        "k": PER_PARTICIPANT_K,
        "mean_same_participant_neighbour_purity": mean,
        "ci_low": low,
        "ci_high": high,
        "median_purity": float(np.median(values)),
        "frames_with_purity_ge_0_75_percent": float(
            100 * np.mean(values >= 0.75)
        ),
    })

purity_by_emotion = (
    pd.DataFrame(emotion_participant_rows)
    .sort_values(
        "mean_same_participant_neighbour_purity",
        ascending=False,
    )
)

purity_by_emotion

	emotion	frames	k	mean_same_participant_neighbour_purity	ci_low	ci_high	median_purity	frames_with_purity_ge_0_75_percent
4	flow	13328	25	0.935105	0.932087	0.937936	1.0	89.503301
5	happiness	13328	25	0.904058	0.900590	0.907269	1.0	83.155762
8	surprise	12361	25	0.896309	0.892886	0.900222	1.0	84.523906
2	disgust	13328	25	0.884256	0.880455	0.887437	1.0	79.891957
1	contempt	13328	25	0.879658	0.875848	0.883224	1.0	80.462185
7	sadness	13328	25	0.877659	0.874312	0.881391	1.0	79.606843
0	anger	13328	25	0.864499	0.860693	0.868475	1.0	77.235894
6	neutral	13328	25	0.857671	0.853607	0.861244	1.0	76.418067
3	fear	13328	25	0.845057	0.840955	0.849831	1.0	73.822029

purity_by_emotion_csv = OUT_DIR / "knn_participant_purity_by_emotion.csv"
purity_by_emotion.to_csv(purity_by_emotion_csv, index=False)
print("Wrote:", purity_by_emotion_csv.resolve())

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/knn_participant_purity_by_emotion.csv

fig_emotion_purity = go.Figure(go.Bar(
    x=purity_by_emotion["emotion"].str.title(),
    y=100 * purity_by_emotion["mean_same_participant_neighbour_purity"],
    error_y=dict(
        type="data",
        symmetric=False,
        array=100 * (
            purity_by_emotion["ci_high"]
            - purity_by_emotion["mean_same_participant_neighbour_purity"]
        ),
        arrayminus=100 * (
            purity_by_emotion["mean_same_participant_neighbour_purity"]
            - purity_by_emotion["ci_low"]
        ),
    ),
    customdata=np.column_stack([
        purity_by_emotion["median_purity"],
        purity_by_emotion["frames_with_purity_ge_0_75_percent"],
    ]),
    hovertemplate=(
        "<b>%{x}</b><br>"
        "Mean participant purity: %{y:.1f}%<br>"
        "Median: %{customdata[0]:.3f}<br>"
        "Frames ≥75% pure: %{customdata[1]:.1f}%"
        "<extra></extra>"
    ),
))

fig_emotion_purity.update_layout(
    title=f"Participant-specific neighbourhood structure by emotion (k={PER_PARTICIPANT_K})",
    xaxis_title="Elicitation state",
    yaxis_title="Same-participant neighbours (%)",
    template="plotly_white",
    width=980,
    height=540,
)

emotion_purity_out = OUT_DIR / "knn_participant_purity_by_emotion.html"
fig_emotion_purity.write_html(
    emotion_purity_out,
    include_plotlyjs="cdn",
)
print("Wrote:", emotion_purity_out.resolve())

fig_emotion_purity.show()

Wrote: /Users/macart/Desktop/jupyter/face/shared_umap_76D_SOT_exports/shared_atlas_outputs/knn_participant_purity_by_emotion.html

Interpretation

Nearest-neighbour purity complements voxel occupancy:

• Voxel occupancy asks how many participants visit a local spatial bin.
• kNN purity asks how strongly each frame's immediate neighbourhood is dominated by its own participant.

A high observed purity relative to the shuffled baseline supports participant-specific local organisation without depending on arbitrary voxel boundaries. Report several values of k rather than selecting only the most favourable one.