In [ ]:
#disable autosave functionality;
#This script is taxing enough, and autosaving tends to push it over the edge
# plus, autosaving seems to zero out the file before restoring it from the backup
# this means an autosave causing a crash will actually delete the file rather than saving it!!!
%autosave 0
In [ ]:
#basic imports
%pylab inline
import time
import pickle
import cPickle
import logging as L
L.basicConfig(level=L.ERROR) # INFO)
import time
import numpy
#import #pylab
import scipy.stats
import os
import pylab
import sklearn
import scipy
import sklearn.linear_model
import re
pylab.rcParams['figure.figsize'] = 10,10 #change the default image size for this session
pylab.ion()
In [ ]:
import mne
In [ ]:
#custom imports
%cd ../scripts
# brian's prototype routines
from protoMEEGutils import *
import protoSpectralWinFFTMapper as specfft
In [ ]:
#### SUBROUTINES ####
# plot the time-scale for ERFs and other epoch figures
def commonPlotProps():
#zeroSample = (abs(epochStart)/float(epochLength)*epochNumTimepoints)
#pylab.plot((0,epochNumTimepoints),(0,0),'k--')
#pylab.ylim((-2.5e-13,2.5e-13)) #((-5e-14,5e-14)) # better way just to get the ones it choose itself?
#pylab.plot((zeroSample,zeroSample),(0,0.01),'k--')
pylab.xticks(numpy.linspace(0,epochNumTimepoints,7),epochStart+(numpy.linspace(0,epochNumTimepoints,7)/samplingRate))
pylab.xlabel('time (s) relative to auditory onset') #+refEvent)
pylab.xlim((62,313))
pylab.show()
pylab.axhline(0, color='k', linestyle='--')
pylab.axvline(125, color='k', linestyle='--')
# adjust R2 down for the artificial inflation you get by increasing the number of explanatory features
def adjustR2(R2, numFeatures, numSamples):
#1/0
#return R2
return R2-(1-R2)*(float(numFeatures)/(numSamples-numFeatures-1))
# normalise (z-scale) the scale of variables (for the explanatory ones, so the magnitude of beta values are comparably interpretable)
def mynormalise(A):
A = scipy.stats.zscore(A)
A[numpy.isnan(A)] = 0
return A
In [ ]:
#change to the data directory to load in the data
%cd ../MEG_data
#Select subjects
(megFileTag1, megFile1) = ('V_TSSS_0.01-50Hz_@125', 'v_hod_allRuns_tsss_audiobookPrepro_stPad1_lp50_resamp125_frac10ICAed.set')#_hp0.010000.set')
(megFileTag2, megFile2) = ('A_TSSS_0.01-50Hz_@125', 'aud_hofd_a_allRuns_tsss_audiobookPrepro_stPad1_lp50_resamp125_frac10ICAed_hp0.010000.set')
(megFileTag3, megFile3) = ('C_TSSS_0.01-50Hz_@125', 'aud_hofd_c_allRuns_tsss_audiobookPrepro_stPad1_lp50_resamp125_frac10ICAed_hp0.010000.set')
## put your on properties in here, as a .tab file of similar format (tab delimited, and field names in a first comment line - should be easy to do in excel...)
#to get the V5.tab:
# python ../scripts/buildtab.py hod_JoeTimes_LoadsaFeaturesV3.tab hod.wsj02to21-comparativized-gcg15-1671-4sm.fullberk.parsed.gcgbadwords > hod_JoeTimes_LoadsaFeaturesV4.tab
# python ../scripts/addsentid.py hod_JoeTimes_LoadsaFeaturesV4.tab > hod_JoeTimes_LoadsaFeaturesV5.tab
tokenPropsFile = 'hod_JoeTimes_LoadsaFeaturesV5.tab' # best yet! have automatic Stanford tagging, several word length and freq measures, and also the 2 and 3 token back-grams
# WHICH CHANNELS TO LOOK AT AS ERFS
#decide which channels to use
channelLabels = ['MEG0111', 'MEG0121', 'MEG0131', 'MEG0211', 'MEG0212', 'MEG0213', 'MEG0341']
# LOAD WORD PROPS
# change dtype to suit the files in your .tab file
tokenProps = scipy.genfromtxt(tokenPropsFile,
delimiter='\t',names=True,
dtype="i4,f4,f4,S50,S50,i2,i2,i2,S10,f4,f4,f4,f4,f4,f4,f4,f4,f4,f4,f4,i1,>i4")
# ... and temporarily save as cpickle archive to satisfy the way I programmed the convenience function loadBookMEGWithAudio (it expects to find the same info in a C-pickle file, and so doesn't need to know about number and type of fields)
tokenPropsPickle = tokenPropsFile+'.cpk'
cPickle.dump(tokenProps, open(tokenPropsPickle, 'wb'))
In [ ]:
triggersOfInterest=['s%d' % i for i in range(1,10)]
refEvent = 'onTime' #,'offTime']
#*# guess an epoch of -0.5 to +1s should be enough #*#
epochStart = -1; # stimulus ref event
epochEnd = +2; #
epochLength = epochEnd-epochStart;
baseline = False #[-1,0]
In [ ]:
# Get the goods on subject 1
(contSignalData, metadata, trackTrials, tokenProps, audioSignal, samplingRate, numChannels) = loadBookMEGWithAudio(megFile1, tokenPropsPickle, triggersOfInterest, epochEnd, epochStart, icaComps=False)
del contSignalData
del trackTrials
In [ ]:
print(metadata.chanlocs[0]._fieldnames)
In [ ]:
#Get big sensor map
%cd ../scripts/
resultType = 'wide'
resultsFile = 'sigresults.dev.full.cpk' # % (resultType)
rsq = cPickle.load(open(resultsFile))
# EXTRACT SENSOR LOCATIONS FROM METADATA
sensorLocations = {}
# don't have the data at hand, but here it would be good to first extract the X/Y co-ordinates of all
for sensor in metadata.chanlocs:
# not sure if the structure is like this - will have to check in metadata returned by the python script
sensorLocations[sensor.labels] = (-sensor.Y, sensor.X)
#sensorLocations[sensor.labels] = (sensor.X, sensor.Y)
# DO A TOPOPLOT FOR EACH BAND AND STAT (three beside eachother, one for each type of sensor)
for stat in ['p']:#, 'max']:
for band in ['alpha']: # add more bands if you want
#sensorPositions = [] # in horizontal plane, looking downwards
#sensorValues = [] # r^2 or whatever other value we want to do map of
pylab.figure()
for sensorType in ['3']:
sensorPositions = [] # in horizontal plane, looking downwards
sensorValues = [] # r^2 or whatever other value we want to do map of
sensorLabels = []
# sensors that end in 1 are magnetometers, those ending in 2 or 3 are gradiometers of different orientations
for sensorLabel in rsq[stat].keys(): # 'ave' or 'max'? not sure
if sensorLabel.endswith(sensorType):
sensorValues = numpy.append(sensorValues,rsq[stat][sensorLabel][band]['syndepth']) # e.g. rsq['ave']['MEG2643']['alpha']
sensorPositions = numpy.append(sensorPositions,sensorLocations[sensorLabel]) # this is already an X-Y tuple
sensorLabels = numpy.append(sensorLabels,sensorLabel)
sensorValues = numpy.array(sensorValues)
sensorPositions = numpy.array(sensorPositions)
sensorLabels = numpy.array(sensorLabels)
#subplotCounter += 1
#pylab.subplot(1,3,subplotCounter)
#print(sensorValues.shape, sensorPositions.reshape(-1,2).shape)
if stat == 'p':
for i,s in enumerate(sensorValues):
if not np.isfinite(s):
sensorValues[i] = 1.0
im, cn = mne.viz.plot_topomap(sensorValues, sensorPositions.reshape(-1,2),
names=sensorLabels,show_names=True,
vmin=0.0,vmax=0.05)
#im, cn = mne.viz.plot_topomap(sensorValues, sensorPositions.reshape(-1,2))
pylab.colorbar(im,fraction=0.046, pad=0.04)
pylab.title('Sensor map')
#pylab.legend(loc='upper right')
#pylab.savefig('graphics/large-post-meg_%s-%s_%s.png' % (resultType,stat,band))
pylab.show()
In [ ]:
#Get blank sensor map
%cd ../scripts/
resultType = 'wide'
resultsFile = 'sigresults.dev.full.cpk' # % (resultType)
rsq = cPickle.load(open(resultsFile))
# EXTRACT SENSOR LOCATIONS FROM METADATA
sensorLocations = {}
# don't have the data at hand, but here it would be good to first extract the X/Y co-ordinates of all
for sensor in metadata.chanlocs:
# not sure if the structure is like this - will have to check in metadata returned by the python script
sensorLocations[sensor.labels] = (-sensor.Y, sensor.X)
#sensorLocations[sensor.labels] = (sensor.X, sensor.Y)
# DO A TOPOPLOT FOR EACH BAND AND STAT (three beside eachother, one for each type of sensor)
for stat in ['p']:#, 'max']:
for factor in ['sentenceSerial','syndepth','bigramLogProbBack_COCA','logFreq_ANC']: # add more bands if you want
#sensorPositions = [] # in horizontal plane, looking downwards
#sensorValues = [] # r^2 or whatever other value we want to do map of
pylab.figure()
subplotCounter = 0
for sensorType in ['1','2','3']:
sensorPositions = [] # in horizontal plane, looking downwards
sensorValues = [] # r^2 or whatever other value we want to do map of
sensorLabels = []
# sensors that end in 1 are magnetometers, those ending in 2 or 3 are gradiometers of different orientations
for sensorLabel in rsq[stat].keys(): # 'ave' or 'max'? not sure
if sensorLabel.endswith(sensorType):
sensorValues = numpy.append(sensorValues,rsq[stat][sensorLabel]['alpha'][factor]) # e.g. rsq['ave']['MEG2643']['alpha']
sensorPositions = numpy.append(sensorPositions,sensorLocations[sensorLabel]) # this is already an X-Y tuple
sensorLabels = numpy.append(sensorLabels,sensorLabel)
sensorValues = numpy.array(sensorValues)
sensorPositions = numpy.array(sensorPositions)
sensorLabels = numpy.array(sensorLabels)
subplotCounter += 1
pylab.subplot(1,3,subplotCounter)
#print(sensorValues.shape, sensorPositions.reshape(-1,2).shape)
if stat == 'p':
for i,s in enumerate(sensorValues):
if not np.isfinite(s):
sensorValues[i] = 1.0
#im, cn = mne.viz.plot_topomap(sensorValues, sensorPositions.reshape(-1,2),
# names=sensorLabels,show_names=True,
# vmin=0.0,vmax=0.05)
im, cn = mne.viz.plot_topomap(sensorValues, sensorPositions.reshape(-1,2),
vmin=0.0,vmax=0.05)
pylab.colorbar(im,fraction=0.046, pad=0.04)
pylab.title('%s map' % (factor))
#pylab.legend(loc='upper right')
#pylab.savefig('graphics/meg_%s-%s_%s.png' % (stat,'alpha',factor))
pylab.show()
In [ ]:
#Get blank sensor map
%cd ../scripts/
resultType = 'wide'
resultsFile = 'sigresults.dev.full.cpk' # % (resultType)
rsq = cPickle.load(open(resultsFile))
# EXTRACT SENSOR LOCATIONS FROM METADATA
sensorLocations = {}
# don't have the data at hand, but here it would be good to first extract the X/Y co-ordinates of all
for sensor in metadata.chanlocs:
# not sure if the structure is like this - will have to check in metadata returned by the python script
sensorLocations[sensor.labels] = (-sensor.Y, sensor.X)
#sensorLocations[sensor.labels] = (sensor.X, sensor.Y)
# DO A TOPOPLOT FOR EACH BAND AND STAT (three beside eachother, one for each type of sensor)
for stat in ['r2']:#, 'max']:
for factor in ['syndepth']: # add more bands if you want
#sensorPositions = [] # in horizontal plane, looking downwards
#sensorValues = [] # r^2 or whatever other value we want to do map of
pylab.figure()
subplotCounter = 0
for sensorType in ['1','2','3']:
sensorPositions = [] # in horizontal plane, looking downwards
sensorValues = [] # r^2 or whatever other value we want to do map of
sensorLabels = []
# sensors that end in 1 are magnetometers, those ending in 2 or 3 are gradiometers of different orientations
for sensorLabel in rsq[stat].keys(): # 'ave' or 'max'? not sure
if sensorLabel.endswith(sensorType):
sensorValues = numpy.append(sensorValues,rsq[stat][sensorLabel]['alpha'][0]) # e.g. rsq['ave']['MEG2643']['alpha']
sensorPositions = numpy.append(sensorPositions,sensorLocations[sensorLabel]) # this is already an X-Y tuple
sensorLabels = numpy.append(sensorLabels,sensorLabel)
sensorValues = numpy.array(sensorValues)
sensorPositions = numpy.array(sensorPositions)
sensorLabels = numpy.array(sensorLabels)
subplotCounter += 1
pylab.subplot(1,3,subplotCounter)
#print(sensorValues.shape, sensorPositions.reshape(-1,2).shape)
if stat == 'p':
for i,s in enumerate(sensorValues):
if not np.isfinite(s):
sensorValues[i] = 1.0
#im, cn = mne.viz.plot_topomap(sensorValues, sensorPositions.reshape(-1,2),
# names=sensorLabels,show_names=True,
# vmin=0.0,vmax=0.05)
im, cn = mne.viz.plot_topomap(sensorValues, sensorPositions.reshape(-1,2),
vmin=0.0)
pylab.colorbar(im,fraction=0.046, pad=0.04)
pylab.title('%s map' % (factor))
#pylab.legend(loc='upper right')
#pylab.savefig('graphics/meg_%s-%s_conf.png' % (stat,'alpha'))
pylab.show()
In [ ]:
%cd ../scripts/
DEV = True
CLUSTER = True
SHORT = False
if CLUSTER:
if DEV:
resultsFile = 'sigresults.dev.cluster.cpk'
else:
resultsFile = 'sigresults.test.cluster.cpk'
else:
if SHORT:
if DEV:
resultsFile = 'sigresults.dev.short.cpk'
else:
resultsFile = 'sigresults.test.short.cpk'
else:
if DEV:
resultsFile = 'sigresults.dev.full.cpk'
else:
resultsFile = 'sigresults.test.full.cpk'
rsq = cPickle.load(open(resultsFile))
print rsq.keys() #print queryable properties
if CLUSTER:
#output cluster coordinates
for cluster in rsq['clustersize']:
print cluster, rsq['clustersize'][cluster]
for stat in ['p']:
for band in ['theta','gamma']:
print 'Significance in', band
if not CLUSTER:
for sensor in rsq[stat]: #['0113','0123','0112','0122','1542','1532','1543','1533']:
#print 'MEG'+sensor, 'p =',rsq[stat]['MEG'+sensor][band]
print 'Sensor:', sensor, 'R2:', rsq['r2'][sensor][band]
print rsq[stat][sensor][band]
else:
#mainly we care about the bottom 6 clusters in a (3,3) cluster setup; the front of the helmet has more air since subj heads rest towards the back
for sensor in [(0,0),(1,0),(2,0),(0,1),(1,1),(2,1)]:
#for sensor in [(1,0)]: check back central cluster in a (x,3) setup
#print 'MEG'+sensor, 'p =',rsq[stat]['MEG'+sensor][band]
print 'Cluster:', sensor, 'R2:', rsq['r2'][sensor][band]
print rsq[stat][sensor][band]
In [ ]: