% Jon Newman

% Files that contain data corresponding protocol 01 (see experimental
% notes)
fids = {
    '../../2013_05_23_A_ExtractB/mat/PI_MUA_FRStep_09';
    '../../2013_05_23_A_ExtractB/mat/PI_MUA_FRStep_10';
    '../../2013_05_23_A_ExtractB/mat/PI_MUA_FRStep_12';
    '../../2013_12_10_A_Extract/mat/P01_01';
    '../../2013_12_10_A_Extract/mat/P01_04';
    '../../2013_12_10_A_Extract/mat/P01_07';
    '../../2013_12_13_A_Extract/mat/P01_02';
    '../../2013_12_13_A_Extract/mat/P01_04';
    '../../2014_01_16_A_Extract/mat/P01_05';
    '../../2014_01_16_A_Extract/mat/P01_06';
    '../../2014_01_16_A_Extract/mat/P01_08';
    '../../2014_04_01_A_Extract/mat/P01_04';%XX
    '../../2014_04_03_A_Extract/mat/Rec4_P01_01'; % MUA; Not in VPM
    '../../2014_04_03_A_Extract/mat/Rec6_P01_01';
    '../../2014_04_03_A_Extract/mat/Rec7_P01_02';
    '../../2014_04_06_A_Extract/mat/Rec1_P01_01';
    '../../2014_04_06_A_Extract/mat/Rec4_P01_01'; % Not in VPM %XX
    '../../2014_04_06_A_Extract/mat/Rec5_P01_01'; % MUA
    '../../2014_04_06_A_Extract/mat/Rec6_P01_01'; % MUA
    '../../2014_04_06_A_Extract/mat/Rec8_P01_02'; 
    '../../2014_04_06_A_Extract/mat/Rec9_P01_01'; 
    '../../2014_04_15_A_Extract/mat/Rec6_P01_02'; 
    '../../2014_04_15_A_Extract/mat/Rec7_P01_01';
    '../../2014_04_15_A_Extract/mat/Rec8_P01_01';
    };

opt_pow_maps = {
    '../../optical-power/2014_01_19_PowerMeasures/volt2power_laser.mat';
    '../../optical-power/2014_01_19_PowerMeasures/vpd2power_led.mat';
    '../../optical-power/2014_04_01_PowerMeasures/vpd2power_led.mat';
    '../../optical-power/2014_04_03_PowerMeasures/vpd2power_led.mat'; 
    };

% Corresponds to fids cell above and provides trials that should be
% excluded from plots and averages (e.g. when we turned the contoller off
% but the system continued to mark trials and continued recording)
tri_exclude = {
    [22 21 20];
    [1 2];
    [1 2 12 13];
    [];
    [];
    [1 2 6 7 8 9 10];
    [10];
    [26 27 28 29];
    [1 22];
    [1 20];
    [22 23];
    [1 2 3 4 17 18 19 20 21 22 23 24 25 26];
    []; 
    [12 13];
    [9 10 11];
    [1 2 3 8];
    [1 2 3 4 23];
    [];
    [];
    [1 2 3 4 11 12 13];
    [1 2 7 8 9 10]
    [];
    [1 2 18];
    [1 2 3 4 5 6 18 19 20];
    };

tau_idx     =  [2 2 2 2 2 2 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2]; % 0 = mixed, 1 = 0.16 sec, 2 = 0.8 sec, 3 = 1.6 sec
ti_idx      =  [2 2 2 2 0 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2]; % 0 = mixed, 1 = 0.1 sec, 2 = 1 sec, 3 = 10 sec
no_units    =  [1 1 1 2 2 1 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]; % no. units controlled
use_led     =  [0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]; % Was the LED being used to stimulate?
pow_map     =  [1 1 1 1 1 1 1 1 2 2 2 3 4 4 4 4 4 4 4 4 4 4 4 4]; % Which light power map to use for the fiber in this experiment, see list above
syms        =  {'s','o','^','d','.','x','v','+','<','>','*','h','p','s','o','^','d','.','x','v','+','<','>','*','h','p'};

% Which fid do you want to create a ts for
fid_idx = 14;

% Waveform sample rate
fs_ds_wave  = 1000;

% ASDR bin
dt = 1; % sec

% Firing rate color map
frcol = hsv(max_fr);

% Within each trial, over what time range should stationary statistics be
% evaluated
stat_bound = [15 30]; % last 25 seconds

% Pre/post time around each trial to use for data display
pre_sec = 10;       % time before each step to show
post_sec = 40;      % time following each step to show

% Firing rate deviation that consitutes control failure
fail_dev = 0.25;
fail_error = 2;


%% Setup figure/axes
% Summarize firing statistics

% Load the selected file
block_loc = fids{fid_idx};
load(block_loc);
load([block_loc '.SS'],'-mat');

% Filter data to get the trials we want to process
these_exclude = tri_exclude{fid_idx};
these_no_units = no_units(fid_idx);
these_syms = syms{fid_idx};
these_use_led = use_led(fid_idx);
these_pow_map = pow_map(fid_idx);

% Timing
pre_samp = ceil(pre_sec * fs_ds_wave);
post_samp = ceil(post_sec * fs_ds_wave);
time = (-pre_samp:post_samp) * (1/fs_ds_wave);

% Extract the data
% Which trials do we want to examine
tri_use = setdiff(1:length(wave.series.ind_trial),these_exclude)

% Storage
T   = (-pre_samp:post_samp)/fs_ds_wave';
TFR = zeros(numel(tri_use),1);
U   = nan(1, numel(-pre_samp:post_samp));
FILT   = nan(1, numel(-pre_samp:post_samp));
SPK = {};
NO_UNI = [];
UNI = {};
wave.series.ind = wave.series.ind(:);

kk = 1;
spk = SS.ReturnClean(); % Because during the extraction process, those units not used for FR clacs were marked as dirty

for i = tri_use
    
    b = wave.series.ind > (wave.series.ind_trial(i) - pre_samp) & wave.series.ind < (wave.series.ind_trial(i) + post_samp);
    
    TFR(kk,1) = wave.series.targ(wave.series.ind_trial(i) + 100);
    
    % Get the correct light power mapping
    switch these_pow_map
        
        case 1
            load(opt_pow_maps{1});
            U(kk,size(U,2)-sum(b)+1:end) = volt2power(10*wave.series.out1(b)');
        case 2
            load(opt_pow_maps{2});
            U(kk,size(U,2)-sum(b)+1:end) = vpd2power(wave.series.LED_pow(b)');
        case 3
            load(opt_pow_maps{3});
            U(kk,size(U,2)-sum(b)+1:end) = vpd2power(wave.series.LED_pow(b)');
        case 4
            load(opt_pow_maps{4});
            U(kk,size(U,2)-sum(b)+1:end) = vpd2power(wave.series.LED_pow(b)');
    end
    
    FILT(kk,size(U,2)-sum(b)+1:end) = wave.series.filt(b)';
    
    b = spk.time > (wave.series.ind_trial(i)/wave.meta.fs_ds_wave - pre_sec - dt) & spk.time < (wave.series.ind_trial(i)/wave.meta.fs_ds_wave + post_sec);
    SPK{kk} = spk.time(b) - wave.series.ind_trial(i)/wave.meta.fs_ds_wave;
    UNI{kk} = spk.unit(b);
    NO_UNI = [NO_UNI; these_no_units];
    kk = kk + 1;
end


%% Make plots
close all
f = figure('Units','centimeters','Position',[5 5 30 20]);

% Target rates to plot
f_star = [8 12 16 20 24];

% Axes limits
s_lim = [12.5 13.5]; % Time to get spikes from
t_lim = [-pre_sec post_sec];
fr_lim = [0 35];
u_lim = [0 15];

% Axes dimensions
zoom = 1;
aw = 1.75;
ah = 1;
sw = 0.15;
sh = 0.15;
wid = aw + sw;
hei = ah + sh;

for f = 1:numel(f_star)
    
    b = TFR == f_star(f);
    
    % Raster
    axes('Units','centimeters','Position',zoom*[2+wid*(f-1)         2+2*hei     aw     ah/6]);
    st = SPK{b};
    st = st(st >= s_lim(1) & st <= s_lim(2));
    for i = 1:numel(st)
       line([st(i) st(i)],[0 1]) 
       hold on
    end
    set(gca,'Xlim',s_lim,'Ylim',[0 1],'TickDir','out','Fontsize',7,'TickLength',[0.02 0.01],'XTick',[],'YTick',[])
    
    % Firing rate
    axes('Units','centimeters','Position',zoom*[2+wid*(f-1)         2+hei       aw     ah]);
    hold on
    plot(T, FILT(b,:),'b-')
    plot([0 30],[f_star(f) f_star(f)],'r')
    set(gca,'Xlim',t_lim,'Ylim',fr_lim,'TickDir','out','Fontsize',7,'TickLength',[0.02 0.01],'XTick',[-10 0 10 20 30 40])
    
    % Control signal
    axes('Units','centimeters','Position',zoom*[2+wid*(f-1)         2           aw     ah]);    
    plot(T, U(b,:),'b-') 
     set(gca,'Xlim',t_lim,'Ylim',u_lim,'TickDir','out','Fontsize',7,'TickLength',[0.02 0.01],'XTick',[-10 0 10 20 30 40])
end


%% Export

export_fig('P01_ts-summary_01','-pdf','-transparent',gcf)
saveas(gcf,'P01_ts-summary_01')

