clear; clc; close all;

% let's simulate some data just to see whether a 4 and 6 category system
% will produce the same values of d' given the same number of random /
% non-random (i.e., correct) responses. (this is just a very coarse way of 
% manipulating signal-to-noise ratio, to verify that out method is behaving as predicted). 
diff = 2;

ncat1 = 4;
corr1 = 1/ncat1*(diff-1);

ncat2 = 5;
corr2 = 1/ncat2*(diff-1);

nt = 10000;
gt1 = randi(ncat1,nt,1);
gt2 = randi(ncat2,nt,1);

resp1 = randi(ncat1,nt,1);
nidx = round(nt*corr1);
addt = randi(ncat1,nidx,1);
resp1 = [resp1;addt];
gt1 = [gt1;addt];

(length(addt)+nt/ncat1) / (nt+length(addt)) 
mean(resp1 == gt1)

resp2 = randi(ncat2,nt,1);
nidx = round(nt*corr2);
addt = randi(ncat2,nidx,1);
resp2 = [resp2;addt];
gt2 = [gt2;addt];

(length(addt)+nt/ncat2) / (nt+length(addt)) 
mean(resp2 == gt2)

% measure dprime
dp1 = mafc_dprime2(resp1,gt1,3);
dp2 = mafc_dprime2(resp2,gt2,3);


% let's check signal-to-noise ratios of 1->inf
ntries = 100;
nests = 1;
ninc = round(nt/ntries);

dp1 = zeros(ntries,nests);
hit1 = zeros(ntries,nests);
fa1 = zeros(ntries,nests);

dp2 = zeros(ntries,nests);
hit2 = zeros(ntries,nests);
fa2 = zeros(ntries,nests);

for i = 1:ntries
    
    for j = 1:nests
        
        % first make all responses match the gt
        resp1 = gt1;
        resp2 = gt2;

        % then we change n trials (increasing per iteration)
        ntochange = ninc*i;
        idxchange = 1:nt;
        idxchange = idxchange(randperm(nt));
        idxchange = idxchange(1:ntochange);
        % and replace the gt-matched resps with random responses
        replaceval = randi(ncat1,ntochange,1);
        resp1(idxchange) = replaceval;
        [dp1(i,j), hit1(i,j),fa1(i,j)] = mafc_dprime(resp1,gt1);

        
        % and replace the gt-matched resps with random responses
        replaceval = randi(ncat2,ntochange,1);
        resp2(idxchange) = replaceval;
        [dp2(i,j), hit2(i,j),fa2(i,j)] = mafc_dprime(resp2,gt2);

        disp([i,j])
    end
end


figure; hold on;
maxd = max([dp1(:);dp2(:)]);
plt1 = plot([0,maxd],[0,maxd],'Color','k','LineStyle','--','LineWidth',3);
plt2 = plot(flipud(dp1),flipud(dp2),'Color','r','LineWidth',3);
set(gca,'LineWidth',3)
set(gca,'FontSize',17)
xlabel('d prime (4 categories)')
ylabel('d prime (6 Categories)')
xlim([0,maxd])
ylim([0,maxd])
box off
axis square
legend([plt1,plt2],'Predicted','Observed')
legend boxoff

