This is an automated email from the git hooks/post-receive script. yoh pushed a commit to annotated tag v0.1 in repository python-mne.
commit bd1860eb7254ceacbc5fa32c5775b8a3d51c17a1 Author: Alexandre Gramfort <[email protected]> Date: Thu Apr 14 18:28:23 2011 -0400 pep8 --- mne/bem_surfaces.py | 69 +++++++++++++++++----------------- mne/cov.py | 30 ++++++++------- mne/epochs.py | 20 +++++----- mne/event.py | 10 ++--- mne/fiff/compensator.py | 50 ++++++++++++------------- mne/fiff/meas_info.py | 2 +- mne/fiff/tag.py | 99 ++++++++++++++++++++++++++----------------------- mne/fiff/write.py | 13 +++---- mne/forward.py | 99 +++++++++++++++++++++++++------------------------ mne/inverse.py | 70 +++++++++++++++++----------------- mne/preprocessing.py | 19 +++++----- mne/transforms.py | 3 +- 12 files changed, 246 insertions(+), 238 deletions(-) diff --git a/mne/bem_surfaces.py b/mne/bem_surfaces.py index 6cec01e..c9fc3db 100755 --- a/mne/bem_surfaces.py +++ b/mne/bem_surfaces.py @@ -14,18 +14,17 @@ from .fiff.tag import find_tag # # These fiff definitions are not needed elsewhere # -FIFFB_BEM = 310 # BEM data -FIFFB_BEM_SURF = 311 # One of the surfaces -# -FIFF_BEM_SURF_ID = 3101 # int surface number -FIFF_BEM_SURF_NAME = 3102 # string surface name -FIFF_BEM_SURF_NNODE = 3103 # int # of nodes on a surface -FIFF_BEM_SURF_NTRI = 3104 # int # number of triangles on a surface -FIFF_BEM_SURF_NODES = 3105 # float surface nodes (nnode,3) -FIFF_BEM_SURF_TRIANGLES = 3106 # int surface triangles (ntri,3) -FIFF_BEM_SURF_NORMALS = 3107 # float surface node normal unit vectors (nnode,3) -FIFF_BEM_COORD_FRAME = 3112 # The coordinate frame of the mode -FIFF_BEM_SIGMA = 3113 # Conductivity of a compartment +FIFFB_BEM = 310 # BEM data +FIFFB_BEM_SURF = 311 # One of the surfaces +FIFF_BEM_SURF_ID = 3101 # int surface number +FIFF_BEM_SURF_NAME = 3102 # string surface name +FIFF_BEM_SURF_NNODE = 3103 # int number of nodes on a surface +FIFF_BEM_SURF_NTRI = 3104 # int number of triangles on a surface +FIFF_BEM_SURF_NODES = 3105 # float surface nodes (nnode,3) +FIFF_BEM_SURF_TRIANGLES = 3106 # int surface triangles (ntri,3) +FIFF_BEM_SURF_NORMALS = 3107 # float surface node normal unit vectors +FIFF_BEM_COORD_FRAME = 3112 # The coordinate frame of the mode +FIFF_BEM_SIGMA = 3113 # Conductivity of a compartment def read_bem_surfaces(fname, add_geom=False): @@ -58,7 +57,7 @@ def read_bem_surfaces(fname, add_geom=False): bem = dir_tree_find(tree, FIFFB_BEM) if bem is None: fid.close() - raise ValueError, 'BEM data not found' + raise ValueError('BEM data not found') bem = bem[0] # @@ -67,7 +66,7 @@ def read_bem_surfaces(fname, add_geom=False): bemsurf = dir_tree_find(bem, FIFFB_BEM_SURF) if bemsurf is None: fid.close() - raise ValueError, 'BEM surface data not found' + raise ValueError('BEM surface data not found') print '\t%d BEM surfaces found' % len(bemsurf) # @@ -117,15 +116,15 @@ def _read_bem_surface(fid, this, def_coord_frame): tag = find_tag(fid, this, FIFF_BEM_SURF_NNODE) if tag is None: fid.close() - raise ValueError, 'Number of vertices not found' + raise ValueError('Number of vertices not found') res = dict() res['np'] = tag.data - tag = find_tag(fid, this,FIFF_BEM_SURF_NTRI) + tag = find_tag(fid, this, FIFF_BEM_SURF_NTRI) if tag is None: fid.close() - raise ValueError, 'Number of triangles not found' + raise ValueError('Number of triangles not found') else: res['ntri'] = tag.data @@ -144,12 +143,12 @@ def _read_bem_surface(fid, this, def_coord_frame): tag = find_tag(fid, this, FIFF_BEM_SURF_NODES) if tag is None: fid.close() - raise ValueError, 'Vertex data not found' + raise ValueError('Vertex data not found') - res['rr'] = tag.data.astype(np.float) # XXX : double because of mayavi bug + res['rr'] = tag.data.astype(np.float) # XXX : double because of mayavi bug if res['rr'].shape[0] != res['np']: fid.close() - raise ValueError, 'Vertex information is incorrect' + raise ValueError('Vertex information is incorrect') tag = find_tag(fid, this, FIFF.FIFF_MNE_SOURCE_SPACE_NORMALS) if tag is None: @@ -158,17 +157,17 @@ def _read_bem_surface(fid, this, def_coord_frame): res['nn'] = tag.data if res['nn'].shape[0] != res['np']: fid.close() - raise ValueError, 'Vertex normal information is incorrect' + raise ValueError('Vertex normal information is incorrect') tag = find_tag(fid, this, FIFF_BEM_SURF_TRIANGLES) if tag is None: fid.close() - raise ValueError, 'Triangulation not found' + raise ValueError('Triangulation not found') - res['tris'] = tag.data - 1 # index start at 0 in Python + res['tris'] = tag.data - 1 # index start at 0 in Python if res['tris'].shape[0] != res['ntri']: fid.close() - raise ValueError, 'Triangulation information is incorrect' + raise ValueError('Triangulation information is incorrect') return res @@ -182,35 +181,35 @@ def _complete_surface_info(this): print '\tCompleting triangulation info...' print 'triangle normals...' this['tri_area'] = np.zeros(this['ntri']) - r1 = this['rr'][this['tris'][:,0],:] - r2 = this['rr'][this['tris'][:,1],:] - r3 = this['rr'][this['tris'][:,2],:] - this['tri_cent'] = (r1+r2+r3) /3.0 - this['tri_nn'] = np.cross((r2-r1), (r3-r1)) + r1 = this['rr'][this['tris'][:, 0], :] + r2 = this['rr'][this['tris'][:, 1], :] + r3 = this['rr'][this['tris'][:, 2], :] + this['tri_cent'] = (r1 + r2 + r3) / 3.0 + this['tri_nn'] = np.cross((r2 - r1), (r3 - r1)) # # Triangle normals and areas # for p in range(this['ntri']): - size = linalg.norm(this['tri_nn'][p,:]) + size = linalg.norm(this['tri_nn'][p, :]) this['tri_area'][p] = size / 2.0 if size > 0.0: - this['tri_nn'][p,:] = this['tri_nn'][p,:] / size + this['tri_nn'][p, :] = this['tri_nn'][p, :] / size # # Accumulate the vertex normals # print 'vertex normals...' this['nn'] = np.zeros((this['np'], 3)) for p in range(this['ntri']): - this['nn'][this['tris'][p,:],:] = this['nn'][this['tris'][p,:],:] \ - + np.kron(np.ones((3, 1)), this['tri_nn'][p,:]) + this['nn'][this['tris'][p, :], :] = this['nn'][this['tris'][p, :], :] \ + + np.kron(np.ones((3, 1)), this['tri_nn'][p, :]) # # Compute the lengths of the vertex normals and scale # print 'normalize...' for p in range(this['np']): - size = linalg.norm(this['nn'][p,:]) + size = linalg.norm(this['nn'][p, :]) if size > 0: - this['nn'][p,:] = this['nn'][p,:] / size + this['nn'][p, :] = this['nn'][p, :] / size print '[done]' return this diff --git a/mne/cov.py b/mne/cov.py index a3b7da4..3843075 100755 --- a/mne/cov.py +++ b/mne/cov.py @@ -25,7 +25,7 @@ from .epochs import Epochs, _is_good def rank(A, tol=1e-8): s = linalg.svd(A, compute_uv=0) - return np.sum(np.where(s > s[0]*tol, 1, 0)) + return np.sum(np.where(s > s[0] * tol, 1, 0)) def _get_whitener(A, rnk, pca, ch_type): @@ -35,10 +35,10 @@ def _get_whitener(A, rnk, pca, ch_type): D = D[I] V = V[:, I] D = 1.0 / D - if not pca: # No PCA case. + if not pca: # No PCA case. print 'Not doing PCA for %s.' % ch_type W = np.sqrt(D)[:, None] * V.T - else: # Rey's approach. MNE has been changed to implement this. + else: # Rey's approach. MNE has been changed to implement this. print 'Setting small %s eigenvalues to zero.' % ch_type D[rnk:] = 0.0 W = np.sqrt(D)[:, None] * V.T @@ -51,8 +51,8 @@ def _get_whitener(A, rnk, pca, ch_type): class Covariance(object): """Noise covariance matrix""" - _kind_to_id = dict(full=1, sparse=2, diagonal=3) # XXX : check - _id_to_kind = {1: 'full', 2: 'sparse', 3: 'diagonal'} # XXX : check + _kind_to_id = dict(full=1, sparse=2, diagonal=3) # XXX : check + _id_to_kind = {1: 'full', 2: 'sparse', 3: 'diagonal'} # XXX : check def __init__(self, fname, kind='full'): self.kind = kind @@ -123,7 +123,7 @@ class Covariance(object): C_idx = [k for k, name in enumerate(self.ch_names) if name in info['ch_names'] and name not in bads] ch_names = [self.ch_names[k] for k in C_idx] - C_noise = self.data[np.ix_(C_idx, C_idx)] # take covariance submatrix + C_noise = self.data[np.ix_(C_idx, C_idx)] # take covariance submatrix # Create the projection operator proj, ncomp, _ = make_projector(info['projs'], ch_names) @@ -160,10 +160,10 @@ class Covariance(object): # estimate noise covariance matrix rank # Loop on all the required data types (MEG MAG, MEG GRAD, EEG) - if has_meg: # Separate rank of MEG + if has_meg: # Separate rank of MEG rank_meg = rank(C_noise[C_ind_meg][:, C_ind_meg]) print 'Rank of MEG part of noise covariance is %d' % rank_meg - if has_eeg: # Separate rank of EEG + if has_eeg: # Separate rank of EEG rank_eeg = rank(C_noise[C_ind_eeg][:, C_ind_eeg]) print 'Rank of EEG part of noise covariance is %d' % rank_eeg @@ -173,7 +173,7 @@ class Covariance(object): # add constant on diagonal C_noise[ind, ind] += reg * np.mean(variances[ind]) - if has_meg and has_eeg: # Sets cross terms to zero + if has_meg and has_eeg: # Sets cross terms to zero C_noise[np.ix_(C_ind_meg, C_ind_eeg)] = 0.0 C_noise[np.ix_(C_ind_eeg, C_ind_meg)] = 0.0 @@ -186,11 +186,11 @@ class Covariance(object): W_eeg = _get_whitener(C_noise[C_ind_eeg][:, C_ind_eeg], rank_eeg, pca, 'EEG') - if has_meg and not has_eeg: # Only MEG case. + if has_meg and not has_eeg: # Only MEG case. W = W_meg - elif has_eeg and not has_meg: # Only EEG case. + elif has_eeg and not has_meg: # Only EEG case. W = W_eeg - elif has_eeg and has_meg: # Bimodal MEG and EEG case. + elif has_eeg and has_meg: # Bimodal MEG and EEG case. # Whitening of MEG and EEG separately, which assumes zero # covariance between MEG and EEG (i.e., a block diagonal noise # covariance). This was recommended by Matti as EEG does not @@ -227,6 +227,7 @@ class Whitener(object): ############################################################################### # IO + def read_cov(fid, node, cov_kind): """Read a noise covariance matrix @@ -299,7 +300,7 @@ def read_cov(fid, node, cov_kind): data = np.zeros((dim, dim)) data[np.tril(np.ones((dim, dim))) > 0] = vals data = data + data.T - data.flat[::dim+1] /= 2.0 + data.flat[::dim + 1] /= 2.0 diagmat = False print '\t%d x %d full covariance (kind = %d) found.' \ % (dim, dim, cov_kind) @@ -338,6 +339,7 @@ def read_cov(fid, node, cov_kind): ############################################################################### # Estimate from data + def _estimate_compute_covariance_from_epochs(epochs, bmin, bmax, reject, flat, keep_sample_mean): """Estimate noise covariance matrix from epochs @@ -453,7 +455,7 @@ def compute_raw_data_covariance(raw, tmin=None, tmax=None, tstep=0.2, print "Artefact detected in [%d, %d]" % (first, last) mu /= n_samples - data -= n_samples * mu[:,None] * mu[None,:] + data -= n_samples * mu[:, None] * mu[None, :] data /= (n_samples - 1.0) print "Number of samples used : %d" % n_samples print '[done]' diff --git a/mne/epochs.py b/mne/epochs.py index 1b0d629..e29af57 100755 --- a/mne/epochs.py +++ b/mne/epochs.py @@ -110,7 +110,7 @@ class Epochs(object): self.ch_names = [raw.info['ch_names'][k] for k in picks] if len(picks) == 0: - raise ValueError, "Picks cannot be empty." + raise ValueError("Picks cannot be empty.") # Set up projection if raw.info['projs'] is None: @@ -155,12 +155,13 @@ class Epochs(object): if n_events > 0: print '%d matching events found' % n_events else: - raise ValueError, 'No desired events found.' + raise ValueError('No desired events found.') # Handle times sfreq = raw.info['sfreq'] - self.times = np.arange(int(round(tmin*sfreq)), int(round(tmax*sfreq)), - dtype=np.float) / sfreq + self.times = np.arange(int(round(tmin * sfreq)), + int(round(tmax * sfreq)), + dtype=np.float) / sfreq # setup epoch rejection self._reject_setup() @@ -184,8 +185,7 @@ class Epochs(object): self.ch_names = self.info['ch_names'] if self.preload: - self._data = self._data[:,idx,:] - + self._data = self._data[:, idx, :] def get_epoch(self, idx): """Load one epoch @@ -207,7 +207,7 @@ class Epochs(object): # Read a data segment first_samp = self.raw.first_samp - start = int(round(event_samp + self.tmin*sfreq)) - first_samp + start = int(round(event_samp + self.tmin * sfreq)) - first_samp stop = start + len(self.times) epoch, _ = self.raw[self.picks, start:stop] @@ -330,7 +330,7 @@ class Epochs(object): evoked.data = data evoked.times = self.times.copy() evoked.comment = self.name - evoked.aspect_kind = np.array([100]) # for standard average file + evoked.aspect_kind = np.array([100]) # for standard average file evoked.nave = n_events evoked.first = - int(np.sum(self.times < 0)) evoked.last = int(np.sum(self.times > 0)) @@ -379,8 +379,8 @@ def _is_good(e, ch_names, channel_type_idx, reject, flat): delta = deltas[idx_min_delta] if delta < thresh: ch_name = ch_names[idx[idx_min_delta]] - print '\tRejecting flat epoch based on %s : %s (%s < %s).' \ - % (name, ch_name, delta, thresh) + print ('\tRejecting flat epoch based on %s : %s (%s < %s).' + % (name, ch_name, delta, thresh)) return False return True diff --git a/mne/event.py b/mne/event.py index 68cbb1c..d5b8aef 100755 --- a/mne/event.py +++ b/mne/event.py @@ -37,7 +37,7 @@ def read_events(filename): if len(events) == 0: fid.close() - raise ValueError, 'Could not find event data' + raise ValueError('Could not find event data') events = events[0] @@ -51,7 +51,7 @@ def read_events(filename): break else: fid.close() - raise ValueError, 'Could not find any events' + raise ValueError('Could not find any events') event_list = event_list.reshape(len(event_list) / 3, 3) return event_list @@ -99,11 +99,11 @@ def find_events(raw, stim_channel='STI 014'): exclude=[]) if len(pick) == 0: raise ValueError('No stim channel found to extract event triggers.') - data, times = raw[pick,:] + data, times = raw[pick, :] data = data.ravel() idx = np.where(np.diff(data.ravel()) > 0)[0] n_events = len(idx) - events_id = data[idx+1].astype(np.int) + events_id = data[idx + 1].astype(np.int) idx += raw.first_samp + 1 events = np.c_[idx, np.zeros_like(idx), events_id] - return events \ No newline at end of file + return events diff --git a/mne/fiff/compensator.py b/mne/fiff/compensator.py index b64dfbc..7c79771 100755 --- a/mne/fiff/compensator.py +++ b/mne/fiff/compensator.py @@ -14,8 +14,8 @@ def get_current_comp(info): if first_comp < 0: first_comp = comp elif comp != first_comp: - raise ValueError, ('Compensation is not set equally on ' - 'all MEG channels') + raise ValueError('Compensation is not set equally on ' + 'all MEG channels') return comp @@ -41,7 +41,7 @@ def get_current_comp(info): # for k in range(len(info['comps'])): # if info['comps'][k]['kind'] == kind: # this_data = info['comps'][k]['data'] -# +# # # Create the preselector # presel = np.zeros((this_data['ncol'], info['nchan'])) # for col, col_name in enumerate(this_data['col_names']): @@ -52,7 +52,7 @@ def get_current_comp(info): # elif len(ind) > 1: # raise ValueError, 'Ambiguous channel %s' % col_name # presel[col, ind] = 1.0 -# +# # # Create the postselector # postsel = np.zeros((info['nchan'], this_data['nrow'])) # for c, ch_name in enumerate(info['ch_names']): @@ -61,10 +61,10 @@ def get_current_comp(info): # raise ValueError, 'Ambiguous channel %s' % ch_name # elif len(ind) == 1: # postsel[c, ind] = 1.0 -# +# # this_comp = np.dot(postsel, np.dot(this_data['data'], presel)) # return this_comp -# +# # return [] @@ -78,17 +78,17 @@ def get_current_comp(info): # % to - desired compensation in the output # % exclude_comp_chs - exclude compensation channels from the output (optional) # % -# +# # % # % Create a compensation matrix to bring the data from one compensation # % state to another # % # """ -# +# # if from_ == to: # comp = np.zeros((info['nchan'], info['nchan'])) # return comp -# +# # if from_ == 0: # C1 = np.zeros((info['nchan'], info['nchan'])) # else: @@ -96,12 +96,12 @@ def get_current_comp(info): # C1 = _make_compensator(info, from_) # except Exception as inst: # raise ValueError, 'Cannot create compensator C1 (%s)' % inst -# +# # if len(C1) == 0: # raise ValueError, ('Desired compensation matrix (kind = %d) not' # ' found' % from_) -# -# +# +# # if to == 0: # C2 = np.zeros((info['nchan'], info['nchan'])) # else: @@ -109,17 +109,17 @@ def get_current_comp(info): # C2 = _make_compensator(info, to) # except Exception as inst: # raise ValueError, 'Cannot create compensator C2 (%s)' % inst -# +# # if len(C2) == 0: # raise ValueError, ('Desired compensation matrix (kind = %d) not ' # 'found' % to) -# -# +# +# # # s_orig = s_from + C1*s_from = (I + C1)*s_from # # s_to = s_orig - C2*s_orig = (I - C2)*s_orig # # s_to = (I - C2)*(I + C1)*s_from = (I + C1 - C2 - C2*C1)*s_from # comp = np.eye(info['nchan']) + C1 - C2 - C2*C1 -# +# # if exclude_comp_chs: # pick = np.zeros((info['nchan'], info['nchan'])) # npick = 0 @@ -127,13 +127,13 @@ def get_current_comp(info): # if chan['kind'] != FIFF.FIFFV_REF_MEG_CH: # npick += 1 # pick[npick] = k -# +# # if npick == 0: # raise ValueError, ('Nothing remains after excluding the ' # 'compensation channels') -# +# # comp = comp[pick[1:npick], :] -# +# # return comp @@ -145,20 +145,20 @@ def get_current_comp(info): # % Apply compensation to the data as desired # % # """ -# +# # newdata = data.copy() # now = get_current_comp(newdata['info']) -# +# # # Are we there already? # if now == to: # print 'Data are already compensated as desired' -# +# # # Make the compensator and apply it to all data sets # comp = make_compensator(newdata['info'], now, to) # for k in range(len(newdata['evoked'])): # newdata['evoked'][k]['epochs'] = np.dot(comp, # newdata['evoked'][k]['epochs']) -# +# # # Update the compensation info in the channel descriptors # newdata['info']['chs'] = set_current_comp(newdata['info']['chs'], to) # return newdata @@ -168,11 +168,11 @@ def get_current_comp(info): # """Set the current compensation value in the channel info structures # """ # new_chs = chs -# +# # lower_half = int('FFFF', 16) # hex2dec('FFFF') # for k in range(len(chs)): # if chs[k].kind == FIFF.FIFFV_MEG_CH: # coil_type = float(chs[k]['coil_type']) & lower_half # new_chs[k]['coil_type'] = int(coil_type | (value << 16)) -# +# # return new_chs diff --git a/mne/fiff/meas_info.py b/mne/fiff/meas_info.py index 52898a6..ddaedd6 100755 --- a/mne/fiff/meas_info.py +++ b/mne/fiff/meas_info.py @@ -201,7 +201,7 @@ def read_meas_info(fid, tree): info['nchan'] = nchan info['sfreq'] = sfreq info['highpass'] = highpass if highpass is not None else 0 - info['lowpass'] = lowpass if lowpass is not None else info['sfreq']/2.0 + info['lowpass'] = lowpass if lowpass is not None else info['sfreq'] / 2.0 # Add the channel information and make a list of channel names # for convenience diff --git a/mne/fiff/tag.py b/mne/fiff/tag.py index d3652ae..e05d1ef 100755 --- a/mne/fiff/tag.py +++ b/mne/fiff/tag.py @@ -65,7 +65,7 @@ class Tag(object): def read_tag_info(fid): """Read Tag info (or header) """ - s = fid.read(4*4) + s = fid.read(4 * 4) tag = Tag(*struct.unpack(">iiii", s)) if tag.next == 0: fid.seek(tag.size, 1) @@ -93,13 +93,13 @@ def read_tag(fid, pos=None): if pos is not None: fid.seek(pos, 0) - s = fid.read(4*4) + s = fid.read(4 * 4) tag = Tag(*struct.unpack(">iIii", s)) # # The magic hexadecimal values # - is_matrix = 4294901760 # ffff0000 + is_matrix = 4294901760 # ffff0000 matrix_coding_dense = 16384 # 4000 matrix_coding_CCS = 16400 # 4010 matrix_coding_RCS = 16416 # 4020 @@ -124,8 +124,8 @@ def read_tag(fid, pos=None): fid.seek(pos, 0) if ndim != 2: - raise ValueError, 'Only two-dimensional matrices are ' \ - 'supported at this time' + raise ValueError('Only two-dimensional matrices are ' + 'supported at this time') matrix_type = data_type & tag.type @@ -142,15 +142,15 @@ def read_tag(fid, pos=None): tag.data = np.fromfile(fid, dtype='>f8', count=dims.prod()).reshape(dims) elif matrix_type == FIFF.FIFFT_COMPLEX_FLOAT: - data = np.fromfile(fid, dtype='>f4', count=2*dims.prod()) + data = np.fromfile(fid, dtype='>f4', count=2 * dims.prod()) # Note: we need the non-conjugate transpose here tag.data = (data[::2] + 1j * data[1::2]).reshape(dims) elif matrix_type == FIFF.FIFFT_COMPLEX_DOUBLE: - data = np.fromfile(fid, dtype='>f8', count=2*dims.prod()) + data = np.fromfile(fid, dtype='>f8', count=2 * dims.prod()) # Note: we need the non-conjugate transpose here tag.data = (data[::2] + 1j * data[1::2]).reshape(dims) else: - raise ValueError, 'Cannot handle matrix of type %d yet' % ( + raise ValueError('Cannot handle matrix of type %d yet' % matrix_type) elif matrix_coding == matrix_coding_CCS or \ @@ -158,13 +158,13 @@ def read_tag(fid, pos=None): from scipy import sparse # Find dimensions and return to the beginning of tag data pos = fid.tell() - fid.seek(tag.size-4, 1) + fid.seek(tag.size - 4, 1) ndim = int(np.fromfile(fid, dtype='>i', count=1)) - fid.seek(-(ndim+2)*4, 1) - dims = np.fromfile(fid, dtype='>i', count=ndim+1) + fid.seek(-(ndim + 2) * 4, 1) + dims = np.fromfile(fid, dtype='>i', count=ndim + 1) if ndim != 2: - raise ValueError, 'Only two-dimensional matrices are ' \ - 'supported at this time' + raise ValueError('Only two-dimensional matrices are ' + 'supported at this time') # Back to where the data start fid.seek(pos, 0) @@ -176,18 +176,18 @@ def read_tag(fid, pos=None): # CCS sparse.csc_matrix() sparse_indices = np.fromfile(fid, dtype='>i4', count=nnz) - sparse_ptrs = np.fromfile(fid, dtype='>i4', count=ncol+1) + sparse_ptrs = np.fromfile(fid, dtype='>i4', count=ncol + 1) tag.data = sparse.csc_matrix((sparse_data, sparse_indices, sparse_ptrs), shape=dims) else: # RCS sparse_indices = np.fromfile(fid, dtype='>i4', count=nnz) - sparse_ptrs = np.fromfile(fid, dtype='>i4', count=nrow+1) + sparse_ptrs = np.fromfile(fid, dtype='>i4', count=nrow + 1) tag.data = sparse.csr_matrix((sparse_data, sparse_indices, sparse_ptrs), shape=dims) else: - raise ValueError, 'Cannot handle other than dense or sparse ' \ - 'matrices yet' + raise ValueError('Cannot handle other than dense or sparse ' + 'matrices yet') else: # All other data types @@ -195,35 +195,37 @@ def read_tag(fid, pos=None): if tag.type == FIFF.FIFFT_BYTE: tag.data = np.fromfile(fid, dtype=">B1", count=tag.size) elif tag.type == FIFF.FIFFT_SHORT: - tag.data = np.fromfile(fid, dtype=">h2", count=tag.size/2) + tag.data = np.fromfile(fid, dtype=">h2", count=tag.size / 2) elif tag.type == FIFF.FIFFT_INT: - tag.data = np.fromfile(fid, dtype=">i4", count=tag.size/4) + tag.data = np.fromfile(fid, dtype=">i4", count=tag.size / 4) elif tag.type == FIFF.FIFFT_USHORT: - tag.data = np.fromfile(fid, dtype=">H2", count=tag.size/2) + tag.data = np.fromfile(fid, dtype=">H2", count=tag.size / 2) elif tag.type == FIFF.FIFFT_UINT: - tag.data = np.fromfile(fid, dtype=">I4", count=tag.size/4) + tag.data = np.fromfile(fid, dtype=">I4", count=tag.size / 4) elif tag.type == FIFF.FIFFT_FLOAT: - tag.data = np.fromfile(fid, dtype=">f4", count=tag.size/4) + tag.data = np.fromfile(fid, dtype=">f4", count=tag.size / 4) elif tag.type == FIFF.FIFFT_DOUBLE: - tag.data = np.fromfile(fid, dtype=">f8", count=tag.size/8) + tag.data = np.fromfile(fid, dtype=">f8", count=tag.size / 8) elif tag.type == FIFF.FIFFT_STRING: tag.data = np.fromfile(fid, dtype=">c", count=tag.size) tag.data = ''.join(tag.data) elif tag.type == FIFF.FIFFT_DAU_PACK16: - tag.data = np.fromfile(fid, dtype=">h2", count=tag.size/2) + tag.data = np.fromfile(fid, dtype=">h2", count=tag.size / 2) elif tag.type == FIFF.FIFFT_COMPLEX_FLOAT: - tag.data = np.fromfile(fid, dtype=">f4", count=tag.size/4) + tag.data = np.fromfile(fid, dtype=">f4", count=tag.size / 4) tag.data = tag.data[::2] + 1j * tag.data[1::2] elif tag.type == FIFF.FIFFT_COMPLEX_DOUBLE: - tag.data = np.fromfile(fid, dtype=">f8", count=tag.size/8) + tag.data = np.fromfile(fid, dtype=">f8", count=tag.size / 8) tag.data = tag.data[::2] + 1j * tag.data[1::2] # # Structures # elif tag.type == FIFF.FIFFT_ID_STRUCT: tag.data = dict() - tag.data['version'] = int(np.fromfile(fid, dtype=">i4", count=1)) - tag.data['version'] = int(np.fromfile(fid, dtype=">i4", count=1)) + tag.data['version'] = int(np.fromfile(fid, dtype=">i4", + count=1)) + tag.data['version'] = int(np.fromfile(fid, dtype=">i4", + count=1)) tag.data['machid'] = np.fromfile(fid, dtype=">i4", count=2) tag.data['secs'] = int(np.fromfile(fid, dtype=">i4", count=1)) tag.data['usecs'] = int(np.fromfile(fid, dtype=">i4", count=1)) @@ -245,22 +247,23 @@ def read_tag(fid, pos=None): # Skip over the inverse transformation # It is easier to just use inverse of trans in Matlab # - fid.seek(12*4, 1) + fid.seek(12 * 4, 1) elif tag.type == FIFF.FIFFT_CH_INFO_STRUCT: - tag.data = Bunch() - tag.data['scanno'] = int(np.fromfile(fid, dtype=">i4", count=1)) - tag.data['logno'] = int(np.fromfile(fid, dtype=">i4", count=1)) - tag.data['kind'] = int(np.fromfile(fid, dtype=">i4", count=1)) - tag.data['range'] = float(np.fromfile(fid, dtype=">f4", count=1)) - tag.data['cal'] = float(np.fromfile(fid, dtype=">f4", count=1)) - tag.data['coil_type'] = int(np.fromfile(fid, dtype=">i4", count=1)) + d = Bunch() + d['scanno'] = int(np.fromfile(fid, dtype=">i4", count=1)) + d['logno'] = int(np.fromfile(fid, dtype=">i4", count=1)) + d['kind'] = int(np.fromfile(fid, dtype=">i4", count=1)) + d['range'] = float(np.fromfile(fid, dtype=">f4", count=1)) + d['cal'] = float(np.fromfile(fid, dtype=">f4", count=1)) + d['coil_type'] = int(np.fromfile(fid, dtype=">i4", count=1)) # # Read the coil coordinate system definition # - tag.data['loc'] = np.fromfile(fid, dtype=">f4", count=12) - tag.data['coil_trans'] = None - tag.data['eeg_loc'] = None - tag.data['coord_frame'] = FIFF.FIFFV_COORD_UNKNOWN + d['loc'] = np.fromfile(fid, dtype=">f4", count=12) + d['coil_trans'] = None + d['eeg_loc'] = None + d['coord_frame'] = FIFF.FIFFV_COORD_UNKNOWN + tag.data = d # # Convert loc into a more useful format # @@ -281,7 +284,8 @@ def read_tag(fid, pos=None): # Unit and exponent # tag.data['unit'] = int(np.fromfile(fid, dtype=">i4", count=1)) - tag.data['unit_mul'] = int(np.fromfile(fid, dtype=">i4", count=1)) + tag.data['unit_mul'] = int(np.fromfile(fid, dtype=">i4", + count=1)) # # Handle the channel name # @@ -295,19 +299,20 @@ def read_tag(fid, pos=None): elif tag.type == FIFF.FIFFT_OLD_PACK: offset = float(np.fromfile(fid, dtype=">f4", count=1)) scale = float(np.fromfile(fid, dtype=">f4", count=1)) - tag.data = np.fromfile(fid, dtype=">h2", count=(tag.size-8)/2) - tag.data = scale*tag.data + offset + tag.data = np.fromfile(fid, dtype=">h2", + count=(tag.size - 8) / 2) + tag.data = scale * tag.data + offset elif tag.type == FIFF.FIFFT_DIR_ENTRY_STRUCT: tag.data = list() - for _ in range(tag.size/16-1): - s = fid.read(4*4) + for _ in range(tag.size / 16 - 1): + s = fid.read(4 * 4) tag.data.append(Tag(*struct.unpack(">iIii", s))) else: - raise ValueError, 'Unimplemented tag data type %s' % tag.type + raise ValueError('Unimplemented tag data type %s' % tag.type) if tag.next != FIFF.FIFFV_NEXT_SEQ: # f.seek(tag.next,0) - fid.seek(tag.next, 1) # XXX : fix? pb when tag.next < 0 + fid.seek(tag.next, 1) # XXX : fix? pb when tag.next < 0 return tag diff --git a/mne/fiff/write.py b/mne/fiff/write.py index 51e0f61..c7cf2d3 100755 --- a/mne/fiff/write.py +++ b/mne/fiff/write.py @@ -93,14 +93,14 @@ def write_id(fid, kind, id_=None): if id_ is None: id_ = dict() id_['version'] = (1 << 16) | 2 - id_['machid'] = 65536 * np.random.rand(2) # Machine id (andom for now) + id_['machid'] = 65536 * np.random.rand(2) # Machine id (andom for now) id_['secs'] = time.time() - id_['usecs'] = 0 # Do not know how we could get this XXX + id_['usecs'] = 0 # Do not know how we could get this XXX FIFFT_ID_STRUCT = 31 FIFFV_NEXT_SEQ = 0 - data_size = 5 * 4 # The id comprises five integers + data_size = 5 * 4 # The id comprises five integers fid.write(np.array(kind, dtype='>i4').tostring()) fid.write(np.array(FIFFT_ID_STRUCT, dtype='>i4').tostring()) fid.write(np.array(data_size, dtype='>i4').tostring()) @@ -179,7 +179,7 @@ def write_coord_trans(fid, trans): # fiff_float_t invmove[3]; /*!< The inverse transform (translation part) */ #} *fiffCoordTrans, fiffCoordTransRec; /*!< Coordinate transformation descriptor */ - data_size = 4*2*12 + 4*2 + data_size = 4 * 2 * 12 + 4 * 2 fid.write(np.array(FIFF_COORD_TRANS, dtype='>i4').tostring()) fid.write(np.array(FIFFT_COORD_TRANS_STRUCT, dtype='>i4').tostring()) fid.write(np.array(data_size, dtype='>i4').tostring()) @@ -216,7 +216,6 @@ def write_ch_info(fid, ch): # fiff_float_t ez[3]; /*!< Coil coordinate system z-axis unit vector */ #} fiffChPosRec,*fiffChPos; /*!< Measurement channel position and coil type */ - #typedef struct _fiffChInfoRec { # fiff_int_t scanNo; /*!< Scanning order # */ # fiff_int_t logNo; /*!< Logical channel # */ @@ -229,7 +228,7 @@ def write_ch_info(fid, ch): # fiff_char_t ch_name[16]; /*!< Descriptive name for the channel */ #} fiffChInfoRec,*fiffChInfo; /*!< Description of one channel */ - data_size = 4*13 + 4*7 + 16 + data_size = 4 * 13 + 4 * 7 + 16 fid.write(np.array(FIFF_CH_INFO, dtype='>i4').tostring()) fid.write(np.array(FIFFT_CH_INFO_STRUCT, dtype='>i4').tostring()) @@ -243,7 +242,7 @@ def write_ch_info(fid, ch): fid.write(np.array(ch['range'], dtype='>f4').tostring()) fid.write(np.array(ch['cal'], dtype='>f4').tostring()) fid.write(np.array(ch['coil_type'], dtype='>i4').tostring()) - fid.write(np.array(ch['loc'], dtype='>f4').tostring()) # writing 12 values + fid.write(np.array(ch['loc'], dtype='>f4').tostring()) # writing 12 values # unit and unit multiplier fid.write(np.array(ch['unit'], dtype='>i4').tostring()) diff --git a/mne/forward.py b/mne/forward.py index 28b08a2..815a6d5 100755 --- a/mne/forward.py +++ b/mne/forward.py @@ -40,23 +40,23 @@ def _block_diag(A, n): """ from scipy import sparse - if not sparse.issparse(A): # then make block sparse + if not sparse.issparse(A): # then make block sparse ma, na = A.shape - bdn = na / int(n) # number of submatrices + bdn = na / int(n) # number of submatrices if na % n > 0: - raise ValueError, 'Width of matrix must be a multiple of n' + raise ValueError('Width of matrix must be a multiple of n') - tmp = np.arange(ma*bdn, dtype=np.int).reshape(bdn, ma) + tmp = np.arange(ma * bdn, dtype=np.int).reshape(bdn, ma) tmp = np.tile(tmp, (1, n)) ii = tmp.ravel() - jj = np.arange(na, dtype=np.int)[None,:] - jj = jj * np.ones(ma, dtype=np.int)[:,None] - jj = jj.T.ravel() # column indices foreach sparse bd + jj = np.arange(na, dtype=np.int)[None, :] + jj = jj * np.ones(ma, dtype=np.int)[:, None] + jj = jj.T.ravel() # column indices foreach sparse bd bd = sparse.coo_matrix((A.T.ravel(), np.c_[ii, jj].T)).tocsc() - else: # already is sparse, unblock it + else: # already is sparse, unblock it import pdb; pdb.set_trace() # [mA,na] = size(A); % matrix always has na columns # % how many entries in the first column? @@ -94,26 +94,26 @@ def _read_one(fid, node): tag = find_tag(fid, node, FIFF.FIFF_MNE_SOURCE_ORIENTATION) if tag is None: fid.close() - raise ValueError, 'Source orientation tag not found' + raise ValueError('Source orientation tag not found') one = dict() one['source_ori'] = tag.data tag = find_tag(fid, node, FIFF.FIFF_MNE_COORD_FRAME) if tag is None: fid.close() - raise ValueError, 'Coordinate frame tag not found' + raise ValueError('Coordinate frame tag not found') one['coord_frame'] = tag.data tag = find_tag(fid, node, FIFF.FIFF_MNE_SOURCE_SPACE_NPOINTS) if tag is None: fid.close() - raise ValueError, 'Number of sources not found' + raise ValueError('Number of sources not found') one['nsource'] = tag.data tag = find_tag(fid, node, FIFF.FIFF_NCHAN) if tag is None: fid.close() - raise ValueError, 'Number of channels not found' + raise ValueError('Number of channels not found') one['nchan'] = tag.data try: @@ -133,17 +133,17 @@ def _read_one(fid, node): if one['sol']['data'].shape[0] != one['nchan'] or \ (one['sol']['data'].shape[1] != one['nsource'] and - one['sol']['data'].shape[1] != 3*one['nsource']): + one['sol']['data'].shape[1] != 3 * one['nsource']): fid.close() - raise ValueError, 'Forward solution matrix has wrong dimensions' + raise ValueError('Forward solution matrix has wrong dimensions') if one['sol_grad'] is not None: if one['sol_grad']['data'].shape[0] != one['nchan'] or \ - (one['sol_grad']['data'].shape[1] != 3*one['nsource'] and - one['sol_grad']['data'].shape[1] != 3*3*one['nsource']): + (one['sol_grad']['data'].shape[1] != 3 * one['nsource'] and + one['sol_grad']['data'].shape[1] != 3 * 3 * one['nsource']): fid.close() - raise ValueError, 'Forward solution gradient matrix has ' \ - 'wrong dimensions' + raise ValueError('Forward solution gradient matrix has ' + 'wrong dimensions') return one @@ -186,26 +186,26 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, fwds = dir_tree_find(tree, FIFF.FIFFB_MNE_FORWARD_SOLUTION) if len(fwds) == 0: fid.close() - raise ValueError, 'No forward solutions in %s' % fname + raise ValueError('No forward solutions in %s' % fname) # Parent MRI data parent_mri = dir_tree_find(tree, FIFF.FIFFB_MNE_PARENT_MRI_FILE) if len(fwds) == 0: fid.close() - raise ValueError, 'No parent MRI information in %s' % fname + raise ValueError('No parent MRI information in %s' % fname) parent_mri = parent_mri[0] # Parent MEG data parent_meg = dir_tree_find(tree, FIFF.FIFFB_MNE_PARENT_MEAS_FILE) if len(parent_meg) == 0: fid.close() - raise ValueError, 'No parent MEG information in %s' % fname + raise ValueError('No parent MEG information in %s' % fname) parent_meg = parent_meg[0] chs = list() for k in range(parent_meg['nent']): kind = parent_meg['directory'][k].kind - pos = parent_meg['directory'][k].pos + pos = parent_meg['directory'][k].pos if kind == FIFF.FIFF_CH_INFO: tag = read_tag(fid, pos) chs.append(tag.data) @@ -214,7 +214,7 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, src = read_source_spaces_from_tree(fid, tree, add_geom=False) except Exception as inst: fid.close() - raise ValueError, 'Could not read the source spaces (%s)' % inst + raise ValueError('Could not read the source spaces (%s)' % inst) for s in src: s['id'] = find_source_space_hemi(s) @@ -228,8 +228,8 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, tag = find_tag(fid, fwds[k], FIFF.FIFF_MNE_INCLUDED_METHODS) if tag is None: fid.close() - raise ValueError, 'Methods not listed for one of the forward ' \ - 'solutions' + raise ValueError('Methods not listed for one of the forward ' + 'solutions') if tag.data == FIFF.FIFFV_MNE_MEG: megnode = fwds[k] @@ -263,7 +263,7 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, megfwd['nsource'] != eegfwd['nsource'] or megfwd['coord_frame'] != eegfwd['coord_frame']): fid.close() - raise ValueError, 'The MEG and EEG forward solutions do not match' + raise ValueError('The MEG and EEG forward solutions do not match') fwd = megfwd fwd['sol']['data'] = np.r_[fwd['sol']['data'], eegfwd['sol']['data']] @@ -293,17 +293,17 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, tag = find_tag(fid, parent_mri, FIFF.FIFF_COORD_TRANS) if tag is None: fid.close() - raise ValueError, 'MRI/head coordinate transformation not found' + raise ValueError('MRI/head coordinate transformation not found') else: - mri_head_t = tag.data; + mri_head_t = tag.data if (mri_head_t['from'] != FIFF.FIFFV_COORD_MRI or mri_head_t['to'] != FIFF.FIFFV_COORD_HEAD): mri_head_t = invert_transform(mri_head_t) if (mri_head_t['from'] != FIFF.FIFFV_COORD_MRI or mri_head_t['to'] != FIFF.FIFFV_COORD_HEAD): fid.close() - raise ValueError, 'MRI/head coordinate transformation not ' \ - 'found' + raise ValueError('MRI/head coordinate transformation not ' + 'found') fid.close() @@ -314,24 +314,24 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, if (fwd['coord_frame'] != FIFF.FIFFV_COORD_MRI and fwd['coord_frame'] != FIFF.FIFFV_COORD_HEAD): - raise ValueError, 'Only forward solutions computed in MRI or head ' \ - 'coordinates are acceptable' + raise ValueError('Only forward solutions computed in MRI or head ' + 'coordinates are acceptable') nuse = 0 for s in src: try: s = transform_source_space_to(s, fwd['coord_frame'], mri_head_t) except Exception as inst: - raise ValueError, 'Could not transform source space (%s)' % inst + raise ValueError('Could not transform source space (%s)' % inst) nuse += s['nuse'] if nuse != fwd['nsource']: - raise ValueError, 'Source spaces do not match the forward solution.' + raise ValueError('Source spaces do not match the forward solution.') print '\tSource spaces transformed to the forward solution ' \ 'coordinate frame' - fwd['src'] = src; + fwd['src'] = src # Handle the source locations and orientations if fwd['source_ori'] == FIFF.FIFFV_MNE_FIXED_ORI or force_fixed == True: @@ -339,8 +339,10 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, fwd['source_rr'] = np.zeros((fwd['nsource'], 3)) fwd['source_nn'] = np.zeros((fwd['nsource'], 3)) for s in src: - fwd['source_rr'][nuse:nuse+s['nuse'],:] = s['rr'][s['vertno'],:] - fwd['source_nn'][nuse:nuse+s['nuse'],:] = s['nn'][s['vertno'],:] + fwd['source_rr'][nuse:nuse + s['nuse'], :] = \ + s['rr'][s['vertno'], :] + fwd['source_nn'][nuse:nuse + s['nuse'], :] = \ + s['nn'][s['vertno'], :] nuse += s['nuse'] # Modify the forward solution for fixed source orientations @@ -354,7 +356,7 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, if fwd['sol_grad'] is not None: fwd['sol_grad']['data'] = np.dot(fwd['sol_grad']['data'], np.kron(fix_rot, np.eye(3))) - fwd['sol_grad']['ncol'] = 3*fwd['nsource'] + fwd['sol_grad']['ncol'] = 3 * fwd['nsource'] print '[done]' else: @@ -365,20 +367,20 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, pp = 0 nuse_total = sum([s['nuse'] for s in src]) fwd['source_rr'] = np.zeros((fwd['nsource'], 3)) - fwd['source_nn'] = np.empty((3*nuse_total, 3), dtype=np.float) + fwd['source_nn'] = np.empty((3 * nuse_total, 3), dtype=np.float) for s in src: - fwd['source_rr'][nuse:nuse + s['nuse'],:] = \ - s['rr'][s['vertno'],:] + fwd['source_rr'][nuse:nuse + s['nuse'], :] = \ + s['rr'][s['vertno'], :] for p in range(s['nuse']): # Project out the surface normal and compute SVD - nn = s['nn'][s['vertno'][p],:].T - nn = nn[:,None] - U, S, _ = linalg.svd(np.eye(3,3) - nn*nn.T) + nn = s['nn'][s['vertno'][p], :].T + nn = nn[:, None] + U, S, _ = linalg.svd(np.eye(3, 3) - nn * nn.T) # Make sure that ez is in the direction of nn - if np.sum(nn*U[:,2]) < 0: + if np.sum(nn * U[:, 2]) < 0: U *= -1 - fwd['source_nn'][pp:pp+3,:] = U.T + fwd['source_nn'][pp:pp + 3, :] = U.T pp += 3 nuse += s['nuse'] @@ -395,8 +397,8 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, nuse = 0 fwd['source_rr'] = np.zeros((fwd['nsource'], 3)) for s in src: - fwd['source_rr'][nuse:nuse+s['nuse'],:] = \ - s['rr'][s['vertno'],:] + fwd['source_rr'][nuse:nuse + s['nuse'], :] = \ + s['rr'][s['vertno'], :] nuse += s['nuse'] fwd['source_nn'] = np.kron(np.ones((fwd['nsource'], 1)), np.eye(3)) @@ -408,4 +410,3 @@ def read_forward_solution(fname, force_fixed=False, surf_ori=False, fwd = pick_channels_forward(fwd, include=include, exclude=exclude) return fwd - diff --git a/mne/inverse.py b/mne/inverse.py index e7a596d..a7a8a2a 100755 --- a/mne/inverse.py +++ b/mne/inverse.py @@ -230,11 +230,11 @@ def read_inverse_operator(fname): # # Some empty fields to be filled in later # - inv['proj'] = [] # This is the projector to apply to the data - inv['whitener'] = [] # This whitens the data - inv['reginv'] = [] # This the diagonal matrix implementing - # regularization and the inverse - inv['noisenorm'] = [] # These are the noise-normalization factors + inv['proj'] = [] # This is the projector to apply to the data + inv['whitener'] = [] # This whitens the data + inv['reginv'] = [] # This the diagonal matrix implementing + # regularization and the inverse + inv['noisenorm'] = [] # These are the noise-normalization factors # nuse = 0 for k in range(len(inv['src'])): @@ -259,6 +259,7 @@ def read_inverse_operator(fname): ############################################################################### # Compute inverse solution + def combine_xyz(vec): """Compute the three Cartesian components of a vector or matrix together @@ -279,7 +280,7 @@ def combine_xyz(vec): raise ValueError('Input must have 3N rows') n, p = vec.shape - return np.sqrt((np.abs(vec).reshape(n/3, 3, p)**2).sum(axis=1)) + return np.sqrt((np.abs(vec).reshape(n / 3, 3, p) ** 2).sum(axis=1)) def prepare_inverse_operator(orig, nave, lambda2, dSPM): @@ -318,14 +319,13 @@ def prepare_inverse_operator(orig, nave, lambda2, dSPM): if inv['eigen_leads_weighted']: inv['eigen_leads']['data'] = sqrt(scale) * inv['eigen_leads']['data'] - print ('\tScaled noise and source covariance from nave = %d to ' 'nave = %d' % (inv['nave'], nave)) inv['nave'] = nave # # Create the diagonal matrix for computing the regularized inverse # - inv['reginv'] = inv['sing'] / (inv['sing']**2 + lambda2) + inv['reginv'] = inv['sing'] / (inv['sing'] ** 2 + lambda2) print '\tCreated the regularized inverter' # # Create the projection operator @@ -376,12 +376,12 @@ def prepare_inverse_operator(orig, nave, lambda2, dSPM): if inv['eigen_leads_weighted']: for k in range(inv['eigen_leads']['nrow']): one = inv['eigen_leads']['data'][k, :] * inv['reginv'] - noise_norm[k] = sqrt(np.sum(one**2)) + noise_norm[k] = sqrt(np.sum(one ** 2)) else: for k in range(inv['eigen_leads']['nrow']): one = sqrt(inv['source_cov']['data'][k]) * \ inv['eigen_leads']['data'][k, :] * inv['reginv'] - noise_norm[k] = sqrt(np.sum(one**2)) + noise_norm[k] = sqrt(np.sum(one ** 2)) # # Compute the final result @@ -395,7 +395,7 @@ def prepare_inverse_operator(orig, nave, lambda2, dSPM): # Even in this case return only one noise-normalization factor # per source location # - noise_norm = combine_xyz(noise_norm[:,None]).ravel() + noise_norm = combine_xyz(noise_norm[:, None]).ravel() inv['noisenorm'] = 1.0 / np.abs(noise_norm) print '[done]' @@ -522,13 +522,13 @@ def _xyz2lf(Lf_xyz, normals): Lf_cortex = np.zeros_like(Lf_xyz) for k in range(n_positions): - lf_normal = np.dot(Lf_xyz[:,k,:], normals[k]) - lf_normal_n = lf_normal[:,None] / linalg.norm(lf_normal) - P = np.eye(n_sensors,n_sensors) - np.dot(lf_normal_n, lf_normal_n.T) - lf_p = np.dot(P, Lf_xyz[:,k,:]) + lf_normal = np.dot(Lf_xyz[:, k, :], normals[k]) + lf_normal_n = lf_normal[:, None] / linalg.norm(lf_normal) + P = np.eye(n_sensors, n_sensors) - np.dot(lf_normal_n, lf_normal_n.T) + lf_p = np.dot(P, Lf_xyz[:, k, :]) U, s, Vh = linalg.svd(lf_p) - Lf_cortex[:,k,0] = lf_normal - Lf_cortex[:,k,1:] = np.c_[U[:,0]*s[0], U[:,1]*s[1]] + Lf_cortex[:, k, 0] = lf_normal + Lf_cortex[:, k, 1:] = np.c_[U[:, 0] * s[0], U[:, 1] * s[1]] Lf_cortex = Lf_cortex.reshape(n_sensors, n_dipoles) return Lf_cortex @@ -596,7 +596,7 @@ def minimum_norm(evoked, forward, whitener, method='dspm', raise ValueError('weight_exp should be a scalar between 0 and 1') # Set regularization parameter based on SNR - lambda2 = 1.0 / snr**2 + lambda2 = 1.0 / snr ** 2 normals = [] for s in forward['src']: @@ -630,14 +630,14 @@ def minimum_norm(evoked, forward, whitener, method='dspm', # compute power if depth: - w = np.sum(gain**2, axis=0) + w = np.sum(gain ** 2, axis=0) w = w.reshape(-1, 3).sum(axis=1) - w = w[:,None] * np.ones((1, n_dip_per_pos)) + w = w[:, None] * np.ones((1, n_dip_per_pos)) w = w.ravel() if orientation is 'fixed': print 'Appying fixed dipole orientations.' - gain = gain * _block_diag(normals.ravel()[None,:], 3).T + gain = gain * _block_diag(normals.ravel()[None, :], 3).T elif orientation is 'free': print 'Using free dipole orientations. No constraints.' elif orientation is 'loose': @@ -666,7 +666,7 @@ def minimum_norm(evoked, forward, whitener, method='dspm', # apply weight limit # Applying weight limit. print 'Applying weight limit.' - weight_limit2 = weight_limit**2 + weight_limit2 = weight_limit ** 2 # limit = min(w(w>min(w) * weight_limit2)); % This is the Matti way. # we do the Rey way (robust to possible weight discontinuity). limit = np.min(w) * weight_limit2 @@ -690,44 +690,44 @@ def minimum_norm(evoked, forward, whitener, method='dspm', # Adjusting Source Covariance matrix to make trace of L*C_J*L' equal # to number of sensors. print 'Adjusting source covariance matrix.' - source_std = np.sqrt(w) # sqrt(C_J) - trclcl = linalg.norm(gain * source_std[None,:], ord='fro') - source_std *= sqrt(rank_noise) / trclcl # correct C_J - gain *= source_std[None,:] + source_std = np.sqrt(w) # sqrt(C_J) + trclcl = linalg.norm(gain * source_std[None, :], ord='fro') + source_std *= sqrt(rank_noise) / trclcl # correct C_J + gain *= source_std[None, :] # Compute SVD. print 'Computing SVD of whitened and weighted lead field matrix.' U, s, Vh = linalg.svd(gain, full_matrices=False) - ss = s / (s**2 + lambda2) + ss = s / (s ** 2 + lambda2) # Compute whitened MNE operator. - Kernel = source_std[:,None] * np.dot(Vh.T, ss[:,None] * U.T) + Kernel = source_std[:, None] * np.dot(Vh.T, ss[:, None] * U.T) # Compute dSPM operator. if method is 'dspm': print 'Computing dSPM inverse operator.' - dspm_diag = np.sum(Kernel**2, axis=1) + dspm_diag = np.sum(Kernel ** 2, axis=1) if n_dip_per_pos == 1: dspm_diag = np.sqrt(dspm_diag) elif n_dip_per_pos in [2, 3]: dspm_diag = dspm_diag.reshape(-1, n_dip_per_pos) dspm_diag = np.sqrt(np.sum(dspm_diag, axis=1)) dspm_diag = (np.ones((1, n_dip_per_pos)) * - dspm_diag[:,None]).ravel() + dspm_diag[:, None]).ravel() - Kernel /= dspm_diag[:,None] + Kernel /= dspm_diag[:, None] # whitened sLORETA imaging kernel elif method is 'sloreta': print 'Computing sLORETA inverse operator.' if n_dip_per_pos == 1: sloreta_diag = np.sqrt(np.sum(Kernel * gain.T, axis=1)) - Kernel /= sloreta_diag[:,None] + Kernel /= sloreta_diag[:, None] elif n_dip_per_pos in [2, 3]: for k in n_positions: - start = k*n_dip_per_pos - stop = (k+1)*n_dip_per_pos - R = np.dot(Kernel[start:stop,:], gain[:,start:stop]) + start = k * n_dip_per_pos + stop = (k + 1) * n_dip_per_pos + R = np.dot(Kernel[start:stop, :], gain[:, start:stop]) SIR = linalg.matfuncs.sqrtm(R, linalg.pinv(R)) Kernel[start:stop] = np.dot(SIR, Kernel[start:stop]) diff --git a/mne/preprocessing.py b/mne/preprocessing.py index d7d31e5..1ad326f 100755 --- a/mne/preprocessing.py +++ b/mne/preprocessing.py @@ -1,23 +1,24 @@ import numpy as np from .fiff.proj import make_projector_info -from .fiff.compensator import get_current_comp #, compensate_to, make_compensator +from .fiff.compensator import get_current_comp +# from .fiff.compensator import compensate_to, make_compensator # XXX # def cancel_noise(data, dest_comp=0): # """Do projection and compensation as needed -# +# # Return the appropriate operators -# +# # [res,proj,comp] = mne_ex_cancel_noise(data,dest_comp) -# +# # res - Data after noise cancellation # proj - The projection operator applied # comp - The compensator which brings uncompensated data to the # desired compensation grade (will be useful in forward # calculations) -# +# # """ # # # # Compensate the data and make a compensator for forward modelling @@ -30,13 +31,13 @@ from .fiff.compensator import get_current_comp #, compensate_to, make_compensato # else: # res = compensate_to(data, dest_comp) # print 'The data are now compensated to grade %d.' % dest_comp -# +# # if dest_comp > 0: # comp = make_compensator(res['info'], 0, dest_comp) # print 'Appropriate forward operator compensator created.' # else: # print 'No forward operator compensator needed.' -# +# # # Do the projection # if data['info']['projs'] is None: # print 'No projector included with these data.' @@ -44,7 +45,7 @@ from .fiff.compensator import get_current_comp #, compensate_to, make_compensato # # Activate the projection items # for k in range(len(res['info']['projs'])): # res['info']['projs'][k]['active'] = True; -# +# # # Create the projector # proj, nproj = make_projector_info(res['info']) # if nproj == 0: @@ -54,5 +55,5 @@ from .fiff.compensator import get_current_comp #, compensate_to, make_compensato # print 'Created an SSP operator (subspace dimension = %d)' % nproj # res['evoked']['epochs'] = np.dot(proj, res['evoked']['epochs']) # print 'Projector applied to the data' -# +# # return res, proj, comp diff --git a/mne/transforms.py b/mne/transforms.py index 5de385d..6e6a895 100755 --- a/mne/transforms.py +++ b/mne/transforms.py @@ -4,6 +4,7 @@ from scipy import linalg # from .fiff import FIFF + def invert_transform(trans): """Invert a transformation between coordinate systems """ @@ -43,7 +44,7 @@ def transform_source_space_to(src, dest, trans): raise ValueError('Cannot transform the source space using this ' 'coordinate transformation') - t = trans['trans'][:3,:] + t = trans['trans'][:3, :] res = src res['coord_frame'] = dest -- Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-mne.git _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
