Aluode/PerceptionLabPortable

1"""Some utility functions for rank estimation."""2 3# Authors: The MNE-Python contributors.4# License: BSD-3-Clause5# Copyright the MNE-Python contributors.6 7import numpy as np8from scipy import linalg9 10from ._fiff.meas_info import Info, _simplify_info11from ._fiff.pick import _picks_by_type, _picks_to_idx, pick_channels_cov, pick_info12from ._fiff.proj import make_projector13from .defaults import _handle_default14from .utils import (15    _apply_scaling_cov,16    _check_on_missing,17    _check_rank,18    _compute_row_norms,19    _on_missing,20    _pl,21    _scaled_array,22    _undo_scaling_cov,23    _validate_type,24    fill_doc,25    logger,26    verbose,27    warn,28)29 30 31@verbose32def estimate_rank(33    data,34    tol="auto",35    return_singular=False,36    norm=True,37    tol_kind="absolute",38    verbose=None,39):40    """Estimate the rank of data.41 42    This function will normalize the rows of the data (typically43    channels or vertices) such that non-zero singular values44    should be close to one.45 46    Parameters47    ----------48    data : array49        Data to estimate the rank of (should be 2-dimensional).50    %(tol_rank)s51    return_singular : bool52        If True, also return the singular values that were used53        to determine the rank.54    norm : bool55        If True, data will be scaled by their estimated row-wise norm.56        Else data are assumed to be scaled. Defaults to True.57    %(tol_kind_rank)s58 59    Returns60    -------61    rank : int62        Estimated rank of the data.63    s : array64        If return_singular is True, the singular values that were65        thresholded to determine the rank are also returned.66    """67    if norm:68        data = data.copy()  # operate on a copy69        norms = _compute_row_norms(data)70        data /= norms[:, np.newaxis]71    s = linalg.svdvals(data)72    rank = _estimate_rank_from_s(s, tol, tol_kind)73    if return_singular is True:74        return rank, s75    else:76        return rank77 78 79def _estimate_rank_from_s(s, tol="auto", tol_kind="absolute"):80    """Estimate the rank of a matrix from its singular values.81 82    Parameters83    ----------84    s : ndarray, shape (..., ndim)85        The singular values of the matrix.86    tol : float | ``'auto'``87        Tolerance for singular values to consider non-zero in calculating the88        rank. Can be 'auto' to use the same thresholding as89        ``scipy.linalg.orth`` (assuming np.float64 datatype) adjusted90        by a factor of 2.91    tol_kind : str92        Can be ``"absolute"`` or ``"relative"``.93 94    Returns95    -------96    rank : ndarray, shape (...)97        The estimated rank.98    """99    s = np.array(s, float)100    max_s = np.amax(s, axis=-1)101    if isinstance(tol, str):102        if tol not in ("auto", "float32"):103            raise ValueError(f'tol must be "auto" or float, got {repr(tol)}')104        # XXX this should be float32 probably due to how we save and105        # load data, but it breaks test_make_inverse_operator (!)106        # The factor of 2 gets test_compute_covariance_auto_reg[None]107        # to pass without breaking minimum norm tests. :(108        # Passing 'float32' is a hack workaround for test_maxfilter_get_rank :(109        if tol == "float32":110            eps = np.finfo(np.float32).eps111        else:112            eps = np.finfo(np.float64).eps113        tol = s.shape[-1] * max_s * eps114        if s.ndim == 1:  # typical115            logger.info(116                "    Using tolerance %0.2g (%0.2g eps * %d dim * %0.2g"117                "  max singular value)",118                tol,119                eps,120                len(s),121                max_s,122            )123    elif not (isinstance(tol, np.ndarray) and tol.dtype.kind == "f"):124        tol = float(tol)125        if tol_kind == "relative":126            tol = tol * max_s127 128    rank = np.sum(s > tol, axis=-1)129    return rank130 131 132def _estimate_rank_raw(133    raw,134    picks=None,135    tol=1e-4,136    scalings="norm",137    with_ref_meg=False,138    tol_kind="absolute",139    on_few_samples="warn",140):141    """Aid the transition away from raw.estimate_rank."""142    if picks is None:143        picks = _picks_to_idx(raw.info, picks, with_ref_meg=with_ref_meg)144    # conveniency wrapper to expose the expert "tol" option + scalings options145    return _estimate_rank_meeg_signals(146        raw[picks][0],147        pick_info(raw.info, picks),148        scalings,149        tol,150        False,151        tol_kind,152        log_ch_type=None,153        on_few_samples=on_few_samples,154    )155 156 157@fill_doc158def _estimate_rank_meeg_signals(159    data,160    info,161    scalings,162    tol="auto",163    return_singular=False,164    tol_kind="absolute",165    log_ch_type=None,166    on_few_samples="warn",167):168    """Estimate rank for M/EEG data.169 170    Parameters171    ----------172    data : np.ndarray of float, shape(n_channels, n_samples)173        The M/EEG signals.174    %(info_not_none)s175    scalings : dict | ``'norm'`` | np.ndarray | None176        The rescaling method to be applied. If dict, it will override the177        following default dict:178 179            dict(mag=1e15, grad=1e13, eeg=1e6)180 181        If ``'norm'`` data will be scaled by channel-wise norms. If array,182        pre-specified norms will be used. If None, no scaling will be applied.183    tol : float | str184        Tolerance. See ``estimate_rank``.185    return_singular : bool186        If True, also return the singular values that were used187        to determine the rank.188    tol_kind : str189        Tolerance kind. See ``estimate_rank``.190    on_few_samples : str191        Can be 'warn' (default), 'ignore', or 'raise' to control behavior when192        there are fewer samples than channels, which can lead to inaccurate rank193        estimates.194 195    Returns196    -------197    rank : int198        Estimated rank of the data.199    s : array200        If return_singular is True, the singular values that were201        thresholded to determine the rank are also returned.202    """203    picks_list = _picks_by_type(info)204    assert data.ndim == 2, data.shape205    n_channels, n_samples = data.shape206    if n_samples < n_channels:207        msg = (208            f"Too few samples ({n_samples=} is less than {n_channels=}), "209            "rank estimate may be unreliable"210        )211        _on_missing(on_few_samples, msg, "on_few_samples")212    with _scaled_array(data, picks_list, scalings):213        out = estimate_rank(214            data,215            tol=tol,216            norm=False,217            return_singular=return_singular,218            tol_kind=tol_kind,219        )220    rank = out[0] if isinstance(out, tuple) else out221    if log_ch_type is None:222        ch_type = " + ".join(list(zip(*picks_list))[0])223    else:224        ch_type = log_ch_type225    logger.info("    Estimated rank (%s): %d", ch_type, rank)226    return out227 228 229@verbose230def _estimate_rank_meeg_cov(231    data,232    info,233    scalings,234    tol="auto",235    return_singular=False,236    *,237    log_ch_type=None,238    on_few_samples="warn",239    verbose=None,240):241    """Estimate rank of M/EEG covariance data, given the covariance.242 243    Parameters244    ----------245    data : np.ndarray of float, shape (n_channels, n_channels)246        The M/EEG covariance.247    %(info_not_none)s248    scalings : dict | 'norm' | np.ndarray | None249        The rescaling method to be applied. If dict, it will override the250        following default dict:251 252            dict(mag=1e12, grad=1e11, eeg=1e5)253 254        If 'norm' data will be scaled by channel-wise norms. If array,255        pre-specified norms will be used. If None, no scaling will be applied.256    tol : float | str257        Tolerance. See ``estimate_rank``.258    return_singular : bool259        If True, also return the singular values that were used260        to determine the rank.261    on_few_samples : str262        Can be 'warn' (default), 'ignore', or 'raise' to control behavior when263        there are fewer samples than channels, which can lead to inaccurate rank264        estimates.265 266    Returns267    -------268    rank : int269        Estimated rank of the data.270    s : array271        If return_singular is True, the singular values that were272        thresholded to determine the rank are also returned.273    """274    picks_list = _picks_by_type(info, exclude=[])275    scalings = _handle_default("scalings_cov_rank", scalings)276    _apply_scaling_cov(data, picks_list, scalings)277    if data.shape[1] < data.shape[0]:278        msg = (279            "You've got fewer samples than channels, your "280            "rank estimate might be inaccurate."281        )282        _on_missing(on_few_samples, msg, "on_few_samples")283    out = estimate_rank(data, tol=tol, norm=False, return_singular=return_singular)284    rank = out[0] if isinstance(out, tuple) else out285    if log_ch_type is None:286        ch_type_ = " + ".join(list(zip(*picks_list))[0])287    else:288        ch_type_ = log_ch_type289    logger.info(f"    Estimated rank ({ch_type_}): {rank}")290    _undo_scaling_cov(data, picks_list, scalings)291    return out292 293 294@verbose295def _get_rank_sss(296    inst, msg="You should use data-based rank estimate instead", verbose=None297):298    """Look up rank from SSS data.299 300    .. note::301        Throws an error if SSS has not been applied.302 303    Parameters304    ----------305    inst : instance of Raw, Epochs or Evoked, or Info306        Any MNE object with an .info attribute307 308    Returns309    -------310    rank : int311        The numerical rank as predicted by the number of SSS312        components.313    """314    # XXX this is too basic for movement compensated data315    # https://github.com/mne-tools/mne-python/issues/4676316    info = inst if isinstance(inst, Info) else inst.info317    del inst318 319    proc_info = info.get("proc_history", [])320    if len(proc_info) > 1:321        logger.info("Found multiple SSS records. Using the first.")322    if (323        len(proc_info) == 0324        or "max_info" not in proc_info[0]325        or "in_order" not in proc_info[0]["max_info"]["sss_info"]326    ):327        raise ValueError(328            f'Could not find Maxfilter information in info["proc_history"]. {msg}'329        )330    proc_info = proc_info[0]331    max_info = proc_info["max_info"]332    inside = max_info["sss_info"]["in_order"]333    nfree = (inside + 1) ** 2 - 1334    nfree -= (335        len(max_info["sss_info"]["components"][:nfree])336        - max_info["sss_info"]["components"][:nfree].sum()337    )338    return nfree339 340 341def _info_rank(info, ch_type, picks, rank):342    if ch_type in ["meg", "mag", "grad"] and rank != "full":343        try:344            return _get_rank_sss(info)345        except ValueError:346            pass347    return len(picks)348 349 350def _compute_rank_int(inst, *args, **kwargs):351    """Wrap compute_rank but yield an int."""352    # XXX eventually we should unify how channel types are handled353    # so that we don't need to do this, or we do it everywhere.354    # Using pca=True in compute_whitener might help.355    return sum(compute_rank(inst, *args, on_few_samples="ignore", **kwargs).values())356 357 358@verbose359def compute_rank(360    inst,361    rank=None,362    scalings=None,363    info=None,364    tol="auto",365    *,366    proj=True,367    tol_kind="absolute",368    on_rank_mismatch="ignore",369    on_few_samples=None,370    verbose=None,371):372    """Compute the rank of data or noise covariance.373 374    This function will normalize the rows of the data (typically375    channels or vertices) such that non-zero singular values376    should be close to one. It operates on :term:`data channels` only.377 378    Parameters379    ----------380    inst : instance of Raw, Epochs, or Covariance381        Raw measurements to compute the rank from or the covariance.382    %(rank_none)s383    scalings : dict | None (default None)384        Defaults to ``dict(mag=1e15, grad=1e13, eeg=1e6)``.385        These defaults will scale different channel types386        to comparable values.387    %(info)s Only necessary if ``inst`` is a :class:`mne.Covariance`388        object (since this does not provide ``inst.info``).389    %(tol_rank)s390    proj : bool391        If True, all projs in ``inst`` and ``info`` will be applied or392        considered when ``rank=None`` or ``rank='info'``.393    %(tol_kind_rank)s394    %(on_rank_mismatch)s395    on_few_samples : str | None396        Can be 'warn', 'ignore', or 'raise' to control behavior when397        there are fewer samples than channels, which can lead to inaccurate rank398        estimates. None (default) means "ignore" if ``inst`` is a399        :class:`mne.Covariance` or ``rank in ("info", "full")``, and "warn" otherwise.400 401        .. versionadded:: 1.11402    %(verbose)s403 404    Returns405    -------406    rank : dict407        Estimated rank of the data for each channel type.408        To get the total rank, you can use ``sum(rank.values())``.409 410    Notes411    -----412    .. versionadded:: 0.18413    """414    return _compute_rank(415        inst=inst,416        rank=rank,417        scalings=scalings,418        info=info,419        tol=tol,420        proj=proj,421        tol_kind=tol_kind,422        on_rank_mismatch=on_rank_mismatch,423        on_few_samples=on_few_samples,424    )425 426 427@verbose428def _compute_rank(429    inst,430    rank=None,431    scalings=None,432    info=None,433    *,434    tol="auto",435    proj=True,436    tol_kind="absolute",437    on_rank_mismatch="ignore",438    on_few_samples=None,439    log_ch_type=None,440    verbose=None,441):442    from .cov import Covariance443    from .epochs import BaseEpochs444    from .io import BaseRaw445 446    rank = _check_rank(rank)447    scalings = _handle_default("scalings_cov_rank", scalings)448    _check_on_missing(on_rank_mismatch, "on_rank_mismatch")449 450    if isinstance(inst, Covariance):451        inst_type = "covariance"452        if info is None:453            raise ValueError("info cannot be None if inst is a Covariance.")454        # Reset bads as it's already taken into account in inst['names']455        info = info.copy()456        info["bads"] = []457        inst = pick_channels_cov(458            inst,459            set(inst["names"]) & set(info["ch_names"]),460            exclude=info["bads"] + inst["bads"],461            ordered=False,462        )463        if info["ch_names"] != inst["names"]:464            info = pick_info(465                info, [info["ch_names"].index(name) for name in inst["names"]]466            )467    else:468        info = inst.info469        inst_type = "data"470    logger.info(f"Computing rank from {inst_type} with rank={repr(rank)}")471 472    _validate_type(rank, (str, dict, None), "rank")473    if isinstance(rank, str):  # string, either 'info' or 'full'474        rank_type = "info"475        info_type = rank476        rank = dict()477    else:  # None or dict478        rank_type = "estimated"479        if rank is None:480            rank = dict()481 482    if on_few_samples is None:483        if inst_type != "covariance" and rank_type == "estimated":484            on_few_samples = "warn"485        else:486            on_few_samples = "ignore"487 488    simple_info = _simplify_info(info)489    picks_list = _picks_by_type(info, meg_combined=True, ref_meg=False, exclude="bads")490    for ch_type, picks in picks_list:491        est_verbose = None492        if ch_type in rank:493            # raise an error of user-supplied rank exceeds number of channels494            if rank[ch_type] > len(picks):495                raise ValueError(496                    f"rank[{repr(ch_type)}]={rank[ch_type]} exceeds the number"497                    f" of channels ({len(picks)})"498                )499            # special case: if whitening a covariance, check the passed rank500            # against the estimated one501            est_verbose = False502            if not (503                on_rank_mismatch != "ignore"504                and rank_type == "estimated"505                and ch_type == "meg"506                and isinstance(inst, Covariance)507                and not inst["diag"]508            ):509                continue510        ch_names = [info["ch_names"][pick] for pick in picks]511        n_chan = len(ch_names)512        if proj:513            proj_op, n_proj, _ = make_projector(info["projs"], ch_names)514        else:515            proj_op, n_proj = None, 0516        if log_ch_type is None:517            ch_type_ = ch_type.upper()518        else:519            ch_type_ = log_ch_type520        if rank_type == "info":521            # use info522            this_rank = _info_rank(info, ch_type, picks, info_type)523            if info_type != "full":524                this_rank -= n_proj525                logger.info(526                    f"    {ch_type_}: rank {this_rank} after "527                    f"{n_proj} projector{_pl(n_proj)} applied to "528                    f"{n_chan} channel{_pl(n_chan)}"529                )530            else:531                logger.info(f"    {ch_type_}: rank {this_rank} from info")532        else:533            # Use empirical estimation534            assert rank_type == "estimated"535            if isinstance(inst, BaseRaw | BaseEpochs):536                if isinstance(inst, BaseRaw):537                    data = inst.get_data(picks, reject_by_annotation="omit")538                else:  # isinstance(inst, BaseEpochs):539                    data = np.concatenate(inst.get_data(picks), axis=1)540                if proj:541                    data = np.dot(proj_op, data)542                this_rank = _estimate_rank_meeg_signals(543                    data,544                    pick_info(simple_info, picks),545                    scalings,546                    tol,547                    False,548                    tol_kind,549                    log_ch_type=log_ch_type,550                    on_few_samples=on_few_samples,551                )552            else:553                assert isinstance(inst, Covariance)554                if inst["diag"]:555                    this_rank = (inst["data"][picks] > 0).sum() - n_proj556                else:557                    data = inst["data"][picks][:, picks]558                    if proj:559                        data = np.dot(np.dot(proj_op, data), proj_op.T)560 561                    this_rank, sing = _estimate_rank_meeg_cov(562                        data,563                        pick_info(simple_info, picks),564                        scalings,565                        tol,566                        return_singular=True,567                        log_ch_type=log_ch_type,568                        on_few_samples=on_few_samples,569                        verbose=est_verbose,570                    )571                    if ch_type in rank:572                        ratio = sing[this_rank - 1] / sing[rank[ch_type] - 1]573                        if ratio > 100:574                            msg = (575                                f"The passed rank[{repr(ch_type)}]="576                                f"{rank[ch_type]} exceeds the estimated rank "577                                f"of the noise covariance ({this_rank}) "578                                f"leading to a potential increase in "579                                f"noise during whitening by a factor "580                                f"of {np.sqrt(ratio):0.1g}. Ensure that the "581                                f"rank correctly corresponds to that of the "582                                f"given noise covariance matrix."583                            )584                            _on_missing(on_rank_mismatch, msg, "on_rank_mismatch")585                        continue586            this_info_rank = _info_rank(info, ch_type, picks, "info")587            logger.info(588                f"    {ch_type_}: rank {this_rank} computed from "589                f"{n_chan} data channel{_pl(n_chan)} with "590                f"{n_proj} projector{_pl(n_proj)}"591            )592            if this_rank > this_info_rank:593                warn(594                    "Something went wrong in the data-driven estimation of the data "595                    "rank as it exceeds the theoretical rank from the info "596                    f"({this_rank} > {this_info_rank}). Consider setting rank "597                    'to "auto" or setting it explicitly as an integer.'598                )599        if ch_type not in rank:600            rank[ch_type] = int(this_rank)601 602    return rank603