{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n\nBasic MAP Estimation\n====================\n\n.. currentmodule:: pydygp.linlatentforcemodels\n\nThis note descibes how to carry out the process of carrying out MAP\nparameter estimation for the MLFM using the Adaptive Gradient matching\napproximation. This uses the :class:`MLFMAdapGrad` object and so our\nfirst step is to import this object.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\nfrom scipy.integrate import odeint\nfrom pydygp.linlatentforcemodels import MLFMAdapGrad\nfrom sklearn.gaussian_process.kernels import RBF\nnp.random.seed(17)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Model Setup\n~~~~~~~~~~~\n\nTo begin we are going to demonstate the model with an ODE on the unit sphere\n\n\\begin{align}S^{2} = \\{ x \\in \\mathbb{R}^3 \\; : \\; \\| x \\| = 1 \\},\\end{align}\n\nwhich is given by the initial value problem\n\n\\begin{align}\\dot{\\mathbf{x}}(t) = \\mathbf{A}(t) \\mathbf{x}(t),\n    \\qquad \\mathbf{x}_0 \\in S^2,\\end{align}\n\nwhere the coefficient matrix, $\\mathbf{A}(t)$, is supported on the Lie\nalgebra $\\mathfrak{so}(3)$. We do this by representing the\n\n\\begin{align}\\mathbf{A}(t) = \\sum_{d=0}^3 \\beta_{0d}\\mathbf{L}_d +\n    \\sum_{r=1}^R g_r(t) \\sum_{d=1}^3 \\beta_{rd}\\mathbf{L}_d,\\end{align}\n\nwhere $\\{\\mathbf{L}_d \\}$ is a basis of the Lie algebra\n$\\mathfrak{so}(3)$. The :class:`so` object returns a tuple\nof basis elements for the Lie algebra, so for our example we will\nbe interested in :code:`so(3)`\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from pydygp.liealgebras import so\nfor d, item in enumerate(so(3)):\n    print(''.join(('\\n', 'L{}'.format(d+1))))\n    print(item)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Simulation\n~~~~~~~~~~\nTo simulate from the model we need to chose the set of coefficients\n$\\beta_{r, d}$. We will consider the model with a single latent\nforcing function, and randomly generate the variables $beta$\n\n:func:`pydygp.linlatentforcemodels.MLFMAdapGrad.sim`\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "g = lambda t: np.exp(-(t-2)**2) * np.cos(t)  # single latent force\nbeta = np.random.randn(2, 3)  \n\nA = [sum(brd*Ld for brd, Ld in zip(br, so(3)))\n     for br in beta]\n\nttd = np.linspace(0., 5., 100)\nx0 = [1., 0., 0.]\nsol = odeint(lambda x, t: (A[0] + g(t)*A[1]).dot(x),\n             x0,\n             ttd)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The MLFM Class\n~~~~~~~~~~~~~~\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mlfm = MLFMAdapGrad(so(3), R=1, lf_kernels=(RBF(), ))\n\nx0 = np.eye(3)\n\n# downsample the dense time vector\ntt = ttd[::10]\nData, _ = mlfm.sim(x0, tt, beta=beta, latent_forces=(g, ), size=3)\n\nfig, ax = plt.subplots()\nax.plot(ttd, sol, '-', alpha=0.3)\nax.plot(tt, Data[0], 'o')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Latent Force Estimation\n~~~~~~~~~~~~~~~~~~~~~~~\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Y = np.column_stack(y.T.ravel() for y in Data)\nres = mlfm.fit(tt, Y, beta0 = beta, beta_is_fixed=True)\n\n# predict the lf using the Laplace approximation\nEg, SDg = mlfm.predict_lf(ttd, return_std=True)\n\n# sphinx_gallery_thumbnail_number = 2\nfig2, ax = plt.subplots()\nax.plot(ttd, g(ttd), 'k-', alpha=0.8)\nax.plot(tt, res.g.T, 'o')\nfor Egr, SDgr in zip(Eg, SDg):\n    ax.fill_between(ttd,\n                    Egr + 2*SDgr, Egr - 2*SDgr,\n                    alpha=0.5)\n\nplt.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.0"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}