% Jon Newman
% Create time series summary figure

% First, produce the correlation figure for the BB controller
function generate_xcorr(SS,bound,dt,maxlag,units,useGPU)

if nargin < 2 || isempty(bound)
    bound = [min(SS.time) max(SS.time)];
end

if nargin < 3 || isempty(bound)
    dt = 0.01; % seconds
end

if nargin < 4 || isempty(bound)
    maxlag = 0.25; % seconds
end

if nargin < 5 || isempty(bound)
    units = unique(SS.unit);
end

if nargin < 6 || isempty(useGPU)
    useGPU = 0;
end


l = ceil(maxlag/dt); % samples in xcorr
t = 0:dt:diff(bound); % normalized time
combos = nchoosek(units,2); % unit pairs
XC =  zeros(length(combos),2*l+1);

for j=1:length(combos)
    
    % Get spike times for units in combo(j) that occured within the time bound
    b1 = SS.time > bound(1) & SS.time <= bound(2) & SS.unit == combos(j,1);
    b2 = SS.time > bound(1) & SS.time <= bound(2) & SS.unit == combos(j,2);
    S1t = hist(SS.time(b1)-bound(1),t);
    S2t = hist(SS.time(b2)-bound(1),t);
    
    if sum(S1t) > 0 && sum(S2t) > 0
        
        if ~useGPU
            [XC(j,:), lags] = xcorr(S1t/dt,S2t/dt,l,'unbiased');
        else
            X1 = gpuArray(S1t/dt);
            X2 = gpuArray(S2t/dt);
            [xc, lags] = xcorr(X1,X2,l,'unbiased');
            XC(j,:) = gather(xc);
        end
    end
end
XC(XC<0) = 0;
uXC = nanmean(XC);


% Plot the results
plot(dt*lags,sqrt(uXC));
xlim([-maxlag maxlag])
ylabel('Firing rate')
xlabel('Lag (sec)')

end









