dpp_example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2, 1, 0]\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "Determinantal Point Process\n",
    "===========================\n",
    "\n",
    "This a pure Python port of the Matlab code of Alex Kulesza. \n",
    "\n",
    "The port was done using the fabulous SMOP package with alterations to ensure\n",
    "it only used scipy/numpy packages. As such various line comments referencing\n",
    "the original `*.m` file is retained in this code\n",
    "\n",
    "If there are any errors in the port please let me know on: chpmn.siu@gmail.com\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "import random\n",
    "import scipy\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "replicate `null`\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "from numpy.linalg import svd\n",
    "\n",
    "\n",
    "def rank(A, atol=1e-13, rtol=0):\n",
    "    \"\"\"Estimate the rank (i.e. the dimension of the nullspace) of a matrix.\n",
    "\n",
    "    The algorithm used by this function is based on the singular value\n",
    "    decomposition of `A`.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    A : ndarray\n",
    "        A should be at most 2-D.  A 1-D array with length n will be treated\n",
    "        as a 2-D with shape (1, n)\n",
    "    atol : float\n",
    "        The absolute tolerance for a zero singular value.  Singular values\n",
    "        smaller than `atol` are considered to be zero.\n",
    "    rtol : float\n",
    "        The relative tolerance.  Singular values less than rtol*smax are\n",
    "        considered to be zero, where smax is the largest singular value.\n",
    "\n",
    "    If both `atol` and `rtol` are positive, the combined tolerance is the\n",
    "    maximum of the two; that is::\n",
    "        tol = max(atol, rtol * smax)\n",
    "    Singular values smaller than `tol` are considered to be zero.\n",
    "\n",
    "    Return value\n",
    "    ------------\n",
    "    r : int\n",
    "        The estimated rank of the matrix.\n",
    "\n",
    "    See also\n",
    "    --------\n",
    "    numpy.linalg.matrix_rank\n",
    "        matrix_rank is basically the same as this function, but it does not\n",
    "        provide the option of the absolute tolerance.\n",
    "    \"\"\"\n",
    "\n",
    "    A = np.atleast_2d(A)\n",
    "    s = svd(A, compute_uv=False)\n",
    "    tol = max(atol, rtol * s[0])\n",
    "    rank = int((s >= tol).sum())\n",
    "    return rank\n",
    "\n",
    "\n",
    "def nullspace(A, atol=1e-13, rtol=0):\n",
    "    \"\"\"Compute an approximate basis for the nullspace of A.\n",
    "\n",
    "    The algorithm used by this function is based on the singular value\n",
    "    decomposition of `A`.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    A : ndarray\n",
    "        A should be at most 2-D.  A 1-D array with length k will be treated\n",
    "        as a 2-D with shape (1, k)\n",
    "    atol : float\n",
    "        The absolute tolerance for a zero singular value.  Singular values\n",
    "        smaller than `atol` are considered to be zero.\n",
    "    rtol : float\n",
    "        The relative tolerance.  Singular values less than rtol*smax are\n",
    "        considered to be zero, where smax is the largest singular value.\n",
    "\n",
    "    If both `atol` and `rtol` are positive, the combined tolerance is the\n",
    "    maximum of the two; that is::\n",
    "        tol = max(atol, rtol * smax)\n",
    "    Singular values smaller than `tol` are considered to be zero.\n",
    "\n",
    "    Return value\n",
    "    ------------\n",
    "    ns : ndarray\n",
    "        If `A` is an array with shape (m, k), then `ns` will be an array\n",
    "        with shape (k, n), where n is the estimated dimension of the\n",
    "        nullspace of `A`.  The columns of `ns` are a basis for the\n",
    "        nullspace; each element in numpy.dot(A, ns) will be approximately\n",
    "        zero.\n",
    "    \"\"\"\n",
    "\n",
    "    A = np.atleast_2d(A)\n",
    "    u, s, vh = svd(A)\n",
    "    tol = max(atol, rtol * s[0])\n",
    "    nnz = (s >= tol).sum()\n",
    "    ns = vh[nnz:].conj().T\n",
    "    return scipy.linalg.orth(ns)\n",
    "\n",
    "\n",
    "def elem_sympoly(lambda_=None,k=None):\n",
    "    # given a vector of lambdas and a maximum size k, determine the value of\n",
    "    # the elementary symmetric polynomials:\n",
    "    #   E(l+1,n+1) = sum_{J \\subseteq 1..n,|J| = l} prod_{i \\in J} lambda(i)\n",
    "    \n",
    "    N= len (lambda_)\n",
    "    # elem_sympoly.m:6\n",
    "    E= np.zeros((k + 1,N + 1))\n",
    "    # elem_sympoly.m:7\n",
    "    E[0,:]=1\n",
    "     # elem_sympoly.m:8\n",
    "    for l_ in range(k):\n",
    "        for n_ in range(N):\n",
    "            l = l_+1\n",
    "            n = n_+1\n",
    "            E[l,n]=E[l,n - 1] + lambda_[n - 1]*E[l-1, n-1]\n",
    "    return E\n",
    "\n",
    "\n",
    "    \n",
    "def sample_k(lambda_=None,k=None):\n",
    "\n",
    "    # pick k lambdas according to p(S) \\propto prod(lambda \\in S)\n",
    "    \n",
    "    # compute elementary symmetric polynomials\n",
    "    E=elem_sympoly(lambda_,k)\n",
    "    \n",
    "    i=len(lambda_)\n",
    "    remaining=k\n",
    "    S = []\n",
    "    while remaining > 0:\n",
    "\n",
    "        # compute marginal of i given that we choose remaining values from 1:i\n",
    "        if i == remaining:\n",
    "            marg=1\n",
    "        else:\n",
    "            marg=np.dot(lambda_[i-1],E[remaining-1,i-1]) / E[remaining,i]\n",
    "        # sample marginal\n",
    "        if np.random.uniform() < marg:\n",
    "            S.append(i-1) # indexing problems\n",
    "            remaining=remaining - 1\n",
    "        i = i - 1\n",
    "    return S\n",
    "\n",
    "\n",
    "def decompose_kernel(M):\n",
    "    \"\"\"\n",
    "    Decomposes the kernel so that dpp function can sample. \n",
    "    based on: https://github.com/javiergonzalezh/dpp/blob/master/dpp/samplers/decompose_kernel.m\n",
    "    \n",
    "    \n",
    "    this function returns:\n",
    "        * M - the original kernel\n",
    "        * V - eigenvectors\n",
    "        * D - diagonals of eigenvalues\n",
    "    \"\"\"\n",
    "    L = {}    \n",
    "    D, V  = np.linalg.eig(M)\n",
    "    L['M'] = M.copy()\n",
    "    L['V'] = np.real(V.copy())\n",
    "    L['D'] = np.real(D.copy())\n",
    "    return L\n",
    "\n",
    "\n",
    "def sample_dpp(L=None,k=None):\n",
    "    \"\"\"\n",
    "    Sample a set from a dpp. L is the (decomposed) kernel, and k is (optionally) \n",
    "    the size of the set to return \n",
    "    \"\"\"    \n",
    "    if k == L['V'].shape[1]: \n",
    "        # error handling\n",
    "        return list(range(k))\n",
    "    if k is None:\n",
    "        # choose eigenvectors randomly\n",
    "        D = np.divide(L['D'], (1+L['D']))\n",
    "        # check this - might just do a random.sample along axis.\n",
    "        v = np.random.randint(0, L['V'].shape[1], random.choice(range(L['V'].shape[1])))\n",
    "        v = np.argwhere(np.random.uniform(size=(len(D), 1) <= D))\n",
    "    else:\n",
    "        v = sample_k(L['D'], k)\n",
    "    \n",
    "    k = len(v)    \n",
    "    V = L['V'][:, v]    \n",
    "\n",
    "    # iterate\n",
    "    y_index = list(range(L['V'].shape[1]))\n",
    "    Y=[]\n",
    "    \n",
    "    for _ in range(k):\n",
    "        # compute probabilities for each item\n",
    "        P=np.sum(np.power(V, 2), 1)\n",
    "        # sample_dpp.m:21\n",
    "        # sample_dpp.m:22\n",
    "        #find(rand <= cumsum(P),1)   \n",
    "        \n",
    "        \"\"\"\n",
    "        This step differs from matlab code in the following way:\n",
    "            \n",
    "        1.  Create vector with valid indices which can be sampled\n",
    "        2.  Normalise the probabilities\n",
    "        3.  Make use of `np.random.choice` to choose (guarentees that it will be a new choice for `Y`)\n",
    "        \"\"\"        \n",
    "        P_index = [(indx, prob) for indx, prob in list(zip(range(len(P)), P)) if indx not in Y]\n",
    "        P_list = [x for x, _ in P_index]\n",
    "        P_norm = np.array([p for _, p in P_index])\n",
    "        P_norm = P_norm/np.sum(P_norm)        \n",
    "        choose_item = np.random.choice(range(len(P_list)) , 1, p=P_norm)[0]\n",
    "        \n",
    "        # add the index into our sampler\n",
    "        Y.append(y_index[choose_item])\n",
    "        if len(Y) == k:\n",
    "            return Y\n",
    "        \n",
    "        # delete item from y_index...\n",
    "        y_index.pop(choose_item)\n",
    "\n",
    "        # update...choose a vector to elinate, lets pick randomly\n",
    "        j = random.choice(range(V.shape[1]))\n",
    "        Vj = V[:, j]\n",
    "        V = np.delete(V, j, axis=1)\n",
    "        \n",
    "        # make sure we do a projection onto Vj, \n",
    "        # https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4524741/\n",
    "        \"\"\"\n",
    "        where Proj⊥Bi Bj is the the projection of Bj on the subspace perpendicular to Bi. \n",
    "        For Bi ≠ 0 and Bj = 0 the projection is ∥Proj⊥Bi Bj∥2 = 0.\n",
    "        \"\"\"\n",
    "        # is orthogonal basis\n",
    "        Vj_basis = nullspace(Vj)\n",
    "        # project onto basis\n",
    "        V = np.apply_along_axis(lambda x: np.dot(x, Vj_basis), 0, V)\n",
    "        V = scipy.linalg.orth(V)\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    # try to do a sample of stuff...\n",
    "    from sklearn.metrics.pairwise import rbf_kernel\n",
    "    from sklearn import datasets\n",
    "    iris = datasets.load_iris()\n",
    "    M = rbf_kernel(iris.data.T)\n",
    "    L = decompose_kernel(M)\n",
    "    indx = sample_dpp(L=L, k=3)\n",
    "    print(indx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import GPy\n",
    "import matplotlib.pyplot as plt\n",
    "from GPyOpt.util.general import multigrid\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Genetate grid\n",
    "Ngrid = 50\n",
    "bounds = [(-2,2),(-2,2)]\n",
    "X = multigrid(bounds, Ngrid)  \n",
    "\n",
    "# Define kernel and kernel matrix\n",
    "kernel = GPy.kern.RBF(len(bounds), variance=1, lengthscale=0.5) \n",
    "L = kernel.K(X)\n",
    "\n",
    "# Number of points of each DPP sample\n",
    "k = 50\n",
    "\n",
    "# Putative inputs\n",
    "set = [25,900, 1655,2125]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x25045b83358>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEICAYAAABWJCMKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGwpJREFUeJzt3X2QXNV55/HvTxKq7Bo5AiSD0KsVy2yAxDZMZBXLJpBg\nG6nwyvYaF9iJjRNKITGbpGI7ISFrx9SSwL44ASNb1mJiWPNiVxkZFSsFI1sOJClRzDhCSMgvE1kT\nSchokIQAQywN8+wf905oRt2a7unb997u+/tUdU133zP3nDnd00/f59x7jiICMzOrnilFN8DMzIrh\nAGBmVlEOAGZmFeUAYGZWUQ4AZmYV5QBgZlZRDgDW9ST9uaSvTPJ3z5K0VdILkn4v67ZNhqSQ9Kai\n22G9zwHAJk3ShZL+UdIRSYck/YOkXyq6XS36I2BzRMyIiFvzrlzSdyRd3aF9L0qDyYvp7RlJD0p6\nx7hyuyW9XFPmy5JOrmnfv6bbnpV0v6Q5nWiv5c8BwCZF0uuBB4HPAacCc4HPAD8tsl2TsBDY0Wij\npKk5tqVTZkbEycBbgIeBdZKuGlfm3WmZ84A+4M9qtl2bbnszMBP4q8432fLgAGCT9WaAiLg3Il6J\niJcj4psRsQ1A0s9J+rakg+k3x7slzRz75fRb5yclbZP0E0lfknS6pI1pOmaTpFPSsmPfZFdJelrS\nfkmfaNQwScvSI5PnJD0h6aIG5b4NXAzcln7DfXP67fcLkjZI+glwsaSflXSXpGFJQ5L+TNKUdB9X\npUc+f5XWt0vSBenzeyQdkPSRBvXfCPynmvpvq9l8iaQfpvtcLUk1v/ebknZKOizpIUkLm3nBIuLH\nEXEL8OfAzWN/w7gy+4CNwLl1th0Cvl5vm3UnBwCbrB8Ar0i6U9LysQ/rGgL+EjgT+HlgPskHT63/\nAryDJJi8m+SD50+B2STvzfE5+YuBJcA7gT+WdMn4RkmaC/w/4L+THJl8Avi6pNnjy0bErwKPkn7D\njYgfpJs+CNwIzAD+nuQo52eBxcCvAB8GPlqzq7cD24DTgHuA+4BfAt4E/DrJB/zJdeq/flz919Zs\nvizdxy8CHwDelf59K9M+el/aT48C947f9wTuB94AnDV+g6T5wArgn+psm0Xymh23zbqTA4BNSkQ8\nD1wIBPB/gGFJ6yWdnm4fjIiHI+KnETEMfJbkw7PW5yLimfRb56PAYxHxTxHxr8A64G3jyn8mIn4S\nEU8CfwNcWadpvw5siIgNETEaEQ8D/SQfas16ICL+ISJGgWPAFcCfRMQLEbEb+N/Ab9SU/1FE/E1E\nvAJ8lSTY3ZD+7d8EjpIEg1bcFBHPRcS/AJuBt6bPXwP8ZUTsjIgR4C+AtzZ7FJB6Ov15as1z35D0\nHEnA+7t0v2NuTbc9AewH/rDFv8VKygHAJi39ELoqIuaRpAXOBP4aIE3n3Cdpn6Tnga8As8bt4pma\n+y/XeTz+W/OemvtDaX3jLQQuT1Mnz6UfXBcCrQxc1tYzCzgpra+27rk1j8e3m4iY6G+ZyI9r7r9U\n8/sLgVtq/rZDJEdbc2neWNlDNc+9JyJmRsTCiPjdiHi5ZtvvpdvmRsSH0oBuPcABwDIREd8Dvsyr\n+eG/IDk6+IWIeD3JN3PV/+2mza+5v4BXv8nW2gP83/QDa+z2uoi4qYV6aqfIfZbkKKD2G/YCYF8L\n+2u2rmbsAX573N/37yLiH1vYx3uBA8D3W6zbeowDgE2KpP8g6eOS5qWP55OkZLakRWYALwJH0rz8\nJzOo9r9J+veSziHJwX+1TpmvAO+W9C5JUyX9jKSLxtrZqjSt8zXgRkkz0lTLH6b1ZOEZkrGFZq0B\n/iTtA9IB6sub+cX0qOxa4NMkKa3RlltrPcUBwCbrBZLBz8fSs2W2ANuBj6fbP0NySuERkkHZ+zOo\n8++AQeBbwP9K8+uvERF7gLGB0mGSb8yfpL33+n8FfgLsIsmR3wPc0cb+at0CvD89o2fC6xAiYh1w\nM3BfmlrbDiyf4NeeS1+jJ0nGQi6PiKzab11MXhDGyk7SIuBHwEnpwKeZZcBHAGZmFeUAYGZWUU4B\nmZlVlI8AzMwqalrRDTiRWbNmxaJFi4puhplZ1xgYGHg2Io6b+qSeUgeARYsW0d/fX3QzzMy6hqSh\niUslnAIyM6soBwAzs4pyADAzqygHADOzinIAMDOrqLYDgKT5kjZLekrSDkm/X6eMJN0qaTBdAvC8\ndus1M7P2ZHEEMAJ8PCLOBpYBH5N09rgyy0mW8lsCrAK+kEG9hRoYOszqzYMMDB3uiXrMrHravg4g\nIvaTLBNHRLwgaSfJikNP1RRbCdwVybwTWyTNlDQn/d2uMzB0mA/dvoWjI6NMnzaFu69exvkLxy+J\n2z31mFk1ZToGkE7b+zbgsXGb5vLaZfb20mAJO0mrJPVL6h8eLufKc1t2HeToyCijAcdGRtmy62BX\n12Nm1ZRZAJB0MvB14A/SBcMnJSLWRkRfRPTNnt3U1cy5W7b4NKZPm8JUwUnTprBs8WldXY+ZVVMm\nU0FIOonkw//uiKi38tM+Xrue6zyyW1M1d+cvPIW7r17Gll0HWbb4tI6lZfKqx3rbwNBhv4esrrYD\ngCQBXwJ2RsRnGxRbD1wr6T6SZQSPdGv+f8z5C0/J5Z8pr3qsN3kcyU4kiyOA/wj8BvCkpK3pc38K\nLACIiDXABpK1SAeBl0gW9DazDqs3juQAYGOyOAvo7wFNUCaAj7Vbl5m1Zmwc6djIqMeR7Dilng7a\nsuVccPV4HMlOxAGgIpwLri6PI1kjnguoInxNgZmN5wBQEb6mwMzGcwqoIpwL7i4er7E8OABUiHPB\n3cHjNZYXp4DMSsbjNZYXBwCzkvF4jeXFKSCzkvF4TTWUYZzHAcCshDxe09vKMs7jFJCZWc7KMs7j\nAGBmlrOyjPM4BWRmlrOyjPM4AJiZFaAM4zxOAZmZVZQDgJlZRWUSACTdIemApO0Ntl8k6Yikrent\nU1nUa2YTGxg6zOrNgwwMHS66KVYyWY0BfBm4DbjrBGUejYjLMqrPzJpQlvPNrZwyOQKIiEeAQ1ns\ny8yyU5bzza2c8hwDuEDSNkkbJZ3TqJCkVZL6JfUPDw/n2Dyz3lOW882tnJSs157BjqRFwIMRcW6d\nba8HRiPiRUkrgFsiYslE++zr64v+/v5M2mdWVWWYc8byI2kgIvqaKZvLdQAR8XzN/Q2SPi9pVkQ8\nm0f9ZlVWhvPNrZxySQFJOkOS0vtL03qdjDQzK1AmRwCS7gUuAmZJ2gt8GjgJICLWAO8HfkfSCPAy\ncEVklXsyM7NJySQARMSVE2y/jeQ0UTMzKwlfCWxmVlEOAGZmFeUAYGZWUQ4AZmYV5QBgZlZRDgBm\nZhXlAGBmVlEOAGZmFeUAYGZWUQ4AZmYV5QBgZlZRDgBmZhXlAGBmVlEOAGZmFeUAYGZWUQ4AZmYV\nlUkAkHSHpAOStjfYLkm3ShqUtE3SeVnUa91jYOgwqzcPMjB0uOimmFkqq0Xhv0yy4tddDbYvB5ak\nt7cDX0h/WgUMDB3mQ7dv4ejIKNOnTeHuq5d5kXKzEsjkCCAiHgEOnaDISuCuSGwBZkqak0XdVn5b\ndh3k6MgoowHHRkbZsutg0U0yM/IbA5gL7Kl5vDd97jiSVknql9Q/PDycS+Oss5YtPo3p06YwVXDS\ntCksW3xa0U0yM7JLAWUmItYCawH6+vqi4OZYBs5feAp3X72MLbsOsmzxaU2nfwaGDrf8O2bWvLwC\nwD5gfs3jeelzVhHnLzylpQ9xjxuYdV5eKaD1wIfTs4GWAUciYn9OdVsX8riBWedlcgQg6V7gImCW\npL3Ap4GTACJiDbABWAEMAi8BH82iXutdY+MGx0ZGPW5g1iGKKG+ava+vL/r7+4tuhhXEYwBmrZM0\nEBF9zZQt3SCw2ZhWxw3MrDWeCsLMrKIcAMzMKsoBwMysohwAzMwqygHAzKyiHADMzCrKAaBLeX59\n62V+f+fD1wF0Ic+TY73M7+/8+AigC3meHOtlfn/nxwGgC3l+fetlfn/nx3MBdSnPk2OdVPT7q+j6\nu5nnAqoAz5NjnVKGHLzf3/lwCsjMXsM5+OpwADCz13AOvjqyWhDmUuAWYCpwe0TcNG77RcADwI/S\np+6PiBuyqNvMsjXZNZyroNfGJtoOAJKmAquBdwB7gcclrY+Ip8YVfTQiLmu3PjPrPOfgj1eGsZGs\nZZECWgoMRsSuiDgK3AeszGC/Zmal0YtjI1kEgLnAnprHe9PnxrtA0jZJGyWd02hnklZJ6pfUPzw8\nnEHzzMza14tjI3mdBvpdYEFEvChpBfANYEm9ghGxFlgLyXUAObXPzLpAkTn4XhwbySIA7APm1zye\nlz73byLi+Zr7GyR9XtKsiHg2g/rNrALKkIPvtbGRLFJAjwNLJL1R0nTgCmB9bQFJZ0hSen9pWm/3\nJ9DMLDe9mIMvWttHABExIula4CGS00DviIgdkq5Jt68B3g/8jqQR4GXgiijzHBRmVjpjOfhjI6M9\nk4MvmucCMrOu0Wvn4XeC5wIys57Uazn4onkqCDOzinIAMDOrKAcAsxx4jVsrI48BmHVYGc5fN6vH\nRwBmHebz162sHADMOqwX55Cx3uAUkFmHlXkOGZ9XX20OAGY5KOP56x6bMKeAzCrKYxPmAGBWUR6b\nMKeAzCrq/IWn8KnLzmHj9v0sP3eO0z8V5ABgVlEDQ4e54cEdHB0Z5fHdhzjrjBkOAhXjFJBZRXkM\nwBwAzCrKYwDWkykgn9tsNrEyX59g+cgkAEi6FLiFZEWw2yPipnHblW5fAbwEXBUR382i7vF8brNZ\n88p4fYLlp+0UkKSpwGpgOXA2cKWks8cVWw4sSW+rgC+0W28jzmuamTUnizGApcBgROyKiKPAfcDK\ncWVWAndFYgswU9KcDOo+jvOaZmbNySIFNBfYU/N4L/D2JsrMBfaP35mkVSRHCSxYsKDlxjivaWbW\nnNINAkfEWmAtJIvCT2YfzmuamU0sixTQPmB+zeN56XOtljEzsxxlEQAeB5ZIeqOk6cAVwPpxZdYD\nH1ZiGXAkIo5L/5iZWX7aTgFFxIika4GHSE4DvSMidki6Jt2+BthAcgroIMlpoB9tt14zM2tPJmMA\nEbGB5EO+9rk1NfcD+FgWdZmZWTY8FYSZWUU5AJiZVZQDgJl1jYGhw6zePMjA0OGim9ITSncdgJlZ\nPZ7nK3s+AjCzruB5vrLnAGBmXcHzfGXPKSAz6wqe5yt7DgBm1jU8z1e2nAIyM6soBwAzs4pyALBc\n+Pzt7uLXqxo8BmAd5/O3u4tfr+rwEYB1nM/f7i5+varDAcA6zudvdxe/XtWhZKbmcurr64v+/v6i\nm2EZGBg67PO3u4hfr+4laSAi+pop6zEAy4XP3+4ufr2qoa0AIOlU4KvAImA38IGIOO60AUm7gReA\nV4CRZqOTmZl1TrtjANcB34qIJcC30seNXBwRb/WHv5lZObQbAFYCd6b37wTe0+b+CuFzns2sitod\nAzg9Ivan938MnN6gXACbJL0CfDEi1jbaoaRVwCqABQsWtNm8ifmcZzOrqgkDgKRNwBl1Nl1f+yAi\nQlKjU4oujIh9kt4APCzpexHxSL2CaXBYC8lZQBO1r131znl2ADCzKpgwAETEJY22SXpG0pyI2C9p\nDnCgwT72pT8PSFoHLAXqBoC8jZ3zfGxk1Oc8m1mltJsCWg98BLgp/fnA+AKSXgdMiYgX0vvvBG5o\ns97MeI5xM6uqdgPATcDXJP0WMAR8AEDSmcDtEbGCZFxgnaSx+u6JiL9ts95M+ZxnM6uitgJARBwE\nfq3O808DK9L7u4C3tFOPmZllz3MBmZlVlAOA9ZRWr+mYzDUgvm6kNe6v8vJcQNYzWr2mYzLXgPi6\nkda4v8rNRwDWM1qdx34y8957rvzWuL/KzQHAekar89hPZt57z5XfGvdXuXk9AOsprc5jP5l573tp\nrvw8/pZe6q9u0Mp6AA4AZhXl/HxvaiUAOAVkVlHOz5sDgFlFOT9vPg3UKsF56ON5HqxileE96QBg\nPc+57sY8D1YxyvKedArIep5z3VY2ZXlPOgBYz3Ou28qmLO9JnwZqlVCGfKtZrU69J1s5DdRjAFYJ\nznVb2ZThPdlWCkjS5ZJ2SBqV1DDiSLpU0vclDUq6rp06zcwsG+2OAWwH3scJ1veVNBVYDSwHzgau\nlHR2m/WamVmb2l0RbCdAutxjI0uBwXRlMCTdB6wEnmqnbjMza08eZwHNBfbUPN6bPleXpFWS+iX1\nDw8PZ9oQL0xRHPe9lY0XA2riCEDSJuCMOpuuj4gHsm5QRKwF1kJyFlBW+y3LhRdV5L63svFiQIkJ\njwAi4pKIOLfOrdkP/33A/JrH89LnclWWCy+qyH1vZePFgBJ5pIAeB5ZIeqOk6cAVwPoc6n2Nslx4\nUUXueysbLwaUaOtCMEnvBT4HzAaeA7ZGxLsknQncHhEr0nIrgL8GpgJ3RMSNzew/6wvBfDFQccra\n92Vtl3Very4G5AVhzJrQizldMy8IY9aEXszpmrXCAcAqqxdzumat8FxAVlleEMWqzgHAKq0ME3KZ\nFcUpIDOzinIAMDOrKAeACum1eUzM2lH0/0PR9YPHACrD57ybvaro/4ei6x/jI4CK8DnvZq8q+v+h\n6PrHOABUhM95N3tV0f8PRdc/xlNBVEg3zGPSKVX+23tJlq9j0e+JMiwK7wBgPa8s+VZrj1/H5ngu\nILMaZcm3Wnv8OmbPAcB6XlnyrdYev47ZcwrIctGr+VbLl1/HibWSAvJ1ANZxZcjdes6f3uDXMVtt\npYAkXS5ph6RRSQ0jjqTdkp6UtFWSv9JXjHO3ZuXU7hHAduB9wBebKHtxRDzbZn3WhcZyt8dGRp27\nNSuRtgJAROwEkJRNa6wned79amiUn3fevjV59ldeYwABbJL0CvDFiFjbqKCkVcAqgAULFuTUPOs0\n5257W6NxnjKM/3STvPtrwjEASZskba9zW9lCPRdGxFuB5cDHJP1yo4IRsTYi+iKib/bs2S1UYWZF\naTTO4/Gf1uTdXxMeAUTEJe1WEhH70p8HJK0DlgKPtLtfMyuHRuM8Hv9pTd79lcl1AJK+A3wiIo47\nw0fS64ApEfFCev9h4IaI+NuJ9uvrAMy6R5ZjAFUeN2j3b89tLiBJ7wU+B8wGngO2RsS7JJ0J3B4R\nKyQtBtalvzINuCcibmxm/w4AZtXjcYP25HYhWESs49UP99rnnwZWpPd3AW9ppx4zq456eXAHgM7w\nXEBmViqe8yc/ngrCrEBVznU34utG8uMAYFYQ57ob83Uj+XAKyKwgPkfeiuYAYFYQ57qtaE4BmRXE\nuW4rmgOAWYGc67YiOQVkZlZRDgBmZhXlAGBmhRgYOszqzYMMDB0uuimvUdZ2dYLHAMwsd2W9BqKs\n7eoUHwGYWe7Keg1EWdvVKQ4AZpa7sl4DUdZ2dUom6wF0SlWmg/Z8MFY2ebwny/q+L2u7mpXbdNDW\nvqrlHK388npPlvUaiLK2qxOcAipY1XKOVn5+T1ZHWwFA0v+U9D1J2yStkzSzQblLJX1f0qCk69qp\ns9dULedo5ef3ZHW0uyTkO4FvR8SIpJsBIuKPx5WZCvwAeAewF3gcuDIinppo/1UZA7jnsX9h4/b9\nLD93Dh98+4Kim2PW9XnwKstzSchv1jzcAry/TrGlwGC6NCSS7gNWAhMGgCoYGDrMDQ/u4OjIKI/v\nPsRZZ8zwP5wVrkp58CrLcgzgN4GNdZ6fC+ypebw3fa4uSask9UvqHx4ezrB55eR8q5kVZcIAIGmT\npO11bitrylwPjAB3t9ugiFgbEX0R0Td79ux2d1d6zreaWVEmTAFFxCUn2i7pKuAy4Nei/oDCPmB+\nzeN56XOG54Q3s+K0NQYg6VLgj4BfiYiXGhR7HFgi6Y0kH/xXAB9sp95e43yrmRWh3TGA24AZwMOS\ntkpaAyDpTEkbACJiBLgWeAjYCXwtIna0Wa+ZmbWp3bOA3tTg+aeBFTWPNwAb2qnLzMyy5SuBzcwq\nygHAzKyiHADMzCrKAcDMrKJKvR6ApGFg6ARFZgHP5tScVrltk+O2TU5Z21bWdkHvtm1hRDR1FW2p\nA8BEJPU3O+lR3ty2yXHbJqesbStru8BtA6eAzMwqywHAzKyiuj0ArC26ASfgtk2O2zY5ZW1bWdsF\nblt3jwGYmdnkdfsRgJmZTZIDgJlZRXVVACjzIvSSLpe0Q9KopIanb0naLenJdPbUXBY8bqFtRfTb\nqZIelvTD9GfdebHz6reJ+kCJW9Pt2ySd16m2TKJtF0k6kvbRVkmfyrFtd0g6IGl7g+2F9FsT7Sqy\nz+ZL2izpqfT/8/frlOlsv0VE19yAdwLT0vs3AzfXKTMV+GdgMTAdeAI4O4e2/TxwFvAdoO8E5XYD\ns3LutwnbVmC//Q/guvT+dfVe07z6rZk+IJnldiMgYBnwWE6vYTNtuwh4MM/3Vk3dvwycB2xvsL2o\nfpuoXUX22RzgvPT+DOAHeb/fuuoIICK+Gcn6ApAsQj+vTrF/W4Q+Io4CY4vQd7ptOyPi+52uZzKa\nbFsh/ZbWcWd6/07gPTnU2UgzfbASuCsSW4CZkuaUpG2FiYhHgEMnKFJIvzXRrsJExP6I+G56/wWS\n9VLGr5fe0X7rqgAwTiaL0BcggE2SBiStKroxNYrqt9MjYn96/8fA6Q3K5dFvzfRBUf3UbL0XpKmC\njZLOyaFdzSrz/2XhfSZpEfA24LFxmzrab20tCNMJkjYBZ9TZdH1EPJCWyWwR+qzb1oQLI2KfpDeQ\nrKT2vfRbShna1hEnalvtg4gISY3OS+5Iv/WY7wILIuJFSSuAbwBLCm5T2RXeZ5JOBr4O/EFEPJ9n\n3aULAFHiRegnaluT+9iX/jwgaR3JoX3bH2QZtK2QfpP0jKQ5EbE/PbQ90GAfHem3cZrpg4710wQm\nrLf2wyMiNkj6vKRZEVGGCc+K6rcTKrrPJJ1E8uF/d0TcX6dIR/utq1JAenUR+v8cTSxCL2k6ySL0\n6/Nq44lIep2kGWP3SQa1656dUICi+m098JH0/keA445Wcuy3ZvpgPfDh9OyMZcCRmhRWJ03YNkln\nSFJ6fynJ//fBHNrWjKL67YSK7LO03i8BOyPisw2Kdbbfihj9nuwNGCTJh21Nb2vS588ENtSUW0Ey\nov7PJCmQPNr2XpL83E+BZ4CHxreN5AyOJ9LbjjK1rcB+Ow34FvBDYBNwapH9Vq8PgGuAa9L7Alan\n25/kBGd8FdC2a9P+eYLkJIkLcmzbvcB+4Fj6XvutMvRbE+0qss8uJBnb2lbzmbYiz37zVBBmZhXV\nVSkgMzPLjgOAmVlFOQCYmVWUA4CZWUU5AJiZVZQDgJlZRTkAmJlV1P8HFJHb61TleK4AAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x2504121d4e0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "# Samples and plot from original and conditional with standard DPPS\n",
    "sample          = sample_dpp(decompose_kernel(L),k)\n",
    "\n",
    "plt.plot(X[sample,0],X[sample,1],'.',)\n",
    "plt.title('Sample from the DPP')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x25042dfb550>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEICAYAAABWJCMKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGoVJREFUeJzt3X+QHOV95/H3RxLYCWBLtmQhJFaCRHECccWGjaIjXE6u\nYBtUvlPsMnf4fGfAVinkjPPjnB9cnLOd1PkK5+5yCUEBqzCFFcsQXyXCqlgEIwxHEpcoVpwsJARB\n1qFIQkZCFgKML2it7/3Rz+JhNbM7s9M73TPP51U1td3Tz/bzzPP0znf7eZ7uVkRgZmb5mVF1AczM\nrBoOAGZmmXIAMDPLlAOAmVmmHADMzDLlAGBmlikHAKsdSZ+R9KUp/u5bJW2X9KKkXy27bFMhKST9\neNXlaIekByWtrroc1hsOAPYqSZdK+qak45K+K+nvJf1s1eXq0G8DD0TEWRFxU68z9xeo9RMHAANA\n0huAvwb+FHgTsBD4feCfqizXFCwGdrXaKGlmD8syrSTNqroM1t8cAGzMTwBExJ0R8YOI+H5EfD0i\ndgBI+jFJ35B0VNJzkjZImj32y5KelvRbknZI+p6kL0iaL+me1B2zRdKclHZJ6hZZI+kZSYck/War\ngklans5Mnpf0LUkrWqT7BvBO4GZJL0n6CUl3SLpF0mZJ3wPeKemNktZLOiJpn6TfkzQj7eOadObz\nP1N+eyVdkt7fL+mwpKtb5P9Z4J835H9zw+bLJD2V9rlWkhp+7yOSdks6JuleSYtb7H+s3j4q6R+B\nb6T3/5ek76Qzt4ckXdjwO3ek/L6W2uFhST/WsP1dkp5Iv3sz0FiuGalu9qXPvV7SG8eV5dpUL8ck\nXSfpZ9Mx8Py4z291FBF++QXwBuAo8EXgCmDOuO0/DrwLeB0wD3gI+OOG7U8DW4H5FGcPh4FHgXcA\nr6f4svp0SrsECOBO4AzgbcAR4LK0/TPAl9LywlSulRT/sLwrrc9r8TkeBFY3rN8BHAd+Pv3+64H1\nwFeBs1JZ/gH4aEp/DTAKXAvMBP4L8I/A2vTZ3w28CJzZTv7pvaA4u5oNDKXPennatgrYA/wUMAv4\nPeCbLfY9Vm/rU739SHr/I+mzvA74Y2D7uM9/FFiW9r8BuCttm5s+yweA04DfSJ99dcN+9wDnA2cC\nfwX8+biy3Jrq9N3A/wPuBt7ScAz8i6qPbb8m+LuvugB+1eeVvoTuAA6kL4JNwPwWaX8J+D8N608D\nH2pY/0vglob1jwN3p+WxL4+fbNj+h8AX0vJn+GEA+J2xL52GtPcCV7co12u+gNPnWd+wPhN4Bbig\n4b1fBh5My9cATzVse1sq6/yG944Cb28n//ReAJc2rH8FuCEt30MKPml9BvAysLjJvsfq7fwJ2nB2\nSvPGhs9/W8P2lcATafnDwNaGbUptPxYA7gf+Q8P2twInKALJWFkWjquXfzPuGPj1qo9rv1q/3AVk\nr4qI3RFxTUQsAn4aOIfiP0pSd85dkg5KegH4EsV/kI2ebVj+fpP1M8el39+wvC/lN95i4MrUpfC8\npOeBS4EFHXy0xnzmUvy3u29c3gsb1seXm4iY7LNM5jsNyy83/P5i4E8aPtt3Kb6IF9Laq59H0kxJ\nN0r6dmqXp9OmxrZplfc5jfuK4lu7sa7O4dR6mkVxljem0za3GnEAsKYi4gmK/x5/Or31Xyn+43tb\nRLwB+Hc09BdP0bkNy0PAM03S7Kc4A5jd8DojIm7sIJ/GW94+R/FfbGM/+xBwsIP9tZtXO/YDvzzu\n8/1IRHyzzTz+LUU30mXAGyn+M4f22uYQDW2QxiUa2+QZTq2nUV77JW99zAHAAJD0k5I+IWlRWj8X\n+CBFvz4UfcwvAcclLQR+q4Rs/7OkH02DltcCf9EkzZeAfynpPem/3ddLWjFWzk5FxA8oumA+K+ms\nNOD6H1M+ZXiWos+8XbcC/2ls4DYNUF/Zwe+fRTFT6yjwoxSBul1fAy6U9H4VM4p+FTi7YfudwG9I\nOk/SmWnffxERox3kYTXmAGBjXgR+Dng4zZbZCuwEPpG2/z5wEcWA6tcoBgS79b8pBhnvB/57RHx9\nfIKI2E/xH+7vUgye7qcIPt0cux8HvgfsBf4O+DJwexf7a/QnwAfSrJhJr0OIiI3A54C7UhfOTopB\n+Hatp+iaOQg8zg8D9qQi4jngSuBGigCyFPj7hiS3A39OMeD/fykGeT/eQdms5pQGa8x6RtISii+U\n0/zfpFl1fAZgZpYpBwAzs0y5C8jMLFM+AzAzy1StbyY1d+7cWLJkSdXFMDPrG9u2bXsuIua1k7bW\nAWDJkiWMjIxUXQwzs74had/kqQruAjIzy5QDgJlZphwAzMwy5QBgZpYpBwAzs0x1HQAknSvpAUmP\nS9ol6deapJGkmyTtSY+Lu6jbfM3MrDtlnAGMAp+IiAuA5cDHJF0wLs0VFHcaXAqsAW4pIV/r0LZ9\nx1j7wB627TtWdVHMrAa6vg4gIg5RPFiCiHhR0m6Kpxk93pBsFcVj+QLYKmm2pAXpd60Htu07xodu\n28oroyc5fdYMNqxezsWL51RdLDOrUKljAOk2v+8AHh63aSGvfdTcAVo88k7SGkkjkkaOHDlSZvGy\ntnXvUV4ZPcnJgBOjJ9m692jVRTKzipUWANITg8YeAv3CVPcTEesiYjgihufNa+tqZmvD8vPfzOmz\nZjBTcNqsGSw//81VF8nMKlbKrSAknUbx5b8hIpo9Keogr33W6CLKewarteHixXPYsHo5W/ceZfn5\nb3b3j5l1HwDSg6S/AOyOiD9qkWwTcL2kuygeO3jc/f+9d/HiOf7iN7NXlXEG8PPAvwcek7Q9vfe7\nwBBARNwKbAZWUjz/9WWKB4CbmVmFypgF9HeAJkkTwMe6zcvMzMrjK4HNzDLlAGA94YvQ6sntkrda\nPxDGBoMvQqsnt4v5DMCmnS9Cqye3izkA2LTzRWj15HYxFRN06ml4eDj8TODBsG3fMV+EVkNul8Ej\naVtEDLeT1mMA1hO+CK2e3C55cxeQmVmmHADMzDLlANCnPH+7ftwm1ok6HC8eA+hDnr9dP24T60Rd\njhefAfQhz9+uH7eJdaIux4sDQB/y/O36cZtYJ+pyvPg6gD7l+dv14zaxTkzX8dLJdQAOAGaWpUEN\n2L4QzMxsAnUZhK2axwDMLDt1GYStWikBQNLtkg5L2tli+wpJxyVtT69PlZGvmdlU1GUQtmpldQHd\nAdwMrJ8gzd9GxHtLys/MbMouXjyHDauXD+QYQCdKCQAR8ZCkJWXsy8ysF3wjvN6OAVwiaYekeyRd\n2MN8zcysiV7NAnoUGIqIlyStBO4GljZLKGkNsAZgaGioR8UzM8tPT84AIuKFiHgpLW8GTpM0t0Xa\ndRExHBHD8+bN60XxzKwCdbgZWu56cgYg6Wzg2YgIScsoAk+e867MzPPwa6KUACDpTmAFMFfSAeDT\nwGkAEXEr8AHgVySNAt8Hroo6X4JsZtOq2Tx8B4DeK2sW0Acn2X4zxTRRM7NX5+GfGD2Z9Tz8qvlW\nEGbWc56HXw8OAGYdGtSbiPWa5+FXzwHArAMevLRB4pvBmXXANxGzQeIAYNYB30TMBom7gMw64MFL\nGyQOAGYd8uClDQp3AZmZZcoBwMwsUw4AJfMNrqaf69isHB4DKJHniE8/17FZeXwGUCLPEZ9+rmOz\n8jgAlMhzxKef69isPKrzXZmHh4djZGSk6mJ0xPeJmX6uY7PWJG2LiOF20noMoGSeIz79XMdm5XAX\nkJlZphwAzMwy5QBgZpapUgKApNslHZa0s8V2SbpJ0h5JOyRdVEa+ZmY2dWWdAdwBXD7B9iuApem1\nBrilpHzNzGyKSgkAEfEQ8N0JkqwC1kdhKzBb0oIy8jYzs6np1RjAQmB/w/qB9N4pJK2RNCJp5MiR\nIz0pnJlZjmo3CBwR6yJiOCKG582bV+q+fROxenK7WJ3kdDz26kKwg8C5DeuL0ns945uI1ZPbxeok\nt+OxV2cAm4APp9lAy4HjEXGoR3kDvolYXbldrE5yOx5LOQOQdCewApgr6QDwaeA0gIi4FdgMrAT2\nAC8D15aRbyfGbiJ2YvSkbyJWI24Xq5Pcjsesbgbnm4jVk9ulM66v6dWr+p2ufDq5GVxWAcCs3+XW\nRz2oprMdOwkAtZsFZGat5dZHPajq0o4OAGZ9xA/EGQx1aUd3AZn1GY8BDIY6jAH4gTBmJenVF7Mf\niDMY6tCODgBmJfDgrPUjjwGYlaAug3pmnXAAMCtBXQb1zDrhLiCzEly8eA4bVi/vu8FZDyjnzQHA\nrCR1GNTrhMctzF1AZpnyuIU5AJhlyuMW5i4gs0z167iFlccBwCxj/TZuYeVyF5CZWaYcAMzMMuUA\nYGaWqVICgKTLJT0paY+kG5psXyHpuKTt6fWpMvI1M6uzbfuOsfaBPWzbd6zqojTV9SCwpJnAWuBd\nwAHgEUmbIuLxcUn/NiLe221+Zmb9oB8utCvjDGAZsCci9kbEK8BdwKoS9mtm1rf64UK7MgLAQmB/\nw/qB9N54l0jaIekeSRe22pmkNZJGJI0cOXKkhOKZmfVeP1xo16vrAB4FhiLiJUkrgbuBpc0SRsQ6\nYB0UTwTrUfnMTuEbpQ226W7ffrjQrowAcBA4t2F9UXrvVRHxQsPyZkl/JmluRDxXQv5mpeuH/lub\nul61b90vtCujC+gRYKmk8ySdDlwFbGpMIOlsSUrLy1K+9esQM0v6of/Wps7tW+j6DCAiRiVdD9wL\nzARuj4hdkq5L228FPgD8iqRR4PvAVVHnp9Fb9sb6b0+Mnqxt/61Nndu3oDp/Dw8PD8fIyEjVxbBM\neQxgsA1q+0raFhHD7aT1zeDMWqh7/611x+3rW0GYmWXLAcDMLFMOAGZmmXIAMDPLlAOAmVmmHADM\nzDLlAGBmlikHADOrnbo/SGVQ+EIwM6sV34ivd3wGYGa14hu19Y4DgJnVSj88SGVQuAvIzGqlHx6k\nMigcAMysdnyjtt5wF5CZWaYcAMzMMjWQAcBziM3MJldKAJB0uaQnJe2RdEOT7ZJ0U9q+Q9JFZeTb\nzNgc4v/x9Sf50G1bHQTMzFroOgBImgmsBa4ALgA+KOmCccmuAJam1xrglm7zbcVziM3M2lPGGcAy\nYE9E7I2IV4C7gFXj0qwC1kdhKzBb0oIS8j6F5xCbmbWnjGmgC4H9DesHgJ9rI81C4FAJ+b+G5xCb\nmbWndtcBSFpD0U3E0NDQlPaRwxzibfuOOciZWVfKCAAHgXMb1hel9zpNA0BErAPWAQwPD0cJ5Rs4\nvlmWmZWhjDGAR4Clks6TdDpwFbBpXJpNwIfTbKDlwPGIKL37Jxce6DazMnR9BhARo5KuB+4FZgK3\nR8QuSdel7bcCm4GVwB7gZeDabvPN2dhA94nRkx7oNrMpU0R9e1mGh4djZGSk6mLUkscAzKwZSdsi\nYridtLUbBLb25DDQbWbTayBvBWFmZpNzADAzy5QDgJlZphwAzMwy5QBgZpYpBwAzs0w5AJiZZcoB\nwLLgp8RZu3I6VnwhmA083zzP2pXbseIzABt4vnmetSu3Y8UBwAaenxJn7crtWPHN4AbMRDeJy/kG\ncjl/dutMvx8rvhlcpibqv8ytb3M83zzP2pXTseIuoAEyUf9lbn2bZjY5B4ABMlH/ZW59m2Y2OY8B\nDBiPAZi9Vm7HvccAMjZR/2VOfZtm4LGvyXTVBSTpTZLuk/RU+tm0ZiU9LekxSdsl+V96M+sJj31N\nrNsxgBuA+yNiKXB/Wm/lnRHx9nZPTczMuuWxr4l12wW0CliRlr8IPAj8Tpf7NJtQbn26NnUXL57D\nhtXLfby00G0AmB8Rh9Lyd4D5LdIFsEXSD4DPR8S6VjuUtAZYAzA0NNRl8WzQuE/XOuWxr9YmDQCS\ntgBnN9n0ycaViAhJraYUXRoRByW9BbhP0hMR8VCzhCk4rINiFtBk5bO8NOvT9R+32dRMGgAi4rJW\n2yQ9K2lBRByStAA43GIfB9PPw5I2AsuApgHAbCJjfbonRk+6T9esS912AW0CrgZuTD+/Oj6BpDOA\nGRHxYlp+N/AHXeZrmXKfrll5ug0ANwJfkfRRYB/wrwEknQPcFhErKcYFNkoay+/LEfE3XeZbWx6g\nnH5l9um6vSxnXQWAiDgK/GKT958BVqblvcDPdJNPv/AAZX9xe1nufC+gEvmik/7i9rLcOQCUyBed\n9Be3l+XON4MrmfuU+4vbywaNbwZXIV900l/cXpYzdwGZmWXKAcDMLFMOABT9wGsf2MO2fceqLorZ\nlPk47i91aK/sxwA8F9wGgY/j/lKX9sr+DMBzwW0Q+DjuL3Vpr+wDgOeC2yDwcdxf6tJevg4AzwW3\nweDjuL9MV3t1ch2AA4CZ2QDpJABk3wVkZpYrBwAzs0w5AJiZZcoBwMwsUw4AZmaZ6ioASLpS0i5J\nJyW1HHWWdLmkJyXtkXRDN3mamVk5uj0D2Am8H3ioVQJJM4G1wBXABcAHJV3QZb5mZtalbp8JvBsg\nPfC9lWXAnvRsYCTdBawCHu8mbzMz604vxgAWAvsb1g+k95qStEbSiKSRI0eOTHvhzMxyNekZgKQt\nwNlNNn0yIr5adoEiYh2wDoorgcvev5mZFSYNABFxWZd5HATObVhflN4zM7MK9aIL6BFgqaTzJJ0O\nXAVs6kG+ZmY2gW6ngb5P0gHgnwFfk3Rvev8cSZsBImIUuB64F9gNfCUidnVXbDMbZHV4WlYOup0F\ntBHY2OT9Z4CVDeubgc3d5GVmeajL07Jy4CuBzaxW6vK0rBw4AJhZrdTlaVk5yP6h8HXgJzmZ/dDF\ni+ewYfVy/030gANAxdzfaXaqixfP8d9BD7gLqGLu7zSzqjgAVMz9nWZWFXcBVcz9nWZWFQeAGnB/\np5lVwV1AZmaZcgAwM8uUA0BGfH+VU7lOLGceA8iErzc4levEcuczgEz4eoNTuU4sdw4AmfD1Bqdy\nnVjuFFHfpy4ODw/HyMhI1cUYGL7n0KlcJzZoJG2LiOF20noMICO+3uBUrhPLmbuAzMwy1e0jIa+U\ntEvSSUktTzkkPS3pMUnbJblPx8ysBrrtAtoJvB/4fBtp3xkRz3WZn5mZlaTbZwLvBpBUTmnMasQD\nxPUzlTZxO7bWq0HgALZI+gHw+YhY1yqhpDXAGoChoaEeFc/stXyRWP1MpU3cjhObdAxA0hZJO5u8\nVnWQz6UR8XbgCuBjkn6hVcKIWBcRwxExPG/evA6yMCuPLxKrn6m0idtxYpOeAUTEZd1mEhEH08/D\nkjYCy4CHut2v2XQZu0jsxOhJXyRWE1NpE7fjxEq5EEzSg8BvRsQpM3wknQHMiIgX0/J9wB9ExN9M\ntt+yLwRzX+Dga9XG7jseDG7HyXVyIVhXAUDS+4A/BeYBzwPbI+I9ks4BbouIlZLOBzamX5kFfDki\nPtvO/ssMAO4LHHyt2thtbznpJAB0dR1ARGyMiEUR8bqImB8R70nvPxMRK9Py3oj4mfS6sN0v/7K5\nL3DwtWpjt71Zc9lcCewbfw2+Vm3stjdrLqubweXWF5ijMscAzPpRz8YAppvvBpo3f2lbFfr9uPPd\nQK3veeDWqpDbcZfNGID1Fw/cWhVyO+4cAKyWPHBrVcjtuPMYQM31e39kN3L+7FPh+ipH1fXYbf4e\nAxgQufVHjuendbUv92OlTFUed71uR3cB1Vhu/ZE2dT5WBkOv29EBoMZy64+0qfOxMhh63Y4eA6i5\nqvsjrX/4WBkMvRwDcAAwMxsgPbsZnJmZ9S8HADOzTDkAmJllygHAzCxTDgBmZplyADAzy1Stp4FK\nOgLsmyDJXOC5HhWnUy7b1LhsU1PXstW1XDC4ZVscEfPaSVjrADAZSSPtznftNZdtaly2qalr2epa\nLnDZwF1AZmbZcgAwM8tUvweAdVUXYAIu29S4bFNT17LVtVzgsvX3GICZmU1dv58BmJnZFDkAmJll\nqq8CgKT/JukJSTskbZQ0u0W6yyU9KWmPpBt6VLYrJe2SdFJSy+lbkp6W9Jik7ZJ6cq/rDspWRb29\nSdJ9kp5KP5veAL1X9TZZHahwU9q+Q9JF01WWKZRthaTjqY62S/pUD8t2u6TDkna22F5JvbVRrirr\n7FxJD0h6PP19/lqTNNNbbxHRNy/g3cCstPw54HNN0swEvg2cD5wOfAu4oAdl+yngrcCDwPAE6Z4G\n5va43iYtW4X19ofADWn5hmZt2qt6a6cOgJXAPYCA5cDDPWrDdsq2AvjrXh5bDXn/AnARsLPF9qrq\nbbJyVVlnC4CL0vJZwD/0+njrqzOAiPh6RIym1a3AoibJlgF7ImJvRLwC3AWs6kHZdkfEk9Odz1S0\nWbZK6i3l8cW0/EXgl3qQZyvt1MEqYH0UtgKzJS2oSdkqExEPAd+dIEkl9dZGuSoTEYci4tG0/CKw\nG1g4Ltm01ltfBYBxPkIRGcdbCOxvWD/AqZVapQC2SNomaU3VhWlQVb3Nj4hDafk7wPwW6XpRb+3U\nQVX11G6+l6SugnskXdiDcrWrzn+XldeZpCXAO4CHx22a1nqbVdaOyiJpC3B2k02fjIivpjSfBEaB\nDXUrWxsujYiDkt4C3CfpifRfSh3KNi0mKlvjSkSEpFbzkqel3gbMo8BQRLwkaSVwN7C04jLVXeV1\nJulM4C+BX4+IF3qZd+0CQERcNtF2SdcA7wV+MVIn2TgHgXMb1hel96a9bG3u42D6eVjSRopT+66/\nyEooWyX1JulZSQsi4lA6tT3cYh/TUm/jtFMH01ZPk5g038Yvj4jYLOnPJM2NiDrc8KyqeptQ1XUm\n6TSKL/8NEfFXTZJMa731VReQpMuB3wb+VUS83CLZI8BSSedJOh24CtjUqzJORNIZks4aW6YY1G46\nO6ECVdXbJuDqtHw1cMrZSg/rrZ062AR8OM3OWA4cb+jCmk6Tlk3S2ZKUlpdR/H0f7UHZ2lFVvU2o\nyjpL+X4B2B0Rf9Qi2fTWWxWj31N9AXso+sO2p9et6f1zgM0N6VZSjKh/m6ILpBdlex9F/9w/Ac8C\n944vG8UMjm+l1646la3CenszcD/wFLAFeFOV9dasDoDrgOvSsoC1aftjTDDjq4KyXZ/q51sUkyQu\n6WHZ7gQOASfSsfbROtRbG+Wqss4upRjb2tHwnbayl/XmW0GYmWWqr7qAzMysPA4AZmaZcgAwM8uU\nA4CZWaYcAMzMMuUAYGaWKQcAM7NM/X80/A/p4ruFwQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x25042dc9588>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import random\n",
    "sample1 = random.sample(range(X.shape[0]), 50)\n",
    "plt.plot(X[sample1,0],X[sample1,1],'.',)\n",
    "plt.title('Sample from the random')\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}