% 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';
    '../../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, but does not matter for simple steps
    '../../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 nan 1 1 1 1]; % no. units controlled
gd_unit     =  [7 7 7 7 7 7 7 7 7 1 6 6 6 6 6 6 6 6 6 8 6 8 7 6];
% use_idx    	=  [0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0]; % use this recording in the plot?
use_idx    	=  [0 0 0 0 0 0 0 0 0 1 1 0 0 1 1 1 1 0 0 0 1 1 1 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        =  {'s','o','^','d','x','v','+','<','>','*','h','p','s','o','^','d','x','v','+','<','>','*','h','p'};

% 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

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

% 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 = 1;
fail_error = 2;

%% Setup figure/axes
% Summarize firing statistics

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

% Axes dimensions
aw = 3;
ah = aw;
as = 1.2;
sw = 1;
ms = 4;

zoom = 1;
a11 = axes('Units','centimeters','Position',zoom*[as         2*as   aw     ah]);
a12 = axes('Units','centimeters','Position',zoom*[2*as+aw    2*as   aw     ah]);
a13 = axes('Units','centimeters','Position',zoom*[3*as+2*aw  2*as   aw     ah]);
a15 = axes('Units','centimeters','Position',zoom*[4*as+3*aw  2*as   aw     ah]);
a14 = axes('Units','centimeters','Position',zoom*[5*as+4*aw  2*as   aw     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_gd_unit = gd_unit(logical(use_idx));
these_syms = unique(syms); % {syms{logical(use_idx)}};
these_use_led = use_led(logical(use_idx));
these_pow_map = pow_map(logical(use_idx));

% Across cells stats
POW = [];
RMS = [];
SETTLE = [];
CV = [];
ACHV = [];
ACHR = [];
SPONTR = [];
POW_CELL = cell(length(these_fids),1);

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(numel(tri_use),1);
    TFRS = zeros(numel(tri_use),1);
    ST = 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 = {};
    
    if isempty(tri_use)
        continue
    end
    
    wave.series.ind = wave.series.ind(:);
    
    kk = 1;
    SS.ResetClean
    
    SS.RemoveUnit(setdiff(unique(SS.unit),these_gd_unit(j)));
    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);
        
        if strcmp('../../2014_04_15_A_Extract/mat/Rec6_P01_02', block_loc) || strcmp('../../2014_04_01_A_Extract/mat/P01_04', block_loc)
            TFR(kk,1) = wave.series.targ(wave.series.ind_trial(i) + 100);
            FILT(kk,size(U,2)-sum(b)+1:end) = wave.series.filt(b)';
            TFRS(kk,1) =  TFR(kk,1)*2;
        else
            TFR(kk,1) = wave.series.targ(wave.series.ind_trial(i) + 100);
            FILT(kk,size(U,2)-sum(b)+1:end) = wave.series.filt(b)';
            TFRS(kk,1) =  TFR(kk,1)*1;
        end
        
        % Get the correct light power mapping
        switch these_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)');
        end
        
        
        
        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)];
        
        
%         % Get the settling time for this step
%         sw = 2.5; %sec
%         dfr_s = 0.25;
%         st_time = 1; % second
%         b = wave.series.ind > wave.series.ind_trial(i) & wave.series.ind < (wave.series.ind_trial(i) + post_samp);
%         tru_fr = wave.series.filt(b);
%         sm_fr = smooth(wave.series.filt(b),round(sw*wave.meta.fs_ds_wave),'lowess');
%         
%         ind_in_bound = find(abs(sm_fr - TFR(kk,1)) < dfr_s); % & t > 2);
%         ind_out_bound = find(abs(sm_fr - TFR(kk,1)) > dfr_s);
%         
%         if  isempty(ind_in_bound)
%             continue;
%         elseif size(ind_in_bound,2) > 1
%             ind_in_bound = ind_in_bound';
%             ind_out_bound = ind_out_bound';
%         end
%         
%         ind_enter_bound = [ind_in_bound(1);  ind_in_bound(find(diff(ind_in_bound) > 1)+1)];
%         ind_exit_bound = [ind_out_bound(1); ind_out_bound(find(diff(ind_out_bound) > 1)+1)];
%         
%         
%         if ind_exit_bound(1) < ind_enter_bound(1)
%             ind_exit_bound(1) = [];
%         end
%         if length(ind_exit_bound) ~= length(ind_enter_bound)
%             ind_exit_bound = [ind_exit_bound; length(sm_fr)];
%         end
%         
%         t_in_bound = ind_exit_bound-ind_enter_bound;
%         ind_settled = ind_enter_bound(find(t_in_bound > st_time*wave.meta.fs_ds_wave,1,'first'));
%         
%         if ~isempty(ind_settled)
%             ST(kk,1) = ind_settled/wave.meta.fs_ds_wave;
%         end
% %         
% %         % Sanity Check
% %         figure
% %         hold on
% %         plot(tru_fr,'y')
% %         plot(sm_fr,'k')
% %         plot(ind_settled,sm_fr(ind_settled),'r*')
% %         close gcf
        
        kk = kk + 1;
    end
    
    % Time-averaged mean/std firing rates
    if ~strcmp(block_loc,'../../2014_04_15_A_Extract/mat/Rec8_P01_01')
        preFR = mean(FILT(:,T < 0),2);
    else
        preFR = nan; % this one did not have down time between steps.
    end
    uFR = mean(FILT(:,T > stat_bound(1) & T < stat_bound(2)),2);
    %sFR = std(FILT(:,T > stat_bound(1) & T < stat_bound(2)),[],2);
    
    % RMS error
    flt = FILT(:,T > stat_bound(1) & T < stat_bound(2));
    rmsError = rms(flt - TFR(:,ones(1,size(flt,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);
    
    % Acheived targets
    b = abs(uFR - TFR) <= fail_dev;
    
    % Mean firing rate plot
    acheived_targs = [];
    acheived_rates = [];
    axes(a11);
    hold on
    SPONTR = [SPONTR nanmean(preFR)];
    plot(0, nanmean(preFR),'k','Marker',these_syms{j},'LineStyle','none','MarkerSize',ms)
    plot(TFR,uFR,'k-');
    for i = 1:length(TFR)
        fr_idx = ceil(TFR(i));
        if abs(uFR(i) - TFR(i)) > fail_dev
            plot(TFR(i),uFR(i),'Marker',these_syms{j},'Color',[0.5 0.5 0.5],'MarkerSize',ms,'LineStyle','none');
        else
            plot(TFR(i),uFR(i),'Marker',these_syms{j},'Color',frcol(fr_idx,:),'MarkerSize',ms,'LineStyle','none');
            acheived_targs = [acheived_targs, TFR(i)];
            acheived_rates = [acheived_rates, uFR(i)];
        end
    end
    
    ACHV = [ACHV, acheived_targs];
    ACHR = [ACHR, acheived_rates];
    
    % Mean control signal plot
    axes(a12);
    hold on
    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',these_syms{j},'Color',[0.5 0.5 0.5],'MarkerSize',ms,'LineStyle','none')
        else
            plot(TFR(i),uU(i),'Marker',these_syms{j},'Color',frcol(fr_idx,:),'MarkerSize',ms,'LineStyle','none')
            POW = [POW uU(i)];
            POW_CELL{j} = [POW_CELL{j}, [TFR(i); uU(i)]];
        end
    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(TFR)
        
        % ISI
        st = SPK{i};
        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(a13);
    hold on
    
    plot(TFR(b),rmsError(b),'k-');
    for i = 1:length(TFR)
        fr_idx = ceil(TFR(i));
        if abs(uFR(i) - TFR(i)) > fail_dev
            %                 plot(TFR(i),rmsError(i),'Marker',these_syms{j},'Color',[0.5 0.5 0.5],'MarkerSize',ms,'LineStyle','none')
        else
            plot(TFR(i),rmsError(i),'Marker',these_syms{j},'Color',frcol(fr_idx,:),'MarkerSize',ms,'LineStyle','none')
            RMS = [RMS rmsError(i)];
        end
    end
    
%     axes(a14);
%     hold on
%     plot(TFR,ST,'k-');
%     for i = 1:length(TFR)
%         fr_idx = ceil(TFR(i));
%         if abs(uFR(i) - TFR(i)) > fail_dev
%             %                 plot(TFR(i),ST(i),'Marker',these_syms{j},'Color',[0.5 0.5 0.5],'MarkerSize',ms,'LineStyle','none')
%         else
%             plot(TFR(i),ST(i),'Marker',these_syms{j},'Color',frcol(fr_idx,:),'MarkerSize',ms,'LineStyle','none')
%             SETTLE = [SETTLE ST(i)];
%         end
%     end
    
    axes(a15);
    hold on
    for i = 1:length(TFR)
        fr_idx = ceil(TFRS(i));
        if abs(uFR(i) - TFR(i)) > fail_dev
            %                 plot(uISI(i),sISI(i),'Marker',these_syms{j},'Color',[0.5 0.5 0.5],'MarkerSize',ms,'LineStyle','none')
        else
            plot(uISI(i),sISI(i),'Marker',these_syms{j},'Color',frcol(fr_idx,:),'MarkerSize',ms,'LineStyle','none')
            CV = [CV sISI(i)/uISI(i)];
        end
    end
    
    
end

% end


%% Plot Options

% Axes limits
max_fr = 50;
min_u = 0.1;
min_u = 0;
max_u = 30;
min_rms = 0;
max_rms = 4;
bound_isi = [0.01 1];
lockcol = [[0,0,0];[1,0,0]];


set(a11,'Xlim',[0 max_fr],'Ylim',[0 max_fr],'TickDir','out','Fontsize',8,'TickLength',[0.02 0.01],'XTick',[0 10 20 30 40 50])
plot(a11,get(a11,'Xlim'),get(a11,'Ylim'),'k:')
xlabel(a11,'Target FR (Hz) ')
ylabel(a11,'Avg. FR (Hz)')

set(a12,'Xlim',[0 max_fr],'Ylim',[min_u max_u],'TickDir','out','Fontsize',8,'TickLength',[0.02 0.01],'YScale','linear','XTick',[0 10 20 30 40 50],'YTick',[0:5:30]) % ,'YTick',[0.01 0.1 1 10 100 1000]
xlabel(a12,'Target FR (Hz) ')
ylabel(a12,'Avg. optical power (mW \cdot mm^{-2})')

set(a13,'Xlim',[0 max_fr],'Ylim',[min_rms max_rms],'TickDir','out','Fontsize',8,'TickLength',[0.02 0.01],'YScale','linear','XTick',[0 10 20 30 40 50])
xlabel(a13,'Target FR (Hz) ')
ylabel(a13,'RMS Error (Hz)')

set(a14,'Xlim',[0 max_fr],'Ylim',[0.1 20],'TickDir','out','Fontsize',8,'TickLength',[0.02 0.01],'YScale','log','XTick',[0 10 20 30 40 50])
xlabel(a14,'Target FR (Hz) ')
ylabel(a14,'Settling time (sec)')

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

%% Firing rate histogram
a16 = axes('Units','centimeters','Position',zoom*[6*as+5*aw  2*as   aw     ah]);

bins = 0:5:50;
c_spont = histc(SPONTR, 0:5:50);
c_clamp = histc(ACHR, 0:5:50);
hold on
plot(a16,[0 c_spont(1:end-1)/max(c_spont(1:end-1))], [ 0 bins(1:end-1) + 2.5])
plot(a16,[0 c_clamp(1:end-1)/max(c_clamp (1:end-1))], [ 0 bins(1:end-1) + 2.5],'r')

set(a16,'Ylim',[0 max_fr],'Xlim',[0 1.2],'TickDir','out','Fontsize',8,'TickLength',[0.02 0.01],'box','off')

%% Export

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

% Save mean optical powers for each target
save('P01_in-vivo_optical_power','POW_CELL');

%% Means

mean(RMS)
std(RMS)

mean(CV)
std(CV)

mean(SETTLE)
std(SETTLE)


% %% Stats
% conf = .95;
% a = 1-conf;
% b1 = norminv(1-a/2);
% b2 = b1/sqrt(2);
% 
% rfac = (1/(norminv(.75)-norminv(.25))) * sqrt(pi/2) * (b1+b2)/2;
% dat = SETTLE;
% n = numel(dat);
% q = quantile(dat,[0.25, 0.5, 0.75]);
% q(2)
% ci_u = q(2) + rfac*(q(3)-q(1))/sqrt(n)
% ci_l = q(2) - rfac*(q(3)-q(1))/sqrt(n)


