adaptive_track_engine(answers,steps,n_reversals,rules)

function [change,reversals] = adaptive_track_engine(answers,steps,n_reversals,rules)
%ADAPTIVE_TRACK_ENGINE determines changes to psychoacoustic adaptive tracks
%   [change,reversals] = adaptive_track_engine(answers,steps,n_reversals,rules)
%
%   answers = logical vector, e.g.,[0 1 1 1 0 1] wrong = 0,correct = 1.
%   steps = vector, e.g., [6 3 1] (dB step or any arbitrary unit).
%   n_reversals = integer vector containing the number of reversals before
%                 switching stages e.g., [3 4 4] meaning use the first step
%                 size for 3 reversals, the second for 4 reversals, etc.
%   rules = integer vector, decrement after this many correct, increment
%                 after this many incorrect, [n_down n_up] e.g., [3 1].
%
%   Based on an adaptive track with steps of 6 dB initial and 3 dB final, 
%   3 reversals for the first stage, 5 for the second, and a 3 down 1 up
%   method, given the listener's answer history so far:
%   answers = [0 1 1 1 1 1 1 1 1 1];
%   [change,reversals] = ADAPTIVE_TRACK_ENGINE(answers,[6 2],[3 5],[3 1])
%
%   change = -6
%   reversals = []
%
%   'change' will tell you if and by how much you need to change your
%   signal on the next trial.
%   'reversals' will be empty until total reversals are exceeded, at which 
%   stage it will contain a logical index of the points in the track when 
%   reversals occurred.
%
%   Later on in the adaptive track:
%   answers = [0 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 0 0 1 1 1 0];
%   [change,reversals] = ADAPTIVE_TRACK_ENGINE(answers,[3 1],[3 5],[3 1])
%
%   change = 2
%   reversals = []
%
%   The code thus required to implement an entire adaptive track is simple:
%   [answers, reversals] = deal([]); %clear out the reversals and answers
%   levels = 50; %this is the starting level
%   while isempty(reversals),
%       play_signal(sum(levels)); %play your stimulus at sum(levels) level
%       answers = [answers;get_answer()]; %get your subject response
%       [levels(end+1),reversals] = adaptive_track_engine(answers,[6 2],[3 4],[3 1]);
%   end
%   levels = cumsum(levels); %this is the vector of levels that were used
%
%   Author: W. Owen Brimijoin, MRC Institute of Hearing Research
%   Edited: 23/03/13 - now returns full reversal index. 

%start input checking:
if nargin~=4,
    error('Incorrect number of input arguments')
end
if length(rules)~=2,
    error('''rules'' invalid: Specify n_down and n_up, e.g., [3 1]')
end
if sum(rem([n_reversals(:);rules(:)],1))~=0 || length(steps)~=length(n_reversals),
    error('''steps'' and ''reversals'' must be equal length vectors with integer elements')
end
if isempty(answers)|| ~isvector(answers) || min(ismember(answers,[0 1]))==0,
    error('''answers'' must be a logical vector')
end
%end input checking

answers = [answers(:),abs(answers(:)-1)];%invert 'answers' for 'wrong'
[changes, reversals] = deal(zeros(1,size(answers,1))); %preallocate change log

%set incorrect answers to -sum of preceding correct and vice versa:
answers(answers(:,1)==0,1) = 1-diff([0;find(answers(:,1)==0)]);
answers(answers(:,2)==0,2) = 1-diff([0;find(answers(:,2)==0)]);
tally = cumsum(answers,1);%determine length of runs of same answers
tally(tally==0) = NaN; %replace zeros with NaN.
changes(mod(tally(:,1),rules(1))==0) = -1;%mod by rules(1) to find decrements (==0)
changes(mod(tally(:,2),rules(2))==0) = 1;%mod by rules(2) to find increments (==0)
change = changes(end); %this is the sign of the needed change, if any

%identify reversals:
changes = [changes;1:length(changes)]; %create trial index
changes(:,changes(1,:)==0) = []; %remove points of no change
changes(:,diff(changes(1,:))==0) = []; %remove points of no difference in change 
reversals(changes(2,:)) = 1; %set reversal points to 1, 
reversals(1) = 0; %the first trial cannot be a reversal; set to 0.

%check if current number of reversals exceeds user specified reversals:
stage = find(cumsum(n_reversals)>sum(reversals)-1,1);
if sum(reversals)-1>=sum(n_reversals),change = 0; %if done, exit.
else change = change*(steps(stage)); reversals = []; %otherwise specify change
end

%the end