Imports


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

Definitions


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

Preprocessing

Input Params


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'))

Trial Params


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]

Epoch Data


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()

Stats


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 [ ]: