% Jon Newman
% Examine the integrated optical delivery for different control types. The
% plot produced by this file should be combined with the same one for the
% in-vitro data.


% 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';
    '../../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
    '../../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];
    []; % last
    [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_idx    	=  [1 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1 1 0 0 1 1 1 1 1]; % use this recording in the plot?
% use_idx    	=  [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]; % use this recording in the plot?
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        =  {'x','.'};

% Waveform sample rate
fs_ds_wave  = 1000;

% Taus we are interested in
Tau = 1;

% Ti's we are interested in
Ti = 2;

% Time bin size
dt = 1; % sec

% Axes limits
max_fr = 60;
min_u = 0.009;
max_u = 1000;
bound_isi = [0.005 10];
lockcol = [[0,0,0];[1,0,0]];

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

%% Setup figure/axes
% Summarize firing statistics

close all
f = figure('Units','centimeters','Position',[5 5 25 10]);

% close all
figure('units','centimeter','position',[5 5 20 20]);
a1 = axes('units','centimeter','position',[5 5 5 5]);
a2 = axes('units','centimeter','position',[12 5 5 5]);

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

% Filter data to get the trials we want to process
these_fids = {fids{logical(use_idx)}};
these_exclude = {tri_exclude{logical(use_idx)}};
these_no_units = no_units(logical(use_idx));
these_use_led = use_led(logical(use_idx));

for j = 1:length(these_fids)
    
    % Get the data
    block_loc = these_fids{j};
    load(block_loc);
    load([block_loc '.SS'],'-mat');
    
    % Which trials do we want to examine
    tri_use = setdiff(1:length(wave.series.ind_trial),these_exclude{j})
    
    % Storage
    T   = (-pre_samp:post_samp)/fs_ds_wave';
    TFR = zeros(1,numel(tri_use));
    U   = nan(1, numel(-pre_samp:post_samp));
    FILT   = nan(1, numel(-pre_samp:post_samp));
    SPK = {};
    NO_UNI = [];
    UNI = {};
    
    if isempty(tri_use)
        continue
    end
    
    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(1,kk) = wave.series.targ(wave.series.ind_trial(i) + 100);
        
        % Get the correct light power mapping
        switch pow_map(j)
            
            case 1
                load(opt_pow_maps{1});
                U(kk,size(U,2)-sum(b)+1:end) = volt2power(10*wave.series.out1(b)');
                sym = syms{1};
            case 2
                load(opt_pow_maps{2});
                U(kk,size(U,2)-sum(b)+1:end) = vpd2power(wave.series.LED_pow(b)');
                sym = syms{2};
            case 3
                load(opt_pow_maps{3});
                U(kk,size(U,2)-sum(b)+1:end) = vpd2power(wave.series.LED_pow(b)');
                sym = syms{2};
            case 4
                load(opt_pow_maps{4});
                U(kk,size(U,2)-sum(b)+1:end) = vpd2power(wave.series.LED_pow(b)');
                sym = syms{2};
            case 5
                load(opt_pow_maps{5});
                U(kk,size(U,2)-sum(b)+1:end) = vpd2power(wave.series.LED_pow(b)');
                sym = syms{2};
        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(j)];
        kk = kk + 1;
    end
    
    % Time-averaged mean/std firing rates
    uFR = mean(FILT(:,T > stat_bound(1) & T < stat_bound(2)),2);
    sFR = std(FILT(:,T > stat_bound(1) & T < stat_bound(2)),[],2);
    
    % Time-averaged mean/std control signal
    uU = mean(U(:,T > stat_bound(1) & T < stat_bound(2)),2);
    uU(uU < 0.01) = 0;
    sU = std(U(:,T > stat_bound(1) & T < stat_bound(2)),[],2);

    % Mean control signal plot
    axes(a1);
    hold on
    b = abs(uFR - TFR') < fail_dev;
    
    plot(TFR(b),uU(b),'k-');
    for i = 1:length(TFR)
        fr_idx = ceil(TFR(i));
        
        if abs(uFR(i) - TFR(i)) < fail_dev
            plot(TFR(i),uU(i),'Marker',sym,'Color',frcol(fr_idx,:),'MarkerSize',ms,'LineStyle','none')
        end
    end
end

set(a1,'TickDir','out','Ylim',[0.0005 2],'Xlim',[0 13])
set(a2,'TickDir','out','Ylim',[0.05 15],'Xlim',[0 13])

% %% Export
% 
% export_fig('P01_stat-summary_01','-pdf','-transparent',gcf)
% saveas(gcf,'P01_stat-summary_01')

