From a0704fc835e35128e4832bb5ff939caee767e77b Mon Sep 17 00:00:00 2001 From: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> Date: Tue, 31 Aug 2021 21:31:39 -0400 Subject: [PATCH 01/12] Add save option to make_dataset Signed-off-by: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> --- datagen.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/datagen.py b/datagen.py index 4e10736..a623045 100644 --- a/datagen.py +++ b/datagen.py @@ -89,7 +89,7 @@ def gen_data_diffeq(f: Callable, projection: Callable, *, t, x0: np.ndarray, dim return ivp["t"], y, projed -def make_dataset(f, x0, num_trajectories, num_dim, begin, end, noise): +def make_dataset(f, x0, num_trajectories, num_dim, begin, end, noise, save=True): xx = [] # states projeds = [] # observations rng = np.random.RandomState(39) @@ -110,8 +110,9 @@ def make_dataset(f, x0, num_trajectories, num_dim, begin, end, noise): ys = projeds us = np.zeros((xx.shape[0], xx.shape[1], 1)) - filename = f"{f.__name__}_{num_trajectories}trajectories_{num_dim}dim_{begin}to{end}_noise{noise}.npz" - np.savez(filename, x = xs, y = ys, u = us) + if save: + filename = f"{f.__name__}_{num_trajectories}trajectories_{num_dim}dim_{begin}to{end}_noise{noise}.npz" + np.savez(filename, x = xs, y = ys, u = us) def generate_lorenz(): From 03a27755cc63c8d6a6af342143c755dae46c3fb7 Mon Sep 17 00:00:00 2001 From: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> Date: Tue, 31 Aug 2021 21:32:28 -0400 Subject: [PATCH 02/12] Add test_make_dataset Signed-off-by: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> --- datagen.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/datagen.py b/datagen.py index a623045..53558f2 100644 --- a/datagen.py +++ b/datagen.py @@ -145,6 +145,19 @@ def generate_vdp(): make_dataset(vdp, x0=np.array([0.1, 0.1]), num_trajectories=num_trajectory, num_dim=num_dim, begin=begin, end=end, noise=noise) +def test_make_dataset(): + params = dict( + num_trajectories=1, + begin=500, + end=20500, + noise=0.05, + save=True + ) + + make_dataset(lorenz, num_dim=3, x0=np.array([0., 1., 1.05]), **params) + make_dataset(vdp, num_dim=2, x0=np.array([0.1, 0.1]) , **params) + + if __name__ == '__main__': if len(sys.argv) == 2 and sys.argv[1] == 'lorenz': generate_lorenz() From bbd388292b33ffb002bdada29c3bfb60c8c6b5f7 Mon Sep 17 00:00:00 2001 From: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> Date: Tue, 31 Aug 2021 22:08:59 -0400 Subject: [PATCH 03/12] Add CPU environment Signed-off-by: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> --- environment-cpu.yml | 15 +++++++++++++++ 1 file changed, 15 insertions(+) create mode 100644 environment-cpu.yml diff --git a/environment-cpu.yml b/environment-cpu.yml new file mode 100644 index 0000000..1ba67aa --- /dev/null +++ b/environment-cpu.yml @@ -0,0 +1,15 @@ +name: neuralflow +channels: + - conda-forge +dependencies: + - python=3.9 + - jax + - jupyter + - jupyterlab + - matplotlib + - numpy + - pip + - pytest + - scipy + - scikit-learn + - tqdm From 1b83577fbc1a2893a1444385d138562626f2361b Mon Sep 17 00:00:00 2001 From: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> Date: Tue, 31 Aug 2021 22:12:32 -0400 Subject: [PATCH 04/12] Change mpl backend to default Signed-off-by: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> --- scripts/run_bubblewrap.py | 1 - 1 file changed, 1 deletion(-) diff --git a/scripts/run_bubblewrap.py b/scripts/run_bubblewrap.py index 549e06a..deca800 100644 --- a/scripts/run_bubblewrap.py +++ b/scripts/run_bubblewrap.py @@ -6,7 +6,6 @@ import numpy as np import matplotlib -matplotlib.use('TkAgg') import matplotlib.pylab as plt from matplotlib.patches import Ellipse from mpl_toolkits.mplot3d import Axes3D From e9096e5d333ff4ee3185ca30815f71e0157f651c Mon Sep 17 00:00:00 2001 From: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> Date: Tue, 31 Aug 2021 22:13:19 -0400 Subject: [PATCH 05/12] Cut run by 10-fold Signed-off-by: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> --- scripts/run_bubblewrap.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/run_bubblewrap.py b/scripts/run_bubblewrap.py index deca800..9dd21d6 100644 --- a/scripts/run_bubblewrap.py +++ b/scripts/run_bubblewrap.py @@ -20,7 +20,7 @@ s = np.load('lorenz_1trajectories_3dim_500to20500_noise0.05.npz') data = s['y'][0] -T = data.shape[0] # should be 20k +T = data.shape[0] // 10 # should be 20k d = data.shape[1] # should be 2 ## Parameters From 0989b69052db36bc0827d58b1e922023536d877e Mon Sep 17 00:00:00 2001 From: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> Date: Tue, 31 Aug 2021 22:17:37 -0400 Subject: [PATCH 06/12] Add workflow. Signed-off-by: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> --- .github/workflows/python.yaml | 45 +++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) create mode 100644 .github/workflows/python.yaml diff --git a/.github/workflows/python.yaml b/.github/workflows/python.yaml new file mode 100644 index 0000000..93ace6d --- /dev/null +++ b/.github/workflows/python.yaml @@ -0,0 +1,45 @@ +name: Run with Conda + +on: + push: + branches: + - main + pull_request: + schedule: + - cron: "2 5 * * *" + +jobs: + test: + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v2 + + - name: Set up Python 3.9 + uses: actions/setup-python@v2 + with: + python-version: 3.9 + + - name: Add conda to system path + run: echo $CONDA/bin >> $GITHUB_PATH + + - name: Install dependencies + run: | + conda config --add channels conda-forge + conda env update --file environment-cpu.yml --name base + + - name: Run Bubblewrap + run: | + conda install jupytext + pytest datagen.py + jupytext --from py --to ipynb scripts/run_bubblewrap.py --execute --run-path . + jupyter nbconvert --to html scripts/run_bubblewrap.ipynb + mkdir public + mv scripts/run_bubblewrap.html public/index.html + + - name: Deploy to GitHub Pages + uses: peaceiris/actions-gh-pages@v3 + if: github.ref == 'refs/heads/main' + with: + github_token: ${{ secrets.GITHUB_TOKEN }} + publish_dir: public From 99234bb6d042a70db81a964431e8843208126fdd Mon Sep 17 00:00:00 2001 From: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> Date: Wed, 1 Sep 2021 20:54:18 -0400 Subject: [PATCH 07/12] Add ephysLoad.m Signed-off-by: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> --- data/neuropixels/ephysLoad.m | 126 +++++++++++++++++++++++++++++++++++ 1 file changed, 126 insertions(+) create mode 100644 data/neuropixels/ephysLoad.m diff --git a/data/neuropixels/ephysLoad.m b/data/neuropixels/ephysLoad.m new file mode 100644 index 0000000..3889f33 --- /dev/null +++ b/data/neuropixels/ephysLoad.m @@ -0,0 +1,126 @@ +ephysroot = '/home/carsen/dm11/data/Spikes/eightprobes/' +% mouse names +mstr = {'Krebs','Waksman','Robbins'}; +% start of spontaneous activity in each mouse +tstart = [3811 3633 3323]; + +% probe location in brain information +% borders tells you region based on depth +% ccfCoords tells you position in microns relative to allen CCF +load(fullfile(ephysroot, 'probeLocations.mat')); + +% possible brain regions (as strings) +areaLabels = {'FrCtx','FrMoCtx','SomMoCtx','SSCtx','V1','V2','RSP',... +'CP','LS','LH','HPF','TH','SC','MB'}; + +%% to plot probe in wire brain, download https://github.com/cortex-lab/allenCCF and npy-matlab +% note that this plot includes every site plotted as a dot - zoom in to see. +addpath(genpath('/github/allenCCF')); +addpath(genpath('/github/npy-matlab')); +plotBrainGrid([], [], [], 0); +hold all; +co = get(gca, 'ColorOrder'); +for imouse = 1:3 + probeColor = co(imouse,:); + for pidx = 1:numel(probeLocations(imouse).probe) + ccfCoords = probeLocations(imouse).probe(pidx).ccfCoords; + + % here we divide by 10 to convert to units of voxels (this atlas is 10um + % voxels, but coordinates are in um) and we swap 3rd with 2nd dimension + % because Allen atlas dimensions are AP/DV/LR, but the wiremesh brain has + % DV as Z, the third dimension (better for view/rotation in matlab). + % So both of these are ultimately quirks of the plotBrainGrid function, + % not of the ccfCoords data + plot3(ccfCoords(:,1)/10, ccfCoords(:,3)/10, ccfCoords(:,2)/10, '.', 'Color', probeColor,'markersize',4) + end +end + +for q = 1:2:360 + view(q, 25); drawnow; +end + +%% +% which mouse +brainLocAll=[]; +for imouse = [1:3] + mouse_name = mstr{imouse}; + + % load spikes + load(fullfile(ephysroot, sprintf('spks/spks%s_Feb18.mat',mouse_name))); + % spks(k).st = spike times (in seconds) + % spks(k).clu = cluster identity of each spike in st (which neuron does spike belong to) + % spks(k).Wheights = height of each cluster on the probe + + % load behavior file + beh = load(fullfile(ephysroot, sprintf('faces/%s_face_proc.mat',mouse_name))); + % movie of mouse face during recording was dimensionality reduced using SVD + % beh.motionSVD = timepoints x components + % beh.motionMask = pixels x pixels x components + % beh.times = times of behavioral frames (in seconds) in same timeframe as spikes + motSVD = beh.motionSVD; + tVid = beh.times; % times of movie frames in spike reference frame + + %% extract spike times and create a matrix neurons x time + stall = zeros(5e3,5500,'uint8'); + ij = 0; + maxt=0; + Wh = []; + iprobe=[]; + brainLoc=[]; + srate = 30; % sampling rate in Hz (how much to bin matrix) + % loop over probes + for k = 1:numel(spks) + clu = spks(k).clu; % cluster ids + st = spks(k).st; % spike times in seconds + st = round(st*srate); % spike times in 30Hz + + % transform spike times into a matrix + % any duplicates of spike times are added together + S = sparse(st, clu, ones(1, numel(clu))); + S = uint8(full(S))'; + % there might be missing numbers (bad spike clusters) + S = S(sort(unique(clu)),:); + + % add these to the big matrix with all probes + stall(ij+[1:size(S,1)],1:size(S,2)) = S; + ij = ij + size(S,1); + maxt = max(maxt, size(S,2)); + + % height of clusters on the probe + % we will use these to determine the brain area + whp = spks(k).Wheights(sort(unique(clu))); + + % borders of brain areas wrt the probe + lowerBorder = probeLocations(imouse).probe(k).borders.lowerBorder; + upperBorder = probeLocations(imouse).probe(k).borders.upperBorder; + acronym = probeLocations(imouse).probe(k).borders.acronym; + loc = zeros(numel(whp),1); + % determine brain area for each cluster based on whp + for j = 1:numel(acronym) + whichArea = find(strcmp(areaLabels, acronym{j})); + loc(whp >= lowerBorder(j) & whp < upperBorder(j)) = whichArea; + end + + % concatenate for all probes + Wh = [Wh; whp]; + brainLoc = [brainLoc; loc]; + iprobe=[iprobe; k * ones(size(S,1),1)]; + end + %% + stall = stall(1:ij, 1:maxt); + + % only use spontaneous time frames + tspont = tstart(imouse)*srate : min(floor(tVid(end)*srate), size(stall,2)-4); + stall = stall(:,tspont); + tspont = tspont / srate; % put back in seconds + + % to put the behavioral data into the spike frames (no time delays) + x = interp1(tVid, motSVD, tspont); + + %%% save the extracted spike matrix and brain locations and faces %%% + %save(fullfile(matroot, sprintf('%swithFaces_KS2.mat',mouse_name)), 'stall','Wh','iprobe',... + %'motSVD','tspont','tVid','srate','brainLoc','areaLabels'); + + brainLocAll=[brainLocAll;brainLoc]; + +end From b42bb6800f8d969bb9eb1918c4e58c5c66aaadd6 Mon Sep 17 00:00:00 2001 From: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> Date: Wed, 1 Sep 2021 20:54:58 -0400 Subject: [PATCH 08/12] Make ephysload usable Signed-off-by: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> --- data/neuropixels/ephysLoad.m | 44 ++++++++++-------------------------- 1 file changed, 12 insertions(+), 32 deletions(-) diff --git a/data/neuropixels/ephysLoad.m b/data/neuropixels/ephysLoad.m index 3889f33..430e8fd 100644 --- a/data/neuropixels/ephysLoad.m +++ b/data/neuropixels/ephysLoad.m @@ -1,6 +1,12 @@ -ephysroot = '/home/carsen/dm11/data/Spikes/eightprobes/' +% Neuropixels data +% Download all the files at https://figshare.com/articles/dataset/Eight-probe_Neuropixels_recordings_during_spontaneous_behaviors/7739750 +% Slightly modified from https://figshare.com/ndownloader/files/14570105 + +ephysroot = '.' % mouse names -mstr = {'Krebs','Waksman','Robbins'}; +% Original: mstr = {'Krebs','Waksman','Robbins'}; +% doing only Waksman +mstr = {'Waksman'}; % start of spontaneous activity in each mouse tstart = [3811 3633 3323]; @@ -13,36 +19,10 @@ areaLabels = {'FrCtx','FrMoCtx','SomMoCtx','SSCtx','V1','V2','RSP',... 'CP','LS','LH','HPF','TH','SC','MB'}; -%% to plot probe in wire brain, download https://github.com/cortex-lab/allenCCF and npy-matlab -% note that this plot includes every site plotted as a dot - zoom in to see. -addpath(genpath('/github/allenCCF')); -addpath(genpath('/github/npy-matlab')); -plotBrainGrid([], [], [], 0); -hold all; -co = get(gca, 'ColorOrder'); -for imouse = 1:3 - probeColor = co(imouse,:); - for pidx = 1:numel(probeLocations(imouse).probe) - ccfCoords = probeLocations(imouse).probe(pidx).ccfCoords; - - % here we divide by 10 to convert to units of voxels (this atlas is 10um - % voxels, but coordinates are in um) and we swap 3rd with 2nd dimension - % because Allen atlas dimensions are AP/DV/LR, but the wiremesh brain has - % DV as Z, the third dimension (better for view/rotation in matlab). - % So both of these are ultimately quirks of the plotBrainGrid function, - % not of the ccfCoords data - plot3(ccfCoords(:,1)/10, ccfCoords(:,3)/10, ccfCoords(:,2)/10, '.', 'Color', probeColor,'markersize',4) - end -end - -for q = 1:2:360 - view(q, 25); drawnow; -end - %% % which mouse brainLocAll=[]; -for imouse = [1:3] +for imouse = [1:length(mstr)] mouse_name = mstr{imouse}; % load spikes @@ -118,9 +98,9 @@ x = interp1(tVid, motSVD, tspont); %%% save the extracted spike matrix and brain locations and faces %%% - %save(fullfile(matroot, sprintf('%swithFaces_KS2.mat',mouse_name)), 'stall','Wh','iprobe',... - %'motSVD','tspont','tVid','srate','brainLoc','areaLabels'); + save(fullfile('.', sprintf('%swithFaces_KS2.mat',mouse_name)), 'stall','Wh','iprobe',... + 'motSVD','tspont','tVid','srate','brainLoc','areaLabels'); brainLocAll=[brainLocAll;brainLoc]; -end +end \ No newline at end of file From 32e78c7b81bf1cd8cd2c9fb753a0dfdc4a591b3d Mon Sep 17 00:00:00 2001 From: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> Date: Wed, 1 Sep 2021 20:57:43 -0400 Subject: [PATCH 09/12] Add dimension reduction Signed-off-by: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> --- .github/workflows/python.yaml | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/.github/workflows/python.yaml b/.github/workflows/python.yaml index 93ace6d..254a83e 100644 --- a/.github/workflows/python.yaml +++ b/.github/workflows/python.yaml @@ -11,6 +11,8 @@ on: jobs: test: runs-on: ubuntu-latest + env: + PATH_NEUROPIXELS: ${{ github.workspace }}/data/neuropixels steps: - uses: actions/checkout@v2 @@ -23,6 +25,37 @@ jobs: - name: Add conda to system path run: echo $CONDA/bin >> $GITHUB_PATH + - name: Cache Data + uses: actions/cache@v2 + with: + path: | + $PATH_NEUROPIXELS/faces.zip + $PATH_NEUROPIXELS/data/spks.zip + key: ${{ runner.os }}-${{ hashFiles('**/neuropixels/*.zip')}} + restore-keys: | + ${{ runner.os }}- + + - name: Download Neuropixels Data + working-directory: ${{ env.PATH_NEUROPIXELS }} + # Neuropixels data see https://figshare.com/articles/dataset/Eight-probe_Neuropixels_recordings_during_spontaneous_behaviors/7739750 + # Only using data from mouse Waksman + run: | + wget -O faces.zip -N https://figshare.com/ndownloader/files/14405048 + wget -O spks.zip -N https://figshare.com/ndownloader/files/14405057 + wget -O probeLocations.mat https://figshare.com/ndownloader/files/14570108 + wget -O probeBorders.mat https://figshare.com/ndownloader/files/14653298 + unzip faces.zip -x faces/Robbins_face_proc.mat faces/Krebs_face_proc.mat + unzip spks.zip -x spks/Robbins_Feb18.mat spks/Krebs_Feb18.mat + + - name: Set up MATLAB + uses: matlab-actions/setup-matlab@v1 + + # Result in WaksmanwithFaces_KS2.mat + - name: Process raw neuropixels data + uses: matlab-actions/run-command@v1 + with: + command: cd('${{ env.PATH_NEUROPIXELS }}'), ephysLoad + - name: Install dependencies run: | conda config --add channels conda-forge @@ -30,6 +63,8 @@ jobs: - name: Run Bubblewrap run: | + python scripts/dimension_reduction_neuropixels.py ${{ env.PATH_NEUROPIXELS }}/WaksmanwithFaces_KS2.mat + conda install jupytext pytest datagen.py jupytext --from py --to ipynb scripts/run_bubblewrap.py --execute --run-path . From ed6ed215af1fb2c1fdf273bbef7e1d74120f0d42 Mon Sep 17 00:00:00 2001 From: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> Date: Wed, 1 Sep 2021 21:40:41 -0400 Subject: [PATCH 10/12] Fix dimension_reduction path Signed-off-by: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> --- scripts/dimension_reduction_neuropixels.py | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/scripts/dimension_reduction_neuropixels.py b/scripts/dimension_reduction_neuropixels.py index b868b3b..53d273f 100644 --- a/scripts/dimension_reduction_neuropixels.py +++ b/scripts/dimension_reduction_neuropixels.py @@ -1,16 +1,23 @@ #%% +import os +import sys + +sys.path.append(os.path.abspath(os.path.dirname("__file__"))) + import time +from pathlib import Path + import numpy as np +import proSVD from scipy import io as sio - from sklearn import random_projection as rp -import proSVD #%% loading processed matlab data # one mouse -dataloc = './' -file = 'WaksmanwithFaces_KS2.mat' -matdict = sio.loadmat(dataloc+file, squeeze_me=True) +assert len(sys.argv) == 2 +file = Path(sys.argv[1]) + +matdict = sio.loadmat(file, squeeze_me=True) spks = matdict['stall'] # truncate so it doesn't take forever @@ -44,6 +51,8 @@ def reduce_sparseRP(X, comps=100, eps=0.1, transformer=None): t = l1 transformer = rp.SparseRandomProjection(n_components=rp_dim) for i in range(num_iters): + if i & 15 == 0: # Check if multiple of 16 + print(f"ssSVD Iteration {i}/{num_iters}") start_rp = time.time() X_rp = reduce_sparseRP(X, transformer=transformer) end_rp = time.time() @@ -75,6 +84,6 @@ def reduce_sparseRP(X, comps=100, eps=0.1, transformer=None): projs[i][:, t:t+l] = curr_basis.T @ curr_neural t += l -np.savez('neuropixel_reduced.npz', ssSVD10=proj_stream_ssSVD) +np.savez(file.parent, ssSVD10=proj_stream_ssSVD) # %% From 0262b41bef6c73978a3a94a4de197760f073a24d Mon Sep 17 00:00:00 2001 From: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> Date: Tue, 7 Sep 2021 18:14:35 -0400 Subject: [PATCH 11/12] Add upload artifact Signed-off-by: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> --- .github/workflows/python.yaml | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/.github/workflows/python.yaml b/.github/workflows/python.yaml index 254a83e..dc7f53c 100644 --- a/.github/workflows/python.yaml +++ b/.github/workflows/python.yaml @@ -1,4 +1,4 @@ -name: Run with Conda +name: Run w/ Synthetic Data on: push: @@ -71,6 +71,12 @@ jobs: jupyter nbconvert --to html scripts/run_bubblewrap.ipynb mkdir public mv scripts/run_bubblewrap.html public/index.html + + - name: Upload Artifact + uses: actions/upload-artifact@v2 + with: + name: Synthetic data run + path: public/index.html - name: Deploy to GitHub Pages uses: peaceiris/actions-gh-pages@v3 From ef736c11f58af3214c612b89441f23e376684207 Mon Sep 17 00:00:00 2001 From: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> Date: Tue, 7 Sep 2021 18:32:49 -0400 Subject: [PATCH 12/12] Deploy in this PR Signed-off-by: Chaichontat Sriworarat <34997334+chaichontat@users.noreply.github.com> --- .github/workflows/python.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/python.yaml b/.github/workflows/python.yaml index dc7f53c..24781c6 100644 --- a/.github/workflows/python.yaml +++ b/.github/workflows/python.yaml @@ -80,7 +80,7 @@ jobs: - name: Deploy to GitHub Pages uses: peaceiris/actions-gh-pages@v3 - if: github.ref == 'refs/heads/main' + if: github.ref == 'refs/heads/main' || github.ref == 'refs/heads/github_actions' with: github_token: ${{ secrets.GITHUB_TOKEN }} publish_dir: public