% Jon Newman
% Examine CL time series quickly

fid = '2012_03_03_Ollie-363_Long-Term-optoBB_00';

% Load simoc file
sim = loadsimoc([fid '.simoc']);

% % Load digital file
% d = loadsimoc([fid '.dig']);

blksec = 3600;
[tend idx] = loadspk([fid '.spk'],'last');
tnow = 0; % current position in file.

while (tnow < tend)


    tnow = tnow + blksec;
    
end

%spk = loadspk([fid '.spk']);
spk = loadspk([fid '.spk'],[0 tend],0);

% For getting average perfomance
t_ignore = 0; %sec
t_pre = 0;
%t_end = (6 + 12)*3600;
t_end = (3 + 12 + 9)*3600;

% Defining plot limits
frmax = 100;
targmax = 1;

% Plot downsample factor
dsfac = 10;

% CTL sig
ctlmax = 5;

% step overlay plots
pre_step = 2*3600; % sec
post_step = 20; % sec

% UNFR dt
dtUNFR = 60; % sec

% spnum
spnum = 5; % number of subplots

%% Process data

% ss obj
ss = SqueakySpk(fid,spk.meta.fs_Hz,1,spk);
ss.ASDR(dtUNFR,[],[],[],[],0);
close all;

% Parse out the step section of the protocol
stepobs = sim.obs(sim.target > 0 & sim.time > t_pre);
steptime = sim.time(sim.target > 0 & sim.time > t_pre);
steps = sim.target(sim.target > 0 & sim.time > t_pre);
usteps = unique(steps);

% Get avg firing rates during steps
ufr = [];

% Save the step sequences for overlay plot
stime = {};
ref = {};
ufr_ts = {};

for i = 1:length(usteps)

%   Find the mean UNFR for the current step
    thistime = steptime(steps == usteps(i));
    thismean = mean(stepobs(steps == usteps(i) & steptime > thistime(1) + t_ignore));
    ufr = [ufr thismean];

%     Put the step snippet into the time series matrix
    stime{i} = sim.time(sim.time > (thistime(1) - pre_step) & sim.time < (thistime(end) + post_step)) - thistime(1);
    ref{i} = sim.target(sim.time > (thistime(1) - pre_step) & sim.time < (thistime(end) + post_step));
    ufr_ts{i} = sim.filt(sim.time > (thistime(1) - pre_step) & sim.time < (thistime(end) + post_step));

end

% bin the control signal the same as your asdr
[u_u_on b_u] = hist(sim.time(sim.output(1,:) > 0),ss.asdr.bin);
[u_u_off b_u] = hist(sim.time,ss.asdr.bin);
ctl_effort_dTsec = u_u_on./u_u_off;

%% Full time series

% Block time
blk_time = t_end/3600; % hours
t1 = 0; % hours

% Time slices of interest (hrs)
tpre = [t1 t1 + blk_time];
snip = 3600*tpre;

close all
figure()
for i = 1:size(snip,1)
    
    b = snip(i,:);
    sim_time = (sim.time(sim.time >=  b(1) & sim.time <=  b(2)) - b(1))/3600; 
    asdr_bin = (ss.asdr.bin(ss.asdr.bin >=  b(1) & ss.asdr.bin <=  b(2))- b(1))/3600;
    xlimit = (b-b(1))/3600;
    
    
    % ASDR, 10 msec bins
    subplot(spnum,3,[1 2 3]);
    plot(sim_time,sim.filt(sim.time >=  b(1) & sim.time <=  b(2)),'k')
    xlim(xlimit)
    set(gca,'Xtick',[]);
    if (i == 1)
        ylabel('UNFR_{10msec}(t)','FontSize',12)
    end
    
    % Rastergram
    subplot(spnum,3,[4 5 6]);
    hold on
    ss.RandScat(b,[]);
    set(gca,'Xtick',[]);
    xlabel('')
    ylabel('Unit #','FontSize',12)
    xlim(xlimit*3600)
    if (i ~= 1)
        ylabel('')
    end
    
    % Control performance, dT bins
    subplot(spnum,3,[10 11 12]);
    unfr = mean(ss.usdr.usdr(:,ss.usdr.units ~= 0),2);
    st = sim.target(sim.time >=  b(1) & sim.time <=  b(2));
    %stairs(sim_time,st,'r','LineWidth',2)
    plot([0 3 3 15 15 24],[0 0 7 7 0 0],'r','LineWidth',2)
    hold on
    %stairs(ss.usdr.bin(ss.usdr.bin >=  b(1) & ss.usdr.bin <=  b(2)),unfr(ss.usdr.bin >=  b(1) & ss.usdr.bin <=  b(2)))
    stairs(ss.usdr.bin/3600,unfr)
    
    xlim(xlimit)
    set(gca,'Xtick',[]);
    if (i == 1)
        legend('ref(t)','UNFR_{60sec}(t)','Location','NorthWest')
        ylabel('Tracking Performance','FontSize',12)
    end
    
    % Control effort, dT bins
    subplot(spnum,3,[13 14 15]);
    stairs(b_u(b_u >=  b(1) & b_u <=  b(2))/3600,100*ctl_effort_dTsec(b_u >=  b(1) & b_u <=  b(2)),'k-')
    ylim([0 100])
    xlim(xlimit)
    ylabel('% effort','FontSize',12)
    xlabel('Time (hrs)','FontSize',12)
    
    
end

%% Zoomed rastergram of pre, during, post control

% Block time
blk_time = 1200; % seconds
t1 = 1.5 - (blk_time/2/3600); % hours
t2 = 9 - (blk_time/2/3600); % hours
t3 = 16.5 - (blk_time/2/3600); % hours

% Time slices of interest (hrs)
tpre = [t1 t1 + blk_time/3600];
tdur = [t2 t2 + blk_time/3600];
tpost = [t3 t3 + blk_time/3600];

%
snip = 3600*[tpre;tdur;tpost];

for i = 1:size(snip,1)
    
    b = snip(i,:);
%     sim_time = (sim.time(sim.time >=  b(1) & sim.time <=  b(2)) - b(1))/3600; 
%     asdr_bin = (ss.asdr.bin(ss.asdr.bin >=  b(1) & ss.asdr.bin <=  b(2))- b(1))/3600;
    xlimit = (b-b(1));
    
    
    subplot(spnum,3,[4 5 6])
    hold on
    plot([b(1) b(2)],[0 0],'r*');
    
    h = subplot(spnum,3,6+i);
    hold on
       ss.RandScat(b,[],0);
    set(gca,'Xtick',[],'Ytick',[]);
    xlabel('')
    xlim(xlimit)
    if (i == 1)
        ylabel('Unit #')
    else
        ylabel('')
    end
    
    % remove space
    p = get(h, 'pos');
    dp = 0.03;
    switch i
        case 1
            p(3) = p(3) + dp;
        case 2
            p(3) = p(3) + dp;
            p(1) = p(1) - dp/2;
        case 3
            p(3) = p(3) + dp;
            p(1) = p(1) - dp;
    end
    
    set(h, 'pos', p);
           
            
   
    
    
end

