ENH: Add support for Info-as-Raw (#5869)

Author: larsoner

doc/whats_new.rst

@@ -39,6 +39,8 @@ Changelog
 
 - Add :func:`mne.labels_to_stc` to facilitate working with label data, by `Eric Larson`_
 
+- Add support for using :class:`mne.Info` and ``duration`` in :func:`mne.simulation.simulate_raw` instead of :class:`mne.io.Raw` by `Eric Larson`_
+
 - Add ``overlap`` argument to :func:`mne.make_fixed_length_events` by `Eric Larson`_
 
 - Add 448-labels subdivided aparc cortical parcellation by `Denis Engemann`_ and `Sheraz Khan`_

mne/simulation/raw.py

@@ -17,7 +17,7 @@
 from ..source_estimate import VolSourceEstimate
 from ..cov import make_ad_hoc_cov, read_cov, Covariance
 from ..bem import fit_sphere_to_headshape, make_sphere_model, read_bem_solution
-from ..io import RawArray, BaseRaw
+from ..io import RawArray, BaseRaw, Info
 from ..chpi import (read_head_pos, head_pos_to_trans_rot_t, _get_hpi_info,
                     _get_hpi_initial_fit)
 from ..io.constants import FIFF
@@ -30,7 +30,8 @@
 from ..source_space import (_ensure_src, _points_outside_surface,
                             _adjust_patch_info)
 from ..source_estimate import _BaseSourceEstimate
-from ..utils import logger, verbose, check_random_state, warn, _pl
+from ..utils import (logger, verbose, check_random_state, warn, _pl,
+                     _validate_type)
 from ..parallel import check_n_jobs
 
 
@@ -65,17 +66,20 @@ def _log_ch(start, info, ch):
 def simulate_raw(raw, stc, trans, src, bem, cov='simple',
                  blink=False, ecg=False, chpi=False, head_pos=None,
                  mindist=1.0, interp='cos2', iir_filter=None, n_jobs=1,
-                 random_state=None, use_cps=True, forward=None, verbose=None):
+                 random_state=None, use_cps=True, forward=None,
+                 duration=None, verbose=None):
     u"""Simulate raw data.
 
     Head movements can optionally be simulated using the ``head_pos``
     parameter.
 
     Parameters
     ----------
-    raw : instance of Raw
+    raw : instance of Raw | instance of Info
         The raw template to use for simulation. The ``info``, ``times``,
-        and potentially ``first_samp`` properties will be used.
+        and ``first_samp`` properties will be used.
+        Can be ``info`` if ``duration`` is also supplied, in which case
+        ``first_samp`` will be set to zero.
     stc : instance of SourceEstimate
         The source estimate to use to simulate data. Must have the same
         sample rate as the raw data.
@@ -136,6 +140,13 @@ def simulate_raw(raw, stc, trans, src, bem, cov='simple',
     forward : instance of Forward | None
         The forward operator to use. If None (default) it will be computed
         using `bem`, `trans`, and `src`.
+
+        .. versionadded:: 0.17
+    duration : float | None
+        The duration to simulate. Can be None to use the duration of ``raw``.
+        Must be supplied if ``raw`` is an instance of :class:`mne.Info`.
+
+        .. versionadded:: 0.18
     verbose : bool, str, int, or None
         If not None, override default verbose level (see :func:`mne.verbose`
         and :ref:`Logging documentation <tut_logging>` for more).
@@ -203,10 +214,22 @@ def simulate_raw(raw, stc, trans, src, bem, cov='simple',
     .. [1] Bentivoglio et al. "Analysis of blink rate patterns in normal
            subjects" Movement Disorders, 1997 Nov;12(6):1028-34.
     """
-    if not isinstance(raw, BaseRaw):
-        raise TypeError('raw should be an instance of Raw')
-    times, info, first_samp = raw.times, raw.info, raw.first_samp
-    raw_verbose = raw.verbose
+    _validate_type(raw, (BaseRaw, Info), 'raw', 'Raw or Info')
+    if duration is not None:
+        duration = float(duration)
+    if isinstance(raw, Info):
+        if duration is None:
+            raise ValueError('duration cannot be None if raw is an instance '
+                             'of Info')
+        info, first_samp = raw, 0
+        raw_verbose = verbose
+    else:
+        info, first_samp = raw.info, raw.first_samp
+        raw_verbose = raw.verbose
+        if duration is None:
+            duration = raw.times[-1] + 1. / info['sfreq']
+    times = np.arange(int(round(info['sfreq'] * duration))) / info['sfreq']
+    del raw
 
     # Check for common flag errors and try to override
     if not isinstance(stc, _BaseSourceEstimate):
@@ -227,7 +250,7 @@ def simulate_raw(raw, stc, trans, src, bem, cov='simple',
             raise ValueError('If forward is not None then trans, src, bem, '
                              'and head_pos must all be None')
         if not np.allclose(forward['info']['dev_head_t']['trans'],
-                           raw.info['dev_head_t']['trans'], atol=1e-6):
+                           info['dev_head_t']['trans'], atol=1e-6):
             raise ValueError('The forward meg<->head transform '
                              'forward["info"]["dev_head_t"] does not match '
                              'the one in raw.info["dev_head_t"]')
@@ -241,7 +264,8 @@ def simulate_raw(raw, stc, trans, src, bem, cov='simple',
         head_pos = head_pos_to_trans_rot_t(head_pos)
     if isinstance(head_pos, tuple):  # can be quats converted to trans, rot, t
         transs, rots, ts = head_pos
-        ts -= raw._first_time  # MF files need reref
+        first_time = first_samp / info['sfreq']
+        ts -= first_time  # MF files need reref
         dev_head_ts = [np.r_[np.c_[r, t[:, np.newaxis]], [[0, 0, 0, 1]]]
                        for r, t in zip(rots, transs)]
         del transs, rots
@@ -267,7 +291,7 @@ def simulate_raw(raw, stc, trans, src, bem, cov='simple',
     dev_head_ts = [{'trans': d, 'to': info['dev_head_t']['to'],
                     'from': info['dev_head_t']['from']}
                    for d in dev_head_ts]
-    offsets = raw.time_as_index(ts)
+    offsets = np.round(ts * info['sfreq']).astype(int)
     offsets = np.concatenate([offsets, [len(times)]])
     assert offsets[-2] != offsets[-1]
     assert np.array_equal(offsets, np.unique(offsets))

mne/simulation/tests/test_raw.py

@@ -99,32 +99,30 @@ def test_simulate_raw_sphere():
     #
     raw_sim = simulate_raw(raw, stc, trans, src, sphere, read_cov(cov_fname),
                            head_pos=head_pos_sim,
-                           blink=True, ecg=True, random_state=seed,
-                           use_cps=True)
+                           blink=True, ecg=True, random_state=seed)
     raw_sim_2 = simulate_raw(raw, stc, trans_fname, src_fname, sphere,
                              cov_fname, head_pos=head_pos_sim,
-                             blink=True, ecg=True, random_state=seed,
-                             use_cps=True)
+                             blink=True, ecg=True, random_state=seed)
     assert_array_equal(raw_sim_2[:][0], raw_sim[:][0])
     std = dict(grad=2e-13, mag=10e-15, eeg=0.1e-6)
     raw_sim = simulate_raw(raw, stc, trans, src, sphere,
                            make_ad_hoc_cov(raw.info, std=std),
                            head_pos=head_pos_sim, blink=True, ecg=True,
-                           random_state=seed, use_cps=True)
+                           random_state=seed)
     raw_sim_2 = simulate_raw(raw, stc, trans_fname, src_fname, sphere,
                              cov=std, head_pos=head_pos_sim, blink=True,
-                             ecg=True, random_state=seed, use_cps=True)
+                             ecg=True, random_state=seed)
     assert_array_equal(raw_sim_2[:][0], raw_sim[:][0])
     sphere_norad = make_sphere_model('auto', None, raw.info)
     raw_meg = raw.copy().pick_types()
     raw_sim = simulate_raw(raw_meg, stc, trans, src, sphere_norad,
                            make_ad_hoc_cov(raw.info, std=None),
                            head_pos=head_pos_sim, blink=True, ecg=True,
-                           random_state=seed, use_cps=True)
+                           random_state=seed)
     raw_sim_2 = simulate_raw(raw_meg, stc, trans_fname, src_fname,
                              sphere_norad,
                              cov='simple', head_pos=head_pos_sim, blink=True,
-                             ecg=True, random_state=seed, use_cps=True)
+                             ecg=True, random_state=seed)
     assert_array_equal(raw_sim_2[:][0], raw_sim[:][0])
     # Test IO on processed data
     tempdir = _TempDir()
@@ -139,12 +137,10 @@ def test_simulate_raw_sphere():
     for ecg, eog in ((True, False), (False, True), (True, True)):
         raw_sim_3 = simulate_raw(raw, stc, trans, src, sphere,
                                  cov=None, head_pos=head_pos_sim,
-                                 blink=eog, ecg=ecg, random_state=seed,
-                                 use_cps=True)
+                                 blink=eog, ecg=ecg, random_state=seed)
         raw_sim_4 = simulate_raw(raw, stc, trans, src, sphere,
                                  cov=None, head_pos=head_pos_sim,
-                                 blink=False, ecg=False, random_state=seed,
-                                 use_cps=True)
+                                 blink=False, ecg=False, random_state=seed)
         picks = np.arange(len(raw.ch_names))
         diff_picks = pick_types(raw.info, meg=False, ecg=ecg, eog=eog)
         these_picks = np.setdiff1d(picks, diff_picks)
@@ -159,68 +155,74 @@ def test_simulate_raw_sphere():
     # make sure it works with EEG-only and MEG-only
     raw_sim_meg = simulate_raw(raw.copy().pick_types(meg=True, eeg=False),
                                stc, trans, src, sphere, cov=None,
-                               ecg=True, blink=True, random_state=seed,
-                               use_cps=True)
+                               ecg=True, blink=True, random_state=seed)
     raw_sim_eeg = simulate_raw(raw.copy().pick_types(meg=False, eeg=True),
                                stc, trans, src, sphere, cov=None,
-                               ecg=True, blink=True, random_state=seed,
-                               use_cps=True)
+                               ecg=True, blink=True, random_state=seed)
     raw_sim_meeg = simulate_raw(raw.copy().pick_types(meg=True, eeg=True),
                                 stc, trans, src, sphere, cov=None,
-                                ecg=True, blink=True, random_state=seed,
-                                use_cps=True)
+                                ecg=True, blink=True, random_state=seed)
     assert_allclose(np.concatenate((raw_sim_meg[:][0], raw_sim_eeg[:][0])),
                     raw_sim_meeg[:][0], rtol=1e-7, atol=1e-20)
     del raw_sim_meg, raw_sim_eeg, raw_sim_meeg
 
+    # check that raw-as-info is supported
+    raw_sim = simulate_raw(raw, stc, trans, src, sphere, cov=None)
+    n_samp = int(round(raw.info['sfreq']))
+    for use_raw in (raw, raw.info):
+        raw_sim_2 = simulate_raw(use_raw, stc, trans, src, sphere, cov=None,
+                                 duration=1.)
+        assert len(raw_sim_2.times) == n_samp
+        assert_allclose(raw_sim[:, :n_samp][0],
+                        raw_sim_2[:, :n_samp][0], rtol=1e-5, atol=1e-30)
+    del raw_sim, raw_sim_2
+
     # check that different interpolations are similar given small movements
     raw_sim = simulate_raw(raw, stc, trans, src, sphere, cov=None,
-                           head_pos=head_pos_sim, interp='linear',
-                           use_cps=True)
+                           head_pos=head_pos_sim, interp='linear')
     raw_sim_hann = simulate_raw(raw, stc, trans, src, sphere, cov=None,
-                                head_pos=head_pos_sim, interp='hann',
-                                use_cps=True)
+                                head_pos=head_pos_sim, interp='hann')
     assert_allclose(raw_sim[:][0], raw_sim_hann[:][0], rtol=1e-1, atol=1e-14)
     del raw_sim, raw_sim_hann
 
     # Make impossible transform (translate up into helmet) and ensure failure
     head_pos_sim_err = deepcopy(head_pos_sim)
     head_pos_sim_err[1.][2, 3] -= 0.1  # z trans upward 10cm
     pytest.raises(RuntimeError, simulate_raw, raw, stc, trans, src, sphere,
-                  ecg=False, blink=False, head_pos=head_pos_sim_err,
-                  use_cps=True)
+                  ecg=False, blink=False, head_pos=head_pos_sim_err)
     pytest.raises(RuntimeError, simulate_raw, raw, stc, trans, src,
                   bem_fname, ecg=False, blink=False,
-                  head_pos=head_pos_sim_err, use_cps=True)
+                  head_pos=head_pos_sim_err)
     # other degenerate conditions
-    pytest.raises(TypeError, simulate_raw, 'foo', stc, trans, src, sphere,
-                  use_cps=True)
-    pytest.raises(TypeError, simulate_raw, raw, 'foo', trans, src, sphere,
-                  use_cps=True)
+    pytest.raises(TypeError, simulate_raw, 'foo', stc, trans, src, sphere)
+    pytest.raises(TypeError, simulate_raw, raw, 'foo', trans, src, sphere)
     pytest.raises(ValueError, simulate_raw, raw, stc.copy().crop(0, 0),
-                  trans, src, sphere, use_cps=True)
+                  trans, src, sphere)
+    with pytest.raises(ValueError, match='duration cannot be None'):
+        simulate_raw(raw.info, stc, trans, src, sphere)
+    with pytest.raises(TypeError, match='must be an instance of Raw or Info'):
+        simulate_raw(0, stc, trans, src, sphere)
     stc_bad = stc.copy()
     stc_bad.tstep += 0.1
-    pytest.raises(ValueError, simulate_raw, raw, stc_bad, trans, src, sphere,
-                  use_cps=True)
+    pytest.raises(ValueError, simulate_raw, raw, stc_bad, trans, src, sphere)
     pytest.raises(TypeError, simulate_raw, raw, stc, trans, src, sphere,
-                  cov=0, use_cps=True)  # wrong covariance type
+                  cov=0)  # wrong covariance type
     pytest.raises(RuntimeError, simulate_raw, raw, stc, trans, src, sphere,
-                  chpi=True, use_cps=True)  # no cHPI info
+                  chpi=True)  # no cHPI info
     pytest.raises(ValueError, simulate_raw, raw, stc, trans, src, sphere,
-                  interp='foo', use_cps=True)
+                  interp='foo')
     pytest.raises(TypeError, simulate_raw, raw, stc, trans, src, sphere,
-                  head_pos=1., use_cps=True)
+                  head_pos=1.)
     pytest.raises(RuntimeError, simulate_raw, raw, stc, trans, src, sphere,
-                  head_pos=pos_fname, use_cps=True)  # ends up with t>t_end
+                  head_pos=pos_fname)  # ends up with t>t_end
     head_pos_sim_err = deepcopy(head_pos_sim)
     head_pos_sim_err[-1.] = head_pos_sim_err[1.]  # negative time
     pytest.raises(RuntimeError, simulate_raw, raw, stc, trans, src, sphere,
-                  head_pos=head_pos_sim_err, use_cps=True)
+                  head_pos=head_pos_sim_err)
     raw_bad = raw.copy()
     raw_bad.info['dig'] = None
     pytest.raises(RuntimeError, simulate_raw, raw_bad, stc, trans, src, sphere,
-                  blink=True, use_cps=True)
+                  blink=True)
 
 
 @pytest.mark.slowtest
@@ -238,10 +240,9 @@ def test_simulate_raw_bem():
     vertices = [s['vertno'] for s in src]
     stc = SourceEstimate(np.eye(sum(len(v) for v in vertices)), vertices,
                          0, 1. / raw.info['sfreq'])
-    raw_sim_sph = simulate_raw(raw, stc, trans, src, sphere, cov=None,
-                               use_cps=True)
+    raw_sim_sph = simulate_raw(raw, stc, trans, src, sphere, cov=None)
     raw_sim_bem = simulate_raw(raw, stc, trans, src, bem_fname, cov=None,
-                               n_jobs=2, use_cps=True)
+                               n_jobs=2)
     # some components (especially radial) might not match that well,
     # so just make sure that most components have high correlation
     assert_array_equal(raw_sim_sph.ch_names, raw_sim_bem.ch_names)
@@ -335,11 +336,10 @@ def test_simulate_raw_chpi():
     stc = _make_stc(raw, src)
     # simulate data with cHPI on
     raw_sim = simulate_raw(raw, stc, None, src, sphere, cov=None, chpi=False,
-                           interp='zero', use_cps=True)
+                           interp='zero')
     # need to trim extra samples off this one
     raw_chpi = simulate_raw(raw, stc, None, src, sphere, cov=None, chpi=True,
-                            head_pos=pos_fname, interp='zero',
-                            use_cps=True)
+                            head_pos=pos_fname, interp='zero')
     # test cHPI indication
     hpi_freqs, hpi_pick, hpi_ons = _get_hpi_info(raw.info)
     assert_allclose(raw_sim[hpi_pick][0], 0.)