%% Simulate DDM with the given parameters

deltat = 0.05;
v = 0.4;
s = 0.2;
nvar = 4;
a = 0.5;
Nt = 1e3;
Nsamples = 1e5;
x = zeros(Nsamples, Nt);
rt = zeros(Nsamples, 1);
choice = zeros(Nsamples, 1);

for i = 1:Nsamples
    x(i, 1) = 0;
    
    for j = 2:Nt
        x(i, j) = x(i, j-1) + v*deltat + sqrt(deltat)*s*randn(1)*sqrt(nvar);
        
        if abs(x(i, j)) > a
            rt(i) = j;
            choice(i) = sign(x(i, j));
            break;
        end
        
        if j == Nt
            % time out
            rt(i) = inf;
            choice(i) = 0;
        end
    end
end

correct_ind = find(choice==1);
correct_ind = correct_ind(1);
incorrect_ind = find(choice==-1);
if isempty(incorrect_ind)
    warning('no incorrect choice')
else
    incorrect_ind = incorrect_ind(1);
end

m = max([rt(correct_ind), rt(incorrect_ind)]);
m = m + 10;

%% Plot two trajectories of the drift-diffusion process (one for each alternative) along with the decision boundaries


figure, plot([1 m], [a a], '-k'), hold on
plot([1 m], [-a -a], '-k')
axis([1 m -0.75 0.75])
plot(1:rt(correct_ind), x(correct_ind, 1:rt(correct_ind)), 'b')
plot(1:rt(incorrect_ind), x(incorrect_ind, 1:rt(incorrect_ind)), 'r')



%% Update a single trace dynamically until one of the processes reaches the threshold. 

figure, plot([1 m], [a a], '-k'), hold on
plot([1 m], [-a -a], '-k')
axis([1 m -0.75 0.75])
xAxis=1:m
for i = 1:m
    if(i<=rt(correct_ind))
        plot(xAxis(1:i), x(correct_ind, 1:i), 'b')
    end
    if(i<=rt(incorrect_ind))
        plot(xAxis(1:i), x(incorrect_ind, 1:i), 'r')
    end
    pause(0.5)
end