Aluode/PerceptionLabPortable

1# Authors: The MNE-Python contributors.2# License: BSD-3-Clause3# Copyright the MNE-Python contributors.4 5import numpy as np6 7from ..fixes import rng_uniform8from ..label import Label9from ..source_estimate import SourceEstimate, VolSourceEstimate10from ..source_space._source_space import _ensure_src11from ..surface import _compute_nearest12from ..utils import (13    _check_option,14    _ensure_events,15    _ensure_int,16    _validate_type,17    check_random_state,18    fill_doc,19    warn,20)21 22 23@fill_doc24def select_source_in_label(25    src,26    label,27    random_state=None,28    location="random",29    subject=None,30    subjects_dir=None,31    surf="sphere",32):33    """Select source positions using a label.34 35    Parameters36    ----------37    src : list of dict38        The source space.39    label : Label40        The label.41    %(random_state)s42    location : str43        The label location to choose. Can be 'random' (default) or 'center'44        to use :func:`mne.Label.center_of_mass` (restricting to vertices45        both in the label and in the source space). Note that for 'center'46        mode the label values are used as weights.47 48        .. versionadded:: 0.1349    subject : str | None50        The subject the label is defined for.51        Only used with ``location='center'``.52 53        .. versionadded:: 0.1354    %(subjects_dir)s55 56        .. versionadded:: 0.1357    surf : str58        The surface to use for Euclidean distance center of mass59        finding. The default here is "sphere", which finds the center60        of mass on the spherical surface to help avoid potential issues61        with cortical folding.62 63        .. versionadded:: 0.1364 65    Returns66    -------67    lh_vertno : list68        Selected source coefficients on the left hemisphere.69    rh_vertno : list70        Selected source coefficients on the right hemisphere.71    """72    lh_vertno = list()73    rh_vertno = list()74    _check_option("location", location, ["random", "center"])75 76    rng = check_random_state(random_state)77    if label.hemi == "lh":78        vertno = lh_vertno79        hemi_idx = 080    else:81        vertno = rh_vertno82        hemi_idx = 183    src_sel = np.intersect1d(src[hemi_idx]["vertno"], label.vertices)84    if location == "random":85        idx = src_sel[rng_uniform(rng)(0, len(src_sel), 1)[0]]86    else:  # 'center'87        idx = label.center_of_mass(88            subject, restrict_vertices=src_sel, subjects_dir=subjects_dir, surf=surf89        )90    vertno.append(idx)91    return lh_vertno, rh_vertno92 93 94@fill_doc95def simulate_sparse_stc(96    src,97    n_dipoles,98    times,99    data_fun=lambda t: 1e-7 * np.sin(20 * np.pi * t),100    labels=None,101    random_state=None,102    location="random",103    subject=None,104    subjects_dir=None,105    surf="sphere",106):107    """Generate sparse (n_dipoles) sources time courses from data_fun.108 109    This function randomly selects ``n_dipoles`` vertices in the whole110    cortex or one single vertex (randomly in or in the center of) each111    label if ``labels is not None``. It uses ``data_fun`` to generate112    waveforms for each vertex.113 114    Parameters115    ----------116    src : instance of SourceSpaces117        The source space.118    n_dipoles : int119        Number of dipoles to simulate.120    times : array121        Time array.122    data_fun : callable123        Function to generate the waveforms. The default is a 100 nAm, 10 Hz124        sinusoid as ``1e-7 * np.sin(20 * pi * t)``. The function should take125        as input the array of time samples in seconds and return an array of126        the same length containing the time courses.127    labels : None | list of Label128        The labels. The default is None, otherwise its size must be n_dipoles.129    %(random_state)s130    location : str131        The label location to choose. Can be ``'random'`` (default) or132        ``'center'`` to use :func:`mne.Label.center_of_mass`. Note that for133        ``'center'`` mode the label values are used as weights.134 135        .. versionadded:: 0.13136    subject : str | None137        The subject the label is defined for.138        Only used with ``location='center'``.139 140        .. versionadded:: 0.13141    %(subjects_dir)s142 143        .. versionadded:: 0.13144    surf : str145        The surface to use for Euclidean distance center of mass146        finding. The default here is "sphere", which finds the center147        of mass on the spherical surface to help avoid potential issues148        with cortical folding.149 150        .. versionadded:: 0.13151 152    Returns153    -------154    stc : SourceEstimate155        The generated source time courses.156 157    See Also158    --------159    simulate_raw160    simulate_evoked161    simulate_stc162 163    Notes164    -----165    .. versionadded:: 0.10.0166    """167    rng = check_random_state(random_state)168    src = _ensure_src(src, verbose=False)169    subject_src = src._subject170    if subject is None:171        subject = subject_src172    elif subject_src is not None and subject != subject_src:173        raise ValueError(174            f"subject argument ({subject}) did not match the source "175            f"space subject_his_id ({subject_src})"176        )177    data = np.zeros((n_dipoles, len(times)))178    for i_dip in range(n_dipoles):179        data[i_dip, :] = data_fun(times)180 181    if labels is None:182        # can be vol or surface source space183        offsets = np.linspace(0, n_dipoles, len(src) + 1).astype(int)184        n_dipoles_ss = np.diff(offsets)185        # don't use .choice b/c not on old numpy186        vs = [187            s["vertno"][np.sort(rng.permutation(np.arange(s["nuse"]))[:n])]188            for n, s in zip(n_dipoles_ss, src)189        ]190        datas = data191    elif n_dipoles > len(labels):192        raise ValueError(193            f"Number of labels ({len(labels)}) smaller than n_dipoles ({n_dipoles:d}) "194            "is not allowed."195        )196    else:197        if n_dipoles != len(labels):198            warn(199                "The number of labels is different from the number of "200                f"dipoles. {min(n_dipoles, len(labels))} dipole(s) will be generated."201            )202        labels = labels[:n_dipoles] if n_dipoles < len(labels) else labels203 204        vertno = [[], []]205        lh_data = [np.empty((0, data.shape[1]))]206        rh_data = [np.empty((0, data.shape[1]))]207        for i, label in enumerate(labels):208            lh_vertno, rh_vertno = select_source_in_label(209                src, label, rng, location, subject, subjects_dir, surf210            )211            vertno[0] += lh_vertno212            vertno[1] += rh_vertno213            if len(lh_vertno) != 0:214                lh_data.append(data[i][np.newaxis])215            elif len(rh_vertno) != 0:216                rh_data.append(data[i][np.newaxis])217            else:218                raise ValueError("No vertno found.")219        vs = [np.array(v) for v in vertno]220        datas = [np.concatenate(d) for d in [lh_data, rh_data]]221        # need to sort each hemi by vertex number222        for ii in range(2):223            order = np.argsort(vs[ii])224            vs[ii] = vs[ii][order]225            if len(order) > 0:  # fix for old numpy226                datas[ii] = datas[ii][order]227        datas = np.concatenate(datas)228 229    tmin, tstep = times[0], np.diff(times[:2])[0]230    assert datas.shape == data.shape231    cls = SourceEstimate if len(vs) == 2 else VolSourceEstimate232    stc = cls(datas, vertices=vs, tmin=tmin, tstep=tstep, subject=subject)233    return stc234 235 236def simulate_stc(237    src, labels, stc_data, tmin, tstep, value_fun=None, allow_overlap=False238):239    """Simulate sources time courses from waveforms and labels.240 241    This function generates a source estimate with extended sources by242    filling the labels with the waveforms given in stc_data.243 244    Parameters245    ----------246    src : instance of SourceSpaces247        The source space.248    labels : list of Label249        The labels.250    stc_data : array, shape (n_labels, n_times)251        The waveforms.252    tmin : float253        The beginning of the timeseries.254    tstep : float255        The time step (1 / sampling frequency).256    value_fun : callable | None257        Function to apply to the label values to obtain the waveform258        scaling for each vertex in the label. If None (default), uniform259        scaling is used.260    allow_overlap : bool261        Allow overlapping labels or not. Default value is False.262 263        .. versionadded:: 0.18264 265    Returns266    -------267    stc : SourceEstimate268        The generated source time courses.269 270    See Also271    --------272    simulate_raw273    simulate_evoked274    simulate_sparse_stc275    """276    if len(labels) != len(stc_data):277        raise ValueError("labels and stc_data must have the same length")278 279    vertno = [[], []]280    stc_data_extended = [[], []]281    hemi_to_ind = {"lh": 0, "rh": 1}282    for i, label in enumerate(labels):283        hemi_ind = hemi_to_ind[label.hemi]284        src_sel = np.intersect1d(src[hemi_ind]["vertno"], label.vertices)285        if len(src_sel) == 0:286            idx = src[hemi_ind]["inuse"].astype("bool")287            xhs = src[hemi_ind]["rr"][idx]288            rr = src[hemi_ind]["rr"][label.vertices]289            closest_src = _compute_nearest(xhs, rr)290            src_sel = src[hemi_ind]["vertno"][np.unique(closest_src)]291 292        if value_fun is not None:293            idx_sel = np.searchsorted(label.vertices, src_sel)294            values_sel = np.array([value_fun(v) for v in label.values[idx_sel]])295 296            data = np.outer(values_sel, stc_data[i])297        else:298            data = np.tile(stc_data[i], (len(src_sel), 1))299        # If overlaps are allowed, deal with them300        if allow_overlap:301            # Search for duplicate vertex indices302            # in the existing vertex matrix vertex.303            duplicates = []304            for src_ind, vertex_ind in enumerate(src_sel):305                ind = np.where(vertex_ind == vertno[hemi_ind])[0]306                if len(ind) > 0:307                    assert len(ind) == 1308                    # Add the new data to the existing one309                    stc_data_extended[hemi_ind][ind[0]] += data[src_ind]310                    duplicates.append(src_ind)311            # Remove the duplicates from both data and selected vertices312            data = np.delete(data, duplicates, axis=0)313            src_sel = list(np.delete(np.array(src_sel), duplicates))314        # Extend the existing list instead of appending it so that we can315        # index its elements316        vertno[hemi_ind].extend(src_sel)317        stc_data_extended[hemi_ind].extend(np.atleast_2d(data))318 319    vertno = [np.array(v) for v in vertno]320    if not allow_overlap:321        for v, hemi in zip(vertno, ("left", "right")):322            d = len(v) - len(np.unique(v))323            if d > 0:324                raise RuntimeError(325                    f"Labels had {d} overlaps in the {hemi} "326                    "hemisphere, they must be non-overlapping"327                )328    # the data is in the order left, right329    data = list()330    for i in range(2):331        if len(stc_data_extended[i]) != 0:332            stc_data_extended[i] = np.vstack(stc_data_extended[i])333            # Order the indices of each hemisphere334            idx = np.argsort(vertno[i])335            data.append(stc_data_extended[i][idx])336            vertno[i] = vertno[i][idx]337 338    stc = SourceEstimate(339        np.concatenate(data),340        vertices=vertno,341        tmin=tmin,342        tstep=tstep,343        subject=src._subject,344    )345    return stc346 347 348class SourceSimulator:349    """Class to generate simulated Source Estimates.350 351    Parameters352    ----------353    src : instance of SourceSpaces354        Source space.355    tstep : float356        Time step between successive samples in data. Default is 0.001 s.357    duration : float | None358        Time interval during which the simulation takes place in seconds.359        If None, it is computed using existing events and waveform lengths.360    first_samp : int361        First sample from which the simulation takes place, as an integer.362        Comparable to the :term:`first_samp` property of `~mne.io.Raw` objects.363        Default is 0.364 365    Attributes366    ----------367    duration : float368        The duration of the simulation in seconds.369    n_times : int370        The number of time samples of the simulation.371    """372 373    def __init__(self, src, tstep=1e-3, duration=None, first_samp=0):374        if duration is not None and duration < tstep:375            raise ValueError("duration must be None or >= tstep.")376        self.first_samp = _ensure_int(first_samp, "first_samp")377        self._src = src378        self._tstep = tstep379        self._labels = []380        self._waveforms = []381        self._events = np.empty((0, 3), dtype=int)382        self._duration = duration  # if not None, sets # samples383        self._last_samples = []384        self._chk_duration = 1000385 386    @property387    def duration(self):388        """Duration of the simulation in same units as tstep."""389        if self._duration is not None:390            return self._duration391        return self.n_times * self._tstep392 393    @property394    def n_times(self):395        """Number of time samples in the simulation."""396        if self._duration is not None:397            return int(self._duration / self._tstep)398        ls = self.first_samp399        if len(self._last_samples) > 0:400            ls = np.max(self._last_samples)401        return ls - self.first_samp + 1  # >= 1402 403    @property404    def last_samp(self):405        return self.first_samp + self.n_times - 1406 407    def add_data(self, label, waveform, events):408        """Add data to the simulation.409 410        Data should be added in the form of a triplet of411        Label (Where) - Waveform(s) (What) - Event(s) (When)412 413        Parameters414        ----------415        label : instance of Label416            The label (as created for example by mne.read_label). If the label417            does not match any sources in the SourceEstimate, a ValueError is418            raised.419        waveform : array, shape (n_times,) or (n_events, n_times) | list420            The waveform(s) describing the activity on the label vertices.421            If list, it must have the same length as events.422        events : array of int, shape (n_events, 3)423            Events associated to the waveform(s) to specify when the activity424            should occur.425        """426        _validate_type(label, Label, "label")427 428        # If it is not a list then make it one429        if not isinstance(waveform, list) and np.ndim(waveform) == 2:430            waveform = list(waveform)431        if not isinstance(waveform, list) and np.ndim(waveform) == 1:432            waveform = [waveform]433        if len(waveform) == 1:434            waveform = waveform * len(events)435        # The length is either equal to the length of events, or 1436        if len(waveform) != len(events):437            raise ValueError(438                "Number of waveforms and events should match or "439                f"there should be a single waveform ({len(waveform)} != {len(events)})."440            )441        events = _ensure_events(events).astype(np.int64)442        # Update the last sample possible based on events + waveforms443        self._labels.extend([label] * len(events))444        self._waveforms.extend(waveform)445        self._events = np.concatenate([self._events, events])446        assert self._events.dtype == np.int64447        # First sample per waveform is the first column of events448        # Last is computed below449        self._last_samples = np.array(450            [self._events[i, 0] + len(w) - 1 for i, w in enumerate(self._waveforms)]451        )452 453    def get_stim_channel(self, start_sample=0, stop_sample=None):454        """Get the stim channel from the provided data.455 456        Returns the stim channel data according to the simulation parameters457        which should be added through the add_data method. If both start_sample458        and stop_sample are not specified, the entire duration is used.459 460        Parameters461        ----------462        start_sample : int463            First sample in chunk. Default is the value of the ``first_samp``464            attribute.465        stop_sample : int | None466            The final sample of the returned stc. If None, then all samples467            from start_sample onward are returned.468 469        Returns470        -------471        stim_data : ndarray of int, shape (n_samples,)472            The stimulation channel data.473        """474        if start_sample is None:475            start_sample = self.first_samp476        if stop_sample is None:477            stop_sample = start_sample + self.n_times - 1478        elif stop_sample < start_sample:479            raise ValueError("Argument start_sample must be >= stop_sample.")480        n_samples = stop_sample - start_sample + 1481 482        # Initialize the stim data array483        stim_data = np.zeros(n_samples, dtype=np.int64)484 485        # Select only events in the time chunk486        stim_ind = np.where(487            np.logical_and(488                self._events[:, 0] >= start_sample, self._events[:, 0] < stop_sample489            )490        )[0]491 492        if len(stim_ind) > 0:493            relative_ind = self._events[stim_ind, 0] - start_sample494            stim_data[relative_ind] = self._events[stim_ind, 2]495 496        return stim_data497 498    def get_stc(self, start_sample=None, stop_sample=None):499        """Simulate a SourceEstimate from the provided data.500 501        Returns a SourceEstimate object constructed according to the simulation502        parameters which should be added through function add_data. If both503        start_sample and stop_sample are not specified, the entire duration is504        used.505 506        Parameters507        ----------508        start_sample : int | None509            First sample in chunk. If ``None`` the value of the ``first_samp``510            attribute is used. Defaults to ``None``.511        stop_sample : int | None512            The final sample of the returned STC. If ``None``, then all samples513            past ``start_sample`` are returned.514 515        Returns516        -------517        stc : SourceEstimate object518            The generated source time courses.519        """520        if len(self._labels) == 0:521            raise ValueError(522                "No simulation parameters were found. Please use "523                "function add_data to add simulation parameters."524            )525        if start_sample is None:526            start_sample = self.first_samp527        if stop_sample is None:528            stop_sample = start_sample + self.n_times - 1529        elif stop_sample < start_sample:530            raise ValueError("start_sample must be >= stop_sample.")531        n_samples = stop_sample - start_sample + 1532 533        # Initialize the stc_data array to span all possible samples534        stc_data = np.zeros((len(self._labels), n_samples))535 536        # Select only the events that fall within the span537        ind = np.where(538            np.logical_and(539                self._last_samples >= start_sample, self._events[:, 0] <= stop_sample540            )541        )[0]542 543        # Loop only over the items that are in the time span544        subset_waveforms = [self._waveforms[i] for i in ind]545        for i, (waveform, event) in enumerate(zip(subset_waveforms, self._events[ind])):546            # We retrieve the first and last sample of each waveform547            # According to the corresponding event548            wf_start = event[0]549            wf_stop = self._last_samples[ind[i]]550 551            # Recover the indices of the event that should be in the chunk552            waveform_ind = np.isin(553                np.arange(wf_start, wf_stop + 1),554                np.arange(start_sample, stop_sample + 1),555            )556 557            # Recover the indices that correspond to the overlap558            stc_ind = np.isin(559                np.arange(start_sample, stop_sample + 1),560                np.arange(wf_start, wf_stop + 1),561            )562 563            # add the resulting waveform chunk to the corresponding label564            stc_data[ind[i]][stc_ind] += waveform[waveform_ind]565 566        start_sample -= self.first_samp  # STC sample ref is 0567        stc = simulate_stc(568            self._src,569            self._labels,570            stc_data,571            start_sample * self._tstep,572            self._tstep,573            allow_overlap=True,574        )575 576        return stc577 578    def __iter__(self):579        """Iterate over 1 second STCs."""580        # Arbitrary chunk size, can be modified later to something else.581        # Loop over chunks of 1 second - or, maximum sample size.582        # Can be modified to a different value.583        last_sample = self.last_samp584        for start_sample in range(self.first_samp, last_sample + 1, self._chk_duration):585            stop_sample = min(start_sample + self._chk_duration - 1, last_sample)586            yield (587                self.get_stc(start_sample, stop_sample),588                self.get_stim_channel(start_sample, stop_sample),589            )590