Fix to sklearn imports. (#3773)

* Fix to sklearn imports.

* Fixes.

* Format fix.

* setup.cfg

* Except.

* Fix to nibabel import.

circle.yml

@@ -34,6 +34,7 @@ dependencies:
     - ls -al /home/ubuntu/miniconda;
     - ls -al /home/ubuntu/miniconda/bin;
     - echo $PATH;
+    - echo $CIRCLE_BRANCH
     - which python;
     - which pip;
     - git clone https://github.com/sphinx-gallery/sphinx-gallery.git;
@@ -43,7 +44,7 @@ dependencies:
 
   override:
     - cd /home/ubuntu/mne-python && python setup.py develop;
-    - if [ "$CIRCLE_BRANCH" == "master" ]; then
+    - if [ "$CIRCLE_BRANCH" == "master" ] || [ "$CIRCLE_BRANCH" == "maint/0.12" ]; then
         mkdir -p ~/mne_data;
         python -c "import mne; mne.datasets._download_all_example_data()";
       fi
@@ -56,13 +57,15 @@ dependencies:
 
 test:
   override:
-    - if [ "$CIRCLE_BRANCH" == "master" ]; then
+    - if [ "$CIRCLE_BRANCH" == "master" ] || [ "$CIRCLE_BRANCH" == "maint/0.12" ]; then
         make test-doc;
       else
         cd doc && make html_dev-noplot;
       fi
     - if [ "$CIRCLE_BRANCH" == "master" ]; then cd doc && make html_dev; fi:
           timeout: 1500
+    - if [ "$CIRCLE_BRANCH" == "maint/0.12" ]; then cd doc && make html_stable; fi:
+          timeout: 1500
 
 general:
   # branches:
@@ -81,3 +84,10 @@ deployment:
       - cd ../mne-tools.github.io && git checkout master && git pull origin master
       - cd doc/_build/html && cp -rf * ~/mne-tools.github.io/dev
       - cd ../mne-tools.github.io && git add -A && git commit -m 'Automated update of dev docs.' && git push origin master
+    #branch: maint/0.12
+    #commands:
+    #  - git config --global user.email "circle@mne.com"
+    #  - git config --global user.name "Circle Ci"
+    #  - cd ../mne-tools.github.io && git checkout master && git pull origin master
+    #  - cd doc/_build/html_stable && cp -rf * ~/mne-tools.github.io/stable
+    #  - cd ../mne-tools.github.io && git add -A && git commit -m 'Automated update of stable docs.' && git push origin master

doc/whats_new.rst

@@ -173,6 +173,8 @@ API
 
     - The ``events`` parameter of :func:`mne.epochs.EpochsArray` is set by default to chronological time-samples and event values to 1, by `Jean-Remi King`_
 
+    - Deprecated ``axis`` parameter from :func:`mne.viz.plot_topomap` (use ``axes`` instead) by `Jaakko Leppakangas`_
+
 Authors
 ~~~~~~~
 

examples/time_frequency/plot_time_frequency_simulated.py

@@ -101,7 +101,7 @@ for width in (0.7, 3.0):
     power = tfr_stockwell(epochs, fmin=fmin, fmax=fmax, width=width)
     power.plot([0], baseline=(0., 0.1), mode='mean',
                title='Sim: Using S transform, width '
-                     '= {:0.1f}'.format(width), show=True)
+                     '= {0:0.1f}'.format(width), show=True)
 
 # #############################################################################
 # Finally, compare to morlet wavelet

mne/channels/channels.py

@@ -220,7 +220,7 @@ class SetChannelsMixin(object):
         n_zero = np.sum(np.sum(np.abs(pos), axis=1) == 0)
         if n_zero > 1:  # XXX some systems have origin (0, 0, 0)
             raise ValueError('Could not extract channel positions for '
-                             '{} channels'.format(n_zero))
+                             '{0} channels'.format(n_zero))
         return pos
 
     def _set_channel_positions(self, pos, names):
@@ -771,7 +771,7 @@ def read_ch_connectivity(fname, picks=None):
                 break
         else:
             raise ValueError('I do not know about this neighbor '
-                             'template: "{}"'.format(fname))
+                             'template: "{0}"'.format(fname))
 
         fname = op.join(templates_dir, fname)
 
@@ -783,8 +783,8 @@ def read_ch_connectivity(fname, picks=None):
     if picks is not None:
         if max(picks) >= len(ch_names):
             raise ValueError('The picks must be compatible with '
-                             'channels. Found a pick ({}) which exceeds '
-                             'the channel range ({})'
+                             'channels. Found a pick ({0}) which exceeds '
+                             'the channel range ({1})'
                              .format(max(picks), len(ch_names)))
     connectivity = _ch_neighbor_connectivity(ch_names, neighbors)
     if picks is not None:

mne/chpi.py

@@ -127,6 +127,9 @@ def read_head_pos(fname):
     _check_fname(fname, must_exist=True, overwrite=True)
     data = np.loadtxt(fname, skiprows=1)  # first line is header, skip it
     data.shape = (-1, 10)  # ensure it's the right size even if empty
+    if np.isnan(data).any():  # make sure we didn't do something dumb
+        raise RuntimeError('positions could not be read properly from %s'
+                           % fname)
     return data
 
 

mne/datasets/utils.py

@@ -167,7 +167,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True,
     path = _get_path(path, key, name)
     # To update the testing or misc dataset, push commits, then make a new
     # release on GitHub. Then update the "releases" variable:
-    releases = dict(testing='0.19', misc='0.1')
+    releases = dict(testing='0.20', misc='0.1')
     # And also update the "hashes['testing']" variable below.
 
     # To update any other dataset, update the data archive itself (upload
@@ -211,7 +211,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True,
         sample='1d5da3a809fded1ef5734444ab5bf857',
         somato='f3e3a8441477bb5bacae1d0c6e0964fb',
         spm='f61041e3f3f2ba0def8a2ca71592cc41',
-        testing='77b2a435d80adb23cbe7e19144e7bc47'
+        testing='625cb9c1c62de03bd2bc6a3ca9a6f2d8'
     )
     folder_origs = dict(  # not listed means None
         misc='mne-misc-data-%s' % releases['misc'],

mne/decoding/tests/test_time_gen.py

@@ -59,13 +59,14 @@ def test_generalization_across_time():
     y_4classes = np.hstack((epochs.events[:7, 2], epochs.events[7:, 2] + 1))
     if check_version('sklearn', '0.18'):
         from sklearn.model_selection import (KFold, StratifiedKFold,
-                                             ShuffleSplit, LeaveOneLabelOut)
+                                             ShuffleSplit, LeaveOneGroupOut)
         cv_shuffle = ShuffleSplit()
-        cv = LeaveOneLabelOut()
+
+        cv = LeaveOneGroupOut()
         # XXX we cannot pass any other parameters than X and y to cv.split
         # so we have to build it before hand
         cv_lolo = [(train, test) for train, test in cv.split(
-                   X=y_4classes, y=y_4classes, labels=y_4classes)]
+                   X=y_4classes, y=y_4classes, groups=y_4classes)]
 
         # With sklearn >= 0.17, `clf` can be identified as a regressor, and
         # the scoring metrics can therefore be automatically assigned.

mne/decoding/time_gen.py

@@ -1571,7 +1571,7 @@ def _set_cv(cv, clf=None, X=None, y=None):
         from sklearn.model_selection import (check_cv, StratifiedKFold, KFold)
         if isinstance(cv, (int, np.int)):
             XFold = StratifiedKFold if is_classifier(clf) else KFold
-            cv = XFold(n_folds=cv)
+            cv = XFold(n_splits=cv)
         cv = check_cv(cv=cv, y=y, classifier=is_classifier(clf))
     else:
         from sklearn.cross_validation import (check_cv, StratifiedKFold, KFold)

mne/io/pick.py

@@ -272,7 +272,7 @@ def pick_types(info, meg=True, eeg=False, stim=False, eog=False, ecg=False,
                   ias, syst, seeg, dipole, gof, bio, ecog):
         if not isinstance(param, bool):
             w = ('Parameters for all channel types (with the exception '
-                 'of "meg" and "ref_meg") must be of type bool, not {}.')
+                 'of "meg" and "ref_meg") must be of type bool, not {0}.')
             raise ValueError(w.format(type(param)))
 
     for k in range(nchan):

mne/preprocessing/_fine_cal.py

@@ -0,0 +1,84 @@
+# -*- coding: utf-8 -*-
+# Authors: Eric Larson <larson.eric.d@gmail.com>
+
+# License: BSD (3-clause)
+
+import numpy as np
+
+from ..utils import check_fname, _check_fname
+
+
+def read_fine_calibration(fname):
+    """Read fine calibration information from a .dat file
+
+    The fine calibration typically includes improved sensor locations,
+    calibration coefficients, and gradiometer imbalance information.
+
+    Parameters
+    ----------
+    fname : str
+        The filename.
+
+    Returns
+    -------
+    calibration : dict
+        Fine calibration information.
+    """
+    # Read new sensor locations
+    _check_fname(fname, overwrite=True, must_exist=True)
+    check_fname(fname, 'cal', ('.dat',))
+    ch_names = list()
+    locs = list()
+    imb_cals = list()
+    with open(fname, 'r') as fid:
+        for line in fid:
+            if line[0] in '#\n':
+                continue
+            vals = line.strip().split()
+            if len(vals) not in [14, 16]:
+                raise RuntimeError('Error parsing fine calibration file, '
+                                   'should have 14 or 16 entries per line '
+                                   'but found %s on line:\n%s'
+                                   % (len(vals), line))
+            # `vals` contains channel number
+            ch_name = vals[0]
+            if len(ch_name) in (3, 4):  # heuristic for Neuromag fix
+                try:
+                    ch_name = int(ch_name)
+                except ValueError:  # something other than e.g. 113 or 2642
+                    pass
+                else:
+                    ch_name = 'MEG' + '%04d' % ch_name
+            ch_names.append(ch_name)
+            # (x, y, z), x-norm 3-vec, y-norm 3-vec, z-norm 3-vec
+            locs.append(np.array([float(x) for x in vals[1:13]]))
+            # and 1 or 3 imbalance terms
+            imb_cals.append([float(x) for x in vals[13:]])
+    locs = np.array(locs)
+    return dict(ch_names=ch_names, locs=locs, imb_cals=imb_cals)
+
+
+def write_fine_calibration(fname, calibration):
+    """Write fine calibration information to a .dat file
+
+    Parameters
+    ----------
+    fname : str
+        The filename to write out.
+    calibration : dict
+        Fine calibration information.
+    """
+    _check_fname(fname, overwrite=True)
+    check_fname(fname, 'cal', ('.dat',))
+
+    with open(fname, 'wb') as cal_file:
+        for ci, chan in enumerate(calibration['ch_names']):
+            # Write string containing 1) channel, 2) loc info, 3) calib info
+            # with field widths (e.g., %.6f) chosen to match how Elekta writes
+            # them out
+            cal_line = np.concatenate([calibration['locs'][ci],
+                                       calibration['imb_cals'][ci]]).round(6)
+            cal_str = str(chan) + ' ' + ' '.join(map(lambda x: "%.6f" % x,
+                                                     cal_line))
+
+            cal_file.write((cal_str + '\n').encode('ASCII'))

mne/preprocessing/maxwell.py

@@ -25,7 +25,7 @@ from ..io.write import _generate_meas_id, _date_now
 from ..io import _loc_to_coil_trans, _BaseRaw
 from ..io.pick import pick_types, pick_info, pick_channels
 from ..utils import verbose, logger, _clean_names, warn, _time_mask
-from ..fixes import _get_args, partial
+from ..fixes import _get_args, partial, in1d
 from ..externals.six import string_types
 from ..channels.channels import _get_T1T2_mag_inds
 
@@ -247,11 +247,12 @@ def maxwell_filter(raw, origin='auto', int_order=8, ext_order=3,
     recon_trans = _check_destination(destination, raw.info, head_frame)
     if st_duration is not None:
         st_duration = float(st_duration)
-        if not 0. < st_duration <= raw.times[-1]:
+        if not 0. < st_duration <= raw.times[-1] + 1. / raw.info['sfreq']:
             raise ValueError('st_duration (%0.1fs) must be between 0 and the '
                              'duration of the data (%0.1fs).'
                              % (st_duration, raw.times[-1]))
         st_correlation = float(st_correlation)
+        st_duration = int(round(st_duration * raw.info['sfreq']))
         if not 0. < st_correlation <= 1:
             raise ValueError('st_correlation must be between 0. and 1.')
     if not isinstance(bad_condition, string_types) or \
@@ -279,7 +280,7 @@ def maxwell_filter(raw, origin='auto', int_order=8, ext_order=3,
         raw, add_channels=add_channels)
     del raw
     _remove_meg_projs(raw_sss)  # remove MEG projectors, they won't apply now
-    info, times = raw_sss.info, raw_sss.times
+    info = raw_sss.info
     meg_picks, mag_picks, grad_picks, good_picks, coil_scale, mag_or_fine = \
         _get_mf_picks(info, int_order, ext_order, ignore_ref, mag_scale=100.)
 
@@ -289,11 +290,11 @@ def maxwell_filter(raw, origin='auto', int_order=8, ext_order=3,
     sss_cal = dict()
     if calibration is not None:
         calibration, sss_cal = _update_sensor_geometry(info, calibration,
-                                                       head_frame, ignore_ref)
+                                                       ignore_ref)
         mag_or_fine.fill(True)  # all channels now have some mag-type data
 
     # Determine/check the origin of the expansion
-    origin = _check_origin(origin, raw_sss.info, coord_frame, disp=True)
+    origin = _check_origin(origin, info, coord_frame, disp=True)
     origin.setflags(write=False)
     n_in, n_out = _get_n_moments([int_order, ext_order])
 
@@ -312,7 +313,9 @@ def maxwell_filter(raw, origin='auto', int_order=8, ext_order=3,
         if len(missing) > 0:
             warn('Not all cross-talk channels in raw:\n%s' % missing)
         ctc_picks = pick_channels(ctc_chs,
-                                  [info['ch_names'][c] for c in good_picks])
+                                  [info['ch_names'][c]
+                                   for c in meg_picks[good_picks]])
+        assert len(ctc_picks) == len(good_picks)  # otherwise we errored
         ctc = sss_ctc['decoupler'][ctc_picks][:, ctc_picks]
         # I have no idea why, but MF transposes this for storage..
         sss_ctc['decoupler'] = sss_ctc['decoupler'].T.tocsc()
@@ -322,8 +325,7 @@ def maxwell_filter(raw, origin='auto', int_order=8, ext_order=3,
     #
     # Translate to destination frame (always use non-fine-cal bases)
     #
-    exp = dict(origin=origin, int_order=int_order, ext_order=0,
-               head_frame=head_frame)
+    exp = dict(origin=origin, int_order=int_order, ext_order=0)
     all_coils = _prep_mf_coils(info, ignore_ref)
     S_recon = _trans_sss_basis(exp, all_coils, recon_trans, coil_scale)
     exp['ext_order'] = ext_order
@@ -341,29 +343,33 @@ def maxwell_filter(raw, origin='auto', int_order=8, ext_order=3,
     # Reconstruct raw file object with spatiotemporal processed data
     max_st = dict()
     if st_duration is not None:
-        max_st.update(job=10, subspcorr=st_correlation, buflen=st_duration)
+        max_st.update(job=10, subspcorr=st_correlation,
+                      buflen=st_duration / info['sfreq'])
         logger.info('    Processing data using tSSS with st_duration=%s'
-                    % st_duration)
+                    % max_st['buflen'])
         st_when = 'before' if st_fixed else 'after'  # relative to movecomp
     else:
-        st_duration = min(raw_sss.times[-1], 10.)  # chunk size
+        # st_duration from here on will act like the chunk size
+        st_duration = max(int(round(10. * info['sfreq'])), 1)
         st_correlation = None
         st_when = 'never'
+    st_duration = min(len(raw_sss.times), st_duration)
     del st_fixed
 
-    # Generate time points to break up data into windows
-    chunk_times = np.arange(times[0], times[-1], st_duration)
-    read_lims = raw_sss.time_as_index(chunk_times)
-    len_last_buf = raw_sss.times[-1] - raw_sss.times[read_lims[-1]]
-    if len_last_buf == st_duration:
+    # Generate time points to break up data into equal-length windows
+    read_lims = np.arange(0, len(raw_sss.times) + 1, st_duration)
+    if len(read_lims) == 1:
         read_lims = np.concatenate([read_lims, [len(raw_sss.times)]])
-    else:
-        # len_last_buf < st_dur so fold it into the previous buffer
+    if read_lims[-1] != len(raw_sss.times):
         read_lims[-1] = len(raw_sss.times)
-        if st_correlation is not None:
+        # len_last_buf < st_dur so fold it into the previous buffer
+        if st_correlation is not None and len(read_lims) > 2:
             logger.info('    Spatiotemporal window did not fit evenly into '
                         'raw object. The final %0.2f seconds were lumped '
-                        'onto the previous window.' % len_last_buf)
+                        'onto the previous window.'
+                        % ((read_lims[-1] - read_lims[-2]) / info['sfreq'],))
+    assert len(read_lims) >= 2
+    assert read_lims[0] == 0 and read_lims[-1] == len(raw_sss.times)
 
     #
     # Do the heavy lifting
@@ -373,7 +379,7 @@ def maxwell_filter(raw, origin='auto', int_order=8, ext_order=3,
     # (and transform pos[1] to times)
     head_pos[1] = raw_sss.time_as_index(head_pos[1], use_rounding=True)
     # Compute the first bit of pos_data for cHPI reporting
-    if raw_sss.info['dev_head_t'] is not None and head_pos[0] is not None:
+    if info['dev_head_t'] is not None and head_pos[0] is not None:
         this_pos_quat = np.concatenate([
             rot_to_quat(info['dev_head_t']['trans'][:3, :3]),
             info['dev_head_t']['trans'][:3, 3],
@@ -387,17 +393,21 @@ def maxwell_filter(raw, origin='auto', int_order=8, ext_order=3,
         grad_picks=grad_picks, mag_picks=mag_picks, good_picks=good_picks,
         mag_or_fine=mag_or_fine, bad_condition=bad_condition)
     S_decomp, pS_decomp, reg_moments, n_use_in = _get_this_decomp_trans(
-        info['dev_head_t'])
+        info['dev_head_t'], t=0.)
     reg_moments_0 = reg_moments.copy()
     # Loop through buffer windows of data
+    n_sig = int(np.floor(np.log10(max(len(read_lims), 0)))) + 1
     pl = 's' if len(read_lims) != 2 else ''
     logger.info('    Processing %s data chunk%s of (at least) %0.1f sec'
-                % (len(read_lims) - 1, pl, st_duration))
-    for start, stop in zip(read_lims[:-1], read_lims[1:]):
-        t_str = '% 8.2f - % 8.2f sec' % tuple(raw_sss.times[[start, stop - 1]])
+                % (len(read_lims) - 1, pl, st_duration / info['sfreq']))
+    for ii, (start, stop) in enumerate(zip(read_lims[:-1], read_lims[1:])):
+        rel_times = raw_sss.times[start:stop]
+        t_str = '%8.3f - %8.3f sec' % tuple(rel_times[[0, -1]])
+        t_str += ('(#%d/%d)'
+                  % (ii + 1, len(read_lims) - 1)).rjust(2 * n_sig + 5)
 
         # Get original data
-        orig_data = raw_sss._data[good_picks, start:stop]
+        orig_data = raw_sss._data[meg_picks[good_picks], start:stop]
         # This could just be np.empty if not st_only, but shouldn't be slow
         # this way so might as well just always take the original data
         out_meg_data = raw_sss._data[meg_picks, start:stop]
@@ -422,7 +432,7 @@ def maxwell_filter(raw, origin='auto', int_order=8, ext_order=3,
                 if avg_trans is not None:
                     # if doing movecomp
                     S_decomp_st, pS_decomp_st, _, n_use_in_st = \
-                        _get_this_decomp_trans(avg_trans, verbose=False)
+                        _get_this_decomp_trans(avg_trans, t=rel_times[0])
                 else:
                     S_decomp_st, pS_decomp_st = S_decomp, pS_decomp
                     n_use_in_st = n_use_in
@@ -446,7 +456,7 @@ def maxwell_filter(raw, origin='auto', int_order=8, ext_order=3,
                 # previous interval)
                 if trans is not None:
                     S_decomp, pS_decomp, reg_moments, n_use_in = \
-                        _get_this_decomp_trans(trans, verbose=False)
+                        _get_this_decomp_trans(trans, t=rel_times[rel_start])
 
                 # Determine multipole moments for this interval
                 mm_in = np.dot(pS_decomp[:n_use_in],
@@ -477,7 +487,7 @@ def maxwell_filter(raw, origin='auto', int_order=8, ext_order=3,
         if st_when == 'after':
             _do_tSSS(out_meg_data, orig_in_data, resid, st_correlation,
                      n_positions, t_str)
-        else:
+        elif st_when == 'never' and head_pos[0] is not None:
             pl = 's' if n_positions > 1 else ''
             logger.info('        Used % 2d head position%s for %s'
                         % (n_positions, pl, t_str))
@@ -620,10 +630,10 @@ def _do_tSSS(clean_data, orig_in_data, resid, st_correlation,
     np.asarray_chkfinite(resid)
     t_proj = _overlap_projector(orig_in_data, resid, st_correlation)
     # Apply projector according to Eq. 12 in [2]_
-    msg = ('        Projecting % 2d intersecting tSSS components '
+    msg = ('        Projecting %2d intersecting tSSS components '
            'for %s' % (t_proj.shape[1], t_str))
     if n_positions > 1:
-        msg += ' (across % 2d positions)' % n_positions
+        msg += ' (across %2d positions)' % n_positions
     logger.info(msg)
     clean_data -= np.dot(np.dot(clean_data, t_proj), t_proj.T)
 
@@ -701,33 +711,23 @@ def _check_pos(pos, head_frame, raw, st_fixed, sfreq):
     return pos
 
 
-@verbose
 def _get_decomp(trans, all_coils, cal, regularize, exp, ignore_ref,
                 coil_scale, grad_picks, mag_picks, good_picks, mag_or_fine,
-                bad_condition, verbose=None):
-    """Helper to get a decomposition matrix"""
+                bad_condition, t):
+    """Helper to get a decomposition matrix and pseudoinverse matrices"""
     #
     # Fine calibration processing (point-like magnetometers and calib. coeffs)
     #
-    S_decomp = _trans_sss_basis(exp, all_coils, trans, coil_scale)
-    if cal is not None:
-        # Compute point-like mags to incorporate gradiometer imbalance
-        cal['grad_cals'] = _sss_basis_point(exp, trans, cal, ignore_ref)
-        # Add point like magnetometer data to bases.
-        S_decomp[grad_picks, :] += cal['grad_cals']
-        # Scale magnetometers by calibration coefficient
-        S_decomp[mag_picks, :] /= cal['mag_cals']
-        # We need to be careful about KIT gradiometers
-    S_decomp = S_decomp[good_picks]
+    S_decomp = _get_s_decomp(exp, all_coils, trans, coil_scale, cal,
+                             ignore_ref, grad_picks, mag_picks, good_picks)
 
     #
     # Regularization
     #
-    reg_moments, n_use_in = _regularize(regularize, exp, S_decomp, mag_or_fine)
-    S_decomp = S_decomp.take(reg_moments, axis=1)
+    S_decomp, pS_decomp, sing, reg_moments, n_use_in = _regularize(
+        regularize, exp, S_decomp, mag_or_fine, t=t)
 
     # Pseudo-inverse of total multipolar moment basis set (Part of Eq. 37)
-    pS_decomp, sing = _col_norm_pinv(S_decomp.copy())
     cond = sing[0] / sing[-1]
     logger.debug('    Decomposition matrix condition: %0.1f' % cond)
     if bad_condition != 'ignore' and cond >= 1000.:
@@ -743,13 +743,31 @@ def _get_decomp(trans, all_coils, cal, regularize, exp, ignore_ref,
     return S_decomp, pS_decomp, reg_moments, n_use_in
 
 
-def _regularize(regularize, exp, S_decomp, mag_or_fine):
+def _get_s_decomp(exp, all_coils, trans, coil_scale, cal, ignore_ref,
+                  grad_picks, mag_picks, good_picks):
+    """Helper to get S_decomp"""
+    S_decomp = _trans_sss_basis(exp, all_coils, trans, coil_scale)
+    if cal is not None:
+        # Compute point-like mags to incorporate gradiometer imbalance
+        grad_cals = _sss_basis_point(exp, trans, cal, ignore_ref)
+        # Add point like magnetometer data to bases.
+        S_decomp[grad_picks, :] += grad_cals
+        # Scale magnetometers by calibration coefficient
+        S_decomp[mag_picks, :] /= cal['mag_cals']
+        # We need to be careful about KIT gradiometers
+    S_decomp = S_decomp[good_picks]
+    return S_decomp
+
+
+@verbose
+def _regularize(regularize, exp, S_decomp, mag_or_fine, t, verbose=None):
     """Regularize a decomposition matrix"""
     # ALWAYS regularize the out components according to norm, since
     # gradiometer-only setups (e.g., KIT) can have zero first-order
     # components
     int_order, ext_order = exp['int_order'], exp['ext_order']
     n_in, n_out = _get_n_moments([int_order, ext_order])
+    t_str = '%8.3f' % t
     if regularize is not None:  # regularize='in'
         logger.info('    Computing regularization')
         in_removes, out_removes = _regularize_in(
@@ -762,11 +780,15 @@ def _regularize(regularize, exp, S_decomp, mag_or_fine):
                                    out_removes)
     n_use_in = len(reg_in_moments)
     n_use_out = len(reg_out_moments)
-    if regularize is not None or n_use_out != n_out:
-        logger.info('        Using %s/%s inside and %s/%s outside harmonic '
-                    'components' % (n_use_in, n_in, n_use_out, n_out))
     reg_moments = np.concatenate((reg_in_moments, reg_out_moments))
-    return reg_moments, n_use_in
+    S_decomp = S_decomp.take(reg_moments, axis=1)
+    pS_decomp, sing = _col_norm_pinv(S_decomp.copy())
+    if regularize is not None or n_use_out != n_out:
+        logger.info('        Using %s/%s harmonic components for %s  '
+                    '(%s/%s in, %s/%s out)'
+                    % (n_use_in + n_use_out, n_in + n_out, t_str,
+                       n_use_in, n_in, n_use_out, n_out))
+    return S_decomp, pS_decomp, sing, reg_moments, n_use_in
 
 
 def _get_mf_picks(info, int_order, ext_order, ignore_ref=False,
@@ -1618,45 +1640,63 @@ def _read_fine_cal(fine_cal):
     return cal_chs, cal_ch_numbers
 
 
-def _update_sensor_geometry(info, fine_cal, head_frame, ignore_ref):
+def _update_sensor_geometry(info, fine_cal, ignore_ref):
     """Helper to replace sensor geometry information and reorder cal_chs"""
+    from ._fine_cal import read_fine_calibration
     logger.info('    Using fine calibration %s' % op.basename(fine_cal))
-    cal_chs, cal_ch_numbers = _read_fine_cal(fine_cal)
+    fine_cal = read_fine_calibration(fine_cal)  # filename -> dict
+    ch_names = _clean_names(info['ch_names'], remove_whitespace=True)
+    info_order = pick_channels(ch_names, fine_cal['ch_names'])
+    meg_picks = pick_types(info, meg=True, exclude=[])
+    if len(set(info_order) - set(meg_picks)) != 0:
+        # this should never happen
+        raise RuntimeError('Found channels in cal file that are not marked '
+                           'as MEG channels in the data file')
+    if len(info_order) != len(meg_picks):
+        raise RuntimeError(
+            'Not all MEG channels found in fine calibration file, missing:\n%s'
+            % sorted(list(set(ch_names[pick] for pick in meg_picks) -
+                          set(fine_cal['ch_names']))))
+    rev_order = np.argsort(info_order)
+    rev_grad = rev_order[in1d(meg_picks,
+                              pick_types(info, meg='grad', exclude=()))]
+    rev_mag = rev_order[in1d(meg_picks,
+                             pick_types(info, meg='mag', exclude=()))]
 
-    # Check that we ended up with correct channels
-    meg_info = pick_info(info, pick_types(info, meg=True, exclude=[]))
-    clean_meg_names = _clean_names(meg_info['ch_names'],
-                                   remove_whitespace=True)
-    cal_names = [c['ch_name'] for c in cal_chs]
-    order = pick_channels(cal_names, clean_meg_names)
-    if meg_info['nchan'] != len(order):
-        raise RuntimeError('Not all MEG channels found in fine calibration '
-                           'file, missing:\n%s'
-                           % sorted(list(set(clean_meg_names) -
-                                         set(cal_names))))
-    # ensure they're ordered like our data
-    cal_chs = [cal_chs[ii] for ii in order]
+    # Determine gradiometer imbalances and magnetometer calibrations
+    grad_imbalances = np.array([fine_cal['imb_cals'][ri] for ri in rev_grad]).T
+    if grad_imbalances.shape[0] not in [1, 3]:
+        raise ValueError('Must have 1 (x) or 3 (x, y, z) point-like ' +
+                         'magnetometers. Currently have %i' %
+                         grad_imbalances.shape[0])
+    mag_cals = np.array([fine_cal['imb_cals'][ri] for ri in rev_mag])
+    del rev_order, rev_grad, rev_mag
+    # Now let's actually construct our point-like adjustment coils for grads
+    grad_coilsets = _get_grad_point_coilsets(
+        info, n_types=len(grad_imbalances), ignore_ref=ignore_ref)
+    calibration = dict(grad_imbalances=grad_imbalances,
+                       grad_coilsets=grad_coilsets, mag_cals=mag_cals)
 
     # Replace sensor locations (and track differences) for fine calibration
-    ang_shift = np.zeros((len(cal_chs), 3))
+    ang_shift = np.zeros((len(fine_cal['ch_names']), 3))
     used = np.zeros(len(info['chs']), bool)
     cal_corrs = list()
-    coil_types = list()
-    grad_picks = pick_types(meg_info, meg='grad')
+    cal_chans = list()
+    grad_picks = pick_types(info, meg='grad', exclude=())
     adjust_logged = False
-    clean_info_names = _clean_names(info['ch_names'], remove_whitespace=True)
-    for ci, cal_ch in enumerate(cal_chs):
-        idx = clean_info_names.index(cal_ch['ch_name'])
-        assert not used[idx]
-        used[idx] = True
-        info_ch = info['chs'][idx]
-        coil_types.append(info_ch['coil_type'])
+    for ci, info_idx in enumerate(info_order):
+        assert ch_names[info_idx] == fine_cal['ch_names'][ci]
+        assert not used[info_idx]
+        used[info_idx] = True
+        info_ch = info['chs'][info_idx]
+        ch_num = int(fine_cal['ch_names'][ci].lstrip('MEG').lstrip('0'))
+        cal_chans.append([ch_num, info_ch['coil_type']])
 
         # Some .dat files might only rotate EZ, so we must check first that
         # EX and EY are orthogonal to EZ. If not, we find the rotation between
         # the original and fine-cal ez, and rotate EX and EY accordingly:
         ch_coil_rot = _loc_to_coil_trans(info_ch['loc'])[:3, :3]
-        cal_loc = cal_ch['loc'].copy()
+        cal_loc = fine_cal['locs'][ci].copy()
         cal_coil_rot = _loc_to_coil_trans(cal_loc)[:3, :3]
         if np.max([np.abs(np.dot(cal_coil_rot[:, ii], cal_coil_rot[:, 2]))
                    for ii in range(2)]) > 1e-6:  # X or Y not orthogonal
@@ -1669,51 +1709,34 @@ def _update_sensor_geometry(info, fine_cal, head_frame, ignore_ref):
             cal_loc[3:] = np.dot(this_trans, ch_coil_rot).T.ravel()
 
         # calculate shift angle
-        v1 = _loc_to_coil_trans(cal_ch['loc'])[:3, :3]
+        v1 = _loc_to_coil_trans(cal_loc)[:3, :3]
         _normalize_vectors(v1)
         v2 = _loc_to_coil_trans(info_ch['loc'])[:3, :3]
         _normalize_vectors(v2)
         ang_shift[ci] = np.sum(v1 * v2, axis=0)
-        if idx in grad_picks:
-            extra = [1., cal_ch['calib_coeff'][0]]
+        if info_idx in grad_picks:
+            extra = [1., fine_cal['imb_cals'][ci][0]]
         else:
-            extra = [cal_ch['calib_coeff'][0], 0.]
+            extra = [fine_cal['imb_cals'][ci][0], 0.]
         cal_corrs.append(np.concatenate([extra, cal_loc]))
         # Adjust channel normal orientations with those from fine calibration
         # Channel positions are not changed
         info_ch['loc'][3:] = cal_loc[3:]
-        assert (info_ch['coord_frame'] == cal_ch['coord_frame'] ==
-                FIFF.FIFFV_COORD_DEVICE)
-    cal_chans = [[sc, ct] for sc, ct in zip(cal_ch_numbers, coil_types)]
+        assert (info_ch['coord_frame'] == FIFF.FIFFV_COORD_DEVICE)
+    assert used[meg_picks].all()
+    assert not used[np.setdiff1d(np.arange(len(used)), meg_picks)].any()
+    # This gets written to the Info struct
     sss_cal = dict(cal_corrs=np.array(cal_corrs),
                    cal_chans=np.array(cal_chans))
 
+    # Log quantification of sensor changes
     # Deal with numerical precision giving absolute vals slightly more than 1.
     np.clip(ang_shift, -1., 1., ang_shift)
     np.rad2deg(np.arccos(ang_shift), ang_shift)  # Convert to degrees
-
-    # Log quantification of sensor changes
     logger.info('        Adjusted coil positions by (μ ± σ): '
                 '%0.1f° ± %0.1f° (max: %0.1f°)' %
                 (np.mean(ang_shift), np.std(ang_shift),
                  np.max(np.abs(ang_shift))))
-
-    # Determine gradiometer imbalances and magnetometer calibrations
-    grad_picks = pick_types(info, meg='grad', exclude=[])
-    mag_picks = pick_types(info, meg='mag', exclude=[])
-    grad_imbalances = np.array([cal_chs[ii]['calib_coeff']
-                                for ii in grad_picks]).T
-    if grad_imbalances.shape[0] not in [1, 3]:
-        raise ValueError('Must have 1 (x) or 3 (x, y, z) point-like ' +
-                         'magnetometers. Currently have %i' %
-                         grad_imbalances.shape[0])
-    mag_cals = np.array([cal_chs[ii]['calib_coeff'] for ii in mag_picks])
-
-    # Now let's actually construct our point-like adjustment coils for grads
-    grad_coilsets = _get_grad_point_coilsets(
-        info, n_types=len(grad_imbalances), ignore_ref=ignore_ref)
-    calibration = dict(grad_imbalances=grad_imbalances,
-                       grad_coilsets=grad_coilsets, mag_cals=mag_cals)
     return calibration, sss_cal
 
 
@@ -1845,10 +1868,10 @@ def _regularize_in(int_order, ext_order, S_decomp, mag_or_fine):
                      'Removing in component %s: l=%s, m=%+0.0f'
                      % (tuple(eigs[ii]) + (eigs[ii, 0] / eigs[ii, 1],
                         ri, degrees[ri], orders[ri])))
-    logger.info('        Resulting information: %0.1f bits/sample '
-                '(%0.1f%% of peak %0.1f)'
-                % (I_tots[lim_idx], 100 * I_tots[lim_idx] / max_info,
-                   max_info))
+    logger.debug('        Resulting information: %0.1f bits/sample '
+                 '(%0.1f%% of peak %0.1f)'
+                 % (I_tots[lim_idx], 100 * I_tots[lim_idx] / max_info,
+                    max_info))
     return in_removes, out_removes
 
 
@@ -1901,6 +1924,7 @@ def _trans_sss_basis(exp, all_coils, trans=None, coil_scale=100.):
     if trans is not None:
         if not isinstance(trans, Transform):
             trans = Transform('meg', 'head', trans)
+        assert not np.isnan(trans['trans']).any()
         all_coils = (apply_trans(trans, all_coils[0]),
                      apply_trans(trans, all_coils[1], move=False),
                      ) + all_coils[2:]

mne/preprocessing/tests/test_fine_cal.py

@@ -0,0 +1,40 @@
+# Author: Mark Wronkiewicz <wronk@uw.edu>
+#
+# License: BSD (3-clause)
+
+import os.path as op
+import warnings
+
+from mne.datasets import testing
+from mne.preprocessing._fine_cal import (read_fine_calibration,
+                                         write_fine_calibration)
+from mne.utils import _TempDir, object_hash, run_tests_if_main
+from nose.tools import assert_equal
+
+warnings.simplefilter('always')  # Always throw warnings
+
+# Define fine calibration filepaths
+data_path = testing.data_path(download=False)
+fine_cal_fname = op.join(data_path, 'SSS', 'sss_cal_3053.dat')
+fine_cal_fname_3d = op.join(data_path, 'SSS', 'sss_cal_3053_3d.dat')
+
+
+@testing.requires_testing_data
+def test_read_write_fine_cal():
+    """Test round trip reading/writing of fine calibration .dat file"""
+    temp_dir = _TempDir()
+    temp_fname = op.join(temp_dir, 'fine_cal_temp.dat')
+
+    for fname in [fine_cal_fname, fine_cal_fname_3d]:
+        # Load fine calibration file
+        fine_cal_dict = read_fine_calibration(fname)
+
+        # Save temp version of fine calibration file
+        write_fine_calibration(temp_fname, fine_cal_dict)
+        fine_cal_dict_reload = read_fine_calibration(temp_fname)
+
+        # Load temp version of fine calibration file and compare hashes
+        assert_equal(object_hash(fine_cal_dict),
+                     object_hash(fine_cal_dict_reload))
+
+run_tests_if_main()

mne/preprocessing/tests/test_maxwell.py

@@ -76,6 +76,20 @@ ctc_mgh_fname = op.join(sss_path, 'ct_sparse_mgh.fif')
 sample_fname = op.join(data_path, 'MEG', 'sample',
                        'sample_audvis_trunc_raw.fif')
 
+triux_path = op.join(data_path, 'SSS', 'TRIUX')
+tri_fname = op.join(triux_path, 'triux_bmlhus_erm_raw.fif')
+tri_sss_fname = op.join(triux_path, 'triux_bmlhus_erm_raw_sss.fif')
+tri_sss_reg_fname = op.join(triux_path, 'triux_bmlhus_erm_regIn_raw_sss.fif')
+tri_sss_st4_fname = op.join(triux_path, 'triux_bmlhus_erm_st4_raw_sss.fif')
+tri_sss_ctc_fname = op.join(triux_path, 'triux_bmlhus_erm_ctc_raw_sss.fif')
+tri_sss_cal_fname = op.join(triux_path, 'triux_bmlhus_erm_cal_raw_sss.fif')
+tri_sss_ctc_cal_fname = op.join(
+    triux_path, 'triux_bmlhus_erm_ctc_cal_raw_sss.fif')
+tri_sss_ctc_cal_reg_in_fname = op.join(
+    triux_path, 'triux_bmlhus_erm_ctc_cal_regIn_raw_sss.fif')
+tri_ctc_fname = op.join(triux_path, 'ct_sparse_BMLHUS.fif')
+tri_cal_fname = op.join(triux_path, 'sss_cal_BMLHUS.dat')
+
 io_dir = op.join(op.dirname(__file__), '..', '..', 'io')
 fname_ctf_raw = op.join(io_dir, 'tests', 'data', 'test_ctf_comp_raw.fif')
 
@@ -629,19 +643,24 @@ def test_regularization():
         assert_meg_snr(raw_sss, sss_reg_in, min_tols[ii], med_tols[ii], msg=rf)
 
         # check components match
-        py_info = raw_sss.info['proc_history'][0]['max_info']['sss_info']
-        assert_true(py_info is not None)
-        assert_true(len(py_info) > 0)
-        mf_info = sss_reg_in.info['proc_history'][0]['max_info']['sss_info']
-        n_in = None
-        for inf in py_info, mf_info:
-            if n_in is None:
-                n_in = _get_n_moments(inf['in_order'])
-            else:
-                assert_equal(n_in, _get_n_moments(inf['in_order']))
-            assert_equal(inf['components'][:n_in].sum(), inf['nfree'])
-        assert_allclose(py_info['nfree'], mf_info['nfree'],
-                        atol=comp_tols[ii], err_msg=rf)
+        _check_reg_match(raw_sss, sss_reg_in, comp_tols[ii])
+
+
+def _check_reg_match(sss_py, sss_mf, comp_tol):
+    """Helper to check regularization"""
+    info_py = sss_py.info['proc_history'][0]['max_info']['sss_info']
+    assert_true(info_py is not None)
+    assert_true(len(info_py) > 0)
+    info_mf = sss_mf.info['proc_history'][0]['max_info']['sss_info']
+    n_in = None
+    for inf in (info_py, info_mf):
+        if n_in is None:
+            n_in = _get_n_moments(inf['in_order'])
+        else:
+            assert_equal(n_in, _get_n_moments(inf['in_order']))
+        assert_equal(inf['components'][:n_in].sum(), inf['nfree'])
+    assert_allclose(info_py['nfree'], info_mf['nfree'],
+                    atol=comp_tol, err_msg=sss_py._filenames[0])
 
 
 @testing.requires_testing_data
@@ -840,7 +859,7 @@ def test_all():
     coord_frames = ('head', 'head', 'meg', 'head')
     ctcs = (ctc_fname, ctc_fname, ctc_fname, ctc_mgh_fname)
     mins = (3.5, 3.5, 1.2, 0.9)
-    meds = (10.9, 10.4, 3.2, 6.)
+    meds = (10.8, 10.4, 3.2, 6.)
     st_durs = (1., 1., 1., None)
     destinations = (None, sample_fname, None, None)
     origins = (mf_head_origin,
@@ -857,4 +876,46 @@ def test_all():
         sss_mf = Raw(sss_fnames[ii])
         assert_meg_snr(sss_py, sss_mf, mins[ii], meds[ii], msg=rf)
 
+
+@slow_test
+@requires_svd_convergence
+@testing.requires_testing_data
+def test_triux():
+    """Test TRIUX system support"""
+    raw = Raw(tri_fname).crop(0, 0.999)
+    raw.fix_mag_coil_types()
+    # standard
+    sss_py = maxwell_filter(raw, coord_frame='meg', regularize=None)
+    assert_meg_snr(sss_py, Raw(tri_sss_fname), 37, 700)
+    # cross-talk
+    sss_py = maxwell_filter(raw, coord_frame='meg', regularize=None,
+                            cross_talk=tri_ctc_fname)
+    assert_meg_snr(sss_py, Raw(tri_sss_ctc_fname), 35, 700)
+    # fine cal
+    sss_py = maxwell_filter(raw, coord_frame='meg', regularize=None,
+                            calibration=tri_cal_fname)
+    assert_meg_snr(sss_py, Raw(tri_sss_cal_fname), 31, 360)
+    # ctc+cal
+    sss_py = maxwell_filter(raw, coord_frame='meg', regularize=None,
+                            calibration=tri_cal_fname,
+                            cross_talk=tri_ctc_fname)
+    assert_meg_snr(sss_py, Raw(tri_sss_ctc_cal_fname), 31, 350)
+    # regularization
+    sss_py = maxwell_filter(raw, coord_frame='meg', regularize='in')
+    sss_mf = Raw(tri_sss_reg_fname)
+    assert_meg_snr(sss_py, sss_mf, 0.6, 9)
+    _check_reg_match(sss_py, sss_mf, 1)
+    # all three
+    sss_py = maxwell_filter(raw, coord_frame='meg', regularize='in',
+                            calibration=tri_cal_fname,
+                            cross_talk=tri_ctc_fname)
+    sss_mf = Raw(tri_sss_ctc_cal_reg_in_fname)
+    assert_meg_snr(sss_py, sss_mf, 0.6, 9)
+    _check_reg_match(sss_py, sss_mf, 1)
+    # tSSS
+    raw = Raw(tri_fname).fix_mag_coil_types()
+    sss_py = maxwell_filter(raw, coord_frame='meg', regularize=None,
+                            st_duration=4., verbose=True)
+    assert_meg_snr(sss_py, Raw(tri_sss_st4_fname), 700., 1600)
+
 run_tests_if_main()

mne/tests/test_source_estimate.py

@@ -149,7 +149,7 @@ def test_volume_stc():
         assert_array_almost_equal(img.get_affine(), t1_img.get_affine(),
                                   decimal=5)
 
-    except ImportError:
+    except Exception:
         print('Save as nifti test skipped, needs NiBabel')
 
 

mne/tests/test_surface.py

@@ -13,7 +13,8 @@ from mne import read_surface, write_surface, decimate_surface
 from mne.surface import (read_morph_map, _compute_nearest,
                          fast_cross_3d, get_head_surf, read_curvature,
                          get_meg_helmet_surf)
-from mne.utils import _TempDir, requires_mayavi, run_tests_if_main, slow_test
+from mne.utils import (_TempDir, requires_mayavi, run_tests_if_main, slow_test,
+                       requires_tvtk)
 from mne.io import read_info
 from mne.transforms import _get_trans
 
@@ -150,6 +151,7 @@ def test_read_curv():
     assert_true(np.logical_or(bin_curv == 0, bin_curv == 1).all())
 
 
+@requires_tvtk
 @requires_mayavi
 def test_decimate_surface():
     """Test triangular surface decimation

mne/time_frequency/tfr.py

@@ -774,10 +774,8 @@ class AverageTFR(ContainsMixin, UpdateChannelsMixin):
         if 'grad' in self:
             types.append('grad')
         fig = figure_nobar()
-        fig.suptitle('{:.2f} s - {:.2f} s, {:.2f} Hz - {:.2f} Hz'.format(tmin,
-                                                                         tmax,
-                                                                         fmin,
-                                                                         fmax),
+        fig.suptitle('{0:.2f} s - {1:.2f} s, '
+                     '{2:.2f} Hz - {3:.2f} Hz'.format(tmin, tmax, fmin, fmax),
                      y=0.04)
         for idx, ch_type in enumerate(types):
             ax = plt.subplot(1, len(types), idx + 1)

mne/utils.py

@@ -682,7 +682,7 @@ def has_nibabel(vox2ras_tkr=False):
                                    'header_class', 0),
                            'get_vox2ras_tkr', None) is not None)
         return out
-    except ImportError:
+    except Exception:
         return False
 
 

setup.cfg

@@ -29,4 +29,4 @@ doctest-fixtures = _fixture
 
 [flake8]
 exclude = __init__.py,*externals*,constants.py,fixes.py
-ignore = E241
+ignore = E241,E305

tutorials/plot_dipole_fit.py

@@ -75,7 +75,7 @@ diff = combine_evoked([evoked, pred_evoked], [1, -1])
 plot_params['colorbar'] = True
 diff.plot_topomap(time_format='Difference', axes=axes[2], **plot_params)
 plt.suptitle('Comparison of measured and predicted fields '
-             'at {:.0f} ms'.format(best_time * 1000.), fontsize=16)
+             'at {0:.0f} ms'.format(best_time * 1000.), fontsize=16)
 
 ###############################################################################
 # Estimate the time course of a single dipole with fixed position and

tutorials/plot_python_intro.py

@@ -34,7 +34,7 @@ print(type(c))
 # If you come from a background of matlab, remember that indexing in python
 # starts from zero:
 a = [1, 2, 3, 4]
-print('This is the zeroth value in the list: {}'.format(a[0]))
+print('This is the zeroth value in the list: {0}'.format(a[0]))
 
 ###############################################################################
 # No need to reinvent the wheel. Scipy and Numpy are battle field tested

tutorials/plot_stats_spatio_temporal_cluster_sensors.py

@@ -70,8 +70,8 @@ print(type(connectivity))  # it's a sparse matrix!
 
 plt.imshow(connectivity.toarray(), cmap='gray', origin='lower',
            interpolation='nearest')
-plt.xlabel('{} Magnetometers'.format(len(ch_names)))
-plt.ylabel('{} Magnetometers'.format(len(ch_names)))
+plt.xlabel('{0} Magnetometers'.format(len(ch_names)))
+plt.ylabel('{0} Magnetometers'.format(len(ch_names)))
 plt.title('Between-sensor adjacency')
 
 ###############################################################################