Fix events sorting bug and clarify some documentation

doc/index.rst

@@ -173,7 +173,7 @@ File Mapping
 ------------
 
 MNE-HCP supports a low level file mapping that allows for quick compilations
-of sets of files for a given subejct and data context.
+of sets of files for a given subject and data context.
 This is done in :func:`hcp.io.file_mapping.get_file_paths`, think of it as a
 file name synthesizer that takes certain data description parameters as inputs
 and lists all corresponding files.
@@ -196,7 +196,7 @@ to be accessed are known in advance.
 Gotchas
 =======
 
-Native coordinates and resulting plotting and processing peculartities
+Native coordinates and resulting plotting and processing peculiarities
 ----------------------------------------------------------------------
 
 The HCP for MEG provides coregistration information for native BTI/4D
@@ -215,7 +215,7 @@ This is not a bug.
 
 2. All channel names and positions are native. Topographic plotting might not
 work as as expected. First of all, the layout file is not recognized. Second,
-the coordinates are not regonized as native ones, eventually rotating and
+the coordinates are not recognized as native ones, eventually rotating and
 distorting the graphical display. To fix this, either a proper layout can be
 computed with :func:`hcp.viz.make_hcp_bti_layout`.
 Or the conversion to MNE can also be
@@ -266,7 +266,7 @@ Convenience functions
 NNE-HCP includes convenience functions that help setting up directory and file layouts
 expected by MNE-Python.
 
-:func:`hcp.make_mne_anatomy` will produce an MNE and Freesurfer compatible directory layout and will create the following outputs by default, mostly using sympbolic links:
+:func:`hcp.make_mne_anatomy` will produce an MNE and Freesurfer compatible directory layout and will create the following outputs by default, mostly using symbolic links:
 
 .. code-block:: bash
 

examples/make_mne_anatomy.py

@@ -6,9 +6,9 @@
 =====================================
 
 The HCP anatomy data needed for MNE based source analysis
-Live in different places and subfolders.
+live in different places and subfolders.
 MNE-HCP has a convenience function that exctracts the required
-files and information and creates a "subjects_dir" as known from
+files and information and creates a "subjects_dir" just as in
 MNE and freesurfer.
 
 The most essential contents are

examples/plot_evoked_data.py

@@ -5,7 +5,7 @@
 Visualize evoked data
 =====================
 
-Here we'll generate some attention for pecularities of visualizing the
+Here we'll generate some attention for peculiarities of visualizing the
 HCP evoked outputs using MNE plotting functions.
 """
 # Author: Denis A. Engemann
@@ -35,9 +35,9 @@
     if not evoked.comment == 'Wrkmem_LM-TIM-face_BT-diff_MODE-mag':
         continue
 ##############################################################################
-# In order to plot topographic patterns we need to transform the sensor
-# positions to MNE coordinates one a copy of the data.
-# We're not using this these transformed data for any source analyses.
+# In order to plot topographic patterns, we need to transform the sensor
+# positions to MNE coordinates using a copy of the data.
+# We're not using these transformed data for any source analyses.
 # These take place in native BTI/4D coordinates.
 
 evoked_viz = evoked.copy()

tutorials/plot_compute_evoked_inverse_solution.py

@@ -5,8 +5,8 @@
 Compute inverse solution for evoked data
 ========================================
 
-Here we'll use our knowledge from the other examples and tutorials
-to compute an inverse solution and apply it on event related fields.
+Here, we'll use our knowledge from the other examples and tutorials
+to compute an inverse solution and apply it to event related fields.
 """
 # Author: Denis A. Engemann
 # License: BSD 3 clause
@@ -17,7 +17,7 @@
 from hcp import preprocessing as preproc
 
 ##############################################################################
-# we assume our data is inside a designated folder under $HOME
+# We assume our data is inside a designated folder under $HOME
 storage_dir = op.expanduser('~/mne-hcp-data')
 hcp_path = op.join(storage_dir, 'HCP')
 recordings_path = op.join(storage_dir, 'hcp-meg')
@@ -27,11 +27,11 @@
 run_index = 0
 
 ##############################################################################
-# We're reading the evoked data.
+# First, read the evoked data.
 # These are the same as in :ref:`tut_plot_evoked`
 
 hcp_evokeds = hcp.read_evokeds(onset='stim', subject=subject,
-                                   data_type=data_type, hcp_path=hcp_path)
+                               data_type=data_type, hcp_path=hcp_path)
 for evoked in hcp_evokeds:
     if not evoked.comment == 'Wrkmem_LM-TIM-face_BT-diff_MODE-mag':
         continue
@@ -43,15 +43,15 @@
 src_outputs = hcp.anatomy.compute_forward_stack(
     subject=subject, subjects_dir=subjects_dir,
     hcp_path=hcp_path, recordings_path=recordings_path,
-    # speed up computations here. Setting `add_dist` to True may improve the
-    # accuracy.
+    # Speed up computations here by setting `add_dist` to False.
+    # Setting `add_dist` to True will slightly improve forward soln. accuracy
     src_params=dict(add_dist=False),
     info_from=dict(data_type=data_type, run_index=run_index))
 
 fwd = src_outputs['fwd']
 
 ##############################################################################
-# Now we can compute the noise covariance. For this purpose we will apply
+# Now, we can compute the noise covariance. For this purpose, we will apply
 # the same filtering as was used for the computations of the ERF in the first
 # place. See also :ref:`tut_reproduce_erf`.
 
@@ -70,13 +70,13 @@
 
 ##############################################################################
 # Note that using the empty room noise covariance will inflate the SNR of the
-# evkoked and renders comparisons  to `baseline` rather uninformative.
+# evoked and renders comparisons  to `baseline` rather uninformative.
 noise_cov = mne.compute_raw_covariance(raw_noise, method='empirical')
 
 
 ##############################################################################
-# Now we assemble the inverse operator, project the data and show the results
-# on the `fsaverage` surface, the freesurfer average brain.
+# Now we assemble the inverse operator, project the data, and show the results
+# on the `fsaverage` surface (i.e., the freesurfer average brain).
 
 inv_op = mne.minimum_norm.make_inverse_operator(
     evoked.info, fwd, noise_cov=noise_cov)

tutorials/plot_compute_forward.py

@@ -24,24 +24,24 @@
 subject = '105923'  # our test subject
 
 ##############################################################################
-# and we assume to have the downloaded data, the MNE/freesurfer style
+# and we assume that we have the downloaded data, the MNE/freesurfer style
 # anatomy directory, and the MNE style MEG directory.
-# these can be obtained from :func:`make_mne_anatomy`.
+# These can be obtained from :func:`make_mne_anatomy`.
 # See also :ref:`tut_make_anatomy`.
 
 ##############################################################################
-# first we read the coregistration.
+# First, we read the coregistration.
 
 head_mri_t = mne.read_trans(
     op.join(recordings_path, subject, '{}-head_mri-trans.fif'.format(
             subject)))
 
 ##############################################################################
-# Now we can setup our source model.
+# Now, we can setup our source model.
 # Note that spacing has to be set to 'all' since no common MNE resampling
 # scheme has been employed in the HCP pipelines.
-# Since this will take very long time to compute and at this point no other
-# decimation scheme is available inside MNE, we will compute the source
+# Since this will take a very long time to compute (and at this point no other
+# decimation scheme is available inside MNE), we will compute the source
 # space on fsaverage, the freesurfer average brain, and morph it onto
 # the subject's native space. With `oct6` we have ~8000 dipole locations.
 
@@ -55,8 +55,9 @@
     src_fsaverage, subject, subjects_dir=subjects_dir)
 
 ##############################################################################
-# For the same reason `ico` has to be set to `None` when computing the bem.
-# The headshape is not computed with MNE and has a none standard configuration.
+# For the same reason, `ico` has to be set to `None` when computing the bem.
+# The headshape is not computed with MNE and has a `None` standard
+# configuration.
 
 bems = mne.make_bem_model(subject, conductivity=(0.3,),
                           subjects_dir=subjects_dir,
@@ -65,7 +66,7 @@
 bem_sol['surfs'][0]['coord_frame'] = 5
 
 ##############################################################################
-# Now we can read the channels that we want to map to the cortical locations.
+# Now, we can read the channels that we want to map to the cortical locations.
 # Then we can compute the forward solution.
 
 info = hcp.read_info(subject=subject, hcp_path=hcp_path, data_type='rest',
@@ -80,7 +81,7 @@
     fwd, projs=None, ch_type='mag', mode='fixed', exclude=[], verbose=None)
 
 ##############################################################################
-# we display sensitivity map on the original surface with little smoothing
+# We display sensitivity map on the original surface with little smoothing
 # and admire the expected curvature-driven sensitivity pattern.
 
 mag_map = mag_map.to_original_src(src_fsaverage, subjects_dir=subjects_dir)

tutorials/plot_reference_correction.py

@@ -3,7 +3,7 @@
 Apply reference channel correction
 ==================================
 
-Apply reference channels and see what happens.
+Apply reference channels and see how it improves noisy channels.
 """
 # Author: Denis A. Engemann
 # License: BSD 3 clause
@@ -30,7 +30,7 @@
 
 
 ###############################################################################
-# Then we define a spectral plotter for convenience
+# Then we define a spectral plotting helper function for convenience
 
 
 def plot_psd(X, label, Fs, NFFT, color=None):
@@ -44,10 +44,10 @@ def plot_psd(X, label, Fs, NFFT, color=None):
 
 
 ###############################################################################
-# Now we read in the data
+# Now, we read in the data.
 #
-# Then we plot the power spectrum of the MEG and reference channels,
-# apply the reference correction and add the resulting cleaned MEG channels
+# After that, we plot the power spectrum of the MEG and reference channels,
+# apply the reference correction, and add the resulting cleaned MEG channels
 # to our comparison.
 
 
@@ -63,12 +63,12 @@ def plot_psd(X, label, Fs, NFFT, color=None):
 # put single channel aside for comparison later
 chan1 = raw[meg_picks[0]][0]
 
-# add some plotting parameter
-decim_fit = 100  # we lean a purely spatial model, we don't need all samples
+# add some plotting parameters
+decim_fit = 100  # we learn a purely spatial model, we don't need all samples
 decim_show = 10  # we can make plotting faster
 n_fft = 2 ** 15  # let's use long windows to see low frequencies
 
-# we put aside the time series for later plotting
+# we put aside the time series for plotting later
 x_meg = raw[meg_picks][0][:, ::decim_show].mean(0)
 x_meg_ref = raw[ref_picks][0][:, ::decim_show].mean(0)
 
@@ -96,11 +96,13 @@ def plot_psd(X, label, Fs, NFFT, color=None):
 plt.show()
 
 ###############################################################################
-# We can see that the ref correction removes low frequencies which is expected
+# We can see that the ref correction removes low frequencies, which is
+# expected.
 
 
 ###############################################################################
-# By comparing single channel time series we can also see the detrending effect
+# By comparing single channel time series, we can also see the detrending
+# effect.
 
 chan1c = raw[meg_picks[0]][0]
 ch_name = raw.ch_names[meg_picks[0]]

tutorials/plot_reproduce_erf.py

@@ -6,10 +6,11 @@
 Computing ERFs from HCP
 =======================
 
-In this tutorial we compare different ways of arriving at event related
-fields (ERF) starting from different HCP outputs. We will first reprocess
-the HCP dat from scratch, then read the preprocessed epochs, finally
-read the ERF files. Subsequently we will compare these outputs.
+In this tutorial, we compare different ways of arriving at the exact same event
+related fields (ERFs) but starting from different points on the HCP output
+pipeline. We will first reprocess the HCP dat from scratch, then read the
+preprocessed epochs, and finally read the fully processed ERF files.
+Subsequently, we will compare these outputs.
 """
 # Author: Denis A. Engemann
 # License: BSD 3 clause
@@ -36,14 +37,16 @@
 # We first reprocess the data from scratch
 #
 # That is, almost from scratch. We're relying on the ICA solutions and
-# data annotations.
+# data annotations from the preprocessed data.
 #
 # In order to arrive at the final ERF we need to pool over two runs.
-# for each run we need to read the raw data, all annotations, apply
+# For each run we need to read the raw data and all annotations, apply
 # the reference sensor compensation, the ICA, bandpass filter, baseline
-# correction and decimation (downsampling)
+# correction and decimation (downsampling).
 
-# these values are looked up from the HCP manual
+# These values were looked up from the HCP manual online.
+# See pg. 133 of the HCP_S500+MEG2 Release Reference Manual:
+# http://www.humanconnectome.org/documentation/S500/HCP_S500+MEG2_Release_Reference_Manual.pdf
 tmin, tmax = -1.5, 2.5
 decim = 4
 event_id = dict(face=1)
@@ -56,15 +59,15 @@
     trial_infos.append(trial_info)
 
 
-# trial_info is a dict
-# it contains a 'comments' vector that maps on the columns of 'codes'
-# 'codes is a matrix with its length corresponding to the number of trials
+# trial_info is a dict that contains a 'comments' vector that maps on the
+# columns of 'codes'
+# 'codes' is a matrix with its length corresponding to the number of trials
 print(trial_info['stim']['comments'][:10])  # which column?
 print(set(trial_info['stim']['codes'][:, 3]))  # check values
 
-# so according to this we need to use the column 7 (index 6)
-# for the time sample and column 4 (index 3) to get the image types
-# with this information we can construct our event vectors
+# So according to this, we need to use the column 7 (index 6)
+# for the time sample, and column 4 (index 3) to get the image types.
+# With this information, we can construct our event vectors.
 
 all_events = list()
 for trial_info in trial_infos:
@@ -80,7 +83,7 @@
 
     all_events.append(events)
 
-# now we can go ahead
+# Now we can go ahead and process the raw data.
 evokeds = list()
 for run_index, events in zip([0, 1], all_events):
 
@@ -106,16 +109,18 @@
     raw.filter(None, 60, method='iir',
                iir_params=dict(order=4, ftype='butter'), n_jobs=1)
 
-    # read ICA and remove EOG ECG
+    # read ICA and remove EOG/ECG artifacts
     # note that the HCP ICA assumes that bad channels have already been removed
     ica_mat = hcp.read_ica(run_index=run_index, **hcp_params)
 
     # We will select the brain ICs only
     exclude = annots['ica']['ecg_eog_ic']
     preproc.apply_ica_hcp(raw, ica_mat=ica_mat, exclude=exclude)
 
+    # Sort events by timestamp
+    events = events[np.argsort(events[:, 0], axis=0)]
+
     # now we can epoch
-    events = np.sort(events, 0)
     epochs = mne.Epochs(raw, events=events[events[:, 2] == 1],
                         event_id=event_id, tmin=tmin, tmax=tmax,
                         reject=None, baseline=baseline, decim=decim,
@@ -128,7 +133,8 @@
     del epochs, raw
 
 ##############################################################################
-# Now we can compute the same ERF based on the preprocessed epochs
+# Now we can compute the same ERF based on the epochs preprocessed by the HCP
+# folks.
 #
 # These are obtained from the 'tmegpreproc' pipeline.
 # Things are pythonized and simplified however, so
@@ -164,12 +170,12 @@
 
 
 ##############################################################################
-# Finally we can read the actual official ERF file
+# Finally, we can read the (completely processed) official ERF file
 #
 # These are obtained from the 'eravg' pipelines.
 # We read the matlab file, MNE-HCP is doing some conversions, and then we
 # search our condition of interest. Here we're looking at the image as onset.
-# and we want the average, not the standard deviation.
+# We want the average, not the standard deviation.
 
 evoked_hcp = None
 hcp_evokeds = hcp.read_evokeds(onset='stim', **hcp_params)