% % Jon Newman
% % Examine CL time series quickly
% 
% fids = {
%     '2012_10_11_Herb-I5387_one-min-clamp_01';
%     '2012_10_11_Herb-I5387_one-min-clamp_02';
%     };
% 
% % Plot parameters
% fr_max = 13;
% fr_min = -1;
% sm_w = 20; %smoothing width
% 
% %% For getting average perfomance
% t_ignore = 15; %sec
% t_pre = 15;
% 
% % 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
% 
% %% Storage
% % 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 = {};
% 
% %% Process data
% k = 1;
% psteps = [];
% for f = 1:length(fids)
%     
%     
%     fid = fids{f};
%     pic = loadpic([fid '.pic']);
%     spk = loadspk([fid '.spk']);
%     
%     t_end = spk.time(end);
%     
%     % 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);
%     
%     if f == 1
%         usteps = usteps(usteps ~= 0);
%     end
%     psteps =[psteps usteps];
%     
%     
%     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{k} = pic.time(b) - thistime(1);
%         ref{k} = pic.targ(b);
%         ufr_ts{k} = pic.filt(b);
%         ctl1_ts{k} = pic.out1(b);
%         ctl2_ts{k} = pic.out2(b);
%         
%         k = k + 1;
%     end
% end
% 
% 
% % Resort
% [psteps ind] = sort(psteps);
% ufr = ufr(ind);
% sfr = sfr(ind) ;
% uctl1 = uctl1(ind);
% sctl1 = sctl1(ind);
% uctl2 = uctl2(ind);
% sctl2 = sctl2(ind);
% 
% % Save the step sequences for overlay plot
% stime = stime(ind);
% ref = ref(ind);
% ufr_ts = ufr_ts(ind);
% ctl1_ts = ctl1_ts(ind);
% ctl2_ts = ctl2_ts(ind);

%% Just the Steps
close all
figure('units','centimeters','position',[10 10 10 10])

% Step overlay
axes('units','centimeters','position',[1 6 9 4])
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
axes('units','centimeters','position',[1 3.5 9 2])
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
axes('units','centimeters','position',[1 1 9 2])
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')

%%
export_fig('1-minute-time-series-summary','-pdf','-transparent',gcf)

