% Jon Newman

% Files that contain data corresponding protocol 01 (see experimental
% notes)
fids = {
    '../../2014_04_15_A_Extract/mat/Rec6_P08_02';
    '../../2014_04_15_A_Extract/mat/Rec7_P08_02';
    '../../2014_04_15_A_Extract/mat/Rec8_P08_02';
    };

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 = {
    [];
    [];
    [];
    };


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

% Waveform sample rate
fs_ds_wave  = 1000;

% Taus we are interested in
Taus = 5/(2*pi);

% Ti's we are interested in
Ti = [0.1 1 10];

% Time bin size
dt = 1; % sec

% Target FR we are interested in
targ_fr = 10;

% Axes limits
max_fr = 100;
max_fr_ts = 100;
max_u = 600;
bound_isi = [0.005 10];
lockcol = [[0,0,0];[1,0,0]];
cellmark = {'*','^'};
frcol = [];

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

% Time series index
ts_idx = 1;

%% Setup figure/axes
% Summarize firing statistics
close all
f = figure('Units','centimeters','Position',[5 5 20 20]);

% Axes dimensions
aw = 3;
ah = 2;
as = 0.33;
sw = 1;
ms = 3;
ab = 2;

a(1,1) = axes('Units','centimeters','Position',[ab         ab+3*as+2*ah   aw     ah/3]);
a(1,2) = axes('Units','centimeters','Position',[ab+1*as+aw    ab+3*as+2*ah   aw     ah/3]);
a(1,3) = axes('Units','centimeters','Position',[ab+2*as+2*aw  ab+3*as+2*ah   aw     ah/3]);
a(2,1) = axes('Units','centimeters','Position',[ab         ab+2*as+ah     aw     ah]);
a(2,2) = axes('Units','centimeters','Position',[ab+1*as+aw    ab+2*as+ah     aw     ah]);
a(2,3) = axes('Units','centimeters','Position',[ab+2*as+2*aw  ab+2*as+ah     aw     ah]);
a(3,1) = axes('Units','centimeters','Position',[ab             ab+as          aw     ah]);
a(3,2) = axes('Units','centimeters','Position',[ab+1*as+aw    ab+as          aw     ah]);
a(3,3) = axes('Units','centimeters','Position',[ab+2*as+2*aw  ab+as          aw     ah]);

a(2,4) = axes('Units','centimeters','Position',[ab+5*as+3*aw  ab+as     ah     (2*ah+2*as+ah/3)]);
% a(3,4) = axes('Units','centimeters','Position',[ab+5*as+3*aw  ab+as          ah     (2*ah+2*as+ah/3)/2 - as/2]);

% 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));
these_pow_map = pow_map(logical(use_idx));

% these_fids = {fids{tau_idx == k & use_idx}};
% these_exclude = {tri_exclude{tau_idx == k & use_idx}};
% these_no_units = no_units(tau_idx == k & 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});
    % tri_use = tri_use(wave.series.targ(wave.series.ind_trial(tri_use)+ 100) == targ_fr);
    
    % 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));
    TAU = [];
    TI = [];
    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(kk,1) = wave.series.targ(wave.series.ind_trial(i) + 100);
        TI(kk,1) = wave.series.Ti(wave.series.ind_trial(i) + 100);
        TAU(kk,1) = wave.series.fr_tau(wave.series.ind_trial(i) + 100);
        
        % 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
        
        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
    
    % Get the average firing rate series
    bins = (-pre_sec -dt):dt:post_sec ;
    FR = zeros(size(U,1),size(bins,2));
    for i = 1:length(SPK)
        these_spks = SPK{i};
        FR(i,:) = histc(these_spks, bins)/NO_UNI(i)/dt;
    end
    bins = bins + dt;
    
    % Time series plots
    
    
    para = unique(TAU);
    
    uUX = [];
    uRX = [];
    
    for p = 1:numel(para)
        
        % Show trial rasters for each value of para
        this_spk = SPK(TAU == para(p));
        
        if j == ts_idx
            axes(a(1,p))
            hold on
            spk.time = [];
            spk.channel = [];
            
            cn = 1;
            for s = 1:numel(this_spk)
                spk.time = [spk.time; this_spk{s}];
                spk.channel = [spk.channel; cn*ones(size(this_spk{s}))];
%                 spk.time = [spk.time; 0.1];
%                 spk.channel = [spk.channel; cn+1];
                
                cn = cn + 1;
%                 plot(this_spk{s},s*ones(size(this_spk{s})) + 0.15*randn(size(this_spk{s})),'b.','MarkerSize',1)
            end
            
            spk.time = spk.time + pre_sec + dt;
            
            ss = SqueakySpk([],SS.fs,1,spk);
            ss.ASDR(dt,[],[],[],[],0);
            ss.csdr.csdr = ss.csdr.csdr/NO_UNI(1);
            ss.PlotCSDR(30,0)
        end
        
        % Get stats for each value of para
        
        % Firing rate
        this_fr = FR(TAU == para(p),:);
        this_targ = TFR(TAU == para(p),:);
        uFR = nanmean(this_fr);
        sFR = nanstd(this_fr);
        %             qFR =   quantile(this_fr,[0.25, 0.5, 0.75]);
        %     sFR =   sFR([1,3],:);
        
        if j == ts_idx
            hold(a(2,p),'on')
            %             plot(a(2,p),bins,this_fr,'color',[0.5 0.5 0.5])
            plot(a(2,p),bins,uFR,'g-')
            plot(a(2,p),[0 30],[TFR(1) TFR(1)],'r')
            plot(a(2,p),bins,uFR-sFR,'b-')
            plot(a(2,p),bins,uFR+sFR,'b-')
        end
        
        % Light power
        this_u = U(TAU == para(p),:);
        uU = nanmean(this_u);
        sU = nanstd(this_u);
        
        if j == ts_idx
            hold(a(3,p),'on')
            %             plot(a(3,p),time,this_u,'-','color',[0.5 0.5 0.5])
            plot(a(3,p),time,uU,'g-')
            plot(a(3,p),time,uU-sU,'b-')
            plot(a(3,p),time,uU+sU,'b-')
        end
        
        % RMS error
        ctl_fr = this_fr(:,bins > 0 & bins < 30);
        this_error = ctl_fr - this_targ(:,ones(1,size(ctl_fr,2)));
        this_rms = rms(this_error,2);
        hold(a(2,4),'on')
        plot(a(2,4),p*ones(size(this_rms))+ 0.1*randn(size(this_rms)),this_rms,'b','Marker',these_syms{j},'LineStyle','none','MarkerSize',ms)
        
        uRX = [uRX mean(this_rms)];
        
%         % Light power
%         uU = nanmean(this_u(:,time > 0 & time < 30),2);
%         hold(a(3,4),'on')
%         plot(a(3,4),p*ones(size(uU))+ 0.1*randn(size(uU)),uU,'b','Marker',these_syms{j},'LineStyle','none','MarkerSize',ms)
%         uUX = [uUX mean(uU)];
        
    end
    
    plot(a(2,4),[1 2 3],uRX,'r-','MarkerSize',ms,'Marker',these_syms{j})
%     plot(a(3,4),[1 2 3],uUX,'r-','MarkerSize',ms)
end



%% Plot options
set(a(1,1),'Xlim',[1 51],'Ylim',[0 9],'TickDir','out','Fontsize',7,'TickLength',[0.01 0.01],'box','off','FontName','Arial','XTickLabel',{},'XTick',[-10:10:40])
set(a(1,2),'Xlim',[1 51],'Ylim',[0 9],'TickDir','out','Fontsize',7,'TickLength',[0.01 0.01],'box','off','FontName','Arial','XTickLabel',{},'YTickLabel',{},'XTick',[-10:10:40])
set(a(1,3),'Xlim',[1 51],'Ylim',[0 9],'TickDir','out','Fontsize',7,'TickLength',[0.01 0.01],'box','off','FontName','Arial','XTickLabel',{},'YTickLabel',{},'XTick',[-10:10:40])
set(a(2,1),'Xlim',[-pre_sec post_sec],'Ylim',[0 30],'TickDir','out','Fontsize',7,'TickLength',[0.01 0.01],'box','off','FontName','Arial','XTickLabel',{},'XTick',[-10:10:40])
set(a(2,2),'Xlim',[-pre_sec post_sec],'Ylim',[0 30],'TickDir','out','Fontsize',7,'TickLength',[0.01 0.01],'box','off','FontName','Arial','XTickLabel',{},'YTickLabel',{},'XTick',[-10:10:40])
set(a(2,3),'Xlim',[-pre_sec post_sec],'Ylim',[0 30],'TickDir','out','Fontsize',7,'TickLength',[0.01 0.01],'box','off','FontName','Arial','XTickLabel',{},'YTickLabel',{},'XTick',[-10:10:40])
set(a(3,1),'Xlim',[-pre_sec post_sec],'Ylim',[0 35],'TickDir','out','Fontsize',7,'TickLength',[0.01 0.01],'box','off','FontName','Arial','XTick',[-10:10:40])
set(a(3,2),'Xlim',[-pre_sec post_sec],'Ylim',[0 35],'TickDir','out','Fontsize',7,'TickLength',[0.01 0.01],'box','off','FontName','Arial','YTickLabel',{},'XTick',[-10:10:40])
set(a(3,3),'Xlim',[-pre_sec post_sec],'Ylim',[0 35],'TickDir','out','Fontsize',7,'TickLength',[0.01 0.01],'box','off','FontName','Arial','YTickLabel',{},'XTick',[-10:10:40])
set(a(2,4),'Xlim',[0.5 3.5],'Ylim',[0 15],'TickDir','out','Fontsize',7,'TickLength',[0.01 0.01],'box','off','FontName','Arial')
% set(a(3,4),'Xlim',[0.5 3.5],'Ylim',[0 27],'TickDir','out','Fontsize',7,'TickLength',[0.01 0.01],'box','off','FontName','Arial')

%% Export

export_fig('pi-stability_in-vivo_tau_01','-pdf','-transparent',gcf)
saveas(gcf,'pi-stability_in-vivo_tau_01')

