import os, re, time, resource
import numpy as np
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt
from multiprocessing import Pool
import fcx

t_start = time.time()
md5 = fcx.checksums()
names = fcx.sessions()
print(f"CRCNS {fcx.DATASET} (doi:{fcx.DOI}), version {fcx.VERSION}; {len(names)} session files")
for n in names:
    print(f"  {md5['data/' + n + '.tar.gz']}  data/{n}.tar.gz")

CRCNS fcx-1 (doi:10.6080/K02N506Q), version unversioned; 27 session files
  76a3a51748d22a4b8f822b59c5367c62  data/20140526_277um.tar.gz
  c9ad3be105258bc5e70b34d50c7964cd  data/20140527_421um.tar.gz
  0342fc85311c394b52a8d29a1290a39f  data/20140528_565um.tar.gz
  5bae671bcfd9f76524b08b9dcbcf3d0e  data/BWRat17_121712.tar.gz
  c744f4a16327b474563f0f615b05e343  data/BWRat17_121912.tar.gz
  f2c6a9a725eb929f1482c48ca0f816ab  data/BWRat18_020513.tar.gz
  908398356c84000b7dce814c3c408dd2  data/BWRat19_032413.tar.gz
  4169ce26875aaaec4ef42c3ab8de402c  data/BWRat19_032513.tar.gz
  65e3bc3f961c6a9b3725e96a2682ed1d  data/BWRat20_101013.tar.gz
  5f1dcb5cc0d2972b6bff9b020a4f26aa  data/BWRat20_101513.tar.gz
  0fc5726ba2ce996996c085f37b02597e  data/BWRat21_121113.tar.gz
  8e2c5d14ee822aaf876cff3b8dff5510  data/BWRat21_121613.tar.gz
  021b1a96327fcc82c8b81d31cd086873  data/BWRat21_121813.tar.gz
  9a65bdc269ec5df500236b4e4b3e79d4  data/Bogey_012615.tar.gz
  4ac889b2421e4c32ccd80c67d96df793  data/Dino_061814_mPFC.tar.gz
  b916aef6b88f570cf91533e4d70c29ed  data/Dino_061914_ACC.tar.gz
  098ce2e0340ef414fad7d56a75fbc946  data/Dino_061914_mPFC.tar.gz
  9225d5b807be4c2b7a0634dab53cc64f  data/Dino_062014_ACC.tar.gz
  3f6d06564687a5be5f53be4f8ec8dbff  data/Dino_062014_mPFC.tar.gz
  9d1bf8b88852e2e7f1e8ea749d2b9199  data/Dino_072114_mPFC.tar.gz
  3bccdfe803cd415593ce8519b6e000dd  data/Dino_072314_mPFC.tar.gz
  d0b8faaa90589c10bebe47492b66b676  data/Dino_072414_mPFC.tar.gz
  4deca5619586d6f7b21110005894e125  data/Rizzo_022615.tar.gz
  acc3d4a3bc2f40172d11aa0f4c344df7  data/Rizzo_022715.tar.gz
  b6b5c4629f528edd74f14af178b32ca8  data/Splinter_020515.tar.gz
  e62faa965caadabedebd6af8beda6b0e  data/Splinter_020915.tar.gz
  9c1e4255a4437a9fe68c3b9e44b221a8  data/Templeton_032415.tar.gz

# Stream every session once (about 90 GB of tar.gz, decompressed in memory and
# reduced as it is read; see fcx.py). Four sessions at a time.
if os.environ.get("FCX_DEV"):
    import pickle
    S = pickle.load(open(os.environ["FCX_DEV"], "rb"))
else:
    with Pool(4) as pool:
        S = pool.map(fcx.load, names, chunksize=1)
S = sorted(S, key=lambda s: s["session"])
tabs = [fcx.seconds_table(s) for s in S]
for t in tabs:
    fcx.add_zband(t)
print(f"loaded {len(S)} sessions in {(time.time() - t_start) / 60:.1f} min; "
      f"peak memory of this process {resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1e6:.1f} GB")

loaded 27 sessions in 8.4 min; peak memory of this process 1.3 GB

rows = []
for s, tab in zip(S, tabs):
    st = s["states"]
    cyc = fcx.cycles(s)
    dur = s["eeg_size"] / 2 / s["meta"]["nChannels"] / s["meta"]["lfp_fs"]
    allspk = np.concatenate([x for x in s["spikes"] if len(x)])
    rows.append(dict(session=s["session"], rat=s["rat"], hours=dur / 3600, units=len(s["spikes"]),
                     pE=len(s["E"]), pI=len(s["I"]), unlabelled=len(s["spikes"]) - len(s["E"]) - len(s["I"]),
                     cycles=len(cyc), wake_min=np.mean([c["wake"][1] - c["wake"][0] for c in cyc]) / 60,
                     sleep_min=np.mean([c["sleep"][1] - c["sleep"][0] for c in cyc]) / 60,
                     packets=len(fcx.state_intervals(s, "SWSPacket")), MA=len(fcx.state_intervals(s, "MA")),
                     REM=len(fcx.state_intervals(s, "REM")), UP=len(s["updown"]["UP"]), DOWN=len(s["updown"]["DN"]),
                     spindles=len(s["spindles"]["times"]), spikes_min_s=allspk.min(), spikes_max_s=allspk.max(),
                     lfp_ch=s["meta"]["goodeegchannel"], spindle_ch=s["spindles"]["channel"],
                     up_ch=s["meta"]["UPstatechannel"], eeg_second_pass=s["second_pass"]))
ov = pd.DataFrame(rows)
pd.set_option("display.width", 250)
pd.set_option("display.max_columns", 40)
print(ov.round(2).to_string())
print("\ntotals:", ov[["units", "pE", "pI", "cycles", "packets", "MA", "REM", "UP", "DOWN", "spindles"]].sum().to_dict())
print("rats:", ov.rat.nunique(), "sessions:", len(ov))
print("sessions per rat:", ov.groupby("rat").size().to_dict())

             session            rat  hours  units   pE  pI  unlabelled  cycles  wake_min  sleep_min  packets  MA  REM    UP   DOWN  spindles  spikes_min_s  spikes_max_s  lfp_ch  spindle_ch  up_ch  eeg_second_pass
0     20140526_277um  JennBuzsaki22   3.77     62   51  11           0       2     17.84      33.31       26  22    6   982   1190       157          0.01      13572.34      33          33     33             True
1     20140527_421um  JennBuzsaki22   4.89     29   24   5           0       1     27.93      39.67        8   4    5   913   1550        92          0.01      17615.25       9           9      9             True
2     20140528_565um  JennBuzsaki22   3.44     33   27   6           0       1     13.15      62.18       15   8   14  1237   1664       158          0.01      12395.73       9           9      9            False
3     BWRat17_121712        BWRat17   1.68     50   43   7           0       1     40.43      26.30        9   6    1   579    995        46          0.04       6059.36      13          13     13            False
4     BWRat17_121912        BWRat17   1.82     35   34   1           0       1     19.73      70.65       33  27    5   885   2269       136          0.02       6534.69      13          13     13             True
5     BWRat18_020513        BWRat18   2.11     40   34   6           0       2     17.52      44.32       17   2   13  1037   3886       172          0.03       7590.45      59          59     59            False
6     BWRat19_032413        BWRat19   1.39     50   37  13           0       1     33.35      40.35       19  20    0  1246   1592        73          0.01       5013.39      54          54     54             True
7     BWRat19_032513        BWRat19   1.09     61   51  10           0       1      7.55      27.12        6   3    9   654    822        56          0.01       3928.79      24          24     24             True
8     BWRat20_101013        BWRat20   3.71     22   17   5           0       1     11.75      54.55       18  13    9  1380   2050        89          0.00      13364.29      22          59     22             True
9     BWRat20_101513        BWRat20   6.39     66   60   6           0       2     10.26      64.37       35  28   27  2917   4282       243          0.01      23011.99      61          61     61             True
10    BWRat21_121113        BWRat21   6.08     32   28   4           0       2     17.93      51.85       28  22   21  1194   2894       124          0.00      21883.19      14           7     14            False
11    BWRat21_121613        BWRat21   5.98     21   20   1           0       1      7.82      32.23        6   1    2    71   1679        66          0.13      21541.79      22          12     22            False
12    BWRat21_121813        BWRat21   6.03     12   11   1           0       3     36.77      71.72       49  17   67    54   2549       147          0.45      21705.63      22          12     22             True
13      Bogey_012615          Bogey   6.74     36   33   3           0       3     20.24      54.30       60  34   18  4022   6364       298          0.01      24263.40      46          46     46            False
14  Dino_061814_mPFC           Dino   5.46    113  104   9           0       2     34.60      41.44       29  23   17  1003   1178       144          0.00      19646.17      69          92     69             True
15   Dino_061914_ACC           Dino   4.87     52   46   6           0       1     10.50      37.25       11   5    5  1030   1376        82          0.06      17519.53      19          19     19            False
16  Dino_061914_mPFC           Dino   4.87     60   55   5           0       1     10.50      37.25       12   5    5   842   1096        62          0.00      17519.54      83          83     83            False
17   Dino_062014_ACC           Dino   6.21     24   23   1           0       3     17.24      59.46       41  20   21  3417   7364       379          0.10      22363.13      18          14     18             True
18  Dino_062014_mPFC           Dino   6.21     36   32   4           0       3     17.33      27.11       17   8   10  1874   2656       159          0.06      22363.20      83          83     83            False
19  Dino_072114_mPFC           Dino   2.21     18   18   0           0       1     55.93      27.18        3   4    5   121    651        57          0.00       7958.64     120         106    120             True
20  Dino_072314_mPFC           Dino   2.87     16   16   0           0       2     19.64      24.45       16   9   11   320   1167       100          0.06      10321.09     120         106    120            False
21  Dino_072414_mPFC           Dino   3.67     46   40   6           0       1     54.95      25.80        6   9    1   348    446        50          0.00      13208.80     120         106    120             True
22      Rizzo_022615          Rizzo  11.27     26   25   1           0       5     24.32      64.74       66  14   52  6449  10322       879          0.03      40555.10      60          17      1            False
23      Rizzo_022715          Rizzo   7.20     93   89   4           0       2     41.81      82.52       44  31   25  4404   5769       306          0.00      25932.10      40          52     17            False
24   Splinter_020515       Splinter   6.62     27   25   2           0       3     25.59      47.19       33  28   13  2046   4691       326          0.02      23836.56      42          42     42             True
25   Splinter_020915       Splinter   6.52     51   42   9           0       2     32.82      88.28       31  19   27  4161   7803       349          0.00      23487.50      43          43     43             True
26  Templeton_032415      Templeton   5.29     10   10   0           0       3     25.19      73.44       63  52   28    44   8534       453          0.10      19057.56      27          27     27            False

totals: {'units': 1121, 'pE': 995, 'pI': 126, 'cycles': 51, 'packets': 701, 'MA': 434, 'REM': 417, 'UP': 43230, 'DOWN': 86839, 'spindles': 5203}
rats: 11 sessions: 27
sessions per rat: {'BWRat17': 2, 'BWRat18': 1, 'BWRat19': 2, 'BWRat20': 2, 'BWRat21': 3, 'Bogey': 1, 'Dino': 8, 'JennBuzsaki22': 3, 'Rizzo': 2, 'Splinter': 2, 'Templeton': 1}

# Every session's scored WAKE, SLEEP and WAKE-SLEEP intervals: whether the 51 cycles are all the files
# hold (the paper reports 54), and cycles whose WAKE starts before the previous SLEEP ends.
tot = dict(WakeSleep=0, Wake=0, Sleep=0, unpaired_Sleep=0, unpaired_Wake=0)
for s in S:
    st = s["states"]
    ws = [(tuple(c[0]), tuple(c[1])) for c in st["WakeSleep"]]
    unp_s = [tuple(x) for x in st["Sleep"] if tuple(x) not in [w[1] for w in ws]]
    unp_w = [tuple(x) for x in st["Wake"] if tuple(x) not in [w[0] for w in ws]]
    for k, v in (("WakeSleep", len(ws)), ("Wake", len(st["Wake"])), ("Sleep", len(st["Sleep"])),
                 ("unpaired_Sleep", len(unp_s)), ("unpaired_Wake", len(unp_w))):
        tot[k] += v
    cyc = fcx.cycles(s)
    ov_ = [f"WAKE {int(d['wake'][0])}-{int(d['wake'][1])} s starts before SLEEP {int(c['sleep'][0])}-{int(c['sleep'][1])} s ends"
           for c, d in zip(cyc, cyc[1:]) if d["wake"][0] < c["sleep"][1]]
    print(f"{s['session']}: WakeSleep {len(ws)}, Wake {len(st['Wake'])}, Sleep {len(st['Sleep'])}, unpaired Sleep {len(unp_s)}, "
          f"unpaired Wake {len(unp_w)}" + (f"; {'; '.join(ov_)}" if ov_ else ""))
print("totals:", tot)
# Drug injections noted in the data set's WatsonSleepHomestasis2016Table.xlsx, against the recording
# files (RecordingFileIntervals of _BasicMetaData.mat) and the cycles.
Sd = {s["session"]: s for s in S}
for row in fcx.overview_table():
    notes = [str(v) for v in row[1:] if isinstance(v, str) and re.search("ketamine|saline", v, re.I)]
    if row and row[0] in Sd and notes:
        s = Sd[row[0]]
        rfi = s["meta"]["RecordingFileIntervals"]
        print(f"{row[0]}: table notes {notes}; {len(rfi)} recording files, the last ending at {rfi[-1, 1]:.0f} s; "
              f"cycles end at {[int(c['sleep'][1]) for c in fcx.cycles(s)]} s")

20140526_277um: WakeSleep 2, Wake 2, Sleep 2, unpaired Sleep 0, unpaired Wake 0
20140527_421um: WakeSleep 1, Wake 1, Sleep 1, unpaired Sleep 0, unpaired Wake 0
20140528_565um: WakeSleep 1, Wake 1, Sleep 1, unpaired Sleep 0, unpaired Wake 0
BWRat17_121712: WakeSleep 1, Wake 1, Sleep 1, unpaired Sleep 0, unpaired Wake 0
BWRat17_121912: WakeSleep 1, Wake 1, Sleep 1, unpaired Sleep 0, unpaired Wake 0
BWRat18_020513: WakeSleep 2, Wake 2, Sleep 2, unpaired Sleep 0, unpaired Wake 0
BWRat19_032413: WakeSleep 1, Wake 1, Sleep 1, unpaired Sleep 0, unpaired Wake 0
BWRat19_032513: WakeSleep 1, Wake 1, Sleep 1, unpaired Sleep 0, unpaired Wake 0
BWRat20_101013: WakeSleep 1, Wake 1, Sleep 1, unpaired Sleep 0, unpaired Wake 0
BWRat20_101513: WakeSleep 2, Wake 2, Sleep 2, unpaired Sleep 0, unpaired Wake 0
BWRat21_121113: WakeSleep 2, Wake 2, Sleep 2, unpaired Sleep 0, unpaired Wake 0
BWRat21_121613: WakeSleep 1, Wake 1, Sleep 1, unpaired Sleep 0, unpaired Wake 0
BWRat21_121813: WakeSleep 3, Wake 3, Sleep 3, unpaired Sleep 0, unpaired Wake 0
Bogey_012615: WakeSleep 3, Wake 3, Sleep 3, unpaired Sleep 0, unpaired Wake 0
Dino_061814_mPFC: WakeSleep 2, Wake 2, Sleep 2, unpaired Sleep 0, unpaired Wake 0
Dino_061914_ACC: WakeSleep 1, Wake 1, Sleep 1, unpaired Sleep 0, unpaired Wake 0
Dino_061914_mPFC: WakeSleep 1, Wake 1, Sleep 1, unpaired Sleep 0, unpaired Wake 0
Dino_062014_ACC: WakeSleep 3, Wake 3, Sleep 3, unpaired Sleep 0, unpaired Wake 0
Dino_062014_mPFC: WakeSleep 3, Wake 3, Sleep 3, unpaired Sleep 0, unpaired Wake 0
Dino_072114_mPFC: WakeSleep 1, Wake 1, Sleep 1, unpaired Sleep 0, unpaired Wake 0
Dino_072314_mPFC: WakeSleep 2, Wake 2, Sleep 2, unpaired Sleep 0, unpaired Wake 0
Dino_072414_mPFC: WakeSleep 1, Wake 1, Sleep 1, unpaired Sleep 0, unpaired Wake 0
Rizzo_022615: WakeSleep 5, Wake 5, Sleep 5, unpaired Sleep 0, unpaired Wake 0
Rizzo_022715: WakeSleep 2, Wake 2, Sleep 2, unpaired Sleep 0, unpaired Wake 0; WAKE 9919-14394 s starts before SLEEP 5655-10010 s ends
Splinter_020515: WakeSleep 3, Wake 3, Sleep 3, unpaired Sleep 0, unpaired Wake 0; WAKE 5335-7210 s starts before SLEEP 3543-5455 s ends
Splinter_020915: WakeSleep 2, Wake 2, Sleep 2, unpaired Sleep 0, unpaired Wake 0
Templeton_032415: WakeSleep 3, Wake 3, Sleep 3, unpaired Sleep 0, unpaired Wake 0
totals: {'WakeSleep': 51, 'Wake': 51, 'Sleep': 51, 'unpaired_Sleep': 0, 'unpaired_Wake': 0}

BWRat20_101513: table notes ['ketamine btw 03&04']; 3 recording files, the last ending at 23348 s; cycles end at [16042, 17896] s
BWRat21_121613: table notes ['ketamine btw 06&07', 'Ketamine']; 6 recording files, the last ending at 21542 s; cycles end at [7781] s
BWRat21_121813: table notes ['saline btw 05&06']; 5 recording files, the last ending at 21706 s; cycles end at [7190, 15911, 21705] s

/home/runner/work/_temp/executed/.venv/lib/python3.13/site-packages/openpyxl/worksheet/_reader.py:329: UserWarning: Unknown extension is not supported and will be removed
  warn(msg)

# Checks of the description against the values.
for s in S:
    md = s["meta"]
    assert md["lfp_fs"] == 1250
    # every pE/pI index points at a stable unit, and no unit is both
    assert set(s["E"]) | set(s["I"]) <= set(range(len(s["spikes"])))
    assert not set(s["E"]) & set(s["I"])
    # spike trains are sorted
    assert all(np.all(np.diff(x) >= 0) for x in s["spikes"])
print("LFP rate 1250 Hz in all sessions; pE/pI indices valid and disjoint; spike times sorted")
# state intervals: whole seconds, positive length
for k in ("SWSPacket", "SWSEpisode", "MA", "REM", "WakeInterruption", "Wake", "Sleep"):
    iv = np.vstack([s["states"][k] for s in S if len(s["states"][k])])
    print(f"{k}: {len(iv)} intervals, non-integer bounds {int(np.sum(iv != np.round(iv)))}, "
          f"zero or negative length {int(np.sum(iv[:, 1] <= iv[:, 0]))}")

# WAKE-SLEEP pairs: WAKE ends where SLEEP starts; WAKE >= 7 min and SLEEP >= 20 min (Table S1)
ws = [(s["session"], c["wake"][1] - c["wake"][0], c["sleep"][1] - c["sleep"][0], c["sleep"][0] - c["wake"][1])
      for s in S for c in fcx.cycles(s)]
ws = pd.DataFrame(ws, columns=["session", "wake_s", "sleep_s", "gap_s"])
print(f"WAKE-SLEEP cycles {len(ws)}; WAKE {ws.wake_s.min():.0f} to {ws.wake_s.max():.0f} s; "
      f"SLEEP {ws.sleep_s.min():.0f} to {ws.sleep_s.max():.0f} s; gap WAKE end to SLEEP start {ws.gap_s.min():.0f} to {ws.gap_s.max():.0f} s")
print("WAKE intervals equal the cycles' WAKE:", all(np.array_equal(s["states"]["Wake"], np.array([c["wake"] for c in fcx.cycles(s)])) for s in S))

# Durations of scored states against their definitions (Table S1, data description)
for key, lo, hi in (("SWSPacket", 20, None), ("MA", None, 40), ("WakeInterruption", 40, 120)):
    d = np.concatenate([np.diff(fcx.state_intervals(s, key), axis=1).ravel() for s in S])
    out = (d < lo).sum() if lo else (d > hi).sum()
    print(f"{key}: n {len(d)}, duration {d.min():.0f} to {d.max():.0f} s; outside definition {out}")
# MAs outside SLEEP
ma_all = sum(len(s["states"]["MA"]) for s in S)
ma_in = sum(len(fcx.state_intervals(s, "MA")) for s in S)
print(f"MA intervals in the files {ma_all}, inside SLEEP of a cycle {ma_in}")

# Units: the spike times span the recording; LFP length vs spikes
print("spike times within the LFP duration in every session:",
      bool(np.all(ov.spikes_max_s <= ov.hours * 3600 + 1)))
# sessions recorded at the same time on two probes (same rat and date)
import re
day = ov.session.str.replace(r"_(ACC|mPFC)$", "", regex=True)
print("sessions recorded at the same time on two probes:",
      [list(g) for _, g in ov.session.groupby(day) if len(g) > 1])
# duplicate spike trains across sessions (identical first 20 spike times)
sig = {}
dups = []
for s in S:
    for i, x in enumerate(s["spikes"]):
        k = tuple(np.round(x[:20], 4))
        if k in sig:
            dups.append((sig[k], (s["session"], i)))
        sig[k] = (s["session"], i)
print("identical spike trains across or within sessions:", dups)
print("motion signal source per session:", {s["session"]: s["motion"]["filttype"] for s in S})

LFP rate 1250 Hz in all sessions; pE/pI indices valid and disjoint; spike times sorted
SWSPacket: 701 intervals, non-integer bounds 0, zero or negative length 0
SWSEpisode: 365 intervals, non-integer bounds 0, zero or negative length 0
MA: 2130 intervals, non-integer bounds 0, zero or negative length 0
REM: 417 intervals, non-integer bounds 1, zero or negative length 0
WakeInterruption: 2139 intervals, non-integer bounds 0, zero or negative length 0
Wake: 51 intervals, non-integer bounds 0, zero or negative length 0
Sleep: 51 intervals, non-integer bounds 1, zero or negative length 0
WAKE-SLEEP cycles 51; WAKE 453 to 4475 s; SLEEP 1239 to 7536 s; gap WAKE end to SLEEP start 0 to 0 s
WAKE intervals equal the cycles' WAKE: True
SWSPacket: n 701, duration 14 to 735 s; outside definition 1
MA: n 434, duration 3 to 39 s; outside definition 0
WakeInterruption: n 498, duration 3 to 117 s; outside definition 345
MA intervals in the files 2130, inside SLEEP of a cycle 434
spike times within the LFP duration in every session: True
sessions recorded at the same time on two probes: [['Dino_061914_ACC', 'Dino_061914_mPFC'], ['Dino_062014_ACC', 'Dino_062014_mPFC']]
identical spike trains across or within sessions: []
motion signal source per session: {'20140526_277um': 'clean', '20140527_421um': 'clean', '20140528_565um': 'clean', 'BWRat17_121712': 'noisybaseline', 'BWRat17_121912': 'noisybaseline', 'BWRat18_020513': 'clean', 'BWRat19_032413': 'clean', 'BWRat19_032513': 'clean', 'BWRat20_101013': 'clean', 'BWRat20_101513': 'clean', 'BWRat21_121113': 'clean', 'BWRat21_121613': 'noisybaseline', 'BWRat21_121813': 'noisybaseline', 'Bogey_012615': 'noisybaseline', 'Dino_061814_mPFC': 'noisybaseline', 'Dino_061914_ACC': 'noisybaseline', 'Dino_061914_mPFC': 'noisybaseline', 'Dino_062014_ACC': 'noisybaseline', 'Dino_062014_mPFC': 'noisybaseline', 'Dino_072114_mPFC': 'clean', 'Dino_072314_mPFC': 'clean', 'Dino_072414_mPFC': 'clean', 'Rizzo_022615': 'noisybaseline', 'Rizzo_022715': 'clean', 'Splinter_020515': 'clean', 'Splinter_020915': 'clean', 'Templeton_032415': 'clean'}

# Clock check (protocol, step 2): the state labels are shifted in time and the
# contrast between nonREM packets and WAKE is computed for three streams, the
# LFP delta power (z-scored), the pE burst fraction (ISIs <= 15 ms over all
# ISIs, from the spike trains) and the incidence of detected UP states. If a
# stream shared the clock of the state table only up to an offset, the
# contrast would peak away from lag 0.
def auc(a, b):
    return stats.mannwhitneyu(a, b).statistic / (len(a) * len(b))

lags = [-600, -120, -30, -10, -3, 0, 3, 10, 30, 120, 600]
clock = []
for s, tab in zip(S, tabs):
    lab = tab["label"]
    T = tab["T"]
    bfrac = tab["burst"][s["E"]].sum(0) / np.maximum(tab["isis"][s["E"]].sum(0), 1)
    for lag in lags:
        sh = np.full(T, -1)
        if lag >= 0:
            sh[lag:] = lab[:T - lag]
        else:
            sh[:lag] = lab[-lag:]
        nr, wk = sh == 2, sh == 1
        if nr.sum() < 100 or wk.sum() < 100:
            continue
        clock.append(dict(session=s["session"], lag=lag,
                          delta=auc(tab["zband"]["delta"][nr], tab["zband"]["delta"][wk]),
                          burst=auc(bfrac[nr], bfrac[wk]),
                          up=tab["upn"][nr].mean() - tab["upn"][wk].mean()))
clock = pd.DataFrame(clock)
best = clock.loc[clock.groupby("session")["delta"].idxmax(), ["session", "lag"]].rename(columns={"lag": "delta_best_lag"})
best["burst_best_lag"] = clock.loc[clock.groupby("session")["burst"].idxmax(), "lag"].values
best["up_best_lag"] = clock.loc[clock.groupby("session")["up"].idxmax(), "lag"].values
best["delta_auc_lag0"] = clock[clock.lag == 0].set_index("session").loc[best.session, "delta"].values
best["burst_auc_lag0"] = clock[clock.lag == 0].set_index("session").loc[best.session, "burst"].values
print(best.round(3).to_string(index=False))
print("\nmean over sessions by lag:")
print(clock.groupby("lag")[["delta", "burst", "up"]].mean().round(3).T.to_string())

         session  delta_best_lag  burst_best_lag  up_best_lag  delta_auc_lag0  burst_auc_lag0
  20140526_277um               0               0            0           0.980           0.701
  20140527_421um               0               3            0           0.994           0.574
  20140528_565um               0               0            0           0.980           0.727
  BWRat17_121712               0              -3            0           0.799           0.692
  BWRat17_121912               3               0            0           0.671           0.707
  BWRat18_020513               0               0            0           0.990           0.638
  BWRat19_032413               0               0            0           0.892           0.724
  BWRat19_032513               0               0            0           0.898           0.805
  BWRat20_101013               0               0            0           1.000           0.724
  BWRat20_101513               3               0            0           0.834           0.722
  BWRat21_121113               0               0            0           0.826           0.724
  BWRat21_121613               0               3            0           0.999           0.622
  BWRat21_121813               0               0            0           0.982           0.660
    Bogey_012615               0               0            0           0.995           0.808
Dino_061814_mPFC               0               3            0           0.983           0.846
 Dino_061914_ACC               0              10            0           0.999           0.734
Dino_061914_mPFC               0               0            0           0.996           0.790
 Dino_062014_ACC               0               3            0           0.857           0.623
Dino_062014_mPFC               0               3            0           0.915           0.775
Dino_072114_mPFC               0              -3           -3           0.979           0.743
Dino_072314_mPFC               0               3            0           0.944           0.696
Dino_072414_mPFC               0              10           -3           0.974           0.648
    Rizzo_022615               0               3            0           0.928           0.583
    Rizzo_022715               0               3            0           0.977           0.722
 Splinter_020515               0               3            0           0.995           0.717
 Splinter_020915               0               0            0           0.995           0.573
Templeton_032415               0               0            0           0.906           0.576

mean over sessions by lag:
lag     -600   -120   -30    -10    -3      0      3      10     30     120    600
delta  0.611  0.773  0.872  0.915  0.933  0.936  0.934  0.913  0.854  0.731  0.543
burst  0.550  0.630  0.662  0.685  0.696  0.698  0.698  0.693  0.673  0.633  0.529
up     0.200  0.283  0.346  0.393  0.411  0.416  0.411  0.395  0.355  0.271  0.037

# Ranges and units of the variables used.
cat = lambda f: np.concatenate([f(s) for s in S])
rng = lambda x: f"{np.nanmin(x):.4g} to {np.nanmax(x):.4g}"
print("spike times (s):", rng(cat(lambda s: np.concatenate([x for x in s['spikes'] if len(x)]))))
print("UP durations (s):", rng(cat(lambda s: np.diff(s['updown']['UP'], axis=1).ravel())),
      "; DOWN durations (s):", rng(cat(lambda s: np.diff(s['updown']['DN'], axis=1).ravel())))
print("spindle duration (s):", rng(cat(lambda s: s['spindles']['duration'])),
      "; frequency (Hz):", rng(cat(lambda s: s['spindles']['frequency'])),
      "; amplitude (V, as the description says):", rng(cat(lambda s: s['spindles']['amplitude'])))
print("motion (the description says z-scored):", rng(cat(lambda s: s['motion']['motion'])),
      "; samples per second of recording:", {s['session']: round(len(s['motion']['motion']) / (s['eeg_size'] / 2 / s['meta']['nChannels'] / 1250), 3) for s in S[:3]}, "...")
print("motion samples / recording seconds, all sessions: ", rng([len(s['motion']['motion']) / (s['eeg_size'] / 2 / s['meta']['nChannels'] / 1250) for s in S]))
print("motion mean and SD per session:", {s['session']: (round(float(np.mean(s['motion']['motion'])), 3), round(float(np.std(s['motion']['motion'])), 3)) for s in S})
for b in ("delta", "sigma", "gamma", "theta"):
    print(f"{b} band power (uV^2 per FFT bin, Hann window, arbitrary scale):", rng(cat(lambda s: s['lfp']['band'][b])))
print("LFP samples / (eeg bytes / 2 / nChannels):", {s['session']: s['lfp']['n_samples'] == s['eeg_size'] // 2 // s['meta']['nChannels'] for s in S[:2]}, "...",
      all(s['lfp']['n_samples'] == s['eeg_size'] // 2 // s['meta']['nChannels'] for s in S))

spike times (s): 0.0011 to 4.056e+04
UP durations (s): 0.201 to 3.998 ; DOWN durations (s): 0.076 to 1.247
spindle duration (s): 0.3 to 2.971 ; frequency (Hz): 2.79 to 26.43 ; amplitude (V, as the description says): 3.688e-06 to 0.000572
motion (the description says z-scored): 0 to 1.412e+05 ; samples per second of recording: {'20140526_277um': 1.0, '20140527_421um': 1.0, '20140528_565um': 1.0} ...
motion samples / recording seconds, all sessions:  0.9998 to 1
motion mean and SD per session: {'20140526_277um': (1085.742, 1032.812), '20140527_421um': (906.913, 1186.519), '20140528_565um': (4076.33, 2490.719), 'BWRat17_121712': (0.777, 0.732), 'BWRat17_121912': (1.403, 1.458), 'BWRat18_020513': (2.405, 2.919), 'BWRat19_032413': (0.615, 0.433), 'BWRat19_032513': (nan, nan), 'BWRat20_101013': (1.034, 1.332), 'BWRat20_101513': (68.596, 362.136), 'BWRat21_121113': (403.124, 1060.778), 'BWRat21_121613': (0.795, 0.382), 'BWRat21_121813': (1.433, 1.751), 'Bogey_012615': (0.788, 0.272), 'Dino_061814_mPFC': (0.734, 0.353), 'Dino_061914_ACC': (0.798, 0.346), 'Dino_061914_mPFC': (0.798, 0.346), 'Dino_062014_ACC': (0.766, 0.354), 'Dino_062014_mPFC': (0.766, 0.354), 'Dino_072114_mPFC': (1036.312, 4259.142), 'Dino_072314_mPFC': (486.36, 2996.26), 'Dino_072414_mPFC': (313.608, 2262.511), 'Rizzo_022615': (0.726, 0.498), 'Rizzo_022715': (0.553, 0.707), 'Splinter_020515': (239.544, 706.952), 'Splinter_020915': (0.746, 0.355), 'Templeton_032415': (0.903, 0.008)}
delta band power (uV^2 per FFT bin, Hann window, arbitrary scale): 8.56e+06 to 3.869e+11
sigma band power (uV^2 per FFT bin, Hann window, arbitrary scale): 2.129e+06 to 4.224e+09
gamma band power (uV^2 per FFT bin, Hann window, arbitrary scale): 1.828e+05 to 2.343e+08
theta band power (uV^2 per FFT bin, Hann window, arbitrary scale): 5.047e+06 to 1.704e+10
LFP samples / (eeg bytes / 2 / nChannels): {'20140526_277um': True, '20140527_421um': True} ... True

counts = dict(pE=int(ov.pE.sum()), pI=int(ov.pI.sum()), cycles=int(ov.cycles.sum()),
              sessions=len(ov), rats=ov.rat.nunique(), mean_pE=round(ov.pE.mean(), 2),
              pE_range=(int(ov.pE.min()), int(ov.pE.max())), mean_pI=round(ov.pI.mean(), 2),
              pI_range=(int(ov.pI.min()), int(ov.pI.max())))
expected = dict(pE=995, pI=126, cycles=54, sessions=27, rats=11, mean_pE=36.9)
for k, v in expected.items():
    print(f"{k:9s} measured {counts[k]!s:8s} paper {v:6}  {'match' if abs(counts[k] - v) <= (0.1 if k == 'mean_pE' else 0) else 'DIFFERS'}")
print("pE range", counts["pE_range"], "(paper 10-104); mean pI", counts["mean_pI"], "range", counts["pI_range"], "(paper 4.67, 0-13)")
print("WAKE-SLEEP cycles per session:", ov.set_index("session").cycles.to_dict())

pE        measured 995      paper    995  match
pI        measured 126      paper    126  match
cycles    measured 51       paper     54  DIFFERS
sessions  measured 27       paper     27  match
rats      measured 11       paper     11  match
mean_pE   measured 36.85    paper   36.9  match
pE range (10, 104) (paper 10-104); mean pI 4.67 range (0, 13) (paper 4.67, 0-13)
WAKE-SLEEP cycles per session: {'20140526_277um': 2, '20140527_421um': 1, '20140528_565um': 1, 'BWRat17_121712': 1, 'BWRat17_121912': 1, 'BWRat18_020513': 2, 'BWRat19_032413': 1, 'BWRat19_032513': 1, 'BWRat20_101013': 1, 'BWRat20_101513': 2, 'BWRat21_121113': 2, 'BWRat21_121613': 1, 'BWRat21_121813': 3, 'Bogey_012615': 3, 'Dino_061814_mPFC': 2, 'Dino_061914_ACC': 1, 'Dino_061914_mPFC': 1, 'Dino_062014_ACC': 3, 'Dino_062014_mPFC': 3, 'Dino_072114_mPFC': 1, 'Dino_072314_mPFC': 2, 'Dino_072414_mPFC': 1, 'Rizzo_022615': 5, 'Rizzo_022715': 2, 'Splinter_020515': 3, 'Splinter_020915': 2, 'Templeton_032415': 3}

ma_rows = []
for s in S:
    for c in fcx.cycles(s):
        iv = fcx.within(s["states"]["MA"], *c["sleep"])
        L = c["sleep"][1] - c["sleep"][0]
        ma_rows.append(dict(session=s["session"], n=len(iv), dur=list(np.diff(iv, axis=1).ravel()),
                            sleep_min=L / 60, frac=np.diff(iv, axis=1).sum() / L))
ma = pd.DataFrame(ma_rows)
alld = np.concatenate([d for d in ma.dur if len(d)])
per_sess = ma.groupby("session").agg(n=("n", "sum"), sleep_min=("sleep_min", "sum"))
per_sess["rate"] = per_sess.n / per_sess.sleep_min
per_sess["pct"] = 100 * ma.assign(t=ma.frac * ma.sleep_min).groupby("session").t.sum() / per_sess.sleep_min
sess_med = [np.median(np.concatenate([d for d in g.dur if len(d)])) for _, g in ma.groupby("session") if g.n.sum()]
print(f"MAs in SLEEP: {len(alld)} (paper 790); all MA intervals in the files: {sum(len(s['states']['MA']) for s in S)}")
print(f"median duration (a) over all MAs {np.median(alld):.1f} s, (b) median of session medians {np.median(sess_med):.1f} s (paper 15.9)")
print(f"incidence per min of SLEEP, per session: median {per_sess.rate.median():.3f}, mean {per_sess.rate.mean():.3f} +- {per_sess.rate.std():.3f} (paper 0.32 +- 0.12)")
print(f"percent of SLEEP time, per session: median {per_sess.pct.median():.2f}, mean {per_sess.pct.mean():.2f} +- {per_sess.pct.std():.2f} (paper 9.1 +- 3.5)")

def ma_frac(s, tab, secs):
    return np.mean(tab["label"][secs] == 3)

ma_tr = {}
for scope in ("episode", "sleep"):
    for nb in (10, 20):
        labels = (2, 3) if scope == "episode" else (1, 2, 3, 4, 5)
        rc = fcx.trend(S, tabs, scope, ma_frac, "c", labels=labels, normalize=False, nbins=nb)
        # pooled: one point per period and bin
        T, V = [], []
        for s, tab in zip(S, tabs):
            for u, tn, secs in fcx.windows(s, tab, scope, "c", labels, nb):
                if len(secs):
                    T.append(tn); V.append(ma_frac(s, tab, secs))
        r, p, n = fcx.pearson(T, V)
        ma_tr[(scope, nb)] = (rc["r"], rc["p"], r, p, n)
        print(f"MA fraction vs normalized {scope} time, {nb} bins: bin means R = {rc['r']:.3f} p = {rc['p']:.2g} (n {rc['n']}); "
              f"pooled period-bins R = {r:.3f} p = {p:.2g} (n {n})")
print("paper: episodes R = -0.85, p < 1e-10; SLEEP not significant")

MAs in SLEEP: 434 (paper 790); all MA intervals in the files: 2130
median duration (a) over all MAs 11.0 s, (b) median of session medians 11.5 s (paper 15.9)
incidence per min of SLEEP, per session: median 0.184, mean 0.185 +- 0.112 (paper 0.32 +- 0.12)
percent of SLEEP time, per session: median 3.70, mean 4.71 +- 3.18 (paper 9.1 +- 3.5)

MA fraction vs normalized episode time, 10 bins: bin means R = -0.970 p = 3.3e-06 (n 10); pooled period-bins R = -0.164 p = 1.1e-22 (n 3528)

MA fraction vs normalized episode time, 20 bins: bin means R = -0.956 p = 4.5e-11 (n 20); pooled period-bins R = -0.147 p = 5.7e-35 (n 7011)

MA fraction vs normalized sleep time, 10 bins: bin means R = -0.495 p = 0.15 (n 10); pooled period-bins R = -0.122 p = 0.0057 (n 510)
MA fraction vs normalized sleep time, 20 bins: bin means R = -0.381 p = 0.098 (n 20); pooled period-bins R = -0.107 p = 0.00061 (n 1020)
paper: episodes R = -0.85, p < 1e-10; SLEEP not significant

# Variant, not the criterion (MA reading b). The files hold WakeInterruption intervals shorter than
# 40 s, which the data description defines as 40-120 s. The Results define an MA as an arousal under
# 40 s and "restricted our classification of these states to only those which occur between nonREM
# periods" (Results, Identification of brain states). First, the states on either side of them
# (within 2 s): P packet, R REM, M MA, W wake interruption, - none.
nb_cnt = {}
for s in S:
    for c in fcx.cycles(s):
        seq = fcx.sleep_sequence(s, c)
        for j, (k, a, b) in enumerate(seq):
            if k == "W" and b - a < 40:
                prev = seq[j - 1][0] if j > 0 and a - seq[j - 1][2] <= 2 else "-"
                nxt = seq[j + 1][0] if j + 1 < len(seq) and seq[j + 1][1] - b <= 2 else "-"
                nb_cnt[prev + "W" + nxt] = nb_cnt.get(prev + "W" + nxt, 0) + 1
print(f"WakeInterruptions < 40 s in SLEEP: {sum(nb_cnt.values())}; by neighbours", dict(sorted(nb_cnt.items(), key=lambda x: -x[1])),
      "; between two packets:", nb_cnt.get("PWP", 0))
# Counted with the MAs they give:
wi_rows = []
for s in S:
    for c in fcx.cycles(s):
        ma_iv = fcx.within(s["states"]["MA"], *c["sleep"])
        wi = fcx.within(s["states"]["WakeInterruption"], *c["sleep"])
        wi = wi[np.diff(wi, axis=1).ravel() < 40]
        both = np.vstack([ma_iv, wi])
        wi_rows.append(dict(session=s["session"], n=len(both), t=np.diff(both, axis=1).sum(),
                            sleep_min=(c["sleep"][1] - c["sleep"][0]) / 60, d=np.diff(both, axis=1).ravel()))
wi_df = pd.DataFrame(wi_rows)
g = wi_df.groupby("session").agg(n=("n", "sum"), t=("t", "sum"), m=("sleep_min", "sum"))
print(f"MA + WakeInterruption < 40 s in SLEEP: n {wi_df.n.sum()}, median duration {np.median(np.concatenate(list(wi_df.d))):.1f} s, "
      f"incidence per session median {np.median(g.n / g.m):.3f}/min, percent of SLEEP median {np.median(100 * g.t / g.m / 60):.2f}")
# The MA trends under this reading (used again for the MA slopes and triplets below)
SB, TB = map(list, zip(*[fcx.short_wi_as_ma(s, t) for s, t in zip(S, tabs)]))
for scope in ("episode", "sleep"):
    for nbin in (10, 20):
        labels = (2, 3) if scope == "episode" else (1, 2, 3, 4, 5)
        rc = fcx.trend(SB, TB, scope, ma_frac, "c", labels=labels, normalize=False, nbins=nbin)
        print(f"MA reading (b), MA fraction vs normalized {scope} time, {nbin} bins: bin means R = {rc['r']:.3f} p = {rc['p']:.2g}")
# Variant: the two Dino ACC sessions, recorded at the same time as the mPFC sessions of the same days, left out
S1_ = [s for s in S if not s["session"].endswith("_ACC")]
T1_ = [t for s, t in zip(S, tabs) if not s["session"].endswith("_ACC")]
for scope in ("episode", "sleep"):
    labels = (2, 3) if scope == "episode" else (1, 2, 3, 4, 5)
    rc = fcx.trend(S1_, T1_, scope, ma_frac, "c", labels=labels, normalize=False, nbins=20)
    print(f"Dino ACC sessions left out, MA fraction vs normalized {scope} time, 20 bins: R = {rc['r']:.3f} p = {rc['p']:.2g}")

WakeInterruptions < 40 s in SLEEP: 345; by neighbours {'RWR': 123, 'RWP': 80, 'RW-': 50, '-WP': 49, 'PWR': 25, 'PW-': 7, '-WR': 6, '-W-': 5} ; between two packets: 0
MA + WakeInterruption < 40 s in SLEEP: n 779, median duration 13.0 s, incidence per session median 0.318/min, percent of SLEEP median 7.85
MA reading (b), MA fraction vs normalized episode time, 10 bins: bin means R = -0.977 p = 1.2e-06

MA reading (b), MA fraction vs normalized episode time, 20 bins: bin means R = -0.960 p = 1.9e-11
MA reading (b), MA fraction vs normalized sleep time, 10 bins: bin means R = 0.172 p = 0.63
MA reading (b), MA fraction vs normalized sleep time, 20 bins: bin means R = 0.022 p = 0.93

Dino ACC sessions left out, MA fraction vs normalized episode time, 20 bins: R = -0.949 p = 1.7e-10
Dino ACC sessions left out, MA fraction vs normalized sleep time, 20 bins: R = -0.390 p = 0.089

pk = np.concatenate([np.diff(fcx.state_intervals(s, "SWSPacket"), axis=1).ravel() for s in S])
print(f"packets in SLEEP of the cycles: {len(pk)}; median {np.median(pk):.1f} s (paper 111 +- 11.1); "
      f"range {pk.min():.0f} to {pk.max():.0f} s, ratio {pk.max() / pk.min():.1f}; 5th-95th percentile {np.percentile(pk, 5):.0f}-{np.percentile(pk, 95):.0f} s")
pk_all = np.concatenate([np.diff(s["states"]["SWSPacket"], axis=1).ravel() for s in S])
print(f"all packets in the files: {len(pk_all)}, median {np.median(pk_all):.1f} s")
print("median per rat:", {r: float(np.median(np.concatenate([np.diff(fcx.state_intervals(s, 'SWSPacket'), axis=1).ravel() for s in S if s['rat'] == r]))) for r in sorted(ov.rat.unique())})

packets in SLEEP of the cycles: 701; median 111.0 s (paper 111 +- 11.1); range 14 to 735 s, ratio 52.5; 5th-95th percentile 28-434 s
all packets in the files: 701, median 111.0 s
median per rat: {'BWRat17': 67.5, 'BWRat18': 245.0, 'BWRat19': 83.0, 'BWRat20': 85.0, 'BWRat21': 126.0, 'Bogey': 95.5, 'Dino': 105.0, 'JennBuzsaki22': 118.0, 'Rizzo': 123.0, 'Splinter': 179.5, 'Templeton': 109.0}

R = {"a": fcx.cell_rates(S, "a"), "b": fcx.cell_rates(S, "b")}
exp_med = {("E", "WAKE"): 0.76, ("E", "nonREM"): 0.69, ("E", "REM"): 0.88,
           ("I", "WAKE"): 5.59, ("I", "nonREM"): 4.69, ("I", "REM"): 4.25}
for rd, df in R.items():
    print(f"reading ({rd}): {len(df[df.type == 'E'])} pE rows, {len(df[df.type == 'I'])} pI rows")
    nok = 0
    for (ty, st), v in exp_med.items():
        x = df[df.type == ty][st].dropna()
        m = x.median()
        ok = abs(m - v) <= 0.1 * v
        nok += ok
        print(f"  {ty} {st:6s} median {m:.3f} Hz (paper {v}; {'within' if ok else 'outside'} 10%), n {len(x)}, missing {df[df.type == ty][st].isna().sum()}")
    for st in ("WAKE", "nonREM", "REM"):
        x = df[df.type == "E"][st].dropna()
        print(f"  pE {st:6s} skew (linear) {stats.skew(x):.2f}; 5-95th pct {np.percentile(x, 5):.3f}-{np.percentile(x, 95):.2f} Hz; "
              f"1-99th pct {np.percentile(x, 1):.3f}-{np.percentile(x, 99):.2f} Hz; min {x.min():.4f} max {x.max():.2f}")
    print(f"  medians within tolerance: {nok} of 6")
# pI REM median: the cells at the middle under each reading, and a check for duplicated rows
for rd, df in R.items():
    x = df[df.type == "I"].dropna(subset=["REM"]).sort_values("REM")
    n = len(x)
    mid = x.iloc[sorted({(n - 1) // 2, n // 2})]
    cols = ["session", "cell"] + (["cycle"] if "cycle" in mid.columns else [])
    print(f"reading ({rd}) pI REM: n {n}, median {x.REM.median():.4f} Hz; middle row(s) {mid[cols].values.tolist()} at {mid.REM.round(4).tolist()} Hz")
print("duplicated (session, cell, cycle) rows in reading (a):", int(R["a"].duplicated(subset=["session", "cell", "cycle"]).sum()))

reading (a): 1787 pE rows, 201 pI rows
  E WAKE   median 0.755 Hz (paper 0.76; within 10%), n 1787, missing 0
  E nonREM median 0.706 Hz (paper 0.69; within 10%), n 1787, missing 0
  E REM    median 0.936 Hz (paper 0.88; within 10%), n 1750, missing 37
  I WAKE   median 6.001 Hz (paper 5.59; within 10%), n 201, missing 0
  I nonREM median 4.565 Hz (paper 4.69; within 10%), n 201, missing 0
  I REM    median 5.724 Hz (paper 4.25; outside 10%), n 188, missing 13
  pE WAKE   skew (linear) 2.70; 5-95th pct 0.042-4.24 Hz; 1-99th pct 0.005-6.78 Hz; min 0.0000 max 13.24
  pE nonREM skew (linear) 2.18; 5-95th pct 0.053-2.59 Hz; 1-99th pct 0.011-4.05 Hz; min 0.0000 max 8.06
  pE REM    skew (linear) 2.65; 5-95th pct 0.048-3.98 Hz; 1-99th pct 0.005-5.96 Hz; min 0.0000 max 14.90
  medians within tolerance: 5 of 6
reading (b): 995 pE rows, 126 pI rows
  E WAKE   median 0.784 Hz (paper 0.76; within 10%), n 995, missing 0
  E nonREM median 0.691 Hz (paper 0.69; within 10%), n 995, missing 0
  E REM    median 0.922 Hz (paper 0.88; within 10%), n 958, missing 37
  I WAKE   median 5.632 Hz (paper 5.59; within 10%), n 126, missing 0
  I nonREM median 4.637 Hz (paper 4.69; within 10%), n 126, missing 0
  I REM    median 5.724 Hz (paper 4.25; outside 10%), n 113, missing 13
  pE WAKE   skew (linear) 2.83; 5-95th pct 0.038-4.28 Hz; 1-99th pct 0.004-7.14 Hz; min 0.0000 max 13.24
  pE nonREM skew (linear) 2.38; 5-95th pct 0.050-2.44 Hz; 1-99th pct 0.014-4.17 Hz; min 0.0000 max 8.06
  pE REM    skew (linear) 2.82; 5-95th pct 0.052-3.86 Hz; 1-99th pct 0.008-6.03 Hz; min 0.0000 max 13.28
  medians within tolerance: 5 of 6
reading (a) pI REM: n 188, median 5.7241 Hz; middle row(s) [['Dino_061814_mPFC', 46, 0], ['Dino_061814_mPFC', 46, 1]] at [5.7056, 5.7426] Hz
reading (b) pI REM: n 113, median 5.7241 Hz; middle row(s) [['Dino_061814_mPFC', 46]] at [5.7241] Hz
duplicated (session, cell, cycle) rows in reading (a): 0

fig, ax = plt.subplots(1, 2, figsize=(9, 3))
for ty, a in zip("EI", ax):
    x = R["b"][R["b"].type == ty]
    for st in ("WAKE", "nonREM", "REM"):
        v = x[st][x[st] > 0]
        a.hist(np.log10(v), bins=40, histtype="step", label=st)
    a.set_xlabel("log10 rate (Hz)"); a.set_title(f"p{ty} cells, reading (b)"); a.legend()
plt.tight_layout(); plt.show()
print("cells with zero rate (not on the log axis):", {st: int((R['b'][st] == 0).sum()) for st in ('WAKE', 'nonREM', 'REM')})

cells with zero rate (not on the log axis): {'WAKE': 4, 'nonREM': 1, 'REM': 5}

for rd, df in R.items():
    e = df[df.type == "E"]
    ls = fcx.log_slope(e.WAKE, e.nonREM)
    lin = fcx.ols_slope(e.WAKE, e.nonREM)
    print(f"reading ({rd}) log10: slope {ls['slope']:.3f} [{ls['lo']:.3f}, {ls['hi']:.3f}], R {ls['r']:.3f}, p {ls['p']:.2g}, n {ls['n']}, dropped {ls['dropped']}")
    print(f"reading ({rd}) linear (M5 b): slope {lin['slope']:.3f} [{lin['lo']:.3f}, {lin['hi']:.3f}], R {lin['r']:.3f}, p {lin['p']:.2g}, n {lin['n']}")
e = R["b"][R["b"].type == "E"]
hi = e.WAKE > e.WAKE.median()
print(f"cells above the median WAKE rate: nonREM < WAKE in {np.mean(e.nonREM[hi] < e.WAKE[hi]):.2f}; "
      f"below: nonREM > WAKE in {np.mean(e.nonREM[~hi] > e.WAKE[~hi]):.2f}")
print("paper: slope 0.685 +- 0.025 (95% CI 0.66-0.71), R 0.82 +- 0.082")

reading (a) log10: slope 0.684 [0.661, 0.708], R 0.803, p 0, n 1780, dropped 7
reading (a) linear (M5 b): slope 0.454 [0.437, 0.471], R 0.781, p 0, n 1787
reading (b) log10: slope 0.687 [0.658, 0.717], R 0.824, p 1.4e-246, n 991, dropped 4
reading (b) linear (M5 b): slope 0.449 [0.428, 0.470], R 0.799, p 1.3e-221, n 995
cells above the median WAKE rate: nonREM < WAKE in 0.82; below: nonREM > WAKE in 0.55
paper: slope 0.685 +- 0.025 (95% CI 0.66-0.71), R 0.82 +- 0.082

for rd, df in R.items():
    e = df[df.type == "E"]
    fr = {st: 100 * np.mean(e[st].dropna() > 2) for st in ("WAKE", "REM", "nonREM")}
    print(f"reading ({rd}): % > 2 Hz WAKE {fr['WAKE']:.1f} (19.2), REM {fr['REM']:.1f} (19.2), nonREM {fr['nonREM']:.1f} (8.8)")
    for other in ("WAKE", "REM"):
        tab2 = [[(e[other].dropna() > 2).sum(), (e[other].dropna() <= 2).sum()],
                [(e["nonREM"].dropna() > 2).sum(), (e["nonREM"].dropna() <= 2).sum()]]
        chi = stats.chi2_contingency(tab2)
        print(f"   chi-square nonREM vs {other}: chi2 {chi[0]:.1f}, p {chi[1]:.2g}, table {tab2}")

reading (a): % > 2 Hz WAKE 18.8 (19.2), REM 20.6 (19.2), nonREM 9.3 (8.8)
   chi-square nonREM vs WAKE: chi2 66.2, p 4.1e-16, table [[np.int64(336), np.int64(1451)], [np.int64(166), np.int64(1621)]]
   chi-square nonREM vs REM: chi2 88.8, p 4.4e-21, table [[np.int64(361), np.int64(1389)], [np.int64(166), np.int64(1621)]]
reading (b): % > 2 Hz WAKE 18.4 (19.2), REM 19.5 (19.2), nonREM 9.1 (8.8)
   chi-square nonREM vs WAKE: chi2 35.0, p 3.2e-09, table [[np.int64(183), np.int64(812)], [np.int64(91), np.int64(904)]]
   chi-square nonREM vs REM: chi2 42.2, p 8.3e-11, table [[np.int64(187), np.int64(771)], [np.int64(91), np.int64(904)]]

def dunn_pair(groups, i, j):
    """Dunn's test between groups i and j after Kruskal-Wallis (two-sided, no correction)."""
    allv = np.concatenate(groups)
    rk = stats.rankdata(allv)
    N = len(allv)
    idx = np.cumsum([0] + [len(g) for g in groups])
    mr = [rk[idx[k]:idx[k + 1]].mean() for k in range(len(groups))]
    _, tc = np.unique(allv, return_counts=True)
    tie = np.sum(tc ** 3 - tc) / (12 * (N - 1))
    se = np.sqrt((N * (N + 1) / 12 - tie) * (1 / len(groups[i]) + 1 / len(groups[j])))
    z = (mr[i] - mr[j]) / se
    return z, 2 * stats.norm.sf(abs(z))

BI = []
for s, tab in zip(S, tabs):
    ivs = {st: fcx.state_intervals(s, k) for st, k in (("WAKE", "WAKE"), ("nonREM", "SWSPacket"), ("REM", "REM"))}
    bins = {st: fcx.runs10(tab, lab) for st, lab in (("WAKE", 1), ("nonREM", 2), ("REM", 4))}
    for i in s["E"]:
        r = dict(session=s["session"], cell=i)
        for st in ivs:
            b, n = fcx.isi_counts(s["spikes"][i], ivs[st])
            r[st + "_a"] = b / n if n else np.nan
            b, n = fcx.burst_spikes(s["spikes"][i], ivs[st])
            r[st + "_c"] = b / n if n else np.nan
            if bins[st]:
                B = np.array([tab["burst"][i, w].sum() for w in bins[st]])
                N = np.array([tab["isis"][i, w].sum() for w in bins[st]])
                r[st + "_b"] = np.mean(B[N > 0] / N[N > 0]) if (N > 0).any() else np.nan
            else:
                r[st + "_b"] = np.nan
        BI.append(r)
BI = pd.DataFrame(BI)
exp_bi = dict(nonREM=11.8, WAKE=4.0, REM=5.5)
for rd in "abc":
    g = [BI[f"{st}_{rd}"].dropna().values for st in ("WAKE", "nonREM", "REM")]
    kw = stats.kruskal(*g)
    z, pd_ = dunn_pair(g, 2, 0)
    rs = stats.mannwhitneyu(g[2], g[0])
    print(f"reading ({rd}): " + ", ".join(f"{st} median {100 * np.median(v):.2f}% (paper {exp_bi[st]}), n {len(v)}"
                                       for st, v in zip(("WAKE", "nonREM", "REM"), g)))
    print(f"   Kruskal-Wallis H {kw.statistic:.1f}, p {kw.pvalue:.2g}; REM vs WAKE Dunn z {z:.2f} p {pd_:.2g}; rank-sum p {rs.pvalue:.2g}")

reading (a): WAKE median 2.43% (paper 4.0), n 988, nonREM median 6.26% (paper 11.8), n 994, REM median 3.56% (paper 5.5), n 948
   Kruskal-Wallis H 209.4, p 3.4e-46; REM vs WAKE Dunn z 4.84 p 1.3e-06; rank-sum p 2.4e-06
reading (b): WAKE median 1.58% (paper 4.0), n 991, nonREM median 4.93% (paper 11.8), n 994, REM median 2.33% (paper 5.5), n 955
   Kruskal-Wallis H 246.1, p 3.6e-54; REM vs WAKE Dunn z 4.81 p 1.5e-06; rank-sum p 2.3e-06
reading (c): WAKE median 4.51% (paper 4.0), n 991, nonREM median 11.47% (paper 11.8), n 994, REM median 6.39% (paper 5.5), n 955
   Kruskal-Wallis H 220.8, p 1.1e-48; REM vs WAKE Dunn z 4.41 p 1e-05; rank-sum p 1.6e-05

ACG = {st: [] for st in ("WAKE", "nonREM", "REM")}
for s in S:
    ivs = {st: fcx.state_intervals(s, k) for st, k in (("WAKE", "WAKE"), ("nonREM", "SWSPacket"), ("REM", "REM"))}
    for i in s["E"]:
        for st in ivs:
            ACG[st].append(fcx.acg(s["spikes"][i], ivs[st]))
lag = np.arange(500) + 0.5   # ms, bin centres
box = np.ones(5) / 5
fits = {st: {"a": [], "b": [], "c": []} for st in ACG}
avg = {}
for st, hs in ACG.items():
    H = np.array(hs)
    tot = H.sum(1)
    keep = tot >= 100
    avg[st] = np.convolve((H[keep] / tot[keep, None]).mean(0), box, mode="same")
    for h in H[keep]:
        sm = np.convolve(h, box, mode="same")
        pk = 1 + int(np.argmax(sm[1:50]))
        for rd, end in (("a", 100), ("b", 500), ("c", 500)):
            x, y = lag[pk:end], h[pk:end]
            ok = y > 0
            if ok.sum() < 5:
                fits[st][rd].append(np.nan); continue
            xx = np.log10(x[ok]) if rd != "c" else x[ok] / 1000
            fits[st][rd].append(np.polyfit(xx, np.log10(y[ok]) if rd != "c" else np.log(y[ok]), 1)[0])
    print(f"{st}: {keep.sum()} of {len(H)} pE cells with >= 100 counts")
for st in avg:
    a = avg[st]
    pk = 1 + int(np.argmax(a[1:]))
    w = (lag >= 100) & (lag <= 150)
    loc = [int(lag[k]) for k in range(101, 160) if a[k] > a[k - 1] and a[k] >= a[k + 1] and a[k] == a[k - 5:k + 6].max()]
    print(f"{st}: average autocorrelogram peak at {lag[pk]:.1f} ms; local maxima (11 ms neighbourhood) in 101-160 ms: {loc}")
exp_fit = dict(nonREM=-1.35, WAKE=-0.83, REM=-0.80)
for rd in "abc":
    g = [np.array(fits[st][rd])[np.isfinite(fits[st][rd])] for st in ("WAKE", "nonREM", "REM")]
    kw = stats.kruskal(*g)
    print(f"reading ({rd}): " + ", ".join(f"{st} median {np.median(v):.3f} (paper {exp_fit[st]})" for st, v in zip(("WAKE", "nonREM", "REM"), g))
          + f"; Kruskal-Wallis p {kw.pvalue:.2g}; most negative: {('WAKE', 'nonREM', 'REM')[int(np.argmin([np.median(v) for v in g]))]}")
# Size of the 125 ms bump: mean of the average autocorrelogram at 115-135 ms over its mean at 80-100 and 150-170 ms
for st in avg:
    a = avg[st]
    bump = a[(lag >= 115) & (lag <= 135)].mean() / np.r_[a[(lag >= 80) & (lag <= 100)], a[(lag >= 150) & (lag <= 170)]].mean()
    print(f"{st}: 115-135 ms over flanks {bump:.4f}")

WAKE: 837 of 995 pE cells with >= 100 counts

nonREM: 864 of 995 pE cells with >= 100 counts
REM: 711 of 995 pE cells with >= 100 counts
WAKE: average autocorrelogram peak at 11.5 ms; local maxima (11 ms neighbourhood) in 101-160 ms: [116, 125, 143]
nonREM: average autocorrelogram peak at 8.5 ms; local maxima (11 ms neighbourhood) in 101-160 ms: [125]
REM: average autocorrelogram peak at 11.5 ms; local maxima (11 ms neighbourhood) in 101-160 ms: [104, 111, 126, 145, 157]
reading (a): WAKE median -0.108 (paper -0.83), nonREM median -0.538 (paper -1.35), REM median -0.169 (paper -0.8); Kruskal-Wallis p 3.3e-123; most negative: nonREM
reading (b): WAKE median -0.118 (paper -0.83), nonREM median -0.336 (paper -1.35), REM median -0.112 (paper -0.8); Kruskal-Wallis p 3.1e-103; most negative: nonREM
reading (c): WAKE median -0.582 (paper -0.83), nonREM median -1.623 (paper -1.35), REM median -0.574 (paper -0.8); Kruskal-Wallis p 1.7e-95; most negative: nonREM
WAKE: 115-135 ms over flanks 1.0270
nonREM: 115-135 ms over flanks 0.9808
REM: 115-135 ms over flanks 1.0255

fig, ax = plt.subplots(figsize=(6, 3))
for st in avg:
    ax.plot(lag, avg[st] / avg[st].max(), label=st)
ax.set_xscale("log"); ax.set_xlabel("lag (ms)"); ax.set_ylabel("normalized to peak"); ax.legend()
plt.tight_layout(); plt.show()
print("all plotted values lie in [0, 1]; lags 0.5 to 499.5 ms")

all plotted values lie in [0, 1]; lags 0.5 to 499.5 ms

SY = []
popsync = []
for s in S:
    E = s["E"]
    for c in fcx.cycles(s):
        a, b = c["wake"][0], c["sleep"][1]
        edges = np.arange(a, b + 1e-9, 0.025)
        M = np.array([np.histogram(s["spikes"][i], edges)[0] > 0 for i in E])
        W = M[:, :-1] | M[:, 1:]                 # 50 ms windows, 25 ms step
        n = W.sum(0)
        popsync.append(n / len(E))
        others = (n[None, :] - W) >= 1
        for k, i in enumerate(E):
            SY.append(dict(session=s["session"], cell=i, cycle=c["k"], nwin=W.shape[1],
                           a=(W[k] & others[k]).sum(), b=W[k].sum()))
SY = pd.DataFrame(SY)
g = SY.groupby(["session", "cell"])[["nwin", "a", "b"]].sum()
ps = np.concatenate(popsync)
print(f"population synchrony per window (fraction of pE cells firing): median {np.median(ps):.3f}, mean {ps.mean():.3f}, skew {stats.skew(ps):.2f}, n windows {len(ps)}")
for rd in "ab":
    f = g[rd] / g.nwin
    print(f"reading ({rd}): per-cell fraction median {f.median():.3f}, skew {stats.skew(f):.2f}; "
          f"cells <= 4% of windows {100 * np.mean(f <= 0.04):.1f}% (paper > 15%), cells >= 50% {100 * np.mean(f >= 0.5):.2f}% (paper < 4.2%), n {len(f)}")

population synchrony per window (fraction of pE cells firing): median 0.039, mean 0.047, skew 1.28, n windows 9340981
reading (a): per-cell fraction median 0.029, skew 2.31; cells <= 4% of windows 62.9% (paper > 15%), cells >= 50% 0.00% (paper < 4.2%), n 995
reading (b): per-cell fraction median 0.036, skew 2.26; cells <= 4% of windows 54.2% (paper > 15%), cells >= 50% 0.00% (paper < 4.2%), n 995

UPd = np.concatenate([np.diff(s["updown"]["UP"], axis=1).ravel() for s in S])
DNd = np.concatenate([np.diff(s["updown"]["DN"], axis=1).ravel() for s in S])
print(f"UP states: n {len(UPd)} (paper 43230), median {np.median(UPd):.3f} s (paper 0.47), skew {stats.skew(UPd):.2f}, range {UPd.min():.3f}-{UPd.max():.3f} s")
print(f"DOWN states: n {len(DNd)} (paper 86839), median {np.median(DNd):.3f} s (paper 0.19), skew {stats.skew(DNd):.2f}, range {DNd.min():.3f}-{DNd.max():.3f} s")
inN = []
for s, tab in zip(S, tabs):
    for key in ("UP", "DN"):
        st = s["updown"][key][:, 0].astype(int)
        inN.append((key, len(st), int(np.sum(tab["label"][np.minimum(st, tab["T"] - 1)] == 2))))
x = pd.DataFrame(inN, columns=["kind", "n", "in_packets"]).groupby("kind").sum()
print("starting inside a scored nonREM packet of a cycle:", (x.in_packets / x.n).round(3).to_dict())
print("per session median UP / DOWN (s):", {s["session"]: (round(float(np.median(np.diff(s['updown']['UP'], axis=1))), 3) if len(s['updown']['UP']) else None,
                                                     round(float(np.median(np.diff(s['updown']['DN'], axis=1))), 3)) for s in S})

UP states: n 43230 (paper 43230), median 0.473 s (paper 0.47), skew 2.56, range 0.201-3.998 s
DOWN states: n 86839 (paper 86839), median 0.188 s (paper 0.19), skew 2.15, range 0.076-1.247 s
starting inside a scored nonREM packet of a cycle: {'DN': 1.0, 'UP': 1.0}
per session median UP / DOWN (s): {'20140526_277um': (0.881, 0.131), '20140527_421um': (0.47, 0.169), '20140528_565um': (0.598, 0.158), 'BWRat17_121712': (0.446, 0.224), 'BWRat17_121912': (0.4, 0.214), 'BWRat18_020513': (0.352, 0.214), 'BWRat19_032413': (0.574, 0.239), 'BWRat19_032513': (0.703, 0.209), 'BWRat20_101013': (0.429, 0.219), 'BWRat20_101513': (0.474, 0.173), 'BWRat21_121113': (0.378, 0.196), 'BWRat21_121613': (0.239, 0.313), 'BWRat21_121813': (0.247, 0.227), 'Bogey_012615': (0.459, 0.232), 'Dino_061814_mPFC': (0.802, 0.206), 'Dino_061914_ACC': (0.534, 0.201), 'Dino_061914_mPFC': (0.679, 0.171), 'Dino_062014_ACC': (0.358, 0.268), 'Dino_062014_mPFC': (0.549, 0.189), 'Dino_072114_mPFC': (0.319, 0.205), 'Dino_072314_mPFC': (0.376, 0.203), 'Dino_072414_mPFC': (0.78, 0.14), 'Rizzo_022615': (0.485, 0.162), 'Rizzo_022715': (0.581, 0.151), 'Splinter_020515': (0.347, 0.421), 'Splinter_020915': (0.461, 0.136), 'Templeton_032415': (0.264, 0.184)}

NB = 10
acc = {k: [] for k in ("E", "I", "sync", "EI", "gamma", "sigma")}
for s, tab in zip(S, tabs):
    up = s["updown"]["UP"]
    if len(up) == 0:
        continue
    zf = {k: fcx.fine_z(s, tab, k) for k in ("gamma", "sigma")}
    edges = up[:, [0]] + (up[:, [1]] - up[:, [0]]) * np.linspace(0, 1, NB + 1)[None, :]
    dt = (up[:, 1] - up[:, 0]) / NB
    def counts(cells):
        C = np.zeros((len(cells), len(up), NB))
        for k, i in enumerate(cells):
            x = s["spikes"][i]
            C[k] = np.diff(np.searchsorted(x, edges), axis=1)
        return C
    CE = counts(s["E"])
    acc["E"].append((CE.sum(0) / len(s["E"]) / dt[:, None]))
    fe = (CE > 0).mean(0)
    acc["sync"].append(fe)
    if len(s["I"]):
        CI = counts(s["I"])
        acc["I"].append(CI.sum(0) / len(s["I"]) / dt[:, None])
        fi = (CI > 0).mean(0)
        with np.errstate(invalid="ignore", divide="ignore"):
            acc["EI"].append(fe / (fe + fi))
    for k in ("gamma", "sigma"):
        z, dt_ = zf[k]
        acc[k].append(fcx.mean_in(z, dt_, edges[:, :-1].ravel(), edges[:, 1:].ravel()).reshape(len(up), NB))
tn = (np.arange(NB) + 0.5) / NB
res1e = {}
for k, v in acc.items():
    m = np.nanmean(np.vstack(v), axis=0)
    r, p, _ = fcx.pearson(tn, m)
    res1e[k] = (m, r, p)
    print(f"{k:6s} per bin: {np.round(m, 3)}; last < mean of first two: {m[-1] < m[:2].mean()}; R vs time {r:.3f} p {p:.2g}; n UP states {len(np.vstack(v))}")
print(f"gamma peak bin {int(np.argmax(res1e['gamma'][0]))}, sigma peak bin {int(np.argmax(res1e['sigma'][0]))} (0 = UP onset)")

E      per bin: [1.823 1.359 1.396 1.4   1.389 1.381 1.366 1.348 1.31  1.619]; last < mean of first two: False; R vs time -0.293 p 0.41; n UP states 43230
I      per bin: [15.881 15.267 15.041 14.758 14.045 13.527 12.734 12.048 11.155 11.631]; last < mean of first two: True; R vs time -0.983 p 3.6e-07; n UP states 42745
sync   per bin: [0.079 0.066 0.067 0.067 0.066 0.066 0.065 0.064 0.063 0.071]; last < mean of first two: True; R vs time -0.447 p 0.2; n UP states 43230
EI     per bin: [0.373 0.294 0.292 0.295 0.303 0.311 0.325 0.337 0.351 0.396]; last < mean of first two: False; R vs time 0.470 p 0.17; n UP states 42745
gamma  per bin: [0.157 0.178 0.205 0.207 0.198 0.193 0.179 0.16  0.134 0.072]; last < mean of first two: True; R vs time -0.621 p 0.056; n UP states 43230
sigma  per bin: [0.314 0.147 0.08  0.054 0.039 0.034 0.022 0.021 0.025 0.02 ]; last < mean of first two: True; R vs time -0.774 p 0.0087; n UP states 43230
gamma peak bin 3, sigma peak bin 0 (0 = UP onset)

rb = R["b"].set_index(["session", "cell"])
LAT = []
for s in S:
    up = s["updown"]["UP"]
    if len(up) == 0:
        continue
    t0 = up[:, 0]
    half_a = np.full(len(t0), -1)
    for c in fcx.cycles(s):
        a, b = c["sleep"]
        m = (t0 >= a) & (t0 < b)
        half_a[m] = (t0[m] >= (a + b) / 2).astype(int)
    T = s["eeg_size"] / 2 / s["meta"]["nChannels"] / 1250
    half_b = (t0 >= T / 2).astype(int)
    for i in s["E"]:
        x = s["spikes"][i]
        j = np.searchsorted(x, t0)
        ok = j < len(x)
        lat = np.full(len(t0), np.nan)
        lat[ok] = x[j[ok]] - t0[ok]
        lat[lat > 0.2] = np.nan
        r = dict(session=s["session"], cell=i)
        for rd, h in (("a", half_a), ("b", half_b)):
            for k in (0, 1):
                v = lat[(h == k)]
                r[f"{rd}{k}"] = np.nanmean(v) if np.isfinite(v).sum() >= 5 else np.nan
        r["all"] = np.nanmean(lat) if np.isfinite(lat).sum() >= 5 else np.nan
        LAT.append(r)
LAT = pd.DataFrame(LAT).set_index(["session", "cell"])
LAT = LAT.join(rb[["WAKE", "nonREM"]])
print("cells with >= 5 UP states with a spike within 200 ms counted per half; n cells", len(LAT))
for rd in "ab":
    r_, p_, n_ = fcx.pearson(LAT[f"{rd}0"], LAT[f"{rd}1"])
    print(f"halves ({rd}): first vs second half latency R {r_:.3f} p {p_:.2g} n {n_}")
for rate in ("WAKE", "nonREM"):
    ok = LAT[rate] > 0
    r_, p_, n_ = fcx.pearson(LAT["all"][ok], np.log10(LAT[rate][ok]))
    print(f"latency (all UP states) vs log10 {rate} rate: R {r_:.3f} p {p_:.2g} n {n_}")
    for rd in "ab":
        for k in (0, 1):
            r_, p_, n_ = fcx.pearson(LAT[f"{rd}{k}"][ok], np.log10(LAT[rate][ok]))
            print(f"   halves ({rd}) half {k}: R {r_:.3f} p {p_:.2g} n {n_}")
print("per session R of latency vs log10 WAKE rate:",
      {sname: round(fcx.pearson(g["all"], np.log10(g["WAKE"].where(g["WAKE"] > 0)))[0], 2) for sname, g in LAT.groupby(level=0)})

cells with >= 5 UP states with a spike within 200 ms counted per half; n cells 995
halves (a): first vs second half latency R 0.733 p 7.9e-158 n 932
halves (b): first vs second half latency R 0.660 p 5.4e-96 n 758
latency (all UP states) vs log10 WAKE rate: R -0.127 p 7.6e-05 n 971
   halves (a) half 0: R -0.125 p 0.00012 n 952
   halves (a) half 1: R -0.163 p 4.9e-07 n 939
   halves (b) half 0: R -0.126 p 0.00029 n 823
   halves (b) half 1: R -0.144 p 1.5e-05 n 899
latency (all UP states) vs log10 nonREM rate: R -0.238 p 6e-14 n 972
   halves (a) half 0: R -0.237 p 1.1e-13 n 953
   halves (a) half 1: R -0.293 p 4.5e-20 n 940
   halves (b) half 0: R -0.234 p 1e-11 n 824
   halves (b) half 1: R -0.261 p 1.9e-15 n 900
per session R of latency vs log10 WAKE rate: {'20140526_277um': -0.25, '20140527_421um': -0.04, '20140528_565um': -0.49, 'BWRat17_121712': -0.17, 'BWRat17_121912': -0.11, 'BWRat18_020513': -0.34, 'BWRat19_032413': 0.14, 'BWRat19_032513': 0.05, 'BWRat20_101013': -0.37, 'BWRat20_101513': -0.24, 'BWRat21_121113': 0.18, 'BWRat21_121613': -0.45, 'BWRat21_121813': 0.06, 'Bogey_012615': -0.25, 'Dino_061814_mPFC': -0.09, 'Dino_061914_ACC': 0.12, 'Dino_061914_mPFC': -0.06, 'Dino_062014_ACC': -0.58, 'Dino_062014_mPFC': 0.39, 'Dino_072114_mPFC': -0.52, 'Dino_072314_mPFC': -0.22, 'Dino_072414_mPFC': -0.35, 'Rizzo_022615': -0.43, 'Rizzo_022715': -0.43, 'Splinter_020515': -0.21, 'Splinter_020915': 0.06, 'Templeton_032415': -0.11}

# M2: sextiles of the rate over all of the session's scored WAKE
WR = fcx.session_wake_rates(S)
sext = fcx.sextiles(WR)
_old = fcx.sextiles(R["b"])
print(f"pE cells whose sextile differs from that of the mean of per-cycle WAKE rates: {int((sext != _old.reindex(sext.index)).sum())} of {len(sext)}")
# M2 reading (b): each cell's sextile in each cycle, from that cycle's WAKE rate (one row per
# cell and cycle, ranked over all of them); sx_b gives it for the cycle holding time t.
ra = R["a"][R["a"].type == "E"].copy()
ra["sext_c"] = np.ceil(ra["WAKE"].rank(method="first") / len(ra) * 6).astype(int)
sext_cyc = ra.set_index(["session", "cycle", "cell"])["sext_c"].to_dict()
cyc_at = {}
for s, tab in zip(S, tabs):
    arr = -np.ones(tab["T"], int)
    for c in fcx.cycles(s):
        arr[int(c["wake"][0]):int(c["wake"][1])] = c["k"]
    for c in fcx.cycles(s):
        arr[int(c["sleep"][0]):int(c["sleep"][1])] = c["k"]
    cyc_at[s["session"]] = arr

def sx_b(s, i, t):
    arr = cyc_at[s["session"]]
    return sext_cyc.get((s["session"], int(arr[min(int(t), len(arr) - 1)]), i), 0)

eb = R["b"][R["b"].type == "E"].copy()
eb["sext_nrem"] = np.ceil(eb["nonREM"].rank(method="first") / len(eb) * 6).astype(int)
sext_n = eb.set_index(["session", "cell"])["sext_nrem"]
pairs = []
for s in S:
    iv = fcx.state_intervals(s, "SWSPacket")
    edges = np.concatenate([np.arange(a, b - 0.1 + 1e-9, 0.1) for a, b in iv])
    # counts in 100 ms bins inside packets (bins that start inside a packet and end by its end)
    C = np.array([np.searchsorted(s["spikes"][i], edges + 0.1) - np.searchsorted(s["spikes"][i], edges) for i in s["E"]], float)
    with np.errstate(invalid="ignore"):
        cc = np.corrcoef(C)
    for a in range(len(s["E"])):
        for b in range(a + 1, len(s["E"])):
            ia, ib = s["E"][a], s["E"][b]
            pairs.append((s["session"], cc[a, b], sext[(s["session"], ia)], sext[(s["session"], ib)],
                          sext_n[(s["session"], ia)], sext_n[(s["session"], ib)]))
pairs = pd.DataFrame(pairs, columns=["session", "r", "g1", "g2", "n1", "n2"]).dropna()
for lab, (c1, c2) in (("WAKE sextiles", ("g1", "g2")), ("nonREM sextiles (b)", ("n1", "n2"))):
    lo = pairs[[c1, c2]].min(axis=1)
    hi = pairs[[c1, c2]].max(axis=1)
    M = np.full((6, 6), np.nan)
    for x in range(1, 7):
        for y in range(x, 7):
            v = pairs.r[(lo == x) & (hi == y)]
            M[x - 1, y - 1] = M[y - 1, x - 1] = v.mean()
    print(f"{lab}: mean pair correlation (rows/cols sextile 1 lowest .. 6 highest)\n{np.round(M, 4)}")
    for x, y in ((6, 6), (1, 6), (1, 1)):
        v = pairs.r[(lo == x) & (hi == y)]
        w = stats.wilcoxon(v)
        print(f"   {x}-{y}: mean {v.mean():.4f}, median {v.median():.4f}, n pairs {len(v)}, signed-rank p {w.pvalue:.2g}")
print("pairs:", len(pairs), "from", pairs.session.nunique(), "sessions")

pE cells whose sextile differs from that of the mean of per-cycle WAKE rates: 65 of 995

WAKE sextiles: mean pair correlation (rows/cols sextile 1 lowest .. 6 highest)
[[0.0062 0.0081 0.0068 0.0087 0.0086 0.0064]
 [0.0081 0.0113 0.0133 0.014  0.0148 0.0148]
 [0.0068 0.0133 0.0134 0.016  0.0171 0.0191]
 [0.0087 0.014  0.016  0.0204 0.0206 0.0202]
 [0.0086 0.0148 0.0171 0.0206 0.0233 0.0269]
 [0.0064 0.0148 0.0191 0.0202 0.0269 0.0328]]
   6-6: mean 0.0328, median 0.0242, n pairs 813, signed-rank p 2.9e-70
   1-6: mean 0.0064, median 0.0027, n pairs 1034, signed-rank p 2e-12
   1-1: mean 0.0062, median 0.0019, n pairs 810, signed-rank p 2.1e-20
nonREM sextiles (b): mean pair correlation (rows/cols sextile 1 lowest .. 6 highest)
[[0.0046 0.008  0.0076 0.0075 0.0089 0.0095]
 [0.008  0.0141 0.014  0.0145 0.0128 0.0153]
 [0.0076 0.014  0.015  0.0161 0.0165 0.0203]
 [0.0075 0.0145 0.0161 0.0186 0.0202 0.0229]
 [0.0089 0.0128 0.0165 0.0202 0.0233 0.0266]
 [0.0095 0.0153 0.0203 0.0229 0.0266 0.0313]]
   6-6: mean 0.0313, median 0.0236, n pairs 811, signed-rank p 1.6e-70
   1-6: mean 0.0095, median 0.0067, n pairs 1066, signed-rank p 8.5e-66
   1-1: mean 0.0046, median 0.0017, n pairs 1024, signed-rank p 1.9e-27
pairs: 24006 from 27 sessions

E_cells = lambda s: s["E"]
I_cells = lambda s: s["I"]

def show_trend(label, scope, measure, readings="abc", labels=(2,), normalize=True, SS=None, TT=None, clocks=None):
    """Trend under each reading; for SLEEP and episodes also under the second clock (ii),
    cumulative time of the used state, stored as reading + "ii"."""
    out = {}
    if clocks is None:
        clocks = ("period", "state") if scope in ("sleep", "episode") else ("period",)
    for ck in clocks:
        for rd in readings:
            r = fcx.trend(S if SS is None else SS, tabs if TT is None else TT, scope, measure, rd,
                          labels=labels, normalize=normalize, clock=ck)
            key = rd if ck == "period" else rd + "ii"
            out[key] = r
            print(f"{label} over {scope}, reading ({rd}){' clock (ii)' if ck == 'state' else ''}: R {r['r']:.3f}, p {r['p']:.2g}, n {r['n']}")
    return out

# Dino_061914 and Dino_062014 were recorded from ACC and mPFC at the same time (two sessions each in
# the files, sharing the state scoring); state- and LFP-level trends are also shown with the ACC copies left out.
S1 = [s for s in S if not s["session"].endswith("_ACC")]
T1 = [t for s, t in zip(S, tabs) if not s["session"].endswith("_ACC")]
def one_copy(label, scope, measure, readings="bc", labels=(2,), normalize=True):
    return show_trend(label + " [Dino ACC sessions left out]", scope, measure, readings, labels, normalize, SS=S1, TT=T1, clocks=("period",))

tr_pop = show_trend("pE mean rate", "sleep", fcx.m_rate(E_cells))
# by rat, reading (b)
for rat in sorted(ov.rat.unique()):
    idx = [k for k, s in enumerate(S) if s["rat"] == rat]
    r = fcx.trend([S[k] for k in idx], [tabs[k] for k in idx], "sleep", fcx.m_rate(E_cells), "b")
    print(f"   rat {rat}: R {r['r']:.3f} p {r['p']:.2g} n {r['n']}")

pE mean rate over sleep, reading (a): R -0.273, p 1.7e-13, n 701

pE mean rate over sleep, reading (b): R -0.154, p 2.5e-57, n 10643

pE mean rate over sleep, reading (c): R -0.732, p 1.5e-09, n 50
pE mean rate over sleep, reading (a) clock (ii): R -0.283, p 2.1e-14, n 701

pE mean rate over sleep, reading (b) clock (ii): R -0.163, p 3.6e-64, n 10643
pE mean rate over sleep, reading (c) clock (ii): R -0.796, p 5.1e-12, n 50
   rat BWRat17: R -0.686 p 6.1e-47 n 328
   rat BWRat18: R -0.410 p 4.7e-18 n 409
   rat BWRat19: R -0.636 p 3e-32 n 272
   rat BWRat20: R -0.172 p 9.3e-06 n 659
   rat BWRat21: R -0.195 p 6.9e-12 n 1221
   rat Bogey: R -0.166 p 3.1e-05 n 626
   rat Dino: R -0.169 p 1.6e-14 n 2045
   rat JennBuzsaki22: R -0.296 p 4.5e-16 n 721

   rat Rizzo: R -0.017 p 0.44 n 2010
   rat Splinter: R -0.049 p 0.063 n 1464
   rat Templeton: R -0.014 p 0.68 n 888

FL = []
for s in S:
    for c in fcx.cycles(s):
        pk = fcx.within(s["states"]["SWSPacket"], *c["sleep"])
        if len(pk) < 2:
            continue
        for i in range(len(s["spikes"])):
            x = s["spikes"][i]
            FL.append(dict(session=s["session"], rat=s["rat"], cycle=c["k"], cell=i,
                           type="E" if i in s["E"] else ("I" if i in s["I"] else "?"),
                           first=fcx.rate_in(x, pk[:1]), last=fcx.rate_in(x, pk[-1:])))
FL = pd.DataFrame(FL)
FLb = FL.groupby(["session", "rat", "cell", "type"], as_index=False)[["first", "last"]].mean()
print("SLEEP periods with >= 2 packets:", FL.groupby(["session", "cycle"]).ngroups, "of", int(ov.cycles.sum()))
for rd, df in (("a", FL), ("b", FLb)):
    e = df[df.type == "E"]
    w = stats.wilcoxon(e["first"], e["last"])
    m = stats.mannwhitneyu(e["first"], e["last"])
    print(f"reading ({rd}): first {e['first'].mean():.3f} +- {e['first'].std():.3f} Hz (paper 1.04 +- 1.20), "
          f"last {e['last'].mean():.3f} +- {e['last'].std():.3f} Hz (paper 0.88 +- 0.91), n {len(e)}; "
          f"signed-rank p {w.pvalue:.2g}, rank-sum p {m.pvalue:.2g}; median change {np.median(e['last'] - e['first']):.3f} Hz")
print("by rat, reading (b), mean first / last:", {r: (round(g['first'].mean(), 2), round(g['last'].mean(), 2)) for r, g in FLb[FLb.type == 'E'].groupby('rat')})

SLEEP periods with >= 2 packets: 51 of 51
reading (a): first 1.027 +- 1.086 Hz (paper 1.04 +- 1.20), last 0.905 +- 0.936 Hz (paper 0.88 +- 0.91), n 1787; signed-rank p 1e-18, rank-sum p 0.0025; median change -0.045 Hz
reading (b): first 1.048 +- 1.121 Hz (paper 1.04 +- 1.20), last 0.883 +- 0.916 Hz (paper 0.88 +- 0.91), n 995; signed-rank p 1.9e-20, rank-sum p 0.0012; median change -0.057 Hz
by rat, reading (b), mean first / last: {'BWRat17': (np.float64(0.8), np.float64(0.47)), 'BWRat18': (np.float64(0.64), np.float64(0.51)), 'BWRat19': (np.float64(1.53), np.float64(0.85)), 'BWRat20': (np.float64(0.83), np.float64(0.8)), 'BWRat21': (np.float64(0.67), np.float64(0.56)), 'Bogey': (np.float64(0.87), np.float64(0.81)), 'Dino': (np.float64(1.05), np.float64(0.94)), 'JennBuzsaki22': (np.float64(1.37), np.float64(1.23)), 'Rizzo': (np.float64(1.32), np.float64(1.26)), 'Splinter': (np.float64(0.66), np.float64(0.63)), 'Templeton': (np.float64(0.59), np.float64(0.52))}

WB = []
for s in S:
    cs = fcx.cycles(s)
    wakes = s["states"]["Wake"]
    for c in cs:
        nxt = wakes[wakes[:, 0] >= c["sleep"][1]]
        if len(nxt) == 0:
            continue
        pre = np.array([[c["wake"][1] - 300, c["wake"][1]]])
        post = np.array([[nxt[0, 0], nxt[0, 0] + 300]])
        for i in range(len(s["spikes"])):
            WB.append(dict(session=s["session"], rat=s["rat"], cycle=c["k"], cell=i, gap=nxt[0, 0] - c["sleep"][1],
                           type="E" if i in s["E"] else ("I" if i in s["I"] else "?"),
                           pre=fcx.rate_in(s["spikes"][i], pre), post=fcx.rate_in(s["spikes"][i], post)))
WB = pd.DataFrame(WB)
WBb = WB.groupby(["session", "rat", "cell", "type"], as_index=False)[["pre", "post"]].mean()
print("cycles with a following WAKE:", WB.groupby(["session", "cycle"]).ngroups, "; gap SLEEP end to next WAKE (s):",
      sorted(WB.groupby(["session", "cycle"]).gap.first().astype(int).tolist()))
for rd, df in (("a", WB), ("b", WBb)):
    e = df[df.type == "E"]
    w = stats.wilcoxon(e.pre, e.post)
    m = stats.mannwhitneyu(e.pre, e.post)
    print(f"reading ({rd}): pre {e.pre.mean():.3f} Hz, post {e.post.mean():.3f} Hz (medians {e.pre.median():.3f}, {e.post.median():.3f}), n {len(e)}, "
          f"cells lower after {np.mean(e.post < e.pre):.2f}; signed-rank p {w.pvalue:.2g}; rank-sum p {m.pvalue:.2g}")

cycles with a following WAKE: 23 ; gap SLEEP end to next WAKE (s): [0, 1, 1, 1, 1, 1, 166, 405, 593, 959, 1291, 1357, 1557, 1581, 1813, 2797, 2814, 4539, 5598, 6280, 7789, 9088, 10344]
reading (a): pre 1.145 Hz, post 1.082 Hz (medians 0.713, 0.703), n 703, cells lower after 0.54; signed-rank p 0.018; rank-sum p 0.95
reading (b): pre 1.069 Hz, post 1.005 Hz (medians 0.660, 0.679), n 494, cells lower after 0.53; signed-rank p 0.048; rank-sum p 0.98

# Variant: only the cycles whose next WAKE starts within 1 s of the end of SLEEP.
imm = WB[WB.gap <= 1]
for rd, df in (("a", imm), ("b", imm.groupby(["session", "rat", "cell", "type"], as_index=False)[["pre", "post"]].mean())):
    e = df[df.type == "E"]
    print(f"immediate next WAKE, reading ({rd}): cycles {imm.groupby(['session', 'cycle']).ngroups}, pre {e.pre.mean():.3f}, post {e.post.mean():.3f} Hz, n {len(e)}, "
          f"signed-rank p {stats.wilcoxon(e.pre, e.post).pvalue:.2g}, rank-sum p {stats.mannwhitneyu(e.pre, e.post).pvalue:.2g}")
# Variant (asked of a reader): the paper's "5-min post-sleep WAKE" read as the first 5 min after SLEEP
# ends, whatever their scored state, for every cycle whose recording continues 5 min past SLEEP.
WA = []
for s, tab in zip(S, tabs):
    for c in fcx.cycles(s):
        a = int(c["sleep"][1])
        if a + 300 > tab["T"] - 1:
            continue
        lab = tab["label"][a:a + 300]
        insleep = np.mean(tab["sleep_id"][a:a + 300] >= 0)
        pre = np.array([[c["wake"][1] - 300, c["wake"][1]]])
        post = np.array([[a, a + 300.0]])
        for i in s["E"]:
            WA.append(dict(session=s["session"], cycle=c["k"], cell=i, f_wake=np.mean(lab == 1), f_sleep=insleep,
                           pre=fcx.rate_in(s["spikes"][i], pre), post=fcx.rate_in(s["spikes"][i], post)))
WA = pd.DataFrame(WA)
WAb = WA.groupby(["session", "cell"], as_index=False)[["pre", "post"]].mean()
cyc_wa = WA.groupby(["session", "cycle"]).agg(f_wake=("f_wake", "first"), f_sleep=("f_sleep", "first"))
print(f"cycles with 5 min of recording after SLEEP: {len(cyc_wa)} of 51; of those 5 min, the fraction scored as a cycle's WAKE: "
      f"mean {cyc_wa.f_wake.mean():.2f} (all of it in {np.sum(cyc_wa.f_wake == 1)} cycles, none in {np.sum(cyc_wa.f_wake == 0)}); "
      f"inside another cycle's SLEEP: mean {cyc_wa.f_sleep.mean():.2f}")
for rd, df in (("a", WA), ("b", WAb)):
    print(f"first 5 min after SLEEP, reading ({rd}): pre {df.pre.mean():.3f}, post {df.post.mean():.3f} Hz, n {len(df)}, cells lower after {np.mean(df.post < df.pre):.2f}, "
          f"signed-rank p {stats.wilcoxon(df.pre, df.post).pvalue:.2g}, rank-sum p {stats.mannwhitneyu(df.pre, df.post).pvalue:.2g}")

immediate next WAKE, reading (a): cycles 6, pre 1.165, post 0.998 Hz, n 201, signed-rank p 0.00026, rank-sum p 0.47
immediate next WAKE, reading (b): cycles 6, pre 1.165, post 0.998 Hz, n 201, signed-rank p 0.00026, rank-sum p 0.47

cycles with 5 min of recording after SLEEP: 47 of 51; of those 5 min, the fraction scored as a cycle's WAKE: mean 0.18 (all of it in 3 cycles, none in 38); inside another cycle's SLEEP: mean 0.00
first 5 min after SLEEP, reading (a): pre 1.264, post 1.141 Hz, n 1690, cells lower after 0.56, signed-rank p 3.9e-10, rank-sum p 0.091
first 5 min after SLEEP, reading (b): pre 1.261, post 1.136 Hz, n 995, cells lower after 0.54, signed-rank p 2.8e-06, rank-sum p 0.17

for rd, df in (("a", FL), ("b", FLb)):
    e = df[df.type == "E"]
    ls = fcx.log_slope(e["first"], e["last"])
    lin = fcx.ols_slope(e["first"], e["last"])
    print(f"reading ({rd}) log10: slope {ls['slope']:.3f} [{ls['lo']:.3f}, {ls['hi']:.3f}], R {ls['r']:.3f}, n {ls['n']}, dropped (zero rate) {ls['dropped']}; "
          f"linear slope {lin['slope']:.3f} [{lin['lo']:.3f}, {lin['hi']:.3f}]")
e = FLb[FLb.type == "E"]
hi = e["first"] > e["first"].median()
print(f"reading (b): cells above the median first-packet rate that decrease {np.mean(e['last'][hi] < e['first'][hi]):.2f}; "
      f"below the median that increase {np.mean(e['last'][~hi] > e['first'][~hi]):.2f}")
print("paper: slope 0.8535 +- 0.0265 (95% CI 0.827-0.880)")

reading (a) log10: slope 0.788 [0.760, 0.817], R 0.797, n 1717, dropped (zero rate) 70; linear slope 0.673 [0.648, 0.698]
reading (b) log10: slope 0.872 [0.840, 0.904], R 0.866, n 972, dropped (zero rate) 23; linear slope 0.647 [0.616, 0.678]
reading (b): cells above the median first-packet rate that decrease 0.72; below the median that increase 0.49
paper: slope 0.8535 +- 0.0265 (95% CI 0.827-0.880)

def sext_cells(k, sx=None, t=None):
    """Cells of sextile k: M2 reading (a) when sx is None, else sx(s, cell, time)."""
    if sx is None:
        return lambda s: np.array([i for i in s["E"] if sext[(s["session"], i)] == k], int)
    return lambda s: np.array([i for i in s["E"] if sx(s, i, t) == k], int)

_sx_cells = {}
def m_rate_sx(k, sx):
    """Mean rate of sextile k's cells in a window, the sextile taken at the window's first second."""
    def f(s, tab, secs):
        key = (id(sx), s["session"], int(cyc_at[s["session"]][min(int(secs[0]), len(cyc_at[s["session"]]) - 1)]), k)
        if key not in _sx_cells:
            _sx_cells[key] = sext_cells(k, sx, secs[0])(s)
        c = _sx_cells[key]
        return tab["cnt"][np.ix_(c, secs)].sum() / (len(c) * len(secs)) if len(c) else np.nan
    return f

tr_sext = {}
for k in range(1, 7):
    tr_sext[k] = show_trend(f"sextile {k} (1 lowest)", "sleep", fcx.m_rate(sext_cells(k)))
# Variant, M2 reading (b): each cell's sextile in each cycle from that cycle's WAKE rate
for k in (1, 6):
    _ = show_trend(f"sextile {k} (1 lowest), M2 reading (b)", "sleep", m_rate_sx(k, sx_b))

sextile 1 (1 lowest) over sleep, reading (a): R 0.092, p 0.022, n 611

sextile 1 (1 lowest) over sleep, reading (b): R 0.046, p 8.5e-06, n 9269

sextile 1 (1 lowest) over sleep, reading (c): R 0.311, p 0.028, n 50

sextile 1 (1 lowest) over sleep, reading (a) clock (ii): R 0.090, p 0.026, n 611

sextile 1 (1 lowest) over sleep, reading (b) clock (ii): R 0.044, p 2.1e-05, n 9269

sextile 1 (1 lowest) over sleep, reading (c) clock (ii): R 0.360, p 0.01, n 50

sextile 2 (1 lowest) over sleep, reading (a): R -0.127, p 0.00078, n 701

sextile 2 (1 lowest) over sleep, reading (b): R -0.101, p 1.2e-25, n 10643

sextile 2 (1 lowest) over sleep, reading (c): R -0.652, p 3e-07, n 50

sextile 2 (1 lowest) over sleep, reading (a) clock (ii): R -0.121, p 0.0013, n 701

tr_cv = show_trend("pE rate CV", "sleep", fcx.m_cv, readings="b")
print("variants, not the criterion:")
_ = show_trend("pE rate CV", "sleep", fcx.m_cv, readings="ac")

pE rate CV over sleep, reading (b): R -0.003, p 0.73, n 10643

pE rate CV over sleep, reading (b) clock (ii): R -0.012, p 0.2, n 10643
variants, not the criterion:
pE rate CV over sleep, reading (a): R -0.101, p 0.0073, n 701
pE rate CV over sleep, reading (c): R -0.025, p 0.86, n 50

pE rate CV over sleep, reading (a) clock (ii): R -0.116, p 0.0022, n 701
pE rate CV over sleep, reading (c) clock (ii): R -0.483, p 0.00039, n 50

sp_rows = []
sp_dur_all = []
for s, tab in zip(S, tabs):
    nr = (tab["label"] == 2) & (tab["sleep_id"] >= 0)
    st = s["spindles"]["times"][:, 0].astype(int)
    m = nr[np.minimum(st, tab["T"] - 1)]
    sp_dur_all.append(s["spindles"]["duration"][m])
    sp_rows.append(dict(session=s["session"], n=int(m.sum()), n_all=len(st), nrem_min=nr.sum() / 60,
                        dur=s["spindles"]["duration"][m].mean(), freq=s["spindles"]["frequency"][m].mean(),
                        amp=1e6 * s["spindles"]["amplitude"][m].mean()))
sp = pd.DataFrame(sp_rows)
sp["inc"] = sp.n / sp.nrem_min
print(f"spindles in the files {sp.n_all.sum()}, starting in nonREM of SLEEP {sp.n.sum()}")
for k, v, u in (("inc", 1.94, "/min"), ("dur", 0.52, "s"), ("freq", 13.93, "Hz"), ("amp", 203.3, "uV")):
    print(f"{k}: mean over sessions {sp[k].mean():.3f} +- {sp[k].std():.3f} {u}, median {sp[k].median():.3f} (paper {v})")
sp_dur_all = np.concatenate(sp_dur_all)
print(f"duration over all {len(sp_dur_all)} spindles: mean {sp_dur_all.mean():.3f} +- {sp_dur_all.std():.3f} s, median {np.median(sp_dur_all):.3f} s")
print("Table S1 incidence 0.049 Hz = 2.94/min; measured", round(sp.inc.mean(), 3), "/min =", round(sp.inc.mean() / 60, 4), "Hz")
print("amplitude: data description gives .data amplitude in V; converted to uV")
tr_spi = {}
tr_spi["inc"] = show_trend("spindle incidence", "sleep", fcx.m_events("spn"))
for k in ("spdur", "spfreq", "spamp"):
    tr_spi[k] = show_trend(f"spindle {k[2:]}", "sleep", fcx.m_mean(k, "spn"))
# Variant: the two Dino ACC sessions left out
_ = one_copy("spindle incidence", "sleep", fcx.m_events("spn"), readings="abc")

spindles in the files 5203, starting in nonREM of SLEEP 5201
inc: mean over sessions 2.852 +- 0.490 /min, median 2.955 (paper 1.94)
dur: mean over sessions 0.616 +- 0.170 s, median 0.581 (paper 0.52)
freq: mean over sessions 13.923 +- 0.180 Hz, median 13.943 (paper 13.93)
amp: mean over sessions 208.401 +- 65.964 uV, median 201.211 (paper 203.3)
duration over all 5201 spindles: mean 0.614 +- 0.320 s, median 0.524 s
Table S1 incidence 0.049 Hz = 2.94/min; measured 2.852 /min = 0.0475 Hz
amplitude: data description gives .data amplitude in V; converted to uV
spindle incidence over sleep, reading (a): R -0.029, p 0.44, n 701

spindle incidence over sleep, reading (b): R -0.021, p 0.03, n 10643
spindle incidence over sleep, reading (c): R -0.064, p 0.66, n 50
spindle incidence over sleep, reading (a) clock (ii): R -0.022, p 0.56, n 701

spindle incidence over sleep, reading (b) clock (ii): R -0.013, p 0.19, n 10643
spindle incidence over sleep, reading (c) clock (ii): R -0.082, p 0.57, n 50
spindle dur over sleep, reading (a): R 0.018, p 0.64, n 663

spindle dur over sleep, reading (b): R 0.030, p 0.068, n 3664
spindle dur over sleep, reading (c): R 0.027, p 0.85, n 50
spindle dur over sleep, reading (a) clock (ii): R 0.026, p 0.51, n 663

spindle dur over sleep, reading (b) clock (ii): R 0.037, p 0.025, n 3664
spindle dur over sleep, reading (c) clock (ii): R 0.133, p 0.36, n 50
spindle freq over sleep, reading (a): R -0.059, p 0.13, n 663

spindle freq over sleep, reading (b): R -0.042, p 0.011, n 3664
spindle freq over sleep, reading (c): R -0.251, p 0.078, n 50
spindle freq over sleep, reading (a) clock (ii): R -0.057, p 0.15, n 663

spindle freq over sleep, reading (b) clock (ii): R -0.042, p 0.011, n 3664
spindle freq over sleep, reading (c) clock (ii): R -0.224, p 0.12, n 50
spindle amp over sleep, reading (a): R 0.077, p 0.047, n 663

spindle amp over sleep, reading (b): R 0.031, p 0.063, n 3664
spindle amp over sleep, reading (c): R 0.188, p 0.19, n 50
spindle amp over sleep, reading (a) clock (ii): R 0.077, p 0.049, n 663

spindle amp over sleep, reading (b) clock (ii): R 0.031, p 0.06, n 3664
spindle amp over sleep, reading (c) clock (ii): R 0.262, p 0.066, n 50
spindle incidence [Dino ACC sessions left out] over sleep, reading (a): R -0.029, p 0.46, n 649
spindle incidence [Dino ACC sessions left out] over sleep, reading (b): R -0.015, p 0.15, n 9736

spindle incidence [Dino ACC sessions left out] over sleep, reading (c): R 0.003, p 0.98, n 50

tr_ud = {}
tr_ud["UP incidence"] = show_trend("UP incidence", "sleep", fcx.m_events("upn"))
tr_ud["DOWN incidence"] = show_trend("DOWN incidence", "sleep", fcx.m_events("dnn"))
tr_ud["DOWN duration"] = show_trend("DOWN duration", "sleep", fcx.m_mean("dnd", "dnn"))
tr_ud["UP duration"] = show_trend("UP duration", "sleep", fcx.m_mean("upd", "upn"))
tr_ud["other"] = show_trend("fraction neither UP nor DOWN", "sleep", fcx.m_other)
want = {"UP incidence": (-1, 1e-6), "DOWN incidence": (-1, 1e-4), "DOWN duration": (1, 1e-4), "UP duration": (0, 0.05), "other": (1, 1e-2)}
for rd in ("a", "b", "c", "aii", "bii", "cii"):
    ok = []
    for k, (sgn, bound) in want.items():
        r = tr_ud[k][rd]
        ok.append((r["p"] >= 0.05) if sgn == 0 else (np.sign(r["r"]) == sgn and r["p"] < bound))
    print(f"reading ({rd}): parts holding {sum(ok)} of 5: {dict(zip(want, ok))}")
# Variant: the two Dino ACC sessions left out, reading (b)
for k, m_ in (("UP incidence", fcx.m_events("upn")), ("DOWN incidence", fcx.m_events("dnn")),
              ("DOWN duration", fcx.m_mean("dnd", "dnn")), ("UP duration", fcx.m_mean("upd", "upn")), ("other", fcx.m_other)):
    _ = one_copy(k, "sleep", m_, readings="b")

UP incidence over sleep, reading (a): R -0.060, p 0.12, n 701
UP incidence over sleep, reading (b): R -0.013, p 0.19, n 10643

UP incidence over sleep, reading (c): R -0.155, p 0.28, n 50
UP incidence over sleep, reading (a) clock (ii): R -0.057, p 0.13, n 701

UP incidence over sleep, reading (b) clock (ii): R -0.010, p 0.3, n 10643

UP incidence over sleep, reading (c) clock (ii): R -0.007, p 0.96, n 50
DOWN incidence over sleep, reading (a): R -0.005, p 0.9, n 701

DOWN incidence over sleep, reading (b): R -0.031, p 0.0015, n 10643
DOWN incidence over sleep, reading (c): R -0.009, p 0.95, n 50
DOWN incidence over sleep, reading (a) clock (ii): R 0.008, p 0.83, n 701

DOWN incidence over sleep, reading (b) clock (ii): R -0.018, p 0.065, n 10643
DOWN incidence over sleep, reading (c) clock (ii): R 0.196, p 0.17, n 50
DOWN duration over sleep, reading (a): R 0.153, p 5.5e-05, n 690

DOWN duration over sleep, reading (b): R 0.117, p 2.4e-32, n 10115
DOWN duration over sleep, reading (c): R 0.664, p 1.4e-07, n 50
DOWN duration over sleep, reading (a) clock (ii): R 0.165, p 1.2e-05, n 690

DOWN duration over sleep, reading (b) clock (ii): R 0.130, p 2.8e-39, n 10115
DOWN duration over sleep, reading (c) clock (ii): R 0.733, p 1.4e-09, n 50
UP duration over sleep, reading (a): R 0.092, p 0.021, n 635

UP duration over sleep, reading (b): R 0.020, p 0.067, n 8534
UP duration over sleep, reading (c): R 0.011, p 0.94, n 50
UP duration over sleep, reading (a) clock (ii): R 0.081, p 0.041, n 635

UP duration over sleep, reading (b) clock (ii): R 0.013, p 0.22, n 8534
UP duration over sleep, reading (c) clock (ii): R -0.176, p 0.22, n 50
fraction neither UP nor DOWN over sleep, reading (a): R 0.014, p 0.71, n 701

tr_up = show_trend("pE within-UP rate", "sleep", fcx.m_uprate(E_cells))
UPFL = []
for s in S:
    up = s["updown"]["UP"]
    if len(up) == 0:
        continue
    for c in fcx.cycles(s):
        pk = fcx.within(s["states"]["SWSPacket"], *c["sleep"])
        if len(pk) < 2:
            continue
        r = []
        for p_ in (pk[0], pk[-1]):
            u = fcx.clip(up, *p_)
            r.append(u)
        if any(np.sum(np.diff(u, axis=1)) <= 0 for u in r):
            continue
        for i in range(len(s["spikes"])):
            x = s["spikes"][i]
            UPFL.append(dict(session=s["session"], cell=i, cycle=c["k"],
                             group=("I" if i in s["I"] else (sext.get((s["session"], i), np.nan) if i in s["E"] else np.nan)),
                             group_b=("I" if i in s["I"] else (sx_b(s, i, c["sleep"][0]) if i in s["E"] else np.nan)),
                             first=fcx.rate_in(x, r[0]), last=fcx.rate_in(x, r[1])))
UPFL = pd.DataFrame(UPFL).dropna(subset=["group"])
UPFL["d"] = UPFL["last"] - UPFL["first"]
g = UPFL.groupby(["session", "cell"]).agg(group=("group", "first"), d=("d", "mean"))
for grp, name in ((6, "fastest sextile"), ("I", "pI"), (1, "slowest sextile")):
    v = g.d[g.group == grp]
    w = stats.wilcoxon(v)
    print(f"{name}: mean change last - first packet {v.mean():.3f} Hz, median {v.median():.3f}, n {len(v)}, signed-rank p {w.pvalue:.2g}")
print("all sextiles:", {k: round(g.d[g.group == k].mean(), 3) for k in range(1, 7)})
# M2 reading (b), per-cycle sextiles: one entry per cell and cycle
for grp, name in ((6, "fastest sextile"), (1, "slowest sextile")):
    v = UPFL.d[UPFL.group_b == grp]
    print(f"M2 reading (b) {name}: mean change {v.mean():.3f} Hz, median {v.median():.3f}, n {len(v)} cell-cycles, signed-rank p {stats.wilcoxon(v).pvalue:.2g}")
# Every group under every combination of M2 reading and M1 unit, so the parts can be read on one unit
for m2, col in (("a", "group"), ("b", "group_b")):
    for m1 in ("a", "b"):
        d = UPFL.dropna(subset=[col])
        if m1 == "b":  # one entry per cell (and sextile, under M2 b), averaged over cycles
            d = d.groupby(["session", "cell", col], as_index=False)["d"].mean()
        parts = []
        for grp, name in ((6, "fastest"), ("I", "pI"), (1, "slowest")):
            v = d.d[d[col] == grp]
            parts.append(f"{name} {v.mean():+.3f} Hz (n {len(v)}, p {stats.wilcoxon(v).pvalue:.2g})")
        print(f"M2 ({m2}), M1 ({m1}): " + "; ".join(parts))

pE within-UP rate over sleep, reading (a): R -0.159, p 5.5e-05, n 635

pE within-UP rate over sleep, reading (b): R -0.087, p 9.4e-16, n 8559

pE within-UP rate over sleep, reading (c): R -0.568, p 1.7e-05, n 50
pE within-UP rate over sleep, reading (a) clock (ii): R -0.159, p 5.5e-05, n 635

pE within-UP rate over sleep, reading (b) clock (ii): R -0.088, p 2.5e-16, n 8559

pE within-UP rate over sleep, reading (c) clock (ii): R -0.673, p 8.5e-08, n 50
fastest sextile: mean change last - first packet -0.637 Hz, median -0.455, n 159, signed-rank p 1.2e-07
pI: mean change last - first packet -0.572 Hz, median -0.092, n 118, signed-rank p 0.045
slowest sextile: mean change last - first packet 0.071 Hz, median 0.000, n 151, signed-rank p 0.27

all sextiles: {1: np.float64(0.071), 2: np.float64(-0.005), 3: np.float64(-0.27), 4: np.float64(0.079), 5: np.float64(-0.181), 6: np.float64(-0.637)}
M2 reading (b) fastest sextile: mean change -0.425 Hz, median -0.374, n 276 cell-cycles, signed-rank p 3.7e-07
M2 reading (b) slowest sextile: mean change 0.123 Hz, median 0.000, n 262 cell-cycles, signed-rank p 0.047
M2 (a), M1 (a): fastest -0.414 Hz (n 279, p 4.9e-07); pI -0.672 Hz (n 182, p 0.047); slowest +0.077 Hz (n 244, p 0.48)
M2 (a), M1 (b): fastest -0.637 Hz (n 159, p 1.2e-07); pI -0.572 Hz (n 118, p 0.045); slowest +0.071 Hz (n 151, p 0.27)
M2 (b), M1 (a): fastest -0.425 Hz (n 276, p 3.7e-07); pI -0.672 Hz (n 182, p 0.047); slowest +0.123 Hz (n 262, p 0.047)
M2 (b), M1 (b): fastest -0.573 Hz (n 181, p 3.3e-07); pI -0.572 Hz (n 118, p 0.045); slowest +0.141 Hz (n 189, p 0.023)

FLS = []
for s in S:
    for c in fcx.cycles(s):
        for key in ("REM", "MA"):
            iv = fcx.within(s["states"][key], *c["sleep"])
            if len(iv) < 2:
                continue
            for i in s["E"]:
                x = s["spikes"][i]
                FLS.append(dict(kind=key, session=s["session"], cell=i, cycle=c["k"],
                                first=fcx.rate_in(x, iv[:1]), last=fcx.rate_in(x, iv[-1:])))
FLS = pd.DataFrame(FLS)
wbE = WB[WB.type == "E"].rename(columns={"pre": "first", "post": "last"}).assign(kind="WAKE")
allfl = pd.concat([FLS, wbE[["kind", "session", "cell", "cycle", "first", "last"]]])
exp_s3a = dict(REM=(0.795, 0.035), MA=(0.63, 0.05), WAKE=(0.83, 0.03))
for rd in "ab":
    for kind, (v, tol) in exp_s3a.items():
        df = allfl[allfl.kind == kind]
        if rd == "b":
            df = df.groupby(["session", "cell"], as_index=False)[["first", "last"]].mean()
        ls = fcx.log_slope(df["first"], df["last"])
        print(f"reading ({rd}) {kind}: slope {ls['slope']:.3f} [{ls['lo']:.3f}, {ls['hi']:.3f}] (paper {v} +- {tol}: "
              f"{'within' if abs(ls['slope'] - v) <= tol else 'outside'}), R {ls['r']:.3f}, p {ls['p']:.2g}, n {ls['n']}, dropped {ls['dropped']}")
# Variant (asked of a reader): WAKE after = the first 5 min after SLEEP ends
for rd, df in (("a", WA), ("b", WAb)):
    ls = fcx.log_slope(df["pre"], df["post"])
    print(f"first 5 min after SLEEP, reading ({rd}) WAKE: slope {ls['slope']:.3f} [{ls['lo']:.3f}, {ls['hi']:.3f}], R {ls['r']:.3f}, p {ls['p']:.2g}, n {ls['n']}")
# MA reading (b): WakeInterruptions shorter than 40 s counted as MAs
FLM = []
for s in SB:
    for c in fcx.cycles(s):
        iv = fcx.within(s["states"]["MA"], *c["sleep"])
        if len(iv) < 2:
            continue
        for i in s["E"]:
            x = s["spikes"][i]
            FLM.append(dict(session=s["session"], cell=i, first=fcx.rate_in(x, iv[:1]), last=fcx.rate_in(x, iv[-1:])))
FLM = pd.DataFrame(FLM)
for rd, df in (("a", FLM), ("b", FLM.groupby(["session", "cell"], as_index=False)[["first", "last"]].mean())):
    ls = fcx.log_slope(df["first"], df["last"])
    print(f"MA reading (b), reading ({rd}): slope {ls['slope']:.3f} [{ls['lo']:.3f}, {ls['hi']:.3f}] (paper 0.63 +- 0.05), n {ls['n']}")

reading (a) REM: slope 0.742 [0.708, 0.776] (paper 0.795 +- 0.035: outside), R 0.751, p 1.6e-260, n 1435, dropped 149
reading (a) MA: slope 0.530 [0.483, 0.577] (paper 0.63 +- 0.05: outside), R 0.539, p 3.3e-92, n 1210, dropped 457
reading (a) WAKE: slope 0.713 [0.663, 0.762] (paper 0.83 +- 0.03: outside), R 0.730, p 5.3e-117, n 697, dropped 6
reading (b) REM: slope 0.807 [0.769, 0.845] (paper 0.795 +- 0.035: within), R 0.823, p 5.7e-203, n 819, dropped 56
reading (b) MA: slope 0.656 [0.601, 0.710] (paper 0.63 +- 0.05: within), R 0.641, p 3.2e-94, n 803, dropped 172
reading (b) WAKE: slope 0.687 [0.627, 0.746] (paper 0.83 +- 0.03: outside), R 0.718, p 8.6e-79, n 489, dropped 5
first 5 min after SLEEP, reading (a) WAKE: slope 0.787 [0.760, 0.815], R 0.808, p 0, n 1672
first 5 min after SLEEP, reading (b) WAKE: slope 0.793 [0.759, 0.828], R 0.821, p 3e-241, n 983
MA reading (b), reading (a): slope 0.509 [0.463, 0.556] (paper 0.63 +- 0.05), n 1301
MA reading (b), reading (b): slope 0.618 [0.563, 0.674] (paper 0.63 +- 0.05), n 839

for rd, df in (("a", WB), ("b", WBb)):
    e = df[df.type == "E"].copy()
    e["sext"] = [sext[(a, b)] for a, b in zip(e.session, e.cell)]
    out = []
    for k in range(1, 7):
        v = e[e.sext == k]
        alt = "greater" if k >= 4 else "less"    # pre > post for the top three; pre < post for the others
        p = stats.wilcoxon(v.pre, v.post, alternative=alt).pvalue
        out.append(f"s{k}: pre {v.pre.mean():.2f} post {v.post.mean():.2f} p({'decrease' if k >= 4 else 'increase'}) {p:.2g} n {len(v)}")
    print(f"reading ({rd}): " + "; ".join(out))
imm_b = imm[imm.type == "E"].groupby(["session", "cell"], as_index=False)[["pre", "post"]].mean()
imm_b["sext"] = [sext[(a, b)] for a, b in zip(imm_b.session, imm_b.cell)]
print("variant, immediate next WAKE only:", "; ".join(
    f"s{k}: p {stats.wilcoxon(v.pre, v.post, alternative='greater' if k >= 4 else 'less').pvalue:.2g} n {len(v)}"
    for k, v in imm_b.groupby("sext")))
# Variant (asked of a reader): WAKE after = the first 5 min after SLEEP ends, all cycles
for rd, df in (("a", WA), ("b", WAb)):
    e = df.copy()
    e["sext"] = [sext[(a, b)] for a, b in zip(e.session, e.cell)]
    print(f"first 5 min after SLEEP, reading ({rd}): " + "; ".join(
        f"s{k}: pre {v.pre.mean():.2f} post {v.post.mean():.2f} p({'decrease' if k >= 4 else 'increase'}) "
        f"{stats.wilcoxon(v.pre, v.post, alternative='greater' if k >= 4 else 'less').pvalue:.2g} n {len(v)}" for k, v in e.groupby("sext")))
# Variant, M2 reading (b): sextile of each cell in each cycle from that cycle's WAKE rate
e = WB[WB.type == "E"].copy()
e["sext"] = [sext_cyc.get((a, c, b)) for a, c, b in zip(e.session, e.cycle, e.cell)]
print("M2 reading (b), per-cycle sextiles, reading (a): " + "; ".join(
    f"s{k}: p({'decrease' if k >= 4 else 'increase'}) {stats.wilcoxon(v.pre, v.post, alternative='greater' if k >= 4 else 'less').pvalue:.2g} n {len(v)}"
    for k, v in e.groupby("sext")))
# Parts holding (top three sextiles decrease, lowest increases, one-tailed p < 0.05)
def held_3e(df, col="sext"):
    return int(sum(stats.wilcoxon(df[df[col] == k].pre, df[df[col] == k].post,
                                  alternative="greater" if k >= 4 else "less").pvalue < 0.05 for k in (6, 5, 4, 1)))
_ea = WB[WB.type == "E"].assign(sext=lambda d: [sext[(a, b)] for a, b in zip(d.session, d.cell)])
_eb = WBb[WBb.type == "E"].assign(sext=lambda d: [sext[(a, b)] for a, b in zip(d.session, d.cell)])
held3e = {"M2 a, M1 a": held_3e(_ea), "M2 a, M1 b": held_3e(_eb), "M2 b, M1 a": held_3e(e)}
print("parts holding of 4:", held3e)

reading (a): s1: pre 0.09 post 0.11 p(increase) 0.15 n 113; s2: pre 0.35 post 0.32 p(increase) 0.83 n 117; s3: pre 0.60 post 0.69 p(increase) 0.26 n 136; s4: pre 0.97 post 1.01 p(decrease) 0.48 n 128; s5: pre 1.62 post 1.51 p(decrease) 0.032 n 97; s6: pre 3.49 post 3.05 p(decrease) 0.0018 n 112
reading (b): s1: pre 0.09 post 0.11 p(increase) 0.1 n 86; s2: pre 0.35 post 0.31 p(increase) 0.78 n 82; s3: pre 0.58 post 0.68 p(increase) 0.17 n 94; s4: pre 0.94 post 1.03 p(decrease) 0.61 n 93; s5: pre 1.56 post 1.48 p(decrease) 0.073 n 73; s6: pre 3.57 post 2.95 p(decrease) 0.00032 n 66
variant, immediate next WAKE only: s1: p 0.7 n 30; s2: p 0.34 n 37; s3: p 0.7 n 41; s4: p 0.11 n 33; s5: p 0.028 n 25; s6: p 0.00037 n 35
first 5 min after SLEEP, reading (a): s1: pre 0.09 post 0.14 p(increase) 0.00096 n 260; s2: pre 0.34 post 0.33 p(increase) 0.98 n 278; s3: pre 0.62 post 0.72 p(increase) 0.074 n 294; s4: pre 1.00 post 0.98 p(decrease) 0.084 n 294; s5: pre 1.61 post 1.43 p(decrease) 1.7e-06 n 271; s6: pre 3.76 post 3.12 p(decrease) 9.7e-13 n 293
first 5 min after SLEEP, reading (b): s1: pre 0.09 post 0.13 p(increase) 0.00027 n 165; s2: pre 0.35 post 0.34 p(increase) 0.96 n 166; s3: pre 0.61 post 0.72 p(increase) 0.089 n 166; s4: pre 0.98 post 0.97 p(decrease) 0.11 n 166; s5: pre 1.57 post 1.40 p(decrease) 0.00036 n 166; s6: pre 3.97 post 3.25 p(decrease) 8.5e-10 n 166
M2 reading (b), per-cycle sextiles, reading (a): s1: p(increase) 2e-05 n 127; s2: p(increase) 0.012 n 123; s3: p(increase) 0.48 n 124; s4: p(decrease) 0.23 n 122; s5: p(decrease) 0.00013 n 91; s6: p(decrease) 3.2e-06 n 116
parts holding of 4: {'M2 a, M1 a': 2, 'M2 a, M1 b': 1, 'M2 b, M1 a': 3}

tr_lfp = {}
for scope in ("episode", "packet"):
    for band in ("delta", "delta14", "sigma", "gamma"):
        tr_lfp[(scope, band)] = show_trend(f"{band} power", scope, fcx.m_band(band),
                                           readings="abc" if scope == "episode" else "bc", normalize=False)
PL = []
for s, tab in zip(S, tabs):
    for a, b in fcx.state_intervals(s, "SWSPacket"):
        PL.append((b - a, tab["band"]["delta"][int(a):int(b)].mean() / tab["nrem_delta_mean"]))
PL = np.array(PL)
r, p, n = fcx.pearson(PL[:, 0], PL[:, 1])
print(f"packet mean delta (over session nonREM mean) vs packet duration: R {r:.3f}, p {p:.2g}, n {n} (paper R 0.31, p < 0.05)")
r2, p2, _ = fcx.pearson(np.log10(PL[:, 0]), PL[:, 1])
print(f"   variant with log10 duration: R {r2:.3f}, p {p2:.2g}")
# Variant: seconds whose 10 s window reaches outside the packet (within 5 s of a packet edge) left out
def edge_free(s, tab, w=5):
    t2 = dict(tab)
    lab = tab["label"].copy()
    for a, b in fcx.state_intervals(s, "SWSPacket"):
        lab[int(a):int(a) + w] = 0
        lab[max(int(b) - w + 1, int(a)):int(b)] = 0
    t2["label"] = lab
    return t2
TE = [edge_free(s, t) for s, t in zip(S, tabs)]
for scope in ("episode", "packet"):
    for band in ("delta", "sigma", "gamma"):
        _ = show_trend(f"{band} power [edge seconds left out]", scope, fcx.m_band(band), readings="bc", normalize=False, TT=TE, clocks=("period",))
# Variant: the two Dino ACC sessions left out
for scope in ("episode", "packet"):
    for band in ("delta", "sigma", "gamma"):
        _ = one_copy(f"{band} power", scope, fcx.m_band(band), normalize=False)

delta power over episode, reading (a): R 0.253, p 1e-11, n 701

delta power over episode, reading (b): R 0.276, p 1.1e-184, n 10643

delta power over episode, reading (c): R 0.619, p 1.7e-06, n 50
delta power over episode, reading (a) clock (ii): R 0.305, p 1.4e-16, n 701

delta power over episode, reading (b) clock (ii): R 0.289, p 2.1e-204, n 10643

delta power over episode, reading (c) clock (ii): R 0.604, p 3.4e-06, n 50

delta14 power over episode, reading (a): R 0.233, p 4.4e-10, n 701
delta14 power over episode, reading (b): R 0.242, p 1.8e-141, n 10643

delta14 power over episode, reading (c): R 0.540, p 5.2e-05, n 50
delta14 power over episode, reading (a) clock (ii): R 0.287, p 9.9e-15, n 701
delta14 power over episode, reading (b) clock (ii): R 0.256, p 1.4e-158, n 10643

delta14 power over episode, reading (c) clock (ii): R 0.520, p 0.00011, n 50
sigma power over episode, reading (a): R 0.179, p 1.9e-06, n 701
sigma power over episode, reading (b): R 0.209, p 4.4e-105, n 10643

sigma power over episode, reading (c): R 0.870, p 2.5e-16, n 50
sigma power over episode, reading (a) clock (ii): R 0.187, p 6.3e-07, n 701
sigma power over episode, reading (b) clock (ii): R 0.218, p 1.1e-114, n 10643

sigma power over episode, reading (c) clock (ii): R 0.879, p 4.9e-17, n 50
gamma power over episode, reading (a): R -0.220, p 3.6e-09, n 701
gamma power over episode, reading (b): R -0.129, p 1e-40, n 10643

def m_ratio_dur(s, tab, secs):
    nu, nd = tab["upn"][secs].sum(), tab["dnn"][secs].sum()
    return (tab["upd"][secs].sum() / nu) / (tab["dnd"][secs].sum() / nd) if nu > 0 and nd > 0 else np.nan

def m_ratio_time(s, tab, secs):
    d = tab["dntime"][secs].sum()
    return tab["uptime"][secs].sum() / d if d > 0 else np.nan

tr_udp = {}
for k, f in (("UP incidence", fcx.m_events("upn")), ("DOWN incidence", fcx.m_events("dnn")),
             ("UP duration", fcx.m_mean("upd", "upn")), ("DOWN duration", fcx.m_mean("dnd", "dnn")),
             ("UP/DOWN ratio (a, durations)", m_ratio_dur), ("UP/DOWN ratio (b, time)", m_ratio_time),
             ("fraction neither", fcx.m_other)):
    tr_udp[k] = show_trend(k, "packet", f, readings="bc")
want = {"UP incidence": (1, 1e-3), "DOWN incidence": (1, 1e-10), "UP duration": (-1, 0.01), "DOWN duration": (1, 1e-10)}
for rd in "bc":
    for ratio in ("UP/DOWN ratio (a, durations)", "UP/DOWN ratio (b, time)"):
        ok = [np.sign(tr_udp[k][rd]["r"]) == sg and tr_udp[k][rd]["p"] < bd for k, (sg, bd) in want.items()]
        ok.append(tr_udp[ratio][rd]["r"] < 0 and tr_udp[ratio][rd]["p"] < 0.05)
        print(f"reading ({rd}) with {ratio}: parts holding {sum(ok)} of 5 {ok}")
# Variant: the two Dino ACC sessions left out
for lab_, m_ in (("UP incidence", fcx.m_events("upn")), ("DOWN incidence", fcx.m_events("dnn")),
                 ("UP duration", fcx.m_mean("upd", "upn")), ("DOWN duration", fcx.m_mean("dnd", "dnn"))):
    _ = one_copy(lab_, "packet", m_, readings="b")

UP incidence over packet, reading (b): R 0.042, p 1.3e-05, n 10643

UP incidence over packet, reading (c): R 0.447, p 0.0011, n 50

DOWN incidence over packet, reading (b): R 0.168, p 3.5e-68, n 10643

DOWN incidence over packet, reading (c): R 0.473, p 0.00052, n 50

UP duration over packet, reading (b): R -0.101, p 9.2e-21, n 8534

UP duration over packet, reading (c): R -0.311, p 0.028, n 50
DOWN duration over packet, reading (b): R 0.190, p 3.6e-83, n 10115

DOWN duration over packet, reading (c): R 0.566, p 1.8e-05, n 50
UP/DOWN ratio (a, durations) over packet, reading (b): R -0.156, p 8.1e-48, n 8529

UP/DOWN ratio (a, durations) over packet, reading (c): R -0.371, p 0.0079, n 50
UP/DOWN ratio (b, time) over packet, reading (b): R -0.035, p 0.00051, n 10120

UP/DOWN ratio (b, time) over packet, reading (c): R 0.204, p 0.16, n 50
fraction neither over packet, reading (b): R -0.176, p 1.1e-74, n 10643

fraction neither over packet, reading (c): R -0.525, p 9e-05, n 50
reading (b) with UP/DOWN ratio (a, durations): parts holding 5 of 5 [True, True, True, True, True]
reading (b) with UP/DOWN ratio (b, time): parts holding 5 of 5 [True, True, True, True, True]
reading (c) with UP/DOWN ratio (a, durations): parts holding 1 of 5 [False, False, False, False, True]
reading (c) with UP/DOWN ratio (b, time): parts holding 0 of 5 [False, False, False, False, False]
UP incidence [Dino ACC sessions left out] over packet, reading (b): R 0.044, p 1.5e-05, n 9736
DOWN incidence [Dino ACC sessions left out] over packet, reading (b): R 0.167, p 6.8e-62, n 9736

tr_spw = {}
for scope in ("episode", "packet"):
    rds = "abc" if scope == "episode" else "bc"
    tr_spw[(scope, "inc")] = show_trend("spindle incidence", scope, fcx.m_events("spn"), readings=rds)
    for k in ("spdur", "spfreq", "spamp"):
        tr_spw[(scope, k)] = show_trend(f"spindle {k[2:]}", scope, fcx.m_mean(k, "spn"), readings=rds)
for scope in ("episode", "packet"):
    pos, dur = [], []
    for s in S:
        st = s["spindles"]["times"][:, 0]
        for a, b, _ in fcx.units(s, scope):
            m = (st >= a) & (st < b)
            pos.append((st[m] - a) / (b - a)); dur.append(s["spindles"]["duration"][m])
    pos, dur = np.concatenate(pos), np.concatenate(dur)
    for lab, f in (("tenth", 0.1), ("third", 1 / 3)):
        print(f"{scope}: duration first {lab} {dur[pos < f].mean():.3f} s (n {np.sum(pos < f)}), last {lab} {dur[pos >= 1 - f].mean():.3f} s (n {np.sum(pos >= 1 - f)})")
print("paper: episodes 0.61 -> 0.67 s, packets 0.58 -> 0.67 s (tolerance 10%)")
# Variant: the two Dino ACC sessions, recorded at the same time as the mPFC sessions of the same days, left out
_ = one_copy("spindle incidence", "episode", fcx.m_events("spn"))
_ = one_copy("spindle incidence", "packet", fcx.m_events("spn"))
_ = one_copy("spindle duration", "episode", fcx.m_mean("spdur", "spn"))
_ = one_copy("spindle duration", "packet", fcx.m_mean("spdur", "spn"))

spindle incidence over episode, reading (a): R 0.012, p 0.75, n 701

spindle incidence over episode, reading (b): R 0.083, p 1.1e-17, n 10643

spindle incidence over episode, reading (c): R 0.669, p 1.1e-07, n 50
spindle incidence over episode, reading (a) clock (ii): R 0.026, p 0.5, n 701

spindle incidence over episode, reading (b) clock (ii): R 0.089, p 5.3e-20, n 10643

spindle incidence over episode, reading (c) clock (ii): R 0.709, p 8.1e-09, n 50
spindle dur over episode, reading (a): R 0.107, p 0.0057, n 663

spindle dur over episode, reading (b): R 0.075, p 5.4e-06, n 3664

spindle dur over episode, reading (c): R 0.468, p 0.00061, n 50
spindle dur over episode, reading (a) clock (ii): R 0.102, p 0.0085, n 663

spindle dur over episode, reading (b) clock (ii): R 0.077, p 3.1e-06, n 3664

spindle dur over episode, reading (c) clock (ii): R 0.520, p 0.00011, n 50
spindle freq over episode, reading (a): R -0.015, p 0.7, n 663

spindle freq over episode, reading (b): R -0.027, p 0.11, n 3664

UPS = []
for s, tab in zip(S, tabs):
    up = s["updown"]["UP"]
    if len(up) == 0:
        continue
    z, dt_ = fcx.fine_z(s, tab, "sigma")
    pw = fcx.mean_in(z, dt_, up[:, 0], up[:, 1])
    for a, b in fcx.state_intervals(s, "SWSPacket"):
        m = (up[:, 0] >= a) & (up[:, 1] <= b)
        pos = (up[m, 0] - a) / (b - a)
        for p_, d, w in zip(pos, np.diff(up[m], axis=1).ravel(), pw[m]):
            UPS.append((s["session"], p_, d, w))
UPS = pd.DataFrame(UPS, columns=["session", "pos", "dur", "sigma"])
for rd, m in (("a", UPS.dur >= 1), ("b", UPS.dur > 0.5)):
    e = UPS[m & (UPS.pos < 1 / 3)].sigma
    l = UPS[m & (UPS.pos >= 2 / 3)].sigma
    t = stats.mannwhitneyu(l, e, alternative="two-sided")
    print(f"reading ({rd}): early n {len(e)} mean z {e.mean():.3f}, late n {len(l)} mean z {l.mean():.3f}; rank-sum p {t.pvalue:.2g}; "
          f"t-test p {stats.ttest_ind(l, e).pvalue:.2g}")
e = UPS[UPS.pos < 1 / 3].sigma; l = UPS[UPS.pos >= 2 / 3].sigma
print(f"all UP states: early {e.mean():.3f}, late {l.mean():.3f}, rank-sum p {stats.mannwhitneyu(l, e).pvalue:.2g}")

reading (a): early n 1866 mean z 0.009, late n 1708 mean z 0.153; rank-sum p 3.7e-34; t-test p 3.8e-37
reading (b): early n 6423 mean z 0.031, late n 6816 mean z 0.151; rank-sum p 2.5e-76; t-test p 3.7e-81
all UP states: early 0.036, late 0.127, rank-sum p 8.5e-80

tr_s4 = {}
for scope in ("episode", "packet"):
    for k in (1, 6):
        tr_s4[(scope, k)] = show_trend(f"sextile {k}", scope, fcx.m_rate(sext_cells(k)), readings="b")
    tr_s4[(scope, "cv")] = show_trend("pE rate CV", scope, fcx.m_cv, readings="b")
    tr_s4[(scope, "burst")] = show_trend("burst ratio", scope, fcx.m_burst, readings="b")
print("other sextiles, reading (b):")
for scope in ("episode", "packet"):
    for k in (2, 3, 4, 5):
        r = fcx.trend(S, tabs, scope, fcx.m_rate(sext_cells(k)), "b")
        print(f"   {scope} sextile {k}: R {r['r']:.3f} p {r['p']:.2g}")
print("variant, reading (c) 50 bins:")
for scope in ("episode", "packet"):
    for k in (1, 6):
        r = fcx.trend(S, tabs, scope, fcx.m_rate(sext_cells(k)), "c")
        print(f"   {scope} sextile {k}: R {r['r']:.3f} p {r['p']:.2g}")

sextile 1 over episode, reading (b): R 0.029, p 0.0048, n 9269

sextile 1 over episode, reading (b) clock (ii): R 0.030, p 0.0038, n 9269

sextile 6 over episode, reading (b): R -0.139, p 1.6e-37, n 8452

sextile 6 over episode, reading (b) clock (ii): R -0.145, p 7.2e-41, n 8452

pE rate CV over episode, reading (b): R -0.103, p 1e-26, n 10643

pE rate CV over episode, reading (b) clock (ii): R -0.109, p 1.1e-29, n 10643

burst ratio over episode, reading (b): R 0.061, p 2.6e-10, n 10643

burst ratio over episode, reading (b) clock (ii): R 0.063, p 8.5e-11, n 10643

sextile 1 over packet, reading (b): R 0.044, p 2.3e-05, n 9269

sextile 6 over packet, reading (b): R -0.137, p 1.4e-36, n 8452

PK = []
for s, tab in zip(S, tabs):
    pk = fcx.state_intervals(s, "SWSPacket")
    if len(pk) < 3:
        continue
    rows = []
    for a, b in pk:
        sec = np.arange(int(a), int(b))
        r = dict(session=s["session"], dur=b - a)
        for k in (1, 6):
            for sfx, sx in (("", None), ("b", sx_b)):   # M2 readings (a) and (b)
                c = sext_cells(k, sx, a)(s)
                r[f"s{k}{sfx}"] = tab["cnt"][np.ix_(c, sec)].sum() / (len(c) * len(sec)) if len(c) else np.nan
        r["delta"] = tab["zband"]["delta"][sec].mean()
        r["gamma"] = tab["zband"]["gamma"][sec].mean()
        r["sigma_last20"] = tab["zband"]["sigma"][sec[-20:]].mean()
        r["sigma_last20_in"] = tab["zband"]["sigma"][sec[-15:-4]].mean()   # windows inside the last 20 s
        r["dn_dur"] = fcx.m_mean("dnd", "dnn")(s, tab, sec)
        r["dn_inc"] = fcx.m_events("dnn")(s, tab, sec)
        r["up_dur"] = fcx.m_mean("upd", "upn")(s, tab, sec)
        rows.append(r)
    df = pd.DataFrame(rows)
    for k in ("s1", "s6", "s1b", "s6b"):
        df[k] = df[k] / df[k].mean()
    for k in ("delta", "gamma", "sigma_last20", "sigma_last20_in", "dn_dur", "dn_inc", "up_dur"):
        df[k] = (df[k] - df[k].mean()) / df[k].std()
    PK.append(df)
PK = pd.concat(PK)
want4 = [("s6", "delta", -1, 1e-10), ("s6", "dn_dur", -1, 1e-4), ("s6", "dn_inc", -1, 1e-10),
         ("s6", "up_dur", 1, 1e-3), ("s6", "gamma", 1, 1e-10), ("s1", "sigma_last20", 1, 1e-10)]
held = 0
for a, b, sg, bd in want4:
    r, p, n = fcx.pearson(PK[a], PK[b])
    ok = np.sign(r) == sg and p < bd
    held += ok
    print(f"{a} vs {b}: R {r:.3f}, p {p:.2g}, n {n}; paper sign {'+' if sg > 0 else '-'}, bound {bd:g}: {'holds' if ok else 'fails'}"
          f"{' (p < 0.05, sign right)' if not ok and np.sign(r) == sg and p < 0.05 else ''}")
print(f"held {held} of 6")
print("other combinations (paper: not significant):")
for a, b in [("s1", "delta"), ("s1", "dn_dur"), ("s1", "dn_inc"), ("s1", "up_dur"), ("s1", "gamma"), ("s6", "sigma_last20")]:
    r, p, n = fcx.pearson(PK[a], PK[b])
    print(f"   {a} vs {b}: R {r:.3f}, p {p:.2g}")
r, p, n = fcx.pearson(PK["s1"], PK["sigma_last20_in"])
print(f"variant, s1 vs sigma from 10 s windows inside the last 20 s only: R {r:.3f}, p {p:.2g}, n {n}")
held_b = 0
for a, b, sg, bd in want4:
    r, p, n = fcx.pearson(PK[a + "b"], PK[b])
    ok = np.sign(r) == sg and p < bd
    held_b += ok
    print(f"M2 reading (b), per-cycle sextiles: {a} vs {b}: R {r:.3f}, p {p:.2g}, n {n}: {'holds' if ok else 'fails'}")
print(f"M2 reading (b): held {held_b} of 6")

s6 vs delta: R -0.200, p 1.6e-06, n 566; paper sign -, bound 1e-10: fails (p < 0.05, sign right)
s6 vs dn_dur: R -0.387, p 1.7e-21, n 560; paper sign -, bound 0.0001: holds
s6 vs dn_inc: R -0.300, p 2.9e-13, n 566; paper sign -, bound 1e-10: holds
s6 vs up_dur: R 0.248, p 2.8e-09, n 560; paper sign +, bound 0.001: holds
s6 vs gamma: R 0.358, p 1.4e-18, n 566; paper sign +, bound 1e-10: holds
s1 vs sigma_last20: R 0.215, p 7.7e-08, n 611; paper sign +, bound 1e-10: fails (p < 0.05, sign right)
held 4 of 6
other combinations (paper: not significant):
   s1 vs delta: R 0.049, p 0.23
   s1 vs dn_dur: R -0.018, p 0.66
   s1 vs dn_inc: R 0.014, p 0.73
   s1 vs up_dur: R 0.110, p 0.01
   s1 vs gamma: R 0.026, p 0.53
   s6 vs sigma_last20: R -0.093, p 0.026
variant, s1 vs sigma from 10 s windows inside the last 20 s only: R 0.199, p 6.9e-07, n 611
M2 reading (b), per-cycle sextiles: s6 vs delta: R -0.223, p 3.2e-08, n 600: fails
M2 reading (b), per-cycle sextiles: s6 vs dn_dur: R -0.355, p 4.4e-19, n 594: holds
M2 reading (b), per-cycle sextiles: s6 vs dn_inc: R -0.273, p 9.9e-12, n 600: holds
M2 reading (b), per-cycle sextiles: s6 vs up_dur: R 0.222, p 5.6e-08, n 587: holds
M2 reading (b), per-cycle sextiles: s6 vs gamma: R 0.322, p 5.9e-16, n 600: holds
M2 reading (b), per-cycle sextiles: s1 vs sigma_last20: R 0.062, p 0.12, n 639: fails
M2 reading (b): held 4 of 6

def fig4e(sx, label):
    groups = {"fastest": lambda s, t: sext_cells(6, sx, t)(s), "pI": lambda s, t: s["I"], "slowest": lambda s, t: sext_cells(1, sx, t)(s)}
    U_i = {g: {0: [], 1: [], 2: []} for g in groups}
    CC = {}   # (group, session, cell, third) -> [spikes, UP time], summed over packets
    for s in S:
        up = s["updown"]["UP"]
        if len(up) == 0:
            continue
        for a, b in fcx.state_intervals(s, "SWSPacket"):
            m = (up[:, 0] >= a) & (up[:, 1] <= b)
            if m.sum() == 0:
                continue
            u = up[m]
            third = np.minimum(((u[:, 0] - a) / (b - a) * 3).astype(int), 2)
            d = u[:, 1] - u[:, 0]
            for g, f in groups.items():
                cells = f(s, a)
                if len(cells) == 0:
                    continue
                C = np.array([fcx.count_in(s["spikes"][i], u) for i in cells])   # cells x UP states
                rate = C.sum(0) / len(cells) / d
                for k in range(3):
                    U_i[g][k] += list(rate[third == k])
                    if (third == k).any():
                        for ci, cell in enumerate(cells):
                            v = CC.setdefault((g, s["session"], cell, k), [0.0, 0.0])
                            v[0] += C[ci, third == k].sum()
                            v[1] += d[third == k].sum()
    # unit (ii): one entry per cell, its within-UP rate in each third over all its packets
    U_ii = {g: {k: [v[0] / v[1] for key, v in CC.items() if key[0] == g and key[3] == k] for k in range(3)} for g in groups}
    want_dir = {"fastest": (-1, 1e-8), "pI": (1, 1e-3), "slowest": (-1, 0.003)}
    for unit, U in (("i", U_i), ("ii", U_ii)):
        for rd, k2 in (("a", 1), ("b", 2)):
            out = []
            held = 0
            for g in groups:
                x0, x1 = np.array(U[g][0]), np.array(U[g][k2])
                p = stats.mannwhitneyu(x1, x0).pvalue
                sg = np.sign(np.mean(x1) - np.mean(x0))
                ok = sg == want_dir[g][0] and p < want_dir[g][1]
                held += ok
                out.append(f"{g}: {np.mean(x0):.2f} -> {np.mean(x1):.2f} Hz (medians {np.median(x0):.2f} -> {np.median(x1):.2f}), p {p:.2g}, n {len(x0)}/{len(x1)}, {'holds' if ok else 'fails'}")
            print(f"{label}, unit ({unit}), reading ({rd}): held {held} of 3; " + "; ".join(out))

fig4e(None, "M2 reading (a)")
fig4e(sx_b, "M2 reading (b), per-cycle sextiles")

M2 reading (a), unit (i), reading (a): held 1 of 3; fastest: 3.43 -> 3.39 Hz (medians 2.71 -> 2.61), p 5.8e-05, n 12679/14950, fails; pI: 13.11 -> 13.65 Hz (medians 10.06 -> 10.33), p 0.00078, n 12899/15208, holds; slowest: 0.29 -> 0.31 Hz (medians 0.00 -> 0.00), p 0.64, n 11045/13029, fails
M2 reading (a), unit (i), reading (b): held 1 of 3; fastest: 3.43 -> 3.42 Hz (medians 2.71 -> 2.66), p 0.0069, n 12679/14395, fails; pI: 13.11 -> 14.01 Hz (medians 10.06 -> 10.44), p 3.1e-05, n 12899/14638, holds; slowest: 0.29 -> 0.33 Hz (medians 0.00 -> 0.00), p 5.3e-07, n 11045/12802, fails
M2 reading (a), unit (ii), reading (a): held 0 of 3; fastest: 3.06 -> 2.96 Hz (medians 2.59 -> 2.44), p 0.4, n 166/166, fails; pI: 8.38 -> 8.68 Hz (medians 7.03 -> 7.30), p 0.8, n 126/126, fails; slowest: 0.27 -> 0.30 Hz (medians 0.20 -> 0.20), p 0.36, n 165/165, fails
M2 reading (a), unit (ii), reading (b): held 0 of 3; fastest: 3.06 -> 2.92 Hz (medians 2.59 -> 2.45), p 0.44, n 166/166, fails; pI: 8.38 -> 9.13 Hz (medians 7.03 -> 7.47), p 0.58, n 126/126, fails; slowest: 0.27 -> 0.32 Hz (medians 0.20 -> 0.21), p 0.11, n 165/165, fails

M2 reading (b), per-cycle sextiles, unit (i), reading (a): held 1 of 3; fastest: 3.16 -> 3.08 Hz (medians 2.74 -> 2.63), p 1.6e-06, n 13012/15314, fails; pI: 13.11 -> 13.65 Hz (medians 10.06 -> 10.33), p 0.00078, n 12899/15208, holds; slowest: 0.29 -> 0.31 Hz (medians 0.00 -> 0.00), p 0.29, n 12027/14178, fails
M2 reading (b), per-cycle sextiles, unit (i), reading (b): held 1 of 3; fastest: 3.16 -> 3.10 Hz (medians 2.74 -> 2.66), p 3.6e-05, n 13012/14752, fails; pI: 13.11 -> 14.01 Hz (medians 10.06 -> 10.44), p 3.1e-05, n 12899/14638, holds; slowest: 0.29 -> 0.34 Hz (medians 0.00 -> 0.00), p 5.5e-10, n 12027/13694, fails
M2 reading (b), per-cycle sextiles, unit (ii), reading (a): held 0 of 3; fastest: 2.98 -> 2.86 Hz (medians 2.48 -> 2.38), p 0.37, n 197/197, fails; pI: 8.38 -> 8.68 Hz (medians 7.03 -> 7.30), p 0.8, n 126/126, fails; slowest: 0.33 -> 0.37 Hz (medians 0.21 -> 0.23), p 0.31, n 209/209, fails
M2 reading (b), per-cycle sextiles, unit (ii), reading (b): held 0 of 3; fastest: 2.98 -> 2.78 Hz (medians 2.48 -> 2.35), p 0.33, n 197/197, fails; pI: 8.38 -> 9.13 Hz (medians 7.03 -> 7.47), p 0.58, n 126/126, fails; slowest: 0.33 -> 0.40 Hz (medians 0.21 -> 0.24), p 0.12, n 209/209, fails

tr_s5 = {}
for m2, mk in (("a", lambda k: fcx.m_rate(sext_cells(k))), ("b", lambda k: m_rate_sx(k, sx_b))):
    for scope, lab in (("wake", (1,)), ("rem", (4,))):
        for k in range(1, 7):
            r = fcx.trend(S, tabs, scope, mk(k), "c", labels=lab)
            tr_s5[(m2, scope, k)] = r
            print(f"M2 reading ({m2}), {scope} sextile {k} (1 lowest): R {r['r']:.3f}, p {r['p']:.2g}, n {r['n']}")
    wake_ok = all(tr_s5[(m2, "wake", k)]["r"] > 0 and tr_s5[(m2, "wake", k)]["p"] < 0.05 for k in (5, 6)) and \
              all(tr_s5[(m2, "wake", k)]["r"] < 0 and tr_s5[(m2, "wake", k)]["p"] < 0.05 for k in (1, 2))
    rem_ok = all(tr_s5[(m2, "rem", k)]["r"] > 0 and tr_s5[(m2, "rem", k)]["p"] < 0.05 for k in (1, 2, 4, 5, 6))
    print(f"M2 reading ({m2}): WAKE part holds: {wake_ok}; REM part holds: {rem_ok}")

M2 reading (a), wake sextile 1 (1 lowest): R -0.286, p 0.044, n 50

M2 reading (a), wake sextile 2 (1 lowest): R 0.603, p 3.5e-06, n 50

M2 reading (a), wake sextile 3 (1 lowest): R 0.153, p 0.29, n 50

M2 reading (a), wake sextile 4 (1 lowest): R 0.281, p 0.048, n 50

M2 reading (a), wake sextile 5 (1 lowest): R 0.471, p 0.00055, n 50

M2 reading (a), wake sextile 6 (1 lowest): R 0.107, p 0.46, n 50

M2 reading (a), rem sextile 1 (1 lowest): R -0.245, p 0.086, n 50

M2 reading (a), rem sextile 2 (1 lowest): R 0.443, p 0.0013, n 50

M2 reading (a), rem sextile 3 (1 lowest): R 0.358, p 0.011, n 50

M2 reading (a), rem sextile 4 (1 lowest): R 0.394, p 0.0047, n 50

wake_rate = WR.set_index(["session", "cell"])["WAKE"].to_dict()

def triplet_changes(pattern, sessions=None, **kw):
    rows = []
    for s in (S if sessions is None else sessions):
        for k, t in enumerate(fcx.triplets(s, pattern, **kw)):
            for i in list(s["E"]) + list(s["I"]):
                x = s["spikes"][i]
                r1 = fcx.rate_in(x, np.array([t["first"]]))
                r2 = fcx.rate_in(x, np.array([t["second"]]))
                rows.append(dict(session=s["session"], rat=s["rat"], trip=k, kind=t["kind"], cell=i,
                                 type="E" if i in s["E"] else "I",
                                 sext=sext.get((s["session"], i), 0) if i in s["E"] else 0,
                                 sext_b=sx_b(s, i, t["first"][0]) if i in s["E"] else 0,
                                 wake=wake_rate[(s["session"], i)],
                                 r1=r1, r2=r2, d=r2 - r1,
                                 dlog=np.log10(r2) - np.log10(r1) if r1 > 0 and r2 > 0 else np.nan,
                                 mid0=t["mid"][0], mid1=t["mid"][1], p0=t["first"][0], p1=t["first"][1],
                                 q0=t["second"][0], q1=t["second"][1]))
    return pd.DataFrame(rows)

def sextile_tests(df, label):
    e = df[df.type == "E"]
    cell = e.groupby(["session", "cell", "sext"]).agg(d=("d", "mean")).reset_index()
    res = {}
    for k in range(1, 7):
        v_b = cell.d[cell.sext == k]
        v_a = e.d[e.sext == k]
        res[k] = dict(mean=v_b.mean(), n=len(v_b), w_b=stats.wilcoxon(v_b).pvalue if len(v_b) > 5 and np.any(v_b != 0) else np.nan,
                      t_b=stats.ttest_1samp(v_b, 0).pvalue if len(v_b) > 2 else np.nan,
                      w_a=stats.wilcoxon(v_a).pvalue if len(v_a) > 5 and np.any(v_a != 0) else np.nan,
                      t_a=stats.ttest_1samp(v_a, 0).pvalue if len(v_a) > 2 else np.nan, n_a=len(v_a))
        print(f"{label} sextile {k} (1 lowest): mean change {v_b.mean():+.3f} Hz, cells {len(v_b)}; signed-rank p (b) {res[k]['w_b']:.2g}, "
              f"t p (b) {res[k]['t_b']:.2g}; per cell-triplet n {len(v_a)}, mean {v_a.mean():+.3f} Hz: signed-rank p (a) {res[k]['w_a']:.2g}, t p (a) {res[k]['t_a']:.2g}")
    return res

def s5_corr(df, label):
    e = df[df.type == "E"].dropna(subset=["dlog"])
    cell = e.groupby(["session", "cell"]).agg(wake=("wake", "first"), dlog=("dlog", "mean"))
    ok = cell.wake > 0
    ra = fcx.pearson(np.log10(cell.wake[ok]), cell.dlog[ok])
    rb = fcx.pearson(cell.wake, cell.dlog)
    print(f"{label} S5 correlation, log rate change vs WAKE rate: (a) log10 WAKE R {ra[0]:.3f} p {ra[1]:.2g} n {ra[2]}; (b) linear R {rb[0]:.3f} p {rb[1]:.2g} n {rb[2]}")
    return ra, rb

def per_triplet(df, cells_sel):
    e = df[cells_sel(df)]
    return e.groupby(["session", "trip"]).agg(d=("d", "mean"), mid0=("mid0", "first"), mid1=("mid1", "first"),
                                              p0=("p0", "first"), p1=("p1", "first"), q0=("q0", "first"), q1=("q1", "first")).reset_index()

TR = {}
for pat in ("MPM", "RPR", "RPR_after", "PMP", "PRP"):
    TR[pat] = triplet_changes(pat)
for pat in ("MPM", "PMP"):
    TR[pat + "_b"] = triplet_changes(pat, sessions=SB)   # MA reading (b), see ma-statistics
TR["PWP_a"] = triplet_changes("PWP", wake_min=None)
TR["PWP_a"] = TR["PWP_a"][TR["PWP_a"].kind == "WAKE"]
TR["PWP_b"] = triplet_changes("PWP", wake_min=40)
for k, v in TR.items():
    print(f"{k}: {v.groupby(['session', 'trip']).ngroups} triplets, {v[v.type == 'E'].groupby(['session', 'cell']).ngroups} pE cells")

MPM: 169 triplets, 941 pE cells
RPR: 25 triplets, 59 pE cells
RPR_after: 80 triplets, 660 pE cells
PMP: 264 triplets, 977 pE cells
PRP: 110 triplets, 789 pE cells
MPM_b: 234 triplets, 961 pE cells
PMP_b: 264 triplets, 977 pE cells
PWP_a: 8 triplets, 315 pE cells
PWP_b: 52 triplets, 866 pE cells

res_mpm = sextile_tests(TR["MPM"], "MA-packet-MA")
s5c = s5_corr(TR["MPM"], "MA-packet-MA")
print("paper: top two sextiles decrease, p < 0.01; S5C R = -0.149, p < 1e-5")
print("MA reading (b):", TR["MPM_b"].groupby(["session", "trip"]).ngroups, "MA-packet-MA triplets")
_ = sextile_tests(TR["MPM_b"], "MA-packet-MA, MA reading (b)")
_ = s5_corr(TR["MPM_b"], "MA-packet-MA, MA reading (b)")
_ = sextile_tests(TR["MPM"].assign(sext=TR["MPM"].sext_b), "MA-packet-MA, M2 reading (b)")

MA-packet-MA sextile 1 (1 lowest): mean change +0.024 Hz, cells 145; signed-rank p (b) 0.48, t p (b) 0.23; per cell-triplet n 1071, mean -0.003 Hz: signed-rank p (a) 0.96, t p (a) 0.77
MA-packet-MA sextile 2 (1 lowest): mean change -0.029 Hz, cells 159; signed-rank p (b) 0.11, t p (b) 0.1; per cell-triplet n 1152, mean -0.023 Hz: signed-rank p (a) 0.27, t p (a) 0.05
MA-packet-MA sextile 3 (1 lowest): mean change -0.035 Hz, cells 157; signed-rank p (b) 0.25, t p (b) 0.36; per cell-triplet n 1060, mean -0.024 Hz: signed-rank p (a) 0.37, t p (a) 0.26
MA-packet-MA sextile 4 (1 lowest): mean change -0.005 Hz, cells 155; signed-rank p (b) 0.37, t p (b) 0.89; per cell-triplet n 1160, mean -0.013 Hz: signed-rank p (a) 0.24, t p (a) 0.59
MA-packet-MA sextile 5 (1 lowest): mean change -0.011 Hz, cells 159; signed-rank p (b) 0.63, t p (b) 0.8; per cell-triplet n 1190, mean -0.031 Hz: signed-rank p (a) 0.24, t p (a) 0.31
MA-packet-MA sextile 6 (1 lowest): mean change -0.010 Hz, cells 166; signed-rank p (b) 0.27, t p (b) 0.91; per cell-triplet n 1129, mean -0.035 Hz: signed-rank p (a) 0.26, t p (a) 0.51
MA-packet-MA S5 correlation, log rate change vs WAKE rate: (a) log10 WAKE R 0.016 p 0.65 n 843; (b) linear R 0.006 p 0.86 n 844
paper: top two sextiles decrease, p < 0.01; S5C R = -0.149, p < 1e-5
MA reading (b): 234 MA-packet-MA triplets
MA-packet-MA, MA reading (b) sextile 1 (1 lowest): mean change +0.047 Hz, cells 158; signed-rank p (b) 0.043, t p (b) 0.093; per cell-triplet n 1532, mean +0.017 Hz: signed-rank p (a) 0.1, t p (a) 0.13
MA-packet-MA, MA reading (b) sextile 2 (1 lowest): mean change -0.016 Hz, cells 161; signed-rank p (b) 0.74, t p (b) 0.32; per cell-triplet n 1580, mean -0.005 Hz: signed-rank p (a) 0.64, t p (a) 0.62
MA-packet-MA, MA reading (b) sextile 3 (1 lowest): mean change -0.010 Hz, cells 158; signed-rank p (b) 0.69, t p (b) 0.8; per cell-triplet n 1489, mean +0.017 Hz: signed-rank p (a) 0.75, t p (a) 0.37
MA-packet-MA, MA reading (b) sextile 4 (1 lowest): mean change -0.016 Hz, cells 158; signed-rank p (b) 0.33, t p (b) 0.58; per cell-triplet n 1644, mean +0.012 Hz: signed-rank p (a) 0.73, t p (a) 0.59
MA-packet-MA, MA reading (b) sextile 5 (1 lowest): mean change +0.014 Hz, cells 160; signed-rank p (b) 0.28, t p (b) 0.68; per cell-triplet n 1652, mean +0.005 Hz: signed-rank p (a) 0.86, t p (a) 0.86
MA-packet-MA, MA reading (b) sextile 6 (1 lowest): mean change -0.012 Hz, cells 166; signed-rank p (b) 0.42, t p (b) 0.9; per cell-triplet n 1502, mean -0.016 Hz: signed-rank p (a) 0.36, t p (a) 0.74
MA-packet-MA, MA reading (b) S5 correlation, log rate change vs WAKE rate: (a) log10 WAKE R -0.054 p 0.11 n 868; (b) linear R -0.034 p 0.32 n 869
MA-packet-MA, M2 reading (b) sextile 1 (1 lowest): mean change +0.021 Hz, cells 175; signed-rank p (b) 0.53, t p (b) 0.22; per cell-triplet n 1109, mean -0.009 Hz: signed-rank p (a) 0.97, t p (a) 0.37
MA-packet-MA, M2 reading (b) sextile 2 (1 lowest): mean change -0.010 Hz, cells 211; signed-rank p (b) 0.29, t p (b) 0.55; per cell-triplet n 1162, mean -0.015 Hz: signed-rank p (a) 0.16, t p (a) 0.22
MA-packet-MA, M2 reading (b) sextile 3 (1 lowest): mean change -0.046 Hz, cells 215; signed-rank p (b) 0.098, t p (b) 0.14; per cell-triplet n 1042, mean -0.019 Hz: signed-rank p (a) 0.39, t p (a) 0.34
MA-packet-MA, M2 reading (b) sextile 4 (1 lowest): mean change -0.063 Hz, cells 210; signed-rank p (b) 0.085, t p (b) 0.12; per cell-triplet n 1092, mean -0.025 Hz: signed-rank p (a) 0.25, t p (a) 0.32
MA-packet-MA, M2 reading (b) sextile 5 (1 lowest): mean change +0.000 Hz, cells 218; signed-rank p (b) 0.7, t p (b) 1; per cell-triplet n 1197, mean -0.003 Hz: signed-rank p (a) 0.52, t p (a) 0.93
MA-packet-MA, M2 reading (b) sextile 6 (1 lowest): mean change -0.040 Hz, cells 191; signed-rank p (b) 0.15, t p (b) 0.65; per cell-triplet n 1160, mean -0.059 Hz: signed-rank p (a) 0.12, t p (a) 0.26

# Adopted reading (assumption): REM n, packet, REM n+1 with any MA or WakeInterruption allowed
# between REM n and the packet, the tolerance M6 gives for the arousal that follows REM.
for key, label in (("RPR_after", "arousal after REM n allowed (adopted)"), ("RPR", "all adjacent")):
    df = TR[key]
    e = df[df.type == "E"]
    print(f"REM-packet-REM, {label}: {df.groupby(['session', 'trip']).ngroups} triplets in {df.session.nunique()} sessions, "
          f"{e.groupby(['session', 'cell']).ngroups} pE cells")
    res = s5_corr(df, f"REM-packet-REM, {label}")
    if key == "RPR_after":
        s5d = res
print("paper S5D R = -0.395, p < 1e-5, about 136 cells")
# Variant: REM, packet, REM with any MA or WakeInterruption between them on either side.
rpr_b = triplet_changes("RPR_bridge")
print("variant with interruptions allowed on both sides:", rpr_b.groupby(["session", "trip"]).ngroups, "triplets,",
      rpr_b[rpr_b.type == "E"].groupby(["session", "cell"]).ngroups, "pE cells")
s5d_b = s5_corr(rpr_b, "REM-packet-REM (both sides)")

REM-packet-REM, arousal after REM n allowed (adopted): 80 triplets in 17 sessions, 660 pE cells
REM-packet-REM, arousal after REM n allowed (adopted) S5 correlation, log rate change vs WAKE rate: (a) log10 WAKE R -0.077 p 0.054 n 627; (b) linear R -0.097 p 0.015 n 627
REM-packet-REM, all adjacent: 25 triplets in 3 sessions, 59 pE cells
REM-packet-REM, all adjacent S5 correlation, log rate change vs WAKE rate: (a) log10 WAKE R 0.019 p 0.89 n 57; (b) linear R -0.055 p 0.69 n 57
paper S5D R = -0.395, p < 1e-5, about 136 cells

variant with interruptions allowed on both sides: 84 triplets, 678 pE cells

REM-packet-REM (both sides) S5 correlation, log rate change vs WAKE rate: (a) log10 WAKE R -0.071 p 0.071 n 646; (b) linear R -0.103 p 0.0089 n 646

tabd = {s["session"]: t for s, t in zip(S, tabs)}
pt = per_triplet(TR["MPM"], lambda d: d.sext == 1)
pt["sigma20"] = [tabd[r.session]["zband"]["sigma"][int(r.mid1) - 20:int(r.mid1)].mean() for r in pt.itertuples()]
r, p, n = fcx.pearson(pt.sigma20, pt.d)
print(f"lowest sextile change vs last-20 s sigma: R {r:.3f}, p {p:.2g}, n {n} triplets (paper R 0.181, p < 0.01)")
pt["sigma20_in"] = [tabd[r.session]["zband"]["sigma"][int(r.mid1) - 15:int(r.mid1) - 4].mean() for r in pt.itertuples()]
r, p, n = fcx.pearson(pt.sigma20_in, pt.d)
print(f"variant, 10 s windows inside the last 20 s only: R {r:.3f}, p {p:.2g}, n {n}")
ptb = per_triplet(TR["MPM_b"], lambda d: d.sext == 1)
ptb["sigma20"] = [tabd[r.session]["zband"]["sigma"][int(r.mid1) - 20:int(r.mid1)].mean() for r in ptb.itertuples()]
r, p, n = fcx.pearson(ptb.sigma20, ptb.d)
print(f"MA reading (b): R {r:.3f}, p {p:.2g}, n {n} triplets")
ptc = per_triplet(TR["MPM"], lambda d: d.sext_b == 1)
ptc["sigma20"] = [tabd[r.session]["zband"]["sigma"][int(r.mid1) - 20:int(r.mid1)].mean() for r in ptc.itertuples()]
r, p, n = fcx.pearson(ptc.sigma20, ptc.d)
print(f"M2 reading (b), per-cycle sextiles: R {r:.3f}, p {p:.2g}, n {n} triplets")

lowest sextile change vs last-20 s sigma: R -0.042, p 0.61, n 150 triplets (paper R 0.181, p < 0.01)
variant, 10 s windows inside the last 20 s only: R -0.018, p 0.83, n 150
MA reading (b): R -0.013, p 0.86, n 209 triplets
M2 reading (b), per-cycle sextiles: R -0.053, p 0.52, n 153 triplets

mot = {s["session"]: s["motion"]["motion"] for s in S}
for rd in ("a", "b"):
    df = TR[f"PWP_{rd}"]
    print(f"reading ({rd}): {df.groupby(['session', 'trip']).ngroups} triplets (paper 24)")
    res = sextile_tests(df, f"packet-WAKE-packet ({rd})")
    pt = per_triplet(df, lambda d: d.sext == 6)
    pt["motion"] = [np.nanmean(mot[r.session][int(r.mid0):int(r.mid1)]) for r in pt.itertuples()]
    r, p, n = fcx.pearson(pt.motion, pt.d)
    print(f"   fastest sextile change vs mean motion: R {r:.3f}, p {p:.2g}, n {n}")
    # variant: motion z-scored per session first (the files are not z-scored, see the data check)
    mz = {k: (v - np.nanmean(v)) / np.nanstd(v) for k, v in mot.items()}
    pt["motion_z"] = [np.nanmean(mz[r_.session][int(r_.mid0):int(r_.mid1)]) for r_ in pt.itertuples()]
    r, p, n = fcx.pearson(pt.motion_z, pt.d)
    print(f"   variant, motion z-scored per session: R {r:.3f}, p {p:.2g}, n {n}")
    s5_corr(df, f"packet-WAKE-packet ({rd})")
# M2 reading (b), per-cycle sextiles
for rd in ("a", "b"):
    df = TR[f"PWP_{rd}"]
    _ = sextile_tests(df.assign(sext=df.sext_b), f"packet-WAKE-packet ({rd}), M2 reading (b)")
    pt = per_triplet(df, lambda d: d.sext_b == 6)
    for nm, mm in (("raw", mot), ("z-scored", mz)):
        r, p, n = fcx.pearson([np.nanmean(mm[r_.session][int(r_.mid0):int(r_.mid1)]) for r_ in pt.itertuples()], pt.d)
        print(f"   M2 reading (b), fastest sextile change vs mean motion ({nm}): R {r:.3f}, p {p:.2g}, n {n}")

reading (a): 8 triplets (paper 24)
packet-WAKE-packet (a) sextile 1 (1 lowest): mean change -0.001 Hz, cells 39; signed-rank p (b) 0.96, t p (b) 0.97; per cell-triplet n 39, mean -0.001 Hz: signed-rank p (a) 0.96, t p (a) 0.97
packet-WAKE-packet (a) sextile 2 (1 lowest): mean change +0.039 Hz, cells 56; signed-rank p (b) 0.21, t p (b) 0.15; per cell-triplet n 56, mean +0.039 Hz: signed-rank p (a) 0.21, t p (a) 0.15
packet-WAKE-packet (a) sextile 3 (1 lowest): mean change +0.002 Hz, cells 62; signed-rank p (b) 0.72, t p (b) 0.95; per cell-triplet n 62, mean +0.002 Hz: signed-rank p (a) 0.72, t p (a) 0.95
packet-WAKE-packet (a) sextile 4 (1 lowest): mean change +0.084 Hz, cells 51; signed-rank p (b) 0.17, t p (b) 0.14; per cell-triplet n 51, mean +0.084 Hz: signed-rank p (a) 0.17, t p (a) 0.14
packet-WAKE-packet (a) sextile 5 (1 lowest): mean change -0.064 Hz, cells 50; signed-rank p (b) 0.85, t p (b) 0.65; per cell-triplet n 50, mean -0.064 Hz: signed-rank p (a) 0.85, t p (a) 0.65
packet-WAKE-packet (a) sextile 6 (1 lowest): mean change +0.092 Hz, cells 57; signed-rank p (b) 0.21, t p (b) 0.35; per cell-triplet n 57, mean +0.092 Hz: signed-rank p (a) 0.21, t p (a) 0.35
   fastest sextile change vs mean motion: R 0.436, p 0.33, n 7
   variant, motion z-scored per session: R 0.869, p 0.011, n 7
packet-WAKE-packet (a) S5 correlation, log rate change vs WAKE rate: (a) log10 WAKE R 0.010 p 0.86 n 302; (b) linear R 0.005 p 0.93 n 302
reading (b): 52 triplets (paper 24)
packet-WAKE-packet (b) sextile 1 (1 lowest): mean change +0.013 Hz, cells 132; signed-rank p (b) 0.33, t p (b) 0.12; per cell-triplet n 306, mean +0.004 Hz: signed-rank p (a) 0.25, t p (a) 0.67
packet-WAKE-packet (b) sextile 2 (1 lowest): mean change +0.002 Hz, cells 147; signed-rank p (b) 0.047, t p (b) 0.94; per cell-triplet n 312, mean +0.013 Hz: signed-rank p (a) 0.059, t p (a) 0.37
packet-WAKE-packet (b) sextile 3 (1 lowest): mean change +0.012 Hz, cells 152; signed-rank p (b) 0.11, t p (b) 0.64; per cell-triplet n 305, mean +0.031 Hz: signed-rank p (a) 0.054, t p (a) 0.073
packet-WAKE-packet (b) sextile 4 (1 lowest): mean change +0.051 Hz, cells 152; signed-rank p (b) 0.068, t p (b) 0.051; per cell-triplet n 322, mean +0.055 Hz: signed-rank p (a) 0.0089, t p (a) 0.0045
packet-WAKE-packet (b) sextile 5 (1 lowest): mean change -0.073 Hz, cells 145; signed-rank p (b) 0.43, t p (b) 0.18; per cell-triplet n 293, mean -0.010 Hz: signed-rank p (a) 0.63, t p (a) 0.77
packet-WAKE-packet (b) sextile 6 (1 lowest): mean change +0.063 Hz, cells 138; signed-rank p (b) 0.15, t p (b) 0.14; per cell-triplet n 291, mean +0.051 Hz: signed-rank p (a) 0.19, t p (a) 0.14
   fastest sextile change vs mean motion: R 0.124, p 0.43, n 43
   variant, motion z-scored per session: R 0.248, p 0.11, n 43
packet-WAKE-packet (b) S5 correlation, log rate change vs WAKE rate: (a) log10 WAKE R 0.011 p 0.75 n 848; (b) linear R -0.002 p 0.95 n 850
packet-WAKE-packet (a), M2 reading (b) sextile 1 (1 lowest): mean change +0.002 Hz, cells 48; signed-rank p (b) 0.74, t p (b) 0.92; per cell-triplet n 48, mean +0.002 Hz: signed-rank p (a) 0.74, t p (a) 0.92
packet-WAKE-packet (a), M2 reading (b) sextile 2 (1 lowest): mean change +0.038 Hz, cells 54; signed-rank p (b) 0.57, t p (b) 0.26; per cell-triplet n 54, mean +0.038 Hz: signed-rank p (a) 0.57, t p (a) 0.26
packet-WAKE-packet (a), M2 reading (b) sextile 3 (1 lowest): mean change +0.055 Hz, cells 60; signed-rank p (b) 0.086, t p (b) 0.16; per cell-triplet n 60, mean +0.055 Hz: signed-rank p (a) 0.086, t p (a) 0.16
packet-WAKE-packet (a), M2 reading (b) sextile 4 (1 lowest): mean change +0.075 Hz, cells 49; signed-rank p (b) 0.36, t p (b) 0.28; per cell-triplet n 49, mean +0.075 Hz: signed-rank p (a) 0.36, t p (a) 0.28
packet-WAKE-packet (a), M2 reading (b) sextile 5 (1 lowest): mean change +0.076 Hz, cells 47; signed-rank p (b) 0.46, t p (b) 0.44; per cell-triplet n 47, mean +0.076 Hz: signed-rank p (a) 0.46, t p (a) 0.44
packet-WAKE-packet (a), M2 reading (b) sextile 6 (1 lowest): mean change -0.070 Hz, cells 57; signed-rank p (b) 0.99, t p (b) 0.59; per cell-triplet n 57, mean -0.070 Hz: signed-rank p (a) 0.99, t p (a) 0.59
   M2 reading (b), fastest sextile change vs mean motion (raw): R 0.451, p 0.31, n 7
   M2 reading (b), fastest sextile change vs mean motion (z-scored): R 0.924, p 0.0029, n 7
packet-WAKE-packet (b), M2 reading (b) sextile 1 (1 lowest): mean change +0.011 Hz, cells 152; signed-rank p (b) 0.28, t p (b) 0.19; per cell-triplet n 324, mean +0.005 Hz: signed-rank p (a) 0.3, t p (a) 0.57
packet-WAKE-packet (b), M2 reading (b) sextile 2 (1 lowest): mean change +0.023 Hz, cells 171; signed-rank p (b) 0.047, t p (b) 0.17; per cell-triplet n 307, mean +0.020 Hz: signed-rank p (a) 0.16, t p (a) 0.12
packet-WAKE-packet (b), M2 reading (b) sextile 3 (1 lowest): mean change +0.044 Hz, cells 178; signed-rank p (b) 0.01, t p (b) 0.032; per cell-triplet n 302, mean +0.041 Hz: signed-rank p (a) 0.006, t p (a) 0.014
packet-WAKE-packet (b), M2 reading (b) sextile 4 (1 lowest): mean change +0.031 Hz, cells 175; signed-rank p (b) 0.014, t p (b) 0.38; per cell-triplet n 308, mean +0.053 Hz: signed-rank p (a) 0.0014, t p (a) 0.037
packet-WAKE-packet (b), M2 reading (b) sextile 5 (1 lowest): mean change -0.007 Hz, cells 168; signed-rank p (b) 0.86, t p (b) 0.83; per cell-triplet n 288, mean +0.016 Hz: signed-rank p (a) 0.54, t p (a) 0.56
packet-WAKE-packet (b), M2 reading (b) sextile 6 (1 lowest): mean change -0.007 Hz, cells 153; signed-rank p (b) 0.68, t p (b) 0.89; per cell-triplet n 300, mean +0.012 Hz: signed-rank p (a) 0.67, t p (a) 0.76
   M2 reading (b), fastest sextile change vs mean motion (raw): R 0.175, p 0.25, n 45
   M2 reading (b), fastest sextile change vs mean motion (z-scored): R 0.282, p 0.061, n 45

res_pmp = sextile_tests(TR["PMP"], "packet-MA-packet")
e = TR["PMP"][TR["PMP"].type == "E"]
cell = e.groupby(["session", "cell"]).agg(sext=("sext", "first"), d=("d", "mean"))
v2 = cell.d[cell.sext <= 2]
print(f"lowest two sextiles together (reading b): mean change {v2.mean():+.3f}, signed-rank p {stats.wilcoxon(v2).pvalue:.2g}, t p {stats.ttest_1samp(v2, 0).pvalue:.2g}, n {len(v2)}")
s5g = s5_corr(TR["PMP"], "packet-MA-packet")
print("paper: S5G R = -0.207, p < 1e-10")
lo_a = e[e.sext == 1].d
print(f"lowest sextile, one entry per cell and triplet (M1 a): mean change {lo_a.mean():+.4f} Hz, median {lo_a.median():+.4f}, "
      f"fraction positive {np.mean(lo_a > 0):.2f}, fraction zero {np.mean(lo_a == 0):.2f}, n {len(lo_a)}")
# MA reading (b): WakeInterruptions shorter than 40 s counted as MAs (see ma-statistics)
print("MA reading (b):", TR["PMP_b"].groupby(["session", "trip"]).ngroups, "packet-MA-packet triplets")
_ = sextile_tests(TR["PMP_b"], "packet-MA-packet, MA reading (b)")
_ = s5_corr(TR["PMP_b"], "packet-MA-packet, MA reading (b)")
# M2 reading (b), per-cycle sextiles
_ = sextile_tests(TR["PMP"].assign(sext=TR["PMP"].sext_b), "packet-MA-packet, M2 reading (b)")

packet-MA-packet sextile 1 (1 lowest): mean change +0.007 Hz, cells 162; signed-rank p (b) 0.36, t p (b) 0.068; per cell-triplet n 1554, mean +0.008 Hz: signed-rank p (a) 0.015, t p (a) 0.016

packet-MA-packet sextile 2 (1 lowest): mean change -0.006 Hz, cells 162; signed-rank p (b) 0.46, t p (b) 0.26; per cell-triplet n 1705, mean -0.001 Hz: signed-rank p (a) 0.68, t p (a) 0.8

packet-MA-packet sextile 3 (1 lowest): mean change -0.001 Hz, cells 164; signed-rank p (b) 0.11, t p (b) 0.9; per cell-triplet n 1628, mean -0.012 Hz: signed-rank p (a) 0.14, t p (a) 0.059
packet-MA-packet sextile 4 (1 lowest): mean change -0.023 Hz, cells 165; signed-rank p (b) 0.019, t p (b) 0.015; per cell-triplet n 1750, mean -0.028 Hz: signed-rank p (a) 0.00024, t p (a) 4.2e-05
packet-MA-packet sextile 5 (1 lowest): mean change -0.051 Hz, cells 163; signed-rank p (b) 8.8e-05, t p (b) 0.00011; per cell-triplet n 1727, mean -0.048 Hz: signed-rank p (a) 1.3e-08, t p (a) 6.5e-09
packet-MA-packet sextile 6 (1 lowest): mean change -0.140 Hz, cells 161; signed-rank p (b) 5.4e-08, t p (b) 1.6e-06; per cell-triplet n 1676, mean -0.080 Hz: signed-rank p (a) 1.9e-08, t p (a) 9.4e-10

lowest two sextiles together (reading b): mean change +0.000, signed-rank p 0.92, t p 0.88, n 324
packet-MA-packet S5 correlation, log rate change vs WAKE rate: (a) log10 WAKE R -0.174 p 4.9e-08 n 968; (b) linear R -0.133 p 3.1e-05 n 970
paper: S5G R = -0.207, p < 1e-10
lowest sextile, one entry per cell and triplet (M1 a): mean change +0.0075 Hz, median +0.0017, fraction positive 0.51, fraction zero 0.03, n 1554
MA reading (b): 264 packet-MA-packet triplets
packet-MA-packet, MA reading (b) sextile 1 (1 lowest): mean change +0.007 Hz, cells 162; signed-rank p (b) 0.36, t p (b) 0.068; per cell-triplet n 1554, mean +0.008 Hz: signed-rank p (a) 0.015, t p (a) 0.016
packet-MA-packet, MA reading (b) sextile 2 (1 lowest): mean change -0.006 Hz, cells 162; signed-rank p (b) 0.46, t p (b) 0.26; per cell-triplet n 1705, mean -0.001 Hz: signed-rank p (a) 0.68, t p (a) 0.8

packet-MA-packet, MA reading (b) sextile 3 (1 lowest): mean change -0.001 Hz, cells 164; signed-rank p (b) 0.11, t p (b) 0.9; per cell-triplet n 1628, mean -0.012 Hz: signed-rank p (a) 0.14, t p (a) 0.059
packet-MA-packet, MA reading (b) sextile 4 (1 lowest): mean change -0.023 Hz, cells 165; signed-rank p (b) 0.019, t p (b) 0.015; per cell-triplet n 1750, mean -0.028 Hz: signed-rank p (a) 0.00024, t p (a) 4.2e-05
packet-MA-packet, MA reading (b) sextile 5 (1 lowest): mean change -0.051 Hz, cells 163; signed-rank p (b) 8.8e-05, t p (b) 0.00011; per cell-triplet n 1727, mean -0.048 Hz: signed-rank p (a) 1.3e-08, t p (a) 6.5e-09
packet-MA-packet, MA reading (b) sextile 6 (1 lowest): mean change -0.140 Hz, cells 161; signed-rank p (b) 5.4e-08, t p (b) 1.6e-06; per cell-triplet n 1676, mean -0.080 Hz: signed-rank p (a) 1.9e-08, t p (a) 9.4e-10
packet-MA-packet, MA reading (b) S5 correlation, log rate change vs WAKE rate: (a) log10 WAKE R -0.174 p 4.9e-08 n 968; (b) linear R -0.133 p 3.1e-05 n 970
packet-MA-packet, M2 reading (b) sextile 1 (1 lowest): mean change +0.012 Hz, cells 202; signed-rank p (b) 0.022, t p (b) 0.011; per cell-triplet n 1646, mean +0.010 Hz: signed-rank p (a) 0.00088, t p (a) 0.0026
packet-MA-packet, M2 reading (b) sextile 2 (1 lowest): mean change +0.009 Hz, cells 219; signed-rank p (b) 0.96, t p (b) 0.22; per cell-triplet n 1685, mean +0.001 Hz: signed-rank p (a) 0.93, t p (a) 0.88
packet-MA-packet, M2 reading (b) sextile 3 (1 lowest): mean change -0.021 Hz, cells 225; signed-rank p (b) 0.058, t p (b) 0.022; per cell-triplet n 1588, mean -0.011 Hz: signed-rank p (a) 0.23, t p (a) 0.055
packet-MA-packet, M2 reading (b) sextile 4 (1 lowest): mean change -0.019 Hz, cells 226; signed-rank p (b) 0.011, t p (b) 0.047; per cell-triplet n 1675, mean -0.031 Hz: signed-rank p (a) 0.00014, t p (a) 1.7e-05
packet-MA-packet, M2 reading (b) sextile 5 (1 lowest): mean change -0.054 Hz, cells 223; signed-rank p (b) 4.1e-06, t p (b) 6.7e-06; per cell-triplet n 1735, mean -0.047 Hz: signed-rank p (a) 4.2e-08, t p (a) 3.1e-08
packet-MA-packet, M2 reading (b) sextile 6 (1 lowest): mean change -0.130 Hz, cells 189; signed-rank p (b) 6.6e-09, t p (b) 2.9e-07; per cell-triplet n 1711, mean -0.085 Hz: signed-rank p (a) 7.3e-10, t p (a) 5.5e-11

res_prp = sextile_tests(TR["PRP"], "packet-REM-packet")
pt = per_triplet(TR["PRP"], lambda d: d.type == "E")
r, p, n = fcx.pearson(pt.mid1 - pt.mid0, pt.d)
print(f"mean change of all pE cells vs REM duration: R {r:.3f}, p {p:.2g}, n {n} (paper R -0.22, p < 0.01)")
pt1 = per_triplet(TR["PRP"], lambda d: d.sext == 1)
pt6 = per_triplet(TR["PRP"], lambda d: d.sext == 6)
for lab, f in (("a, theta ratio", lambda tb, a, b: np.mean(tb["band"]["theta"][a:b] / tb["band"]["theta_wide"][a:b])),
               ("b, theta power", lambda tb, a, b: np.mean(tb["zband"]["theta"][a:b]))):
    for nm, ptx in (("slowest", pt1), ("fastest", pt6)):
        th = [f(tabd[r_.session], int(r_.mid0), int(r_.mid1)) for r_ in ptx.itertuples()]
        r, p, n = fcx.pearson(th, ptx.d)
        print(f"theta ({lab}) vs {nm} sextile change: R {r:.3f}, p {p:.2g}, n {n}")
s5f = s5_corr(TR["PRP"], "packet-REM-packet")
print("paper: S5F R = -0.038, p = 0.24 (no rate dependence)")
# M2 reading (b), per-cycle sextiles
_ = sextile_tests(TR["PRP"].assign(sext=TR["PRP"].sext_b), "packet-REM-packet, M2 reading (b)")
pt1b = per_triplet(TR["PRP"], lambda d: d.sext_b == 1)
for lab, f in (("a, theta ratio", lambda tb, a, b: np.mean(tb["band"]["theta"][a:b] / tb["band"]["theta_wide"][a:b])),
               ("b, theta power", lambda tb, a, b: np.mean(tb["zband"]["theta"][a:b]))):
    th = [f(tabd[r_.session], int(r_.mid0), int(r_.mid1)) for r_ in pt1b.itertuples()]
    r, p, n = fcx.pearson(th, pt1b.d)
    print(f"M2 reading (b), theta ({lab}) vs slowest sextile change: R {r:.3f}, p {p:.2g}, n {n}")

packet-REM-packet sextile 1 (1 lowest): mean change -0.016 Hz, cells 132; signed-rank p (b) 0.11, t p (b) 0.15; per cell-triplet n 390, mean -0.006 Hz: signed-rank p (a) 0.069, t p (a) 0.41
packet-REM-packet sextile 2 (1 lowest): mean change -0.041 Hz, cells 139; signed-rank p (b) 0.00056, t p (b) 0.0028; per cell-triplet n 616, mean -0.022 Hz: signed-rank p (a) 1.4e-05, t p (a) 0.062
packet-REM-packet sextile 3 (1 lowest): mean change -0.012 Hz, cells 137; signed-rank p (b) 0.72, t p (b) 0.57; per cell-triplet n 678, mean -0.002 Hz: signed-rank p (a) 0.27, t p (a) 0.88
packet-REM-packet sextile 4 (1 lowest): mean change +0.006 Hz, cells 140; signed-rank p (b) 0.5, t p (b) 0.78; per cell-triplet n 674, mean +0.020 Hz: signed-rank p (a) 0.18, t p (a) 0.085
packet-REM-packet sextile 5 (1 lowest): mean change +0.036 Hz, cells 128; signed-rank p (b) 0.23, t p (b) 0.33; per cell-triplet n 644, mean +0.029 Hz: signed-rank p (a) 0.26, t p (a) 0.1
packet-REM-packet sextile 6 (1 lowest): mean change +0.164 Hz, cells 113; signed-rank p (b) 0.065, t p (b) 0.037; per cell-triplet n 711, mean +0.061 Hz: signed-rank p (a) 0.071, t p (a) 0.0066
mean change of all pE cells vs REM duration: R -0.019, p 0.85, n 110 (paper R -0.22, p < 0.01)
theta (a, theta ratio) vs slowest sextile change: R 0.143, p 0.18, n 89
theta (a, theta ratio) vs fastest sextile change: R -0.005, p 0.96, n 86
theta (b, theta power) vs slowest sextile change: R 0.068, p 0.53, n 89
theta (b, theta power) vs fastest sextile change: R -0.179, p 0.1, n 86
packet-REM-packet S5 correlation, log rate change vs WAKE rate: (a) log10 WAKE R 0.068 p 0.06 n 771; (b) linear R 0.122 p 0.00069 n 771
paper: S5F R = -0.038, p = 0.24 (no rate dependence)
packet-REM-packet, M2 reading (b) sextile 1 (1 lowest): mean change +0.002 Hz, cells 164; signed-rank p (b) 0.46, t p (b) 0.85; per cell-triplet n 459, mean +0.004 Hz: signed-rank p (a) 0.3, t p (a) 0.63
packet-REM-packet, M2 reading (b) sextile 2 (1 lowest): mean change -0.012 Hz, cells 193; signed-rank p (b) 0.01, t p (b) 0.42; per cell-triplet n 609, mean -0.009 Hz: signed-rank p (a) 0.0033, t p (a) 0.31
packet-REM-packet, M2 reading (b) sextile 3 (1 lowest): mean change -0.010 Hz, cells 182; signed-rank p (b) 0.51, t p (b) 0.63; per cell-triplet n 597, mean +0.003 Hz: signed-rank p (a) 0.46, t p (a) 0.81
packet-REM-packet, M2 reading (b) sextile 4 (1 lowest): mean change -0.005 Hz, cells 189; signed-rank p (b) 0.66, t p (b) 0.85; per cell-triplet n 655, mean +0.008 Hz: signed-rank p (a) 0.63, t p (a) 0.61
packet-REM-packet, M2 reading (b) sextile 5 (1 lowest): mean change +0.022 Hz, cells 187; signed-rank p (b) 0.83, t p (b) 0.46; per cell-triplet n 663, mean +0.014 Hz: signed-rank p (a) 0.83, t p (a) 0.44
packet-REM-packet, M2 reading (b) sextile 6 (1 lowest): mean change +0.131 Hz, cells 136; signed-rank p (b) 0.16, t p (b) 0.049; per cell-triplet n 730, mean +0.064 Hz: signed-rank p (a) 0.031, t p (a) 0.0038
M2 reading (b), theta (a, theta ratio) vs slowest sextile change: R 0.083, p 0.41, n 99
M2 reading (b), theta (b, theta power) vs slowest sextile change: R 0.088, p 0.39, n 99

DP = []
for s, tab in zip(S, tabs):
    for c in fcx.cycles(s):
        pk = fcx.within(s["states"]["SWSPacket"], *c["sleep"])
        for (a1, b1), (a2, b2) in zip(pk, pk[1:]):
            r = dict(session=s["session"], dd=tab["zband"]["delta"][int(a2):int(b2)].mean() - tab["zband"]["delta"][int(a1):int(b1)].mean(),
                     # variant: only seconds whose 10 s window lies inside the packet
                     dd_in=tab["zband"]["delta"][int(a2) + 5:int(b2) - 4].mean() - tab["zband"]["delta"][int(a1) + 5:int(b1) - 4].mean())
            for k in (1, 6):
                for sfx, sx in (("", None), ("b", sx_b)):   # M2 readings (a) and (b)
                    cl = sext_cells(k, sx, a1)(s)
                    if len(cl):
                        r1 = tab["cnt"][cl, int(a1):int(b1)].sum() / (len(cl) * (b1 - a1))
                        r2 = tab["cnt"][cl, int(a2):int(b2)].sum() / (len(cl) * (b2 - a2))
                        r[f"d{k}{sfx}"] = r2 - r1
            DP.append(r)
DP = pd.DataFrame(DP)
for k in (1, 6):
    r, p, n = fcx.pearson(DP.dd, DP[f"d{k}"])
    print(f"sextile {k} (1 lowest) rate change vs delta change: R {r:.3f}, p {p:.2g}, n {n} packet pairs")
for k in (1, 6):
    r, p, n = fcx.pearson(DP.dd_in, DP[f"d{k}"])
    print(f"variant, delta from 10 s windows inside the packets only: sextile {k} R {r:.3f}, p {p:.2g}, n {n}")
for k in (1, 6):
    r, p, n = fcx.pearson(DP.dd, DP[f"d{k}b"])
    print(f"M2 reading (b), per-cycle sextiles: sextile {k} R {r:.3f}, p {p:.2g}, n {n}")

sextile 1 (1 lowest) rate change vs delta change: R 0.003, p 0.95, n 567 packet pairs
sextile 6 (1 lowest) rate change vs delta change: R -0.285, p 2.8e-11, n 524 packet pairs
variant, delta from 10 s windows inside the packets only: sextile 1 R -0.008, p 0.86, n 567
variant, delta from 10 s windows inside the packets only: sextile 6 R -0.256, p 2.7e-09, n 524
M2 reading (b), per-cycle sextiles: sextile 1 R -0.011, p 0.78, n 593
M2 reading (b), per-cycle sextiles: sextile 6 R -0.394, p 5.2e-22, n 555

for rd, df in (("a", FL), ("b", FLb)):
    i_ = df[df.type == "I"]
    ls = fcx.log_slope(i_["first"], i_["last"])
    print(f"reading ({rd}) pI log10: slope {ls['slope']:.3f} [{ls['lo']:.3f}, {ls['hi']:.3f}], includes 1: {ls['lo'] <= 1 <= ls['hi']}, "
          f"within 1.0005 +- 0.0045: {abs(ls['slope'] - 1.0005) <= 0.0045}; R {ls['r']:.3f}, n {ls['n']}, dropped {ls['dropped']}")

reading (a) pI log10: slope 0.836 [0.766, 0.906], includes 1: False, within 1.0005 +- 0.0045: False; R 0.863, n 192, dropped 9
reading (b) pI log10: slope 0.888 [0.802, 0.973], includes 1: False, within 1.0005 +- 0.0045: False; R 0.882, n 123, dropped 3

e = FLb[FLb.type == "E"]
ls = fcx.log_slope(e["first"], e["last"])
print(f"fig3d (reading b): slope {ls['slope']:.3f} [{ls['lo']:.3f}, {ls['hi']:.3f}] against 0.8535 +- 0.0265")
print(f"cv-decreases-over-sleep (reading b): R {tr_cv['b']['r']:.3f}, p {tr_cv['b']['p']:.2g} (paper R -0.05, p < 0.01)")
print(f"fig3e-wake-before-after-sextiles, parts holding of 4: {held3e}")

fig3d (reading b): slope 0.872 [0.840, 0.904] against 0.8535 +- 0.0265
cv-decreases-over-sleep (reading b): R -0.003, p 0.73 (paper R -0.05, p < 0.01)
fig3e-wake-before-after-sextiles, parts holding of 4: {'M2 a, M1 a': 2, 'M2 a, M1 b': 1, 'M2 b, M1 a': 3}

Claim	Tested
table1-dataset-counts	tested
ma-statistics	tested
packet-duration	tested
fig2a-pe-rate-distribution	tested
fig2b-wake-nrem-slope	tested
fig2a-fraction-above-2hz	tested
burst-index-by-state	tested
fig2c-autocorrelogram	tested
synchrony-skew	tested
updown-durations	tested (deposited detections)
fig1e-up-state-dynamics	tested
fig2de-fast-before-slow	tested
fig2f-rate-group-correlations	tested
fig3ab-population-rate-declines	tested
first-last-packet-mean-rate	tested
post-sleep-wake-lower	tested
fig3d-first-last-packet-slope	tested
fig3b-sextile-changes	tested
cv-decreases-over-sleep	tested
fig3c-spindles-over-sleep	tested (deposited detections)
fig3c-updown-over-sleep	tested
fig3f-within-up-changes	tested
figs3a-first-last-rem-ma-slopes	tested
fig3e-wake-before-after-sextiles	tested
fig4-lfp-within-nrem	tested
fig4-updown-within-packets	tested
fig4-spindles-within-nrem	tested
fig4e-late-packet-sigma	tested
fig4-sextiles-within-nrem	tested
fig4cd-packet-rate-lfp-correlations	tested
fig4e-within-up-sextiles	tested
figs5-within-wake-rem	tested
fig5b-ma-nrem-ma	tested
figs5d-rem-nrem-rem	tested
fig5c-sigma-ma-change	tested
fig5b-wake-between-packets	tested
fig5b-ma-between-packets	tested
fig5b-rem-between-packets	tested
fig5c-delta-change-slow-cells	tested
figs2d-interneurons-unchanged	tested
sleep-homogenizes-rates	tested

Claim	Outcome
table1-dataset-counts	partially reproduced
ma-statistics	partially reproduced
packet-duration	reproduced
fig2a-pe-rate-distribution	not reproduced
fig2b-wake-nrem-slope	reproduced
fig2a-fraction-above-2hz	reproduced
burst-index-by-state	partially reproduced
fig2c-autocorrelogram	partially reproduced
synchrony-skew	reproduced
updown-durations	reproduced
fig1e-up-state-dynamics	not reproduced
fig2de-fast-before-slow	reproduced
fig2f-rate-group-correlations	not reproduced
fig3ab-population-rate-declines	reproduced
first-last-packet-mean-rate	reproduced
post-sleep-wake-lower	partially reproduced
fig3d-first-last-packet-slope	reproduced
fig3b-sextile-changes	reproduced
cv-decreases-over-sleep	not reproduced
fig3c-spindles-over-sleep	partially reproduced (assumed)
fig3c-updown-over-sleep	not reproduced
fig3f-within-up-changes	reproduced
figs3a-first-last-rem-ma-slopes	partially reproduced
fig3e-wake-before-after-sextiles	partially reproduced
fig4-lfp-within-nrem	partially reproduced
fig4-updown-within-packets	reproduced
fig4-spindles-within-nrem	reproduced
fig4e-late-packet-sigma	reproduced
fig4-sextiles-within-nrem	reproduced (assumed)
fig4cd-packet-rate-lfp-correlations	partially reproduced
fig4e-within-up-sextiles	not reproduced
figs5-within-wake-rem	not reproduced
fig5b-ma-nrem-ma	not reproduced (assumed)
figs5d-rem-nrem-rem	partially reproduced (assumed)
fig5c-sigma-ma-change	not reproduced (assumed)
fig5b-wake-between-packets	partially reproduced (assumed)
fig5b-ma-between-packets	partially reproduced (assumed)
fig5b-rem-between-packets	not reproduced (assumed)
fig5c-delta-change-slow-cells	not reproduced
figs2d-interneurons-unchanged	not reproduced
sleep-homogenizes-rates	partially reproduced

Reproducing Watson et al. (2016) from CRCNS fcx-1¶

Data check¶

table1-dataset-counts¶

ma-statistics¶

packet-duration¶

fig2a-pe-rate-distribution¶

fig2b-wake-nrem-slope¶

fig2a-fraction-above-2hz¶

burst-index-by-state¶

fig2c-autocorrelogram¶

synchrony-skew¶

updown-durations¶

fig1e-up-state-dynamics¶

fig2de-fast-before-slow¶

fig2f-rate-group-correlations¶

fig3ab-population-rate-declines¶

first-last-packet-mean-rate¶

post-sleep-wake-lower¶

fig3d-first-last-packet-slope¶

fig3b-sextile-changes¶

cv-decreases-over-sleep¶

fig3c-spindles-over-sleep¶

fig3c-updown-over-sleep¶

fig3f-within-up-changes¶

figs3a-first-last-rem-ma-slopes¶

fig3e-wake-before-after-sextiles¶

fig4-lfp-within-nrem¶

fig4-updown-within-packets¶

fig4-spindles-within-nrem¶

fig4e-late-packet-sigma¶

fig4-sextiles-within-nrem¶

fig4cd-packet-rate-lfp-correlations¶

fig4e-within-up-sextiles¶

figs5-within-wake-rem¶

Triplets (Figure 5 and S5)¶

fig5b-ma-nrem-ma¶

figs5d-rem-nrem-rem¶

fig5c-sigma-ma-change¶

fig5b-wake-between-packets¶

fig5b-ma-between-packets¶

fig5b-rem-between-packets¶

fig5c-delta-change-slow-cells¶

figs2d-interneurons-unchanged¶

sleep-homogenizes-rates¶

Summary¶