s03_compart_data.htb 10.8 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
%% SCRIPT DESCRIPTION
% This script is loads the results of a preliminary ICA and uses the
% standard deviation of the components to flag and remove data time that
% does not fit a well behaved model. This will help make the next
% components stronger. The log likelihood of the model is also flagged for
% reference. Sections marked for high SD often have low log likelihoods
% as well. It then prepares the required parameter files for the parrallel amicas ABC.
%
%% From Config          key_strswap         Description
%-----------------------------------------------------------------------------------------------------------------------------
%  data_path =          [in_path]           Path to input and output data files assuming cd = work_path
%  recur_sec =          [recur_sec]         Recurrence (sec) for artifact detection epoching (e.g. .3)
%  limit_sec =          [limit_sec]         Limits (sec) for artifact detection epoching (e.g. [-.3 .3])
14
15
16
17
18
19
%  sd_t_meth =          [sd_t_meth]         Method used for flagging epochs (e.g. 'q' (quantiles), or 'na' (default))
%  sd_t_vals  =         [sd_t_vals]         Percentage trim for confidence intervals during epoch standard deviation criteria (e.g. [.3 .7])
%  sd_t_o =             [sd_t_o]            z threshold for flagging epochs during standard deviation criteria (e.g. 6)
%  sd_t_f_meth =        [sd_t_f_meth]       Fixed method used for flagging epochs (should be 'fixed')
%  sd_t_f_vals  =       [sd_t_f_vals]       Percentage trim for confidence intervals during epoch standard deviation criteria (e.g. [.3 .7], leave empty for 'fixed')
%  sd_t_f_o =           [sd_t_f_o]          z threshold for flagging epochs during fixed standard deviation criteria (e.g. .2)
20
%  out_path =           [out_path]          Relative path to output data files assuming cd = work_path
21
22
23
%  min_gap_ms =         [min_gap_ms]        Minimum time (ms) to allow between periods marked for rejection 
%  amica_param_file =   [amica_param_file]  template amicadefs.param file to modify
%  amica_threads_s04 =  [amica_threads_s04] number of threads to use for running s04a, s04b, s04c amica scripts (Default: 8)
24
25
26
27
28
29
30
31
32
33

%% LOAD SA DATASET
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loads the batched EEG datafile and the AMICA output files. Since ICA is
% done on a channel subset the components need to be returned to the
% correct indicies.

% Load the EEG datafile
logging_log('INFO', 'Loading set file: [batch_dfn,_,-1]_sa.set...');
tic;
34
EEG = pop_loadset('filepath','[in_path]','filename','[batch_dfn,_,-1]_sa.set');
35
36
37
38
39
40
41
42
43
44
EEG = eeg_checkset( EEG );
logging_log('INFO', 'TIME TO: LOAD DATASET...');
toc

% Prepare returning components
EEG.icachansind = marks_label2index(EEG.marks.chan_info,{'manual','rank'},'indices','invert','on');
icatimeind = find(~EEG.marks.time_info(1).flags);

% Load AMICA Model 1
tic;
45
EEG.amica(1).models=loadmodout15('[out_path]/[batch_dfn,_,-1]_amicaout_init');
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
EEG.icaweights=EEG.amica(1).models.W;
EEG.icasphere=EEG.amica(1).models.S(1:EEG.amica(1).models.num_pcs,:);
EEG=eeg_checkset(EEG);
logging_log('INFO', 'TIME TO: LOAD ICA INFO...');
toc
%% PASS EEGDATA THROUGH THE AMICA MATRIX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Uses the ICA weights and the channel data to recreate the ICA components.

tic;
% In place of EEG checkset
tmpdata = EEG.data(EEG.icachansind,:);
tmpindices = find(~sum(isnan(tmpdata)));
EEG.icaact = zeros(size(EEG.icaweights,1), size(tmpdata,2)); EEG.icaact(:) = NaN;
EEG.icaact(:,tmpindices) = (EEG.icaweights*EEG.icasphere)*tmpdata(:,tmpindices);
clear tmpdata;

%% LOG LIKELIHOOD FLAGS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Creates marks based off of the log likelihood of the ICA model to the data.
% Create initial log likelihood marks

logl_init_flags=zeros(size(EEG.data(1,:)));
logl_init_flag_inds=marks_label2index(EEG.marks.time_info,{'manual'},'indices','invert','on');
logl_init_flags(logl_init_flag_inds)=(EEG.amica(1).models.Lt-min(EEG.amica(1).models.Lt))/(max(EEG.amica(1).models.Lt)-min(EEG.amica(1).models.Lt));

EEG.marks=marks_add_label(EEG.marks,'time_info', ...
	{'logl_init',[1,0.45,0],logl_init_flags});
logging_log('INFO', 'TIME TO: RECREATE ICA COMPONENTS...');
toc

% Window the continous Data
logging_log('INFO', 'Windowing the continous data...');
tic;
EEG=marks_continuous2epochs(EEG,'recurrence',[recur_sec],'limits',[limit_sec]);
logging_log('INFO', 'TIME TO: WINDOW THE CONTINUOUS DATA...');
toc
%% CALCULATE ICAACT SD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculates the IC standard Deviation by epoch window. Flags windows with
% too much standard deviation.

% Calculate IC sd by window
logging_log('INFO', 'Calculating the IC sd array for window criteria...');
tic;
epoch_inds=marks_label2index(EEG.marks.time_info,{'manual'},'indexes','invert','on');
92
[EEG,icaact_sd1_t]=chan_variance(EEG,'data_field','icaact', ...
93
94
95
96
97
		'epoch_inds',epoch_inds, ...
		'plot_figs','off');

% Create the windowing sd criteria
logging_log('INFO', 'Assessing window icaact sd distributions...')
98
[~,flag_t_inds]=marks_array2flags(icaact_sd1_t, ...
99
    'flag_dim','col', ...
100
    'init_method','[sd_t_meth]', ...
101
102
    'init_vals',[sd_t_vals], ...
    'init_crit',[sd_t_o], ...
103
104
105
    'flag_method','[sd_t_f_meth]', ...
    'flag_vals',[ [sd_t_f_vals] ], ...
    'flag_crit',[sd_t_f_o], ...
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
    'plot_figs','off');

% Update marks rejection structure
logging_log('INFO', 'Updating latflaginfo structure...');
icsd_epoch_flags=zeros(size(EEG.data(1,:,:)));
icsd_epoch_flags(1,:,epoch_inds(flag_t_inds))=1;
icsd_epoch_flags=padflags(EEG,icsd_epoch_flags,1,'value',.5);
EEG.marks=marks_add_label(EEG.marks,'time_info', ...
	{'ic_sd1',[.4,.6,.6],icsd_epoch_flags});
logging_log('INFO', 'TIME TO: UPDATE REJECTION STRUCTURE WITH IC SD MARKS...');
toc

% Concat the windowed data back to continuous
logging_log('INFO', 'Concatenating windowed data...');
tic;
EEG=marks_epochs2continuous(EEG);
EEG=eeg_checkset(EEG,'eventconsistency');
logging_log('INFO', 'TIME TO: CONCATENATE THE WINDOWED DATA INTO CONTINUOUS DATA...');
toc
%% FLAG GAPS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Flags small gaps between rejection flags so that there are not as many
% short segments to be analyzed

% Mark any flag gaps < 2 sec.
tic;
EEG=pop_marks_flag_gap(EEG,{'manual','mark_gap','ic_sd1'}, ...
                            [min_gap_ms],'mark_gap',[.8,.8,.8],'offsets',[0 0],'ref_point','both');

%%COMBINE MARKS STRUCTURE INTO MANUAL FLAGS...
EEG=pop_marks_merge_labels(EEG,'time_info',{'manual','ic_sd1','mark_gap'},'target_label','manual');
logging_log('INFO', 'TIME TO: FLAG GAPS AND COMBINE MARKS STRUCTURE INTO MANUAL FLAGS...');
toc
%% SAVE compart_data.set FILE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A saved copy is made before the file is purged to its ready for ICA
% version.

tic;
logging_log('INFO', 'Saving file: [batch_dfn,_,-1]_compart_data.set...');
146
EEG = pop_saveset( EEG, 'filename','[out_path]/[batch_dfn,_,-1]_compart_data.set');
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
logging_log('INFO', 'TIME TO: SAVE compart_data FILE...');
toc;
%% SAVE compart_data_purge FILE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% All EEG parts associated with ICA need to be cleared as models seem to
% adapt based on previous data. Old and New flagged channels/time need to
% be purged again in preparation of ICA.

% Clearing ICA Data
EEG.icaact = [];
EEG.icawinv = [];
EEG.icasphere = [];
EEG.icaweights = [];
EEG.icachansind = [];
EEG.amica = []; 
logging_log('INFO', 'ICA data Reset');

% Removing Flagged Channels
tic;
sprintf('%s','Purging flagged channels...\n');
EEG=pop_marks_select_data(EEG,'channel marks',[],'labels',{'manual','rank'},'remove','on');
EEG=pop_marks_select_data(EEG,'time marks',[],'labels',{'manual'},'remove','on');
EEG=eeg_checkset(EEG);
logging_log('INFO', 'TIME TO: PURGE DATA...');
toc;

% Saving File to enter next ICAs
logging_log('INFO', 'Saving the Param files');
sprintf('%s','Saving file: [batch_dfn,_,-1]_compart_data_purge.set...\n');
176
EEG = pop_saveset( EEG, 'filename','[out_path]/[batch_dfn,_,-1]_compart_data_purge.set');
177
178
179

% Save diagnostic arrays
try OCTAVE_VERSION;
180
    save('-mat7-binary','[out_path]/[batch_dfn,_,-1]_icaact_sd1_t.mat','icaact_sd1_t');
181
catch
182
    save('[out_path]/[batch_dfn,_,-1]_icaact_sd1_t.mat','icaact_sd1_t');
183
184
end

185
186
187
188
189
190
191
192
193
%% SAVE AMICA PARAMETER FILES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A parameter file need to be made for each of the ICAs running after this
% script. In This case AMICA A, B and C are run.

tic;
% A
% Saving AMICA Parameter File
swapstr=file_strswap('[amica_param_file]', ...
194
195
			'[repstr_fdt_fname]','[out_path]/[batch_dfn,_,-1]_compart_data_purge.fdt', ...
			'[repstr_outpath]','[out_path]/[batch_dfn,_,-1]_amicaout_A', ...
196
197
			'[repstr_nbchan]', num2str(EEG.nbchan), ...
			'[repstr_pnts]', sprintf('%12.0f',EEG.pnts), ...
198
            '[repstr_tpp]', '[amica_threads_s04]', ...
199
			'[repstr_pca]', num2str(EEG.nbchan));
200
fidparam=fopen('[out_path]/[batch_dfn,_,-1]_A.param','w');
201
202
203
204
205
206
fprintf(fidparam,'%s',swapstr);
fclose(fidparam);

% B
% Saving AMICA Parameter File
swapstr=file_strswap('[amica_param_file]', ...
207
208
			'[repstr_fdt_fname]','[out_path]/[batch_dfn,_,-1]_compart_data_purge.fdt', ...
			'[repstr_outpath]','[out_path]/[batch_dfn,_,-1]_amicaout_B', ...
209
210
			'[repstr_nbchan]', num2str(EEG.nbchan), ...
			'[repstr_pnts]', sprintf('%12.0f',EEG.pnts), ...
211
            '[repstr_tpp]', '[amica_threads_s04]', ...
212
			'[repstr_pca]', num2str(EEG.nbchan));
213
fidparam=fopen('[out_path]/[batch_dfn,_,-1]_B.param','w');
214
215
216
217
218
219
fprintf(fidparam,'%s',swapstr);
fclose(fidparam);

% C
% Saving AMICA Parameter File
swapstr=file_strswap('[amica_param_file]', ...
220
221
			'[repstr_fdt_fname]','[out_path]/[batch_dfn,_,-1]_compart_data_purge.fdt', ...
			'[repstr_outpath]','[out_path]/[batch_dfn,_,-1]_amicaout_C', ...
222
223
			'[repstr_nbchan]', num2str(EEG.nbchan), ...
			'[repstr_pnts]', sprintf('%12.0f',EEG.pnts), ...
224
            '[repstr_tpp]', '[amica_threads_s04]', ...
225
			'[repstr_pca]', num2str(EEG.nbchan));
226
fidparam=fopen('[out_path]/[batch_dfn,_,-1]_C.param','w');
227
228
229
230
231
232
233
fprintf(fidparam,'%s',swapstr);
fclose(fidparam);

%% Create the Amica Output Folder
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Needed for MPI jobs on cluster. Otherwise race effect will occur.

234
235
236
[~] = rmdir('[out_path]/[batch_dfn,_,-1]_amicaout_A','s');
[~] = rmdir('[out_path]/[batch_dfn,_,-1]_amicaout_B','s');
[~] = rmdir('[out_path]/[batch_dfn,_,-1]_amicaout_C','s');
237
238


239
240
241
mkdir '[out_path]/[batch_dfn,_,-1]_amicaout_A';
mkdir '[out_path]/[batch_dfn,_,-1]_amicaout_B';
mkdir '[out_path]/[batch_dfn,_,-1]_amicaout_C';
242
243
244
245
246
247
248
249

logging_log('INFO', 'TIME TO: CREATE AMICA PARAMETER FILES...');
toc

% FINISH
logging_log('INFO', '***********COMPART DATA COMPLETE!**********');
logging_log('INFO', 'Scheduler: [scheduler]');
print_chan_sample(EEG);