In [5]:
from matplotlib.mlab import PCA as mlabPCA
import h5py
from pylab import*
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d import proj3d
import matplotlib.animation as animation
params = { 'figure.figsize' : [10, 10],
'figure.dpi:' : 150,
'figure.autolayout:' : True,
'text.usetex' : True,
'font.size' : 11,
'font.family' : 'lmodern',
'text.latex.unicode': True,
'legend.fontsize' : 'medium',
'legend.handlelength' : 2.2,
'lines.linewidth' : 0.7,
}
plt.rcParams.update(params)
In [6]:
class data:
def __init__(self,animal, n_ignoredDays):
self.animal = animal
self.valueMats = []
self.trialTypeMats = []
self.trialIndicesVecs = []
self.trialTypeVecs = []
self.n_days = len(animal)
self.n_ignoredDays = n_ignoredDays
self.dayIndex = zeros(self.n_days, dtype=numpy.int)
self.extractData()
self.mergedValueMats = self.valueMats[0]
self.mergeValueMats()
self.mapTrialIndexToType()
self.pcaMatrices = []
self.pcaMatrix = []
def mapTrialIndexToType(self):
trialTypes = array([1,2,3,4,5,6], dtype=int)
for i, matrix in enumerate(self.valueMats):
trialTypeVec = dot(self.trialTypeMats[i], trialTypes)
trial = trialTypeVec[array(self.trialIndicesVecs[i]-1, dtype = int)]
self.trialTypeVecs.append(trial)
def extractData(self):
for dayId, day in enumerate(self.animal):
self.valueMats.append(self.animal.get(day + "/valueMatrix"))
self.trialTypeMats.append(self.animal.get(day + "/trialTypeMat"))
self.trialIndicesVecs.append(self.animal.get(day + "/trialIndices"))
self.removeNan()
if dayId < self.n_days - self.n_ignoredDays:
self.dayIndex[dayId+1] = int(self.dayIndex[dayId] + self.valueMats[-1].shape[0])
def removeNan(self):
for index, matrix in enumerate(self.valueMats):
badData = array(unique(where(isnan(matrix))[0]))
matrix = delete(matrix, badData,0)
self.trialIndicesVecs[index] = delete(self.trialIndicesVecs[index], badData,0)
self.valueMats[index] = matrix
def mergeValueMats(self):
for i in range(1, self.n_days-self.n_ignoredDays):
self.mergedValueMats = concatenate((self.mergedValueMats, self.valueMats[i]))
def pcaOnEachDay(self):
for matrix in self.valueMats:
matrix = array(matrix)
self.pcaMatrices.append(mlabPCA(matrix))
def pcaOnAllDays(self):
self.pcaMatrix = mlabPCA(self.mergedValueMats)
In [7]:
def axSettings(ax):
ax.legend()
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')
ax.set_zlabel('PC3')
ax.set_xlim([-40., 40.])
ax.set_ylim([-40., 40.])
ax.set_zlim([-10., 10.])
return ax
def plotAllDays(pcaMatrix, dayIndex,nPoints = 10, n_ignoredDays = 1):
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
day = 0
if dayIndex == None:
for matrix in pcaMatrix:
day +=1
ax.plot(matrix.Y[::nPoints,0], matrix.Y[::nPoints,1],matrix.Y[::nPoints,2],
'o', markersize=7, alpha=0.1, label = 'day: ' + str(day))
ax = axSettings(ax)
else:
for i in range(len(dayIndex)-n_ignoredDays):
day +=1
ax.plot(pcaMatrix.Y[dayIndex[i]:dayIndex[i+1]:nPoints,0],
pcaMatrix.Y[dayIndex[i]:dayIndex[i+1]:nPoints,1],
pcaMatrix.Y[dayIndex[i]:dayIndex[i+1]:nPoints,2],
'o', markersize=7, alpha=0.1, label = 'day: ' + str(day))
ax = axSettings(ax)
def plotEachDay(pcaMatrices, dayIndex, nPoints = 10, n_ignoredDays = 1):
fig = plt.figure()
n_cols = 5
if dayIndex == None:
n_rows = ceil(len(pcaMatrices)/5.)
for day, matrix in enumerate(pcaMatrices):
ax = fig.add_subplot(n_rows,n_cols,day+1, projection='3d')
ax.plot(matrix.Y[:,0],matrix.Y[:,1],matrix.Y[:,2],
'o', markersize=7, alpha=0.5, label = 'day: ' + str(day+1))
ax = axSettings(ax)
else:
n_rows = ceil(len(dayIndex)/5.)
day = 0
for i in range(len(dayIndex)- n_ignoredDays):
day +=1
ax= fig.add_subplot(n_rows,n_cols,day, projection='3d')
ax.plot(pcaMatrices.Y[dayIndex[i]:dayIndex[i+1]:20,0],
pcaMatrices.Y[dayIndex[i]:dayIndex[i+1]:20,1],
pcaMatrices.Y[dayIndex[i]:dayIndex[i+1]:20,2],
'o', markersize=7, alpha=0.5, label = 'day: ' + str(day))
ax = axSettings(ax)
def weigthPlot(matrix):
print matrix.Wt.shape
imshow(abs(matrix.Wt), interpolation="nearest", origin = "lower")
def plotTrialType(pcaMatrices, trialType, nPoints=10, day = None):
fig = plt.figure()
if day == None:
n_cols = 5
n_rows = ceil(len(pcaMatrices)/5.)
for day, matrix in enumerate(pcaMatrices):
ax = fig.add_subplot(n_rows,n_cols,day+1, projection='3d')
ax.scatter(matrix.Y[::nPoints,0],matrix.Y[::nPoints,1],matrix.Y[::nPoints,2],
marker = 'o', alpha=0.5, c = array(trialType[day])[::nPoints])
ax = axSettings(ax)
else:
ax = fig.add_subplot(111, projection='3d')
ax.scatter(pcaMatrices[day].Y[::nPoints,0],pcaMatrices[day].Y[::nPoints,1],pcaMatrices[day].Y[::nPoints,2],
marker = 'o', alpha=0.5, c = array(trialType[day])[::nPoints])
ax = axSettings(ax)
def plotTrialIndex(pcaMatrices, trialIndex, nPoints = 10, day = None):
fig = plt.figure()
if day == None:
n_cols = 5
n_rows = ceil(len(pcaMatrices)/5.)
for day, matrix in enumerate(pcaMatrices):
ax = fig.add_subplot(n_rows,n_cols,day+1, projection='3d')
ax.scatter(matrix.Y[::nPoints,0],
matrix.Y[::nPoints,1],
matrix.Y[::nPoints,2],
marker = 'o', alpha=0.5, c = array(trialIndex[day])[::nPoints])
ax = axSettings(ax)
else:
ax = fig.add_subplot(111, projection='3d')
ax.scatter(pcaMatrices[day].Y[::nPoints,0],
pcaMatrices[day].Y[::nPoints,1],
pcaMatrices[day].Y[::nPoints,2],
marker = 'o', alpha=0.5, c = array(trialIndex[day])[::nPoints])
ax = axSettings(ax)
In [8]:
h5_file = h5py.File("data/koendata.h5", 'r')
animal_A_h5 = h5_file.get("/an0074")
animal_B_h5 = h5_file.get("/an0075")
In [9]:
animal_A = data(animal_A_h5, n_ignoredDays = 1)
animal_A.pcaOnEachDay()
animal_A.pcaOnAllDays()
animal_B = data(animal_B_h5, n_ignoredDays = 2)
animal_B.pcaOnEachDay()
animal_B.pcaOnAllDays()
In [12]:
%pylab qt
plotTrialType(animal_B.pcaMatrices, animal_B.trialTypeVecs,nPoints = 20)
plotTrialType(animal_A.pcaMatrices, animal_A.trialTypeVecs,nPoints = 20)
In [11]:
%pylab inline
plotEachDay(animal_B.pcaMatrix,animal_B.dayIndex, n_ignoredDays=2)
plotEachDay(animal_A.pcaMatrix,animal_A.dayIndex, n_ignoredDays=1)
plotEachDay(animal_B.pcaMatrices,None)
plotEachDay(animal_A.pcaMatrices,None)
plotAllDays(animal_B.pcaMatrix,animal_B.dayIndex,n_ignoredDays=2)
plotAllDays(animal_A.pcaMatrix,animal_A.dayIndex,n_ignoredDays=1)
plotAllDays(animal_B.pcaMatrices,None)
plotAllDays(animal_A.pcaMatrices,None)
In [11]: