notebook.community

Edit and run



In [246]:

    
clear all



In [247]:

    
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy import io
import scipy.io as sio
%matplotlib inline 
import pylab
import csv
from Tkinter import Tk
from tkFileDialog import askopenfilename
from tkFileDialog import askdirectory
import nibabel as nb
from scipy import io
from nifti import NiftiImage
import nibabel as nb
from scipy.interpolate import interp1d



In [248]:

    
# from http://stackoverflow.com/questions/3579568/choosing-a-file-in-python-with-simple-dialog
from Tkinter import Tk
from tkFileDialog import askopenfilename

Tk().withdraw() # we don't want a full GUI, so keep the root window from appearing
filename = askopenfilename() # show an "Open" dialog box and return the path to the selected file
print(filename)









    



/home/sophie/Downloads/100106ss2on500regcregcU10sMpsfkf147Smith0_4TS.mat



In [249]:

    
Ua=sio.loadmat(filename)
DT=Ua['TSo']
DT.shape









    Out[249]:





(5426, 147)



In [250]:

    
# from http://stackoverflow.com/questions/3579568/choosing-a-file-in-python-with-simple-dialog
from Tkinter import Tk
from tkFileDialog import askopenfilename

Tk().withdraw() # we don't want a full GUI, so keep the root window from appearing
filename2 = askopenfilename() # show an "Open" dialog box and return the path to the selected file
print(filename2)









    



/home/sophie/Downloads/100106ss2on500regcregcU10sMpsfkf147Smith0_4IC.nii



In [251]:

    
img1 = nb.load(filename2)
data = img1.get_data()
S=data.shape
S









    Out[251]:





(185, 116, 10, 147)



In [252]:

    
Tk().withdraw() # we don't want a full GUI, so keep the root window from appearing
filename = askopenfilename() # show an "Open" dialog box and return the path to the selected file
print(filename)









    



/home/sophie/Downloads/100106ss1TimeFluoOn.mat



In [253]:

    
#Old version from csv
#with open(filename) as inputfile:
#    results = list(csv.reader(inputfile))
#Raw_timea=[float(results[i][2]) for i in (range(1,len(results)))]
#Raw_time=np.asarray(Raw_timea[0:len(results)])
#Raw_time.shape
#on=126
#off=6253
#Time_fluo=Raw_time[on:off]-Raw_time[on]

Ua=sio.loadmat(filename)
Time_fluo=Ua['TimeFluoOn']
Time_fluo.shape









    Out[253]:





(1, 5927)



In [254]:

    
Time_fluoICA=Time_fluo[:,500:5926]



In [255]:

    
Time_fluoICA.shape









    Out[255]:





(1, 5426)



In [256]:

    
#if stimulus
#Tstim=np.array(range(1,100*(int(Tmax)+1)))
#Tstim=Tstim/100
#Odor=np.zeros(len(Tstim*100))
#UV=np.zeros(len(Tstim*100))
#for i in range(0,7):
#   UV[10*100+i*600:100*10+i*600+300]=1
#    
#for i in range(0,5):    
#    Odor[(10+8*6+15)*100+i*600:(10+8*6+15)*100+i*600+200]=1
#plt.plot(Tstim,UV)
#plt.plot(Tstim,Odor)

Open behavior file



In [257]:

    
Tk().withdraw() # we don't want a full GUI, so keep the root window from appearing
filename = askopenfilename() # show an "Open" dialog box and return the path to the selected file
print(filename)
B=sio.loadmat(filename)









    



/home/sophie/Downloads/100106BehaviorAll.mat



In [258]:

    
B









    Out[258]:





{'T': array([[  0.00000000e+00,   1.00000000e-02,   2.00000000e-02, ...,
           1.17437000e+02,   1.17447000e+02,   1.17457000e+02]]),
 '__globals__': [],
 '__header__': 'MATLAB 5.0 MAT-file, Platform: PCWIN64, Created on: Wed Jun 22 20:59:02 2016',
 '__version__': '1.0',
 'behavior': array([[0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        ..., 
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0],
        [0, 0, 0, ..., 0, 0, 0]], dtype=uint8),
 'left': array([[ 0.33898797,  0.34203277,  0.53269537, ...,  0.17775845,
          0.37721715,  0.53259871]]),
 'right': array([[ 0.53327534,  0.37528394,  0.50041081, ...,  0.34899232,
          0.47793727,  0.40075395]]),
 'straight': array([[ 0.69513315,  0.60277415,  0.43878981, ...,  0.19008264,
          0.35194046,  0.67212798]])}



In [259]:

    
b=B['behavior']



In [260]:

    
FrontGroom=b[0,:]
RearGroom=b[1,:]
Rest=b[3,:]
Tbehavior=B['T']
Walk=b[2,:]
Weird=b[4,:]



In [261]:

    
Weird.shape









    Out[261]:





(11593,)



In [262]:

    
#Threshold ball movements



In [263]:

    
left=B['left'].T
right=B['right'].T
straight=B['straight'].T



In [264]:

    
plt.plot(left)









    Out[264]:





[<matplotlib.lines.Line2D at 0xa5a0250>]



In [265]:

    
max(left)









    Out[265]:





array([ 167.66497511])



In [266]:

    
#Fill gaps (1 if i 1 and i+4 1)



In [267]:

    
List=['left','right','straight','FrontGroom','RearGroom','Rest','Walk','Weird']



In [268]:

    
dataset=eval(List[0])
for i in range(len(List)):
    del dataset
    dataset=eval(List[i])
    for j in range(len(dataset)-3):
        if dataset[j]==1 and dataset[j+3]==1:
            dataset[j+1]=1
            dataset[j+2]=1
    exec(List[i]+'=dataset')



In [269]:

    
Tb=np.squeeze(Tbehavior)



In [270]:

    
Tfi=Time_fluoICA.T
Tf=Tfi[range(len(Time_fluoICA.T)-58)]
f=interp1d(Tb,FrontGroom)
FrontGroomint=np.squeeze(f(Tf))
f=interp1d(Tb,np.squeeze(RearGroom))
RearGroomint=f(Tf)
f=interp1d(Tb,np.squeeze(Rest))
Restint=f(Tf)



In [271]:

    
Tball=np.zeros(len(left))
for j in range(len(left)):
    Tball[j]=(Tb[j]+Tb[j+1])/2



In [272]:

    
f=interp1d(Tball,left.T)
leftint=np.squeeze(f(Tf))
f=interp1d(Tball,np.squeeze(right))
rightint=f(Tf)
f=interp1d(Tball,np.squeeze(straight))
straightint=f(Tf)



In [273]:

    
leftintbin=leftint
lowvals=leftint<50
leftintbin[lowvals]=0
highvals=leftint>50
leftintbin[highvals]=1

rightintbin=rightint
lowvals=rightint<50
rightintbin[lowvals]=0
highvals=rightint>50
rightintbin[highvals]=1

straightintbin=straightint
lowvals=straightint<50
straightintbin[lowvals]=0
highvals=straightint>50
straightintbin[highvals]=1



In [274]:

    
#Align times (verify what frames have been removed)



In [275]:

    
#interpolate for ball movements



In [276]:

    
#check where the first image is compared to extrapolated full light
Tb=Tb
max(Tb)









    Out[276]:





117.45699999999998



In [277]:

    
Tb.shape









    Out[277]:





(11593,)



In [278]:

    
max(np.squeeze(Time_fluoICA))









    Out[278]:





118.66205986795977



In [279]:

    
#falls between +/- 5ms. note that the offset for the video is about 6 ms later thanlast image



In [280]:

    
#Do regressions with delay



In [281]:

    
#Regression integrator straight



In [282]:

    
#scikit learn GUI



In [283]:

    
Demean=np.zeros(S)
Dmaps=np.zeros(S)
Dvar=np.zeros(S)
Var=np.zeros(S[3])
D2=np.zeros([S[0],S[1],5,S[3]])
Tvar=np.zeros(S[3])



In [284]:

    
for i in range(S[3]):
    Demean[:,:,:,i]=data[:,:,:,i]-np.mean(np.mean(np.mean(data[:,:,:,i],0),0),0)



In [285]:

    
for i in range(S[3]):
    Dsq=np.reshape(Demean[:,:,:,i],S[0]*S[1]*S[2])
    Var[i]=np.sqrt(np.var(Dsq))
    Dvar=Demean[:,:,:,i]/Var[i]
    Dmaps[:,:,:,i]=Dvar-2.5
    Tvar[i]=np.var(DT[i,:])
Dmaps[Dmaps<0]=0



In [286]:

    
Order=np.argsort(Var)[::-1]
datao=data[:,:,:,Order[:]]
Dmapso=Dmaps[:,:,:,Order[:]]
Varor=Var[Order]



In [287]:

    
my_cmap=plt.cm.jet
my_cmap.set_bad(alpha=0)
Good_ICs=np.zeros(S[3])
Label_ICs=[]
pylab.rcParams['figure.figsize'] = (13, 2.5)



In [288]:

    
Dtemp=np.mean(data,3)



In [289]:

    
Dtemp.shape









    Out[289]:





(185, 116, 10)



In [290]:

    
%%javascript
IPython.OutputArea.auto_scroll_threshold =4000;



In [291]:

    
if S[2]>5:
    Nstack=5
    Int100=[(i+1)*100/Nstack for i in range(Nstack)]
    Percs=np.percentile(range(S[2]),Int100)
    Indices=np.split(range(S[2]),Percs)
    D1=np.zeros([S[0],S[1],Nstack])
    Dmean=np.squeeze(data[:,:,range(Nstack),2])
    for i in range(Nstack):
        Vmean=np.mean(Dtemp[:,:,Indices[i]],2)
        Dmean[:,:,i]=Vmean
else:
    Nstack=S[2]
    D1=np.zeros([S[0],S[1],S[2]])
    Dmean=data[:,:,range(S[2])]  
    Dmean=np.squeeze(Dtemp[:,:,:])



In [292]:

    
plt.imshow(Dmean[:,:,1],cmap=plt.cm.gray)









    Out[292]:





<matplotlib.image.AxesImage at 0x7d36fd0>



In [293]:

    
Tb.shape









    Out[293]:





(11593,)



In [294]:

    
S









    Out[294]:





(185, 116, 10, 147)



In [295]:

    
RearGroomint.shape









    Out[295]:





(5368, 1)



In [296]:

    
Tf.shape









    Out[296]:





(5368, 1)



In [297]:

    
for j in range(S[3]):

    if S[2]>5:
        for i in range(Nstack):
            V=Dmaps[:,:,Indices[i],Order[j]]
            D1[:,:,i]=np.max(V,2)
        D2[:,:,:,j]=D1
        D1[D1==0]=np.nan
           
    else:
        for i in range(S[2]):
            V=Dmaps[:,:,i,Order[j]]
            D1[:,:,i]=V 
 

    print(j)
    for i in range(Nstack):
        plt.subplot(1,5,i+1)
        plt.imshow(Dmean[:,:,i],cmap=plt.cm.gray)
        plt.imshow(D1[:,:,i], cmap=my_cmap,interpolation='none')
        frame1 = plt.gca()
        frame1.axes.get_xaxis().set_visible(False)
        frame1.axes.get_yaxis().set_visible(False)
        
    plt.show()
    
    plt.plot(Tf,straightintbin*1.5)   
    plt.plot(Tf,leftintbin*1.5)
    plt.plot(Tf,rightintbin*1.5)  
    plt.plot(Tf,FrontGroomint.T)
    plt.plot(Tf,RearGroomint) 
    plt.plot(Time_fluoICA.T,DT[:,Order[j]]+2)
    plt.show()
    
    Label_ICs.append(raw_input())
    if Label_ICs[j]=='':
        Good_ICs[j]=0
    else:
        Good_ICs[j]=1



In [806]:

    
Tk().withdraw() 
filenamet = askopenfilename() 
print(filenamet)
nimt=NiftiImage(filenamet)
Dtemp=np.squeeze(nimt.data.T)
Dtemp.shape









    



()






    



---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-806-d6de7db42465> in <module>()
      2 filenamet = askopenfilename()
      3 print(filenamet)
----> 4 nimt=NiftiImage(filenamet)
      5 Dtemp=np.squeeze(nimt.data.T)
      6 Dtemp.shape

/usr/lib/pymodules/python2.7/nifti/image.pyc in __init__(self, source, header, load, **kwargs)
     76         """
     77         # setup all nifti header related stuff
---> 78         NiftiFormat.__init__(self, source, header, **kwargs)
     79 
     80         # where the data will go to

/usr/lib/pymodules/python2.7/nifti/format.pyc in __init__(self, source, header, loadmeta)
     83             raise ValueError, \
     84                   "Unsupported source type. Only NumPy arrays and filename " \
---> 85                   + "string are supported."
     86 
     87 

ValueError: Unsupported source type. Only NumPy arrays and filename string are supported.



In [298]:

    
Dmaps.shape









    Out[298]:





(185, 116, 10, 147)



In [299]:

    
set(Label_ICs)









    Out[299]:





{'',
 'A',
 'AMMC',
 'C',
 'DCB',
 'EB',
 'FB',
 'KC',
 'L',
 'LH',
 'MB',
 'PB',
 'PN',
 'PS',
 'SOG',
 'i',
 'o'}



In [300]:

    
if len(Label_ICs)<S[3]:
    for j in range(S[3]-len(Label_ICs)):
      Label_ICs.append('')



In [301]:

    
Dict={'a':0,'AL':0,'Co':1,'CO':1,'O':1,'o':1,'OG':2,'AVLP':3,'A':0,'G':3,'PN':3,'CA':4,'L':4,'LH':4,'lh':4,'SMP':5,'SLP':5,'l':6,'KC':7,'kc':7,'M':7,'MB':7,'mb':7,'MBON':7,'FB':9,'EB':10,'eb':10,'C':10,'c':10,'PB':9,'pb':9,'AMMC':12,'ammc':12,'V':5,'S':13,'SOG':13,'PS':13,'s':13,'DCB':14,'I':12,'i':12,'D':5,'mvt':15,'':15,'0':15}



In [302]:

    
Translated=[Dict[X] for X in Label_ICs]



In [303]:

    
G=Good_ICs.tolist();



In [304]:

    
len(Good_ICs)









    Out[304]:





147



In [305]:

    
G.count(1)









    Out[305]:





61



In [306]:

    
where_are_NaNs = np.isnan(D2)

D2[where_are_NaNs] = 0



In [369]:

    
if S[2]>5:
    Final_map=np.zeros([S[0],S[1],5,3])
    Fmaps=np.zeros([S[0],S[1],5,3])
else:
    Final_map=npLabel_ICs[Neworder[i]] 
    Fmaps=np.zeros([S[0],S[1],3])    
C=np.zeros([S[3],3])



In [370]:

    
D2.shape









    Out[370]:





(185, 116, 5, 147)



In [371]:

    
for j in range(S[3]):    
    #if Good_ICs[j]:
    if Translated[j]==12 or Translated[j]==13:
        C[j,:]=np.squeeze(np.random.rand(3,1))
        for k in range(3):
            M=np.max(np.squeeze(np.reshape(D2[:,:,:,j],S[0]*S[1]*5)))
            #M[M==0]=1
  #          Fmaps[:,:,:,k]=0.7*Dmaps[:,:,:,j]*C[j,k]/np.max(np.squeeze(np.reshape(Dmaps[:,:,:,j]*np.max(C[j,:],S[0]*S[1]*5))
            Fmaps[:,:,:,k]=0.7*D2[:,:,:,j]*C[j,k]/(M*np.max(C[j,:]))
        Final_map=Final_map+Fmaps



In [373]:

    
pylab.rcParams['figure.figsize'] = (18, 5)

#Final_map[Final_map<0.1]=np.NaN

if S[2]>5:
    N=Nstack
else:
    N=S[2]
for i in range(N):
    #if Good_ICs[j]:
        plt.subplot(1,N,i+1)
        plt.imshow(Dmean[:,:,i],cmap=plt.cm.gray)
        plt.imshow(Final_map[:,:,i],interpolation='none')
        #plt.imshow(Final_map[:,:,i]) 
        frame1 = plt.gca()
        frame1.axes.get_xaxis().set_visible(False)
        frame1.axes.get_yaxis().set_visible(False)



In [375]:

    
pylab.rcParams['figure.figsize'] = (18, 10)
h=0.5
i=0

for j in range(S[3]):
    #if Good_ICs[j]:
    if Translated[j]==12 or Translated[j]==13:

        #plt.plot(Time_fluoICA.T,DT[:,Order[j]]+2)
        plt.plot(Time_fluoICA.T,(DT[:,Order[j]]+h*i+2),color=C[j,:])
        plt.plot(Tf,straightintbin*1.5)   
        plt.plot(Tf,leftintbin*1.5)
        plt.plot(Tf,rightintbin*1.5)  
        plt.plot(Tf,FrontGroomint.T)
        plt.plot(Tf,RearGroomint) 
        i=i+1
#plt.xlim([10,17])
plt.show()

Compare significant components for turning left, turning right, and going straight



In [327]:

    
from sklearn import linear_model
import scipy

Extend linearregression class to get p values for the betas (from http://stackoverflow.com/questions/27928275/find-p-value-significance-in-scikit-learn-linearregression)



In [328]:

    
class LinearRegression(linear_model.LinearRegression):

    def __init__(self,*args,**kwargs):
        # *args is the list of arguments that might go into the LinearRegression object
        # that we don't know about and don't want to have to deal with. Similarly, **kwargs
        # is a dictionary of key words and values that might also need to go into the orginal
        # LinearRegression object. We put *args and **kwargs so that we don't have to look
        # these up and write them down explicitly here. Nice and easy.

        if not "fit_intercept" in kwargs:
            kwargs['fit_intercept'] = False

        super(LinearRegression,self).__init__(*args,**kwargs)

    # Adding in t-statistics for the coefficients.
    def fit(self,x,y):
        # This takes in numpy arrays (not matrices). Also assumes you are leaving out the column
        # of constants.

        # Not totally sure what 'super' does here and why you redefine self...
        self = super(LinearRegression, self).fit(x,y)
        n, k = x.shape
        yHat = np.matrix(self.predict(x)).T

        # Change X and Y into numpy matricies. x also has a column of ones added to it.
        x = np.hstack((np.ones((n,1)),np.matrix(x)))
        y = np.matrix(y).T

        # Degrees of freedom.
        df = float(n-k-1)

        # Sample variance.     
        sse = np.sum(np.square(yHat - y),axis=0)
        self.sampleVariance = sse/df

        # Sample variance for x.
        self.sampleVarianceX = x.T*x

        # Covariance Matrix = [(s^2)(X'X)^-1]^0.5. (sqrtm = matrix square root.  ugly)
        self.covarianceMatrix = scipy.linalg.sqrtm(self.sampleVariance[0,0]*self.sampleVarianceX.I)

        # Standard erros for the difference coefficients: the diagonal elements of the covariance matrix.
        self.se = self.covarianceMatrix.diagonal()[1:]

        # T statistic for each beta.
        self.betasTStat = np.zeros(len(self.se))
        for i in xrange(len(self.se)):
        #    self.betasTStat[i] = self.coef_[0,i]/self.se[i]  sophie: I get an error here so change the line to the one below:
             self.betasTStat[i] = self.coef_[i]/self.se[i]

        # P-value for each beta. This is a two sided t-test, since the betas can be 
        # positive or negative.
        self.betasPValue = 1 - scipy.stats.t.cdf(abs(self.betasTStat),df)



In [329]:

    
clf2=LinearRegression()



In [330]:

    
DT.shape









    Out[330]:





(5426, 147)



In [331]:

    
X=np.squeeze(np.array([np.squeeze(straightint),np.squeeze(rightint),np.squeeze(leftint)]))
Nreg=X.shape[0]
Sig=np.zeros([S[3],Nreg])
beta=np.zeros([S[3],Nreg])

for j in range(S[3]):
    if Good_ICs[j]:
        Out=DT[range(len(Tf)),j]    
        clf2.fit(X.T,Out)
        if clf2.betasPValue[2]<(0.005/(len(Good_ICs)*3))and clf2.coef_[2]>0:
            Sig[j,2]=1
            beta[j,2]=clf2.betasPValue[2]
        if clf2.betasPValue[1]<(0.005/(len(Good_ICs)*3))and clf2.coef_[1]>0:
            Sig[j,1]=1 
            beta[j,1]=clf2.betasPValue[1]
        if clf2.betasPValue[0]<(0.005/(len(Good_ICs)*3))and clf2.coef_[0]>0:
            Sig[j,0]=1
            beta[j,0]=clf2.betasPValue[0]



In [332]:

    
Nreg









    Out[332]:





3



In [333]:

    
plt.plot(Sig)









    Out[333]:





[<matplotlib.lines.Line2D at 0x65a7a50>,
 <matplotlib.lines.Line2D at 0x65a78d0>,
 <matplotlib.lines.Line2D at 0x65a74d0>]



In [334]:

    
del Final_map
if S[2]>5:
    Final_map=np.zeros([S[0],S[1],5,3])
    Fmaps=np.zeros([S[0],S[1],5,3])
else:
    Final_map=npLabel_ICs[Neworder[i]] 
    Fmaps=np.zeros([S[0],S[1],3])    
C=np.zeros([S[3],3])



In [335]:

    
for j in range(S[3]):
    if ((Sig[j,2] or Sig[j,1] or Sig[j,0]) and Good_ICs[j]):
        C[j,:]=Sig[j,:]
        for k in range(3):
            M=np.max(np.squeeze(np.reshape(D2[:,:,:,j],S[0]*S[1]*5)))
            #M[M==0]=1
  #          Fmaps[:,:,:,k]=0.7*Dmaps[:,:,:,j]*C[j,k]/np.max(np.squeeze(np.reshape(Dmaps[:,:,:,j]*np.max(C[j,:],S[0]*S[1]*5))
            Fmaps[:,:,:,k]=0.7*D2[:,:,:,j]*C[j,k]/M
        Final_map=Final_map+Fmaps



In [336]:

    
pylab.rcParams['figure.figsize'] = (14, 5)

#Final_map[Final_map<0.1]=np.NaN

if S[2]>5:
    N=Nstack
else:
    N=S[2]
for i in range(N):
        plt.subplot(1,N,i+1)
        #plt.imshow(Dmean[:,:,i],cmap=plt.cm.gray)
        plt.imshow(Final_map[:,:,i],interpolation='none')
        #plt.imshow(Final_map[:,:,i]) 
        frame1 = plt.gca()
        frame1.axes.get_xaxis().set_visible(False)
        frame1.axes.get_yaxis().set_visible(False)



In [337]:

    
pylab.rcParams['figure.figsize'] = (7, 10)
h=1
i=0

for j in range(S[3]):
    if (Sig[j,2] or Sig[j,1] or Sig[j,0]) and Good_ICs[j]:
        plt.plot(Tball,straightbin*1.5)   
        plt.plot(Tball,leftbin*1.5)
        plt.plot(Tball,rightbin*1.5)  
        plt.plot(Tb,FrontGroom.T)
        plt.plot(Tb,RearGroom.T) 
        #plt.plot(Time_fluoICA.T,DT[:,Order[j]]+2)
        plt.plot(Time_fluoICA.T,(DT[:,Order[j]]+h*i+2),color=C[j,:]/2)
        i=i+1

plt.show()









    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-337-569f8aecaf00> in <module>()
      5 for j in range(S[3]):
      6     if (Sig[j,2] or Sig[j,1] or Sig[j,0]) and Good_ICs[j]:
----> 7         plt.plot(Tball,straightbin*1.5)
      8         plt.plot(Tball,leftbin*1.5)
      9         plt.plot(Tball,rightbin*1.5)

NameError: name 'straightbin' is not defined

Compare significant components for walking, grooming, and resting



In [338]:

    
Groom=FrontGroomint+np.squeeze(RearGroomint)
Groom=np.ceil(Groom/2)



In [339]:

    
Walk=np.squeeze(straightint)+np.squeeze(rightint)
Walk=np.ceil(Walk/2)
Walk=Walk+leftint
Walk=np.ceil(Walk/2)



In [340]:

    
len(Good_ICs)









    Out[340]:





147



In [341]:

    
X=np.squeeze(np.array([np.squeeze(Walk),np.squeeze(Groom),np.squeeze(Restint)]))
Nreg=X.shape[0]
Sig=np.zeros([S[3],Nreg])

for j in range(S[3]):
    if Good_ICs[j]==1:
        Out=DT[range(len(Tf)),j]    
        clf2.fit(X.T,Out)
        if clf2.betasPValue[2]<(0.005/(len(Good_ICs)*3)) and clf2.coef_[2]>0:
            Sig[j,2]=1
        if clf2.betasPValue[1]<(0.005/(len(Good_ICs)*3))and clf2.coef_[1]>0:
            Sig[j,1]=1  
        if clf2.betasPValue[0]<(0.005/(len(Good_ICs)*3))and clf2.coef_[0]>0:
            Sig[j,0]=1



In [342]:

    
plt.plot(Sig)









    Out[342]:





[<matplotlib.lines.Line2D at 0xc6b5690>,
 <matplotlib.lines.Line2D at 0xc6b5910>,
 <matplotlib.lines.Line2D at 0xc6b5b50>]



In [343]:

    
del Final_map
if S[2]>5:
    Final_map=np.zeros([S[0],S[1],5,3])
    Fmaps=np.zeros([S[0],S[1],5,3])
else:
    Final_map=npLabel_ICs[Neworder[i]] 
    Fmaps=np.zeros([S[0],S[1],3])    
C=np.zeros([S[3],3])



In [344]:

    
for j in range(S[3]):    
    if (Sig[j,2] or Sig[j,1] or Sig[j,0]) and Good_ICs[j]:
#        C[j,:]=np.squeeze(np.random.rand(3,1))
        C[j,:]=Sig[j,:]
        C[j,2]=0
        for k in range(3):
            M=np.max(np.squeeze(np.reshape(D2[:,:,:,j],S[0]*S[1]*5)))
            #M[M==0]=1
  #          Fmaps[:,:,:,k]=0.7*Dmaps[:,:,:,j]*C[j,k]/np.max(np.squeeze(np.reshape(Dmaps[:,:,:,j]*np.max(C[j,:],S[0]*S[1]*5))
            Fmaps[:,:,:,k]=0.7*D2[:,:,:,j]*C[j,k]/M
        Final_map=Final_map+Fmaps



In [345]:

    
pylab.rcParams['figure.figsize'] = (14, 5)

#Final_map[Final_map<0.1]=np.NaN

if S[2]>5:
    N=Nstack
else:
    N=S[2]
for i in range(N):
        plt.subplot(1,N,i+1)
        plt.imshow(Dmean[:,:,i],cmap=plt.cm.gray)
        plt.imshow(Final_map[:,:,i],interpolation='none')
        #plt.imshow(Final_map[:,:,i]) 
        frame1 = plt.gca()
        frame1.axes.get_xaxis().set_visible(False)
        frame1.axes.get_yaxis().set_visible(False)



In [346]:

    
pylab.rcParams['figure.figsize'] = (7, 10)
h=1
i=0

for j in range(S[3]):
    if (Sig[j,2] or Sig[j,1] or Sig[j,0]) and Good_ICs[j]:
        plt.plot(Tball,straightbin*1.5)   
        plt.plot(Tball,leftbin*1.5)
        plt.plot(Tball,rightbin*1.5)  
        plt.plot(Tb,FrontGroom.T)
        plt.plot(Tb,RearGroom.T) 
        #plt.plot(Time_fluoICA.T,DT[:,Order[j]]+2)
        plt.plot(Time_fluoICA.T,(DT[:,Order[j]]+h*i+2),color=C[j,:]/2)
        i=i+1

plt.show()









    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-346-569f8aecaf00> in <module>()
      5 for j in range(S[3]):
      6     if (Sig[j,2] or Sig[j,1] or Sig[j,0]) and Good_ICs[j]:
----> 7         plt.plot(Tball,straightbin*1.5)
      8         plt.plot(Tball,leftbin*1.5)
      9         plt.plot(Tball,rightbin*1.5)

NameError: name 'straightbin' is not defined

Compare significant components for window just before walking and just before stoping



In [347]:

    
plt.plot(range(600,800),Walk[range(600,800)])









    Out[347]:





[<matplotlib.lines.Line2D at 0xc1f0590>]

Detect onsets and offsets



In [352]:

    
OnSteps=[]
OffSteps=[]
PreStart=np.zeros(Walk.shape)
PreStop=np.zeros(Walk.shape)

for j in range(20,len(Walk)-20):
    if (not(any(Walk[range(j-10,j)])) and all(Walk[range(j,j+10)])):
        OnSteps.append(j)
        PreStart[range(j-20,j)]=1
    elif (all(Walk[range(j-10,j)]) and not(any(Walk[range(j,j+10)]))):
        OffSteps.append(j)
        PreStop[range(j-20,j)]=1



In [353]:

    
plt.plot(PreStart)
plt.plot(PreStop)
plt.plot(Walk)
plt.show()



In [354]:

    
X=np.squeeze(np.array([np.squeeze(PreStart),np.squeeze(PreStop)]))
Nreg=X.shape[0]
Sig=np.zeros([S[3],Nreg])

for j in range(S[3]):
    if Good_ICs[j]==1:
        Out=DT[range(len(Tf)),j]    
        clf2.fit(X.T,Out)
        if clf2.betasPValue[1]<(0.0005/(len(Good_ICs)*2))and clf2.coef_[1]>0:
            Sig[j,1]=1  
        if clf2.betasPValue[0]<(0.0005/(len(Good_ICs)*2))and clf2.coef_[0]>0:
            Sig[j,0]=1



In [355]:

    
del Final_map
if S[2]>5:
    Final_map=np.zeros([S[0],S[1],5,3])
    Fmaps=np.zeros([S[0],S[1],5,3])
else:
    Final_map=npLabel_ICs[Neworder[i]] 
    Fmaps=np.zeros([S[0],S[1],3])    
C=np.zeros([S[3],3])



In [356]:

    
for j in range(S[3]):    
    if (Sig[j,1] or Sig[j,0]) and Good_ICs[j]:
#        C[j,:]=np.squeeze(np.random.rand(3,1))
        C[j,range(2)]=Sig[j,:]
        for k in range(3):
            M=np.max(np.squeeze(np.reshape(D2[:,:,:,j],S[0]*S[1]*5)))
            #M[M==0]=1
  #          Fmaps[:,:,:,k]=0.7*Dmaps[:,:,:,j]*C[j,k]/np.max(np.squeeze(np.reshape(Dmaps[:,:,:,j]*np.max(C[j,:],S[0]*S[1]*5))
            Fmaps[:,:,:,k]=0.7*D2[:,:,:,j]*C[j,k]/M
        Final_map=Final_map+Fmaps



In [357]:

    
pylab.rcParams['figure.figsize'] = (14, 5)

#Final_map[Final_map<0.1]=np.NaN

if S[2]>5:
    N=Nstack
else:
    N=S[2]
for i in range(N):
        plt.subplot(1,N,i+1)
        plt.imshow(Dmean[:,:,i],cmap=plt.cm.gray)
        plt.imshow(Final_map[:,:,i],interpolation='none')
        #plt.imshow(Final_map[:,:,i]) 
        frame1 = plt.gca()
        frame1.axes.get_xaxis().set_visible(False)
        frame1.axes.get_yaxis().set_visible(False)

Open data and plot average in each category



In [380]:

    
# from http://stackoverflow.com/questions/3579568/choosing-a-file-in-python-with-simple-dialog
from Tkinter import Tk
from tkFileDialog import askopenfilename

Tk().withdraw() # we don't want a full GUI, so keep the root window from appearing
filename2 = askopenfilename() # show an "Open" dialog box and return the path to the selected file
print(filename2)









    



/home/sophie/Downloads/100106ss2on500regcregcU10sMpsfkf.nii



In [389]:

    
img1 = nb.load(filename2)
InitData = np.array(img1.get_data())
S=InitData.shape
S









    Out[389]:





(185, 116, 10, 5426)



In [390]:

    
InitDatap=np.zeros(S)
for i in range(S[0]):
    for j in range(S[1]):
        for k in range(S[2]):
            InitDatap[i,j,k,:]=InitData[i,j,k,:]-np.min(InitData[i,j,k,:])



In [393]:

    
StraightAv=np.zeros([S[0],S[1],S[2]])
for i in range(S[3]-200):
    StraightAv=StraightAv+InitDatap[:,:,:,i]*straightintbin[i]
LeftAv=np.zeros([S[0],S[1],S[2]])*len(OffSteps)
for i in range(S[3]-200):
    LeftAv=LeftAv+InitDatap[:,:,:,i]*leftintbin[i]
rightAv=np.zeros([S[0],S[1],S[2]])
for i in range(S[3]-200):
    rightAv=StraightAv+InitDatap[:,:,:,i]*rightintbin[i]
SLR=np.zeros([S[0],S[1],S[2],3])
SLR[:,:,:,0]=StraightAv/(S[3]-2)
SLR[:,:,:,1]=LeftAv/(S[3]-2)
SLR[:,:,:,2]=rightAv/(S[3]-2)



In [394]:

    
S2=StraightAv.shape



In [395]:

    
N=S2[2]/2
for i in range(N):
    plt.subplot(1,N,i+1)
    plt.imshow(1+SLR[:,:,i*2,:],interpolation='none')
    #plt.imshow(Final_map[:,:,i]) 
    frame1 = plt.gca()
    frame1.axes.get_xaxis().set_visible(False)
    frame1.axes.get_yaxis().set_visible(False)



In [397]:

    
WalkAv=np.mean(SLR,3)*3
GroomAv=np.zeros([S[0],S[1],S[2]])
for i in range(S[3]-200):
    GroomAv=GroomAv+InitDatap[:,:,:,i]*(FrontGroomint[i]+RearGroomint[i])

WG=np.zeros([S[0],S[1],S[2],3])
WG[:,:,:,0]=WalkAv
WG[:,:,:,1]=GroomAv/(S[3]-2)
WG[:,:,:,2]=0



In [399]:

    
N=S2[2]/2
for i in range(N):
    plt.subplot(1,N,i+1)
    plt.imshow(WG[:,:,i*2,:],interpolation='none')
    #plt.imshow(Final_map[:,:,i]) 
    frame1 = plt.gca()
    frame1.axes.get_xaxis().set_visible(False)
    frame1.axes.get_yaxis().set_visible(False)

Cut data around transitions



In [400]:

    
OnSteps









    Out[400]:





[854, 1625, 1714, 2343, 4113, 4148]



In [401]:

    
OffSteps









    Out[401]:





[398, 1733, 4123]



In [402]:

    
len(OnSteps)









    Out[402]:





6



In [403]:

    
DataOnset=np.zeros([S[0],S[1],S[2],100,len(OnSteps)])
for i in range(len(OnSteps)):
    j=OnSteps[i]
    DataOnset[:,:,:,:,i]=InitData[:,:,:,range(j-50,j+50)]



In [404]:

    
DataOffset=np.zeros([S[0],S[1],S[2],100,len(OffSteps)])
for i in range(len(OffSteps)):
    j=OffSteps[i]
    DataOffset[:,:,:,:,i]=InitData[:,:,:,range(j-50,j+50)]



In [405]:

    
DataOnsetAv=np.mean(DataOnset,4)
DataOffsetAv=np.mean(DataOffset,4)



In [413]:

    
D3=np.transpose(DataOffsetAv[:,:,:,:],(3,2,1,0))
nim = NiftiImage(D3)
nim.save('/home/sophie/100106Off.nii')



In [414]:

    
D3=np.transpose(DataOnsetAv[:,:,:,:],(3,2,1,0))
nim = NiftiImage(D3)
nim.save('/home/sophie/100106On.nii')



In [415]:

    
S









    Out[415]:





(185, 116, 10, 5426)



In [419]:

    
DataOnsetcat=np.zeros([S[0],S[1]*len(OnSteps),S[2],100])
for i in range(len(OnSteps)):
    j=OnSteps[i]
    DataOnsetcat[:,range(S[1]*i,S[1]*(i+1)),:,:]=InitData[:,:,:,range(j-50,j+50)]
    
DataOffsetcat=np.zeros([S[0],S[1]*len(OnSteps),S[2],100])
for i in range(len(OffSteps)):
    j=OffSteps[i]
    DataOffsetcat[:,range(S[1]*i,S[1]*(i+1)),:,:]=InitData[:,:,:,range(j-50,j+50)]



In [420]:

    
D3=np.transpose(DataOnsetcat[:,:,:,:],(3,2,1,0))
nim = NiftiImage(D3)
nim.save('/home/sophie/100106Oncat.nii')
D3=np.transpose(DataOffsetcat[:,:,:,:],(3,2,1,0))
nim = NiftiImage(D3)
nim.save('/home/sophie/100106Offcat.nii')



In [270]:

    
List1=[(Translated[i],i) for i in range(S[3])]
Newlist=sorted(List1, key=lambda List1: List1[0])



In [271]:

    
Neworder=[Newlist[i][1] for i in range(S[3])  if Newlist[i][0] != 15]



In [272]:

    
[Label_ICs[Neworder[i]] for i in range(len(Neworder))]









    Out[272]:





['A',
 'A',
 'A',
 'A',
 'a',
 'A',
 'A',
 'A',
 'a',
 'a',
 'A',
 'A',
 'o',
 'o',
 'o',
 'o',
 'o',
 'o',
 'o',
 'o',
 'o',
 'O',
 'o',
 'o',
 'o',
 'o',
 'o',
 'o',
 'o',
 'OG',
 'PN',
 'PN',
 'LH',
 'LH',
 'LH',
 'LH',
 'LH',
 'LH',
 'LH',
 'SLP',
 'SMP',
 'SMP',
 'SMP',
 'SMP',
 'SLP',
 'SMP',
 'SMP',
 'SLP',
 'SMP',
 'SMP',
 'KC',
 'PB',
 'PB',
 'PB',
 'PB',
 'PB',
 'PB',
 'PB',
 'FB',
 'PB',
 'FB',
 'FB',
 'EB',
 'AMMC',
 'AMMC',
 'AMMC',
 'AMMC',
 'AMMC',
 'AMMC',
 'SOG',
 'PS',
 'PS',
 'PS',
 'PS',
 'PS',
 'DCB']



In [273]:

    
NewDT=DT[:,Order[Neworder[:]]].T
for j in range(len(Neworder)):
    A=NewDT[:,j]
    V=np.sqrt(np.var(A))
    NewDT[:,j]=A/V



In [274]:

    
C1=np.zeros([6,3])
C1[0][:]=(1,0,0)
C1[1][:]=(0,1,0)
C1[2][:]=(0,0,1)
C1[3][:]=(0.8,0.8,0)
C1[4][:]=(0,1,1)
C1[5][:]=(1,0,1)
S1=DT.shape



In [275]:

    
h=0.5
pylab.rcParams['figure.figsize'] = (7, 14)
i=0
t=[j*0.005 for j in range(S1[0])]

for j in range(len(Neworder)):
        plt.plot(t,NewDT[j,:]+h*j,color=C1[i%6][:])
        i=i+1
plt.ylim([-3,h*j+5])

frame1=plt.gca()
frame1.axes.get_yaxis().set_ticks([])
matplotlib.rcParams.update({'font.size': 18})
plt.show()



In [276]:

    
Newmaps=Dmaps[:,:,:,Order[Neworder[:]]]



In [277]:

    
L=len(set([Translated[Neworder[i]] for i in range(len(Neworder))]))
Regionmaps=np.zeros([S[0],S[1],L,3])
Datasort=np.zeros([S[0],S[1],S[2],L,3])


Mapsordered=datao[:,:,:,Neworder[:]]
mapfilename=foldername+'mapsordered.nii'
Mapsordered2=np.transpose(Mapsordered,(3,2,1,0))
#nimap = NiftiImage(Mapsordered2)
#nimap.save(mapfilename)

DMapsordered=Dmapso[:,:,:,Neworder[:]]
mapfilename=foldername+'zscoredmapsordered.nii'
DMapsordered2=np.transpose(DMapsordered,(3,2,1,0))
#nimap = NiftiImage(DMapsordered2)
#nimap.save(mapfilename)



In [278]:

    
j=0
i=0
k=Translated[Neworder[i]]
m=0
Regionname.append(Label_ICs[Neworder[i]])
for i in range(len(Neworder)):
    #C2=C1[i%6][:]
    for l in range(3):
        M=np.max(np.squeeze(np.reshape(Newmaps[:,:,:,i],S[0]*S[1]*S[2])))
        Regionmaps[:,:,j,l]=Regionmaps[:,:,j,l]+0.7*np.max(DMapsordered[:,:,:,i],2)*C1[i%6][l]/M
        Datasort[:,:,:,j,l]=Datasort[:,:,:,j,l]+Dmaps[:,:,:,Order[Neworder[i]]]*C1[i%6][l] 
    i=i+1
    m=m+1
    if i<len(Neworder):
        k1=Translated[Neworder[i]]
        
    if k1 != k:
        j=j+1
        k=k1
        m=0
        Regionname.append(Label_ICs[Neworder[i]])

Datasort.shape









    Out[278]:





(178, 102, 34, 12, 3)



In [279]:

    
Regionname









    Out[279]:





['AL',
 'AL',
 'o',
 'PN',
 'LH',
 'D',
 'MBON',
 'FB',
 'EB',
 'PB',
 'AMMC',
 'PS',
 '',
 'AL',
 'o',
 'PN',
 'LH',
 'D',
 'MBON',
 'FB',
 'EB',
 'PB',
 'AMMC',
 'PS',
 'DCB',
 'AL',
 'o',
 'PN',
 'LH',
 'D',
 'MBON',
 'FB',
 'EB',
 'PB',
 'AMMC',
 'PS',
 'DCB',
 'AL',
 'o',
 'PN',
 'LH',
 'D',
 'MBON',
 'FB',
 'EB',
 'PB',
 'AMMC',
 'PS',
 'DCB',
 'A',
 'o',
 'OG',
 'PN',
 'LH',
 'SLP',
 'KC',
 'PB',
 'EB',
 'AMMC',
 'SOG',
 'DCB']



In [281]:

    
DMapscolor=np.sum(Datasort, axis=3)
DMapscolor.shape
mapfilename=foldername+'zscoredmapscolor.nii'
DMapscolor2=np.transpose(DMapscolor,(3,2,1,0))
DMapscolor2.shape
nimap = NiftiImage(DMapscolor2)
nimap.save(mapfilename)
Dmeanproj=np.max(Dmean,2)
Regionmaps[Regionmaps==0]=np.nan



In [282]:

    
pylab.rcParams['figure.figsize'] = (7, 14)
import scipy
from scipy import ndimage
#from skimage import color

Dmeanprojrgb=np.dstack((Dmeanproj.T,Dmeanproj.T,Dmeanproj.T))



In [284]:

    
j=0
i=0
k=Translated[Neworder[i]]
m=0
Regionname.append(Label_ICs[Neworder[i]])
for i in range(len(Neworder)):
    #C2=C1[i%6][:]
    for l in range(3):
        M=np.max(np.squeeze(np.reshape(Newmaps[:,:,:,i],S[0]*S[1]*S[2])))
        Regionmaps[:,:,j,l]=Regionmaps[:,:,j,l]+0.7*np.max(DMapsordered[:,:,:,i],2)*C1[i%6][l]/M
        Datasort[:,:,:,j,l]=Datasort[:,:,:,j,l]+Dmaps[:,:,:,Order[Neworder[i]]]*C1[i%6][l] 
    i=i+1
    m=m+1
    if i<len(Neworder):
        k1=Translated[Neworder[i]]
        
    if k1 != k:
        j=j+1
        k=k1
        m=0
        Regionname.append(Label_ICs[Neworder[i]])



In [286]:

    
Regionmaps=np.zeros([S[0],S[1],L,3])
Datasort=np.zeros([S[0],S[1],S[2],L,3])



In [287]:

    
j=0
i=0
k=Translated[Neworder[i]]
m=0
Regionname.append(Label_ICs[Neworder[i]])
for i in range(len(Neworder)):
    #C2=C1[i%6][:]
    for l in range(3):
        M=np.max(np.squeeze(np.reshape(Newmaps[:,:,:,i],S[0]*S[1]*S[2])))
        Regionmaps[:,:,j,l]=Regionmaps[:,:,j,l]+0.7*np.max(DMapsordered[:,:,:,i],2)*C1[i%6][l]/M
        Datasort[:,:,:,j,l]=Datasort[:,:,:,j,l]+Dmaps[:,:,:,Order[Neworder[i]]]*C1[i%6][l] 
    i=i+1
    m=m+1
    if i<len(Neworder):
        k1=Translated[Neworder[i]]
        
    if k1 != k:
        j=j+1
        k=k1
        m=0
        Regionname.append(Label_ICs[Neworder[i]])



In [288]:

    
pylab.rcParams['figure.figsize'] = (7, 14)
import scipy
from scipy import ndimage

for i in range(L):
        plt.subplot(L,1,L-i)
        Rotated_Plot = ndimage.rotate(Regionmaps[:,:,i], -90)
        IM=plt.imshow(Rotated_Plot) 
        frame1 = plt.gca()
        frame1.axes.get_xaxis().set_visible(False)
        frame1.axes.get_yaxis().set_visible(False)



In [156]:

    
import scipy.io as sio
sio.savemat(foldername+'NewDT.mat',{'NewDT':NewDT})



In [157]:

    
names=set([Label_ICs[Neworder[i]] for i in range(len(Neworder))])  
for i in range(L):
    regionfilename=foldername+Regionname[i]+'.nii'
    D3=np.transpose(Datasort[:,:,:,i,:],(2,3,1,0))
    nim = NiftiImage(D3)
    nim.save(regionfilename)



In [ ]: