% Jon Newman
% Examine CL time series quickly

fid = '2012_03_02_Ollie363_5minSteps_00';
sim = loadsimoc([fid '.simoc']);
%spk = loadspk([fid '.spk']);

%% For getting average perfomance
t_ignore = 15; %sec
t_pre = 60;
%t_end = spk.time(end);
t_end = loadspk([fid '.spk'],'last');


% Defining plot limits
frmax = 15;

% Plot downsample factor
dsfac = 1;

% CTL sig
ctlmax = 5;

% step overlay plots
pre_step = 30; % sec
post_step = 60; % sec

% UNFR dt
dtUNFR = 10; % sec

% first step ratio in the logarithm progression 
first_step_ratio = 1.0;

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

%% 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 = sim.obs(sim.target > 0 & sim.time > t_pre);
stepfilt = sim.filt(sim.target > 0 & sim.time > t_pre);
steptime = sim.time(sim.target > 0 & sim.time > t_pre);
stepctl = sim.output(1,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 = [];
sfr = [];
uctl = [];
sctl = [];

% Save the step sequences for overlay plot
stime = {};
ref = {};
ufr_ts = {};
ctl_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(stepctl(steps == usteps(i) & steptime > thistime(1) + t_ignore));
    thiss = std(stepctl(steps == usteps(i) & steptime > thistime(1) + t_ignore));
    uctl = [uctl thismean];
    sctl = [sctl thiss];
    
    % 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));
    ctl_ts{i} = sim.output(1, sim.time > (thistime(1) - pre_step) & sim.time < (thistime(end) + post_step));
    
end

unfr_avg = steps(1)/first_step_ratio;

%%
spnum = 5;
xlimit = [0 t_end];

close all
figure()
subplot(spnum,2,[1 2]);
hold on
ss.RandScat([],[],0);
xlim(xlimit)
%
subplot(spnum,2,[3 4]);

plot(sim.time(1:dsfac:end),sim.filt(1:dsfac:end),'b','LineWidth',0.5)
hold on
hold on
plot(sim.time,sim.target,'k','LineWidth',1)
legend('ref(t)','UNFR(t)','Location','NorthWest')
% plot(sim.time,sf,'b')
ylim([0 frmax])
xlim(xlimit)
ylabel('Pulse Variables')


subplot(spnum,2,[5 6]);
hold on
plot(sim.time(1:dsfac:end),sim.output(5:8,(1:dsfac:end))')
legend('Freq(t)','PW(t)','V_{I}(t)','Location','NorthWest')
% plot(sim.time,sf)
xlim(xlimit)
ylabel('Pulse Variables')

subplot(spnum,2,[7 8]);
hold on
plot(sim.time(1:dsfac:end),sim.output(1:4,(1:dsfac:end))')
legend('U(t)','P','I','D','Location','NorthWest')
% plot(sim.time,sf)
ylim([-ctlmax ctlmax])
xlim(xlimit)
ylabel('FB states')
xlabel('Time (sec)')

% Step overlay
subplot(spnum,2,9)
hold on
col = hsv(length(stime));
for i = 1:length(stime)
    
    plot(stime{i},ref{i},'k')
    plot(stime{i},ufr_ts{i},'Color',col(i,:))
end
% plot([0 0], [0.01 max(ref{end})],'k');
st = stime{1};
r = ref{1};
% plot([max(st(r > 0)) max(st(r > 0))], [0.01 max(ref{end})],'k');
xlim([-pre_step max(stime{1}) + post_step])
ylim([0.01 frmax])
ylabel('UNFR (Hz/Unit)')
xlabel('Time (sec)')

% Steady State mean tracking
subplot(spnum,2,10)
plot(usteps,usteps,'r-s')
hold on
plot(usteps,ufr,'k-o')
grid on
axis([0 20 0 20])
legend('<ref_N>','<UNFR>','Location','NorthWest')
ylabel('Hz/Unit')
xlabel('Hz/Unit')


%% Just the Steps
figure()

% Step overlay
subplot(2,3,[1 2])
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},100),'Color',col(i,:))
end
st = stime{1};
r = ref{1};
xlim([-pre_step ceil(max(stime{1}))])
ylim([0 frmax])
ylabel('UNFR (Hz/Unit)')
xlabel('Time (sec)')

% Average performance
subplot(2,3,3)
plot(usteps,usteps,'r-s')
hold on
errorbar(usteps,ufr,sfr,'k-o')
grid on
axis([0 20 0 20])
legend('<ref_N>','<UNFR>','Location','NorthWest')
ylabel('Hz/Unit')
xlabel('Hz/Unit')

% Control Signal
subplot(2,3,[4 5])
hold on
col = hsv(length(stime));
for i = 1:length(stime)
    plot(stime{i},smooth(ctl_ts{i},100),'Color',col(i,:))
end
st = stime{1};
r = ref{1};
xlim([-pre_step ceil(max(stime{1}))])
ylim([0 1.1])
ylabel('UNFR (Hz/Unit)')
xlabel('Time (sec)')


% Average performance
subplot(2,3,6)
errorbar(usteps,uctl,sctl,'k-o')
grid on
axis([0 20 0 1.1 ])
ylabel('<U(t)>')
xlabel('Hz/Unit')

