/
ica.py
2719 lines (2343 loc) · 113 KB
/
ica.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
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
146
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
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# -*- coding: utf-8 -*-
#
# Authors: Denis A. Engemann <denis.engemann@gmail.com>
# Alexandre Gramfort <alexandre.gramfort@inria.fr>
# Juergen Dammers <j.dammers@fz-juelich.de>
#
# License: BSD (3-clause)
from inspect import isfunction
from collections import namedtuple
from copy import deepcopy
from numbers import Integral
from time import time
import os
import json
import numpy as np
from scipy import linalg
from .ecg import (qrs_detector, _get_ecg_channel_index, _make_ecg,
create_ecg_epochs)
from .eog import _find_eog_events, _get_eog_channel_index
from .infomax_ import infomax
from ..cov import compute_whitener
from .. import Covariance, Evoked
from ..io.pick import (pick_types, pick_channels, pick_info,
_picks_to_idx, _get_channel_types, _DATA_CH_TYPES_SPLIT)
from ..io.write import (write_double_matrix, write_string,
write_name_list, write_int, start_block,
end_block)
from ..io.tree import dir_tree_find
from ..io.open import fiff_open
from ..io.tag import read_tag
from ..io.meas_info import write_meas_info, read_meas_info
from ..io.constants import FIFF
from ..io.base import BaseRaw
from ..io.eeglab.eeglab import _get_info, _check_load_mat
from ..epochs import BaseEpochs
from ..viz import (plot_ica_components, plot_ica_scores,
plot_ica_sources, plot_ica_overlay)
from ..viz.ica import plot_ica_properties
from ..viz.topomap import _plot_corrmap
from ..channels.channels import _contains_ch_type, ContainsMixin
from ..io.write import start_file, end_file, write_id
from ..utils import (check_version, logger, check_fname, verbose,
_reject_data_segments, check_random_state, _validate_type,
compute_corr, _get_inst_data, _ensure_int,
copy_function_doc_to_method_doc, _pl, warn, Bunch,
_check_preload, _check_compensation_grade, fill_doc,
_check_option, _PCA, deprecated)
from ..utils.check import _check_all_same_channel_names
from ..fixes import _get_args, _safe_svd
from ..filter import filter_data
from .bads import _find_outliers
from .ctps_ import ctps
from ..io.pick import pick_channels_regexp
__all__ = ('ICA', 'ica_find_ecg_events', 'ica_find_eog_events',
'get_score_funcs', 'read_ica', 'run_ica', 'read_ica_eeglab')
def _make_xy_sfunc(func, ndim_output=False):
"""Aux function."""
if ndim_output:
def sfunc(x, y):
return np.array([func(a, y.ravel()) for a in x])[:, 0]
else:
def sfunc(x, y):
return np.array([func(a, y.ravel()) for a in x])
sfunc.__name__ = '.'.join(['score_func', func.__module__, func.__name__])
sfunc.__doc__ = func.__doc__
return sfunc
# makes score funcs attr accessible for users
def get_score_funcs():
"""Get the score functions.
Returns
-------
score_funcs : dict
The score functions.
"""
from scipy import stats
from scipy.spatial import distance
score_funcs = Bunch()
xy_arg_dist_funcs = [(n, f) for n, f in vars(distance).items()
if isfunction(f) and not n.startswith('_')]
xy_arg_stats_funcs = [(n, f) for n, f in vars(stats).items()
if isfunction(f) and not n.startswith('_')]
score_funcs.update({n: _make_xy_sfunc(f)
for n, f in xy_arg_dist_funcs
if _get_args(f) == ['u', 'v']})
score_funcs.update({n: _make_xy_sfunc(f, ndim_output=True)
for n, f in xy_arg_stats_funcs
if _get_args(f) == ['x', 'y']})
return score_funcs
def _check_for_unsupported_ica_channels(picks, info, allow_ref_meg=False):
"""Check for channels in picks that are not considered valid channels.
Accepted channels are the data channels
('seeg','ecog','eeg', 'hbo', 'hbr', 'mag', and 'grad'), 'eog' and 'ref_meg'
This prevents the program from crashing without
feedback when a bad channel is provided to ICA whitening.
"""
types = _DATA_CH_TYPES_SPLIT + ('eog',)
types += ('ref_meg',) if allow_ref_meg else ()
chs = _get_channel_types(info, picks, unique=True, only_data_chs=False)
check = all([ch in types for ch in chs])
if not check:
raise ValueError('Invalid channel type%s passed for ICA: %s.'
'Only the following types are supported: %s'
% (_pl(chs), chs, types))
@fill_doc
class ICA(ContainsMixin):
u"""M/EEG signal decomposition using Independent Component Analysis (ICA).
This object estimates independent components from :class:`mne.io.Raw`,
:class:`mne.Epochs`, or :class:`mne.Evoked` objects. Components can
optionally be removed (for artifact repair) prior to signal reconstruction.
.. warning:: ICA is sensitive to low-frequency drifts and therefore
requires the data to be high-pass filtered prior to fitting.
Typically, a cutoff frequency of 1 Hz is recommended.
Parameters
----------
n_components : int | float | None
Number of principal components (from the pre-whitening PCA step) that
are passed to the ICA algorithm during fitting. If :class:`int`, must
not be larger than ``max_pca_components``. If :class:`float` between 0
and 1, the number of components with cumulative explained variance less
than ``n_components`` will be used. If ``None``, ``max_pca_components``
will be used. Defaults to ``None``; the actual number used when
executing the :meth:`ICA.fit` method will be stored in the attribute
``n_components_`` (note the trailing underscore).
max_pca_components : int | None
Number of principal components (from the pre-whitening PCA step) that
are retained for later use (i.e., for signal reconstruction in
:meth:`ICA.apply`; see the ``n_pca_components`` parameter). Use this
parameter to reduce the dimensionality of the input data via PCA before
any further processing is performed. If ``None``, no dimensionality
reduction occurs and ``max_pca_components`` will equal the number of
channels in the :class:`mne.io.Raw`, :class:`mne.Epochs`, or
:class:`mne.Evoked` object passed to :meth:`ICA.fit`. Defaults to
``None``.
n_pca_components : int | float | None
Total number of components (ICA + PCA) used for signal reconstruction
in :meth:`ICA.apply`. At minimum, at least ``n_components`` will be
used (unless modified by ``ICA.include`` or ``ICA.exclude``). If
``n_pca_components > n_components``, additional PCA components will be
incorporated. If :class:`float` between 0 and 1, the number is chosen
as the number of *PCA* components with cumulative explained variance
less than ``n_components`` (without accounting for ``ICA.include`` or
``ICA.exclude``). If :class:`int` or :class:`float`, ``n_components_ ≤
n_pca_components ≤ max_pca_components`` must hold. If ``None``,
``max_pca_components`` will be used. Defaults to ``None``.
noise_cov : None | instance of Covariance
Noise covariance used for pre-whitening. If None (default), channels
are scaled to unit variance ("z-standardized") prior to the whitening
by PCA.
%(random_state)s
As estimation can be non-deterministic it can be useful to fix the
random state to have reproducible results.
method : {'fastica', 'infomax', 'picard'}
The ICA method to use in the fit method. Use the fit_params argument to
set additional parameters. Specifically, if you want Extended Infomax,
set method='infomax' and fit_params=dict(extended=True) (this also
works for method='picard'). Defaults to 'fastica'. For reference, see
[1]_, [2]_, [3]_ and [4]_.
fit_params : dict | None
Additional parameters passed to the ICA estimator as specified by
`method`.
max_iter : int
Maximum number of iterations during fit. Defaults to 200. The actual
number of iterations it took :meth:`ICA.fit` to complete will be stored
in the ``n_iter_`` attribute.
allow_ref_meg : bool
Allow ICA on MEG reference channels. Defaults to False.
.. versionadded:: 0.18
%(verbose)s
Attributes
----------
current_fit : str
Flag informing about which data type (raw or epochs) was used for the
fit.
ch_names : list-like
Channel names resulting from initial picking.
n_components_ : int
If fit, the actual number of PCA components used for ICA decomposition.
pre_whitener_ : ndarray, shape (n_channels, 1) or (n_channels, n_channels)
If fit, array used to pre-whiten the data prior to PCA.
pca_components_ : ndarray, shape (`max_pca_components`, n_channels)
If fit, the PCA components.
pca_mean_ : ndarray, shape (n_channels,)
If fit, the mean vector used to center the data before doing the PCA.
pca_explained_variance_ : ndarray, shape (`max_pca_components`,)
If fit, the variance explained by each PCA component.
mixing_matrix_ : ndarray, shape (`n_components_`, `n_components_`)
If fit, the whitened mixing matrix to go back from ICA space to PCA
space.
It is, in combination with the `pca_components_`, used by
:meth:`ICA.apply` and :meth:`ICA.get_components` to re-mix/project
a subset of the ICA components into the observed channel space.
The former method also removes the pre-whitening (z-scaling) and the
de-meaning.
unmixing_matrix_ : ndarray, shape (`n_components_`, `n_components_`)
If fit, the whitened matrix to go from PCA space to ICA space.
Used, in combination with the `pca_components_`, by the methods
:meth:`ICA.get_sources` and :meth:`ICA.apply` to unmix the observed
data.
exclude : array-like of int
List or np.array of sources indices to exclude when re-mixing the data
in the :meth:`ICA.apply` method, i.e. artifactual ICA components.
The components identified manually and by the various automatic
artifact detection methods should be (manually) appended
(e.g. ``ica.exclude.extend(eog_inds)``).
(There is also an `exclude` parameter in the :meth:`ICA.apply` method.)
To scrap all marked components, set this attribute to an empty list.
info : None | instance of Info
The measurement info copied from the object fitted.
n_samples_ : int
The number of samples used on fit.
labels_ : dict
A dictionary of independent component indices, grouped by types of
independent components. This attribute is set by some of the artifact
detection functions.
n_iter_ : int
If fit, the number of iterations required to complete ICA.
Notes
-----
A trailing ``_`` in an attribute name signifies that the attribute was
added to the object during fitting, consistent with standard scikit-learn
practice.
Prior to fitting and applying the ICA, data is whitened (de-correlated and
scaled to unit variance, also called sphering transformation) by means of
a Principal Component Analysis (PCA). In addition to the whitening, this
step introduces the option to reduce the dimensionality of the data, both
prior to fitting the ICA (with the ``max_pca_components`` parameter) and
prior to reconstructing the sensor signals (with the ``n_pca_components``
parameter). In this way, we separate the question of how many ICA
components to estimate from the question of how much to reduce the
dimensionality of the signal. For example: by setting high values for
``max_pca_components`` and ``n_pca_components``, relatively little
dimensionality reduction will occur when the signal is reconstructed,
regardless of the value of ``n_components`` (the number of ICA components
estimated).
.. note:: Commonly used for reasons of i) computational efficiency and
ii) additional noise reduction, it is a matter of current debate
whether pre-ICA dimensionality reduction could decrease the
reliability and stability of the ICA, at least for EEG data and
especially during preprocessing [5]_. (But see also [6]_ for a
possibly confounding effect of the different whitening/sphering
methods used in this paper (ZCA vs. PCA).)
On the other hand, for rank-deficient data such as EEG data after
average reference or interpolation, it is recommended to reduce
the dimensionality (by 1 for average reference and 1 for each
interpolated channel) for optimal ICA performance (see the
`EEGLAB wiki <eeglab_wiki_>`_).
If you are migrating from EEGLAB and intend to reduce dimensionality via
PCA, similarly to EEGLAB's ``runica(..., 'pca', n)`` functionality, simply
pass ``max_pca_components=n``, while leaving ``n_components`` and
``n_pca_components`` at their respective default values.
Caveat! If supplying a noise covariance, keep track of the projections
available in the cov or in the raw object. For example, if you are
interested in EOG or ECG artifacts, EOG and ECG projections should be
temporally removed before fitting ICA, for example::
>> projs, raw.info['projs'] = raw.info['projs'], []
>> ica.fit(raw)
>> raw.info['projs'] = projs
Methods currently implemented are FastICA (default), Infomax, and Picard.
Standard Infomax can be quite sensitive to differences in floating point
arithmetic. Extended Infomax seems to be more stable in this respect,
enhancing reproducibility and stability of results; use Extended Infomax
via ``method='infomax', fit_params=dict(extended=True)``. Allowed entries
in ``fit_params`` are determined by the various algorithm implementations:
see :class:`~sklearn.decomposition.FastICA`, :func:`~picard.picard`,
:func:`~mne.preprocessing.infomax`.
Reducing the tolerance (set in `fit_params`) speeds up estimation at the
cost of consistency of the obtained results. It is difficult to directly
compare tolerance levels between Infomax and Picard, but for Picard and
FastICA a good rule of thumb is ``tol_fastica == tol_picard ** 2``.
.. _eeglab_wiki: https://sccn.ucsd.edu/wiki/Chapter_09:_Decomposing_Data_Using_ICA#Issue:_ICA_returns_near-identical_components_with_opposite_polarities
References
----------
.. [1] Hyvärinen, A., 1999. Fast and robust fixed-point algorithms for
independent component analysis. IEEE transactions on Neural
Networks, 10(3), pp.626-634.
.. [2] Bell, A.J., Sejnowski, T.J., 1995. An information-maximization
approach to blind separation and blind deconvolution. Neural
computation, 7(6), pp.1129-1159.
.. [3] Lee, T.W., Girolami, M., Sejnowski, T.J., 1999. Independent
component analysis using an extended infomax algorithm for mixed
subgaussian and supergaussian sources. Neural computation, 11(2),
pp.417-441.
.. [4] Ablin P, Cardoso J, Gramfort A, 2018. Faster Independent Component
Analysis by Preconditioning With Hessian Approximations.
IEEE Transactions on Signal Processing 66:4040–4049
.. [5] Artoni, F., Delorme, A., und Makeig, S, 2018. Applying Dimension
Reduction to EEG Data by Principal Component Analysis Reduces the
Quality of Its Subsequent Independent Component Decomposition.
NeuroImage 175, pp.176–187.
.. [6] Montoya-Martínez, J., Cardoso, J.-F., Gramfort, A, 2017. Caveats
with stochastic gradient and maximum likelihood based ICA for EEG.
LVA-ICA International Conference, Feb 2017, Grenoble, France.
`〈hal-01451432〉 <hal-01451432_>`_
.. _hal-01451432: https://hal.archives-ouvertes.fr/hal-01451432/document
""" # noqa: E501
@verbose
def __init__(self, n_components=None, max_pca_components=None,
n_pca_components=None, noise_cov=None, random_state=None,
method='fastica', fit_params=None, max_iter=200,
allow_ref_meg=False, verbose=None): # noqa: D102
_validate_type(method, str, 'method')
if method != 'imported_eeglab': # internal use only
_check_option('method', method, ['fastica', 'infomax', 'picard'])
if method == 'fastica' and not check_version('sklearn'):
raise RuntimeError(
'The scikit-learn package is required for FastICA.')
self.noise_cov = noise_cov
if (n_components is not None and
max_pca_components is not None and
n_components > max_pca_components):
raise ValueError('n_components must be smaller than '
'max_pca_components')
if isinstance(n_components, float) \
and not 0 < n_components <= 1:
raise ValueError('Selecting ICA components by explained variance '
'needs values between 0.0 and 1.0 ')
self.current_fit = 'unfitted'
self.verbose = verbose
self.n_components = n_components
self.max_pca_components = max_pca_components
self.n_pca_components = n_pca_components
self.ch_names = None
self.random_state = random_state
if fit_params is None:
fit_params = {}
fit_params = deepcopy(fit_params) # avoid side effects
if method == 'fastica':
update = {'algorithm': 'parallel', 'fun': 'logcosh',
'fun_args': None}
fit_params.update({k: v for k, v in update.items() if k
not in fit_params})
elif method == 'infomax':
# extended=True is default in underlying function, but we want
# default False here unless user specified True:
fit_params.setdefault('extended', False)
fit_params.setdefault('max_iter', max_iter)
self.max_iter = max_iter
self.fit_params = fit_params
self.exclude = []
self.info = None
self.method = method
self.labels_ = dict()
self.allow_ref_meg = allow_ref_meg
def __repr__(self):
"""ICA fit information."""
if self.current_fit == 'unfitted':
s = 'no'
elif self.current_fit == 'raw':
s = 'raw data'
else:
s = 'epochs'
s += ' decomposition, '
s += 'fit (%s): %s samples, ' % (self.method,
str(getattr(self, 'n_samples_', '')))
s += ('%s components' % str(self.n_components_) if
hasattr(self, 'n_components_') else
'no dimension reduction')
if self.info is not None:
ch_fit = ['"%s"' % c for c in _DATA_CH_TYPES_SPLIT if c in self]
s += ', channels used: {}'.format('; '.join(ch_fit))
if self.exclude:
s += ', %i sources marked for exclusion' % len(self.exclude)
return '<ICA | %s>' % s
@verbose
def fit(self, inst, picks=None, start=None, stop=None, decim=None,
reject=None, flat=None, tstep=2.0, reject_by_annotation=True,
verbose=None):
"""Run the ICA decomposition on raw data.
Caveat! If supplying a noise covariance keep track of the projections
available in the cov, the raw or the epochs object. For example,
if you are interested in EOG or ECG artifacts, EOG and ECG projections
should be temporally removed before fitting the ICA.
Parameters
----------
inst : instance of Raw, Epochs or Evoked
Raw measurements to be decomposed.
%(picks_good_data_noref)s
This selection remains throughout the initialized ICA solution.
start : int | float | None
First sample to include. If float, data will be interpreted as
time in seconds. If None, data will be used from the first sample.
stop : int | float | None
Last sample to not include. If float, data will be interpreted as
time in seconds. If None, data will be used to the last sample.
decim : int | None
Increment for selecting each nth time slice. If None, all samples
within ``start`` and ``stop`` are used.
reject : dict | None
Rejection parameters based on peak-to-peak amplitude.
Valid keys are 'grad', 'mag', 'eeg', 'seeg', 'ecog', 'eog', 'ecg',
'hbo', 'hbr'.
If reject is None then no rejection is done. Example::
reject = dict(grad=4000e-13, # T / m (gradiometers)
mag=4e-12, # T (magnetometers)
eeg=40e-6, # V (EEG channels)
eog=250e-6 # V (EOG channels)
)
It only applies if `inst` is of type Raw.
flat : dict | None
Rejection parameters based on flatness of signal.
Valid keys are 'grad', 'mag', 'eeg', 'seeg', 'ecog', 'eog', 'ecg',
'hbo', 'hbr'.
Values are floats that set the minimum acceptable peak-to-peak
amplitude. If flat is None then no rejection is done.
It only applies if `inst` is of type Raw.
tstep : float
Length of data chunks for artifact rejection in seconds.
It only applies if `inst` is of type Raw.
reject_by_annotation : bool
Whether to omit bad segments from the data before fitting. If True,
annotated segments with a description that starts with 'bad' are
omitted. Has no effect if ``inst`` is an Epochs or Evoked object.
Defaults to True.
.. versionadded:: 0.14.0
%(verbose_meth)s
Returns
-------
self : instance of ICA
Returns the modified instance.
"""
_validate_type(inst, (BaseRaw, BaseEpochs), 'inst', 'Raw or Epochs')
picks = _picks_to_idx(inst.info, picks, allow_empty=False,
with_ref_meg=self.allow_ref_meg)
_check_for_unsupported_ica_channels(
picks, inst.info, allow_ref_meg=self.allow_ref_meg)
t_start = time()
if isinstance(inst, BaseRaw):
self._fit_raw(inst, picks, start, stop, decim, reject, flat,
tstep, reject_by_annotation, verbose)
elif isinstance(inst, BaseEpochs):
self._fit_epochs(inst, picks, decim, verbose)
# sort ICA components by explained variance
var = _ica_explained_variance(self, inst)
var_ord = var.argsort()[::-1]
_sort_components(self, var_ord, copy=False)
t_stop = time()
logger.info("Fitting ICA took {:.1f}s.".format(t_stop - t_start))
return self
def _reset(self):
"""Aux method."""
del self.pre_whitener_
del self.unmixing_matrix_
del self.mixing_matrix_
del self.n_components_
del self.n_samples_
del self.pca_components_
del self.pca_explained_variance_
del self.pca_mean_
del self.n_iter_
if hasattr(self, 'drop_inds_'):
del self.drop_inds_
if hasattr(self, 'reject_'):
del self.reject_
def _fit_raw(self, raw, picks, start, stop, decim, reject, flat, tstep,
reject_by_annotation, verbose):
"""Aux method."""
if self.current_fit != 'unfitted':
self._reset()
logger.info('Fitting ICA to data using %i channels '
'(please be patient, this may take a while)' % len(picks))
if self.max_pca_components is None:
self.max_pca_components = len(picks)
logger.info('Inferring max_pca_components from picks')
self.info = pick_info(raw.info, picks)
if self.info['comps']:
self.info['comps'] = []
self.ch_names = self.info['ch_names']
start, stop = _check_start_stop(raw, start, stop)
reject_by_annotation = 'omit' if reject_by_annotation else None
# this will be a copy
data = raw.get_data(picks, start, stop, reject_by_annotation)
# this will be a view
if decim is not None:
data = data[:, ::decim]
# this will make a copy
if (reject is not None) or (flat is not None):
self.reject_ = reject
data, self.drop_inds_ = _reject_data_segments(data, reject, flat,
decim, self.info,
tstep)
self.n_samples_ = data.shape[1]
# this may operate inplace or make a copy
data, self.pre_whitener_ = self._pre_whiten(data, raw.info, picks)
self._fit(data, self.max_pca_components, 'raw')
return self
def _fit_epochs(self, epochs, picks, decim, verbose):
"""Aux method."""
if self.current_fit != 'unfitted':
self._reset()
if epochs.events.size == 0:
raise RuntimeError('Tried to fit ICA with epochs, but none were '
'found: epochs.events is "{}".'
.format(epochs.events))
logger.info('Fitting ICA to data using %i channels '
'(please be patient, this may take a while)' % len(picks))
# filter out all the channels the raw wouldn't have initialized
self.info = pick_info(epochs.info, picks)
if self.info['comps']:
self.info['comps'] = []
self.ch_names = self.info['ch_names']
if self.max_pca_components is None:
self.max_pca_components = len(picks)
logger.info('Inferring max_pca_components from picks')
# this should be a copy (picks a list of int)
data = epochs.get_data()[:, picks]
# this will be a view
if decim is not None:
data = data[:, :, ::decim]
self.n_samples_ = np.prod(data[:, 0, :].shape)
# This will make at least one copy (one from hstack, maybe one
# more from _pre_whiten)
data, self.pre_whitener_ = \
self._pre_whiten(np.hstack(data), epochs.info, picks)
self._fit(data, self.max_pca_components, 'epochs')
return self
def _pre_whiten(self, data, info, picks):
"""Aux function."""
has_pre_whitener = hasattr(self, 'pre_whitener_')
if not has_pre_whitener and self.noise_cov is None:
# use standardization as whitener
# Scale (z-score) the data by channel type
info = pick_info(info, picks)
pre_whitener = np.empty([len(data), 1])
for ch_type in _DATA_CH_TYPES_SPLIT + ('eog', "ref_meg"):
if _contains_ch_type(info, ch_type):
if ch_type == 'seeg':
this_picks = pick_types(info, meg=False, seeg=True)
elif ch_type == 'ecog':
this_picks = pick_types(info, meg=False, ecog=True)
elif ch_type == 'eeg':
this_picks = pick_types(info, meg=False, eeg=True)
elif ch_type in ('mag', 'grad'):
this_picks = pick_types(info, meg=ch_type)
elif ch_type == 'eog':
this_picks = pick_types(info, meg=False, eog=True)
elif ch_type in ('hbo', 'hbr'):
this_picks = pick_types(info, meg=False, fnirs=ch_type)
elif ch_type == 'ref_meg':
this_picks = pick_types(info, meg=False, ref_meg=True)
else:
raise RuntimeError('Should not be reached.'
'Unsupported channel {}'
.format(ch_type))
pre_whitener[this_picks] = np.std(data[this_picks])
data /= pre_whitener
elif not has_pre_whitener and self.noise_cov is not None:
pre_whitener, _ = compute_whitener(self.noise_cov, info, picks)
assert data.shape[0] == pre_whitener.shape[1]
data = np.dot(pre_whitener, data)
elif has_pre_whitener and self.noise_cov is None:
data /= self.pre_whitener_
pre_whitener = self.pre_whitener_
else:
data = np.dot(self.pre_whitener_, data)
pre_whitener = self.pre_whitener_
return data, pre_whitener
def _fit(self, data, max_pca_components, fit_type):
"""Aux function."""
random_state = check_random_state(self.random_state)
pca = _PCA(n_components=max_pca_components, whiten=True)
n_channels, n_samples = data.shape
data = pca.fit_transform(data.T)
assert data.shape == (n_samples, max_pca_components or n_channels)
if isinstance(self.n_components, float):
self.n_components_ = np.sum(
pca.explained_variance_ratio_.cumsum() <= self.n_components)
if self.n_components_ < 1:
raise RuntimeError('One PCA component captures most of the '
'explained variance, your threshold resu'
'lts in 0 components. You should select '
'a higher value.')
msg = 'Selecting by explained variance'
else:
if self.n_components is not None: # normal n case
self.n_components_ = _ensure_int(self.n_components)
msg = 'Selecting by number'
else: # None case
self.n_components_ = len(pca.components_)
msg = 'Selecting all PCA components'
logger.info('%s: %s components' % (msg, self.n_components_))
# the things to store for PCA
self.pca_mean_ = pca.mean_
self.pca_components_ = pca.components_
self.pca_explained_variance_ = pca.explained_variance_
del pca
# update number of components
self._update_ica_names()
if self.n_pca_components is not None:
if self.n_pca_components > len(self.pca_components_):
self.n_pca_components = len(self.pca_components_)
# take care of ICA
sel = slice(0, self.n_components_)
if self.method == 'fastica':
from sklearn.decomposition import FastICA
ica = FastICA(whiten=False, random_state=random_state,
**self.fit_params)
ica.fit(data[:, sel])
self.unmixing_matrix_ = ica.components_
self.n_iter_ = ica.n_iter_
elif self.method in ('infomax', 'extended-infomax'):
unmixing_matrix, n_iter = infomax(data[:, sel],
random_state=random_state,
return_n_iter=True,
**self.fit_params)
self.unmixing_matrix_ = unmixing_matrix
self.n_iter_ = n_iter
del unmixing_matrix, n_iter
elif self.method == 'picard':
from picard import picard
_, W, _, n_iter = picard(data[:, sel].T, whiten=False,
return_n_iter=True,
random_state=random_state,
**self.fit_params)
self.unmixing_matrix_ = W
self.n_iter_ = n_iter + 1 # picard() starts counting at 0
del _, n_iter
assert self.unmixing_matrix_.shape == (self.n_components_,) * 2
self.unmixing_matrix_ /= np.sqrt(
self.pca_explained_variance_[sel])[None, :] # whitening
self._update_mixing_matrix()
self.current_fit = fit_type
def _update_mixing_matrix(self):
self.mixing_matrix_ = linalg.pinv(self.unmixing_matrix_)
def _update_ica_names(self):
"""Update ICA names when n_components_ is set."""
self._ica_names = ['ICA%03d' % ii for ii in range(self.n_components_)]
def _transform(self, data):
"""Compute sources from data (operates inplace)."""
if self.pca_mean_ is not None:
data -= self.pca_mean_[:, None]
# Apply first PCA
pca_data = np.dot(self.pca_components_[:self.n_components_], data)
# Apply unmixing to low dimension PCA
sources = np.dot(self.unmixing_matrix_, pca_data)
return sources
def _transform_raw(self, raw, start, stop, reject_by_annotation=False):
"""Transform raw data."""
if not hasattr(self, 'mixing_matrix_'):
raise RuntimeError('No fit available. Please fit ICA.')
start, stop = _check_start_stop(raw, start, stop)
picks = pick_types(raw.info, include=self.ch_names, exclude='bads',
meg=False, ref_meg=False)
if len(picks) != len(self.ch_names):
raise RuntimeError('Raw doesn\'t match fitted data: %i channels '
'fitted but %i channels supplied. \nPlease '
'provide Raw compatible with '
'ica.ch_names' % (len(self.ch_names),
len(picks)))
if reject_by_annotation:
data = raw.get_data(picks, start, stop, 'omit')
else:
data = raw[picks, start:stop][0]
data, _ = self._pre_whiten(data, raw.info, picks)
return self._transform(data)
def _transform_epochs(self, epochs, concatenate):
"""Aux method."""
if not hasattr(self, 'mixing_matrix_'):
raise RuntimeError('No fit available. Please fit ICA.')
picks = pick_types(epochs.info, include=self.ch_names, exclude='bads',
meg=False, ref_meg=False)
# special case where epochs come picked but fit was 'unpicked'.
if len(picks) != len(self.ch_names):
raise RuntimeError('Epochs don\'t match fitted data: %i channels '
'fitted but %i channels supplied. \nPlease '
'provide Epochs compatible with '
'ica.ch_names' % (len(self.ch_names),
len(picks)))
data = np.hstack(epochs.get_data()[:, picks])
data, _ = self._pre_whiten(data, epochs.info, picks)
sources = self._transform(data)
if not concatenate:
# Put the data back in 3D
sources = np.array(np.split(sources, len(epochs.events), 1))
return sources
def _transform_evoked(self, evoked):
"""Aux method."""
if not hasattr(self, 'mixing_matrix_'):
raise RuntimeError('No fit available. Please fit ICA.')
picks = pick_types(evoked.info, include=self.ch_names, exclude='bads',
meg=False, ref_meg=False)
if len(picks) != len(self.ch_names):
raise RuntimeError('Evoked doesn\'t match fitted data: %i channels'
' fitted but %i channels supplied. \nPlease '
'provide Evoked compatible with '
'ica.ch_names' % (len(self.ch_names),
len(picks)))
data, _ = self._pre_whiten(evoked.data[picks], evoked.info, picks)
sources = self._transform(data)
return sources
def get_components(self):
"""Get ICA topomap for components as numpy arrays.
Returns
-------
components : array, shape (n_channels, n_components)
The ICA components (maps).
"""
return np.dot(self.mixing_matrix_[:, :self.n_components_].T,
self.pca_components_[:self.n_components_]).T
def get_sources(self, inst, add_channels=None, start=None, stop=None):
"""Estimate sources given the unmixing matrix.
This method will return the sources in the container format passed.
Typical usecases:
1. pass Raw object to use `raw.plot` for ICA sources
2. pass Epochs object to compute trial-based statistics in ICA space
3. pass Evoked object to investigate time-locking in ICA space
Parameters
----------
inst : instance of Raw, Epochs or Evoked
Object to compute sources from and to represent sources in.
add_channels : None | list of str
Additional channels to be added. Useful to e.g. compare sources
with some reference. Defaults to None.
start : int | float | None
First sample to include. If float, data will be interpreted as
time in seconds. If None, the entire data will be used.
stop : int | float | None
Last sample to not include. If float, data will be interpreted as
time in seconds. If None, the entire data will be used.
Returns
-------
sources : instance of Raw, Epochs or Evoked
The ICA sources time series.
"""
if isinstance(inst, BaseRaw):
_check_compensation_grade(self.info, inst.info, 'ICA', 'Raw',
ch_names=self.ch_names)
sources = self._sources_as_raw(inst, add_channels, start, stop)
elif isinstance(inst, BaseEpochs):
_check_compensation_grade(self.info, inst.info, 'ICA', 'Epochs',
ch_names=self.ch_names)
sources = self._sources_as_epochs(inst, add_channels, False)
elif isinstance(inst, Evoked):
_check_compensation_grade(self.info, inst.info, 'ICA', 'Evoked',
ch_names=self.ch_names)
sources = self._sources_as_evoked(inst, add_channels)
else:
raise ValueError('Data input must be of Raw, Epochs or Evoked '
'type')
return sources
def _sources_as_raw(self, raw, add_channels, start, stop):
"""Aux method."""
# merge copied instance and picked data with sources
sources = self._transform_raw(raw, start=start, stop=stop)
if raw.preload: # get data and temporarily delete
data = raw._data
del raw._data
out = raw.copy() # copy and reappend
if raw.preload:
raw._data = data
# populate copied raw.
start, stop = _check_start_stop(raw, start, stop)
if add_channels is not None:
raw_picked = raw.copy().pick_channels(add_channels)
data_, times_ = raw_picked[:, start:stop]
data_ = np.r_[sources, data_]
else:
data_ = sources
_, times_ = raw[0, start:stop]
out._data = data_
out._times = times_
out._filenames = [None]
out.preload = True
# update first and last samples
out._first_samps = np.array([raw.first_samp +
(start if start else 0)])
out._last_samps = np.array([out.first_samp + stop
if stop else raw.last_samp])
out._projector = None
self._export_info(out.info, raw, add_channels)
out._update_times()
return out
def _sources_as_epochs(self, epochs, add_channels, concatenate):
"""Aux method."""
out = epochs.copy()
sources = self._transform_epochs(epochs, concatenate)
if add_channels is not None:
picks = [epochs.ch_names.index(k) for k in add_channels]
else:
picks = []
out._data = np.concatenate([sources, epochs.get_data()[:, picks]],
axis=1) if len(picks) > 0 else sources
self._export_info(out.info, epochs, add_channels)
out.preload = True
out._raw = None
out._projector = None
return out
def _sources_as_evoked(self, evoked, add_channels):
"""Aux method."""
if add_channels is not None:
picks = [evoked.ch_names.index(k) for k in add_channels]
else:
picks = []
sources = self._transform_evoked(evoked)
if len(picks) > 1:
data = np.r_[sources, evoked.data[picks]]
else:
data = sources
out = evoked.copy()
out.data = data
self._export_info(out.info, evoked, add_channels)
return out
def _export_info(self, info, container, add_channels):
"""Aux method."""
# set channel names and info
ch_names = []
ch_info = info['chs'] = []
for ii, name in enumerate(self._ica_names):
ch_names.append(name)
ch_info.append(dict(
ch_name=name, cal=1, logno=ii + 1,
coil_type=FIFF.FIFFV_COIL_NONE, kind=FIFF.FIFFV_MISC_CH,
coord_Frame=FIFF.FIFFV_COORD_UNKNOWN, unit=FIFF.FIFF_UNIT_NONE,
loc=np.array([0., 0., 0., 1.] * 3, dtype='f4'),
range=1.0, scanno=ii + 1, unit_mul=0))
if add_channels is not None:
# re-append additionally picked ch_names
ch_names += add_channels
# re-append additionally picked ch_info
ch_info += [k for k in container.info['chs'] if k['ch_name'] in
add_channels]
info['bads'] = [ch_names[k] for k in self.exclude]
info['projs'] = [] # make sure projections are removed.
info._update_redundant()
info._check_consistency()
@verbose
def score_sources(self, inst, target=None, score_func='pearsonr',
start=None, stop=None, l_freq=None, h_freq=None,
reject_by_annotation=True, verbose=None):
"""Assign score to components based on statistic or metric.
Parameters
----------
inst : instance of Raw, Epochs or Evoked
The object to reconstruct the sources from.
target : array-like | str | None
Signal to which the sources shall be compared. It has to be of
the same shape as the sources. If str, a routine will try to find
a matching channel name. If None, a score
function expecting only one input-array argument must be used,
for instance, scipy.stats.skew (default).
score_func : callable | str
Callable taking as arguments either two input arrays
(e.g. Pearson correlation) or one input
array (e. g. skewness) and returns a float. For convenience the
most common score_funcs are available via string labels:
Currently, all distance metrics from scipy.spatial and All
functions from scipy.stats taking compatible input arguments are
supported. These function have been modified to support iteration
over the rows of a 2D array.
start : int | float | None
First sample to include. If float, data will be interpreted as
time in seconds. If None, data will be used from the first sample.
stop : int | float | None
Last sample to not include. If float, data will be interpreted as
time in seconds. If None, data will be used to the last sample.
l_freq : float
Low pass frequency.
h_freq : float
High pass frequency.
reject_by_annotation : bool
If True, data annotated as bad will be omitted. Defaults to True.
.. versionadded:: 0.14.0
%(verbose_meth)s
Returns
-------
scores : ndarray
Scores for each source as returned from score_func.
"""
if isinstance(inst, BaseRaw):
_check_compensation_grade(self.info, inst.info, 'ICA', 'Raw',
ch_names=self.ch_names)
sources = self._transform_raw(inst, start, stop,