Aluode/PerceptionLabPortable

1# Authors: The MNE-Python contributors.2# License: BSD-3-Clause3# Copyright the MNE-Python contributors.4 5import numbers6 7import numpy as np8from scipy.stats import pearsonr9from sklearn.base import (10    BaseEstimator,11    MetaEstimatorMixin,12    clone,13    is_regressor,14)15from sklearn.exceptions import NotFittedError16from sklearn.metrics import r2_score17 18from ..utils import _validate_type, fill_doc, pinv19from ._fixes import _check_n_features_3d, validate_data20from .base import _check_estimator, get_coef21from .time_delaying_ridge import TimeDelayingRidge22 23 24@fill_doc25class ReceptiveField(MetaEstimatorMixin, BaseEstimator):26    """Fit a receptive field model.27 28    This allows you to fit an encoding model (stimulus to brain) or a decoding29    model (brain to stimulus) using time-lagged input features (for example, a30    spectro- or spatio-temporal receptive field, or STRF)31    :footcite:`TheunissenEtAl2001,WillmoreSmyth2003,CrosseEtAl2016,HoldgrafEtAl2016`.32 33    Parameters34    ----------35    tmin : float36        The starting lag, in seconds (or samples if ``sfreq`` == 1).37    tmax : float38        The ending lag, in seconds (or samples if ``sfreq`` == 1).39        Must be >= tmin.40    sfreq : float41        The sampling frequency used to convert times into samples.42    feature_names : array, shape (n_features,) | None43        Names for input features to the model. If None, feature names will44        be auto-generated from the shape of input data after running `fit`.45    estimator : instance of sklearn.base.BaseEstimator | float | None46        The model used in fitting inputs and outputs. This can be any47        scikit-learn-style model that contains a fit and predict method. If a48        float is passed, it will be interpreted as the ``alpha`` parameter49        to be passed to a Ridge regression model. If `None`, then a Ridge50        regression model with an alpha of 0 will be used.51    fit_intercept : bool | None52        If True (default), the sample mean is removed before fitting.53        If ``estimator`` is a :class:`sklearn.base.BaseEstimator`,54        this must be None or match ``estimator.fit_intercept``.55    scoring : ['r2', 'corrcoef']56        Defines how predictions will be scored. Currently must be one of57        'r2' (coefficient of determination) or 'corrcoef' (the correlation58        coefficient).59    patterns : bool60        If True, inverse coefficients will be computed upon fitting using the61        covariance matrix of the inputs, and the cross-covariance of the62        inputs/outputs, according to :footcite:`HaufeEtAl2014`. Defaults to63        False.64    n_jobs : int | str65        Number of jobs to run in parallel. Can be 'cuda' if CuPy66        is installed properly and ``estimator is None``.67 68        .. versionadded:: 0.1869    edge_correction : bool70        If True (default), correct the autocorrelation coefficients for71        non-zero delays for the fact that fewer samples are available.72        Disabling this speeds up performance at the cost of accuracy73        depending on the relationship between epoch length and model74        duration. Only used if ``estimator`` is float or None.75 76        .. versionadded:: 0.1877 78    Attributes79    ----------80    coef_ : array, shape ([n_outputs, ]n_features, n_delays)81        The coefficients from the model fit, reshaped for easy visualization.82        During :meth:`mne.decoding.ReceptiveField.fit`, if ``y`` has one83        dimension (time), the ``n_outputs`` dimension here is omitted.84    patterns_ : array, shape ([n_outputs, ]n_features, n_delays)85        If fit, the inverted coefficients from the model.86    delays_ : array, shape (n_delays,), dtype int87        The delays used to fit the model, in indices. To return the delays88        in seconds, use ``self.delays_ / self.sfreq``89    valid_samples_ : slice90        The rows to keep during model fitting after removing rows with91        missing values due to time delaying. This can be used to get an92        output equivalent to using :func:`numpy.convolve` or93        :func:`numpy.correlate` with ``mode='valid'``.94 95    See Also96    --------97    mne.decoding.TimeDelayingRidge98 99    Notes100    -----101    For a causal system, the encoding model will have significant102    non-zero values only at positive lags. In other words, lags point103    backward in time relative to the input, so positive lags correspond104    to previous input time samples, while negative lags correspond to105    future input time samples.106 107    References108    ----------109    .. footbibliography::110    """  # noqa E501111 112    def __init__(113        self,114        tmin,115        tmax,116        sfreq,117        feature_names=None,118        estimator=None,119        fit_intercept=None,120        scoring="r2",121        patterns=False,122        n_jobs=None,123        edge_correction=True,124    ):125        self.tmin = tmin126        self.tmax = tmax127        self.sfreq = sfreq128        self.feature_names = feature_names129        self.estimator = estimator130        self.fit_intercept = fit_intercept131        self.scoring = scoring132        self.patterns = patterns133        self.n_jobs = n_jobs134        self.edge_correction = edge_correction135 136    def __repr__(self):  # noqa: D105137        s = f"tmin, tmax : ({self.tmin:.3f}, {self.tmax:.3f}), "138        estimator = self.estimator139        if not isinstance(estimator, str):140            estimator = type(self.estimator)141        s += f"estimator : {estimator}, "142        if hasattr(self, "coef_"):143            if self.feature_names is not None:144                feats = self.feature_names145                if len(feats) == 1:146                    s += f"feature: {feats[0]}, "147                else:148                    s += f"features : [{feats[0]}, ..., {feats[-1]}], "149            s += "fit: True"150        else:151            s += "fit: False"152        if hasattr(self, "scores_"):153            s += f"scored ({self.scoring})"154        return f"<ReceptiveField | {s}>"155 156    def __sklearn_tags__(self):157        """..."""158        from sklearn.utils import RegressorTags159 160        tags = super().__sklearn_tags__()161        tags.estimator_type = "regressor"162        tags.regressor_tags = RegressorTags()163        tags.input_tags.three_d_array = True164        tags.target_tags.one_d_labels = True165        tags.target_tags.multi_output = True166        tags.target_tags.required = True167        return tags168 169    def _delay_and_reshape(self, X, y=None):170        """Delay and reshape the variables."""171        if not isinstance(self.estimator_, TimeDelayingRidge):172            # X is now shape (n_times, n_epochs, n_feats, n_delays)173            X = _delay_time_series(174                X,175                self.tmin,176                self.tmax,177                self.sfreq_,178                fill_mean=self.fit_intercept_,179            )180            X = _reshape_for_est(X)181            # Concat times + epochs182            if y is not None:183                y = y.reshape(-1, y.shape[-1], order="F")184        return X, y185 186    def _check_data(self, X, y=None, reset=False):187        if reset:188            X, y = validate_data(189                self,190                X=X,191                y=y,192                reset=reset,193                validate_separately=(  # to take care of 3D y194                    dict(allow_nd=True, ensure_2d=False),195                    dict(allow_nd=True, ensure_2d=False),196                ),197            )198        else:199            X = validate_data(self, X=X, allow_nd=True, ensure_2d=False, reset=reset)200        _check_n_features_3d(self, X, reset)201        return X, y202 203    def _validate_params(self, X):204        if self.scoring not in _SCORERS.keys():205            raise ValueError(206                f"scoring must be one of {sorted(_SCORERS.keys())}, got {self.scoring}"207            )208        self.sfreq_ = float(self.sfreq)209        if self.tmin > self.tmax:210            raise ValueError(f"tmin ({self.tmin}) must be at most tmax ({self.tmax})")211 212    def fit(self, X, y):213        """Fit a receptive field model.214 215        Parameters216        ----------217        X : array, shape (n_times[, n_epochs], n_features)218            The input features for the model.219        y : array, shape (n_times[, n_epochs][, n_outputs])220            The output features for the model.221 222        Returns223        -------224        self : instance225            The instance so you can chain operations.226        """227        X, y = self._check_data(X, y, reset=True)228        self._validate_params(X)229        X, y, _, self._y_dim = self._check_dimensions(X, y)230 231        # Initialize delays232        self.delays_ = _times_to_delays(self.tmin, self.tmax, self.sfreq_)233 234        # Define the slice that we should use in the middle235        self.valid_samples_ = _delays_to_slice(self.delays_)236 237        if self.estimator is None or isinstance(self.estimator, numbers.Real):238            alpha = self.estimator if self.estimator is not None else 0.0239            if self.fit_intercept is None:240                self.fit_intercept_ = True241            else:242                self.fit_intercept_ = self.fit_intercept243            estimator = TimeDelayingRidge(244                self.tmin,245                self.tmax,246                self.sfreq_,247                alpha=alpha,248                fit_intercept=self.fit_intercept_,249                n_jobs=self.n_jobs,250                edge_correction=self.edge_correction,251            )252        elif is_regressor(self.estimator):253            estimator = clone(self.estimator)254            if (255                self.fit_intercept is not None256                and estimator.fit_intercept != self.fit_intercept257            ):258                raise ValueError(259                    f"Estimator fit_intercept ({estimator.fit_intercept}) != "260                    f"initialization fit_intercept ({self.fit_intercept}), initialize "261                    "ReceptiveField with the same fit_intercept value or use "262                    "fit_intercept=None"263                )264            self.fit_intercept_ = estimator.fit_intercept265        else:266            raise ValueError(267                "`estimator` must be a float or an instance of `BaseEstimator`, got "268                f"type {self.estimator}."269            )270        self.estimator_ = estimator271        del estimator272        _check_estimator(self.estimator_)273 274        # Create input features275        n_times, n_epochs, n_feats = X.shape276        n_outputs = y.shape[-1]277        n_delays = len(self.delays_)278 279        # Update feature names if we have none280        if (self.feature_names is not None) and (len(self.feature_names) != n_feats):281            raise ValueError(282                f"n_features in X does not match feature names ({n_feats} != "283                f"{len(self.feature_names)})"284            )285 286        # Create input features287        X, y = self._delay_and_reshape(X, y)288 289        self.estimator_.fit(X, y)290        coef = get_coef(self.estimator_, "coef_")  # (n_targets, n_features)291        shape = [n_feats, n_delays]292        if self._y_dim > 1:293            shape.insert(0, -1)294        self.coef_ = coef.reshape(shape)295 296        # Inverse-transform model weights297        if self.patterns:298            n_total_samples = n_times * n_epochs299            if n_total_samples < 2:300                raise ValueError(301                    "Cannot compute patterns with only one sample; "302                    f"got n_samples = {n_total_samples}."303                )304            if isinstance(self.estimator_, TimeDelayingRidge):305                cov_ = self.estimator_.cov_ / float(n_times * n_epochs - 1)306                y = y.reshape(-1, y.shape[-1], order="F")307            else:308                X = X - X.mean(0, keepdims=True)309                cov_ = np.cov(X.T)310            del X311 312            # Inverse output covariance313            if y.ndim == 2 and y.shape[1] != 1:314                y = y - y.mean(0, keepdims=True)315                inv_Y = pinv(np.cov(y.T))316            else:317                inv_Y = 1.0 / float(n_times * n_epochs - 1)318            del y319 320            # Inverse coef according to Haufe's method321            # patterns has shape (n_feats * n_delays, n_outputs)322            coef = np.reshape(self.coef_, (n_feats * n_delays, n_outputs))323            patterns = cov_.dot(coef.dot(inv_Y))324            self.patterns_ = patterns.reshape(shape)325 326        return self327 328    def predict(self, X):329        """Generate predictions with a receptive field.330 331        Parameters332        ----------333        X : array, shape (n_times[, n_epochs], n_channels)334            The input features for the model.335 336        Returns337        -------338        y_pred : array, shape (n_times[, n_epochs][, n_outputs])339            The output predictions. "Note that valid samples (those340            unaffected by edge artifacts during the time delaying step) can341            be obtained using ``y_pred[rf.valid_samples_]``.342        """343        if not hasattr(self, "delays_"):344            raise NotFittedError("Estimator has not been fit yet.")345 346        X, _ = self._check_data(X)347        X, _, X_dim = self._check_dimensions(X, None, predict=True)[:3]348 349        del _350        # convert to sklearn and back351        pred_shape = X.shape[:-1]352        if self._y_dim > 1:353            pred_shape = pred_shape + (self.coef_.shape[0],)354        X, _ = self._delay_and_reshape(X)355        y_pred = self.estimator_.predict(X)356        y_pred = y_pred.reshape(pred_shape, order="F")357        shape = list(y_pred.shape)358        if X_dim <= 2:359            shape.pop(1)  # epochs360            extra = 0361        else:362            extra = 1363        shape = shape[: self._y_dim + extra]364        y_pred.shape = shape365        return y_pred366 367    def score(self, X, y):368        """Score predictions generated with a receptive field.369 370        This calls ``self.predict``, then masks the output of this371        and ``y` with ``self.valid_samples_``. Finally, it passes372        this to a :mod:`sklearn.metrics` scorer.373 374        Parameters375        ----------376        X : array, shape (n_times[, n_epochs], n_channels)377            The input features for the model.378        y : array, shape (n_times[, n_epochs][, n_outputs])379            Used for scikit-learn compatibility.380 381        Returns382        -------383        scores : list of float, shape (n_outputs,)384            The scores estimated by the model for each output (e.g. mean385            R2 of ``predict(X)``).386        """387        # Create our scoring object388        scorer_ = _SCORERS[self.scoring]389 390        # Generate predictions, then reshape so we can mask time391        X, y = self._check_dimensions(X, y, predict=True)[:2]392        n_times, n_epochs, n_outputs = y.shape393        y_pred = self.predict(X)394        y_pred = y_pred[self.valid_samples_]395        y = y[self.valid_samples_]396 397        # Re-vectorize and call scorer398        y = y.reshape([-1, n_outputs], order="F")399        y_pred = y_pred.reshape([-1, n_outputs], order="F")400        assert y.shape == y_pred.shape401        scores = scorer_(y, y_pred, multioutput="raw_values")402        return scores403 404    def _check_dimensions(self, X, y, predict=False):405        _validate_type(X, "array-like", "X")406        _validate_type(y, ("array-like", None), "y")407        X_dim = X.ndim408        y_dim = y.ndim if y is not None else 0409        if X_dim == 2:410            # Ensure we have a 3D input by adding singleton epochs dimension411            X = X[:, np.newaxis, :]412            if y is not None:413                if y_dim == 1:414                    y = y[:, np.newaxis, np.newaxis]  # epochs, outputs415                elif y_dim == 2:416                    y = y[:, np.newaxis, :]  # epochs417                else:418                    raise ValueError(419                        "y must be shape (n_times[, n_epochs][,n_outputs], got "420                        f"{y.shape}"421                    )422        elif X.ndim == 3:423            if y is not None:424                if y.ndim == 2:425                    y = y[:, :, np.newaxis]  # Add an outputs dim426                elif y.ndim != 3:427                    raise ValueError(428                        "If X has 3 dimensions, y must have 2 or 3 dimensions"429                    )430        else:431            raise ValueError(432                "X must be shape (n_times[, n_epochs], n_features), "433                f"got {X.shape}. Reshape your data to 2D or 3D "434                "(e.g., array.reshape(-1, 1) for a single feature, "435                "or array.reshape(1, -1) for a single sample)."436            )437        if y is not None:438            if X.shape[0] != y.shape[0]:439                raise ValueError(440                    f"X and y do not have the same n_times\n{X.shape[0]} != "441                    f"{y.shape[0]}"442                )443            if X.shape[1] != y.shape[1]:444                raise ValueError(445                    f"X and y do not have the same n_epochs\n{X.shape[1]} != "446                    f"{y.shape[1]}"447                )448            if predict and y.shape[-1] not in (len(self.estimator_.coef_), 1):449                raise ValueError(450                    "Number of outputs does not match estimator coefficients dimensions"451                )452        return X, y, X_dim, y_dim453 454 455def _delay_time_series(X, tmin, tmax, sfreq, fill_mean=False):456    """Return a time-lagged input time series.457 458    Parameters459    ----------460    X : array, shape (n_times[, n_epochs], n_features)461        The time series to delay. Must be 2D or 3D.462    tmin : int | float463        The starting lag.464    tmax : int | float465        The ending lag.466        Must be >= tmin.467    sfreq : int | float468        The sampling frequency of the series. Defaults to 1.0.469    fill_mean : bool470        If True, the fill value will be the mean along the time dimension471        of the feature, and each cropped and delayed segment of data472        will be shifted to have the same mean value (ensuring that mean473        subtraction works properly). If False, the fill value will be zero.474 475    Returns476    -------477    delayed : array, shape(n_times[, n_epochs][, n_features], n_delays)478        The delayed data. It has the same shape as X, with an extra dimension479        appended to the end.480 481    Examples482    --------483    >>> tmin, tmax = -0.1, 0.2484    >>> sfreq = 10.485    >>> x = np.arange(1, 6)486    >>> x_del = _delay_time_series(x, tmin, tmax, sfreq)487    >>> print(x_del)  # doctest:+SKIP488    [[2. 1. 0. 0.]489     [3. 2. 1. 0.]490     [4. 3. 2. 1.]491     [5. 4. 3. 2.]492     [0. 5. 4. 3.]]493    """494    _check_delayer_params(tmin, tmax, sfreq)495    delays = _times_to_delays(tmin, tmax, sfreq)496    # Iterate through indices and append497    delayed = np.zeros(X.shape + (len(delays),))498    if fill_mean:499        mean_value = X.mean(axis=0)500        if X.ndim == 3:501            mean_value = np.mean(mean_value, axis=0)502        delayed[:] = mean_value[:, np.newaxis]503    for ii, ix_delay in enumerate(delays):504        # Create zeros to populate w/ delays505        if ix_delay < 0:506            out = delayed[:ix_delay, ..., ii]507            use_X = X[-ix_delay:]508        elif ix_delay > 0:509            out = delayed[ix_delay:, ..., ii]510            use_X = X[:-ix_delay]511        else:  # == 0512            out = delayed[..., ii]513            use_X = X514        out[:] = use_X515        if fill_mean:516            out[:] += mean_value - use_X.mean(axis=0)517    return delayed518 519 520def _times_to_delays(tmin, tmax, sfreq):521    """Convert a tmin/tmax in seconds to delays."""522    # Convert seconds to samples523    delays = np.arange(int(np.round(tmin * sfreq)), int(np.round(tmax * sfreq) + 1))524    return delays525 526 527def _delays_to_slice(delays):528    """Find the slice to be taken in order to remove missing values."""529    # Negative values == cut off rows at the end530    min_delay = None if delays[-1] <= 0 else delays[-1]531    # Positive values == cut off rows at the end532    max_delay = None if delays[0] >= 0 else delays[0]533    return slice(min_delay, max_delay)534 535 536def _check_delayer_params(tmin, tmax, sfreq):537    """Check delayer input parameters. For future custom delay support."""538    _validate_type(sfreq, "numeric", "`sfreq`")539 540    for tlim in (tmin, tmax):541        _validate_type(tlim, "numeric", "tmin/tmax")542    if not tmin <= tmax:543        raise ValueError("tmin must be <= tmax")544 545 546def _reshape_for_est(X_del):547    """Convert X_del to a sklearn-compatible shape."""548    n_times, n_epochs, n_feats, n_delays = X_del.shape549    X_del = X_del.reshape(n_times, n_epochs, -1)  # concatenate feats550    X_del = X_del.reshape(n_times * n_epochs, -1, order="F")551    return X_del552 553 554# Create a correlation scikit-learn-style scorer555def _corr_score(y_true, y, multioutput=None):556    assert multioutput == "raw_values"557    for this_y in (y_true, y):558        if this_y.ndim != 2:559            raise ValueError(560                f"inputs must be shape (samples, outputs), got {this_y.shape}"561            )562    return np.array([pearsonr(y_true[:, ii], y[:, ii])[0] for ii in range(y.shape[-1])])563 564 565def _r2_score(y_true, y, multioutput=None):566    return r2_score(y_true, y, multioutput=multioutput)567 568 569_SCORERS = {"r2": _r2_score, "corrcoef": _corr_score}570