[MRG] Rank refactor part2 (#5870)

* start moving things to rank.py

* more

* more

* more

* more

* fix

* FIX: Minor fixes

Author: agramfort

mne/cov.py

@@ -14,17 +14,18 @@
 import numpy as np
 from scipy import linalg, sparse
 
+from .io.meas_info import _simplify_info
 from .io.write import start_file, end_file
 from .io.proj import (make_projector, _proj_equal, activate_proj,
-                      _needs_eeg_average_ref_proj, _check_projs,
+                      _check_projs, _needs_eeg_average_ref_proj,
                       _has_eeg_average_ref_proj)
 from .io import fiff_open
 from .io.pick import (pick_types, pick_channels_cov, pick_channels, pick_info,
                       _picks_by_type, _pick_data_channels,
                       _DATA_CH_TYPES_SPLIT)
 
 from .io.constants import FIFF
-from .io.meas_info import read_bad_channels, _simplify_info, create_info
+from .io.meas_info import read_bad_channels, create_info
 from .io.proj import _read_proj, _write_proj
 from .io.tag import find_tag
 from .io.tree import dir_tree_find
@@ -33,10 +34,13 @@
 from .defaults import _handle_default
 from .epochs import Epochs
 from .event import make_fixed_length_events
-from .utils import (check_fname, logger, verbose, estimate_rank,
-                    _compute_row_norms, check_version, _time_mask, warn,
-                    copy_function_doc_to_method_doc, _pl)
+from .utils import (check_fname, logger, verbose,
+                    check_version, _time_mask, warn,
+                    copy_function_doc_to_method_doc, _pl,
+                    _undo_scaling_cov, _undo_scaling_array,
+                    _apply_scaling_array)
 from . import viz
+from .rank import _estimate_rank_meeg_cov
 
 from .fixes import BaseEstimator, EmpiricalCovariance, _logdet
 
@@ -629,7 +633,7 @@ def compute_covariance(epochs, keep_sample_mean=True, tmin=None, tmax=None,
 
              ['shrunk', 'diagonal_fixed', 'empirical', 'factor_analysis']
 
-        ``'factor_analysis'`` is removed when `rank` is not 'full'.
+        ``'factor_analysis'`` is removed when ``rank`` is not 'full'.
         The ``'auto'`` mode is not recommended if there are many
         segments of data, since computation can take a long time.
 
@@ -1990,175 +1994,3 @@ def _write_cov(fid, cov):
 
     #   Done!
     end_block(fid, FIFF.FIFFB_MNE_COV)
-
-
-def _apply_scaling_array(data, picks_list, scalings):
-    """Scale data type-dependently for estimation."""
-    scalings = _check_scaling_inputs(data, picks_list, scalings)
-    if isinstance(scalings, dict):
-        picks_dict = dict(picks_list)
-        scalings = [(picks_dict[k], v) for k, v in scalings.items()
-                    if k in picks_dict]
-        for idx, scaling in scalings:
-            data[idx, :] *= scaling  # F - order
-    else:
-        data *= scalings[:, np.newaxis]  # F - order
-
-
-def _invert_scalings(scalings):
-    if isinstance(scalings, dict):
-        scalings = dict((k, 1. / v) for k, v in scalings.items())
-    elif isinstance(scalings, np.ndarray):
-        scalings = 1. / scalings
-    return scalings
-
-
-def _undo_scaling_array(data, picks_list, scalings):
-    scalings = _invert_scalings(_check_scaling_inputs(data, picks_list,
-                                                      scalings))
-    return _apply_scaling_array(data, picks_list, scalings)
-
-
-def _apply_scaling_cov(data, picks_list, scalings):
-    """Scale resulting data after estimation."""
-    scalings = _check_scaling_inputs(data, picks_list, scalings)
-    scales = None
-    if isinstance(scalings, dict):
-        n_channels = len(data)
-        covinds = list(zip(*picks_list))[1]
-        assert len(data) == sum(len(k) for k in covinds)
-        assert list(sorted(np.concatenate(covinds))) == list(range(len(data)))
-        scales = np.zeros(n_channels)
-        for ch_t, idx in picks_list:
-            scales[idx] = scalings[ch_t]
-    elif isinstance(scalings, np.ndarray):
-        if len(scalings) != len(data):
-            raise ValueError('Scaling factors and data are of incompatible '
-                             'shape')
-        scales = scalings
-    elif scalings is None:
-        pass
-    else:
-        raise RuntimeError('Arff...')
-    if scales is not None:
-        assert np.sum(scales == 0.) == 0
-        data *= (scales[None, :] * scales[:, None])
-
-
-def _undo_scaling_cov(data, picks_list, scalings):
-    scalings = _invert_scalings(_check_scaling_inputs(data, picks_list,
-                                                      scalings))
-    return _apply_scaling_cov(data, picks_list, scalings)
-
-
-def _check_scaling_inputs(data, picks_list, scalings):
-    """Aux function."""
-    rescale_dict_ = dict(mag=1e15, grad=1e13, eeg=1e6)
-
-    scalings_ = None
-    if isinstance(scalings, str) and scalings == 'norm':
-        scalings_ = 1. / _compute_row_norms(data)
-    elif isinstance(scalings, dict):
-        rescale_dict_.update(scalings)
-        scalings_ = rescale_dict_
-    elif isinstance(scalings, np.ndarray):
-        scalings_ = scalings
-    elif scalings is None:
-        pass
-    else:
-        raise NotImplementedError("No way! That's not a rescaling "
-                                  'option: %s' % scalings)
-    return scalings_
-
-
-def _estimate_rank_meeg_signals(data, info, scalings, tol='auto',
-                                return_singular=False):
-    """Estimate rank for M/EEG data.
-
-    Parameters
-    ----------
-    data : np.ndarray of float, shape(n_channels, n_samples)
-        The M/EEG signals.
-    info : Info
-        The measurement info.
-    scalings : dict | 'norm' | np.ndarray | None
-        The rescaling method to be applied. If dict, it will override the
-        following default dict:
-
-            dict(mag=1e15, grad=1e13, eeg=1e6)
-
-        If 'norm' data will be scaled by channel-wise norms. If array,
-        pre-specified norms will be used. If None, no scaling will be applied.
-    tol : float | str
-        Tolerance. See ``estimate_rank``.
-    return_singular : bool
-        If True, also return the singular values that were used
-        to determine the rank.
-
-    Returns
-    -------
-    rank : int
-        Estimated rank of the data.
-    s : array
-        If return_singular is True, the singular values that were
-        thresholded to determine the rank are also returned.
-    """
-    picks_list = _picks_by_type(info)
-    _apply_scaling_array(data, picks_list, scalings)
-    if data.shape[1] < data.shape[0]:
-        ValueError("You've got fewer samples than channels, your "
-                   "rank estimate might be inaccurate.")
-    out = estimate_rank(data, tol=tol, norm=False,
-                        return_singular=return_singular)
-    rank = out[0] if isinstance(out, tuple) else out
-    ch_type = ' + '.join(list(zip(*picks_list))[0])
-    logger.info('estimated rank (%s): %d' % (ch_type, rank))
-    _undo_scaling_array(data, picks_list, scalings)
-    return out
-
-
-def _estimate_rank_meeg_cov(data, info, scalings, tol='auto',
-                            return_singular=False):
-    """Estimate rank of M/EEG covariance data, given the covariance.
-
-    Parameters
-    ----------
-    data : np.ndarray of float, shape (n_channels, n_channels)
-        The M/EEG covariance.
-    info : Info
-        The measurement info.
-    scalings : dict | 'norm' | np.ndarray | None
-        The rescaling method to be applied. If dict, it will override the
-        following default dict:
-
-            dict(mag=1e12, grad=1e11, eeg=1e5)
-
-        If 'norm' data will be scaled by channel-wise norms. If array,
-        pre-specified norms will be used. If None, no scaling will be applied.
-    tol : float | str
-        Tolerance. See ``estimate_rank``.
-    return_singular : bool
-        If True, also return the singular values that were used
-        to determine the rank.
-
-    Returns
-    -------
-    rank : int
-        Estimated rank of the data.
-    s : array
-        If return_singular is True, the singular values that were
-        thresholded to determine the rank are also returned.
-    """
-    picks_list = _picks_by_type(info)
-    scalings = _handle_default('scalings_cov_rank', scalings)
-    _apply_scaling_cov(data, picks_list, scalings)
-    if data.shape[1] < data.shape[0]:
-        ValueError("You've got fewer samples than channels, your "
-                   "rank estimate might be inaccurate.")
-    out = estimate_rank(data, tol=tol, norm=False,
-                        return_singular=return_singular)
-    rank = out[0] if isinstance(out, tuple) else out
-    ch_type = ' + '.join(list(zip(*picks_list))[0])
-    logger.info('estimated rank (%s): %d' % (ch_type, rank))
-    _undo_scaling_cov(data, picks_list, scalings)
-    return out

mne/io/base.py

@@ -1959,7 +1959,7 @@ def estimate_rank(self, tstart=0.0, tstop=30.0, tol=1e-4,
 
         Bad channels will be excluded from calculations.
         """
-        from ..cov import _estimate_rank_meeg_signals
+        from ..rank import _estimate_rank_meeg_signals
 
         start = max(0, self.time_as_index(tstart)[0])
         if tstop is None:

mne/io/fiff/tests/test_raw_fiff.py

@@ -6,7 +6,6 @@
 
 from copy import deepcopy
 from functools import partial
-import itertools as itt
 import os.path as op
 import pickle
 import sys
@@ -26,8 +25,6 @@
                  pick_info)
 from mne.utils import (_TempDir, requires_pandas, object_diff,
                        requires_mne, run_subprocess, run_tests_if_main)
-from mne.io.proc_history import _get_rank_sss
-from mne.io.pick import _picks_by_type
 from mne.annotations import Annotations
 
 testing_path = testing.data_path(download=False)
@@ -202,37 +199,6 @@ def test_copy_append():
     assert_equal(data.shape[1], 2 * raw._data.shape[1])
 
 
-@pytest.mark.slowtest
-@testing.requires_testing_data
-def test_rank_estimation():
-    """Test raw rank estimation."""
-    iter_tests = itt.product(
-        [fif_fname, hp_fif_fname],  # sss
-        ['norm', dict(mag=1e11, grad=1e9, eeg=1e5)]
-    )
-    for fname, scalings in iter_tests:
-        raw = read_raw_fif(fname).crop(0, 4.).load_data()
-        (_, picks_meg), (_, picks_eeg) = _picks_by_type(raw.info,
-                                                        meg_combined=True)
-        n_meg = len(picks_meg)
-        n_eeg = len(picks_eeg)
-
-        if len(raw.info['proc_history']) == 0:
-            expected_rank = n_meg + n_eeg
-        else:
-            expected_rank = _get_rank_sss(raw.info) + n_eeg
-        assert_array_equal(raw.estimate_rank(scalings=scalings), expected_rank)
-        assert_array_equal(raw.estimate_rank(picks=picks_eeg,
-                                             scalings=scalings), n_eeg)
-        if 'sss' in fname:
-            raw.add_proj(compute_proj_raw(raw))
-        raw.apply_proj()
-        n_proj = len(raw.info['projs'])
-        assert_array_equal(raw.estimate_rank(tstart=0, tstop=3.,
-                                             scalings=scalings),
-                           expected_rank - (0 if 'sss' in fname else n_proj))
-
-
 @testing.requires_testing_data
 def test_output_formats():
     """Test saving and loading raw data using multiple formats."""

mne/io/proc_history.py

@@ -15,7 +15,7 @@
                     write_float_sparse, write_id)
 from .tag import find_tag
 from .constants import FIFF
-from ..utils import warn, logger
+from ..utils import warn
 
 _proc_keys = ['parent_file_id', 'block_id', 'parent_block_id',
               'date', 'experimenter', 'creator']
@@ -295,52 +295,3 @@ def _write_maxfilter_record(fid, record):
             if key in sss_cal:
                 writer(fid, id_, sss_cal[key])
         end_block(fid, FIFF.FIFFB_SSS_CAL)
-
-
-def _get_sss_rank(sss):
-    """Get SSS rank."""
-    inside = sss['sss_info']['in_order']
-    nfree = (inside + 1) ** 2 - 1
-    nfree -= (len(sss['sss_info']['components'][:nfree]) -
-              sss['sss_info']['components'][:nfree].sum())
-    return nfree
-
-
-def _get_rank_sss(inst):
-    """Look up rank from SSS data.
-
-    .. note::
-        Throws an error if SSS has not been applied.
-
-    Parameters
-    ----------
-    inst : instance of Raw, Epochs or Evoked, or Info
-        Any MNE object with an .info attribute
-
-    Returns
-    -------
-    rank : int
-        The numerical rank as predicted by the number of SSS
-        components.
-    """
-    from .meas_info import Info
-    info = inst if isinstance(inst, Info) else inst.info
-    del inst
-
-    max_infos = list()
-    for proc_info in info.get('proc_history', list()):
-        max_info = proc_info.get('max_info')
-        if max_info is not None:
-            if len(max_info) > 0:
-                max_infos.append(max_info)
-            elif len(max_info) > 1:
-                logger.info('found multiple SSS records. Using the first.')
-            elif len(max_info) == 0:
-                raise ValueError(
-                    'Did not find any SSS record. You should use data-based '
-                    'rank estimate instead')
-    if len(max_infos) > 0:
-        max_info = max_infos[0]
-    else:
-        raise ValueError('There is no `max_info` here. Sorry.')
-    return _get_sss_rank(max_info)

mne/io/tests/test_proc_history.py

@@ -7,9 +7,8 @@
 import numpy as np
 from numpy.testing import assert_array_equal
 
-from mne.io import read_info, read_raw_fif
+from mne.io import read_info
 from mne.io.constants import FIFF
-from mne.io.proc_history import _get_rank_sss
 
 base_dir = op.join(op.dirname(__file__), 'data')
 raw_fname = op.join(base_dir, 'test_chpi_raw_sss.fif')
@@ -37,12 +36,3 @@ def test_maxfilter_io():
     assert mf['sss_cal']['cal_chans'].shape == (306, 2)
     vv_coils = [v for k, v in FIFF.items() if 'FIFFV_COIL_VV' in k]
     assert all(k in vv_coils for k in set(mf['sss_cal']['cal_chans'][:, 1]))
-
-
-def test_maxfilter_get_rank():
-    """Test maxfilter rank lookup."""
-    raw = read_raw_fif(raw_fname)
-    mf = raw.info['proc_history'][0]['max_info']
-    rank1 = mf['sss_info']['nfree']
-    rank2 = _get_rank_sss(raw)
-    assert rank1 == rank2

mne/preprocessing/tests/test_maxwell.py

@@ -18,12 +18,12 @@
 from mne.cov import _estimate_rank_meeg_cov
 from mne.datasets import testing
 from mne.forward import use_coil_def
-from mne.io import (read_raw_fif, proc_history, read_info, read_raw_bti,
-                    read_raw_kit, BaseRaw)
+from mne.io import read_raw_fif, read_info, read_raw_bti, read_raw_kit, BaseRaw
 from mne.preprocessing.maxwell import (
     maxwell_filter, _get_n_moments, _sss_basis_basic, _sh_complex_to_real,
     _sh_real_to_complex, _sh_negate, _bases_complex_to_real, _trans_sss_basis,
     _bases_real_to_complex, _prep_mf_coils)
+from mne.rank import _get_sss_rank
 from mne.tests.common import assert_meg_snr
 from mne.utils import (_TempDir, run_tests_if_main, catch_logging,
                        requires_version, object_diff, buggy_mkl_svd)
@@ -439,8 +439,7 @@ def test_basic():
 
     # Check against SSS functions from proc_history
     sss_info = raw_sss.info['proc_history'][0]['max_info']
-    assert_equal(_get_n_moments(int_order),
-                 proc_history._get_sss_rank(sss_info))
+    assert_equal(_get_n_moments(int_order), _get_sss_rank(sss_info))
 
     # Degenerate cases
     pytest.raises(ValueError, maxwell_filter, raw, coord_frame='foo')