cmu_mocap() example mostly working except some fiddling with axes for visualization. Also changes to naming of scaling and offset parameters in GP.py and deal with the case where the scale parameter is zero.

This commit is contained in:
Neil Lawrence 2013-04-27 10:39:55 +01:00
parent d7ac1d025b
commit ac842d51e6
6 changed files with 202 additions and 98 deletions

View file

@ -2,19 +2,30 @@ import numpy as np
from likelihood import likelihood
class Gaussian(likelihood):
"""
Likelihood class for doing Expectation propagation
:param Y: observed output (Nx1 numpy.darray)
..Note:: Y values allowed depend on the likelihood_function used
:param variance :
:param normalize: whether to normalize the data before computing (predictions will be in original scales)
:type normalize: False|True
"""
def __init__(self,data,variance=1.,normalize=False):
self.is_heteroscedastic = False
self.Nparams = 1
self.Z = 0. # a correction factor which accounts for the approximation made
N, self.D = data.shape
#normaliztion
#normalization
if normalize:
self._mean = data.mean(0)[None,:]
self._std = data.std(0)[None,:]
self._bias = data.mean(0)[None,:]
self._scale = data.std(0)[None,:]
# Don't scale outputs which have zero variance to zero.
self._scale[np.nonzero(self._scale==0.)] = 1.0e-3
else:
self._mean = np.zeros((1,self.D))
self._std = np.ones((1,self.D))
self._bias = np.zeros((1,self.D))
self._scale = np.ones((1,self.D))
self.set_data(data)
@ -24,7 +35,7 @@ class Gaussian(likelihood):
self.data = data
self.N,D = data.shape
assert D == self.D
self.Y = (self.data - self._mean)/self._std
self.Y = (self.data - self._bias)/self._scale
if D > self.N:
self.YYT = np.dot(self.Y,self.Y.T)
self.trYYT = np.trace(self.YYT)
@ -47,19 +58,19 @@ class Gaussian(likelihood):
"""
Un-normalize the prediction and add the likelihood variance, then return the 5%, 95% interval
"""
mean = mu*self._std + self._mean
mean = mu*self._scale + self._bias
if full_cov:
if self.D >1:
raise NotImplementedError, "TODO"
#Note. for D>1, we need to re-normalise all the outputs independently.
# This will mess up computations of diag(true_var), below.
#note that the upper, lower quantiles should be the same shape as mean
true_var = (var + np.eye(var.shape[0])*self._variance)*self._std**2
_5pc = mean + - 2.*np.sqrt(np.diag(true_var))
true_var = (var + np.eye(var.shape[0])*self._variance)*self._scale**2
_5pc = mean - 2.*np.sqrt(np.diag(true_var))
_95pc = mean + 2.*np.sqrt(np.diag(true_var))
else:
true_var = (var + self._variance)*self._std**2
_5pc = mean + - 2.*np.sqrt(true_var)
true_var = (var + self._variance)*self._scale**2
_5pc = mean - 2.*np.sqrt(true_var)
_95pc = mean + 2.*np.sqrt(true_var)
return mean, true_var, _5pc, _95pc

View file

@ -19,8 +19,6 @@ class GP(model):
:parm likelihood: a GPy likelihood
:param normalize_X: whether to normalize the input data before computing (predictions will be in original scales)
:type normalize_X: False|True
:param normalize_Y: whether to normalize the input data before computing (predictions will be in original scales)
:type normalize_Y: False|True
:param Xslices: how the X,Y data co-vary in the kernel (i.e. which "outputs" they correspond to). See (link:slicing)
:rtype: model object
:param epsilon_ep: convergence criterion for the Expectation Propagation algorithm, defaults to 0.1

View file

@ -24,12 +24,12 @@ class GPLVM(GP):
:type init: 'PCA'|'random'
"""
def __init__(self, Y, Q, init='PCA', X = None, kernel=None, **kwargs):
def __init__(self, Y, Q, init='PCA', X = None, kernel=None, normalize_Y=False, **kwargs):
if X is None:
X = self.initialise_latent(init, Q, Y)
if kernel is None:
kernel = kern.rbf(Q) + kern.bias(Q)
likelihood = Gaussian(Y)
likelihood = Gaussian(Y, normalize=normalize_Y)
GP.__init__(self, X, likelihood, kernel, **kwargs)
def initialise_latent(self, init, Q, Y):

View file

@ -225,49 +225,92 @@ def creep_data():
X = all_data[:, features].copy()
return {'X': X, 'y' : y}
def cmu_35_walk_jog():
def cmu_mocap_49_balance():
"""Load CMU subject 49's one legged balancing motion that was used by Alvarez, Luengo and Lawrence at AISTATS 2009."""
train_motions = ['18', '19']
test_motions = ['20']
data = cmu_mocap('49', train_motions, test_motions, sample_every=4)
data['info'] = "One legged balancing motions from CMU data base subject 49. As used in Alvarez, Luengo and Lawrence at AISTATS 2009. It consists of " + data['info']
return data
skel = GPy.util.mocap.acclaim_skeleton(os.path.join(data_path, 'mocap', 'cmu', '35', '35.asf'))
examples = ['01', '02', '03', '04', '05', '06',
def cmu_mocap_35_walk_jog():
"""Load CMU subject 35's walking and jogging motions, the same data that was used by Taylor, Roweis and Hinton at NIPS 2007. but without their preprocessing. Also used by Lawrence at AISTATS 2007."""
train_motions = ['01', '02', '03', '04', '05', '06',
'07', '08', '09', '10', '11', '12',
'13', '14', '15', '16', '17', '19',
'20', '21', '22', '23', '24', '25',
'26', '28', '30', '31', '32', '33', '34']
test_examples = ['18', '29']
# Label differently for each sequence
exlbls = np.eye(31)
testexlbls = np.eye(2)
tot_length = 0
tot_test_length = 0
tY = []
tlbls = []
for i in range(len(examples)):
tmpchan = skel.load_channels(os.path.join(data_path, 'mocap', 'cmu', '35', '35_' + examples[i] + '.amc'))
tY.append(tmpchan[::4, :])
tlbls.append(np.tile(exlbls[i, :], (tY[i].shape[0], 1)))
tot_length += tY[i].shape[0]
Y = np.zeros((tot_length, tY[0].shape[1]))
lbls = np.zeros((tot_length, tlbls[0].shape[1]))
endInd = 0
for i in range(len(tY)):
startInd = endInd
endInd += tY[i].shape[0]
Y[startInd:endInd, :] = tY[i]
lbls[startInd:endInd, :] = tlbls[i]
tYtest = []
tlblstest = []
for i in range(len(test_examples)):
tmpchan = skel.load_channels(os.path.join(data_path, 'mocap', 'cmu', '35', '35_' + test_examples[i] + '.amc'))
tYtest.append(tmpchan[::4, :])
tlblstest.append(np.tile(testexlbls[i, :], (tYtest[i].shape[0], 1)))
tot_test_length += tYtest[i].shape[0]
test_motions = ['18', '29']
data = cmu_mocap('35', train_motions, test_motions, sample_every=4)
data['info'] = "Walk and jog data from CMU data base subject 35. As used in Tayor, Roweis and Hinton at NIPS 2007, but without their pre-processing (i.e. as used by Lawrence at AISTATS 2007). It consists of " + data['info']
return data
Ytest = np.zeros((tot_test_length, tYtest[0].shape[1]))
lblstest = np.zeros((tot_test_length, tlblstest[0].shape[1]))
endInd = 0
for i in range(len(tYtest)):
startInd = endInd
endInd += tYtest[i].shape[0]
Ytest[startInd:endInd, :] = tYtest[i]
lblstest[startInd:endInd, :] = tlblstest[i]
return {'Y': Y, 'lbls' : lbls, 'Ytest': Ytest, 'lblstest' : lblstest, 'info': "Walk and jog data from CMU data base subject 35."}
def cmu_mocap(subject, train_motions, test_motions=[], sample_every=4):
"""Load a given subject's training and test motions from the CMU motion capture data."""
# Load in subject skeleton.
subject_dir = os.path.join(data_path, 'mocap', 'cmu', subject)
skel = GPy.util.mocap.acclaim_skeleton(os.path.join(subject_dir, subject + '.asf'))
# Set up labels for each sequence
exlbls = np.eye(len(train_motions))
# Load sequences
tot_length = 0
temp_Y = []
temp_lbls = []
for i in range(len(train_motions)):
temp_chan = skel.load_channels(os.path.join(subject_dir, subject + '_' + train_motions[i] + '.amc'))
temp_Y.append(temp_chan[::sample_every, :])
temp_lbls.append(np.tile(exlbls[i, :], (temp_Y[i].shape[0], 1)))
tot_length += temp_Y[i].shape[0]
Y = np.zeros((tot_length, temp_Y[0].shape[1]))
lbls = np.zeros((tot_length, temp_lbls[0].shape[1]))
end_ind = 0
for i in range(len(temp_Y)):
start_ind = end_ind
end_ind += temp_Y[i].shape[0]
Y[start_ind:end_ind, :] = temp_Y[i]
lbls[start_ind:end_ind, :] = temp_lbls[i]
if len(test_motions)>0:
temp_Ytest = []
temp_lblstest = []
testexlbls = np.eye(len(test_motions))
tot_test_length = 0
for i in range(len(test_motions)):
temp_chan = skel.load_channels(os.path.join(subject_dir, subject + '_' + test_motions[i] + '.amc'))
temp_Ytest.append(temp_chan[::sample_every, :])
temp_lblstest.append(np.tile(testexlbls[i, :], (temp_Ytest[i].shape[0], 1)))
tot_test_length += temp_Ytest[i].shape[0]
# Load test data
Ytest = np.zeros((tot_test_length, temp_Ytest[0].shape[1]))
lblstest = np.zeros((tot_test_length, temp_lblstest[0].shape[1]))
end_ind = 0
for i in range(len(temp_Ytest)):
start_ind = end_ind
end_ind += temp_Ytest[i].shape[0]
Ytest[start_ind:end_ind, :] = temp_Ytest[i]
lblstest[start_ind:end_ind, :] = temp_lblstest[i]
else:
Ytest = None
lblstest = None
info = 'Subject: ' + subject + '. Training motions: '
for motion in train_motions:
info += motion + ', '
info = info[:-2]
if len(test_motions)>0:
info += '. Test motions: '
for motion in test_motions:
info += motion + ', '
info = info[:-2] + '.'
else:
info += '.'
if sample_every != 1:
info += ' Data is sub-sampled to every ' + str(sample_every) + ' frames.'
return {'Y': Y, 'lbls' : lbls, 'Ytest': Ytest, 'lblstest' : lblstest, 'info': info, 'skel': skel}

View file

@ -157,6 +157,13 @@ class skeleton(tree):
def __init__(self):
tree.__init__(self)
def connection_matrix(self):
connection = np.zeros((len(self.vertices), len(self.vertices)), dtype=bool)
for i in range(len(self.vertices)):
for j in range(len(self.vertices[i].children)):
connection[i, self.vertices[i].children[j]] = True
return connection
def to_xyz(self, channels):
raise NotImplementedError, "this needs to be implemented to use the skeleton class"
@ -557,6 +564,7 @@ class acclaim_skeleton(skeleton):
lin = self.read_line(fid)
else:
raise Error, 'Unrecognised file format'
self.finalize()
def read_units(self, fid):
"""Read units from an acclaim skeleton file stream."""

View file

@ -184,71 +184,115 @@ class image_show(data_show):
#if self.invert:
# self.vals = -self.vals
class stick_show(data_show):
"""Show a three dimensional point cloud as a figure. Connect elements of the figure together using the matrix connect."""
class mocap_data_show(data_show):
"""Base class for visualizing motion capture data."""
def __init__(self, vals, axes=None, connect=None):
if axes==None:
fig = plt.figure()
axes = fig.add_subplot(111, projection='3d')
data_show.__init__(self, vals, axes)
self.vals = vals.reshape((3, vals.shape[1]/3)).T
self.x_lim = np.array([self.vals[:, 0].min(), self.vals[:, 0].max()])
self.y_lim = np.array([self.vals[:, 1].min(), self.vals[:, 1].max()])
self.z_lim = np.array([self.vals[:, 2].min(), self.vals[:, 2].max()])
self.points_handle = self.axes.scatter(self.vals[:, 0], self.vals[:, 1], self.vals[:, 2])
self.axes.set_xlim(self.x_lim)
self.axes.set_ylim(self.y_lim)
self.axes.set_zlim(self.z_lim)
self.axes.set_aspect(1)
self.axes.autoscale(enable=False)
self.connect = connect
if not self.connect==None:
x = []
y = []
z = []
self.I, self.J = np.nonzero(self.connect)
for i in range(len(self.I)):
x.append(self.vals[self.I[i], 0])
x.append(self.vals[self.J[i], 0])
x.append(np.NaN)
y.append(self.vals[self.I[i], 1])
y.append(self.vals[self.J[i], 1])
y.append(np.NaN)
z.append(self.vals[self.I[i], 2])
z.append(self.vals[self.J[i], 2])
z.append(np.NaN)
self.line_handle = self.axes.plot(np.array(x), np.array(y), np.array(z), 'b-')
self.process_values(vals)
self.initialize_axes()
self.draw_vertices()
self.finalize_axes()
self.draw_edges()
self.axes.figure.canvas.draw()
def modify(self, vals):
self.points_handle.remove()
self.line_handle[0].remove()
self.vals = vals.reshape((3, vals.shape[1]/3)).T
def draw_vertices(self):
self.points_handle = self.axes.scatter(self.vals[:, 0], self.vals[:, 1], self.vals[:, 2])
self.axes.set_xlim(self.x_lim)
self.axes.set_ylim(self.y_lim)
self.axes.set_zlim(self.z_lim)
def draw_edges(self):
self.line_handle = []
if not self.connect==None:
x = []
y = []
z = []
self.I, self.J = np.nonzero(self.connect)
for i in range(len(self.I)):
x.append(self.vals[self.I[i], 0])
x.append(self.vals[self.J[i], 0])
for i, j in zip(self.I, self.J):
x.append(self.vals[i, 0])
x.append(self.vals[j, 0])
x.append(np.NaN)
y.append(self.vals[self.I[i], 1])
y.append(self.vals[self.J[i], 1])
y.append(self.vals[i, 1])
y.append(self.vals[j, 1])
y.append(np.NaN)
z.append(self.vals[self.I[i], 2])
z.append(self.vals[self.J[i], 2])
z.append(self.vals[i, 2])
z.append(self.vals[j, 2])
z.append(np.NaN)
self.line_handle = self.axes.plot(np.array(x), np.array(y), np.array(z), 'b-')
def modify(self, vals):
self.process_values(vals)
self.initialize_axes_modify()
self.draw_vertices()
self.finalize_axes_modify()
self.draw_edges()
self.axes.figure.canvas.draw()
def process_values(self, vals):
raise NotImplementedError, "this needs to be implemented to use the data_show class"
def initialize_axes(self):
"""Set up the axes with the right limits and scaling."""
self.x_lim = np.array([self.vals[:, 0].min(), self.vals[:, 0].max()])
self.y_lim = np.array([self.vals[:, 1].min(), self.vals[:, 1].max()])
self.z_lim = np.array([self.vals[:, 2].min(), self.vals[:, 2].max()])
def initialize_axes_modify(self):
self.points_handle.remove()
self.line_handle[0].remove()
def finalize_axes(self):
self.axes.set_xlim(self.x_lim)
self.axes.set_ylim(self.y_lim)
self.axes.set_zlim(self.z_lim)
self.axes.set_aspect(1)
self.axes.autoscale(enable=False)
def finalize_axes_modify(self):
self.axes.set_xlim(self.x_lim)
self.axes.set_ylim(self.y_lim)
self.axes.set_zlim(self.z_lim)
class stick_show(mocap_data_show):
"""Show a three dimensional point cloud as a figure. Connect elements of the figure together using the matrix connect."""
def __init__(self, vals, axes=None, connect=None):
mocap_data_show.__init__(self, vals, axes, connect)
def process_values(self, vals):
self.vals = vals.reshape((3, vals.shape[1]/3)).T
class skeleton_show(mocap_data_show):
"""data_show class for visualizing motion capture data encoded as a skeleton with angles."""
def __init__(self, vals, skel, padding=0, axes=None):
self.skel = skel
self.padding = padding
connect = skel.connection_matrix()
mocap_data_show.__init__(self, vals, axes, connect)
def process_values(self, vals):
if self.padding>0:
channels = np.zeros((vals.shape[0], vals.shape[1]+self.padding))
channels[:, 0:vals.shape[0]] = vals
else:
channels = vals
vals_mat = self.skel.to_xyz(channels.flatten())
self.vals = vals_mat
# Flip the Y and Z axes
self.vals[:, 0] = vals_mat[:, 0]
self.vals[:, 1] = vals_mat[:, 2]
self.vals[:, 2] = vals_mat[:, 1]
def wrap_around(vals, lim, connect):
quot = lim[1] - lim[0]
vals = rem(vals, quot)+lim[0]
nVals = floor(vals/quot)
for i in range(connect.shape[0]):
for j in find(connect[i, :]):
if nVals[i] != nVals[j]:
connect[i, j] = False
return vals, connect