% TRYME.m
%
% This archive contains Matlab code for pattern jitter.
%
% This mfile demonstrates the use of PatternJitter.
%
% Typing TRYME at the Matlab prompt should create/modify
% figures 1 and 2 and display rasters of an original spike train
% (top horizontal row of dots in each figure) and several 
% independent pattern jittered versions (subsequent horizontal
% rows of dots).
%
% Figure 1 uses discretized pattern jitter (discretized to 1 ms
% bins).  This is the recommended usage.  
%
% Figure 2 uses continuous pattern jitter.  Continuous pattern
% jitter can be more efficient for very fine discretization, but 
% it is also more prone to numerical instabilities.  Numerical 
% instabilities will usually trigger an error, but there are no 
% guarantees.
%
% 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.
%
% All uses should include a citation to
%
% MT Harrison & S Geman. (2009) "A Rate and History-Preserving Resampling
%  Algorithm for Neural Spike Trains." Neural Computation 21(5): 1244-1258.
%

help TRYME

%---------------------------
% SET-UP
%---------------------------

% create a random spike train (one 1-second trial, 25 spikes)
t = sort(rand(25,1)); 

% set up the pattern jitter parameters
L = .05; % jitter length (50 ms)
R = .02; % history length (20 ms)
N = 15; % number of jitter resamples

%---------------------------
% DISCRETIZED PATTERN JITTER
%---------------------------

% discretize to 1 ms bins

t1000 = ceil(t*1000);
L1000 = round(L*1000);
R1000 = round(R*1000);

% get the pattern jitter resamples
% used fixed partition jitter, fix the times of the first and last spike

s1000 = PatternJitter(N,t1000,L1000,R1000,'fixed','integer',[1 inf]);

% plot the rasters
figure(1)
plot([t1000 s1000]'/1000,repmat([1.1*N ;(1:N)'],1,numel(t)),'.k')

%---------------------------
% CONTINUOUS PATTERN JITTER
%---------------------------

% get the pattern jitter resamples
% used fixed partition jitter, fix the times of the first and last spike

s = PatternJitter(N,t,L,R,'fixed','real',[1 inf]);

% plot the rasters
figure(2)
plot([t s]',repmat([N+1.5 ;(1:N)'],1,numel(t)),'.k')

