% Jon Newman
% Examine CL time series quickly

fid = '2012_07_02_Xeko-20709_one-min-clamp_07';
sim = loadsimoc([fid '.simoc']);
spk = loadspk([fid '.spk']);

%% For getting average perfomance
t_ignore = 30; %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
b = sim.output(6,:) ~= 0 | sim.output(5,:) ~= 0;
stepobs = sim.obs(b) ;
stepfilt = sim.filt(b);
steptime = sim.time(b);
stepctl1 = sim.output(6,b);
stepctl2 = sim.output(5,b);
steps = sim.target(b);
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 = {};
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(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 = sim.time > (thistime(1) - pre_step) & sim.time < (thistime(end) + post_step);
    stime{i} = sim.time(b) - thistime(1);
    ref{i} = sim.target(b);
    ufr_ts{i} = sim.filt(b);
    ctl1_ts{i} = sim.output(6,b);
    ctl2_ts{i} = sim.output(5,b);
    
end


avg_perf.steps = usteps;
avg_perf.ufr = ufr;
avg_perf.sfr = sfr;
avg_perf.uctl1 = uctl1;
avg_perf.sctl1 = sctl1;
avg_perf.uctl2 = uctl2;
avg_perf.sctl2 = sctl2;

ts.t = stime;
ts.ref = ref;
ts.fr = ufr_ts;
ts.ctl1 = ctl1_ts;
ts.ctl2 = ctl2_ts;


% Spont data estimate window
spont_t = [1200 1400];
avg_perf.u_spont_fr = numel(ss.time(ss.unit ~= 0 & ss.time > spont_t(1) & ss.time <= spont_t(2)))/range(spont_t)/numel(unique(ss.unit(ss.unit~=0)));


save([fid '.mat'],'avg_perf','ts')


