% Run independent components analysis (ICA) on MEG data,
% classify components related to blinks, eye-movements, 
% and pulse artefacts, then remove artefact component(s) 
% related to, e.g., blinks.
%
% This script assumes you've run the session 2 script
%   cbu_meeg_masterclass2
% or at least that you have a directory for the project
% directory in your imaging space.
%

%--Set Parameters------------------------------------------

% Working directory 
wd = '/home/ms02/MEGMasterclass/Marie/'; %% !!!CHANGE TO YOUR OWN !!!
cd(wd)

% Extra bit for MEG masterclass - copy ICAed files from here:
%ica_sourcedir = '/imaging/jt03/public/demo/';
ica_sourcedir = '/imaging/jt03/public/demo/icadata/';

% Add relevant directories to your path:
addpath /imaging/local/spm/spm5/cbu_updates/
addpath /imaging/local/meg_misc/  
addpath /imaging/local/spm_eeglab/  % ica + other EEGLAB functions
addpath /imaging/jt03/public/demo/  % functions not yet uploaded to meg_misc!

% Processing flags & variables for ICA
do_runica    = 0; % do ICA (takes time!)?
do_copyica   = 1; % just copy ICAed files instead!
mem_flag     = 0; % flag if run into memory problems!

icaflag      = [1 1 1 0];  % don't run on grms
ica_Npc      = 32;  % Number of ICs to extract (PCA run first)
ica_samppct  = 1;  % Proportion of data to give to ICA
ica_class    = {'blink','horiz','pulse'}; % Artefact(s) to classify
ica_chans    = {'veog','heog','ecg'};     % Reference channels for correlations
ica_rem      = {'blink'};                 % Artefact(s) to remove 

%%%%%%%%%% Different types of sensor
typ          = {'mags' 'grds' 'eeg' 'grms'};
Ntyp         = length(typ);
eegflag      = [0 0 1 0];  % Binary flag for which typ is EEG

%%%%%%%%%%% subject-specific parameters
dosubs       = 1:8;
Nsub         = length(dosubs);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%--COPY ICAed DATA-----------------------------------------
% Extra bit for MEG masterclass. Instead of running ICA on 
% all 8 subjects, just copy ICAed files:

if do_copyica
	for sub=1:Nsub
		swd = sprintf('%ssub%d',wd,sub);
		try   cd(swd)
    catch mkdir(swd)
          cd(swd)
    end
    
		ses=1;
	
		% Assign filenames:
		basefile=sprintf('dsub%d_ses%d_trans',sub,ses);
		icafile=sprintf('ica_%s',basefile);
		[p fstem ext]=fileparts(basefile);
		
		for typn = 1:Ntyp
      if icaflag(typn)
        % Copy ica*.mat file to local directory:
        fn=sprintf('%s-%s.mat',icafile,typ{typn});
        if exist(fn)~=2
          disp(sprintf('Copying file %s...',fn))
          [cpfail tmp]=unix(sprintf('! cp %s/%s ./',ica_sourcedir,fn));
          if cpfail
            error('Copy failed!')
          end
          disp(sprintf('Done.'))
        end
        % Ensure .dat file is present (copy over if not):
        D=spm_eeg_ldata(fn);
        D.path=swd;
        D.fnamedat=sprintf('%s-%s.dat',fstem,typ{typn});
        if exist(D.fnamedat)~=2
          disp(sprintf('Copying .dat file ...'))
          [cpfail tmp]=unix(sprintf('! cp %s/%s ./',ica_sourcedir,D.fnamedat));
          if cpfail
            error('Copy failed!')
          end
          disp(sprintf('Done.'))
        end
        % Remove any previous classifications:
        D.ica.class=[];
        save(D.fname,'D');

        clear D
      end % icaflag
		end % typn
	end % sub
end % do_copyica


% End of extra bit
%----------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%--Processing begins here!---------------------------------

for sub = 1:Nsub;
	subnum    = dosubs(sub);
	Nses(sub)=1;

	cd(wd)
	swd = sprintf('sub%d',subnum);
	try cd(swd); catch eval(sprintf('!mkdir %s',swd)); cd(swd); end

	for ses = 1:Nses(sub)

		% We are using the files that have been imported to SPM,
		% downsampled, and split (separated into mags and grds)
		% Thus, basefile for subject 1 is 'dsub1_ses1_trans'
		basefile=sprintf('dsub%d_ses%d_trans',sub,ses);

		for typn = 1:Ntyp


%%%%% Artefact compensation using ICA

			if icaflag(typn)

%--Run ICA--------------------------------------------------
%
% ICA takes somewhere between 5-10 minutes on each file,
%  longer if pca=65 or if using non-split mags+grds file.
%  I would recommend running this bit of code in a batch over
%  subjects, say overnight, then running the next bit 
%  (classify and remove artefacts) in the morning!

			if do_runica % If...else...end is Extra bit for MEG masterclass
	
				clear S; S.D = sprintf('%s-%s.mat',basefile,typ{typn});
				S.samppct = ica_samppct;
				S.args    = {'extended',1,...
										 'maxsteps',800,...
										 'pca',ica_Npc};
				if ~mem_flag
					D = spm_eeglab_runica(S);
				elseif mem_flag
					save('tmp.mat','S');
					! matlab -nodesktop -nosplash -r "spm eeg; load tmp.mat; spm_eeglab_runica(S); exit"
					try delete('tmp.mat'); end
					D.fname=['ica_' S.D];
				end %mem_flag
				
			else % Extra bit for MEG masterclass
					
					D.fname=sprintf('ica_%s-%s.mat',basefile,typ{typn});
					
			end % Extra bit for MEG masterclass

					
%--Classify and Remove Artefact ICs-------------------------
%
% View raw data on blink- and pulse-sensitive channels,
%  compute IC activations, correlate with physiological 
%  signals, plot corr Z-scores, plot topographies,
%  classify ICs as blink, horizontal eye mov'ts, pulse,
%  and remove selected artefact components.

				clear S; S.D = D.fname;
				S.mode     = {'both'};  % classify + remove
				S.artlabs  = ica_class;
				S.chanlabs = ica_chans;
				S.remlabs  = ica_rem;
				S.remfname = ['rem_' S.D];

				if ~mem_flag
					D = meg_ica_artefact(S);
				elseif mem_flag
					save('tmp.mat','S');
					! matlab -nodesktop -nosplash -r "spm eeg; load tmp.mat; meg_ica_artefact(S); exit"
					try delete('tmp.mat'); end
					D.fname=S.remfname;
				end %mem_flag

%%%%% Done!

			end % if icaflag
			
			% The rest of pre-processing filter, epoch, average, etc.
			%  would go here.
	
			
		end % typn
	end % ses
end % sub