import warnings

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats

import vmhvl

warnings.filterwarnings("ignore")
plt.rcParams.update({"figure.dpi": 100, "axes.spines.top": False, "axes.spines.right": False})
pd.set_option("display.width", 160)
pd.set_option("display.max_columns", 30)

versions = ["0.240814.1849", "0.260531.0005", "0.260624.0330"]
tables = {v: vmhvl.session_table(v) for v in versions}
for v, t in tables.items():
    print(f"{v}: {len(t)} files, {t.subject.nunique()} subject labels, {t['size'].sum() / 1e6:.0f} MB")
same = tables["0.260531.0005"][["path", "asset_id"]].equals(tables["0.260624.0330"][["path", "asset_id"]])
print("0.260531.0005 and 0.260624.0330 hold the same files with the same asset ids:", same)
old = tables["0.240814.1849"]
print("0.240814.1849 holds:", old.path.tolist())
print("  the same paths in 0.260624.0330 have the same asset ids:",
      old.set_index("path").asset_id.equals(tables["0.260624.0330"].set_index("path").asset_id.loc[old.path]))
print("  and the same sizes:",
      old.set_index("path")["size"].equals(tables["0.260624.0330"].set_index("path")["size"].loc[old.path]))
for p in old.path:
    a, b = vmhvl.digest(p, "0.240814.1849"), vmhvl.digest(p, "0.260624.0330")
    print(f"  {p}: SHA-256 {a[:12]}... in 0.240814.1849 and {b[:12]}... in 0.260624.0330; identical {a == b}")

0.240814.1849: 2 files, 2 subject labels, 16 MB
0.260531.0005: 107 files, 26 subject labels, 414 MB
0.260624.0330: 107 files, 26 subject labels, 414 MB
0.260531.0005 and 0.260624.0330 hold the same files with the same asset ids: True
0.240814.1849 holds: ['sub-m541/sub-m541_behavior.nwb', 'sub-m542/sub-m542_behavior.nwb']
  the same paths in 0.260624.0330 have the same asset ids: False
  and the same sizes: True

  sub-m541/sub-m541_behavior.nwb: SHA-256 354c77d3bc10... in 0.240814.1849 and 354c77d3bc10... in 0.260624.0330; identical True

  sub-m542/sub-m542_behavior.nwb: SHA-256 a104d039f845... in 0.240814.1849 and a104d039f845... in 0.260624.0330; identical True

files = tables[vmhvl.VERSION]
sessions_t = files[files.session.notna()].copy()
annotations = vmhvl.load_many(sessions_t.path.tolist(), neural=False)
sessions_t["n_frames"] = [annotations[p].n_frames for p in sessions_t.path]
sessions_t["n_cells"] = [annotations[p].n_cells for p in sessions_t.path]
sessions_t["kind"] = [annotations[p].kind for p in sessions_t.path]
sessions_t["female_annotated"] = [annotations[p].has_female() for p in sessions_t.path]
sessions_t["receptive"] = [annotations[p].receptive() for p in sessions_t.path]

receptive_paths = sessions_t[sessions_t.experiment == "receptive"].path.tolist()
oestrogen_paths = sessions_t[sessions_t.experiment == "oestrogen-injected (7 days)"].path.tolist()
cycle_paths = sessions_t[(sessions_t.experiment == "natural cycle (3 days)")
                         & sessions_t.subject.isin(["m541", "543"])].path.tolist()
duplicate_pairs = [
    ("sub-m764_ses-WT_3d_mating_d1_ecephys+behavior.nwb", "sub-m764_ses-ovx_primed_receptive_mating_ecephys+behavior.nwb"),
    ("sub-m776_ses-ovx_primed_receptive_mating_ecephys+behavior.nwb", "sub-m776_ses-wireless_opto_ecephys+behavior.nwb"),
    ("sub-m551_ses-OVX_3d_mating_d3_ecephys+behavior.nwb", "sub-m551_ses-ovx_primed_receptive_mating_ecephys+behavior.nwb"),
]
m438_days = [f"sub-m438_ses-WT_3d_mating_d{d}_ecephys+behavior.nwb" for d in (1, 2, 3)]
same_frames = [f"sub-{s}_ses-{e}_ecephys+behavior.nwb" for s in ("m438", "m541")
               for e in ("OVX_3d_mating_d3", "ovx_primed_receptive_mating")]
neural_paths = sorted(set(receptive_paths + oestrogen_paths + cycle_paths + m438_days + same_frames + [p for pair in duplicate_pairs for p in pair]))
data = vmhvl.load_many(neural_paths, neural=True)
weights = {s: vmhvl.load_weights(s) for s in ("m541", "m542")}
print(f"annotations read from {len(annotations)} files, traces from {len(data)} files")

annotations read from 105 files, traces from 76 files

# The behaviour raster marks frames from an annotation's start up to its stop, not including it.
test = pd.DataFrame({"channel": ["Ch2", "Ch2"], "behavior": ["x", "x"], "start": [1.0, 3.02], "stop": [1.5, 3.03]})
assert vmhvl.raster(test, ["x"], 50).nonzero()[0].tolist() == [10, 11, 12, 13, 14, 30]
print("Files by experiment (sessions, distinct subject labels, cells per session, frames per session):")
summary = sessions_t.groupby("experiment").agg(
    sessions=("path", "size"), subjects=("subject", "nunique"),
    cells_min=("n_cells", "min"), cells_max=("n_cells", "max"),
    frames_min=("n_frames", "min"), frames_max=("n_frames", "max"),
    female_annotated=("female_annotated", "sum"), receptive=("receptive", "sum"))
print(summary.to_string())
print()
print("Subjects per experiment:")
for exp, g in sessions_t.groupby("experiment"):
    print(f"  {exp}: {', '.join(sorted(g.subject.unique()))}")
print()
print("Trace type stored in each file, by experiment (Fluorescence is raw, DfOverF is described as z-scored):")
print(pd.crosstab(sessions_t.experiment, sessions_t.kind).to_string())
print()
print("Stored traces:")
for kind, g in sessions_t.groupby("kind"):
    s = annotations[g.path.iloc[0]]
    print(f"  {kind}: {len(g)} files, unit '{s.unit}', description '{s.description[:60]}...'")
for kind in ("Fluorescence", "DfOverF"):
    ps = [p for p in neural_paths if data[p].kind == kind]
    allv = np.concatenate([data[p].traces.ravel() for p in ps])
    means = np.concatenate([data[p].traces.mean(0) for p in ps])
    sds = np.concatenate([data[p].traces.std(0) for p in ps])
    print(f"  {kind} (loaded {len(ps)} files): values {allv.min():.1f} to {allv.max():.1f}; "
          f"per-cell mean {means.min():.2f} to {means.max():.2f}; per-cell SD {sds.min():.2f} to {sds.max():.2f}; "
          f"NaN {np.isnan(allv).sum()}")
dts = np.concatenate([np.diff(data[p].timestamps) for p in neural_paths])
starts = [data[p].timestamps[0] for p in neural_paths]
print(f"  timestamps: first {min(starts):.3f} to {max(starts):.3f} s; steps {dts.min():.4f} to {dts.max():.4f} s "
      f"(rate {1 / np.median(dts):.2f} Hz)")
print()
b_all = pd.concat([a.behaviour.assign(path=p) for p, a in annotations.items()])
on_grid = np.abs(b_all[["start", "stop"]].to_numpy() * 30 - np.round(b_all[["start", "stop"]].to_numpy() * 30))
print(f"Annotations: {len(b_all)} intervals; stop before start in {(b_all.stop < b_all.start).sum()}; "
      f"zero-length {(b_all.stop == b_all.start).sum()}; times off the 30-Hz grid by more than 1 ms: {(on_grid > 0.03).sum()}")
ends = pd.Series({p: annotations[p].behaviour.stop.max() - annotations[p].n_frames / vmhvl.RATE for p in annotations})
print(f"  last annotation minus end of imaging: {ends.min():.1f} to {ends.max():.1f} s")
print(f"  sessions whose imaging runs more than 60 s past the last annotation: "
      f"{', '.join(p.split('_ecephys')[0] for p in ends[ends < -60].index)}")
entry = pd.Series({p: annotations[p].interaction()[0] for p in annotations})
print(f"  male entry (into_male_cage) at {entry.min():.0f} to {entry.max():.0f} s after the start of imaging")
print("  labels, Ch1 (male):", sorted(b_all[b_all.channel == "Ch1"].behavior.unique()))
print("  labels, Ch2 (female):", sorted(b_all[b_all.channel == "Ch2"].behavior.unique()))
overlap = []
for p, a in annotations.items():
    f = a.behaviour[(a.behaviour.channel == "Ch2") & (a.behaviour.behavior != "into_male_cage")].sort_values("start")
    overlap.append(np.clip(f.stop.to_numpy()[:-1] - f.start.to_numpy()[1:], 0, None).sum())
print(f"  overlap between consecutive female annotations, summed over all sessions: {sum(overlap):.1f} s")
print()
print("Descriptive metadata, distinct values over the session files:")
meta = pd.DataFrame([a.meta for a in annotations.values()])
for col in meta.columns:
    print(f"  {col}: {sorted(meta[col].astype(str).unique())}")

Files by experiment (sessions, distinct subject labels, cells per session, frames per session):
                             sessions  subjects  cells_min  cells_max  frames_min  frames_max  female_annotated  receptive
experiment                                                                                                                
male inhibition                     8         8         43        128        7386       18006                 0          0
natural cycle (3 days)             21         7         66        176        5327       15644                 3          3
oestrogen-injected (7 days)        42         6         83        249        4360        6133                42         31
ovariectomised, primed             10        10         56        184        4216       12404                10         10
ovariectomy (3 days)                9         3         64        184        4216       10191                 0          0
receptive                          15        15         16        207        4511       15769                15         15

Subjects per experiment:
  male inhibition: m764, m774, m775, m776, m858, m918, m931, m932
  natural cycle (3 days): 543, m438, m541, m542, m551, m610, m764
  oestrogen-injected (7 days): m918, m928, m929, m930, m931, m940
  ovariectomised, primed: m438, m449, m541, m551, m764, m774, m775, m776, m858, m859
  ovariectomy (3 days): m438, m541, m551
  receptive: m438, m540, m541, m542, m543, m546, m551, m553, m857, m858, m859, m864, m928, m930, m940

Trace type stored in each file, by experiment (Fluorescence is raw, DfOverF is described as z-scored):
kind                         DfOverF  Fluorescence
experiment                                        
male inhibition                    0             8
natural cycle (3 days)            21             0
oestrogen-injected (7 days)       42             0
ovariectomised, primed             0            10
ovariectomy (3 days)               9             0
receptive                          0            15

Stored traces:
  DfOverF: 72 files, unit 'n.a.', description 'Z-scored fluorescence traces. Shape: (frames=10527, neurons=...'
  Fluorescence: 33 files, unit 'lumens', description 'Raw fluorescence traces. Shape: (frames=4216, neurons=64). S...'
  Fluorescence (loaded 21 files): values -5476.8 to 8110.5; per-cell mean 1.82 to 1417.48; per-cell SD 12.02 to 1493.29; NaN 0

  DfOverF (loaded 55 files): values -6.6 to 16.1; per-cell mean -1.61 to 2.35; per-cell SD 0.10 to 1.78; NaN 0

  timestamps: first 0.000 to 0.000 s; steps 0.1000 to 0.1000 s (rate 10.00 Hz)

Annotations: 17983 intervals; stop before start in 0; zero-length 0; times off the 30-Hz grid by more than 1 ms: 0
  last annotation minus end of imaging: -426.3 to 0.0 s
  sessions whose imaging runs more than 60 s past the last annotation: sub-m540_ses-wt_proestrus_receptive_mating, sub-m541_ses-wt_proestrus_receptive_mating, sub-m542_ses-wt_proestrus_receptive_mating, sub-m543_ses-wt_proestrus_receptive_mating, sub-m546_ses-wt_proestrus_receptive_mating, sub-m551_ses-wt_proestrus_receptive_mating, sub-m553_ses-wt_proestrus_receptive_mating, sub-m764_ses-WT_3d_mating_d3, sub-m764_ses-wireless_opto, sub-m775_ses-ovx_primed_receptive_mating, sub-m858_ses-wireless_opto, sub-m859_ses-ovx_primed_receptive_mating, sub-m928_ses-WT_7d_mating_d2, sub-m928_ses-WT_7d_mating_d5, sub-m931_ses-WT_7d_mating_d4, sub-m931_ses-wireless_opto, sub-m932_ses-wireless_opto
  male entry (into_male_cage) at 123 to 233 s after the start of imaging
  labels, Ch1 (male): ['ejaculation', 'grooming', 'into_male_cage', 'intromission', 'laser_on', 'male_sniff', 'mating_inhibition', 'mount', 'other', 'sniff']
  labels, Ch2 (female): ['approach', 'check_genital', 'dart', 'dig', 'eat', 'into_male_cage', 'intromission', 'kick', 'lordose', 'other', 'rear', 'sniff', 'stay', 'top_up', 'turn', 'wiggle']
  overlap between consecutive female annotations, summed over all sessions: 23.1 s

Descriptive metadata, distinct values over the session files:
  device: ['Two-photon calcium imaging microscope']
  excitation_nm: ['920.0']
  indicator: ['GCaMP']
  location: ['VMH']
  experimenter: ["('Your Name',)"]
  lab: ['Your Lab Name']
  institution: ['Your Institution']
  session_start: ['2026-05-30']
  image_mask: ['(1, 1)']

def match(a, b):
    """Median over cells of the best correlation of each cell of a with any cell of b."""
    n = min(len(a), len(b))
    za, zb = stats.zscore(a[:n], axis=0), stats.zscore(b[:n], axis=0)
    return np.median(np.max(za.T @ zb / n, axis=1))

for p, q in duplicate_pairs:
    a, b = data[p], data[q]
    same_shape = a.traces.shape == b.traces.shape
    r = np.median([np.corrcoef(a.traces[:, i], b.traces[:, i])[0, 1] for i in range(a.n_cells)]) if same_shape else np.nan
    ca = a.behaviour[a.behaviour.channel == "Ch1"].behavior.value_counts().to_dict()
    cb = b.behaviour[b.behaviour.channel == "Ch1"].behavior.value_counts().to_dict()
    print(f"{p.split('_ecephys')[0]} ({a.kind}) and {q.split('_ecephys')[0]} ({b.kind}): shape {a.traces.shape} / {b.traces.shape}; "
          f"median correlation of cell i with cell i {r:.4f}; male annotations identical {ca == cb}")
print()
for s, w in weights.items():
    print(f"per-mouse file sub-{s}: subject_id '{w['subject_id']}', traces {w['traces'].shape}, weight columns {w['columns']}, "
          f"session_description '{w['session_description']}', comments '{w['comments']}'")
    print(f"    its traces: per-cell mean {w['traces'].mean(0).min():.2f} to {w['traces'].mean(0).max():.2f}, "
          f"SD {w['traces'].std(0).min():.2f} to {w['traces'].std(0).max():.2f}")
    for p in cycle_paths:
        d = data[p]
        if d.traces.shape[1] == w["traces"].shape[1]:
            print(f"    {p.split('_ecephys')[0]}: {d.traces.shape}, best-match correlation {match(w['traces'], d.traces):.2f}")
        else:
            print(f"    {p.split('_ecephys')[0]}: {d.traces.shape}, cell counts differ")
m542_days = sessions_t[(sessions_t.subject == "m542") & sessions_t.experiment.str.startswith("natural")]
print(f"    the natural-cycle sessions of m542 have {sorted(set(m542_days.n_cells))} cells")
print()
print("Contents of the per-mouse files under analysis and processing (behaviour epochs summarised):")
for s in weights:
    objects = vmhvl.list_objects(f"sub-{s}/sub-{s}_behavior.nwb")
    epochs = [o for o in objects if "BehavioralEpochs/" in o[0]]
    print(f"  sub-{s}: {len(epochs)} objects under processing/behavior/BehavioralEpochs (annotation series); the rest:")
    for name, shape, kind in objects:
        if "BehavioralEpochs/" not in name:
            print(f"    {name} {shape if shape is not None else ''} {kind}")

sub-m764_ses-WT_3d_mating_d1 (DfOverF) and sub-m764_ses-ovx_primed_receptive_mating (Fluorescence): shape (9077, 138) / (9077, 138); median correlation of cell i with cell i 1.0000; male annotations identical True
sub-m776_ses-ovx_primed_receptive_mating (Fluorescence) and sub-m776_ses-wireless_opto (Fluorescence): shape (12404, 75) / (12404, 75); median correlation of cell i with cell i 1.0000; male annotations identical False
sub-m551_ses-OVX_3d_mating_d3 (DfOverF) and sub-m551_ses-ovx_primed_receptive_mating (Fluorescence): shape (10191, 141) / (10191, 141); median correlation of cell i with cell i 1.0000; male annotations identical False

per-mouse file sub-m541: subject_id 'm541', traces (6451, 176), weight columns ['latent factor 1', 'latent factor 2'], session_description 'resident-intruder assay', comments 'GCaMP6m was used. Units are DF/F.'
    its traces: per-cell mean -0.00 to 0.00, SD 1.00 to 1.00
    sub-543_ses-WT_3d_mating_d1: (10527, 72), cell counts differ
    sub-543_ses-WT_3d_mating_d2: (5327, 72), cell counts differ
    sub-543_ses-WT_3d_mating_d3: (6063, 72), cell counts differ
    sub-m541_ses-WT_3d_mating_d1: (9415, 176), best-match correlation 0.54

    sub-m541_ses-WT_3d_mating_d2: (6415, 176), best-match correlation 0.54

    sub-m541_ses-WT_3d_mating_d3: (6451, 176), best-match correlation 0.86
per-mouse file sub-m542: subject_id 'm542', traces (10527, 72), weight columns ['latent factor 1', 'latent factor 2'], session_description 'resident-intruder assay', comments 'GCaMP6m was used. Units are DF/F.'
    its traces: per-cell mean -0.00 to 0.00, SD 1.00 to 1.00
    sub-543_ses-WT_3d_mating_d1: (10527, 72), best-match correlation 0.80
    sub-543_ses-WT_3d_mating_d2: (5327, 72), best-match correlation 0.49
    sub-543_ses-WT_3d_mating_d3: (6063, 72), best-match correlation 0.42
    sub-m541_ses-WT_3d_mating_d1: (9415, 176), cell counts differ
    sub-m541_ses-WT_3d_mating_d2: (6415, 176), cell counts differ
    sub-m541_ses-WT_3d_mating_d3: (6451, 176), cell counts differ
    the natural-cycle sessions of m542 have [129] cells

Contents of the per-mouse files under analysis and processing (behaviour epochs summarised):

  sub-m541: 42 objects under processing/behavior/BehavioralEpochs (annotation series); the rest:
    analysis  
    analysis/rSLDS  DynamicTable
    analysis/rSLDS/id (176,) ElementIdentifiers
    analysis/rSLDS/latent factor 1 (176,) VectorData
    analysis/rSLDS/latent factor 2 (176,) VectorData
    processing  
    processing/behavior  ProcessingModule
    processing/behavior/BehavioralEpochs  BehavioralEpochs
    processing/ophys  ProcessingModule
    processing/ophys/NeuralTrace  TimeSeries
    processing/ophys/NeuralTrace/data (6451, 176) 
    processing/ophys/NeuralTrace/starting_time ()

  sub-m542: 39 objects under processing/behavior/BehavioralEpochs (annotation series); the rest:
    analysis  
    analysis/rSLDS  DynamicTable
    analysis/rSLDS/id (72,) ElementIdentifiers
    analysis/rSLDS/latent factor 1 (72,) VectorData
    analysis/rSLDS/latent factor 2 (72,) VectorData
    processing  
    processing/behavior  ProcessingModule
    processing/behavior/BehavioralEpochs  BehavioralEpochs
    processing/ophys  ProcessingModule
    processing/ophys/NeuralTrace  TimeSeries
    processing/ophys/NeuralTrace/data (10527, 72) 
    processing/ophys/NeuralTrace/starting_time ()

order_files = {"m541": "sub-m541_ses-WT_3d_mating_d3_ecephys+behavior.nwb",
               "m542": "sub-543_ses-WT_3d_mating_d1_ecephys+behavior.nwb"}
for s, p in order_files.items():
    A, B = weights[s]["traces"], data[p].traces
    n = min(len(A), len(B))
    za, zb = stats.zscore(A[:n], axis=0), stats.zscore(B[:n], axis=0)
    corr = za.T @ zb / n
    same = np.diag(corr)
    print(f"sub-{s} against {p.split('_ecephys')[0]} ({A.shape} and {B.shape}):")
    print(f"  best match of cell i is cell i for {(corr.argmax(1) == np.arange(len(corr))).sum()} of {len(corr)} cells; "
          f"correlation of cell i with cell i: median {np.median(same):.3f}, range {same.min():.2f} to {same.max():.2f}")
    by_shift = {}
    for lag in (-20, -10, -5, -2, 0, 2, 5, 10, 20):
        x, y = (A[lag:n], B[:n - lag]) if lag >= 0 else (A[:n + lag], B[-lag:n])
        by_shift[lag] = round(float(np.median([np.corrcoef(x[:, i], y[:, i])[0, 1] for i in range(A.shape[1])])), 3)
    print(f"  median same-index correlation by shift of the per-mouse traces (frames): {by_shift}")
    sm = {k: round(float(np.median([np.corrcoef(A[:, i], vmhvl.smooth(B, k)[:, i])[0, 1] for i in range(A.shape[1])])), 3)
          for k in (3, 9, 21, 31)}
    print(f"  the same after smoothing the stored traces with a moving average of k frames: {sm}")
    W = weights[s]["weights"]
    for k in range(W.shape[1]):
        pa, pb = A @ W[:, k], B @ W[:, k]
        print(f"  factor {k + 1}: correlation of the projection from the per-mouse traces with that from the stored traces "
              f"{np.corrcoef(pa, pb)[0, 1]:.3f}; ACHW up to 300 s {vmhvl.achw(pa[:, None], max_lag=3000)[0]:.1f} s "
              f"and {vmhvl.achw(pb[:, None], max_lag=3000)[0]:.1f} s")

sub-m541 against sub-m541_ses-WT_3d_mating_d3 ((6451, 176) and (6451, 176)):
  best match of cell i is cell i for 171 of 176 cells; correlation of cell i with cell i: median 0.862, range 0.33 to 0.96
  median same-index correlation by shift of the per-mouse traces (frames): {-20: 0.831, -10: 0.854, -5: 0.861, -2: 0.863, 0: 0.862, 2: 0.862, 5: 0.858, 10: 0.851, 20: 0.828}

  the same after smoothing the stored traces with a moving average of k frames: {3: 0.883, 9: 0.901, 21: 0.922, 31: 0.932}
  factor 1: correlation of the projection from the per-mouse traces with that from the stored traces 0.960; ACHW up to 300 s 62.7 s and 56.8 s
  factor 2: correlation of the projection from the per-mouse traces with that from the stored traces 0.957; ACHW up to 300 s 115.3 s and 104.8 s
sub-m542 against sub-543_ses-WT_3d_mating_d1 ((10527, 72) and (10527, 72)):
  best match of cell i is cell i for 71 of 72 cells; correlation of cell i with cell i: median 0.802, range 0.43 to 0.91
  median same-index correlation by shift of the per-mouse traces (frames): {-20: 0.788, -10: 0.8, -5: 0.802, -2: 0.802, 0: 0.802, 2: 0.802, 5: 0.801, 10: 0.798, 20: 0.782}

  the same after smoothing the stored traces with a moving average of k frames: {3: 0.829, 9: 0.853, 21: 0.881, 31: 0.901}
  factor 1: correlation of the projection from the per-mouse traces with that from the stored traces 0.951; ACHW up to 300 s 168.7 s and 128.9 s
  factor 2: correlation of the projection from the per-mouse traces with that from the stored traces 0.944; ACHW up to 300 s 40.3 s and 32.1 s

def male_onsets(s):
    b = s.behaviour
    return np.sort(b[(b.channel == "Ch1") & b.behavior.isin(vmhvl.MALE_SELF)].start.to_numpy())

receptive_day = {}
cycle = sessions_t[sessions_t.experiment == "natural cycle (3 days)"]
for subject, g in cycle.groupby("subject"):
    rec = sessions_t[(sessions_t.experiment == "receptive") & (sessions_t.subject.str.lstrip("m") == subject.lstrip("m"))]
    if rec.empty:
        print(f"{subject}: no receptive session deposited")
        continue
    r = annotations[rec.path.iloc[0]]
    for _, row in g.sort_values("day").iterrows():
        d = annotations[row.path]
        a, b = male_onsets(d), male_onsets(r)
        same_n = len(a) == len(b)
        offset = np.median(b - a) if same_n else np.nan
        spread = np.ptp(b - a) if same_n else np.nan
        line = f"{subject} day {row.day}: {len(a)} male annotations against {len(b)}"
        if same_n:
            line += f", offset {offset:.2f} s (range {spread:.2f} s)"
        if row.path in data and rec.path.iloc[0] in data:
            line += f", trace best-match correlation {match(data[row.path].traces, data[rec.path.iloc[0]].traces):.2f}"
        print(line)
        if same_n and spread < 0.5:
            receptive_day[subject] = row.day
print("receptive day by matching:", receptive_day)

543 day 1.0: 63 male annotations against 63, offset 0.07 s (range 0.03 s), trace best-match correlation 0.84
543 day 2.0: 93 male annotations against 63, trace best-match correlation 0.44
543 day 3.0: 61 male annotations against 63, trace best-match correlation 0.37
m438 day 1.0: 60 male annotations against 76, trace best-match correlation 0.78
m438 day 2.0: 79 male annotations against 76, trace best-match correlation 0.44
m438 day 3.0: 96 male annotations against 76, trace best-match correlation 0.33
m541 day 1.0: 120 male annotations against 31, trace best-match correlation 0.39
m541 day 2.0: 83 male annotations against 31, trace best-match correlation 0.44
m541 day 3.0: 31 male annotations against 31, offset 0.17 s (range 0.00 s), trace best-match correlation 0.81
m542 day 1.0: 83 male annotations against 82
m542 day 2.0: 82 male annotations against 82, offset 0.17 s (range 0.03 s)
m542 day 3.0: 123 male annotations against 82
m551 day 1.0: 100 male annotations against 100, offset 0.10 s (range 0.03 s)
m551 day 2.0: 75 male annotations against 100
m551 day 3.0: 102 male annotations against 100
m610: no receptive session deposited
m764: no receptive session deposited
receptive day by matching: {'543': 1.0, 'm541': 3.0, 'm542': 2.0, 'm551': 1.0}

rng_clock = np.random.default_rng(1097_00)
shifts = [-60, -30, -10, -3, 0, 3, 10, 30, 60]
clock = []
for p in receptive_paths + oestrogen_paths + cycle_paths:
    s = data[p]
    t0, t1 = s.interaction()
    i0, i1 = int(round(t0 * vmhvl.RATE)), min(int(round(t1 * vmhvl.RATE)), s.n_frames)
    y = vmhvl.bin_frames(vmhvl.zscore(s.traces)[i0:i1], 10)
    lab = vmhvl.bin_frames(vmhvl.raster(s.behaviour, vmhvl.MALE_SELF, s.n_frames, "Ch1")[i0:i1].astype(float), 10) >= 0.5
    acc = vmhvl.shifted_decoding(y, lab, shifts, rng_clock)
    clock.append({"session": p.split("_ecephys")[0], "positives": int(lab.sum()), "bins": len(lab), **dict(zip(shifts, acc))})
clock = pd.DataFrame(clock).set_index("session")
best = clock[shifts].idxmax(axis=1)
near = clock[[-3, 0, 3]].max(axis=1)
far = clock[[-60, 60]].max(axis=1)
print(f"{len(clock)} sessions. Shift of best accuracy: {best.value_counts().sort_index().to_dict()}")
print(f"best accuracy within 3 s of zero: median {near.median():.2f} (range {near.min():.2f} to {near.max():.2f}); "
      f"best at +-60 s: median {far.median():.2f}")
clock_fail = clock.index[(near - far) < 0.05].tolist()
print(f"sessions where accuracy near zero exceeds that at +-60 s by less than 0.05: {', '.join(clock_fail) or 'none'}")
fig, ax = plt.subplots(figsize=(6, 3.2))
for _, r in clock.iterrows():
    ax.plot(shifts, r[shifts].to_numpy(), color="0.6", lw=0.7)
ax.plot(shifts, clock[shifts].median(), color="k", lw=2, label="median")
ax.set_xlabel("shift of the labels (s)"); ax.set_ylabel("balanced accuracy"); ax.axhline(0.5, ls="--", color="k", lw=0.8)
ax.legend(frameon=False); plt.tight_layout()
print(f"plotted accuracies range from {clock[shifts].min().min():.2f} to {clock[shifts].max().max():.2f}")

63 sessions. Shift of best accuracy: {-10: 1, -3: 1, 0: 26, 3: 35}
best accuracy within 3 s of zero: median 0.79 (range 0.59 to 0.90); best at +-60 s: median 0.56
sessions where accuracy near zero exceeds that at +-60 s by less than 0.05: sub-m928_ses-WT_7d_mating_d2, sub-m940_ses-WT_7d_mating_d7
plotted accuracies range from 0.35 to 0.90

def length(p):
    t0, t1 = annotations[p].interaction()
    return t1 - t0


def describe(p, q):
    a, b = annotations[p], annotations[q]
    trace = f"; trace best-match correlation {match(data[p].traces, data[q].traces):.2f}" if p in data and q in data else ""
    if p in data and q in data and data[p].traces.shape == data[q].traces.shape:
        trace += (", median correlation of cell i with cell i "
                  f"{np.median([np.corrcoef(data[p].traces[:, i], data[q].traces[:, i])[0, 1] for i in range(a.n_cells)]):.3f}")
    return (f"frames {a.n_frames} / {b.n_frames}, cells {a.n_cells} / {b.n_cells}, male onsets {len(male_onsets(a))} / "
            f"{len(male_onsets(b))}, interaction {length(p):.1f} / {length(q):.1f} s{trace}")


dup_of = {}
for _, row in sessions_t[sessions_t.experiment.str.contains("days")].iterrows():
    others = sessions_t[(sessions_t.subject.str.lstrip("m") == row.subject.lstrip("m"))
                        & sessions_t.experiment.isin(["receptive", "ovariectomised, primed"])]
    for _, o in others.iterrows():
        matched = receptive_day.get(row.subject) == row.day and o.experiment == "receptive" and row.experiment.startswith("natural")
        if matched or o.n_frames == row.n_frames:
            dup_of[row.path] = o.path
print(f"{len(dup_of)} days of 3-day or 7-day series duplicate a separately deposited session:")
for p, q in dup_of.items():
    print(f"  {p.split('_ecephys')[0]} = {q.split('_ecephys')[0]}: {describe(p, q)}")
p438 = "sub-m438_ses-wt_proestrus_receptive_mating_ecephys+behavior.nwb"
print("  for comparison, m438's other natural-cycle days against its receptive session:")
for p in m438_days[1:]:
    print(f"    {p.split('_ecephys')[0]}: {describe(p, p438)}")
print("Closest 7-day session to the receptive session, by frame count, for the mice that have both:")
for subject in ("m928", "m930", "m940"):
    rec = sessions_t[(sessions_t.subject == subject) & (sessions_t.experiment == "receptive")].path.iloc[0]
    days = sessions_t[(sessions_t.subject == subject) & sessions_t.experiment.str.startswith("oestrogen")]
    for p in days.iloc[(days.n_frames - annotations[rec].n_frames).abs().argsort()[:2]].path:
        print(f"  {p.split('_ecephys')[0]} against the receptive session: {describe(p, rec)}")


def onset_alignment(p, q, tol=0.2):
    """Largest share of p's male onsets that fall within tol s of one of q's under one constant offset."""
    a, b = male_onsets(annotations[p]), male_onsets(annotations[q])
    offsets = (b[None, :] - a[:, None]).ravel()
    offsets = np.unique(np.round(offsets[np.abs(offsets) < 600], 1))
    best = max(offsets, key=lambda o: np.sum(np.min(np.abs(a[:, None] + o - b[None, :]), 1) < tol))
    return np.mean(np.min(np.abs(a[:, None] + best - b[None, :]), 1) < tol), best


def shifted_match(p, q, max_shift=100, step=5):
    """Median over cells of p of the best correlation with any cell of q, at the frame shift that maximises it."""
    x, y = data[p].traces, data[q].traces
    best = (-1.0, 0)
    for k in range(-max_shift, max_shift + 1, step):
        u, v = (x[k:], y) if k >= 0 else (x, y[-k:])
        n = min(len(u), len(v))
        zu, zv = stats.zscore(u[:n], axis=0), stats.zscore(v[:n], axis=0)
        best = max(best, (float(np.median(np.max(zu.T @ zv / n, axis=1))), k))
    return best


print("Alignment checks, not set out in advance. For each pair: the share of male onsets matched within 0.2 s under\n"
      "one constant offset, and the trace best-match correlation at the best shift of up to 10 s in 0.5-s steps.")
p930 = sessions_t[(sessions_t.subject == "m930") & (sessions_t.experiment == "receptive")].path.iloc[0]
pairs = [(p, dup_of[p]) for p in dup_of if p in data and dup_of[p] in data]
pairs += [(p, p438) for p in m438_days[1:]]
pairs += [(p, p930) for p in sessions_t[(sessions_t.subject == "m930") & sessions_t.experiment.str.startswith("oestrogen")]
          .sort_values("day").path]
for p, q in pairs:
    share, offset = onset_alignment(p, q)
    r, k = shifted_match(p, q) if p in data and q in data else (np.nan, 0)
    print(f"  {p.split('_ecephys')[0]} against {q.split('_ecephys')[0]}: onsets matched {share:.2f} (offset {offset:.1f} s); "
          f"best-match correlation {r:.2f} at shift {k / vmhvl.RATE:+.1f} s")

distinct = sessions_t[~sessions_t.path.isin(dup_of)]
trial_sets = {
    "all annotated": distinct[distinct.female_annotated].path.tolist(),
    "receptive": distinct[distinct.female_annotated & distinct.receptive].path.tolist(),
    "all distinct, male annotation": distinct[distinct.experiment != "male inhibition"].path.tolist(),
}
for name, ps in trial_sets.items():
    subj = sessions_t.set_index("path").loc[ps, "subject"].str.lstrip("m")
    print(f"{name}: {len(ps)} trials from {subj.nunique()} mice")


def clipped(iv, t0, t1):
    return np.clip(iv, t0, t1)


def duration(s, labels, channel):
    t0, t1 = s.interaction()
    iv = clipped(vmhvl.intervals(s.behaviour, labels, channel), t0, t1)
    return float((iv[:, 1] - iv[:, 0]).sum())

9 days of 3-day or 7-day series duplicate a separately deposited session:
  sub-543_ses-WT_3d_mating_d1 = sub-m543_ses-wt_proestrus_receptive_mating: frames 10527 / 13682, cells 72 / 68, male onsets 63 / 63, interaction 915.3 / 915.2 s; trace best-match correlation 0.84
  sub-m438_ses-OVX_3d_mating_d3 = sub-m438_ses-ovx_primed_receptive_mating: frames 4216 / 4216, cells 64 / 64, male onsets 45 / 45, interaction 290.7 / 274.6 s; trace best-match correlation 1.00, median correlation of cell i with cell i 1.000
  sub-m438_ses-WT_3d_mating_d1 = sub-m438_ses-wt_proestrus_receptive_mating: frames 9484 / 9484, cells 66 / 64, male onsets 60 / 76, interaction 808.4 / 808.5 s; trace best-match correlation 0.78
  sub-m541_ses-OVX_3d_mating_d3 = sub-m541_ses-ovx_primed_receptive_mating: frames 6287 / 6287, cells 184 / 184, male onsets 54 / 57, interaction 498.7 / 498.1 s; trace best-match correlation 1.00, median correlation of cell i with cell i 1.000
  sub-m541_ses-WT_3d_mating_d3 = sub-m541_ses-wt_proestrus_receptive_mating: frames 6451 / 10455, cells 176 / 155, male onsets 31 / 31, interaction 507.6 / 503.3 s; trace best-match correlation 0.81
  sub-m542_ses-WT_3d_mating_d2 = sub-m542_ses-wt_proestrus_receptive_mating: frames 6613 / 9970, cells 129 / 124, male onsets 82 / 82, interaction 534.4 / 533.1 s
  sub-m551_ses-OVX_3d_mating_d3 = sub-m551_ses-ovx_primed_receptive_mating: frames 10191 / 10191, cells 141 / 141, male onsets 113 / 113, interaction 890.9 / 889.9 s; trace best-match correlation 1.00, median correlation of cell i with cell i 1.000
  sub-m551_ses-WT_3d_mating_d1 = sub-m551_ses-wt_proestrus_receptive_mating: frames 11510 / 15769, cells 171 / 122, male onsets 100 / 100, interaction 1020.5 / 1020.4 s
  sub-m764_ses-WT_3d_mating_d1 = sub-m764_ses-ovx_primed_receptive_mating: frames 9077 / 9077, cells 138 / 138, male onsets 72 / 72, interaction 763.7 / 763.7 s; trace best-match correlation 1.00, median correlation of cell i with cell i 1.000
  for comparison, m438's other natural-cycle days against its receptive session:
    sub-m438_ses-WT_3d_mating_d2: frames 7361 / 9484, cells 66 / 64, male onsets 79 / 76, interaction 606.7 / 808.5 s; trace best-match correlation 0.44
    sub-m438_ses-WT_3d_mating_d3: frames 9528 / 9484, cells 66 / 64, male onsets 96 / 76, interaction 821.4 / 808.5 s; trace best-match correlation 0.33
Closest 7-day session to the receptive session, by frame count, for the mice that have both:
  sub-m928_ses-WT_7d_mating_d5 against the receptive session: frames 5111 / 5398, cells 179 / 92, male onsets 70 / 74, interaction 295.9 / 413.1 s; trace best-match correlation 0.47
  sub-m928_ses-WT_7d_mating_d3 against the receptive session: frames 5037 / 5398, cells 179 / 92, male onsets 52 / 74, interaction 364.0 / 413.1 s; trace best-match correlation 0.42
  sub-m930_ses-WT_7d_mating_d3 against the receptive session: frames 4520 / 4511, cells 83 / 91, male onsets 55 / 65, interaction 322.7 / 322.3 s; trace best-match correlation 0.61
  sub-m930_ses-WT_7d_mating_d4 against the receptive session: frames 4924 / 4511, cells 83 / 91, male onsets 74 / 65, interaction 364.6 / 322.3 s; trace best-match correlation 0.73
  sub-m940_ses-WT_7d_mating_d5 against the receptive session: frames 6133 / 6160, cells 249 / 207, male onsets 82 / 44, interaction 486.2 / 487.9 s; trace best-match correlation 0.48

  sub-m940_ses-WT_7d_mating_d6 against the receptive session: frames 5479 / 6160, cells 249 / 207, male onsets 85 / 44, interaction 415.5 / 487.9 s; trace best-match correlation 0.50
Alignment checks, not set out in advance. For each pair: the share of male onsets matched within 0.2 s under
one constant offset, and the trace best-match correlation at the best shift of up to 10 s in 0.5-s steps.

  sub-543_ses-WT_3d_mating_d1 against sub-m543_ses-wt_proestrus_receptive_mating: onsets matched 1.00 (offset 0.0 s); best-match correlation 0.84 at shift +0.0 s
  sub-m438_ses-OVX_3d_mating_d3 against sub-m438_ses-ovx_primed_receptive_mating: onsets matched 1.00 (offset 0.1 s); best-match correlation 1.00 at shift +0.0 s

  sub-m438_ses-WT_3d_mating_d1 against sub-m438_ses-wt_proestrus_receptive_mating: onsets matched 0.65 (offset -0.1 s); best-match correlation 0.78 at shift +0.0 s

  sub-m541_ses-OVX_3d_mating_d3 against sub-m541_ses-ovx_primed_receptive_mating: onsets matched 1.00 (offset 0.1 s); best-match correlation 1.00 at shift +0.0 s

  sub-m541_ses-WT_3d_mating_d3 against sub-m541_ses-wt_proestrus_receptive_mating: onsets matched 1.00 (offset 0.2 s); best-match correlation 0.81 at shift +0.0 s

  sub-m551_ses-OVX_3d_mating_d3 against sub-m551_ses-ovx_primed_receptive_mating: onsets matched 1.00 (offset 0.1 s); best-match correlation 1.00 at shift +0.0 s

  sub-m764_ses-WT_3d_mating_d1 against sub-m764_ses-ovx_primed_receptive_mating: onsets matched 1.00 (offset -0.1 s); best-match correlation 1.00 at shift +0.0 s

  sub-m438_ses-WT_3d_mating_d2 against sub-m438_ses-wt_proestrus_receptive_mating: onsets matched 0.13 (offset -18.3 s); best-match correlation 0.49 at shift -9.0 s

  sub-m438_ses-WT_3d_mating_d3 against sub-m438_ses-wt_proestrus_receptive_mating: onsets matched 0.10 (offset -3.8 s); best-match correlation 0.38 at shift -8.5 s

BIN, SPAN = 20, 600
n_bins = SPAN // BIN


def group_curves(paths, labels):
    out = np.full((len(paths), n_bins), np.nan)
    for i, p in enumerate(paths):
        s = annotations[p]
        t0, t1 = s.interaction()
        frames = int(np.ceil(t1 * vmhvl.RATE)) + 1
        r = vmhvl.raster(s.behaviour, labels, frames, "Ch2")
        for k in range(n_bins):
            a, b = t0 + k * BIN, t0 + (k + 1) * BIN
            if b <= t1:
                out[i, k] = r[int(round(a * vmhvl.RATE)):int(round(b * vmhvl.RATE))].mean()
    return out


def judge_escalation(acc, res):
    counts = np.sum(~np.isnan(acc), 0)
    keep = counts >= 5
    ma, mr = np.nanmean(acc, 0)[keep], np.nanmean(res, 0)[keep]
    k = len(ma)
    rho, p = stats.spearmanr(np.arange(k), ma)
    third = k // 3
    first, last = ma[:third].mean(), ma[-third:].mean()
    run = np.convolve(mr, np.ones(3) / 3, "same")
    run[0], run[-1] = mr[:2].mean(), mr[-2:].mean()
    peak = int(np.argmax(run))
    accept_ok = rho > 0 and p < 0.05 and last > first
    resist_ok = 0 < peak < 2 * k / 3 and mr[-third:].mean() < run[peak]
    return dict(bins_kept=k, accept_spearman=round(rho, 2), accept_p=p, accept_first_third=round(first, 3),
                accept_last_third=round(last, 3), resist_peak_bin=peak, resist_peak=round(run[peak], 3),
                resist_first_bin=round(mr[0], 3), resist_last_third=round(mr[-third:].mean(), 3),
                accepting_escalates=accept_ok, resistance_rises_then_falls=resist_ok), ma, mr, keep


fig, axes = plt.subplots(1, 2, figsize=(10, 3.4), sharey=True)
rows = {}
for ax, name in zip(axes, ["all annotated", "receptive"]):
    ps = trial_sets[name]
    acc, res, app = (group_curves(ps, g) for g in (vmhvl.ACCEPT, vmhvl.RESIST, vmhvl.APPETITIVE))
    verdict, ma, mr, keep = judge_escalation(acc, res)
    rows[name] = verdict
    t = (np.arange(n_bins)[keep] + 0.5) * BIN
    for curve, label, color in ((acc, "accept", "C3"), (res, "resistance", "C0"), (app, "appetitive", "C2")):
        m = np.nanmean(curve, 0)[keep]
        se = np.nanstd(curve, 0)[keep] / np.sqrt(np.sum(~np.isnan(curve), 0)[keep])
        ax.plot(t, m, color=color, label=label); ax.fill_between(t, m - se, m + se, color=color, alpha=0.2)
    ax.set_title(f"{name} ({len(ps)} trials)", fontsize=10); ax.set_xlabel("time from male entry (s)")
    print(f"{name}: trials in each kept bin {np.sum(~np.isnan(acc), 0)[keep].tolist()}")
    print(f"{name}: trials per bin from {np.sum(~np.isnan(acc), 0)[keep].min()} to {np.sum(~np.isnan(acc), 0)[keep].max()}; "
          f"plotted means range {min(np.nanmean(c, 0)[keep].min() for c in (acc, res, app)):.3f} to "
          f"{max(np.nanmean(c, 0)[keep].max() for c in (acc, res, app)):.3f}")
axes[0].set_ylabel("probability"); axes[0].legend(frameon=False); plt.tight_layout()
pd.DataFrame(rows).T

all annotated: trials in each kept bin [69, 69, 69, 69, 69, 69, 68, 68, 68, 68, 68, 67, 67, 66, 63, 57, 54, 53, 45, 33, 28, 26, 23, 19, 15, 14, 12, 11, 11, 11]

all annotated: trials per bin from 11 to 69; plotted means range 0.005 to 0.255
receptive: trials in each kept bin [58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 58, 57, 57, 56, 54, 52, 51, 50, 43, 31, 27, 25, 22, 19, 15, 14, 12, 11, 11, 11]
receptive: trials per bin from 11 to 58; plotted means range 0.005 to 0.255

ps = trial_sets["receptive"]
acc = group_curves(ps, vmhvl.ACCEPT)
subj = sessions_t.set_index("path").loc[ps, "subject"].str.lstrip("m").to_numpy()
by_mouse = []
for m in np.unique(subj):
    curve = np.nanmean(acc[subj == m], 0)
    ok = ~np.isnan(curve)
    by_mouse.append(dict(mouse=m, trials=int((subj == m).sum()), bins=int(ok.sum()),
                         spearman=round(stats.spearmanr(np.arange(ok.sum()), curve[ok])[0], 2)))
by_mouse = pd.DataFrame(by_mouse)
print(f"{(by_mouse.spearman > 0).sum()} of {len(by_mouse)} mice have a positive correlation")
print(by_mouse.set_index("mouse").T.to_string())
for name in ("all annotated", "receptive"):
    acc = group_curves(trial_sets[name], vmhvl.ACCEPT)
    rhos = []
    for row in acc:
        ok = ~np.isnan(row)
        if ok.sum() >= 5 and np.ptp(row[ok]) > 0:
            rhos.append(stats.spearmanr(np.arange(ok.sum()), row[ok])[0])
    rhos = np.array(rhos)
    print(f"{name}, within trials ({len(rhos)} trials with at least 5 bins and some change in accepting behaviour; "
          f"{len(acc) - len(rhos)} left out): positive in {(rhos > 0).sum()}, median Spearman {np.median(rhos):.2f}, "
          f"sign test p = {stats.binomtest(int((rhos > 0).sum()), len(rhos)).pvalue:.2g}")
    full = ~np.isnan(acc[:, :15]).any(1)
    m = acc[full, :15].mean(0)
    rho, p = stats.spearmanr(np.arange(15), m)
    print(f"{name}, first 15 bins, the {full.sum()} trials that cover them all: Spearman {rho:.2f} (p = {p:.2g}); "
          f"first third {m[:5].mean():.3f}, last third {m[-5:].mean():.3f}")

12 of 23 mice have a positive correlation
mouse       438    449    540    541    542    543    546    551   553    764    774    775    776    857    858    859    864    918    928    929   930    931    940
trials     2.00   1.00   1.00   2.00   1.00   1.00   1.00   2.00   1.0   3.00   1.00   1.00   1.00   1.00   2.00   2.00   1.00   5.00   5.00   7.00   7.0   5.00   5.00
bins      30.00  23.00  24.00  25.00  26.00  30.00  30.00  30.00  27.0  30.00  15.00  22.00  30.00  30.00  30.00  30.00  23.00  20.00  20.00  22.00  22.0  21.00  24.00
spearman  -0.31   0.39  -0.21  -0.27   0.18   0.12   0.02   0.52   0.0   0.13   0.23   0.12  -0.38  -0.28   0.17   0.19  -0.13   0.25  -0.08  -0.29   0.1  -0.18  -0.26
all annotated, within trials (58 trials with at least 5 bins and some change in accepting behaviour; 11 left out): positive in 33, median Spearman 0.01, sign test p = 0.36
all annotated, first 15 bins, the 63 trials that cover them all: Spearman 0.63 (p = 0.012); first third 0.080, last third 0.151
receptive, within trials (58 trials with at least 5 bins and some change in accepting behaviour; 0 left out): positive in 33, median Spearman 0.01, sign test p = 0.36
receptive, first 15 bins, the 54 trials that cover them all: Spearman 0.63 (p = 0.012); first third 0.093, last third 0.176

def behaviour_times(paths):
    rows = []
    for p in paths:
        s = annotations[p]
        t0, t1 = s.interaction()
        rows.append(dict(path=p, subject=sessions_t.set_index("path").subject[p].lstrip("m"),
                         interaction=t1 - t0,
                         responsive=duration(s, vmhvl.RESPONSIVE, "Ch2"),
                         female_self=duration(s, vmhvl.FEMALE_SELF, "Ch2"),
                         male_self=duration(s, vmhvl.MALE_SELF, "Ch1"),
                         accept=duration(s, vmhvl.ACCEPT, "Ch2")))
    return pd.DataFrame(rows)

times = {name: behaviour_times(trial_sets[name]) for name in ("all annotated", "receptive")}
rows = []
for name, t in times.items():
    frac = t.responsive / (t.responsive + t.female_self)
    rows.append({"trials": name, "n": len(t), "mean of per-trial fractions": round(frac.mean(), 3),
                 "SEM": round(frac.sem(), 3), "pooled fraction": round(t.responsive.sum() / (t.responsive + t.female_self).sum(), 3),
                 "range": f"{frac.min():.2f} to {frac.max():.2f}"})
print("expected 0.839, tolerance 0.084 (0.755 to 0.923)")
pd.DataFrame(rows).set_index("trials")

expected 0.839, tolerance 0.084 (0.755 to 0.923)

t = times["receptive"]
t["fraction"] = t.responsive / (t.responsive + t.female_self)
print("Receptive trials, by mouse:")
t.groupby("subject").fraction.agg(["size", "mean", "min", "max"]).round(3).T

Receptive trials, by mouse:

rows = []
for name, t in times.items():
    ok = t.female_self > 0
    ratio = t.male_self[ok] / t.female_self[ok]
    rows.append({"trials": name, "n": len(t), "trials without female self-initiated time": int((~ok).sum()),
                 "mean of per-trial ratios": round(ratio.mean(), 2), "SEM": round(ratio.sem(), 2),
                 "ratio of summed times": round(t.male_self.sum() / t.female_self.sum(), 2),
                 "median per-trial ratio": round(ratio.median(), 2),
                 "range": f"{ratio.min():.1f} to {ratio.max():.1f}"})
print("expected 11.3, tolerance 1.13 (10.17 to 12.43)")
pd.DataFrame(rows).set_index("trials")

expected 11.3, tolerance 1.13 (10.17 to 12.43)

t = times["receptive"]
t = t[t.female_self > 0].assign(ratio=lambda d: d.male_self / d.female_self)
print("Receptive trials, by mouse (per-trial ratio):")
t.groupby("subject").ratio.agg(["size", "mean", "min", "max"]).round(1).T

Receptive trials, by mouse (per-trial ratio):

def copulation(p, gap):
    s = annotations[p]
    return vmhvl.bouts(vmhvl.intervals(s.behaviour, vmhvl.COPULATION, "Ch1"), gap)

rows = []
for name, ps in trial_sets.items():
    for gap_name, gap in (("0.1 s", 0.1), ("1 s", 1.0), ("no merging", None)):
        b = [copulation(p, gap) for p in ps]
        ibi = np.concatenate([x[1:, 0] - x[:-1, 1] for x in b if len(x) > 1])
        dur = np.concatenate([x[:, 1] - x[:, 0] for x in b])
        rows.append({"trials": name, "merge gaps up to": gap_name, "trials n": len(ps),
                     "bouts": len(dur), "intervals": len(ibi), "mean interval (s)": round(ibi.mean(), 2),
                     "median interval (s)": round(np.median(ibi), 2), "mean bout (s)": round(dur.mean(), 2),
                     "intervals below 0": int((ibi < 0).sum())})
print("expected mean interval 13.7 s, tolerance 1.37 s (12.33 to 15.07)")
pd.DataFrame(rows).set_index(["trials", "merge gaps up to"])

expected mean interval 13.7 s, tolerance 1.37 s (12.33 to 15.07)

by_mouse = []
for p in trial_sets["receptive"]:
    x = copulation(p, 0.1)
    by_mouse.append(dict(subject=sessions_t.set_index("path").subject[p].lstrip("m"),
                         intervals=len(x) - 1 if len(x) else 0,
                         total=(x[1:, 0] - x[:-1, 1]).sum() if len(x) > 1 else 0.0))
by_mouse = pd.DataFrame(by_mouse).groupby("subject").sum()
by_mouse["mean interval (s)"] = (by_mouse.total / by_mouse.intervals).round(1)
print("Receptive trials, 0.1-s merging, by mouse:")
by_mouse[["intervals", "mean interval (s)"]].T

Receptive trials, 0.1-s merging, by mouse:

def social_fractions(ps, gap=0.1):
    in_bout, in_ibi, tb, sb, ti, si = [], [], 0, 0, 0, 0
    for p in ps:
        s = annotations[p]
        frames = int(np.ceil(s.behaviour.stop.max() * vmhvl.RATE)) + 2
        social = vmhvl.raster(s.behaviour, vmhvl.SOCIAL, frames, "Ch2")
        b = copulation(p, gap)
        for (a, z) in b:
            seg = social[int(round(a * vmhvl.RATE)):max(int(round(z * vmhvl.RATE)), int(round(a * vmhvl.RATE)) + 1)]
            in_bout.append(seg.mean()); tb += len(seg); sb += seg.sum()
        for (a, z) in zip(b[:-1, 1], b[1:, 0]):
            seg = social[int(round(a * vmhvl.RATE)):int(round(z * vmhvl.RATE))]
            if len(seg):
                in_ibi.append(seg.mean()); ti += len(seg); si += seg.sum()
    return np.array(in_bout), np.array(in_ibi), sb / tb, si / ti

rows = []
for name in ("all annotated", "receptive"):
    fb, fi, pb, pi = social_fractions(trial_sets[name])
    rows.append({"trials": name, "bouts": len(fb), "intervals": len(fi),
                 "bout, mean of fractions": round(fb.mean(), 3), "bout, pooled": round(pb, 3),
                 "interval, mean of fractions": round(fi.mean(), 3), "interval, pooled": round(pi, 3)})
print("expected 0.62 (0.558 to 0.682) in bouts and 0.23 (0.207 to 0.253) in intervals")
pd.DataFrame(rows).set_index("trials")

expected 0.62 (0.558 to 0.682) in bouts and 0.23 (0.207 to 0.253) in intervals

m930_days = sessions_t[(sessions_t.subject == "m930") & sessions_t.experiment.str.startswith("oestrogen")]
m930_close = m930_days[m930_days.day.astype(int).isin([3, 4])].sort_values("day").path.tolist()
rows = []
for name, drop in (("all annotated", []), ("without day 3", m930_close[:1]), ("without day 4", m930_close[1:]),
                   ("without days 3 and 4", m930_close)):
    ps = [p for p in trial_sets["all annotated"] if p not in drop]
    t = behaviour_times(ps)
    ok = t.female_self > 0
    b = [copulation(p, 0.1) for p in ps]
    ibi = np.concatenate([x[1:, 0] - x[:-1, 1] for x in b if len(x) > 1])
    fb, fi, _, _ = social_fractions(ps)
    rows.append({"trials": name, "n": len(ps),
                 "1c responsive fraction (0.755 to 0.923)": round((t.responsive / (t.responsive + t.female_self)).mean(), 3),
                 "1d mean ratio (10.17 to 12.43)": round((t.male_self[ok] / t.female_self[ok]).mean(), 2),
                 "1e mean interval, s (12.33 to 15.07)": round(ibi.mean(), 2),
                 "1g bout (0.558 to 0.682)": round(fb.mean(), 3), "1g interval (0.207 to 0.253)": round(fi.mean(), 3)})
print("m930 days left out:", [p.split("_ecephys")[0] for p in m930_close])
pd.DataFrame(rows).set_index("trials")

m930 days left out: ['sub-m930_ses-WT_7d_mating_d3', 'sub-m930_ses-WT_7d_mating_d4']

rng_cp = np.random.default_rng(1097_02)
cp_behaviours = [(b, "Ch2") for b in vmhvl.FEMALE_MATING] + [(b, "Ch1") for b in vmhvl.MALE_SELF]
cp_rows = []
for p in receptive_paths:
    s = data[p]
    t0, t1 = s.interaction()
    i0, i1 = int(round(t0 * vmhvl.RATE)), min(int(round(t1 * vmhvl.RATE)), s.n_frames)
    y = vmhvl.bin_frames(vmhvl.zscore(s.traces)[i0:i1], 10)
    for b, ch in cp_behaviours:
        lab = vmhvl.bin_frames(vmhvl.raster(s.behaviour, [b], s.n_frames, ch)[i0:i1].astype(float), 10) >= 0.5
        if lab.sum() < 10 or (~lab).sum() < 10:
            continue
        cp, mu, sd = vmhvl.choice_probability(y, lab, rng_cp)
        tuned = ((cp > 0.7) | (cp < 0.3)) & (np.abs(cp - mu) > 2 * sd)
        cp_rows.append(dict(subject=s.subject, behaviour=b, bins=int(lab.sum()), cells=s.n_cells,
                            tuned=int(tuned.sum()), tuned_up=int((tuned & (cp > 0.5)).sum()),
                            beyond_2sd_only=int((np.abs(cp - mu) > 2 * sd).sum())))
cp_rows = pd.DataFrame(cp_rows)
mice = cp_rows.groupby("behaviour").subject.nunique()
analysed = mice[mice >= 5].index
cp_table = cp_rows[cp_rows.behaviour.isin(analysed)].groupby("behaviour").agg(
    mice=("subject", "nunique"), cells=("cells", "sum"), tuned=("tuned", "sum"), tuned_up=("tuned_up", "sum"),
    beyond_2sd=("beyond_2sd_only", "sum"))
cp_table["tuned %"] = (100 * cp_table.tuned / cp_table.cells).round(1)
cp_table["beyond 2 SD of shuffles, any CP %"] = (100 * cp_table.beyond_2sd / cp_table.cells).round(1)
print(f"behaviours with at least 10 bins in at least 5 mice: {len(analysed)}; left out: "
      f"{', '.join(f'{b} ({n} mice)' for b, n in mice[mice < 5].items()) or 'none'}")
print(f"median tuned % over behaviours: {cp_table['tuned %'].median():.1f}; range {cp_table['tuned %'].min():.1f} to {cp_table['tuned %'].max():.1f}")
cp_table.sort_values("tuned %")

behaviours with at least 10 bins in at least 5 mice: 11; left out: turn (3 mice)
median tuned % over behaviours: 22.0; range 6.9 to 37.4

per_mouse = cp_rows[cp_rows.behaviour.isin(analysed)].assign(pct=lambda d: 100 * d.tuned / d.cells)
print("Tuned % by mouse (rows) and behaviour (columns); blank where the mouse had fewer than 10 bins:")
per_mouse.pivot(index="subject", columns="behaviour", values="pct").round(0)

Tuned % by mouse (rows) and behaviour (columns); blank where the mouse had fewer than 10 bins:

rng_cp_pairs = np.random.default_rng(1097_06)
pair_rows = []
for p in receptive_paths:
    s = data[p]
    t0, t1 = s.interaction()
    i0, i1 = int(round(t0 * vmhvl.RATE)), min(int(round(t1 * vmhvl.RATE)), s.n_frames)
    y = vmhvl.bin_frames(vmhvl.zscore(s.traces)[i0:i1], 10)
    for group, channel, names in (("female", "Ch2", vmhvl.FEMALE_MATING), ("male", "Ch1", vmhvl.MALE_SELF)):
        labs = {b: vmhvl.bin_frames(vmhvl.raster(s.behaviour, [b], s.n_frames, channel)[i0:i1].astype(float), 10) >= 0.5
                for b in names}
        for i, a in enumerate(names):
            for b in names[i + 1:]:
                la, lb = labs[a] & ~labs[b], labs[b] & ~labs[a]
                if la.sum() < 10 or lb.sum() < 10:
                    continue
                keep = la | lb
                cp, mu, sd = vmhvl.choice_probability(y[keep], la[keep], rng_cp_pairs)
                tuned = ((cp > 0.7) | (cp < 0.3)) & (np.abs(cp - mu) > 2 * sd)
                pair_rows.append(dict(subject=s.subject, group=group, pair=f"{a} vs {b}", cells=s.n_cells, tuned=int(tuned.sum())))
pair_rows = pd.DataFrame(pair_rows)
pairs = pair_rows.groupby(["group", "pair"]).agg(mice=("subject", "nunique"), cells=("cells", "sum"), tuned=("tuned", "sum"))
print(f"pairs with at least 10 bins of each behaviour in some mouse: {len(pairs)}; in at least 5 mice: {(pairs.mice >= 5).sum()}")
pairs = pairs[pairs.mice >= 5]
pairs["tuned %"] = (100 * pairs.tuned / pairs.cells).round(1)
print(f"tuned % over the kept pairs: median {pairs['tuned %'].median():.1f}, range {pairs['tuned %'].min():.1f} to "
      f"{pairs['tuned %'].max():.1f}; pairs within 1 to 15%: {pairs['tuned %'].between(1, 15).sum()} of {len(pairs)}")
pairs.sort_values("tuned %")

pairs with at least 10 bins of each behaviour in some mouse: 39; in at least 5 mice: 20
tuned % over the kept pairs: median 39.3, range 15.4 to 57.6; pairs within 1 to 15%: 0 of 20

lags = np.arange(0, 101, 10)
rng_glm = np.random.default_rng(1097_03)


def glm_designs(s):
    t0, t1 = s.interaction()
    i0, i1 = int(round(t0 * vmhvl.RATE)), min(int(round(t1 * vmhvl.RATE)), s.n_frames)
    Y = vmhvl.zscore(s.traces)[i0:i1]
    fem = np.column_stack([vmhvl.raster(s.behaviour, [b], s.n_frames, "Ch2")[i0:i1] for b in vmhvl.FEMALE_MATING])
    mal = np.column_stack([vmhvl.raster(s.behaviour, [b], s.n_frames, "Ch1")[i0:i1] for b in vmhvl.MALE_SELF])
    Xf, Xm = vmhvl.lagged(fem.astype(float), lags), vmhvl.lagged(mal.astype(float), lags)
    Xb = np.hstack([Xf, Xm])
    return Y, {"female": (Xf, None), "male": (Xm, None), "both": (Xb, None),
               "coupling": (np.hstack([Xb, Y]), Xb.shape[1] + np.arange(Y.shape[1]))}


glm, glm_random = [], []
for p in receptive_paths:
    s = data[p]
    Y, designs = glm_designs(s)
    for rng, out in ((None, glm), (rng_glm, glm_random)):
        r2 = {}
        for name, (X, exclude) in designs.items():
            r2[name], pred, fold_of = vmhvl.cv_ridge(X, Y, exclude=exclude, rng=rng, return_pred=True)
            r2[f"{name}, within folds"] = vmhvl.within_fold_r2(Y, pred, fold_of)
        pred = vmhvl.intercept_predictions(Y, fold_of)
        r2["intercept only"] = vmhvl.pooled_r2(Y, pred)
        r2["intercept only, within folds"] = vmhvl.within_fold_r2(Y, pred, fold_of)
        for i in range(Y.shape[1]):
            out.append(dict(subject=s.subject, cell=i, **{k: v[i] for k, v in r2.items()}))
glm, glm_random = pd.DataFrame(glm), pd.DataFrame(glm_random)
models = ["female", "male", "both", "coupling"]
expected = {"female": 0.08, "male": 0.14, "both": 0.15, "coupling": 0.46}
summary = pd.DataFrame({
    "expected": expected,
    "mean over cells": glm[models].mean(),
    "mean of mouse means": glm.groupby("subject")[models].mean().mean(),
    "median over cells": glm[models].median(),
    "cells below 0": (glm[models] < 0).sum(),
    "variant, random folds, mean over cells": glm_random[models].mean(),
    "variant, random folds, mean of mouse means": glm_random.groupby("subject")[models].mean().mean(),
}).round(3)
print(f"{len(glm)} cells from {glm.subject.nunique()} mice; cvR2 range over cells: "
      + "; ".join(f"{m} {glm[m].min():.2f} to {glm[m].max():.2f}" for m in models))
summary

1335 cells from 15 mice; cvR2 range over cells: female -1.78 to 0.50; male -0.58 to 0.60; both -5.15 to 0.63; coupling -1.27 to 0.97

rows = {}
for folds, g in (("contiguous", glm), ("random", glm_random)):
    rows[(folds, "intercept only")] = {"pooled R2, mean over cells": g["intercept only"].mean(),
                                       "within-fold R2, mean over cells": g["intercept only, within folds"].mean(),
                                       "within-fold R2, median over cells": g["intercept only, within folds"].median(),
                                       "cells above the baseline": np.nan}
    for m in models:
        rows[(folds, m)] = {"pooled R2, mean over cells": g[m].mean(),
                            "within-fold R2, mean over cells": g[f"{m}, within folds"].mean(),
                            "within-fold R2, median over cells": g[f"{m}, within folds"].median(),
                            "gain over intercept only, mean over cells": (g[m] - g["intercept only"]).mean(),
                            "gain, mean of mouse means": (g[m] - g["intercept only"]).groupby(g.subject).mean().mean(),
                            "cells above the baseline": int((g[m] > g["intercept only"]).sum())}
print(f"{len(glm)} cells; expected from the paper: female 0.08, male 0.14, both 0.15, coupling 0.46")
pd.DataFrame(rows).T.round(3)

1335 cells; expected from the paper: female 0.08, male 0.14, both 0.15, coupling 0.46

glm_inner = []
for p in receptive_paths:
    s = data[p]
    Y, designs = glm_designs(s)
    r2 = {name: vmhvl.cv_ridge(X, Y, exclude=exclude, inner_folds=5) for name, (X, exclude) in designs.items()}
    for i in range(Y.shape[1]):
        glm_inner.append(dict(subject=s.subject, cell=i, **{k: v[i] for k, v in r2.items()}))
glm_inner = pd.DataFrame(glm_inner)
assert (glm_inner[["subject", "cell"]].to_numpy() == glm[["subject", "cell"]].to_numpy()).all()
gain = glm_inner[models].sub(glm["intercept only"], axis=0)
print(f"inner CV, both minus male: mean over cells {(glm_inner.both - glm_inner.male).mean():+.4f}, "
      f"mean of mouse means {(glm_inner.groupby('subject').both.mean() - glm_inner.groupby('subject').male.mean()).mean():+.4f}")
pd.DataFrame({
    "expected": expected,
    "last fifth, mean over cells": glm[models].mean(),
    "inner CV, mean over cells": glm_inner[models].mean(),
    "inner CV, mean of mouse means": glm_inner.groupby("subject")[models].mean().mean(),
    "inner CV, gain over intercept only": gain.mean(),
    "inner CV, gain, mean of mouse means": gain.groupby(glm_inner.subject).mean().mean(),
    "inner CV, cells below 0": (glm_inner[models] < 0).sum(),
}).round(3).assign(**{"cells where both is at least male, last fifth / inner CV":
                      [np.nan, np.nan, f"{(glm.both >= glm.male).sum()} / {(glm_inner.both >= glm_inner.male).sum()}", np.nan]})

inner CV, both minus male: mean over cells +0.0003, mean of mouse means +0.0014

print("Mean cvR2 by mouse:")
glm.groupby("subject")[models].mean().round(3).T

Mean cvR2 by mouse:

fig, ax = plt.subplots(figsize=(6, 3.2))
for i, m in enumerate(models):
    ax.violinplot(np.clip(glm[m], -0.2, 1), positions=[i], showmedians=True)
ax.set_xticks(range(4), models); ax.set_ylabel("cvR2 (clipped at -0.2)"); plt.tight_layout()
print(f"values below -0.2 clipped: " + ", ".join(f"{m} {(glm[m] < -0.2).sum()}" for m in models))

values below -0.2 clipped: female 166, male 99, both 160, coupling 17

achw_rows = []
for p in receptive_paths:
    s = data[p]
    t0, t1 = s.interaction()
    i0, i1 = int(round(t0 * vmhvl.RATE)), min(int(round(t1 * vmhvl.RATE)), s.n_frames)
    windows = (("whole recording", (0, s.n_frames)), ("interaction", (i0, i1)),
               ("variant, start of imaging to end of annotation", (0, i1)))
    for window, (a, b) in windows:
        for smoothed in (True, False):
            y = vmhvl.prepare(s.traces[a:b], smoothed=smoothed)
            for first in (False, True):
                h = vmhvl.achw(y, first_crossing=first)
                achw_rows.append(dict(subject=s.subject, window=window, smoothed=smoothed,
                                      rule="first crossing" if first else "closest to 0.5 (code)",
                                      cells=len(h), above_25=(h > 25).sum(), sum=h.sum(), at_cap=(h >= 50).sum()))
achw_rows = pd.DataFrame(achw_rows)
achw_table = achw_rows.groupby(["window", "smoothed", "rule"])[["cells", "above_25", "sum", "at_cap"]].sum()
achw_table["fraction above 25 s"] = (achw_table.above_25 / achw_table.cells).round(3)
achw_table["mean ACHW (s)"] = (achw_table["sum"] / achw_table.cells).round(1)
print("expected fraction 0.30 (0.27 to 0.33) and mean 20 s (18 to 22); ACHW is capped at 50 s (500 lags)")
achw_table[["cells", "fraction above 25 s", "mean ACHW (s)", "at_cap"]]

expected fraction 0.30 (0.27 to 0.33) and mean 20 s (18 to 22); ACHW is capped at 50 s (500 lags)

main = achw_rows[(achw_rows.window == "whole recording") & achw_rows.smoothed & (achw_rows.rule == "closest to 0.5 (code)")]
print("Whole recording, smoothed, code's rule, by mouse:")
main.assign(fraction=lambda d: (d.above_25 / d.cells).round(2), mean=lambda d: (d["sum"] / d.cells).round(1)
            ).set_index("subject")[["cells", "fraction", "mean"]].T

Whole recording, smoothed, code's rule, by mouse:

cycle_mice = {"m541": ("m541", "sub-m541_ses-WT_3d_mating_d{}_ecephys+behavior.nwb"),
              "543": ("m542", "sub-543_ses-WT_3d_mating_d{}_ecephys+behavior.nwb")}
proj = {}
for mouse, (wfile, pattern) in cycle_mice.items():
    rday = int(receptive_day[mouse])
    days = {d: data[pattern.format(d)] for d in (1, 2, 3)}
    W = weights[wfile]["weights"]
    rec = days[rday]
    t_entry = int(round(rec.interaction()[0] * vmhvl.RATE))
    capped = [vmhvl.achw((rec.traces @ W[:, k])[:, None])[0] for k in range(W.shape[1])]
    # The 500-lag cap of the authors' ACHW can tie; the factor is chosen with lags up to 300 s.
    halfwidths = [vmhvl.achw((rec.traces @ W[:, k])[:, None], max_lag=3000)[0] for k in range(W.shape[1])]
    chosen = int(np.argmax(halfwidths))
    print(f"mouse {mouse} (weights from sub-{wfile}): ACHW of the receptive-day projection, capped at 50 s, "
          f"factor 1 {capped[0]:.1f} s and factor 2 {capped[1]:.1f} s; up to 300 s, {halfwidths[0]:.1f} s and "
          f"{halfwidths[1]:.1f} s; integration factor {chosen + 1}")
    for k in (chosen, 1 - chosen):
        x_rec = rec.traces @ W[:, k]
        sign = 1 if x_rec[t_entry:].mean() > x_rec[:t_entry].mean() else -1
        w = sign * W[:, k]
        xs = {d: s.traces @ w for d, s in days.items()}
        lo, hi = min(x.min() for x in xs.values()), max(x.max() for x in xs.values())
        rlo, rhi = xs[rday].min(), xs[rday].max()
        rows = []
        for d, s in days.items():
            mount = vmhvl.raster(s.behaviour, ["mount"], s.n_frames, "Ch1")
            t0, t1 = s.interaction()
            e = int(round(t0 * vmhvl.RATE))
            end = min(int(round(t1 * vmhvl.RATE)), s.n_frames)
            x = xs[d]
            rows.append({"day": d, "receptive": d == rday, "mount frames": int(mount.sum()),
                         "intromission bouts": int((s.behaviour.behavior == "intromission").sum()),
                         "mounting, min-max over 3 days": round(((x[mount] - lo) / (hi - lo)).mean(), 3),
                         "mounting, receptive-day range": round(((x[mount] - rlo) / (rhi - rlo)).mean(), 3),
                         "ramp: last 100 s minus 60 s before entry (SD of receptive day)":
                             round((x[end - 1000:end].mean() - x[e - 600:e].mean()) / xs[rday].std(), 2)})
        table = pd.DataFrame(rows).set_index("day")
        if k == chosen:
            proj[mouse] = dict(factor=k + 1, sign=sign, xs=xs, rday=rday, lo=lo, hi=hi, table=table)
            print(f"  factor {k + 1} (the integration factor), sign {sign:+d}:")
        else:
            print(f"  variant, factor {k + 1}, sign {sign:+d}:")
        print(table.to_string())
    print()

mouse m541 (weights from sub-m541): ACHW of the receptive-day projection, capped at 50 s, factor 1 50.0 s and factor 2 50.0 s; up to 300 s, 56.8 s and 104.8 s; integration factor 2
  factor 2 (the integration factor), sign +1:
     receptive  mount frames  intromission bouts  mounting, min-max over 3 days  mounting, receptive-day range  ramp: last 100 s minus 60 s before entry (SD of receptive day)
day                                                                                                                                                                           
1        False          1460                   0                          0.231                          0.187                                                            0.03
2        False           251                   0                          0.191                          0.144                                                            0.51
3         True           176                   7                          0.614                          0.592                                                            2.54
  variant, factor 1, sign +1:
     receptive  mount frames  intromission bouts  mounting, min-max over 3 days  mounting, receptive-day range  ramp: last 100 s minus 60 s before entry (SD of receptive day)
day                                                                                                                                                                           
1        False          1460                   0                          0.483                          0.483                                                            0.63
2        False           251                   0                          0.488                          0.488                                                           -0.08
3         True           176                   7                          0.619                          0.619                                                            1.34

mouse 543 (weights from sub-m542): ACHW of the receptive-day projection, capped at 50 s, factor 1 50.0 s and factor 2 32.1 s; up to 300 s, 128.9 s and 32.1 s; integration factor 1
  factor 1 (the integration factor), sign +1:
     receptive  mount frames  intromission bouts  mounting, min-max over 3 days  mounting, receptive-day range  ramp: last 100 s minus 60 s before entry (SD of receptive day)
day                                                                                                                                                                           
1         True           337                   4                          0.606                          0.606                                                            2.05
2        False           876                   0                          0.643                          0.643                                                            0.18
3        False           722                   1                          0.626                          0.626                                                           -0.35
  variant, factor 2, sign -1:
     receptive  mount frames  intromission bouts  mounting, min-max over 3 days  mounting, receptive-day range  ramp: last 100 s minus 60 s before entry (SD of receptive day)
day                                                                                                                                                                           
1         True           337                   4                          0.655                          0.655                                                            2.05
2        False           876                   0                          0.538                          0.538                                                            1.07
3        False           722                   1                          0.505                          0.505                                                            1.03

fig, axes = plt.subplots(1, 2, figsize=(10, 3.2), sharey=True)
for ax, (mouse, pr) in zip(axes, proj.items()):
    for d, x in pr["xs"].items():
        s = data[cycle_mice[mouse][1].format(d)]
        e = int(round(s.interaction()[0] * vmhvl.RATE))
        t = (np.arange(len(x)) - e) / vmhvl.RATE
        ax.plot(t, (x - pr["lo"]) / (pr["hi"] - pr["lo"]), lw=0.8,
                label=f"day {d}" + (" (receptive)" if d == pr["rday"] else ""))
    ax.set_title(f"mouse {mouse}", fontsize=10); ax.set_xlabel("time from male entry (s)"); ax.legend(frameon=False, fontsize=8)
axes[0].set_ylabel("integration factor, min-max over 3 days"); plt.tight_layout()
print("all plotted values lie in 0 to 1 by construction")

all plotted values lie in 0 to 1 by construction

ed7 = []
for mouse, (wfile, pattern) in cycle_mice.items():
    pr = proj[mouse]
    w = weights[wfile]["weights"][:, pr["factor"] - 1]
    weighted = np.abs(w) > np.std(w, ddof=1)
    h = {d: vmhvl.achw(vmhvl.prepare(data[pattern.format(d)].traces)) for d in (1, 2, 3)}
    unrec = [d for d in (1, 2, 3) if d != pr["rday"]]
    for group, mask in (("weighted", weighted), ("other", ~weighted)):
        hr = h[pr["rday"]][mask]
        comparisons = {"unreceptive pooled": np.concatenate([h[d][mask] for d in unrec])}
        comparisons.update({f"unreceptive day {d}": h[d][mask] for d in unrec})
        for name, hu in comparisons.items():
            ks = stats.ks_2samp(hr, hu)
            ed7.append({"mouse": mouse, "cells": group, "n": int(mask.sum()), "against": name,
                        "receptive mean (s)": round(hr.mean(), 1), "receptive SEM": round(stats.sem(hr), 1),
                        "unreceptive mean (s)": round(hu.mean(), 1), "unreceptive SEM": round(stats.sem(hu), 1),
                        "KS statistic": round(ks.statistic, 2), "KS p": ks.pvalue})
ed7 = pd.DataFrame(ed7).set_index(["mouse", "cells", "against"])
print("expected: weighted cells 25.2 s on receptive against 16.1 s on unreceptive days; other cells unchanged")
ed7

expected: weighted cells 25.2 s on receptive against 16.1 s on unreceptive days; other cells unchanged

print("Variants: the one-SD threshold of the authors' code replaced by the top 10% and the top 25% of |weight|,\n"
      "and the one-SD threshold on the other factor; unreceptive days pooled.")
rows = []
for mouse, (wfile, pattern) in cycle_mice.items():
    pr = proj[mouse]
    W = np.abs(weights[wfile]["weights"])
    h = {d: vmhvl.achw(vmhvl.prepare(data[pattern.format(d)].traces)) for d in (1, 2, 3)}
    unrec = [d for d in (1, 2, 3) if d != pr["rday"]]
    k, other = pr["factor"] - 1, 2 - pr["factor"]
    sd = np.std(weights[wfile]["weights"], axis=0, ddof=1)
    on_k, on_other = W[:, k] > sd[k], W[:, other] > sd[other]
    print(f"mouse {mouse}: {on_k.sum()} cells above 1 SD on factor {k + 1}, {on_other.sum()} on factor {other + 1}, "
          f"{(on_k & on_other).sum()} on both")
    selections = {f"factor {k + 1}, top 10%": W[:, k] > np.quantile(W[:, k], 0.90),
                  f"factor {k + 1}, top 25%": W[:, k] > np.quantile(W[:, k], 0.75),
                  f"factor {other + 1}, above 1 SD": W[:, other] > np.std(weights[wfile]["weights"][:, other], ddof=1)}
    for name, mask in selections.items():
        for group, m in (("weighted", mask), ("other", ~mask)):
            hr, hu = h[pr["rday"]][m], np.concatenate([h[d][m] for d in unrec])
            rows.append({"mouse": mouse, "selection": name, "cells": group, "n": int(m.sum()),
                         "receptive mean": round(hr.mean(), 1), "unreceptive mean": round(hu.mean(), 1),
                         "KS p": stats.ks_2samp(hr, hu).pvalue})
pd.DataFrame(rows).set_index(["mouse", "selection", "cells"])

Variants: the one-SD threshold of the authors' code replaced by the top 10% and the top 25% of |weight|,
and the one-SD threshold on the other factor; unreceptive days pooled.
mouse m541: 44 cells above 1 SD on factor 2, 44 on factor 1, 23 on both

mouse 543: 14 cells above 1 SD on factor 1, 20 on factor 2, 7 on both

ed10 = []
for p in oestrogen_paths:
    s = data[p]
    t0, t1 = s.interaction()
    e, end = int(round(t0 * vmhvl.RATE)), min(int(round(t1 * vmhvl.RATE)), s.n_frames)
    m = s.traces.mean(1)
    base = m[max(0, e - 600):e].mean()
    contact = s.behaviour[(s.behaviour.channel == "Ch1") & s.behaviour.behavior.isin(vmhvl.MALE_SELF)].start.min()
    c = int(round(contact * vmhvl.RATE))
    ed10.append({"session": p.split("_ecephys")[0], "mouse": s.subject,
                 "accepting %": 100 * duration(s, vmhvl.ACCEPT, "Ch2") / (t1 - t0),
                 "AUC entry to end, baseline-subtracted": (m[e:end] - base).sum() / vmhvl.RATE,
                 "AUC entry to end, raw": m[e:end].sum() / vmhvl.RATE,
                 "AUC 200 s from first contact, baseline-subtracted": (m[c:c + 2000] - base).sum() / vmhvl.RATE})
ed10 = pd.DataFrame(ed10)
rows = []
for col in ed10.columns[3:]:
    r = stats.linregress(ed10["accepting %"], ed10[col])
    rows.append({"reading": col, "trials": len(ed10), "r": round(r.rvalue, 3), "r2": round(r.rvalue ** 2, 3),
                 "p": r.pvalue, "slope": r.slope})
print(f"{len(ed10)} trials from {ed10.mouse.nunique()} mice; accepting % from {ed10['accepting %'].min():.1f} to "
      f"{ed10['accepting %'].max():.1f}; {int((ed10['accepting %'] == 0).sum())} trials without accepting behaviour")
keep = ~ed10.session.isin(clock_fail)
for col in ed10.columns[3:]:
    r = stats.linregress(ed10.loc[keep, "accepting %"], ed10.loc[keep, col])
    rows.append({"reading": f"{col}, without {', '.join(clock_fail) or 'none'}", "trials": int(keep.sum()),
                 "r": round(r.rvalue, 3), "r2": round(r.rvalue ** 2, 3), "p": r.pvalue, "slope": r.slope})
pd.DataFrame(rows).set_index("reading")

42 trials from 6 mice; accepting % from 0.0 to 29.2; 11 trials without accepting behaviour

fig, ax = plt.subplots(figsize=(5, 3.4))
col = "AUC entry to end, baseline-subtracted"
for i, (mouse, g) in enumerate(ed10.groupby("mouse")):
    ax.scatter(g["accepting %"], g[col], s=18, label=mouse, color=f"C{i}")
ax.set_xlabel("accepting behaviour (% of interaction)"); ax.set_ylabel("population mean AUC (z x s)")
ax.legend(frameon=False, fontsize=8); plt.tight_layout()
print(f"plotted AUC range {ed10[col].min():.0f} to {ed10[col].max():.0f}")
print("By mouse, Spearman correlation of accepting % with this AUC:")
print(ed10.groupby("mouse").apply(lambda g: round(stats.spearmanr(g["accepting %"], g[col])[0], 2)).to_dict())

plotted AUC range -13 to 156
By mouse, Spearman correlation of accepting % with this AUC:
{'m918': -0.67, 'm928': -0.19, 'm929': 0.25, 'm930': 0.04, 'm931': 0.52, 'm940': -0.52}

print("Male-inhibition sessions (Figure 4c,d) and their inhibition epochs:")
for p in sessions_t[sessions_t.experiment == "male inhibition"].path:
    b = annotations[p].behaviour
    inh = b[b.behavior.isin(["laser_on", "mating_inhibition"])]
    print(f"  {p.split('_ecephys')[0]}: " + "; ".join(f"{r.behavior} {r.start:.0f} to {r.stop:.0f} s ({r.stop - r.start:.0f} s)" for r in inh.itertuples()))
print()
ej = b_all[b_all.behavior == "ejaculation"].path.unique()
print(f"Sessions with an annotated ejaculation (Figure 4g,h): {len(ej)}: {', '.join(p.split('_ecephys')[0] for p in ej)}")
print()
print("Labels that could mark a pencil-cup period (Extended Data Figure 3a-c):",
      [l for l in b_all.behavior.unique() if "cup" in l.lower() or "pencil" in l.lower()] or "none")
print("Labels that could mark a 10-s optogenetic inhibition of VMHvl (Figure 3h-k):",
      [l for l in b_all.behavior.unique() if "led" in l.lower() or "nvoke" in l.lower()] or "none")

Male-inhibition sessions (Figure 4c,d) and their inhibition epochs:
  sub-m764_ses-wireless_opto: laser_on 383 to 696 s (313 s)
  sub-m774_ses-wireless_opto: laser_on 815 to 1122 s (307 s)
  sub-m775_ses-wireless_opto: laser_on 488 to 794 s (305 s)
  sub-m776_ses-wireless_opto: laser_on 368 to 669 s (301 s)
  sub-m858_ses-wireless_opto: mating_inhibition 230 to 581 s (351 s)
  sub-m918_ses-wireless_opto: mating_inhibition 512 to 739 s (227 s); mating_inhibition 755 to 933 s (178 s)
  sub-m931_ses-wireless_opto: mating_inhibition 376 to 722 s (346 s); mating_inhibition 874 to 980 s (106 s)
  sub-m932_ses-wireless_opto: mating_inhibition 243 to 422 s (179 s)

Sessions with an annotated ejaculation (Figure 4g,h): 12: sub-m764_ses-WT_3d_mating_d1, sub-m764_ses-WT_3d_mating_d3, sub-m764_ses-ovx_primed_receptive_mating, sub-m764_ses-wireless_opto, sub-m775_ses-ovx_primed_receptive_mating, sub-m858_ses-wireless_opto, sub-m859_ses-ovx_primed_receptive_mating, sub-m928_ses-WT_7d_mating_d2, sub-m928_ses-WT_7d_mating_d5, sub-m931_ses-WT_7d_mating_d4, sub-m932_ses-wireless_opto, sub-m940_ses-wt_proestrus_receptive_mating

Labels that could mark a pencil-cup period (Extended Data Figure 3a-c): none
Labels that could mark a 10-s optogenetic inhibition of VMHvl (Figure 3h-k): none

tails = pd.Series({p: annotations[p].n_frames / vmhvl.RATE - annotations[p].behaviour.stop.max() for p in annotations})
long_tails = tails[tails > 60].sort_values(ascending=False)
print(f"{len(long_tails)} sessions with more than 60 s of imaging after the last annotation:")
for p, v in long_tails.items():
    print(f"  {p.split('_ecephys')[0]} ({sessions_t.set_index('path').experiment[p]}): {v:.0f} s, {annotations[p].n_cells} cells")
rows = []
for p in receptive_paths:
    s = data[p]
    t0, t1 = s.interaction()
    i0, i1 = int(round(t0 * vmhvl.RATE)), min(int(round(t1 * vmhvl.RATE)), s.n_frames)
    n_tail = s.n_frames - i1
    if n_tail < 1800:
        continue
    h_int = vmhvl.achw(vmhvl.prepare(s.traces[i0:i1]))
    h_end = vmhvl.achw(vmhvl.prepare(s.traces[max(i0, i1 - n_tail):i1]))
    h_tail = vmhvl.achw(vmhvl.prepare(s.traces[i1:]))
    rows.append({"mouse": s.subject, "tail (s)": n_tail / vmhvl.RATE, "cells": s.n_cells,
                 "interaction, mean ACHW": h_int.mean(), "same-length end of interaction": h_end.mean(),
                 "tail": h_tail.mean(), "tail above 25 s": (h_tail > 25).mean(), "interaction above 25 s": (h_int > 25).mean(),
                 "Wilcoxon p, tail against same-length end": stats.wilcoxon(h_end, h_tail).pvalue})
tail_table = pd.DataFrame(rows).set_index("mouse")
print(f"{len(tail_table)} receptive sessions, {tail_table.cells.sum()} cells")
tail_table.round(3)

17 sessions with more than 60 s of imaging after the last annotation:
  sub-m551_ses-wt_proestrus_receptive_mating (receptive): 426 s, 122 cells
  sub-m541_ses-wt_proestrus_receptive_mating (receptive): 405 s, 155 cells
  sub-m542_ses-wt_proestrus_receptive_mating (receptive): 337 s, 124 cells
  sub-m764_ses-wireless_opto (male inhibition): 317 s, 124 cells
  sub-m543_ses-wt_proestrus_receptive_mating (receptive): 316 s, 68 cells
  sub-m540_ses-wt_proestrus_receptive_mating (receptive): 307 s, 39 cells
  sub-m546_ses-wt_proestrus_receptive_mating (receptive): 302 s, 20 cells
  sub-m858_ses-wireless_opto (male inhibition): 273 s, 128 cells
  sub-m553_ses-wt_proestrus_receptive_mating (receptive): 272 s, 16 cells
  sub-m932_ses-wireless_opto (male inhibition): 246 s, 82 cells
  sub-m775_ses-ovx_primed_receptive_mating (ovariectomised, primed): 233 s, 56 cells
  sub-m928_ses-WT_7d_mating_d2 (oestrogen-injected (7 days)): 216 s, 179 cells
  sub-m764_ses-WT_3d_mating_d3 (natural cycle (3 days)): 171 s, 138 cells
  sub-m931_ses-WT_7d_mating_d4 (oestrogen-injected (7 days)): 131 s, 128 cells
  sub-m859_ses-ovx_primed_receptive_mating (ovariectomised, primed): 119 s, 120 cells
  sub-m928_ses-WT_7d_mating_d5 (oestrogen-injected (7 days)): 88 s, 179 cells
  sub-m931_ses-wireless_opto (male inhibition): 61 s, 82 cells

7 receptive sessions, 544 cells

for mouse, (wfile, pattern) in cycle_mice.items():
    pr = proj[mouse]
    s = data[pattern.format(pr["rday"])]
    t0, t1 = s.interaction()
    i0, i1 = int(round(t0 * vmhvl.RATE)), min(int(round(t1 * vmhvl.RATE)), s.n_frames)
    w = weights[wfile]["weights"][:, pr["factor"] - 1]
    weighted = np.abs(w) > np.std(w, ddof=1)
    y = vmhvl.prepare(s.traces)[:i1]
    ramp = np.zeros(i1)
    ramp[i0:] = np.linspace(0, 1, i1 - i0)
    r = np.array([np.corrcoef(y[:, i], ramp)[0, 1] for i in range(y.shape[1])])
    for group, m in (("weighted", weighted), ("other", ~weighted)):
        print(f"mouse {mouse}, {group} cells ({m.sum()}): r2 above 0.5 in {(r[m] ** 2 > 0.5).sum()} "
              f"({(r[m] ** 2 > 0.5).mean():.0%}); r from {r[m].min():.2f} to {r[m].max():.2f}, median {np.median(r[m]):.2f}")

mouse m541, weighted cells (44): r2 above 0.5 in 0 (0%); r from -0.60 to 0.71, median 0.19
mouse m541, other cells (132): r2 above 0.5 in 0 (0%); r from -0.52 to 0.54, median -0.07
mouse 543, weighted cells (14): r2 above 0.5 in 1 (7%); r from -0.76 to 0.59, median -0.34
mouse 543, other cells (58): r2 above 0.5 in 0 (0%); r from -0.57 to 0.69, median -0.05

Claim	Paper	Here
`fig1b-accepting-escalates`	Accepting behaviour rises over mating; resistance rises then falls	Tested
`fig1c-responsive-dominates`	83.9% of female mating time is male-responsive	Tested
`fig1d-male-drives-mating`	Males spend 11.3 times more time in self-initiated mating behaviour	Tested
`fig1e-interbout-interval`	Mean interbout interval 13.7 s	Tested
`fig1g-social-in-bouts-and-intervals`	Social behaviour for 62% of bouts and 23% of intervals	Tested
`fig2c-few-tuned-cells`	About 2 to 13% of cells tuned to a behaviour	Tested
`fig2e-glm-behaviour-explains-little`	Behaviour GLMs explain 8 to 15% of variance	Tested
`fig2f-glm-coupling`	Adding the other cells raises it to 46%	Tested
`fig2i-persistent-cells`	30% of cells have ACHW above 25 s; mean 20 s	Tested
`ed3b-pencil-cup`	ACHW shorter with the male behind a pencil cup	Blocked: the period is not marked
`fig3a-rslds-fit`	rSLDS explains 64% of variance	Blocked: needs rSLDS fits
`fig3c-integration-dimension`	One dimension with a 111-s time constant against 25 s	Blocked: needs rSLDS fits
`ed3s-integrator-model`	Male behaviour drives an integrator with a time constant above 200 s	Blocked: needs rSLDS fits
`ed4b-single-cell-ramps`	56% of integration cells ramp	Blocked: needs rSLDS fits
`fig3k-perturbation-relaxation`	Activity relaxes back to the attractor after inhibition	Blocked: the sessions are not in the dandiset
`fig4b-persistence-during-intervals`	Integration dimension as high in intervals as in bouts	Blocked: needs rSLDS fits
`fig4d-persistence-male-inhibition`	It stays high while the male is inhibited	Blocked: needs rSLDS fits
`fig4h-peak-before-ejaculation`	It peaks before ejaculation	Blocked: needs rSLDS fits
`fig5c-attractor-receptive-days`	Line attractor on receptive days only	Blocked: needs rSLDS fits
`ed6e-projection-unreceptive-days`	The receptive-day dimension stays low on unreceptive days	Tested with the authors' weights, 2 mice
`ed7-achw-across-cycle`	Weighted cells are less persistent on unreceptive days	Tested with the authors' weights, 2 mice
`ed9j-ovx-hormone-priming`	The attractor returns with hormone priming	Blocked: needs rSLDS fits
`fig5k-receptivity-correlation`	Integration dimension tracks receptivity (r2 0.62)	Blocked: needs rSLDS fits
`ed10c-population-mean-uncorrelated`	The population mean does not track receptivity	Tested

subject	438	449	540	541	542	543	546	551	553	764	774	775	776	857	858	859	864	918	928	929	930	931	940
size	2.000	1.000	1.00	2.000	1.000	1.000	1.000	2.000	1.000	3.000	1.000	1.000	1.000	1.000	2.000	2.000	1.000	5.000	5.000	7.000	7.000	5.000	5.000
mean	0.786	0.842	0.94	0.834	0.646	0.797	0.703	0.815	0.762	0.842	0.816	0.792	0.862	0.912	0.923	0.889	0.786	0.838	0.848	0.882	0.758	0.799	0.856
min	0.764	0.842	0.94	0.757	0.646	0.797	0.703	0.744	0.762	0.785	0.816	0.792	0.862	0.912	0.916	0.852	0.786	0.696	0.792	0.736	0.661	0.717	0.745
max	0.808	0.842	0.94	0.911	0.646	0.797	0.703	0.885	0.762	0.884	0.816	0.792	0.862	0.912	0.929	0.926	0.786	0.913	0.945	0.960	0.866	0.851	0.940

		trials n	bouts	intervals	mean interval (s)	median interval (s)	mean bout (s)	intervals below 0
trials	merge gaps up to
all annotated	0.1 s	69	1540	1471	14.00	10.07	5.22	0
	1 s	69	1394	1325	15.48	11.23	5.84	0
	no merging	69	2180	2111	9.76	4.57	3.69	3
receptive	0.1 s	58	1278	1220	14.86	11.05	5.81	0
	1 s	58	1178	1120	16.13	11.90	6.36	0
	no merging	58	1893	1835	9.88	4.57	3.92	3
all distinct, male annotation	0.1 s	88	2070	1982	14.61	10.40	4.60	0
	1 s	88	1860	1772	16.27	11.72	5.19	0
	no merging	88	2764	2676	10.82	5.37	3.45	3

subject	438	449	540	541	542	543	546	551	553	764	774	775	776	857	858	859	864	918	928	929	930	931	940
intervals	38.0	29.0	20.0	18.0	18.0	8.0	38.0	53.0	10.0	79.0	17.0	14.0	33.0	36.0	48.0	31.0	24.0	119.0	101.0	155.0	136.0	108.0	87.0
mean interval (s)	19.8	12.8	14.5	32.6	21.4	73.0	15.5	20.8	30.9	12.8	10.3	21.0	27.6	10.7	14.6	20.5	15.7	11.1	12.0	11.8	12.6	12.5	14.3

	bouts	intervals	bout, mean of fractions	bout, pooled	interval, mean of fractions	interval, pooled
trials
all annotated	1540	1471	0.633	0.741	0.218	0.114
receptive	1278	1220	0.593	0.748	0.174	0.101

Reproducing Liu, Nair et al. (2024) from dandiset 001097¶

What can and cannot be tested¶

Versions of the dandiset¶

Load the data and check it against its description¶

Duplicated recordings and the per-mouse files¶

Are the weights in the same cell order as the session files?¶

Which natural-cycle day is the receptive one¶

Trial sets for the behaviour claims¶

Figure 1b: accepting behaviour escalates over the interaction¶

Figure 1c: most female mating behaviour responds to the male¶

Figure 1d: the male drives the interaction¶

Figure 1e: copulation bouts and interbout intervals¶

Figure 2c: few cells are tuned to any one behaviour¶

Figure 2e and 2f: behaviour explains little of a cell's activity; other cells explain more¶

Figure 2i: many cells are persistent¶

Extended Data Figure 6e and Figure 5h: the receptive-day dimension on unreceptive days¶

Extended Data Figure 7: single-cell persistence across the cycle¶

Extended Data Figure 10c: the population mean does not track receptivity¶

Claims that cannot be tested here¶

Summary¶

How this differs from the paper's analysis¶

What this attempt found about the dandiset¶

	bins_kept	accept_spearman	accept_p	accept_first_third	accept_last_third	resist_peak_bin	resist_peak	resist_first_bin	resist_last_third	accepting_escalates	resistance_rises_then_falls
all annotated	30	0.74	0.000003	0.106	0.192	2	0.143	0.102	0.066	True	True
receptive	30	0.59	0.000626	0.126	0.194	2	0.114	0.076	0.065	True	True

	n	mean of per-trial fractions	SEM	pooled fraction	range
trials
all annotated	69	0.841	0.010	0.844	0.65 to 0.99
receptive	58	0.829	0.011	0.835	0.65 to 0.96

	n	trials without female self-initiated time	mean of per-trial ratios	SEM	ratio of summed times	median per-trial ratio	range
trials
all annotated	69	0	11.75	1.40	7.51	9.24	2.2 to 79.3
receptive	58	0	9.31	0.74	7.12	8.20	2.2 to 28.5

subject	438	449	540	541	542	543	546	551	553	764	774	775	776	857	858	859	864	918	928	929	930	931	940
size	2.0	1.0	1.0	2.0	1.0	1.0	1.0	2.0	1.0	3.0	1.0	1.0	1.0	1.0	2.0	2.0	1.0	5.0	5.0	7.0	7.0	5.0	5.0
mean	4.7	8.2	16.2	7.0	2.2	5.4	4.9	7.3	4.0	7.4	5.4	7.4	9.9	11.1	15.2	11.7	4.8	8.6	12.0	14.0	6.7	8.3	11.5
min	3.7	8.2	16.2	3.9	2.2	5.4	4.9	3.2	4.0	5.7	5.4	7.4	9.9	11.1	12.8	8.2	4.8	3.7	5.7	3.8	2.7	4.4	5.1
max	5.7	8.2	16.2	10.2	2.2	5.4	4.9	11.4	4.0	9.6	5.4	7.4	9.9	11.1	17.6	15.2	4.8	12.9	28.5	25.6	10.2	10.2	25.8

	mice	cells	tuned	tuned_up	beyond_2sd	tuned %	beyond 2 SD of shuffles, any CP %
behaviour
approach	6	394	27	22	84	6.9	21.3
mount	14	1128	109	64	480	9.7	42.6
stay	14	1128	138	80	479	12.2	42.5
lordose	15	1335	250	113	833	18.7	62.4
male_sniff	15	1335	285	210	871	21.3	65.2
top_up	8	623	137	74	243	22.0	39.0
intromission	15	1335	369	146	1001	27.6	75.0
check_genital	8	690	201	84	345	29.1	50.0
sniff	8	676	216	154	341	32.0	50.4
wiggle	11	925	319	121	604	34.5	65.3
dart	5	454	170	100	207	37.4	45.6

behaviour	approach	check_genital	dart	intromission	lordose	male_sniff	mount	sniff	stay	top_up	wiggle
subject
m438	NaN	41.0	44.0	50.0	42.0	31.0	12.0	30.0	30.0	9.0	48.0
m540	NaN	NaN	33.0	36.0	23.0	10.0	18.0	NaN	8.0	18.0	NaN
m541	NaN	25.0	NaN	26.0	21.0	37.0	11.0	NaN	5.0	NaN	NaN
m542	4.0	NaN	11.0	35.0	40.0	0.0	31.0	6.0	41.0	27.0	NaN
m543	7.0	10.0	NaN	31.0	15.0	31.0	15.0	NaN	7.0	NaN	66.0
m546	0.0	NaN	20.0	50.0	35.0	20.0	0.0	15.0	5.0	NaN	50.0
m551	NaN	22.0	NaN	11.0	6.0	34.0	2.0	52.0	2.0	9.0	37.0
m553	12.0	62.0	NaN	62.0	50.0	38.0	38.0	NaN	50.0	6.0	62.0
m857	14.0	NaN	NaN	15.0	20.0	0.0	0.0	1.0	24.0	22.0	11.0
m858	NaN	21.0	NaN	10.0	5.0	2.0	2.0	NaN	2.0	NaN	17.0
m859	5.0	NaN	NaN	32.0	17.0	28.0	1.0	43.0	5.0	30.0	38.0
m864	NaN	NaN	NaN	29.0	9.0	4.0	9.0	43.0	1.0	NaN	35.0
m928	NaN	17.0	NaN	25.0	14.0	17.0	8.0	NaN	11.0	39.0	NaN
m930	NaN	65.0	NaN	26.0	0.0	37.0	2.0	52.0	7.0	NaN	34.0
m940	NaN	NaN	54.0	30.0	21.0	25.0	NaN	NaN	NaN	NaN	29.0

		mice	cells	tuned	tuned %
group	pair
female	lordose vs stay	14	1128	174	15.4
male	male_sniff vs mount	14	1128	191	16.9
female	stay vs top_up	8	623	127	20.4
male	mount vs intromission	14	1128	236	20.9
female	lordose vs wiggle	11	925	236	25.5
	lordose vs top_up	8	623	160	25.7
	approach vs stay	6	394	129	32.7
male	male_sniff vs intromission	15	1335	473	35.4
female	lordose vs dart	5	454	171	37.7
	sniff vs top_up	5	476	184	38.7
	stay vs check_genital	8	690	275	39.9
	sniff vs stay	8	676	272	40.2
	lordose vs check_genital	8	690	293	42.5
	wiggle vs stay	10	718	313	43.6
	approach vs lordose	6	394	173	43.9
	sniff vs lordose	8	676	305	45.1
	wiggle vs top_up	5	368	168	45.7
	sniff vs wiggle	7	552	272	49.3
	approach vs wiggle	5	270	147	54.4
	wiggle vs check_genital	6	443	255	57.6

	expected	mean over cells	mean of mouse means	median over cells	cells below 0	variant, random folds, mean over cells	variant, random folds, mean of mouse means
female	0.08	-0.077	-0.050	-0.071	1027	0.229	0.238
male	0.14	-0.013	0.008	-0.035	829	0.207	0.220
both	0.15	-0.057	-0.017	-0.048	889	0.322	0.334
coupling	0.46	0.450	0.434	0.480	73	0.814	0.778

		pooled R2, mean over cells	within-fold R2, mean over cells	within-fold R2, median over cells	cells above the baseline	gain over intercept only, mean over cells	gain, mean of mouse means
contiguous	intercept only	-0.086	-1.400	-0.698	NaN	NaN	NaN
	female	-0.077	-1.418	-0.711	703.0	0.009	0.031
	male	-0.013	-1.208	-0.610	856.0	0.073	0.089
	both	-0.057	-1.409	-0.643	801.0	0.029	0.064
	coupling	0.450	-0.008	0.105	1284.0	0.536	0.515
random	intercept only	-0.000	-0.002	-0.002	NaN	NaN	NaN
	female	0.229	0.227	0.209	1335.0	0.229	0.238
	male	0.207	0.205	0.179	1335.0	0.208	0.220
	both	0.322	0.320	0.297	1335.0	0.322	0.334
	coupling	0.814	0.813	0.848	1335.0	0.815	0.778

subject	m438	m540	m541	m542	m543	m546	m551	m553	m857	m858	m859	m864	m928	m930	m940
female	0.023	-0.059	-0.078	0.013	-0.078	0.047	-0.088	0.094	-0.030	-0.100	-0.039	-0.076	-0.085	-0.104	-0.192
male	0.098	0.052	-0.044	0.039	-0.020	0.076	-0.033	0.188	-0.025	-0.067	0.013	0.004	-0.035	-0.104	-0.019
both	0.075	0.010	-0.063	0.036	-0.061	0.122	-0.020	0.181	-0.040	-0.076	0.026	-0.013	-0.054	-0.137	-0.237
coupling	0.437	0.396	0.444	0.443	0.313	0.334	0.446	0.411	0.420	0.571	0.492	0.409	0.374	0.524	0.500

subject	m438	m540	m541	m542	m543	m546	m551	m553	m857	m858	m859	m864	m928	m930	m940
cells	64.00	39.00	155.00	124.00	68.00	20.0	122.00	16.0	79.00	82.00	87.00	89.00	92.00	91.00	207.0
fraction	0.02	0.28	0.22	0.16	0.13	0.1	0.33	0.0	0.08	0.17	0.17	0.03	0.07	0.48	0.3
mean	6.80	20.20	17.70	14.10	13.00	8.8	20.60	9.7	11.40	15.50	14.50	9.80	9.00	24.90	18.3

			n	receptive mean (s)	receptive SEM	unreceptive mean (s)	unreceptive SEM	KS statistic	KS p
mouse	cells	against
m541	weighted	unreceptive pooled	44	28.4	2.1	17.2	1.3	0.43	0.000024
		unreceptive day 1	44	28.4	2.1	15.6	1.7	0.48	0.000066
		unreceptive day 2	44	28.4	2.1	18.7	1.9	0.39	0.002557
	other	unreceptive pooled	132	14.9	1.1	12.1	0.6	0.11	0.233400
		unreceptive day 1	132	14.9	1.1	12.1	1.0	0.14	0.171996
		unreceptive day 2	132	14.9	1.1	12.1	0.9	0.11	0.449216
543	weighted	unreceptive pooled	14	26.7	3.4	9.1	1.5	0.68	0.000179
		unreceptive day 2	14	26.7	3.4	10.2	2.3	0.64	0.004900
		unreceptive day 3	14	26.7	3.4	8.0	1.9	0.79	0.000163
	other	unreceptive pooled	58	12.8	1.3	10.8	1.0	0.28	0.005017
		unreceptive day 2	58	12.8	1.3	10.0	1.2	0.29	0.013315
		unreceptive day 3	58	12.8	1.3	11.6	1.5	0.26	0.040917

	tail (s)	cells	interaction, mean ACHW	same-length end of interaction	tail	tail above 25 s	interaction above 25 s	Wilcoxon p, tail against same-length end
mouse
m540	307.0	39	9.267	4.195	5.785	0.026	0.077	0.106
m541	404.9	155	9.602	7.679	6.652	0.013	0.045	0.474
m542	337.1	124	9.070	5.595	6.703	0.016	0.048	0.158
m543	315.9	68	11.353	4.644	5.338	0.000	0.103	0.247
m546	302.2	20	7.310	3.155	3.960	0.000	0.100	0.236
m551	426.3	122	17.877	8.288	11.492	0.057	0.270	0.000
m553	271.9	16	8.056	5.625	4.862	0.000	0.000	0.274

			cells	fraction above 25 s	mean ACHW (s)	at_cap
window	smoothed	rule
interaction	False	closest to 0.5 (code)	1335	0.100	9.9	19
	False	first crossing	1335	0.100	9.9	23
	True	closest to 0.5 (code)	1335	0.124	12.0	23
	True	first crossing	1335	0.121	11.9	29
variant, start of imaging to end of annotation	False	closest to 0.5 (code)	1335	0.158	12.7	23
	False	first crossing	1335	0.155	12.7	32
	True	closest to 0.5 (code)	1335	0.197	15.0	36
	True	first crossing	1335	0.195	15.0	43
whole recording	False	closest to 0.5 (code)	1335	0.166	13.1	31
	False	first crossing	1335	0.163	13.1	45
	True	closest to 0.5 (code)	1335	0.200	15.5	51
	True	first crossing	1335	0.199	15.5	60

			n	receptive mean	unreceptive mean	KS p
mouse	selection	cells
m541	factor 2, top 10%	weighted	18	32.9	17.5	0.001796
	factor 2, top 10%	other	158	16.6	12.9	0.026897
	factor 2, top 25%	weighted	44	28.4	17.2	0.000024
	factor 2, top 25%	other	132	14.9	12.1	0.233400
	factor 1, above 1 SD	weighted	44	30.8	17.5	0.000007
	factor 1, above 1 SD	other	132	14.1	12.0	0.073436
543	factor 1, top 10%	weighted	8	27.0	11.8	0.125566
	factor 1, top 10%	other	64	14.1	10.3	0.002081
	factor 1, top 25%	weighted	18	24.7	8.5	0.000053
	factor 1, top 25%	other	54	12.5	11.2	0.028808
	factor 2, above 1 SD	weighted	20	21.0	12.9	0.069840
	factor 2, above 1 SD	other	52	13.5	9.6	0.001606

	trials	r	r2	p	slope
reading
AUC entry to end, baseline-subtracted	42	0.121	0.015	0.445170	0.872317
AUC entry to end, raw	42	0.054	0.003	0.732813	0.236698
AUC 200 s from first contact, baseline-subtracted	42	0.164	0.027	0.298175	0.804148
AUC entry to end, baseline-subtracted, without sub-m928_ses-WT_7d_mating_d2, sub-m940_ses-WT_7d_mating_d7	40	0.104	0.011	0.524530	0.737282
AUC entry to end, raw, without sub-m928_ses-WT_7d_mating_d2, sub-m940_ses-WT_7d_mating_d7	40	0.073	0.005	0.654207	0.319013
AUC 200 s from first contact, baseline-subtracted, without sub-m928_ses-WT_7d_mating_d2, sub-m940_ses-WT_7d_mating_d7	40	0.157	0.025	0.331796	0.763578

Claim	Outcome	Note
`fig1b-accepting-escalates`	Reproduced	Accepting rises (Spearman 0.74); resistance peaks early
`fig1c-responsive-dominates`	Reproduced	84.1% against 83.9%
`fig1d-male-drives-mating`	Reproduced	11.75 against 11.3, on one of six readings
`fig1e-interbout-interval`	Reproduced	14.0 s against 13.7 s
`fig1g-social-in-bouts-and-intervals`	Reproduced	63.3% and 21.8% against 62% and 23%, on one of four readings
`fig2c-few-tuned-cells`	Not reproduced	7 to 37% tuned, median 22%; more with pairs of behaviours
`fig2e-glm-behaviour-explains-little`	Not reproduced	The ordering fails with the verified design (-0.077, -0.013, -0.057); a cross-validated penalty would make it partially reproduced, proposed and not verified; a weak test
`fig2f-glm-coupling`	Reproduced	0.477 (0.450 with the last-fifth penalty) against 0.46, with the fold design of Figure 2e
`fig2i-persistent-cells`	Not reproduced	20% above 25 s, mean 15.5 s
`ed3b-pencil-cup`	Blocked: metadata insufficient	The pencil-cup period is not marked; an exploratory look at the unannotated tail is shown
`fig3a-rslds-fit`, `fig3c-integration-dimension`, `ed3s-integrator-model`, `ed4b-single-cell-ramps`, `fig4b-persistence-during-intervals`, `fig4d-persistence-male-inhibition`, `fig4h-peak-before-ejaculation`, `fig5c-attractor-receptive-days`, `ed9j-ovx-hormone-priming`, `fig5k-receptivity-correlation`	Blocked: other	The rSLDS cannot be built or fit within the registry's environment and time limit
`fig3k-perturbation-relaxation`	Blocked: data missing	The VMHvl inhibition sessions are not in the dandiset
`ed6e-projection-unreceptive-days`	Partially reproduced	Holds in m541, not in 543; on the other factor it would hold in both
`ed7-achw-across-cycle`	Partially reproduced	Weighted cells hold in both mice; other cells change in 543; m541's other factor shows the same change
`ed10c-population-mean-uncorrelated`	Reproduced	r2 at most 0.03

Reproducing Liu, Nair et al. (2024) from dandiset 001097¶

What can and cannot be tested¶

Versions of the dandiset¶

Load the data and check it against its description¶

Duplicated recordings and the per-mouse files¶

Are the weights in the same cell order as the session files?¶

Which natural-cycle day is the receptive one¶

Do the traces and the annotations share a clock?¶

Trial sets for the behaviour claims¶

Figure 1b: accepting behaviour escalates over the interaction¶

Figure 1c: most female mating behaviour responds to the male¶

Figure 1d: the male drives the interaction¶

Figure 1e: copulation bouts and interbout intervals¶

Figure 1g: females stay social during bouts and after them¶

Figure 2c: few cells are tuned to any one behaviour¶

Figure 2e and 2f: behaviour explains little of a cell's activity; other cells explain more¶

Figure 2i: many cells are persistent¶

Extended Data Figure 6e and Figure 5h: the receptive-day dimension on unreceptive days¶

Extended Data Figure 7: single-cell persistence across the cycle¶

Extended Data Figure 10c: the population mean does not track receptivity¶

Claims that cannot be tested here¶

Summary¶

How this differs from the paper's analysis¶

What this attempt found about the dandiset¶