From 0e74639dfbf22a67b53fc497d4667f8a29c19d69 Mon Sep 17 00:00:00 2001 From: cowenowner Date: Mon, 9 Nov 2020 11:04:35 -0700 Subject: [PATCH 1/6] minor changes. initialized a variable and defined some values before a loop to slightly increase efficiency. --- .gitignore | 3 + Programs_and_data/Main_assemblies_detection.m | 54 ++-- Programs_and_data/TestPair.m | 240 +++++++++--------- 3 files changed, 151 insertions(+), 146 deletions(-) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..144333c --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ + +Programs_and_data/TestPair.asv +Programs_and_data/Main_assemblies_detection.asv diff --git a/Programs_and_data/Main_assemblies_detection.m b/Programs_and_data/Main_assemblies_detection.m index ec686a4..bd2eee4 100644 --- a/Programs_and_data/Main_assemblies_detection.m +++ b/Programs_and_data/Main_assemblies_detection.m @@ -1,34 +1,34 @@ function [assembly]=Main_assemblies_detection(spM,MaxLags,BinSizes,ref_lag,alph,No_th,O_th,bytelimit) -% this function returns cell assemblies detected in spM spike matrix binned -% at a temporal resolution specified in 'BinSizes' vector and testing for all +% this function returns cell assemblies detected in spM spike matrix binned +% at a temporal resolution specified in 'BinSizes' vector and testing for all % lags between '-MaxLags(i)' and 'MaxLags(i)' % % USAGE: [assembly]=Main_assemblies_detection(spM, MaxLags, BinSizes, ref_lag, alph, Dc, No_th, O_th, bytelimit) % % ARGUMENTS: -% spM := matrix with population spike trains; each row is the spike train (time stamps, not binned) relative to a unit. +% spM := matrix with population spike trains; each row is the spike train (time stamps, not binned) relative to a unit. % BinSizes:= vector of bin sizes to be tested; % MaxLags:= vector of maximal lags to be tested. For a binning dimension of BinSizes(i) the program will test all pairs configurations with a time shift between -MaxLags(i) and MaxLags(i); % (optional) ref_lag := reference lag. Default value 2 % (optional) alph := alpha level. Default value 0.05 % (optional) No_th := minimal number of occurrences required for an assembly (all assemblies, even if significant, with fewer occurrences than No_th are discarded). Default value 0. % (optional) O_th := maximal assembly order (the algorithm will return assemblies of composed by maximum O_th elements). -% (optional) bytelimit := maximal size (in bytes) allocated for all assembly structures detected with a bin dimension. When the size limit is reached the algorithm stops adding new units. +% (optional) bytelimit := maximal size (in bytes) allocated for all assembly structures detected with a bin dimension. When the size limit is reached the algorithm stops adding new units. % % RETURNS: % assembly - structure containing assembly information: % assembly.parameters - parameters used to run Main_assemblies_detection -% assembly.bin{i} contains information about assemblies detected with -% 'BinSizes(i)' bin size tested for all lags between +% assembly.bin{i} contains information about assemblies detected with +% 'BinSizes(i)' bin size tested for all lags between % '-MaxLags(i)' and 'MaxLags(i)' -% +% % assembly.bin{i}.bin_edges - bin edges (common to all assemblies in assembly.bin{i}) -% assembly.bin{i}.n{j} information about the j-th assembly detected with BinSizes(i) bin size +% assembly.bin{i}.n{j} information about the j-th assembly detected with BinSizes(i) bin size % elements: vector of units taking part to the assembly (unit order correspond to the agglomeration order) % lag: vector of time lags. '.lag(z)' is the activation delay between .elements(1) and .elements(z+1) % pr: vector of pvalues. '.pr(z)' is the pvalue of the statistical test between performed adding .elements(z+1) to the structure .elements(1:z) % Time: assembly activation time. It reports how many times the complete assembly activates in that bin. .Time always refers to the activation of the first listed assembly element (.elements(1)), that doesn't necessarily corresponds to the first unit firing. -% Noccurrences: number of assembly occurrence. '.Noccurrences(z)' is the occurrence number of the structure composed by the units .elements(1:z+1) +% Noccurrences: number of assembly occurrence. '.Noccurrences(z)' is the occurrence number of the structure composed by the units .elements(1:z+1) % % % @@ -38,40 +38,40 @@ % last update 11/01/2016 % last update 13/02/2016 TestPair.m update -if nargin<4 || isempty(ref_lag), ref_lag=2; end; -if nargin<5 || isempty(alph), alph=0.05; end; -if nargin<6 || isempty(No_th), No_th=0; end; % no limitation on the number of assembly occurrences -if nargin<7 || isempty(O_th), O_th=Inf; end; % no limitation on the assembly order (=number of elements in the assembly) -if nargin<8 || isempty(bytelimit), bytelimit=Inf; end; % no limitation on assembly dimension +if nargin<4 || isempty(ref_lag), ref_lag=2; end +if nargin<5 || isempty(alph), alph=0.05; end +if nargin<6 || isempty(No_th), No_th=0; end % no limitation on the number of assembly occurrences +if nargin<7 || isempty(O_th), O_th=Inf; end % no limitation on the assembly order (=number of elements in the assembly) +if nargin<8 || isempty(bytelimit), bytelimit=Inf; end % no limitation on assembly dimension nneu=size(spM,1); % number of units -assemblybin=cell(1,length(BinSizes)); +assemblybin=cell(1,length(BinSizes)); Dc=100; %length (in # bins) of the segments in which the spike train is divided to compute #abba variance (parameter k). +start_t = nanmin(spM(:)); +end_t = nanmax(spM(:)); -parfor gg=1:length(BinSizes) - - +for gg=1:length(BinSizes) % Cowen: parfor does not really save much time. int=BinSizes(gg); maxlag=MaxLags(gg); - fprintf('%d - testing: bin size=%f sec; max tested lag=%d \n', gg, int, maxlag); + fprintf('%d - testing: bin size=%f sec; max tested lag=%d \n', gg, int, maxlag); %%%%%%%%%%%%%%%%%%%%% Binning %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % the spike train is binned with temporal resolution 'int' - tb=min(spM(:)):int:max(spM(:)); + tb=start_t:int:end_t; binM=zeros(nneu,length(tb)-1,'uint8'); for n=1:nneu - [ binM(n,:),~] = histcounts(spM(n,:),tb); - end + binM(n,:) = histcounts(spM(n,:),tb); + end if size(binM,2)-MaxLags(gg)<100 fprintf('Warning: testing bin size=%f. The time series is too short, consider taking a longer portion of spike train or diminish the bin size to be tested \n', int); else %%%%%%%%%%%%%%%%%%%%%%% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - [assemblybin{gg}]=FindAssemblies_recursive(binM,maxlag,alph,gg,Dc,No_th,O_th,bytelimit,ref_lag); % it returns assemblies detected at specific temporal resolution 'int' + [assemblybin{gg}]=FindAssemblies_recursive(binM,maxlag,alph,gg,Dc,No_th,O_th,bytelimit,ref_lag); % it returns assemblies detected at specific temporal resolution 'int' if ~isempty(assemblybin{gg}) assemblybin{gg}.bin_edges=tb; end - fprintf('%d - testing done\n',gg); + fprintf('%d - testing done\n',gg); fname = sprintf('assembly%d.mat', gg); parsave(fname,assemblybin{gg}) end @@ -84,8 +84,8 @@ assembly.parameters.bytelimit=bytelimit; fprintf('\n'); -end +end function parsave(fname,aus) - save(fname,'aus') -end +save(fname,'aus') +end diff --git a/Programs_and_data/TestPair.m b/Programs_and_data/TestPair.m index f107020..40733d8 100644 --- a/Programs_and_data/TestPair.m +++ b/Programs_and_data/TestPair.m @@ -4,10 +4,10 @@ % ensemble := structure with the previously formed assembly and its spike train % spikeTrain2 := spike train of the new unit to be tested for significance (candidate to be a new assembly member) -% n2 := new unit tested +% n2 := new unit tested % maxlag := maximum lag to be tested % Dc := size (in bins) of chunks in which I divide the spike trains to compute the variance (to reduce non stationarity effects on variance estimation) -% reference_lag := lag of reference; if zero or negative reference lag=-l +% reference_lag := lag of reference; if zero or negative reference lag=-l % % % @@ -15,54 +15,56 @@ % for information please contact eleonora.russo@zi-mannheim.de; daniel.durstewitz@zi-mannheim.de. % % last update 11/01/2016 -% last update 13/02/2016 added correction for continuity approximation (ExpAB limit and Yates correction) +% last update 13/02/2016 added correction for continuity approximation (ExpAB limit and Yates correction) +% Cowen - minor improvements 9/11/2020 - couple=[ensemble.Time-min(ensemble.Time(:)) spikeTrain2-min(spikeTrain2(:))]; % spike train pair I am going to test - - nu=2; - ntp=size(couple,1); %%%% trial length - - % I divide in parallel trials with 0/1 elements - maxrate=max(max(couple)); - Zaa=cell(1,maxrate); - for i=1:maxrate - Zaa{i}=zeros(size(couple), 'uint8'); - Zaa{i}(couple>=i)=1; - ExpABi(i)=prod(sum(Zaa{i},1))/size(couple,1); - end - - % % % decide which is the lag with most coincidences (l_:=best lag) - ctAB=nan(1,maxlag+1); - ctAB_=nan(1,maxlag+1); - for l=0:maxlag - trAB=[couple(1:end-maxlag,1), couple(l+1:end-maxlag+l,2)]; - trBA=[couple(l+1:end-maxlag+l,1), couple(1:end-maxlag,2)]; - ctAB(l+1)=sum(min(trAB')); - ctAB_(l+1)=sum(min(trBA')); - end - - if reference_lag<=0 - aus=[ctAB; ctAB_]; - [a,b]=max(aus(:)); - [I,J] = ind2sub(size(aus),b); - l_=(I==1)*(J-1)-(I==2)*(J-1); %% I select l_ +couple=[ensemble.Time-min(ensemble.Time(:)) spikeTrain2-min(spikeTrain2(:))]; % spike train pair I am going to test + +nu=2; +ntp=size(couple,1); %%%% trial length + +% I divide in parallel trials with 0/1 elements +maxrate=max(max(couple)); +Zaa=cell(1,maxrate); +ExpABi = nan(maxrate,1); % Cowen - added this as sometimes maxrate is nothing. +for i=1:maxrate + Zaa{i}=zeros(size(couple), 'uint8'); + Zaa{i}(couple>=i)=1; + ExpABi(i)=prod(sum(Zaa{i},1))/size(couple,1); +end + +% % % decide which is the lag with most coincidences (l_:=best lag) +ctAB=nan(1,maxlag+1); +ctAB_=nan(1,maxlag+1); +for l=0:maxlag + trAB=[couple(1:end-maxlag,1), couple(l+1:end-maxlag+l,2)]; + trBA=[couple(l+1:end-maxlag+l,1), couple(1:end-maxlag,2)]; + ctAB(l+1)=sum(min(trAB,[],2)); + ctAB_(l+1)=sum(min(trBA,[],2)); +end + +if reference_lag<=0 + aus=[ctAB; ctAB_]; + [a,b]=max(aus(:)); + [I,J] = ind2sub(size(aus),b); + l_=(I==1)*(J-1)-(I==2)*(J-1); %% I select l_ +else + Hab_l=[ctAB_(end:-1:2), ctAB]; + [a,b]=max(Hab_l(:)); + lags=-maxlag:maxlag; + l_=lags(b); + Hab=a; + if l_<0 + l_ref=l_+reference_lag; + Hab_ref=Hab_l(lags==l_ref); else - Hab_l=[ctAB_(end:-1:2), ctAB]; - [a,b]=max(Hab_l(:)); - lags=-maxlag:maxlag; - l_=lags(b); - Hab=a; - if l_<0 - l_ref=l_+reference_lag; - Hab_ref=Hab_l(find(lags==l_ref)); - else - l_ref=l_-reference_lag; - Hab_ref=Hab_l(find(lags==l_ref)); - end + l_ref=l_-reference_lag; + Hab_ref=Hab_l(lags==l_ref); end - +end + + - %% HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH @@ -75,112 +77,112 @@ assemD.Time=[]; assemD.Noccurrences=[ensemble.Noccurrences 0]; else - + % % % construct the activation time series for the couple len=size(couple,1); %%%% trial length Time=uint8(zeros(len,1)); %%%% activation vector - if reference_lag<=0 + if reference_lag<=0 if l_==0 - for i=1:maxrate + for i=1:maxrate sA=Zaa{i}(:,1); sB=Zaa{i}(:,2); Time(1:len)=Time(1:len)+sA(1:end).*sB(1:end); end TPrMTot=[0 ctAB(1); ctAB_(3) 0]; % matrix with #AB and #BA - elseif l_>0 - for i=1:maxrate + elseif l_>0 + for i=1:maxrate sA=Zaa{i}(:,1); - sB=Zaa{i}(:,2); + sB=Zaa{i}(:,2); Time(1:len-l_)=Time(1:len-l_)+sA(1:end-l_).*sB(l_+1:end); end - TPrMTot=[0 ctAB(J); ctAB_(J) 0]; % matrix with #AB and #BA + TPrMTot=[0 ctAB(J); ctAB_(J) 0]; % matrix with #AB and #BA else - for i=1:maxrate + for i=1:maxrate sA=Zaa{i}(:,1); - sB=Zaa{i}(:,2); + sB=Zaa{i}(:,2); Time(-l_+1:end)=Time(-l_+1:end)+sA(-l_+1:end).*sB(1:end+l_); end - TPrMTot=[0 ctAB(J); ctAB_(J) 0]; % matrix with #AB and #BA + TPrMTot=[0 ctAB(J); ctAB_(J) 0]; % matrix with #AB and #BA end - - else - + + else + if l_==0 - for i=1:maxrate + for i=1:maxrate sA=Zaa{i}(:,1); sB=Zaa{i}(:,2); Time(1:len)=Time(1:len)+sA(1:end).*sB(1:end); - end - elseif l_>0 - for i=1:maxrate + end + elseif l_>0 + for i=1:maxrate sA=Zaa{i}(:,1); - sB=Zaa{i}(:,2); + sB=Zaa{i}(:,2); Time(1:len-l_)=Time(1:len-l_)+sA(1:end-l_).*sB(l_+1:end); end - + else - for i=1:maxrate + for i=1:maxrate sA=Zaa{i}(:,1); - sB=Zaa{i}(:,2); + sB=Zaa{i}(:,2); Time(-l_+1:end)=Time(-l_+1:end)+sA(-l_+1:end).*sB(1:end+l_); end end - TPrMTot=[0 Hab; Hab_ref 0]; % matrix with #AB and #BA + TPrMTot=[0 Hab; Hab_ref 0]; % matrix with #AB and #BA end -%% --------------------------------------------------------------------% + %% --------------------------------------------------------------------% % % % HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH % I cut the spike train in stationary segments %%%%%%%%%%%%%%%%%%%%% chunking %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nch=ceil((size(couple,1)-maxlag)/Dc); - Dc=floor((size(couple,1)-maxlag)/nch); %% new chunk size, this is to have all chunks of rougly the same size + Dc=floor((size(couple,1)-maxlag)/nch); %% new chunk size, this is to have all chunks of rougly the same size chunked=cell(nch,1); % in order to take into account the same time series part that I used % for MargPr I cut the time series according to l_: couple_cut=nan((size(couple,1)-maxlag),2); - + if l_==0 - couple_cut=couple(1:end-maxlag,:); - elseif l_>0 + couple_cut=couple(1:end-maxlag,:); + elseif l_>0 couple_cut(:,1)=couple(1:end-maxlag,1); couple_cut(:,2)=couple(l_+1:end-maxlag+l_,2); else couple_cut(:,1)=couple(1-l_:end-maxlag-l_,1); couple_cut(:,2)=couple(1:end-maxlag,2); end - - + + for iii=1:nch-1 chunked{iii}=couple_cut((1+Dc*(iii-1)):Dc*iii,:); end chunked{nch}=couple_cut((1+Dc*(nch-1)):end,:); % last chunk can be of slightly different size - -%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - + + %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + MargPr_t=cell(nch,maxrate); %% number of spikes in each parallel process in each unit maxrate_t=nan(nch,1); ch_nn=nan(nch,1); - for iii=1:nch - couple_t=chunked{iii}; + for iii=1:nch + couple_t=chunked{iii}; maxrate_t(iii)=max(max(couple_t)); - ch_nn(iii)=size(chunked{iii},1); + ch_nn(iii)=size(chunked{iii},1); Zaa_t=cell(1,maxrate_t(iii)); for i=1:maxrate_t(iii) Zaa_t{i}=zeros(size(couple_t), 'uint8'); Zaa_t{i}(couple_t>=i)=1; end - - for i=1:maxrate_t(iii) + + for i=1:maxrate_t(iii) sA=Zaa_t{i}(:,1); - sB=Zaa_t{i}(:,2); - MargPr_t{iii,i}=[sum(sA); sum(sB)]; - end + sB=Zaa_t{i}(:,2); + MargPr_t{iii,i}=[sum(sA); sum(sB)]; + end end %%%%%%%%%%%%%%%%%%%%% chunks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - + n=ntp-maxlag; Mx0=cell(nch,maxrate); @@ -190,41 +192,41 @@ covX=cell(1,nch); varX=cell(1,nch); varXtot=zeros(2); - for iii=1:nch + for iii=1:nch maxrate_t=max(chunked{iii}(:)); - ch_n=ch_nn(iii); + ch_n=ch_nn(iii); % % % HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH - % % % evaluation of #AB - for i=1:maxrate_t + % % % evaluation of #AB + for i=1:maxrate_t if ~isempty(MargPr_t{iii,i}) - Mx0{iii,i}= MargPr_t{iii,i} *ones(1,2); + Mx0{iii,i}= MargPr_t{iii,i} *ones(1,2); end - end - + end + % % variance & covariance varT{iii}=zeros(nu); covABAB{iii}=cell(maxrate_t,maxrate_t); - for i=1:maxrate_t - Mx0{iii,i}= MargPr_t{iii,i} *ones(1,2); - covABAB{iii}{i,i}=(Mx0{iii,i}.*Mx0{iii,i}'./ch_n).*(ch_n-Mx0{iii,i}).*(ch_n-Mx0{iii,i}')./(ch_n*(ch_n-1)); - varT{iii}=varT{iii}+covABAB{iii}{i,i}; - for j=(i+1):maxrate_t - covABAB{iii}{i,j}=2*(Mx0{iii,j}.*Mx0{iii,j}'./ch_n).*(ch_n-Mx0{iii,i}).*(ch_n-Mx0{iii,i}')./(ch_n*(ch_n-1)); - varT{iii}=varT{iii}+covABAB{iii}{i,j}; + for i=1:maxrate_t + Mx0{iii,i}= MargPr_t{iii,i} *ones(1,2); + covABAB{iii}{i,i}=(Mx0{iii,i}.*Mx0{iii,i}'./ch_n).*(ch_n-Mx0{iii,i}).*(ch_n-Mx0{iii,i}')./(ch_n*(ch_n-1)); + varT{iii}=varT{iii}+covABAB{iii}{i,i}; + for j=(i+1):maxrate_t + covABAB{iii}{i,j}=2*(Mx0{iii,j}.*Mx0{iii,j}'./ch_n).*(ch_n-Mx0{iii,i}).*(ch_n-Mx0{iii,i}')./(ch_n*(ch_n-1)); + varT{iii}=varT{iii}+covABAB{iii}{i,j}; end end %HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH - % % evaluation of X=#AB-#BA + % % evaluation of X=#AB-#BA covX{iii}=zeros(nu); covABBA{iii}=cell(maxrate_t,maxrate_t); - - for i=1:maxrate_t + + for i=1:maxrate_t covABBA{iii}{i,i}=(Mx0{iii,i}.*Mx0{iii,i}'./ch_n).*(ch_n-Mx0{iii,i}).*(ch_n-Mx0{iii,i}')./(ch_n*(ch_n-1)^2); - covX{iii}=covX{iii}+covABBA{iii}{i,i}; - for j=(i+1):maxrate_t - covABBA{iii}{i,j}=2*(Mx0{iii,j}.*Mx0{iii,j}'./ch_n).*(ch_n-Mx0{iii,i}).*(ch_n-Mx0{iii,i}')./(ch_n*(ch_n-1)^2); - covX{iii}=covX{iii}+covABBA{iii}{i,j}; + covX{iii}=covX{iii}+covABBA{iii}{i,i}; + for j=(i+1):maxrate_t + covABBA{iii}{i,j}=2*(Mx0{iii,j}.*Mx0{iii,j}'./ch_n).*(ch_n-Mx0{iii,i}).*(ch_n-Mx0{iii,i}')./(ch_n*(ch_n-1)^2); + covX{iii}=covX{iii}+covABBA{iii}{i,j}; end end @@ -236,31 +238,31 @@ %HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH X=TPrMTot-TPrMTot'; if abs(X(1,2))>0 - X=abs(TPrMTot-TPrMTot')-0.5; %Yates correction + X=abs(TPrMTot-TPrMTot')-0.5; %Yates correction end if varXtot(1,2)==0 % if variance is zero - prF=1; + prF=1; else F=X.^2./varXtot; - prF=fcdf(F(1,2),1,n,'upper'); -% prF=fcdf(F(1,2),1,2*double(maxrate)*n-1,'upper'); + prF=fcdf(F(1,2),1,n,'upper'); + % prF=fcdf(F(1,2),1,2*double(maxrate)*n-1,'upper'); end - - -%% + + + %% %HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH %HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH - %All information about the assembly and test are returned + %All information about the assembly and test are returned assemD.elements=[ensemble.elements n2]; assemD.lag=[ensemble.lag, l_]; assemD.pr=[ensemble.pr prF]; assemD.Time=Time; assemD.Noccurrences=[ensemble.Noccurrences sum(Time)]; - - + + end %HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH %HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH From eed1755836920f1e98833a35df46499b8d9d2bbf Mon Sep 17 00:00:00 2001 From: cowenowner Date: Mon, 9 Nov 2020 11:08:40 -0700 Subject: [PATCH 2/6] Update .gitignore --- .gitignore | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/.gitignore b/.gitignore index 144333c..9affd4a 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,2 @@ +*.asv -Programs_and_data/TestPair.asv -Programs_and_data/Main_assemblies_detection.asv From 6d446a49089bc3d3ee91a75435de0476518378fe Mon Sep 17 00:00:00 2001 From: cowenowner Date: Mon, 9 Nov 2020 13:03:05 -0700 Subject: [PATCH 3/6] Update Main_assemblies_detection.m --- Programs_and_data/Main_assemblies_detection.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Programs_and_data/Main_assemblies_detection.m b/Programs_and_data/Main_assemblies_detection.m index bd2eee4..c30ce40 100644 --- a/Programs_and_data/Main_assemblies_detection.m +++ b/Programs_and_data/Main_assemblies_detection.m @@ -50,7 +50,7 @@ start_t = nanmin(spM(:)); end_t = nanmax(spM(:)); -for gg=1:length(BinSizes) % Cowen: parfor does not really save much time. +parfor gg=1:length(BinSizes) % Cowen: parfor does not really save much time - well - guess that depends on the number of binsizes explored.. int=BinSizes(gg); maxlag=MaxLags(gg); fprintf('%d - testing: bin size=%f sec; max tested lag=%d \n', gg, int, maxlag); From 30cf39edeb04ebaeb03bb1440815dac7d6f3a3df Mon Sep 17 00:00:00 2001 From: cowenowner Date: Tue, 17 Nov 2020 11:07:38 -0700 Subject: [PATCH 4/6] mods to allow custom binning --- .../Assembly_activity_function.m | 6 +- Programs_and_data/Main_assemblies_detection.m | 32 +++-- .../assembly_assignment_matrix.m | 117 +++++++++--------- 3 files changed, 85 insertions(+), 70 deletions(-) diff --git a/Programs_and_data/Assembly_activity_function.m b/Programs_and_data/Assembly_activity_function.m index 5fb8deb..f9ef762 100644 --- a/Programs_and_data/Assembly_activity_function.m +++ b/Programs_and_data/Assembly_activity_function.m @@ -35,6 +35,8 @@ % % last update 11/01/2016 % +% Cowen: more work to do here to get this to work with custion binnned +% spike time matrices. if nargin<5 || isempty(lagChoice), lagChoice='duration'; end; @@ -45,8 +47,8 @@ assembly_activity=cell(length(as),1); for n=1:length(as) - - tb=assembly.bin{find(BinSizes==as{n}.bin)}.bin_edges; + ix = find(BinSizes==as{n}.bin); + tb=assembly.bin{ix}.bin_edges; switch act_count case 'full' diff --git a/Programs_and_data/Main_assemblies_detection.m b/Programs_and_data/Main_assemblies_detection.m index c30ce40..31364e7 100644 --- a/Programs_and_data/Main_assemblies_detection.m +++ b/Programs_and_data/Main_assemblies_detection.m @@ -37,30 +37,42 @@ % % last update 11/01/2016 % last update 13/02/2016 TestPair.m update +% +% Cowen 2020: Added the option of passing in a cell array of pre-binned data +% instead of the spM which is binned in this function. This allows different binning approaches to be used +% (e.g., by phase or variable bin sizes - that you can't do with +% histcounts). if nargin<4 || isempty(ref_lag), ref_lag=2; end if nargin<5 || isempty(alph), alph=0.05; end if nargin<6 || isempty(No_th), No_th=0; end % no limitation on the number of assembly occurrences if nargin<7 || isempty(O_th), O_th=Inf; end % no limitation on the assembly order (=number of elements in the assembly) if nargin<8 || isempty(bytelimit), bytelimit=Inf; end % no limitation on assembly dimension - -nneu=size(spM,1); % number of units +if ~iscell(spM) + nneu=size(spM,1); % number of units + start_t = nanmin(spM(:)); + end_t = nanmax(spM(:)); +else + nneu=size(spM{1},1); % number of units +end assemblybin=cell(1,length(BinSizes)); Dc=100; %length (in # bins) of the segments in which the spike train is divided to compute #abba variance (parameter k). -start_t = nanmin(spM(:)); -end_t = nanmax(spM(:)); -parfor gg=1:length(BinSizes) % Cowen: parfor does not really save much time - well - guess that depends on the number of binsizes explored.. +for gg=1:length(BinSizes) % Cowen: parfor does not really save much time - well - guess that depends on the number of binsizes explored.. int=BinSizes(gg); maxlag=MaxLags(gg); fprintf('%d - testing: bin size=%f sec; max tested lag=%d \n', gg, int, maxlag); %%%%%%%%%%%%%%%%%%%%% Binning %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % the spike train is binned with temporal resolution 'int' - tb=start_t:int:end_t; - binM=zeros(nneu,length(tb)-1,'uint8'); - - for n=1:nneu - binM(n,:) = histcounts(spM(n,:),tb); + if ~iscell(spM) + tb=start_t:int:end_t; + binM=zeros(nneu,length(tb)-1,'uint8'); + for n=1:nneu + binM(n,:) = histcounts(spM(n,:),tb); + end + else + tb=0:size(spM{gg},2); + binM = spM{gg}; end if size(binM,2)-MaxLags(gg)<100 diff --git a/Programs_and_data/assembly_assignment_matrix.m b/Programs_and_data/assembly_assignment_matrix.m index 68d7779..eb692e9 100644 --- a/Programs_and_data/assembly_assignment_matrix.m +++ b/Programs_and_data/assembly_assignment_matrix.m @@ -1,4 +1,4 @@ -function [Amatrix,Binvector,Unit_order,As_order]=assembly_assignment_matrix(As_across_bins,nneu,BinSizes,display) +function [Amatrix,Binvector,Unit_order,As_order]=assembly_assignment_matrix(As_across_bins,nneu,BinSizes,display,plot_it) % Display the Assembly assignment matrix. % USAGE: [Amatrix,Binvector,Unit_order]=assembly_assignment_matrix(As_across_bins,nneu,BinSizes,display) % @@ -7,15 +7,15 @@ % nneu - # of units % BinSizes - list of investigated bin sizes % display - Possible display options: -% 1) display='raw' -> no rearrangement made, assemblies and units displayed as +% 1) display='raw' -> no rearrangement made, assemblies and units displayed as % from the detection algorithm -% 2) display='ordunit' -> Assembly order is left unchanged, units rearranged +% 2) display='ordunit' -> Assembly order is left unchanged, units rearranged % in order to separate assembly units from units not taking part to any assembly. % 3) display='clustered' -> Assembly and unit order changed in order to % display close to one other assemblies with similar elements % % RETURNS: -% Amatrix - ASSEMBLY ASSIGNMENT MATRIX: each column is an assembly +% Amatrix - ASSEMBLY ASSIGNMENT MATRIX: each column is an assembly % each raw a unit. If unit i takes part in assembly j, Amatrix(i,j) displays % the firing delay (lag) of i with respect to the first assembly unit. % Binvector - bin size. Binvector(j) is the bin size at which @@ -30,7 +30,11 @@ % % last update 11/01/2016 %% -if nargin<4 , display='raw'; end; +if nargin < 5 + plot_it = true; +end + +if nargin<4 , display='raw'; end if nargin == 4 if ~strcmp(display,'raw') && ~strcmp(display,'ordunit') && ~strcmp(display,'clustered') fprintf('''raw'', ''ordunit'', ''clustered'' are the only acceptable style specifications\n') @@ -38,9 +42,6 @@ end end - - - nAss_final=size(As_across_bins,2); AAT=nan(nneu,nAss_final); Binvector=nan(1,nAss_final); @@ -51,7 +52,7 @@ end switch display - + case 'raw' Unit_order=1:nneu; As_order=1:length(As_across_bins); @@ -59,8 +60,8 @@ case 'ordunit' aus=all(isnan(AAT),2); idx_nan=find(aus==1); - idx_activ=find(aus==0); - AAT=[AAT(~~sum(~isnan(AAT),2),:);AAT(~sum(~isnan(AAT),2),:)]; + idx_activ=find(aus==0); + AAT=[AAT(~~sum(~isnan(AAT),2),:);AAT(~sum(~isnan(AAT),2),:)]; Unit_order=[idx_activ;idx_nan]; As_order=1:length(As_across_bins); @@ -70,16 +71,16 @@ idx_activ=find(aus==0); aus=~all(isnan(AAT),2); AAT_activ=AAT(aus,:); - - %%% order on the base of units co-occurrence + + %%% order on the base of units co-occurrence A01=~isnan(AAT_activ); M_assemb=zeros(size(AAT_activ,1)); - + for n=1:size(A01,2) aus=find(A01(:,n)==1); M_assemb(aus,aus)=M_assemb(aus,aus)+1; end - + M_assemb(M_assemb==0)=0.0001; d_assemb=ones(size(AAT_activ,1))./M_assemb; d_assemb=d_assemb-diag(diag(d_assemb)); @@ -92,13 +93,13 @@ Z = squareform(D); Q2=linkage(Z,'average'); [~,~,perm2]=dendrogram(Q2,0); - + AAT=AAT_activ(perm,perm2); AAT=[AAT; nan(length(idx_nan),size(AAT,2))]; Unit_order=[perm1;idx_nan]; Binvector=Binvector(perm2); % I apply the same order on the bin size vector As_order=perm2; - + end Amatrix=AAT; @@ -107,49 +108,49 @@ binmat(end+1,end+1)=0; - - -subplot(2,1,1) - -h=pcolor(AAT); -set(h, 'EdgeColor', [0.8 0.8 0.8]); -caxis([0, max(AAT(:))]) -colormap jet -ylabel('Unit #') -set(gca, 'XTick', []); -hcb=colorbar; -hcb.Label.String = 'Time lag \it l \rm (# bins)'; -yy = 1:size(AAT,1); -set(gca,'YTick',yy(1)+0.5:yy(end),'Yticklabel',Unit_order) -set(gcf, 'Color', [1,1,1]); - - - - - - -subplot(2,1,2) - -pcolor(binmat) -yy = 1:size(AAT,2); -set(gca,'XTick',yy(1)+0.5:2:yy(end)+1+0.05,'Xticklabel',As_order(1:2:end)) -xlabel('Assembly #') -set(gca, 'YTick', []); -hC = colorbar; -if length(BinSizes)>1 - caxis([log10(min(BinSizes)), log10(max(BinSizes))]); -else - caxis([log10(BinSizes-0.001), log10(BinSizes+0.001)]); +if plot_it + + subplot(2,1,1) + + h=pcolor(AAT); + set(h, 'EdgeColor', [0.8 0.8 0.8]); + caxis([0, max(AAT(:))]) + colormap jet + ylabel('Unit #') + set(gca, 'XTick', []); + hcb=colorbar; + hcb.Label.String = 'Time lag \it l \rm (# bins)'; + yy = 1:size(AAT,1); + set(gca,'YTick',yy(1)+0.5:yy(end),'Yticklabel',Unit_order) + set(gcf, 'Color', [1,1,1]); + + + + + + + subplot(2,1,2) + + pcolor(binmat) + yy = 1:size(AAT,2); + set(gca,'XTick',yy(1)+0.5:2:yy(end)+1+0.05,'Xticklabel',As_order(1:2:end)) + xlabel('Assembly #') + set(gca, 'YTick', []); + hC = colorbar; + if length(BinSizes)>1 + caxis([log10(min(BinSizes)), log10(max(BinSizes))]); + else + caxis([log10(BinSizes-0.001), log10(BinSizes+0.001)]); + end + L=[0.001,0.002,0.003,0.004,0.005, 0.006,0.007,0.008,0.009,0.01,0.02,0.03,... + 0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,... + 1,2,3,4,5,6,7,8,9,10]; + set(hC,'Ytick',log10(L),'YTicklabel',L); + set(hC,'Location','southoutside') + hC.Label.String = 'Temporal precision \Delta (sec)'; + end -L=[0.001,0.002,0.003,0.004,0.005, 0.006,0.007,0.008,0.009,0.01,0.02,0.03,... - 0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,... - 1,2,3,4,5,6,7,8,9,10]; -set(hC,'Ytick',log10(L),'YTicklabel',L); -set(hC,'Location','southoutside') -hC.Label.String = 'Temporal precision \Delta (sec)'; - -end From 7771c9187b51e82dddce3a69a9ad1d155ef65607 Mon Sep 17 00:00:00 2001 From: cowenowner Date: Thu, 19 Nov 2020 14:53:12 -0700 Subject: [PATCH 5/6] updates --- .../Assembly_activity_function.m | 187 ++++++++++-------- Programs_and_data/Main_assemblies_detection.m | 6 +- 2 files changed, 104 insertions(+), 89 deletions(-) diff --git a/Programs_and_data/Assembly_activity_function.m b/Programs_and_data/Assembly_activity_function.m index f9ef762..3e000d9 100644 --- a/Programs_and_data/Assembly_activity_function.m +++ b/Programs_and_data/Assembly_activity_function.m @@ -4,14 +4,14 @@ % ARGUMENTS: % As_across_bins := structure containing assembly information output of ASSEMBLIES_ACROSS_BINS; % assembly := structure containing assembly information output of Main_assemblies_detection.m -% spM := matrix with population spike trains; each row is the spike train (time stamps, not binned) relative to a unit. +% spM := matrix with population spike trains; each row is the spike train (time stamps, not binned) relative to a unit. % BinSizes := vector of bin sizes to be tested; -% -% (optional) lagChoice - The user can chose if to consider only the assembly activation (marked at -% the occurrence of the first assembly unit spike) or the activity of the whole +% +% (optional) lagChoice - The user can chose if to consider only the assembly activation (marked at +% the occurrence of the first assembly unit spike) or the activity of the whole % assembly (from the spike of the first assembly unit to the spike of the last assembly unit) -% lagChoice='beginning' ---> only assembly activation -% lagChoice='duration' ---> activity stays up during the whole assembly duration +% lagChoice='beginning' ---> only assembly activation +% lagChoice='duration' ---> activity stays up during the whole assembly duration % (optional) act_count - options on how to count activity @@ -25,8 +25,8 @@ % RETURNS: -% assembly_activity - assembly activation profile across recording time -% according to different possible definitions. +% assembly_activity - assembly activation profile across recording time +% according to different possible definitions. % % % @@ -39,64 +39,79 @@ % spike time matrices. -if nargin<5 || isempty(lagChoice), lagChoice='duration'; end; -if nargin<6 || isempty(act_count), act_count='full'; end; +if nargin<5 || isempty(lagChoice), lagChoice='duration'; end +if nargin<6 || isempty(act_count), act_count='full'; end as=As_across_bins; assembly_activity=cell(length(as),1); - + for n=1:length(as) - ix = find(BinSizes==as{n}.bin); - tb=assembly.bin{ix}.bin_edges; - - switch act_count - case 'full' - [activity]=DetectAssembly(spM ,as{n},tb); % only full assembly occurrences are considered - case 'partial' - [activity]=DetectAssembly_additive(spM ,as{n},tb); % partial activations are also took into account - case 'combined' - [activity]=DetectAssembly_additive_weighted(spM ,as{n},tb); % partial activations are also took into account but weighting differently full assembly contributions from partial assembly contributions. - end - ta=tb(1:end-1)+as{n}.bin/2; %tb are the binning edges, I want the center - WOO=[ta',activity']; - - - - switch lagChoice - case 'beginning' - assembly_activity{n}=WOO; - case 'duration' - - alag=As_across_bins{n}.lag(end); - activity_lag=activity; - for j=1:alag + ix = find(BinSizes==as{n}.bin); + tb=assembly.bin{ix}.bin_edges; + + switch act_count + case 'full' + if iscell(spM) + [activity]=DetectAssembly(spM{ix} ,as{n}); % only full assembly occurrences are considered + else + [activity]=DetectAssembly(spM ,as{n},tb); % only full assembly occurrences are considered + end + case 'partial' + if iscell(spM) + [activity]=DetectAssembly_additive(spM{ix} ,as{n}); % only full assembly occurrences are considered + else + [activity]=DetectAssembly_additive(spM ,as{n},tb); % partial activations are also took into account + end + case 'combined' + if iscell(spM) + [activity]=DetectAssembly_additive_weighted(spM{ix} ,as{n}); % only full assembly occurrences are considered + else + [activity]=DetectAssembly_additive_weighted(spM ,as{n},tb); % partial activations are also took into account but weighting differently full assembly contributions from partial assembly contributions. + end + end + ta=tb(1:end-1)+as{n}.bin/2; %tb are the binning edges, I want the center + WOO=[ta',activity']; + + + + switch lagChoice + case 'beginning' + assembly_activity{n}=WOO; + case 'duration' + + alag=As_across_bins{n}.lag(end); + activity_lag=activity; + for j=1:alag activity_lag=[activity_lag; [zeros(1,j), activity(1:end-j)]]; - end - VOO=[ta',sum(activity_lag,1)']; - - assembly_activity{n}=VOO; - end - + end + VOO=[ta',sum(activity_lag,1)']; + + assembly_activity{n}=VOO; + end + end end - - - - function [activity]=DetectAssembly(spM, Assembly, tb) % only occurrences of full assemblies are considered. The 'activity' vector % returns the number of times the complete set of assembly units is firing -% in a bin. +% in a bin. eelements=Assembly.elements; alag=Assembly.lag; -Sbin_assembly=zeros(length(eelements),length(tb)-1); -for i=1:length(eelements) - [Sbin_assembly(i,:),~] = histcounts(spM(eelements(i),:),tb); +if nargin < 3 + Sbin_assembly=zeros(length(eelements),size(spM,2)); + for i=1:length(eelements) + Sbin_assembly(i,:) = spM(eelements(i),:); + end +else + Sbin_assembly=zeros(length(eelements),length(tb)-1); + for i=1:length(eelements) + [Sbin_assembly(i,:),~] = histcounts(spM(eelements(i),:),tb); + end end Sbin_shifted=nan(size(Sbin_assembly)); @@ -104,14 +119,14 @@ if alag(e1)==0 Sbin_shifted(e1,:)=Sbin_assembly(e1,:); elseif alag(e1)>0 - Sbin_shifted(e1,1:end-alag(e1))=Sbin_assembly(e1,1+alag(e1):end); + Sbin_shifted(e1,1:end-alag(e1))=Sbin_assembly(e1,1+alag(e1):end); else Sbin_shifted(e1,1-alag(e1):end)=Sbin_assembly(e1,1:end+alag(e1)); end end aus=sum(Sbin_shifted,1); Sbin_shifted(:,isnan(aus))=0; % cut at the borders -activity=min(Sbin_shifted,[],1); +activity=min(Sbin_shifted,[],1); end @@ -119,16 +134,24 @@ eelements=Assembly.elements; alag=Assembly.lag; -Sbin_assembly=zeros(length(eelements),length(tb)-1); -for i=1:length(eelements) - [Sbin_assembly(i,:),~] = histcounts(spM(eelements(i),:),tb); + +if nargin < 3 + Sbin_assembly=zeros(length(eelements),size(spM,2)); + for i=1:length(eelements) + Sbin_assembly(i,:) = spM(eelements(i),:); + end +else + Sbin_assembly=zeros(length(eelements),length(tb)-1); + for i=1:length(eelements) + [Sbin_assembly(i,:),~] = histcounts(spM(eelements(i),:),tb); + end end -maxrate=max(Sbin_assembly(:)); +maxrate=max(Sbin_assembly(:)); Sbin_p=cell(1,maxrate); for i=1:maxrate Sbin_p{i}=zeros(size(Sbin_assembly)); - Sbin_p{i}(find(Sbin_assembly>=i))=1; + Sbin_p{i}(Sbin_assembly>=i)=1; end activity=zeros(maxrate,size(Sbin_assembly,2)); @@ -138,16 +161,16 @@ if alag(e1)==0 Sbin_shifted(e1,:)=Sbin_p{i}(e1,:); elseif alag(e1)>0 - Sbin_shifted(e1,1:end-alag(e1))=Sbin_p{i}(e1,1+alag(e1):end); + Sbin_shifted(e1,1:end-alag(e1))=Sbin_p{i}(e1,1+alag(e1):end); else Sbin_shifted(e1,1-alag(e1):end)=Sbin_p{i}(e1,1:end+alag(e1)); end end - activity(i,:)=sum(Sbin_shifted,1); - activity(i,find(activity(i,:)==1))=0; + activity(i,:)=sum(Sbin_shifted,1); + activity(i,activity(i,:)==1)=0; end - -activationtot=sum(activity,1); + +activationtot=sum(activity,1); activationtot(isnan(activationtot))=0; activationtot=activationtot/length(eelements); end @@ -158,16 +181,23 @@ eelements=Assembly.elements; alag=Assembly.lag; -Sbin_assembly=zeros(length(eelements),length(tb)-1); -for i=1:length(eelements) - [Sbin_assembly(i,:),~] = histcounts(spM(eelements(i),:),tb); +if nargin < 3 + Sbin_assembly=zeros(length(eelements),size(spM,2)); + for i=1:length(eelements) + Sbin_assembly(i,:) = spM(eelements(i),:); + end +else + Sbin_assembly=zeros(length(eelements),length(tb)-1); + for i=1:length(eelements) + [Sbin_assembly(i,:),~] = histcounts(spM(eelements(i),:),tb); + end end -maxrate=max(Sbin_assembly(:)); +maxrate=max(Sbin_assembly(:)); Sbin_p=cell(1,maxrate); for i=1:maxrate Sbin_p{i}=zeros(size(Sbin_assembly)); - Sbin_p{i}(find(Sbin_assembly>=i))=1; + Sbin_p{i}(Sbin_assembly>=i)=1; end @@ -178,35 +208,20 @@ if alag(e1)==0 Sbin_shifted(e1,:)=Sbin_p{i}(e1,:); elseif alag(e1)>0 - Sbin_shifted(e1,1:end-alag(e1))=Sbin_p{i}(e1,1+alag(e1):end); + Sbin_shifted(e1,1:end-alag(e1))=Sbin_p{i}(e1,1+alag(e1):end); else Sbin_shifted(e1,1-alag(e1):end)=Sbin_p{i}(e1,1:end+alag(e1)); end end - activity(i,:)=sum(Sbin_shifted,1); + activity(i,:)=sum(Sbin_shifted,1); activity(i,:)= activity(i,:).^5; - activity(i,find(activity(i,:)==1))=0; + activity(i,activity(i,:)==1)=0; end - -activationtot=sum(activity,1); + +activationtot=sum(activity,1); activationtot(isnan(activationtot))=0; activationtot=activationtot/length(eelements)^5; end - - - - - - - - - - - - - - - diff --git a/Programs_and_data/Main_assemblies_detection.m b/Programs_and_data/Main_assemblies_detection.m index 31364e7..da1ea04 100644 --- a/Programs_and_data/Main_assemblies_detection.m +++ b/Programs_and_data/Main_assemblies_detection.m @@ -83,9 +83,9 @@ if ~isempty(assemblybin{gg}) assemblybin{gg}.bin_edges=tb; end - fprintf('%d - testing done\n',gg); - fname = sprintf('assembly%d.mat', gg); - parsave(fname,assemblybin{gg}) + % fprintf('%d - testing done\n',gg); + % fname = sprintf('assembly%d.mat', gg); + % parsave(fname,assemblybin{gg}) end end assembly.bin=assemblybin; From 0b007dd58fc6fc26916765694f195989c07d13a1 Mon Sep 17 00:00:00 2001 From: cowenowner Date: Sun, 29 Nov 2020 10:08:36 -0700 Subject: [PATCH 6/6] Update FindAssemblies_recursive.m --- Programs_and_data/FindAssemblies_recursive.m | 224 +++++++++---------- 1 file changed, 112 insertions(+), 112 deletions(-) diff --git a/Programs_and_data/FindAssemblies_recursive.m b/Programs_and_data/FindAssemblies_recursive.m index 3479070..423d04d 100644 --- a/Programs_and_data/FindAssemblies_recursive.m +++ b/Programs_and_data/FindAssemblies_recursive.m @@ -9,9 +9,9 @@ % for information please contact eleonora.russo@zi-mannheim.de; daniel.durstewitz@zi-mannheim.de. % % last update 11/01/2016 -%% zero order +%% zero order nu=size(binM,1); - + for w1=1:nu assembly_in.n{w1}.elements=w1; assembly_in.n{w1}.lag=[]; @@ -27,36 +27,36 @@ clear assembly_out ANfo -alpha=alph*2/(nu*(nu-1)*(2*maxlag+1)); %% bonferroni correction +alpha=alph*2/(nu*(nu-1)*(2*maxlag+1)); %% bonferroni correction n_as=size(ANin,1); nu=size(ANin,2); prANout=ones(nu,nu); % max possible dimension, then I cut -ANfo=zeros(nu,nu); +ANfo=zeros(nu,nu); assembly_out.n=cell(1,n_as*nu); nns=1; -for w1=1:(nu-1) +for w1=1:(nu-1) for w2=w1+1:nu spikeTrain2=binM(w2,:)'; n2=w2; - true=0; + true=0; [assemD]=TestPair(assembly_in.n{w1},spikeTrain2,n2,maxlag,Dc,reference_lag); - if assemD.pr(end)No_th - assembly_out.n{nns}=assemD; - prANout(w1,w2)=assemD.pr; - ANfo(w1,w2)=1; - true=1; + if assemD.pr(end)No_th + assembly_out.n{nns}=assemD; + prANout(w1,w2)=assemD.pr; + ANfo(w1,w2)=1; + true=1; end - - if true + + if true nns=nns+1; end clear spikeTrain2 - end - + end + end assembly_out.n(nns:end)=[]; @@ -72,8 +72,8 @@ if isempty(assembly.n) assembly_output=[]; end - - + + fname = sprintf('Assembly_O%d_b%d.mat', 1,gg); % save(fname,'-v6','assembly'); save(fname,'-v7.3','assembly'); @@ -88,14 +88,14 @@ % fname = sprintf('Assembly_O%d.mat', O); % assembly=load(fname); while Oincrement && O<(O_th-1) % this loop stop when no more units can be added to the existing assemblies or if we force a maximal assembly size (O_th:= max O_th+1 elements in the assembly ) - -% fname = sprintf('Assembly_O%d.mat', O); -% assembly=load(fname); + + % fname = sprintf('Assembly_O%d.mat', O); + % assembly=load(fname); Oincrement=0; - n_as=length(assembly.n); % number of groups previously found - assembly_out.n=cell(1,n_as*nu); % max possible dimension, then I cut - + n_as=length(assembly.n); % number of groups previously found + assembly_out.n=cell(1,n_as*nu); % max possible dimension, then I cut + nns=1; for w1=1:n_as % it runs over the existing assemblies w1_elements=assembly.n{w1}.elements; @@ -104,89 +104,89 @@ w2_to_test=unique(w2_to_test); alpha=alph/(length(w2_to_test)*n_as*(2*maxlag+1)); %% bonferroni correction only for the test that I actually perform - + for ww2=1:length(w2_to_test) - + w2=w2_to_test(ww2); spikeTrain2=binM(w2,:)'; -% n2=w2; + % n2=w2; true=0; - + [assemD]=TestPair(assembly.n{w1},spikeTrain2,w2,maxlag,Dc,reference_lag); - - if assemD.pr(end)No_th - assembly_out.n{nns}=assemD; - ANfo(w1,w2)=1; - true=1; - Oincrement=1; + + if assemD.pr(end)No_th + assembly_out.n{nns}=assemD; + ANfo(w1,w2)=1; + true=1; + Oincrement=1; end - - if true + + if true nns=nns+1; end clear spikeTrain2 assemD - end + end end assembly_out.n(nns:end)=[]; - - if nns>1 - O=O+1; - assembly=assembly_out; - clear assembly_out -%%% pruning step 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% between two assemblies with the same unit set arranged into different configurations I choose the most significant one - - na=length(assembly.n); % number assemblies - nelement=O+1; % number elements for assembly - selection=nan(na,nelement+1+1); - assembly_final.n=cell(1,na); %max possible dimension - nns=1; - for i=1:na - elem=sort(assembly.n{i}.elements); - [ism,indx]=ismember(elem,selection(:,1:nelement),'rows'); - if ~ism - assembly_final.n{nns}=assembly.n{i}; - selection(nns,1:nelement)=elem; - selection(nns,nelement+1)=assembly.n{i}.pr(end); - selection(nns,nelement+2)=i; - nns=nns+1; - else - if selection(indx,nelement+1)>assembly.n{i}.pr(end) - assembly_final.n{indx}=assembly.n{i}; - selection(indx,nelement+1)=assembly.n{i}.pr(end); - selection(indx,nelement+2)=i; + if nns>1 + O=O+1; + assembly=assembly_out; + clear assembly_out + + %%% pruning step 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % between two assemblies with the same unit set arranged into different configurations I choose the most significant one + + na=length(assembly.n); % number assemblies + nelement=O+1; % number elements for assembly + selection=nan(na,nelement+1+1); + assembly_final.n=cell(1,na); %max possible dimension + nns=1; + for i=1:na + elem=sort(assembly.n{i}.elements); + [ism,indx]=ismember(elem,selection(:,1:nelement),'rows'); + if ~ism + assembly_final.n{nns}=assembly.n{i}; + selection(nns,1:nelement)=elem; + selection(nns,nelement+1)=assembly.n{i}.pr(end); + selection(nns,nelement+2)=i; + nns=nns+1; + else + if selection(indx,nelement+1)>assembly.n{i}.pr(end) + assembly_final.n{indx}=assembly.n{i}; + selection(indx,nelement+1)=assembly.n{i}.pr(end); + selection(indx,nelement+2)=i; + end end end + assembly_final.n(nns:end)=[]; + assembly=assembly_final; + clear assembly_final + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + end - assembly_final.n(nns:end)=[]; - assembly=assembly_final; - clear assembly_final + clear assemS2 assemS1 + -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - + fname = sprintf('Assembly_O%d_b%d.mat', O,gg); + save(fname,'-v6','assembly'); + % fname = sprintf('Assembly_2O%d.mat', O); + % save(fname,'-v6','assembly'); + bytesize=whos('assembly'); + if bytesize.bytes>bytelimit + fprintf('The algorithm has been interrupted because assembly structures reached a global size of %d bytes, this limit can be changed in size or removed with the ''bytelimit'' option\n',bytelimit); + O=O_th; end -clear assemS2 assemS1 - - -fname = sprintf('Assembly_O%d_b%d.mat', O,gg); -save(fname,'-v6','assembly'); -% fname = sprintf('Assembly_2O%d.mat', O); -% save(fname,'-v6','assembly'); - -bytesize=whos('assembly'); -if bytesize.bytes>bytelimit - fprintf('The algorithm has been interrupted because assembly structures reached a global size of %d bytes, this limit can be changed in size or removed with the ''bytelimit'' option\n',bytelimit); - O=O_th; end -end - + maxOrder=O; - + %% % pruning step 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% I remove assemblies whom elements are already ALL included in a bigger @@ -195,55 +195,55 @@ nns=1; for o=1:maxOrder - - fname = sprintf('Assembly_O%d_b%d.mat', o,gg); - load(fname); - minor=assembly; + + fname = sprintf('Assembly_O%d_b%d.mat', o,gg); + load(fname); + minor=assembly; clear assembly -% minor=assembly.order{o}; + % minor=assembly.order{o}; no=length(minor.n); % number assemblies selection_o=ones(no,1); - + for O=maxOrder:-1:o+1 fname = sprintf('Assembly_O%d_b%d.mat', O,gg); load(fname); major=assembly; - clear assembly -% major=assembly.order{O}; + clear assembly + % major=assembly.order{O}; nO=length(major.n); % number assemblies index_elemo=find(selection_o==1); - for i=1:sum(selection_o) + for i=1:sum(selection_o) elemo=minor.n{index_elemo(i)}.elements; - + for j=1:nO elemO=major.n{j}.elements; if ismember(elemo, elemO) - selection_o(index_elemo(i),end)=0; - j=nO; - end + selection_o(index_elemo(i),end)=0; + j=nO; + end end end if ~selection_o - O=0; + O=0; end - + end -%%% give as output just the maximal groups - -index_elemo=find(selection_o==1); -for i=1:sum(selection_o) -assembly_output.n{nns}=minor.n{index_elemo(i)}; -nns=nns+1; -end - - -recycle('off') -fname = sprintf('Assembly_O%d_b%d.mat', o,gg); -delete(fname) -recycle('on') + %%% give as output just the maximal groups + + index_elemo=find(selection_o==1); + for i=1:sum(selection_o) + assembly_output.n{nns}=minor.n{index_elemo(i)}; + nns=nns+1; + end + + + recycle('off') + fname = sprintf('Assembly_O%d_b%d.mat', o,gg); + delete(fname) + recycle('on') end @@ -251,7 +251,7 @@ end - +