%% LATER demo: rise-to-threshold with Gaussian rate -> recinormal RTs
% Inspired by Fig. 4 in Noorani & Carpenter (2016)  [oai_citation:2‡Noorani, Carpenter - 2016 - The LATER model of reaction time and decision.pdf](sediment://file_00000000a3b0720cbd9174ef9accd180)

clear; 
close all; 
clc;

col = cbrewer('qual','Paired',6);
%% Parameters (edit these)
N       = 2000;        % number of trials
S0      = 0;           % starting level
Delta   = 1;           % threshold range (ST - S0)
ST      = S0 + Delta;  % threshold

mu_r    = 4.0;         % mean rate (units of S per second)
sigma_r = 0.8;         % SD of rate

r_lo = max(1e-6, norminv(0.025, mu_r, sigma_r));
r_hi = max(1e-6, norminv(0.975, mu_r, sigma_r));

r_lo = max(1e-6, norminv(0.025, mu_r, sigma_r));
r_hi = max(1e-6, norminv(0.975, mu_r, sigma_r));

tmax    = 0.8;         % seconds (for plotting decision trajectories)
dt      = 0.001;       % seconds (1 ms resolution for plotting)

fixedDelay = 0.00;     % optional non-decision delay (seconds), set e.g. 0.05

rng(1);                % reproducible

t		= -0.2:dt:tmax;

%% stimulus 
stimulus = zeros(size(t));
stimulus(t>0) = 1;

%% Average & 95% CI
    % r = r(idx(k));
S_mu = max(S0,min(ST, S0 + mu_r .* t));
S_lo = max(S0,min(ST, S0 + r_lo .* t));
S_hi = max(S0,min(ST, S0 + r_hi .* t));

% S_mu = max(S0,S0 + mu_r .* t);
% S_lo = max(S0,S0 + r_lo .* t);
% S_hi = max(S0,S0 + r_hi .* t);

%% Draw random rates and enforce r>0 (so every trial reaches threshold)
% LATER permits r<=0 (omissions), but for a clean RT histogram you may prefer truncation.
r = mu_r + sigma_r.*randn(N,1);

% Option A (default): truncate at small positive value
rmin = 1e-3;
r(r < rmin) = rmin;

% Option B: if you prefer to *drop* nonresponses instead, uncomment:
% r = r(r > 0);
% N = numel(r);

%% Compute reaction times
T_decision	= Delta ./ r;             % seconds
RT			= fixedDelay + T_decision;        % include optional fixed delay

%% Build example decision trajectories for a subset of trials
nShow	= 20;
idx		= round(linspace(1, N, nShow));

S		= nan(nShow, numel(t));
tCross	= nan(nShow,1);

for k = 1:nShow
    rk = r(idx(k));
    Sk = S0 + rk.*t;                 % linear rise from stimulus onset
	Sk(Sk<S0)=S0;
    % clip at threshold for plotting
    crossIdx = find(Sk >= ST, 1, 'first');
    if isempty(crossIdx)
        crossIdx = numel(t);
    end
    tCross(k) = t(crossIdx);
    Sk(crossIdx:end) = ST;
    S(k,:) = Sk;
end



%% Plot (Fig. 4-like): decision signals + RT histogram
figure(1);
clf

% --- Left: decision signals
subplot(3,1,2);

plot([-250 tmax*1000],[S0 S0],'k-');
text(-255,S0,'S0','Color','k','FontSize',24,'HorizontalAlignment','right','VerticalAlignment','middle');
hold on
plot([-250 tmax*1000],[ST ST],'k-');
text(-255,ST,'S_T','Color','k','FontSize',24,'HorizontalAlignment','right','VerticalAlignment','middle');

% plot(t*1000, S', 'LineWidth', 1,'Color',[0.7 0.7 0.7]); hold on; % time in ms
x_ms = t*1000;
patch([x_ms, fliplr(x_ms)], [S_lo, fliplr(S_hi)], 0.5*[1 1 1], ...
    'EdgeColor','none', 'FaceAlpha',1,'FaceColor',col(1,:));

hold on
plot(t*1000, S_mu, 'k-','LineWidth', 3,'Color',col(2,:)); hold on; % time in ms



% yline(ST,'k--','Threshold');
% xline(0,'k-','Stimulus onset');
xlabel('time (ms)');
ylabel('decision signal S');
% title('LATER: linear rise with Gaussian rate r');
xlim([min(t)*1000-100 tmax*1000]);
ylim([S0-0.1, ST+0.3]);

% Add a few markers where threshold is hit (decision time only)
% scatter(tCross*1000, ST*ones(size(tCross)), 18, 'k', 'filled', 'MarkerFaceAlpha', 0.6);
% text(0.02*tmax*1000, ST-0.08, sprintf('\\Delta = %.3g,  r \\sim N(%.3g, %.3g^2)',Delta,mu_r,sigma_r));
% nicegraph;
% axis normal;
text(-100,0.2,'decision signal','Color',col(2,:),'FontSize',24,'HorizontalAlignment','center');

xline(0,'k--');
axis off;




% --- Right: RT histogram (skewed)
% edges = linspace(0, max(RT)*1000, 40);
% histogram(RT*1000, edges, 'EdgeColor','none'); % ms
% xlabel('Reaction time (ms)');
% ylabel('Count');
% % title('Resulting RT distribution (recinormal, skewed)');
% box off;

% % Overlay a quick sanity check: promptness should be ~normal
% % (1/RT is not strictly normal if fixedDelay>0; use decision component)
% promptness = 1 ./ T_decision; % = r/Delta
% mu_p = mean(promptness);
% sd_p = std(promptness);
% yl = ylim;
% text(0.05*max(RT)*1000, 0.92*yl(2), sprintf('Promptness 1/T ~ N(%.3g, %.3g^2)',mu_p,sd_p));


% RT in ms
subplot(311)
RTms = RT(:) * 1000;

[xi, f] = ksdensity(RTms, 'Support','positive');   % xi: ms, f: density (1/ms)
xi = xi./max(xi);
hold on;
patch([zeros(size(f)) fliplr(f)], [xi fliplr(xi)], 0.5*[1 1 1], ...
    'EdgeColor','none', 'FaceAlpha',0.35, 'FaceColor',col(3,:));
plot(f, xi, 'k-', 'LineWidth', 1.5,'Color',col(4,:));
plot([-200,median(RTms)],[0 0],'LineWidth', 2,'Color',col(4,:));
plot([median(RTms) median(RTms)],[0 0.5],'LineWidth', 2,'Color',col(4,:));
plot([median(RTms) tmax*1000],[0.5 0.5],'LineWidth', 2,'Color',col(4,:));

xlabel('Reaction time (ms)');
ylabel('Density (1/ms)');
% title('RT distribution (KDE)');
box off;
xlim([min(t)*1000-100 tmax*1000]);
ylim([S0-0.3, ST+0.1]);
text(-100,0.2,'response','Color',col(4,:),'FontSize',24,'HorizontalAlignment','center');
xline(0,'k--');
plot([0 median(RTms)],[-0.2 -0.2],'LineWidth', 1,'Color',col(4,:));
text(median(RTms)/2,-0.1,'latency','Color',col(4,:),'FontSize',24,'HorizontalAlignment','center');

axis off;

% nicegraph;
% axis normal;

%% Plot stimulus
subplot(3,1,3);
plot(t*1000,stimulus,'k-','LineWidth',2,'Color',col(6,:))
ylim([S0-0.1, ST+0.1]);
box off;
xlim([min(t)*1000-100 tmax*1000]);
text(-100,0.2,'stimulus','Color',col(6,:),'FontSize',24,'HorizontalAlignment','center');
xline(0,'k--');
axis off;

%% Optional: reciprobit-style diagnostic (promptness normality)
% Uncomment if you want a quick check:
% figure('Color','w');
% z = sort(promptness);
% p = ((1:N) - 0.5) ./ N;
% % probit transform: inverse normal CDF
% y = sqrt(2)*erfinv(2*p - 1);
% plot(z, y, '.'); grid on;
% xlabel('Promptness 1/T (s^{-1})');
% ylabel('Probit (Z)');
% title('Reciprobit idea: if promptness is normal, points fall on a line');



% optional, save figure
fname = fullfile('..','..','images',[mfilename '_ori.svg']);
savegraph(fname,'svg');