custom gp.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate Dummy Data #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X = np.linspace(0,3*np.pi,31)[:,None]\n",
    "Y = np.sin(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup of `GPy` #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# GPy normalises Y, does this need to be done?\n",
    "\n",
    "kern       = kern.RBF(X.shape[1])\n",
    "likelihood = likelihoods.Gaussian(variance=noise_var)\n",
    "\n",
    "inference_method = exact_gaussian_inference.ExactGaussianInference()\n",
    "\n",
    "signal_variance = 1**2\n",
    "lengthscale     = 1**2\n",
    "noise_variance  = 1**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parameters Changed #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def _symmetrify_numpy(A, upper=False):\n",
    "    triu = np.triu_indices_from(A,k=1)\n",
    "    if upper:\n",
    "        A.T[triu] = A[triu]\n",
    "    else:\n",
    "        A[triu] = A.T[triu]\n",
    "\n",
    "def symmetrify(A, upper=False):\n",
    "    \"\"\"\n",
    "    Take the square matrix A and make it symmetrical by copting elements from\n",
    "    the lower half to the upper\n",
    "\n",
    "    works IN PLACE.\n",
    "\n",
    "    note: tries to use cython, falls back to a slower numpy version\n",
    "    \"\"\"\n",
    "    #if config.getboolean('cython', 'working'):\n",
    "    #    _symmetrify_cython(A, upper)\n",
    "    #else:\n",
    "        _symmetrify_numpy(A, upper)\n",
    "\n",
    "def force_F_ordered(A):\n",
    "    \"\"\"\n",
    "    return a F ordered version of A, assuming A is triangular\n",
    "    \"\"\"\n",
    "    if A.flags['F_CONTIGUOUS']:\n",
    "        return A\n",
    "    print(\"why are your arrays not F order?\")\n",
    "    return np.asfortranarray(A)\n",
    "\n",
    "def dpotri(A, lower=1):\n",
    "    \"\"\"\n",
    "    Wrapper for lapack dpotri function\n",
    "\n",
    "    DPOTRI - compute the inverse of a real symmetric positive\n",
    "      definite matrix A using the Cholesky factorization A =\n",
    "      U**T*U or A = L*L**T computed by DPOTRF\n",
    "\n",
    "    :param A: Matrix A\n",
    "    :param lower: is matrix lower (true) or upper (false)\n",
    "    :returns: A inverse\n",
    "\n",
    "    \"\"\"\n",
    "    if _fix_dpotri_scipy_bug:\n",
    "        assert lower==1, \"scipy linalg behaviour is very weird. please use lower, fortran ordered arrays\"\n",
    "        lower = 0\n",
    "\n",
    "    A = force_F_ordered(A)\n",
    "    R, info = lapack.dpotri(A, lower=lower) #needs to be zero here, seems to be a scipy bug\n",
    "\n",
    "    symmetrify(R)\n",
    "    return R, info\n",
    "\n",
    "def dtrtri(L):\n",
    "    \"\"\"\n",
    "    Inverts a Cholesky lower triangular matrix\n",
    "\n",
    "    :param L: lower triangular matrix\n",
    "    :rtype: inverse of L\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    L = force_F_ordered(L)\n",
    "    return lapack.dtrtri(L, lower=1)[0]\n",
    "\n",
    "def dpotrs(A, B, lower=1):\n",
    "    \"\"\"\n",
    "    Wrapper for lapack dpotrs function\n",
    "    :param A: Matrix A\n",
    "    :param B: Matrix B\n",
    "    :param lower: is matrix lower (true) or upper (false)\n",
    "    :returns:\n",
    "    \"\"\"\n",
    "    A = force_F_ordered(A)\n",
    "    return lapack.dpotrs(A, B, lower=lower)\n",
    "\n",
    "def jitchol(A, maxtries=5):\n",
    "    A = np.ascontiguousarray(A)\n",
    "    L, info = lapack.dpotrf(A, lower=1)\n",
    "    if info == 0:\n",
    "        return L\n",
    "    else:\n",
    "        diagA = np.diag(A)\n",
    "        if np.any(diagA <= 0.):\n",
    "            raise linalg.LinAlgError(\"not pd: non-positive diagonal elements\")\n",
    "        jitter = diagA.mean() * 1e-6\n",
    "        num_tries = 1\n",
    "        while num_tries <= maxtries and np.isfinite(jitter):\n",
    "            try:\n",
    "                L = linalg.cholesky(A + np.eye(A.shape[0]) * jitter, lower=True)\n",
    "                return L\n",
    "            except:\n",
    "                jitter *= 10\n",
    "            finally:\n",
    "                num_tries += 1\n",
    "        raise linalg.LinAlgError(\"not positive definite, even with jitter.\")\n",
    "    import traceback\n",
    "    try: raise\n",
    "    except:\n",
    "        logging.warning('\\n'.join(['Added jitter of {:.10e}'.format(jitter),\n",
    "            '  in '+traceback.format_list(traceback.extract_stack(limit=2)[-2:-1])[0][2:]]))\n",
    "    return L\n",
    "\n",
    "def pdinv(A, *args):\n",
    "    \"\"\"\n",
    "    :param A: A DxD pd numpy array\n",
    "\n",
    "    :rval Ai: the inverse of A\n",
    "    :rtype Ai: np.ndarray\n",
    "    :rval L: the Cholesky decomposition of A\n",
    "    :rtype L: np.ndarray\n",
    "    :rval Li: the Cholesky decomposition of Ai\n",
    "    :rtype Li: np.ndarray\n",
    "    :rval logdet: the log of the determinant of A\n",
    "    :rtype logdet: float64\n",
    "\n",
    "    \"\"\"\n",
    "    L = jitchol(A, *args)\n",
    "    logdet = 2.*np.sum(np.log(np.diag(L)))\n",
    "    Li = dtrtri(L)\n",
    "    Ai, _ = dpotri(L, lower=1)\n",
    "    # Ai = np.tril(Ai) + np.tril(Ai,-1).T\n",
    "    symmetrify(Ai)\n",
    "\n",
    "    return Ai, L, Li, logdet\n",
    "\n",
    "# GPy uses tdot_blas instead of np.dot (suggests it is good for \"large\" matricies)\n",
    "def unscaled_dist(X, X2=None):\n",
    "    \"\"\"\n",
    "    Compute the Euclidean distance between each row of X and X2, or between\n",
    "    each pair of rows of X if X2 is None.\n",
    "    \"\"\"\n",
    "    if X2 is None:\n",
    "        Xsq = np.sum(np.square(X),1)\n",
    "        r2 = -2.*np.dot(X, X.T) + (Xsq[:,None] + Xsq[None,:])\n",
    "        util.diag.view(r2)[:,]= 0. # force diagnoal to be zero: sometime numerically a little negative\n",
    "        r2 = np.clip(r2, 0, np.inf)\n",
    "        return np.sqrt(r2)\n",
    "    else:\n",
    "        X1sq = np.sum(np.square(X),1)\n",
    "        X2sq = np.sum(np.square(X2),1)\n",
    "        r2 = -2.*np.dot(X, X2.T) + X1sq[:,None] + X2sq[None,:]\n",
    "        r2 = np.clip(r2, 0, np.inf)\n",
    "\n",
    "# self.posterior, self._log_marginal_likelihood, self.grad_dict = self.inference_method.inference(self.kern, self.X, \n",
    "#     self.likelihood, self.Y_normalized, self.mean_function, self.Y_metadata)\n",
    "# m = 0\n",
    "# YYT_factor = self.get_YYTfactor(Y-m)\n",
    "\n",
    "YYT_factor = Y\n",
    "\n",
    "# K = kern.K(X)\n",
    "# r = self._scaled_dist(X, X2)\n",
    "# r = self._unscaled_dist(X, X2)/self.lengthscale\n",
    "r = unscaled_dist(X) / lengthscale:\n",
    "# K = self.K_of_r(r)\n",
    "K = signal_variance * np.exp(-0.5 * r**2)\n",
    "\n",
    "Ky = K.copy()\n",
    "diag.add(Ky, noise_variance+1e-8)\n",
    "Wi, LW, LWi, W_logdet = pdinv(Ky)\n",
    "\n",
    "alpha, _ = dpotrs(LW, YYT_factor, lower=1)\n",
    "\n",
    "log_marginal =  0.5*(-Y.size * log_2_pi - Y.shape[1] * W_logdet - np.sum(alpha * YYT_factor))\n",
    "\n",
    "dL_dK = 0.5 * (np.dot(alpha,alpha.T) - Y.shape[1] * Wi)\n",
    "\n",
    "dL_dthetaL = likelihood.exact_inference_gradients(np.diag(dL_dK),Y_metadata)\n",
    "\n",
    "return Posterior(woodbury_chol=LW, woodbury_vector=alpha, K=K), log_marginal, {'dL_dK':dL_dK, 'dL_dthetaL':dL_dthetaL, 'dL_dm':alpha}\n",
    "\n",
    "\n",
    "\n",
    "# self.likelihood.update_gradients(self.grad_dict['dL_dthetaL'])\n",
    "self.variance.gradient = grad\n",
    "\n",
    "# self.kern.update_gradients_full(self.grad_dict['dL_dK'], self.X)\n",
    "self.variance.gradient = np.sum(self.K(X, X2)* dL_dK)/self.variance\n",
    "dL_dr = self.dK_dr_via_X(X, X2) * dL_dK\n",
    "r = self._scaled_dist(X, X2)\n",
    "self.lengthscale.gradient = -np.sum(dL_dr*r)/self.lengthscale"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Evaluate Likelihood and Derivatives #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# call: m.log_likelihood()\n",
    "self._log_marginal_likelihood\n",
    "# call: m.gradient\n",
    "self.gradient"
   ]
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}