% Jon Newman
% Examine CL time series quickly

fid = '2012_10_19_Herb-I5386_one-min-clamp_10';
pic = loadpic([fid '.pic']);
spk = loadspk([fid '.spk']);

% Plot parameters
fr_max = 15;
fr_min = -1;
sm_w = 20; %smoothing width

%% For getting average perfomance
t_ignore = 15; %sec
t_pre = 15;
t_end = spk.time(end);

% Plot downsample factor
dsfac = 1;

% CTL sig
ctlmax = 5;

% step overlay plots
pre_step = 15; % sec
post_step = 10; % sec

% UNFR dt
dtUNFR = 10; % sec


%% Process data
% Smoothed the output
% sf = smooth(sim.filt,100);

% 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 = pic.meas(pic.clon & pic.time > t_pre);
stepfilt = pic.filt(pic.clon & pic.time > t_pre);
steptime = pic.time(pic.clon & pic.time > t_pre);
stepctl1 = pic.out1(pic.clon & pic.time > t_pre);
stepctl2 = pic.out2(pic.clon & pic.time > t_pre);
steps = pic.targ(pic.clon & pic.time > t_pre);
usteps = unique(steps);

% Get avg firing rates during steps
ufr = [];
sfr = [];
uctl1 = [];
sctl1 = [];
uctl2 = [];
sctl2 = [];

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



for i = 1:length(usteps)
    
    % Find the mean UNFR for the current step
    thistime = steptime(steps == usteps(i));
    thismean = mean(stepfilt(steps == usteps(i) & steptime > thistime(1) + t_ignore));
    thiss = std(stepfilt(steps == usteps(i) & steptime > thistime(1) + t_ignore));
    ufr = [ufr thismean];
    sfr = [sfr thiss];
    
    thismean = mean(stepctl1(steps == usteps(i) & steptime > thistime(1) + t_ignore));
    thiss = std(stepctl1(steps == usteps(i) & steptime > thistime(1) + t_ignore));
    uctl1 = [uctl1 thismean];
    sctl1 = [sctl1 thiss];
    
    thismean = mean(stepctl2(steps == usteps(i) & steptime > thistime(1) + t_ignore));
    thiss = std(stepctl2(steps == usteps(i) & steptime > thistime(1) + t_ignore));
    uctl2 = [uctl2 thismean];
    sctl2 = [sctl2 thiss];
    
    % Put the step snippet into the time series matrix
    b = pic.time > (thistime(1) - pre_step) & pic.time < (thistime(end) + post_step);
    stime{i} = pic.time(b) - thistime(1);
    ref{i} = pic.targ(b);
    ufr_ts{i} = pic.filt(b);
    ctl1_ts{i} = pic.out1(b);
    ctl2_ts{i} = pic.out2(b);
    
end

unfr_avg = steps(1);


%% Just the Steps
close all
figure()

% Step overlay
subplot(4,3,[1 2 4 5])
hold on
col = hsv(length(stime));
for i = 1:length(stime)
    
    rt = ref{i};
    tt = stime{i};
%     tt = tt(rt ~=0);
%     rt = rt(rt ~=0);
    plot(tt,rt,'k')
    plot(stime{i},smooth(ufr_ts{i},sm_w),'Color',col(i,:))
end
st = stime{1};
r = ref{1};
xlim([-pre_step ceil(max(stime{1}))])
ylim([fr_min fr_max+4])
ylabel('UNFR (Hz/Unit)')
set(gca,'XTick',[],'TickDir','out')

% Control Signal
subplot(4,3,[7 8])
hold on
col = hsv(length(stime));
for i = 1:length(stime)
    plot(stime{i},smooth(ctl1_ts{i},sm_w),'Color',col(i,:))
end
st = stime{1};
r = ref{1};
xlim([-pre_step ceil(max(stime{1}))])
ylim([-0.1 1.1])
ylabel('UNFR (Hz/Unit)')
set(gca,'XTick',[],'TickDir','out')

% Control Signal
subplot(4,3,[10 11])
hold on
col = hsv(length(stime));
for i = 1:length(stime)
    plot(stime{i},smooth(ctl2_ts{i},1),'Color',col(i,:))
end
st = stime{1};
r = ref{1};
xlim([-pre_step ceil(max(stime{1}))])
ylim([-0.1 1.1])
ylabel('UNFR (Hz/Unit)')
xlabel('Time (sec)')
set(gca,'TickDir','out')

% Average performance
subplot(4,3,[3 6])
plot(usteps,usteps,'r-s')
hold on
errorbar(usteps,ufr,sfr,'k-o')

axis([fr_min fr_max fr_min fr_max])
ylabel('Hz/Unit')
xlabel('Hz/Unit')
set(gca,'XTick',[],'TickDir','out')


subplot(4,3,[9 12])
hold on
errorbar(usteps,uctl1,sctl1,'b-o')
errorbar(usteps,uctl2,sctl2,'y-o')

axis([fr_min fr_max -0.1 1.1 ])
ylabel('<U(t)>')
xlabel('Hz/Unit')
set(gca,'TickDir','out')

%%
export_fig(fid,'-pdf','-transparent',gcf)

