initial

.flake8

@@ -0,0 +1,11 @@
+[flake8]
+extend-ignore = E124,E128,E301,E302,E305,E402,E501,E261,W504
+# E124: closing bracket does not match visual indentation
+# E128: continuation line under-indented for visual indent
+# E301: expected 1 blank line, found 0
+# E302: expected 2 blank lines, found 1
+# E305: expected 2 blank lines after class or function definition, found 1
+# E402: module level import not at top of file
+# E501: line too long (82 > 79 characters)
+# E261: at least two spaces before inline comment
+# W504: line break after binary operator

.github/workflows/gcr-spikeforestxyz.yml

@@ -0,0 +1,36 @@
+name: Build spikeforestxyz image and push to GCR
+
+on:
+  push:
+    branches:
+      - main
+    paths:
+      - "dendro_apps/spikeforestxyz/**"
+  workflow_dispatch:
+
+jobs:
+  publish-docker-image:
+    runs-on: ubuntu-latest
+
+    permissions:
+      contents: read
+      packages: write
+
+    steps:
+      - name: Checkout code
+        uses: actions/checkout@v3
+
+      - name: Login to GitHub Container Registry
+        uses: docker/login-action@v1
+        with:
+          registry: ghcr.io
+          username: ${{ github.actor }}
+          password: ${{ secrets.GITHUB_TOKEN }}
+
+      - name: Build the Docker image
+        run: |
+          cd dendro_apps/spikeforestxyz && \
+          NAME="spikeforestxyz" && \
+          docker buildx build --push \
+          -t ghcr.io/magland/$NAME:latest \
+          -f Dockerfile .

.gitignore

@@ -0,0 +1,34 @@
+tmp
+
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+dist/
+build/
+*.egg-info/
+
+# Virtual environments
+venv/
+env/
+.env/
+
+# IDE specific files
+.idea/
+.vscode/
+
+# Compiled Python files
+*.pyc
+
+# Logs and databases
+*.log
+*.sqlite3
+
+# OS generated files
+.DS_Store
+Thumbs.db

README.md

@@ -0,0 +1,4 @@
+# SpikeForestXYZ
+
+https://dendro.vercel.app/project/9e302504?tab=project-home
+

dendro_apps/make_spec.sh

@@ -0,0 +1,5 @@
+#!/bin/bash
+
+set -ex
+
+dendro make-app-spec-file --app-dir spikeforestxyz --spec-output-file spikeforestxyz/spec.json

dendro_apps/spikeforestxyz/Dockerfile

@@ -0,0 +1,22 @@
+FROM python:3.9-slim
+
+# Install dendro
+RUN pip install dendro==0.2.15
+
+# Install kachery-cloud
+RUN pip install kachery-cloud==0.4.9
+
+# Install spikeinterface
+RUN pip install spikeinterface==0.100.6
+
+# Install lindi
+RUN pip install lindi==0.3.4
+
+# Install sortingview
+RUN pip install sortingview==0.13.3
+
+# Copy files into the container
+RUN mkdir /app
+COPY *.py /app/
+COPY recording_summary/*.py /app/recording_summary/
+COPY recording_summary/helpers/*.py /app/recording_summary/helpers/

dendro_apps/spikeforestxyz/main.py

@@ -0,0 +1,18 @@
+#!/usr/bin/env python
+
+
+from dendro.sdk import App
+from recording_summary.recording_summary import RecordingSummaryProcessor
+
+app = App(
+    name="spikeforestxyz",
+    description="Processors for SpikeForestXYZ",
+    app_image="ghcr.io/magland/spikeforestxyz:latest",
+    app_executable="/app/main.py",
+)
+
+
+app.add_processor(RecordingSummaryProcessor)
+
+if __name__ == "__main__":
+    app.run()

dendro_apps/spikeforestxyz/recording_summary/create_recording_summary.py

@@ -0,0 +1,111 @@
+import time
+import uuid
+import shutil
+import numpy as np
+import h5py
+import lindi
+import kachery_cloud as kcl
+from .helpers.compute_correlogram_data import compute_correlogram_data
+# from .nwbextractors import NwbRecordingExtractor, NwbSortingExtractor
+
+
+def create_recording_summary(
+    nwb_lindi_fname: str
+):
+    staging_area = lindi.StagingArea.create(dir=nwb_lindi_fname + '.d')
+    f = lindi.LindiH5pyFile.from_lindi_file(
+        nwb_lindi_fname,
+        mode='r+',
+        staging_area=staging_area,
+        local_cache=lindi.LocalCache()
+    )
+    # rec = NwbRecordingExtractor(h5py_file=f)  # type: ignore
+    # sorting_true = NwbSortingExtractor(h5py_file=f)  # type: ignore
+
+    # Load the spike times from the units group
+    units_group = f['/units']
+    assert isinstance(units_group, h5py.Group)
+    print('Loading spike times')
+    spike_times = units_group['spike_times'][()]  # type: ignore
+    spike_times_index = units_group['spike_times_index'][()]  # type: ignore
+    num_units = len(spike_times_index)
+    total_num_spikes = len(spike_times)
+    print(f'Loaded {num_units} units with {total_num_spikes} total spikes')
+
+    # Compute autocorrelograms for all the units
+    print('Computing autocorrelograms')
+    auto_correlograms = []
+    p = 0
+    timer = time.time()
+    for i in range(num_units):
+        spike_train = spike_times[p:spike_times_index[i]]
+        elapsed = time.time() - timer
+        if elapsed > 2:
+            print(f'Computing autocorrelogram for unit {i + 1} of {num_units} ({len(spike_train)} spikes)')
+            timer = time.time()
+        r = compute_correlogram_data(
+            spike_train_1=spike_train,
+            spike_train_2=None,
+            window_size_msec=100,
+            bin_size_msec=1
+        )
+        bin_edges_sec = r['bin_edges_sec']
+        bin_counts = r['bin_counts']
+        auto_correlograms.append({
+            'bin_edges_sec': bin_edges_sec,
+            'bin_counts': bin_counts
+        })
+        p = spike_times_index[i]
+    autocorrelograms_array = np.zeros(
+        (num_units, len(auto_correlograms[0]['bin_counts'])),
+        dtype=np.uint32
+    )
+    for i, ac in enumerate(auto_correlograms):
+        autocorrelograms_array[i, :] = ac['bin_counts']
+    bin_edges_array = np.zeros(
+        (num_units, len(auto_correlograms[0]['bin_edges_sec'])),
+        dtype=np.float32
+    )
+    for i, ac in enumerate(auto_correlograms):
+        bin_edges_array[i, :] = ac['bin_edges_sec']
+
+    # Create a new dataset in the units group to store the autocorrelograms
+    print('Writing autocorrelograms to output file')
+    ds = units_group.create_dataset('acg', data=autocorrelograms_array)
+    ds.attrs['description'] = 'the autocorrelogram for each spike unit'
+    ds.attrs['namespace'] = 'hdmf-common'
+    ds.attrs['neurodata_type'] = 'VectorData'
+    ds.attrs['object_id'] = str(uuid.uuid4())
+
+    ds = units_group.create_dataset('acg_bin_edges', data=bin_edges_array)
+    ds.attrs['description'] = 'the bin edges in seconds for the autocorrelogram for each spike unit'
+    ds.attrs['namespace'] = 'hdmf-common'
+    ds.attrs['neurodata_type'] = 'VectorData'
+    ds.attrs['object_id'] = str(uuid.uuid4())
+
+    # Update the colnames attribute of the units group
+    colnames = units_group.attrs['colnames']
+    assert isinstance(colnames, np.ndarray)
+    colnames = colnames.tolist()
+    colnames.append('acg')
+    colnames.append('acg_bin_edges')
+    units_group.attrs['colnames'] = colnames
+
+    f.flush()  # write changes to the file
+
+    def on_store_blob(filename: str):
+        url = kcl.store_file(filename)
+        return url
+
+    def on_store_main(filename: str):
+        shutil.copyfile(filename, nwb_lindi_fname)
+        return nwb_lindi_fname
+
+    staging_store = f.staging_store
+    assert staging_store is not None
+    print('Uploading supporting files')
+    f.upload(
+        on_upload_blob=on_store_blob,
+        on_upload_main=on_store_main
+    )
+

dendro_apps/spikeforestxyz/recording_summary/helpers/compute_correlogram_data.py

@@ -0,0 +1,96 @@
+from typing import Union
+import numpy as np
+
+
+def compute_correlogram_data(
+    *,
+    spike_train_1: np.ndarray,
+    spike_train_2: Union[np.ndarray, None] = None,
+    window_size_msec: float = 100,
+    bin_size_msec: float = 1
+):
+    times1 = spike_train_1
+    num_bins = int(window_size_msec / bin_size_msec)
+    if num_bins % 2 == 0:
+        num_bins = num_bins - 1  # odd number of bins
+    num_bins_half = int((num_bins + 1) / 2)
+    bin_edges_msec = np.array(
+        (np.arange(num_bins + 1) - num_bins / 2) * bin_size_msec, dtype=np.float32
+    )
+    bin_counts = np.zeros((num_bins,), dtype=np.int32)
+    if spike_train_2 is None:
+        # autocorrelogram
+        offset = 1
+        while True:
+            if offset >= len(times1):
+                break
+            deltas_msec = (times1[offset:] - times1[:-offset]) * 1000
+            deltas_msec = deltas_msec[deltas_msec <= bin_edges_msec[-1]]
+            if len(deltas_msec) == 0:
+                break
+            for i in range(num_bins_half):
+                start_msec = bin_edges_msec[num_bins_half - 1 + i]
+                end_msec = bin_edges_msec[num_bins_half + i]
+                ct = len(
+                    deltas_msec[(start_msec <= deltas_msec) & (deltas_msec < end_msec)]
+                )
+                bin_counts[num_bins_half - 1 + i] += ct
+                bin_counts[num_bins_half - 1 - i] += ct
+            offset = offset + 1
+    else:
+        # cross-correlogram
+        times2 = spike_train_2
+        all_times = np.concatenate((times1, times2))
+        all_labels = np.concatenate(
+            (1 * np.ones(times1.shape), 2 * np.ones(times2.shape))
+        )
+        sort_inds = np.argsort(all_times)
+        all_times = all_times[sort_inds]
+        all_labels = all_labels[sort_inds]
+        offset = 1
+        while True:
+            if offset >= len(all_times):
+                break
+            deltas_msec = (all_times[offset:] - all_times[:-offset]) * 1000
+
+            deltas12_msec = deltas_msec[
+                (all_labels[offset:] == 2) & (all_labels[:-offset] == 1)
+            ]
+            deltas21_msec = deltas_msec[
+                (all_labels[offset:] == 1) & (all_labels[:-offset] == 2)
+            ]
+            deltas11_msec = deltas_msec[
+                (all_labels[offset:] == 1) & (all_labels[:-offset] == 1)
+            ]
+            deltas22_msec = deltas_msec[
+                (all_labels[offset:] == 2) & (all_labels[:-offset] == 2)
+            ]
+
+            deltas12_msec = deltas12_msec[deltas12_msec <= bin_edges_msec[-1]]
+            deltas21_msec = deltas21_msec[deltas21_msec <= bin_edges_msec[-1]]
+            deltas11_msec = deltas11_msec[deltas11_msec <= bin_edges_msec[-1]]
+            deltas22_msec = deltas22_msec[deltas22_msec <= bin_edges_msec[-1]]
+
+            if len(deltas12_msec) + len(deltas21_msec) + len(deltas11_msec) + len(deltas22_msec) == 0:
+                break
+
+            for i in range(num_bins_half):
+                start_msec = bin_edges_msec[num_bins_half - 1 + i]
+                end_msec = bin_edges_msec[num_bins_half + i]
+                ct12 = len(
+                    deltas12_msec[
+                        (start_msec <= deltas12_msec) & (deltas12_msec < end_msec)
+                    ]
+                )
+                ct21 = len(
+                    deltas21_msec[
+                        (start_msec <= deltas21_msec) & (deltas21_msec < end_msec)
+                    ]
+                )
+                bin_counts[num_bins_half - 1 + i] += ct12
+                bin_counts[num_bins_half - 1 - i] += ct21
+            offset = offset + 1
+    return {
+        "bin_edges_sec": (bin_edges_msec / 1000).astype(np.float32),
+        "bin_counts": bin_counts.astype(np.int32),
+    }