% Jon Newman

fid = '2013_12_22_Smiley-24729_01';
load('2013_12_22_Smiley-24729_01_parse.SS','-mat');
load('2013_12_22_Smiley-24729_01_oc_parse.mat');
oc.time = oc.time - t_pre;

%%
% For the correlation and sync measures, only use the most reliable unit on
% each electrode.
SScop = SS.Copy;
uc = unique(SScop.channel);
bu = [0];
for i = 1:length(uc)
    uu = unique(SScop.unit(SS.channel == uc(i) & SScop.unit~=0));
    c = [];
    for j = 1:length(uu)
        c = [c; numel(SScop.unit(SS.unit == uu(j)))];
    end
    
    [~, mi] = max(c);
    bu = [bu, setdiff(uu, uu(mi))'];
    
end
SScop.RemoveUnit(bu);
uu = unique(SScop.unit);
idx = 1:length(uu);

% XCorr parameters
dt = 0.01;
stat_window = [-150 150];
stat_time = [0.5 2.5 5.5]*3600; 

maxlag = 0.5;
CORR.dt = dt;
CORR.tcorr = diff(stat_window);
CORR.maxlag = maxlag;

XC = nan(numel(stat_time),length(uu),length(uu),numel(-maxlag:dt:maxlag)-1);
combos = nchoosek(uu,2);

 for i = 1:numel(stat_time)
    
        
    bound = [stat_time(i) stat_time(i)] + stat_window;
    T = SScop.time(SScop.time > bound(1) & SScop.time < bound(2));
    U = SScop.unit(SScop.time > bound(1) & SScop.time < bound(2));
    
    ncomb = size(combos,1);
    for j = 1:ncomb
        
        j1 = idx(uu == combos(j,1));
        j2 = idx(uu == combos(j,2));
        
        T1 = T(U == combos(j,1));
        T2 = T(U == combos(j,2));
        
        if ~isempty(T1) && ~isempty(T2)
            [XC(i,j1,j2,:), bins] = spk_xcorr(T1,T2,dt,maxlag,0);
        end
        
        disp([num2str(j) ' of ' num2str(ncomb)])
    end
end

CORR.XC = XC;
CORR.bins = bins;
save('2013_12_22_Smiley-24729_01_CORR','CORR')

%% Mean firing levels
bounds = [[0 1];[1 4];[4 7]]*3600;
bounds_oc = bounds - 3600;
FR = zeros(size(bounds,1),1);
uFR = zeros(size(bounds,1),1);
sFR = zeros(size(bounds,1),1);

for i = 1:size(bounds,1)

    bound = bounds(i,:);
    T = SS.time(SS.time > bound(1) & SS.time < bound(2) & SS.unit~=0);
    FR(i) = numel(T)/max(uu)/diff(bound);
    
    bound = bounds_oc(i,:);
    f = oc.filt(oc.time >= bound(1) & oc.time <= bound(2));
    
   uFR(i) =  mean(f);
   sFR(i) =  std(f);

end


% Figures
close all
ah = 3;
aw = 3;
as = 1;

f = figure('Units','centimeters','Position',[2 2 20 20]);

axes('Units','centimeters','Position',[2 2+as+aw aw ah]);
hold on
% errorbar([1 2 3],uFR,sFR,'o')
bar(uFR)
plot([0 4],[0.75 0.75],'r-')
set(gca,'TickDir','out','Xlim',[0 4],'Ylim',[0 1.5])

axes('Units','centimeters','Position',[2 2 aw ah]);
hold all

col = [[0 0 1];
        [1 0 0];
        [0 1 0]];
for i = 1:size(CORR.XC,1)
    
    bins = CORR.bins;
    xc = squeeze(CORR.XC(i,:,:,:));
    xc = squeeze(reshape(xc,[],size(xc,2)^2,size(xc,3)));
    plot(bins,nanmean(xc),'Color',col(i,:))
end
set(gca,'TickDir','out')

%%
saveas(gcf, fid)
export_fig([fid '_stats'],'-pdf','-transparent',gcf)


