function  stats = PLM_stats(stim,resp,act,inp,warnp)

%PLM_STATS -- Aggregate statistics, perceptual-learning data.
%
% stats = PLM_stats(stim,resp)
% stats = PLM_stats(stim,resp,act,inp)
%
% STIM is is a matrix describing the trial sequence. It has N_TRIALS rows
% and N_RUNS pages. It must have at least 3 columns (it usually has 5 as
% generated by BLOCK_FROM_CACHE):
%   1. Gabor orientation: 1=-10 or 2=+10 deg, see PLEXP1_PARAMS
%   2. Gabor peak contrast: 1=0.23, 2=0.15 or 3=0.10, see PLEXP1_PARAMS
%   3. noise orientation: 1=-15 or 2=+15 deg, see PLEXP1_PARAMS
% RESP is a matrix of 1s and 2s. N_TRIALS rows and N_RUNS columns.
% ACT and INP are optional and have the same format as RESP. They 
% correspond to the second and third columns of PLMODEL1's OUTPUT.
% (Singleton dimensions are squeezed automatically, see example below.)
% When WARNP=0, no warning messages are prnted. Default=1.
%
% The result is a structure STATS with the following fields:
%  - group -- 1="first block has context L", 2="1st context R", [Nx1]
%  - zP_congr -- z-probability correct for the congruent stimuli,
%          N_BLOCKS=32 by DIFFICULTY=3 values per run, [32x3xN]
%  - zP_incon -- z-probability correct for the incongruent stim,[32x3xN]
%  - dprime := zP_congr+zP_incon, [32x3xN]
%  - accZP -- zP collapsed across blocks, [Nx6], 1:3=congr, 4:6=incongr
% % - seqP -- sequential: 1=P(Rt=R|St-1=L,Rt-1=L),2=P(Rt=R|St-1=R,Rt-1=L),
% %           3=P(Rt=R|St-1=L,Rt-1=2),4=P(Rt=R|St-1=R,Rt-1=R), [Nx4]
%  - mnA_congr -- mean activation (aligned with context), congruent stim.
%  - mnA_incon -- mean activation (aligned against ctxt), incongr,[32x3xN]
%  - Aprime := mnA_congr+mnA_incon, [32x3xN]
%  - accA -- A collapsed across blocks, [Nx6], 1:3=congr, 4:6=incongr
% % - seqA -- A as a function of S(t-1) and R(t-1), [Nx4], analog seqP
%  - mnI_congr -- mean bottom-up input (align with ctxt), congr, [32x3xN]
%  - mnI_incon -- mean bottom-up input (align against ctext), incongruent
%  - Iprime := mnI_congr+mnI_incon, [32x3xN]
%  - accI -- I collapsed across blocks, [Nx6], 1:3=congr, 4:6=incongr
% % - seqI -- I as a function of S(t-1) and R(t-1), [Nx4], analog seqP
%  - contextc -- zP averaged over 5 context-switch transitions, Nx(10+10)
%  - contexti -- zP, incongr. stim, context-switch transitions, Nx(10+10)
%  - otherc -- zP,   congruent stim, remaining 7+19 transitions,Nx(10+10)
%  - otheri -- zP, incongruent stim, remaining 7+19 transitions,Nx(10+10)
%
% Example (2 runs on the same stimulus sequence):
%  o1=PLModel1(descr,Mparams);o2=PLModel1(descr,Mparams);o=cat(3,o1,o2);
%  stats=PLM_stats(descr,o(:,1,:),o(:,2,:),o(:,3,:)) ;
%  sumstats = PLM_summary(stats) ; plot_PLM_dynamics(sumstats)
%
% See also PLM_SUMMARY, PLMODEL1, PLESTATS_BY_SBJ, MERGE_PLM_STATS.

% Original coding by Alexander Petrov, UC Irvine.
% $Revision: 1.0 $  $Date: 2003/03/05 10:15 $

%-- Check and unify the stimulus format
[N_trials,k,N_runs_idx] = size(stim) ;
%if (k<3) error('STIM must have at least 3 columns: Gabor,Cpeak,Noise.');end

resp = squeeze(resp) ;    % useful when resp=output(:,1,:), [9600x1xN]
[k,N_runs] = size(resp) ;
if (N_runs==N_runs_idx)
   N_runs_idx = [1:N_runs] ;
elseif (N_runs_idx==1)    % shorthand
   N_runs_idx = repmat(1,1,N_runs) ;
else
   error('STIM and RESP must involve the same number of runs.') ;
end

if (nargin==2)
   act = zeros(0,N_runs) ;  inp = zeros(0,N_runs) ;
else
   act = squeeze(act) ; inp = squeeze(inp) ;    % out(:,2,:),out(:,3,:)
end
if (nargin<5)  warnp = 1 ; end

%--  Analyze each run separately
stats = do_one(stim(:,1:3,1),resp(:,1),act(:,1),inp(:,1),warnp) ;
for k = 2:N_runs
   new_stats = do_one(stim(:,1:3,N_runs_idx(k)),...
                      resp(:,k),act(:,k),inp(:,k),warnp) ;
   stats = merge_PLM_stats(stats,new_stats) ;
end

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  stats = do_one(stim,resp,act,inp,warnp)

N_trials = size(stim,1) ;
BLOCK_SIZE = 300 ;                % hardwired, change if necessary
N_blocks = N_trials/BLOCK_SIZE ;
if (N_blocks~=fix(N_blocks))
   error('N_TRIALS is not a multiple of BLOCK_SIZE.') ;
end

%-- Unpack the arguments
Gabor_orient = stim(:,1) ;       % 1="-10", 2="+10"      [Nx1]
Cpeak = stim(:,2) ;              % 1=.23, 2=.15, 3=.10   [Nx1]
Noise_orient = stim(:,3) ;       % 1="-15", 2="+15"      [Nx1]

stim_congr = 2-(Gabor_orient==Noise_orient) ;  % 1="congr", 2="incon"
correct = (resp==Gabor_orient) ;               % 0="incorr",1="correct"

%-- Partition into blocks, identify contexts
bl = floor([0:N_trials-1]'./BLOCK_SIZE)+1 ;

if (warnp)
  [counts,vals] = xtab2(Noise_orient,bl) ;
  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
  %context{1} = find(counts(1,:)>0) ;  % negative background-noise orient
  %context{2} = find(counts(2,:)>0) ;  % positive background-noise orient
end

stats.group = Noise_orient(1) ;      % trial 1 representing block 1

%-- Probability correct by Block, Difficulty, and Congruence
[X,vals] = xtab(bl,Cpeak,stim_congr,correct) ;  % [32x3x2x2]
counts = squeeze(sum(X,4)) ;                    % [32x3x2]
P = X(:,:,:,2) ./ counts ;     % Probab. correct, [32x3x2]
P = max(0.01,min(P,0.99)) ;    % Avoid extreme frequencies

zP_congr = norminv(squeeze(P(:,:,1))) ;   % [32x3]
zP_incon = norminv(squeeze(P(:,:,2))) ;   % [32x3]
%zP_congr = (squeeze(P(:,:,1))) ;
%zP_incon = (squeeze(P(:,:,2))) ;

%-- Export the z-transformed probabilities
stats.zP_congr = zP_congr ;
stats.zP_incon = zP_incon ;
stats.dprime = zP_congr + zP_incon ;
stats.accZP = [mean(zP_congr) , mean(zP_incon)] ;    % [ccc,iii]

%-- Handle the mean output activations and mean bottom-up inputs
if (~isempty(act))
  %- Reverse the values in the negative context
  contextL = find(Noise_orient==1);
  act(contextL) = -act(contextL) ;
  inp(contextL) = -inp(contextL) ;
  %- Split by Block, Cpeak, and Stim_congr
  [X,vals] = split([act inp],[stim_congr,Cpeak,bl]) ;  % {192x1}
  M = aggreg(X,'mean') ;
  congr_idx = 1:3*N_blocks ;
  incon_idx = congr_idx + 3*N_blocks ;
  mnA_congr = +reshape(M(congr_idx,1),N_blocks,3) ;    % [32x3]
  mnA_incon = -reshape(M(incon_idx,1),N_blocks,3) ;    % [32x3]
  mnI_congr = +reshape(M(congr_idx,2),N_blocks,3) ;    % [32x3]
  mnI_incon = -reshape(M(incon_idx,2),N_blocks,3) ;    % [32x3]
  %- Export the means for each condition
  stats.mnA_congr = mnA_congr ;
  stats.mnA_incon = mnA_incon ;
  stats.Aprime = mnA_congr + mnA_incon ;
  stats.accA = [mean(mnA_congr) , mean(mnA_incon)] ;   % [ccc,iii]
  stats.mnI_congr = mnI_congr ;
  stats.mnI_incon = mnI_incon ;
  stats.Iprime = mnI_congr + mnI_incon ;
  stats.accI = [mean(mnI_congr) , mean(mnI_incon)] ;   % [ccc,iii]
end

%-- Call PLM_TRANSITIONS for fine-grain dynamics around transition points
trans = PLM_transitions(stim_congr,correct) ;
stats.contextc = trans.contextc' ;
stats.contexti = trans.contexti' ;
stats.otherc = (7/26).*trans.nightc' + (19/26).*trans.blockc' ;
stats.otheri = (7/26).*trans.nighti' + (19/26).*trans.blocki' ;

%-- Return STATS
%%%% End of subfunction DO_ONE

%%%%% End of PLM_STATS.M
