<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 jPCA
Churchland, Mark M., et al. "Neural population dynamics during reaching." Nature 487.7405 (2012)
Step1:  The data can be loaded with pandas 
Step2:  It's the responses of theta-modulated thalamic neurons to hippocampal sharp-waves ripples 
Step3:  There are 767 neurons here with 201 time bins. The order of the data matrix is 
Step4:  First step is the classical PCA to reduce the dimensionality of the dataset.
Similar to the original article, I reduced it to 6 dimensions
Step5:  We can thus work on the 6 first components of the PCA
Step6:  We can plot the 6 components 
Step7:  Now we can compute $\dot{X}$ using the function written below 
Step8:  The function derivative is called for each component 
Step9:  Next step is to build the H mapping using this function 
Step10:  $\tilde{X}$ is the block diagonal matrix
Step11:  We can put $\dot{X}$ in columns 
Step12:  Multiply $\tilde{X}$ by $H$ 
Step13:  and solve $(\tilde{X}.H).k = \dot{X}$
Step14:  Do $m = H.k$ to get $M_{skew}$
Step15:  Construct the two vectors for projection with $M_{skew}$
Step16:  and get the jpc vectors as $X_r = X.u$
Step17:  We can now look at the two jpc components 
Step18:  We can now project the data on rX to find the swr angle 
Step19:  We can now represent the sharp-waves phase for all neurons as 
<ASSISTANT_TASK:>
Python Code:
import numpy as np
from sklearn.decomposition import PCA
import pandas as pd
from pylab import *
Explanation: jPCA
Churchland, Mark M., et al. "Neural population dynamics during reaching." Nature 487.7405 (2012): 51.
End of explanation
data = pd.read_hdf("swr_modth.h5")
Explanation: The data can be loaded with pandas :
End of explanation
figure()
plot(data)
xlabel("Time lag (ms)")
ylabel("Modulation (z-scored)")
show()
Explanation: It's the responses of theta-modulated thalamic neurons to hippocampal sharp-waves ripples :
End of explanation
print(data.shape)
Explanation: There are 767 neurons here with 201 time bins. The order of the data matrix is :
End of explanation
n = 6
pca = PCA(n_components = n)
new_data = pca.fit_transform(data.values.T) # data needs to be inverted here depending of how you do the PCA
Explanation: First step is the classical PCA to reduce the dimensionality of the dataset.
Similar to the original article, I reduced it to 6 dimensions
End of explanation
X = pca.components_.transpose()
times = data.index.values
Explanation: We can thus work on the 6 first components of the PCA
End of explanation
figure()
plot(times, X)
xlabel("Time lag (ms)")
ylabel("pc")
show()
Explanation: We can plot the 6 components :
End of explanation
def derivative(x, f):
    ''' 
        Compute the derivative of a time serie
        Used for jPCA
    '''
    from scipy.stats import linregress
    fish = np.zeros(len(f))
    slopes_ = []
    tmpf = np.hstack((f[0],f,f[-1])) # not circular
    binsize = x[1]-x[0]	
    tmpx = np.hstack((np.array([x[0]-binsize]),x,np.array([x[-1]+binsize])))	
    # plot(tmpx, tmpf, 'o')
    # plot(x, f, '+')
    for i in range(len(f)):
        slope, intercept, r_value, p_value, std_err = linregress(tmpx[i:i+3], tmpf[i:i+3])
        slopes_.append(slope)	
        # plot(tmpx[i:i+3], tmpx[i:i+3]*slope+intercept, '-')
    return np.array(slopes_)/binsize
Explanation: Now we can compute $\dot{X}$ using the function written below :
End of explanation
dX = np.zeros_like(X)
for i in range(n):
    dX[:,i] = derivative(times, X[:,i])
Explanation: The function derivative is called for each component :
End of explanation
def buildHMap(n, ):
    '''
        build the H mapping for a given n
        used for the jPCA
    '''
    from scipy.sparse import lil_matrix
    M = np.zeros((n,n), dtype = np.int)
    M[np.triu_indices(n,1)] = np.arange(1,int(n*(n-1)/2)+1)
    M = M - M.transpose()
    m = np.vstack(M.reshape(n*n))
    k = np.vstack(M[np.triu_indices(n,1)]).astype('int')
    H = lil_matrix( (len(m), len(k)), dtype = np.float16)
    H = np.zeros( (len(m), len(k) ))
    # first column 
    for i in k.flatten():
        # positive
        H[np.where(m == i)[0][0],i-1] = 1.0
        # negative
        H[np.where(m == -i)[0][0],i-1] = -1.0
    return H
H = buildHMap(n)
Explanation: Next step is to build the H mapping using this function :
End of explanation
Xtilde = np.zeros( (X.shape[0]*X.shape[1], X.shape[1]*X.shape[1]) )
for i, j in zip( (np.arange(0,n**2,n) ), np.arange(0, n*X.shape[0], X.shape[0]) ):
    Xtilde[j:j+X.shape[0],i:i+X.shape[1]] = X
Explanation: $\tilde{X}$ is the block diagonal matrix:
End of explanation
dXv = np.vstack(dX.transpose().reshape(X.shape[0]*X.shape[1]))
Explanation: We can put $\dot{X}$ in columns :
End of explanation
XtH = np.dot(Xtilde, H)
Explanation: Multiply $\tilde{X}$ by $H$ :
End of explanation
k, residuals, rank, s = np.linalg.lstsq(XtH, dXv, rcond = None)
Explanation: and solve $(\tilde{X}.H).k = \dot{X}$
End of explanation
m = np.dot(H, k)
Mskew = m.reshape(n,n).transpose()
Explanation: Do $m = H.k$ to get $M_{skew}$:
End of explanation
evalues, evectors = np.linalg.eig(Mskew)
index = np.argsort(np.array([np.linalg.norm(i) for i in evalues]).reshape(int(n/2),2)[:,0])
evectors = evectors.transpose().reshape(int(n/2),2,n)
u = np.vstack([np.real(evectors[index[-1]][0] + evectors[index[-1]][1]),
                np.imag(evectors[index[-1]][0] - evectors[index[-1]][1])]).transpose()
Explanation: Construct the two vectors for projection with $M_{skew}$:
End of explanation
rX = np.dot(X, u)
Explanation: and get the jpc vectors as $X_r = X.u$
End of explanation
figure(figsize=(15, 5))
subplot(121)
plot(times, rX)
xlabel("Time lag (ms)")
subplot(122)
plot(rX[:,0], rX[:,1])
show()
Explanation: We can now look at the two jpc components :
End of explanation
score = np.dot(data.values.T, rX)
phi = np.mod(np.arctan2(score[:,1], score[:,0]), 2*np.pi)
Explanation: We can now project the data on rX to find the swr angle :
End of explanation
figure(figsize = (10,10))
scatter(score[:,0], score[:,1], c = phi)
scatter(np.cos(phi)*np.max(score), np.sin(phi)*np.max(score), c = phi)
show()
Explanation: We can now represent the sharp-waves phase for all neurons as :
End of explanation
<END_TASK>