function [pat,patstart,patn,patlen,tndx] = patjitPattern(t,r);
% function [pat,patstart,patn,patlen,tndx] = patjitPattern(t,r)
%
% extracts the maximal r-patterns from the (sorted) spike train t
% (see PatternJitter.m)
%
% pat is a cell array of patterns (a pattern is a list of offsets from
% the initial spike time, pat{k} is column vector) 
%
% patstart(k) is the time of the first spike in pat{k} 
%
% patn(k) is the number of spikes in pat{k}
%
% patlen(k) = pat{k}(patn(k))
%
% tndx(j) = k such that spike t(j) is in pattern k

% Copyright (C) 2008 Matthew T Harrison
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    See <http://www.gnu.org/licenses/> for a copy of the license.

% sorting
t = t(:);
sortflag = false;
if ~issorted(t)
    [t,sortndx] = sort(t);
    sortflag = true;
end

n = numel(t);

% initialization
pat = cell(n,1);
patstart = zeros(n,1);
patn = zeros(n,1);
patlen = zeros(n,1);
tndx = zeros(n,1);

if n == 0
    return
end

% pattern extraction
kk = 1;
patk = 1;
tndx(1) = patk;
for k = 2:n
    if t(k)-t(k-1) <= r
        % this spike goes with the last pattern
        tndx(k) = patk;
        continue
    end
    % this spike starts a new pattern
    % remember the previous pattern
    pat{patk} = t(kk:k-1)-t(kk);
    patstart(patk) = t(kk);
    patn(patk) = k-kk; 
    patlen(patk) = t(k-1)-t(kk);
    % start a new pattern here
    kk = k;
    patk = patk + 1;
    tndx(k) = patk;
end
% get the final pattern
pat{patk} = t(kk:n) - t(kk);
patstart(patk) = t(kk);
patn(patk) = n+1-kk;
patlen(patk) = t(n)-t(kk);

% resizing
pat = pat(1:patk);
patstart = patstart(1:patk);
patn = patn(1:patk);
patlen = patlen(1:patk);

% unsorting
if sortflag
    tndx(sortndx) = tndx;
end
