% Jon Newman

% Files that contain data corresponding protocol 09 (see experimental
% notes)

fids = {
    '../../2014_01_16_A_Extract/mat/P09_02';
    '../../2014_04_01_A_Extract/mat/P09_03'; % Unsorted
    '../../2014_04_01_A_Extract/mat/P09_05'; % MUA
    '../../2014_04_03_A_Extract/mat/Rec3_P09_01'; % Did not have successful clamp
    '../../2014_04_03_A_Extract/mat/Rec4_P09_01'; % Not in VPM
    '../../2014_04_03_A_Extract/mat/Rec5_P09_01';
    '../../2014_04_03_A_Extract/mat/Rec6_P09_01'; % Not in VPM
    '../../2014_04_03_A_Extract/mat/Rec7_P09_01'; % Not in VPM; forepaw area instead
    '../../2014_04_03_A_Extract/mat/Rec8_P09_01'; % Not in VPM
    '../../2014_04_06_A_Extract/mat/Rec4_P09_01'; % Not in VPM
    '../../2014_04_06_A_Extract/mat/Rec5_P09_01'; % Not in VPM; MUA?
    '../../2014_04_06_A_Extract/mat/Rec6_P09_01';
    '../../2014_04_06_A_Extract/mat/Rec9_P09_01';
    '../../2014_04_15_A_Extract/mat/Rec6_P09_01';
    '../../2014_04_15_A_Extract/mat/Rec7_P09_01';
    '../../2014_04_15_A_Extract/mat/Rec8_P09_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 = {
    [10 11];
    [];
    [];
    [];
    [];
    [26];
    [];
    [1];
    [];
    [];
    [1];
    [21 22 23];
    [21];
    [2];
    [];
    [19]
    };

no_units    =  [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]; % no. units controlled
use_idx    	=  [1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1]; % use this recording in the plot?
use_led     =  [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     =  [2 3 3 4 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'};

% Waveform sample rate
fs_ds_wave  = 1000;

% Time bin size
dt = 1; % sec

% TS example index
ts_idx = 8; %1,8

% Axes limits
max_fr = 90;
min_u = 0.009;
max_pow = 30;
bound_isi = [0.005 10];
lockcol = [[0,0,0];[1,0,0]];

% Firing rate color map[0 100 255]/255
% frcol = ((5:max_pow+5)'./max_pow)*([0 100 255]/255);
frcol = gray(40);
frcol = frcol(5:5+max_pow,:);
frcol = hsv(max_pow);

% Within each trial, over what time range should stationary statistics be
% evaluated
stat_bound = [15 30]; % last 15 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

%% Setup figure/axes
% Summarize firing statistics

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

% Axes dimensions
wts = 7.5;
aw = 3;
ah = aw;
as = 1.5;
sw = 1;
ms = 4;

a10 = axes('Units','centimeters','Position',[as        as+ah+as   wts     ah]);
a11 = axes('Units','centimeters','Position',[as    as   aw     ah]);
a12 = axes('Units','centimeters','Position',[2*as+aw   as   aw     ah]);
adumb = axes('Units','centimeters','Position',[10        as+ah+as   1     ah]);

% 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_syms = {syms{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';
    TU = 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);
        
        % 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)');
            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)');
            case 5
                load(opt_pow_maps{5});
                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(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);
    
    
    % Time series plots
    if j == ts_idx
        axes(a10);
        hold on
        for i = 1:length(uU)
            u_idx = ceil(uU(i));
            plot(T,(FILT(i,:)),'Color',frcol(u_idx,:));
            % plot([0 30],[TFR(i) TFR(i)],'k-');
        end
    end
    
    
    % Mean firing rate vs light power
    axes(a11);
    hold on
    plot(uU,uFR,'k-')
    %     for i = 1:length(uU)
    %         u_idx = ceil(uU(i));
    %         plot(uU(i),uFR(i),'Marker',these_syms{j},'Color',frcol(u_idx,:),'MarkerSize',ms,'LineStyle','none');
    %     end
    
    % ISI stats
    uu = unique(UNI{1});
    for k = 1:these_no_units(j)
        
        ISI = cell(1,numel(tri_use));
        uISI = zeros(1,numel(tri_use));
        sISI = zeros(1,numel(tri_use));
        
        for i = 1:length(uU)
            % ISI
            st = SPK{i};
            st = st(UNI{i} == uu(k));
            st = st(st > stat_bound(1) & st < stat_bound(2));
            ISI{i} = diff(st);
            uISI(i) = mean(ISI{i});
            sISI(i) = std(ISI{i});
        end
        
        axes(a12);
        hold on
        for i = 1:length(uU)
            u_idx = ceil(uU(i));
            plot(uISI(i),sISI(i),'Marker',these_syms{j},'Color',[0 0 0],'MarkerSize',ms,'LineStyle','none')
            hold on
        end
    end
end


%% Plot Options
set(a10,'Xlim',[-pre_sec post_sec],'Ylim',[0 60],'TickDir','out','Fontsize',8,'TickLength',[0.02 0.01],'YScale','linear')
xlabel(a10,'Time (sec)')
ylabel(a10,'FR (Hz)')

colormap(adumb,frcol);
colorbar('peer',adumb,'Location','EastOutside')

set(a11,'Xlim',[0 25],'Ylim',[0 max_fr],'TickDir','out','Fontsize',8,'TickLength',[0.02 0.01],'YScale','linear')
xlabel(a11,'Avg. power (mW \cdot mm^{-2})')
ylabel(a11,'Avg. FR (Hz)')

set(a12,'XScale','log','YScale','log','Fontsize',8)
xlabel(a12,'Avg. ISI','Fontsize',8)
ylabel(a12,'S.D. ISI','Fontsize',8)
y_lim = get(gca,'Ylim');
set(a12,'Ylim',bound_isi,'Xlim',bound_isi,'TickDir','out','Fontsize',8,'TickLength',[0.02 0.01],'box','off','XTick',[0.01 0.1 1 20],'YTick',[0.01 0.1 1 10])
plot(get(a12,'Xlim'),get(a12,'Ylim'),'k:')

%% Export

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

