%model_RTsplit_demo Simulated population activity as a function of RT
%
% This script computes population activity aligned on the go and and
% on saccade onset for responses generated by model1DR_2pop_demo or
% model1DR_3pop_demo. It compares slow (long-RT) vs fast (short-RT)
% responses over a range of RTs.
%
% This script should be run immediately after model1DR_2pop_demo or
% model1DR_3pop_demo, as it uses arrays generated by those scripts,
% with one trial per row:
%
%         RT --> reaction times
%        cho --> choices, in (1, correct) or out (0, error)
%   allrg_in --> cue driven neural responses aligned on go
%  allrg_out --> bias driven neural responses aligned on go
%   allrs_in --> cue driven neural responses aligned on saccade
%  allrs_out --> bias driven neural responses aligned on saccade
%      
% plus similar arrays with _in_cnt or _in_var extensions.
%  
% The parameter which_in determines which IN responses to analyze,
% those of the constant (baseline-independent) component, those of the
% variable (baseline-dependent) component, or both together. Only the
% third option is available after running model1DR_2pop_demo.

% Emilio Salinas, November 2017

%-----------
% parameters
%-----------

% which IN data: 'all', 'cnt', or 'var'
%which_in = 'var';   
%which_in = 'cnt';     
which_in = 'all';      

Nbin = 20;           % number of RT bins for splitting the data
qRT_ovlp = 0.10;     % {0--0.5}, half ovelap in RT quantiles

% data arrays
Ngo = length(ttg);
Nsac = length(tts);
Ntq = NaN(Nbin,1);
RTq = NaN(Nbin,1);
RTq_lo = NaN(Nbin,1);
RTq_hi = NaN(Nbin,1);

% RT bin centers for splitting the data
qRT_bin = [qRT_ovlp:(1-2*qRT_ovlp)/(Nbin-1):(1-qRT_ovlp)];   

% target output arrays: averages over RT quantiles
rpop_cin_g = NaN(Nbin,Ngo);
rpop_cin_s = NaN(Nbin,Nsac);
rpop_cout_g = NaN(Nbin,Ngo);
rpop_cout_s = NaN(Nbin,Nsac);

Ntq_c = zeros(Nbin,1);
Ntq_e = zeros(Nbin,1);

% report which data are being analyzed
if strcmpi(which_in, 'all')
	disp('Complete responses considered (r_cnt, r_var combined)')
elseif strcmpi(which_in, 'cnt')
	disp('Only constant-build-up responses considered (r_cnt)')
elseif strcmpi(which_in, 'var')
	disp('Only variable-build-up responses considered (r_var)')
else
	error('Invalid which_in value ({''all'', ''cnt'', ''var''})')
end

%----------------------------------
% sort population data into RT bins
%----------------------------------
%
% firing rate vs time (i.e., spike density functions) for each of the
% Nbin RT bins in each condition
%
iic = (cho > 0);
iie = (cho < 0);
for j=1:Nbin
	% split current dataset into RT quantiles
	qlo = qRT_bin(j) - qRT_ovlp;
	qhi = qRT_bin(j) + qRT_ovlp;
	if qlo < 0
		qlo = 0;
	end
	if qhi > 1
		qhi = 1;
	end
	q = quantile(RT, [qlo qhi]);
	ii = (RT >= q(1)) & (RT <= q(2)); 

	% save RT value representative (~ center) of current bin
	RTq(j) = median(RT(ii));

	% save edges of RT quantiles
	RTq_lo(j) = q(1);
	RTq_hi(j) = q(2);

	% number of trials in each quantile
	Ntq(j) = sum(ii);

	% find correct and error responses at given RT range
	iic1 = (ii & iic);
	iie1 = (ii & iie);
	nc1 = sum(iic1);
	ne1 = sum(iie1);

	% create average traces aligned on go
	if nc1
		Ntq_c(j) = nc1;
		sqn = sqrt(nc1);
		if strcmpi(which_in, 'all')
			rpop_cin_g(j,:) = mean(allrg_in(iic1,:));
		elseif strcmpi(which_in, 'cnt')
			rpop_cin_g(j,:) = mean(allrg_in_cnt(iic1,:));
		elseif strcmpi(which_in, 'var')
			rpop_cin_g(j,:) = mean(allrg_in_var(iic1,:));
		end
		rpop_cout_g(j,:) = mean(allrg_out(iic1,:));
	end

	% create average traces aligned on saccade
	if nc1
		sqn = sqrt(nc1);
		if strcmpi(which_in, 'all')
			rpop_cin_s(j,:) = mean(allrs_in(iic1,:));
		elseif strcmpi(which_in, 'cnt')
			rpop_cin_s(j,:) = mean(allrs_in_cnt(iic1,:));
		elseif strcmpi(which_in, 'var')
			rpop_cin_s(j,:) = mean(allrs_in_var(iic1,:));
		end
		rpop_cout_s(j,:) = mean(allrs_out(iic1,:));
	end

	% just to track progress
	disp(['quantile ' num2str(j) ', RTq = ' num2str(RTq(j)) ...
		  ' ms, nc1 = ' num2str(nc1) ', ne1 = ' num2str(ne1)])
end

%-----------------
% plot the results
%-----------------

% graph parameters
xlim1 = [-200 400];  % for data aligned on go
xlim2 = [-400 200];  % for data aligned on saccade

fnt1 = 8;
line1 = 1;
msz1 = 5;
tickfac = 2.0;
xtick = [-1000:200:1000];

% condition labels
cond_str = {['Correct T (in) activity' char(10) 'target at T, saccade to T']; ...   
            ['Correct D (out) activity' char(10) 'target at T, saccade to T']; ...
            ['Incorrect T (in) activity' char(10) 'target at T, saccade to D']; ...
            ['Incorrect D (out) activity' char(10) 'target at T, saccade to D']; ...
           };

% colors for the Nbin curves sorted by RT
Nbin = size(rpop_cin_g,1);        
if Nbin == 20
    % with 20 bins, use colormap used in the paper
	clm = [     0         0         0; ...
		   0.1080    0.1080    0.1088; ...
		   0.2242    0.2242    0.2260; ...
		   0.3322    0.3322    0.3348; ...
		   0.4484    0.4484    0.4520; ...
		   0.5564    0.5564    0.5608; ...
		   0.6726    0.6726    0.6779; ...
		   0.7806    0.7806    0.7867; ...
		   0.8885    0.8885    0.8955; ...
		   0.9882    0.9479    0.7507; ...
		   0.9822    0.9106    0.5608; ...
		   0.9756    0.8705    0.3563; ...
		   0.9696    0.8332    0.1665; ...
		   0.9664    0.7636    0.1263; ...
		   0.9636    0.6823    0.1123; ...
		   0.9702    0.5142    0.0845; ...
		   0.9856    0.2476    0.0407; ...
		   1.0000         0         0; ...
		   0.7300         0         0; ...
		   0.6000         0         0];
else
    % otherwise, use a built-in colormap
	map = colormap('jet');
	foo = size(map,1);
	ii = round([1:(foo-1)/(Nbin-1):foo]);
	clm = map(ii,:);
end

% reorder the color map entries for Go-aligned data
foo = round(Nbin/2);
jmap_g=foo;
for j=1:foo
	if foo+j <= Nbin
		jmap_g = [jmap_g foo+j];
	end
	if foo-j >= 1
		jmap_g = [jmap_g foo-j];
	end
end

% reorder the color map entries for Sac-aligned data
foo = round(Nbin/2);
jmap_s=foo;
for j=1:foo
	if foo+j <= Nbin
		jmap_s = [jmap_s foo+j];
	end
	if foo-j >= 1
		jmap_s = [jmap_s foo-j];
	end
end

ymaxx_g = max(max(rpop_cin_g));
ymaxx_s = max(max(rpop_cin_s));
ymaxx = max([ymaxx_g ymaxx_s]);
ylim1 = [0 ymaxx*1.10];

dy = diff(ylim1);
ymark = ylim1(1) - 0.07*dy;
ylab = ymark - 0.08*dy;

%
% initialize figure to white background
%
clf
whitebg('white') 
set(gcf,'color', 'white')

%
% plot fast-to-slow responses in 2 conditions (IN/OUT correct)
%
for jj=1:2
    switch jj
      case 1
        h1 = subplot(6,6,[8 9]);
        h2 = subplot(6,6,[10 11]);
        rpop_g = rpop_cin_g;
        rpop_s = rpop_cin_s;
      case 2
        h1 = subplot(6,6,[14 15]);
        h2 = subplot(6,6,[16 17]);
        rpop_g = rpop_cout_g;
        rpop_s = rpop_cout_s;
      otherwise
        error('Invalid condition ({1--2})')
    end
    mssg = cond_str{jj};

    % fast-to-slow responses aligned on go
	axes(h1)
    ticklen = get(h1, 'TickLength');
    ticklen = ticklen*tickfac;
	set(h1, 'Xtick', xtick, 'XtickLabel', [], 'TickDir', 'out', ...
        'Ycolor', 'w', 'FontSize', fnt1, 'TickLength', ticklen)

	hold on
	xlim(xlim1);
	ylim(ylim1);
    % firing rate traces
	for j=jmap_g
		plot(ttg, rpop_g(j,:), '-', 'linewidth', line1, 'color', clm(j,:))
	end
	plot([0 0], ylim, 'k-')
    plot(0, ymark, '^k', 'markersize', msz1, 'markerfacecolor', 'k', 'clipping', 'off')
    text(0, ylab, 'Go', 'VerticalAlignment', 'cap', ...
         'HorizontalAlignment', 'center', 'Fontsize', fnt1)

    % calibration bars
    plot((xlim1(1)-30)+[-20 0 0 -20], [0 0 1000 1000], 'k-', 'clipping', 'off')
    text(xlim1(1)-70, 0, '0', 'HorizontalAlignment', 'right', 'Fontsize', fnt1)
    text(xlim1(1)-70, 1000, '1000', 'HorizontalAlignment', 'right', 'Fontsize', fnt1)
    if jj ==  1
        ybar1 = ylim1(1) - diff(ylim1)*0.12;
        plot(xlim1(2)+[-200 0], ybar1*[1 1], ...
             'k-', 'linewidth', 2.5, 'clipping', 'off')
        ybar1 = ybar1 - diff(ylim1)*0.06;
        text(xlim1(2)-200/2, ybar1, '200 ms', 'FontSize', fnt1, ...
             'HorizontalAlignment', 'center', 'VerticalAlignment', 'top')
    end

    % fast-to-slow responses aligned on saccade
	axes(h2)
    ticklen = get(h2, 'TickLength');
    ticklen = ticklen*tickfac;
	set(h2, 'Xtick', xtick, 'XtickLabel', [], 'TickDir', 'out', 'Ycolor', 'w', ...
        'YtickLabel', [], 'FontSize', fnt1, 'TickLength', ticklen)

	hold on
	xlim(xlim2);
	ylim(ylim1);
    % firing rate traces
	for j=jmap_s
		plot(tts, rpop_s(j,:), '-', 'linewidth', line1, 'color', clm(j,:))
	end
	plot([0 0], ylim, 'k-')
    plot(0, ymark, '^k', 'markersize', msz1, 'markerfacecolor', 'k', 'clipping', 'off')
    text(0, ylab, 'Saccade', 'VerticalAlignment', 'cap', ...
         'HorizontalAlignment', 'center', 'Fontsize', fnt1)
    % condition label
	text(xlim2(2)-50, ylim1(2), mssg, 'Fontsize', fnt1, ...
         'HorizontalAlignment', 'left', 'VerticalAlignment', 'top')
end




