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

import sleepstates as ss

pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
pd.set_option("display.max_rows", 200)

# Streams every usable file: units, states, and the LFP power of the reference channel in
# Non-REM epochs.
sessions, excluded = ss.load_all(workers=6)
edges, sext = ss.wake_sextiles(sessions)
print(len(sessions), "sessions loaded; excluded:", excluded)

21 sessions loaded; excluded: {'sub-Dino/sub-Dino_ses-Dino-072114_ecephys.nwb': 'spike times are whole seconds divided by 20000'}

import collections
import h5py, remfile

urls = dict(ss.asset_urls())
rows = []
for s in sessions:
    dur = s.stops - s.starts
    spk_all = np.concatenate([x for x in s.spikes if len(x)])
    rows.append(dict(
        subject=s.subject, session=s.session_id,
        pE=int((s.cell_type == "excitatory").sum()), pI=int((s.cell_type == "inhibitory").sum()),
        unknown=int((s.cell_type == "unknown").sum()),
        states_min=round((s.stops[-1] - s.starts[0]) / 60, 1),
        spikes_first=round(spk_all.min(), 2), spikes_last_min=round(spk_all.max() / 60, 1),
        lfp_min=round(s.lfp_samples / ss.LFP_RATE / 60, 1),
        WAKE_min=round(dur[s.state == "WAKE"].sum() / 60, 1), MA_n=int((s.state == "MA").sum()),
        NREM_min=round(dur[s.state == "NREM"].sum() / 60, 1), REM_min=round(dur[s.state == "REM"].sum() / 60, 1),
        epochs=len(s.labels), lfp_packets=len(s.lfp)))
tab = pd.DataFrame(rows)
print(tab.to_string())
print("\nTotals over the 21 usable sessions:", tab[["pE", "pI", "unknown", "MA_n"]].sum().to_dict())
print("Subjects:", tab.subject.nunique(), "| sessions per subject:", tab.groupby("subject").size().to_dict())
print("pE cells per subject:", tab.groupby("subject").pE.sum().to_dict())
print("State labels in the files:", sorted(set(np.concatenate([s.labels for s in sessions]))))

# Units, ranges and timing of what is used.
for s in sessions:
    dur = s.stops - s.starts
    assert np.all(np.diff(s.starts) > 0) and np.all(dur >= 0)
    assert all(np.all(np.diff(x) >= 0) for x in s.spikes)
gaps = np.concatenate([s.starts[1:] - s.stops[:-1] for s in sessions])
allstart = np.concatenate([s.starts for s in sessions])
print("\nState table: start times are whole seconds:", bool(np.all(allstart == np.round(allstart))),
      "| gaps between consecutive epochs (s), values and counts:",
      dict(zip(*[x.tolist() for x in np.unique(gaps, return_counts=True)])))
print("State durations (s): min", min((s.stops - s.starts).min() for s in sessions),
      "max", max((s.stops - s.starts).max() for s in sessions))
print("Spikes after the end of the state table (fraction of spikes), per session:",
      {s.session_id: round(float(np.mean(np.concatenate(s.spikes) > s.stops[-1])), 3) for s in sessions})
isi_viol = []
for s in sessions:
    for u in s.pe:
        x = s.spikes[u]
        if len(x) > 1:
            isi_viol.append(np.mean(np.diff(x) < 0.001))
isi_viol = np.array(isi_viol)
print("pE units: fraction of intervals under 1 ms, median %.4f, 95th percentile %.4f, max %.4f"
      % (np.median(isi_viol), np.percentile(isi_viol, 95), isi_viol.max()))
print("pE spike counts per unit: min", min(len(s.spikes[u]) for s in sessions for u in s.pe),
      "median", int(np.median([len(s.spikes[u]) for s in sessions for u in s.pe])))
print("LFP reference column per session:", {s.session_id: s.lfp_channel for s in sessions})

# The excluded file: spike times are whole seconds divided by 20000.
p = list(excluded)[0]
f = h5py.File(remfile.File(urls[p]), "r")
t = f["units/spike_times"][()]
print("\n%s: %d spikes, %d distinct values, max %.4f s; max x 20000 = %.0f; times x 20000 all integers: %s;"
      " LFP lasts %.0f s, states end at %.0f s" % (
          p, len(t), len(np.unique(t)), t.max(), t.max() * 20000, bool(np.allclose(t * 20000, np.round(t * 20000), atol=1e-6)),
          f["processing/ecephys/LFP/lfp/data"].shape[0] / 1250, f["processing/behavior/states/stop_time"][()].max()))
print("LFP data attributes: unit %r, conversion %r" % (
    f["processing/ecephys/LFP/lfp/data"].attrs["unit"], f["processing/ecephys/LFP/lfp/data"].attrs["conversion"]))
el = f["general/extracellular_ephys/electrodes"]
print("electrode location values:", set(ss._str(el["location"][()])),
      "| electrode group locations:", sorted({f["general/extracellular_ephys"][g].attrs.get("location", "")
                                             for g in f["general/extracellular_ephys"] if g.startswith("shank")}))
print("subject age, weight:", f["general/subject/age"][()], f["general/subject/weight"][()])

# Recordings deposited more than once: identical spike trains or state tables across files.
sig_units = collections.defaultdict(list)
for s in sessions:
    for u, x in enumerate(s.spikes):
        if len(x) >= 5:
            sig_units[(len(x), tuple(np.round(x[:5], 5)))].append((s.session_id, u))
dups = [v for v in sig_units.values() if len({a for a, _ in v}) > 1]
print("\nUnits with identical spike trains in different files:", dups if dups else "none")
st_sig = collections.defaultdict(list)
for s in sessions:
    st_sig[(len(s.starts), tuple(s.starts[:10]), tuple(s.stops[-3:]))].append(s.session_id)
print("Identical state tables in different files:", [v for v in st_sig.values() if len(v) > 1] or "none")
same_counts = collections.Counter(tuple(sorted(len(x) for x in s.spikes)) for s in sessions)
print("Files with identical sets of per-unit spike counts:", [k for k, v in same_counts.items() if v > 1] or "none")

      subject           session   pE  pI  unknown  states_min  spikes_first  spikes_last_min  lfp_min  WAKE_min  MA_n  NREM_min  REM_min  epochs  lfp_packets
0     BWRat17    BWRat17_121712   44   6        0       100.8          0.04            101.0    101.0      72.2     6      24.2      1.8      33           13
1     BWRat17    BWRat17_121912   29   0        0       108.7          0.02            108.9    108.9      37.0    16      61.1      3.6      71           31
2     BWRat18    BWRat18_020513   35   5        0       126.4          0.03            126.5    126.5      42.8     2      73.2      8.9      53           19
3     BWRat19    BWRat19_032413   37  12        0        83.4          0.01             83.6     83.6      43.4     9      35.7      0.2      34           16
4     BWRat19    BWRat19_032513   50   9        0        65.3          0.01             65.5     65.5      19.6     9      37.0      4.8      42           18
5     BWRat20    BWRat20_101013   13   5        0       222.6          0.00            222.7    222.7     136.4     9      69.4     12.8      62           26
6     BWRat20    BWRat20_101513   58   3        0       383.4          0.01            383.5    383.5     126.4    47     201.8     38.5     298          107
7     BWRat21    BWRat21_121113   28   4        0       364.6          0.00            364.7    364.7     119.9    17     171.0     64.1     213           67
8     BWRat21    BWRat21_121613   20   1       23       358.9          0.13            359.0    359.0     144.5    20     152.4     51.6     183           59
9     BWRat21    BWRat21_121813   11   1        0       361.6          0.45            361.8    361.8     143.7    11     155.5     56.4     148           49
10      Bogey      Bogey_012615   32   4        0       404.2          0.01            404.4    404.4     146.2    42     214.4     22.3     218           91
11       Dino       Dino_061814  105   9        0       327.3          0.00            327.4    327.4     148.2    23     142.3     27.6     161           61
12       Dino       Dino_061914   55   5        0       291.8          0.00            292.0    331.7     140.8    13     125.1     19.7     106           42
13       Dino       Dino_062014   32   4        0       372.6          0.06            372.7    372.7     168.2    20     162.8     33.2     146           60
14       Dino       Dino_072314   24   0        0       171.8          0.03            172.0    172.0      71.0     7      73.9     23.4      67           24
15       Dino       Dino_072414   47   6        0       220.0          0.00            220.1    220.1      89.3     6     106.6     21.0      85           26
16      Rizzo      Rizzo_022615   25   1        0       675.8          0.03            675.9    675.9     260.6    37     359.6     38.0     287          110
17      Rizzo      Rizzo_022715   89   4        0       432.0          0.00            432.2    432.2     215.0    20     179.5     29.1     154           61
18   Splinter   Splinter_020515   26   1        0       397.1          0.02            397.3    397.3     131.2    17     221.4     36.9     144           53
19   Splinter   Splinter_020915   42   9        0       391.3          0.00            391.5    391.5     140.8    22     215.7     25.3     127           51
20  Templeton  Templeton_032415   10   0        0       317.4          0.10            317.6    317.6     117.4    38     146.4     37.0     202           85

Totals over the 21 usable sessions: {'pE': 812, 'pI': 89, 'unknown': 23, 'MA_n': 391}
Subjects: 10 | sessions per subject: {np.str_('BWRat17'): 2, np.str_('BWRat18'): 1, np.str_('BWRat19'): 2, np.str_('BWRat20'): 2, np.str_('BWRat21'): 3, np.str_('Bogey'): 1, np.str_('Dino'): 5, np.str_('Rizzo'): 2, np.str_('Splinter'): 2, np.str_('Templeton'): 1}
pE cells per subject: {np.str_('BWRat17'): 73, np.str_('BWRat18'): 35, np.str_('BWRat19'): 87, np.str_('BWRat20'): 71, np.str_('BWRat21'): 59, np.str_('Bogey'): 32, np.str_('Dino'): 263, np.str_('Rizzo'): 114, np.str_('Splinter'): 68, np.str_('Templeton'): 10}
State labels in the files: [np.str_('Awake'), np.str_('Non-REM'), np.str_('REM')]

State table: start times are whole seconds: True | gaps between consecutive epochs (s), values and counts: {1.0: 2813}
State durations (s): min 6.0 max 3952.0
Spikes after the end of the state table (fraction of spikes), per session: {np.str_('BWRat17_121712'): 0.001, np.str_('BWRat17_121912'): 0.001, np.str_('BWRat18_020513'): 0.0, np.str_('BWRat19_032413'): 0.001, np.str_('BWRat19_032513'): 0.002, np.str_('BWRat20_101013'): 0.0, np.str_('BWRat20_101513'): 0.0, np.str_('BWRat21_121113'): 0.0, np.str_('BWRat21_121613'): 0.0, np.str_('BWRat21_121813'): 0.0, np.str_('Bogey_012615'): 0.0, np.str_('Dino_061814'): 0.0, np.str_('Dino_061914'): 0.0, np.str_('Dino_062014'): 0.0, np.str_('Dino_072314'): 0.0, np.str_('Dino_072414'): 0.0, np.str_('Rizzo_022615'): 0.0, np.str_('Rizzo_022715'): 0.0, np.str_('Splinter_020515'): 0.0, np.str_('Splinter_020915'): 0.0, np.str_('Templeton_032415'): 0.001}
pE units: fraction of intervals under 1 ms, median 0.0000, 95th percentile 0.0003, max 0.0018
pE spike counts per unit: min 0 median 10443
LFP reference column per session: {np.str_('BWRat17_121712'): 13, np.str_('BWRat17_121912'): 13, np.str_('BWRat18_020513'): 59, np.str_('BWRat19_032413'): 54, np.str_('BWRat19_032513'): 24, np.str_('BWRat20_101013'): 59, np.str_('BWRat20_101513'): 61, np.str_('BWRat21_121113'): 7, np.str_('BWRat21_121613'): 12, np.str_('BWRat21_121813'): 12, np.str_('Bogey_012615'): 46, np.str_('Dino_061814'): 92, np.str_('Dino_061914'): 83, np.str_('Dino_062014'): 83, np.str_('Dino_072314'): 106, np.str_('Dino_072414'): 106, np.str_('Rizzo_022615'): 17, np.str_('Rizzo_022715'): 52, np.str_('Splinter_020515'): 42, np.str_('Splinter_020915'): 43, np.str_('Templeton_032415'): 27}

sub-Dino/sub-Dino_ses-Dino-072114_ecephys.nwb: 215233 spikes, 8451 distinct values, max 0.4252 s; max x 20000 = 8504; times x 20000 all integers: True; LFP lasts 7959 s, states end at 7953 s
LFP data attributes: unit 'volts', conversion np.float64(nan)

electrode location values: {np.str_('unknown')} | electrode group locations: ['mPFC']
subject age, weight: b'3-7 months' b'250-500g'

Units with identical spike trains in different files: none
Identical state tables in different files: none
Files with identical sets of per-unit spike counts: none

# LFP: the lag-1 sample autocorrelation of the reference channel, from the raw samples kept at
# load time. LFP at 1250 Hz is smooth (near 1); scrambled samples are not.
lag1 = {s.session_id: float(np.mean([np.corrcoef(x[:-1].astype(float), x[1:].astype(float))[0, 1] for _, x in s.lfp_raw]))
        for s in sessions}
print("LFP lag-1 autocorrelation per session:", {k: round(v, 3) for k, v in lag1.items()})
lfp_bad = [s.session_id for s in sessions if lag1[s.session_id] < 0.8]
print("Sessions whose LFP is left out of the LFP analyses (lag-1 autocorrelation under 0.8):", lfp_bad)
for s in sessions:
    if s.session_id in lfp_bad:
        print("  %s: LFP has %d columns for %d electrodes and lasts %.1f min; states end at %.1f min" % (
            s.session_id, len(h5py.File(remfile.File(urls[s.path]), "r")["processing/ecephys/LFP/lfp/electrodes"][()]),
            len(h5py.File(remfile.File(urls[s.path]), "r")["general/extracellular_ephys/electrodes/id"][()]),
            s.lfp_samples / ss.LFP_RATE / 60, s.stops[-1] / 60))
        s.lfp = {}

LFP lag-1 autocorrelation per session: {np.str_('BWRat17_121712'): 0.998, np.str_('BWRat17_121912'): 0.998, np.str_('BWRat18_020513'): 0.999, np.str_('BWRat19_032413'): 0.998, np.str_('BWRat19_032513'): 0.998, np.str_('BWRat20_101013'): 0.995, np.str_('BWRat20_101513'): 0.993, np.str_('BWRat21_121113'): 0.996, np.str_('BWRat21_121613'): 0.997, np.str_('BWRat21_121813'): 0.996, np.str_('Bogey_012615'): 0.997, np.str_('Dino_061814'): 0.993, np.str_('Dino_061914'): 0.28, np.str_('Dino_062014'): 0.998, np.str_('Dino_072314'): 0.996, np.str_('Dino_072414'): 0.995, np.str_('Rizzo_022615'): 0.997, np.str_('Rizzo_022715'): 0.989, np.str_('Splinter_020515'): 0.999, np.str_('Splinter_020915'): 0.999, np.str_('Templeton_032415'): 0.999}
Sessions whose LFP is left out of the LFP analyses (lag-1 autocorrelation under 0.8): [np.str_('Dino_061914')]

  Dino_061914: LFP has 125 columns for 128 electrodes and lasts 331.7 min; states end at 291.9 min

shifts = np.array([-300, -120, -60, -20, -5, 0, 5, 20, 60, 120, 300])
clock = []
for s in sessions:
    pop = np.sort(np.concatenate([s.spikes[u] for u in s.pe]))
    edges50 = np.arange(0, s.stops[-1], 0.05)
    silent = np.histogram(pop, edges50)[0] == 0
    centres = edges50[:-1] + 0.025
    row = {"session": s.session_id, "pE": len(s.pe)}
    for d in shifts:
        k = np.searchsorted(s.starts + d, centres, side="right") - 1
        ok = (k >= 0) & (centres < (s.stops + d)[np.clip(k, 0, None)])
        lab = np.where(ok, s.state[np.clip(k, 0, None)], "none")
        row[int(d)] = silent[lab == "NREM"].mean() - silent[lab == "WAKE"].mean()
    # LFP against spikes: cross-correlation of LFP and pE spike count in 8 ms bins (10 samples)
    lags = np.arange(-125, 126)
    acc = np.zeros(len(lags)); nseg = 0
    for t0, x in s.lfp_raw:
        nb = len(x) // 10
        lf = x[:nb * 10].reshape(nb, 10).mean(axis=1)
        sp = np.diff(np.searchsorted(pop, t0 + np.arange(nb + 1) * 0.008))
        lf = (lf - lf.mean()) / lf.std(); sp = (sp - sp.mean()) / (sp.std() or 1)
        acc += np.array([np.mean(lf[max(0, -L):nb - max(0, L)] * sp[max(0, L):nb - max(0, -L)]) for L in lags]); nseg += 1
    cc = acc / nseg
    j = np.argmax(np.abs(cc))
    row.update({"lfp seconds": round(sum(len(x) for _, x in s.lfp_raw) / ss.LFP_RATE), "xcorr extreme lag ms": lags[j] * 8,
                "xcorr extreme": round(cc[j], 3), "xcorr at -1 s": round(cc[0], 3), "xcorr at +1 s": round(cc[-1], 3)})
    clock.append(row)
clock = pd.DataFrame(clock)
best = clock[[int(d) for d in shifts]].idxmax(axis=1)
clock["best shift"] = best
clock["labels consistent"] = (best.abs() <= 60) & (clock[0] > clock[[-300, 300]].max(axis=1))
clock["lfp consistent"] = (clock["xcorr extreme lag ms"].abs() <= 50) & (clock["xcorr extreme"].abs() > 3 * clock[["xcorr at -1 s", "xcorr at +1 s"]].abs().max(axis=1))
print("Contrast of silent-bin fraction, NREM minus WAKE, by shift of the state table (s), and the LFP-spike cross-correlation:")
print(clock.round(3).to_string())
print("\nState labels consistent with the spikes (largest contrast within 60 s of zero shift and above the contrast at +-300 s):",
      int(clock["labels consistent"].sum()), "of", len(clock), "| not established:", list(clock.session[~clock["labels consistent"]]))
print("LFP consistent with the spikes (largest |r| within 50 ms of zero lag and over 3 times |r| at +-1 s):",
      int(clock["lfp consistent"].sum()), "of", len(clock), "| not established:", list(clock.session[~clock["lfp consistent"]]))
print("Shift with the largest contrast, counts:", best.value_counts().sort_index().to_dict())

Contrast of silent-bin fraction, NREM minus WAKE, by shift of the state table (s), and the LFP-spike cross-correlation:
             session   pE   -300   -120    -60    -20     -5      0      5     20     60    120    300  lfp seconds  xcorr extreme lag ms  xcorr extreme  xcorr at -1 s  xcorr at +1 s best shift  labels consistent  lfp consistent
0     BWRat17_121712   44  0.057  0.114  0.128  0.159  0.168  0.167  0.161  0.144  0.125  0.110  0.092          310                     0         -0.150          0.014         -0.010         -5               True            True
1     BWRat17_121912   29  0.050  0.054  0.075  0.100  0.106  0.103  0.102  0.086  0.074  0.064  0.055          334                     0         -0.057          0.021         -0.016         -5               True           False
2     BWRat18_020513   35  0.108  0.131  0.150  0.159  0.159  0.156  0.161  0.172  0.167  0.160  0.115          494                   -40         -0.118          0.006          0.003         20               True            True
3     BWRat19_032413   37  0.127  0.166  0.179  0.230  0.257  0.257  0.251  0.232  0.198  0.191  0.153          316                   -16         -0.169          0.003         -0.026          0               True            True
4     BWRat19_032513   50  0.044  0.099  0.128  0.156  0.169  0.172  0.170  0.160  0.129  0.105  0.008          322                   -16         -0.189         -0.013          0.022          0               True            True
5     BWRat20_101013   13  0.159  0.227  0.259  0.291  0.308  0.308  0.307  0.297  0.264  0.228  0.150          710                     8         -0.142          0.003         -0.006          0               True            True
6     BWRat20_101513   58  0.034  0.076  0.100  0.134  0.151  0.154  0.149  0.133  0.102  0.076  0.029          304                     0         -0.113         -0.002          0.019          0               True            True
7     BWRat21_121113   28  0.014  0.086  0.134  0.165  0.158  0.153  0.148  0.131  0.077  0.039 -0.004          374                   -64         -0.094         -0.000         -0.000        -20               True           False
8     BWRat21_121613   20 -0.019 -0.011  0.007  0.029  0.034  0.032  0.027  0.014  0.004 -0.005 -0.021          391                   -88         -0.042          0.022          0.002         -5               True           False
9     BWRat21_121813   11 -0.071 -0.062 -0.058 -0.065 -0.072 -0.071 -0.067 -0.064 -0.068 -0.077 -0.077          665                   -72         -0.034          0.015         -0.000        -60              False           False
10      Bogey_012615   32  0.070  0.100  0.124  0.171  0.182  0.177  0.170  0.148  0.108  0.085  0.054          505                   -64         -0.165         -0.001         -0.031         -5               True           False
11       Dino_061814  105  0.030  0.065  0.095  0.129  0.141  0.141  0.138  0.127  0.104  0.076  0.034          355                    16         -0.073          0.005         -0.012          0               True            True
12       Dino_061914   55  0.059  0.087  0.107  0.126  0.134  0.135  0.133  0.126  0.113  0.095  0.053          936                   128         -0.011         -0.002          0.004          0               True           False
13       Dino_062014   32  0.044  0.091  0.110  0.113  0.110  0.108  0.105  0.107  0.098  0.079  0.053          613                     0         -0.044         -0.012          0.008        -20               True            True
14       Dino_072314   24  0.007  0.047  0.068  0.065  0.064  0.066  0.068  0.065  0.055  0.034 -0.005          341                   -32         -0.076         -0.012         -0.010        -60               True            True
15       Dino_072414   47  0.030  0.068  0.087  0.100  0.104  0.104  0.102  0.098  0.087  0.072  0.039          425                   -56         -0.099          0.005         -0.006         -5               True           False
16      Rizzo_022615   25  0.026  0.078  0.119  0.155  0.159  0.156  0.151  0.138  0.110  0.078  0.031          556                    -8         -0.162         -0.001         -0.004         -5               True            True
17      Rizzo_022715   89  0.065  0.090  0.105  0.120  0.125  0.124  0.121  0.115  0.105  0.093  0.066          333                     0         -0.233          0.050         -0.016         -5               True            True
18   Splinter_020515   26  0.052  0.114  0.142  0.147  0.134  0.129  0.125  0.121  0.102  0.078  0.043          390                   -40         -0.259          0.004         -0.000        -20               True            True
19   Splinter_020915   42  0.144  0.174  0.192  0.205  0.200  0.200  0.201  0.198  0.191  0.184  0.162          300                   -32         -0.101         -0.004          0.006        -20               True            True
20  Templeton_032415   10  0.015  0.054  0.075  0.074  0.054  0.051  0.051  0.046  0.022  0.017 -0.007          455                     0         -0.032         -0.001         -0.001        -60               True            True

State labels consistent with the spikes (largest contrast within 60 s of zero shift and above the contrast at +-300 s): 20 of 21 | not established: [np.str_('BWRat21_121813')]
LFP consistent with the spikes (largest |r| within 50 ms of zero lag and over 3 times |r| at +-1 s): 14 of 21 | not established: [np.str_('BWRat17_121912'), np.str_('BWRat21_121113'), np.str_('BWRat21_121613'), np.str_('BWRat21_121813'), np.str_('Bogey_012615'), np.str_('Dino_061914'), np.str_('Dino_072414')]
Shift with the largest contrast, counts: {-60: 3, -20: 4, -5: 7, 0: 6, 20: 1}

n_files = len(sessions) + len(excluded)
dep = dict(rats=tab.subject.nunique(), sessions=n_files)
# the excluded file's units are counted from its own table
ct_ex = collections.Counter(ss._str(f["units/cell_type"][()]))
st_ex = f["processing/behavior/states"]
dur_ex = (st_ex["stop_time"][()].max() - st_ex["start_time"][()].min()) / 60
durs = list(tab.states_min) + [dur_ex]
dep.update(mean_duration_min=np.mean(durs), sd_duration_min=np.std(durs, ddof=1),
           pE=int(tab.pE.sum() + ct_ex["excitatory"]), pI=int(tab.pI.sum() + ct_ex["inhibitory"]),
           unknown=int(tab.unknown.sum() + ct_ex["unknown"]))
print("As deposited (22 files):", {k: (round(v, 1) if isinstance(v, float) else v) for k, v in dep.items()})
print("Usable (21 sessions): rats %d, sessions %d, pE %d, pI %d" % (tab.subject.nunique(), len(tab), tab.pE.sum(), tab.pI.sum()))
print("Mean spike span of the usable sessions: %.1f min" % tab.spikes_last_min.mean())
exp = dict(rats=(11, 1.1), sessions=(27, 2.7), mean_duration_min=(294, 137), pE=(995, 99.5), pI=(126, 12.6))
ok = {k: abs(dep[k] - v) <= t for k, (v, t) in exp.items()}
print("Within tolerance:", ok, "->", sum(ok.values()), "of 5")
print("pE cells per session: mean %.1f, range %d to %d (paper: mean 36.9, range 10 to 104)" % (
    tab.pE.mean(), tab.pE.min(), tab.pE.max()))
print("pI cells per session: mean %.2f, range %d to %d (paper: mean 4.67, range 0 to 13)" % (
    tab.pI.mean(), tab.pI.min(), tab.pI.max()))

As deposited (22 files): {'rats': 10, 'sessions': 22, 'mean_duration_min': np.float64(286.8), 'sd_duration_min': np.float64(150.2), 'pE': 826, 'pI': 96, 'unknown': 23}
Usable (21 sessions): rats 10, sessions 21, pE 812, pI 89
Mean spike span of the usable sessions: 294.3 min
Within tolerance: {'rats': True, 'sessions': False, 'mean_duration_min': np.True_, 'pE': False, 'pI': False} -> 2 of 5
pE cells per session: mean 38.7, range 10 to 105 (paper: mean 36.9, range 10 to 104)
pI cells per session: mean 4.24, range 0 to 12 (paper: mean 4.67, range 0 to 13)

def cell_table(wake_states):
    out = []
    for s in sessions:
        for cls, units in (("pE", s.pe), ("pI", s.pi)):
            r = {lab: ss.state_rates(s, units, sts) for lab, sts in
                 (("WAKE", wake_states), ("NREM", ["NREM"]), ("REM", ["REM"]))}
            for k, u in enumerate(units):
                out.append(dict(subject=s.subject, session=s.session_id, unit=u, cls=cls,
                                WAKE=r["WAKE"][k], NREM=r["NREM"][k], REM=r["REM"][k]))
    return pd.DataFrame(out)

cells_a = cell_table(["WAKE"])
cells_b = cell_table(["WAKE", "MA"])
expected = {("pE", "WAKE"): 0.76, ("pE", "NREM"): 0.69, ("pE", "REM"): 0.88,
            ("pI", "WAKE"): 5.59, ("pI", "NREM"): 4.69, ("pI", "REM"): 4.25}
for name, cells in (("(a) WAKE without MA", cells_a), ("(b) WAKE with MA", cells_b)):
    med = cells.groupby("cls")[["WAKE", "NREM", "REM"]].median()
    sd = cells.groupby("cls")[["WAKE", "NREM", "REM"]].std()
    n_in = sum(abs(med.loc[c, s_] - v) <= 0.1 * v for (c, s_), v in expected.items())
    print(name, "| medians (Hz):\n", med.round(3).to_string(), "\nSD:\n", sd.round(2).to_string(),
          "\nwithin 10%:", n_in, "of 6\n")
print("n cells:", cells_a.cls.value_counts().to_dict())
print("Per subject, pE medians (reading a):")
print(cells_a[cells_a.cls == "pE"].groupby("subject")[["WAKE", "NREM", "REM"]].median().round(3).to_string())
print("\nPer subject, pI medians (reading a):")
print(cells_a[cells_a.cls == "pI"].groupby("subject")[["WAKE", "NREM", "REM"]].agg(["median", "size"]).round(2).to_string())

(a) WAKE without MA | medians (Hz):
       WAKE   NREM    REM
cls                     
pE   0.755  0.662  0.825
pI   3.917  2.730  3.062 
SD:
      WAKE  NREM   REM
cls                  
pE   1.97  1.33  2.35
pI   6.89  5.07  8.39 
within 10%: 3 of 6

(b) WAKE with MA | medians (Hz):
       WAKE   NREM    REM
cls                     
pE   0.748  0.662  0.825
pI   3.912  2.730  3.062 
SD:
      WAKE  NREM   REM
cls                  
pE   1.95  1.33  2.35
pI   6.85  5.07  8.39 
within 10%: 3 of 6

n cells: {'pE': 812, 'pI': 89}
Per subject, pE medians (reading a):
            WAKE   NREM    REM
subject                       
BWRat17    0.573  0.445  0.620
BWRat18    0.606  0.444  0.208
BWRat19    1.119  0.867  0.855
BWRat20    0.489  0.420  0.580
BWRat21    0.370  0.482  0.617
Bogey      0.746  0.698  0.813
Dino       0.754  0.738  1.034
Rizzo      1.281  1.051  1.364
Splinter   0.517  0.367  0.630
Templeton  0.477  0.512  1.192

Per subject, pI medians (reading a):
           WAKE        NREM         REM     
         median size median size median size
subject                                     
BWRat17    3.54    6   2.76    6   2.05    6
BWRat18    3.23    5   0.21    5   0.07    5
BWRat19    4.63   21   2.76   21   3.59   21
BWRat20    2.79    8   1.20    8   1.70    8
BWRat21    3.46    6   3.18    6   3.84    6
Bogey      3.31    4   2.21    4   3.09    4
Dino       5.41   24   3.90   24   4.49   24
Rizzo     13.37    5  13.01    5  13.18    5
Splinter   3.61   10   2.81   10   2.57   10

pe_a = cells_a[cells_a.cls == "pE"]
for state in ("WAKE", "NREM", "REM"):
    r = pe_a[state].values
    lo, hi = np.percentile(r, [2.5, 97.5])
    pos = r[r > 0]
    print("%s: skewness %.2f, skewness of log10 (rates > 0) %.2f, 2.5th pct %.4f Hz, 97.5th pct %.3f Hz, "
          "min %.4f, max %.2f, zero rates %d of %d" % (state, st.skew(r), st.skew(np.log10(pos)), lo, hi,
                                                        r.min(), r.max(), (r == 0).sum(), len(r)))
r = pe_a.WAKE.values
lo, hi = np.percentile(r, [2.5, 97.5])
print("\nWAKE: skewness > 1:", st.skew(r) > 1, "| 2.5th pct within 0.025-0.1 Hz:", 0.025 <= lo <= 0.1,
      "| 97.5th pct within 2.5-10 Hz:", 2.5 <= hi <= 10)
fig, ax = plt.subplots(figsize=(5, 3.5))
for state, c in (("WAKE", "k"), ("NREM", "tab:blue"), ("REM", "tab:red")):
    x = np.sort(pe_a[state].values)
    x = x[x > 0]
    ax.plot(x, np.arange(1, len(x) + 1) / len(x), color=c, label=state)
for e in edges[1:-1]:
    ax.axvline(e, color="0.7", lw=0.5)
ax.set_xscale("log"); ax.set_xlabel("rate (Hz)"); ax.set_ylabel("cumulative fraction of pE cells"); ax.legend()
ax.set_title("pE rates by state (zero rates not drawn)")
plt.show()
print("Rates drawn: from %.4f to %.2f Hz; zero rates not drawn: WAKE %d, NREM %d, REM %d" % (
    min(pe_a[s_][pe_a[s_] > 0].min() for s_ in ("WAKE", "NREM", "REM")),
    max(pe_a[s_].max() for s_ in ("WAKE", "NREM", "REM")),
    (pe_a.WAKE == 0).sum(), (pe_a.NREM == 0).sum(), (pe_a.REM == 0).sum()))
print("WAKE-rate sextile edges (Hz):", np.round(edges, 4))

WAKE: skewness 4.57, skewness of log10 (rates > 0) -0.74, 2.5th pct 0.0124 Hz, 97.5th pct 7.165 Hz, min 0.0000, max 20.86, zero rates 13 of 812
NREM: skewness 6.58, skewness of log10 (rates > 0) -0.97, 2.5th pct 0.0169 Hz, 97.5th pct 3.887 Hz, min 0.0000, max 20.24, zero rates 13 of 812
REM: skewness 6.29, skewness of log10 (rates > 0) -0.73, 2.5th pct 0.0018 Hz, 97.5th pct 5.448 Hz, min 0.0000, max 30.00, zero rates 20 of 812

WAKE: skewness > 1: True | 2.5th pct within 0.025-0.1 Hz: False | 97.5th pct within 2.5-10 Hz: True

Rates drawn: from 0.0009 to 30.00 Hz; zero rates not drawn: WAKE 13, NREM 13, REM 20
WAKE-rate sextile edges (Hz): [ 0.      0.1781  0.4247  0.7545  1.1668  2.0164 20.8578]

frac = {s_: 100 * np.mean(pe_a[s_] > 2) for s_ in ("WAKE", "REM", "NREM")}
print("Percent of pE cells above 2 Hz:", {k: round(v, 2) for k, v in frac.items()})
for other in ("WAKE", "REM"):
    a = (pe_a.NREM > 2).sum(); b = (pe_a[other] > 2).sum(); n = len(pe_a)
    chi2, p, _, _ = st.chi2_contingency([[a, n - a], [b, n - b]])
    print("NREM against %s: chi2 = %.1f, p = %.2g" % (other, chi2, p))
print("Within tolerance:", {k: abs(frac[k] - v) <= 0.1 * v for k, v in (("WAKE", 19.2), ("REM", 19.2), ("NREM", 8.8))})
print("\nPer subject, percent above 2 Hz:")
print((pe_a.groupby("subject")[["WAKE", "REM", "NREM"]].agg(lambda x: 100 * np.mean(x > 2))).round(1).to_string())

Percent of pE cells above 2 Hz: {'WAKE': np.float64(17.0), 'REM': np.float64(18.72), 'NREM': np.float64(9.11)}
NREM against WAKE: chi2 = 21.5, p = 3.5e-06
NREM against REM: chi2 = 30.5, p = 3.4e-08
Within tolerance: {'WAKE': np.False_, 'REM': np.True_, 'NREM': np.True_}

Per subject, percent above 2 Hz:
           WAKE   REM  NREM
subject                    
BWRat17     4.1   4.1   1.4
BWRat18     2.9   2.9   2.9
BWRat19    37.9  31.0  14.9
BWRat20    19.7  25.4  12.7
BWRat21     0.0  10.2   0.0
Bogey      12.5  12.5   6.2
Dino       14.4  20.2   8.0
Rizzo      31.6  30.7  19.3
Splinter   13.2   7.4   7.4
Templeton   0.0   0.0   0.0

def sd_ratio(x, y, groups, seed, n_boot=2000):
    """SD of y over SD of x (the symmetric slope), the Pitman-Morgan p for equal variances of paired
    values, and a 95% bootstrap interval resampling the groups (sessions or cycles)."""
    rng = np.random.default_rng(seed)
    x, y = np.asarray(x, float), np.asarray(y, float)
    _, inv = np.unique(np.asarray(groups), return_inverse=True)
    idx = [np.where(inv == g)[0] for g in range(inv.max() + 1)]
    boot = []
    for _ in range(n_boot):
        sel = np.concatenate([idx[j] for j in rng.integers(0, len(idx), len(idx))])
        boot.append(np.std(y[sel], ddof=1) / np.std(x[sel], ddof=1))
    return np.std(y, ddof=1) / np.std(x, ddof=1), st.pearsonr(x + y, x - y).pvalue, *np.percentile(boot, [2.5, 97.5])

def ols(x, y):
    res = st.linregress(x, y)
    t = st.t.ppf(0.975, len(x) - 2)
    return res.slope, res.slope - t * res.stderr, res.slope + t * res.stderr, res.rvalue, res.pvalue, len(x)

keep = (pe_a.WAKE > 0) & (pe_a.NREM > 0)
fits = {"(a) log10": ols(np.log10(pe_a.WAKE[keep]), np.log10(pe_a.NREM[keep])),
        "(b) Hz": ols(pe_a.WAKE, pe_a.NREM),
        "(a) reverse, WAKE on NREM": ols(np.log10(pe_a.NREM[keep]), np.log10(pe_a.WAKE[keep])),
        "(b) reverse": ols(pe_a.NREM, pe_a.WAKE)}
for k, (b, lo, hi, r, p, n) in fits.items():
    print("%s: slope %.3f (95%% CI %.3f to %.3f), R = %.3f, p = %.2g, n = %d" % (k, b, lo, hi, r, p, n))
print("Dropped for zero rate in the log reading:", (~keep).sum())
for k in ("(a) log10", "(b) Hz"):
    b, lo, hi, r, p, n = fits[k]
    print(k, "| R within 0.738-0.902:", 0.738 <= r <= 0.902, "| slope within 0.594-0.781:", 0.594 <= b <= 0.781,
          "| CI excludes 1:", hi < 1)
print("\nPer subject, reading (a):")
for subj, g in pe_a[keep].groupby("subject"):
    b, lo, hi, r, p, n = ols(np.log10(g.WAKE), np.log10(g.NREM))
    print("  %s: slope %.2f (%.2f to %.2f), R %.2f, n %d" % (subj, b, lo, hi, r, n))
fig, ax = plt.subplots(figsize=(4, 4))
ax.loglog(pe_a.WAKE[keep], pe_a.NREM[keep], ".", ms=3)
ax.plot([1e-3, 30], [1e-3, 30], "k--", lw=0.8)
ax.set_xlabel("WAKE rate (Hz)"); ax.set_ylabel("nonREM rate (Hz)"); ax.set_title("pE cells")
plt.show()
print("Plotted: WAKE %.4f to %.2f Hz, NREM %.4f to %.2f Hz; %d cells with a zero rate not drawn" % (
    pe_a.WAKE[keep].min(), pe_a.WAKE[keep].max(), pe_a.NREM[keep].min(), pe_a.NREM[keep].max(), (~keep).sum()))
# Control for regression to the mean: each cell's WAKE epochs are split alternately into two halves.
# The slope of one half on the other shows what measurement noise alone does to a slope, and the
# reliability of the whole WAKE rate (Spearman-Brown) gives the slope corrected for attenuation.
half = []
for s in sessions:
    w_ = np.where(s.state == "WAKE")[0]
    half += list(zip(ss.rates(s.spikes, s.pe, s.starts[w_[0::2]], s.stops[w_[0::2]]),
                     ss.rates(s.spikes, s.pe, s.starts[w_[1::2]], s.stops[w_[1::2]])))
half = np.array(half)
kh = (half > 0).all(axis=1)
b, lo, hi, r_h, p, n = ols(np.log10(half[kh, 0]), np.log10(half[kh, 1]))
rel = 2 * r_h / (1 + r_h)
print("\nSplit-half control, alternate WAKE epochs: slope of one half on the other %.3f (95%% CI %.3f to %.3f), "
      "R = %.3f, n = %d (%d cells with a zero rate in a half left out)" % (b, lo, hi, r_h, n, (~kh).sum()))
print("Reliability of the log10 WAKE rate (Spearman-Brown) %.3f; nonREM-on-WAKE slope corrected for attenuation %.3f "
      "(%.3f to %.3f, the reliability's own uncertainty not counted)" % (rel, *(v / rel for v in fits["(a) log10"][:3])))
# Regression to the mean comes from any correlation below 1, not only from measurement noise. The
# ratio of the SDs (nonREM over WAKE, log10) is the symmetric slope: below 1 only if the spread of
# rates is narrower in nonREM.
v = sd_ratio(np.log10(pe_a.WAKE[keep]), np.log10(pe_a.NREM[keep]), pe_a.session[keep], seed=21)
print("SD of log10 rates, nonREM over WAKE: %.3f (Pitman-Morgan p = %.2g; 95%% interval resampling sessions %.3f to %.3f)" % v)

(a) log10: slope 0.757 (95% CI 0.727 to 0.788), R = 0.862, p = 9.1e-238, n = 799
(b) Hz: slope 0.583 (95% CI 0.560 to 0.607), R = 0.865, p = 1.8e-244, n = 812
(a) reverse, WAKE on NREM: slope 0.982 (95% CI 0.942 to 1.022), R = 0.862, p = 9.1e-238, n = 799
(b) reverse: slope 1.282 (95% CI 1.231 to 1.333), R = 0.865, p = 1.8e-244, n = 812
Dropped for zero rate in the log reading: 13
(a) log10 | R within 0.738-0.902: True | slope within 0.594-0.781: True | CI excludes 1: True
(b) Hz | R within 0.738-0.902: True | slope within 0.594-0.781: False | CI excludes 1: True

Per subject, reading (a):
  BWRat17: slope 0.82 (0.70 to 0.94), R 0.85, n 73
  BWRat18: slope 0.90 (0.66 to 1.14), R 0.80, n 35
  BWRat19: slope 0.49 (0.42 to 0.57), R 0.81, n 87
  BWRat20: slope 0.91 (0.85 to 0.96), R 0.97, n 71
  BWRat21: slope 0.88 (0.75 to 1.02), R 0.90, n 46
  Bogey: slope 0.41 (0.26 to 0.55), R 0.73, n 32
  Dino: slope 0.80 (0.74 to 0.86), R 0.87, n 263
  Rizzo: slope 0.79 (0.74 to 0.84), R 0.95, n 114
  Splinter: slope 0.69 (0.54 to 0.84), R 0.75, n 68
  Templeton: slope 0.49 (0.15 to 0.83), R 0.77, n 10

Plotted: WAKE 0.0009 to 20.86 Hz, NREM 0.0016 to 20.24 Hz; 13 cells with a zero rate not drawn

Split-half control, alternate WAKE epochs: slope of one half on the other 0.936 (95% CI 0.909 to 0.963), R = 0.924, n = 797 (15 cells with a zero rate in a half left out)
Reliability of the log10 WAKE rate (Spearman-Brown) 0.961; nonREM-on-WAKE slope corrected for attenuation 0.788 (0.756 to 0.821, the reliability's own uncertainty not counted)

SD of log10 rates, nonREM over WAKE: 0.878 (Pitman-Morgan p = 9.8e-13; 95% interval resampling sessions 0.777 to 0.971)

bi = []
for s in sessions:
    for u in s.pe:
        row = dict(subject=s.subject, session=s.session_id, unit=u)
        for state in ("WAKE", "NREM", "REM"):
            m = s.state == state
            for rd in ("a", "b"):
                row[f"{state}_{rd}"] = ss.burst_index(s.spikes[u], s.starts[m], s.stops[m], rd)
        bi.append(row)
bi = pd.DataFrame(bi)
for rd in ("a", "b"):
    cols = [f"{x}_{rd}" for x in ("NREM", "WAKE", "REM")]
    med = 100 * bi[cols].median()
    h, p = st.kruskal(*[bi[c].dropna() for c in cols])
    print("Reading (%s): medians %% NREM %.2f, WAKE %.2f, REM %.2f; n %s; Kruskal-Wallis H = %.1f, p = %.2g" % (
        rd, *med.values, bi[cols].notna().sum().to_dict(), h, p))
    full = bi[cols].dropna()
    print("   Friedman test on the %d cells with all three values (paired): p = %.2g" % (
        len(full), st.friedmanchisquare(*[full[c] for c in cols]).pvalue))
    print("   within tolerance:", [abs(m_ - v) <= 0.1 * v for m_, v in zip(med.values, (11.8, 4.0, 5.5))],
          "| order NREM > REM > WAKE:", med.iloc[0] > med.iloc[2] > med.iloc[1])
    pair = st.wilcoxon(bi[f"REM_{rd}"], bi[f"WAKE_{rd}"], nan_policy="omit")
    print("   REM against WAKE, paired Wilcoxon p = %.2g" % pair.pvalue)
print("Cells left out for fewer than 10 intervals:", bi[[f"{x}_a" for x in ("NREM", "WAKE", "REM")]].isna().sum().to_dict())
print("\nPer subject, median burst index % (reading a):")
print((100 * bi.groupby("subject")[["NREM_a", "WAKE_a", "REM_a"]].median()).round(2).to_string())

Reading (a): medians % NREM 6.35, WAKE 2.64, REM 3.65; n {'NREM_a': 796, 'WAKE_a': 796, 'REM_a': 755}; Kruskal-Wallis H = 162.3, p = 5.7e-36
   Friedman test on the 753 cells with all three values (paired): p = 3.9e-151
   within tolerance: [np.False_, np.False_, np.False_] | order NREM > REM > WAKE: True
   REM against WAKE, paired Wilcoxon p = 1e-32
Reading (b): medians % NREM 11.67, WAKE 4.95, REM 6.71; n {'NREM_b': 796, 'WAKE_b': 796, 'REM_b': 755}; Kruskal-Wallis H = 171.0, p = 7.3e-38
   Friedman test on the 753 cells with all three values (paired): p = 1e-151
   within tolerance: [np.True_, np.False_, np.False_] | order NREM > REM > WAKE: True
   REM against WAKE, paired Wilcoxon p = 1.2e-31
Cells left out for fewer than 10 intervals: {'NREM_a': 16, 'WAKE_a': 16, 'REM_a': 57}

Per subject, median burst index % (reading a):
           NREM_a  WAKE_a  REM_a
subject                         
BWRat17      4.72    2.25   2.77
BWRat18      7.32    2.37   1.86
BWRat19      8.01    2.48   4.73
BWRat20      2.60    1.11   1.76
BWRat21      3.99    1.48   2.04
Bogey        6.51    1.44   2.52
Dino         5.06    2.07   3.18
Rizzo       10.49    6.67   7.00
Splinter     3.99    2.48   2.16
Templeton   11.26    5.40   6.90

# The helper on a synthetic train with a spike every 10 ms: every pair must fall in the 10, 20, ...,
# 290 ms bins (lags are rounded to 1 ns before binning, so float error cannot move a pair across an edge).
h_syn, _ = ss.autocorrelogram(np.arange(0, 10, 0.01), [0.0], [10.0])
print("Synthetic 10 ms train, non-empty bins (ms):", np.nonzero(h_syn)[0].tolist())
assert np.nonzero(h_syn)[0].tolist() == list(range(10, 300, 10))
acg = {state: [] for state in ("NREM", "WAKE", "REM")}
acg_sess = {state: [] for state in acg}
for s in sessions:
    for u in s.pe:
        for state in acg:
            m = s.state == state
            h, n = ss.autocorrelogram(s.spikes[u], s.starts[m], s.stops[m])
            if n >= 100:
                acg[state].append(h / n)
                acg_sess[state].append(s.session_id)
lags = np.arange(300) + 0.5
res = {}
for state, hs in acg.items():
    mean = np.mean(hs, axis=0)
    early = np.argmax(mean[2:50]) + 2
    sm = np.convolve(mean, np.ones(21) / 21, mode="same")
    w = (lags >= 60) & (lags <= 250)
    k = np.where(w)[0][np.argmax(sm[w])]
    res[state] = (lags[early], lags[k], len(hs), k)
    print("%s: cells %d; peak 2-50 ms at %.1f ms; smoothed maximum 60-250 ms at %.1f ms (%s the window edges)" % (
        state, len(hs), lags[early], lags[k], "inside" if 63 <= lags[k] <= 247 else "at"))
for state, hs in acg.items():
    mean = np.mean(hs, axis=0)
    hp = np.convolve(mean, np.ones(21) / 21, mode="same") / np.convolve(mean, np.ones(101) / 101, mode="same")
    w = (lags >= 80) & (lags <= 200)
    print("Variant, %s: 21 ms over 101 ms moving average, maximum 80-200 ms at %.1f ms (ratio %.3f)" % (
        state, lags[w][np.argmax(hp[w])], hp[w].max()))
print("\nNREM peak within 7-9 ms:", 7 <= res["NREM"][0] <= 9,
      "| theta peak within 112.5-137.5 ms in REM:", 112.5 <= res["REM"][1] <= 137.5,
      "and WAKE:", 112.5 <= res["WAKE"][1] <= 137.5)
fig, ax = plt.subplots(1, 2, figsize=(9, 3))
for state, c in (("WAKE", "k"), ("NREM", "tab:blue"), ("REM", "tab:red")):
    mean = np.mean(acg[state], axis=0)
    ax[0].plot(lags[:50], mean[:50], color=c, label=state)
    ax[1].plot(lags[30:], np.convolve(mean, np.ones(21) / 21, mode="same")[30:] / mean[30:].mean(), color=c)
ax[0].set_xlabel("lag (ms)"); ax[0].set_ylabel("count per spike per ms"); ax[0].legend()
ax[1].set_xlabel("lag (ms)"); ax[1].set_ylabel("smoothed, relative to mean 30-300 ms")
plt.show()
for state in acg:
    mean = np.mean(acg[state], axis=0)
    print("%s mean ACG range 0-50 ms %.4f to %.4f; smoothed relative 30-300 ms %.3f to %.3f" % (
        state, mean[:50].min(), mean[:50].max(),
        *(lambda v: (v.min(), v.max()))((np.convolve(mean, np.ones(21) / 21, mode="same")[30:] / mean[30:].mean())[10:-10])))
print("\nPer subject, NREM early peak (ms) and theta-window maximum in REM and WAKE (ms):")
for subj in sorted(tab.subject.unique()):
    out = []
    for state in ("NREM", "REM", "WAKE"):
        hs = []
        for s in sessions:
            if s.subject != subj:
                continue
            for u in s.pe:
                m = s.state == state
                h, n = ss.autocorrelogram(s.spikes[u], s.starts[m], s.stops[m])
                if n >= 100:
                    hs.append(h / n)
        if not hs:
            out.append(None); continue
        mean = np.mean(hs, axis=0)
        sm = np.convolve(mean, np.ones(21) / 21, mode="same")
        w = (lags >= 60) & (lags <= 250)
        out.append(lags[np.argmax(mean[2:50]) + 2] if state == "NREM" else lags[np.where(w)[0][np.argmax(sm[w])]])
    print("  %s: NREM %s, REM %s, WAKE %s" % (subj, *out))

Synthetic 10 ms train, non-empty bins (ms): [10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250, 260, 270, 280, 290]

NREM: cells 787; peak 2-50 ms at 8.5 ms; smoothed maximum 60-250 ms at 60.5 ms (at the window edges)
WAKE: cells 781; peak 2-50 ms at 14.5 ms; smoothed maximum 60-250 ms at 61.5 ms (at the window edges)
REM: cells 685; peak 2-50 ms at 13.5 ms; smoothed maximum 60-250 ms at 60.5 ms (at the window edges)
Variant, NREM: 21 ms over 101 ms moving average, maximum 80-200 ms at 167.5 ms (ratio 0.995)
Variant, WAKE: 21 ms over 101 ms moving average, maximum 80-200 ms at 118.5 ms (ratio 1.016)
Variant, REM: 21 ms over 101 ms moving average, maximum 80-200 ms at 142.5 ms (ratio 1.017)

NREM peak within 7-9 ms: True | theta peak within 112.5-137.5 ms in REM: False and WAKE: False

NREM mean ACG range 0-50 ms 0.0001 to 0.0078; smoothed relative 30-300 ms 0.742 to 1.569
WAKE mean ACG range 0-50 ms 0.0000 to 0.0043; smoothed relative 30-300 ms 0.886 to 1.097
REM mean ACG range 0-50 ms 0.0001 to 0.0058; smoothed relative 30-300 ms 0.878 to 1.188

Per subject, NREM early peak (ms) and theta-window maximum in REM and WAKE (ms):
  BWRat17: NREM 7.5, REM 62.5, WAKE 60.5
  BWRat18: NREM 10.5, REM 84.5, WAKE 99.5

  BWRat19: NREM 4.5, REM 77.5, WAKE 88.5

  BWRat20: NREM 8.5, REM 60.5, WAKE 61.5
  BWRat21: NREM 11.5, REM 138.5, WAKE 140.5
  Bogey: NREM 8.5, REM 135.5, WAKE 138.5

  Dino: NREM 9.5, REM 125.5, WAKE 131.5

  Rizzo: NREM 5.5, REM 60.5, WAKE 60.5
  Splinter: NREM 8.5, REM 107.5, WAKE 60.5
  Templeton: NREM 6.5, REM 60.5, WAKE 60.5

# sextiles of the whole-session rate, reading (b)
whole = {s.path: np.array([len(s.spikes[u]) / (s.stops[-1] - s.starts[0]) for u in s.pe]) for s in sessions}
edges_b = np.quantile(np.concatenate(list(whole.values())), np.linspace(0, 1, 7))
groups_b = {p: np.searchsorted(edges_b[1:-1], r, side="right") + 1 for p, r in whole.items()}
pairs = []
n_pairs_all = 0
for s in sessions:
    m = s.state == "NREM"
    binned = []
    for a, b in zip(s.starts[m], s.stops[m]):
        e = np.arange(a, b, 0.1)
        if len(e) > 1:
            binned.append(np.array([np.diff(np.searchsorted(s.spikes[u], e)) for u in s.pe]))
    X = np.concatenate(binned, axis=1)
    with np.errstate(invalid="ignore", divide="ignore"):
        C = np.corrcoef(X)
    n_pairs_all += len(s.pe) * (len(s.pe) - 1) // 2
    ga, gb = sext[s.path][1], groups_b[s.path]
    iu = np.triu_indices(len(s.pe), 1)
    for i, j in zip(*iu):
        if np.isfinite(C[i, j]):
            pairs.append(dict(subject=s.subject, session=s.session_id, r=C[i, j],
                              a=tuple(sorted((ga[i], ga[j]))), b=tuple(sorted((gb[i], gb[j])))))
pairs = pd.DataFrame(pairs)
print("pairs with a defined correlation %d of %d (cells without Non-REM spikes have none)" % (len(pairs), n_pairs_all))
print("pairs:", len(pairs), "| overall mean r %.4f, fraction negative %.3f" % (pairs.r.mean(), np.mean(pairs.r < 0)))
for rd in ("a", "b"):
    g = pairs.groupby(rd).r.agg(["mean", "size"])
    mat = np.full((6, 6), np.nan)
    for (i, j), v in g["mean"].items():
        mat[i - 1, j - 1] = mat[j - 1, i - 1] = v
    print("\nReading (%s), mean r by group pair (rows and columns: groups 1 to 6):" % rd)
    print(pd.DataFrame(mat, index=range(1, 7), columns=range(1, 7)).round(4).to_string())
    r16 = pairs[pairs[rd] == (1, 6)].r
    t = st.ttest_1samp(r16, 0)
    top = g["mean"].idxmax()
    print("(6,6) mean %.4f, largest pair is %s; (1,6) mean %.4f, n %d, t = %.2f, p = %.2g" % (
        g["mean"][(6, 6)], top, r16.mean(), len(r16), t.statistic, t.pvalue))
print("\nPer subject, reading (a): mean r of (6,6) and (1,6) pairs")
print(pairs.groupby(["subject", "a"]).r.mean().unstack()[[(6, 6), (1, 6)]].round(4).to_string())

pairs with a defined correlation 20764 of 20933 (cells without Non-REM spikes have none)
pairs: 20764 | overall mean r 0.0156, fraction negative 0.228

Reading (a), mean r by group pair (rows and columns: groups 1 to 6):
        1       2       3       4       5       6
1  0.0061  0.0081  0.0076  0.0076  0.0073  0.0089
2  0.0081  0.0131  0.0119  0.0119  0.0122  0.0142
3  0.0076  0.0119  0.0188  0.0167  0.0145  0.0196
4  0.0076  0.0119  0.0167  0.0181  0.0185  0.0206
5  0.0073  0.0122  0.0145  0.0185  0.0213  0.0249
6  0.0089  0.0142  0.0196  0.0206  0.0249  0.0424
(6,6) mean 0.0424, largest pair is (np.int64(6), np.int64(6)); (1,6) mean 0.0089, n 867, t = 10.85, p = 8.7e-26

Reading (b), mean r by group pair (rows and columns: groups 1 to 6):
        1       2       3       4       5       6
1  0.0057  0.0073  0.0062  0.0080  0.0083  0.0087
2  0.0073  0.0125  0.0133  0.0118  0.0119  0.0148
3  0.0062  0.0133  0.0145  0.0155  0.0152  0.0186
4  0.0080  0.0118  0.0155  0.0256  0.0197  0.0220
5  0.0083  0.0119  0.0152  0.0197  0.0221  0.0258
6  0.0087  0.0148  0.0186  0.0220  0.0258  0.0401
(6,6) mean 0.0401, largest pair is (np.int64(6), np.int64(6)); (1,6) mean 0.0087, n 896, t = 13.45, p = 1e-37

Per subject, reading (a): mean r of (6,6) and (1,6) pairs
a          (6, 6)  (1, 6)
subject                  
BWRat17    0.0558  0.0067
BWRat18       NaN -0.0021
BWRat19    0.0551  0.0066
BWRat20    0.0663  0.0144
BWRat21       NaN     NaN
Bogey      0.0606  0.0073
Dino       0.0319  0.0074
Rizzo      0.0322  0.0105
Splinter   0.0818  0.0154
Templeton     NaN     NaN

rows = []
all_ma = []
ep_time = ep_ma_n = 0
for s in sessions:
    dur = s.stops - s.starts
    all_ma.append(dur[s.state == "MA"])
    periods = ss.sleep_periods(s.starts, s.stops, s.labels)
    n = T = tma = 0
    durs_in = []
    for p_ in periods:
        rr = np.arange(p_["first"], p_["last"] + 1)
        ma = rr[s.state[rr] == "MA"]
        n += len(ma); tma += dur[ma].sum(); T += p_["stop"] - p_["start"]; durs_in += list(dur[ma])
    for a, b in ss.nrem_episodes(s.state):
        ep_time += s.stops[b] - s.starts[a]
        ep_ma_n += int(np.sum(s.state[a:b + 1] == "MA"))
    rows.append(dict(subject=s.subject, session=s.session_id, periods=len(periods), MA_in_sleep=n,
                     sleep_min=T / 60, freq=n / (T / 60) if T else np.nan, frac=100 * tma / T if T else np.nan,
                     durs=durs_in))
ma = pd.DataFrame(rows)
in_durs = np.concatenate([np.array(x, float) for x in ma.durs])
all_ma = np.concatenate(all_ma)
print(ma.drop(columns="durs").round(3).to_string())
a = (np.median(in_durs), ma.freq.mean(), ma.frac.mean())
b = (np.median(all_ma), ep_ma_n / (ep_time / 60), 100 * all_ma.sum() / ep_time)
for name, v in (("(a) inside SLEEP, per-session means", a), ("(b) all MAs, pooled, per nonREM episode time", b)):
    ok = [abs(v[0] - 15.9) <= 3.3, abs(v[1] - 0.32) <= 0.12, abs(v[2] - 9.1) <= 3.5]
    print("%s: median duration %.1f s, frequency %.3f per min, fraction %.2f%%; within tolerance %s" % (name, *v, ok))
print("SD over sessions (a): frequency %.3f, fraction %.2f" % (ma.freq.std(), ma.frac.std()))
print("MAs: %d in all, %d inside SLEEP periods (paper: 790 in 27 sessions)" % (len(all_ma), len(in_durs)))
print("MA durations (s): min %.0f, max %.0f; the state table has 1 s resolution" % (all_ma.min(), all_ma.max()))
print("Awake epochs under 40 s not between two Non-REM epochs (left as WAKE):",
      sum(int(np.sum((s.labels == "Awake") & (s.stops - s.starts < 40) & (s.state == "WAKE"))) for s in sessions))

      subject           session  periods  MA_in_sleep  sleep_min   freq   frac
0     BWRat17    BWRat17_121712        1            3     26.300  0.114  4.119
1     BWRat17    BWRat17_121912        1           12     70.667  0.170  6.108
2     BWRat18    BWRat18_020513        1            1     28.100  0.036  0.593
3     BWRat19    BWRat19_032413        1            8     40.567  0.197  7.642
4     BWRat19    BWRat19_032513        1            4     24.483  0.163  4.697
5     BWRat20    BWRat20_101013        0            0      0.000    NaN    NaN
6     BWRat20    BWRat20_101513        5           35    256.233  0.137  3.181
7     BWRat21    BWRat21_121113        6           11    221.750  0.050  1.841
8     BWRat21    BWRat21_121613        3           10    149.933  0.067  2.034
9     BWRat21    BWRat21_121813        4            7    185.750  0.038  1.355
10      Bogey      Bogey_012615        4           38    256.417  0.148  6.181
11       Dino       Dino_061814        4           15    151.267  0.099  2.666
12       Dino       Dino_061914        2            6     60.367  0.099  4.307
13       Dino       Dino_062014        5           14    186.300  0.075  2.308
14       Dino       Dino_072314        3            4     88.883  0.045  1.669
15       Dino       Dino_072414        2            4     70.433  0.057  1.514
16      Rizzo      Rizzo_022615       11           25    372.267  0.067  2.306
17      Rizzo      Rizzo_022715        3           14    201.467  0.069  1.911
18   Splinter   Splinter_020515        5           12    211.767  0.057  1.850
19   Splinter   Splinter_020915        4           17    246.717  0.069  2.229
20  Templeton  Templeton_032415        3           22    131.367  0.167  5.950
(a) inside SLEEP, per-session means: median duration 19.0 s, frequency 0.096 per min, fraction 3.22%; within tolerance [np.True_, np.False_, np.False_]
(b) all MAs, pooled, per nonREM episode time: median duration 19.0 s, frequency 0.127 per min, fraction 4.23%; within tolerance [np.True_, np.False_, np.False_]
SD over sessions (a): frequency 0.051, fraction 1.97
MAs: 391 in all, 262 inside SLEEP periods (paper: 790 in 27 sessions)
MA durations (s): min 7, max 39; the state table has 1 s resolution
Awake epochs under 40 s not between two Non-REM epochs (left as WAKE): 345

# Awake epochs under 40 s by their neighbours in the state table, and MA counts under looser rules,
# counted inside SLEEP periods (per minute of SLEEP and percent of SLEEP time, per-session means as in
# reading (a)). The first rule is the one used throughout.
def neighbour(s, j):
    return s.labels[j] if 0 <= j < len(s.labels) else "edge"

short_n = collections.Counter()
short_d = collections.defaultdict(list)
for s in sessions:
    dur = s.stops - s.starts
    for i in np.where((s.labels == "Awake") & (dur < ss.MA_MAX))[0]:
        key = (neighbour(s, i - 1), neighbour(s, i + 1))
        short_n[key] += 1
        short_d[key].append(dur[i])
print("Awake epochs under 40 s by previous and next label: count, median duration")
for key, v in short_n.most_common():
    print("  %s, Awake, %s: %d, median %.0f s" % (key[0], key[1], v, np.median(short_d[key])))

both = lambda s, i: neighbour(s, i - 1) == "Non-REM" and neighbour(s, i + 1) == "Non-REM"
one = lambda s, i: "Non-REM" in (neighbour(s, i - 1), neighbour(s, i + 1))
rules = {"Non-REM on both sides, under 40 s (used here)": lambda s, i, d: d < 40 and both(s, i),
         "Non-REM on at least one side, under 40 s": lambda s, i, d: d < 40 and one(s, i),
         "any neighbours, under 40 s": lambda s, i, d: d < 40,
         "Non-REM on both sides, under 60 s": lambda s, i, d: d < 60 and both(s, i),
         "Non-REM on at least one side, under 60 s": lambda s, i, d: d < 60 and one(s, i),
         "every Awake epoch inside SLEEP (all are 120 s or shorter)": lambda s, i, d: True}
print("\nMAs inside SLEEP under each rule (paper: 0.32 per min, 9.1% of SLEEP):")
for name, rule in rules.items():
    freq, frac, n_all = [], [], 0
    for s in sessions:
        dur = s.stops - s.starts
        n = T = tma = 0
        for p_ in ss.sleep_periods(s.starts, s.stops, s.labels):
            for i in range(p_["first"], p_["last"] + 1):
                if s.labels[i] == "Awake" and rule(s, i, dur[i]):
                    n += 1
                    tma += dur[i]
            T += p_["stop"] - p_["start"]
        if T:
            freq.append(n / (T / 60)); frac.append(100 * tma / T); n_all += n
    print("  %s: %d, %.3f per min, %.2f%% of SLEEP" % (name, n_all, np.mean(freq), np.mean(frac)))

Awake epochs under 40 s by previous and next label: count, median duration
  Non-REM, Awake, Non-REM: 391, median 19 s
  REM, Awake, Non-REM: 186, median 23 s
  REM, Awake, REM: 120, median 15 s
  Non-REM, Awake, REM: 36, median 14 s
  edge, Awake, Non-REM: 2, median 14 s
  REM, Awake, edge: 1, median 22 s

MAs inside SLEEP under each rule (paper: 0.32 per min, 9.1% of SLEEP):
  Non-REM on both sides, under 40 s (used here): 262, 0.096 per min, 3.22% of SLEEP
  Non-REM on at least one side, under 40 s: 421, 0.145 per min, 5.11% of SLEEP
  any neighbours, under 40 s: 518, 0.175 per min, 5.97% of SLEEP
  Non-REM on both sides, under 60 s: 309, 0.110 per min, 4.36% of SLEEP
  Non-REM on at least one side, under 60 s: 524, 0.177 per min, 7.63% of SLEEP
  every Awake epoch inside SLEEP (all are 120 s or shorter): 732, 0.245 per min, 13.66% of SLEEP

fr, fr_sess = [], []
for s in sessions:
    for a, b in ss.nrem_episodes(s.state):
        rr = np.arange(a, b + 1)
        mas = [(s.starts[i], s.stops[i]) for i in rr if s.state[i] == "MA"]
        if not mas:
            continue
        bins = ss.binned_intervals(s.starts[a], s.stops[b], 10, mas)
        width = (s.stops[b] - s.starts[a]) / 10
        fr_sess.append(s.session_id)
        fr.append([np.sum(e - st_) / width for st_, e in bins])
fr = np.array(fr)
x = np.tile(np.arange(10) + 0.5, len(fr))
ra = st.pearsonr(np.arange(10) + 0.5, fr.mean(axis=0))
rb = st.pearsonr(x, fr.ravel())
print("Episodes with at least one MA:", len(fr))
print("Mean MA fraction per bin:", np.round(fr.mean(axis=0), 4))
print("(a) bin means: R = %.3f, p = %.2g | (b) pooled: R = %.3f, p = %.2g" % (ra[0], ra[1], rb[0], rb[1]))
print("Within -0.935 to -0.765:", {"a": -0.935 <= ra[0] <= -0.765, "b": -0.935 <= rb[0] <= -0.765})
# the first and last bins start and end in a packet by construction
print("Variant without the first and last bins: (a) R = %.3f, p = %.2g" % st.pearsonr(np.arange(1, 9) + 0.5, fr.mean(axis=0)[1:-1]))

Episodes with at least one MA: 243
Mean MA fraction per bin: [0.0674 0.1117 0.131  0.1506 0.1439 0.1566 0.1228 0.083  0.0651 0.01  ]
(a) bin means: R = -0.451, p = 0.19 | (b) pooled: R = -0.087, p = 1.7e-05
Within -0.935 to -0.765: {'a': np.False_, 'b': np.False_}
Variant without the first and last bins: (a) R = -0.571, p = 0.14

pa, pb = [], []
for s in sessions:
    dur = s.stops - s.starts
    pa.append(dur[s.state == "NREM"])
    for p_ in ss.sleep_periods(s.starts, s.stops, s.labels):
        rr = np.arange(p_["first"], p_["last"] + 1)
        pb.append(dur[rr][s.state[rr] == "NREM"])
pa, pb = np.concatenate(pa), np.concatenate(pb)
for name, v in (("(a) all", pa), ("(b) inside SLEEP", pb)):
    print("%s: n %d, median %.1f s, 5th to 95th pct %.0f to %.0f s, min %.0f, max %.0f; within 99.9-122.1: %s" % (
        name, len(v), np.median(v), *np.percentile(v, [5, 95]), v.min(), v.max(), 99.9 <= np.median(v) <= 122.1))
print("Per subject median (a):", {subj: float(np.median(np.concatenate([(s.stops - s.starts)[s.state == "NREM"]
      for s in sessions if s.subject == subj]))) for subj in sorted(tab.subject.unique())})

(a) all: n 1069, median 108.0 s, 5th to 95th pct 12 to 484 s, min 6, max 1199; within 99.9-122.1: True
(b) inside SLEEP: n 685, median 128.0 s, 5th to 95th pct 16 to 524 s, min 6, max 1199; within 99.9-122.1: False
Per subject median (a): {np.str_('BWRat17'): 90.0, np.str_('BWRat18'): 219.0, np.str_('BWRat19'): 78.5, np.str_('BWRat20'): 83.0, np.str_('BWRat21'): 123.0, np.str_('Bogey'): 99.0, np.str_('Dino'): 115.0, np.str_('Rizzo'): 113.0, np.str_('Splinter'): 176.5, np.str_('Templeton'): 76.0}

cyc = {s.path: ss.cycles(s) for s in sessions}
cyc_open = {s.path: ss.cycles(s, require_wake_end=False) for s in sessions}
nc = pd.Series({s.session_id: len(cyc[s.path]) for s in sessions})
nco = pd.Series({s.session_id: len(cyc_open[s.path]) for s in sessions})
print(pd.DataFrame({"cycles": nc, "cycles incl. SLEEP ended by the recording": nco}).to_string())
print("Total %d, mean per session %.2f, sessions with none: %d | variant: total %d, mean %.2f, none: %d" % (
    nc.sum(), nc.mean(), (nc == 0).sum(), nco.sum(), nco.mean(), (nco == 0).sum()))
print("Within 1.8-2.2:", 1.8 <= nc.mean() <= 2.2)
pw = [s.stops[p_["pre"]] - s.starts[p_["pre"]] for s in sessions for p_ in cyc[s.path]]
sl = [(p_["stop"] - p_["start"]) / 60 for s in sessions for p_ in cyc[s.path]]
print("Pre-sleep WAKE (min): %.1f to %.1f; SLEEP duration (min): median %.1f, range %.1f to %.1f" % (
    min(pw) / 60, max(pw) / 60, np.median(sl), min(sl), max(sl)))
print("\nSessions without a cycle: each sleep segment between Awake epochs over 120 s (leading and trailing Awake"
      " trimmed), the Awake before it, and why it is not a cycle; segments under 1 min not listed")
for s in sessions:
    if cyc[s.path]:
        continue
    dur = s.stops - s.starts
    long_wake = np.where((s.labels == "Awake") & (dur > ss.SLEEP_BREAK))[0]
    cuts = [-1] + list(long_wake) + [len(s.labels)]
    print("  %s: %.0f min of states, %d Awake epochs over 120 s" % (s.session_id, (s.stops[-1] - s.starts[0]) / 60, len(long_wake)))
    for lo_, hi_ in zip(cuts[:-1], cuts[1:]):
        rr = [r for r in range(lo_ + 1, hi_) if s.labels[r] != "Awake"]
        if not rr or s.stops[rr[-1]] - s.starts[rr[0]] < 60:
            continue
        seg = (s.stops[rr[-1]] - s.starts[rr[0]]) / 60
        pre_ = dur[lo_] / 60 if lo_ >= 0 else None
        why = ([] if seg >= 20 else ["shorter than 20 min"]) + ([] if hi_ < len(s.labels) else ["ended by the recording"]) + \
              (["no Awake before"] if pre_ is None else [] if pre_ >= 7 else ["Awake before under 7 min"])
        print("    %.1f to %.1f min: %.1f min, Awake before %s; %s" % (
            s.starts[rr[0]] / 60, s.stops[rr[-1]] / 60, seg, "none" if pre_ is None else "%.1f min" % pre_,
            ", ".join(why) or "a cycle"))

                  cycles  cycles incl. SLEEP ended by the recording
BWRat17_121712         1                                          1
BWRat17_121912         1                                          1
BWRat18_020513         0                                          1
BWRat19_032413         0                                          0
BWRat19_032513         1                                          1
BWRat20_101013         0                                          1
BWRat20_101513         1                                          1
BWRat21_121113         2                                          2
BWRat21_121613         1                                          1
BWRat21_121813         2                                          3
Bogey_012615           2                                          2
Dino_061814            1                                          1
Dino_061914            1                                          1
Dino_062014            2                                          2
Dino_072314            2                                          2
Dino_072414            2                                          3
Rizzo_022615           4                                          4
Rizzo_022715           1                                          1
Splinter_020515        3                                          3
Splinter_020915        2                                          2
Templeton_032415       1                                          1
Total 30, mean per session 1.43, sessions with none: 3 | variant: total 34, mean 1.62, none: 1
Within 1.8-2.2: False
Pre-sleep WAKE (min): 7.5 to 55.4; SLEEP duration (min): median 35.2, range 20.7 to 125.6

Sessions without a cycle: each sleep segment between Awake epochs over 120 s (leading and trailing Awake trimmed), the Awake before it, and why it is not a cycle; segments under 1 min not listed
  BWRat18_020513: 126 min of states, 3 Awake epochs over 120 s
    21.4 to 49.5 min: 28.1 min, Awake before 3.6 min; Awake before under 7 min
    60.6 to 126.4 min: 65.8 min, Awake before 11.1 min; ended by the recording
  BWRat19_032413: 83 min of states, 3 Awake epochs over 120 s
    23.8 to 28.4 min: 4.6 min, Awake before 23.7 min; shorter than 20 min
    33.5 to 74.0 min: 40.6 min, Awake before 5.0 min; Awake before under 7 min
  BWRat20_101013: 223 min of states, 6 Awake epochs over 120 s
    111.0 to 122.2 min: 11.2 min, Awake before 31.8 min; shorter than 20 min
    129.7 to 140.3 min: 10.6 min, Awake before 7.5 min; shorter than 20 min
    151.8 to 221.3 min: 69.5 min, Awake before 11.4 min; ended by the recording

def group_mean(rmat, units_mask):
    """Mean over the cells in the mask, cells with no spikes over the whole interval dropped,
    divided by its mean over the points; None when undefined."""
    sub = rmat[units_mask]
    sub = sub[np.nansum(sub, axis=1) > 0]
    if len(sub) == 0:
        return None
    m = np.nanmean(sub, axis=0)
    mu = np.nanmean(m)
    return m / mu if mu > 0 else None

def sleep_points(nbins=10, min_nrem=10.0):
    """Per cycle and reading: x, rate matrix (pE x points), session path."""
    out = {"a": [], "b": []}
    for s in sessions:
        for p_ in cyc[s.path]:
            rr = np.arange(p_["first"], p_["last"] + 1)
            pk = rr[s.state[rr] == "NREM"]
            T = p_["stop"] - p_["start"]
            xa = ((s.starts[pk] + s.stops[pk]) / 2 - p_["start"]) / T
            ra = np.array([[ss.counts(s.spikes[u], s.starts[i], s.stops[i]) / (s.stops[i] - s.starts[i]) for i in pk] for u in s.pe])
            out["a"].append((xa, ra, s))
            bins = ss.binned_intervals(p_["start"], p_["stop"], nbins, list(zip(s.starts[pk], s.stops[pk])))
            rb, dur = ss.bin_rates(s, s.pe, bins)
            ok = dur >= min_nrem
            out["b"].append(((np.arange(nbins) + 0.5)[ok] / nbins, rb[:, ok], s))
    return out

pts = sleep_points()

def pooled_r(pts_r, groups=None):
    """Pooled Pearson R of normalized group means against x. groups None = all pE cells."""
    xs, ys = [], []
    for x, rmat, s in pts_r:
        mask = np.ones(len(s.pe), bool) if groups is None else np.isin(sext[s.path][1], groups)
        v = group_mean(rmat, mask)
        if v is not None:
            xs.append(x); ys.append(v)
    x, y = np.concatenate(xs), np.concatenate(ys)
    ok = np.isfinite(y)
    r, p = st.pearsonr(x[ok], y[ok])
    return r, p, ok.sum()

for rd in ("a", "b"):
    r, p, n = pooled_r(pts[rd])
    print("Reading (%s): R = %.3f, p = %.2g, points %d; within -0.121 to -0.099: %s" % (rd, r, p, n, -0.121 <= r <= -0.099))
print("\nPer subject, reading (a):")
for subj in sorted(tab.subject.unique()):
    sub = [t_ for t_ in pts["a"] if t_[2].subject == subj]
    if sub:
        r, p, n = pooled_r(sub)
        print("  %s: R = %.3f, p = %.2g, points %d, cycles %d" % (subj, r, p, n, len(sub)))

Reading (a): R = -0.267, p = 4.3e-06, points 287; within -0.121 to -0.099: False
Reading (b): R = -0.285, p = 7.6e-07, points 291; within -0.121 to -0.099: False

Per subject, reading (a):
  BWRat17: R = -0.712, p = 7.1e-06, points 31, cycles 2
  BWRat19: R = -0.811, p = 0.05, points 6, cycles 1
  BWRat20: R = -0.381, p = 0.11, points 19, cycles 1
  BWRat21: R = -0.155, p = 0.31, points 44, cycles 5
  Bogey: R = -0.340, p = 0.037, points 38, cycles 2
  Dino: R = -0.210, p = 0.14, points 52, cycles 8
  Rizzo: R = -0.321, p = 0.06, points 35, cycles 5
  Splinter: R = -0.132, p = 0.4, points 43, cycles 5
  Templeton: R = 0.031, p = 0.9, points 19, cycles 1

fl = []
for s in sessions:
    for c, p_ in enumerate(cyc[s.path]):
        rr = np.arange(p_["first"], p_["last"] + 1)
        pk = rr[s.state[rr] == "NREM"]
        r1 = ss.rates(s.spikes, s.pe, s.starts[pk[0]], s.stops[pk[0]])
        r2 = ss.rates(s.spikes, s.pe, s.starts[pk[-1]], s.stops[pk[-1]])
        for k, u in enumerate(s.pe):
            fl.append(dict(subject=s.subject, session=s.session_id, cycle=c, unit=u, first=r1[k], last=r2[k],
                           sextile=sext[s.path][1][k], first_dur=s.stops[pk[0]] - s.starts[pk[0]],
                           last_dur=s.stops[pk[-1]] - s.starts[pk[-1]],
                           first_half=ss.rates(s.spikes, [u], s.starts[pk[0]], (s.starts[pk[0]] + s.stops[pk[0]]) / 2)[0],
                           second_half=ss.rates(s.spikes, [u], (s.starts[pk[0]] + s.stops[pk[0]]) / 2, s.stops[pk[0]])[0]))
fl = pd.DataFrame(fl)
w = st.wilcoxon(fl["first"], fl["last"])
print("cell-cycle pairs %d (cycles %d); first %.4f +- %.3f Hz, last %.4f +- %.3f Hz; Wilcoxon p = %.2g; last lower: %s" % (
    len(fl), fl.groupby(["session", "cycle"]).ngroups, fl["first"].mean(), fl["first"].std(), fl["last"].mean(),
    fl["last"].std(), w.pvalue, fl["last"].mean() < fl["first"].mean()))
print("Within tolerance: first", abs(fl["first"].mean() - 1.04) <= 0.104, "last", abs(fl["last"].mean() - 0.88) <= 0.088)
print("Medians: first %.3f, last %.3f; first packet duration median %.0f s, last %.0f s" % (
    fl["first"].median(), fl["last"].median(), fl.groupby(["session", "cycle"]).first_dur.first().median(),
    fl.groupby(["session", "cycle"]).last_dur.first().median()))
print("\nPer subject:")
for subj, g in fl.groupby("subject"):
    print("  %s: first %.3f, last %.3f, n %d, Wilcoxon p = %.2g" % (subj, g["first"].mean(), g["last"].mean(), len(g),
                                                                    st.wilcoxon(g["first"], g["last"]).pvalue))
pcyc = fl.groupby(["session", "cycle"])[["first", "last"]].mean()
print("\nOne value per cycle (mean over its pE cells): first %.3f, last %.3f Hz; last lower in %d of %d cycles, "
      "Wilcoxon over cycles p = %.2g" % (pcyc["first"].mean(), pcyc["last"].mean(), (pcyc["last"] < pcyc["first"]).sum(),
                                        len(pcyc), st.wilcoxon(pcyc["first"], pcyc["last"]).pvalue))

cell-cycle pairs 1070 (cycles 30); first 1.0915 +- 1.483 Hz, last 0.9679 +- 1.241 Hz; Wilcoxon p = 3.9e-06; last lower: True
Within tolerance: first True last True
Medians: first 0.692, last 0.657; first packet duration median 114 s, last 146 s

Per subject:
  BWRat17: first 0.687, last 0.505, n 73, Wilcoxon p = 0.025
  BWRat19: first 2.172, last 1.470, n 50, Wilcoxon p = 0.0094
  BWRat20: first 1.098, last 0.984, n 58, Wilcoxon p = 0.11
  BWRat21: first 0.689, last 0.641, n 98, Wilcoxon p = 0.35
  Bogey: first 0.877, last 0.794, n 64, Wilcoxon p = 0.8
  Dino: first 1.157, last 1.066, n 366, Wilcoxon p = 2.3e-05
  Rizzo: first 1.488, last 1.320, n 189, Wilcoxon p = 0.63
  Splinter: first 0.690, last 0.663, n 162, Wilcoxon p = 0.28
  Templeton: first 0.540, last 0.741, n 10, Wilcoxon p = 0.084

One value per cycle (mean over its pE cells): first 1.040, last 0.949 Hz; last lower in 22 of 30 cycles, Wilcoxon over cycles p = 0.015

k2 = (fl["first"] > 0) & (fl["last"] > 0)
fits = {"(a) log10, last on first": ols(np.log10(fl["first"][k2]), np.log10(fl["last"][k2])),
        "(b) Hz, last on first": ols(fl["first"], fl["last"]),
        "(a) reverse, first on last": ols(np.log10(fl["last"][k2]), np.log10(fl["first"][k2]))}
k3 = (fl.first_half > 0) & (fl.second_half > 0)
fits["check: second half on first half of the first packet, log10"] = ols(np.log10(fl.first_half[k3]), np.log10(fl.second_half[k3]))
for k, (b, lo, hi, r, p, n) in fits.items():
    print("%s: slope %.3f (95%% CI %.3f to %.3f), R = %.3f, n = %d" % (k, b, lo, hi, r, n))
print("Pairs dropped for a zero rate (a): %d" % (~k2).sum())
for k in ("(a) log10, last on first", "(b) Hz, last on first"):
    b, lo, hi = fits[k][:3]
    print(k, "| slope within 0.744-0.968:", 0.744 <= b <= 0.968, "| CI excludes 1:", hi < 1)
r_h = fits["check: second half on first half of the first packet, log10"][3]
rel = 2 * r_h / (1 + r_h)
print("Reliability of the log10 first-packet rate (Spearman-Brown, from its halves) %.3f; last-on-first slope corrected "
      "for attenuation %.3f (%.3f to %.3f, the reliability's own uncertainty not counted)" % (
          rel, *(v / rel for v in fits["(a) log10, last on first"][:3])))
# Spread of log10 rates in the last packet against the first: if fast cells slow down and slow cells
# speed up beyond regression to the mean, the SD is smaller in the last packet.
v = sd_ratio(np.log10(fl["first"][k2]), np.log10(fl["last"][k2]), (fl.session + "/" + fl.cycle.astype(str))[k2], seed=52)
print("SD of log10 rates, last over first packet: %.3f (Pitman-Morgan p = %.2g; 95%% interval resampling cycles %.3f to %.3f)" % v)
print("SD of log10 rates: first packet %.3f, last %.3f" % (np.log10(fl["first"][k2]).std(), np.log10(fl["last"][k2]).std()))

(a) log10, last on first: slope 0.801 (95% CI 0.768 to 0.834), R = 0.832, n = 1021
(b) Hz, last on first: slope 0.723 (95% CI 0.697 to 0.748), R = 0.863, n = 1070
(a) reverse, first on last: slope 0.864 (95% CI 0.829 to 0.900), R = 0.832, n = 1021
check: second half on first half of the first packet, log10: slope 0.855 (95% CI 0.829 to 0.881), R = 0.898, n = 987
Pairs dropped for a zero rate (a): 49
(a) log10, last on first | slope within 0.744-0.968: True | CI excludes 1: True
(b) Hz, last on first | slope within 0.744-0.968: False | CI excludes 1: True
Reliability of the log10 first-packet rate (Spearman-Brown, from its halves) 0.946; last-on-first slope corrected for attenuation 0.846 (0.812 to 0.881, the reliability's own uncertainty not counted)

SD of log10 rates, last over first packet: 0.963 (Pitman-Morgan p = 0.028; 95% interval resampling cycles 0.892 to 1.029)
SD of log10 rates: first packet 0.532, last 0.512

ba = []
for s in sessions:
    for c, p_ in enumerate(cyc[s.path]):
        pre, post = p_["pre"], p_["post"]
        r1 = ss.rates(s.spikes, s.pe, s.stops[pre] - 300, s.stops[pre])
        r2 = ss.rates(s.spikes, s.pe, s.starts[post], min(s.starts[post] + 300, s.stops[post]))
        for k, u in enumerate(s.pe):
            ba.append(dict(subject=s.subject, session=s.session_id, cycle=c, unit=u, pre=r1[k], post=r2[k],
                           sextile=sext[s.path][1][k], post_len=min(300, s.stops[post] - s.starts[post])))
ba = pd.DataFrame(ba)
print("Post-sleep WAKE windows shorter than 300 s: %d of %d cycles" % (
    (ba.groupby(["session", "cycle"]).post_len.first() < 300).sum(), ba.groupby(["session", "cycle"]).ngroups))
for rd, d in (("(a) cell-cycle pairs", ba),
              ("(b) per-cell means", ba.groupby(["session", "unit"]).agg(pre=("pre", "mean"), post=("post", "mean"),
                                                                         sextile=("sextile", "first")).reset_index())):
    w = st.wilcoxon(d.pre, d.post)
    print("\n%s: n %d, mean pre %.3f, post %.3f Hz, median change %.4f, Wilcoxon p = %.2g" % (
        rd, len(d), d.pre.mean(), d.post.mean(), (d.post - d.pre).median(), w.pvalue))
    for g in range(1, 7):
        x = d[d.sextile == g]
        alt = "greater" if g == 1 else "less"
        pv = st.wilcoxon(x.post, x.pre, alternative=alt).pvalue
        print("  sextile %d: n %d, mean pre %.3f, post %.3f, mean change %+.3f, one-tailed (%s) p = %.2g" % (
            g, len(x), x.pre.mean(), x.post.mean(), (x.post - x.pre).mean(), "increase" if g == 1 else "decrease", pv))
print("\nPer subject, (a): mean change post - pre (Hz) for all pE, sextile 1 and sextile 6")
print(ba.assign(ch=ba.post - ba.pre).groupby("subject").apply(
    lambda g: pd.Series({"all": g.ch.mean(), "s1": g.ch[g.sextile == 1].mean(), "s6": g.ch[g.sextile == 6].mean(), "n": len(g)}),
    include_groups=False).round(3).to_string())
# Variant with windows of equal length: the pre-sleep window is cut to the length of the post-sleep one.
ba2 = []
for s in sessions:
    for c, p_ in enumerate(cyc[s.path]):
        pre, post = p_["pre"], p_["post"]
        L = min(300, s.stops[post] - s.starts[post])
        r1 = ss.rates(s.spikes, s.pe, s.stops[pre] - L, s.stops[pre])
        r2 = ss.rates(s.spikes, s.pe, s.starts[post], s.starts[post] + L)
        for k, u in enumerate(s.pe):
            ba2.append(dict(session=s.session_id, cycle=c, unit=u, pre=r1[k], post=r2[k], sextile=sext[s.path][1][k], L=L))
ba2 = pd.DataFrame(ba2)

def before_after(d, label):
    out = []
    for g in range(1, 7):
        x = d[d.sextile == g]
        out.append("s%d %+.3f (p %.2g)" % (g, (x.post - x.pre).mean(),
                                           st.wilcoxon(x.post, x.pre, alternative="greater" if g == 1 else "less").pvalue))
    print("%s: n %d, cycles %d, mean pre %.3f, post %.3f Hz, p = %.2g\n   %s" % (
        label, len(d), d.groupby(["session", "cycle"]).ngroups, d.pre.mean(), d.post.mean(), st.wilcoxon(d.pre, d.post).pvalue,
        ", ".join(out)))

print("\nVariants (mean change per sextile in Hz, one-tailed p):")
before_after(ba2, "pre window cut to the length of the post window")
print("   window lengths (s): median %.0f, range %.0f to %.0f" % (
    ba2.groupby(["session", "cycle"]).L.first().median(), ba2.L.min(), ba2.L.max()))
before_after(ba[ba.post_len >= 300], "only cycles with 300 s of WAKE after SLEEP")
print("\nOne value per cycle (mean change over its cells), and per session; Wilcoxon over these (two-sided):")
chg = ba.assign(ch=ba.post - ba.pre)
for label, g in [("all pE", None)] + [("sextile %d" % g, g) for g in range(1, 7)]:
    x = chg if g is None else chg[chg.sextile == g]
    cm, sm = x.groupby(["session", "cycle"]).ch.mean(), x.groupby("session").ch.mean()
    print("  %s: per cycle mean %+.3f, lower in %d of %d, p = %.2g | per session lower in %d of %d, p = %.2g" % (
        label, cm.mean(), (cm < 0).sum(), len(cm), st.wilcoxon(cm).pvalue, (sm < 0).sum(), len(sm), st.wilcoxon(sm).pvalue))

Post-sleep WAKE windows shorter than 300 s: 23 of 30 cycles

(a) cell-cycle pairs: n 1070, mean pre 1.341, post 1.180 Hz, median change -0.0394, Wilcoxon p = 4.1e-09
  sextile 1: n 162, mean pre 0.088, post 0.086, mean change -0.002, one-tailed (increase) p = 0.3
  sextile 2: n 183, mean pre 0.303, post 0.292, mean change -0.011, one-tailed (decrease) p = 0.0043
  sextile 3: n 183, mean pre 0.608, post 0.537, mean change -0.071, one-tailed (decrease) p = 0.018
  sextile 4: n 173, mean pre 1.007, post 0.888, mean change -0.120, one-tailed (decrease) p = 0.03
  sextile 5: n 189, mean pre 1.633, post 1.444, mean change -0.189, one-tailed (decrease) p = 0.00063
  sextile 6: n 180, mean pre 4.284, post 3.727, mean change -0.557, one-tailed (decrease) p = 2.4e-05

(b) per-cell means: n 727, mean pre 1.318, post 1.131 Hz, median change -0.0463, Wilcoxon p = 1.4e-09
  sextile 1: n 119, mean pre 0.081, post 0.080, mean change -0.001, one-tailed (increase) p = 0.36
  sextile 2: n 123, mean pre 0.301, post 0.282, mean change -0.019, one-tailed (decrease) p = 0.0094
  sextile 3: n 125, mean pre 0.617, post 0.537, mean change -0.080, one-tailed (decrease) p = 0.0087
  sextile 4: n 120, mean pre 1.021, post 0.855, mean change -0.167, one-tailed (decrease) p = 0.0051
  sextile 5: n 128, mean pre 1.626, post 1.463, mean change -0.163, one-tailed (decrease) p = 0.0031
  sextile 6: n 112, mean pre 4.496, post 3.761, mean change -0.735, one-tailed (decrease) p = 2.3e-05

Per subject, (a): mean change post - pre (Hz) for all pE, sextile 1 and sextile 6
             all     s1     s6      n
subject                              
BWRat17   -0.223  0.029 -1.296   73.0
BWRat19   -0.857  0.003 -1.896   50.0
BWRat20    0.003  0.001 -0.462   58.0
BWRat21   -0.044  0.018    NaN   98.0
Bogey     -0.092  0.003 -0.918   64.0
Dino      -0.142 -0.021 -0.549  366.0
Rizzo     -0.198 -0.010 -0.379  189.0
Splinter  -0.072 -0.009  0.245  162.0
Templeton -0.150    NaN    NaN   10.0

Variants (mean change per sextile in Hz, one-tailed p):

pre window cut to the length of the post window: n 1070, cycles 30, mean pre 1.340, post 1.180 Hz, p = 1.7e-09
   s1 +0.000 (p 0.18), s2 -0.017 (p 0.014), s3 -0.075 (p 0.014), s4 -0.103 (p 0.036), s5 -0.199 (p 0.00063), s6 -0.545 (p 1.8e-05)
   window lengths (s): median 194, range 122 to 300
only cycles with 300 s of WAKE after SLEEP: n 254, cycles 7, mean pre 1.505, post 1.428 Hz, p = 0.055
   s1 +0.033 (p 0.081), s2 -0.004 (p 0.099), s3 -0.021 (p 0.51), s4 -0.026 (p 0.42), s5 -0.306 (p 0.0034), s6 -0.099 (p 0.27)

One value per cycle (mean change over its cells), and per session; Wilcoxon over these (two-sided):
  all pE: per cycle mean -0.126, lower in 20 of 30, p = 0.016 | per session lower in 14 of 18, p = 0.0034
  sextile 1: per cycle mean +0.002, lower in 11 of 29, p = 0.47 | per session lower in 6 of 17, p = 0.71
  sextile 2: per cycle mean +0.014, lower in 15 of 29, p = 0.75 | per session lower in 9 of 17, p = 0.31
  sextile 3: per cycle mean -0.071, lower in 19 of 30, p = 0.23 | per session lower in 11 of 18, p = 0.26
  sextile 4: per cycle mean +0.056, lower in 16 of 30, p = 0.87 | per session lower in 11 of 18, p = 0.39
  sextile 5: per cycle mean -0.157, lower in 19 of 28, p = 0.12 | per session lower in 12 of 17, p = 0.35
  sextile 6: per cycle mean -0.491, lower in 17 of 23, p = 0.048 | per session lower in 12 of 13, p = 0.01

for rd in ("a", "b"):
    print("Reading (%s):" % rd)
    for g in range(1, 7):
        r, p, n = pooled_r(pts[rd], [g])
        print("  sextile %d: R = %+.3f, p = %.2g, points %d" % (g, r, p, n))
print("\nPer subject, reading (a), sextiles 1 and 6:")
for subj in sorted(tab.subject.unique()):
    sub = [t_ for t_ in pts["a"] if t_[2].subject == subj]
    if sub:
        out = []
        for g in (1, 6):
            try:
                out.append("%+.3f (p %.2g)" % pooled_r(sub, [g])[:2])
            except ValueError:
                out.append("n/a")
        print("  %s: sextile 1 %s, sextile 6 %s" % (subj, *out))

Reading (a):
  sextile 1: R = +0.163, p = 0.0074, points 268
  sextile 2: R = +0.143, p = 0.017, points 281
  sextile 3: R = -0.052, p = 0.38, points 287
  sextile 4: R = -0.040, p = 0.5, points 287
  sextile 5: R = -0.341, p = 1.2e-08, points 265
  sextile 6: R = -0.447, p = 3.1e-11, points 200
Reading (b):
  sextile 1: R = +0.094, p = 0.11, points 281
  sextile 2: R = +0.070, p = 0.24, points 281
  sextile 3: R = -0.121, p = 0.04, points 291
  sextile 4: R = -0.004, p = 0.95, points 291
  sextile 5: R = -0.258, p = 1.7e-05, points 271
  sextile 6: R = -0.430, p = 2e-11, points 223

Per subject, reading (a), sextiles 1 and 6:
  BWRat17: sextile 1 +0.391 (p 0.03), sextile 6 -0.812 (p 0.026)
  BWRat19: sextile 1 +0.799 (p 0.057), sextile 6 -0.852 (p 0.031)
  BWRat20: sextile 1 +0.046 (p 0.85), sextile 6 -0.313 (p 0.19)
  BWRat21: sextile 1 +0.181 (p 0.24), sextile 6 n/a
  Bogey: sextile 1 +0.195 (p 0.24), sextile 6 -0.656 (p 7.8e-06)
  Dino: sextile 1 +0.140 (p 0.32), sextile 6 -0.209 (p 0.14)
  Rizzo: sextile 1 -0.058 (p 0.74), sextile 6 -0.257 (p 0.14)
  Splinter: sextile 1 +0.218 (p 0.16), sextile 6 -0.477 (p 0.0012)
  Templeton: sextile 1 n/a, sextile 6 n/a

def cv_sleep(pts_r):
    xs, ys = [], []
    for x, rmat, s in pts_r:
        mu = rmat.mean(axis=0); sd = rmat.std(axis=0)
        ok = mu > 0
        xs.append(x[ok]); ys.append(sd[ok] / mu[ok])
    return tuple(st.pearsonr(np.concatenate(xs), np.concatenate(ys))) + (len(np.concatenate(xs)),)

for rd in ("a", "b"):
    print("Reading (%s): R = %.3f, p = %.2g, points %d" % (rd, *cv_sleep(pts[rd])))

Reading (a): R = -0.026, p = 0.66, points 287
Reading (b): R = -0.014, p = 0.81, points 291

def intervals_of(kind, min_dur=20.0, sess=None):
    """(session, span start, span stop, included intervals) for packets, episodes, WAKE or REM epochs,
    over the sessions in `sess` (all by default)."""
    out = []
    for s in (sessions if sess is None else sess):
        if kind in ("packet", "WAKE", "REM"):
            lab = "NREM" if kind == "packet" else kind
            for i in np.where(s.state == lab)[0]:
                if s.stops[i] - s.starts[i] >= min_dur:
                    out.append((s, s.starts[i], s.stops[i], None, [i]))
        else:
            for a, b in ss.nrem_episodes(s.state):
                if s.stops[b] - s.starts[a] >= min_dur:
                    rows_ = [i for i in range(a, b + 1) if s.state[i] == "NREM"]
                    out.append((s, s.starts[a], s.stops[b], [(s.starts[i], s.stops[i]) for i in rows_], rows_))
    return out

def course(kind, nbins=10, min_dur=20.0, sess=None):
    """Per interval: bin centres, pE rate matrix, burst counts per bin (short, all), session."""
    res = []
    for s, a, b, inc, rows_ in intervals_of(kind, min_dur, sess):
        bins = ss.binned_intervals(a, b, nbins, inc)
        rmat, dur = ss.bin_rates(s, s.pe, bins)
        ok = dur > 0
        # burst: intervals within the same packet (or epoch), later spike binned
        pk_s = np.array([s.starts[i] for i in rows_]); pk_e = np.array([s.stops[i] for i in rows_])
        short = np.zeros(nbins); alln = np.zeros(nbins)
        for u in s.pe:
            isi, t2, _ = ss.isis_in(s.spikes[u], pk_s, pk_e)
            k = np.clip(((t2 - a) / (b - a) * nbins).astype(int), 0, nbins - 1)
            alln += np.bincount(k, minlength=nbins)
            short += np.bincount(k[isi < ss.BURST_ISI], minlength=nbins)
        x = (np.arange(nbins) + 0.5) / nbins
        res.append((x[ok], rmat[:, ok], s, short[ok], alln[ok]))
    return res

_cache = {}
def course_c(kind, nbins=10, min_dur=20.0):
    if (kind, nbins, min_dur) not in _cache:
        _cache[(kind, nbins, min_dur)] = course(kind, nbins, min_dur)
    return _cache[(kind, nbins, min_dur)]

courses = {k: course_c(k) for k in ("episode", "packet")}
print({k: len(v) for k, v in courses.items()}, "intervals")
print("Packets of at least 20 s by subject:", collections.Counter(c[2].subject for c in courses["packet"]))

{'episode': 622, 'packet': 952} intervals
Packets of at least 20 s by subject: Counter({np.str_('Dino'): 185, np.str_('BWRat21'): 167, np.str_('Rizzo'): 159, np.str_('BWRat20'): 108, np.str_('Splinter'): 92, np.str_('Bogey'): 86, np.str_('Templeton'): 72, np.str_('BWRat17'): 37, np.str_('BWRat19'): 27, np.str_('BWRat18'): 19})

def sextile_table(cs):
    out = {}
    for g in range(1, 7):
        out[g] = pooled_r([(x, r, s) for x, r, s, *_ in cs], [g])
    out["all pE"] = pooled_r([(x, r, s) for x, r, s, *_ in cs], None)
    return out

for kind in ("episode", "packet"):
    t = sextile_table(courses[kind])
    print("%s: " % kind + "; ".join("%s R = %+.3f (p %.2g)" % (("sextile %d" % g) if g != "all pE" else g, v[0], v[1])
                                     for g, v in t.items()))
print("\nVariants (sextile 1 and sextile 6 R, p):")
for kind in ("episode", "packet"):
    for nb, md in ((5, 20.0), (20, 20.0), (10, 60.0)):
        t = sextile_table(course_c(kind, nb, md))
        print("  %s, %d bins, min %d s: s1 %+.3f (p %.2g), s6 %+.3f (p %.2g)" % (kind, nb, md, *t[1][:2], *t[6][:2]))
print("\nPer subject, packets, sextile 1 and 6:")
for subj in sorted(tab.subject.unique()):
    cs = [c for c in courses["packet"] if c[2].subject == subj]
    out = []
    for g in (1, 6):
        try:
            out.append("%+.3f (p %.2g, n %d)" % pooled_r([(x, r, s) for x, r, s, *_ in cs], [g]))
        except ValueError:
            out.append("n/a")
    print("  %s: s1 %s, s6 %s" % (subj, *out))
fig, ax = plt.subplots(1, 2, figsize=(9, 3), sharey=True)
for j, kind in enumerate(("episode", "packet")):
    for g, c in zip(range(1, 7), plt.cm.viridis(np.linspace(0, 1, 6))):
        vals = []
        for x, r, s, *_ in courses[kind]:
            if len(x) == 10:
                v = group_mean(r, sext[s.path][1] == g)
                if v is not None:
                    vals.append(v)
        m = np.nanmean(vals, axis=0)
        ax[j].plot(np.arange(10) + 0.5, m, color=c, label=f"sextile {g}")
        print("%s, sextile %d: mean normalized rate per bin from %.3f to %.3f" % (kind, g, np.nanmin(m), np.nanmax(m)))
    ax[j].set_title(kind); ax[j].set_xlabel("normalized time (bin)")
ax[0].set_ylabel("rate / interval mean"); ax[1].legend(fontsize=7)
plt.show()
print("\nPer session, R of sextiles 1 and 6 over normalized time (pooled within the session); Wilcoxon over sessions:")
for kind in ("episode", "packet"):
    for g in (1, 6):
        rs = []
        for sid in sorted({c[2].session_id for c in courses[kind]}):
            try:
                rs.append(pooled_r([(x, r, s) for x, r, s, *_ in courses[kind] if s.session_id == sid], [g])[0])
            except ValueError:
                pass
        rs = np.array(rs)
        rs = rs[np.isfinite(rs)]
        print("  %s, sextile %d: %d sessions, median R %+.3f, positive in %d, Wilcoxon p = %.2g" % (
            kind, g, len(rs), np.median(rs), (rs > 0).sum(), st.wilcoxon(rs).pvalue))

episode: sextile 1 R = +0.060 (p 6.3e-06); sextile 2 R = +0.094 (p 6.4e-13); sextile 3 R = +0.047 (p 0.00034); sextile 4 R = +0.029 (p 0.025); sextile 5 R = -0.010 (p 0.47); sextile 6 R = -0.083 (p 3.7e-08); all pE R = +0.003 (p 0.81)

packet: sextile 1 R = +0.061 (p 2e-08); sextile 2 R = +0.082 (p 1e-14); sextile 3 R = +0.045 (p 1.8e-05); sextile 4 R = +0.015 (p 0.16); sextile 5 R = +0.022 (p 0.04); sextile 6 R = -0.055 (p 4.4e-06); all pE R = +0.046 (p 8.3e-06)

Variants (sextile 1 and sextile 6 R, p):

  episode, 5 bins, min 20 s: s1 +0.090 (p 1.6e-06), s6 -0.126 (p 2.3e-09)

  episode, 20 bins, min 20 s: s1 +0.043 (p 5.8e-06), s6 -0.061 (p 1.1e-08)

  episode, 10 bins, min 60 s: s1 +0.057 (p 7.8e-05), s6 -0.156 (p 3.2e-21)

  packet, 5 bins, min 20 s: s1 +0.088 (p 8.4e-09), s6 -0.083 (p 1.1e-06)

  packet, 20 bins, min 20 s: s1 +0.043 (p 2.6e-08), s6 -0.039 (p 3.8e-06)

  packet, 10 bins, min 60 s: s1 +0.076 (p 8e-10), s6 -0.113 (p 2.9e-16)

Per subject, packets, sextile 1 and 6:
  BWRat17: s1 +0.137 (p 0.0086, n 370), s6 -0.382 (p 8.8e-05, n 100)
  BWRat18: s1 -0.059 (p 0.42, n 190), s6 -0.351 (p 6.7e-07, n 190)
  BWRat19: s1 +0.241 (p 6.5e-05, n 270), s6 +0.040 (p 0.52, n 270)
  BWRat20: s1 -0.004 (p 0.91, n 860), s6 +0.050 (p 0.1, n 1080)
  BWRat21: s1 +0.033 (p 0.18, n 1670), s6 n/a
  Bogey: s1 +0.173 (p 3.3e-07, n 860), s6 -0.067 (p 0.048, n 860)
  Dino: s1 +0.041 (p 0.078, n 1840), s6 -0.079 (p 0.00063, n 1850)
  Rizzo: s1 -0.010 (p 0.69, n 1570), s6 +0.121 (p 1.2e-06, n 1590)
  Splinter: s1 +0.250 (p 1.5e-14, n 920), s6 -0.090 (p 0.0064, n 920)
  Templeton: s1 n/a, s6 n/a
episode, sextile 1: mean normalized rate per bin from 0.924 to 1.112
episode, sextile 2: mean normalized rate per bin from 0.951 to 1.108
episode, sextile 3: mean normalized rate per bin from 0.961 to 1.101
episode, sextile 4: mean normalized rate per bin from 0.958 to 1.152
episode, sextile 5: mean normalized rate per bin from 0.967 to 1.074
episode, sextile 6: mean normalized rate per bin from 0.956 to 1.062
packet, sextile 1: mean normalized rate per bin from 0.917 to 1.069

packet, sextile 2: mean normalized rate per bin from 0.933 to 1.057
packet, sextile 3: mean normalized rate per bin from 0.953 to 1.061
packet, sextile 4: mean normalized rate per bin from 0.945 to 1.083
packet, sextile 5: mean normalized rate per bin from 0.959 to 1.074
packet, sextile 6: mean normalized rate per bin from 0.960 to 1.064

Per session, R of sextiles 1 and 6 over normalized time (pooled within the session); Wilcoxon over sessions:
  episode, sextile 1: 19 sessions, median R +0.064, positive in 11, Wilcoxon p = 0.049

def cv_r(cs):
    xs, ys = [], []
    for x, r, s, *_ in cs:
        mu, sd = r.mean(axis=0), r.std(axis=0)
        ok = mu > 0
        xs.append(x[ok]); ys.append(sd[ok] / mu[ok])
    return st.pearsonr(np.concatenate(xs), np.concatenate(ys)) + (len(np.concatenate(xs)),)

for kind in ("episode", "packet"):
    r, p, n = cv_r(courses[kind])
    print("%s: R = %.3f, p = %.2g, points %d" % (kind, r, p, n))
    for nb, md in ((5, 20.0), (20, 20.0), (10, 60.0)):
        print("   variant %d bins, min %d s: R = %.3f, p = %.2g, points %d" % (nb, md, *cv_r(course_c(kind, nb, md))))

episode: R = -0.005, p = 0.71, points 6151
   variant 5 bins, min 20 s: R = -0.007, p = 0.72, points 3095
   variant 20 bins, min 20 s: R = -0.007, p = 0.45, points 12177
   variant 10 bins, min 60 s: R = -0.011, p = 0.43, points 5161
packet: R = 0.011, p = 0.28, points 9519
   variant 5 bins, min 20 s: R = 0.008, p = 0.58, points 4760

   variant 20 bins, min 20 s: R = 0.008, p = 0.25, points 19020
   variant 10 bins, min 60 s: R = -0.001, p = 0.91, points 7270

def burst_r(cs):
    xs, ys = [], []
    for x, r, s, short, alln in cs:
        ok = alln > 0
        xs.append(x[ok]); ys.append(short[ok] / alln[ok])
    return st.pearsonr(np.concatenate(xs), np.concatenate(ys)) + (len(np.concatenate(xs)),)

for kind in ("episode", "packet"):
    r, p, n = burst_r(courses[kind])
    print("%s: R = %.3f, p = %.2g, points %d" % (kind, r, p, n))
    for nb, md in ((5, 20.0), (20, 20.0), (10, 60.0)):
        print("   variant %d bins, min %d s: R = %.3f, p = %.2g, points %d" % (nb, md, *burst_r(course_c(kind, nb, md))))

episode: R = 0.029, p = 0.022, points 6151
   variant 5 bins, min 20 s: R = 0.033, p = 0.069, points 3095
   variant 20 bins, min 20 s: R = 0.026, p = 0.0038, points 12174
   variant 10 bins, min 60 s: R = 0.047, p = 0.0008, points 5161
packet: R = 0.024, p = 0.021, points 9519

   variant 5 bins, min 20 s: R = 0.030, p = 0.042, points 4760
   variant 20 bins, min 20 s: R = 0.020, p = 0.0061, points 19012
   variant 10 bins, min 60 s: R = 0.034, p = 0.0038, points 7270

def lfp_course(kind, nbins=10, min_dur=20.0, sess=None):
    out = {b: ([], []) for b in ss.BANDS}
    for s, a, b, inc, rows_ in intervals_of(kind, min_dur, sess):
        rows_ = [i for i in rows_ if i in s.lfp]
        if not rows_:
            continue
        t = np.concatenate([s.lfp[i][0] for i in rows_])
        k = np.clip(((t - a) / (b - a) * nbins).astype(int), 0, nbins - 1)
        for band in ss.BANDS:
            v = np.concatenate([s.lfp[i][1][band] for i in rows_])
            m = np.array([v[k == j].mean() if np.any(k == j) else np.nan for j in range(nbins)])
            ok = np.isfinite(m)
            if ok.sum() < 2:
                continue
            out[band][0].append(((np.arange(nbins) + 0.5) / nbins)[ok]); out[band][1].append(m[ok] - m[ok].mean())
    return {band: st.pearsonr(np.concatenate(x), np.concatenate(y)) + (len(np.concatenate(x)),) for band, (x, y) in out.items()}

for kind in ("episode", "packet"):
    res = lfp_course(kind)
    print("%s: " % kind + "; ".join("%s R = %+.3f (p %.2g, n %d)" % (b, *v) for b, v in res.items()))
    for nb, md in ((5, 20.0), (20, 20.0), (10, 60.0)):
        res = lfp_course(kind, nb, md)
        print("   variant %d bins, min %d s: " % (nb, md) + "; ".join("%s %+.3f (p %.2g)" % (b, v[0], v[1]) for b, v in res.items()))
print("\nPer subject, packets:")
for subj in sorted(tab.subject.unique()):
    res = lfp_course("packet", sess=[s for s in sessions if s.subject == subj])
    print("  %s: " % subj + "; ".join("%s %+.3f (p %.2g)" % (b, v[0], v[1]) for b, v in res.items()))

episode: delta R = +0.069 (p 1.1e-07, n 5880); sigma R = +0.135 (p 3.1e-25, n 5880); gamma R = +0.062 (p 1.7e-06, n 5880)

   variant 5 bins, min 20 s: delta +0.104 (p 1.6e-08); sigma +0.193 (p 3.7e-26); gamma +0.065 (p 0.00036)

   variant 20 bins, min 20 s: delta +0.049 (p 1.4e-07); sigma +0.101 (p 6e-28); gamma +0.054 (p 5.4e-09)

   variant 10 bins, min 60 s: delta +0.135 (p 1.6e-21); sigma +0.249 (p 8.4e-71); gamma +0.009 (p 0.51)

packet: delta R = -0.005 (p 0.61, n 9150); sigma R = +0.104 (p 1.6e-23, n 9150); gamma R = +0.182 (p 4.5e-69, n 9150)

   variant 5 bins, min 20 s: delta +0.009 (p 0.53); sigma +0.154 (p 8.6e-26); gamma +0.213 (p 3.5e-48)

   variant 20 bins, min 20 s: delta -0.009 (p 0.2); sigma +0.075 (p 1.7e-24); gamma +0.152 (p 1.5e-95)
   variant 10 bins, min 60 s: delta +0.022 (p 0.063); sigma +0.214 (p 1.4e-73); gamma +0.156 (p 2.9e-39)

Per subject, packets:
  BWRat17: delta +0.054 (p 0.3); sigma +0.116 (p 0.025); gamma +0.237 (p 4e-06)
  BWRat18: delta +0.317 (p 8.2e-06); sigma +0.417 (p 2.3e-09); gamma -0.066 (p 0.37)
  BWRat19: delta +0.266 (p 9.4e-06); sigma +0.244 (p 5.2e-05); gamma +0.108 (p 0.075)
  BWRat20: delta +0.335 (p 8.8e-30); sigma +0.373 (p 6e-37); gamma -0.105 (p 0.00054)

  BWRat21: delta -0.004 (p 0.88); sigma +0.160 (p 4.9e-11); gamma +0.283 (p 4.8e-32)
  Bogey: delta -0.143 (p 2.5e-05); sigma +0.003 (p 0.93); gamma +0.440 (p 5.4e-42)
  Dino: delta +0.036 (p 0.17); sigma +0.046 (p 0.079); gamma +0.209 (p 3.9e-16)
  Rizzo: delta -0.110 (p 1.2e-05); sigma +0.004 (p 0.86); gamma +0.133 (p 9.5e-08)
  Splinter: delta -0.092 (p 0.005); sigma +0.069 (p 0.035); gamma +0.055 (p 0.098)
  Templeton: delta -0.091 (p 0.015); sigma +0.154 (p 3.2e-05); gamma +0.214 (p 6.4e-09)

pk = []
for s in sessions:
    if not s.lfp:
        continue
    allv = {b: np.concatenate([10 ** s.lfp[i][1][b] for i in s.lfp]) for b in ss.BANDS}
    for i, (t, pw) in s.lfp.items():
        if s.stops[i] - s.starts[i] < 20:
            continue
        rr = {}
        for g in (1, 6):
            units = s.pe[sext[s.path][1] == g]
            rr[g] = np.mean(ss.rates(s.spikes, units, s.starts[i], s.stops[i])) if len(units) else np.nan
        last = t >= s.stops[i] - 20
        pk.append(dict(subject=s.subject, session=s.session_id, row=i, dur=s.stops[i] - s.starts[i],
                       delta_norm=np.mean(10 ** pw["delta"]) / allv["delta"].mean(),
                       delta=pw["delta"].mean(), gamma=pw["gamma"].mean(), sigma_last20=pw["sigma"][last].mean(),
                       s1=rr[1], s6=rr[6]))
pk = pd.DataFrame(pk)
r, p = st.pearsonr(pk.dur, pk.delta_norm)
print("Packets of at least 20 s with LFP: %d; R = %.3f, p = %.2g; within 0.279-0.341: %s" % (len(pk), r, p, 0.279 <= r <= 0.341))
print("Spearman rho = %.3f; R with log10 duration = %.3f" % (st.spearmanr(pk.dur, pk.delta_norm)[0], st.pearsonr(np.log10(pk.dur), pk.delta_norm)[0]))
print("Normalized delta power: %.2f to %.2f; durations %.0f to %.0f s" % (pk.delta_norm.min(), pk.delta_norm.max(), pk.dur.min(), pk.dur.max()))
print("Per subject R:", {subj: round(st.pearsonr(g.dur, g.delta_norm)[0], 3) for subj, g in pk.groupby("subject")})

Packets of at least 20 s with LFP: 915; R = 0.334, p = 2.8e-25; within 0.279-0.341: True
Spearman rho = 0.555; R with log10 duration = 0.385
Normalized delta power: 0.12 to 3.81; durations 20 to 1199 s
Per subject R: {np.str_('BWRat17'): np.float64(0.348), np.str_('BWRat18'): np.float64(0.742), np.str_('BWRat19'): np.float64(0.532), np.str_('BWRat20'): np.float64(0.33), np.str_('BWRat21'): np.float64(0.558), np.str_('Bogey'): np.float64(0.509), np.str_('Dino'): np.float64(0.459), np.str_('Rizzo'): np.float64(0.644), np.str_('Splinter'): np.float64(0.524), np.str_('Templeton'): np.float64(-0.101)}

z = lambda v: (v - v.mean()) / v.std() if v.std() > 0 else v * np.nan
pz = pk.copy()
for c in ("delta", "gamma", "sigma_last20", "s1", "s6"):
    pz[c] = pk.groupby("session")[c].transform(z)
out = {}
for name, a, b, exp_, tol in (("s6 vs delta", "s6", "delta", -0.28, 0.028), ("s6 vs gamma", "s6", "gamma", 0.44, 0.044),
                              ("s1 vs sigma last 20 s", "s1", "sigma_last20", 0.25, 0.025)):
    d = pz[[a, b]].dropna()
    r, p = st.pearsonr(d[a], d[b])
    out[name] = (r, p)
    print("%s: R = %.3f, p = %.2g, n %d; within tolerance %s; sign as stated %s" % (
        name, r, p, len(d), abs(r - exp_) <= tol, np.sign(r) == np.sign(exp_)))
print("\nPer subject:")
for subj, g in pz.groupby("subject"):
    vals = []
    for a, b in (("s6", "delta"), ("s6", "gamma"), ("s1", "sigma_last20")):
        d = g[[a, b]].dropna()
        vals.append(st.pearsonr(d[a], d[b])[0] if len(d) > 2 else np.nan)
    print("  %s: s6-delta %.3f, s6-gamma %.3f, s1-sigma %.3f, packets %d" % (subj, *vals, len(g)))

s6 vs delta: R = -0.416, p = 1.3e-28, n 649; within tolerance False; sign as stated True
s6 vs gamma: R = 0.570, p = 4.2e-57, n 649; within tolerance False; sign as stated True
s1 vs sigma last 20 s: R = 0.104, p = 0.0029, n 821; within tolerance False; sign as stated True

Per subject:
  BWRat17: s6-delta 0.615, s6-gamma 0.703, s1-sigma -0.127, packets 37
  BWRat18: s6-delta -0.097, s6-gamma 0.086, s1-sigma 0.005, packets 19
  BWRat19: s6-delta -0.761, s6-gamma 0.443, s1-sigma 0.797, packets 27
  BWRat20: s6-delta -0.489, s6-gamma 0.590, s1-sigma -0.165, packets 108
  BWRat21: s6-delta nan, s6-gamma nan, s1-sigma -0.144, packets 167
  Bogey: s6-delta -0.241, s6-gamma 0.740, s1-sigma 0.328, packets 86
  Dino: s6-delta -0.493, s6-gamma 0.401, s1-sigma 0.254, packets 148
  Rizzo: s6-delta -0.664, s6-gamma 0.606, s1-sigma 0.041, packets 159
  Splinter: s6-delta -0.013, s6-gamma 0.709, s1-sigma 0.384, packets 92
  Templeton: s6-delta nan, s6-gamma nan, s1-sigma nan, packets 72

def triplets(a_state, b_state, states=None):
    """Triplets over all sessions; `states` maps a session path to another state labelling."""
    out = []
    for s in sessions:
        lab = s.state if states is None else states[s.path]
        for i in range(1, len(lab) - 1):
            if lab[i - 1] == a_state and lab[i] == b_state and lab[i + 1] == a_state:
                r1 = ss.rates(s.spikes, s.pe, s.starts[i - 1], s.stops[i - 1])
                r2 = ss.rates(s.spikes, s.pe, s.starts[i + 1], s.stops[i + 1])
                for k, u in enumerate(s.pe):
                    out.append(dict(subject=s.subject, session=s.session_id, row=i, unit=u, sextile=sext[s.path][1][k],
                                    pre=r1[k], post=r2[k], ch=r2[k] - r1[k], mid_dur=s.stops[i] - s.starts[i]))
    return pd.DataFrame(out)

def wp(v, **kw):
    """Wilcoxon signed-rank p, or NaN when there is nothing to test."""
    v = np.asarray(v, float)
    return st.wilcoxon(v, **kw).pvalue if len(v) and np.any(v != 0) else np.nan

def triplet_report(d, alpha):
    print("triplets %d, cell-triplet pairs %d" % (d.groupby(["session", "row"]).ngroups, len(d)))
    res = {}
    for g in range(1, 7):
        x = d[d.sextile == g]
        p = st.wilcoxon(x.ch).pvalue if (x.ch != 0).any() else 1.0
        cm = x.groupby(["session", "unit"]).ch.mean()
        pc = st.wilcoxon(cm).pvalue if (cm != 0).any() else 1.0
        res[g] = (x.ch.mean(), p)
        print("  sextile %d: n %d, mean change %+.4f Hz (pre mean %.3f), Wilcoxon p = %.2g%s | per-cell means: n %d, mean %+.4f, p = %.2g" % (
            g, len(x), x.ch.mean(), x.pre.mean(), p, " *" if p < alpha else "", len(cm), cm.mean(), pc))
    print("  one value per triplet (mean over its cells), and per session; Wilcoxon over these:")
    for g in range(1, 7):
        x = d[d.sextile == g]
        t_, s_ = x.groupby(["session", "row"]).ch.mean(), x.groupby("session").ch.mean()
        print("    sextile %d: triplets %d, mean %+.4f, p = %.2g | sessions %d, mean %+.4f, p = %.2g, lower in %d" % (
            g, len(t_), t_.mean(), wp(t_), len(s_), s_.mean(), wp(s_), (s_ < 0).sum()))
    return res

def triplet_by_subject(d, gs):
    print("  per subject, mean change (Hz) for sextiles %s:" % (gs,))
    print(d[d.sextile.isin(gs)].groupby(["subject", "sextile"]).ch.mean().unstack().round(4).to_string())

trip = {k: triplets(*k) for k in (("MA", "NREM"), ("NREM", "WAKE"), ("NREM", "MA"), ("NREM", "REM"))}

res_ma = triplet_report(trip[("MA", "NREM")], 0.01)
triplet_by_subject(trip[("MA", "NREM")], [5, 6])

triplets 148, cell-triplet pairs 6073
  sextile 1: n 1137, mean change +0.0062 Hz (pre mean 0.079), Wilcoxon p = 0.32 | per-cell means: n 128, mean +0.0103, p = 0.35
  sextile 2: n 977, mean change -0.0076 Hz (pre mean 0.278), Wilcoxon p = 0.68 | per-cell means: n 130, mean -0.0175, p = 0.37
  sextile 3: n 988, mean change -0.0240 Hz (pre mean 0.551), Wilcoxon p = 0.039 | per-cell means: n 126, mean -0.0277, p = 0.14
  sextile 4: n 1001, mean change +0.0029 Hz (pre mean 0.744), Wilcoxon p = 0.7 | per-cell means: n 126, mean -0.0177, p = 0.31
  sextile 5: n 1027, mean change -0.0107 Hz (pre mean 1.177), Wilcoxon p = 0.24 | per-cell means: n 132, mean -0.0045, p = 0.79
  sextile 6: n 943, mean change -0.1037 Hz (pre mean 3.325), Wilcoxon p = 0.077 | per-cell means: n 135, mean -0.1197, p = 0.029
  one value per triplet (mean over its cells), and per session; Wilcoxon over these:
    sextile 1: triplets 129, mean +0.0043, p = 0.19 | sessions 18, mean +0.0142, p = 0.18, lower in 7
    sextile 2: triplets 140, mean -0.0047, p = 1 | sessions 19, mean -0.0118, p = 0.54, lower in 9
    sextile 3: triplets 145, mean -0.0415, p = 0.089 | sessions 19, mean -0.0748, p = 0.1, lower in 12
    sextile 4: triplets 145, mean -0.0086, p = 0.26 | sessions 19, mean -0.0283, p = 0.4, lower in 11
    sextile 5: triplets 146, mean -0.0220, p = 0.29 | sessions 19, mean -0.0196, p = 0.95, lower in 8
    sextile 6: triplets 110, mean -0.0864, p = 0.12 | sessions 15, mean -0.0946, p = 0.17, lower in 10
  per subject, mean change (Hz) for sextiles [5, 6]:
sextile         5       6
subject                  
BWRat17    0.1158 -0.0487
BWRat19   -0.1234 -0.0844
BWRat20    0.0839 -0.2476
BWRat21    0.0455     NaN
Bogey     -0.0232 -0.1632
Dino      -0.0625 -0.0056
Rizzo      0.0191 -0.1611
Splinter   0.0258  0.1326
Templeton -0.1159     NaN

res_w = triplet_report(trip[("NREM", "WAKE")], 0.05)
d = trip[("NREM", "WAKE")]
print("WAKE durations in these triplets (s): median %.0f, range %.0f to %.0f" % (
    d.groupby(["session", "row"]).mid_dur.first().median(), d.mid_dur.min(), d.mid_dur.max()))
triplet_by_subject(d, [1, 6])
print("\nTriplets and the change of sextile 6 by the minimum duration of the intervening WAKE:")
for mn in (40, 60, 120, 180, 300, 420, 600, 900):
    x = d[d.mid_dur >= mn]
    x6 = x[x.sextile == 6]
    print("  WAKE of at least %d s: triplets %d, sextile 6 n %d, mean change %+.4f Hz, Wilcoxon p = %.2g" % (
        mn, x.groupby(["session", "row"]).ngroups, len(x6), x6.ch.mean(), wp(x6.ch)))

triplets 265, cell-triplet pairs 10431
  sextile 1: n 1773, mean change +0.0049 Hz (pre mean 0.106), Wilcoxon p = 0.0063 * | per-cell means: n 136, mean +0.0083, p = 0.013
  sextile 2: n 1709, mean change +0.0079 Hz (pre mean 0.323), Wilcoxon p = 0.11 | per-cell means: n 135, mean +0.0193, p = 0.072
  sextile 3: n 1640, mean change +0.0114 Hz (pre mean 0.634), Wilcoxon p = 0.015 * | per-cell means: n 135, mean +0.0157, p = 0.055
  sextile 4: n 1749, mean change +0.0115 Hz (pre mean 0.844), Wilcoxon p = 0.21 | per-cell means: n 135, mean +0.0068, p = 0.15
  sextile 5: n 1817, mean change +0.0240 Hz (pre mean 1.275), Wilcoxon p = 0.37 | per-cell means: n 135, mean +0.0506, p = 0.25
  sextile 6: n 1743, mean change -0.0509 Hz (pre mean 3.030), Wilcoxon p = 0.092 | per-cell means: n 136, mean -0.0941, p = 0.17
  one value per triplet (mean over its cells), and per session; Wilcoxon over these:
    sextile 1: triplets 241, mean +0.0030, p = 0.53 | sessions 19, mean +0.0074, p = 0.13, lower in 6
    sextile 2: triplets 252, mean +0.0138, p = 0.22 | sessions 20, mean +0.0286, p = 0.0056, lower in 6
    sextile 3: triplets 259, mean +0.0143, p = 0.16 | sessions 20, mean +0.0193, p = 0.083, lower in 7
    sextile 4: triplets 259, mean +0.0143, p = 0.75 | sessions 20, mean +0.0114, p = 0.81, lower in 10
    sextile 5: triplets 258, mean +0.0312, p = 0.2 | sessions 20, mean +0.0666, p = 0.097, lower in 6
    sextile 6: triplets 203, mean -0.0474, p = 0.27 | sessions 16, mean -0.0798, p = 0.13, lower in 11
WAKE durations in these triplets (s): median 102, range 40 to 3945
  per subject, mean change (Hz) for sextiles [1, 6]:
sextile        1       6
subject                 
BWRat17  -0.0004 -0.0530
BWRat18   0.0025 -0.0591
BWRat19   0.0394 -0.0698
BWRat20   0.0019 -0.1636
BWRat21   0.0003     NaN
Bogey     0.0082 -0.0589
Dino      0.0021 -0.1149
Rizzo    -0.0034 -0.0111
Splinter  0.0215  0.2449

Triplets and the change of sextile 6 by the minimum duration of the intervening WAKE:
  WAKE of at least 40 s: triplets 265, sextile 6 n 1743, mean change -0.0509 Hz, Wilcoxon p = 0.092
  WAKE of at least 60 s: triplets 193, sextile 6 n 1415, mean change -0.0277 Hz, Wilcoxon p = 0.43
  WAKE of at least 120 s: triplets 120, sextile 6 n 937, mean change +0.0163 Hz, Wilcoxon p = 0.59
  WAKE of at least 180 s: triplets 93, sextile 6 n 800, mean change +0.0497 Hz, Wilcoxon p = 0.32
  WAKE of at least 300 s: triplets 61, sextile 6 n 526, mean change +0.0063 Hz, Wilcoxon p = 0.61
  WAKE of at least 420 s: triplets 47, sextile 6 n 422, mean change -0.0338 Hz, Wilcoxon p = 0.8
  WAKE of at least 600 s: triplets 36, sextile 6 n 319, mean change -0.0642 Hz, Wilcoxon p = 0.87
  WAKE of at least 900 s: triplets 23, sextile 6 n 184, mean change -0.1854 Hz, Wilcoxon p = 0.5

res_m = triplet_report(trip[("NREM", "MA")], 0.05)
triplet_by_subject(trip[("NREM", "MA")], [1, 6])

triplets 391, cell-triplet pairs 15436
  sextile 1: n 2846, mean change +0.0034 Hz (pre mean 0.110), Wilcoxon p = 0.082 | per-cell means: n 136, mean +0.0020, p = 0.76
  sextile 2: n 2533, mean change +0.0099 Hz (pre mean 0.325), Wilcoxon p = 0.00019 * | per-cell means: n 135, mean +0.0232, p = 0.0027
  sextile 3: n 2458, mean change +0.0016 Hz (pre mean 0.633), Wilcoxon p = 1 | per-cell means: n 135, mean -0.0074, p = 0.74
  sextile 4: n 2533, mean change -0.0114 Hz (pre mean 0.809), Wilcoxon p = 0.48 | per-cell means: n 135, mean -0.0255, p = 0.0086
  sextile 5: n 2631, mean change -0.0302 Hz (pre mean 1.222), Wilcoxon p = 0.0023 * | per-cell means: n 135, mean -0.0522, p = 0.00016
  sextile 6: n 2435, mean change -0.0725 Hz (pre mean 2.909), Wilcoxon p = 6.8e-05 * | per-cell means: n 136, mean -0.1095, p = 2.5e-07
  one value per triplet (mean over its cells), and per session; Wilcoxon over these:
    sextile 1: triplets 344, mean +0.0064, p = 0.037 | sessions 19, mean +0.0047, p = 0.17, lower in 8
    sextile 2: triplets 371, mean +0.0126, p = 0.004 | sessions 20, mean +0.0288, p = 0.024, lower in 5
    sextile 3: triplets 382, mean +0.0016, p = 0.98 | sessions 20, mean -0.0083, p = 0.81, lower in 9
    sextile 4: triplets 382, mean -0.0163, p = 0.33 | sessions 20, mean -0.0209, p = 0.28, lower in 11
    sextile 5: triplets 380, mean -0.0327, p = 0.003 | sessions 20, mean -0.0704, p = 0.0094, lower in 16
    sextile 6: triplets 289, mean -0.0785, p = 0.00095 | sessions 16, mean -0.1117, p = 0.00015, lower in 14
  per subject, mean change (Hz) for sextiles [1, 6]:
sextile        1       6
subject                 
BWRat17   0.0094 -0.0925
BWRat18  -0.0682 -0.0655
BWRat19  -0.0017 -0.1663
BWRat20  -0.0016 -0.0722
BWRat21  -0.0001     NaN
Bogey     0.0105  0.0010
Dino      0.0104 -0.1161
Rizzo     0.0056 -0.0160
Splinter  0.0112 -0.1387

res_r = triplet_report(trip[("NREM", "REM")], 0.05)
triplet_by_subject(trip[("NREM", "REM")], [1, 6])
print("REM epochs followed by Non-REM: %d of %d REM epochs" % (
    sum(int(np.sum((s.state[:-1] == "REM") & (s.state[1:] == "NREM"))) for s in sessions),
    sum(int(np.sum(s.state == "REM")) for s in sessions)))

triplets 17, cell-triplet pairs 465
  sextile 1: n 80, mean change +0.0032 Hz (pre mean 0.124), Wilcoxon p = 0.51 | per-cell means: n 38, mean +0.0075, p = 0.11
  sextile 2: n 82, mean change +0.0039 Hz (pre mean 0.430), Wilcoxon p = 0.85 | per-cell means: n 44, mean -0.0044, p = 0.75
  sextile 3: n 101, mean change -0.0315 Hz (pre mean 0.588), Wilcoxon p = 0.14 | per-cell means: n 43, mean -0.0577, p = 0.12
  sextile 4: n 74, mean change -0.1146 Hz (pre mean 0.947), Wilcoxon p = 0.022 * | per-cell means: n 43, mean -0.1830, p = 0.0093
  sextile 5: n 47, mean change -0.1885 Hz (pre mean 1.609), Wilcoxon p = 0.011 * | per-cell means: n 39, mean -0.2005, p = 0.0082
  sextile 6: n 81, mean change -0.2656 Hz (pre mean 2.677), Wilcoxon p = 0.19 | per-cell means: n 57, mean -0.3989, p = 0.084
  one value per triplet (mean over its cells), and per session; Wilcoxon over these:
    sextile 1: triplets 15, mean +0.0007, p = 0.98 | sessions 6, mean +0.0125, p = 0.84, lower in 2
    sextile 2: triplets 17, mean +0.0306, p = 0.4 | sessions 8, mean +0.0150, p = 0.95, lower in 5
    sextile 3: triplets 16, mean -0.0433, p = 0.13 | sessions 7, mean -0.0834, p = 0.11, lower in 5
    sextile 4: triplets 16, mean -0.0765, p = 0.044 | sessions 7, mean -0.1082, p = 0.016, lower in 7
    sextile 5: triplets 13, mean -0.0762, p = 0.45 | sessions 7, mean -0.0410, p = 0.38, lower in 5
    sextile 6: triplets 12, mean -0.0096, p = 0.47 | sessions 6, mean -0.0702, p = 0.56, lower in 2
  per subject, mean change (Hz) for sextiles [1, 6]:
sextile        1       6
subject                 
BWRat20   0.0069  0.2135
BWRat21  -0.0208     NaN
Dino     -0.0182  0.1931
Rizzo     0.0617 -0.7548
Splinter  0.0086 -0.0201
REM epochs followed by Non-REM: 49 of 534 REM epochs

d = trip[("NREM", "REM")]
per = d.groupby(["session", "row"]).agg(pre=("pre", "mean"), post=("post", "mean"), dur=("mid_dur", "first"),
                                        subject=("subject", "first"))
per = per[per.pre > 0]
ra = st.pearsonr(per.dur, (per.post - per.pre) / per.pre)
dc = d[d.pre > 0]
rb = st.pearsonr(dc.mid_dur, (dc.post - dc.pre) / dc.pre)
print("(a) per triplet: R = %.3f, p = %.2g, n %d | (b) per cell and triplet: R = %.3f, p = %.2g, n %d" % (
    ra[0], ra[1], len(per), rb[0], rb[1], len(dc)))
print("Within -0.242 to -0.198:", {"a": -0.242 <= ra[0] <= -0.198, "b": -0.242 <= rb[0] <= -0.198})
print("REM durations %.0f to %.0f s; fractional change (a) %.2f to %.2f" % (
    per.dur.min(), per.dur.max(), ((per.post - per.pre) / per.pre).min(), ((per.post - per.pre) / per.pre).max()))
print("Spearman (a): rho = %.3f, p = %.2g" % st.spearmanr(per.dur, (per.post - per.pre) / per.pre))

(a) per triplet: R = -0.047, p = 0.86, n 17 | (b) per cell and triplet: R = -0.061, p = 0.21, n 436
Within -0.242 to -0.198: {'a': np.False_, 'b': np.False_}
REM durations 8 to 157 s; fractional change (a) -0.31 to 0.53
Spearman (a): rho = -0.190, p = 0.46

sig_all = {}
for s in sessions:
    if not s.lfp:
        continue
    v = {i: s.lfp[i][1]["sigma"][s.lfp[i][0] >= s.stops[i] - 20].mean() for i in s.lfp}
    vals = np.array(list(v.values()))
    sig_all[s.path] = {i: (x - np.nanmean(vals)) / np.nanstd(vals) for i, x in v.items()}
d = trip[("MA", "NREM")]
d1 = d[d.sextile == 1].groupby(["session", "row"]).agg(ch=("ch", "mean")).reset_index()
path_of = {s.session_id: s.path for s in sessions}
d1["sigma"] = [sig_all.get(path_of[a], {}).get(b, np.nan) for a, b in zip(d1.session, d1.row)]
d1 = d1.dropna()
r, p = st.pearsonr(d1.sigma, d1.ch)
print("MA-packet-MA triplets with sextile 1 cells and LFP: %d; R = %.3f, p = %.2g; within 0.163-0.199: %s" % (len(d1), r, p, 0.163 <= r <= 0.199))
print("Spearman rho = %.3f, p = %.2g" % st.spearmanr(d1.sigma, d1.ch))

MA-packet-MA triplets with sextile 1 cells and LFP: 125; R = 0.033, p = 0.71; within 0.163-0.199: False
Spearman rho = 0.023, p = 0.79

dd = []
for sess, g in pk.groupby("session", sort=False):
    g = g.sort_values("row")
    dd.append(pd.DataFrame({"subject": g.subject.values[1:], "session": sess, "d_s1": np.diff(g.s1.values), "d_delta": np.diff(g.delta_norm.values),
                            "d_s6": np.diff(g.s6.values)}))
dd = pd.concat(dd).dropna()
r, p = st.pearsonr(dd.d_delta, dd.d_s1)
print("Consecutive packet pairs (both at least 20 s): %d; sextile 1: R = %.3f, p = %.2g" % (len(dd), r, p))
print("Sextile 6 for comparison: R = %.3f, p = %.2g" % st.pearsonr(dd.d_delta, dd.d_s6))
print("Per subject, sextile 1 R:", {subj: round(st.pearsonr(g.d_delta, g.d_s1)[0], 3) for subj, g in dd.groupby("subject") if len(g) > 2})
rs = []
for sess, g in dd.groupby("session"):
    if len(g) >= 5:
        rs.append(st.pearsonr(g.d_delta, g.d_s1)[0])
rs = np.array(rs)
print("Per session (at least 5 packet pairs): %d sessions, median R %+.3f, positive in %d, Wilcoxon over sessions p = %.2g" % (
    len(rs), np.median(rs), (rs > 0).sum(), st.wilcoxon(rs).pvalue))

Consecutive packet pairs (both at least 20 s): 613; sextile 1: R = 0.183, p = 5.3e-06

Sextile 6 for comparison: R = -0.368, p = 4.7e-21
Per subject, sextile 1 R: {np.str_('BWRat17'): np.float64(-0.359), np.str_('BWRat18'): np.float64(-0.444), np.str_('BWRat19'): np.float64(0.747), np.str_('BWRat20'): np.float64(-0.266), np.str_('Bogey'): np.float64(0.516), np.str_('Dino'): np.float64(0.305), np.str_('Rizzo'): np.float64(0.094), np.str_('Splinter'): np.float64(-0.069)}
Per session (at least 5 packet pairs): 14 sessions, median R +0.025, positive in 7, Wilcoxon over sessions p = 0.43

wk = course("WAKE", 10, 40.0)
rm = course("REM", 10, 20.0)
tw, trm = sextile_table(wk), sextile_table(rm)
print("WAKE epochs %d, REM epochs %d" % (len(wk), len(rm)))
for name, t in (("WAKE", tw), ("REM", trm)):
    print(name + ": " + "; ".join("%s R = %+.3f (p %.2g)" % (g, v[0], v[1]) for g, v in t.items()))
wake_ok = all(tw[g][0] > 0 and tw[g][1] < 0.05 for g in (5, 6)) and all(tw[g][0] < 0 and tw[g][1] < 0.05 for g in (1, 2))
rem_ok = all(trm[g][0] > 0 and trm[g][1] < 0.05 for g in (1, 2, 3, 5, 6))
print("WAKE part holds:", wake_ok, "| REM part holds:", rem_ok)

WAKE epochs 495, REM epochs 388
WAKE: 1 R = +0.021 (p 0.16); 2 R = +0.067 (p 5.2e-06); 3 R = +0.028 (p 0.055); 4 R = +0.009 (p 0.55); 5 R = -0.122 (p 3.5e-17); 6 R = -0.120 (p 4.9e-13); all pE R = -0.131 (p 1.6e-20)
REM: 1 R = +0.014 (p 0.41); 2 R = -0.024 (p 0.15); 3 R = +0.015 (p 0.38); 4 R = +0.028 (p 0.088); 5 R = +0.011 (p 0.52); 6 R = -0.005 (p 0.82); all pE R = +0.037 (p 0.021)
WAKE part holds: False | REM part holds: False

drop_labels = sorted(str(x) for x in clock.session[~clock["labels consistent"]])
drop_lfp = sorted(str(x) for x in clock.session[~clock["lfp consistent"]])
all_ids = set(tab.session)
subsets = {"all 21": all_ids, "without " + ", ".join(drop_labels): all_ids - set(drop_labels),
           "without %d (LFP)" % len(drop_lfp): all_ids - set(drop_lfp)}
print("Labels not confirmed:", drop_labels, "| LFP not confirmed:", drop_lfp)
acg_lags = np.arange(300) + 0.5

def headlines(ids):
    S = [s for s in sessions if s.session_id in ids]
    out = []

    def add(claim, name, fn):
        try:
            v = float(fn())
        except (ValueError, KeyError, IndexError, TypeError):
            v = np.nan
        out.append((claim, name, v))

    c = cells_a[cells_a.session.isin(ids)]
    pe = c[c.cls == "pE"]
    add("dataset-composition", "pE cells", lambda: len(pe))
    add("dataset-composition", "pI cells", lambda: (c.cls == "pI").sum())
    for cls in ("pE", "pI"):
        for state in ("WAKE", "NREM", "REM"):
            add("fig2a-median-rates", "%s %s median (Hz)" % (cls, state), lambda: c[c.cls == cls][state].median())
    add("fig2a-skewed-wide-distribution", "WAKE skewness", lambda: st.skew(pe.WAKE))
    add("fig2a-skewed-wide-distribution", "WAKE 2.5th percentile (Hz)", lambda: np.percentile(pe.WAKE, 2.5))
    add("fig2a-skewed-wide-distribution", "WAKE 97.5th percentile (Hz)", lambda: np.percentile(pe.WAKE, 97.5))
    for state in ("WAKE", "REM", "NREM"):
        add("fig2a-fraction-above-2hz", "%% above 2 Hz, %s" % state, lambda: 100 * np.mean(pe[state] > 2))
    for other in ("WAKE", "REM"):
        a_, b_, n_ = (pe.NREM > 2).sum(), (pe[other] > 2).sum(), len(pe)
        add("fig2a-fraction-above-2hz", "chi-square p, NREM against %s" % other,
            lambda: st.chi2_contingency([[a_, n_ - a_], [b_, n_ - b_]])[1])
    k = (pe.WAKE > 0) & (pe.NREM > 0)
    f2b = ols(np.log10(pe.WAKE[k]), np.log10(pe.NREM[k]))
    add("fig2b-wake-nrem-slope", "R, log10", lambda: f2b[3])
    add("fig2b-wake-nrem-slope", "slope, log10", lambda: f2b[0])
    add("fig2b-wake-nrem-slope", "slope, upper 95% bound", lambda: f2b[2])
    bb = bi[bi.session.isin(ids)]
    for rd in ("a", "b"):
        for state in ("NREM", "WAKE", "REM"):
            add("fig2-burst-index-by-state", "median %% %s (%s)" % (state, rd), lambda: 100 * bb["%s_%s" % (state, rd)].median())
    add("fig2-burst-index-by-state", "Kruskal-Wallis p (a)",
        lambda: st.kruskal(*[bb[x + "_a"].dropna() for x in ("NREM", "WAKE", "REM")]).pvalue)
    for state in ("NREM", "REM", "WAKE"):
        mean = np.mean([h for h, sid in zip(acg[state], acg_sess[state]) if sid in ids], axis=0)
        sm = np.convolve(mean, np.ones(21) / 21, mode="same")
        w = (acg_lags >= 60) & (acg_lags <= 250)
        if state == "NREM":
            add("fig2c-acg-peaks", "NREM peak 2-50 ms (ms)", lambda: acg_lags[np.argmax(mean[2:50]) + 2])
        else:
            add("fig2c-acg-peaks", "%s smoothed maximum 60-250 ms (ms)" % state, lambda: acg_lags[w][np.argmax(sm[w])])
    g = pairs[pairs.session.isin(ids)].groupby("a").r.mean()
    add("fig2f-rate-group-correlations", "(6,6) mean r", lambda: g[(6, 6)])
    add("fig2f-rate-group-correlations", "(6,6) the largest (1 if so)", lambda: g.idxmax() == (6, 6))
    add("fig2f-rate-group-correlations", "(1,6) mean r", lambda: g[(1, 6)])
    m = ma[ma.session.isin(ids)]
    add("ma-statistics", "median duration (s)", lambda: np.median(np.concatenate([np.array(x, float) for x in m.durs])))
    add("ma-statistics", "per min of SLEEP", lambda: m.freq.mean())
    add("ma-statistics", "% of SLEEP", lambda: m.frac.mean())
    F = fr[np.isin(fr_sess, list(ids))]
    add("ma-fraction-over-nrem-episodes", "R of bin means", lambda: st.pearsonr(np.arange(10) + 0.5, F.mean(axis=0))[0])
    add("ma-fraction-over-nrem-episodes", "R pooled", lambda: st.pearsonr(np.tile(np.arange(10) + 0.5, len(F)), F.ravel())[0])
    add("packet-duration", "median (s)", lambda: np.median(np.concatenate([(s.stops - s.starts)[s.state == "NREM"] for s in S])))
    add("wake-sleep-cycles", "cycles per session", lambda: nc[nc.index.isin(ids)].mean())
    P = {rd: [t_ for t_ in pts[rd] if t_[2].session_id in ids] for rd in ("a", "b")}
    for rd in ("a", "b"):
        add("fig3b-population-decline", "R (%s)" % rd, lambda: pooled_r(P[rd])[0])
    f_ = fl[fl.session.isin(ids)]
    add("fig3-first-last-packet", "first mean (Hz)", lambda: f_["first"].mean())
    add("fig3-first-last-packet", "last mean (Hz)", lambda: f_["last"].mean())
    add("fig3-first-last-packet", "Wilcoxon p", lambda: st.wilcoxon(f_["first"], f_["last"]).pvalue)
    k = (f_["first"] > 0) & (f_["last"] > 0)
    f3d = ols(np.log10(f_["first"][k]), np.log10(f_["last"][k]))
    add("fig3d-first-last-slope", "slope, log10", lambda: f3d[0])
    add("fig3d-first-last-slope", "slope, upper 95% bound", lambda: f3d[2])
    b_ = ba[ba.session.isin(ids)]
    add("fig3e-wake-before-after", "Wilcoxon p, all pE", lambda: st.wilcoxon(b_.pre, b_.post).pvalue)
    for gg in (1, 4, 5, 6):
        x = b_[b_.sextile == gg]
        add("fig3e-wake-before-after", "sextile %d one-tailed p" % gg,
            lambda: st.wilcoxon(x.post, x.pre, alternative="greater" if gg == 1 else "less").pvalue)
    for gg in (1, 6):
        add("fig3b-sextiles-over-sleep", "sextile %d R (a)" % gg, lambda: pooled_r(P["a"], [gg])[0])
        add("fig3b-sextiles-over-sleep", "sextile %d p (a)" % gg, lambda: pooled_r(P["a"], [gg])[1])
    add("fig3b-cv-over-sleep", "R (a)", lambda: cv_sleep(P["a"])[0])
    add("fig3b-cv-over-sleep", "p (a)", lambda: cv_sleep(P["a"])[1])
    for kind in ("episode", "packet"):
        cs = [c_ for c_ in courses[kind] if c_[2].session_id in ids]
        try:
            t6 = sextile_table(cs)
        except ValueError:
            t6 = {}
        for gg in (1, 6):
            add("fig4iv-sextiles-within-nrem", "%s, sextile %d R" % (kind, gg), lambda: t6[gg][0])
            add("fig4iv-sextiles-within-nrem", "%s, sextile %d p" % (kind, gg), lambda: t6[gg][1])
        add("fig4v-cv-within-nrem", "%s R" % kind, lambda: cv_r(cs)[0])
        add("fig4v-cv-within-nrem", "%s p" % kind, lambda: cv_r(cs)[1])
        add("fig4vi-burst-within-nrem", "%s R" % kind, lambda: burst_r(cs)[0])
        add("fig4vi-burst-within-nrem", "%s p" % kind, lambda: burst_r(cs)[1])
        lc = lfp_course(kind, sess=S)
        for band in ss.BANDS:
            add("fig4i-band-power-within-nrem", "%s, %s R" % (kind, band), lambda: lc[band][0])
            add("fig4i-band-power-within-nrem", "%s, %s p" % (kind, band), lambda: lc[band][1])
    q = pk[pk.session.isin(ids)]
    add("fig4-delta-packet-length", "R", lambda: st.pearsonr(q.dur, q.delta_norm)[0])
    add("fig4-delta-packet-length", "p", lambda: st.pearsonr(q.dur, q.delta_norm)[1])
    qz = pz[pz.session.isin(ids)]
    for a_, b2, nm in (("s6", "delta", "sextile 6 against delta R"), ("s6", "gamma", "sextile 6 against gamma R"),
                       ("s1", "sigma_last20", "sextile 1 against last-20-s sigma R")):
        add("fig4cd-packet-rate-lfp", nm, lambda: st.pearsonr(*qz[[a_, b2]].dropna().values.T)[0])
    for key, claim, gs in ((("MA", "NREM"), "fig5b-nrem-between-ma", (5, 6)), (("NREM", "WAKE"), "fig5b-wake-between-packets", (6,)),
                           (("NREM", "MA"), "fig5b-ma-between-packets", (1, 6)), (("NREM", "REM"), "fig5b-rem-between-packets", ())):
        dt = trip[key][trip[key].session.isin(ids)]
        add(claim, "triplets", lambda: dt.groupby(["session", "row"]).ngroups)
        for gg in gs:
            x = dt[dt.sextile == gg].ch
            add(claim, "sextile %d mean change (Hz)" % gg, lambda: x.mean())
            add(claim, "sextile %d p" % gg, lambda: wp(x))
        if key == ("NREM", "REM"):
            add(claim, "sextiles lower at p < 0.05", lambda: sum(
                dt[dt.sextile == g_].ch.mean() < 0 and wp(dt[dt.sextile == g_].ch) < 0.05 for g_ in range(1, 7)))
    dr = trip[("NREM", "REM")][trip[("NREM", "REM")].session.isin(ids)]
    per_ = dr.groupby(["session", "row"]).agg(pre=("pre", "mean"), post=("post", "mean"), dur=("mid_dur", "first"))
    per_ = per_[per_.pre > 0]
    add("fig5c-rem-duration", "R per triplet", lambda: st.pearsonr(per_.dur, (per_.post - per_.pre) / per_.pre)[0])
    add("fig5c-rem-duration", "p per triplet", lambda: st.pearsonr(per_.dur, (per_.post - per_.pre) / per_.pre)[1])
    e1 = d1[d1.session.isin(ids)]
    add("fig5c-sigma-ma", "R", lambda: st.pearsonr(e1.sigma, e1.ch)[0])
    add("fig5c-sigma-ma", "p", lambda: st.pearsonr(e1.sigma, e1.ch)[1])
    e2 = dd[dd.session.isin(ids)]
    add("fig5c-delta-change", "R", lambda: st.pearsonr(e2.d_delta, e2.d_s1)[0])
    add("fig5c-delta-change", "p", lambda: st.pearsonr(e2.d_delta, e2.d_s1)[1])
    for name_, cs_all in (("WAKE", wk), ("REM", rm)):
        try:
            tt = sextile_table([c_ for c_ in cs_all if c_[2].session_id in ids])
        except ValueError:
            tt = {}
        if name_ == "WAKE":
            for gg in (1, 2, 5, 6):
                add("figs5-within-wake-rem", "WAKE sextile %d R" % gg, lambda: tt[gg][0])
                add("figs5-within-wake-rem", "WAKE sextile %d p" % gg, lambda: tt[gg][1])
        else:
            add("figs5-within-wake-rem", "REM sextiles higher at p < 0.05",
                lambda: sum(tt[g_][0] > 0 and tt[g_][1] < 0.05 for g_ in range(1, 7)))
    return out

tables = {name: headlines(ids) for name, ids in subsets.items()}
first = next(iter(tables.values()))
sens = pd.DataFrame({name: [v for _, _, v in rows_] for name, rows_ in tables.items()},
                    index=pd.MultiIndex.from_tuples([(c_, q_) for c_, q_, _ in first], names=["claim", "quantity"]))
pos = {c_: i for i, c_ in enumerate(dict.fromkeys(c_ for c_, _, _ in first))}
sens = sens.iloc[sorted(range(len(first)), key=lambda i: pos[first[i][0]])]
print(sens.to_string(float_format=lambda v: "%.4g" % v))

Labels not confirmed: ['BWRat21_121813'] | LFP not confirmed: ['BWRat17_121912', 'BWRat21_121113', 'BWRat21_121613', 'BWRat21_121813', 'Bogey_012615', 'Dino_061914', 'Dino_072414']

                                                                       all 21  without BWRat21_121813  without 7 (LFP)
claim                          quantity                                                                               
dataset-composition            pE cells                                   812                     801              590
                               pI cells                                    89                      88               68
fig2a-median-rates             pE WAKE median (Hz)                     0.7545                  0.7552           0.8101
                               pE NREM median (Hz)                      0.662                  0.6595           0.6639
                               pE REM median (Hz)                      0.8246                  0.8199           0.8926
                               pI WAKE median (Hz)                      3.917                   3.957            4.369
                               pI NREM median (Hz)                       2.73                   2.747            2.906
                               pI REM median (Hz)                       3.062                   3.078            3.121
fig2a-skewed-wide-distribution WAKE skewness                            4.575                   4.545            3.567
                               WAKE 2.5th percentile (Hz)              0.0124                 0.01226          0.02249
                               WAKE 97.5th percentile (Hz)              7.165                   7.266            7.538
fig2a-fraction-above-2hz       % above 2 Hz, WAKE                          17                   17.23            20.17
                               % above 2 Hz, REM                        18.72                   18.73            20.34
                               % above 2 Hz, NREM                       9.113                   9.238            10.85
                               chi-square p, NREM against WAKE      3.479e-06               3.399e-06        1.407e-05
                               chi-square p, NREM against REM       3.381e-08               6.549e-08        1.018e-05
fig2b-wake-nrem-slope          R, log10                                0.8623                   0.863           0.8693
                               slope, log10                            0.7574                  0.7562           0.7804
                               slope, upper 95% bound                  0.7884                  0.7872           0.8164
fig2-burst-index-by-state      median % NREM (a)                        6.353                   6.378            6.477
                               median % WAKE (a)                        2.645                   2.682             2.78
                               median % REM (a)                         3.649                   3.689             3.96
                               median % NREM (b)                        11.67                    11.7            11.99
                               median % WAKE (b)                        4.953                   4.999            5.278
                               median % REM (b)                         6.707                   6.754            7.271
                               Kruskal-Wallis p (a)                  5.74e-36               3.544e-35        7.821e-23
fig2c-acg-peaks                NREM peak 2-50 ms (ms)                     8.5                     8.5              6.5
                               REM smoothed maximum 60-250 ms (ms)       60.5                    60.5             60.5
                               WAKE smoothed maximum 60-250 ms (ms)      61.5                    60.5             61.5
fig2f-rate-group-correlations  (6,6) mean r                            0.0424                  0.0424          0.04278
                               (6,6) the largest (1 if so)                  1                       1                1
                               (1,6) mean r                          0.008865                0.008865         0.008945
ma-statistics                  median duration (s)                         19                      19               17
                               per min of SLEEP                        0.0962                 0.09928          0.09968
                               % of SLEEP                               3.223                   3.321            3.163
ma-fraction-over-nrem-episodes R of bin means                         -0.4514                 -0.4255          -0.5174
                               R pooled                              -0.08713                 -0.0833          -0.0976
packet-duration                median (s)                                 108                     106            106.5
wake-sleep-cycles              cycles per session                       1.429                     1.4            1.357
fig3b-population-decline       R (a)                                  -0.2675                 -0.2867          -0.1972
                               R (b)                                  -0.2852                 -0.2793          -0.2173
fig3-first-last-packet         first mean (Hz)                          1.091                   1.098            1.134
                               last mean (Hz)                          0.9679                   0.971            1.006
                               Wilcoxon p                           3.918e-06               2.222e-06         0.001708
fig3d-first-last-slope         slope, log10                            0.8007                  0.8027           0.8002
                               slope, upper 95% bound                  0.8336                  0.8355           0.8381
fig3e-wake-before-after        Wilcoxon p, all pE                   4.079e-09               9.593e-10        2.447e-09
                               sextile 1 one-tailed p                  0.3039                   0.283           0.6236
                               sextile 4 one-tailed p                 0.03003                  0.0114          0.02323
                               sextile 5 one-tailed p               0.0006272               0.0006272         0.002218
                               sextile 6 one-tailed p               2.411e-05               2.411e-05        0.0001242
fig3b-sextiles-over-sleep      sextile 1 R (a)                         0.1633                  0.1464           0.1689
                               sextile 1 p (a)                       0.007391                 0.02162          0.04605
                               sextile 6 R (a)                        -0.4475                 -0.4475          -0.3673
                               sextile 6 p (a)                      3.066e-11               3.066e-11        8.042e-06
fig3b-cv-over-sleep            R (a)                                 -0.02627                -0.03836         -0.07173
                               p (a)                                   0.6576                  0.5341           0.3689
fig4iv-sextiles-within-nrem    episode, sextile 1 R                   0.06033                  0.0678            0.066
                               episode, sextile 1 p                 6.266e-06               9.033e-07        0.0001252
                               episode, sextile 6 R                  -0.08264                -0.08264         -0.09016
                               episode, sextile 6 p                  3.74e-08                3.74e-08        8.055e-08
                               packet, sextile 1 R                    0.06062                  0.0664          0.06133
                               packet, sextile 1 p                  2.032e-08               2.315e-09        1.063e-05
                               packet, sextile 6 R                   -0.05539                -0.05539          -0.0478
                               packet, sextile 6 p                  4.421e-06               4.421e-06        0.0004413
fig4v-cv-within-nrem           episode R                            -0.004703               -0.005099        -0.003879
                               episode p                               0.7123                  0.6981            0.808
                               packet R                               0.01108                 0.01145          0.01501
                               packet p                                0.2797                  0.2763           0.2403
fig4vi-burst-within-nrem       episode R                              0.02924                 0.03312          0.03612
                               episode p                              0.02182                 0.01172          0.02361
                               packet R                                0.0237                 0.02688          0.03441
                               packet p                               0.02076                 0.01055         0.007093
fig4i-band-power-within-nrem   episode, delta R                       0.06909                 0.07376          0.06534
                               episode, delta p                     1.139e-07               4.098e-08        4.248e-05
                               episode, sigma R                        0.1348                  0.1301           0.1364
                               episode, sigma p                     3.114e-25               2.851e-22        9.578e-18
                               episode, gamma R                       0.06233                 0.04448         0.004732
                               episode, gamma p                     1.723e-06               0.0009465           0.7671
                               packet, delta R                      -0.005331               -0.002496          0.01425
                               packet, delta p                         0.6101                  0.8162           0.2651
                               packet, sigma R                         0.1042                  0.1027           0.1054
                               packet, sigma p                      1.569e-23               8.461e-22        1.351e-16
                               packet, gamma R                         0.1821                  0.1708           0.1165
                               packet, gamma p                      4.454e-69               7.593e-58        6.019e-20
fig4-delta-packet-length       R                                        0.334                  0.3263           0.2965
                               p                                    2.825e-25               5.544e-23        6.897e-14
fig4cd-packet-rate-lfp         sextile 6 against delta R              -0.4165                 -0.4165          -0.4282
                               sextile 6 against gamma R               0.5697                  0.5697           0.5326
                               sextile 1 against last-20-s sigma R     0.1038                  0.1128           0.1323
fig5b-nrem-between-ma          triplets                                   148                     146              105
                               sextile 5 mean change (Hz)            -0.01074                -0.01074         -0.02236
                               sextile 5 p                             0.2362                  0.2362           0.1599
                               sextile 6 mean change (Hz)             -0.1037                 -0.1037         -0.09586
                               sextile 6 p                            0.07661                 0.07661           0.1473
fig5b-wake-between-packets     triplets                                   265                     258              172
                               sextile 6 mean change (Hz)            -0.05092                -0.05092         -0.03679
                               sextile 6 p                            0.09198                 0.09198           0.1984
fig5b-ma-between-packets       triplets                                   391                     380              266
                               sextile 1 mean change (Hz)             0.00339                0.003336         0.001232
                               sextile 1 p                            0.08198                 0.09945           0.4859
                               sextile 6 mean change (Hz)            -0.07253                -0.07253         -0.06828
                               sextile 6 p                           6.76e-05                6.76e-05        0.0009182
fig5b-rem-between-packets      triplets                                    17                      13               13
                               sextiles lower at p < 0.05                   2                       2                2
fig5c-rem-duration             R per triplet                         -0.04663                 -0.4519          -0.4519
                               p per triplet                           0.8589                  0.1211           0.1211
fig5c-sigma-ma                 R                                      0.03319                 0.03445       -0.0005352
                               p                                       0.7133                  0.7052           0.9961
fig5c-delta-change             R                                       0.1828                  0.1828           0.1136
                               p                                    5.261e-06               5.261e-06          0.01052
figs5-within-wake-rem          WAKE sextile 1 R                       0.02084                  0.0278          0.03253
                               WAKE sextile 1 p                        0.1623                 0.06837          0.09163
                               WAKE sextile 2 R                       0.06682                 0.06029          0.08578
                               WAKE sextile 2 p                     5.225e-06               5.818e-05        1.481e-06
                               WAKE sextile 5 R                       -0.1219                 -0.1219         -0.06031
                               WAKE sextile 5 p                     3.504e-17               3.504e-17         0.000722
                               WAKE sextile 6 R                       -0.1202                 -0.1202          -0.1227
                               WAKE sextile 6 p                     4.941e-13               4.941e-13        6.292e-11
                               REM sextiles higher at p < 0.05              0                       1                0

loose = {s.path: ss.classify_states(s.starts, s.stops, s.labels, ma_max=60.0, both_sides=False) for s in sessions}
print("MAs: %d under the rule used here, %d under the looser rule" % (
    sum(int(np.sum(s.state == "MA")) for s in sessions), sum(int(np.sum(v == "MA")) for v in loose.values())))
# ma-fraction-over-nrem-episodes
fr2 = []
for s in sessions:
    lab = loose[s.path]
    for a, b in ss.nrem_episodes(lab):
        mas = [(s.starts[i], s.stops[i]) for i in range(a, b + 1) if lab[i] == "MA"]
        if mas:
            width = (s.stops[b] - s.starts[a]) / 10
            fr2.append([np.sum(e - st_) / width for st_, e in ss.binned_intervals(s.starts[a], s.stops[b], 10, mas)])
fr2 = np.array(fr2)
print("\nma-fraction-over-nrem-episodes: %d episodes with an MA; bin means R = %.3f (p = %.2g); pooled R = %.3f (p = %.2g)" % (
    len(fr2), *st.pearsonr(np.arange(10) + 0.5, fr2.mean(axis=0)), *st.pearsonr(np.tile(np.arange(10) + 0.5, len(fr2)), fr2.ravel())))
print("   mean MA fraction per bin:", np.round(fr2.mean(axis=0), 4))
trip_l = {k: triplets(*k, states=loose) for k in (("MA", "NREM"), ("NREM", "MA"))}
for key, claim, gs in ((("MA", "NREM"), "fig5b-nrem-between-ma", (5, 6)), (("NREM", "MA"), "fig5b-ma-between-packets", (1, 5, 6))):
    dt = trip_l[key]
    print("\n%s: triplets %d (rule used here: %d)" % (claim, dt.groupby(["session", "row"]).ngroups,
                                                    trip[key].groupby(["session", "row"]).ngroups))
    for g in gs:
        x = dt[dt.sextile == g]
        t_ = x.groupby(["session", "row"]).ch.mean()
        print("   sextile %d: mean change %+.4f Hz, Wilcoxon p = %.2g | one value per triplet p = %.2g, per session p = %.2g" % (
            g, x.ch.mean(), wp(x.ch), wp(t_), wp(x.groupby("session").ch.mean())))
# fig5c-sigma-ma
e = trip_l[("MA", "NREM")]
e = e[e.sextile == 1].groupby(["session", "row"]).agg(ch=("ch", "mean")).reset_index()
e["sigma"] = [sig_all.get(path_of[a], {}).get(b, np.nan) for a, b in zip(e.session, e.row)]
e = e.dropna()
print("\nfig5c-sigma-ma: %d triplets with sextile 1 cells and LFP; R = %.3f, p = %.2g" % (len(e), *st.pearsonr(e.sigma, e.ch)))

MAs: 391 under the rule used here, 757 under the looser rule

ma-fraction-over-nrem-episodes: 249 episodes with an MA; bin means R = -0.602 (p = 0.066); pooled R = -0.111 (p = 2.8e-08)
   mean MA fraction per bin: [0.0947 0.152  0.1513 0.1494 0.1406 0.1603 0.1422 0.0934 0.0785 0.0138]

fig5b-nrem-between-ma: triplets 343 (rule used here: 148)
   sextile 5: mean change -0.0000 Hz, Wilcoxon p = 0.14 | one value per triplet p = 0.84, per session p = 0.52
   sextile 6: mean change +0.0286 Hz, Wilcoxon p = 0.92 | one value per triplet p = 0.4, per session p = 0.71

fig5b-ma-between-packets: triplets 463 (rule used here: 391)
   sextile 1: mean change +0.0048 Hz, Wilcoxon p = 0.016 | one value per triplet p = 0.029, per session p = 0.11
   sextile 5: mean change -0.0262 Hz, Wilcoxon p = 0.0041 | one value per triplet p = 0.019, per session p = 0.024
   sextile 6: mean change -0.0819 Hz, Wilcoxon p = 5.3e-06 | one value per triplet p = 0.00013, per session p = 3.1e-05

fig5c-sigma-ma: 293 triplets with sextile 1 cells and LFP; R = 0.110, p = 0.06

Claim	Tested or blocked
dataset-composition	tested
fig2a-median-rates	tested
fig2a-skewed-wide-distribution	tested
fig2a-fraction-above-2hz	tested
fig2b-wake-nrem-slope	tested
fig2-burst-index-by-state	tested
fig2c-acg-peaks	tested
fig2c-acg-decay	blocked: methods underspecified
fig2f-rate-group-correlations	tested
ma-statistics	tested
ma-fraction-over-nrem-episodes	tested
packet-duration	tested
wake-sleep-cycles	tested
fig3b-population-decline	tested
fig3-first-last-packet	tested
fig3d-first-last-slope	tested
fig3e-wake-before-after	tested
fig3b-sextiles-over-sleep	tested
fig3b-cv-over-sleep	tested
fig4iv-sextiles-within-nrem	tested
fig4v-cv-within-nrem	tested
fig4vi-burst-within-nrem	tested
fig4i-band-power-within-nrem	tested
fig4-delta-packet-length	tested
fig4cd-packet-rate-lfp	tested
fig5b-nrem-between-ma	tested
fig5b-wake-between-packets	tested
fig5b-ma-between-packets	tested
fig5b-rem-between-packets	tested
fig5c-rem-duration	tested
fig5c-sigma-ma	tested
fig5c-delta-change	tested
figs5-within-wake-rem	tested
fig5c-wake-movement	blocked: data missing
updown-durations	blocked: methods underspecified
fig3c-slow-oscillation-spindles-over-sleep	blocked: methods underspecified
fig4ii-iii-updown-spindles-within-packets	blocked: methods underspecified
fig2de-up-state-sequences	blocked: methods underspecified
fig3f-4e-within-up-rates	blocked: methods underspecified

Claim	Outcome
dataset-composition	not reproduced
fig2a-median-rates	partially reproduced
fig2a-skewed-wide-distribution	partially reproduced
fig2a-fraction-above-2hz	not reproduced
fig2b-wake-nrem-slope	reproduced
fig2-burst-index-by-state	partially reproduced
fig2c-acg-peaks	partially reproduced
fig2c-acg-decay	blocked (methods underspecified)
fig2f-rate-group-correlations	partially reproduced
ma-statistics	not reproduced
ma-fraction-over-nrem-episodes	partially reproduced
packet-duration	reproduced
wake-sleep-cycles	not reproduced
fig3b-population-decline	partially reproduced
fig3-first-last-packet	reproduced
fig3d-first-last-slope	reproduced
fig3e-wake-before-after	partially reproduced
fig3b-sextiles-over-sleep	reproduced
fig3b-cv-over-sleep	not reproduced
fig4iv-sextiles-within-nrem	reproduced
fig4v-cv-within-nrem	not reproduced
fig4vi-burst-within-nrem	not reproduced
fig4i-band-power-within-nrem	not reproduced
fig4-delta-packet-length	reproduced
fig4cd-packet-rate-lfp	partially reproduced
fig5b-nrem-between-ma	not reproduced
fig5b-wake-between-packets	not reproduced
fig5b-ma-between-packets	not reproduced
fig5b-rem-between-packets	not reproduced
fig5c-rem-duration	not reproduced
fig5c-sigma-ma	not reproduced
fig5c-delta-change	reproduced
figs5-within-wake-rem	not reproduced
fig5c-wake-movement	blocked (data missing)
updown-durations	blocked (methods underspecified)
fig3c-slow-oscillation-spindles-over-sleep	blocked (methods underspecified)
fig4ii-iii-updown-spindles-within-packets	blocked (methods underspecified)
fig2de-up-state-sequences	blocked (methods underspecified)
fig3f-4e-within-up-rates	blocked (methods underspecified)

Reproducing Watson et al. (2016) from dandiset 000041¶

Data check¶

dataset-composition¶

fig2a-median-rates¶

fig2a-skewed-wide-distribution¶

fig2a-fraction-above-2hz¶

fig2b-wake-nrem-slope¶

fig2-burst-index-by-state¶

fig2c-acg-peaks¶

fig2f-rate-group-correlations¶

ma-statistics¶

ma-fraction-over-nrem-episodes¶

packet-duration¶

wake-sleep-cycles¶

fig3b-population-decline¶

fig3-first-last-packet¶

fig3d-first-last-slope¶

fig3e-wake-before-after¶

fig3b-sextiles-over-sleep¶

fig3b-cv-over-sleep¶

Time courses within nonREM packets and episodes¶

fig4iv-sextiles-within-nrem¶

fig4v-cv-within-nrem¶

fig4vi-burst-within-nrem¶

fig4i-band-power-within-nrem¶

fig4-delta-packet-length¶

fig4cd-packet-rate-lfp¶

State triplets¶

fig5b-nrem-between-ma¶

fig5b-wake-between-packets¶

fig5b-ma-between-packets¶

fig5b-rem-between-packets¶

fig5c-rem-duration¶

fig5c-sigma-ma¶

fig5c-delta-change¶

figs5-within-wake-rem¶

Sessions whose clocks could not be confirmed¶

Microarousal-based claims under a looser MA rule¶

Summary¶

Reproducing Watson et al. (2016) from dandiset 000041¶

Data check¶

Do the spikes, the states and the LFP share one clock?¶

dataset-composition¶

fig2a-median-rates¶

fig2a-skewed-wide-distribution¶

fig2a-fraction-above-2hz¶

fig2b-wake-nrem-slope¶

fig2-burst-index-by-state¶

fig2c-acg-peaks¶

fig2f-rate-group-correlations¶

ma-statistics¶

ma-fraction-over-nrem-episodes¶

packet-duration¶

wake-sleep-cycles¶

fig3b-population-decline¶

fig3-first-last-packet¶

fig3d-first-last-slope¶

fig3e-wake-before-after¶

fig3b-sextiles-over-sleep¶

fig3b-cv-over-sleep¶

Time courses within nonREM packets and episodes¶

fig4iv-sextiles-within-nrem¶

fig4v-cv-within-nrem¶

fig4vi-burst-within-nrem¶

fig4i-band-power-within-nrem¶

fig4-delta-packet-length¶

fig4cd-packet-rate-lfp¶

State triplets¶

fig5b-nrem-between-ma¶

fig5b-wake-between-packets¶

fig5b-ma-between-packets¶

fig5b-rem-between-packets¶

fig5c-rem-duration¶

fig5c-sigma-ma¶

fig5c-delta-change¶

figs5-within-wake-rem¶

Sessions whose clocks could not be confirmed¶

Microarousal-based claims under a looser MA rule¶

Summary¶