function [s,SS] = PatternJitter(N,t,L,r,Lmethod,Gmethod,tfix,t0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% function [s,SamplerStructure] = PatternJitter(N,t,L,r,Lmethod,Gmethod,tfix,t0)
% function s = PatternJitter(N,SamplerStructure,t0)
%
% DESCRIPTION OF JITTER:
%
% Segments the (sorted) spike train into maximal r-patterns.  An r-pattern
% is a sequence of consecutive spikes such that any two consecutive spikes 
% in the sequence are separated by (temporal) distance r or less.  A 
% maximal r-pattern of a spike train t is an r-pattern in t that is not 
% imbedded in a larger r-pattern in t.  That is, the beginning and end of 
% the pattern are separated by more than r from the spikes preceeding and 
% following the pattern, respectively.  A single spike is a maximal 
% r-pattern if there are no other spikes within distance r of it.
%
% "Independently and uniformly" jitters the maximal r-patterns in t over
% windows of length L up to the hard constraint that the set of maximal 
% r-patterns is not changed by jittering.  The jittered spike train is
% stored in s(:,1).  This is repeated N times with the kth jittered spike
% train stored in s(:,k).  Jitter windows for a pattern are defined by 
% the time of its first spike.
%
% Note: r = -inf corresponds to SpikeJitter.  The ouput may not be sorted.
%
% DESCRIPTION OF INPUT:
%
% N is the the number of MonteCarlo trials.
% t is a sequence of increasing (non-decreasing) spike times.
% L is a non-negative scalar (duration of jitter window).
% r is a non-negative scalar (pattern size) or -inf.
% Lmethod is one of {'centered','fixed'}, 'fixed' is default
%   'centered': the jitter window for t(k) is t(k) + [-L/2,L/2].
%   'fixed': the jitter window for t(k) is floor(t(k)/L)*L + [0,L].
% Gmethod is one of {'real','integer'}, 'integer' is default
%   'real': no changes from above
%   'integer': t and jittered spikes are constrained to the integers
%       L is assumed to be an integer >= 1
%       r is assumed to be an integer
%       the jitter positions are
%       'centered': t(k)-floor(L/2) + {0,1,...,L-1}
%       'fixed': floor(t(k)/L)*L + {0,1,...,L-1}
%
% OPTIONAL ARGUMENTS:
%
% tfix is a list of indices into t (or a logical array the same size as t).  
% The spike times t(tfix) are not allowed to jitter.  This can also cause 
% entire patterns to remain fixed.  Use an index of inf to get the number of
% spikes in t.
%
% t0 is the spike time preceding the first pattern in t. 
% It can be a scalar or a vector of length N (so that it is different for
% each sample).  The initial spike time in the jittered spike train will
% be constrained to be greater than t0+r.  t0 defaults to -inf so that
% this constraint has no effect.
%
% The optional ouput argument SamplerStructure contains a collection of
% precomputed data structures for multiple calls to PatternJitter using
% the same input arguments (although N and t0 are allowed to vary).
% SamplerStructure(1:k) can be used to sample from (the marginal
% distribution of) the first k patterns.  
% [Note: SamplerStructure(j:k) can also be used to sample from patterns j
% to k, but this will only be the marginal distribution if the spike times 
% in t0 correspond to samples from the marginal distribution of the first 
% spike in pattern j-1.]
%
% REMARK: It is much more efficient to call PatternJitter(N,...) than
% to call PatternJitter(1,...) N different times, even (and especially) if 
% the SamplerStructure has been precomputed.  The only exception is when
% N (and the number of spikes) is so large that memory becomes a problem.
%
%
% Uses: 
% patjitPattern.m
% patjitSample.m (patjitRNG.m, patjitRNGint.m, mypolyval.m, polyroot.m)
% patjitDP.m (mypolyval.m) 
% patjitDPint.m
%
% Tested on Matlab 7.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Matt Harrison 
% 4/03/2006

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% input checking
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if numel(N) ~= 1 || N ~= round(N) || N < 0 || isinf(N)
    error('bad N')
end
if isstruct(t)
    if nargin > 3
        error('too many inputs for structure usage')
    end
    if nargin == 2
        L = -inf;
    end
    s = patjitSample(N,t,L);
    if nargout > 1
        SS = t;
    end
    return
end
if numel(L) ~= 1 || L < 0
    error('bad L')
end
if numel(r) ~= 1 
    error('bad r')
end

if nargin < 8 || isempty(t0), t0 = -inf; end
if numel(t0) ~= 1 && numel(t0) ~= N
    error('t0 is incorrectly sized')
end

if nargin < 7, tfix = []; end
if ~islogical(tfix)
    tfix(isinf(tfix))=numel(t);
    tfixtmp = tfix;
    tfix = false(size(t));
    if ~isempty(tfix)
        tfix(tfixtmp) = true;
    end
elseif numel(tfix) ~= numel(t)
    error('tfix is incorrectly sized')
end
if L == 0
    tfix(:) = true;
end

t = t(:);
tfix = tfix(:);
sortflag = false;
if ~issorted(t)
    [t,sortndx] = sort(t); 
    tfix = tfix(sortndx);
    sortflag = true;
    if nargout > 1
        warning('The output SamplerStructure will not unsort the spike train.')
    end
end
fixany = any(tfix);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PATTERN EXTRACTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[pat,patstart,patn,patlen,tndx] = patjitPattern(t,r);
n = numel(pat);

% get the appropriate r's
patr = [r; patlen(1:n-1)+r];

% find any patterns to fix 
patfix = false(n,1);
if fixany
    patfix(tndx(tfix)) = true;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine the methods and 
% compute the appropriate intervals.
% More input checking.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 6 || isempty(Gmethod), Gmethod = 'integer'; end
if nargin < 5 || isempty(Lmethod), Lmethod = 'fixed'; end

switch lower(Lmethod)

    case 'centered'

        switch lower(Gmethod)

            case 'real'
                Greal = true;
                u = patstart - L/2;
            case 'integer'
                Greal = false;
                if L ~= round(L) || L < 1, error('bad L'), end
                if any(t ~= round(t)) || any(t0 ~= round(t0)), error('non integer t'), end
                if r ~= round(r), error('bad r'), end
                u = round(patstart - floor(L/2));
            otherwise
                error('unknown Gmethod')
        end
        
    case 'fixed'

        switch lower(Gmethod)
            
            case 'real'
                Greal = true;
                if L > 0
                    u = floor(patstart/L)*L;
                else
                    u = patstart;
                end
            case 'integer'
                Greal = false;
                if L ~= round(L) || L < 1, error('bad L'), end
                if any(t ~= round(t)) || any(t0 ~= round(t0)), error('non integer t'), end
                if r ~= round(r), error('bad r'), end
                u = round(floor(patstart/L)*L);
            otherwise
                error('unknown Gmethod')
        end
        
    otherwise
        error('unknown Lmethod')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CREATE THE SAMPLER STRUCTURE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if Greal
    SS = struct('r',cell(n,1),'alpha',cell(n,1),'c',cell(n,1),...
        'x',cell(n,1),'d',cell(n,1),'m',cell(n,1),'n',cell(n,1),...
        'beta',cell(n,1),'pat',cell(n,1),'patn',cell(n,1),'Greal',cell(n,1));
else
    SS = struct('r',cell(n,1),'alpha',cell(n,1),'n',cell(n,1),...
        'beta',cell(n,1),'pat',cell(n,1),'patn',cell(n,1),'Greal',cell(n,1));
end

for k = 1:n
    SS(k).r = patr(k);
    SS(k).pat = pat{k};
    SS(k).patn = patn(k);
    SS(k).Greal = Greal;
end

if n > 0
    if Greal
        % initialize
        if patfix(n)
            SS(n).alpha = patstart(n);
            SS(n).n = 0;
        else
            [SS(n).alpha,SS(n).c,SS(n).x,SS(n).d,SS(n).m,SS(n).n,SS(n).beta] ...
                = patjitDP(u(n),u(n)+L,-inf,inf);
        end

        % dynamic programming to compute the rest of the structure
        for k = n-1:-1:1
            if patfix(k)
                SS(k).alpha = patstart(k);
                SS(k).n = 0;
            else
                [SS(k).alpha,SS(k).c,SS(k).x,...
                    SS(k).d,SS(k).m,SS(k).n,...
                    SS(k).beta] = patjitDP(u(k),u(k)+L,SS(k+1));
            end
        end

    else

        % initialize
        if patfix(n)
            SS(n).alpha = patstart(n);
            SS(n).n = 1;
            SS(n).beta = 1;
        else
            [SS(n).alpha,SS(n).n,SS(n).beta] ...
                = patjitDPint(u(n),u(n)+L-1,-inf,inf,1,1);
        end

        % dynamic programming to compute the rest of the structure
        for k = n-1:-1:1
            if patfix(k)
                SS(k).alpha = patstart(k);
                SS(k).n = 1;
                SS(k).beta = 1;
            else
                [SS(k).alpha,SS(k).n,SS(k).beta] ...
                    = patjitDPint(u(k),u(k)+L-1,SS(k+1));
            end
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SAMPLE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s = patjitSample(N,SS,t0);

% unsort if necessary
if sortflag && N > 0
    s(sortndx,:) = s;
end
