%model1DR_3pop_demo Three populations rising to threshold in the 1DR task
%
% This is exactly like model1DR_2pop_demo except that the cue-driven
% activity is split into two parts (or subpopulations) having constant
% (or, more precisely baseline-independent) versus variable (i.e.,
% baseline-dependent) build-up rates.  
%
% The extensions _cnt and _var identify the two components here (as
% opposed to the superscripts U and C used in the manuscript, for
% 'uncoupled' and 'coupled' components). They can be mixed in
% different proportions via the parameter pmix; their means can be
% offset with parameter omix_* (which correspond to alpha in the
% manuscript); and the random noise can be assigned to them in
% different proportions via parameter nmix_*.
%
% The sensitivity to RT of the simulated responses can be seen by
% running model_RTsplit_demo right after running this script. See its
% header for further details.
%
% A couple of things to try:
%
%  (1) To enhance the RT sensitivity, especially for the
%  fast-preferring neurons, compare the results generated by 
%  run model_RTsplit_demo with which_in = 'var' in two conditions,
%  after running this script with
%
%      omix_io = -3.16
%      omix_oi = 0
% 
%  (the default below) versus after running this script with
%
%      omix_io = -3.16 - 0.5;
%      omix_oi = 0 - 0.5;
% 
%  This is equivalent to changing the alpha parameter in the manuscript.
%
%  (2) To see how strongly the intrinsic noise affects the RT
%  sensitivity in the high-reward condition, set which_condition = 1.
%  Then compare the results generated by model_RTplit_demo after
%  running this script with 
%  
%     nmix_oi = 0
%     nmix_io = 0 
%
%  (the default below) versus after running this script with 
%
%     nmix_oi = 1
%     nmix_io = 1 
%
% This parameter change shifts the intrinsic randomness entirely from
% one component to the other. You'll note that the fast/slow
% preferences reverse.

% Emilio Salinas, November 2017

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

which_condition = 2;   % {1--3} --> hi in/lo out, lo in/hi out, ADR
do_traces = 0;         % {0--4} show single-trial neural responses
smix = 0.0;            % {0--1} --> proportion of _var slowdown 
pmix = 0.0;            % {0--0.5} --> proportions of _cnt and _var

Nit = 10000;       % number of simulated trials
pst_tau = 120;     % time constant for postsaccadic tail (ms)
pst_ydip = 200;    % asymptote level for postsaccadic tail
wbin = 20;         % bin width for RT histograms
Tmax = 1000;       % max trial length
RTmax = 520;       % max RT for calculating RT statistics
prec = 3;          % precision for numerical reports

%
% race parameters that are fixed across 1DR/ADR conditions
%
% scale factor for threshold
th = 1000;                  
% threshold coefficients
th_c = [1.185 1.2 -1.2 0.723];    

mTi1 = 40;        % mean interrupt start
mTi2 = 155;       % mean interrupt end
sdTi = 2.0*0;     % SD of Ti1 and Ti2
vTi_fac = 0.38;   % attenuation factor during interrupt

a_in = 35;        % mean afferent delay of cue
y_in = 0.10*0;    % SD of afferent delay of cue
a_out = 50;       % mean afferent delay of internal plan
y_out = 0.05*0;   % SD of afferent delay of internal plan
rho_bb = -0.5;    % correlation between baseline fluctuations

lose_out_fac = 0;     % gain drop when bias-driven plan loses
win_out_fac = 1.0;    % gain boost when bias-driven plan wins
lose_in_fac = 1.0;    % gain drop when cue-driven plan loses 
win_in_c1 = 5.5;      % neutral value of cue-driven build-up rate
win_in_c2 = 2.6;      % gain boost/drop when cue-driven plan wins

% build-up rate parameters when baseline_in > baseline_out 
% (initial condition with no conflict)
% here _out gain is irrelevant because activity doesn't rise
gc_in_c1 = 6.16;   
gc_in_c2 = 2.50;  
gc_in_c3 = 0;  
SDgc_in_c1 = 0.55;
omix_io = -3.16;   
nmix_io = 0.0;

% build-up rate parameters when baseline_in < baseline_out 
% (initial condition with conflict)
gc_in_c4 = 3.00;
gc_in_c5 = 23.25;  
gc_in_c6 = 1.30;  
SDgc_in_c4 = 0.30;
gc_out_c1 = 1.40; 
gc_out_c2 = 1.70;             
SDgc_out_c1 = 0;
omix_oi = 0;   
nmix_oi = 0.0;

% 
% race parameters that vary across 1DR/ADR conditions;
% i.e., mainly the two baseline levels
% 
if which_condition == 1
    % for congruent (high-reward) condition
	p_in = 0.28; 
	p_out = 0.28; 
	b_in = th*0.340; 
	b_out = th*0.16; 
elseif which_condition == 2
    % for incongruent (low-reward) condition
	p_in = 0.28; 
	p_out = 0.28; 
	b_in = th*0.16; 
	b_out = th*0.34; 
elseif which_condition == 3
    % for ADR condition
	p_in = 0.14; 
	p_out = 0.14; 
	b_in = th*0.20; 
	b_out = th*0.20; 
    % to capture the ADR condition just right, a small adjustment
    % is necessary in this parameter
    gc_in_c1 = 5.95;
else
    error('Invalid which_condition value ({1,2,3})')
end

% report condition to screen (useful while waiting)
if which_condition == 1
	disp('Simulating congruent (high-reward) condition')
elseif which_condition == 2
	disp('Simulating incongruent (low-reward) condition')
elseif which_condition == 3
	disp('Simulating ADR condition')
end

% minimum firing rate values allowed before saccade
rmin1 = 0;
rmin2 = 0;

%------  ----------------------------------------------------
% compute baselines, gain, threshold, and afferent delays for 
% each trial
%------------------------------------------------------------
rr = mvnrnd([0 0], [1 rho_bb; rho_bb 1], Nit);
csi_in = rr(:,1);
csi_out = rr(:,2);
eps_in = randn(Nit,1);
eps_out = randn(Nit,1);
phi_in = randn(Nit,1);
phi_out = randn(Nit,1);

% afferent delays
A_in = a_in*(1 + y_in*phi_in);
A_in = round(A_in);
A_out = a_out*(1 + y_out*phi_out);
A_out = round(A_out);

% baseline levels
B_in = b_in*(1 + p_in*csi_in);
B_in(B_in <= 0) = 0;
B_out = b_out*(1 + p_out*csi_out);
B_out(B_out <= 0) = 0;

% gain factors (i.e., build-up rates), which increase with the
% baseline level of the corresponding motor plan and decrease with the
% baseline level of the competing plan
noise_out = SDgc_out_c1*randn(Nit,1);
G_out = gc_out_c1 + noise_out + gc_out_c2*(B_out - B_in)/th;
G_out(G_out < 0) = 0;
% case 1: when baseline_in > baseline_out
noise_in = SDgc_in_c1*randn(Nit,1);
n_cnt = (1 - nmix_io)*noise_in;
n_var = nmix_io*noise_in;
G_in_cnt1_io = 2*(omix_io + n_cnt + gc_in_c1)./(1 + gc_in_c3*B_out/th);
G_in_var1_io = 2*(-omix_io + n_var + gc_in_c2*B_in/th)./(1 + gc_in_c3*B_out/th);
% case 2: when baseline_out > baseline_in
noise_in = SDgc_in_c4*randn(Nit,1);
n_cnt = (1 - nmix_oi)*noise_in;
n_var = nmix_oi*noise_in;
G_in_cnt1_oi = 2*(omix_oi + n_cnt + gc_in_c4)./(1 + gc_in_c6*B_out/th);
G_in_var1_oi = 2*(-omix_oi + n_var + gc_in_c5*B_in/th)./(1 + gc_in_c6*B_out/th);
% sort cases 1 and 2 and assign appropriate build-up rates
ii = find(B_in >= B_out);
G_in_cnt1 = G_in_cnt1_oi;
G_in_cnt1(ii) = G_in_cnt1_io(ii);
G_in_var1 = G_in_var1_oi;
G_in_var1(ii) = G_in_var1_io(ii);
G_in_cnt1(G_in_cnt1 < 0) = 0;
G_in_var1(G_in_var1 < 0) = 0;
% mix constant and variable terms in desired proportions
G_in_cnt = (1 - pmix)*G_in_cnt1 + pmix*G_in_var1;
G_in_var = (1 - pmix)*G_in_var1 + pmix*G_in_cnt1;

% sanity check
G_in = (G_in_cnt + G_in_var)/2;
if any(G_in_cnt < 0) | any(G_in_var < 0) | any(G_in < 0) | any(G_out < 0) 
    error('Build-up rates (G_*) cannot be negative')
end

%
% threshold level in each trial; fluctuations depend on the baseline
% levels
%
th1 = th_c(1)*th + th_c(2)*B_in + th_c(3)*B_out;

% the threshold can't be lower than a certain level
ii = find(th1 < th_c(4)*th);
th1(ii) = th_c(4)*th;

% interrupt times
if mTi2 > mTi1
    dTi = round(sdTi*randn(Nit,1));
    Ti1 = mTi1 + dTi;
    Ti2 = mTi2 + dTi;
else
    Ti1 = zeros(Nit,1);
    Ti2 = zeros(Nit,1);
end

%----------------------------------------------
% generate neural responses (rise to threshold)
%----------------------------------------------
% time bases for data aligned on go and saccade
ttg = [-200:Tmax];
Nttg = length(ttg);
igo = find(ttg == 0);
tts = [(-Tmax-200):200];
Ntts = length(tts);
iiz = find(tts == 0);
Ntail1 = Ntts - iiz;
absmix = abs(smix);

% for mean response traces (correct and error) aligned on go
rg_cin_cnt = zeros(1,Nttg);
rg_cin_var = zeros(1,Nttg);
rg_cin = zeros(1,Nttg);
rg_cout = zeros(1,Nttg);
rg_ein_cnt = zeros(1,Nttg);
rg_ein_var = zeros(1,Nttg);
rg_ein = zeros(1,Nttg);
rg_eout = zeros(1,Nttg);

% for mean response traces (correct and error) aligned on saccade
rs_cin_cnt = zeros(1,Ntts);
rs_cin_var = zeros(1,Ntts);
rs_cin = zeros(1,Ntts);
rs_cout = zeros(1,Ntts);
rs_ein_cnt = zeros(1,Ntts);
rs_ein_var = zeros(1,Ntts);
rs_ein = zeros(1,Ntts);
rs_eout = zeros(1,Ntts);

% for temporary responses
r1_cnt = NaN(1,Nttg);
r1_var = NaN(1,Nttg);
r1 = NaN(1,Nttg);
r2 = NaN(1,Nttg);
r3_cnt = NaN(1,Ntts);
r3_var = NaN(1,Ntts);
r3 = NaN(1,Ntts);
r4 = NaN(1,Ntts);

% to save individual neural trajectories
allrg_in_cnt = NaN(Nit,Nttg);
allrg_in_var = NaN(Nit,Nttg);
allrg_in = NaN(Nit,Nttg);
allrg_out = NaN(Nit,Nttg);
allrs_in_cnt = NaN(Nit,Ntts);
allrs_in_var = NaN(Nit,Ntts);
allrs_in = NaN(Nit,Ntts);
allrs_out = NaN(Nit,Ntts);

% trial outcomes
rwin = NaN(Nit,1);    % winner firing level @ saccade
rlose = NaN(Nit,1);   % loser firing level @ saccade
cho = NaN(Nit,1);     % choice 
RT = NaN(Nit,1);      % RT  

% falling tail for postsaccadic responses
ii = [0:1200];
fall1 = exp(-ii/pst_tau);

% in case individual trial graphs are requested
if do_traces
   clf
   ylo = -100;
   yhi = ceil((mean(th1) + 2*std(th1))/50)*50;
   clg = 0.5*[1 1 1];
end

if do_traces 
	if do_traces == 1
		disp('showing neural trajectories in valid trials')
	elseif do_traces == 2
		disp('showing neural trajectories in error trials')
	elseif do_traces == 3
		disp('showing neural trajectories in correct trials')
	elseif do_traces == 4
		disp('showing neural trajectories in no-movement trials')
	else
		error('Unknown do_traces option ({0--4})')
	end
	disp('press ctrl-C to stop...')
end

%
% generate neural trajectories over time
%
% loop over trials
for j=1:Nit
    % initial conditions
    rt1_cnt = B_in(j);
    rt1_var = B_in(j);
    rt1 = (rt1_cnt + rt1_var)/2;
    rt2 = B_out(j);
    r1_cnt(1:igo) = rt1_cnt;
    r1_var(1:igo) = rt1_var;
    r2(1:igo) = rt2;
    v1_cnt = G_in_cnt(j);
    v1_var = G_in_var(j);
    v2 = G_out(j);
    vv1_cnt = 0;
    vv1_var = 0;
    vv2 = 0;
    won_flag = 0;
    % loop over time, for each trial
    for t=1:Tmax
        % responses allowed to rise only after afferent delay
        if t < A_in(j)
            vv1_cnt = 0;
            vv1_var = 0;
        end
        if t < A_out(j)
            vv2 = 0;
			intrpt = false;
        end

        % by default, IN response rises after afferent delay
        if t >= A_in(j)
             vv1_cnt = v1_cnt;
             vv1_var = v1_var;
             v1 = (v1_cnt + v1_var)/2;
        end

        % OUT response also rises after afferent delay...
        if t >= A_out(j)
			% ...but not at full rate during the cue-driven interruption 
			if (t > Ti1(j)) && (t <= Ti2(j)) && (won_flag == 0)
				vv2 = v2*vTi_fac;
				intrpt = true;
			else
				vv2 = v2;
				intrpt = false;
			end
        end

		% the plan that first goes past the other, wins
        if won_flag == 0
            % rt1 exceeds rt2 and wins the race
            if (t >= A_in(j)) && (rt1 > rt2)
				% adjust the target build-up rate of combined cue-driven 
                % population
                v1_avg = win_in_c1 + win_in_c2*(v1 - win_in_c1);
				% if the combined rate must decrease, reduce the variable 
				% and constant components in desired proportions, from 
				% attenuating only the _var component (smix=1) to
				% attenuating both equally (smix=0) to attenuating the
				% _cnt component only (smix=-1)
                if smix >= 0
					v1_var_1 = 2*win_in_c1 + win_in_c2*(v1_var - 2*win_in_c1) ...
							   + v1_cnt*(win_in_c2 - 1);
					v1_cnt_1 = v1_cnt;
                else
					v1_cnt_1 = 2*win_in_c1 + win_in_c2*(v1_cnt - 2*win_in_c1) ...
							   + v1_var*(win_in_c2 - 1);
					v1_var_1 = v1_var;
                end
				v1_var_2 = win_in_c1 + win_in_c2*(v1_var - win_in_c1);
				v1_cnt_2 = win_in_c1 + win_in_c2*(v1_cnt - win_in_c1);
				v1_var = v1_var_1*absmix + v1_var_2*(1 - absmix);
				v1_cnt = v1_cnt_1*absmix + v1_cnt_2*(1 - absmix);
				% the constant component, v1_cnt, may need to change a bit even 
                % for smax=1, if strong reduction of the combined build-up rate 
                % is necessary; the complementary case should be rare
				if v1_var < 0
					v1_var = 0;
					v1_cnt = 2*v1_avg; 
                elseif v1_cnt < 0
					v1_cnt = 0;
					v1_var = 2*v1_avg; 
				end
				% recompute combined rate, which should always be equal to 
                % v1_avg above (in case a check is needed)
                v1 = (v1_cnt + v1_var)/2;
                % adjust build-up rate of loser
			    v2 = v2*lose_out_fac;
                vv2 = v2;
                won_flag = 1;
            end

            % rt2 exceeds rt1 and wins the race
            if (t >= A_out(j)) && (intrpt == 0) && (rt2 > rt1)
                % adjust build-up rate of winner
                v2 = v2*win_out_fac;
                vv2 = v2;
                % adjust build-up rate of loser
                v1_cnt = v1_cnt*lose_in_fac;
                v1_var = v1_var*lose_in_fac;
                v1 = (v1_cnt + v1_var)/2;
                vv1_cnt = v1_cnt;
                vv1_var = v1_var;
                won_flag = 2;
            end
        end

        % update trajectories (assumes delta_t = 1 ms)
        rt1_cnt = rt1_cnt + vv1_cnt;  
        rt1_var = rt1_var + vv1_var;  
        rt1 = (rt1_cnt + rt1_var)/2;  
        rt2 = rt2 + vv2;  

        % prevent the loser from overtaking the winner once the race
        % has been decided
        if (won_flag == 1) && (rt2 > 0.98*rt1)
            % this should never happen
            rt2 = 0.98*rt1;
        end
        if (won_flag == 2) && (rt1 > 0.98*rt2)
            % this happens often
            dr = rt1 - 0.98*rt2;
            rt1_cnt = rt1_cnt - dr;
            rt1_var = rt1_var - dr;
            rt1 = (rt1_cnt + rt1_var)/2;  
        end

        % prevent firing rates from dipping below minimum (0)
        if rt1 < rmin1
            rt1_cnt = rmin1/2;
            rt1_var = rmin1/2;
            rt1 =  (rt1_cnt + rt1_var)/2;  
        end
        if rt2 < rmin2
            rt2 = rmin2;
        end

        % save trajectories
        r1_cnt(igo+t) = rt1_cnt;  
        r1_var(igo+t) = rt1_var;  
        r2(igo+t) = rt2;  

		% during the interrupt, the bias-driven plan cannot trigger a
		% saccade; if that's about to happen, threshold increases a
		% bit (this does not normally happen)
        if (rt2 > th1(j)) && (intrpt == 1)
            th1(j) = rt2 + 50 + 100*rand;
        end

        % stop if threshold is crossed
        if (rt1 > th1(j)) || (rt2 > th1(j)) 
            break
        end
    end
    % set full cue-driven activity 
    r1 = (r1_cnt + r1_var)/2;

    % set RT and choice
    if t < Tmax 
        RTj = t;
        % cho == NaN means no movement
	    if rt1 > rt2
			cho(j) = 1;
		elseif rt1 < rt2
			cho(j) = -1;
		else
			cho(j) = -1 + 2*(rand < 0.5);
		end
    else
        % RT == NaN means no movement
        RTj = NaN;
    end
    RT(j) = RTj;
    isac = igo + RTj;
        
	% append postsaccadic tails
	ii = find(ttg >= RTj);
	Ntail = length(ii);
	if Ntail
		r1_cnt(ii) = (r1_cnt(isac) - pst_ydip)*fall1(1:Ntail) + pst_ydip;
		r1_var(ii) = (r1_var(isac) - pst_ydip)*fall1(1:Ntail) + pst_ydip;
		r1(ii) = (r1(isac) - pst_ydip)*fall1(1:Ntail) + pst_ydip;
		r2(ii) = (r2(isac) - pst_ydip)*fall1(1:Ntail) + pst_ydip;
	end

	% generate saccade-aligned responses
	r3_cnt = r3_cnt*NaN;
	r3_var = r3_var*NaN;
	r3 = r3*NaN;
	r4 = r4*NaN;
    if t < Tmax 
	    ilo = iiz-isac+1;
        % inward
    	r3_cnt(1:ilo) = r1_cnt(1);
    	r3_cnt(ilo:iiz) = r1_cnt(1:isac);
    	r3_cnt(iiz+1:end) = (r3_cnt(iiz) - pst_ydip)*fall1(1:Ntail1) + pst_ydip;
    	r3_var(1:ilo) = r1_var(1);
    	r3_var(ilo:iiz) = r1_var(1:isac);
    	r3_var(iiz+1:end) = (r3_var(iiz) - pst_ydip)*fall1(1:Ntail1) + pst_ydip;
    	r3(1:ilo) = r1(1);
    	r3(ilo:iiz) = r1(1:isac);
    	r3(iiz+1:end) = (r3(iiz) - pst_ydip)*fall1(1:Ntail1) + pst_ydip;
        % outward
	    r4(1:ilo) = r2(1);
    	r4(ilo:iiz) = r2(1:isac);
    	r4(iiz+1:end) = (r4(iiz) - pst_ydip)*fall1(1:Ntail1) + pst_ydip;
    end

	% record firing levels at saccade, for a later sanity check
    if t < Tmax 
	    rwin(j) = max([r1(isac) r2(isac)]);
	    rlose(j) = min([r1(isac) r2(isac)]);
    end

	% show single-trial responses, if requested
	if ((do_traces == 1) && ~isnan(cho(j))) || ...
	   ((do_traces == 2) && (cho(j) < 0)) || ...
	   ((do_traces == 3) && (cho(j) > 0)) || ...
	   ((do_traces == 4) && (isnan(cho(j))))

        % responses aligned on go
		subplot(2,1,1)
        cl1 = get(gca, 'Xcolor');
		cla
		hold on
		if (mTi2 > mTi1) && (Ti2(j) > Ti1(j))
			rectangle('Position', [Ti1(j) ylo Ti2(j)-Ti1(j) yhi-ylo], ...
					  'facecolor', clg, 'edgecolor', clg)
		end
		plot([-200 600], th1(j)*[1 1], 'color', cl1)
		plot([0 0], [ylo yhi], '-', 'color', cl1)
		plot(RTj*[1 1], [ylo yhi], '--', 'color', cl1)
		plot(ttg,r1_cnt,'b-', ttg,r1_var,'c-', 'linewidth', 1)
		plot(ttg,r1,'r-', ttg,r2,'g-', 'linewidth', 1)
		xlim([-200 600])
		ylim([ylo yhi])
		xlabel('Time from go (ms)')
		ylabel('Neural activity')
		mssg = ['Trial ' num2str(j) ', RT = ' num2str(RTj) ' ms'];
		if cho(j) > 0
			mssg = [mssg ', correct'];
		elseif cho(j) < 0
            if RTj <= 1
			    mssg = [mssg ', fixation break'];
            else
			    mssg = [mssg ', incorrect'];
            end
		else
			mssg = [mssg ', no movement in required time'];
		end
		title(mssg)

        % responses aligned on saccade
		subplot(2,1,2)
		cla
		hold on
		plot([-600 200], th1(j)*[1 1], '-', 'color', cl1)
		plot([0 0], [ylo yhi], '-', 'color', cl1)
		plot(-RTj*[1 1], [ylo yhi], '--', 'color', cl1)
		plot(tts,r3_cnt,'b-', tts,r3_var,'c-', 'linewidth', 1)
		plot(tts,r3,'r-', tts,r4,'g-', 'linewidth', 1)
		xlim([-600 200])
		ylim([ylo yhi])
		ylabel('Neural activity')
		xlabel('Time from saccade (ms)')

		% report some useful data for debugging/checking
        %
        % trial number
		disp(['j = ' num2str(j)])
        % afferent delays
		disp(['  [A_in A_out] = ['  num2str(A_in(j),3) '  ' num2str(A_out(j),3) ']'])
        % interruption interval
		disp(['  [Ti1 Ti2] = [' num2str(Ti1(j),3) '  ' num2str(Ti2(j),3) ']'])
        % firing levels and threshold at threshold crossing
        if ~isnan(cho(j))
		    disp(['  [r1 r2 th1] = [' num2str(r1(isac),4) '  ' ...
		    	  num2str(r2(isac),4) '  ' num2str(th1(j),4) ']'])
        else
		    disp(['  [r1 r2 th1] = [xxx  xxx  ' num2str(th1(j),4) ']'])
        end

		pause
	end

	% sort in/out responses by correct/incorrect 
	%
	% note that *_ein and *_eout indicate the responses into and out
	% of the RF during errors, for which the saccade was *away*; so
	% during errors, eout is the winner and ein the loser
    if RTj > 1
		if cho(j) > 0
			rg_cin_cnt = rg_cin_cnt + r1_cnt;
			rg_cin_var = rg_cin_var + r1_var;
			rg_cin = rg_cin + r1;
			rg_cout = rg_cout + r2;
			rs_cin_cnt = rs_cin_cnt + r3_cnt;
			rs_cin_var = rs_cin_var + r3_var;
			rs_cin = rs_cin + r3;
			rs_cout = rs_cout + r4;
		elseif cho(j) < 0
			rg_ein_cnt = rg_ein_cnt + r1_cnt;
			rg_ein_var = rg_ein_var + r1_var;
			rg_ein = rg_ein + r1;
			rg_eout = rg_eout + r2;
			rs_ein_cnt = rs_ein_cnt + r3_cnt;
			rs_ein_var = rs_ein_var + r3_var;
			rs_ein = rs_ein + r3;
			rs_eout = rs_eout + r4;
		else
			% cho == NaN means no movement; 
			% RT == 1 means fixation break; 
            % exclude those cases from neural averages
		end
    end

	% save all individual response traces (for later analyses)
	allrg_in_cnt(j,:) = r1_cnt;
	allrg_in_var(j,:) = r1_var;
	allrg_in(j,:) = r1;
	allrg_out(j,:) = r2;
	allrs_in_cnt(j,:) = r3_cnt;
	allrs_in_var(j,:) = r3_var;
	allrs_in(j,:) = r3;
	allrs_out(j,:) = r4;

end

%---------------------------------------------
% compile statistics after all trials have run
%---------------------------------------------

% fractions correct, incorrect, fix-breaks, and no-movement
nnm = sum(isnan(cho));
fnm = nnm/Nit;
nbk = sum(RT <= 1);
fbk = nbk/Nit;
iok = (RT > 1);
nok = sum(iok);
nc = sum(cho(iok) > 0);
ne = sum(cho(iok) < 0);
nce = nc + ne;
if nce ~= nok
    % sanity check
    error('Mismatch in number of valid trials')
end
fc = nc/nce;
fe = 1 - fc;

% RT statistics for correct and incorrect choices 
%
% Note that RTmax is a separate cutoff for 'valid' RTs, to mimic the
% response time window (i.e., trial selection criterion) used during 
% the experiment
RTc = RT((cho > 0) & iok); 
RTe = RT((cho < 0) & iok); 

mRTc = mean(RTc(RTc <= RTmax));
mdRTc = median(RTc(RTc <= RTmax));
sdRTc = std(RTc(RTc <= RTmax));
skRTc = skewness(RTc(RTc <= RTmax));
mRTe = mean(RTe(RTe <= RTmax));
mdRTe = median(RTe(RTe <= RTmax));
sdRTe = std(RTe(RTe <= RTmax));
skRTe = skewness(RTe(RTe <= RTmax));

% report results
disp(['no movement trials: ' num2str(nnm) ' (' num2str(fnm, prec) ')'])
disp(['fixation breaks: ' num2str(nbk) ' (' num2str(fbk, prec) ')'])
disp(['competition between inward (cue-driven) and outward (bias-driven) activity'])

disp(['fraction correct (IN): ' num2str(fc, prec) ...
	  ' (' num2str(nc) '/' num2str(nce) ')'])
disp(['   valid RTs: ' num2str(mRTc, prec) ' +/- ' ...
	  num2str(sdRTc, prec) ',  sk = ' num2str(skRTc, prec) ...
	  ',  median = ' num2str(mdRTc, prec)])

disp(['fraction incorrect (OUT): ' num2str(fe, prec) ...
	  ' (' num2str(ne) '/' num2str(nce) ')'])
disp(['   valid RTs: ' num2str(mRTe, prec) ' +/- ' ...
	  num2str(sdRTe, prec) ',  sk = ' num2str(skRTe, prec) ...
	  ',  median = ' num2str(mdRTe, prec)])

%
% mean neural trajectories
%

% mean go-aligned responses
rg_cin_cnt = rg_cin_cnt/nc;
rg_cin_var = rg_cin_var/nc;
rg_cin = rg_cin/nc;
rg_cout = rg_cout/nc;
rg_ein_cnt = rg_ein_cnt/ne;
rg_ein_var = rg_ein_var/ne;
rg_ein = rg_ein/ne;
rg_eout = rg_eout/ne;

% mean saccade-aligned responses
rs_cin_cnt = rs_cin_cnt/nc;
rs_cin_var = rs_cin_var/nc;
rs_cin = rs_cin/nc;
rs_cout = rs_cout/nc;
rs_ein_cnt = rs_ein_cnt/ne;
rs_ein_var = rs_ein_var/ne;
rs_ein = rs_ein/ne;
rs_eout = rs_eout/ne;

% go and saccade values
iig = find(ttg == 0);
Rp_cin = rs_cin(iiz);
Rp_cout = rs_cout(iiz);
Rp_ein = rs_ein(iiz);
Rp_eout = rs_eout(iiz);
Rg_cin = rg_cin(iig);
Rg_cout = rg_cout(iig);
Rg_eout = rg_eout(iig);
Rg_ein = rg_ein(iig);

% normalize go and saccade values
%Rnorm = Rp_cin;
Rnorm = th;
Rp_cin = Rp_cin/Rnorm;
Rp_cout = Rp_cout/Rnorm;
Rp_ein = Rp_ein/Rnorm;
Rp_eout = Rp_eout/Rnorm;
Rg_cin = Rg_cin/Rnorm;
Rg_cout = Rg_cout/Rnorm;
Rg_eout = Rg_eout/Rnorm;
Rg_ein = Rg_ein/Rnorm;

% sanity check: winner response at saccade must be very close to
% threshold
th_err = abs(th1 - rwin);
th_err(isnan(cho)) = 0;
th_err(RT <= 1) = 0;
emax = max(th_err);

% report results
disp(['check: max(|th1 - rwin|) = ' num2str(emax,5)])
disp(['normalized go and saccade responses:'])
disp(['     correct in: ' num2str(Rg_cin, prec) ',  ' ...
						  num2str(Rp_cin, prec)])
disp(['    correct out: ' num2str(Rg_cout, prec) ',  ' ...
						  num2str(Rp_cout, prec)])
disp(['   incorrect in: ' num2str(Rg_ein, prec) ',  ' ...
						  num2str(Rp_ein, prec)])
disp(['  incorrect out: ' num2str(Rg_eout, prec) ',  ' ...
						  num2str(Rp_eout, prec)])

%----------------------
% compute RT histograms
%----------------------

% bin the data into 1 ms time bins and then smooth using a filter of
% width fwin
RTbin = [0:1:1000];
fwin = ones(wbin,1)/wbin;
hc = conv(hist(RTc, RTbin), fwin, 'same');
he = conv(hist(RTe, RTbin), fwin, 'same');

% normalize RT histograms
%
% {0-2} normalize RT histograms to sum (1) or max (2)
do_norm = 2;       
if do_norm == 1
	hnc = hc/sum(hc);
	hne = he/sum(he);
elseif do_norm == 2
	hnc = hc/max(hc);
	hne = he/max(he);
else
	hnc = hc;
	hne = he;
end

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

clf

% RT distributions
subplot(3,2,[1 2])
xlim([0 600])
hold on
hb = bar(RTbin, hnc, 'r');
set(hb, 'EdgeColor', 'r')
plot(RTbin, hne, 'c.-')
set(gca,'TickDir', 'out')
xlabel('RT (ms)')
title('Correct/incorrect RTs')

% mean neural responses
ylo = 0;
yhi = 1300;
cl1 = get(gca, 'Xcolor');
cl2 = 0.5*[1 1 1];

subplot(3,2,3)
xlim([-200 600])
ylim([ylo yhi])
hold on
plot(ttg, rg_cout, 'g.-')
plot(ttg, rg_cin, 'r.-')
plot(ttg, rg_cin_cnt, 'b-')
plot(ttg, rg_cin_var, 'c-')
plot(ttg, (rg_cin_cnt+rg_cin_var)/2, '-', 'color', cl2)
plot([0 0], ylim, '-', 'color', cl1)
plot(xlim, th*[1 1], ':', 'color', cl1)
set(gca,'TickDir', 'out')
xlabel('Time from go (ms)')
title('IN/OUT trajectories in correct trials')

subplot(3,2,4)
xlim([-600 200])
ylim([ylo yhi])
hold on
plot(tts, rs_cout, 'g.-')
plot(tts, rs_cin, 'r.-')
plot(tts, rs_cin_cnt, 'b-')
plot(tts, rs_cin_var, 'c-')
plot(tts, (rs_cin_cnt+rs_cin_var)/2, '-', 'color', cl2)
plot([0 0], ylim, '-', 'color', cl1)
plot(xlim, th*[1 1], ':', 'color', cl1)
set(gca,'TickDir', 'out')
xlabel('Time from saccade (ms)')

subplot(3,2,5)
xlim([-200 600])
ylim([ylo yhi])
hold on
plot(ttg, rg_eout, '.-', 'color', 0.7*[0 1 0])
plot(ttg, rg_ein, '.-', 'color', 0.7*[1 0 0])
plot(ttg, rg_ein_cnt, 'b-')
plot(ttg, rg_ein_var, 'c-')
plot(ttg, (rg_ein_cnt+rg_ein_var)/2, '-', 'color', cl2)
plot([0 0], ylim, '-', 'color', cl1)
plot(xlim, th*[1 1], ':', 'color', cl1)
set(gca,'TickDir', 'out')
xlabel('Time from go (ms)')
title('IN/OUT trajectories in incorrect trials')

subplot(3,2,6)
xlim([-600 200])
ylim([ylo yhi])
hold on
plot(tts, rs_eout, '.-', 'color', 0.7*[0 1 0])
plot(tts, rs_ein, '.-', 'color', 0.7*[1 0 0])
plot(tts, rs_ein_cnt, 'b-')
plot(tts, rs_ein_var, 'c-')
plot(tts, (rs_ein_cnt+rs_ein_var)/2, 'w-', 'color', cl2)
plot([0 0], ylim, '-', 'color', cl1)
plot(xlim, th*[1 1], ':', 'color', cl1)
set(gca,'TickDir', 'out')
xlabel('Time from saccade (ms)')

disp(' ')



