function  S = PLEdata_by_sbj(ext_p)

%PLEDATA_BY_SBJ -- Organize PL data in a struct array with 1 elt per sbj.
%
% S = PLEdata_by_sbj
% S = PLEdata_by_sbj(0)
%
% PLEDATA_BY_SBJ calls IMPORT_PLE_DATA which reads the raw data file
% named 'raw_data.dat'. This creates 19 variables (column-vectors).
% It then rearragnes these vars into a structure array S organized by sbj.
% Each structure has the following fields:
%    1. sbj
%    2. group
%    3. data -- a matrix with one row per trial and 10 columns:
%       1=block [1:32], 2=noise_degr [-15,+15], 3=Cpeak [.10 .15 .23],
%       4=offset [-111,+111], 5=Gabor_degr [-10,+10], 6=congr [0,1]
%       7=resp [0,1,2], 8=correct [0,1], 9=RT [msec], 10=exitcode [0,1,2]
%    4. miss_idx -- a column-vector of indices in DATA
%    5. context -- a cell array of vectors of block indices
%    6. resp_freq -- by block x Cpeak x Gabor_degr x resp, 32x3x2x2
%    7. P_correct -- by block x Cpeak x Gabor_degr, 32x3x2
%    8. P_congr -- P(correct) for "congruent" vs. "incongruent" stimuli;
%                  by block x Cpeak x congr, 32x3x2
%    9. tmean_RT -- by block x Cpeak x Gabor_degr, 32x3x2, 10% trimmed
%   10. congr_RT -- by block x Cpeak x congr, 32x3x2, 10% trimmed
%    ---- The following fields are not computed if 0 is passed as arg.----
%    ---- Use this feature when the STATS Toolbox is not available.   ----
%   11. dprime -- by block x Cpeak  (pHit=P(R=2|S=2),pFA=P(R=2|S=1)), 32x3
%   12. log_beta -- criterion log-values by block x Cpeak, 32x3
%   13. zP_correct := 2*norminv(P_correct) , 32x3x2
%   14. zP_congr   := 2*norminv(P_congr)   , 32x3x2
%
%  See also PLEXP1, IMPORT_PLE_DATA, PLE_SEQUENTIAL.

% Original coding by Alex Petrov, Ohio State University
% $Revision: 1.1.1 $  $Date: 2008-12-12 $
%
% Part of the Hebbian Perceptual Learning Model version 1.1 for MATLAB.
% http://alexpetrov.com/proj/plearn/
% Copyright (c) Alexander Petrov 1999-2008, http://alexpetrov.com
% Please read the LICENSE and NO WARRANTY statement in PercLearn_license.m

%-- Do extended analyses by default.
if (nargin < 1)  ext_p = 1 ; end

%-- Read the raw data file
import_PLE_data ;    % This creates 19 column-vectors in the workspace.

%-- Determine the unique subject numbers
u_sbj = unique(sbj) ;
N_sbj = length(u_sbj) ;

%-- Pack and analyze the data for each individual subject
for k = 1:N_sbj
   curr_sbj = u_sbj(k) ;
   idx = find(sbj==curr_sbj) ;
   fprintf(1,'Processing %d trials for participant %2d...',...
           length(idx),curr_sbj) ;
   S(k) = do_one(ext_p,curr_sbj,group(idx),session(idx),block(idx),...
                 noise_degr(idx),Gabor_Cpeak(idx),offset(idx),...
                 Gabor_degr(idx),resp(idx),correctp(idx),...
                 RT(idx),exitcode(idx)) ;
   fprintf(1,' Done.\n') ;
end

%--- Return S
%%%% End of the main function


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  S = do_one(ext_p,sbj,group,session,block,...
                     noise_degr,Gabor_Cpeak,offset,Gabor_degr,...
                     resp,correctp,RT,exitcode) ;

%-- Check that the group is consistent throughout
gr = min(group) ;
if (gr~=max(group))
   warning(sprintf('\n  min(group)=%d while max(group)=%d.',...
           gr,max(group))) ;
end

%-- Identify the valid trials
% Exitcodes (see DO_PLE_TRIAL): 0='OK',1='bad key',2='missed video frame'
ok = find(exitcode==0) ;
miss_idx = find(exitcode~=0) ;

%-- Check that each session has four blocks of 300 trials each
[counts,vals] = xtab2(session(ok),block(ok)) ;
if (length(vals{2})~=4)
   warning('\n Nonstandard number of blocks per session.') ;
elseif (any(counts(:)~=300))
   warning('\n Nonstandard number of trials per block.');
end
bl = block + 4*(session-1) ;   % the main time variable

%-- Check that all trials within a block have the same NOISE_DEGR
[counts,vals] = xtab2(noise_degr(ok),bl(ok)) ;
if (length(vals{1})~=2)
   warning('\n Confusing noise orientations.') ;
elseif  (any(prod(counts)~=0))  % one of the two frequencies should be 0
   warning('\n Block(s) with mixed noise orientations.');
end

%-- Indices of blocks with one context versus the other
context{1} = find(counts(1,:)>0) ;   % negative background-noise orient
context{2} = find(counts(2,:)>0) ;   % positive background-noise orient

%-- Probability correct by block and difficulty level
corr_freq = xtab(bl(ok),Gabor_Cpeak(ok),Gabor_degr(ok),correctp(ok)) ;
counts = squeeze(sum(corr_freq,4)) ;    % number of cases per cell
if (any(counts(:)~=50))  warning('\n Imbalanced design.') ; end
P_correct = corr_freq(:,:,:,2) ./ counts ;

%-- Congruent vs. incongruent stimuli
congruence = ((noise_degr.*Gabor_degr)>0) ;  % e.g. (-15,-10) v. (-15,+10)
congr_corr = xtab(bl(ok),Gabor_Cpeak(ok),congruence(ok),correctp(ok)) ;
counts = squeeze(sum(congr_corr,4)) ;
P_congr = congr_corr(:,:,:,2) ./ counts ;

%-- Sanity check: stimulus 1 must be congruent with context 1, etc.
mismatch = (P_correct(context{1},:,[1,2])~=P_congr(context{1},:,[1,0]+1)) ;
if (any(mismatch(:)))
   warning('\n Problems with "congruence" within Context 1.') ;
end
mismatch = (P_correct(context{2},:,[1,2])~=P_congr(context{2},:,[0,1]+1)) ;
if (any(mismatch(:)))
   warning('\n Problems with "congruence" within Context 2.') ;
end

%-- Cross-tabulation of RESPonse
resp_freq = xtab(bl(ok),Gabor_Cpeak(ok),Gabor_degr(ok),resp(ok)) ;
counts = squeeze(sum(resp_freq,4)) ;    % number of cases per cell

%-- Trimmed mean response times
RTs_by_cell = split(RT(ok),[Gabor_degr(ok),Gabor_Cpeak(ok),bl(ok)]) ;
RTs_by_cell = cat(2,RTs_by_cell{:}) ;    % 50x192 = 50x(2x3x32)
N_blocks = size(RTs_by_cell,2)/6 ;
sorted_RTs = sort(RTs_by_cell) ;
tmean_RT = mean(sorted_RTs(1:45,:)) ;    % ignore 5 longest RTs per cell
tmean_RT = reshape(tmean_RT,[N_blocks,3,2]) ;

%-- RTs as a function of congruence
congr_RT = zeros(size(tmean_RT)) ;
congr_RT(context{1},:,[1,2]) = tmean_RT(context{1},:,[1,0]+1) ;
congr_RT(context{2},:,[1,2]) = tmean_RT(context{2},:,[0,1]+1) ;

%----- Extended analyses: done only when EXT_P is not zero.
%-- Transform probabilities to Z-values under the normal integral.
%   Zero-frequency samples are set to 0.5 to avoid division by zero
Pmin = 0.5/50 ; Pmax = 1-Pmin ;
if ext_p
   zP_correct = max(Pmin,min(P_correct(:),Pmax)) ;
   zP_correct = reshape(2*norminv(zP_correct),size(P_correct)) ;
   zP_congr = max(Pmin,min(P_congr(:),Pmax)) ;
   zP_congr = reshape(2*norminv(zP_congr),size(P_congr)) ;
end

%-- D-prime conventions:
%   Gabor_degr < 0 --> "noise" ; Gabor_degr > 0 --> "stimulus+noise"
%   resp = 1 = "L" --> "No"    ; resp = 2 = "R" --> "Yes"
if ext_p
   pFA  = resp_freq(:,:,1,2) ./ counts(:,:,1) ;  % false alarms
   pHit = resp_freq(:,:,2,2) ./ counts(:,:,2) ;  % hits
   %- Zero-frequency samples are set to 0.5 to avoid division by zero
   pFA  = max(Pmin,min(pFA,Pmax)) ;
   pHit = max(Pmin,min(pHit,Pmax)) ;
   [dpr,beta] = dprime(pHit,pFA) ;
   %- Sanity check: d'=z(Hit)-z(FA)=z(Hit)+z(CR)=mean(zP_correct)
   mismatch = abs(dpr - squeeze(mean(zP_correct,3))) ;
   if (any(mismatch > 1e-8))
      warning('DPRIME and zP_CORRECT are not consitent with each other.') ;
   end
end

%-- Pack everything into a big structure and return
S.sbj = sbj ;
S.group = gr ;
S.data = [bl,noise_degr,Gabor_Cpeak,offset,Gabor_degr,congruence,...
          resp,correctp,RT,exitcode] ;
S.miss_idx = miss_idx ;
S.context = context ;
S.resp_freq = resp_freq ;
S.P_correct = P_correct ;
S.P_congr = P_congr ;
S.tmean_RT = tmean_RT ;
S.congr_RT = congr_RT ;
if ext_p   % extended analysis, STATS Toolbox available
   S.dprime = dpr ;
   S.log_beta = log(beta) ;
   S.zP_correct = zP_correct ;
   S.zP_congr = zP_congr ;
end

%--- Return S
%%%% End of subfunction DO_ONE

%%%%% End of PLEDATA_BY_SBJ.M
