function  D = import_PLE2_data(fname,fid)

%IMPORT_PLE2_DATA --  Import the raw data of the PercLearn Experiment 2.
%
% D = import_PLE2_data(fname)
% D = import_PLE2_data(fname,verbose)
%
% FNAME is the name of a raw data file in ASCII format. If not supplied,
% the default fname='raw_data.dat'. It should contain 21 columns as 
% documented in PLEXP2.
% Set VERBOSE to 0 to suppress messages, 1(=default) to have them
% printed in the console window, or pass a fileid of an open ASCII file.
%
% IMPORT_PLE2_DATA reads the data file and returns a structure array 
% D with one element per subject and the following fields:
%   1. sbj
%   2. group
%   3. stim -- a matrix with one row per trial and 7 columns:
%       1=Gabor_orient : 1=-10 or 2=+10 deg, see MAKE_PLE2_STIM
%       2=difficulty   : 1=0.23, 2=0.15, 3=0.10; 4=0.50=demo
%       3=noise_orient : 1=-15 or 2=+15 deg, see PLEXP2_PARAMS
%       4=position     : 1=-111=above fixation, 2=+111=below fixation
%       5=NC_fnumber   : 001-036, see CACHE_PLE2_BACKGR
%       6=NC_idx       : 1-300, see CACHE_PLE2_BACKGR & MAKE_PLE2_STIM
%       7=LUT_idx : 0-255, see the 'sanity check' in DO_PLE2_TRIAL
%   4. data -- a matrix with one row per trial and 10 columns:
%       1=session(day)
%       2=block        : 0=demo, 1-4
%       3=Gabor_degr   : -10 or +10 degrees    , cf stim(:,1) above
%       4=Gabor_Cpeak  : [0.23 0.15 0.10, 0.50], cf stim(:,2) above
%       5=noise_degr   : +15 or -15 degrees    , cf stim(:,3) above
%       6=offset (dy)  : -111 or +111 pixels   , cf stim(:,4) above
%       7=response     : 1="Left", 2="Right", 0=bad_key
%       8=correct      : 0=incorrect, 1=correct
%       9=RT           : in milliseconds from the onset of the stimulus
%      10=exit_code    : 0=successful, 1=bad_key, 2=missed_frame
%   5. miss_idx -- a column-vector of indices in STIM and DATA
%   6. demo_idx
%   7. skip_idx -- =union(miss_idx,demo_idx)
%
%  Example 1 (Assuming all data are in a single file 'raw_data.dat'):
%    D=import_PLE2_data ; st=PLE2_stats(D) ; plot_PLE2_stats(D) ;
%
%  Example 2 (Read data for subjects 1,2,3, and 20):
%    Ss = [1,2,3,20] ; N = length(Ss) ;
%    for k=1:N D(k)=import_PLE2_data(sprintf('sbj%03d.dat',Ss(k)),0);end
%    st=PLE2_stats(D);for k=1:4 subplot(2,2,k);plot(st(k).dprime,'.-');end
%
%  Example 3 (Invoke MS-DOS COPY command to append the data files):
%    !copy C:\work\PLExp2\data\sbj0??.dat C:\work\PLExp2\raw_data.dat
%    D = import_PLE2_data('C:\work\PLExp2\raw_data.dat')
%
%  See also PLE2_STATS, PLEXP2, PLEXP2A MAKE_PLE2_STIM.

%  Original coding by Alexander Petrov, UC Irvine.
%  $Revision: 1.1 $  $Date: 2003/09/08 17:55 $

% Format of the raw data file---see PLEXP2 and its fubfun RECORD_DATA
%    1. subject number
%    2. subject group
%    3. session (day)
%    4. block
%    5. overall trial number
%    ---
%    6. time stamp
%    7. noise orientation (1=negative="Left", 2=positive="Right")
%    8. actual noise orientation (in degrees)
%    9. NCache file number (typically b/n 001:036, 000=newly_generated)
%   10. NCache index (1:300)
%   11. Gabor orientation (1=negative, 2=positive)
%   12. actual Gabor orientation (in degrees)
%   13. difficulty level (max Gabor contrast, 1=0.23,2=0.15,3=0.10,4=demo)
%   14. actual Gabor peak contrast (0.23,0.15,0.10)
%   15. position (1=above fixation, 2=below fixation)
%   16. actual position (in pixels relative to fixation at center screen)
%   17. sanity check -- actual LUT index of STIM(P.Gabor_argmax)
%   ---
%   18. response (1="Left", 2="Right", 0=bad_key)
%   19. correct (1 if response==Gabor_orientation, 0 otherwise)
%   20. response time (in msec from the onset of the stimulus)
%   21. exit code (0='successful trial', 1='bad key', 2='missed frame')
%

if (nargin < 2)  fid = 1 ; end    % verbose messages go to stdout

%%%%%%%%%%%%%%%%%  Read the raw ASCII file

if (nargin < 1)  fname = 'raw_data.dat' ; end   % default data fname

fprintf(fid,'\nLoading ASCII file %s.  Please wait...\n',fname) ;
raw_data = load(fname,'-ascii') ;
[rows cols] = size(raw_data) ;
fprintf(fid,'\n%d rows and %d columns loaded.\n',rows,cols) ;

%-- Do a (superficial) consistency check
if (find(isnan(raw_data)))
  warning ('NaNs found in raw_data. Check the raw ASCII file.') ;
else
  fprintf(fid,'The RAW_DATA matrix contains no missing values.\n') ;
end;


%%%%%%%%%%%%%%%%%  Slice the big matrix into variables
fprintf(fid,'\nSlicing the big matrix into variables, columnwise:\n') ;

%-- Description of the subject and session: columns 1 through 5
sbj   = raw_data(:,1) ;     % subject number
group = raw_data(:,2) ;     % group number
session = raw_data(:,3) ;   % session (day)
block = raw_data(:,4) ;     % block
T = raw_data(:,5) ;         % overall trial number
fprintf(fid,'  Extracted variables SBJ, GROUP, SESSION, BLOCK, and T.\n');

%-- Absolute time since the beginning of experiment: column 6
clock = raw_data(:,6) ;
fprintf(fid,'  Extracted variable  CLOCK.\n') ;

%-- Description of the background noise: columns 7-10
noise_orient = raw_data(:,7) ;      % 1=negative, 2=positive
noise_degr = raw_data(:,8) ;        % degrees relative to vertical
fprintf(fid,'  Extracted variables NOISE_ORIENT and NOISE_DEGR.\n') ;
NCache_fnumber = raw_data(:,9) ;    % 001-036, see the CD-ROM
NCache_idx = raw_data(:,10) ;       % NOISE_CACHE(:,:,NCache_idx)
fprintf(fid,'  Extracted variables NCACHE_FNUMBER and NCACHE_IDX.\n') ;

%-- Description of the Gabor target: columns 11-14
Gabor_orient = raw_data(:,11) ;     % 1=negative, 2=positive
Gabor_degr = raw_data(:,12) ;       % degrees relative to vertical
fprintf(fid,'  Extracted variables GABOR_ORIENT and GABOR_DEGR.\n') ;
difficulty = raw_data(:,13) ;       % difficulty level: 1,2,3,(4)
Gabor_Cpeak = raw_data(:,14) ;      % peak contrast of the Gabor patch
fprintf(fid,'  Extracted variables DIFFICULTY and GABOR_CPEAK.\n') ;

%-- Retinal position and "sanity check" value for one pixel: cols 15-17
position = raw_data(:,15) ;         % 1=above fixation, 2=below fix.
offset = raw_data(:,16) ;           % pixels above or below fixation cross
LUT_idx = raw_data(:,17) ;          % actual STIM(P.Gabor_argmax)
fprintf(fid,'  Extracted variables POSITION, OFFSET, and LUT_IDX.\n') ;

%-- Responses obtained from the subject: columns 18-21
resp = raw_data(:,18) ;     % 1='negative', 2='positive', 0='bad key'
corr = raw_data(:,19) ;     % correct = (resp==Gabor_orient)
RT   = raw_data(:,20) ;     % response time [msec] from stimulus onset
exitcode = raw_data(:,21) ; % 0='successful', 1='bad key', 2='missed frame'
fprintf(fid,'  Extracted variables RESP, CORR, RT, and EXITCODE.\n') ;

%-- Clear the big raw matrix
clear raw_data rows cols ;
fprintf(fid,'Cleared the big matrix RAW_DATA.\n\n') ;

% The following 21 variables now exist in the local workspace:
% sbj, group, session, block, T, clock,
% noise_orient, noise_degr, NCache_fnumber, NCache_idx, 
% Gabor_orient, Gabor_degr, difficulty, Gabor_Cpeak, 
% position, offset, LUT_idx, resp, corr, RT, exitcode.


%%%%%%%%%%%%%%%%%  Partition data by subject and pack into a structure

%-- 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(fid,'Processing %5d trials for participant %2d... ',...
           length(idx),curr_sbj) ;
   D(k) = do_one(curr_sbj,group(idx),session(idx),block(idx),T(idx),...
                 noise_orient(idx),noise_degr(idx),...
                 NCache_fnumber(idx),NCache_idx(idx),...
                 Gabor_orient(idx),Gabor_degr(idx),...
                 difficulty(idx),Gabor_Cpeak(idx),...
                 position(idx),offset(idx),LUT_idx(idx),...
                 resp(idx),corr(idx),RT(idx),exitcode(idx)) ;
   fprintf(fid,'Done.\n') ;
end

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  D = do_one(sbj,group,session,block,T,...
                 noise_orient,noise_degr,  NCache_fnumber,NCache_idx,...
                 Gabor_orient,Gabor_degr,  difficulty,Gabor_Cpeak,...
                 position,offset,LUT_idx,  resp,corr,RT,exitcode)

%-- Define some constants describing the experimental design
Ph1_sessions = 9 ;   % 9 sessions [1:9]   in phase 1, see PLExp2
Ph2_sessions = 2 ;   % 2 sessions [10:11] in phase 2, see PLExp2a
Ph12_sessions = Ph1_sessions + Ph2_sessions ;
Ph3_sessions = 1 ;   % 1 session  [12]    in phase 3, see PLExp2b
Ph12_trials_per_block = 300 ;   % 12 stimuli * 25 replications
Ph3_trials_per_block = 312 ;   % 24 stimuli * 13 replications

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

%-- Identify the valid trials and the demo trials
% Exitcodes (see DO_PLE2_TRIAL): 0='OK',1='bad key',2='missed video frame'
miss_idx = find(exitcode~=0) ;
demo_idx = find(block==0) ;
skip_idx = union(miss_idx,demo_idx) ;
ok = [1:length(block)]' ; ok(skip_idx) = [] ;   % valid experim. trials

%-- Partition sessions by phases
phase1 = intersect(ok,find(session<=Ph1_sessions)) ;
phase2 = intersect(ok,find(session>Ph1_sessions & session<=Ph12_sessions));
phase12= intersect(ok,find(session<=Ph12_sessions)) ;      % [1:9,10:11]
phase3 = intersect(ok,find(session>Ph12_sessions & ...
                           session<=Ph3_sessions+Ph12_sessions)) ;  % 12

%-- Check that each session has four blocks of 300 (312) trials each
[counts,vals] = xtab2(session(phase12),block(phase12)) ;
if (length(vals{2})~=4)
   warning(' Nonstandard number of blocks per session.') ;
elseif (any(counts(:)~=Ph12_trials_per_block))
   warning(' Nonstandard number of trials per block (phase 1 or 2).');
end

if (~isempty(phase3))     % this subject has not reached that far yet
  [counts,vals] = xtab2(session(phase3),block(phase3)) ;
  if (length(vals{2})~=4)
     warning(' Nonstandard number of blocks per session.') ;
  elseif (any(counts(:)~=Ph3_trials_per_block))
     warning(' Nonstandard number of trials per block (phase 3).');
  end
end

%-- Check that all trials within a block have the same NOISE_DEGR
bl = block + 4*(session-1) ;   % the main time variable

[counts,vals] = xtab2(noise_degr(phase1),bl(phase1)) ;
if (length(vals{1})~=2)        % 1 and 2
   warning(' Confusing noise orientations (phase 1).') ;
elseif  (any(prod(counts)~=0)) % one of the two frequencies should be 0
   warning(' Block(s) with mixed noise orientations.');
end

if (~isempty(phase2))          % this sbj has not reached that far yet
  [counts,vals] = xtab2(noise_degr(phase2),bl(phase2)) ;
  if (length(vals{1})>2)       % should be [3,4] or [3] or [4]
     warning(' Confusing noise orientations (phase 2).') ;
  elseif  (length(vals{1})==2 & any(prod(counts)~=0))  
     warning(' Block(s) with mixed noise orientations.');
  end
end


%-- Check for one-to-one correspondence b/n blocks and NCACHE_FNUMBERs
[counts,vals] = xtab2(NCache_fnumber(phase1),bl(phase1)) ;
if ((length(vals{1})~=length(vals{2})) | ...
    (min(max(counts))~=Ph12_trials_per_block))
   warning(' Imperfect correspondence b/n blocks & NOISE_CACHE files.');
   %xtab2(NCache_fnumber(ok),bl(ok))
end

%-- Arrange the STIM and DATA matrices
stim = [Gabor_orient,difficulty,noise_orient,...
        position,NCache_fnumber,NCache_idx,LUT_idx] ;
data = [session,block,Gabor_degr,Gabor_Cpeak,noise_degr,...
        offset,resp,corr,RT,exitcode];

%-- Pack everything into a big structure and return
D.sbj = sbj ;
D.group = group(1) ;
D.stim = stim ;
D.data = data ;
D.miss_idx = miss_idx ;
D.demo_idx = demo_idx ;
D.skip_idx = skip_idx ;

%--- Return D
%%%% End of subfunction DO_ONE

%%%%% End of IMPORT_PLE2_DATA.M
