Aluode/PerceptionLabPortable
0
1# Authors: The MNE-Python contributors.2# License: BSD-3-Clause3# Copyright the MNE-Python contributors.4 5import contextlib6import copy7import os.path as op8from types import GeneratorType9 10import numpy as np11from scipy import sparse12from scipy.spatial.distance import cdist, pdist13 14from ._fiff.constants import FIFF15from ._fiff.meas_info import Info16from ._fiff.pick import _picks_to_idx, pick_types17from ._freesurfer import _get_atlas_values, _get_mri_info_data, read_freesurfer_lut18from .baseline import rescale19from .cov import Covariance20from .evoked import _get_peak21from .filter import FilterMixin, _check_fun, resample22from .fixes import _eye_array, _safe_svd23from .parallel import parallel_func24from .source_space._source_space import (25 SourceSpaces,26 _check_volume_labels,27 _ensure_src,28 _ensure_src_subject,29 _get_morph_src_reordering,30 _get_src_nn,31 get_decimated_surfaces,32)33from .surface import _get_ico_surface, _project_onto_surface, mesh_edges, read_surface34from .transforms import _get_trans, apply_trans35from .utils import (36 TimeMixin,37 _build_data_frame,38 _check_fname,39 _check_option,40 _check_pandas_index_arguments,41 _check_pandas_installed,42 _check_preload,43 _check_src_normal,44 _check_stc_units,45 _check_subject,46 _check_time_format,47 _convert_times,48 _ensure_int,49 _import_h5io_funcs,50 _import_nibabel,51 _path_like,52 _pl,53 _time_mask,54 _validate_type,55 copy_function_doc_to_method_doc,56 fill_doc,57 get_subjects_dir,58 logger,59 object_size,60 sizeof_fmt,61 verbose,62 warn,63)64from .viz import (65 plot_source_estimates,66 plot_vector_source_estimates,67 plot_volume_source_estimates,68)69 70 71def _read_stc(filename):72 """Aux Function."""73 with open(filename, "rb") as fid:74 buf = fid.read()75 76 stc = dict()77 offset = 078 num_bytes = 479 80 # read tmin in ms81 stc["tmin"] = (82 float(np.frombuffer(buf, dtype=">f4", count=1, offset=offset).item()) / 1000.083 )84 offset += num_bytes85 86 # read sampling rate in ms87 stc["tstep"] = (88 float(np.frombuffer(buf, dtype=">f4", count=1, offset=offset).item()) / 1000.089 )90 offset += num_bytes91 92 # read number of vertices/sources93 vertices_n = int(np.frombuffer(buf, dtype=">u4", count=1, offset=offset).item())94 offset += num_bytes95 96 # read the source vector97 stc["vertices"] = np.frombuffer(buf, dtype=">u4", count=vertices_n, offset=offset)98 offset += num_bytes * vertices_n99 100 # read the number of timepts101 data_n = int(np.frombuffer(buf, dtype=">u4", count=1, offset=offset).item())102 offset += num_bytes103 104 if (105 vertices_n106 and ( # vertices_n can be 0 (empty stc)107 (len(buf) // 4 - 4 - vertices_n) % (data_n * vertices_n)108 )109 != 0110 ):111 raise ValueError("incorrect stc file size")112 113 # read the data matrix114 stc["data"] = np.frombuffer(115 buf, dtype=">f4", count=vertices_n * data_n, offset=offset116 )117 stc["data"] = stc["data"].reshape([data_n, vertices_n]).T118 119 return stc120 121 122def _write_stc(filename, tmin, tstep, vertices, data):123 """Write an STC file.124 125 Parameters126 ----------127 filename : path-like128 The name of the STC file.129 tmin : float130 The first time point of the data in seconds.131 tstep : float132 Time between frames in seconds.133 vertices : array of integers134 Vertex indices (0 based).135 data : 2D array136 The data matrix (nvert * ntime).137 """138 with open(filename, "wb") as fid:139 # write start time in ms140 fid.write(np.array(1000 * tmin, dtype=">f4").tobytes())141 # write sampling rate in ms142 fid.write(np.array(1000 * tstep, dtype=">f4").tobytes())143 # write number of vertices144 fid.write(np.array(vertices.shape[0], dtype=">u4").tobytes())145 # write the vertex indices146 fid.write(np.array(vertices, dtype=">u4").tobytes())147 # write the number of timepts148 fid.write(np.array(data.shape[1], dtype=">u4").tobytes())149 # write the data150 fid.write(np.array(data.T, dtype=">f4").tobytes())151 152 153def _read_3(fid):154 """Read 3 byte integer from file."""155 data = np.fromfile(fid, dtype=np.uint8, count=3).astype(np.int32)156 157 out = np.left_shift(data[0], 16) + np.left_shift(data[1], 8) + data[2]158 159 return out160 161 162def _read_w(filename):163 """Read a w file.164 165 w files contain activations or source reconstructions for a single time166 point.167 168 Parameters169 ----------170 filename : path-like171 The name of the w file.172 173 Returns174 -------175 data: dict176 The w structure. It has the following keys:177 vertices vertex indices (0 based)178 data The data matrix (nvert long)179 """180 with open(filename, "rb", buffering=0) as fid: # buffering=0 for np bug181 # skip first 2 bytes182 fid.read(2)183 184 # read number of vertices/sources (3 byte integer)185 vertices_n = int(_read_3(fid))186 187 vertices = np.zeros((vertices_n), dtype=np.int32)188 data = np.zeros((vertices_n), dtype=np.float32)189 190 # read the vertices and data191 for i in range(vertices_n):192 vertices[i] = _read_3(fid)193 data[i] = np.fromfile(fid, dtype=">f4", count=1).item()194 195 w = dict()196 w["vertices"] = vertices197 w["data"] = data198 199 return w200 201 202def _write_3(fid, val):203 """Write 3 byte integer to file."""204 f_bytes = np.zeros((3), dtype=np.uint8)205 f_bytes[0] = (val >> 16) & 255206 f_bytes[1] = (val >> 8) & 255207 f_bytes[2] = val & 255208 fid.write(f_bytes.tobytes())209 210 211def _write_w(filename, vertices, data):212 """Write a w file.213 214 w files contain activations or source reconstructions for a single time215 point.216 217 Parameters218 ----------219 filename: path-like220 The name of the w file.221 vertices: array of int222 Vertex indices (0 based).223 data: 1D array224 The data array (nvert).225 """226 assert len(vertices) == len(data)227 228 with open(filename, "wb") as fid:229 # write 2 zero bytes230 fid.write(np.zeros((2), dtype=np.uint8).tobytes())231 232 # write number of vertices/sources (3 byte integer)233 vertices_n = len(vertices)234 _write_3(fid, vertices_n)235 236 # write the vertices and data237 for i in range(vertices_n):238 _write_3(fid, vertices[i])239 # XXX: without float() endianness is wrong, not sure why240 fid.write(np.array(float(data[i]), dtype=">f4").tobytes())241 242 243def read_source_estimate(fname, subject=None):244 """Read a source estimate object.245 246 Parameters247 ----------248 fname : path-like249 Path to (a) source-estimate file(s).250 subject : str | None251 Name of the subject the source estimate(s) is (are) from.252 It is good practice to set this attribute to avoid combining253 incompatible labels and SourceEstimates (e.g., ones from other254 subjects). Note that due to file specification limitations, the255 subject name isn't saved to or loaded from files written to disk.256 257 Returns258 -------259 stc : SourceEstimate | VectorSourceEstimate | VolSourceEstimate | MixedSourceEstimate260 The source estimate object loaded from file.261 262 Notes263 -----264 - for volume source estimates, ``fname`` should provide the path to a265 single file named ``'*-vl.stc``` or ``'*-vol.stc'``266 - for surface source estimates, ``fname`` should either provide the267 path to the file corresponding to a single hemisphere (``'*-lh.stc'``,268 ``'*-rh.stc'``) or only specify the asterisk part in these patterns. In269 any case, the function expects files for both hemisphere with names270 following this pattern.271 - for vector surface source estimates, only HDF5 files are supported.272 - for mixed source estimates, only HDF5 files are supported.273 - for single time point ``.w`` files, ``fname`` should follow the same274 pattern as for surface estimates, except that files are named275 ``'*-lh.w'`` and ``'*-rh.w'``.276 """ # noqa: E501277 fname_arg = fname278 279 # expand `~` without checking whether the file actually exists – we'll280 # take care of that later, as it's complicated by the different suffixes281 # STC files can have282 fname = str(_check_fname(fname=fname, overwrite="read", must_exist=False))283 284 # make sure corresponding file(s) can be found285 ftype = None286 if op.exists(fname):287 if fname.endswith(("-vl.stc", "-vol.stc", "-vl.w", "-vol.w")):288 ftype = "volume"289 elif fname.endswith(".stc"):290 ftype = "surface"291 if fname.endswith(("-lh.stc", "-rh.stc")):292 fname = fname[:-7]293 else:294 err = (295 f"Invalid .stc filename: {fname!r}; needs to end with "296 "hemisphere tag ('...-lh.stc' or '...-rh.stc')"297 )298 raise OSError(err)299 elif fname.endswith(".w"):300 ftype = "w"301 if fname.endswith(("-lh.w", "-rh.w")):302 fname = fname[:-5]303 else:304 err = (305 f"Invalid .w filename: {fname!r}; needs to end with "306 "hemisphere tag ('...-lh.w' or '...-rh.w')"307 )308 raise OSError(err)309 elif fname.endswith(".h5"):310 ftype = "h5"311 fname = fname[:-3]312 else:313 raise RuntimeError(f"Unknown extension for file {fname_arg}")314 315 if ftype != "volume":316 stc_exist = [op.exists(f) for f in [fname + "-rh.stc", fname + "-lh.stc"]]317 w_exist = [op.exists(f) for f in [fname + "-rh.w", fname + "-lh.w"]]318 if all(stc_exist) and ftype != "w":319 ftype = "surface"320 elif all(w_exist):321 ftype = "w"322 elif op.exists(fname + ".h5"):323 ftype = "h5"324 elif op.exists(fname + "-stc.h5"):325 ftype = "h5"326 fname += "-stc"327 elif any(stc_exist) or any(w_exist):328 raise OSError(f"Hemisphere missing for {fname_arg!r}")329 else:330 raise OSError(f"SourceEstimate File(s) not found for: {fname_arg!r}")331 332 # read the files333 if ftype == "volume": # volume source space334 if fname.endswith(".stc"):335 kwargs = _read_stc(fname)336 elif fname.endswith(".w"):337 kwargs = _read_w(fname)338 kwargs["data"] = kwargs["data"][:, np.newaxis]339 kwargs["tmin"] = 0.0340 kwargs["tstep"] = 0.0341 else:342 raise OSError("Volume source estimate must end with .stc or .w")343 kwargs["vertices"] = [kwargs["vertices"]]344 elif ftype == "surface": # stc file with surface source spaces345 lh = _read_stc(fname + "-lh.stc")346 rh = _read_stc(fname + "-rh.stc")347 assert lh["tmin"] == rh["tmin"]348 assert lh["tstep"] == rh["tstep"]349 kwargs = lh.copy()350 kwargs["data"] = np.r_[lh["data"], rh["data"]]351 kwargs["vertices"] = [lh["vertices"], rh["vertices"]]352 elif ftype == "w": # w file with surface source spaces353 lh = _read_w(fname + "-lh.w")354 rh = _read_w(fname + "-rh.w")355 kwargs = lh.copy()356 kwargs["data"] = np.atleast_2d(np.r_[lh["data"], rh["data"]]).T357 kwargs["vertices"] = [lh["vertices"], rh["vertices"]]358 # w files only have a single time point359 kwargs["tmin"] = 0.0360 kwargs["tstep"] = 1.0361 ftype = "surface"362 elif ftype == "h5":363 read_hdf5, _ = _import_h5io_funcs()364 kwargs = read_hdf5(fname + ".h5", title="mnepython")365 ftype = kwargs.pop("src_type", "surface")366 if isinstance(kwargs["vertices"], np.ndarray):367 kwargs["vertices"] = [kwargs["vertices"]]368 369 if ftype != "volume":370 # Make sure the vertices are ordered371 vertices = kwargs["vertices"]372 if any(np.any(np.diff(v.astype(int)) <= 0) for v in vertices):373 sidx = [np.argsort(verts) for verts in vertices]374 vertices = [verts[idx] for verts, idx in zip(vertices, sidx)]375 data = kwargs["data"][np.r_[sidx[0], len(sidx[0]) + sidx[1]]]376 kwargs["vertices"] = vertices377 kwargs["data"] = data378 379 if "subject" not in kwargs:380 kwargs["subject"] = subject381 if subject is not None and subject != kwargs["subject"]:382 raise RuntimeError(383 f'provided subject name "{subject}" does not match '384 f'subject name from the file "{kwargs["subject"]}'385 )386 387 if ftype in ("volume", "discrete"):388 klass = VolVectorSourceEstimate389 elif ftype == "mixed":390 klass = MixedVectorSourceEstimate391 else:392 assert ftype == "surface"393 klass = VectorSourceEstimate394 if kwargs["data"].ndim < 3:395 klass = klass._scalar_class396 return klass(**kwargs)397 398 399def _get_src_type(src, vertices, warn_text=None):400 src_type = None401 if src is None:402 if warn_text is None:403 warn("src should not be None for a robust guess of stc type.")404 else:405 warn(warn_text)406 if isinstance(vertices, list) and len(vertices) == 2:407 src_type = "surface"408 elif (409 isinstance(vertices, np.ndarray)410 or isinstance(vertices, list)411 and len(vertices) == 1412 ):413 src_type = "volume"414 elif isinstance(vertices, list) and len(vertices) > 2:415 src_type = "mixed"416 else:417 src_type = src.kind418 assert src_type in ("surface", "volume", "mixed", "discrete")419 return src_type420 421 422def _make_stc(423 data,424 vertices,425 src_type=None,426 tmin=None,427 tstep=None,428 subject=None,429 vector=False,430 source_nn=None,431 warn_text=None,432):433 """Generate a surface, vector-surface, volume or mixed source estimate."""434 435 def guess_src_type():436 return _get_src_type(src=None, vertices=vertices, warn_text=warn_text)437 438 src_type = guess_src_type() if src_type is None else src_type439 440 if vector and src_type == "surface" and source_nn is None:441 raise RuntimeError("No source vectors supplied.")442 443 # infer Klass from src_type444 if src_type == "surface":445 Klass = VectorSourceEstimate if vector else SourceEstimate446 elif src_type in ("volume", "discrete"):447 Klass = VolVectorSourceEstimate if vector else VolSourceEstimate448 elif src_type == "mixed":449 Klass = MixedVectorSourceEstimate if vector else MixedSourceEstimate450 else:451 raise ValueError(452 "vertices has to be either a list with one or more arrays or an array"453 )454 455 # Rotate back for vector source estimates456 if vector:457 n_vertices = sum(len(v) for v in vertices)458 assert data.shape[0] in (n_vertices, n_vertices * 3)459 if len(data) == n_vertices:460 assert src_type == "surface" # should only be possible for this461 assert source_nn.shape == (n_vertices, 3)462 data = data[:, np.newaxis] * source_nn[:, :, np.newaxis]463 else:464 data = data.reshape((-1, 3, data.shape[-1]))465 assert source_nn.shape in ((n_vertices, 3, 3), (n_vertices * 3, 3))466 # This will be an identity transform for volumes, but let's keep467 # the code simple and general and just do the matrix mult468 data = np.matmul(469 np.transpose(source_nn.reshape(n_vertices, 3, 3), axes=[0, 2, 1]), data470 )471 472 return Klass(data=data, vertices=vertices, tmin=tmin, tstep=tstep, subject=subject)473 474 475def _verify_source_estimate_compat(a, b):476 """Make sure two SourceEstimates are compatible for arith. operations."""477 compat = False478 if type(a) is not type(b):479 raise ValueError(f"Cannot combine {type(a)} and {type(b)}.")480 if len(a.vertices) == len(b.vertices):481 if all(np.array_equal(av, vv) for av, vv in zip(a.vertices, b.vertices)):482 compat = True483 if not compat:484 raise ValueError(485 "Cannot combine source estimates that do not have "486 "the same vertices. Consider using stc.expand()."487 )488 if a.subject != b.subject:489 raise ValueError(490 "source estimates do not have the same subject "491 f"names, {repr(a.subject)} and {repr(b.subject)}"492 )493 494 495class _BaseSourceEstimate(TimeMixin, FilterMixin):496 _data_ndim = 2497 498 @verbose499 def __init__(self, data, vertices, tmin, tstep, subject=None, verbose=None):500 assert hasattr(self, "_data_ndim"), self.__class__.__name__501 assert hasattr(self, "_src_type"), self.__class__.__name__502 assert hasattr(self, "_src_count"), self.__class__.__name__503 kernel, sens_data = None, None504 if isinstance(data, tuple):505 if len(data) != 2:506 raise ValueError("If data is a tuple it has to be length 2")507 kernel, sens_data = data508 data = None509 if kernel.shape[1] != sens_data.shape[0]:510 raise ValueError(511 f"kernel ({kernel.shape}) and sens_data ({sens_data.shape}) "512 "have invalid dimensions"513 )514 if sens_data.ndim != 2:515 raise ValueError(516 "The sensor data must have 2 dimensions, got {sens_data.ndim}"517 )518 519 _validate_type(vertices, list, "vertices")520 if self._src_count is not None:521 if len(vertices) != self._src_count:522 raise ValueError(523 f"vertices must be a list with {self._src_count} entries, "524 f"got {len(vertices)}."525 )526 vertices = [np.array(v, np.int64) for v in vertices] # makes copy527 if any(np.any(np.diff(v) <= 0) for v in vertices):528 raise ValueError("Vertices must be ordered in increasing order.")529 530 n_src = sum([len(v) for v in vertices])531 532 # safeguard the user against doing something silly533 if data is not None:534 if data.ndim not in (self._data_ndim, self._data_ndim - 1):535 raise ValueError(536 f"Data (shape {data.shape}) must have {self._data_ndim} "537 f"dimensions for {self.__class__.__name__}"538 )539 if data.shape[0] != n_src:540 raise ValueError(541 f"Number of vertices ({n_src}) and stc.data.shape[0] "542 f"({data.shape[0]}) must match"543 )544 if self._data_ndim == 3:545 if data.shape[1] != 3:546 raise ValueError(547 "Data for VectorSourceEstimate must have "548 f"shape[1] == 3, got shape {data.shape}"549 )550 if data.ndim == self._data_ndim - 1: # allow upbroadcasting551 data = data[..., np.newaxis]552 553 self._data = data554 self._tmin = tmin555 self._tstep = tstep556 self.vertices = vertices557 self._kernel = kernel558 self._sens_data = sens_data559 self._kernel_removed = False560 self._times = None561 self._update_times()562 self.subject = _check_subject(None, subject, raise_error=False)563 564 def __repr__(self): # noqa: D105565 s = f"{sum(len(v) for v in self.vertices)} vertices"566 if self.subject is not None:567 s += f", subject : {self.subject}"568 s += ", tmin : %s (ms)" % (1e3 * self.tmin)569 s += ", tmax : %s (ms)" % (1e3 * self.times[-1])570 s += ", tstep : %s (ms)" % (1e3 * self.tstep)571 s += f", data shape : {self.shape}"572 sz = sum(object_size(x) for x in (self.vertices + [self.data]))573 s += f", ~{sizeof_fmt(sz)}"574 return f"<{type(self).__name__} | {s}>"575 576 @fill_doc577 def get_peak(578 self, tmin=None, tmax=None, mode="abs", vert_as_index=False, time_as_index=False579 ):580 """Get location and latency of peak amplitude.581 582 Parameters583 ----------584 %(get_peak_parameters)s585 586 Returns587 -------588 pos : int589 The vertex exhibiting the maximum response, either ID or index.590 latency : float591 The latency in seconds.592 """593 stc = self.magnitude() if self._data_ndim == 3 else self594 if self._n_vertices == 0:595 raise RuntimeError("Cannot find peaks with no vertices")596 vert_idx, time_idx, _ = _get_peak(stc.data, self.times, tmin, tmax, mode)597 if not vert_as_index:598 vert_idx = np.concatenate(self.vertices)[vert_idx]599 if not time_as_index:600 time_idx = self.times[time_idx]601 return vert_idx, time_idx602 603 @verbose604 def extract_label_time_course(605 self, labels, src, mode="auto", allow_empty=False, verbose=None606 ):607 """Extract label time courses for lists of labels.608 609 This function will extract one time course for each label. The way the610 time courses are extracted depends on the mode parameter.611 612 Parameters613 ----------614 %(labels_eltc)s615 %(src_eltc)s616 %(mode_eltc)s617 %(allow_empty_eltc)s618 %(verbose)s619 620 Returns621 -------622 %(label_tc_el_returns)s623 624 See Also625 --------626 extract_label_time_course : Extract time courses for multiple STCs.627 628 Notes629 -----630 %(eltc_mode_notes)s631 """632 return extract_label_time_course(633 self,634 labels,635 src,636 mode=mode,637 return_generator=False,638 allow_empty=allow_empty,639 verbose=verbose,640 )641 642 @verbose643 def apply_function(644 self, fun, picks=None, dtype=None, n_jobs=None, verbose=None, **kwargs645 ):646 """Apply a function to a subset of vertices.647 648 %(applyfun_summary_stc)s649 650 Parameters651 ----------652 %(fun_applyfun_stc)s653 %(picks_all)s654 %(dtype_applyfun)s655 %(n_jobs)s Ignored if ``vertice_wise=False`` as the workload656 is split across vertices.657 %(verbose)s658 %(kwargs_fun)s659 660 Returns661 -------662 self : instance of SourceEstimate663 The SourceEstimate object with transformed data.664 """665 _check_preload(self, "source_estimate.apply_function")666 picks = _picks_to_idx(len(self._data), picks, exclude=(), with_ref_meg=False)667 668 if not callable(fun):669 raise ValueError("fun needs to be a function")670 671 data_in = self._data672 if dtype is not None and dtype != self._data.dtype:673 self._data = self._data.astype(dtype)674 675 # check the dimension of the source estimate data676 _check_option("source_estimate.ndim", self._data.ndim, [2, 3])677 678 parallel, p_fun, n_jobs = parallel_func(_check_fun, n_jobs)679 if n_jobs == 1:680 # modify data inplace to save memory681 for idx in picks:682 self._data[idx, :] = _check_fun(fun, data_in[idx, :], **kwargs)683 else:684 # use parallel function685 data_picks_new = parallel(686 p_fun(fun, data_in[p, :], **kwargs) for p in picks687 )688 for pp, p in enumerate(picks):689 self._data[p, :] = data_picks_new[pp]690 691 return self692 693 @verbose694 def apply_baseline(self, baseline=(None, 0), *, verbose=None):695 """Baseline correct source estimate data.696 697 Parameters698 ----------699 %(baseline_stc)s700 Defaults to ``(None, 0)``, i.e. beginning of the the data until701 time point zero.702 %(verbose)s703 704 Returns705 -------706 stc : instance of SourceEstimate707 The baseline-corrected source estimate object.708 709 Notes710 -----711 Baseline correction can be done multiple times.712 """713 self.data = rescale(self.data, self.times, baseline, copy=False)714 return self715 716 @verbose717 def save(self, fname, ftype="h5", *, overwrite=False, verbose=None):718 """Save the full source estimate to an HDF5 file.719 720 Parameters721 ----------722 fname : path-like723 The file name to write the source estimate to, should end in724 ``'-stc.h5'``.725 ftype : str726 File format to use. Currently, the only allowed values is ``"h5"``.727 %(overwrite)s728 729 .. versionadded:: 1.0730 %(verbose)s731 """732 fname = _check_fname(fname=fname, overwrite=True) # check below733 if ftype != "h5":734 raise ValueError(735 f"{self.__class__.__name__} objects can only be written as HDF5 files."736 )737 _, write_hdf5 = _import_h5io_funcs()738 if fname.suffix != ".h5":739 fname = fname.with_name(f"{fname.name}-stc.h5")740 fname = _check_fname(fname=fname, overwrite=overwrite)741 write_hdf5(742 fname,743 dict(744 vertices=self.vertices,745 data=self.data,746 tmin=self.tmin,747 tstep=self.tstep,748 subject=self.subject,749 src_type=self._src_type,750 ),751 title="mnepython",752 overwrite=True,753 )754 755 @copy_function_doc_to_method_doc(plot_source_estimates)756 def plot(757 self,758 subject=None,759 surface="inflated",760 hemi="lh",761 colormap="auto",762 time_label="auto",763 smoothing_steps=10,764 transparent=True,765 alpha=1.0,766 time_viewer="auto",767 *,768 subjects_dir=None,769 figure=None,770 views="auto",771 colorbar=True,772 clim="auto",773 cortex="classic",774 size=800,775 background="black",776 foreground=None,777 initial_time=None,778 time_unit="s",779 backend="auto",780 spacing="oct6",781 title=None,782 show_traces="auto",783 src=None,784 volume_options=1.0,785 view_layout="vertical",786 add_data_kwargs=None,787 brain_kwargs=None,788 verbose=None,789 ):790 brain = plot_source_estimates(791 self,792 subject,793 surface=surface,794 hemi=hemi,795 colormap=colormap,796 time_label=time_label,797 smoothing_steps=smoothing_steps,798 transparent=transparent,799 alpha=alpha,800 time_viewer=time_viewer,801 subjects_dir=subjects_dir,802 figure=figure,803 views=views,804 colorbar=colorbar,805 clim=clim,806 cortex=cortex,807 size=size,808 background=background,809 foreground=foreground,810 initial_time=initial_time,811 time_unit=time_unit,812 backend=backend,813 spacing=spacing,814 title=title,815 show_traces=show_traces,816 src=src,817 volume_options=volume_options,818 view_layout=view_layout,819 add_data_kwargs=add_data_kwargs,820 brain_kwargs=brain_kwargs,821 verbose=verbose,822 )823 return brain824 825 @property826 def sfreq(self):827 """Sample rate of the data."""828 return 1.0 / self.tstep829 830 @property831 def _n_vertices(self):832 return sum(len(v) for v in self.vertices)833 834 def _remove_kernel_sens_data_(self):835 """Remove kernel and sensor space data and compute self._data."""836 if self._kernel is not None or self._sens_data is not None:837 self._kernel_removed = True838 self._data = np.dot(self._kernel, self._sens_data)839 self._kernel = None840 self._sens_data = None841 842 @fill_doc843 def crop(self, tmin=None, tmax=None, include_tmax=True):844 """Restrict SourceEstimate to a time interval.845 846 Parameters847 ----------848 tmin : float | None849 The first time point in seconds. If None the first present is used.850 tmax : float | None851 The last time point in seconds. If None the last present is used.852 %(include_tmax)s853 854 Returns855 -------856 stc : instance of SourceEstimate857 The cropped source estimate.858 """859 mask = _time_mask(860 self.times, tmin, tmax, sfreq=self.sfreq, include_tmax=include_tmax861 )862 self.tmin = self.times[np.where(mask)[0][0]]863 if self._kernel is not None and self._sens_data is not None:864 self._sens_data = self._sens_data[..., mask]865 else:866 self.data = self.data[..., mask]867 868 return self # return self for chaining methods869 870 @verbose871 def resample(872 self,873 sfreq,874 *,875 npad=100,876 method="fft",877 window="auto",878 pad="auto",879 n_jobs=None,880 verbose=None,881 ):882 """Resample data.883 884 If appropriate, an anti-aliasing filter is applied before resampling.885 See :ref:`resampling-and-decimating` for more information.886 887 Parameters888 ----------889 sfreq : float890 New sample rate to use.891 npad : int | str892 Amount to pad the start and end of the data.893 Can also be "auto" to use a padding that will result in894 a power-of-two size (can be much faster).895 %(method_resample)s896 897 .. versionadded:: 1.7898 %(window_resample)s899 900 .. versionadded:: 1.7901 %(pad_resample_auto)s902 903 .. versionadded:: 1.7904 %(n_jobs)s905 %(verbose)s906 907 Returns908 -------909 stc : instance of SourceEstimate910 The resampled source estimate.911 912 Notes913 -----914 For some data, it may be more accurate to use npad=0 to reduce915 artifacts. This is dataset dependent -- check your data!916 917 Note that the sample rate of the original data is inferred from tstep.918 """919 from .filter import _check_resamp_noop920 921 o_sfreq = 1.0 / self.tstep922 if _check_resamp_noop(sfreq, o_sfreq):923 return self924 925 # resampling in sensor instead of source space gives a somewhat926 # different result, so we don't allow it927 self._remove_kernel_sens_data_()928 929 data = self.data930 if data.dtype == np.float32:931 data = data.astype(np.float64)932 self.data = resample(933 data, sfreq, o_sfreq, npad=npad, window=window, n_jobs=n_jobs, method=method934 )935 936 # adjust indirectly affected variables937 self.tstep = 1.0 / sfreq938 return self939 940 @property941 def data(self):942 """Numpy array of source estimate data."""943 if self._data is None:944 # compute the solution the first time the data is accessed and945 # remove the kernel and sensor data946 self._remove_kernel_sens_data_()947 return self._data948 949 @data.setter950 def data(self, value):951 value = np.asarray(value)952 if self._data is not None and value.ndim != self._data.ndim:953 raise ValueError(f"Data array should have {self._data.ndim} dimensions.")954 n_verts = sum(len(v) for v in self.vertices)955 if value.shape[0] != n_verts:956 raise ValueError(957 "The first dimension of the data array must match the number of "958 f"vertices ({value.shape[0]} != {n_verts})."959 )960 self._data = value961 self._update_times()962 963 @property964 def shape(self):965 """Shape of the data."""966 if self._data is not None:967 return self._data.shape968 return (self._kernel.shape[0], self._sens_data.shape[1])969 970 @property971 def tmin(self):972 """The first timestamp."""973 return self._tmin974 975 @tmin.setter976 def tmin(self, value):977 self._tmin = float(value)978 self._update_times()979 980 @property981 def tstep(self):982 """The change in time between two consecutive samples (1 / sfreq)."""983 return self._tstep984 985 @tstep.setter986 def tstep(self, value):987 if value <= 0:988 raise ValueError(".tstep must be greater than 0.")989 self._tstep = float(value)990 self._update_times()991 992 @property993 def times(self):994 """A timestamp for each sample."""995 return self._times996 997 @times.setter998 def times(self, value):999 raise ValueError(1000 "You cannot write to the .times attribute directly. "1001 "This property automatically updates whenever "1002 ".tmin, .tstep or .data changes."1003 )1004 1005 def _update_times(self):1006 """Update the times attribute after changing tmin, tmax, or tstep."""1007 self._times = self.tmin + (self.tstep * np.arange(self.shape[-1]))1008 self._times.flags.writeable = False1009 1010 def __add__(self, a):1011 """Add source estimates."""1012 stc = self.copy()1013 stc += a1014 return stc1015 1016 def __iadd__(self, a): # noqa: D1051017 self._remove_kernel_sens_data_()1018 if isinstance(a, _BaseSourceEstimate):1019 _verify_source_estimate_compat(self, a)1020 self.data += a.data1021 else:1022 self.data += a1023 return self1024 1025 def mean(self):1026 """Make a summary stc file with mean over time points.1027 1028 Returns1029 -------1030 stc : SourceEstimate | VectorSourceEstimate1031 The modified stc.1032 """1033 out = self.sum()1034 out /= len(self.times)1035 return out1036 1037 def sum(self):1038 """Make a summary stc file with sum over time points.1039 1040 Returns1041 -------1042 stc : SourceEstimate | VectorSourceEstimate1043 The modified stc.1044 """1045 data = self.data1046 tmax = self.tmin + self.tstep * data.shape[-1]1047 tmin = (self.tmin + tmax) / 2.01048 tstep = tmax - self.tmin1049 sum_stc = self.__class__(1050 self.data.sum(axis=-1, keepdims=True),1051 vertices=self.vertices,1052 tmin=tmin,1053 tstep=tstep,1054 subject=self.subject,1055 )1056 return sum_stc1057 1058 def __sub__(self, a):1059 """Subtract source estimates."""1060 stc = self.copy()1061 stc -= a1062 return stc1063 1064 def __isub__(self, a): # noqa: D1051065 self._remove_kernel_sens_data_()1066 if isinstance(a, _BaseSourceEstimate):1067 _verify_source_estimate_compat(self, a)1068 self.data -= a.data1069 else:1070 self.data -= a1071 return self1072 1073 def __truediv__(self, a): # noqa: D1051074 return self.__div__(a)1075 1076 def __div__(self, a): # noqa: D1051077 """Divide source estimates."""1078 stc = self.copy()1079 stc /= a1080 return stc1081 1082 def __itruediv__(self, a): # noqa: D1051083 return self.__idiv__(a)1084 1085 def __idiv__(self, a): # noqa: D1051086 self._remove_kernel_sens_data_()1087 if isinstance(a, _BaseSourceEstimate):1088 _verify_source_estimate_compat(self, a)1089 self.data /= a.data1090 else:1091 self.data /= a1092 return self1093 1094 def __mul__(self, a):1095 """Multiply source estimates."""1096 stc = self.copy()1097 stc *= a1098 return stc1099 1100 def __imul__(self, a): # noqa: D1051101 self._remove_kernel_sens_data_()1102 if isinstance(a, _BaseSourceEstimate):1103 _verify_source_estimate_compat(self, a)1104 self.data *= a.data1105 else:1106 self.data *= a1107 return self1108 1109 def __pow__(self, a): # noqa: D1051110 stc = self.copy()1111 stc **= a1112 return stc1113 1114 def __ipow__(self, a): # noqa: D1051115 self._remove_kernel_sens_data_()1116 self.data **= a1117 return self1118 1119 def __radd__(self, a): # noqa: D1051120 return self + a1121 1122 def __rsub__(self, a): # noqa: D1051123 return self - a1124 1125 def __rmul__(self, a): # noqa: D1051126 return self * a1127 1128 def __rdiv__(self, a): # noqa: D1051129 return self / a1130 1131 def __neg__(self): # noqa: D1051132 """Negate the source estimate."""1133 stc = self.copy()1134 stc._remove_kernel_sens_data_()1135 stc.data *= -11136 return stc1137 1138 def __pos__(self): # noqa: D1051139 return self1140 1141 def __abs__(self):1142 """Compute the absolute value of the data.1143 1144 Returns1145 -------1146 stc : instance of _BaseSourceEstimate1147 A version of the source estimate, where the data attribute is set1148 to abs(self.data).1149 """1150 stc = self.copy()1151 stc._remove_kernel_sens_data_()1152 stc._data = abs(stc._data)1153 return stc1154 1155 def sqrt(self):1156 """Take the square root.1157 1158 Returns1159 -------1160 stc : instance of SourceEstimate1161 A copy of the SourceEstimate with sqrt(data).1162 """1163 return self ** (0.5)1164 1165 def copy(self):1166 """Return copy of source estimate instance.1167 1168 Returns1169 -------1170 stc : instance of SourceEstimate1171 A copy of the source estimate.1172 """1173 return copy.deepcopy(self)1174 1175 def bin(self, width, tstart=None, tstop=None, func=np.mean):1176 """Return a source estimate object with data summarized over time bins.1177 1178 Time bins of ``width`` seconds. This method is intended for1179 visualization only. No filter is applied to the data before binning,1180 making the method inappropriate as a tool for downsampling data.1181 1182 Parameters1183 ----------1184 width : scalar1185 Width of the individual bins in seconds.1186 tstart : scalar | None1187 Time point where the first bin starts. The default is the first1188 time point of the stc.1189 tstop : scalar | None1190 Last possible time point contained in a bin (if the last bin would1191 be shorter than width it is dropped). The default is the last time1192 point of the stc.1193 func : callable1194 Function that is applied to summarize the data. Needs to accept a1195 numpy.array as first input and an ``axis`` keyword argument.1196 1197 Returns1198 -------1199 stc : SourceEstimate | VectorSourceEstimate1200 The binned source estimate.