function [traj, infStates] = tapas_ph_binary(r, p, varargin)
% Calculates the trajectories of v under the Pearce-Hall learning model
%
% This function can be called in two ways:
% 
% (1) tapas_ph_binary(r, p)
%   
%     where r is the structure generated by tapas_fitModel and p is the parameter vector in native space;
%
% (2) tapas_ph_binary(r, ptrans, 'trans')
% 
%     where r is the structure generated by tapas_fitModel, ptrans is the parameter vector in
%     transformed space, and 'trans' is a flag indicating this.
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2015 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Transform paramaters back to their native space if needed
if ~isempty(varargin) && strcmp(varargin{1},'trans');
    p = tapas_ph_binary_transp(r, p);
end

% Unpack parameters
v_0  = p(1);
al_0 = p(2);
S    = p(3);

% Add dummy "zeroth" trial
u = [0; r.u(:,1)];
n = length(u);

% Initialize updated quantities: value v and associability al
v  = NaN(n,1);
al = NaN(n,1);
da = NaN(n,1);

% Prior
v(1)  = v_0;
al(1) = al_0;
da(1) = al_0;

% Pass through value update loop
for k = 2:1:n
    if not(ismember(k-1, r.ign))
        
        %%%%%%%%%%%%%%%%%%%%%%
        % Effect of input u(k)
        %%%%%%%%%%%%%%%%%%%%%%
        
        % Prediction error
        da(k) = u(k)-v(k-1);
        
        % Learning rate
        al(k) = abs(da(k-1));
        
        % Value
        v(k) = v(k-1)+S*al(k)*da(k);
    else
        da(k) = 0;
        al(k) = al(k-1);
        v(k)  = v(k-1);
    end
end

% Predicted value
vhat = v;
vhat(end) = [];

% Remove representation priors
v(1)  = [];
al(1) = [];
da(1) = [];

% Create result data structure
traj = struct;

traj.v     = v;
traj.vhat  = vhat;
traj.al    = al;
traj.da    = da;

% Create matrix (in this case: vector) needed by observation model
infStates = NaN(n-1,3);
infStates(:,1) = traj.vhat;
infStates(:,2) = traj.v;
infStates(:,3) = traj.al;

return;