From: <?R???�ö???p?~?u?~?? Microsoft Internet Explorer 5>
Date: Wed, 8 Jul 2009 03:21:26 +0300

% 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;