From a25c4a7d6df24be27647e336e1db23b2801fc10d Mon Sep 17 00:00:00 2001 From: "srinivas.gs" Date: Fri, 16 Dec 2016 14:18:29 -0500 Subject: [PATCH] usability fixes; cleanup --- .DS_Store | Bin 6148 -> 8196 bytes @spikesort/removeDoublets.m | 4 +- README.md | 4 +- build_number | 2 +- console.log | 33 -- docs/README.md | 1 - obsolete/mergeData.m | 154 ------ obsolete/removeArtifactsUsingTemplate.m | 48 -- obsolete/ssdm_2DFracAmpISI.m | 81 --- obsolete/ssdm_2DnormPCA.m | 13 - old_spikesort.m | 671 ------------------------ preCachetSNE.m | 14 - preCachetSNE_engine.m | 103 ---- pref.m | 2 - 14 files changed, 6 insertions(+), 1124 deletions(-) delete mode 100644 console.log delete mode 100644 obsolete/mergeData.m delete mode 100644 obsolete/removeArtifactsUsingTemplate.m delete mode 100644 obsolete/ssdm_2DFracAmpISI.m delete mode 100644 obsolete/ssdm_2DnormPCA.m delete mode 100644 old_spikesort.m delete mode 100644 preCachetSNE.m delete mode 100644 preCachetSNE_engine.m diff --git a/.DS_Store b/.DS_Store index 5008ddfcf53c02e82d7eee2e57c38e5672ef89f6..057cab3a913e069c382ac78919ab6279f8946e47 100644 GIT binary patch literal 8196 zcmeHMPfrs;6n_H+cR}b%`6mZmB}5OHv_?TR1|$RjB?Bx!=-3&?1y(r*sRM_g1V9YobkI=7Iw0`T1;zzdIR>t{Ca)fda3!L| zAlx1GIx|O%3#@YF?m)OZ5b$*4OXZ zidV{7P8RhavL(y({AlzOX=)C&w4PBErCsSzR#NLuCFL7_*-IGiQt->n)2xuwj$`e~ zMt0t$m5kOs?%2Lz+SIKxUw&ZPd8d+a%68rh)=3)vI{K=rs;C)lU}qmw-ds!HEb?9TLRacg_$@$QpXJVGL}1cVzRetlH9n<>|}9J_#e<{N@M zUSs6ItE1ul-=?9T(SR&mAB+8xgcm&MBT0B!A7De^t+)${U_gQ0A=q`RlLulIaLi)4-5M!{w2Y}~ONH({0OQo%LzsgjX1 z?ZSdjANeVBn`X4G2Gcp3f4D!D-YTI@_x?n%+>}eH9jrIIw&qcPp-x_0_An_=`?NT1 zuQ|B>CZ^rll5aY;H&0y;am6yS)S9&uj_*6gMZ$^Z#Xq$U#h_O+DT7znGib?mFySz~rpj)4mf u%!JJcw5KPGg_gAhnpCgF}!R pBnTAa1`@77th%xAJM(0I8AV3M$)+;eJWLRCKt?lcj^~-f3;-`E4{!hg diff --git a/@spikesort/removeDoublets.m b/@spikesort/removeDoublets.m index 65156d0..690975f 100644 --- a/@spikesort/removeDoublets.m +++ b/@spikesort/removeDoublets.m @@ -33,7 +33,7 @@ function removeDoublets(s) end end -if s.pref.ssDebug +if s.verbosity cprintf('green','\n[INFO]') cprintf('text', [' B2A doublet resolution.' oval(length(B2A)) ' spikes swapped']) @@ -65,7 +65,7 @@ function removeDoublets(s) B = sort(unique([B(:); A2B_cand(:)])); A = setdiff(A,A2B_cand); -if s.pref.ssDebug +if s.verbosity cprintf('green','\n[INFO]') cprintf('text', [' B2A doublet resolution.' oval(length(A2B_cand)) ' spikes swapped']) end diff --git a/README.md b/README.md index 83e8a37..fc714c3 100644 --- a/README.md +++ b/README.md @@ -57,7 +57,7 @@ brew install tag This implementation of t-SNE is the the fastest I know of, and gets faster the more cores you throw at it. Follow the installation instructions [here](https://github.com/DmitryUlyanov/Multicore-TSNE#install). macOS users, [see this](https://github.com/DmitryUlyanov/Multicore-TSNE/issues/1#issuecomment-262938483) for installation troubleshooting. -You then need to tell MATLAB where `Multicore-TSNE is installed`. To do this, +You then need to tell MATLAB where `Multicore-TSNE` is installed. To do this, 1. first determine where it is installed (`pip show MulticoreTSNE`) 2. Copy and paste this code to your MATLAB `startup.m` (use `edit startup.m`) @@ -79,6 +79,8 @@ setenv('PATH', path1); ## Architecture +![](images/arch.png) + `spikesort` is built around a plugin architecture for the three most important things it does: * Data handling diff --git a/build_number b/build_number index dda8f88..f520b75 100644 --- a/build_number +++ b/build_number @@ -1 +1 @@ - 247 + 248 diff --git a/console.log b/console.log deleted file mode 100644 index 2ec0da6..0000000 --- a/console.log +++ /dev/null @@ -1,33 +0,0 @@ -30-Nov-2016 17:36:18-spikesort-Loading file:/local-data/gain-flip/pentyl-acetate-2-butanone/ab2//2016_07_12_CS_F2_pentyl-acetate_2-butanone_ab2_S15.mat -13-Dec-2016 17:07:38-spikesort-Loading file:/Volumes/sgs_data/Dropbox (emonetlab)/users/srinivas_gs/data/DA-paper/data-for-paper/paired-pulses//2016_12_13_paired_pulses_2ac_CS_ab3_F1_S1_1.mat -13-Dec-2016 17:08:14-spikesort-Loading file:/Volumes/sgs_data/Dropbox (emonetlab)/users/srinivas_gs/data/DA-paper/data-for-paper/paired-pulses//2016_12_13_paired_pulses_2ac_CS_ab3_F1_S1_1.mat -13-Dec-2016 17:08:37-spikesort-Loading file:/Volumes/sgs_data/Dropbox (emonetlab)/users/srinivas_gs/data/DA-paper/data-for-paper/paired-pulses//2016_12_13_paired_pulses_2ac_CS_ab3_F1_S1_1.mat -13-Dec-2016 17:09:34-spikesort-Loading file:/Volumes/sgs_data/Dropbox (emonetlab)/users/srinivas_gs/data/DA-paper/data-for-paper/paired-pulses//2016_12_13_paired_pulses_2ac_CS_ab3_F1_S2_b0.mat -13-Dec-2016 17:11:15-spikesort-Loading file:/Volumes/sgs_data/Dropbox (emonetlab)/users/srinivas_gs/data/DA-paper/data-for-paper/paired-pulses//2016_12_13_paired_pulses_2ac_CS_ab3_F1_S3_b0.mat -13-Dec-2016 17:11:28-spikesort-Loading file:/Volumes/sgs_data/Dropbox (emonetlab)/users/srinivas_gs/data/DA-paper/data-for-paper/paired-pulses//2016_12_13_paired_pulses_2ac_CS_ab3_F1_S4_b0.mat -13-Dec-2016 17:11:35-spikesort-Loading file:/Volumes/sgs_data/Dropbox (emonetlab)/users/srinivas_gs/data/DA-paper/data-for-paper/paired-pulses//2016_12_13_paired_pulses_2ac_CS_ab3_F1_s1.mat -13-Dec-2016 17:11:54-spikesort-Loading file:/Volumes/sgs_data/Dropbox (emonetlab)/users/srinivas_gs/data/DA-paper/data-for-paper/paired-pulses//2016_12_13_paired_pulses_2ac_CS_ab3_F2_S5_b0.mat -13-Dec-2016 17:12:07-spikesort-Rejecting point: Y exceeded -13-Dec-2016 17:13:16-spikesort-Made a raster plot. -13-Dec-2016 17:13:20-spikesort-Made a firing rate plot. -13-Dec-2016 17:14:13-spikesort-Rejecting point: Y exceeded -13-Dec-2016 17:14:20-spikesort-Made a raster plot. -13-Dec-2016 17:14:24-spikesort-Made a firing rate plot. -13-Dec-2016 17:14:58-spikesort-Made a firing rate plot. -14-Dec-2016 09:09:54-spikesort-Loading file:/Volumes/sgs_data/Dropbox (emonetlab)/users/srinivas_gs/data/DA-paper/data-for-paper/paired-pulses//2016_12_13_paired_pulses_2ac_CS_ab3_F1_S1_1.mat -14-Dec-2016 09:10:22-spikesort-Loading file:/Volumes/sgs_data/Dropbox (emonetlab)/users/srinivas_gs/data/DA-paper/data-for-paper/paired-pulses//2016_12_13_paired_pulses_2ac_CS_ab3_F1_S2_b0.mat -14-Dec-2016 09:10:30-spikesort-Rejecting point: Y exceeded -14-Dec-2016 09:10:56-spikesort-Loading file:/Volumes/sgs_data/Dropbox (emonetlab)/users/srinivas_gs/data/DA-paper/data-for-paper/paired-pulses//2016_12_13_paired_pulses_2ac_CS_ab3_F1_S3_b0.mat -14-Dec-2016 09:11:02-spikesort-Loading file:/Volumes/sgs_data/Dropbox (emonetlab)/users/srinivas_gs/data/DA-paper/data-for-paper/paired-pulses//2016_12_13_paired_pulses_2ac_CS_ab3_F1_S4_b0.mat -14-Dec-2016 09:11:05-spikesort-Loading file:/Volumes/sgs_data/Dropbox (emonetlab)/users/srinivas_gs/data/DA-paper/data-for-paper/paired-pulses//2016_12_13_paired_pulses_2ac_CS_ab3_F1_s1.mat -14-Dec-2016 09:11:12-spikesort-Loading file:/Volumes/sgs_data/Dropbox (emonetlab)/users/srinivas_gs/data/DA-paper/data-for-paper/paired-pulses//2016_12_13_paired_pulses_2ac_CS_ab3_F2_S5_b0.mat -14-Dec-2016 09:11:16-spikesort-Loading file:/Volumes/sgs_data/Dropbox (emonetlab)/users/srinivas_gs/data/DA-paper/data-for-paper/paired-pulses//2016_12_13_paired_pulses_2ac_CS_ab3_F2_S5_b1.mat -14-Dec-2016 09:11:18-spikesort-Loading file:/Volumes/sgs_data/Dropbox (emonetlab)/users/srinivas_gs/data/DA-paper/data-for-paper/paired-pulses//2016_12_13_paired_pulses_2ac_CS_ab3_F2_S5_b0.mat -14-Dec-2016 09:49:51-spikesort-Loading file:/Volumes/sgs_data/Dropbox (emonetlab)/users/srinivas_gs/data/DA-paper/data-for-paper/paired-pulses/ok//2016_12_13_paired_pulses_2ac_CS_ab3_F1_S2_b0.mat -14-Dec-2016 09:54:22-spikesort-Loading file:/Volumes/sgs_data/Dropbox (emonetlab)/users/srinivas_gs/data/DA-paper/data-for-paper/paired-pulses/ok//2016_12_13_paired_pulses_2ac_CS_ab3_F1_S2_b0.mat -14-Dec-2016 09:54:23-spikesort-No need to autosort -14-Dec-2016 09:54:29-spikesort-Loading file:/Volumes/sgs_data/Dropbox (emonetlab)/users/srinivas_gs/data/DA-paper/data-for-paper/paired-pulses/ok//2016_12_13_paired_pulses_2ac_CS_ab3_F2_S5_b0.mat -14-Dec-2016 09:54:31-spikesort-No need to autosort -14-Dec-2016 09:54:33-spikesort-No need to autosort -14-Dec-2016 09:54:34-spikesort-No need to autosort -14-Dec-2016 09:56:27-spikesort-Made a firing rate plot. diff --git a/docs/README.md b/docs/README.md index 0cd457b..fc7c193 100644 --- a/docs/README.md +++ b/docs/README.md @@ -15,7 +15,6 @@ deltat = 1e-4; % what is the time step of the data? %% ~~~~~~~~~~~~~~~~~ GENERAL ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -ssDebug = false; % should spikesort run in debug mode? useFastBandPass = true; % use a fast, FFT-based bandPass? %% ~~~~~~~~~~~~~~~~~ DISPLAY ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ diff --git a/obsolete/mergeData.m b/obsolete/mergeData.m deleted file mode 100644 index 7a19d14..0000000 --- a/obsolete/mergeData.m +++ /dev/null @@ -1,154 +0,0 @@ -% mergeData.m -% merges data from different files into one file -% use this wisely -% this will only merge the following variables: -% -% data -% spikes -% -% WARNING: NOTHING ELSE WILL BE MERGED. -% -% created by Srinivas Gorur-Shandilya at 7:57 , 24 June 2015. Contact me at http://srinivas.gs/contact/ -% -% This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. -% To view a copy of this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/. - -function [] = mergeData() - -allfiles = dir('*.mat'); -if length(allfiles) < 2 - error('WTF am I supposed to merge?') -elseif length(allfiles) > 10 - error('Are you sure?') -end - -disp('Merging the following files:') -disp({allfiles.name}') -% load the first one to get a sense of what this is like -load(allfiles(1).name) - -% reorder fields because of MATLAB stupidity -data = orderfields(data); -ControlParadigm = orderfields(ControlParadigm); - -hash = dataHash(ControlParadigm); - -merged_data = data; -merged_spikes = spikes; -merged_ControlParadigm = ControlParadigm; - -for i = 2:length(allfiles) - load(allfiles(i).name) - - data = orderfields(data); - ControlParadigm = orderfields(ControlParadigm); - - if ~strcmp(dataHash(fieldnames(data)),dataHash(fieldnames(merged_data))) - error('data that I just loaded has variables that I did not expect.') - end - - if ~strcmp(dataHash(fieldnames(spikes)),dataHash(fieldnames(merged_spikes))) - warning('spikes that I just loaded has variables that I did not expect.') - end - - if strcmp(hash,dataHash(ControlParadigm)) - haz_data = find(Kontroller_ntrials(data)); - for j = haz_data - if j > length(merged_data) - error('not coded #39') - else - % merge data - fn = fieldnames(data); - for k = 1:length(fn) - eval(strcat('merged_data(j).',fn{k},'= [merged_data(j).',fn{k} ,' ; data(j).',fn{k},'];')) - end - - % merge spikes - fn = fieldnames(spikes); - for k = 1:length(fn) - try - eval(strcat('merged_spikes(j).',fn{k},'= [merged_spikes(j).',fn{k} ,' ; spikes(j).',fn{k},'];')) - catch er - disp(er.message) - end - end - - end - end - else - % need to match paradigm by paradigm, and make up new ones if needed. - disp('ControlParadigm mismatch. Attempting to fit as best as I can...') - haz_data = find(Kontroller_ntrials(data)); - for j = haz_data - this_hash = dataHash(ControlParadigm(j)); - - % check if this belongs somewhere in the merged master - m = []; - for k = 1:length(merged_ControlParadigm) - if strcmp(dataHash(merged_ControlParadigm(k)),this_hash) - m = k; - disp('Matched') - disp(ControlParadigm(j).Name) - disp('to:') - disp(merged_ControlParadigm(k).Name) - end - end - - - - if isempty(m) - disp('Looks like a new paradigm...') - merged_ControlParadigm(end+1) = ControlParadigm(j); - m = length(merged_ControlParadigm); - end - - % merge data - fn = fieldnames(data); - if m > length(merged_data) - for k = 1:length(fn) - eval((strcat('merged_data(m).',fn{k},'= data(j).',fn{k},';'))) - end - else - for k = 1:length(fn) - eval(strcat('merged_data(m).',fn{k},'= [merged_data(m).',fn{k} ,' ; data(j).',fn{k},'];')) - end - end - - % merge spikes - fn = fieldnames(spikes); - if m > length(merged_spikes) - for k = 1:length(fn) - try - eval((strcat('merged_spikes(m).',fn{k},'= spikes(j).',fn{k},';'))) - catch - warning('Could not merge spikes...') - end - end - else - for k = 1:length(fn) - eval(strcat('merged_spikes(m).',fn{k},'= [merged_spikes(m).',fn{k} ,' ; spikes(j).',fn{k},'];')) - end - end - - - end - - end - -end - -data = merged_data; -spikes = merged_spikes; -ControlParadigm = merged_ControlParadigm; - -l = Inf; -use_this = 1; -for i = 1:length(allfiles) - if length(allfiles(i).name) < l - l = min([l length(allfiles(i).name)]); - use_this = i; - end - -end - -save([allfiles(use_this).name(1:end-4) '_merged.mat'],'ControlParadigm','data','spikes','timestamps','metadata','SamplingRate','OutputChannelNames') \ No newline at end of file diff --git a/obsolete/removeArtifactsUsingTemplate.m b/obsolete/removeArtifactsUsingTemplate.m deleted file mode 100644 index b4f1673..0000000 --- a/obsolete/removeArtifactsUsingTemplate.m +++ /dev/null @@ -1,48 +0,0 @@ -% removeArtifactsUsingTemplate.m -% -% created by Srinivas Gorur-Shandilya at 9:21 , 09 February 2016. Contact me at http://srinivas.gs/contact/ -% -% This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. -% To view a copy of this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/. - -function [VC] = removeArtifactsUsingTemplate(V,this_control,pref) - -if pref.ssDebug - disp('removing artifacts using template match...') -end - -VC = V; - -% read the template file -load('template.mat','template') - -% figure out the control signal -c = this_control(template.control_signal_channels,:); -assert(width(c) == 1, 'too many control signals specified. I expected only one control channel') -c = c(:); - -% find ons and offs and build templates -on_transitions = find(diff(c)==1); -off_transitions = find(diff(c)==-1); - -after = length(template.on_template) - 1; -if isempty(on_transitions) -else - if pref.ssDebug - disp('control signal has transitions, using them to find artifacts...') - end - % trim some edge cases - on_transitions(find(on_transitions+after>(length(V)-1))) = []; - off_transitions(find(off_transitions+after>(length(V)-1))) = []; - - for ti = 1:length(on_transitions) - if pref.use_on_template - VC(on_transitions(ti):on_transitions(ti)+after) = VC(on_transitions(ti):on_transitions(ti)+after) - template.on_template'; - end - if pref.use_off_template - VC(off_transitions(ti):off_transitions(ti)+after) = VC(off_transitions(ti):off_transitions(ti)+after) - template.off_template'; - end - end -end - -keyboard diff --git a/obsolete/ssdm_2DFracAmpISI.m b/obsolete/ssdm_2DFracAmpISI.m deleted file mode 100644 index 0c1f67c..0000000 --- a/obsolete/ssdm_2DFracAmpISI.m +++ /dev/null @@ -1,81 +0,0 @@ -% ssdm_2DFracAmpISI.m -% -% this is a plugin for spikesort.m -% reduces spikes to a amplitude, measured from the minimum to preceding maximum. -% part of the spikesort package -% https://github.com/sg-s/spikesort -% created by Srinivas Gorur-Shandilya at 10:20 , 09 April 2014. Contact me at http://srinivas.gs/contact/ -% -% This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. -% To view a copy of this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/. -function R = ssdm_2DFracAmpISI(V,deltat,loc,ax,ax2) - - -wb = waitbar(0.2,'Computing Fractional amplitudes...'); - -h = (40*1e-4)/deltat; % deltat in seconds -% 1D - find total spike amplitude for each -R = zeros*loc; -loc_max = 0*loc; -[R(1),loc_max(1)] = max(V(loc(1)-h:loc(1)) - V(loc(1))); -loc_max(1) = loc(1) + loc_max(1) - h; -for i = 2:length(loc) - before = max([loc(i)-h loc(i-1)]); - [R(i),loc_max(i)] = max(V(before:loc(i)) - V(loc(i))); - loc_max(i) = loc_max(i) + before; -end - - -time = deltat:deltat:(deltat*length(V)); -% figure, hold on -% plot(time,V) -% scatter(time(loc),V(loc)) -% scatter(time(loc_max),V(loc_max)) - - - -waitbar(0.4,wb,'Building spike envelopes...'); -% build an upper and lower envelope - -upper_envelope = interp1(time(loc_max),V(loc_max),time); -upper_envelope(1:find(~isnan(upper_envelope),1,'first')) = upper_envelope(find(~isnan(upper_envelope),1,'first')); -upper_envelope((find(~isnan(upper_envelope),1,'last')):end) = upper_envelope(find(~isnan(upper_envelope),1,'last')-1); -lower_envelope = interp1(time(loc),V(loc),time); -lower_envelope(1:find(~isnan(lower_envelope),1,'first')) = lower_envelope(find(~isnan(lower_envelope),1,'first')); -lower_envelope((find(~isnan(lower_envelope),1,'last')):end) = lower_envelope(find(~isnan(lower_envelope),1,'last')-1); -% plot(ax,time,lower_envelope,'g') -% plot(ax,time,upper_envelope,'r') - -waitbar(0.6,wb,'Estimating spike density...'); -% build a time-varying estimate of ISI -t_prev_spike = [0 diff(loc)]; -t_next_spike = [diff(loc) 0]; -t_closest_spike = (t_next_spike + t_prev_spike)/2; -t_closest_spike(1) = t_next_spike(1); -t_closest_spike(end) = t_prev_spike(end); -isi = interp1(time(loc),t_closest_spike,time); -isi(1:find(~isnan(isi),1,'first')) = isi(find(~isnan(isi),1,'first')); -isi((find(~isnan(isi),1,'last')):end) = isi(find(~isnan(isi),1,'last')-1); -isi = filtfilt(ones(1,h)/h,1,isi); clear t_closest_spike t_prev_spike t_next_spike - - -waitbar(0.8,wb,'Filtering...'); -% filter the envelopes in a time-dependant manner -upper_envelope2 = upper_envelope; -lower_envelope2 = lower_envelope; -scaling_factor = 2; -for i = loc - before = floor(max([1 i-isi(i)*scaling_factor])); - after = floor(min([length(isi) i+isi(i)*scaling_factor])); - upper_envelope2(i) = max(upper_envelope(before:after)); - lower_envelope2(i) = min(lower_envelope(before:after)); - -end -clear upper_envelope lower_envelope -envelope_amplitude = upper_envelope2- lower_envelope2; - -R = R./envelope_amplitude(loc); -close(wb) - -R = [R; isi(loc)]; - diff --git a/obsolete/ssdm_2DnormPCA.m b/obsolete/ssdm_2DnormPCA.m deleted file mode 100644 index 8131098..0000000 --- a/obsolete/ssdm_2DnormPCA.m +++ /dev/null @@ -1,13 +0,0 @@ -% ssdm_normV -% -% created by Srinivas Gorur-Shandilya at 3:27 , 15 April 2015. Contact me at http://srinivas.gs/contact/ -% -% This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. -% To view a copy of this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/. - -function R = ssdm_2DnormPCA(V_snippets) -for i = 1:size(V_snippets,2) - V_snippets(:,i) = V_snippets(:,i)/min(V_snippets(:,i)); -end -[~,R]=princomp(V_snippets'); -R = R(:,1:2)'; diff --git a/old_spikesort.m b/old_spikesort.m deleted file mode 100644 index 106c8e4..0000000 --- a/old_spikesort.m +++ /dev/null @@ -1,671 +0,0 @@ -% old spikesort - - -%% begin subfunctions -% all subfunctions here are listed alphabetically - - function autosortCallback(~,~) - if autosort_control.Value == 1 - autosort_control.FontWeight = 'bold'; - autosort_control.String = 'AUTOSORT ON'; - else - autosort_control.FontWeight = 'normal'; - autosort_control.String = 'Autosort Off'; - end - end - - - - function [A,B,N] = autosort() - % automatically sorts data, if possible. - reduceDimensionsCallback; - [A,B,N]=findCluster; - - end - - function [] = matchTemplate(~,~) - % template match - nchannels = length(get(handles.valve_channel,'Value')); - plot_these = get(handles.valve_channel,'Value'); - control_signal = ControlParadigm(s.this_paradigm).Outputs(plot_these,:); - - if size(control_signal,2) > size(control_signal,1) - control_signal = control_signal'; - end - - % make the template object - template.control_signal_channels = plot_these; - - figure, hold on - for i = 1:width(control_signal) - temp = control_signal(:,i); - % find ons and offs and build templates - on_transitions = find(diff(temp)==1); - off_transitions = find(diff(temp)==-1); - - after = round(pref.template_width); - if isnan(after) || after < 11 - after = 50; - end - - if isempty(on_transitions) - else - % trim some edge cases - on_transitions(find(on_transitions+after>(length(V)-1))) = []; - off_transitions(find(off_transitions+after>(length(V)-1))) = []; - - on_template = zeros(after+1,length(on_transitions)); - off_template = zeros(after+1,length(on_transitions)); - for ti = 1:length(on_transitions) - snippet = V(on_transitions(ti):on_transitions(ti)+after); - snippet = snippet - snippet(1); - on_template(:,ti) = snippet; - - snippet = V(off_transitions(ti):off_transitions(ti)+after); - snippet = snippet - snippet(1); - off_template(:,ti) = snippet(:); - end - off_template = (mean(off_template,2)); - on_template = (mean(on_template,2)); - - subplot(width(control_signal),2,2*(i-1)+1); - plot(on_template) - title('On Template') - - subplot(width(control_signal),2,2*(i-1)+2); - plot(off_template) - title('Off Template') - - template(i).on_template = on_template; - template(i).off_template = off_template; - - end - - save('template.mat','template'); - end - end - - - function [] = removeArtifacts(~,~) - if strcmp(get(handles.remove_artifacts_menu,'Checked'),'off') - set(handles.remove_artifacts_menu,'Checked','on'); - else - set(handles.remove_artifacts_menu,'Checked','off'); - end - plotResp; - end - - - - - - - function discard(~,~) - - - if get(discard_control,'Value') == 0 - % reset discard - if isfield(spikes,'discard') - spikes(s.this_paradigm).discard(ThisTrial) = 0; - end - set(discard_control,'String','Discard','FontWeight','normal') - else - set(discard_control,'String','Discarded!','FontWeight','bold') - - % need to reset spikes - if length(spikes) >= s.this_paradigm - if width(spikes(s.this_paradigm).A) >= ThisTrial - spikes(s.this_paradigm).A(ThisTrial,:) = 0; - spikes(s.this_paradigm).B(ThisTrial,:) = 0; - spikes(s.this_paradigm).amplitudes_A(ThisTrial,:) = 0; - spikes(s.this_paradigm).amplitudes_B(ThisTrial,:) = 0; - - else - % all cool - end - else - % should have no problem - end - - % mark as discarded - spikes(s.this_paradigm).discard(ThisTrial) = 1; - save(strcat(path_name,file_name),'spikes','-append') - - end - - - % update screen - plotResp; - end - - - - - function [A,B,N] = findCluster(~,~) - % cluster based on the method - methodname = get(cluster_control,'String'); - method = get(cluster_control,'Value'); - methodname = strcat('sscm_',methodname{method}); - req_arg = argInNames(methodname); % find out what arguments the external method needs - % start constructing the eval string - es = strcat('[A,B,N]=',methodname,'('); - for ri = 1:length(req_arg) - es = strcat(es,req_arg{ri},','); - end - clear ri - es = es(1:end-1); - es = strcat(es,');'); - try - eval(es); - catch exc - ms = strcat(methodname, ' ran into an error: ', exc.message); - msgbox(ms,'spikesort'); - disp('The full stack is:') - for ei = 1:length(exc.stack) - disp(exc.stack(ei)) - end - return - end - clear es - - - - % mark them - set(handles.ax1_A_spikes,'XData',time(A),'YData',V(A),'Marker','o','MarkerSize',pref.marker_size,'Parent',handles.ax1,'MarkerEdgeColor','r','LineStyle','none'); - set(handles.ax1_B_spikes,'XData',time(B),'YData',V(B),'Marker','o','MarkerSize',pref.marker_size,'Parent',handles.ax1,'MarkerEdgeColor','b','LineStyle','none'); - set(handles.ax1_all_spikes,'XData',NaN,'YData',NaN); - - % remove noise spikes from the loc vector - loc = setdiff(loc,N); - - % save them - try - spikes(s.this_paradigm).A(ThisTrial,:) = sparse(1,length(time)); - spikes(s.this_paradigm).B(ThisTrial,:) = sparse(1,length(time)); - spikes(s.this_paradigm).N(ThisTrial,:) = sparse(1,length(time)); - spikes(s.this_paradigm).amplitudes_A(ThisTrial,:) = sparse(1,length(time)); - spikes(s.this_paradigm).amplitudes_B(ThisTrial,:) = sparse(1,length(time)); - catch - spikes(s.this_paradigm).A = sparse(ThisTrial,length(time)); - spikes(s.this_paradigm).B = sparse(ThisTrial,length(time)); - spikes(s.this_paradigm).N = sparse(ThisTrial,length(time)); - spikes(s.this_paradigm).amplitudes_A = sparse(ThisTrial,length(time)); - spikes(s.this_paradigm).amplitudes_B = sparse(ThisTrial,length(time)); - - end - - - % also save spike amplitudes - try - spikes(s.this_paradigm).amplitudes_A(ThisTrial,A) = ssdm_1DAmplitudes(V,A); - catch - end - try - spikes(s.this_paradigm).amplitudes_B(ThisTrial,B) = ssdm_1DAmplitudes(V,B); - catch - end - - % save them - save(strcat(path_name,file_name),'spikes','-append') - - end - - - - - - function generateSummary(~,~) - allfiles = dir(strcat(path_name,'*.mat')); - if any(find(strcmp('cached.mat',{allfiles.name}))) - allfiles(find(strcmp('cached.mat',{allfiles.name}))) = []; - end - if any(find(strcmp('cached_log.mat',{allfiles.name}))) - allfiles(find(strcmp('cached_log.mat',{allfiles.name}))) = []; - end - summary_string = ''; - fileID = fopen('summary.log','w'); - for i = 1:length(allfiles) - summary_string = strcat(summary_string,'\n', allfiles(i).name); - temp = load(allfiles(i).name,'metadata'); - metadata = temp.metadata; - if size(metadata.spikesort_comment,1) > 1 - metadata.spikesort_comment = metadata.spikesort_comment(1,:); - end - if isfield(metadata,'spikesort_comment') - summary_string = strcat(summary_string,'\t\t', metadata.spikesort_comment); - else - % no comment on this file - summary_string = strcat(summary_string,'\t\t', 'no comment'); - end - - end - - fprintf(fileID,summary_string); - fclose(fileID); - end - - function jump(src,~) - % get the digital channels - digital_channels = get(handles.valve_channel,'Value'); - - % find out where we are - xl= floor(get(handles.ax1,'XLim')/pref.deltat); - - - if src == jump_fwd - next_on = Inf; - - % find the next digital channel switch in any channel - for i = 1:length(digital_channels) - this_channel = ControlParadigm(s.this_paradigm).Outputs(digital_channels(i),:); - [ons] = computeOnsOffs(this_channel); - ons(onsxl(1)-1) = []; - prev_on = max([prev_on(:); ons(:)]); - end - if ~isinf(-prev_on) - set(handles.ax1,'Xlim',[time(prev_on) time(prev_on+diff(xl))]); - end - else - error('Unknown source of call to jump'); - end - end - - function loadFileCallback(src,~) - if strcmp(get(src,'String'),'Load File') - [file_name,path_name,filter_index] = uigetfile({'*.mat';'*.kontroller'}); - if ~file_name - return - end - elseif strcmp(get(src,'String'),'<') - if isempty(file_name) - return - else - % first save what we had before - save(strcat(path_name,file_name),'spikes','-append') - - - if filter_index == 1 - allfiles = dir(strcat(path_name,'*.mat')); - else - allfiles = dir(strcat(path_name,'*.kontroller')); - end - thisfile = find(strcmp(file_name,{allfiles.name}))-1; - if thisfile < 1 - file_name = allfiles(end).name; - else - file_name = allfiles(thisfile).name; - end - - end - else - if isempty(file_name) - return - else - % first save what we had before - save(strcat(path_name,file_name),'spikes','-append') - - if filter_index == 1 - allfiles = dir(strcat(path_name,'*.mat')); - else - allfiles = dir(strcat(path_name,'*.kontroller')); - end - thisfile = find(strcmp(file_name,{allfiles.name}))+1; - if thisfile > length(allfiles) - file_name = allfiles(1).name; - else - file_name = allfiles(thisfile).name; - end - - end - end - - % reset some pushbuttons and other things - set(discard_control,'Value',0) - s.this_paradigm = 1; - ThisTrial = 1; - temp = []; - clear spikes - spikes.A = 0; - spikes.B = 0; - spikes.artifacts = 0; - R = 0; % this holds the dimensionality reduced data - V = 0; % holds the current trace - Vf = 0; % filtered V - V_snippets = []; - time = 0; - loc =0; % holds current spike times - - console(strcat('Loading file:',path_name,'/',file_name)) - - temp=load(strcat(path_name,file_name),'-mat'); - try - try - delete(handles.load_waitbar) - catch - end - handles.load_waitbar = waitbar(0.2, 'Loading data...'); - ControlParadigm = temp.ControlParadigm; - data = temp.data; - SamplingRate = temp.SamplingRate; - OutputChannelNames = temp.OutputChannelNames; - try - metadata = temp.metadata; - timestamps = temp.timestamps; - catch - end - if isfield(temp,'spikes') - spikes = temp.spikes; - end - clear temp - - - - waitbar(0.3,handles.load_waitbar,'Updating listboxes...') - % update control signal listboxes with OutputChannelNames - set(handles.valve_channel,'String',OutputChannelNames) - - % update stimulus listbox with all input channel names - fl = fieldnames(data); - - % also add all the control signals - try - set(stim_channel,'String',[fl(:); OutputChannelNames]); - catch - set(stim_channel,'String',[fl(:) OutputChannelNames]); - end - - % update response listbox with all the input channel names - set(resp_channel,'String',fl); - - % some sanity checks - if length(data) > length(ControlParadigm) - error('Something is wrong with this file: more data than control paradigms.') - end - - % find out which paradigms have data - n = structureElementLength(data); - - % only show the paradigms with data - temp = {ControlParadigm.Name}; - set(paradigm_chooser,'String',temp(find(n)),'Value',1); - - - % go to the first paradigm with data. - s.this_paradigm = find(n); - s.this_paradigm = s.this_paradigm(1); - - - n = n(s.this_paradigm); - if n - temp ={}; - for i = 1:n - temp{i} = strcat('Trial-',mat2str(i)); - end - set(trial_chooser,'String',temp); - ThisTrial = 1; - set(trial_chooser,'String',temp); - else - set(trial_chooser,'String','No data'); - ThisTrial = NaN; - end - - waitbar(0.4,handles.load_waitbar,'Guessing control signals...') - % automatically default to picking the digital signals as the control signals - digital_channels = zeros(1,length(OutputChannelNames)); - for i = 1:length(ControlParadigm) - for j = 1:width(ControlParadigm(i).Outputs) - uv = (unique(ControlParadigm(i).Outputs(j,:))); - if length(uv) == 2 && sum(uv) == 1 - digital_channels(j) = 1; - end - - end - end - digital_channels = find(digital_channels); - set(handles.valve_channel,'Value',digital_channels); - - - waitbar(0.5,handles.load_waitbar,'Guessing stimulus and response...') - temp = find(strcmp('PID', fl)); - if ~isempty(temp) - set(stim_channel,'Value',temp); - end - temp = find(strcmp('voltage', fl)); - if ~isempty(temp) - set(resp_channel,'Value',temp); - - end - - set(handles.main_fig,'Name',strcat(versionname,'--',file_name)) - - % enable all controls - waitbar(.7,handles.load_waitbar,'Enabling UI...') - set(handles.sine_control,'Enable','on'); - set(autosort_control,'Enable','on'); - set(redo_control,'Enable','on'); - set(findmode,'Enable','on'); - set(filtermode,'Enable','on'); - set(cluster_control,'Enable','on'); - set(prev_trial,'Enable','on'); - set(next_trial,'Enable','on'); - set(prev_paradigm,'Enable','on'); - set(next_paradigm,'Enable','on'); - set(trial_chooser,'Enable','on'); - set(paradigm_chooser,'Enable','on'); - set(discard_control,'Enable','on'); - set(metadata_text_control,'Enable','on') - - % check for amplitudes - waitbar(.7,handles.load_waitbar,'Checking to see amplitude data exists...') - % check if we have spike_amplitude data - if length(spikes) - for i = 1:length(spikes) - for j = 1:width(spikes(i).A) - haz_data = 0; - if length(spikes(i).A(j,:)) > 2 - if isfield(spikes,'discard') - if length(spikes(i).discard) < j - haz_data = 1; - else - if ~spikes(i).discard(j) - haz_data = 1; - end - end - else - haz_data = 1; - end - end - if haz_data - recompute = 0; - if isfield(spikes,'amplitudes_A') - if width(spikes(i).amplitudes_A) < j - recompute = 1; - spikes(i).amplitudes_A = []; - spikes(i).amplitudes_B = []; - elseif length(spikes(i).amplitudes_A(j,:)) < length(spikes(i).A(j,:)) - spikes(i).amplitudes_A = []; - spikes(i).amplitudes_B = []; - recompute = 1; - - end - end - if recompute - A = spikes(i).A(j,:); - - spikes(i).amplitudes_A(j,:) = sparse(1,length(A)); - spikes(i).amplitudes_B(j,:) = sparse(1,length(A)); - V = data(i).voltage(j,:); - spikes(i).amplitudes_A(j,find(A)) = ssdm_1DAmplitudes(V,pref.deltat,find(A),pref.invert_V); - B = spikes(i).B(j,:); - spikes(i).amplitudes_B(j,find(B)) = ssdm_1DAmplitudes(V,pref.deltat,find(B),pref.invert_V); - end - end - - end - end - end - - % check to see if metadata exists - try - set(metadata_text_control,'String',metadata.spikesort_comment) - catch - set(metadata_text_control,'String','') - end - - % check to see if this file is tagged. - if isunix - clear es - es{1} = 'tag -l '; - es{2} = strcat(path_name,file_name); - [~,temp] = unix(strjoin(es)); - set(tag_control,'String',temp(strfind(temp,'.mat')+5:end-1)); - end - - % clean up - close(handles.load_waitbar) - - plotStim; - plotResp(@loadFileCallback); - catch err - if strcmp(get(src,'String'),'>') - loadFileCallback(src) - elseif strcmp(get(src,'String'),'<') - loadFileCallback(src) - else - warning('Something went wrong with loading the file. The error was:') - disp(err) - disp('The error was here:') - for ei = 1:length(err.stack) - disp(err.stack(ei)) - end - end - - end - end - - - - function resetZoom(~,~) - set(handles.ax1,'XLim',[min(time) max(time)]); - - temp = sort(V); - handles.ax1.YLim(1) = mean(temp(1:floor(length(temp)*.05)))*3; - temp = sort(V,'descend'); - handles.ax1.YLim(2) = mean(temp(1:floor(length(temp)*.05)))*3; - - end - - - function markAllCallback(~,~) - % get view - xmin = get(handles.ax1,'XLim'); - xmin = xmin/pref.deltat; - xmax = xmin(2); xmin=xmin(1); - - % get mode - if get(mode_B2A,'Value') - % add to A spikes - spikes(s.this_paradigm).A(ThisTrial,loc(loc>xmin & locxmin & locxmin & locxmin & locxmin & locxmin & locxmin & loc length(V) - xl(2) = length(V); - end - - % check if we already have some discard information stored in spikes - if length(spikes) < s.this_paradigm - spikes(s.this_paradigm).use_trace_fragment = ones(1,length(V)); - else - if isfield(spikes,'use_trace_fragment') - if width(spikes(s.this_paradigm).use_trace_fragment) < ThisTrial - spikes(s.this_paradigm).use_trace_fragment(ThisTrial,:) = ones(1,length(V)); - else - - end - else - spikes(s.this_paradigm).use_trace_fragment(ThisTrial,:) = ones(1,length(V)); - end - end - - - if strcmp(get(src,'String'),'Discard View') - spikes(s.this_paradigm).use_trace_fragment(ThisTrial,xl(1):xl(2)) = 0; - % disp('Discarding view for trial #') - % disp(ThisTrial) - % disp('Discarding data from:') - % disp(xl*pref.deltat) - elseif strcmp(get(src,'String'),'Retain View') - spikes(s.this_paradigm).use_trace_fragment(ThisTrial,xl(1):xl(2)) = 1; - else - error('modifyTraceDiscard ran into an error because I was called by a function that I did not expect. I am meant to be called only by the discard view or the retain view pushbuttons.') - end - - plotResp(@modifyTraceDiscard); - end - - - - - function plotValve(~,~) - % get the channels to plot - handles.valve_channels = get(handles.valve_channel,'Value'); - c = jet(length(handles.valve_channels)); - for i = 1:length(handles.valve_channels) - this_valve = ControlParadigm(s.this_paradigm).Outputs(handles.valve_channels(i),:); - end - plotStim; - end - - - - - function templateMatch(~,~) - plotResp(@templateMatch); - end - - function updateMetadata(src,~) - metadata.spikesort_comment = get(src,'String'); - save(strcat(path_name,file_name),'metadata','-append') - end - - - - diff --git a/preCachetSNE.m b/preCachetSNE.m deleted file mode 100644 index 6e68d9b..0000000 --- a/preCachetSNE.m +++ /dev/null @@ -1,14 +0,0 @@ -function [] = preCachetSNE() - -if license('test', 'Distrib_Computing_Toolbox') - if ~length(gcp('nocreate')) - else - delete(gcp) - end - p = parpool(1); - f = parfeval(p,@preCachetSNE_engine,0); - return -else - preCachetSNE_engine; - return -end diff --git a/preCachetSNE_engine.m b/preCachetSNE_engine.m deleted file mode 100644 index 873164f..0000000 --- a/preCachetSNE_engine.m +++ /dev/null @@ -1,103 +0,0 @@ -% preCachetSNE_engine.m -% this function attempts to pre-calculate t-SNE embeddings for all the data, so that the actual process of spike sorting is faster and less annoying -% -% usage: -% cd /folder/with/data/from/kontroller -% preCachetSNE('invert_V',true,'variable_name','voltage') -% created by Srinivas Gorur-Shandilya at 9:35 , 11 September 2015. Contact me at http://srinivas.gs/contact/ -% -% This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. -% To view a copy of this license, visit http://creativecommons.org/licenses/by-nc-sa/4.0/. - -function [] = preCachetSNE_engine() - -variable_name = 'voltage'; - -% use pref.m to change how this function behaves. - -pref = readPref(fileparts(which(mfilename))); -dbm = ['The hash of the pref. file I read is ' dataHash(pref)]; -system(['echo "' dbm '" >> spikesort.log']); - -% add src to path -% add src folder to path -addpath([fileparts(which(mfilename)) oss 'src']) - - -allfiles = dir('*.mat'); - -for i = 1:length(allfiles) - thisfile = allfiles(i).name; - dbm = ['Working on: ' thisfile]; - system(['echo "' dbm '" >> spikesort.log']); - if ~strcmp(thisfile,'consolidated_data.mat') && ~strcmp(thisfile,'cached.mat') && ~strcmp(thisfile,'cached_log.mat') - load(thisfile) - for j = 1:length(data) - this_control = ControlParadigm(j).Outputs; - if eval(['~isempty(data(j).' variable_name ')']) - this_data = eval(['(data(j).' variable_name ')']); - - for k = 1:width(this_data) - try - V = this_data(k,:); - - % use templates to remove artifacts - if exist([pwd oss 'template.mat'],'file') - disp('Template exists...') - if pref.use_on_template || pref.use_off_template - V = removeArtifactsUsingTemplate(V,this_control,pref); - end - end - - - lc = 1/pref.band_pass(1); - lc = floor(lc/pref.deltat); - hc = 1/pref.band_pass(2); - hc = floor(hc/pref.deltat); - - if pref.useFastBandPass - [V,Vf] = fastBandPass(V,lc,hc); - else - [V,Vf] = bandPass(V,lc,hc); - end - - % find spikes - loc = findSpikes(V); - - % take snippets for each putative spike - V_snippets = NaN(pref.t_before+pref.t_after,length(loc)); - if loc(1) < pref.t_before+1 - loc(1) = []; - V_snippets(:,1) = []; - end - if loc(end) + pref.t_after+1 > length(V) - loc(end) = []; - V_snippets(:,end) = []; - end - for l = 1:length(loc) - V_snippets(:,l) = V(loc(l)-pref.t_before+1:loc(l)+pref.t_after); - end - - if pref.ssDebug - disp('We have these many V_snippets') - disp(length(V_snippets)) - end - - if exist('fast_tsne','file') - % run the fast tSNE algorithm on this - dbm = ['starting fast_tsne @ ' datestr(now) '.working on paradigm: ' oval(j) ' , trial: ' oval(k)]; - system(['echo "' dbm '" >> spikesort.log']); - dbm = ['The hash of this data is: ' dataHash(V_snippets)]; - system(['echo "' dbm '" >> spikesort.log']); - fast_tsne(V_snippets,2,10,60,.5); - else - error('You need bhtsne to use this.') - end - catch err - disp(err) - end - end - end - end - end -end diff --git a/pref.m b/pref.m index 844a6b1..b8fd36c 100644 --- a/pref.m +++ b/pref.m @@ -15,7 +15,6 @@ %% ~~~~~~~~~~~~~~~~~ GENERAL ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -ssDebug = true; % should spikesort run in debug mode? useFastBandPass = false; % use a fast, FFT-based bandPass? %% ~~~~~~~~~~~~~~~~~ DISPLAY ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -82,4 +81,3 @@ multicore_tsne_path = '~/anaconda3/bin'; % change this to the path where MultiCoreTSNE is installed -