% Jon Newman
% Examine CL time series quickly

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

%% First file
fid = '2012_10_22_Herb-I5386_one-min-clamp_09';
pic = loadpic([fid '.pic']);
spk = loadspk([fid '.spk']);

ss = SqueakySpk(fid,spk.meta.fs_Hz,1,spk);

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

% 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);


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

for i = 1:length(usteps)

        step = [step usteps(i)];
        
        % 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


%% Second file
fid = '2012_10_22_Herb-I5386_one-min-clamp_10';
pic = loadpic([fid '.pic']);
spk = loadspk([fid '.spk']);

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

% 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);



k = length(stime);

for i = 1:length(usteps)
    if usteps(i) == 2

        % 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(step == 2) = thismean;
        sfr(step == 2) = thiss;
        
        b = pic.time > (thistime(1) - pre_step) & pic.time < (thistime(end) + post_step);
        stime{k+1} = pic.time(b) - thistime(1);
        ref{k+1} = pic.targ(b);
        ufr_ts{k+1} = pic.filt(b);
        ctl1_ts{k+1} = pic.out1(b);
        ctl2_ts{k+1} = pic.out2(b);
        
        k = k+1;
 
    end
    
end

%%
[step ind] = sort(step);
avg_perf.steps = step;
avg_perf.ufr = ufr(ind);
avg_perf.sfr = sfr(ind);
avg_perf.uctl1 = uctl1(ind);
avg_perf.sctl1 = sctl1(ind);
avg_perf.uctl2 = uctl2(ind);
avg_perf.sctl2 = sctl2(ind);


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 1900];
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')