diff --git a/.circleci/config.yml b/.circleci/config.yml index dace4cb7c61..96052d67bd3 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -113,115 +113,13 @@ jobs: python setup.py develop --user mkdir -p ~/mne_data touch pattern.txt; - if [ "$CIRCLE_BRANCH" == "master" ] || [[ $(cat gitlog.txt) == *"[circle full]"* ]]; then - echo "Doing a full dev build"; - echo html_dev > build.txt; - python -c "import mne; mne.datasets._download_all_example_data()"; - elif [ "$CIRCLE_BRANCH" == "maint/0.21" ]; then - echo "Doing a full stable build"; - echo html_stable > build.txt; - python -c "import mne; mne.datasets._download_all_example_data()"; - else - echo "Doing a partial build"; - if ! git remote -v | grep upstream ; then git remote add upstream git://github.com/mne-tools/mne-python.git; fi - git fetch upstream - FNAMES=$(git diff --name-only $(git merge-base $CIRCLE_BRANCH upstream/master) $CIRCLE_BRANCH); - if [[ $(cat gitlog.txt) == *"[circle front]"* ]]; then - FNAMES="tutorials/source-modeling/plot_mne_dspm_source_localization.py tutorials/machine-learning/plot_receptive_field.py examples/connectivity/plot_mne_inverse_label_connectivity.py tutorials/machine-learning/plot_sensors_decoding.py tutorials/stats-source-space/plot_stats_cluster_spatio_temporal.py tutorials/evoked/plot_20_visualize_evoked.py "${FNAMES}; - python -c "import mne; print(mne.datasets.testing.data_path(update_path=True))"; - fi; - echo FNAMES="$FNAMES"; - for FNAME in $FNAMES; do - if [[ `expr match $FNAME "\(tutorials\|examples\)/.*plot_.*\.py"` ]] ; then - echo "Checking example $FNAME ..."; - PATTERN=`basename $FNAME`"\\|"$PATTERN; - if [[ $(cat $FNAME | grep -x ".*datasets.*sample.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.sample.data_path(update_path=True))"; - fi; - if [[ $(cat $FNAME | grep -x ".*datasets.*fetch_fsaverage.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.fetch_fsaverage(verbose=True))"; - fi; - if [[ $(cat $FNAME | grep -x ".*datasets.*spm_face.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.spm_face.data_path(update_path=True))"; - fi; - if [[ $(cat $FNAME | grep -x ".*datasets.*somato.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.somato.data_path(update_path=True))"; - fi; - if [[ $(cat $FNAME | grep -x ".*datasets.*eegbci.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.eegbci.load_data(1, [3, 6, 10, 14], update_path=True))"; - python -c "import mne; print(mne.datasets.eegbci.load_data(2, [3], update_path=True))"; - python -c "import mne; print(mne.datasets.eegbci.load_data(3, [3], update_path=True))"; - python -c "import mne; print(mne.datasets.eegbci.load_data(4, [3], update_path=True))"; - fi; - if [[ $(cat $FNAME | grep -x ".*datasets.*sleep_physionet.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.sleep_physionet.age.fetch_data([0, 1], recording=[1], update_path=True))"; - fi; - if [[ $(cat $FNAME | grep -x ".*datasets.*hf_sef.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.hf_sef.data_path(update_path=True))"; - fi; - if [[ $(cat $FNAME | grep -x ".*brainstorm.*bst_auditory.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.brainstorm.bst_auditory.data_path(update_path=True))" --accept-brainstorm-license; - fi; - if [[ $(cat $FNAME | grep -x ".*brainstorm.*bst_resting.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.brainstorm.bst_resting.data_path(update_path=True))" --accept-brainstorm-license; - fi; - if [[ $(cat $FNAME | grep -x ".*brainstorm.*bst_raw.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.brainstorm.bst_raw.data_path(update_path=True))" --accept-brainstorm-license; - fi; - if [[ $(cat $FNAME | grep -x ".*brainstorm.*bst_phantom_ctf.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.brainstorm.bst_phantom_ctf.data_path(update_path=True))" --accept-brainstorm-license; - fi; - if [[ $(cat $FNAME | grep -x ".*brainstorm.*bst_phantom_elekta.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.brainstorm.bst_phantom_elekta.data_path(update_path=True))" --accept-brainstorm-license; - fi; - if [[ $(cat $FNAME | grep -x ".*datasets.*hcp_mmp_parcellation.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.sample.data_path(update_path=True))"; - python -c "import mne; print(mne.datasets.fetch_hcp_mmp_parcellation(subjects_dir=mne.datasets.sample.data_path() + '/subjects'))" --accept-hcpmmp-license; - fi; - if [[ $(cat $FNAME | grep -x ".*datasets.*misc.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.misc.data_path(update_path=True))"; - fi; - if [[ $(cat $FNAME | grep -x ".*datasets.*testing.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.testing.data_path(update_path=True))"; - fi; - if [[ $(cat $FNAME | grep -x ".*datasets.*kiloword.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.kiloword.data_path(update_path=True))"; - fi; - if [[ $(cat $FNAME | grep -x ".*datasets.*mtrf.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.mtrf.data_path(update_path=True))"; - fi; - if [[ $(cat $FNAME | grep -x ".*datasets.*fieldtrip_cmc.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.fieldtrip_cmc.data_path(update_path=True))"; - fi; - if [[ $(cat $FNAME | grep -x ".*datasets.*multimodal.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.multimodal.data_path(update_path=True))"; - fi; - if [[ $(cat $FNAME | grep -x ".*datasets.*fnirs_motor.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.fnirs_motor.data_path(update_path=True))"; - fi; - if [[ $(cat $FNAME | grep -x ".*datasets.*opm.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.opm.data_path(update_path=True))"; - fi; - if [[ $(cat $FNAME | grep -x ".*datasets.*phantom_4dbti.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.phantom_4dbti.data_path(update_path=True))"; - fi; - if [[ $(cat $FNAME | grep -x ".*datasets.*limo.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.limo.data_path(subject=1, update_path=True))"; - fi; - if [[ $(cat $FNAME | grep -x ".*datasets.*refmeg_noise.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.refmeg_noise.data_path(update_path=True))"; - fi; - fi; - done; - echo PATTERN="$PATTERN"; - if [[ $PATTERN ]]; then - PATTERN="\(${PATTERN::-2}\)"; - echo html_dev-pattern > build.txt; - else - echo html_dev-noplot > build.txt; - fi; - fi; - echo "$PATTERN" > pattern.txt; + ./tools/circleci_download.sh + + - run: + name: Get data (again) + when: on_fail + command: | + ./tools/circleci_download.sh - run: name: Verify build type @@ -234,7 +132,7 @@ jobs: - run: name: make test-doc command: | - if [[ $(cat gitlog.txt) == *"[circle front]"* ]] || [[ $(cat build.txt) == "html_dev" ]] || [[ $(cat build.txt) == "html_stable" ]]; then + if [[ $(cat gitlog.txt) == *"[circle front]"* ]] || [[ $(cat build.txt) == "html_dev-memory" ]] || [[ $(cat build.txt) == "html_stable-memory" ]]; then make test-doc; mkdir -p doc/_build/test-results/test-doc; cp junit-results.xml doc/_build/test-results/test-doc/junit.xml; @@ -245,6 +143,16 @@ jobs: command: | cd doc; PATTERN=$(cat ../pattern.txt) make $(cat ../build.txt); + - run: + name: Show profiling output + when: always + command: | + if compgen -G "doc/*.dat" > /dev/null; then + mkdir -p doc/generated + mprof plot doc/*.dat --output doc/generated/memory.png + else + echo "No profile data found in doc/" + fi - run: name: Sanity check system state command: | @@ -254,7 +162,7 @@ jobs: - run: name: Reduce artifact upload time command: | - if grep -q html_dev-pattern build.txt || grep -q html_dev-noplot build.txt; then + if grep -q html_dev-pattern-memory build.txt || grep -q html_dev-noplot build.txt; then tar czf doc/_build/html/_downloads.tgz doc/_build/html/_downloads rm -Rf doc/_build/html/_downloads rm -f doc/auto_*/*/*.pickle diff --git a/.gitignore b/.gitignore index 8f00a1f3864..dcb21f31197 100644 --- a/.gitignore +++ b/.gitignore @@ -61,6 +61,8 @@ pip-log.txt tags doc/coverages doc/samples +doc/*.dat +doc/fil-result cover *.html diff --git a/doc/Makefile b/doc/Makefile index dbe11d7fdeb..34e315f1e73 100644 --- a/doc/Makefile +++ b/doc/Makefile @@ -5,6 +5,7 @@ SPHINXOPTS = -nWT --keep-going SPHINXBUILD = sphinx-build PAPER = +MPROF = SG_STAMP_STARTS=true mprof run --python sphinx # Internal variables. PAPEROPT_a4 = -D latex_paper_size=a4 @@ -40,13 +41,28 @@ html_stable: @echo @echo "Build finished. The HTML pages are in _build/html_stable." +html_stable-memory: + $(MPROF) -b html $(ALLSPHINXOPTS) _build/html_stable + @echo + @echo "Build finished. The HTML pages are in _build/html_stable." + html_dev: BUILD_DEV_HTML=1 $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) _build/html @echo @echo "Build finished. The HTML pages are in _build/html" +html_dev-memory: + BUILD_DEV_HTML=1 $(MPROF) -b html $(ALLSPHINXOPTS) _build/html + @echo + @echo "Build finished. The HTML pages are in _build/html" + html_dev-pattern: - BUILD_DEV_HTML=1 $(SPHINXBUILD) -D plot_gallery=1 -D sphinx_gallery_conf.filename_pattern=$(PATTERN) -b html $(ALLSPHINXOPTS) _build/html + BUILD_DEV_HTML=1 $(SPHINXBUILD) -D sphinx_gallery_conf.filename_pattern=$(PATTERN) -D sphinx_gallery_conf.run_stale_examples=True -b html $(ALLSPHINXOPTS) _build/html + @echo + @echo "Build finished. The HTML pages are in _build/html" + +html_dev-pattern-memory: + BUILD_DEV_HTML=1 $(MPROF) -D sphinx_gallery_conf.filename_pattern=$(PATTERN) -D sphinx_gallery_conf.run_stale_examples=True -b html $(ALLSPHINXOPTS) _build/html @echo @echo "Build finished. The HTML pages are in _build/html" diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc index 8c0e6f51696..c563dc23eaa 100644 --- a/doc/changes/latest.inc +++ b/doc/changes/latest.inc @@ -23,6 +23,8 @@ Enhancements - Add ``proj`` argument to :func:`mne.make_fixed_length_epochs` by `Eric Larson`_ (:gh:`8351`) +- Reduce memory usage of volume source spaces by `Eric Larson`_ (:gh:`8379`) + - Speed up :class:`mne.decoding.TimeDelayingRidge` with edge correction using Numba by `Eric Larson`_ (:gh:`8323`) Bugs @@ -65,7 +67,7 @@ Bugs - Pass ``rank`` everyhwere in forward preparation for source imaging. This bug affected sparse solvers when using maxfilter data, by `Alex Gramfort`_ (:gh:`8368`) -- Fix bug in :func:`mne.viz.plot_alignment` where ECoG and sEEG channels were not plotted and fNIRS channels were always plotted in the head coordinate frame by `Eric Larson`_ +- Fix bug in :func:`mne.viz.plot_alignment` where ECoG and sEEG channels were not plotted and fNIRS channels were always plotted in the head coordinate frame by `Eric Larson`_ (:gh:`8393`) API changes ~~~~~~~~~~~ @@ -79,3 +81,5 @@ API changes - The ``n_pca_components`` argument of :class:`~mne.preprocessing.ICA` has been deprecated, use ``n_pca_components`` in :meth:`~mne.preprocessing.ICA.apply` by `Eric Larson`_ (:gh:`8356`) - The ``trans`` argument of :func:`mne.extract_label_time_course` is deprecated and will be removed in 0.23 as it is no longer necessary by `Eric Larson`_ + +- Update the ``backend`` parameter of :func:`mne.viz.plot_source_estimates` to integrate ``pyvista`` by `Guillaume Favelier`_ diff --git a/doc/conf.py b/doc/conf.py index 874536e3762..def92d114ca 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -14,10 +14,10 @@ from datetime import date from distutils.version import LooseVersion -import gc import os import os.path as op import sys +import time import warnings import sphinx_gallery @@ -25,7 +25,9 @@ from numpydoc import docscrape import matplotlib import mne -from mne.utils import linkcode_resolve # noqa, analysis:ignore +from mne.viz import Brain +from mne.utils import (linkcode_resolve, # noqa, analysis:ignore + _assert_no_instances, sizeof_fmt) if LooseVersion(sphinx_gallery.__version__) < LooseVersion('0.2'): raise ImportError('Must have at least version 0.2 of sphinx-gallery, got ' @@ -370,8 +372,8 @@ ############################################################################## # sphinx-gallery -examples_dirs = ['../examples', '../tutorials'] -gallery_dirs = ['auto_examples', 'auto_tutorials'] +examples_dirs = ['../tutorials', '../examples'] +gallery_dirs = ['auto_tutorials', 'auto_examples'] os.environ['_MNE_BUILDING_DOC'] = 'true' scrapers = ('matplotlib',) @@ -441,17 +443,33 @@ def setup(app): class Resetter(object): """Simple class to make the str(obj) static for Sphinx build env hash.""" + def __init__(self): + self.t0 = time.time() + def __repr__(self): return '<%s>' % (self.__class__.__name__,) def __call__(self, gallery_conf, fname): import matplotlib.pyplot as plt + try: + from pyvista import Plotter + except ImportError: + Plotter = None reset_warnings(gallery_conf, fname) # in case users have interactive mode turned on in matplotlibrc, # turn it off here (otherwise the build can be very slow) plt.ioff() - gc.collect() plt.rcParams['animation.embed_limit'] = 30. + _assert_no_instances(Brain, 'running') # calls gc.collect() + if Plotter is not None: + _assert_no_instances(Plotter, 'running') + # This will overwrite some Sphinx printing but it's useful + # for memory timestamps + if os.getenv('SG_STAMP_STARTS', '').lower() == 'true': + import psutil + process = psutil.Process(os.getpid()) + mem = sizeof_fmt(process.memory_info().rss) + print(f'{time.time() - self.t0:6.1f} s : {mem}'.ljust(22)) def reset_warnings(gallery_conf, fname): diff --git a/examples/connectivity/plot_mixed_source_space_connectivity.py b/examples/connectivity/plot_mixed_source_space_connectivity.py index b1be975130e..c3a724de57b 100644 --- a/examples/connectivity/plot_mixed_source_space_connectivity.py +++ b/examples/connectivity/plot_mixed_source_space_connectivity.py @@ -179,7 +179,7 @@ # # By default matplotlib does not save using the facecolor, even though this was # set when the figure was generated. If not set via savefig, the labels, title, -# and legend will be cut off from the output png file. - -# fname_fig = data_path + '/MEG/sample/plot_mixed_connect.png' -# plt.savefig(fname_fig, facecolor='black') +# and legend will be cut off from the output png file:: +# +# >>> fname_fig = data_path + '/MEG/sample/plot_mixed_connect.png' +# >>> plt.savefig(fname_fig, facecolor='black') diff --git a/examples/decoding/plot_receptive_field_mtrf.py b/examples/decoding/plot_receptive_field_mtrf.py index af29e51bd7c..a45ca91b3c8 100644 --- a/examples/decoding/plot_receptive_field_mtrf.py +++ b/examples/decoding/plot_receptive_field_mtrf.py @@ -150,7 +150,6 @@ ax.set(title="Topomap of model coefficients\nfor delay %s" % time_plot) mne.viz.tight_layout() - ############################################################################### # Create and fit a stimulus reconstruction model # ---------------------------------------------- @@ -259,8 +258,6 @@ % (time_plot[0], time_plot[1])) mne.viz.tight_layout() -plt.show() - ############################################################################### # References # ---------- diff --git a/examples/inverse/plot_evoked_ers_source_power.py b/examples/inverse/plot_evoked_ers_source_power.py index a92bdab2252..2e1346c7768 100644 --- a/examples/inverse/plot_evoked_ers_source_power.py +++ b/examples/inverse/plot_evoked_ers_source_power.py @@ -33,7 +33,8 @@ raw_fname = op.join(data_path, 'sub-{}'.format(subject), 'meg', 'sub-{}_task-{}_meg.fif'.format(subject, task)) -raw = mne.io.read_raw_fif(raw_fname) +# crop to 5 minutes to save memory +raw = mne.io.read_raw_fif(raw_fname).crop(0, 300) # We are interested in the beta band (12-30 Hz) raw.load_data().filter(12, 30) @@ -45,7 +46,7 @@ # Read epochs events = mne.find_events(raw) epochs = mne.Epochs(raw, events, event_id=1, tmin=-1.5, tmax=2, picks=picks, - preload=True) + preload=True, decim=3) # Read forward operator and point to freesurfer subject directory fname_fwd = op.join(data_path, 'derivatives', 'sub-{}'.format(subject), @@ -86,7 +87,8 @@ def _gen_dics(active_win, baseline_win, epochs): tmax=baseline_win[1], decim=20) csd_ers = csd_morlet(epochs, freqs, tmin=active_win[0], tmax=active_win[1], decim=20) - filters = make_dics(epochs.info, fwd, csd.mean(), pick_ori='max-power') + filters = make_dics(epochs.info, fwd, csd.mean(), pick_ori='max-power', + reduce_rank=True) stc_base, freqs = apply_dics_csd(csd_baseline.mean(), filters) stc_act, freqs = apply_dics_csd(csd_ers.mean(), filters) stc_act /= stc_base @@ -122,11 +124,22 @@ def _gen_mne(active_cov, baseline_cov, common_cov, fwd, info, method='dSPM'): ############################################################################### # Plot source estimates # --------------------- +# DICS: -brains = list() -for method, stc in zip(['DICS', 'LCMV', 'dSPM'], - [stc_dics, stc_lcmv, stc_dspm]): - title = '%s source power in the 12-30 Hz frequency band' % method - brains.append(stc.plot( - hemi='rh', subjects_dir=subjects_dir, subject=subject, - time_label=title)) +brain_dics = stc_dics.plot( + hemi='rh', subjects_dir=subjects_dir, subject=subject, + time_label='DICS source power in the 12-30 Hz frequency band') + +############################################################################### +# LCMV: + +brain_lcmv = stc_lcmv.plot( + hemi='rh', subjects_dir=subjects_dir, subject=subject, + time_label='LCMV source power in the 12-30 Hz frequency band') + +############################################################################### +# dSPM: + +brain_dspm = stc_dspm.plot( + hemi='rh', subjects_dir=subjects_dir, subject=subject, + time_label='dSPM source power in the 12-30 Hz frequency band') diff --git a/examples/inverse/plot_mixed_source_space_inverse.py b/examples/inverse/plot_mixed_source_space_inverse.py index e45f7ade2f4..f9d7a6a56ce 100644 --- a/examples/inverse/plot_mixed_source_space_inverse.py +++ b/examples/inverse/plot_mixed_source_space_inverse.py @@ -1,7 +1,7 @@ """ -=================================================================== -Compute MNE inverse solution on evoked data in a mixed source space -=================================================================== +===================================================================== +Compute MNE inverse solution on evoked data with a mixed source space +===================================================================== Create a mixed source space and compute an MNE inverse solution on an evoked dataset. @@ -71,16 +71,16 @@ # Generate the mixed source space src += vol_src +print(f"The source space contains {len(src)} spaces and " + f"{sum(s['nuse'] for s in src)} vertices") -# Visualize the source space. -src.plot(subjects_dir=subjects_dir) +############################################################################### +# View the source space +# --------------------- -n = sum(src[i]['nuse'] for i in range(len(src))) -print('the src space contains %d spaces and %d points' % (len(src), n)) +src.plot(subjects_dir=subjects_dir) ############################################################################### -# Viewing the source space -# ------------------------ # We could write the mixed source space with:: # # >>> write_source_spaces(fname_mixed_src, src, overwrite=True) @@ -98,13 +98,12 @@ fname_evoked, fname_trans, src, fname_bem, mindist=5.0, # ignore sources<=5mm from innerskull meg=True, eeg=False, n_jobs=1) +del src # save memory leadfield = fwd['sol']['data'] print("Leadfield size : %d sensors x %d dipoles" % leadfield.shape) - -src_fwd = fwd['src'] -n = sum(src_fwd[i]['nuse'] for i in range(len(src_fwd))) -print('the fwd src space contains %d spaces and %d points' % (len(src_fwd), n)) +print(f"The fwd source space contains {len(fwd['src'])} spaces and " + f"{sum(s['nuse'] for s in fwd['src'])} vertices") # Load data condition = 'Left Auditory' @@ -124,6 +123,7 @@ inverse_operator = make_inverse_operator( evoked.info, fwd, noise_cov, depth=None, loose=loose, verbose=True) +del fwd stc = apply_inverse(evoked, inverse_operator, lambda2, inv_method, pick_ori=None) @@ -169,7 +169,7 @@ # plot the times series of 2 labels fig, axes = plt.subplots(1) axes.plot(1e3 * stc.times, label_ts[0][0, :], 'k', label='bankssts-lh') -axes.plot(1e3 * stc.times, label_ts[0][71, :].T, 'r', label='Brain-stem') +axes.plot(1e3 * stc.times, label_ts[0][-1, :].T, 'r', label='Brain-stem') axes.set(xlabel='Time (ms)', ylabel='MNE current (nAm)') axes.legend() mne.viz.tight_layout() diff --git a/examples/inverse/plot_psf_ctf_vertices.py b/examples/inverse/plot_psf_ctf_vertices.py index 02e78e01260..26db8d8e8b3 100644 --- a/examples/inverse/plot_psf_ctf_vertices.py +++ b/examples/inverse/plot_psf_ctf_vertices.py @@ -26,8 +26,8 @@ # read forward solution forward = mne.read_forward_solution(fname_fwd) # forward operator with fixed source orientations -forward = mne.convert_forward_solution(forward, surf_ori=True, - force_fixed=True) +mne.convert_forward_solution(forward, surf_ori=True, + force_fixed=True, copy=False) # noise covariance matrix noise_cov = mne.read_cov(fname_cov) @@ -58,8 +58,9 @@ stc_ctf = get_cross_talk(rm_lor, forward['src'], sources, norm=True) ############################################################################## -# Visualise +# Visualize # --------- +# PSF: # Which vertex corresponds to selected source vertno_lh = forward['src'][0]['vertno'] @@ -69,8 +70,7 @@ vert_max_psf = vertno_lh[stc_psf.data.argmax()] vert_max_ctf = vertno_lh[stc_ctf.data.argmax()] -brain_psf = stc_psf.plot('sample', 'inflated', 'lh', subjects_dir=subjects_dir, - figure=1) +brain_psf = stc_psf.plot('sample', 'inflated', 'lh', subjects_dir=subjects_dir) brain_psf.show_view('ventral') brain_psf.add_text(0.1, 0.9, 'sLORETA PSF', 'title', font_size=16) @@ -82,8 +82,10 @@ brain_psf.add_foci(vert_max_psf, coords_as_verts=True, scale_factor=1., hemi='lh', color='black') -brain_ctf = stc_ctf.plot('sample', 'inflated', 'lh', subjects_dir=subjects_dir, - figure=2) +############################################################################### +# CTF: + +brain_ctf = stc_ctf.plot('sample', 'inflated', 'lh', subjects_dir=subjects_dir) brain_ctf.add_text(0.1, 0.9, 'sLORETA CTF', 'title', font_size=16) brain_ctf.show_view('ventral') brain_ctf.add_foci(verttrue, coords_as_verts=True, scale_factor=1., hemi='lh', @@ -93,5 +95,7 @@ brain_ctf.add_foci(vert_max_ctf, coords_as_verts=True, scale_factor=1., hemi='lh', color='black') -print('The green spheres indicate the true source location, and the black ' - 'spheres the maximum of the distribution.') + +############################################################################### +# The green spheres indicate the true source location, and the black +# spheres the maximum of the distribution. diff --git a/examples/inverse/plot_psf_ctf_vertices_lcmv.py b/examples/inverse/plot_psf_ctf_vertices_lcmv.py index 6f4cfb6ebff..8c274e39a12 100644 --- a/examples/inverse/plot_psf_ctf_vertices_lcmv.py +++ b/examples/inverse/plot_psf_ctf_vertices_lcmv.py @@ -27,7 +27,7 @@ raw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif' # Read raw data -raw = mne.io.read_raw_fif(raw_fname, preload=True) +raw = mne.io.read_raw_fif(raw_fname) # only pick good EEG/MEG sensors raw.info['bads'] += ['EEG 053'] # bads + 1 more @@ -42,6 +42,7 @@ tmin, tmax = -.2, .25 # epoch duration epochs = mne.Epochs(raw, events, event_id=event_id, tmin=tmin, tmax=tmax, picks=picks, baseline=(-.2, 0.), preload=True) +del raw # covariance matrix for pre-stimulus interval tmin, tmax = -.2, 0. @@ -52,20 +53,18 @@ tmin, tmax = 0.05, .25 cov_post = mne.compute_covariance(epochs, tmin=tmin, tmax=tmax, method='empirical') +info = epochs.info +del epochs # read forward solution forward = mne.read_forward_solution(fname_fwd) # use forward operator with fixed source orientations -forward = mne.convert_forward_solution(forward, surf_ori=True, - force_fixed=True) +mne.convert_forward_solution(forward, surf_ori=True, + force_fixed=True, copy=False) # read noise covariance matrix noise_cov = mne.read_cov(fname_cov) -# get valid measurement info -raw = raw.pick_types(meg=True, eeg=True, exclude='bads') -info = raw.info - # regularize noise covariance (we used 'empirical' above) noise_cov = mne.cov.regularize(noise_cov, info, mag=0.1, grad=0.1, eeg=0.1, rank='info') @@ -104,12 +103,13 @@ stc_pre = get_cross_talk(rm_pre, forward['src'], sources, norm=True) stc_post = get_cross_talk(rm_post, forward['src'], sources, norm=True) +verttrue = [forward['src'][0]['vertno'][sources[0]]] # pick one vertex +del forward ############################################################################## -# Visualise +# Visualize # --------- -vertno_lh = forward['src'][0]['vertno'] # vertex of selected source -verttrue = [vertno_lh[sources[0]]] # pick one vertex +# Pre: brain_pre = stc_pre.plot('sample', 'inflated', 'lh', subjects_dir=subjects_dir, figure=1, clim=dict(kind='value', lims=(0, .2, .4))) @@ -117,6 +117,13 @@ brain_pre.add_text(0.1, 0.9, 'LCMV beamformer with pre-stimulus\ndata ' 'covariance matrix', 'title', font_size=16) +# mark true source location for CTFs +brain_pre.add_foci(verttrue, coords_as_verts=True, scale_factor=1., hemi='lh', + color='green') + +############################################################################### +# Post: + brain_post = stc_post.plot('sample', 'inflated', 'lh', subjects_dir=subjects_dir, figure=2, clim=dict(kind='value', lims=(0, .2, .4))) @@ -124,10 +131,6 @@ brain_post.add_text(0.1, 0.9, 'LCMV beamformer with post-stimulus\ndata ' 'covariance matrix', 'title', font_size=16) -# mark true source location for CTFs -brain_pre.add_foci(verttrue, coords_as_verts=True, scale_factor=1., hemi='lh', - color='green') - brain_post.add_foci(verttrue, coords_as_verts=True, scale_factor=1., hemi='lh', color='green') diff --git a/examples/inverse/plot_resolution_metrics.py b/examples/inverse/plot_resolution_metrics.py index c3cca34e5d4..5ce714c442f 100644 --- a/examples/inverse/plot_resolution_metrics.py +++ b/examples/inverse/plot_resolution_metrics.py @@ -73,22 +73,26 @@ function='psf', metric='peak_err') sd_dspm_psf = resolution_metrics(rm_dspm, inverse_operator['src'], function='psf', metric='sd_ext') -del rm_dspm +del rm_dspm, forward ############################################################################### # Visualize results # ----------------- -# Visualise peak localisation error (PLE) across the whole cortex for PSF +# Visualise peak localisation error (PLE) across the whole cortex for MNE PSF: brain_ple_mne = ple_mne_psf.plot('sample', 'inflated', 'lh', subjects_dir=subjects_dir, figure=1, clim=dict(kind='value', lims=(0, 2, 4))) brain_ple_mne.add_text(0.1, 0.9, 'PLE MNE', 'title', font_size=16) +############################################################################### +# And dSPM: + brain_ple_dspm = ple_dspm_psf.plot('sample', 'inflated', 'lh', subjects_dir=subjects_dir, figure=2, clim=dict(kind='value', lims=(0, 2, 4))) brain_ple_dspm.add_text(0.1, 0.9, 'PLE dSPM', 'title', font_size=16) +############################################################################### # Subtract the two distributions and plot this difference diff_ple = ple_mne_psf - ple_dspm_psf @@ -102,19 +106,25 @@ # color) than MNE in deeper brain areas, but higher error (blue color) in more # superficial areas. # -# Next we'll visualise spatial deviation (SD) across the whole cortex for PSF: +# Next we'll visualise spatial deviation (SD) across the whole cortex for MNE +# PSF: brain_sd_mne = sd_mne_psf.plot('sample', 'inflated', 'lh', subjects_dir=subjects_dir, figure=4, clim=dict(kind='value', lims=(0, 2, 4))) brain_sd_mne.add_text(0.1, 0.9, 'SD MNE', 'title', font_size=16) +############################################################################### +# And dSPM: + brain_sd_dspm = sd_dspm_psf.plot('sample', 'inflated', 'lh', subjects_dir=subjects_dir, figure=5, clim=dict(kind='value', lims=(0, 2, 4))) brain_sd_dspm.add_text(0.1, 0.9, 'SD dSPM', 'title', font_size=16) -# Subtract the two distributions and plot this difference +############################################################################### +# Subtract the two distributions and plot this difference: + diff_sd = sd_mne_psf - sd_dspm_psf brain_sd_diff = diff_sd.plot('sample', 'inflated', 'lh', diff --git a/examples/inverse/plot_resolution_metrics_eegmeg.py b/examples/inverse/plot_resolution_metrics_eegmeg.py index cc98bec782f..3a7c90050af 100644 --- a/examples/inverse/plot_resolution_metrics_eegmeg.py +++ b/examples/inverse/plot_resolution_metrics_eegmeg.py @@ -94,13 +94,18 @@ brain_ple_emeg.add_text(0.1, 0.9, 'PLE PSF EMEG', 'title', font_size=16) +############################################################################### +# For MEG only: + brain_ple_meg = ple_psf_meg.plot('sample', 'inflated', 'lh', subjects_dir=subjects_dir, figure=2, clim=dict(kind='value', lims=(0, 2, 4))) brain_ple_meg.add_text(0.1, 0.9, 'PLE PSF MEG', 'title', font_size=16) -# Subtract the two distributions and plot this difference +############################################################################### +# Subtract the two distributions and plot this difference: + diff_ple = ple_psf_emeg - ple_psf_meg brain_ple_diff = diff_ple.plot('sample', 'inflated', 'lh', @@ -121,13 +126,18 @@ brain_sd_emeg.add_text(0.1, 0.9, 'SD PSF EMEG', 'title', font_size=16) +############################################################################### +# For MEG only: + brain_sd_meg = sd_psf_meg.plot('sample', 'inflated', 'lh', subjects_dir=subjects_dir, figure=5, clim=dict(kind='value', lims=(0, 2, 4))) brain_sd_meg.add_text(0.1, 0.9, 'SD PSF MEG', 'title', font_size=16) -# Subtract the two distributions and plot this difference +############################################################################### +# Subtract the two distributions and plot this difference: + diff_sd = sd_psf_emeg - sd_psf_meg brain_sd_diff = diff_sd.plot('sample', 'inflated', 'lh', diff --git a/examples/inverse/plot_vector_mne_solution.py b/examples/inverse/plot_vector_mne_solution.py index 46b4a7285e5..0a8983480fa 100644 --- a/examples/inverse/plot_vector_mne_solution.py +++ b/examples/inverse/plot_vector_mne_solution.py @@ -73,6 +73,10 @@ brain_max = stc_max.plot( initial_time=peak_time, hemi='lh', subjects_dir=subjects_dir, time_label='Max power') + +############################################################################### +# The normal is very similar: + brain_normal = stc.project('normal', inv['src'])[0].plot( initial_time=peak_time, hemi='lh', subjects_dir=subjects_dir, time_label='Normal') diff --git a/examples/preprocessing/plot_find_ref_artifacts.py b/examples/preprocessing/plot_find_ref_artifacts.py index 8fc714753c6..d72de7204e1 100644 --- a/examples/preprocessing/plot_find_ref_artifacts.py +++ b/examples/preprocessing/plot_find_ref_artifacts.py @@ -41,9 +41,10 @@ data_path = refmeg_noise.data_path() ############################################################################### -# Read raw data +# Read raw data, cropping to 5 minutes to save memory + raw_fname = data_path + '/sample_reference_MEG_noise-raw.fif' -raw = io.read_raw_fif(raw_fname, preload=True) +raw = io.read_raw_fif(raw_fname).crop(300, 600).load_data() ############################################################################### # Note that even though standard noise removal has already @@ -66,12 +67,14 @@ raw_tog = raw.copy() ica_kwargs = dict( method='picard', - fit_params=dict(tol=1e-3), # use a high tol here for speed + fit_params=dict(tol=1e-4), # use a high tol here for speed ) all_picks = mne.pick_types(raw_tog.info, meg=True, ref_meg=True) ica_tog = ICA(n_components=60, allow_ref_meg=True, **ica_kwargs) ica_tog.fit(raw_tog, picks=all_picks) -bad_comps, scores = ica_tog.find_bads_ref(raw_tog, threshold=2.5) +# low threshold (2.0) here because of cropped data, entire recording can use +# a higher threshold (2.5) +bad_comps, scores = ica_tog.find_bads_ref(raw_tog, threshold=2.0) # Plot scores with bad components marked. ica_tog.plot_scores(scores, bad_comps) diff --git a/examples/time_frequency/plot_source_power_spectrum_opm.py b/examples/time_frequency/plot_source_power_spectrum_opm.py index 65fb46c62bf..1ace1f65d4a 100644 --- a/examples/time_frequency/plot_source_power_spectrum_opm.py +++ b/examples/time_frequency/plot_source_power_spectrum_opm.py @@ -122,20 +122,30 @@ del src[0]['dist'], src[1]['dist'] bem = mne.read_bem_solution(bem_fname) fwd = dict() -trans = dict(vv=vv_trans_fname, opm=opm_trans_fname) -# check alignment and generate forward + +# check alignment and generate forward for VectorView +kwargs = dict(azimuth=0, elevation=90, distance=0.6, focalpoint=(0., 0., 0.)) +fig = mne.viz.plot_alignment( + raws['vv'].info, trans=vv_trans_fname, subject=subject, + subjects_dir=subjects_dir, dig=True, coord_frame='mri', + surfaces=('head', 'white')) +mne.viz.set_3d_view(figure=fig, **kwargs) +fwd['vv'] = mne.make_forward_solution( + raws['vv'].info, vv_trans_fname, src, bem, eeg=False, verbose=True) + +############################################################################## +# And for OPM: + with mne.use_coil_def(opm_coil_def_fname): - for kind in kinds: - dig = True if kind == 'vv' else False - fig = mne.viz.plot_alignment( - raws[kind].info, trans=trans[kind], subject=subject, - subjects_dir=subjects_dir, dig=dig, coord_frame='mri', - surfaces=('head', 'white')) - mne.viz.set_3d_view(figure=fig, azimuth=0, elevation=90, - distance=0.6, focalpoint=(0., 0., 0.)) - fwd[kind] = mne.make_forward_solution( - raws[kind].info, trans[kind], src, bem, eeg=False, verbose=True) -del trans, src, bem + fig = mne.viz.plot_alignment( + raws['opm'].info, trans=opm_trans_fname, subject=subject, + subjects_dir=subjects_dir, dig=False, coord_frame='mri', + surfaces=('head', 'white')) + mne.viz.set_3d_view(figure=fig, **kwargs) + fwd['opm'] = mne.make_forward_solution( + raws['opm'].info, opm_trans_fname, src, bem, eeg=False, verbose=True) + +del src, bem ############################################################################## # Compute and apply inverse to PSD estimated using multitaper + Welch. @@ -187,6 +197,7 @@ def plot_band(kind, band): brain = stcs[kind][band].plot( subject=subject, subjects_dir=subjects_dir, views='cau', hemi='both', time_label=title, title=title, colormap='inferno', + time_viewer=False, show_traces=False, clim=dict(kind='percent', lims=(70, 85, 99)), smoothing_steps=10) brain.show_view(dict(azimuth=0, elevation=0), roll=0) return fig, brain @@ -204,9 +215,12 @@ def plot_band(kind, band): # Beta # ---- # Here we also show OPM data, which shows a profile similar to the VectorView -# data beneath the sensors. +# data beneath the sensors. VectorView first: fig_beta, brain_beta = plot_band('vv', 'beta') + +############################################################################### +# Then OPM: fig_beta_opm, brain_beta_opm = plot_band('opm', 'beta') ############################################################################### diff --git a/examples/time_frequency/plot_time_frequency_erds.py b/examples/time_frequency/plot_time_frequency_erds.py index 67311c20f70..5d18f309ae5 100644 --- a/examples/time_frequency/plot_time_frequency_erds.py +++ b/examples/time_frequency/plot_time_frequency_erds.py @@ -112,7 +112,7 @@ ax.set_title(epochs.ch_names[ch], fontsize=10) ax.axvline(0, linewidth=1, color="black", linestyle=":") # event - if not ax.is_first_col(): + if ch != 0: ax.set_ylabel("") ax.set_yticklabels("") fig.colorbar(axes[0].images[-1], cax=axes[-1]) diff --git a/examples/visualization/plot_3d_to_2d.py b/examples/visualization/plot_3d_to_2d.py index 6d0a1da4bbe..56276bbe394 100644 --- a/examples/visualization/plot_3d_to_2d.py +++ b/examples/visualization/plot_3d_to_2d.py @@ -44,7 +44,7 @@ # Load data # --------- # -# First we'll load a sample ECoG dataset which we'll use for generating +# First we will load a sample ECoG dataset which we'll use for generating # a 2D snapshot. mat = loadmat(path_data) diff --git a/examples/visualization/plot_meg_sensors.py b/examples/visualization/plot_meg_sensors.py index 037f6b9f24a..01414f8cf97 100644 --- a/examples/visualization/plot_meg_sensors.py +++ b/examples/visualization/plot_meg_sensors.py @@ -6,6 +6,10 @@ ====================================== Show sensor layouts of different MEG systems. + +.. contents:: + :local: + :depth: 1 """ # Author: Eric Larson # @@ -21,27 +25,51 @@ print(__doc__) +############################################################################### +# Neuromag +# -------- + +kwargs = dict(eeg=False, coord_frame='meg', show_axes=True, verbose=True) + +raw = read_raw_fif(sample.data_path() + '/MEG/sample/sample_audvis_raw.fif') +fig = plot_alignment(raw.info, meg=('helmet', 'sensors'), **kwargs) +set_3d_title(figure=fig, title='Neuromag') + +############################################################################### +# CTF +# --- + +raw = read_raw_ctf(spm_face.data_path() + + '/MEG/spm/SPM_CTF_MEG_example_faces1_3D.ds') +fig = plot_alignment(raw.info, meg=('helmet', 'sensors', 'ref'), **kwargs) +set_3d_title(figure=fig, title='CTF 275') + +############################################################################### +# BTi +# --- + bti_path = op.abspath(op.dirname(mne.__file__)) + '/io/bti/tests/data/' +raw = read_raw_bti(op.join(bti_path, 'test_pdf_linux'), + op.join(bti_path, 'test_config_linux'), + op.join(bti_path, 'test_hs_linux')) +fig = plot_alignment(raw.info, meg=('helmet', 'sensors', 'ref'), **kwargs) +set_3d_title(figure=fig, title='Magnes 3600wh') + +############################################################################### +# KIT +# --- + kit_path = op.abspath(op.dirname(mne.__file__)) + '/io/kit/tests/data/' -raws = { - 'Neuromag': read_raw_fif(sample.data_path() + - '/MEG/sample/sample_audvis_raw.fif'), - 'CTF 275': read_raw_ctf(spm_face.data_path() + - '/MEG/spm/SPM_CTF_MEG_example_faces1_3D.ds'), - 'Magnes 3600wh': read_raw_bti(op.join(bti_path, 'test_pdf_linux'), - op.join(bti_path, 'test_config_linux'), - op.join(bti_path, 'test_hs_linux')), - 'KIT': read_raw_kit(op.join(kit_path, 'test.sqd')), - 'Artemis123': read_raw_artemis123(op.join( - testing.data_path(), 'ARTEMIS123', - 'Artemis_Data_2017-04-14-10h-38m-59s_Phantom_1k_HPI_1s.bin')), -} - -for system, raw in sorted(raws.items()): - meg = ['helmet', 'sensors'] - # We don't have coil definitions for KIT refs, so exclude them - if system != 'KIT': - meg.append('ref') - fig = plot_alignment(raw.info, eeg=False, meg=('helmet', 'sensors'), - coord_frame='meg', show_axes=True, verbose=True) - set_3d_title(figure=fig, title=system) +raw = read_raw_kit(op.join(kit_path, 'test.sqd')) +fig = plot_alignment(raw.info, meg=('helmet', 'sensors'), **kwargs) +set_3d_title(figure=fig, title='KIT') + +############################################################################### +# Artemis123 +# ---------- + +raw = read_raw_artemis123(op.join( + testing.data_path(), 'ARTEMIS123', + 'Artemis_Data_2017-04-14-10h-38m-59s_Phantom_1k_HPI_1s.bin')) +fig = plot_alignment(raw.info, meg=('helmet', 'sensors', 'ref'), **kwargs) +set_3d_title(figure=fig, title='Artemis123') diff --git a/mne/bem.py b/mne/bem.py index 4aed0b673bf..274232f9712 100644 --- a/mne/bem.py +++ b/mne/bem.py @@ -197,8 +197,8 @@ def _fwd_bem_lin_pot_coeff(surfs): # if sidx1 == sidx2 and (tri == j).any(): # continue # Otherwise do the hard job - coeffs = _lin_pot_coeff(surf1['rr'], tri_rr[k], tri_nn[k], - tri_area[k]) + coeffs = _lin_pot_coeff(fros=surf1['rr'], tri_rr=tri_rr[k], + tri_nn=tri_nn[k], tri_area=tri_area[k]) coeffs[skip_idx] = 0. submat[:, tri] -= coeffs if si_1 == si_2: @@ -233,7 +233,7 @@ def _fwd_bem_multi_solution(solids, gamma, nps): def _fwd_bem_homog_solution(solids, nps): """Make a homogeneous solution.""" - return _fwd_bem_multi_solution(solids, None, nps) + return _fwd_bem_multi_solution(solids, gamma=None, nps=nps) def _fwd_bem_ip_modify_solution(solution, ip_solution, ip_mult, n_tri): @@ -273,33 +273,33 @@ def _check_complete_surface(surf, copy=False, incomplete='raise', extra=''): return surf -def _fwd_bem_linear_collocation_solution(m): +def _fwd_bem_linear_collocation_solution(bem): """Compute the linear collocation potential solution.""" # first, add surface geometries - for surf in m['surfs']: + for surf in bem['surfs']: _check_complete_surface(surf) logger.info('Computing the linear collocation solution...') logger.info(' Matrix coefficients...') - coeff = _fwd_bem_lin_pot_coeff(m['surfs']) - m['nsol'] = len(coeff) + coeff = _fwd_bem_lin_pot_coeff(bem['surfs']) + bem['nsol'] = len(coeff) logger.info(" Inverting the coefficient matrix...") - nps = [surf['np'] for surf in m['surfs']] - m['solution'] = _fwd_bem_multi_solution(coeff, m['gamma'], nps) - if len(m['surfs']) == 3: - ip_mult = m['sigma'][1] / m['sigma'][2] + nps = [surf['np'] for surf in bem['surfs']] + bem['solution'] = _fwd_bem_multi_solution(coeff, bem['gamma'], nps) + if len(bem['surfs']) == 3: + ip_mult = bem['sigma'][1] / bem['sigma'][2] if ip_mult <= FWD.BEM_IP_APPROACH_LIMIT: logger.info('IP approach required...') logger.info(' Matrix coefficients (homog)...') - coeff = _fwd_bem_lin_pot_coeff([m['surfs'][-1]]) + coeff = _fwd_bem_lin_pot_coeff([bem['surfs'][-1]]) logger.info(' Inverting the coefficient matrix (homog)...') ip_solution = _fwd_bem_homog_solution(coeff, - [m['surfs'][-1]['np']]) + [bem['surfs'][-1]['np']]) logger.info(' Modify the original solution to incorporate ' 'IP approach...') - _fwd_bem_ip_modify_solution(m['solution'], ip_solution, ip_mult, + _fwd_bem_ip_modify_solution(bem['solution'], ip_solution, ip_mult, nps) - m['bem_method'] = FWD.BEM_LINEAR_COLL + bem['bem_method'] = FWD.BEM_LINEAR_COLL logger.info("Solution ready.") diff --git a/mne/conftest.py b/mne/conftest.py index aff52d4bfff..85bad3432c9 100644 --- a/mne/conftest.py +++ b/mne/conftest.py @@ -28,7 +28,7 @@ import numpy as np import mne from mne.datasets import testing -from mne.utils import _pl +from mne.utils import _pl, _assert_no_instances test_path = testing.data_path(download=False) s_path = op.join(test_path, 'MEG', 'sample') @@ -471,6 +471,39 @@ def download_is_error(monkeypatch): monkeypatch.setattr(mne.utils.fetching, '_get_http', _fail) +@pytest.fixture() +def brain_gc(request): + """Ensure that brain can be properly garbage collected.""" + keys = ('renderer_interactive', 'renderer') + assert set(request.fixturenames) & set(keys) != set() + for key in keys: + if key in request.fixturenames: + is_pv = request.getfixturevalue(key)._get_3d_backend() == 'pyvista' + break + if not is_pv: + yield + return + from mne.viz import Brain + _assert_no_instances(Brain, 'before') + ignore = set(id(o) for o in gc.get_objects()) + yield + _assert_no_instances(Brain, 'after') + # We only check VTK for PyVista -- Mayavi/PySurfer is not as strict + objs = gc.get_objects() + bad = list() + for o in objs: + try: + name = o.__class__.__name__ + except Exception: # old Python, probably + pass + else: + if name.startswith('vtk') and id(o) not in ignore: + bad.append(name) + del o + del objs, ignore, Brain + assert len(bad) == 0, 'VTK objects linger:\n' + '\n'.join(bad) + + def pytest_sessionfinish(session, exitstatus): """Handle the end of the session.""" n = session.config.option.durations diff --git a/mne/fixes.py b/mne/fixes.py index 868eb40e73d..f89b73da407 100644 --- a/mne/fixes.py +++ b/mne/fixes.py @@ -916,12 +916,21 @@ def _crop_colorbar(cbar, cbar_vmin, cbar_vmax): # matplotlib >= 3.2.0 no longer normalizes axes between 0 and 1 # See https://matplotlib.org/3.2.1/api/prev_api_changes/api_changes_3.2.0.html + # _outline was removed in + # https://github.com/matplotlib/matplotlib/commit/03a542e875eba091a027046d5ec652daa8be6863 + # so we use the code from there if LooseVersion(matplotlib.__version__) >= LooseVersion("3.2.0"): cbar.ax.set_ylim(cbar_vmin, cbar_vmax) X, _ = cbar._mesh() - new_X = np.array([X[0], X[-1]]) - new_Y = np.array([[cbar_vmin, cbar_vmin], [cbar_vmax, cbar_vmax]]) - xy = cbar._outline(new_X, new_Y) + X = np.array([X[0], X[-1]]) + Y = np.array([[cbar_vmin, cbar_vmin], [cbar_vmax, cbar_vmax]]) + N = X.shape[0] + ii = [0, 1, N - 2, N - 1, 2 * N - 1, 2 * N - 2, N + 1, N, 0] + x = X.T.reshape(-1)[ii] + y = Y.T.reshape(-1)[ii] + xy = (np.column_stack([y, x]) + if cbar.orientation == 'horizontal' else + np.column_stack([x, y])) cbar.outline.set_xy(xy) else: cbar.ax.set_ylim(cbar.norm(cbar_vmin), cbar.norm(cbar_vmax)) diff --git a/mne/forward/_compute_forward.py b/mne/forward/_compute_forward.py index a2cbb745a63..e8c09c5e489 100644 --- a/mne/forward/_compute_forward.py +++ b/mne/forward/_compute_forward.py @@ -181,7 +181,7 @@ def _bem_specify_coils(bem, coils, coord_frame, mults, n_jobs): Parameters ---------- - bem : dict + bem : instance of ConductorModel BEM information coils : list of dict, len(n_MEG_sensors) MEG sensor information dicts @@ -231,7 +231,7 @@ def _bem_specify_els(bem, els, mults): Parameters ---------- - bem : dict + bem : instance of ConductorModel BEM information els : list of dict, len(n_EEG_sensors) List of EEG sensor information dicts @@ -732,7 +732,7 @@ def _prep_field_computation(rr, bem, fwd_data, n_jobs, verbose=None): ---------- rr : ndarray, shape (n_dipoles, 3) 3D dipole source positions in head coordinates - bem : dict + bem : instance of ConductorModel Boundary Element Model information fwd_data : dict Dict containing sensor information. Gets updated here with BEM and @@ -893,7 +893,7 @@ def _compute_forwards(rr, bem, coils_list, ccoils_list, infos, coil_types, ---------- rr : ndarray, shape (n_sources, 3) 3D dipole in head coordinates - bem : dict + bem : instance of ConductorModel Boundary Element Model information for all surfaces coils_list : list List of MEG and/or EEG sensor information dicts diff --git a/mne/minimum_norm/spatial_resolution.py b/mne/minimum_norm/spatial_resolution.py index 8e372aad0c0..99704a695e1 100644 --- a/mne/minimum_norm/spatial_resolution.py +++ b/mne/minimum_norm/spatial_resolution.py @@ -146,7 +146,8 @@ def _localisation_error(resmat, src, function, metric): resmat = _rectify_resolution_matrix(resmat) locations = _get_src_locations(src) # locs used in forw. and inv. operator locations = 100. * locations # convert to cm (more common) - resmat = np.absolute(resmat) # only use absolute values + # we want to use absolute values, but doing abs() mases a copy and this + # can be quite expensive in memory. So let's just use abs() in place below. # The code below will operate on columns, so transpose if you want CTFs if function == 'ctf': @@ -154,16 +155,16 @@ def _localisation_error(resmat, src, function, metric): # Euclidean distance between true location and maximum if metric == 'peak_err': - resmax = resmat.argmax(axis=0) # find indices of maxima along columns + resmax = [abs(col).argmax() for col in resmat.T] # max inds along cols maxloc = locations[resmax, :] # locations of maxima diffloc = locations - maxloc # diff btw true locs and maxima locs - locerr = np.sqrt(np.sum(diffloc ** 2, 1)) # Euclidean distance + locerr = np.linalg.norm(diffloc, axis=1) # Euclidean distance # centre of gravity elif metric == 'cog_err': locerr = np.empty(locations.shape[0]) # initialise result array for ii, rr in enumerate(locations): - resvec = resmat[:, ii].T # corresponding column of resmat + resvec = abs(resmat[:, ii].T) # corresponding column of resmat cog = resvec.dot(locations) / np.sum(resvec) # centre of gravity locerr[ii] = np.sqrt(np.sum((rr - cog) ** 2)) # Euclidean distance @@ -202,14 +203,12 @@ def _spatial_extent(resmat, src, function, metric, threshold=0.5): """ locations = _get_src_locations(src) # locs used in forw. and inv. operator locations = 100. * locations # convert to cm (more common) - resmat = np.absolute(resmat) # only use absolute values # The code below will operate on columns, so transpose if you want CTFs if function == 'ctf': resmat = resmat.T - resmax = resmat.argmax(axis=0) # find indices of maxima along rows - width = np.empty(len(resmax)) # initialise output array + width = np.empty(resmat.shape[1]) # initialise output array # spatial deviation as in Molins et al. if metric == 'sd_ext': @@ -217,16 +216,16 @@ def _spatial_extent(resmat, src, function, metric, threshold=0.5): diffloc = locations - locations[ii, :] # locs w/r/t true source locerr = np.sum(diffloc**2, 1) # squared Eucl dists to true source - resvec = resmat[:, ii]**2 # pick current row + resvec = abs(resmat[:, ii]) ** 2 # pick current row # spatial deviation (Molins et al, NI 2008, eq. 12) width[ii] = np.sqrt(np.sum(np.multiply(locerr, resvec)) / np.sum(resvec)) # maximum radius to 50% of max amplitude elif metric == 'maxrad_ext': - maxamp = resmat.max(axis=0) # peak ampl. per location across columns - for ii, amps in enumerate(maxamp): # for all locations - resvec = resmat[:, ii] # pick current column + for ii, resvec in enumerate(resmat.T): # iterate over columns + resvec = abs(resvec) # operate on absolute values + amps = resvec.max() # indices of elements with values larger than fraction threshold # of peak amplitude thresh_idx = np.where(resvec > threshold * amps) @@ -269,17 +268,17 @@ def _relative_amplitude(resmat, src, function, metric): if function == 'ctf': resmat = resmat.T - resmat = np.absolute(resmat) # only use absolute values - # Ratio between amplitude at peak and global peak maximum if metric == 'peak_amp': - maxamps = resmat.max(axis=0) # maximum amplitudes per column + # maximum amplitudes per column + maxamps = np.array([abs(col).max() for col in resmat.T]) maxmaxamps = maxamps.max() # global absolute maximum relamp = maxamps / maxmaxamps # ratio between sums of absolute amplitudes elif metric == 'sum_amp': - sumamps = np.sum(resmat, axis=0) # sum of amplitudes per column + # sum of amplitudes per column + sumamps = np.array([abs(col).sum() for col in resmat.T]) sumampsmax = sumamps.max() # maximum of summed amplitudes relamp = sumamps / sumampsmax diff --git a/mne/morph.py b/mne/morph.py index f4fe49bd2a1..723e667e7fb 100644 --- a/mne/morph.py +++ b/mne/morph.py @@ -711,7 +711,7 @@ def _get_src_data(src, mri_resolution=True): shape = use_src['shape'] start = 0 if src_kind == 'volume' else 2 for si, s in enumerate(src_t[start:], start): - if s['interpolator'] is None: + if s.get('interpolator', None) is None: if mri_resolution: raise RuntimeError( 'MRI interpolator not present in src[%d], ' diff --git a/mne/preprocessing/tests/test_maxwell.py b/mne/preprocessing/tests/test_maxwell.py index f00b93d9fa7..ccdbe972934 100644 --- a/mne/preprocessing/tests/test_maxwell.py +++ b/mne/preprocessing/tests/test_maxwell.py @@ -28,7 +28,7 @@ _sh_real_to_complex, _sh_negate, _bases_complex_to_real, _trans_sss_basis, _bases_real_to_complex, _prep_mf_coils, find_bad_channels_maxwell) from mne.rank import _get_rank_sss, _compute_rank_int -from mne.utils import (assert_meg_snr, run_tests_if_main, catch_logging, +from mne.utils import (assert_meg_snr, catch_logging, object_diff, buggy_mkl_svd, use_log_level) data_path = testing.data_path(download=False) @@ -1339,14 +1339,11 @@ def test_find_bad_channels_maxwell(fname, bads, annot, add_ch, ignore_ref, # Check "flat" scores. scores_flat = got_scores['scores_flat'] limits_flat = got_scores['limits_flat'] - # The following essentially is just this: - # n_segments_below_limit = (scores_flat < limits_flat).sum(-1) - # made to work with NaN's in the scores. - n_segments_below_limit = np.less( - scores_flat, limits_flat, - where=np.equal(np.isnan(scores_flat), False), - out=np.full_like(scores_flat, fill_value=False)).sum(-1) - + # Deal with NaN's in the scores (can't use np.less directly because of + # https://github.com/numpy/numpy/issues/17198) + scores_flat[np.isnan(scores_flat)] = np.inf + limits_flat[np.isnan(limits_flat)] = -np.inf + n_segments_below_limit = (scores_flat < limits_flat).sum(-1) ch_idx = np.where(n_segments_below_limit >= min(min_count, len(got_scores['bins']))) flats = set(got_scores['ch_names'][ch_idx]) @@ -1355,18 +1352,10 @@ def test_find_bad_channels_maxwell(fname, bads, annot, add_ch, ignore_ref, # Check "noisy" scores. scores_noisy = got_scores['scores_noisy'] limits_noisy = got_scores['limits_noisy'] - # The following essentially is just this: - # n_segments_beyond_limit = (scores_noisy > limits_noisy).sum(-1) - # made to work with NaN's in the scores. - n_segments_beyond_limit = np.greater( - scores_noisy, limits_noisy, - where=np.equal(np.isnan(scores_noisy), False), - out=np.full_like(scores_noisy, fill_value=False)).sum(-1) - + scores_noisy[np.isnan(scores_noisy)] = -np.inf + limits_noisy[np.isnan(limits_noisy)] = np.inf + n_segments_beyond_limit = (scores_noisy > limits_noisy).sum(-1) ch_idx = np.where(n_segments_beyond_limit >= min(min_count, len(got_scores['bins']))) bads = set(got_scores['ch_names'][ch_idx]) assert bads == set(want_bads) - - -run_tests_if_main() diff --git a/mne/simulation/raw.py b/mne/simulation/raw.py index b5d2c18daf4..4d07829f17e 100644 --- a/mne/simulation/raw.py +++ b/mne/simulation/raw.py @@ -321,10 +321,10 @@ def simulate_raw(info, stc=None, trans=None, src=None, bem=None, head_pos=None, logger.info(' %d STC iteration%s provided' % (n + 1, _pl(n + 1))) break + del fwd else: raise RuntimeError('Maximum number of STC iterations (%d) ' 'exceeded' % (n,)) - del fwd raw_data = np.concatenate(raw_datas, axis=-1) raw = RawArray(raw_data, info, first_samp=first_samp, verbose=False) raw.set_annotations(raw.annotations) diff --git a/mne/source_space.py b/mne/source_space.py index 25f9c7105dd..e2d44835e7b 100644 --- a/mne/source_space.py +++ b/mne/source_space.py @@ -305,6 +305,22 @@ def copy(self): """ return deepcopy(self) + def __deepcopy__(self, memodict): + """Make a deepcopy.""" + # don't copy read-only views (saves a ton of mem for split-vol src) + info = deepcopy(self.info, memodict) + ss = list() + for s in self: + for key in ('rr', 'nn'): + if key in s: + arr = s[key] + id_ = id(arr) + if id_ not in memodict: + if not arr.flags.writeable: + memodict[id_] = arr + ss.append(deepcopy(s, memodict)) + return SourceSpaces(ss, info) + def save(self, fname, overwrite=False): """Save the source spaces to a fif file. @@ -753,6 +769,8 @@ def _read_one_source_space(fid, this): tag = find_tag(fid, mri, FIFF.FIFF_MNE_SOURCE_SPACE_INTERPOLATOR) if tag is not None: res['interpolator'] = tag.data + if tag.data.data.size == 0: + del res['interpolator'] else: logger.info("Interpolation matrix for MRI not found.") @@ -1139,16 +1157,20 @@ def _write_one_source_space(fid, this, verbose=None): if mri_volume_name is not None: write_string(fid, FIFF.FIFF_MNE_FILE_NAME, mri_volume_name) + mri_width, mri_height, mri_depth, nvox = _src_vol_dims(this) + interpolator = this.get('interpolator') + if interpolator is None: + interpolator = sparse.csr_matrix((nvox, this['np'])) write_float_sparse_rcs(fid, FIFF.FIFF_MNE_SOURCE_SPACE_INTERPOLATOR, - this['interpolator']) + interpolator) if 'mri_file' in this and this['mri_file'] is not None: write_string(fid, FIFF.FIFF_MNE_SOURCE_SPACE_MRI_FILE, this['mri_file']) - write_int(fid, FIFF.FIFF_MRI_WIDTH, this['mri_width']) - write_int(fid, FIFF.FIFF_MRI_HEIGHT, this['mri_height']) - write_int(fid, FIFF.FIFF_MRI_DEPTH, this['mri_depth']) + write_int(fid, FIFF.FIFF_MRI_WIDTH, mri_width) + write_int(fid, FIFF.FIFF_MRI_HEIGHT, mri_height) + write_int(fid, FIFF.FIFF_MRI_DEPTH, mri_depth) end_block(fid, FIFF.FIFFB_MNE_PARENT_MRI_FILE) @@ -1819,7 +1841,8 @@ def setup_volume_source_space(subject=None, pos=5.0, mri=None, # Compute an interpolation matrix to show data in MRI_VOXEL coord frame if mri is not None: - _add_interpolator(sp, add_interpolator) + if add_interpolator: + _add_interpolator(sp) elif sp[0]['type'] == 'vol': # If there is no interpolator, it's actually a discrete source space sp[0]['type'] = 'discrete' @@ -2030,7 +2053,7 @@ def _make_volume_source_space(surf, grid, exclude, mindist, mri=None, x, y, z = rr[2].ravel(), rr[1].ravel(), rr[0].ravel() rr = np.array([x * grid, y * grid, z * grid]).T sp = dict(np=npts, nn=np.zeros((npts, 3)), rr=rr, - inuse=np.ones(npts, int), type='vol', nuse=npts, + inuse=np.ones(npts, bool), type='vol', nuse=npts, coord_frame=FIFF.FIFFV_COORD_MRI, id=-1, shape=ns) sp['nn'][:, 2] = 1.0 assert sp['rr'].shape[0] == npts @@ -2070,8 +2093,15 @@ def _make_volume_source_space(surf, grid, exclude, mindist, mri=None, 'do_neighbors is True') sps = list() orig_sp = sp + # reduce the sizes when we deepcopy for volume_label, id_ in volume_labels.items(): - sp = deepcopy(orig_sp) + # this saves us some memory + memodict = dict() + for key in ('rr', 'nn'): + if key in orig_sp: + arr = orig_sp[key] + memodict[id(arr)] = arr + sp = deepcopy(orig_sp, memodict) good = _get_atlas_values(vol_info, sp['rr'][sp['vertno']]) == id_ n_good = good.sum() logger.info(' Selected %d voxel%s from %s' @@ -2083,6 +2113,7 @@ def _make_volume_source_space(surf, grid, exclude, mindist, mri=None, sp['seg_name'] = volume_label sp['mri_file'] = mri sps.append(sp) + del orig_sp assert len(sps) == len(volume_labels) # This will undo some of the work above, but the calculations are # pretty trivial so allow it @@ -2249,18 +2280,16 @@ def _get_mgz_header(fname): return header -def _add_interpolator(sp, add_interpolator): +def _src_vol_dims(s): + w, h, d = [s[f'mri_{key}'] for key in ('width', 'height', 'depth')] + return w, h, d, np.prod([w, h, d]) + + +def _add_interpolator(sp): """Compute a sparse matrix to interpolate the data into an MRI volume.""" # extract transformation information from mri s = sp[0] - mri_width = s['mri_width'] - mri_height = s['mri_height'] - mri_depth = s['mri_depth'] - nvox = mri_width * mri_height * mri_depth - if not add_interpolator: - for s in sp: - s['interpolator'] = sparse.csr_matrix((nvox, s['np'])) - return + mri_width, mri_height, mri_depth, nvox = _src_vol_dims(s) # # Convert MRI voxels from destination (MRI volume) to source (volume @@ -2801,7 +2830,7 @@ def _get_vertex_map_nn(fro_src, subject_from, subject_to, hemi, subjects_dir, reg_to['neighbor_tri'] = _triangle_neighbors(reg_to['tris'], reg_to['np']) - morph_inuse = np.zeros(len(reg_to['rr']), bool) + morph_inuse = np.zeros(len(reg_to['rr']), int) best = np.zeros(fro_src['np'], int) ones = _compute_nearest(reg_to['rr'], reg_fro['rr'][fro_src['vertno']]) for v, one in zip(fro_src['vertno'], ones): @@ -2986,6 +3015,8 @@ def _compare_source_spaces(src0, src1, mode='exact', nearest=True, assert_equal(s0[name], s1[name], name) for name in ['interpolator']: if name in s0 or name in s1: + assert name in s0, f'{name} in s1 but not s0' + assert name in s1, f'{name} in s1 but not s0' diffs = (s0['interpolator'] - s1['interpolator']).data if len(diffs) > 0 and 'nointerp' not in mode: # 5% diff --git a/mne/tests/test_morph.py b/mne/tests/test_morph.py index 529481830a9..70327eecece 100644 --- a/mne/tests/test_morph.py +++ b/mne/tests/test_morph.py @@ -488,7 +488,7 @@ def test_volume_source_morph_round_trip( src['fsaverage'] = mne.read_source_spaces(fname_fs_vol) src['fsaverage'][0].update( vol_dims=np.array([23, 29, 25]), seg_name='brain') - _add_interpolator(src['fsaverage'], True) + _add_interpolator(src['fsaverage']) assert src['fsaverage'][0]['nuse'] == 6379 src_to, src_from = src[subject_to], src[subject_from] del src diff --git a/mne/utils/__init__.py b/mne/utils/__init__.py index 2b1965dcf88..f3bd8e6ff62 100644 --- a/mne/utils/__init__.py +++ b/mne/utils/__init__.py @@ -32,7 +32,8 @@ ClosingStringIO) from .misc import (run_subprocess, _pl, _clean_names, pformat, _file_like, _explain_exception, _get_argvalues, sizeof_fmt, - running_subprocess, _DefaultEventParser) + running_subprocess, _DefaultEventParser, + _assert_no_instances) from .progressbar import ProgressBar from ._testing import (run_tests_if_main, run_command_if_main, requires_sklearn, diff --git a/mne/utils/misc.py b/mne/utils/misc.py index f6b21e7c368..af5b429048f 100644 --- a/mne/utils/misc.py +++ b/mne/utils/misc.py @@ -6,6 +6,7 @@ from contextlib import contextmanager import fnmatch +import gc import inspect from math import log import os @@ -319,3 +320,32 @@ def _file_like(obj): # but this might be more robust to file-like objects not properly # inheriting from these classes: return all(callable(getattr(obj, name, None)) for name in ('read', 'seek')) + + +def _assert_no_instances(cls, when=''): + __tracebackhide__ = True + n = 0 + ref = list() + gc.collect() + objs = gc.get_objects() + for obj in objs: + try: + check = isinstance(obj, cls) + except Exception: # such as a weakref + check = False + if check: + rr = gc.get_referrers(obj) + count = 0 + for r in rr: + if r is not objs and \ + r is not globals() and \ + r is not locals() and \ + not inspect.isframe(r): + ref.append( + f'{r.__class__.__name__}: ' + + repr(r)[:100].replace('\n', ' ')) + count += 1 + del r + del rr + n += count > 0 + assert n == 0, f'{n} {when}:\n' + '\n'.join(ref) diff --git a/mne/utils/numerics.py b/mne/utils/numerics.py index 415e21932e4..ad077b08549 100644 --- a/mne/utils/numerics.py +++ b/mne/utils/numerics.py @@ -651,7 +651,7 @@ def object_hash(x, h=None): return int(h.hexdigest(), 16) -def object_size(x): +def object_size(x, memo=None): """Estimate the size of a reasonable python object. Parameters @@ -660,6 +660,8 @@ def object_size(x): Object to approximate the size of. Can be anything comprised of nested versions of: {dict, list, tuple, ndarray, str, bytes, float, int, None}. + memo : dict | None + The memodict. Returns ------- @@ -668,28 +670,36 @@ def object_size(x): """ # Note: this will not process object arrays properly (since those only) # hold references + if memo is None: + memo = dict() + id_ = id(x) + if id_ in memo: + return 0 # do not add already existing ones if isinstance(x, (bytes, str, int, float, type(None))): size = sys.getsizeof(x) elif isinstance(x, np.ndarray): # On newer versions of NumPy, just doing sys.getsizeof(x) works, # but on older ones you always get something small :( - size = sys.getsizeof(np.array([])) + x.nbytes + size = sys.getsizeof(np.array([])) + if x.base is None or id(x.base) not in memo: + size += x.nbytes elif isinstance(x, np.generic): size = x.nbytes elif isinstance(x, dict): size = sys.getsizeof(x) for key, value in x.items(): - size += object_size(key) - size += object_size(value) + size += object_size(key, memo) + size += object_size(value, memo) elif isinstance(x, (list, tuple)): - size = sys.getsizeof(x) + sum(object_size(xx) for xx in x) + size = sys.getsizeof(x) + sum(object_size(xx, memo) for xx in x) elif isinstance(x, datetime): - size = object_size(_dt_to_stamp(x)) + size = object_size(_dt_to_stamp(x), memo) elif sparse.isspmatrix_csc(x) or sparse.isspmatrix_csr(x): size = sum(sys.getsizeof(xx) for xx in [x, x.data, x.indices, x.indptr]) else: raise RuntimeError('unsupported type: %s (%s)' % (type(x), x)) + memo[id_] = size return size diff --git a/mne/utils/tests/test_numerics.py b/mne/utils/tests/test_numerics.py index d63328f6e4f..3e8f2b27e3a 100644 --- a/mne/utils/tests/test_numerics.py +++ b/mne/utils/tests/test_numerics.py @@ -304,6 +304,21 @@ def test_object_size(): size = object_size(obj) assert lower < size < upper, \ '%s < %s < %s:\n%s' % (lower, size, upper, obj) + # views work properly + x = dict(a=1) + assert object_size(x) < 1000 + x['a'] = np.ones(100000, float) + nb = x['a'].nbytes + sz = object_size(x) + assert nb < sz < nb * 1.01 + x['b'] = x['a'] + sz = object_size(x) + assert nb < sz < nb * 1.01 + x['b'] = x['a'].view() + x['b'].flags.writeable = False + assert x['a'].flags.writeable + sz = object_size(x) + assert nb < sz < nb * 1.01 def test_object_diff_with_nan(): diff --git a/mne/viz/_3d.py b/mne/viz/_3d.py index 800723236b5..fdca94a1560 100644 --- a/mne/viz/_3d.py +++ b/mne/viz/_3d.py @@ -1701,9 +1701,9 @@ def plot_source_estimates(stc, subject=None, surface='inflated', hemi='lh', time_unit : 's' | 'ms' Whether time is represented in seconds ("s", default) or milliseconds ("ms"). - backend : 'auto' | 'mayavi' | 'matplotlib' + backend : 'auto' | 'mayavi' | 'pyvista' | 'matplotlib' Which backend to use. If ``'auto'`` (default), tries to plot with - mayavi, but resorts to matplotlib if mayavi is not available. + pyvista, but resorts to matplotlib if no 3d backend is available. .. versionadded:: 0.15.0 spacing : str @@ -1747,11 +1747,15 @@ def plot_source_estimates(stc, subject=None, surface='inflated', hemi='lh', subjects_dir = get_subjects_dir(subjects_dir=subjects_dir, raise_error=True) subject = _check_subject(stc.subject, subject, True) - _check_option('backend', backend, ['auto', 'matplotlib', 'mayavi']) + _check_option('backend', backend, + ['auto', 'matplotlib', 'mayavi', 'pyvista']) plot_mpl = backend == 'matplotlib' if not plot_mpl: try: - set_3d_backend(_get_3d_backend()) + if backend == 'auto': + set_3d_backend(_get_3d_backend()) + else: + set_3d_backend(backend) except (ImportError, ModuleNotFoundError): if backend == 'auto': warn('No 3D backend found. Resorting to matplotlib 3d.') diff --git a/mne/viz/_brain/_brain.py b/mne/viz/_brain/_brain.py index 1fed9c9d711..04fcaccb955 100644 --- a/mne/viz/_brain/_brain.py +++ b/mne/viz/_brain/_brain.py @@ -267,6 +267,12 @@ def __init__(self, subject_id, hemi, surf, title=None, bgcolor=background, shape=shape, fig=figure) + + if _get_3d_backend() == "pyvista": + self.plotter = self._renderer.plotter + self.window = self.plotter.app_window + self.window.signal_close.connect(self._clean) + for h in self._hemis: # Initialize a Surface object as the geometry geo = Surface(subject_id, h, surf, subjects_dir, offset, @@ -303,6 +309,7 @@ def __init__(self, subject_id, hemi, surf, title=None, actor, mesh = mesh_data.actor, mesh_data self._hemi_meshes[h] = mesh self._hemi_actors[h] = actor + del mesh_data, actor, mesh else: self._renderer.polydata( self._hemi_meshes[h], @@ -375,20 +382,18 @@ def setup_time_viewer(self, time_viewer=True, show_traces=True): self.slider_tube_color = (0.69803922, 0.70196078, 0.70980392) # Direct access parameters: - self.plotter = self._renderer.plotter + self._iren = self._renderer.plotter.iren self.main_menu = self.plotter.main_menu - self.window = self.plotter.app_window self.tool_bar = self.window.addToolBar("toolbar") self.status_bar = self.window.statusBar() self.interactor = self.plotter.interactor - self.window.signal_close.connect(self._clean) # Derived parameters: self.playback_speed = self.default_playback_speed_value _validate_type(show_traces, (bool, str, 'numeric'), 'show_traces') self.interactor_fraction = 0.25 if isinstance(show_traces, str): - assert 'show_traces' == 'separate' # should be guaranteed earlier + assert show_traces == 'separate' # should be guaranteed earlier self.show_traces = True self.separate_canvas = True else: @@ -425,23 +430,46 @@ def setup_time_viewer(self, time_viewer=True, show_traces=True): def _clean(self): # resolve the reference cycle self.clear_points() + for h in self._hemis: + # clear init actors + actor = self._hemi_actors[h] + if actor is not None: + mapper = actor.GetMapper() + mapper.SetLookupTable(None) + actor.SetMapper(None) + # clear data actors + hemi_data = self._data.get(h) + if hemi_data is not None: + if hemi_data['actors'] is not None: + for actor in hemi_data['actors']: + mapper = actor.GetMapper() + mapper.SetLookupTable(None) + actor.SetMapper(None) self._clear_callbacks() - self.actions.clear() - self.sliders.clear() - self.reps = None - self.plotter = None - self.main_menu = None - self.window = None - self.tool_bar = None - self.status_bar = None - self.interactor = None - if self.mpl_canvas is not None: + if getattr(self, 'mpl_canvas', None) is not None: self.mpl_canvas.clear() - self.mpl_canvas = None - self.time_actor = None - self.picked_renderer = None - for key in list(self.act_data_smooth.keys()): - self.act_data_smooth[key] = None + if getattr(self, 'act_data_smooth', None) is not None: + for key in list(self.act_data_smooth.keys()): + self.act_data_smooth[key] = None + # XXX this should be done in PyVista + for renderer in self.plotter.renderers: + renderer.RemoveAllLights() + # app_window cannot be set to None because it is used in __del__ + for key in ('lighting', 'interactor', '_RenderWindow'): + setattr(self.plotter, key, None) + # Qt LeaveEvent requires _Iren so we use _FakeIren instead of None + # to resolve the ref to vtkGenericRenderWindowInteractor + self.plotter._Iren = _FakeIren() + if getattr(self.plotter, 'scalar_bar', None) is not None: + self.plotter.scalar_bar = None + if getattr(self.plotter, 'picker', None) is not None: + self.plotter.picker = None + # XXX end PyVista + for key in ('reps', 'plotter', 'main_menu', 'window', 'tool_bar', + 'status_bar', 'interactor', 'mpl_canvas', 'time_actor', + 'picked_renderer', 'act_data_smooth', '_iren', + 'actions', 'sliders', 'geo', '_hemi_actors', '_data'): + setattr(self, key, None) @contextlib.contextmanager def ensure_minimum_sizes(self): @@ -922,6 +950,9 @@ def _load_icons(self): self.icons["visibility_on"] = QIcon(":/visibility_on.svg") self.icons["visibility_off"] = QIcon(":/visibility_off.svg") + def _save_movie_noname(self): + return self.save_movie(None) + def _configure_tool_bar(self): self.actions["screenshot"] = self.tool_bar.addAction( self.icons["screenshot"], @@ -931,7 +962,7 @@ def _configure_tool_bar(self): self.actions["movie"] = self.tool_bar.addAction( self.icons["movie"], "Save movie...", - partial(self.save_movie, filename=None) + self._save_movie_noname, ) self.actions["visibility"] = self.tool_bar.addAction( self.icons["visibility_on"], @@ -1208,6 +1239,8 @@ def remove_point(self, mesh): def clear_points(self): """Clear the picked points.""" + if not hasattr(self, '_spheres'): + return for sphere in list(self._spheres): # will remove itself, so copy self.remove_point(sphere) assert sum(len(v) for v in self.picked_points.values()) == 0 @@ -1308,6 +1341,9 @@ def help(self): ) def _clear_callbacks(self): + from ..backends._pyvista import _remove_picking_callback + if not hasattr(self, 'callbacks'): + return for callback in self.callbacks.values(): if callback is not None: if hasattr(callback, "plotter"): @@ -1317,6 +1353,8 @@ def _clear_callbacks(self): if hasattr(callback, "slider_rep"): callback.slider_rep = None self.callbacks.clear() + if self.show_traces: + _remove_picking_callback(self._iren, self.plotter.picker) @property def interaction(self): @@ -2923,3 +2961,17 @@ def _get_range(brain): def _normalize(point, shape): return (point[0] / shape[1], point[1] / shape[0]) + + +class _FakeIren(): + def EnterEvent(self): + pass + + def MouseMoveEvent(self): + pass + + def LeaveEvent(self): + pass + + def SetEventInformation(self, *args, **kwargs): + pass diff --git a/mne/viz/_brain/tests/test_brain.py b/mne/viz/_brain/tests/test_brain.py index 5ebd6b3f59f..b2907a27a42 100644 --- a/mne/viz/_brain/tests/test_brain.py +++ b/mne/viz/_brain/tests/test_brain.py @@ -93,7 +93,14 @@ def GetPosition(self): @testing.requires_testing_data -def test_brain_init(renderer, tmpdir, pixel_ratio): +def test_brain_gc(renderer, brain_gc): + """Test that a minimal version of Brain gets GC'ed.""" + brain = Brain('fsaverage', 'both', 'inflated', subjects_dir=subjects_dir) + brain.close() + + +@testing.requires_testing_data +def test_brain_init(renderer, tmpdir, pixel_ratio, brain_gc): """Test initialization of the Brain instance.""" from mne.label import read_label hemi = 'lh' @@ -105,14 +112,15 @@ def test_brain_init(renderer, tmpdir, pixel_ratio): kwargs = dict(subject_id=subject_id, subjects_dir=subjects_dir) with pytest.raises(ValueError, match='"size" parameter must be'): Brain(hemi=hemi, surf=surf, size=[1, 2, 3], **kwargs) + with pytest.raises(KeyError): + Brain(hemi='foo', surf=surf, **kwargs) with pytest.raises(TypeError, match='figure'): Brain(hemi=hemi, surf=surf, figure='foo', **kwargs) with pytest.raises(TypeError, match='interaction'): Brain(hemi=hemi, surf=surf, interaction=0, **kwargs) with pytest.raises(ValueError, match='interaction'): Brain(hemi=hemi, surf=surf, interaction='foo', **kwargs) - with pytest.raises(KeyError): - Brain(hemi='foo', surf=surf, **kwargs) + renderer.backend._close_all() brain = Brain(hemi=hemi, surf=surf, size=size, title=title, cortex=cortex, units='m', **kwargs) @@ -203,6 +211,8 @@ def test_brain_init(renderer, tmpdir, pixel_ratio): # add annotation annots = ['aparc', path.join(subjects_dir, 'fsaverage', 'label', 'lh.PALS_B12_Lobes.annot')] + brain.close() + borders = [True, 2] alphas = [1, 0.5] colors = [None, 'r'] @@ -229,7 +239,7 @@ def test_brain_init(renderer, tmpdir, pixel_ratio): @testing.requires_testing_data @pytest.mark.slowtest -def test_brain_save_movie(tmpdir, renderer): +def test_brain_save_movie(tmpdir, renderer, brain_gc): """Test saving a movie of a Brain instance.""" if renderer._get_3d_backend() == "mayavi": pytest.skip('Save movie only supported on PyVista') @@ -249,7 +259,7 @@ def test_brain_save_movie(tmpdir, renderer): @testing.requires_testing_data @pytest.mark.slowtest -def test_brain_time_viewer(renderer_interactive, pixel_ratio): +def test_brain_time_viewer(renderer_interactive, pixel_ratio, brain_gc): """Test time viewer primitives.""" if renderer_interactive._get_3d_backend() != 'pyvista': pytest.skip('TimeViewer tests only supported on PyVista') @@ -291,6 +301,7 @@ def test_brain_time_viewer(renderer_interactive, pixel_ratio): img = brain.screenshot(mode='rgb') want_shape = np.array([300 * pixel_ratio, 300 * pixel_ratio, 3]) assert_allclose(img.shape, want_shape) + brain.close() @testing.requires_testing_data @@ -306,7 +317,8 @@ def test_brain_time_viewer(renderer_interactive, pixel_ratio): pytest.param('mixed', marks=pytest.mark.slowtest), ]) @pytest.mark.slowtest -def test_brain_traces(renderer_interactive, hemi, src, tmpdir): +def test_brain_traces(renderer_interactive, hemi, src, tmpdir, + brain_gc): """Test brain traces.""" if renderer_interactive._get_3d_backend() != 'pyvista': pytest.skip('Only PyVista supports traces') @@ -399,6 +411,7 @@ def test_brain_traces(renderer_interactive, hemi, src, tmpdir): # only test one condition to save time if not (hemi == 'rh' and src == 'surface' and check_version('sphinx_gallery')): + brain.close() return fnames = [str(tmpdir.join(f'temp_{ii}.png')) for ii in range(2)] block_vars = dict(image_path_iterator=iter(fnames), @@ -412,6 +425,7 @@ def test_brain_traces(renderer_interactive, hemi, src, tmpdir): gallery_conf = dict(src_dir=str(tmpdir), compress_images=[]) scraper = _BrainScraper() rst = scraper(block, block_vars, gallery_conf) + assert brain.plotter is None # closed gif_0 = fnames[0][:-3] + 'gif' for fname in (gif_0, fnames[1]): assert path.basename(fname) in rst @@ -424,7 +438,7 @@ def test_brain_traces(renderer_interactive, hemi, src, tmpdir): @testing.requires_testing_data @pytest.mark.slowtest -def test_brain_linkviewer(renderer_interactive): +def test_brain_linkviewer(renderer_interactive, brain_gc): """Test _LinkViewer primitives.""" if renderer_interactive._get_3d_backend() != 'pyvista': pytest.skip('Linkviewer only supported on PyVista') @@ -455,9 +469,13 @@ def test_brain_linkviewer(renderer_interactive): link_viewer.set_fmax(1) link_viewer.set_playback_speed(value=0.1) link_viewer.toggle_playback() + del link_viewer + brain1.close() + brain2.close() + brain_data.close() -def test_brain_colormap(): +def test_calculate_lut(): """Test brain's colormap functions.""" colormap = "coolwarm" alpha = 1.0 diff --git a/mne/viz/backends/_pyvista.py b/mne/viz/backends/_pyvista.py index 15fc22a170b..d22ee9cd679 100644 --- a/mne/viz/backends/_pyvista.py +++ b/mne/viz/backends/_pyvista.py @@ -754,6 +754,7 @@ def _get_camera_direction(focalpoint, position): def _set_3d_view(figure, azimuth, elevation, focalpoint, distance, roll=None): position = np.array(figure.plotter.camera_position[0]) + figure.plotter.reset_camera() if focalpoint is None: focalpoint = np.array(figure.plotter.camera_position[1]) r, theta, phi, fp = _get_camera_direction(focalpoint, position) @@ -931,6 +932,13 @@ def _update_picking_callback(plotter, plotter.picker = picker +def _remove_picking_callback(interactor, picker): + interactor.RemoveObservers(vtk.vtkCommand.RenderEvent) + interactor.RemoveObservers(vtk.vtkCommand.LeftButtonPressEvent) + interactor.RemoveObservers(vtk.vtkCommand.EndInteractionEvent) + picker.RemoveObservers(vtk.vtkCommand.EndPickEvent) + + def _arrow_glyph(grid, factor): glyph = vtk.vtkGlyphSource2D() glyph.SetGlyphTypeToArrow() diff --git a/mne/viz/evoked.py b/mne/viz/evoked.py index f6c2510a103..d1bf0e573a5 100644 --- a/mne/viz/evoked.py +++ b/mne/viz/evoked.py @@ -223,7 +223,7 @@ def _plot_evoked(evoked, picks, exclude, unit, show, ylim, proj, xlim, hline, raise ValueError("If `group_by` is a dict, `axes` must be " "a dict of axes or None.") _validate_if_list_of_axes(list(axes.values())) - remove_xlabels = any([ax.is_last_row() for ax in axes.values()]) + remove_xlabels = any([_is_last_row(ax) for ax in axes.values()]) for sel in group_by: # ... we loop over selections if sel not in axes: raise ValueError(sel + " present in `group_by`, but not " @@ -243,7 +243,7 @@ def _plot_evoked(evoked, picks, exclude, unit, show, ylim, proj, xlim, hline, mask_alpha=mask_alpha, time_unit=time_unit, show_names=show_names, sphere=sphere) - if remove_xlabels and not ax.is_last_row(): + if remove_xlabels and not _is_last_row(ax): ax.set_xticklabels([]) ax.set_xlabel("") ims = [ax.images[0] for ax in axes.values()] @@ -369,6 +369,13 @@ def _plot_evoked(evoked, picks, exclude, unit, show, ylim, proj, xlim, hline, return fig +def _is_last_row(ax): + try: + return ax.get_subplotspec().is_last_row() + except AttributeError: # XXX old mpl + return ax.is_last_row() + + def _plot_lines(data, info, picks, fig, axes, spatial_colors, unit, units, scalings, hline, gfp, types, zorder, xlim, ylim, times, bad_ch_idx, titles, ch_types_used, selectable, psd, diff --git a/mne/viz/tests/test_3d.py b/mne/viz/tests/test_3d.py index 35224d518ed..e367a1a7d25 100644 --- a/mne/viz/tests/test_3d.py +++ b/mne/viz/tests/test_3d.py @@ -121,10 +121,10 @@ def test_plot_sparse_source_estimates(renderer_interactive): stc = SourceEstimate(stc_data, vertices, 1, 1) colormap = 'mne_analyze' - plot_source_estimates(stc, 'sample', colormap=colormap, - background=(1, 1, 0), - subjects_dir=subjects_dir, colorbar=True, - clim='auto') + brain = plot_source_estimates( + stc, 'sample', colormap=colormap, background=(1, 1, 0), + subjects_dir=subjects_dir, colorbar=True, clim='auto') + brain.close() pytest.raises(TypeError, plot_source_estimates, stc, 'sample', figure='foo', hemi='both', clim='auto', subjects_dir=subjects_dir) @@ -592,6 +592,7 @@ def test_plot_source_estimates(renderer_interactive, all_src_types_inv_evoked, ) if pick_ori != 'vector': kwargs['surface'] = 'white' + kwargs['backend'] = renderer_interactive._get_3d_backend() # Mayavi can't handle non-surface if kind != 'surface' and not is_pyvista: with pytest.raises(RuntimeError, match='PyVista'): diff --git a/mne/viz/tests/test_ica.py b/mne/viz/tests/test_ica.py index 355d23e6f79..a10ef0e3e41 100644 --- a/mne/viz/tests/test_ica.py +++ b/mne/viz/tests/test_ica.py @@ -297,6 +297,13 @@ def test_plot_ica_overlay(): plt.close('all') +def _get_geometry(fig): + try: + return fig.axes[0].get_subplotspec().get_geometry() # pragma: no cover + except AttributeError: # MPL < 3.4 (probably) + return fig.axes[0].get_geometry() # pragma: no cover + + @requires_sklearn def test_plot_ica_scores(): """Test plotting of ICA scores.""" @@ -323,20 +330,20 @@ def test_plot_ica_scores(): # check setting number of columns fig = ica.plot_scores([[0.3, 0.2], [0.3, 0.2], [0.3, 0.2]], axhline=[0.1, -0.1]) - assert 2 == fig.get_axes()[0].get_geometry()[1] + assert 2 == _get_geometry(fig)[1] fig = ica.plot_scores([[0.3, 0.2], [0.3, 0.2]], axhline=[0.1, -0.1], n_cols=1) - assert 1 == fig.get_axes()[0].get_geometry()[1] + assert 1 == _get_geometry(fig)[1] # only use 1 column (even though 2 were requested) fig = ica.plot_scores([0.3, 0.2], axhline=[0.1, -0.1], n_cols=2) - assert 1 == fig.get_axes()[0].get_geometry()[1] + assert 1 == _get_geometry(fig)[1] - pytest.raises( - ValueError, - ica.plot_scores, - [0.3, 0.2], axhline=[0.1, -0.1], labels=['one', 'one-too-many']) - pytest.raises(ValueError, ica.plot_scores, [0.2]) + with pytest.raises(ValueError, match='Need as many'): + ica.plot_scores([0.3, 0.2], axhline=[0.1, -0.1], + labels=['one', 'one-too-many']) + with pytest.raises(ValueError, match='The length of'): + ica.plot_scores([0.2]) plt.close('all') diff --git a/mne/viz/topomap.py b/mne/viz/topomap.py index 623d390db01..6ae2c60fd35 100644 --- a/mne/viz/topomap.py +++ b/mne/viz/topomap.py @@ -2550,7 +2550,9 @@ def plot_arrowmap(data, info_from, info_to=None, scale=3e-10, vmin=None, dyy = -dx.data axes.quiver(x, y, dxx, dyy, scale=scale, color='k', lw=1, clip_on=False) axes.figure.canvas.draw_idle() - tight_layout(fig=fig) + with warnings.catch_warnings(record=True): + warnings.simplefilter('ignore') + tight_layout(fig=fig) plt_show(show) return fig diff --git a/tools/circleci_download.sh b/tools/circleci_download.sh new file mode 100755 index 00000000000..b1f3ec14b60 --- /dev/null +++ b/tools/circleci_download.sh @@ -0,0 +1,111 @@ +#!/bin/bash -ef + +if [ "$CIRCLE_BRANCH" == "master" ] || [[ $(cat gitlog.txt) == *"[circle full]"* ]]; then + echo "Doing a full dev build"; + echo html_dev-memory > build.txt; + python -c "import mne; mne.datasets._download_all_example_data()"; + elif [ "$CIRCLE_BRANCH" == "maint/0.21" ]; then + echo "Doing a full stable build"; + echo html_stable-memory > build.txt; + python -c "import mne; mne.datasets._download_all_example_data()"; +else + echo "Doing a partial build"; + if ! git remote -v | grep upstream ; then git remote add upstream git://github.com/mne-tools/mne-python.git; fi + git fetch upstream + FNAMES=$(git diff --name-only $(git merge-base $CIRCLE_BRANCH upstream/master) $CIRCLE_BRANCH); + if [[ $(cat gitlog.txt) == *"[circle front]"* ]]; then + FNAMES="tutorials/source-modeling/plot_mne_dspm_source_localization.py tutorials/machine-learning/plot_receptive_field.py examples/connectivity/plot_mne_inverse_label_connectivity.py tutorials/machine-learning/plot_sensors_decoding.py tutorials/stats-source-space/plot_stats_cluster_spatio_temporal.py tutorials/evoked/plot_20_visualize_evoked.py "${FNAMES}; + python -c "import mne; print(mne.datasets.testing.data_path(update_path=True))"; + fi; + echo FNAMES="$FNAMES"; + for FNAME in $FNAMES; do + if [[ `expr match $FNAME "\(tutorials\|examples\)/.*plot_.*\.py"` ]] ; then + echo "Checking example $FNAME ..."; + PATTERN=`basename $FNAME`"\\|"$PATTERN; + if [[ $(cat $FNAME | grep -x ".*datasets.*sample.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.sample.data_path(update_path=True))"; + fi; + if [[ $(cat $FNAME | grep -x ".*datasets.*fetch_fsaverage.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.fetch_fsaverage(verbose=True))"; + fi; + if [[ $(cat $FNAME | grep -x ".*datasets.*spm_face.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.spm_face.data_path(update_path=True))"; + fi; + if [[ $(cat $FNAME | grep -x ".*datasets.*somato.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.somato.data_path(update_path=True))"; + fi; + if [[ $(cat $FNAME | grep -x ".*datasets.*eegbci.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.eegbci.load_data(1, [3, 6, 10, 14], update_path=True))"; + python -c "import mne; print(mne.datasets.eegbci.load_data(2, [3], update_path=True))"; + python -c "import mne; print(mne.datasets.eegbci.load_data(3, [3], update_path=True))"; + python -c "import mne; print(mne.datasets.eegbci.load_data(4, [3], update_path=True))"; + fi; + if [[ $(cat $FNAME | grep -x ".*datasets.*sleep_physionet.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.sleep_physionet.age.fetch_data([0, 1], recording=[1], update_path=True))"; + fi; + if [[ $(cat $FNAME | grep -x ".*datasets.*hf_sef.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.hf_sef.data_path(update_path=True))"; + fi; + if [[ $(cat $FNAME | grep -x ".*brainstorm.*bst_auditory.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.brainstorm.bst_auditory.data_path(update_path=True))" --accept-brainstorm-license; + fi; + if [[ $(cat $FNAME | grep -x ".*brainstorm.*bst_resting.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.brainstorm.bst_resting.data_path(update_path=True))" --accept-brainstorm-license; + fi; + if [[ $(cat $FNAME | grep -x ".*brainstorm.*bst_raw.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.brainstorm.bst_raw.data_path(update_path=True))" --accept-brainstorm-license; + fi; + if [[ $(cat $FNAME | grep -x ".*brainstorm.*bst_phantom_ctf.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.brainstorm.bst_phantom_ctf.data_path(update_path=True))" --accept-brainstorm-license; + fi; + if [[ $(cat $FNAME | grep -x ".*brainstorm.*bst_phantom_elekta.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.brainstorm.bst_phantom_elekta.data_path(update_path=True))" --accept-brainstorm-license; + fi; + if [[ $(cat $FNAME | grep -x ".*datasets.*hcp_mmp_parcellation.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.sample.data_path(update_path=True))"; + python -c "import mne; print(mne.datasets.fetch_hcp_mmp_parcellation(subjects_dir=mne.datasets.sample.data_path() + '/subjects'))" --accept-hcpmmp-license; + fi; + if [[ $(cat $FNAME | grep -x ".*datasets.*misc.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.misc.data_path(update_path=True))"; + fi; + if [[ $(cat $FNAME | grep -x ".*datasets.*testing.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.testing.data_path(update_path=True))"; + fi; + if [[ $(cat $FNAME | grep -x ".*datasets.*kiloword.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.kiloword.data_path(update_path=True))"; + fi; + if [[ $(cat $FNAME | grep -x ".*datasets.*mtrf.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.mtrf.data_path(update_path=True))"; + fi; + if [[ $(cat $FNAME | grep -x ".*datasets.*fieldtrip_cmc.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.fieldtrip_cmc.data_path(update_path=True))"; + fi; + if [[ $(cat $FNAME | grep -x ".*datasets.*multimodal.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.multimodal.data_path(update_path=True))"; + fi; + if [[ $(cat $FNAME | grep -x ".*datasets.*fnirs_motor.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.fnirs_motor.data_path(update_path=True))"; + fi; + if [[ $(cat $FNAME | grep -x ".*datasets.*opm.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.opm.data_path(update_path=True))"; + fi; + if [[ $(cat $FNAME | grep -x ".*datasets.*phantom_4dbti.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.phantom_4dbti.data_path(update_path=True))"; + fi; + if [[ $(cat $FNAME | grep -x ".*datasets.*limo.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.limo.data_path(subject=1, update_path=True))"; + fi; + if [[ $(cat $FNAME | grep -x ".*datasets.*refmeg_noise.*" | wc -l) -gt 0 ]]; then + python -c "import mne; print(mne.datasets.refmeg_noise.data_path(update_path=True))"; + fi; + fi; + done; + echo PATTERN="$PATTERN"; + if [[ $PATTERN ]]; then + PATTERN="\(${PATTERN::-2}\)"; + echo html_dev-pattern-memory > build.txt; + else + echo html_dev-noplot > build.txt; + fi; +fi; +echo "$PATTERN" > pattern.txt; \ No newline at end of file diff --git a/tutorials/misc/plot_ecog.py b/tutorials/misc/plot_ecog.py index 9f40bf58afb..0566e14fc8d 100644 --- a/tutorials/misc/plot_ecog.py +++ b/tutorials/misc/plot_ecog.py @@ -166,7 +166,7 @@ def animate(i, activity): # create the figure and apply the animation of the # gamma frequency band activity -fig, ax = plt.subplots(figsize=(10, 10)) +fig, ax = plt.subplots(figsize=(5, 5)) ax.imshow(im) ax.set_axis_off() paths = ax.scatter(*xy_pts.T, c=np.zeros(len(xy_pts)), s=200, @@ -176,11 +176,12 @@ def animate(i, activity): size='large') # avoid edge artifacts and decimate, showing just a short chunk -show_power = gamma_power_t[:, 100:150] +sl = slice(100, 150) +show_power = gamma_power_t[:, sl] anim = animation.FuncAnimation(fig, animate, init_func=init, fargs=(show_power,), frames=show_power.shape[1], - interval=200, blit=True) + interval=100, blit=True) ############################################################################### # Alternatively, we can project the sensor data to the nearest locations on @@ -188,13 +189,13 @@ def animate(i, activity): # sphinx_gallery_thumbnail_number = 4 -evoked = mne.EvokedArray(gamma_power_t, raw.info) +evoked = mne.EvokedArray( + gamma_power_t[:, sl], raw.info, tmin=raw.times[sl][0]) stc = mne.stc_near_sensors(evoked, trans, subject, subjects_dir=subjects_dir) clim = dict(kind='value', lims=[vmin * 0.9, vmin, vmax]) brain = stc.plot(surface='pial', hemi='both', initial_time=0.68, colormap='viridis', clim=clim, views='parietal', - subjects_dir=subjects_dir, size=(600, 600)) + subjects_dir=subjects_dir, size=(500, 500)) # You can save a movie like the one on our documentation website with: -# brain.save_movie(time_dilation=20, tmin=0.62, tmax=0.72, -# interpolation='linear', framerate=5, +# brain.save_movie(time_dilation=50, interpolation='linear', framerate=10, # time_viewer=True) diff --git a/tutorials/simulation/plot_dics.py b/tutorials/simulation/plot_dics.py index 3e1eef4650c..387c8460a0c 100644 --- a/tutorials/simulation/plot_dics.py +++ b/tutorials/simulation/plot_dics.py @@ -7,9 +7,9 @@ locations on the cortex. These signals will be sinusoids, so we'll be looking at oscillatory activity (as opposed to evoked activity). -We'll use dynamic imaging of coherent sources (DICS) [1]_ to map out -spectral power along the cortex. Let's see if we can find our two simulated -sources. +We'll use dynamic imaging of coherent sources (DICS) :footcite:`GrossEtAl2001` +to map out spectral power along the cortex. Let's see if we can find our two +simulated sources. """ # Author: Marijn van Vliet # @@ -237,10 +237,11 @@ def coh_signal_gen(): # consensus regarding the best approach. This is why we will demonstrate two # approaches here: # -# 1. The approach as described in [2]_, which first normalizes the forward -# solution and computes a vector beamformer. -# 2. The scalar beamforming approach based on [3]_, which uses weight -# normalization instead of normalizing the forward solution. +# 1. The approach as described in :footcite:`vanVlietEtAl2018`, which first +# normalizes the forward solution and computes a vector beamformer. +# 2. The scalar beamforming approach based on +# :footcite:`SekiharaNagarajan2008`, which uses weight normalization +# instead of normalizing the forward solution. # Estimate the cross-spectral density (CSD) matrix on the trial containing the # signal. @@ -264,22 +265,31 @@ def coh_signal_gen(): power_approach1, f = apply_dics_csd(csd_signal, filters_approach1) power_approach2, f = apply_dics_csd(csd_signal, filters_approach2) -# Plot the DICS power maps for both approaches. -brains = list() -for approach, power in enumerate([power_approach1, power_approach2], 1): - title = 'DICS power map, approach %d' % approach - brain = power.plot('sample', subjects_dir=subjects_dir, hemi='both', - figure=approach + 1, size=600, time_label=title, - title=title) +############################################################################### +# Plot the DICS power maps for both approaches, starting with the first: + +def plot_approach(power, n): + """Plot the results on a brain.""" + title = 'DICS power map, approach %d' % n + brain = power_approach1.plot( + 'sample', subjects_dir=subjects_dir, hemi='both', + size=600, time_label=title, title=title) # Indicate the true locations of the source activity on the plot. brain.add_foci(vertices[0][0], coords_as_verts=True, hemi='lh', color='b') brain.add_foci(vertices[1][0], coords_as_verts=True, hemi='rh', color='b') - # Rotate the view and add a title. brain.show_view(view={'azimuth': 0, 'elevation': 0, 'distance': 550, - 'focalpoint': [0, 0, 0]}) - brains.append(brain) + 'focalpoint': [0, 0, 0]}) + return brain + + +brain1 = plot_approach(power_approach1, 1) + +############################################################################### +# Now the second: + +brain2 = plot_approach(power_approach2, 2) ############################################################################### # Excellent! All methods found our two simulated sources. Of course, with a @@ -290,12 +300,4 @@ def coh_signal_gen(): # # References # ---------- -# .. [1] Gross et al. (2001). Dynamic imaging of coherent sources: Studying -# neural interactions in the human brain. Proceedings of the National -# Academy of Sciences, 98(2), 694-699. -# https://doi.org/10.1073/pnas.98.2.694 -# .. [2] van Vliet, et al. (2018) Analysis of functional connectivity and -# oscillatory power using DICS: from raw MEG data to group-level -# statistics in Python. bioRxiv, 245530. https://doi.org/10.1101/245530 -# .. [3] Sekihara & Nagarajan. Adaptive spatial filters for electromagnetic -# brain imaging (2008) Springer Science & Business Media +# .. footbibliography:: diff --git a/tutorials/source-modeling/plot_dipole_orientations.py b/tutorials/source-modeling/plot_dipole_orientations.py index a49aa3a57d8..4387830516d 100644 --- a/tutorials/source-modeling/plot_dipole_orientations.py +++ b/tutorials/source-modeling/plot_dipole_orientations.py @@ -63,16 +63,14 @@ coord_frame = 'mri' # Plot the cortex -fig = mne.viz.plot_alignment(subject=subject, subjects_dir=subjects_dir, - trans=trans, surfaces='white', - coord_frame=coord_frame, fig=fig) +mne.viz.plot_alignment( + subject=subject, subjects_dir=subjects_dir, trans=trans, surfaces='white', + coord_frame=coord_frame, fig=fig) # Mark the position of the dipoles with small red dots -fig = mne.viz.plot_dipole_locations(dipoles=dipoles, trans=trans, - mode='sphere', subject=subject, - subjects_dir=subjects_dir, - coord_frame=coord_frame, - scale=7e-4, fig=fig) +mne.viz.plot_dipole_locations( + dipoles=dipoles, trans=trans, mode='sphere', subject=subject, + subjects_dir=subjects_dir, coord_frame=coord_frame, scale=7e-4, fig=fig) mne.viz.set_3d_view(figure=fig, azimuth=180, distance=0.25) @@ -96,16 +94,14 @@ fig = mne.viz.create_3d_figure(size=(600, 400)) # Plot the cortex -fig = mne.viz.plot_alignment(subject=subject, subjects_dir=subjects_dir, - trans=trans, - surfaces='white', coord_frame='head', fig=fig) +mne.viz.plot_alignment( + subject=subject, subjects_dir=subjects_dir, trans=trans, + surfaces='white', coord_frame='head', fig=fig) # Show the dipoles as arrows pointing along the surface normal -fig = mne.viz.plot_dipole_locations(dipoles=dipoles, trans=trans, - mode='arrow', subject=subject, - subjects_dir=subjects_dir, - coord_frame='head', - scale=7e-4, fig=fig) +mne.viz.plot_dipole_locations( + dipoles=dipoles, trans=trans, mode='arrow', subject=subject, + subjects_dir=subjects_dir, coord_frame='head', scale=7e-4, fig=fig) mne.viz.set_3d_view(figure=fig, azimuth=180, distance=0.1) @@ -147,14 +143,14 @@ fig = mne.viz.create_3d_figure(size=(600, 400)) # Plot the cortex -fig = mne.viz.plot_alignment(subject=subject, subjects_dir=subjects_dir, - trans=trans, - surfaces='white', coord_frame='head', fig=fig) +mne.viz.plot_alignment( + subject=subject, subjects_dir=subjects_dir, trans=trans, + surfaces='white', coord_frame='head', fig=fig) # Show the three dipoles defined at each location in the source space -fig = mne.viz.plot_alignment(subject=subject, subjects_dir=subjects_dir, - trans=trans, fwd=fwd, - surfaces='white', coord_frame='head', fig=fig) +mne.viz.plot_alignment( + subject=subject, subjects_dir=subjects_dir, trans=trans, fwd=fwd, + surfaces='white', coord_frame='head', fig=fig) mne.viz.set_3d_view(figure=fig, azimuth=180, distance=0.1) diff --git a/tutorials/source-modeling/plot_eeg_mri_coords.py b/tutorials/source-modeling/plot_eeg_mri_coords.py index e4b9f3a9f96..af3a73ed90e 100644 --- a/tutorials/source-modeling/plot_eeg_mri_coords.py +++ b/tutorials/source-modeling/plot_eeg_mri_coords.py @@ -138,8 +138,9 @@ # New we have all of the components we need to compute a forward solution, # but first we should sanity check that everything is well aligned: -plot_alignment(evoked.info, trans=trans, show_axes=True, surfaces='head-dense', - subject='sample', subjects_dir=subjects_dir) +fig = plot_alignment( + evoked.info, trans=trans, show_axes=True, surfaces='head-dense', + subject='sample', subjects_dir=subjects_dir) ############################################################################## # Now we can actually compute the forward: @@ -153,4 +154,4 @@ inv = mne.minimum_norm.make_inverse_operator( evoked.info, fwd, cov, verbose=True) stc = mne.minimum_norm.apply_inverse(evoked, inv) -stc.plot(subjects_dir=subjects_dir, initial_time=0.1) +brain = stc.plot(subjects_dir=subjects_dir, initial_time=0.1) diff --git a/tutorials/source-modeling/plot_eeg_no_mri.py b/tutorials/source-modeling/plot_eeg_no_mri.py index d04fc5da7a8..d6359990675 100644 --- a/tutorials/source-modeling/plot_eeg_no_mri.py +++ b/tutorials/source-modeling/plot_eeg_no_mri.py @@ -78,5 +78,5 @@ # for illustration purposes use fwd to compute the sensitivity map eeg_map = mne.sensitivity_map(fwd, ch_type='eeg', mode='fixed') -eeg_map.plot(time_label='EEG sensitivity', subjects_dir=subjects_dir, - clim=dict(lims=[5, 50, 100])) +brain = eeg_map.plot(time_label='EEG sensitivity', subjects_dir=subjects_dir, + clim=dict(lims=[5, 50, 100])) diff --git a/tutorials/source-modeling/plot_visualize_stc.py b/tutorials/source-modeling/plot_visualize_stc.py index e446b07b781..09c3fa867d8 100644 --- a/tutorials/source-modeling/plot_visualize_stc.py +++ b/tutorials/source-modeling/plot_visualize_stc.py @@ -136,8 +136,7 @@ fname_aseg = op.join(subjects_dir, 'sample', 'mri', 'aparc.a2009s+aseg.mgz') label_names = mne.get_volume_labels_from_aseg(fname_aseg) -label_tc = stc.extract_label_time_course( - fname_aseg, src=src, trans=mri_head_t) +label_tc = stc.extract_label_time_course(fname_aseg, src=src) lidx, tidx = np.unravel_index(np.argmax(label_tc), label_tc.shape) fig, ax = plt.subplots(1)