function  seq = PLM_sequential(stim,resp,block)

%PLM_SEQUENTIAL -- Frequencies as a fun. of previous stimulus and response.
%
% seq = PLM_sequential(stim,resp)
% seq = PLM_sequential(stim,resp,block)
%
% STIM must be a matrix with one row per trial and at least three 
% columns (usually generated by BLOCK_FROM_CACHE or IMPORT_PLE2_DATA).
% Only columns 1 and 3 are relevant---they specify the orientation of 
% the Gabor target (1=Left,2=Right) and the background noise, respvly.
% RESP is a column vector (1="Left", 2="Right", 0=bad_trial).
% BLOCK is a column vector of unique block numbers. If not supplied,
% it defaults to reshape(repmat([1:N_trials/300],300,1),N_trials,1).
%
% The first trial in each block and any trials with RESP=0 are discarded.
% The remaining trials are grouped by the three IVs below and the
% frequencies of giving four kinds of responses (DVs) are counted.
% Note that all contrast levels are averaged out. Thus, the result 
% SEQ is a matrix with 16=2x2x2x2 rows (conditions) and 4+4+2 columns 
% (4 IVs + 4 DVs + 2 probabilities) as follows:
%  1. curr_Gabor -- orientation of the current target (1=Left,2=Right)
%  2. backgr -- background noise orientation (1=Left,2=Right)
%  3. prev_Gabor -- orientation of the previous Gabor target (1=L,2=R)
%  4. prev_resp -- previous response (1="Left",2="Right")
% ----
%  5. N_total -- total number of cases in the condition
%  6. N_Left -- number of cases with current response = "Left"
%  7. N_Right -- number of cases w curr resp = "Right" (==N_total-N_Left)
%  8. N_correct -- number of correct responses
% ----
%  9. P_Right = N_Right ./ N_total
% 10. delta_P -- P_Right with curr_Gabor and backgr partialled out
%
%  See also XTAB, PLM_STATS, PLE2_STATS, PLE_SEQUENLR.

%  Original coding by Alexander Petrov, UC Irvine.
%  $Revision: 1.0 $  $Date: 2003/08/24 19:29 $

N_trials = size(resp,1) ;

%-- Default BLOCK if necessary
DEFAULT_BLOCK_SIZE = 300 ;
if (nargin<3)
   if (rem(N_trials,DEFAULT_BLOCK_SIZE)~=0)
      error('Non-standard sequence. Please supply BLOCK explicitly.') ;
   end
   block = reshape(repmat([1:N_trials/300],300,1),N_trials,1) ;
end

%-- Arrange all relevant data in a big matrix
data = zeros(N_trials,6) ;
data(:,1) = [NaN;block(1:N_trials-1)] ;   % previous block
data(:,2) = stim(:,3) ;                   % current background B(t)=B(t-1)
data(:,3) = [NaN;stim(1:N_trials-1,1)] ;  % previous stimulus  S(t-1)
data(:,4) = stim(:,1) ;                   % current stimulus   S(t)
data(:,5) = [0;resp(1:N_trials-1,1)] ;    % previous response  R(t-1)
data(:,6) = resp ;                        % current response   R(t)

%-- Discard invalid or incomplete trials
idx = find(block==data(:,1)) ; data = data(idx,:) ;  % curr_bl==prev_bl
idx = find(data(:,5)~=0) ; data = data(idx,:) ;      % R(t-1) is valid
idx = find(data(:,6)~=0) ; data = data(idx,:) ;      % R(t) is valid too

%-- Cross-tabulate CURR_RESP on the four independent variables
[N_LR,vals] = xtab2(data(:,[4,2,3,5]),data(:,6)) ;
if (length(vals{2})~=2)   % vals{2} must be [1 2], not [-9 1 2] or whatnot
  error('Invalid response value(s): [',num2str(vals{2}),']') ;
end
IVs = vals{1} ;             % [curr_stim,bkgr,prev_stim,prev_resp]
N_total = sum(N_LR,2) ;     % total number of cases per condition

%-- Correct responses
N_correct = zeros(16,1) ;
idx = find(IVs(:,1)==1) ; N_correct(idx) = N_LR(idx,1) ;  % R(t)=S(t)=1
idx = find(IVs(:,1)==2) ; N_correct(idx) = N_LR(idx,2) ;  % R(t)=S(t)=2

counts = [N_total,N_LR,N_correct] ;

%-- Probability of responding "Right" on the current trial
P = N_LR(:,2) ./ N_total ;      % [16x1]

%-- Partial out the current Gabor and background orientation
mP = mean(reshape(P,4,4)) ;           % [4x1]
mP = reshape(repmat(mP,4,1),16,1) ;
P1 = P - mP ;

%-- Pack and return
seq = [IVs, counts, P P1] ;

%--- Return SEQ
%%%% End of PLM_SEQUENTIAL.M
