notebook.community

Edit and run



In [1]:

    
clear all



In [2]:

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

    
# 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/100106ss2on250regcregcmedU10sMpsfkf250Smith0_4TS.mat



In [104]:

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









    Out[104]:





(5676, 250)



In [105]:

    
# 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/100106ss2on250regcregcmedU10sMpsfkf250Smith0_4IC.nii



In [106]:

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









    Out[106]:





(185, 116, 10, 250)



In [107]:

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

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





(1, 5927)



In [109]:

    
Time_fluoICA=Time_fluo[:,250:5926]



In [110]:

    
Time_fluoICA.shape









    Out[110]:





(1, 5676)



In [111]:

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

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

    
B









    Out[18]:





{'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 [21]:

    
b=B['behavior']



In [22]:

    
b.shape









    Out[22]:





(7, 11593)



In [23]:

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



In [24]:

    
Weird.shape









    Out[24]:





(11593,)



In [433]:

    
#Threshold ball movements



In [25]:

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



In [26]:

    
plt.plot(left)









    Out[26]:





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



In [206]:

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



In [30]:

    
Tb=np.squeeze(Tbehavior)



In [29]:

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



In [212]:

    
#interpolate for ball movements



In [33]:

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



In [53]:

    
max(Tb)









    Out[53]:





117.45699999999998



In [57]:

    
max(Time_fluoICA.T)









    Out[57]:





array([ 118.66205987])



In [40]:

    
plt.plot(Tb,'g')
plt.plot(Time_fluoICA.T,'r')









    Out[40]:





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



In [97]:

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



In [72]:

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

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

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

    
#Align times (verify what frames have been removed)



In [79]:

    
Tb.shape









    Out[79]:





(11593,)



In [80]:

    
rightintbin.shape









    Out[80]:





(5600, 1)



In [81]:

    
plt.plot(Tf,rightintbin)









    Out[81]:





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



In [82]:

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









    Out[82]:





117.46299999999998

Process the maps



In [112]:

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

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



In [114]:

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

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



In [116]:

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

    
Dtemp=np.mean(data,3)



In [118]:

    
Dtemp.shape









    Out[118]:





(185, 116, 10)



In [119]:

    
%%javascript
IPython.OutputArea.auto_scroll_threshold =4000;



In [120]:

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

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









    Out[121]:





<matplotlib.image.AxesImage at 0x15ba2510>



In [122]:

    
Tb.shape









    Out[122]:





(11593,)



In [123]:

    
FrontGroom.shape









    Out[123]:





(11593,)



In [ ]:

    
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 
            straightbin

    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.T) 
    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 [67]:

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

    
Dmaps.shape









    Out[342]:





(178, 102, 34, 600)



In [343]:

    
set(Label_ICs)









    Out[343]:





{'',
 'A',
 'AMMC',
 'DCB',
 'EB',
 'FB',
 'KC',
 'LH',
 'MBON',
 'O',
 'OG',
 'PB',
 'PN',
 'PS',
 'SLP',
 'SMP',
 'SOG',
 'a',
 'o'}



In [344]:

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



In [345]:

    
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,'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 [346]:

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



In [347]:

    
G=Good_ICs.tolist();



In [348]:

    
len(Good_ICs)









    Out[348]:





600



In [349]:

    
G.count(1)









    Out[349]:





76



In [350]:

    
where_are_NaNs = np.isnan(D2)

D2[where_are_NaNs] = 0



In [669]:

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

    
D2.shape









    Out[670]:





(178, 102, 5, 600)



In [671]:

    
for j in range(S[3]):    
    if Good_ICs[j]:
        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 [672]:

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

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

for j in range(S[3]):
    if Translated[j]==5:
        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,:])
        i=i+1

plt.show()

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



In [459]:

    
from sklearn import linear_model

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

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

    
clf2=LinearRegression()



In [577]:

    
X.shape









    Out[577]:





(3, 5327)



In [654]:

    
X=np.squeeze(np.array([np.squeeze(straightint),np.squeeze(rightint),np.squeeze(leftint)]))
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.05/(len(Good_ICs)*3)):
            Sig[j,2]=1
        if clf2.betasPValue[1]<(0.05/(len(Good_ICs)*3)):
            Sig[j,1]=1  
        if clf2.betasPValue[0]<(0.05/(len(Good_ICs)*3)):
            Sig[j,0]=1



In [666]:

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

    
for j in range(S[3]):    
    if (Sig[j,2] or Sig[j,1] or Sig[j,0]) and Good_ICs[j]:
#    if (Sig[j,1] or Sig[j,0]) and Good_ICs[j]:
#        C[j,:]=np.squeeze(np.random.rand(3,1))
        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 [665]:

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

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

Compare significant components for walking, grooming, and resting



In [628]:

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



In [640]:

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



In [673]:

    
len(Good_ICs)









    Out[673]:





600



In [641]:

    
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.05/(len(Good_ICs)*3)):
            Sig[j,2]=1
        if clf2.betasPValue[1]<(0.05/(len(Good_ICs)*3)):
            Sig[j,1]=1  
        if clf2.betasPValue[0]<(0.05/(len(Good_ICs)*3)):
            Sig[j,0]=1



In [649]:

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

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

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

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

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



In [694]:

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









    Out[694]:





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

Detect onsets and offsets



In [698]:

    
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-20,j)])) and all(Walk[range(j+1,j+20)])):
        OnSteps.append(j)
        PreStart[range(j-20,j)]=1
    elif (all(Walk[range(j-20,j)]) and not(any(Walk[range(j+1,j+20)]))):
        OffSteps.append(j)
        PreStop[range(j-20,j)]=1



In [702]:

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



In [704]:

    
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.05/(len(Good_ICs)*2)):
            Sig[j,1]=1  
        if clf2.betasPValue[0]<(0.05/(len(Good_ICs)*2)):
            Sig[j,0]=1



In [705]:

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

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

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