function Data = simulDDM_valShiftNormApprox(varargin)
% function simulDDM_valShiftNormApprox([params])
% Determines the probability of boundary crossing for a DDM process given
% particular parameters, using several approximations to achieve a
% reasonable balance between accuracy and speed of estimation (the irony of
% which tradeoff is deeply appreciated)
% This approximation includes a discrete transition probability matrix, 
% in both time and diffusion space (see simPrecision and numSteps
% variables). The model also allows to estimate a single "kink" in the
% drift rate, assuming that prior to some point t the evidence takes one
% value, and after it, a different value.
%
% Author: Cendri Hutcherson
% Date: 04/05/2015
% Last updated (mostly comments and cleaning): 05/01/2019
%
% Arguments:
% params - a structure with fields that specify core variables of the DDM 
% and levels of precision desired for the approximation. If params is not 
% given, or does not contain a fieldname for a given variable, the default
% value will be used
%
% initDrift: scalar specifying the initial rate of value accumulation
%            default = .08
% finalDrift: scalar specifying the final rate of value accumulation
%            default = .08 (i.e. no change)
% timeShift: scalar identifying (in seconds) when the drift begins to 
%            shift from initDrift to finalDrift
% rateShift: scalar from 0 to 1000 specifying how quickly the evidence
%            shifts from intialDrift to finalDrift (allows for gradual
%            transition from one to the other instead of an abrupt shift)
% barrier: decision threshold
%          default = .15;
% startPt: bias in starting point
%          default = 0;
% startPtVar: variability in starting point
%          default = 0;
% nonDec: non-decision time (sum of sensory and motor delays, in secs)
%         default = .3;
% maxRT: maximum response time limit (in secs)
%        default = 4 seconds
% s: noise parameter, by convention set to .1 (all other parameter values
%    will scale to this noise term)
% simPrecision: discrete time step (in secs) to estimate state transitions
%               default = .001
% summaryPrecision: time step for reporting probabilities (in secs)
%                   helps if you want to compress variables for storage;                 
%                   default = .001 (no compression)
% numSteps: bin size for estimate diffusion movement (this is the heart of
%           approximation that allows for faster computation)
%
% nSim: variable that determines how many simulated individual choices and
% RTs are generated (default = 0, no simulations), useful for generating
% simulated data to compare against the observed.

% Output:
% Data structure with following fields:
% pCross1: prob of a no response in each time-bin (uses summaryPrecision)
% pCross2: prob of a yes response in each time-bin (uses summaryPrecision)
% pNonResp: probability of no response within given time limit
% simData: if requested via the use of "nSim", a data structure consisting
% of nSim realizations of choices and RTs from the generated probability
% distributions
%
% Author: Cendri Hutcherson
%
% With thanks to Gabriela Tavares, who wrote a version of the discrete bin 
% estimation in Python, from which this script is adapted
%
% If you use this script, please cite this paper:
%    Harris, A., Clithero, J., & Hutcherson, C. (2018). Accounting for 
%    taste: A multi-attribute neurocomputational model explains the neural 
%    dynamics of choices for self and others. Journal of Neuroscience, 38, 
%    7952-7968.



%% ----------------- INITIALIZE PARAMETERS OF MODEL ---------------------%%

% set defaults for evidence inputs to DDM (e.g. integrated value of the option)
initDrift = .08; % value of the option initially
finalDrift = .08; % value of the option after a transition to final value
timeShift = 0; % time at which value begins to transition from initial to final
rateShift = 1000; % speed at which value transitions from initial to final. Given in units of drift/sec. 1000 is essentially immediately

% set defaults for other parameters of DDM
barrier = .15;  % bounds on DDM process (defines -barrier and +barrier)
startPt = 0; % starting point
startPtVar = 0; % starting point variability
nonDec = .3; % non-decision time for each attribute
maxRT = 4;

s = .1; % within-trial drift parameter, serving as scaling parameter that by convention is set to .1

% set defaults for estimation
simPrecision = .001; % simulation time precision in seconds
summaryPrecision = .001; % precision for summary - helps to condense for data storage

numSteps = 50; % defines granularity of states

nSim = 0; 

% overwrite defaults for any specified parameters (contained in params
% structure)
if ~isempty(varargin)
    paramNames = fieldnames(varargin{1});
    params = varargin{1};
    
    for i = 1:length(paramNames)
        eval([paramNames{i} '= params.' paramNames{i} ';']);
    end
end

%% -------------- SET UP THE DDM  -------------- %%
nTimeSteps = floor(maxRT/simPrecision);

binLowEdge = [-Inf linspace(-barrier, barrier, numSteps + (mod(numSteps,2) == 1))]; % define bins for each drift state
binHiEdge = [linspace(-barrier, barrier, numSteps + (mod(numSteps,2) == 1)) Inf];
CoMstates = (binLowEdge + binHiEdge)/2; % define the center of each state-bin

% create distance matrices for the low and high end of each bin from each
% bin center
allDiffsLo = repmat(binLowEdge',1,length(CoMstates)) - repmat(CoMstates,length(CoMstates),1);
allDiffsHi = repmat(binHiEdge',1,length(CoMstates)) - repmat(CoMstates,length(CoMstates),1);

pStates = zeros(length(CoMstates),1); % # of states + 2 absorbing states (boundary crosses)
pStates(binLowEdge <= (startPt + startPtVar) & binHiEdge > (startPt - startPtVar)) = ...
    1/sum(binLowEdge <= startPt + startPtVar & binHiEdge > startPt - startPtVar); % start centered at startPt with uniform distribution startPtVar
cumCrossUp = zeros(nTimeSteps,1);
cumCrossDown = zeros(nTimeSteps,1);

% estimate initial transition probabilities
pDiffsLo = normcdf(allDiffsLo,initDrift*simPrecision,s*sqrt(simPrecision));
pDiffsHi = normcdf(allDiffsHi,initDrift*simPrecision,s*sqrt(simPrecision));
pTransit = pDiffsHi - pDiffsLo; % creates transition matrix from state to state
pTransit(isnan(pTransit)) = 1;

absorbingHi = binHiEdge > barrier; % state beyond upper barrier
absorbingLo = binLowEdge < -barrier; % state beyond lower barrier
absorbMat = diag(absorbingLo + absorbingHi);
pTransit(:,absorbingHi | absorbingLo) = absorbMat(:,absorbingHi | absorbingLo);

% estimate final transition probabilities
pDiffsFLo = normcdf(allDiffsLo,finalDrift*simPrecision,s*sqrt(simPrecision));
pDiffsFHi = normcdf(allDiffsHi,finalDrift*simPrecision,s*sqrt(simPrecision));
pTransitF = pDiffsFHi - pDiffsFLo; % creates transition matrix from state to state
pTransitF(isnan(pTransitF)) = 1;

pTransitF(:,absorbingHi | absorbingLo) = absorbMat(:,absorbingHi | absorbingLo);

% calculate number of timesteps over which transition from initDrift to
% finalDrift will take place
timeStepsNeededToChange = max((abs(initDrift - finalDrift)/(rateShift*simPrecision)),1);
mixingRate = 1/timeStepsNeededToChange;

%% -------------- RUN THE DDM  -------------- %%
for t = floor(nonDec/simPrecision):nTimeSteps
    % estimate cumulative fraction of DDMs in absorbing states
    cumCrossUp(t) = absorbingHi*pStates;
    cumCrossDown(t) = absorbingLo*pStates;

    if t < timeShift/simPrecision + floor(nonDec/simPrecision) % if in initial evidence phase
        pTransitAtT = pTransit;
    else
        timePastShift = t - timeShift/simPrecision - floor(nonDec/simPrecision);
        if timePastShift < timeStepsNeededToChange
            % average two transition matrices proportional to time past
            % value shift
            pTransitAtT = (1 - timePastShift*mixingRate)*pTransit + (timePastShift*mixingRate)*pTransitF;
        else
            % otherwise, use final value
            pTransitAtT = pTransitF;
        end
    end
    
    % compute movements of probability mass from current location
    pStates = pTransitAtT*pStates;
    if(sum(pStates(~absorbingLo & ~absorbingHi)) < .0001) % quit if sufficiently low probability not to have terminated
        break;
    end
end

% compute cumulative probabilities of boundary crossing at each time pt
cumCrossUp(t+1:end) = cumCrossUp(t);
cumCrossDown(t+1:end) = cumCrossDown(t);

% Use this to compute hte instantaneous probability of crossing at any
% given time pt
pCross1 = [0;diff(cumCrossDown)];
pCross2 = [0;diff(cumCrossUp)];

% if summarizing the distribution at a lower precision than estimated (e.g.
% for purposes of saving memory storage)
condenseFactor = summaryPrecision/simPrecision;
if condenseFactor > 1
    pCross1 = reshape(pCross1,condenseFactor,length(pCross1)/condenseFactor);
    pCross1 = sum(pCross1,1);
    pCross2 = reshape(pCross2,condenseFactor,length(pCross2)/condenseFactor);
    pCross2 = sum(pCross2,1);
end

% create the data structure to be returned by the function
Data.pCross1 = pCross1;
Data.pCross2 = pCross2;
Data.pNonResp = 1 - sum(Data.pCross1) - sum(Data.pCross2);

% simulate nSim individual choices and RTs, if requested
if nSim > 0
    rt = NaN*zeros(nSim,1);
    resp = NaN*zeros(nSim,1);

    temp = rand(nSim,1);
    resp = zeros(nSim,1);
    rt = NaN*zeros(nSim,1);
    resp(temp <= sum(Data.pCross1)) = -1;
    resp(temp > sum(Data.pCross1)) = 1;
    resp(temp > (1 - Data.pNonResp)) = NaN;

    pRTNo = cumsum(Data.pCross1);
    pRTYes = cumsum(Data.pCross2) + sum(Data.pCross1);
    for j = 1:nSim
        if ~isnan(resp(j))
            if resp(j) == -1
                rt(j) = find(pRTNo > temp(j),1)/1000;
            else
                rt(j) = find(pRTYes > temp(j),1)/1000;
            end
        end

    end
    Data.simData.resp = resp;
    Data.simData.rt = rt;
end