diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 847740f..8ee19ec 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -8,6 +8,9 @@ on: pull_request: branches: [ "master" ] +env: + BUILD_ENV: ci + permissions: contents: read @@ -25,7 +28,7 @@ jobs: - name: Install dependencies run: | python -m pip install --upgrade pip - pip install -e .[docs] + pip install -e .[docs,plot] - name: Sphinx build run: | diff --git a/.readthedocs.yaml b/.readthedocs.yaml index e8a5093..3c68071 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -24,3 +24,4 @@ python: path: . extra_requirements: - docs + - plot diff --git a/CHANGELOG.md b/CHANGELOG.md index cfe85b0..09f2fda 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,8 +18,8 @@ noise ([#58](https://github.com/ctrltz/meegsim/pull/58)) - Allow specifying standard deviation via a SourceEstimate object ([#67](https://github.com/ctrltz/meegsim/pull/67)) - A method for setting phase-phase coupling by adding noise to the shifted copy of input waveform ([#71](https://github.com/ctrltz/meegsim/pull/71)) - Function to convert the sources to mne.Label ([#73](https://github.com/ctrltz/meegsim/pull/73)) -- Quick list-like access to the simulated sources ([#82](https://github.com/ctrltz/meegsim/pull/82)) -- Control over the amplitude envelope of the coupled waveform ([#87](https://github.com/ctrltz/meegsim/pull/87)) +- Quick dict-like access to the simulated sources ([#82](https://github.com/ctrltz/meegsim/pull/82)) +- Partial control over the amplitude envelope of the coupled waveform: same as input or randomly generated ([#87](https://github.com/ctrltz/meegsim/pull/87)) ### Changed diff --git a/CITATION.cff b/CITATION.cff new file mode 100644 index 0000000..52547cf --- /dev/null +++ b/CITATION.cff @@ -0,0 +1,37 @@ +# This CITATION.cff file was generated with cffinit. +# Visit https://bit.ly/cffinit to generate yours today! + +cff-version: 1.2.0 +title: MEEGsim +message: >- + If you use this software, please cite it using the + metadata from this file. +type: software +authors: + - given-names: Nikolai + family-names: Kapralov + email: kapralov@cbs.mpg.de + affiliation: >- + Max Planck Institute for Human Cognitive and Brain + Sciences, Leipzig, Germany + orcid: 'https://orcid.org/0000-0002-5659-7307' + - given-names: Alina + family-names: Studenova + email: studenova@cbs.mpg.de + affiliation: >- + Max Planck Institute for Human Cognitive and Brain + Sciences, Leipzig, Germany + orcid: 'https://orcid.org/0000-0003-0821-9966' + - given-names: Mina + family-names: Jamshidi Idaji + orcid: 'https://orcid.org/0000-0003-1593-3201' +repository-code: 'https://github.com/ctrltz/meegsim' +url: 'https://meegsim.readthedocs.io/en/stable/' +keywords: + - simulation + - MEG + - EEG + - connectivity +license: BSD-3-Clause +version: 0.0.1 +date-released: '2024-10-31' diff --git a/Makefile b/Makefile index baaf39a..0a0a4e3 100644 --- a/Makefile +++ b/Makefile @@ -7,11 +7,21 @@ SPHINXOPTS ?= SPHINXBUILD ?= sphinx-build SOURCEDIR = docs BUILDDIR = _build +EXAMPLESDIR = docs/auto_examples +AUTODOCDIR = docs/generated # Put it first so that "make" without argument is like "make help". help: @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) +clean: + rm -rf $(BUILDDIR)/* + rm -rf $(EXAMPLESDIR)/* + rm -rf $(AUTODOCDIR)/* + +collect: html + ./collect_example_stubs.sh + show: python -m webbrowser -t "$(BUILDDIR)/html/index.html" diff --git a/README.md b/README.md index af52c2f..e3e2f8e 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,7 @@ # MEEGsim +[![DOI](https://zenodo.org/badge/832185753.svg)](https://doi.org/10.5281/zenodo.15106042) + ## Overview **MEEGsim** is a Python package that provides template waveforms for simulating M/EEG data with known ground truth source activity. In addition, it simplifies the manipulation of relevant simulation parameters (e.g., signal-to-noise ratio and source connectivity). As a result, the users can focus on _what_ to simulate, not on _how_ to implement the simulation. The package is compatible with MNE-Python and re-uses the forward and inverse modeling functionality provided by MNE. diff --git a/collect_example_stubs.sh b/collect_example_stubs.sh new file mode 100755 index 0000000..c3a1fa6 --- /dev/null +++ b/collect_example_stubs.sh @@ -0,0 +1,15 @@ +# NOTE: brain plots cannot be rendered with ReadTheDocs, therefore +# we build the docs locally and add generated brain images in correct places +# to show in the online version + +STATIC_THUMB=docs/_static/example_stubs/thumb +STATIC_IMAGES=docs/_static/example_stubs/images + +# examples/basics/02_plot_brain_configuration.py +cp docs/auto_examples/basics/images/thumb/sphx_glr_02_plot_brain_configuration_thumb.png $STATIC_THUMB +cp docs/auto_examples/basics/images/sphx_glr_02_plot_brain_configuration_001.png $STATIC_IMAGES + +# examples/building_blocks/04_plot_std.py +cp docs/auto_examples/building_blocks/images/thumb/sphx_glr_04_plot_brain_std_thumb.png $STATIC_THUMB +cp docs/auto_examples/building_blocks/images/sphx_glr_04_plot_brain_std_001.png $STATIC_IMAGES +cp docs/auto_examples/building_blocks/images/sphx_glr_04_plot_brain_std_002.png $STATIC_IMAGES diff --git a/docs/_static/example_stubs/images/sphx_glr_02_plot_brain_configuration_001.png b/docs/_static/example_stubs/images/sphx_glr_02_plot_brain_configuration_001.png new file mode 100644 index 0000000..69d0336 Binary files /dev/null and b/docs/_static/example_stubs/images/sphx_glr_02_plot_brain_configuration_001.png differ diff --git a/docs/_static/example_stubs/images/sphx_glr_04_plot_brain_std_001.png b/docs/_static/example_stubs/images/sphx_glr_04_plot_brain_std_001.png new file mode 100644 index 0000000..f570f63 Binary files /dev/null and b/docs/_static/example_stubs/images/sphx_glr_04_plot_brain_std_001.png differ diff --git a/docs/_static/example_stubs/images/sphx_glr_04_plot_brain_std_002.png b/docs/_static/example_stubs/images/sphx_glr_04_plot_brain_std_002.png new file mode 100644 index 0000000..a54cebc Binary files /dev/null and b/docs/_static/example_stubs/images/sphx_glr_04_plot_brain_std_002.png differ diff --git a/docs/_static/example_stubs/thumb/sphx_glr_02_plot_brain_configuration_thumb.png b/docs/_static/example_stubs/thumb/sphx_glr_02_plot_brain_configuration_thumb.png new file mode 100644 index 0000000..73301f0 Binary files /dev/null and b/docs/_static/example_stubs/thumb/sphx_glr_02_plot_brain_configuration_thumb.png differ diff --git a/docs/_static/example_stubs/thumb/sphx_glr_04_plot_brain_std_thumb.png b/docs/_static/example_stubs/thumb/sphx_glr_04_plot_brain_std_thumb.png new file mode 100644 index 0000000..a5d1cbe Binary files /dev/null and b/docs/_static/example_stubs/thumb/sphx_glr_04_plot_brain_std_thumb.png differ diff --git a/docs/_templates/autosummary/class.rst b/docs/_templates/autosummary/class.rst index b3eedfe..0b397d4 100644 --- a/docs/_templates/autosummary/class.rst +++ b/docs/_templates/autosummary/class.rst @@ -6,3 +6,8 @@ :special-members: __contains__,__getitem__,__iter__,__len__,__add__,__sub__,__mul__,__div__,__neg__ :members: :inherited-members: + +.. _sphx_glr_backreferences_{{ fullname }}: + +.. minigallery:: {{ fullname }} + :add-heading: diff --git a/docs/_templates/autosummary/function.rst b/docs/_templates/autosummary/function.rst new file mode 100644 index 0000000..f911483 --- /dev/null +++ b/docs/_templates/autosummary/function.rst @@ -0,0 +1,10 @@ +{{ fullname | escape | underline }} + +.. currentmodule:: {{ module }} + +.. autofunction:: {{ objname }} + +.. _sphx_glr_backref_{{ fullname }}: + +.. minigallery:: {{ fullname }} + :add-heading: diff --git a/docs/changelog/devel.md b/docs/changelog/devel.md index 4046631..0877b09 100644 --- a/docs/changelog/devel.md +++ b/docs/changelog/devel.md @@ -1 +1,23 @@ # [Unreleased] + +## Added + +- A desired level of white noise can be added in sensor space to model measurement +noise ([#58](https://github.com/ctrltz/meegsim/pull/58)) +- A possibility to plot the source configuration ([#59](https://github.com/ctrltz/meegsim/pull/59)) +- Adjustment of global (all signal vs. all noise sources) SNR ([#64](https://github.com/ctrltz/meegsim/pull/64)) +- Adjustment of the standard deviation of source activity ([#66](https://github.com/ctrltz/meegsim/pull/66)) +- Allow specifying standard deviation via a SourceEstimate object ([#67](https://github.com/ctrltz/meegsim/pull/67)) +- A method for setting phase-phase coupling by adding noise to the shifted copy of input waveform ([#71](https://github.com/ctrltz/meegsim/pull/71)) +- Function to convert the sources to mne.Label ([#73](https://github.com/ctrltz/meegsim/pull/73)) +- Quick dict-like access to the simulated sources ([#82](https://github.com/ctrltz/meegsim/pull/82)) +- Partial control over the amplitude envelope of the coupled waveform: same as input or randomly generated ([#87](https://github.com/ctrltz/meegsim/pull/87)) + +## Changed + +- Reworked normalization of source activity: by default, all source time courses are scaled to make their standard deviation equal to 1 nAm ([#66](https://github.com/ctrltz/meegsim/pull/66)) +- Improved performance when adjusting the SNR for a large number of sources ([#68](https://github.com/ctrltz/meegsim/pull/68)) + +## Fixed + +- Fixed a bug that caused different sources to have the same location and/or waveform when random state was explicitly provided ([#76](https://github.com/ctrltz/meegsim/pull/76)) diff --git a/docs/conf.py b/docs/conf.py index 8981281..afeea9d 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -7,6 +7,7 @@ # https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information +import pyvista import os import sys @@ -128,8 +129,9 @@ def linkcode_resolve(domain, info): # Intersphinx -intersphinx_mapping = get_intersphinx_mapping( - packages={"matplotlib", "mne", "numpy", "python"} +intersphinx_mapping = {"neurodsp": ("https://neurodsp-tools.github.io/neurodsp/", None)} +intersphinx_mapping.update( + get_intersphinx_mapping(packages={"matplotlib", "mne", "numpy", "python"}) ) @@ -138,8 +140,30 @@ def linkcode_resolve(domain, info): bibtex_default_style = "unsrt" +# Enabling PyVista scraper +pyvista.BUILDING_GALLERY = True +pyvista.OFF_SCREEN = False + + # Sphinx Gallery sphinx_gallery_conf = { + "backreferences_dir": "generated", + "doc_module": ("meegsim",), "examples_dirs": "../examples", + "filename_pattern": "/\\d{2}_plot_", "gallery_dirs": "auto_examples", + "image_scrapers": ("matplotlib", "pyvista"), + "subsection_order": [ + "../examples/basics", + "../examples/building_blocks", + "../examples/advanced", + ], + "within_subsection_order": "FileNameSortKey", } + +# Disable brain plots for CIs +if os.environ.get("BUILD_ENV", "local") == "ci": + sphinx_gallery_conf["filename_pattern"] = "/\\d{2}_plot_(?!brain)" + +print(os.environ.get("BUILD_ENV", "local")) +print(sphinx_gallery_conf) diff --git a/docs/development/index.rst b/docs/development/index.rst index 99e05ac..d8ce0c3 100644 --- a/docs/development/index.rst +++ b/docs/development/index.rst @@ -5,5 +5,4 @@ Development :maxdepth: 2 setup - whats_new internals/index diff --git a/docs/references.bib b/docs/references.bib index a568821..a5a7b83 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -1,3 +1,14 @@ +@article{HaufeEwald2019, + title = {{A} {S}imulation {F}ramework for {B}enchmarking {EEG}-{B}ased {B}rain {C}onnectivity {E}stimation {M}ethodologies}, + author = {Haufe, Stefan and Ewald, Arne}, + journal = {{B}rain {T}opography}, + journal-short={{B}rain {T}opogr}, + volume = {32}, + pages = {625-642}, + year = 2019, + doi = "10.1007/s10548-016-0498-y" +} + @article{Nikulin2011, title = "A novel method for reliable and fast extraction of neuronal {EEG/MEG} oscillations on the basis of spatio-spectral diff --git a/docs/user_guide/advanced/customization.rst b/docs/user_guide/advanced/customization.rst index 92d51ac..69985c5 100644 --- a/docs/user_guide/advanced/customization.rst +++ b/docs/user_guide/advanced/customization.rst @@ -3,8 +3,9 @@ Customization ============= If the built-in functions for location, waveform or coupling do not satisfy the -needs of your project, you are always welcome use a custom function instead. -In this section, we provide a short description of requirements for such functions: +needs of your project, you are always welcome use a custom function instead (either +created by you or re-used from another toolbox). In this section, we provide a short +description of requirements for such functions: * which arguments should be accepted (but not necessarily used) @@ -13,8 +14,8 @@ In this section, we provide a short description of requirements for such functio In addition, we provide minimal examples of custom functions for each case. .. note:: - When adding sources, we always try to execute the provided functions with 0 - as ``random_state`` to making the debugging a bit easier. + When adding sources, we always call the provided functions with 0 + as ``random_state`` to making the debugging a bit easier in case errors occur. Location ======== @@ -59,7 +60,12 @@ The waveform function should accept: * keyword argument ``random_state`` -The result is expected to be an array with shape ``(n_series, n_times)``. +The result is expected to be an array with shape ``(n_series, n_times)``. Below, +we show two examples: the first relies on an own custom function, the second shows +how to adapt the function from another toolbox to be used with MEEGsim. + +Own custom function +------------------- The function below returns white noise, and it produces different results every time unless ``random_state`` is fixed: @@ -83,6 +89,63 @@ are not required in this case): waveform=my_white_noise ) +Function from another package +----------------------------- + +For this example, we use a +`function `_ +from the NeuroDSP package that allows simulating bursty oscillations (currently not possible with our toolbox). +First, we need to create a wrapper function to adapt the input and output formats to match +the built-in functions of MEEGsim: + +.. code-block:: python + + from neurodsp.sim import sim_bursty_oscillation + from neurodsp.sim.multi import sim_multiple + from meegsim.utils import normalize_variance + + def bursty_osc(n_series, times, **kwargs): + # Convert MEEGsim input to NeuroDSP input + tstep = (times[1] - times[0]) + n_seconds = times.max() + tstep + fs = 1.0 / tstep + + params = dict(n_seconds=n_seconds, fs=fs) + params.update(kwargs) + params.pop("random_state") # is not accepted by NeuroDSP function + + sims = sim_multiple(sim_bursty_oscillation, params, n_sims=n_series) + + return normalize_variance(sims.signals) + +.. note:: + We use ``**kwargs`` in the example above to forward all additional arguments to the + simulation function from the NeuroDSP package. This way, the names and meaning of + each argument remains the same. + +Once adapted, the function can be used similar to other built-in functions when +adding sources: + +.. code-block:: python + + # src should be loaded before + sim = SourceSimulator(src) + + sim.add_point_sources( + location=[(0, 123), (1, 456)], + waveform=bursty_osc, + waveform_params=dict( # NeuroDSP parameters + freq=20, + burst_def='durations', + burst_params={'n_cycles_burst' : 3, 'n_cycles_off' : 3} + ), + ... # snr / std / names + ) + +However, it is important to keep in mind that coupling methods might also need to be +adapted in order to preserve any special features of the simulated time series (e.g., +presence of bursts). + Coupling ======== diff --git a/docs/user_guide/advanced/index.rst b/docs/user_guide/advanced/index.rst index 975d91d..5498a9f 100644 --- a/docs/user_guide/advanced/index.rst +++ b/docs/user_guide/advanced/index.rst @@ -9,4 +9,6 @@ for the toolbox. :maxdepth: 2 patches + standard_deviation + local_snr customization diff --git a/docs/user_guide/advanced/local_snr.rst b/docs/user_guide/advanced/local_snr.rst new file mode 100644 index 0000000..f072536 --- /dev/null +++ b/docs/user_guide/advanced/local_snr.rst @@ -0,0 +1,68 @@ +======================= +Local adjustment of SNR +======================= + +To enable the *local* adjustment of SNR, you need to enable the corresponding ``snr_mode`` +when initializing the ``SourceSimulator`` class: + +.. code-block:: python + + from meegsim.simulate import SourceSimulator + + # src should be loaded or created beforehand + sim = SourceSimulator(src, snr_mode="local") + +For the *local* adjustment of SNR, we calculate the mean variance of each point or patch source +across all sensors and adjust it relative to the mean variance of all noise sources +(mean over sensors but summed over noise sources). The calculation of variance is +performed after filtering both time series (signal and noise) in the frequency +band of interest. + +.. note:: + For the adjustment of SNR, you always need to add noise sources to the + simulation. + +By default, no adjustment of *local* SNR is performed. To enable it, you need to specify +the value of SNR using the ``snr`` argument and provide the limits of the frequency +band in ``snr_params`` when adding the point or patch sources: + +.. code-block:: python + + import numpy as np + + # add some noise sources + sim.add_noise_sources( + location=select_random, + location_params=dict(n=10) + ) + + # now add point sources with adjustment of SNR + sim.add_point_sources( + location=[(0, 123), (0, 456), (1, 789)], + waveform=np.ones((3, 1000)), + snr=5, + snr_params=dict(fmin=8, fmax=12), + ... + ) + +It is also possible to specify SNR for each of the sources separately by providing +one value for each source: + +.. code-block:: python + + import numpy as np + + # add some noise sources first + sim.add_noise_sources( + location=select_random, + location_params=dict(n=10) + ) + + # now add point sources with adjustment of SNR + sim.add_point_sources( + location=[(0, 123), (0, 456), (1, 789)], + waveform=np.ones((3, 1000)), + snr=[1, 2.5, 5], + snr_params=dict(fmin=8, fmax=12), + ... + ) diff --git a/docs/user_guide/advanced/patches.rst b/docs/user_guide/advanced/patches.rst index a454df8..9d201fb 100644 --- a/docs/user_guide/advanced/patches.rst +++ b/docs/user_guide/advanced/patches.rst @@ -12,7 +12,8 @@ Specifying all vertices belonging to the patch (default) By default, you are expected to provide indices of all vertices that belong to the patch in the second element of ``location`` tuples. -For example, below we define a patch that contains three vertices: +For example, below we define a patch in a left hemisphere (in the case of a surface +source space) that contains three vertices: .. code-block:: python @@ -21,8 +22,8 @@ For example, below we define a patch that contains three vertices: ... ) -And here we define two patches (with 3 and 4 vertices, respectively) in -the same call: +And here we define two patches (with 3 and 4 vertices in left and right hemispheres, +respectively) in the same call: .. code-block:: python @@ -48,7 +49,7 @@ Example 1 - single patch: sim.add_patch_sources( location=[(0, 123)], - extent=15 # in mm + extents=15 # in mm ) Example 2 - several patches of different size: @@ -57,5 +58,5 @@ Example 2 - several patches of different size: sim.add_patch_sources( location=[(0, 123), (1, 456)], - extent=[15, 30] # in mm + extents=[15, 30] # in mm ) diff --git a/docs/user_guide/advanced/standard_deviation.rst b/docs/user_guide/advanced/standard_deviation.rst new file mode 100644 index 0000000..1fe6a33 --- /dev/null +++ b/docs/user_guide/advanced/standard_deviation.rst @@ -0,0 +1,37 @@ +===================================== +Standard deviation of source activity +===================================== + +By default, all generated waveforms (but not the ones explicitly provided by the user) +are normalized to have a standard deviation (SD) of 1 nAm. If necessary, this value can be +modified using the ``base_std`` argument when initializing the ``SourceSimulator``: + +.. code-block:: python + + from meegsim.simulate import SourceSimulator + + sim = SourceSimulator(src, base_std=1e-8) # 10 nAm + +Normalization also means that the SD of activity will be the same for +all added sources by default. However, in real data this does not always hold. +For example, alpha power is usually higher in parieto-occipital areas. Therefore, we +also provide an option to configure the SD of activity for each +source when adding them: + +.. code-block:: python + + sim.add_point_sources( + location=..., + waveform=..., + std=5, # 5 times stronger than by default + ) + +.. note:: + :doc:`Local adjustment of SNR ` will have priority + over the provided value of SD. The adjustment of the SD is mainly designed to be + used in combination with the global adjustment of SNR (the default option). + +The SD can also be set for all candidate source locations at once +using an instance of :py:class:`mne.SourceEstimate`. If the source locations are +generated randomly, the SD will be adjusted according to the +actual location of the source. diff --git a/docs/user_guide/get_started/configuration.rst b/docs/user_guide/get_started/configuration.rst index 4238381..37ac40f 100644 --- a/docs/user_guide/get_started/configuration.rst +++ b/docs/user_guide/get_started/configuration.rst @@ -25,7 +25,7 @@ The result of this function call is a :class:`SourceConfiguration` object that contains all simulated sources and their waveforms (with desired SNR and coupling). Now you can use the :meth:`SourceConfiguration.to_stc()` and -:meth:`SourceConfiguration.to_raw()` to obtain source time courses and +:meth:`SourceConfiguration.to_raw()` methods to obtain source time courses and sensor-space data, respectively. The projection to sensor space requires a forward model (:class:`mne.Forward`) and an :class:`mne.Info` object describing the sensor layout: @@ -36,6 +36,41 @@ the sensor layout: raw = sc.to_raw(fwd, info) +Global adjustment of the SNR +---------------------------- + +If necessary, the global SNR can be adjusted when simulating the data. To achieve this, +you need to specify the desired value of SNR (``snr_global``) and the limits of +the frequency band that should be considered during the adjustment (``fmin`` and +``fmax`` keys in the ``snr_params`` dictionary). As a result, the mean sensor-space +variance of all point and patch sources will be adjusted relative to the mean variance +of all noise sources to achieve the desired SNR: + +.. code-block:: python + + sfreq = 250 # in Hz + duration = 30 # in seconds + sc = sim.simulate( + sfreq, + duration, + fwd=fwd, # Forward model is required for the adjustment of SNR + snr_global=3.0, + snr_params=dict(fmin=8, fmax=12) + ) + +Adding sensor noise +------------------- + +When projecting the simulated source activity to sensor space, it is also possible +to add a desired amount of sensor noise (currently modeled as white noise): + +.. code-block:: python + + raw = sc.to_raw(fwd, info, sensor_noise_level=0.01) + +The ``sensor_noise_level`` controls the ratio of noise to total power in sensor space. +In the example above, 1% of total sensor-space power will originate from noise. + Reproducibility =============== diff --git a/docs/user_guide/get_started/example.rst b/docs/user_guide/get_started/example.rst deleted file mode 100644 index 808ce85..0000000 --- a/docs/user_guide/get_started/example.rst +++ /dev/null @@ -1,56 +0,0 @@ -A minimal example -================= - -Below you can find an example script that contains all the ideas showcased -in the previous sections of the overview. It may serve as a good starting point -for your own simulation. - -.. code-block:: python - - import numpy as np - - from meegsim.coupling import ppc_von_mises - from meegsim.location import select_random - from meegsim.simulate import SourceSimulator - from meegsim.waveform import narrowband_oscillation - - - # Here you need to load the prerequisites: fwd, src, and info - - # Simulation parameters - sfreq = 250 - duration = 120 - - # Initialize - sim = SourceSimulator(src) - - # Add 500 noise sources with random locations - sim.add_noise_sources( - location=select_random, - location_params=dict(n=500) - ) - - # Add two point sources with fixed locations - # (vertex indices are chosen arbitrarily) - sim.add_point_sources( - location=[(0, 123), (1, 456)], - waveform=narrowband_oscillation, - waveform_params=dict(fmin=8, fmax=12), - snr=[2, 5], - snr_params=dict(fmin=8, fmax=12), - names=['s1', 's2'] - ) - - # Set the coupling between point sources - sim.set_coupling( - ('s1', 's2'), - method=ppc_von_mises, - kappa=1, phase_lag=np.pi/2, - fmin=8, fmax=12 - ) - - # Obtain the data - sc = sim.simulate(sfreq, duration, fwd, random_state=0) - - stc = sc.to_stc() - raw = sc.to_raw(fwd, info) diff --git a/docs/user_guide/get_started/index.rst b/docs/user_guide/get_started/index.rst index 9150a6e..3f08578 100644 --- a/docs/user_guide/get_started/index.rst +++ b/docs/user_guide/get_started/index.rst @@ -2,14 +2,16 @@ Getting started with MEEGsim ============================ -Welcome to MEEGsim! In this overview, you will learn how to simulate an M/EEG dataset -with our toolbox. In particular, we will cover the following aspects: +Welcome to MEEGsim! In this overview tutorial, you will learn how to simulate an +M/EEG dataset using building blocks provided by our toolbox. In particular, we will +cover the following aspects: * prerequisites for the simulation and general workflow * different types of sources and how to add them to the simulation -* how to set up location and waveform of any source -* how to adjust the signal-to-noise ratio (SNR) of added sources -* how to set up phase coupling between waveforms of source activity +* how to set up location and waveform of a source +* how to adjust the signal-to-noise ratio (SNR) and/or standard deviation (SD) of activity for added sources +* how to set up phase coupling between waveforms of different sources +* how to obtain sensor-space data for the simulated source configuration .. toctree:: :maxdepth: 2 @@ -17,4 +19,4 @@ with our toolbox. In particular, we will cover the following aspects: workflow simulation configuration - example + next_steps diff --git a/docs/user_guide/get_started/next_steps.rst b/docs/user_guide/get_started/next_steps.rst new file mode 100644 index 0000000..c9ba236 --- /dev/null +++ b/docs/user_guide/get_started/next_steps.rst @@ -0,0 +1,29 @@ +========== +Next steps +========== + +That was it for the general introduction to MEEGsim! Follow the links from the list +below to learn more about the functionality of the package: + + * :doc:`Advanced topics `: visit this page to learn more + about features that were left of the introductory tutorial: patch sources, local + adjustment of SNR, etc. + + * :doc:`Examples <../../auto_examples/index>`: here you can find various examples + that illustrate the nuances of built-in components + (e.g., :doc:`coupling parameters <../../auto_examples/building_blocks/01_plot_coupling>`), + provide a :doc:`backbone <../../auto_examples/basics/00_plot_minimal_example>` + for your own simulations or help with debugging (e.g., + :doc:`plotting <../../auto_examples/basics/02_plot_brain_configuration>` the simulated configuration). + + * Obviously, the toolbox does not cover the needs of every possible simulation. + Visit the :doc:`Customization ` page to learn how to extend + the provided functionality to fit your use case, or how to integrate input from + other packages providing methods for simulating data, such as + `NeuroDSP `_. + + * If you ran into problems while trying to use our package, feel free to describe your + problem in a `new discussion `_ on GitHub. + + * Found a bug? Have an idea for a new feature? Feel free to create an issue + `here `_. diff --git a/docs/user_guide/get_started/simulation.rst b/docs/user_guide/get_started/simulation.rst index aebdef0..dc480e8 100644 --- a/docs/user_guide/get_started/simulation.rst +++ b/docs/user_guide/get_started/simulation.rst @@ -29,15 +29,15 @@ using different methods: * Point-like sources (:meth:`SourceSimulator.add_point_sources`): the time course of activity is assigned to one vertex. -* Cortical patches (:meth:`SourceSimulator.add_patch_sources`): the time course of +* Patch-like sources (:meth:`SourceSimulator.add_patch_sources`): the time course of activity is assigned to a set of vertices. .. note:: Currently, all vertices belonging to the patch have identical waveforms. * Noise sources (:meth:`SourceSimulator.add_noise_sources`): these sources are - automatically considered as `noise` when adjusting the SNR, while point-like - or patch sources are treated as `signal`. + automatically considered as `noise` when adjusting the SNR, while point- + or patch-like sources are treated as `signal`. .. note:: Currently, only point-like noise sources are supported. @@ -55,10 +55,10 @@ object: sim.add_noise_sources(...) The parameters in brackets can be used to configure the location of the source, -the waveform of source activity, the SNR and, for patch sources only, the extent -of the source. Below we discuss the parameters relevant for all sources in more -detail. For more information about patch sources, please visit -:doc:`this page `. +the waveform of source activity and its standard deviation, the sensor-space SNR and, +for patch sources only, the extent of the source. Below we discuss the parameters +relevant for all sources in more detail. For more information about patch sources, +please visit :doc:`this page `. Source location --------------- @@ -73,7 +73,7 @@ In the first case, you need to provide a list of 2-element tuples, where the first element specifies the index of the source space, while the second contains the index of the vertex (`vertno`). If we consider the typical surface source space with two hemispheres, the call below will add two sources to the left -hemisphere (123 and 456) and one source to the right hemisphere (789). +hemisphere (vertices 123 and 456) and one source to the right hemisphere (vertex 789). .. code-block:: python @@ -90,8 +90,8 @@ hemisphere (123 and 456) and one source to the right hemisphere (789). In the second case, you need to provide a function that accepts ``src`` as the first argument and returns the indices of selected vertices in the same format as -described above. If the function relies on additional arguments, these can be -provided in the ``location_params`` argument. In the example below, +described above. If the function relies on additional arguments, their values can be +provided in the ``location_params`` dictionary. In the example below, we use a built-in function :meth:`select_random` to place 10 sources in random locations: @@ -105,8 +105,9 @@ locations: ... ) -As one of the benefits of specifying location as a function, the simulated source -configurations may differ between subsequent calls of the same ``sim`` object, +If location is specified as a function, the simulated source +configurations may differ between subsequent calls of the +:meth:`~meegsim.simulate.SourceSimulator.simulate()` method, simplifying the simulation of multiple datasets with the same underlying idea. .. note:: @@ -121,7 +122,7 @@ Similar to location, activity waveforms can also be defined manually or through a function. In the first case, the toolbox expects an array with one time series per added source, and the number of waveforms should match the number of defined locations. We can now add some activity to the point sources described above -(here we use constant time series for simplicity): +(using constant time series for simplicity): .. code-block:: python @@ -169,65 +170,30 @@ Signal-to-noise ratio (SNR) --------------------------- By default, the waveforms provided by the user are saved as is, while the built-in -waveforms are normalized to have the same variance. In practice, it is often useful -to set up a specific SNR for each source, for example, to test how different analysis -methods perform depending on the SNR of target sources. - -We use the approach from :cite:p:`Nikulin2011` for adjusting the sensor-space SNR -of sources. Namely, we calculate the mean variance of each point or patch source -across all sensors and adjust it relative to the mean variance of all noise sources -(mean over sensors but summed over noise sources). The calculation of variance is -performed after filtering both time series (signal and noise) in the frequency -band of interest. +waveforms are normalized to have the same standard deviation (1 nAm by default). +In practice, it is often useful to set up a specific SNR for each source separately +or for all sources combined. This is relevant, for example, when testing how different +analysis methods perform depending on the SNR of target sources. + +We provide two approaches for adjusting the sensor-space SNR of added sources. The +first approach was used in :cite:`Nikulin2011`, where the SNR was adjusted separately for +each point or patch source relative to the total power of all noise sources. +The second approach was introduced in :cite:`HaufeEwald2019`: in this case, the total power +of all point and patch sources is adjusted relative to the total power of all noise sources +to obtain the desired SNR. We refer to these approaches as *local* and *global* +adjustment of the SNR, respectively. .. note:: - For the adjustment of SNR, you always need to add noise sources to the - simulation. - -By default, no adjustment of SNR is performed. To enable it, you need to specify -the value of SNR using the ``snr`` argument and provide the limits of the frequency -band in ``snr_params`` when adding the point or patch sources: - -.. code-block:: python - - import numpy as np - - # add some noise sources first - sim.add_noise_sources( - location=select_random, - location_params=dict(n=10) - ) - - # now add point sources with adjustment of SNR - sim.add_point_sources( - location=[(0, 123), (0, 456), (1, 789)], - waveform=np.ones((3, 1000)), - snr=5, - snr_params=dict(fmin=8, fmax=12), - ... - ) + By default, the *global* adjustment of SNR is performed. -It is also possible to specify SNR for each of the sources separately by providing -one value for each source: +The *global* adjustment of the SNR is discussed later in the tutorial as it is +performed only at the final simulation step. +The *local* adjustment of the SNR can be performed when adding the sources to +the simulation and is described :doc:`here `. -.. code-block:: python - - import numpy as np - - # add some noise sources first - sim.add_noise_sources( - location=select_random, - location_params=dict(n=10) - ) - - # now add point sources with adjustment of SNR - sim.add_point_sources( - location=[(0, 123), (0, 456), (1, 789)], - waveform=np.ones((3, 1000)), - snr=[1, 2.5, 5], - snr_params=dict(fmin=8, fmax=12), - ... - ) +.. note:: + For any adjustment of the SNR, you need to add noise sources to the + simulation. Configuring coupling between sources ==================================== @@ -235,7 +201,7 @@ Configuring coupling between sources With the toolbox, we aim to provide a convenient interface for the generation of source waveforms with desired coupling. To set the coupling between sources, you only need to specify the names of sources to be coupled and the coupling -parameters. The waveforms will be then generated automatically to obtain the +parameters. The waveforms will then be generated automatically to obtain the desired coupling. Source names @@ -272,10 +238,10 @@ coupling parameters as shown below: .. code-block:: python - from meegsim.coupling import constant_phase_shift + from meegsim.coupling import ppc_constant_phase_shift sim.set_coupling(('source', 'sink'), - method=constant_phase_shift, phase_lag=np.pi/3 + method=ppc_constant_phase_shift, phase_lag=np.pi/3 ) .. currentmodule:: meegsim.coupling @@ -284,10 +250,10 @@ coupling parameters as shown below: Find more details on the parameters of built-in coupling methods :doc:`here `. -In the example above, we used the :meth:`constant_phase_shift` coupling method. +In the example above, we used the :meth:`ppc_constant_phase_shift` coupling method. As the name suggests, it generates time series with a constant phase shift relative to each other. To have more control over the strength of coupling, -you can use the :meth:`ppc_von_mises` method that set probabilistic phase shifts +you can use the :meth:`ppc_von_mises` method that sets probabilistic phase shifts based on the von Mises distribution. For more details about both approaches, see chapter 3 of :cite:`JamshidiIdaji2022_PhDThesis`. diff --git a/docs/user_guide/how_to_cite.rst b/docs/user_guide/how_to_cite.rst new file mode 100644 index 0000000..434e8ed --- /dev/null +++ b/docs/user_guide/how_to_cite.rst @@ -0,0 +1,46 @@ +How to cite +=========== + +If you used MEEGsim in your project, please consider citing the software using the +DOI provided by Zenodo: https://doi.org/10.5281/zenodo.15106042 + +.. code-block:: + + @software{MEEGsim, + author = {Kapralov, Nikolai and Studenova, Alina and Jamshidi Idaji, Mina}, + license = {BSD-3-Clause}, + title = {{MEEGsim}}, + doi = {https://doi.org/10.5281/zenodo.15106042}, + version = {0.0.2}, + year = {2025} + } + +In addition, please `cite MNE-Python `_ +since MEEGsim re-uses lots of important functionality +provided by this toolbox. + +Finally, check the section below to find reference papers for methods behind MEEGsim. + +Methods +------- + +Below we also provide reference papers for some of the approaches and functions provided by MEEGsim. +By citing these references, you can guide the reader to the description of underlying methods: + ++------------------------------------------------------+------------------------------------------+ +| **Approach / Function** | **Reference** | ++------------------------------------------------------+------------------------------------------+ +| Global adjustment of the SNR | :footcite:`HaufeEwald2019` | ++------------------------------------------------------+------------------------------------------+ +| Local adjustment of the SNR | :footcite:`Nikulin2011` | ++------------------------------------------------------+------------------------------------------+ +| :py:func:`meegsim.coupling.ppc_constant_phase_shift` | :footcite:`JamshidiIdaji2022_PhDThesis` | ++------------------------------------------------------+------------------------------------------+ +| :py:func:`meegsim.coupling.ppc_von_mises` | :footcite:`JamshidiIdaji2022_PhDThesis` | ++------------------------------------------------------+------------------------------------------+ + + +References +---------- + +.. footbibliography:: diff --git a/docs/user_guide/index.rst b/docs/user_guide/index.rst index 50472c0..3a23ccd 100644 --- a/docs/user_guide/index.rst +++ b/docs/user_guide/index.rst @@ -6,5 +6,7 @@ User Guide install get_started/index - ../auto_examples/index advanced/index + ../auto_examples/index + whats_new + how_to_cite diff --git a/docs/development/whats_new.rst b/docs/user_guide/whats_new.rst similarity index 100% rename from docs/development/whats_new.rst rename to docs/user_guide/whats_new.rst diff --git a/examples/GALLERY_HEADER.rst b/examples/GALLERY_HEADER.rst index 78fe4b4..952b599 100644 --- a/examples/GALLERY_HEADER.rst +++ b/examples/GALLERY_HEADER.rst @@ -1,4 +1,6 @@ +======== Examples ======== -This gallery illustrates the bulding blocks provided by meegsim and highlights the practical nuances. +In this gallery, you will find examples showcasing the functionality provided +by MEEGsim and how to use it for your own simulations. diff --git a/examples/advanced/00_plot_customization.py b/examples/advanced/00_plot_customization.py new file mode 100644 index 0000000..3046c00 --- /dev/null +++ b/examples/advanced/00_plot_customization.py @@ -0,0 +1,101 @@ +""" +Using functions from other packages +----------------------------------- + +This example shows how to adapt functions from other packages to +use them with MEEGsim. +""" + +import matplotlib.pyplot as plt +import mne +import numpy as np + +from neurodsp.sim import sim_bursty_oscillation +from neurodsp.sim.multi import sim_multiple +from mne.datasets import sample + +from meegsim.location import select_random +from meegsim.simulate import SourceSimulator +from meegsim.utils import normalize_variance + +# %% +# For this example, we use the :func:`~neurodsp.sim.sim_bursty_oscillation` function +# from the `NeuroDSP `_ package. To make +# the function compatible with MEEGsim, we need to wrap it in another function and +# adapt the input and output parameters: + + +def bursty_osc(n_series, times, **kwargs): + # Convert MEEGsim input to NeuroDSP input + tstep = times[1] - times[0] + n_seconds = times.max() + tstep + fs = 1.0 / tstep + + params = dict(n_seconds=n_seconds, fs=fs) + params.update(kwargs) + + # Random state is not accepted by NeuroDSP function, use it to set + # random starting phase for some variability + seed = params.pop("random_state") + phase = np.random.default_rng(seed).random() + params["phase"] = phase + + sims = sim_multiple(sim_bursty_oscillation, params, n_sims=n_series) + + return normalize_variance(sims.signals) + + +# %% +# Load the source space and a head model for EEG channels: + +# Paths +data_path = sample.data_path() / "MEG" / "sample" +fwd_path = data_path / "sample_audvis-meg-eeg-oct-6-fwd.fif" +raw_path = data_path / "sample_audvis_raw.fif" + +# Load the prerequisites: fwd, src, and info +fwd = mne.read_forward_solution(fwd_path) +raw = mne.io.read_raw(raw_path) +src = fwd["src"] +info = raw.info + +# Pick EEG channels only +eeg_idx = mne.pick_types(info, eeg=True) +info_eeg = mne.pick_info(info, eeg_idx) +fwd_eeg = fwd.pick_channels(info_eeg.ch_names) + +# %% +# Now we can use the adapted function in a simulation, providing the required +# parameters via the ``waveform_params`` dictionary: + +sim = SourceSimulator(src) + +sim.add_point_sources( + location=select_random, + location_params=dict(n=2), + waveform=bursty_osc, + waveform_params=dict( # NeuroDSP parameters + freq=20, + burst_def="durations", + burst_params={"n_cycles_burst": 3, "n_cycles_off": 3}, + ), + names=["lh", "rh"], +) + +sc = sim.simulate(sfreq=250, duration=60, fwd=fwd, random_state=123) + +# %% +# We can check that the simulation of a bursty oscillation actually worked: + +n_samples_to_plot = 500 +fig, ax = plt.subplots() +ax.plot(sc.times[:n_samples_to_plot], sc["lh"].waveform[:n_samples_to_plot]) +ax.set_xlabel("Time (s)") +ax.set_ylabel("Amplitude (nAm)") + +# %% +# +# .. note:: +# Coupling methods provided by MEEGsim might not preserve specific properties +# of time courses generated via functions from other packages. (e.g., presence +# of bursts in this example). diff --git a/examples/advanced/GALLERY_HEADER.rst b/examples/advanced/GALLERY_HEADER.rst new file mode 100644 index 0000000..07357b7 --- /dev/null +++ b/examples/advanced/GALLERY_HEADER.rst @@ -0,0 +1,4 @@ +Advanced +======== + +These examples illustrate advanced features provided by MEEGsim. diff --git a/examples/basics/00_plot_minimal_example.py b/examples/basics/00_plot_minimal_example.py new file mode 100644 index 0000000..04d3213 --- /dev/null +++ b/examples/basics/00_plot_minimal_example.py @@ -0,0 +1,112 @@ +""" +=============== +Minimal example +=============== + +Below you can find an example script that contains all the ideas showcased +in the Getting Started tutorial. It may serve as a good starting point +for your own simulation. +""" + +import numpy as np +import mne + +from mne.datasets import sample + +from meegsim.coupling import ppc_von_mises +from meegsim.location import select_random +from meegsim.simulate import SourceSimulator +from meegsim.waveform import narrowband_oscillation + +# %% +# First, we need to load all the prerequisites for our simulation: +# +# * ``src`` - the :class:`mne.SourceSpaces` object that describes all candidate +# source locations +# * ``fwd`` - the :class:`mne.Forward` object that describes the forward model +# * ``info`` - the :class:`mne.Info` object that describes the channel layout +# +# In this simulation, we only use the EEG channels. + +# Paths +data_path = sample.data_path() / "MEG" / "sample" +fwd_path = data_path / "sample_audvis-meg-eeg-oct-6-fwd.fif" +raw_path = data_path / "sample_audvis_raw.fif" + +# Load the prerequisites: fwd, src, and info +fwd = mne.read_forward_solution(fwd_path) +fwd = mne.convert_forward_solution(fwd, force_fixed=True) +raw = mne.io.read_raw(raw_path) +src = fwd["src"] +info = raw.info + +# Pick EEG channels only +eeg_idx = mne.pick_types(info, eeg=True) +info_eeg = mne.pick_info(info, eeg_idx) +fwd_eeg = fwd.pick_channels(info_eeg.ch_names) + +# %% +# Below we define the parameters of the simulation and the simulation itself. +# In this case, we place 50 noise (1/f) sources randomly and add two phase-coupled +# sources of alpha (8-12 Hz) activity in arbitrary locations: + +# Simulation parameters +sfreq = 100 # in Hz +duration = 60 # in seconds + +# Initialize +sim = SourceSimulator(src) + +# Add 50 noise sources with random locations +sim.add_noise_sources(location=select_random, location_params=dict(n=50)) + +# Add two point sources with fixed locations, vertex indices are chosen +# arbitrarily to have one source in each hemisphere +lh_vertno = src[0]["vertno"][0] +rh_vertno = src[1]["vertno"][0] +sim.add_point_sources( + location=[(0, lh_vertno), (1, rh_vertno)], + waveform=narrowband_oscillation, + waveform_params=dict(fmin=8, fmax=12), + names=["s1", "s2"], +) + +# Set the coupling between point sources +sim.set_coupling( + ("s1", "s2"), method=ppc_von_mises, kappa=1, phase_lag=np.pi / 2, fmin=8, fmax=12 +) + +# %% +# Now let's simulate the configuration with a desired level of global SNR: + +sc = sim.simulate( + sfreq, + duration, + fwd=fwd_eeg, + snr_global=3, + snr_params=dict(fmin=8, fmax=12), + random_state=0, +) + +# %% +# We can double-check that the defined sources were successfully added +# (for noise sources, we only print the total count): + +print(f"Point sources: {sc._sources}") +print(f"The number of noise sources: {len(sc._noise_sources)}") + +# %% +# Finally, we can obtain the simulated source activity and/or project it to +# sensor space while adding 1% of sensor noise: + +stc = sc.to_stc() +raw = sc.to_raw(fwd_eeg, info_eeg, sensor_noise_level=0.01) + +# %% +# We can now plot the power spectra of simulated sensor-space data to verify +# that it has a mixture of 1/f and alpha activity as defined in the simulation: + +spec = raw.compute_psd( + method="welch", n_fft=2 * sfreq, n_overlap=sfreq, n_per_seg=2 * sfreq +) +spec.plot() diff --git a/examples/basics/01_plot_random_state.py b/examples/basics/01_plot_random_state.py new file mode 100644 index 0000000..3b5c787 --- /dev/null +++ b/examples/basics/01_plot_random_state.py @@ -0,0 +1,101 @@ +""" +================================ +Random state and reproducibility +================================ + +This example showcases how reproducible source configurations +can be achieved by fixing the random state when simulating. +""" + +import matplotlib.pyplot as plt +import mne +import numpy as np + +from mne.datasets import sample + +from meegsim.coupling import ppc_shifted_copy_with_noise +from meegsim.location import select_random +from meegsim.waveform import narrowband_oscillation +from meegsim.simulate import SourceSimulator + + +# %% +# First, we load the head model and associated source space: + +# Paths +data_path = sample.data_path() / "MEG" / "sample" +fwd_path = data_path / "sample_audvis-meg-eeg-oct-6-fwd.fif" +raw_path = data_path / "sample_audvis_raw.fif" + +# Load the prerequisites: fwd, src, and info +fwd = mne.read_forward_solution(fwd_path) +fwd = mne.convert_forward_solution(fwd, force_fixed=True) +raw = mne.io.read_raw(raw_path) +src = fwd["src"] +info = raw.info + +# Pick EEG channels only +eeg_idx = mne.pick_types(info, eeg=True) +info_eeg = mne.pick_info(info, eeg_idx) +fwd_eeg = fwd.pick_channels(info_eeg.ch_names) + +# %% +# For this example, we place some sources in random locations and set up +# one connectivity edge. Unless the ``random_state`` is fixed, all randomly +# generated components (``location``, ``waveform``, ``coupling``) will differ +# between simulated configurations. With a fixed ``random_state``, the results become +# reproducible. + +sim = SourceSimulator(src) + +# Select some vertices randomly +n_sources = 3 +sim.add_point_sources( + location=select_random, + waveform=narrowband_oscillation, + location_params=dict(n=n_sources), + waveform_params=dict(fmin=8, fmax=12), + names=[str(i) for i in range(n_sources)], +) +sim.set_coupling( + ("0", "1"), + method=ppc_shifted_copy_with_noise, + phase_lag=np.pi / 2, + coh=0.5, + fmin=8, + fmax=12, + band_limited=False, +) + +# %% +# We simulate three configurations, the first and the last one of them +# have the same ``random_state``: + +sfreq = 250 +duration = 60 +sc1 = sim.simulate(sfreq=sfreq, duration=duration, random_state=123) +sc2 = sim.simulate(sfreq=sfreq, duration=duration, random_state=456) +sc3 = sim.simulate(sfreq=sfreq, duration=duration, random_state=123) + +# %% +# First, we can check the locations (``vertno``) of the simulated sources: + +for i, sc in enumerate([sc1, sc2, sc3]): + print(f"Configuration {i+1}: {[int(s.vertno) for s in sc._sources.values()]}") + +# %% +# For the source with name ``"1"``, we additionally plot the waveform in +# all three configurations: + +n_samples_to_plot = 1000 +fig, axes = plt.subplots(nrows=3, figsize=(8, 6), layout="constrained") +for i, (ax, sc) in enumerate(zip(axes, [sc1, sc2, sc3])): + waveform = np.squeeze(sc["1"].waveform) + ax.plot(sc.times[:n_samples_to_plot], waveform[:n_samples_to_plot]) + ax.set_xlabel("Time (s)") + ax.set_ylabel("Amplitude (nAm)") + ax.set_title(f"Configuration {i+1} | random_state={sc.random_state}") + +# %% +# As expected, both locations and waveforms of the simulated sources are the +# same for configurations 1 and 3 but different for configuration 2. diff --git a/examples/basics/02_plot_brain_configuration.py b/examples/basics/02_plot_brain_configuration.py new file mode 100644 index 0000000..6fed9ca --- /dev/null +++ b/examples/basics/02_plot_brain_configuration.py @@ -0,0 +1,109 @@ +""" +========================== +Plotting the configuration +========================== + +This example illustrates how to plot the simulated source configuration. +""" + +# sphinx_gallery_thumbnail_path = '_static/example_stubs/thumb/sphx_glr_02_plot_brain_configuration_thumb.png' + +import mne + +from mne.datasets import sample + +from meegsim.location import select_random +from meegsim.simulate import SourceSimulator +from meegsim.waveform import narrowband_oscillation + +# %% +# First, we load all the prerequisites for our simulation and restrict to the EEG +# channels only + +# Paths +subjects_dir = sample.data_path() / "subjects" +data_path = sample.data_path() / "MEG" / "sample" +fwd_path = data_path / "sample_audvis-meg-eeg-oct-6-fwd.fif" +raw_path = data_path / "sample_audvis_raw.fif" + +print(subjects_dir) + +# Load the prerequisites: fwd, src, and info +fwd = mne.read_forward_solution(fwd_path) +raw = mne.io.read_raw(raw_path) +src = fwd["src"] +info = raw.info + +# Pick EEG channels only +eeg_idx = mne.pick_types(info, eeg=True) +info_eeg = mne.pick_info(info, eeg_idx) +fwd_eeg = fwd.pick_channels(info_eeg.ch_names) + +# %% +# Below we define the simulation itself. In this case, we place 50 noise (1/f) +# sources and add a couple of point and patch sources for demonstration purposes. +# All sources are placed randomly. + +# Simulation parameters +sfreq = 100 # in Hz +duration = 60 # in seconds + +# Initialize +sim = SourceSimulator(src) + +# Add 50 noise sources with random locations +sim.add_noise_sources(location=select_random, location_params=dict(n=50)) + +# Add point sources +sim.add_point_sources( + location=select_random, + location_params=dict(n=3), + waveform=narrowband_oscillation, + waveform_params=dict(fmin=8, fmax=12), +) + +# Add patch sources +sim.add_patch_sources( + location=select_random, + location_params=dict(n=3), + waveform=narrowband_oscillation, + waveform_params=dict(fmin=8, fmax=12), + extents=[10, 20, 50], + subject="sample", + subjects_dir=subjects_dir, +) + +# %% +# Now we simulate the configuration with an arbitrary level of global SNR: + +sc = sim.simulate( + sfreq, + duration, + fwd=fwd_eeg, + snr_global=3, + snr_params=dict(fmin=8, fmax=12), + random_state=1234, +) + +# %% +# We can now plot the source configuration using the dedicated method +# :py:meth:`~meegsim.configuration.SourceConfiguration.plot()` of the +# :py:class:`~meegsim.configuration.SourceConfiguration` class. The +# method returns a :py:class:`~mne.viz.Brain` object, which can be +# used to plot additional information, e.g., parcellation of interest. + +brain = sc.plot( + scale_factors=dict(point=1.25), + subject="sample", + subjects_dir=subjects_dir, + size=(1000, 800), + background="black", + hemi="split", + views=["lat", "med"], +) +brain.add_annotation("aparc") + +# %% +# .. image:: ../../_static/example_stubs/images/sphx_glr_02_plot_brain_configuration_001.png +# :alt: brain plot configuration +# :width: 600 diff --git a/examples/basics/GALLERY_HEADER.rst b/examples/basics/GALLERY_HEADER.rst new file mode 100644 index 0000000..53a6f53 --- /dev/null +++ b/examples/basics/GALLERY_HEADER.rst @@ -0,0 +1,4 @@ +Basics +====== + +These examples illustrate some basic aspects of using MEEGsim for simulations. diff --git a/examples/building_blocks/00_plot_kappa_von_mises.py b/examples/building_blocks/00_plot_kappa_von_mises.py new file mode 100644 index 0000000..a3611f2 --- /dev/null +++ b/examples/building_blocks/00_plot_kappa_von_mises.py @@ -0,0 +1,77 @@ +""" +Phase-phase coupling based on von Mises distribution +==================================================== + +.. currentmodule:: meegsim.coupling + +This example illustrates how the choice of kappa affects the obtained phase-phase +coupling if the :func:`ppc_von_mises` method based on the von Mises distribution is used. +""" + +import numpy as np +from harmoni.extratools import compute_plv +from matplotlib import pyplot as plt + +from meegsim.coupling import ppc_von_mises +from meegsim.waveform import narrowband_oscillation + +# %% +# To illustrate the effect, we consider a range of values from 0.01 to 10 for `kappa`, +# performing several simulations for each value. In addition, we vary the length of +# the simulated data: + +lengths = [10, 30, 60] +n_kappas = 11 +n_runs = 25 +kappas = np.logspace(-2, 1, n_kappas) +sfreq = 250 +fmin = 8 +fmax = 12 + +# %% +# For each considered setting (kappa, data length), we randomly generate a narrowband +# alpha (8-12 Hz) oscillation and a coupled oscillation with the mean phase lag +# of :math:`{\pi}/2`. We then estimate the phase-locking value (PLV) of the simulated +# oscillations and plot the mean PLV (black line) as well as its standard deviation +# (shaded region: +-1.96SD): + +fig, axes = plt.subplots(ncols=len(lengths), figsize=(9, 4), layout="constrained") +for ax, length in zip(axes, lengths): + times = np.arange(0, length, 1 / sfreq) + waveform = narrowband_oscillation(1, times, fmin=fmin, fmax=fmax) + + phase_lag = np.pi / 4 + plv = np.zeros((n_runs, n_kappas)) + for i_run in range(n_runs): + for i_kappa, kappa in enumerate(kappas): + result = ppc_von_mises( + waveform, sfreq, phase_lag, kappa=kappa, fmin=fmin, fmax=fmax + ) + cplv = compute_plv(waveform, result, m=1, n=1, plv_type="complex") + plv[i_run, i_kappa] = np.abs(cplv)[0][0] + + plv_mean = np.mean(plv, axis=0) + plv_sd = np.std(plv, axis=0) + ax.plot(kappas, plv_mean, c="k", linewidth=1.5) + ax.fill_between( + kappas, plv_mean - 1.96 * plv_sd, plv_mean + 1.96 * plv_sd, color="grey" + ) + ax.set_xscale("log") + ax.set_xticks([0.01, 0.1, 1, 10]) + ax.set_xticklabels([0.01, 0.1, 1, 10]) + ax.set_ylim([0, 1]) + ax.set_xlabel("Kappa") + ax.set_ylabel("PLV") + ax.set_title(f"duration = {length} s") + +# %% +# As shown in the plots above, kappa is monotoneously related to the resulting PLV +# between coupled time series. This allows for flexible control of coupling, and +# these plots can be used as reference when picking a suitable value of kappa. +# For lower values of kappa and shorter recordings, the estimated connectivity +# has a non-zero noise floor (mean does not reach 0) and becomes less stable +# (the SD becomes larger). +# +# If the connectivity measure depends on both amplitude and phase (e.g., coherence), +# its estimated value will also depend on the correspondence between amplitude +# envelopes of the time series. diff --git a/examples/plot_coupling.py b/examples/building_blocks/01_plot_coupling.py similarity index 97% rename from examples/plot_coupling.py rename to examples/building_blocks/01_plot_coupling.py index 9b4d22c..3d18ec3 100644 --- a/examples/plot_coupling.py +++ b/examples/building_blocks/01_plot_coupling.py @@ -1,4 +1,5 @@ """ +================================================== Phase-phase coupling using shifted copy with noise ================================================== @@ -61,6 +62,7 @@ coh=target_coh, fmin=fmin, fmax=fmax, + band_limited=False, random_state=seed, ) # x is leading y @@ -108,4 +110,3 @@ ax_lag.set_xlabel("Target coherence") ax_lag.set_ylabel("Obtained phase lag (degrees)") fig.tight_layout() -plt.show() diff --git a/examples/building_blocks/02_plot_global_snr.py b/examples/building_blocks/02_plot_global_snr.py new file mode 100644 index 0000000..80a8563 --- /dev/null +++ b/examples/building_blocks/02_plot_global_snr.py @@ -0,0 +1,82 @@ +""" +Adjustment of global SNR +======================== + +This example shows how the global SNR can be adjusted. +""" + +import mne +import matplotlib.pyplot as plt + +from mne.datasets import sample + +from meegsim.location import select_random +from meegsim.simulate import SourceSimulator +from meegsim.waveform import narrowband_oscillation + +# %% +# First, we load the head model and associated source space: + +# Paths +data_path = sample.data_path() / "MEG" / "sample" +fwd_path = data_path / "sample_audvis-meg-eeg-oct-6-fwd.fif" +raw_path = data_path / "sample_audvis_raw.fif" + +# Load the prerequisites: fwd, src, and info +fwd = mne.read_forward_solution(fwd_path) +fwd = mne.convert_forward_solution(fwd, force_fixed=True) +raw = mne.io.read_raw(raw_path) +src = fwd["src"] +info = raw.info + +# Pick EEG channels only +eeg_idx = mne.pick_types(info, eeg=True) +info_eeg = mne.pick_info(info, eeg_idx) +fwd_eeg = fwd.pick_channels(info_eeg.ch_names) + +# %% +# We simulate the same configuration (100 noise sources and 3 point sources) +# several times with different levels of SNR. As shown in the picture below, +# the average alpha power increases relative to the 1/f level with higher SNR: + +# Simulation parameters +sfreq = 250 +duration = 60 +seed = 123 + +fig, axes = plt.subplots(ncols=3, figsize=(8, 3)) +snr_values = [1, 5, 10] + +for i_snr, target_snr in enumerate(snr_values): + sim = SourceSimulator(src) + + # Select some vertices randomly + sim.add_point_sources( + location=select_random, + waveform=narrowband_oscillation, + location_params=dict(n=3), + waveform_params=dict(fmin=8, fmax=12), + names=["s1", "s2", "s3"], + ) + + sim.add_noise_sources(location=select_random, location_params=dict(n=100)) + + sc = sim.simulate( + sfreq, + duration, + fwd=fwd, + snr_global=target_snr, + snr_params=dict(fmin=8, fmax=12), + random_state=seed, + ) + raw = sc.to_raw(fwd, info) + + spec = raw.compute_psd(fmax=40, n_fft=sfreq, n_overlap=sfreq // 2, n_per_seg=sfreq) + spec.plot(average=True, dB=False, axes=axes[i_snr], amplitude=False) + + axes[i_snr].set_title(f"SNR={target_snr}") + axes[i_snr].set_xlabel("Frequency (Hz)") + axes[i_snr].set_ylabel("PSD (uV^2/Hz)") + axes[i_snr].set_ylim([0, 0.125]) + +fig.tight_layout() diff --git a/examples/snr_adjustment.py b/examples/building_blocks/03_plot_local_snr.py similarity index 56% rename from examples/snr_adjustment.py rename to examples/building_blocks/03_plot_local_snr.py index 399f4dd..bedced7 100644 --- a/examples/snr_adjustment.py +++ b/examples/building_blocks/03_plot_local_snr.py @@ -1,46 +1,54 @@ """ -Testing the adjustment of SNR -============================= +Adjustment of local SNR +======================== +This example shows how the local SNR can be adjusted. """ import mne import matplotlib.pyplot as plt -from pathlib import Path +from mne.datasets import sample from meegsim.location import select_random from meegsim.simulate import SourceSimulator from meegsim.waveform import narrowband_oscillation +# %% +# First, we load the head model and associated source space: -# Load the head model -fs_dir = Path("~/mne_data/MNE-fsaverage-data/fsaverage/") -fwd_path = fs_dir / "bem_copy" / "fsaverage-oct6-fwd.fif" -src_path = fs_dir / "bem_copy" / "fsaverage-oct6-src.fif" -src = mne.read_source_spaces(src_path) +# Paths +data_path = sample.data_path() / "MEG" / "sample" +fwd_path = data_path / "sample_audvis-meg-eeg-oct-6-fwd.fif" +raw_path = data_path / "sample_audvis_raw.fif" + +# Load the prerequisites: fwd, src, and info fwd = mne.read_forward_solution(fwd_path) +fwd = mne.convert_forward_solution(fwd, force_fixed=True) +raw = mne.io.read_raw(raw_path) +src = fwd["src"] +info = raw.info + +# Pick EEG channels only +eeg_idx = mne.pick_types(info, eeg=True) +info_eeg = mne.pick_info(info, eeg_idx) +fwd_eeg = fwd.pick_channels(info_eeg.ch_names) + +# %% +# We simulate the same configuration (100 noise sources and 3 point sources) +# several times with different levels of SNR. As shown in the picture below, +# the average alpha power increases relative to the 1/f level with higher SNR: # Simulation parameters sfreq = 250 duration = 60 seed = 123 -# Channel info -montage = mne.channels.make_standard_montage("standard_1020") -ch_names = [ch for ch in montage.ch_names if ch not in ["O9", "O10"]] -info = mne.create_info(ch_names, sfreq, ch_types="eeg") -info.set_montage("standard_1020") - -# Adapt fwd to the info (could be done by our structure in principle) -fwd = mne.convert_forward_solution(fwd, force_fixed=True) -fwd = mne.pick_channels_forward(fwd, info.ch_names, ordered=True) - fig, axes = plt.subplots(ncols=3, figsize=(8, 3)) -snr_values = [1, 5, 10] +snr_values = [1, 2.5, 5] for i_snr, target_snr in enumerate(snr_values): - sim = SourceSimulator(src) + sim = SourceSimulator(src, snr_mode="local") # Select some vertices randomly sim.add_point_sources( @@ -53,7 +61,7 @@ names=["s1", "s2", "s3"], ) - sim.add_noise_sources(location=select_random, location_params=dict(n=10)) + sim.add_noise_sources(location=select_random, location_params=dict(n=100)) sc = sim.simulate(sfreq, duration, fwd=fwd, random_state=seed) raw = sc.to_raw(fwd, info) @@ -64,7 +72,6 @@ axes[i_snr].set_title(f"SNR={target_snr}") axes[i_snr].set_xlabel("Frequency (Hz)") axes[i_snr].set_ylabel("PSD (uV^2/Hz)") - axes[i_snr].set_ylim([0, 1.25]) + axes[i_snr].set_ylim([0, 0.125]) fig.tight_layout() -plt.show() diff --git a/examples/building_blocks/04_plot_brain_std.py b/examples/building_blocks/04_plot_brain_std.py new file mode 100644 index 0000000..2b18bf0 --- /dev/null +++ b/examples/building_blocks/04_plot_brain_std.py @@ -0,0 +1,136 @@ +""" +Controlling the standard deviation of activity +============================================== + +This example illustrates how the standard deviation (SD) of source activity can be +manipulated. +""" + +# sphinx_gallery_thumbnail_path = '_static/example_stubs/thumb/sphx_glr_04_plot_brain_std_thumb.png' + +import matplotlib.pyplot as plt +import mne +import numpy as np + +from mne.datasets import sample + +from meegsim.location import select_random +from meegsim.simulate import SourceSimulator +from meegsim.waveform import narrowband_oscillation + +# %% +# First, we load the head model and associated source space: + +# Paths +subjects_dir = sample.data_path() / "subjects" +data_path = sample.data_path() / "MEG" / "sample" +fwd_path = data_path / "sample_audvis-meg-eeg-oct-6-fwd.fif" +raw_path = data_path / "sample_audvis_raw.fif" + +# Load the prerequisites: fwd, src, and info +fwd = mne.read_forward_solution(fwd_path) +fwd = mne.convert_forward_solution(fwd, force_fixed=True) +raw = mne.io.read_raw(raw_path) +src = fwd["src"] +info = raw.info + +# Pick EEG channels only +eeg_idx = mne.pick_types(info, eeg=True) +info_eeg = mne.pick_info(info, eeg_idx) +fwd_eeg = fwd.pick_channels(info_eeg.ch_names) + +# %% +# Simulation parameters are listed below: + +sfreq = 250 +duration = 60 +seed = 1234 +target_snr = 4 +fmin = 8 +fmax = 12 + +# %% +# To illustrate the effect, we set the SD based on the y-position of the sources, +# with higher SDs for parieto-occipital areas. By wrapping the SD values in a +# :class:`mne.SourceEstimate` object, we can set the SD for whole sources at once +# even if there actual positions are generated randomly. In this case, however, we +# pick the source locations manually to show the effect better with one frontal +# and one occipital source: + +ypos = np.hstack([s["rr"][s["inuse"] > 0, 1] for s in src]) +std = 1 - 8 * ypos +vertno = [s["vertno"] for s in src] +std_stc = mne.SourceEstimate( + data=std, vertices=vertno, tmin=0, tstep=0.01, subject="sample" +) +source_vertno = [126371, 10957] # frontal & occipital + +# The resulting spatial distribution of SD along with chosen locations for patch +# sources (white borders) are shown below: + +patches = mne.grow_labels( + subject="sample", + seeds=source_vertno, + extents=10, + hemis=[0, 0], + subjects_dir=subjects_dir, +) +brain = std_stc.plot( + subject="sample", + subjects_dir=subjects_dir, + hemi="lh", + views="lat", + clim=dict(kind="value", lims=[0, 1, 2]), + transparent=False, + background="white", +) +for patch in patches: + brain.add_label(patch, color="white", borders=True) + +# %% +# .. image:: ../../_static/example_stubs/images/sphx_glr_04_plot_brain_std_001.png +# :alt: custom SD with source locations +# :width: 600 + +# %% +# Below, we create two identical simulations except for the standard deviation +# (``std`` argument of :meth:`~meegsim.simulate.SourceSimulator.add_patch_sources` call). +# We then illustrate the difference in the topomap of alpha power betweeen two cases: +# equal SD for all sources (``std=1.0``) vs. custom variance based on ``std_stc``. +# As expected, occipital source dominates the topomap when the custom SD is used: + +fig, axes = plt.subplots(ncols=2, figsize=(8, 4)) +for ax, std, case in zip(axes, [1.0, std_stc], ["equal", "unequal"]): + sim = SourceSimulator(src) + + sim.add_noise_sources(location=select_random, location_params=dict(n=10)) + + # Use manually selected vertices, put all sources in the left hemisphere + sim.add_patch_sources( + location=[(0, v) for v in source_vertno], + waveform=narrowband_oscillation, + location_params=dict(n=3), + waveform_params=dict(fmin=fmin, fmax=fmax), + std=std, + extents=10, + subject="sample", + subjects_dir=subjects_dir, + ) + + sc = sim.simulate( + sfreq, + duration, + fwd=fwd, + snr_global=target_snr, + snr_params=dict(fmin=fmin, fmax=fmax), + random_state=seed, + ) + raw = sc.to_raw(fwd, info, sensor_noise_level=0.05) + + spec = raw.compute_psd(n_fft=sfreq, n_overlap=sfreq // 2, n_per_seg=sfreq) + spec.plot_topomap(bands={case: (8, 12)}, axes=ax) + +# %% +# .. image:: ../../_static/example_stubs/images/sphx_glr_04_plot_brain_std_002.png +# :alt: topomaps of alpha power in both cases +# :width: 600 diff --git a/examples/sensor_space_noise.py b/examples/building_blocks/05_plot_sensor_space_noise.py similarity index 60% rename from examples/sensor_space_noise.py rename to examples/building_blocks/05_plot_sensor_space_noise.py index 6770469..77e37f9 100644 --- a/examples/sensor_space_noise.py +++ b/examples/building_blocks/05_plot_sensor_space_noise.py @@ -1,41 +1,43 @@ """ -Testing the sensor space noise -============================== +Adding sensor space noise +========================= """ import mne import matplotlib.pyplot as plt -from pathlib import Path +from mne.datasets import sample from meegsim.location import select_random from meegsim.simulate import SourceSimulator from meegsim.waveform import narrowband_oscillation +# %% +# First, we load the head model and associated source space: -# Load the head model -fs_dir = Path("~/mne_data/MNE-fsaverage-data/fsaverage/") -fwd_path = fs_dir / "bem_copy" / "fsaverage-oct6-fwd.fif" -src_path = fs_dir / "bem_copy" / "fsaverage-oct6-src.fif" -src = mne.read_source_spaces(src_path) +# Paths +data_path = sample.data_path() / "MEG" / "sample" +fwd_path = data_path / "sample_audvis-meg-eeg-oct-6-fwd.fif" +raw_path = data_path / "sample_audvis_raw.fif" + +# Load the prerequisites: fwd, src, and info fwd = mne.read_forward_solution(fwd_path) +fwd = mne.convert_forward_solution(fwd, force_fixed=True) +raw = mne.io.read_raw(raw_path) +src = fwd["src"] +info = raw.info + +# Pick EEG channels only +eeg_idx = mne.pick_types(info, eeg=True) +info_eeg = mne.pick_info(info, eeg_idx) +fwd_eeg = fwd.pick_channels(info_eeg.ch_names) # Simulation parameters sfreq = 250 duration = 60 seed = 123 -# Channel info -montage = mne.channels.make_standard_montage("standard_1020") -ch_names = [ch for ch in montage.ch_names if ch not in ["O9", "O10"]] -info = mne.create_info(ch_names, sfreq, ch_types="eeg") -info.set_montage("standard_1020") - -# Adapt fwd to the info (could be done by our structure in principle) -fwd = mne.convert_forward_solution(fwd, force_fixed=True) -fwd = mne.pick_channels_forward(fwd, info.ch_names, ordered=True) - sim = SourceSimulator(src) # Select some vertices randomly @@ -54,7 +56,7 @@ sc = sim.simulate(sfreq, duration, fwd=fwd, random_state=seed) raw = sc.to_raw(fwd, info) -noise_levels = [0, 0.05, 0.1, 0.25, 0.5, 1] +noise_levels = [0.05, 0.25, 0.5] n_levels = len(noise_levels) fig, axes = plt.subplots(ncols=n_levels, figsize=(3 * n_levels, 3)) @@ -62,10 +64,9 @@ raw = sc.to_raw(fwd, info, sensor_noise_level=noise_level) spec = raw.compute_psd(fmax=60, n_fft=sfreq, n_overlap=sfreq // 2, n_per_seg=sfreq) - spec.plot(axes=axes[i_level], amplitude=False, sphere="eeglab") + spec.plot(axes=axes[i_level], amplitude=False) axes[i_level].set_title(f"{noise_level=}") axes[i_level].set_xlabel("Frequency (Hz)") fig.tight_layout() -plt.show() diff --git a/examples/building_blocks/GALLERY_HEADER.rst b/examples/building_blocks/GALLERY_HEADER.rst new file mode 100644 index 0000000..cd0a841 --- /dev/null +++ b/examples/building_blocks/GALLERY_HEADER.rst @@ -0,0 +1,5 @@ +Building blocks +=============== + +These examples illustrate the building blocks provided by MEEGsim and highlights +the practical nuances related to their usage. diff --git a/examples/dummy.py b/examples/dummy.py deleted file mode 100644 index 015dff5..0000000 --- a/examples/dummy.py +++ /dev/null @@ -1,115 +0,0 @@ -""" -Testing the configuration structure -=================================== - -""" - -import json -import matplotlib.pyplot as plt -import mne -import numpy as np - -from pathlib import Path - -from meegsim.location import select_random -from meegsim.simulate import SourceSimulator -from meegsim.waveform import narrowband_oscillation - - -def to_json(sources): - return json.dumps({k: str(s) for k, s in sources.items()}, indent=4) - - -def data2stc(data, src): - vertno = [s["vertno"] for s in src] - return mne.SourceEstimate( - data=data, vertices=vertno, tmin=0, tstep=0.01, subject="fsaverage" - ) - - -def extents_from_areas_cm2(areas_cm2): - return list(np.sqrt(np.array(areas_cm2) * 100 / np.pi)) - - -# Load the head model -fs_dir = Path("~/mne_data/MNE-fsaverage-data/fsaverage/") -fwd_path = fs_dir / "bem_copy" / "fsaverage-oct6-fwd.fif" -src_path = fs_dir / "bem_copy" / "fsaverage-oct6-src.fif" -src = mne.read_source_spaces(src_path) -fwd = mne.read_forward_solution(fwd_path) - -# Simulation parameters -sfreq = 250 -duration = 60 -seed = 123 -target_snr = 20 - -# Channel info -montage = mne.channels.make_standard_montage("standard_1020") -ch_names = [ - ch for ch in montage.ch_names if ch not in ["O9", "O10", "T3", "T4", "T5", "T6"] -] -info = mne.create_info(ch_names, sfreq, ch_types="eeg") -info.set_montage("standard_1020") - -# Adapt fwd to the info (could be done by our structure in principle) -fwd = mne.convert_forward_solution(fwd, force_fixed=True) -fwd = mne.pick_channels_forward(fwd, info.ch_names, ordered=True) - -# Create a dummy stc for std based on the y-position of the sources -ypos = np.hstack([1 - 8 * np.abs(s["rr"][s["inuse"] > 0, 1]) for s in src]) -std_stc = data2stc(ypos, src) -std_stc.plot( - subject="fsaverage", - hemi="split", - views=["lat", "med"], - clim=dict(kind="value", lims=[0, 1, 2]), - transparent=False, - background="white", -) - -sim = SourceSimulator(src) - -sim.add_noise_sources(location=select_random, location_params=dict(n=10)) - -# Select some vertices randomly -sim.add_patch_sources( - location=select_random, - waveform=narrowband_oscillation, - location_params=dict(n=3), - waveform_params=dict(fmin=8, fmax=12), - std=std_stc, - extents=extents_from_areas_cm2([2, 4, 8]), -) - -sc = sim.simulate( - sfreq, - duration, - fwd=fwd, - snr_global=4, - snr_params=dict(fmin=8, fmax=12), - random_state=seed, -) -stc = sc.to_stc() -raw = sc.to_raw(fwd, info, sensor_noise_level=0.05) - -source_std = np.std(stc.data, axis=1) -lim = np.max(source_std) -std_stc_est = mne.SourceEstimate(source_std, stc.vertices, tmin=0, tstep=0.01) -std_stc_est.plot( - subject="fsaverage", - hemi="split", - views=["lat", "med"], - clim=dict(kind="value", lims=[0, lim / 2, lim]), - colormap="Reds", - time_viewer=False, - transparent=False, - background="white", -) - -sc.plot(subject="fsaverage", hemi="split", views=["lat", "med"]) - -spec = raw.compute_psd(n_fft=sfreq, n_overlap=sfreq // 2, n_per_seg=sfreq) -spec.plot_topomap(bands={"alpha": (8, 12)}, sphere="eeglab") -spec.plot(sphere="eeglab") -plt.show(block=True) diff --git a/examples/random_state.py b/examples/random_state.py deleted file mode 100644 index f186ed6..0000000 --- a/examples/random_state.py +++ /dev/null @@ -1,86 +0,0 @@ -""" -Random state and reproducibility --------------------------------- - -This example will showcase how reproducible source configurations -can be achieved by fixing the random state when simulating. -""" - -import matplotlib.pyplot as plt -import mne -import numpy as np - -from matplotlib.gridspec import GridSpec -from pathlib import Path - -from meegsim.coupling import ppc_shifted_copy_with_noise -from meegsim.location import select_random -from meegsim.waveform import narrowband_oscillation -from meegsim.simulate import SourceSimulator - - -# Load the head model -fs_dir = Path("~/mne_data/MNE-fsaverage-data/fsaverage/") -fwd_path = fs_dir / "bem_copy" / "fsaverage-oct6-fwd.fif" -src_path = fs_dir / "bem_copy" / "fsaverage-oct6-src.fif" -src = mne.read_source_spaces(src_path) -fwd = mne.read_forward_solution(fwd_path) - - -# Simulation parameters -sfreq = 250 -duration = 60 -seed = 123 -n_sources = 3 - -# Channel info -montage = mne.channels.make_standard_montage("standard_1020") -ch_names = [ - ch for ch in montage.ch_names if ch not in ["O9", "O10", "T3", "T4", "T5", "T6"] -] -info = mne.create_info(ch_names, sfreq, ch_types="eeg") -info.set_montage("standard_1020") - -# Adapt fwd to the info (could be done by our structure in principle) -fwd = mne.convert_forward_solution(fwd, force_fixed=True) -fwd = mne.pick_channels_forward(fwd, info.ch_names, ordered=True) - - -sim = SourceSimulator(src) - -# Select some vertices randomly -for i in range(3): - sim.add_point_sources( - location=select_random, - waveform=narrowband_oscillation, - location_params=dict(n=1), - waveform_params=dict(fmin=8, fmax=12), - names=[str(i)], - ) -sim.set_coupling( - ("0", "1"), - method=ppc_shifted_copy_with_noise, - phase_lag=np.pi / 2, - coh=0.5, - fmin=8, - fmax=12, -) - - -sc = sim.simulate(sfreq=sfreq, duration=duration, random_state=1234) -brain = sc.plot(subject="fsaverage", hemi="split", views=["lat", "med"]) -screenshot = brain.screenshot() -brain.close() - -print([s.vertno for s in sc._sources.values()]) - -fig = plt.figure(figsize=(12, 6)) -gs = GridSpec(3, 2, figure=fig) -ax_brain = fig.add_subplot(gs[:, 0]) -ax_brain.imshow(screenshot) -ax_brain.axis("off") -for i in range(3): - s = sc._sources[str(i)] - ax = fig.add_subplot(gs[i, 1]) - ax.plot(sc.times[:1000], np.squeeze(s.waveform)[:1000]) -fig.tight_layout() diff --git a/examples/select_kappa.py b/examples/select_kappa.py deleted file mode 100644 index 4c0d0bb..0000000 --- a/examples/select_kappa.py +++ /dev/null @@ -1,38 +0,0 @@ -""" -Phase-phase coupling based on von Mises distribution -==================================================== - -This example illustrates how the choice of kappa affects the obtained phase-phase -coupling if the method based on von Mises distribution is used. -""" - -import numpy as np -from harmoni.extratools import compute_plv -from matplotlib import pyplot as plt -from meegsim.coupling import ppc_von_mises -from meegsim.utils import get_sfreq - - -lenghts = [0.5, 1, 10, 60] -n_kappas = 100 -n_runs = 1000 -kappas = np.logspace(-5, 3, n_kappas) -fs = 1000 -for length in lenghts: - print(length) - times = np.arange(0, length, 1 / fs) - waveform = np.sin(2 * np.pi * 10 * times) - phase_lag = 0 - plv = np.zeros((n_runs, n_kappas)) - for run in range(n_runs): - for ikappa, kappa in enumerate(kappas): - result = ppc_von_mises(waveform, get_sfreq(times), phase_lag, kappa=kappa) - cplv = compute_plv(waveform, result, m=1, n=1, plv_type="complex") - plv[run, ikappa] = np.abs(cplv)[0][0] - - plt.plot(np.log10(kappas), plv.T, c="grey") - plt.plot(np.log10(kappas), np.mean(plv, axis=0), c="k", linewidth=3) - plt.xlabel("Kappa, logscale") - plt.ylabel("plv") - plt.savefig("kappa-to-plv_length" + str(length) + "sec.png") - plt.close() diff --git a/pyproject.toml b/pyproject.toml index 98032ab..01d72e2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -48,6 +48,7 @@ docs = [ "intersphinx-registry", "harmoni", "myst-parser", + "neurodsp", "numpydoc", "pydata-sphinx-theme", "sphinx", diff --git a/src/meegsim/coupling.py b/src/meegsim/coupling.py index 003d1c9..ae7f4e6 100644 --- a/src/meegsim/coupling.py +++ b/src/meegsim/coupling.py @@ -321,7 +321,7 @@ def ppc_shifted_copy_with_noise( multiple simulations. For every individual output, coherence and phase lag might deviate from the desired values: The lower the requested coherence is, the higher is the variance of the output. For more information, see :doc:`this example - `. + `. """ check_numeric("coherence", coh, bounds=(0, 1), allow_none=False) snr = _get_required_snr(coh, band_limited) diff --git a/src/meegsim/simulate.py b/src/meegsim/simulate.py index 135baef..7cd55c9 100644 --- a/src/meegsim/simulate.py +++ b/src/meegsim/simulate.py @@ -1,5 +1,6 @@ import networkx as nx import numpy as np +import warnings from ._check import check_coupling, check_option, check_numeric_array, check_snr_params from .configuration import SourceConfiguration @@ -152,8 +153,15 @@ def add_point_sources( self._sources.extend(point_sg.names) # Check if SNR should be adjusted - if point_sg.snr is not None and self.snr_mode == "local": - self.is_local_snr_adjusted = True + if point_sg.snr is not None: + if self.snr_mode == "local": + self.is_local_snr_adjusted = True + else: + warnings.warn( + "Ignoring the provided value of local SNR since global adjustment " + "is enabled. To enable the local adjustment, set snr_mode to 'local' " + "when initializing the SourceSimulator." + ) # Return the names of newly added sources return point_sg.names @@ -236,7 +244,7 @@ def add_patch_sources( subject_dir : str, optional Path to the directory with FreeSurfer output, only used when growing patch sources from the central vertex. If None (default), it is resolved automatically - by MNE-Python. + by MNE-Python. Provide the path explicitly if errors arise. names : list, optional A list of names for each source. If not specified, the names will be autogenerated using the format 'auto-sgN-sM', where N is the index @@ -271,8 +279,15 @@ def add_patch_sources( self._sources.extend(patch_sg.names) # Check if SNR should be adjusted - if patch_sg.snr is not None and self.snr_mode == "local": - self.is_local_snr_adjusted = True + if patch_sg.snr is not None: + if self.snr_mode == "local": + self.is_local_snr_adjusted = True + else: + warnings.warn( + "Ignoring the provided value of local SNR since global adjustment " + "is enabled. To enable the local adjustment, set snr_mode to 'local' " + "when initializing the SourceSimulator." + ) # Return the names of newly added sources return patch_sg.names @@ -472,6 +487,12 @@ def simulate( # Check the forward model and auto-fill info if needed is_global_snr_adjusted = self.snr_mode == "global" and snr_global is not None is_local_snr_adjusted = self.snr_mode == "local" and self.is_local_snr_adjusted + if snr_global is not None and self.snr_mode == "local": + warnings.warn( + "Ignoring the provided value of global SNR since local adjustment " + "is enabled. To enable the global adjustment, set snr_mode to 'global' " + "when initializing the SourceSimulator." + ) if (is_global_snr_adjusted or is_local_snr_adjusted) and fwd is None: raise ValueError("A forward model is required for the adjustment of SNR.")