function  stats = PLE2_stats(D)

%PLE2_STATS -- Aggregate statistics for Perceptual Learning Experiment 2.
%
% stats = PLE2_stats(D)
%
% D is a 'raw data' structure array imported by IMPORT_PLE2_DATA.
% The result is a structure array STATS with the following fields:
%  - sbj -- subject identificaiton number, [1x1]
%  - group -- 1="first block has context L", 2="1st context R", [1x1]
%  - zP_congr -- z-probability correct for the congruent stimuli, [36x3]
%  - zP_incon -- z-probability correct for the incongruent stim., [36x3]
%  - dprime := zP_congr+zP_incon, [36x3]
%  - accZP -- zP collapsed across blocks, [1x6], 1:3=congr, 4:6=incongr
%  - RT_congr -- 10%-trimmed mean response time, congr. stim, [36x3]
%  - RT_incon -- 10%-trimmed mean response time, incongruent, [36x3]
%  - RT = (RT_conr+RT_incon)/2, [36x3]
%  - accRT -- RT collapsed across blocks, [1x6], 1:3=congr, 4:6=incongr
%  - seq_data -- 16x10 matrix produced by PLM_SEQUENTIAL
%  - prSprR -- P("R") conditional on <S(t-1),R(t-1)>=[<LL> <LR> <RL> <RR>]
%  - seq_eff -- [maineff_prevS , maineff_prevR]
%  % contextc -- zP averaged over 5 context-switch transitions, (10+10)
%  % contexti -- zP, incongr. stim, context-switch transitions, (10+10)
%  % otherc -- zP,   congruent stim, remaining 7+19 transitions,(10+10)
%  % otheri -- zP, incongruent stim, remaining 7+19 transitions,(10+10)
%
%  Example (Plot partial d' curves as data become available):
%    D=import_PLE2_data ; st=PLE2_stats(D) ; plot(st(1).dprime,'.-') ;
%
%  Example (Aggregate into a single struct. All STATS must be same size.)
%    st = PLE2_stats(import_PLE2_data('raw_data.dat')) ; N = length(st);
%    S=st(1);for k=2:N S=merge_PLM_stats(S,st(k));end; sst=PLM_summary(S);
%
%  See also IMPORT_PLE2_DATA,PLM_STATS,PLM_SUMMARY,MERGE_PLM_STATS.

%  Original coding by Alexander Petrov, UC Irvine.
%  $Revision: 1.0 $  $Date: 2003/08/27 10:25 $

%--  Analyze each subject separately and combine the results
for k = 1:length(D)
   ph1 = do_phase1(D(k)) ;
   ph2 = do_phase2(D(k)) ;
   ph3 = do_phase3(D(k)) ;
   st.sbj = ph1.sbj ; st.group = ph1.group ;
   st.zP_congr = [ph1.zP_congr ; ph2.zP_congr ; ph3.zP_congr ] ;
   st.zP_incon = [ph1.zP_incon ; ph2.zP_incon ; ph3.zP_incon ] ;
   st.dprime = [ph1.dprime ; ph2.dprime ; ph3.dprime ] ;
   st.accZP1 = ph1.accZP1 ;
   st.accZP2 = ph2.accZP2 ;
   st.accZP3 = ph3.accZP3 ;
   st.RT_congr = [ph1.RT_congr ; ph2.RT_congr ; ph3.RT_congr ] ;
   st.RT_incon = [ph1.RT_incon ; ph2.RT_incon ; ph3.RT_incon ] ;
   st.RT = [ph1.RT ; ph2.RT ; ph3.RT ] ;
   st.accRT1 = ph1.accRT1 ;
   st.accRT2 = ph2.accRT2 ;
   st.accRT3 = ph3.accRT3 ;
   st.seq_data = ph1.seq_data ;
   st.prSprR = ph1.prSprR ;
   st.deltaP = ph1.deltaP ;
   st.seq_eff = ph1.seq_eff ;
   stats(k) = st ;
end

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  st = do_phase1(D) ;

st.sbj = D.sbj ;
st.group = D.group ;

%-- Phase 1 comprises the first 9 sessions. See PLExp2.
Ph1_sessions = 9 ;
session = D.data(:,1) ;
block = D.data(:,2) ;
exitcode = D.data(:,10) ;
ph1= find((session<=Ph1_sessions) & (block>0) ) ;  % bad trials included
ok = find((session<=Ph1_sessions) & (block>0) & (exitcode==0)) ;

%-- Unpack the relevant variables from the raw data
Gabor_degr = D.data(ok,3) ;     % -10 or +10 degrees
noise_degr = D.data(ok,5) ;     % -15 or +15 degrees
difficulty = D.stim(ok,2) ;     % 1=0.23, 2=0.15, 3=0.10

%resp = D.data(ok,7) ;           % 1="Left", 2="Right"
correct = D.data(ok,8) ;        % 0="incorr",1="correct"
RT = D.data(ok,9) ;             % msec from the stimulus onset

%-- Derivative variables
congruence = 2-((noise_degr.*Gabor_degr)>0) ; % 1=congruent, 2=incongr
bl = (session(ok)-1).*4 + block(ok) ;

%-- Probability correct by Block, Difficulty, and Congruence
[X,vals] = xtab(bl,difficulty,congruence,correct) ;  % [36x3x2x2]
N_blocks = size(X,1) ;       % 36 = 9 days * 4 blocks/day in a full set
counts = squeeze(sum(X,4)) ;                         % [36x3x2]
if (any(counts(:)~=50))  warning('\n Imbalanced design.') ; end
P = X(:,:,:,2) ./ counts ;       % Probab. correct, [36x3x2]
P = max(0.01,min(P,0.99)) ;      % Avoid extreme frequencies

%-- Z-transform all probabilities and split by Congruence
zP_congr = norminv(squeeze(P(:,:,1))) ;      % [36x3]
zP_incon = norminv(squeeze(P(:,:,2))) ;      % [36x3]

%-- Export the z-transformed probabilities
st.zP_congr = zP_congr ;
st.zP_incon = zP_incon ;
st.dprime = zP_congr + zP_incon ;            % d' = zP(Hit) + zP(CR)
st.accZP1 = [mean(zP_congr) , mean(zP_incon)] ;    % [ccc,iii]

%-- Group RTs in "cells" defined by Block, Difficulty, and Congruence
RTs_by_cell = split(RT,[congruence,difficulty,bl]) ;
RTs_by_cell = cat(2,RTs_by_cell{:}) ;    % [50x216] = 50x(2*3*36)

%-- Ignore the 5 longest response times in each cell
sorted_RTs = sort(RTs_by_cell) ;
tmean_RT = mean(sorted_RTs(1:45,:)) ;                % [1x216]
tmean_RT = reshape(tmean_RT,[N_blocks,3,2]) ;        % [36x3x2]

%-- Export as separate fields
st.RT_congr = tmean_RT(:,:,1) ;
st.RT_incon = tmean_RT(:,:,2) ;
st.RT = mean(tmean_RT,3) ;    % = (RT_congr+RT_incon)/2
st.accRT1 = [mean(st.RT_congr) , mean(st.RT_incon)] ; % [ccc,iii]

%-- Sequential analysis
stim = D.stim(ph1,1:3) ;
resp = D.data(ph1,7) ;
bl1 = (session(ph1)-1).*4 + block(ph1) ;

seq = PLM_sequential(stim,resp,bl1) ;  % <-- main work
condP = reshape(seq(:,9),4,4) ;        % 4 non-lagged x 4 lagged vars
prSprR = mean(condP') ;   % <S(t-1),R(t-1)> = [<LL> <LR> <RL> <RR>]
maineff_prevS = (prSprR(3)-prSprR(1)+prSprR(4)-prSprR(2))/2 ;
maineff_prevR = (prSprR(2)-prSprR(1)+prSprR(4)-prSprR(3))/2 ;

st.seq_data = seq ;
st.prSprR = prSprR ;
st.deltaP = mean(reshape(seq(:,10),4,4)') ;
st.seq_eff = [maineff_prevS maineff_prevR] ;

%--- Return STATS
%%%% End of subfunction DO_PHASE1     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  st = do_phase2(D) ;

st.sbj = D.sbj ;
st.group = D.group ;

%-- Phase 2 consists of sessions 10 and 11. See PLExp2a.
Ph1_sessions = 9 ;
Ph2_sessions = 11 ;
session = D.data(:,1) ;
block = D.data(:,2) ;
exitcode = D.data(:,10) ;
ok = find((session>Ph1_sessions)&(session<=Ph2_sessions) & ...
          (block>0) & (exitcode==0)) ;

%-- Abort if this subject hasn't gotten so far yet...
if (isempty(ok))
   st.zP_congr = zeros(0,3) ; st.zP_incon = zeros(0,3) ;
   st.dprime = zeros(0,3) ; st.accZP2 = zeros(0,6) ;
   st.RT_congr = zeros(0,3) ; st.RT_incon = zeros(0,3) ;
   st.RT = zeros(0,3) ; st.accRT2 = zeros(0,6) ;
   return
end

%-- Unpack the relevant variables from the raw data
Gabor_degr = D.data(ok,3) ;     % 080 or 100 degrees
noise_degr = D.data(ok,5) ;     % 075 or 105 degrees
difficulty = D.stim(ok,2) ;     % 1=0.23, 2=0.15, 3=0.10

%resp = D.data(ok,7) ;           % 1="Left", 2="Right"
correct = D.data(ok,8) ;        % 0="incorr",1="correct"
RT = D.data(ok,9) ;             % msec from the stimulus onset

%-- Derivative variables
congruence = 2-(((noise_degr-90).*(Gabor_degr-90))>0) ; % 1=congr,2=inc
bl = (session(ok)-1).*4 + block(ok) ;

%-- Probability correct by Block, Difficulty, and Congruence
[X,vals] = xtab(bl,difficulty,congruence,correct) ;  % [36x3x2x2]
N_blocks = size(X,1) ;       % 36 = 9 days * 4 blocks/day in a full set
counts = squeeze(sum(X,4)) ;                         % [36x3x2]
if (any(counts(:)~=50))  warning('\n Imbalanced design.') ; end
P = X(:,:,:,2) ./ counts ;       % Probab. correct, [36x3x2]
P = max(0.01,min(P,0.99)) ;      % Avoid extreme frequencies

%-- Z-transform all probabilities and split by Congruence
zP_congr = norminv(squeeze(P(:,:,1))) ;      % [36x3]
zP_incon = norminv(squeeze(P(:,:,2))) ;      % [36x3]

%-- Export the z-transformed probabilities
st.zP_congr = zP_congr ;
st.zP_incon = zP_incon ;
st.dprime = zP_congr + zP_incon ;            % d' = zP(Hit) + zP(CR)
st.accZP2 = [mean(zP_congr) , mean(zP_incon)] ;    % [ccc,iii]

%-- Group RTs in "cells" defined by Block, Difficulty, and Congruence
RTs_by_cell = split(RT,[congruence,difficulty,bl]) ;
RTs_by_cell = cat(2,RTs_by_cell{:}) ;    % [50x216] = 50x(2*3*36)

%-- Ignore the 5 longest response times in each cell
sorted_RTs = sort(RTs_by_cell) ;
tmean_RT = mean(sorted_RTs(1:45,:)) ;                % [1x216]
tmean_RT = reshape(tmean_RT,[N_blocks,3,2]) ;        % [36x3x2]

%-- Export as separate fields
st.RT_congr = tmean_RT(:,:,1) ;
st.RT_incon = tmean_RT(:,:,2) ;
st.RT = mean(tmean_RT,3) ;    % = (RT_congr+RT_incon)/2
st.accRT2 = [mean(st.RT_congr) , mean(st.RT_incon)] ; % [ccc,iii]

%--- Return STATS
%%%% End of subfunction DO_PHASE2


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  st = do_phase3(D) ;

st.sbj = D.sbj ;
st.group = D.group ;

%-- Phase 3 consists of session 12 only. There is no demo. See PLExp2b.
Ph3_session = 12 ;
session = D.data(:,1) ;
exitcode = D.data(:,10) ;
ok = find((session==Ph3_session) & (exitcode==0)) ;

%-- Abort if this subject hasn't gotten so far yet...
if (isempty(ok))
   st.zP_congr = zeros(0,3) ; st.zP_incon = zeros(0,3) ;
   st.dprime = zeros(0,3) ; st.accZP3 = zeros(0,6) ;
   st.RT_congr = zeros(0,3) ; st.RT_incon = zeros(0,3) ;
   st.RT = zeros(0,3) ; st.accRT3 = zeros(0,6) ;
   return
end

%-- Unpack the relevant variables from the raw data
Gabor_degr = D.data(ok,3) ;     % -10 or +10 degrees
noise_degr = D.data(ok,5) ;     % -15 or +15 degrees
difficulty = D.stim(ok,2) ;     % 1=0.23, 2=0.15, 3=0.10

%resp = D.data(ok,7) ;           % 1="Left", 2="Right"
correct = D.data(ok,8) ;        % 0="incorr",1="correct"
RT = D.data(ok,9) ;             % msec from the stimulus onset

%-- Derivative variables
congruence = 2-((noise_degr.*Gabor_degr)>0) ; % 1=congruent, 2=incongr
bl = (session(ok)-1).*4 + D.data(ok,2) ;

%-- Probability correct by Block, Difficulty, and Congruence
[X,vals] = xtab(bl,difficulty,congruence,correct) ;  % [36x3x2x2]
N_blocks = size(X,1) ;       % 36 = 9 days * 4 blocks/day in a full set
counts = squeeze(sum(X,4)) ;                         % [36x3x2]
if (any(counts(:)~=52))  warning('\n Imbalanced design.') ; end
P = X(:,:,:,2) ./ counts ;       % Probab. correct, [36x3x2]
P = max(0.01,min(P,0.99)) ;      % Avoid extreme frequencies

%-- Z-transform all probabilities and split by Congruence
zP_congr = norminv(squeeze(P(:,:,1))) ;      % [36x3]
zP_incon = norminv(squeeze(P(:,:,2))) ;      % [36x3]

%-- Export the z-transformed probabilities
st.zP_congr = zP_congr ;
st.zP_incon = zP_incon ;
st.dprime = zP_congr + zP_incon ;            % d' = zP(Hit) + zP(CR)
st.accZP3 = [mean(zP_congr) , mean(zP_incon)] ;    % [ccc,iii]

%-- Group RTs in "cells" defined by Block, Difficulty, and Congruence
RTs_by_cell = split(RT,[congruence,difficulty,bl]) ;
RTs_by_cell = cat(2,RTs_by_cell{:}) ;    % [50x216] = 50x(2*3*36)

%-- Ignore the 5 longest response times in each cell
sorted_RTs = sort(RTs_by_cell) ;
tmean_RT = mean(sorted_RTs(1:45,:)) ;                % [1x216]
tmean_RT = reshape(tmean_RT,[N_blocks,3,2]) ;        % [36x3x2]

%-- Export as separate fields
st.RT_congr = tmean_RT(:,:,1) ;
st.RT_incon = tmean_RT(:,:,2) ;
st.RT = mean(tmean_RT,3) ;    % = (RT_congr+RT_incon)/2
st.accRT3 = [mean(st.RT_congr) , mean(st.RT_incon)] ; % [ccc,iii]

%--- Return STATS
%%%% End of subfunction DO_PHASE3

%%%%% End of PLE2_STATS.M
