# MOCAP Modelling with Second Order Differential Equation¶

## MOCAP Data¶

Motion capture data consists of x,y,z point clouds of a person, animal or character moving around. It is either captured directly from a special multiple camera set up, or it can be created by an animator in software. It often takes the form of a high dimensional time series.

The CMU Mocap database provides a public resource of motion capture data recordings of humans participating in every day activities. The pods software provides facilities for loading in this data for modelling.

In [1]:
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import matplotlib.pyplot as plt
import numpy as np
import pods

# Load in subject 35, activity 01 from the CMU mocap data base.
data = pods.datasets.cmu_mocap('35',['01'])

# Shift the data to start from 0 and scale it by its standard deviation
data2 = data['Y']
data2 = (data2 - data2[0, :])/data2.std(axis=0)
data2.shape

/Users/neil/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/pytz/__init__.py:29: UserWarning: Module pods was already imported from /Users/neil/sods/ods/pods/__init__.pyc, but /Users/neil/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages is being added to sys.path
from pkg_resources import resource_stream

Out[1]:
(90, 62)

The motion we have loaded in is of subject 35 walking. It consists of 62 channels and 90 points in the time series. As an example we will consider removing some of the channels of motion and predicting what they should be using the Gaussian process. First we create the data set by extracting the relevant channels.

In [2]:
index = np.empty([data2.size, 1], dtype = int)
Y = np.empty([data2.size, 1])
t = np.empty([data2.size, 1])
ti = np.arange(0, data2.shape[0])/30.
d = data2.shape[1]
start = 0
for k in range(d):
t[start:start+data2.shape[0], 0] = ti
Y[start:start+data2.shape[0], 0] = data2[:, k]
index[start:start+data2.shape[0], 0] = k*np.ones(ti.size, dtype = int)
start = start + data2.shape[0]

X = np.hstack((t, index))


Whilst there are only 90 time points, there are 62 channels in the time series. Since we are consdering a multiple output process this gives a GP with around five thousand data points. We use the sparse GP regression to allow us to fit the model.

In [3]:
#Inducing inputs definition
Mq = 30 #Number of inducing points per latent function
q = np.int_(data2.shape[1]/9) #Number of latent functions
indexz = np.empty([Mq*q, 1], dtype = int)
z = np.empty([Mq*q, 1])
start = 0
zi = np.linspace(ti.min(), ti.max(), Mq)
for i in range(q):
z[start:start+Mq] = zi.reshape((zi.size, 1))
indexz[start:start+Mq] = i*np.ones((zi.size, 1), dtype = int)
start = start + Mq
indexz += d
Z = np.hstack((z, indexz))


Now, We erase some frames from channel 8 to show how the ODE2 kernel is capable of recovering this samples.

In [4]:
cut = np.where((t[:, 0] > 0.5) & (t[:, 0] < 1.5) & (index[:, 0] == 8))
Y_train = np.delete(Y, cut, 0)
X_train = np.delete(X, cut, 0)

ind = np.where(X_train[:, 1] == 8) #Change the number in this line to change the output
plt.plot(X_train[ind[0], 0], Y_train[ind[0]], 'x')
ind = np.where(X[:, 1] == 8)
plt.plot(X[cut[0], 0], Y[cut[0]], 'o')

Out[4]:
[<matplotlib.lines.Line2D at 0x10db48cd0>]

Green circles from previous plots are the missing samples. To predict them we require to define the kernel and type of inference for the GP. Kernel is created using the second order differential equation covariance function. We train the model based on SparseGPRegression (model optimization can take several minutes)

In [13]:
import GPy
K = GPy.kern.EQ_ODE2(2, output_dim=d, rank=q)
m = GPy.models.SparseGPRegression(X=X_train, Y=Y_train, Z=Z, kernel=K)
m.kern.lengthscale = 0.1
m.inducing_inputs.fix()
#m.optimize('bfgs', max_iters=100)

Out[13]:
array([  0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,
13,  14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24,  25,
26,  27,  28,  29,  30,  31,  32,  33,  34,  35,  36,  37,  38,
39,  40,  41,  42,  43,  44,  45,  46,  47,  48,  49,  50,  51,
52,  53,  54,  55,  56,  57,  58,  59,  60,  61,  62,  63,  64,
65,  66,  67,  68,  69,  70,  71,  72,  73,  74,  75,  76,  77,
78,  79,  80,  81,  82,  83,  84,  85,  86,  87,  88,  89,  90,
91,  92,  93,  94,  95,  96,  97,  98,  99, 100, 101, 102, 103,
104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116,
117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129,
130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142,
143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155,
156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168,
169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181,
182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194,
195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207,
208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220,
221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233,
234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246,
247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259,
260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272,
273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285,
286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298,
299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311,
312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324,
325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337,
338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350,
351, 352, 353, 354, 355, 356, 357, 358, 359])

Next, we show the prediction for channel 8.

In [14]:
Yp_mean = m.predict(X)
Yp = Yp_mean[0][:, 0]
ind = np.where(index == 8) #Change the number in this line to change the output
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
plt.plot(X[ind[0], 0], Yp[ind[0]], X[ind[0], 0], Y[ind[0]])

Out[14]:
[<matplotlib.lines.Line2D at 0x10fd79410>,
<matplotlib.lines.Line2D at 0x10fd79690>]
In [ ]:
m.plot_f(fixed_inputs=[[1, 2]])

In [ ]:
X

In [ ]:
from IPython.display import display
display(m)

In [ ]:
m.eq_ode2.W

In [ ]: