Encoding model with all three cameras: Davatolhagh et al. (2025), dandiset 001712¶
A variant of the encoding model of reproduce_davatolhagh_2025.ipynb, for the claims
fig4c-full-cvr2 and fig4cd-movement-dominates. The paper's movement predictors include the video
components of the cameras; the claims file fixed a design with the body camera's components only,
for runtime, and that design decides the outcomes in the main notebook. Here the top 200 components
of the video and of the video motion energy of the left and right cameras are added, as in the
paper. Everything else follows the main notebook and wfprior.py. Streaming and decoding the three
videos takes most of the run, about 15 minutes.
Paper: doi:10.1101/2025.09.12.675910. Dandiset:
DANDI:001712, version draft.
In [1]:
import time
import warnings
import numpy as np
import pandas as pd
import wfprior as wf
pd.set_option("display.width", 200)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
info = wf.draft_info()
print(f"dandiset {wf.DANDISET_ID} version {wf.DANDISET_VERSION}, modified {info['modified']}")
diff = wf.check_assets(info)
print("Assets differ from those this entry was written against:" if diff else
"The assets read match those this entry was written against.", *diff, sep="\n")
f = wf.open_nwb(wf.PROCESSED)
trials = wf.load_trials(f)
svd = wf.load_svd(f)
frame_t = svd["t"]
movement = wf.load_movement(f, frame_t)
vids = {}
for cam in ("Body", "Left", "Right"):
t0 = time.time()
vids[cam] = wf.video_components(f, frame_t, cam)
print(f"{cam} camera decoded in {time.time() - t0:.0f} s")
dandiset 001712 version draft, modified 2026-05-28 10:04:28.686333+00:00 The assets read match those this entry was written against.
Body camera decoded in 84 s
Left camera decoded in 423 s
Right camera decoded in 288 s
Data check¶
In [2]:
print(f"Trials: {len(trials)}; choice trials {int((trials.mouse_wheel_choice != 'none').sum())}; imaging frames {len(frame_t)}; "
f"temporal components {svd['svt'].shape[1]}")
for cam, v in vids.items():
print(f"{cam} camera: {v['n_video_frames']} decoded frames, {v['n_timestamps']} timestamps, {v['pixels']} pixels after binning "
f"({wf.VIDEO_BIN[cam]} x {wf.VIDEO_BIN[cam]} blocks); top {wf.N_VIDEO_COMPONENTS} components explain {v['video_explained']:.3f} "
f"of video and {v['motion_explained']:.3f} of motion-energy variance; component range {v['video'].min():.0f} to {v['video'].max():.0f}")
print(movement.describe().T[["mean", "min", "max"]].round(3).to_string())
Trials: 913; choice trials 911; imaging frames 124886; temporal components 200
Body camera: 119964 decoded frames, 119964 timestamps, 1280 pixels after binning (16 x 16 blocks); top 200 components explain 0.996 of video and 0.969 of motion-energy variance; component range -1056 to 833
Left camera: 239237 decoded frames, 239237 timestamps, 1280 pixels after binning (32 x 32 blocks); top 200 components explain 0.979 of video and 0.956 of motion-energy variance; component range -540 to 665
Right camera: 600397 decoded frames, 600397 timestamps, 1280 pixels after binning (16 x 16 blocks); top 200 components explain 0.989 of video and 0.957 of motion-energy variance; component range -378 to 561
mean min max
wheel_velocity -0.182 -7.026 6.478
paw_left_leftcam 343.419 0.183 8266.333
paw_right_leftcam 472.724 0.156 10671.467
paw_left_rightcam 298.151 0.684 4750.300
paw_right_rightcam 224.384 0.005 5833.166
motion_energy_left 7.642 1.505 117.863
motion_energy_right 3.500 0.764 47.627
motion_energy_body 2.420 0.348 67.980
Every camera has as many decoded frames as timestamps, so frames and times pair one to one. The movement variables are those of the main notebook.
fig4c-full-cvr2 and fig4cd-movement-dominates, with all three cameras¶
In [3]:
t0 = time.time()
ok = (trials.mouse_wheel_choice != "none").to_numpy()
ak = wf.fit_action_kernel(trials)
rows_frame, rows_tp, rows_trial = wf.segment_frames(trials, frame_t, ok)
points = wf.STIM_WINDOW[1] - wf.STIM_WINDOW[0] + wf.FEEDBACK_WINDOW[1] - wf.FEEDBACK_WINDOW[0]
mov = np.hstack([movement.to_numpy()] + [vids[c][k] for c in ("Body", "Left", "Right") for k in ("video", "motion")])
X, index, analog = wf.build_design(trials, rows_trial, rows_tp, points, ak["prior"], mov, rows_frame)
Y = svd["svt"][rows_frame]
task_cols = np.concatenate([index[k] for k in ("choice", "stim_left", "stim_right", "outcome", "true_prior", "subjective_prior")])
mov_cols = index["movement_columns"]
icp = index["intercept"]
u_flat = svd["u"][:, ::2, ::2].reshape(len(svd["u"]), -1).T.astype(np.float64)
print(f"Rows {X.shape[0]}, columns {X.shape[1]}: intercept {len(icp)}, task {len(task_cols)}, movement {len(mov_cols)}")
test_folds = wf.folds(rows_tp, np.random.default_rng(404)) # the folds of the main notebook
fits = {name: wf.fit_model(X, Y, analog, test_folds, cols, np.random.default_rng(404 + i), u_flat)
for i, (name, cols) in enumerate((("full", []), ("task shuffled", task_cols), ("movement shuffled", mov_cols)))}
fits["intercept only"] = wf.fit_model(X[:, icp], Y, [], test_folds, [], np.random.default_rng(403), u_flat)
# Variant as in the main notebook: each frame once per trial, and folds that keep each trial's rows together.
_, first = np.unique(np.stack([rows_trial, rows_frame], 1), axis=0, return_index=True)
keep_rows = np.sort(first)
trial_test = wf.trial_folds(rows_trial[keep_rows], np.random.default_rng(430))
Xg, Yg = X[keep_rows], Y[keep_rows]
for i, (name, cols) in enumerate((("full", []), ("task shuffled", task_cols), ("movement shuffled", mov_cols))):
fits[f"trial folds: {name}"] = wf.fit_model(Xg, Yg, analog, trial_test, cols, np.random.default_rng(431 + i), u_flat)
fits["trial folds: intercept only"] = wf.fit_model(Xg[:, icp], Yg, [], trial_test, [], np.random.default_rng(434), u_flat)
readings = {"a pooled pixel space": "pooled", "b mean over pixels": "pixel_mean", "c mean over components": "component_mean"}
table = pd.DataFrame({name: {r: fit[k] for r, k in readings.items()} for name, fit in fits.items()}).T
print(f"cvR2 by model and reading, mean over {wf.N_FOLDS} folds ({time.time() - t0:.0f} s):")
print(table.round(4).to_string())
lo, hi = 0.573 - 0.087, 0.573 + 0.087
mlo, mhi = 0.3685 - 0.0715, 0.3685 + 0.0715
rows = []
for r, k in readings.items():
full, icpt = table.loc["full", r], table.loc["intercept only", r]
dm, dt_ = full - table.loc["movement shuffled", r], full - table.loc["task shuffled", r]
fm = fits["full"]["folds"][k] - fits["movement shuffled"]["folds"][k]
ft = fits["full"]["folds"][k] - fits["task shuffled"]["folds"][k]
rows.append({"reading": r, "full": full, "full minus intercept": full - icpt, "full in 0.486-0.660": lo <= full <= hi,
"movement dR2": dm, "task dR2": dt_, "movement in 0.297-0.440": mlo <= dm <= mhi,
"movement folds min": fm.min(), "task folds max": ft.max()})
print(pd.DataFrame(rows).set_index("reading").round(4).to_string())
print("Variant, trial folds without repeated frames:")
for r in readings:
full = table.loc["trial folds: full", r]
print(f" {r}: full {full:.4f} (minus intercept {full - table.loc['trial folds: intercept only', r]:.4f}), "
f"movement dR2 {full - table.loc['trial folds: movement shuffled', r]:.4f}, task dR2 {full - table.loc['trial folds: task shuffled', r]:.4f}")
al = fits["full"]["alphas"]
print(f"Full model penalties: median {np.median(al):g}, at grid minimum {(al == wf.ALPHAS[0]).mean():.3f}, at grid maximum {(al == wf.ALPHAS[-1]).mean():.3f}")
Rows 63770, columns 1698: intercept 70, task 420, movement 1208
cvR2 by model and reading, mean over 10 folds (229 s):
a pooled pixel space b mean over pixels c mean over components
full 0.5740 0.5251 0.2075
task shuffled 0.5660 0.5172 0.2003
movement shuffled 0.3580 0.2929 0.0466
intercept only 0.3481 0.2886 0.0486
trial folds: full 0.4722 0.4050 0.0461
trial folds: task shuffled 0.4640 0.3976 0.0457
trial folds: movement shuffled 0.3630 0.2905 0.0383
trial folds: intercept only 0.3541 0.2888 0.0447
full full minus intercept full in 0.486-0.660 movement dR2 task dR2 movement in 0.297-0.440 movement folds min task folds max
reading
a pooled pixel space 0.5740 0.2259 True 0.2160 0.0080 False 0.2089 0.0094
b mean over pixels 0.5251 0.2364 True 0.2322 0.0079 False 0.2241 0.0092
c mean over components 0.2075 0.1589 False 0.1608 0.0072 False 0.1586 0.0076
Variant, trial folds without repeated frames:
a pooled pixel space: full 0.4722 (minus intercept 0.1181), movement dR2 0.1092, task dR2 0.0082
b mean over pixels: full 0.4050 (minus intercept 0.1163), movement dR2 0.1145, task dR2 0.0075
c mean over components: full 0.0461 (minus intercept 0.0014), movement dR2 0.0078, task dR2 0.0004
Full model penalties: median 56.2341, at grid minimum 0.000, at grid maximum 0.000
Result: variant; it does not decide an outcome. With the video components of all three cameras, the full model's cvR2 is 0.574 pooled in pixel space and 0.525 as the mean over pixels, both within 0.486 to 0.660, and 0.208 as the mean over components. Shuffling the movement predictors lowers it by 0.216, 0.232 and 0.161, below 0.297 to 0.440 under every reading, and shuffling the task predictors by 0.008, 0.008 and 0.007. With folds that keep each trial's frames together and each frame once per trial, the full model's cvR2 is 0.472, 0.405 and 0.046, below the range under every reading, and the movement drop is 0.109, 0.115 and 0.008 against 0.008, 0.008 and 0.000 for task. The side cameras bring the full model into the range only under the chiCa folds, which share frames of a trial between training and test rows.
In [4]:
print(f"Notebook ran in {(time.time() - started) / 60:.1f} min")
Notebook ran in 17.4 min