% Finds polychronous groups in the workspace generated by spnet.m
% and saved in the matlab file fname, e.g., fname = '3600.mat';
% Created by Eugene M.Izhikevich.           November 21, 2005
% Modified on April 2, 2008 based on suggestions of Petra Vertes (UK).

% Main idea: for each mother neuron, consider various combinations of 
% pre-synatic (anchor) neurons and see whether any activity of a silent
% network could emerge if these anchors are fired. 

fname = '18000.mat';

global a d N D pp s ppre dpre post pre delay T

load(fname);

%uncomment the command below to test the algorithm for shuffled (randomized e->e) synapses
%e2e = find(s>=0 & post<Ne); s(e2e) = s(e2e(randperm(length(e2e))));
            
groups={};          % the list of all polychronous groups
anchor_width=3;     % the number of anchor neurons, from which a group starts
min_group_path=7;   % discard all groups having shorter paths from the anchor neurons 
T=150;              % the max length of a group to be considered;
                    % longer groups will be cut at t=T

% Make necessary initializations to speed-up simulations.

% This matrix provides the delay values for each synapse.
for i=1:N
    for j=1:D
        delay(i, delays{i,j})=j;
    end;
end;

%This cell element tells what the presynaptic delay is; 
for i=1:N
    dpre{i}=delay( pre{i} );
end;

%This cell element tells where to put PSPs in the matrix I (N by 1000)
for i=1:N
    pp{i}=post(i,:)+N*(delay(i,:)-1);
end;

%This cell element tells who (indexes) the pre-synaptic neurons are; 
for i=1:N
    ppre{i}=mod( pre{i}, N);
end;

sm_threshold = 0.95*sm;     % discard all weak exc->exc synapses
s(find(post<Ne & s>0 & s<=sm_threshold))=0;

for i=1:Ne
    
    anchors=1:anchor_width;                     % initial choice of anchor neurons
    strong_pre=find(s(pre{i})>sm_threshold);    % candidates for anchor neurons
    if length(strong_pre) >= anchor_width       % must be enough candidates
     while 1        % will get out of the loop via the 'break' command below
         
        gr=polygroup( ppre{i}(strong_pre(anchors)), D-dpre{i}(strong_pre(anchors)) ); 
        
        %Calculate the longest path from the first to the last spike
        fired_path=sparse(N,1);        % the path length of the firing (from the anchor neurons)       
        for j=1:length(gr.gr(:,2))
            fired_path( gr.gr(j,4) ) = max( fired_path( gr.gr(j,4) ), 1+fired_path( gr.gr(j,2) ));
        end;
        
        longest_path = max(fired_path);
        
        if longest_path>=min_group_path 
            
            gr.longest_path = longest_path(1,1); % the path is a cell
            
            % How many times were the spikes from the anchor neurons useful?
            % (sometimes an anchor neuron does not participate in any
            % firing, because the mother neuron does its job; such groups
            % should be excluded. They are found when the mother neuron is
            % an anchor neuron for some other neuron).
            useful = zeros(1,anchor_width);
            anch = ppre{i}(strong_pre(anchors));
            for j=1:anchor_width
                useful(j) = length( find(gr.gr(:,2) == anch(j) ) );
            end;
       
            
            
            if all(useful>=2)
                
                groups{end+1}=gr;           % add found group to the list
                disp([num2str(round(100*i/Ne)) '%: groups=' num2str(length(groups)) ', size=' num2str(size(gr.firings,1)) ', path_length=' num2str(gr.longest_path)])   % display of the current status
            
                plot(gr.firings(:,1),gr.firings(:,2),'o');
                hold on;
                for j=1:size(gr.gr,1)
                    plot(gr.gr(j,[1 3 5]),gr.gr(j,[2 4 4]),'.-');
                end;
                axis([0 T 0 N]);
                hold off
                drawnow;
            end;
        end

        % Now, get a different combination of the anchor neurons
        k=anchor_width;
        while k>0 & anchors(k)==length(strong_pre)-(anchor_width-k)
            k=k-1;
        end;
        
        if k==0, break, end;    % exhausted all possibilities
        
        anchors(k)=anchors(k)+1;
        for j=k+1:anchor_width
            anchors(j)=anchors(j-1)+1;
        end;
        
        pause(0); % to avoid feezing when no groups are found for long time
        
     end;
    end;
end;