% Jon Newman
% Create time series summary figure

%% Load SS
fid = '';
spk = loadspk([fid '.spk']);
ss = SqueakySpk(fid,spk.meta.fs_Hz,1,spk);


% sort_bound = [0 3*3600]; % This time bound indicates how sorted units will be organized with a rastergram

% Parameters
dt = 0.1; % Seconds, defines spike time binning granularity
data_bound = [0 300]; %seconds to to recording start time, this is the window that the cross correlation function will be calculated from
maxlag = 0.5; % seconds, The maximum lag to calculate the xcorr for
nunit = 10; % number to calculate the xcorr for



%% First Raster zoom in
close all
figure('units','centimeter','outerposition',[5 5 10 10]); % fill A4 paper width

    

    l = ceil(maxlag/dt);
    t = data_bound(1):dt:data_bound(2);
    XC =  zeros(1,2*l+1);
    
    num_units = max(SS1.SS.unit);
    r =  randperm(num_units);
    u = 1:MaxUnit;
    combos = nchoosek(u(1:nunit),2);
    XC =  zeros(length(combos),2*l+1);
    for j=1:length(combos)
        b1 = SS1.SS.time > rb(i,1) & SS1.SS.time <= rb(i,2) & SS1.SS.unit == combos(j,1);
        b2 = SS1.SS.time > rb(i,1) & SS1.SS.time <= rb(i,2) & SS1.SS.unit == combos(j,2);
        S1t = hist(SS1.SS.time(b1)-rb(i,1),t);
        S2t = hist(SS1.SS.time(b2)-rb(i,1),t);
        if sum(S1t) > 0 && sum(S2t) > 0
            [XC(j,:) lags] = xcorr(S1t,S2t,l);
            XC(j,:) = XC(j,:)/(sqrt(var(S1t)*var(S2t)));
        end
    end
    
    uXC = mean(XC);
    plot(tau*lags,tau*uXC);
    ylim([-10 100])
    
    

    
    ip = get(gca, 'Position');
    op = get(gca, 'OuterPosition');
    set(gca, 'Position', [op(1) ip(2) op(3) ip(4)]);
    if i == 1
        w = op(3);
        h = ip(4);
    else
        set(gca, 'Position', [op(1) ip(2) w h]);
    end
end
%%
% export_fig '2012-01-12_lollipop-10803_DW-MEStim-combined_pearsons_xcorr_FINAL' -pdf -transparent









