function thresh = getThreshold_1
%function [thresh results pBest] = getThreshold(k)
%
% Winning entry - uses 3-down 1-up staircase and a psychometric function
% derived from signal detection theory.


nTrials = 60;

I = .5;             %starting intensity
stepSize = .2;
firstChange = 5;    %number of steps before changing stepSize
newStepSize = .1;

correctInARow = 0;
results.intensity = zeros(1,nTrials);
results.response = zeros(1,nTrials);
for i = 1:nTrials
    if i == firstChange
        stepSize = newStepSize;
    end
    
    results.intensity(i) = I;
    results.response(i) = idealObserver(I);
%     results.response(i) = myIdealObserver(I,k);

    if results.response(i) == 1
        correctInARow = correctInARow+1;
        if correctInARow == 3
            I = I - stepSize;
            correctInARow = 0;
        end
    else
        correctInARow = 0;
        I = I + stepSize;
    end
    %limit floor and ceiling
    I = max([.1, min([I .9])]);
end

%assume HOV and that noise mean = 0
pInit.shutup = 1;
pInit.noiseMean = 0;
pInit.theta = 5;        %where theta = k / sigma and the sole free param
% pInit.lambda = 1;

SDTfun_h = @SDTfun;

pBest = fit('fitPsychometricFunction',pInit,{'theta'},results,SDTfun_h); %1 free param

%calculate d' for a thresh of about ~80% correct
dprimeThresh = inverseNormalCDF((1/2)^(1/3));
%calculate the threshold analytically based on best fit params
thresh = 2*normcdf(sqrt(2)*dprimeThresh / pBest.theta) - 1;

end

%% SDT function

function pC = SDTfun(p,x)
%calculate the response using SDT

dprimes = p.theta*inverseNormalCDF((x+1)/2);
pC = normcdf(dprimes / sqrt(2));

%Hack! bound free params by exploding error when they are outside the bounds
if p.theta > 9.221 || p.theta < .70444
    pC = zeros*pC;
end

%these numbers exploit the prior knowledge that Geoff will not test with an
%observer outside of .1 and .9

% thetaFloor = sqrt(2)*inverseNormalCDF((1/2)^(1/3)) / (inverseNormalCDF((.1+1)/2));   % 9.2209
% thetaCeiling = sqrt(2)*inverseNormalCDF((1/2)^(1/3)) / (inverseNormalCDF((.9+1)/2)); % .70444 BUT at < 1.5 the slope drops bellow 1 which is a thresh of about .56016%
%
% %unlock the ceiling
% pC = pC*(p.lambda) + (1-p.lambda)/2;
% if p.lambda > 1 || p.lambda < .95
%     pC = zeros*pC;
% end

end

%% my ideal observer

function resp = myIdealObserver(I,k)
%use my own idealObserver to test params

if rand(1) < .02;
    %lapses, suck
    resp = round(rand(1));
else
    signalMean = k*inverseNormalCDF((I+1)/2);
    internalRespSignal= randn(1,length(signalMean))+signalMean;
    internalRespNoise= randn(1,length(signalMean));
    resp = internalRespSignal>internalRespNoise;
end

end