Demo 04 — parameter sweeps

How the detection parameters were chosen, using demo 03's hippocampus job.

Two 3x3 sweeps:

• corrmap — cell_filt_xy_um x intensity_thresh. The lateral scale Suite3D looks for, and how much of the movie survives the per-voxel noise gate. Together these largely set the diameter of what gets detected.
• segmentation — segmentation_spatial_filt x vox_snr_thresh. Frame smoothing before peak detection, and the fraction of a voxel's variance the ROI's trace must explain for that voxel to join it.

Run ../03-hippocampus/run_pipeline.py first.

%matplotlib inline
# inline PNGs, never ipywidgets -- these pages go on the public web.

import sys, os, shutil
sys.path.insert(0, os.path.abspath('..'))

import numpy as np
import matplotlib.pyplot as plt

from suite3d.job import Job
from suite3d.curation import collate_sweep_results

Job.setup_sweep requires the job's current value of each swept parameter to appear in that parameter's list, so each list brackets the default.

JOB_ROOT = os.environ.get('SUITE3D_DEMO_OUT', './results')  # demo 03's --out-dir
JOB_ID   = 'demo-hippocampus'

job = Job(JOB_ROOT, JOB_ID, create=False)
job_dir = job.dirs['job_dir']

CORRMAP_SWEEP = {'cell_filt_xy_um': [0.5, 1.5, 2.5],   # default 1.5
                 'intensity_thresh': [1, 3, 5]}        # default 3
SEG_SWEEP     = {'segmentation_spatial_filt': [0, 2, 4],   # default 2
                 'vox_snr_thresh': [0.05, 0.10, 0.20]}     # default 0.10

The trap

Job.save_params() rewrites the job's root params.npy on every mutation, and a sweep mutates it many times. Snapshot it first, and restore it in a finally: block, or an interrupted sweep leaves your job pinned to the last cell's parameters.

root   = os.path.join(job_dir, 'params.npy')
backup = os.path.join(job_dir, 'params.npy.sweep_backup')
if not os.path.exists(backup):
    shutil.copy2(root, backup)

try:
    job.sweep_corrmap(CORRMAP_SWEEP, sweep_name='corrmap')
finally:
    shutil.copy2(backup, root)
    print('restored root params.npy')

Looking at the corrmap sweep

sweep_dir = os.path.join(job_dir, 'sweeps', 'corrmap')
summary = np.load(os.path.join(sweep_dir, 'sweep_summary.npy'),
                  allow_pickle=True).item()

# sweep_corrmap() stores each run's correlation map under the key 'corrmap'
# (sweep_segmentation() uses 'stats'). `results` comes back as an object array
# with one axis per swept parameter, in the order they appear in the sweep dict.
results, sweep_params = collate_sweep_results(summary, result_key='corrmap')
print(results.shape, sweep_params.shape)

# One panel per (cell_filt_xy_um, intensity_thresh) cell, at a single plane.
Z = 2
vals_a = CORRMAP_SWEEP['cell_filt_xy_um']
vals_b = CORRMAP_SWEEP['intensity_thresh']
fig, axs = plt.subplots(len(vals_a), len(vals_b),
                        figsize=(3.2 * len(vals_b), 3.2 * len(vals_a)))
for i, a in enumerate(vals_a):
    for j, b in enumerate(vals_b):
        v = np.asarray(results[i, j])
        img = v[Z] if v.ndim == 3 else v
        axs[i, j].imshow(img, cmap='magma', vmax=np.percentile(img, 99.5))
        axs[i, j].set_title(f'cell_filt={a}, int_thresh={b}', fontsize=8)
        axs[i, j].axis('off')
fig.suptitle(f'Correlation map, z={Z}')
plt.tight_layout(); plt.show()

Smaller cell_filt_xy_um sharpens small structures and amplifies noise; larger blurs neighbouring somata together. intensity_thresh trades sensitivity to dim cells against background.

The segmentation sweep

This one re-runs detection for each of its 9 cells, so it is much more expensive than the corrmap sweep. Same snapshot/restore discipline.

try:
    job.sweep_segmentation(SEG_SWEEP, sweep_name='seg', all_combinations=True)
finally:
    shutil.copy2(backup, root)
    print('restored root params.npy')

seg_dir = os.path.join(job_dir, 'sweeps', 'seg')
seg_summary = np.load(os.path.join(seg_dir, 'sweep_summary.npy'),
                      allow_pickle=True).item()
seg_results, seg_params = collate_sweep_results(seg_summary, result_key='stats')

counts = np.vectorize(lambda s: len(s), otypes=[int])(seg_results)
print('ROIs detected per (segmentation_spatial_filt, vox_snr_thresh):')
print(counts)

Opening the sweep in napari

The sweep viewer lays the grid out interactively -- scrub planes, compare cells, and see which ROIs appear and vanish.

# From a terminal (recommended -- napari wants its own event loop):
#   python open_sweep_napari.py --sweep-dir ./results/s3d-demo-hippocampus/sweeps/corrmap
#
# Or inline:
#   import napari
#   from suite3d.curation import SweepUI
#   ui = SweepUI(sweep_dir)
#   ui.load_outputs()
#   ui.create_ui()
#   napari.run()