{ "metadata": { "name": "", "signature": "sha256:09e41d07a6ce90010de9160f19393792ca5f47869101292192431caffd90d63b" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "DRAFT\n", "===\n", "\n", "# MOCAP Modelling with Second Order Differential Equation\n", "\n", "## Gaussian Process Winter School, Genova, Italy\n", "\n", "### 21st January 2015\n", "\n", "### Cristian Guarnizo and Neil D. Lawrence\n", "\n", "## MOCAP Data\n", "\n", "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.\n", "\n", "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." ] }, { "cell_type": "code", "collapsed": false, "input": [ "%matplotlib inline\n", "%config InlineBackend.figure_format = 'svg'\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pods\n", "\n", "# Load in subject 35, activity 01 from the CMU mocap data base.\n", "data = pods.datasets.cmu_mocap('35',['01'])\n", "\n", "# Shift the data to start from 0 and scale it by its standard deviation\n", "data2 = data['Y']\n", "data2 = (data2 - data2[0, :])/data2.std(axis=0)\n", "data2.shape" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stderr", "text": [ "/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\n", " from pkg_resources import resource_stream\n" ] }, { "metadata": {}, "output_type": "pyout", "prompt_number": 1, "text": [ "(90, 62)" ] } ], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "collapsed": false, "input": [ "index = np.empty([data2.size, 1], dtype = int)\n", "Y = np.empty([data2.size, 1])\n", "t = np.empty([data2.size, 1])\n", "ti = np.arange(0, data2.shape[0])/30.\n", "d = data2.shape[1]\n", "start = 0\n", "for k in range(d):\n", " t[start:start+data2.shape[0], 0] = ti\n", " Y[start:start+data2.shape[0], 0] = data2[:, k]\n", " index[start:start+data2.shape[0], 0] = k*np.ones(ti.size, dtype = int)\n", " start = start + data2.shape[0]\n", "\n", "X = np.hstack((t, index))" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "collapsed": false, "input": [ "#Inducing inputs definition\n", "Mq = 30 #Number of inducing points per latent function\n", "q = np.int_(data2.shape[1]/9) #Number of latent functions\n", "indexz = np.empty([Mq*q, 1], dtype = int)\n", "z = np.empty([Mq*q, 1])\n", "start = 0\n", "zi = np.linspace(ti.min(), ti.max(), Mq)\n", "for i in range(q):\n", " z[start:start+Mq] = zi.reshape((zi.size, 1))\n", " indexz[start:start+Mq] = i*np.ones((zi.size, 1), dtype = int)\n", " start = start + Mq\n", "indexz += d\n", "Z = np.hstack((z, indexz))" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, We erase some frames from channel 8 to show how the ODE2 kernel is capable of recovering this samples." ] }, { "cell_type": "code", "collapsed": false, "input": [ "cut = np.where((t[:, 0] > 0.5) & (t[:, 0] < 1.5) & (index[:, 0] == 8))\n", "Y_train = np.delete(Y, cut, 0)\n", "X_train = np.delete(X, cut, 0)\n", "\n", "ind = np.where(X_train[:, 1] == 8) #Change the number in this line to change the output\n", "plt.plot(X_train[ind[0], 0], Y_train[ind[0]], 'x')\n", "ind = np.where(X[:, 1] == 8)\n", "plt.plot(X[cut[0], 0], Y[cut[0]], 'o')" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 4, "text": [ "[]" ] }, { "metadata": {}, "output_type": "display_data", "svg": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text": [ "" ] } ], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": {}, "source": [ "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)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "import GPy\n", "K = GPy.kern.EQ_ODE2(2, output_dim=d, rank=q)\n", "m = GPy.models.SparseGPRegression(X=X_train, Y=Y_train, Z=Z, kernel=K)\n", "m.kern.lengthscale = 0.1\n", "m.inducing_inputs.fix()\n", "#m.optimize('bfgs', max_iters=100)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 13, "text": [ "array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,\n", " 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25,\n", " 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,\n", " 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51,\n", " 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64,\n", " 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77,\n", " 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90,\n", " 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,\n", " 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116,\n", " 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129,\n", " 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142,\n", " 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155,\n", " 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168,\n", " 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181,\n", " 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194,\n", " 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207,\n", " 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220,\n", " 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233,\n", " 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246,\n", " 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259,\n", " 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272,\n", " 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285,\n", " 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298,\n", " 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311,\n", " 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324,\n", " 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337,\n", " 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350,\n", " 351, 352, 353, 354, 355, 356, 357, 358, 359])" ] } ], "prompt_number": 13 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we show the prediction for channel 8." ] }, { "cell_type": "code", "collapsed": false, "input": [ "Yp_mean = m.predict(X)\n", "Yp = Yp_mean[0][:, 0]\n", "ind = np.where(index == 8) #Change the number in this line to change the output\n", "%matplotlib inline\n", "%config InlineBackend.figure_format = 'svg'\n", "plt.plot(X[ind[0], 0], Yp[ind[0]], X[ind[0], 0], Y[ind[0]])" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 14, "text": [ "[,\n", " ]" ] }, { "metadata": {}, "output_type": "display_data", "svg": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text": [ "" ] } ], "prompt_number": 14 }, { "cell_type": "code", "collapsed": false, "input": [ "m.plot_f(fixed_inputs=[[1, 2]])" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "X" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "from IPython.display import display\n", "display(m)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "m.eq_ode2.W" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }