diff --git a/DetPoisson_Python_w_DPPy.ipynb b/DetPoisson_Python_w_DPPy.ipynb new file mode 100644 index 0000000..a887311 --- /dev/null +++ b/DetPoisson_Python_w_DPPy.ipynb @@ -0,0 +1,279 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# DetPoisson [GitHub repo](https://github.com/hpaulkeeler/DetPoisson_Python) dedicated to [article](https://arxiv.org/pdf/1810.08672.pdf) from Błaszczyszyn & Keeler" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### The authors mention [DPPy](https://github.com/guilgautier/DPPy) as a way to sample DPPs\n", + "#### Here is a more Pythonic implementation that uses [DPPy](https://github.com/guilgautier/DPPy) " + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from sklearn.metrics import euclidean_distances # to compute |x-y|^2\n", + "\n", + "def kernel(name='RBF', data=np.random.rand(2,2), **kwargs):\n", + " \"\"\"Compute kernel matrix associated to data points\"\"\"\n", + "\n", + " sigma = kwargs.get('sigma', 1.0) # Scale parameter of kernel\n", + " alpha = kwargs.get('alpha', 1.0) # Exponent parameter for Cauchy\n", + " \n", + " dist_sq = euclidean_distances(data, data)\n", + " \n", + " if name == 'RBF':\n", + " return np.exp(-dist_sq/sigma**2)\n", + " \n", + " elif name == 'Cauchy':\n", + " return 1/(1+dist_sq/sigma**2)**(alpha+1/2)\n", + " \n", + " else:\n", + " raise ValueError('Invalid kernel name, choose `RBF` or `Cauchy`')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Clone DPPy and create access its source code" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "fatal: destination path 'DPPy' already exists and is not an empty directory.\r\n" + ] + } + ], + "source": [ + "#!rm -rf DPPy\n", + "!git clone https://github.com/guilgautier/DPPy.git\n", + "\n", + "import os\n", + "import sys\n", + "sys.path.insert(0, os.path.abspath('DPPy')) # Access DPPy source code\n", + "\n", + "from dppy.finite_dpps import FiniteDPP" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Sample from 2D Poisson Process" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAAEICAYAAADMa/SXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xl8VOXZ//HPxRLWAIpgAQmbaIyggMiiKAg8GGWtUgpWRKHi0xY3QKtFBLTVuhQfioqoRSoo4gY/cAFUNhGQTRZlUZGgQRRECJE1kOv3x4wxYZGgk5zJzPf9euXFzD1n5nxzAly5z7nPfZu7IyIiEq+KBR1AREQkSCqEIiIS11QIRUQkrqkQiohIXFMhFBGRuKZCKCIicU2FUKKKmf3NzJ6N9Lb5+Cw3szMj8VkiUrSY7iOUgmJm1wODgHrAbmAKcLe77woy17GYmQP13f3zoLOISOFSj1AKhJkNAh4C7gAqAi2AWsA7ZpZwnPeUKLyEBSMWvgeReKNCKBFnZhWAEcDN7j7D3bPcPQ3oAdQGrg1vN9zMXjWziWa2G7g+3DYx12ddZ2abzWyHmQ01szQza5/r/RPDj2uHT2/2MbMvzew7MxuS63OamdkiM9tlZlvN7PHjFeRjfD9zzexBM1tiZrvN7P+Z2alH7LefmX0JzA63dzGzT8L7m2tm5+T6vJpm9rqZbQ9/X4/neq2vma0zs51mNtPMaoXbzcweM7Nt4QxrzKxB+LUrzWytmWWa2RYzG3zSPzSROKZCKAXhIqA08HruRnf/AXgL+J9czV2BV4FKwAu5tzezFOBJ4A9ANUI9yxon2Hcr4GygHXBvrgJ0GLgdOA1oGX79zyfxPV0H9A3nOAT8+4jXWwPnAJeb2VnAJOA2oAqh73m6mSWYWXHgDWAzoV8KagAvhb/frsDfgKvC73s//DkAHYBLgbMIHYcewI7wa/8BbnL3RKAB4WIsIvmjQigF4TTgO3c/dIzXtoZf/9Eid5/q7tnuvu+IbbsD0919gbsfBO4FTnRRe4S773P3VcAq4HwAd1/u7ovd/VC4dzqWUPHKrwnu/rG77wGGAj3CRe1Hw919T/h7+D3wpru/4+5ZwKNAGUK/IDQDqgN3hLff7+4Lwp/xv8CD7r4ufOweABqFe4VZQCKQTOja/jp33xp+XxaQYmYV3H2nu684ie9LJO6pEEpB+A447TjXy6qFX//RVz/zOdVzv+7ue/mpF3Q83+R6vBcoD2BmZ5nZG2b2Tfg07APkLcgnkjvnZqDkEe/P/Xr18DY/5s4Ov14DqAlsPs4vCbWAUeHTqbuA7wEDarj7bOBx4Algm5k9HT4FDXA1cCWw2czmmVnLk/i+ROKeCqEUhEXAAUKn+HKYWXngCuC9XM0/18PbCpyR6/1lgMq/MNMYYD2hkaEVCJ2CtJN4f81cj5MI9cJyF/Tc38fXhIoaELq+F37/FkIFMek4vyR8RegUZ6VcX2XcfSGAu//b3S8AUgidIr0j3L7U3bsCVYGpwMsn8X2JxD0VQok4d88gNFhmtJmlmllJM6tN6D/odGBCPj/qVaCzmV0UHtgynJMrXrklErqF4wczSwb+dJLvv9bMUsysLHAf8Kq7Hz7Oti8DHc2snZmVJHQLyQFgIbCEUIH/p5mVM7PSZnZx+H1PAXeb2bkAZlbRzH4XfnyhmTUPf94eYD+QHb7u+Aczqxg+DbsbyD7J700krqkQSoFw94cJ9boeJfSf84eEejzt3P1APj/jE+BmQoNJtgI/ANsIFZWTNRi4BsgEngEmn+T7JwDjCZ16LQ3ccrwN3X0DoZGxown1GjsDnd39YLh4dgbOBL4k9IvB78Pvm0LolpOXwqdvPybUgwaoEM69k9Bp1x3AI+HXegNp4ff8L6HBRSKST7qhXoqM8KnVXYROb24qxP3OBSa6e0RmsRGR6KIeoUQ1M+tsZmXNrByh3uUaIC3YVCISS1QIJdp1JTT45GugPtDTdRpDRCJIp0ZFRCSuqUcoIiJxLbAJgk877TSvXbt2ULsXESmSli9f/p27Vwk6RywJrBDWrl2bZcuWBbV7EZEiycw2n3grORk6NSoiInFNhVBEROKaCqGIiMQ1FUIREYlrKoQiIhLXTlgIzWycmW0zs4+P87qZ2b/N7HMzW21mTSIfU0REpGDkp0c4Hkj9mdevIDT1VX2gP6F130RERIqEE95H6O7zw2vJHU9X4Pnw/I+LzaySmVVz960Ryigikm+HDx/mnXfeYcOGDaSkpNCuXTuKFdNVIDm+SNxQX4PQOnM/Sg+3HVUIzaw/oV4jSUlJEdi1iMhPduzYweWXXw5Ay5YtGT9+PKVKleLtt9/mlFNOCThdbDOzw4RWhykBrAP6uPven9l+obtfVFj5fk6h/prk7k+7e1N3b1qlimYIEpHIuuOOO7joootYunQpo0ePZsWKFZx//vkMGTIk6GhRI7ys2UAze8/MppnZ78zMIvDR+9y9kbs3AA4SWiT6uKKlCEJkCuEWoGau52eE20SkkK1bt45bbrmFrl27cv/997N9+/agIxUad2fy5MkMHTqUH/9fNzPuueceXnrppYDTRQczKwW8A7QGRgIvAUOBhyO8q/eBM8P7HGhmH4e/bsuV5Yfwn9XMbL6ZrQxvc4mZFTez8eHna8zs9vC2jcxscXhg5hQzOyXcPtfMHjKzJWb2qZldcjJhI1EIpwHXhUePtgAydH1QpPC9++67XHrppZx66qn06dOHL7/8kqZNm/LVV1+d+M0x4tChQyQkJORpK1WqFIcOHQooUdTpARwAurn7m+7+IqGieIOZ1YrEDsysBKFBlGvM7ALgBqA50AK40cwaH/GWa4CZ7t4IOB9YCTQCarh7A3dvCDwX3vZ54K/ufh6h07DDcn1OCXdvBtx2RPsJ5ef2iUnAIuBsM0s3s35m9r9m9mO39y3gC+Bz4BngzycTQER+PXfn1ltvZfz48QwfPpyrrrqKZ555hmuuuYYHHngg6HiFwszo2rUro0aNytP+f//3f3Tr1i2QTDt37uThhx/m6quv5tZbb2Xt2rWB5MilNfBy7sWt3X0n8B5wUr2oYyhjZiuBZcCXwH+AVsAUd9/j7j8Arx9jP0sJFeLhQEN3zyRUU+qa2WgzSwV2m1lFoJK7zwu/77/Apbk+5/Xwn8uB2icTPD+jRnud4HUH/nIyOxWRyPrmm2/49ttvufLKK/O09+7dmy5dugSUqvA9+uijtGnThmXLlnHRRRcxf/58PvvsM+bOnVvoWb799ltatWpFixYt+P3vf8/atWtp06YN48ePP+rnVIi2AXWO0V4H+PZXfva+cK8uR34uPYbvTLgU6AiMN7OR7v68mZ0PXE7oWmMP4PYTfNSB8J+HOcmBoIEtwyQikVO+fHkOHjzI7t27qVixYk771q1b42q0ZFJSEmvWrGHy5MmsX7+eHj160KNHD8qWLVvoWR5++GFSU1MZPXp0Ttsll1zCTTfdRGpqalC3dDwHLDSz/+fuC82sGKFCcwowuwD29z6h4vZPwIDfAr1zbxA+JZvu7s+Er2E2MbO3gIPu/pqZbQAmunuGme00s0vc/f3w58wjAlQIRWJAYmIinTt35q677mL06NGUKFGCnTt3MmTIEPr16xd0vEJVrlw5+vbtG3QM3nnnHcaNG5enrW3bthw4cIC0tDTq1q1b6Jnc/TMz6we8ambfARWBHUBHdz9cAPtbYWbjgSXhpmfd/aMjNmsD3GFmWcAPwHWEbsF7LlyoAe4O/9kHeMrMyhI6fXpDJHKqEIrEiCeeeIKePXtSu3ZtUlJSWLZsGX369KF///5BR4tLlSpV4ptvvsnTtm/fPjIzM0lMTAwoFbj7NDN7m9CAlD3AutzXDH/F55Y/TvtIQiNUj7m9u/+X0PW+Ix01Xae7ryQ06ObI9ja5Hn9HpK8RikjRUKlSJWbMmMHatWvZvHkzjRo1olq1akHHilt9+/Zl+PDhXHTRRZx66qkcPnyYYcOG0bp1a4K+j9rdswgNUhFUCEViTkpKCikpKUHHiAoHDhwgMzOTypUr52vgRiT16dOHTz75hHr16tG8eXM2bNhAzZo1ee211wo1h5yYJuATkZhz8OBBBg8ezG9+8xvq16/P2WefXegFyMx45JFHWLt2LTfffDNTp05l/vz5gfcG5WjqEYpIzBk4cCBpaWl8/PHHVK9enfnz53PNNddQuXJl2rRpU6hZqlWrRseOHQt1n3JyLALXSH+Rpk2b+rJlywLZt4jEroyMDGrVqsXGjRupXLlyTvu4ceOYPn06U6ZMCTDdr2dmy929adA5YolOjYpITNm6dStVq1bNUwQBGjVqxKZNmwJKJdFMhVBEYkqtWrX4/vvv2bhxY572WbNm0bjxkdNciugaoYjEmDJlynDXXXfRpUsXRo4cSXJyMlOnTmXkyJGBTLUm0U+FsAjIzs4mIyODxMRESpTQj0zkRAYNGsTpp5/Ovffey5YtW2jZsiXvvvuubiuRY9Kp0Sj3/PPPU69ePWrXrk21atX4+9//TnZ2dtCxRKKamdG7d28+/PBD0tPTeeWVVzjvvPOCjiVRSoUwik2ZMoVhw4YxefJkMjIyWLRoEW+//Tb/+Mc/go4mIhIzdPtEFGvVqhV33nlnnmV0Nm7cSIsWLfj6668pWbJkgOlEJAi6fSLy1COMYl988QVNmuSdd7ZevXocOnSIzMzMgFKJiMQWFcIo1qhRI9599908bcuWLaNChQpUqlQpoFQiIrFFQxCj2D333EO3bt0A6NChAytWrODWW29l+PDhQS3qKSISc3SNMMp98MEH/P3vf2flypXUqVOHQYMGcfXVVwcdS0QComuEkaceYZS7+OKLefvtt4OOISK5/PDDD7zxxhtkZmbSoUMHatWqFXQk+RV0fk1E5CTMmzePOnXqMHHiRN5//32aNGnC/fffH3Qs+RXUIxQRyaf9+/fz+9//nkmTJtG+fXsAtm3bRvPmzWnTpg2XXHJJwAnll1CPUEQkn2bPns3ZZ5+dUwQBqlatyp///GdefPHFAJPJr6FCKCKSTwcOHKBs2bJHtZcrV44DBw4EkEgiQYVQRCSf2rZty+LFi1m3bl1O2/79+3nmmWfyzAAlRYuuEYqI5FPFihX597//TevWrbnuuuuoXLkyEydOpFGjRiqERZh6hCIiJ6F3794sWLCAcuXKsX37dkaNGsXEiRM1yUURphvqRUSOwd2ZNWsW06ZNo2TJkvTq1YvmzZsHHUs31BcA/QojInIEd+fGG29k4MCB1K1blypVqtC9e3f++c9/Bh1NCoCuEYqIHGHevHm8//77rFixgnLlygHQt29fGjRowDXXXENSUlLACSWS1CMUETnCm2++Se/evXOKIEC1atXo1KkTM2bMCDCZFIR8FUIzSzWzDWb2uZnddYzXk8xsjpl9ZGarzezKyEcVESkcZcuWJSMj46j2jIyMY95HKEXbCQuhmRUHngCuAFKAXmaWcsRm9wAvu3tjoCfwZKSDiogUll69evHf//6Xzz77LKdtwYIFLFiwQLdJxKD8XCNsBnzu7l8AmNlLQFdgba5tHKgQflwR+DqSIUVEClNycjL//Oc/adasGa1bt2bfvn0sX76cF198kQoVKpz4A6RIyU8hrAF8let5OnDkGOLhwCwzuxkoB7TnGMysP9Af0MVmEYlqffv2pWvXrrzzzjskJCRw+eWX57lmKLEjUqNGewHj3f1fZtYSmGBmDdw9O/dG7v408DSE7iOM0L5FRApE5cqV6dmzZ9AxpIDlZ7DMFqBmrudnhNty6we8DODui4DSwGmRCCgiIlKQ8lMIlwL1zayOmSUQGgwz7YhtvgTaAZjZOYQK4fZIBhURESkIJzw16u6HzGwAMBMoDoxz90/M7D5gmbtPAwYBz5jZ7YQGzlzvQc3dJiIRsWXLFqZMmUJWVhZdunShXr16QUcSKRD5uo/Q3d9y97PcvZ67/yPcdm+4COLua939Ync/390bufusggwtIgVr4sSJNGzYkBUrVrBhwwZatGjBI488EnQskQKhSbdFJI9vv/2W5ORkFi5cyDnnnAPA119/zQUXXMCsWbNo2LBhwAnjmybdjjxNsSYieUybNo0rr7wypwgCVK9eneuuu45XXnklwGQiBUOFUETyyM7OPubaesWLFyc7O/sY7xAp2lQIRSSPzp078+abb7Jx48actu3bt/P8889z9dVXB5hMpGBoGSYRyaN69eo89NBDNG/enJ49e1KqVCkmTZpE//79ady4cdDxRCJOhVBEjnLjjTfSrl07XnnlFbKyspg1axYNGjQIOpZIgdCoURGRIkSjRiNP1whFRCSuqRCKiEhcUyEUEZG4pkIoIiJxTaNGRSSquDvz5s1j9erV1K9fnw4dOlC8ePGgY0kMU49QRKJGZmYmbdq0YcCAAXz66acMGzaMJk2a8O233wYdTWKYeoQiEjVGjBhBUlISc+bMoVixYrg7d955J7fffjsvvvhi0PEkRqkQikSZzMxMJk6cyEcffUTdunW54YYbOP3004OOVSgmTZqUUwQBzIwhQ4ZQvXp1srKyKFmyZMAJJRbp1KhIFNm6dStNmjThnXfeoUmTJmzcuJHzzjuPlStXBh2tUGRlZVGqVKk8bQkJCRw+fBit9S0FRT1CkSgyYsQIunbtyqOPPprT1qJFC2699VbmzZsXYLKCkZ6ezltvvUXp0qXp0qUL3bp1Y9SoUYwcOTJnmyeffJIOHTqQkJAQYFKJZZpiTQKVkZHBnDlzKFWqFG3btj2qNxBvzjjjDObOncuZZ56Z05aVlUWlSpXYunUrFSpUCDBdZD322GPcf//9dOrUiT179jBnzhxGjRrFgw8+SFJSEpdddhlLly5lyZIlzJ49O88xiWeaYi3y1COUwEycOJGbb76Z5s2bs3fvXm644QZefvllLr300qCjBaZs2bLs3r07T9vevXsxs5i6PrZy5UoeeeQRVq9ezRlnnAHAkiVLSE1NZd26dcycOZM1a9bQoUMHnnvuORITEwNOLLFMhVAKzZYtW5g1axZly5bl7LPP5vbbb+eDDz4gJSUFgHfffZfu3buzadMmypUrF3DaYFx77bUMGzaM1157jYSEBNydESNG0KlTJ8qUKRN0vIh56aWX6NevX04RBGjWrBmXXHIJ7733Htddd12A6STeaLCMFIqRI0fSsGFDZs2axfPPP89FF11EmzZtcoogQPv27WnWrBnTp08PMGnkpKWlMXjwYFJTU7nlllv49NNPT/iev/71ryQkJFCvXj169+5Nw4YN+eCDDxg9enQhJC48WVlZlC5d+qj2MmXKcPDgwQASSTxTIZQCt2zZMv71r3+xevVqJk2axJtvvkmnTp1466232LNnT55tK1eufNSpwaLo448/pnnz5hQrVoxbbrmFSpUqcfHFF7N48eKffV+pUqV47bXXeOONN2jfvj1jxoxh8eLFVKlSpZCSF44uXbowfvz4PD/rTZs2MXPmTK644ooAk0lccvdAvi644AKX+DBw4EAfMWJEnrZZs2Z5uXLlfNKkSTlt3377rVeuXNnT0tIKO2LEde3a1UeNGpWnbcKECX7JJZcElCi6ZGdn+4ABA7xu3bo+fPhwv+OOO7xq1ao+ZsyYoKNFPWCZB/T/dqx+adQosG/fPqZNm8b27dtp3bo1DRs2DDpSTPnLX/5C3bp1GTRoUE6bu5OUlETx4sW588472bt3L0888QR9+/Zl6NChAaaNjAoVKrBp0yYqV66c05aVlUWZMmXYt29fVA982bFjBwkJCQU+QMXdWbhwIdOmTaN06dL07NmTc845p0D3GQs0ajTy4v7U6MqVKznzzDN57rnnWLNmDampqdx4441kZ2cHHS1mdO7cmXHjxrF3796ctk2bNvHDDz9w3333sWzZMj7//HPGjx8fE0UQQqd409PT87Rt3bqVxMRESpSIzjFqS5YsoXnz5tStW5dq1arRvXt3tm3bVmD7MzMuvvhiHnroIUaMGKEiKMEJqisaDadGs7Oz/ZxzzvEXXnghp+2HH37wCy+80CdOnBhgstiSnZ3t119/vScnJ/uDDz7od999t1etWtWfeuqpoKMVmH/84x/etm1bz8jIcHf3PXv2eLdu3Xzw4MEBJzu29PR0r1Klir/wwgt+6NAh3717tw8ePNibNWvm2dnZQceTXNCp0Yh/xXWPcPXq1WRlZdGrV6+ctnLlyjFo0CAmTZoUYLLYYmaMGzeO0aNH56wi8N5773HTTTcFnKzg3HnnnZx11lnUqlWLNm3akJSURPny5bn//vuDjnZM48aNo0ePHlxzzTUUL16cxMREHn74YTIzM1m4cGHQ8UQKVHSeoykkP85raGZ52kuXLk1WVlZAqWKTmdG+fXvat29/wm137tzJY489xttvv01iYiJ9+vThuuuuO+rnFM1KlCjBmDFjGDp0KBs2bKBevXokJSUFHeu40tLSaNGiRZ42M+P8888nLS2Niy++OKBkIgUvrnuEjRs3JjMzk/feey+n7dChQzz++ON069YtwGTxa8+ePbRu3Zovv/ySxx57jNtuu43Ro0czcODAoKP9ItWrV+eyyy6L6iII0KRJE2bMmJGnbf/+/cydO5fGjRsHlEqkcMT9qNF3332XXr168dvf/pbatWvz2muvUbVqVaZOnRr3814GYcyYMcyYMYOpU6fm9AB37dpFvXr1WLFiBbVq1Qo4YWzKzMzkggsu4Morr+Smm25i165dDBs2jCpVqvDCCy8EHU9y0ajRyIvrHiGEZjNZtWoV9evXZ+fOnYwYMYI33nhDRTAgixYtolu3bnlOg1aqVInWrVuzZMmSAJPFtsTERN5//30OHz5Mx44d6d+/P+3bt2f8+PFBRxMpcHF9jfBH1atX54477gg6hgDVqlU7aioyd2fDhg1Uq1YtoFTx4fTTT2f06NExN52byInkq0doZqlmtsHMPjezu46zTQ8zW2tmn5jZi5GNKfHij3/8I88++yyzZ8/G3Tl48CB///vfSUhI0IANESkQJ+wRmllx4Angf4B0YKmZTXP3tbm2qQ/cDVzs7jvNrGpBBZbYVr9+fSZMmMCNN96Iu5OZmcn555/P9OnTi9SoUREpOvJzarQZ8Lm7fwFgZi8BXYG1uba5EXjC3XcCuHvBTUchMS81NZXPPvuMTz/9lPLly+dZqkdEJNLyUwhrAF/lep4OND9im7MAzOwDoDgw3N1nHLENZtYf6A9E/XByCVaxYsVITk4OOoaIxIFIjRotAdQH2gC9gGfMrNKRG7n70+7e1N2bxtqyMiIiUjTlpxBuAWrmen5GuC23dGCau2e5+ybgU0KFUUREJKrlpxAuBeqbWR0zSwB6AtOO2GYqod4gZnYaoVOlX0Qwp4iISIE4YSF090PAAGAmsA542d0/MbP7zKxLeLOZwA4zWwvMAe5w9x0FFfqX+PTTTxkwYACXXXYZf/7zn1m/fn3QkUREJArk6xqhu7/l7me5ez13/0e47V53nxZ+7O4+0N1T3L2hu79UkKFP1rJly2jVqhWnnXYa99xzD7/5zW+49NJL+fDDD4OOJiIiAYuLmWX+9re/8eCDD9KvXz8A2rVrR+3atbn77ruZPXt2wOlERCRIcTHpdqlSpdixYwfly5fPaTtw4ABly5bl0KFDulFbRIoMTbodeXEx6XaVKlXYvHlznra0tDROO+00FUERkTgXF4Xwpptu4tZbb2XXrl0AZGRkcOutt8b0CukiQXJ3lixZwpQpU0hPTw86jsjPiotrhHfffTdbt26lTp06JCcns379enr06MHQoUODjiYR8vnnnzNp0iT27dtHp06daNmypXr7Adm6dStdu3YlIyODs88+m379+tGnTx9Gjhypn4lEpbjoEZYoUYInn3ySDRs28Oijj7J+/XrGjh1LyZIlg44mETB+/HhatmzJzp07KVmyJL1792bAgAEEdf073vXt25cOHTqwfv16pk2bxhdffMGCBQv473//G3Q0kWOKi8EyErt27NjBmWeeyYcffshZZ50FhFZbb9q0KU8++STt2rULOGF8+eabb0hJSWHr1q15FreePn06//rXv5g7d25w4WKEBstEXlz0CCV2zZw5kzZt2uQUQQittt6vXz+mTp0aYLL4lJmZSWJiIgkJCXnaq1Spwu7duwNKJfLzVAilSCtZsiQHDx7MeZ6RkcGoUaN46aWXWL58OV9//XWA6eJPvXr1KFmy5FH3544bN47U1NSAUon8PBVCKdJSU1NZvHgxH374IVu2bKFx48bMmTOHtLQ0qlevTuPGjVm+fHnQMeNGsWLFePzxx+nZsydDhgxhwoQJdO/enQULFjBo0KCg44kckwqhFGmJiYlMmDCBK6+8kksuuYTy5cuzYMECBg8ezKuvvsrDDz/MbbfdFnTMuJKamsqCBQs4ePAgb731Vs50hpUrVw46msgxabCMxISdO3dSt25dbrvtNvr06UPt2rUBOHToEBUrVmTr1q1UqFAh2JAiEaDBMpGnHqHEhFNOOYWqVavSsWPHnCII8MMPPwAcNXhDRORHKoQSM3r37s3QoUM5cOAAANnZ2QwdOpTf/va3lC5dOuB0IhKtVAglZtx5551UrFiROnXq0LNnT5KTk1m5ciX//ve/g44mIlFM1wgl5qxdu5ZVq1ZRt25dmjVrpmm9JKboGmHkxcVcoxJfUlJSSElJCTqGiBQROjUqIiJxTYVQRETimk6NisgJbdq0icmTJ7N//346depE06a6RCWxQz1CEflZzz//PBdeeCHp6ekcOHCAq666ioEDB2qZK4kZGjUqIsf14zJXixYtIjk5GQhNbN6kSROee+45Lr300oATxh+NGo089QhF5LhmzJjBZZddllMEASpWrMgf//hHXn/99QCTiUSOCqGIHFfx4sU5dOjQUe1ZWVkUL148gEQikadCKCLHdcUVV7BgwQJWrFiR07Zt2zaeffZZevToEWAykcjRqFEROa6KFSsybtw42rdvzxVXXEFiYiKvv/46N998M82bNw86nkhEqBCKyM/q1q0bF198MVOmTGHfvn0sXLiQM888M+hYIhGjQigiJ1SlShX69+8fdAyRAqFrhCIiEteNIGgLAAAPT0lEQVSKdCHcunUrq1evzll/TkRE5GQVyUK4a9curr76as4991x69uxJzZo1efbZZ9m2bRubN2+O2xkvZs6cSYsWLUhISOCss87iqaeeittjISKSX0WyEPbt25cqVaqQnp7O2rVrmTx5MjfffDP16tWjRYsWNGzYkPnz5wcds1DNnTuXPn36cNddd7F7926ef/55xowZw2OPPRZ0NBGRqJavKdbMLBUYBRQHnnX3fx5nu6uBV4EL3f1n50/7pVOsbdmyhfPOO4/09HTKlClDdnY2TZo0oVatWgBMnTqV6dOn069fP5YvX05SUtJJ76MoSk1N5dprr+Xaa6/NaVu/fj1t2rThq6++omTJkgGmE5FI0RRrkXfCHqGZFQeeAK4AUoBeZnbUqqdmlgjcCnwY6ZC5bdu2jRo1alCmTBkA3n//fdydQYMGsW3bNsyMLl26cM011zBu3LiCjBJV1q5dS6tWrfK0JScnk52dzY4dOwJKJSIS/fJzarQZ8Lm7f+HuB4GXgK7H2O5+4CFgfwTzHSU5OZmtW7eyYcMGANLT00lJSeH111/PUwjOPfdc0tPTCzJKVDnnnHNYuHBhnrZPP/0UM6Ny5coBpRIRiX75uY+wBvBVrufpQJ4pJcysCVDT3d80szuO90Fm1h/oD/ziU5ZlypThvvvuIzU1leHDh5OQkMC0adM45ZRTWLJkCQDuzrRp0+jSpcsv2kdRdOedd9K7d28qVqzI5ZdfzqpVq7jxxhsZNGiQTouKiPyMX31DvZkVA0YC159oW3d/GngaQtcIf+k+//SnP1GnTh3GjBnDN998Q926dSlbtiyffPIJW7ZsYezYsWzevJk//OEPv3QXRU67du145plnGDZsGJ07d6ZWrVoMHDiQAQMGBB1NRCSqnXCwjJm1BIa7++Xh53cDuPuD4ecVgY3AD+G3/Ab4HujycwNmIrke4eHDhxk7diwTJkxgz549dOzYkTvuuINTTz01Ip9f1Lg7ZhZ0DBEpABosE3n5KYQlgE+BdsAWYClwjbt/cpzt5wKDC2rUqIhIPFMhjLwTDpZx90PAAGAmsA542d0/MbP7zCx+LsKJiEhMytc1Qnd/C3jriLZ7j7Ntm18fS0REpHBo9YkY4+4sXryYdevWkZycTMuWLXW9UETkZxTJKdbk2DIzM2nXrh19+vRh/vz5XH/99bRt25bdu3cHHU1EJGqpEMaQIUOGkJSUxPr16xk/fjzr16+nTp063H333UFHExGJWvmaa7QgaNRo5J166qmsWrWKmjVr5rRt2bKFc889l127dgWYTEQiRaNGI089whiyb98+ypcvn6ctMTGRffv2aTkmEZHjUCGMIZ07d2bMmDF52p588kk6deqkATMiIsehUaMx5KGHHqJ169asWbOGVq1a8cEHH/D+++8zd+7coKOJiEQt9QhjSJ06dVi9ejUtW7bk448/pkWLFqxZs4Z69eoFHU1EJGqpRxhjKlWqxC233BJ0DBGRIkM9QhERiWsqhDEoMzOTbdu2aaSoiEg+qBDGkO+//55evXpRvXp1kpOTady4MfPmzQs6lohIVFMhjCHdu3encuXKbNmyhe+++47hw4fzu9/9jo0bNwYdTUQkaqkQxohVq1bxxRdfMGrUKCpUqECxYsXo1q0b/fr1Y+zYsUHHExGJWiqEMSItLY0GDRpQvHjxPO3nnXceaWlpwYQSESkCVAhjROPGjVm0aBGZmZl52mfMmEHTppqWUETkeFQIY0RSUhI9evSgY8eOzJs3j3Xr1vHXv/6V+fPn88c//vGkP2/btm0MHjyY8847j1atWvGf//yH7OzsAkguIhIsFcIY8vjjj9O9e3duu+02OnXqRGZmJgsWLODUU089qc/JyMigVatW7Nu3j3HjxjFkyBDGjh3LoEGDCii5iEhwtAyTHGXkyJEsXbqUSZMm5bTt2rWLunXrsnr1as4444wA04nENy3DFHnqEcpRlixZQufOnfO0VapUiYsuuojly5cHlEpEpGCoEMpRatSowdq1a/O0ZWdns379evUGRSTmqBDKUfr378/YsWOZNWsW7s6+ffv429/+xumnn06TJk2CjiciElFafUKOcvbZZ/PCCy8wYMAA9uzZw549e7j00kt5/fXXtcCviMQcFUI5pg4dOrB+/Xo2b95M+fLlqVKlStCRREQKhAqhHFexYsWoU6dO0DFERAqUrhGKiEhcUyEUEZG4pkIoIiJxTYVQRETimgqhiIjENRVCERGJa/kqhGaWamYbzOxzM7vrGK8PNLO1ZrbazN4zs1qRjyoiIhJ5JyyEZlYceAK4AkgBeplZyhGbfQQ0dffzgFeBhyMdVEREpCDkp0fYDPjc3b9w94PAS0DX3Bu4+xx33xt+uhjQzMwiIlIk5KcQ1gC+yvU8Pdx2PP2At39NKBERkcIS0SnWzOxaoCnQ+jiv9wf6AyQlJUVy1yIiIr9IfgrhFqBmrudnhNvyMLP2wBCgtbsfONYHufvTwNMQWqH+pNOKRCF3Z8aMGbzyyitkZ2dz1VVX0blzZ63UIVJE5OfU6FKgvpnVMbMEoCcwLfcGZtYYGAt0cfdtkY8pEr0GDhzIwIEDadKkCc2aNWPIkCHcdNNNuOt3PZGi4IQ9Qnc/ZGYDgJlAcWCcu39iZvcBy9x9GvAIUB54Jfxb8Jfu3qUAc4tEhdWrV/Pyyy+zdu1aKlasCECfPn1o0KABS5YsoXnz5gEnFJETydc1Qnd/C3jriLZ7cz1uH+FcIkXCrFmzuPrqq3OKIEC5cuXo2bMnM2bMUCEUKQI0s4zIr5CYmMh33313VPv27dupUKFCAIlE5GSpEIr8Ct27d2fmzJksWLAgp23p0qVMmTKFnj17BphMRPJLK9SL/AqVK1dm0qRJXHXVVSQnJ1O8eHHWrFnDc889R7Vq1YKOJyL5oEIo8it16NCBL7/8kjlz5pCdnc1ll11G2bJlg44lIvmkQigSAaVLl+aKK64IOoaI/AK6RigiMeOrr77i3XffJS0tLegoUoSoEIpIkZeVlUW/fv1o1KgRDzzwAE2bNuUPf/gD+/fvDzqaFAEqhCJS5D344INs2bKFzZs3M3v2bL766iv27t3Lvffee+I3S9yzoKaBatq0qS9btiyQfYtIbKlZsyYzZszg3HPPzWlLS0ujSZMm7NixI6bmfTWz5e7eNOgcsUQ9QhEp8r7//ntq1Mi7Olz16tXJyMggOzs7oFRSVKgQikiR165dOyZOnJin7YUXXqB169YUL148oFRSVOj2CREp8h544AHatWtHWloarVu3ZtGiRTz77LO8+eabQUeTIkA9QhEp8n5c7aNYsWI89dRT7N+/n0WLFnHhhRcGHU2KAA2WEREpQjRYJvLUIxQRkbimQigiInFNhVBEROKaCqGIiMQ1FUIREYlrKoQiIhLXVAhFRCSuqRCKiEhcUyEUEZG4prlGJaqtW7eODz/8kOrVq9OuXTtNoCwiEaceoUSlw4cP07dvX9q2bct7773HkCFDaNCgAZs3bw46mojEGPUIJSo9/fTTfPbZZ2zcuJGyZcsC8NBDD3HDDTcwe/bsgNOJSCxRj1Ci0gsvvMA999yTUwQBbr/9dlatWsWWLVsCTCYisUaFUKLSvn37SExMzNNWsmRJypQpw/79+wNKJSKxSIVQolLHjh0ZM2YMuZcJmz59OhUqVKBu3boBJhORWKNrhBKVBg0aRNu2benQoQNdu3Zl/fr1TJ48mddffx0zCzqeiMQQFUKJShUrVuSDDz5g8uTJLF68mBo1avDRRx9xxhlnBB1NRGKMVqgXESlCtEJ95OXrGqGZpZrZBjP73MzuOsbrpcxscvj1D82sdqSDioiIFIQTFkIzKw48AVwBpAC9zCzliM36ATvd/UzgMeChSAcVEREpCPnpETYDPnf3L9z9IPAS0PWIbboC/w0/fhVoZxrRICIiRUB+CmEN4Ktcz9PDbcfcxt0PARlA5SM/yMz6m9kyM1u2ffv2X5ZYREQkggr1PkJ3f9rdm7p70ypVqhTmrkVERI4pP4VwC1Az1/Mzwm3H3MbMSgAVgR2RCCgiIlKQ8lMIlwL1zayOmSUAPYFpR2wzDegTftwdmO1B3ZchIiJyEk54Q727HzKzAcBMoDgwzt0/MbP7gGXuPg34DzDBzD4HvidULEVERKJeYDfUm9l2ID+Ly50GfFfAcYoKHYuf6Fj8RMfiJ/FwLGq5uwZZRFBghTC/zGyZZlEI0bH4iY7FT3QsfqJjIb+EVp8QEZG4pkIoIiJxrSgUwqeDDhBFdCx+omPxEx2Ln+hYyEmL+muEIiIiBako9AhFREQKjAqhiIjEtagphFrzMCQfx2Ggma01s9Vm9p6Z1QoiZ2E40bHItd3VZuZmFrPD5vNzLMysR/jvxidm9mJhZyws+fg3kmRmc8zso/C/kyuDyClFiLsH/kVoxpqNQF0gAVgFpByxzZ+Bp8KPewKTg84d0HG4DCgbfvynWDwO+T0W4e0SgfnAYqBp0LkD/HtRH/gIOCX8vGrQuQM8Fk8Dfwo/TgHSgs6tr+j+ipYeodY8DDnhcXD3Oe6+N/x0MaFJ0GNRfv5OANxPaCHo/YUZrpDl51jcCDzh7jsB3H1bIWcsLPk5Fg5UCD+uCHxdiPmkCIqWQhixNQ+LuPwch9z6AW8XaKLgnPBYmFkToKa7v1mYwQKQn78XZwFnmdkHZrbYzFILLV3hys+xGA5ca2bpwFvAzYUTTYqqE066LdHJzK4FmgKtg84SBDMrBowErg84SrQoQej0aBtCZwnmm1lDd98VaKpg9ALGu/u/zKwloQUBGrh7dtDBJDpFS49Qax6G5Oc4YGbtgSFAF3c/UEjZCtuJjkUi0ACYa2ZpQAtgWowOmMnP34t0YJq7Z7n7JuBTQoUx1uTnWPQDXgZw90VAaUKTcYscU7QUQq15GHLC42BmjYGxhIpgrF4HghMcC3fPcPfT3L22u9cmdL20i7svCyZugcrPv4+phHqDmNlphE6VflGYIQtJfo7Fl0A7ADM7h1Ah3F6oKaVIiYpCGL7m9+Oah+uAlz285qGZdQlv9h+gcnjNw4HAcYfTF1X5PA6PAOWBV8xspZkd+Z9ATMjnsYgL+TwWM4EdZrYWmAPc4e6xdsYkv8diEHCjma0CJgHXx+AvzRJBmmJNRETiWlT0CEVERIKiQigiInFNhVBEROKaCqGIiMQ1FUIREYlrKoQiIhLXVAhFRCSu/X9dyTjFlh+bywAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Simulation window parameters\n", + "x_min, x_max = 0, 1; x_delta = x_max-x_min\n", + "y_min, y_max = 0, 1; y_delta = y_max-y_min\n", + "area = x_delta*y_delta\n", + "\n", + "# Sample from 2D Poisson process\n", + "lmbda = 30\n", + "nb_pts = np.random.poisson(lam=lmbda*area, size=1)[0]\n", + "poiss_pts = np.array([x_min, y_min]) + np.random.rand(nb_pts, 2)\n", + "\n", + "# Display\n", + "fig, ax = plt.subplots(1,1)\n", + "plt.title('Original process')\n", + "# Poisson point process\n", + "ax.scatter(poiss_pts[:, 0], poiss_pts[:, 1], \n", + " label='Poisson',\n", + " marker='o', edgecolor='k', facecolor='none')\n", + "plt.legend(bbox_to_anchor=(1.0, 1.0), frameon=False)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Build kernel and create DPPy object" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# Build kernel\n", + "kernel_name = 'RBF' # 'RBF','Cauchy'\n", + "sigma, alpha = 2.0, 1.0\n", + "L = lmbda*kernel(name=kernel_name, data=poiss_pts, **{'sigma':sigma, 'alpha':alpha})\n", + "\n", + "# Sampling from DPP(L)\n", + "DPP = FiniteDPP(\"marginal\", **{\"L\":L})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Sample from DPP(L)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[15, 5, 25, 29, 30, 20, 11, 4, 9, 18, 28, 1, 2, 0, 7, 26, 3, 23]" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Sample from DPP(L)\n", + "DPP.sample_exact()\n", + "dpp_sample = DPP.list_of_samples[-1]\n", + "dpp_sample" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Display the original process and its DPP-thinned version" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdoAAAEICAYAAAD4EjWLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3XmcT3X///HHaywz9goJYbhIiVREksZydYVEC0oLdWn7XV2l5VralHYtX3V1pX2TcimSpRBZS4hKCS2KshURlX1mXr8/zqGPMWM+mJkzn5nn/Xabm5n3OZ9znp/j8/m8zvt9zuccc3dEREQkfyRFHUBERKQoU6EVERHJRyq0IiIi+UiFVkREJB+p0IqIiOQjFVoREZF8VCwKrZndambP5/W8cSzLzax+XiyruNjXNjOzCWbWp6AzHQgze8DMro86RywzSzazL82s6n4+brqZXZ7DtNpm9ruZlciblPuVa7mZ/bmg15slQ8K8JiU6CVdozexSM1toZlvM7Ecze8rMDtnXY9z9fnfP9oPiYOaVguXundx9COx+HXxwIMsJH5sRFojfzWyZmb1kZkfFzJMaFv1d8yw3s5tjpruZbQ6nrTKzQbuKTVjIegPPhH+3NbPpux4Xs4zpZrYtXMbPZjbKzKrHTB9gZjtjMiwxs/Niprc1s8yY6b+b2biYxw6IXbe7bwdeBHY/j2y2zQAzezXebenuP7h7eXfPiPcxRUnsa1IkJwlVaM3sJuBB4J9AJeBkoA4w2cxK5/CYkgWXMH8UhedQCM129/IEr6M/A1uBj82scZb5Dgnn6wXcYWYdY6Y1Dad1AC4ErgjbLwXGu/vWOHL8PVxGfaA88EiW6a+Hhaw8cD3wqplVi5m+etf08OesXNY3DOhjZslxZCv29N6TvJAwhdbMKgJ3Ade6+0R33+nuy4GeQCpwcTjfADMbaWavmtmvwKVZ99LNrLeZfW9m682sf+wQVOy8Mb2aPmb2Q9jruC1mOS3MbLaZbTSzNWb2RE4FP5vnMz0cXvzIzH41szFmdliW9fY1sx+AqWF7VzNbFK5vupkdE7O8WmGPaF34vJ6ImfbXsDf0i5m9a2Z1wnYzs0fNbG2YYeGuQmNmnc1ssZn9FvbY/pHD8/iTmU0N1/mzmb1mMSMM4bb9h5l9bmabzOx1M0uJmf7PcNutNrO/xrHNLg+f99NAq7AXt3F/Msdy9wx3/9bd/wbMAAbkMN9sYBGQtRDj7l8C78dM6xQuK27uvhEYDRy/j3neBX4D/rQ/y86yjJXALwQ7qXsIdyJuBc4Pt+tnMZPrmNmscNtOMrMq4WN2vVZLhn9PN7N7cpk3p/dTkpndbGbfhq+nN3a9J8Lpl8S8b28jB2bW0oLRrhIxbeeY2ee5rSe7956ZpVjwebI+fO/Ns3Bnx2KG1cPl3h5mXGtmr5hZpXieuxRtCVNogVOAFGBUbKO7/w6MB06Pae4GjAQOAV6Lnd/MGgFPAhcB1Ql6NDVzWfepQEOCnssdMQUuA7gBqAK0Cqf/bT+eU2/gr2GOdODxLNPTgGOAMywY1vwfQa+mKsFzHmdmpcMPlLeB7wl2OmoCw8Pn243gw/Pc8HHvh8sB+AtwGnAUwXboCawPp70AXOXuFQgKyNQcnoMBDwA1wqy12LtY9QQ6AnWB4wh6fLs+2P9B8H/XgKBnmSt3XwJcTdgrdfddhT3ezDkZBbTJ2hjukLQGjgU+zWZ6o/Bxu6Y1Ab6KyTvd3duGv1t2KzazygT/R0tzmG5mdiZQGlic2xNx9wHhz+51x1gCNM3mMROB+/mjFx07z4XAZcDhYYZ97cTkNm9O76drgbMJXvc1CHYIBsPubfwUcEk4rTJwZA7PfS6wGWifJdOw3NYTY/d7D+hD8P6oFa73aoIRkKwuDX/aAfUIRiieyDJPTs9dirBEKrRVgJ/dPT2baWvC6bvMdvfR7p6ZzfBdd2Ccu3/g7juAO4DcLvh8l7tvdffPgM8IP6Tc/WN3n+Pu6WHv+hmCN2i8hrr7F+6+GegP9LQ9TyoZ4O6bw+dwPvCOu092950EQ4xlCHZAWhB8YPwznH+bu+86fnk18IC7Lwm33f3A8Rb0ancCFYCjAQvnWRM+bifQyMwquvsv7v5Jdk/A3ZeGmba7+zpgUDbb4HF3X+3uG4Bx/NFr6wm8FLMNBuzHtstOXJn3YTVwWJa2n4ENwPPAze4+JWbaJ2b2C8Fzeh54KWw/hKDnGY/HzWxTuJ4qBEUgVs+wx/47MBa4P+z97lIj7GXt+ukZxzp/CzPuj5fc/evwtfgG++h5xzFvtu8ngtfqbe6+MjyePADoHvaWuwNvu/vMcFp/IHMfGf5HMNyPmVUAOvPHDua+1rNL7HtvJ0GBrR+OgHzs7r9ms86LgEHu/l3YAbgFuCDLcnN67lKEJVKh/RmoYtkfM6keTt9lxT6WUyN2urtv4Y9eXE5+jPl9C8GeKmZ2lJm9HQ5T/UpQxKpkt4AcxOb8HiiV5fGx02uE8+zKnRlOr0mwp/19DjshdYD/7PogJigaBtR096kEe9yDgbVm9qwFQ/QA5xF8OH1vZjPMrFV2T8DMqpnZ8HCo9lfgVfbeBtluP7L8X8Q+vwMUV+Z9qEmwfWJVcfdD3f0Yd8864nBiOO1P7n57+H8CQQ+pQpzrvM7dKxH09A9l717aG+5+iLuXIxgy7m1mV8VMXx1O3/XzRhzrrABszHWuPeX0f3gg8+Y0vQ7wVsxrdQnBqFE19n7fbmbf79thwLkWHIs+F/jE3Xe9vva1nl1iX5dDgXeB4RYc4njIzEpls8493qPh7yWzLHd/tqMUEYlUaGcD2wneNLuZWXmCY2KxPY199VDXEPNhZmZlCPZWD8RTwJdAA3evSDBEm+3QYA5qxfxem2DPOXaHIfZ5rCb4gACCocTw8asIPhRq57ATsoJgODX2w7iMu38I4O6Pu3szoBHBEPI/w/Z57t6NYPhvNEHPJDv3hzmbhNvgYuLfBmvYexvEa6//4/3InJNzCIbWD9bnBNsybu6+ELgXGBz+32Y3z3JgApDbCU+5OYagN5Xtag5y2QdjBdApy2s1xd1XkeW1YmZl2cf71t0XExS6Tuw5bJzbenYvImZZO939LndvRDCC1IXgsE9We7xHCV7P6cBPcT17KbISptC6+yaCk6H+a2YdzayUmaUSfJiuJNjrjMdI4CwzO8WCE5cGsH/FMVYF4FfgdzM7Gvh/+/n4i82sUfihcTcw0nP+msQbwJlm1iHcm76JYMfjQ+Ajgg+igWZWLjx5o3X4uKeBW8zsWAAzq2RmPcLfTwpPHClFcExrG5AZHve9yMwqhcPUv5LzMF0FgmHNTWZWk7BQx+kNgpPVdm2DO/fjsT8BR4b/h+xn5t3MrISZ1TWz/wJtCV5jB2s8+3cIYZchBL2frtlNNLMjCY51LzrQYOH/0WHAnBxm+QlINbMoPhueBu6zP07WqxqeYwDB+7aLmZ0a/p/fTe6fX8OAfgTnIYyIcz17MbN2ZtYkPKzzK8EOcXavrf8BN4Svp/L8cbw7u5EmKUYSptACuPtDBL3GRwhe8HMJ9k47hMda4lnGIoLjYMMJitPvwFqCorW//kGwt/wb8Bzw+n4+fijwMsFwUgpwXU4zuvtXBL3F/xL0es8CznL3HWFxPovgKyI/EOx4nB8+7i2Cr0QND4d2vyDYyweoGOb+hWDvfz3wcDjtEmB5+JirCY4/Zecu4ERgE/AOWU5W2xd3nwA8RnDS0lL27+SlqQQF50cz2zUKEG9mCM9YJngdTSfYFieFPcuD9QrQORwtiVt4zsB/CI4/7rLrDODfgXnALA5uZ+BCYMg+3i+7CtJ6M9vfY9wH6z8Ex6EnmdlvBDsDLWH3+/YaguK5huA1uzKX5f2PYIdnqrvHjhTluJ4cHEFQ6H8lGGaeQfY79i+G7TOBZQQ7rlmPuUsxZF7Mb/we7nluJBj+XVaA650OvOrueXIVKilczOx+YK27PxZ1ll3C45WfAae5+9qo84gUF8Xyy9hmdhbBMV0j6B0vBJZHmUmKFne/NeoMWYW92KOjziFS3CTU0HEe6kZw4sJqgu9vXuDFvWsvIiL5otgPHYuIiOSn4tqjFRERKRCRHaOtUqWKp6amRrV6EZGE9PHHH//s7vt1q0OJVmSFNjU1lfnz50e1ehGRhGRmB3sFNSlgGjoWERHJRyq0IiIi+UiFVkREJB+p0IqIiOQjFVoREZF8lGuhNbMXzWytmX2Rw3Qzs8fNbKmZfW5mJ+Z9TBERkcQUT4/2ZYJbc+WkE8FlDBsAVxLco1VERESI43u07j4zvO9rTroBr4TXCp5jZoeYWXV3X5NHGUVE4paRkcHkyZP58ssvOfbYY+nQoQNJSTpKJtHJiwtW1CS4J+wuK8O2vQqtmV1J0Ouldu3aebBqEZE//Pzzz5xxxhkkJSXRqlUrhgwZQnJyMhMmTODQQw+NOl6RZmYZBHdCK0lw394+7r5lH/N/6O6nFFS+KBXobp67P+vuzd29edWquoKYiOStf/3rX5x66ql89NFHPP7443zyySc0bdqUW28tdHctjIyZlTWzG81sipmNNbMeZmZ5sOit7n68uzcGdgBX72vm4lJkIW8K7SqgVszfR4ZtIlLAFi9ezLXXXkvXrl25++67Wbu2+Nzf3d15/fXXuf3229lVN8yM22+/nddffz3idIWDmSUDk4E0YBAwHOgPPJTHq3ofqB+u80Yz+yL8uT4my+/hv9XNbKaZLQjnaWNmJczs5fDvhWZ2Qzjv8WY2Jzzx9i0zOzRsn25mD5rZR2b2tZm1yePnc1DyotCOBXqHZx+fDGzS8VmRgjd58mTatm3LuHE38d13L7By5UqaN2/OihUrcn9wEZGenk7p0qX3aEtOTiY9PT2iRIVOT2A7cLa7v+PuwwiK7mVmVicvVmBmJQlOkl1oZs2Ay4CWwMnAFWZ2QpaHXAi86+7HA02BBcDxQE13b+zuTYCXwnlfAf7t7scRDFPfGbOcku7eArg+S3vk4vl6z/+A2UBDM1tpZn3N7Goz2zUsMB74DlgKPAf8Ld/Siki23J1+/foxZMgQUlNTqVKlKs8++ywXX3wx9913X9TxCoSZ0a1bNx577LE92h999FHOPvvsSDJt2LCBgQMHcu6553LttdfyxRfZfkuyIKUBb3jMjcjd/RdgCnCwvcAyZrYAmA/8ALwAnAq85e6b3f13YFQ265lHUOgHAE3c/TeCmlLPzP5rZh2BX82sEnCIu88IHzcEOC1mOaPCfz8GUg/yueSpeM467pXLdAeuybNEIrLffvzxR7755jkefPAUZoQfQ23bwpYtt7FhQ9NIsxWkRx55hHbt2vHxxx/TqlUrZs6cybfffsv06dMLPMtPP/1E69atad26NV999QyffrqZ119vwUsvvcSZZ55Z4HlCa4G62bTXBX46yGVvDXulu8Vz6Df8ZstpwJnAy2Y2yN1fMbOmwBkEx3p7Ajfksqjt4b8ZRHhnuuzonHeRIqBcuXK4Z5KenrFH+/btO4rV2ba1a9fm888/55xzzmHjxo1ccMEFfPbZZ9SoUaPAszz00EN07tyZIUOGULVqVerUSWX48OH069ePzMzMAs8Tegn4q5mdAmBmSWb2N+BQYGo+rO994OzwBKxywDlh227hkPVP7v4c8DxwoplVAZLc/U3gduBEd98E/BJz/PUSYAYJoFBVfRE5MBUrVuT885+mQoXXKFHiScySeOutX+jUqROXXXZ51PEKVLly5bjsssuijsHkyZMpVWoWbduye5Th7rvbsWLFUJYtW8af/vSnAs/k7t+YWV9gpJn9DFQC1gNnunvGvh99QOv7xMxeBj4Km55390+zzNYW+KeZ7QR+B3oTfEX0JTPb1Rm8Jfy3D/C0mZUlGF6O/j86Diq0IkXE4MGD6dWrF3PmzKFs2XLUq9eWv/71r1xxxRVRRyuWDjnkEH75ZccebZmZmWRkZFCxYsWIUoG7jzWzCQQnHG0GlsQesz2I5ZbPoX0QwRnO2c7v7kMIjrdmtdflfN19AcFJVVnb28b8/jOJdoxWRBLDIYccwoQJE/jyyy/5/vvvadp0CUcccUTUsYqtvn378vjjf+HNNydx3nmVcXeaNfsnhx32HVWrjo40m7vvJDgJSQqACq1IEXP00Udz9NFHRx2jUNi6dSu//vorVatWLfDLMPbu3ZvFixdTv359zKazZcsWzD5hxIgRBZpDoqeToUSkyNm+fTv9+vWjevXqHHvssdSvX5833nijQDOYGQ8++CBLlizhtddWM29eeaZPn46uilf8qEcrIkXO9ddfz6pVq/jyyy+pVq0as2bN4oILLqBq1aq0a9euQLMcccQRdOrUqUDXKYWL5cEx8APSvHlznz9/fiTrFpGia+PGjaSmpvLdd99x2GGH7W5/6aWXGDNmDKNHR3t89GCZ2cfu3jzqHBI/DR2LSJGyZs0aqlWrtkeRBWjatCnLly+PJpQUayq0IlKk1K1blw0bNrB06dI92idOnEizZs0iSlX0mVlGeGOARWb2mZndtOt7sGbW1sw2hdOXmNmd+2rPZtkvm9mycL5PzKxVLlnuNrM/5/2zPDA6RisiRUpKSgq33norXbt25eGHH+boo49mzJgxPPbYY8ycOTPqeEXZ7kswmtnhwDCgIn9c4P99d+8SXiFqgZmNy6nd3T/JZvn/dPeRZvYX4BnguJyCuPsdefWk8oJ6tAkgIyODtWvXsmPHjtxnFhFuuOEGbr/9du677z46dOjARx99xNSpU/W1pxhmTDdjen4s293XAlcCf896r1t330xw4f/68bRnYyZ/3IIvp9vmvWxm3cPfB5rZ4nCeR8K2HuEt+D4zs5lhW4qZvRTelu9TM2sXtl9qZqPMbKKZfWNm+31LQRXaQu6FF14gNTWVmjW/oVy5j7jzzjvJyMjzK6WJFDkXXnghH374IcuXL2f48OE0btw46kjFirt/B5QADo9tN7PKBFd3WhRPezbOIrhFHuz7tnm7lnkOcGw4z73hpDuAM9y9KdA1bLsmiO1NgF7AEDNLCacdD5wPNAHON7PYe7DnSoW2EBs5ciQPPPAAY8aMoXXr1jRr1pxp06Zxzz33RB1NRBJUTE82DUjLz55tFm3M7FNgEjDQ3Rfl0p7Vw+Ft+K4E+sZx2zyATcA24AUzOxfYErbPIrhT0BUEOwMQ3NLvVQB3/xL4HjgqnDbF3Te5+zZgMbBf9+5VoS3EBg0aRJkyc7nxxhOZMQPmzk1hx453ue++09m5c2fU8UREcmRm9QhuWbc2bHrf3U9w92bu/nTMrHu1h0O4C8xsfMx8/3T34939dHeP68a+7p4OtABGAl2AiWH71QR3BaoFfBz2fPdle8zv+30bPp0MVYgtX76cOnX2vE53SkoZ3J3ffvttr68viIjkxp22EPRsY//OS2ZWFXgaeMLdPZ770sZy91zvyuPum8zsFzNr4+7vk81t88ysPFDW3ceb2SyCO/5gZn9y97nAXDPrRFBw3wcuAqaa2VFAbeArsrm5wf5SoS3ETjzxRM4551X69u1L27ZB28CBc+jV6xIOOeTbSLOJiGRRJhzaLQWkA0PJ5q49eSy32+ZVAMaEx1oNuDFsf9jMGoRtU4DPgC+Bp8xsYZj/Unffvr87CdnRlaEKsTlz5tC1a1fuu+8+XnjhEn7//Tc2b27BXXfdRe/evaOOJyIR0JWhEo+O0RZiJ598MuPGjePtt99m9eqjqFy5O4899piKrIhIAtHQcSHXsmVLxowZE3UMEYmxceNGRo8eze+//84ZZ5xBgwYNoo4khZh6tCIi++G9996jfv36jBs3js8++4zWrVvTv3//qGNJIaYerYhInLZu3cqFF17IqFGjOO204CubAwcOpEWLFnTo0IG2u85aFImhHq2ISJymTp3Kscceu7vIAlSuXJm//e1vDB8+PMJkUpip0IqIxGnnzp2ULl16r/bk5GTS09MjSCSJQIVWRCRO7du3Z968eSxcuHB325YtW3j22Wc5++yzI0wmhZmO0YqIxKlixYo89dRTtGvXjl69elG5cmWGDRvGqaeeyplnnhl1PCmk1KMVEdkP559/PvPnz6d69eps376dF154gRdeeIG8uIKQFE26MpSISDYyMzMZO3YsY8eOpVSpUvTq1atQnFWsK0MlHg0di4hk4e707t2bRYsWsXnz22RmZjJ5chq9e/dmwIABUceTBKOhYxGRLKZMmcInn3zC7NmzqVGjJkceWYu5c+fyxBNPsGzZsqjjSYJRoRURyWLixIls2zaBjh1TmDEDZsyAHj2qUrLkB0yaNCnqeJJg4iq0ZtbRzL4ys6VmdnM202ub2TQz+9TMPjezznkfVUSkYFSoUIGdO3fu1b5z504qVqwYQSJJZLkWWjMrAQwGOgGNgF5m1ijLbLcDb7j7CcAFwJN5HVREpKBcdNFFbNvWisGDF5GWBmlp0L//FJKS2nPWWWdFHU8STDwnQ7UAlrr7rjvTDwe6AYtj5nFg125eJWB1XoYUESlI9evX57HHHqNNmzbANDIzM7nooosYOXIk5cuXjzqeJJh4Cm1NYEXM3yuBllnmGQBMMrNrgXLAn7NbkJldCVwJULt27f3NKiJSYC666CLOOusspk2bRunSpWnf/nuSk5OjjiUJKK++3tMLeNnd/8/MWgFDzayxu2fGzuTuzwLPQvA92jxat4hIvqhYsSLdunWLOoYkuHhOhloF1Ir5+8iwLVZf4A0Ad58NpABV8iKgiIhIIoun0M4DGphZXTMrTXCy09gs8/wAdAAws2MICu26vAwqIiKSiHIdOnb3dDP7O/AuUAJ40d0XmdndwHx3HwvcBDxnZjcQnBh1qUd1bUcRyRPffvstI0aMYOfOnXTr1o3jjjsu6kgiCSmu79G6+3h3P8rd/+Tu94Vtd4RFFndf7O6t3b2pux/v7vpGt0gCe/bZZ2nZsiWrV69m06ZNdOzYUZceFDlAuqmAiOxh9erVNG7cmPnz51OvXj0A1q1bx/HHH8+4ceM48cQTI05YvOmmAolHl2AUkT2MGzeOLl267C6yAFWrVqVPnz6MGjUqwmQiiUmFVkT2kJSUREZGxl7tGRkZJCXpI0Nkf+ldIyJ76Nq1KxMmTGDJkiW721avXs2QIUPo0aNHhMlEEpPuRysie6hWrRr//e9/ad26NV27diU5OZlRo0bx73//myZNmkQdTyThqNCKyF4uuugi2rdvz6hRo9i5cydz587d45itiMRPhVZEslW9enWuueaaqGOIJDwdoxUREclHKrQiIiL5SIVWREQkH6nQioiI5COdDCUihUp6ejrvvPMOn3/+OQ0aNOCcc87RDdcloanQikihsWHDBk4//XSSk5NZtepVfv31N/r3b8zUqVOpVatW7gsQKYQ0dCwihcYdd9zBSSedxKxZs6hbtx5Nmzbl4osv5vrrr486msgBU49WpJD5+eefee6551iwYAH16tXjqquuIjU1NepYBeLNN9+kTp3vaNfOmDEjaMvIuI1Zsz5k586dlCpVKtqAIgdAPVqRQuT777/nxBNP5JtvvmHhwsd59dXLadGiBXPmzIk6mogcIPVoRQqRAQMGcOmll3L33XfTti0cfng1rrrqP9xwww3Mnj076nh5bvHixbz99tskJyfTvXt3unfvzs6dN/DUU0/Rrp0B0LbtvVSr9gWlSp0WcVqRA6NCK5H68ccfmThxIqVLl6ZLly5UrFgx6kiRmjRpErVqLWXmTHYPnbpfwEcf1WTTpk1UqlQp2oB56M477+SZZ56hZ8+ebN68mbvuuouHHnqIp59+mlNOOYXVq1/j119/Y/XqYUyZMiXquCIHTIVWIjN48GBuv/12SpWaRWZmBtdeW5ehQ4fSuXPnqKNFpmLFiuzcuZOUlDK729LT0wErUl9xmTt3Li+99BILFy6katWqANxwww20adOGr776itmzZ7Nw4TAaNGjA2WcvLFLPXYofFVopMF9++SXjx4+nbNmyNGrUiLvvvptPP/2USy9NBeDBB8fTuXNnvvvuuyLVc9sfl112GZMmncuYMWM488xyuDsNG17DUUftJCWlTdTx8szIkSPp27fv7iIL0LhxY9q3b8+7777LJZdcQrdu3SJMKJJ3dDKUFIg77riDtLQ0li5dygcffMAZZ5xBZuYULr00lRkzgmHSf/+7JWYzePvtt6OOmycWLVpE3759admyJZdccgmffPJJro+58cYbqV27NqmpqSxcuJC5c+ewfPlyHnvssQJIXLDMLNs2d48gjUj+UY9W8t2sWbN45ZVXWLRoEVWqVAGCD9RhwzaQkZFO7MuwRIkktm3bFlHSvDNv3jw6d+7MjTfeyOWXX86cOXPo2LEjw4cPp3379jk+rmTJkrz44ossW7aMzz//nLp163LccZMKMHnBOO+88zj//PO55pprqFy5MhCcGDVlyhSeeeaZiNOJ5C0VWsl3I0aM4Morr9xdZAGuuOIKRo3qxOWXP0mJEn0AGDp0BU2bnkqnTl9EFTXP9O/fn4EDB9K3b18AWrVqRZ06dbjllluYO3duro+vW7cudevWze+YkTn55JPp3bs3TZo04fzzz2fz5s28+eabPPHEE7sLr0hRYVEN0zRv3tznz58fybqz2rBhA6+//jrr1q0jLS2N0047LdthLTkw1113HTVr1uTf//737jZ3p27dumzbto0yZeaSmZnBjh2tufnmm+nXr1+EafNGhQoVWLFiBYcccsjutoyMDJKTk9m6dWuhvfBCZmYmS5cupWzZshx55JH5vr4vvviCcePGkZKSQo8ePQpknYnOzD529+ZR55D4FftjtLNmzaJhw4a8//77bN26lauuuooePXqEZ3pKXjj33HN5/vnn2bRp0+62xYsXs2nTJp5//nnOOec/XHjhc7z77rtFosgCHHHEEXzzzTd7tC1btoxDDz2UkiUL50DS1KlTOfrooznuuPXUrfs97du354cffsjXdTZu3JhbbrmFG264QUVWii53j+SnWbNmHrWMjAyvV6+ejx07dnfb9u3b/bTTTvPnn38+wmRFS2Zmpl9//fVeq1Yt/9e//uVXX321H3bYYT5s2LCoo+WbQYMFiVULAAAYuUlEQVQGeatWrXzNmjXu7v7zzz/76aef7v379484WfaWLVvmVapU8QkTJnhaWqafdlqG33fffd64cWPPyMiIOp7EAOZ7RJ/b+jmwn2Ldo12wYMHuCyXsUrp0afr168fIkSMjTFa0mBmPPvoob731FuXKlaNevXp8+umn9OrVK+po+aZfv36kpaVxzDHHcNxxx1G/fn2OPvpo7rjjjqijZevFF18kJWU2Awd2ZMYMY+bMJCZNuoVly15i5syZUccTSWiFcwyrgJgZmZmZe7VnZmbqGG0+aNasGc2aNct1vjVr1nDfffcxfvx4KlSoQJ8+fbjuuusK7ZBrdpKSknjggQe4+eabWbZsGbVr1+awww6LOlaOVq1aRZkyZbO0GmXLlmXVqq8iySRSVBTrHm3Tpk3JzMxkzJgxu9u2b9/OY489Ro8ePSJMVnxt2rSJNm3akJyczDvvvMOTTz7J22+/zZVXXhl1tANSqVIljj/++EJdZCE4K7p+/cuZNs1JS4O0NBg37jfc02jZsmXU8UQSWrE/63jOnDl069aNNm3akJqayujRoznppJN49dVXKVGiRNTxip1HH32UuXPnMnz48N1tW7ZsITU1lQ8//JD69etHmK7o2rJlC61bt6ZJkyYsXPhf0tPTKV++C8cffzxPPfVU1PEkhs46TjzFukcLwff5vv76azp16kTlypV55ZVXGDZsmIpsRD755BM6duy4R1vZsmVp06YNn376aUSpir6yZcsybdo06tatS4kSHTjssHPp27cvgwcPjjqaSMIr9j1aKVz69+/P5s2bGTRo0O62zMxMGjZsyGuvvUaLFi0iTCcSPfVoE09cPVoz62hmX5nZUjO7OYd5eprZYjNbZGbD8jamFBeXX345Q4cOZcSIEWRmZvLbb79x0003Ua1aNU466aSo44mI7LdcC62ZlQAGA52ARkAvM2uUZZ4GwC1Aa3c/Frg+H7JKMVCnTh3GjBnDgw8+SOXKlalRowarV69m9OjROhNcRBJSPN+XaAEsdffvAMxsONANWBwzzxXAYHf/BcDd1+Z1UCk+TjnlFObNm8e6desoU6YMFSpUiDqSiMgBi2fouCawIubvlWFbrKOAo8xslpnNMbOOZMPMrjSz+WY2f926dQeWWIoFM+Pwww9XkRWRhJdXZx2XBBoAbYFewHNmdkjWmdz9WXdv7u7NY2/4LCIiUlTFU2hXAbVi/j4ybIu1Ehjr7jvdfRnwNUHhFRERKdbiKbTzgAZmVtfMSgMXAGOzzDOaoDeLmVUhGEr+Lg9zioiIJKRcC627pwN/B94FlgBvuPsiM7vbzLqGs70LrDezxcA04J/uvj6/Qh+IefPm0bNnTxo3bkz37t2ZM2dO1JFERKQYiOsYrbuPd/ej3P1P7n5f2HaHu48Nf3d3v9HdG7l7E3cfvu8lFqzp06dz5pln0rZtW8qUmcPnnz9Ot27deO+996KOJiIiRVyxuARj//79+e9//8vf/vY3ypUrT40aNXjyySfp379/1NFERKSIS5z7jh2Ejz76iMGDu/PUUzBjRtDmfg5z5lTG3XUhBBERyTfFokdbs2ZNNm/eskfb5s1bKF06WUVWRETyVbEotMFNw//MK6/8QFoanHzydkqX/gv33vtB1NFEiqSdO3cyduxYnn76aRYsWBB1HJFIFYuh4+uuu46NGzdywgknsGXLO6Snp3PzzR246aaboo4meWT27Nm8+uqrbN26lS5dutCtWzfd6jAiS5cupWPHjqxbN4KyZRtSqtTptGnThiFDhlCyZLH4yBHZQ7Ho0SYlJTFgwABWrFjB559XZuPGE7jnnntISioWT7/Ie+SRR+jZsye1a9emefPm3H///VxwwQVkZmZGHa1Yuuyyy7j22ms54YQTaNiwIV9//TWrV6/mmWeeiTqaSCR0P1pJaKtXr+bYY49l0aJF1KhRA4AdO3bQokUL7r33Xrp06RJxwuJlxYoV1Kv3A6eccgozZwbnP6SlwS+/bKB8+bOYNWtWxAkTn+5Hm3jUpZOENmXKFE4//fTdRRagdOnS9OnThwkTJkSYrHjasWMHSUlG1nMMk5KS2L59ezShRCKmAyaS0MqXL8+GDRt2/71y5UqeeOIJRowYQaVKlVi8eDGNGjXaxxIkL9WrV48GDbpx3XV3A+cCMG2ac8EFV9G4cbdow4lERD1aSWgdO3bkiy++YMKECXzzzTecdNJJ/PTTT3z//RB++ul12rZty7Rp06KOWWyYGc899xxXX301X331FStWrCAtLY0ffviBG264Iep4IpHQMVpJeB9++CHdu3cnPT2dww47jJ9++onDD19M9erV6dfvLe655x4+/vhjfWe6AP30008MHTqUVatWcfLJJ3PuuedSqlSpqGMVCTpGm3hUaKVI2L59O9WrV6dKlS84/PCqzJoVfKinpTnvv/8BGzYcR6VKlSJOKXLwVGgTj47RSpGQnJxMjRo1KFWqHCVL/tFz2rFjJ2ZQpkyZCNOJSHGmY7RSZFxxxRWUK3cmb731C2lp0KZNJnXrXsYVVwyjdOnSUccTkWJKPVopMq699lqWL19OvXr1cJ/Kli1bqFZtB4888kjU0USkGNMxWilyfvrpJ5YsWUJqaiqpqalRxxHJUzpGm3jUo5Uip1q1alSrVi3qGCIigI7RioiI5CsVWhERkXykoWMRydXs2bMZNmwY27Zto0uXLpx11lm6+5VInFRoRWSfHnnkEf7zn/9QqtQHlChRkvnzuzB8+HBee+01FVuROOhdIiI5Wr16Nffffz9z586ldu061KxZk9mzZ7No0SImTZoUdTyRhKBCKyI5mjx5MiVLfsCFF9ZgxgyYMQM6dkxh8+Z3eOedd6KOJ5IQVGhFJEflypUjPX3nXu3p6TspX758BIlEEo+O0YpIjjp37kypUnW57bZXgdMBeP75pbRu3ZqLL54abTiRBKFCKyI5Klu2LCNHjqRHjx5s2zaREiVKcNJJp/Hwww9z7LHHRh1PJCGo0IrIPrVp04bly5fz3nvvsW3bNtq3/5bDDjss6lgiCUOFVkRylZKSQpcuXaKOIZKQdDKUiIhIPkrYQuvuzJ8/n7fffpu1a9dGHafQyMjI4Mcff2T79u1RRxERERK00K5cuZIWLVpw4YUX8sQTT9CwYUPuuOMO5s+fz/vvv8+2bduijhiJF198kdTUVGrX/pby5edx6623kp6eHnUsEZFiLSEL7SWXXEKXLl346quvmDhxImPGjOGhhx6idesdnHlmWWrXrs2IESOijlmgRo8ezb333svo0aM55ZTWNG/enDlz5nDnnXdGHU1EpFiLq9CaWUcz+8rMlprZzfuY7zwzczPLt5sSL1++nMWLF3PrrbdiZqSnp9OnTx/69OlD+fLlOfHEZkycOJFrrrmGL7/8Mr9iFDqDBg2ibNmPuOmmZsyYAXPmpLB163geeqiThpFFRCKUa6E1sxLAYKAT0AjoZWaNspmvAtAPmJvXIWNt2rSJypUrU6pUKQDee+891q9/k7lzB7Jhw3HMmAE33ngiKSlzePnll/MzSqGyfPnyva7Uk5ycAhgbN26MJpSIiMTVo20BLHX379x9BzAc6JbNfPcADwL5eoC0UaNGbNq0iY8++giADRs2kJKSzI8//rjHfCkpyaxfvz4/oxQqzZo149JLX2b6dEhLC37+7/8+pnr1XlSpUiXqeCIixVY836OtCayI+Xsl0DJ2BjM7Eajl7u+Y2T9zWpCZXQlcCVC7du39TwuUKlWKxx9/nLPOOou///3vHHrooWzceAL16tWjfv2FlCxZiqlTM0lLu4B27f7fAa0jEd1222106tQJd2fbtkv5/fff6dGjB3fddRclSpSIOp6ISLF10CdDmVkSMAi4Kbd53f1Zd2/u7s2rVq16wOs877zzmDx5MmvXrmXGjBm0a9eOzMxM1q37mXXr1tG5c2dKlixJ9+7dD3gdiaZ58+ZMmDCByZMns2bN0VSu3J1BgwZx2WWXRR1NRKRYM3ff9wxmrYAB7n5G+PctAO7+QPh3JeBb4PfwIUcAG4Cu7j4/p+U2b97c58/PcfJ+cXfeeecdhg4dyubNm+ncuTN//etfSUlJyZPli4gUFmb2sbvn2wmnkvfiKbQlga+BDsAqYB5wobsvymH+6cA/9lVkIW8LrYhIcaFCm3hyHTp293Tg78C7wBLgDXdfZGZ3m1nX/A4oIiKSyOK6qYC7jwfGZ2m7I4d52x58LBERkaJBd+8pYrZt28aoUaNYsmQJxxxzDOeee66OVYuIRCghL8Eo2VuzZg3HH388L730EiVLluTll1+madOmrF69OupoIiLFVq4nQ+UXnQyV93r37k316tV58MEHd7fdcsstrFy5kqFDh0aYTETyik6GSjwqtEVI+fLlWbZsGbHfUV6/fj21atVi8+bNmFmE6UQkL6jQJh4NHRchSUlJZGRk7NGWkZFBUpL+m0VEoqJP4CKkZ8+e3H///ewapXB37r//fnr27KnerIhIRHTWcREycOBATj/9dE455RROPfVUPvjgA7Zu3crkyZOjjiYiUmzpGG0Rk56ezoQJE1i8eDHHHHPM7us+i0jRoGO0iUeFVkQkgajQJh4doxUREclHKrRFjLuzaNEiZs+ezbZt26KOIyJS7KnQFiHLli2jZcuWNGv2K3/5Sylq167NsGHDoo4lIlKs6SyZIsLd6datG5deeilly54MGP/5z2Q6duzIMcccwwknnBB1RBGRYkk92iJizpw5fPvt84wdewMzZhgzZkC/fk1JSZnDCy+8EHU8EZFiS4W2iFi3bh3JySnAnhemSElJ4eeff44mlIiIqNAWFa1atcI9jf/9bw1paZCWBtOmOX/6U1/S0tKijiciUmyp0BYRVatW5R//+AdpaWn8+OMa1q9fT48ePfjxxx/p06fPfi9vyZIlnH322ZQtW5YjjjiCm2++WWcxi4gcABXaIuS2227j0UcfpUmT66hevRetWrVi2rRplC1bdr+Ws2bNGtq3b0/btm1ZvXo177//Pl999RW9e/fOp+QiIkWXrgwle7nrrrtYu3YtgwcP3t22fft26tSpw8yZMznqqKMiTCdSvOnKUIlHPVrZy+LFizn11FP3aEtOTqZFixYsXrw4olQiIolJhVb20rBhQz788MM92nbs2MG8efM4+uijI0olIpKYVGhlL1dddRUjRozgySefZMuWLXz//fdccsklnHzyySq0IiL7SYVW9lKzZk0mT57M22+/TcWKFTnhhBOoWbMmr732WtTRREQSji7BKNlq0qQJ48ePx90xs9wfICIi2VKPVvZJRVZE5OCo0IqIiOQjFVoREZF8pEIrIiKSj1RoRURE8pEKrYiISD5SoRUREclHcRVaM+toZl+Z2VIzuzmb6Tea2WIz+9zMpphZnbyPKiIiknhyLbRmVgIYDHQCGgG9zKxRltk+BZq7+3HASOChvA4qIiKSiOLp0bYAlrr7d+6+AxgOdIudwd2nufuW8M85wJF5G1NERCQxxVNoawIrYv5eGbblpC8wIbsJZnalmc03s/nr1q2LP6WIiEiCytOToczsYqA58HB20939WXdv7u7Nq1atmperFhERKZTiKbSrgFoxfx8Ztu3BzP4M3AZ0dffteRNPpPDbsWMHgwYNolWrVrRs2ZIHH3yQrVu3Rh1LRAqJeO7eMw9oYGZ1CQrsBcCFsTOY2QnAM0BHd1+b5ylFCil3p3v37uzYsYMtW8ZjBrNm9WHy5MlMmjSJpCR9g06kuMv1U8Dd04G/A+8CS4A33H2Rmd1tZl3D2R4GygMjzGyBmY3Nt8Qihcj777/P119/zbhx4zj00EM55JBDeeutt1i/fj2TJk2KOp6IFAJx3Y/W3ccD47O03RHz+5/zOJdIQpg9ezabN7/D6aeXYsaMoK1DhxL88ssoZs9+mY4dO0YbUEQip3EtkYNQs2ZNtmzZslf7li1bqFlzXyfni0hxoUIrchDOPfdcypTpxMUXP89ppzlpac7ll79KUlIHzj///KjjiUghENfQsYhkr2zZsrz77rtcdtllLFgQXDBt06b/Y8KECVSqVCnidCJSGKjQihykY489lo8++ojly5fj7qSmfoKZRR1LRAoJFVqRPJKamhp1BBEphFRoRaRIyMzMZMqUKSxatIijjjqKM844gxIlSkQdS0SFVkQS38aNG+nUqRPbtm1j3bo32LTpVxo0OIlJkyZRpUqVqONJMaezjkUk4d122200adKETz75hPr1G9Cs2Ymcdtpp/OMf/4g6mgjm7pGsuHnz5j5//vxI1i0iRUuVKlVo0GAlyckpuy8ccsopO5kzZzbp6W2K1MlpZvaxuzePOofETz1aEUl46enpmO35cZaUZLg7UXUmRHbRMVoRSXhnn302RxxxBwMHDqRt26CtQ4cHOOKIz0lKSos0m4gKrYgkvAceeIC2bduyYMECfvjhGX799Vd++OElpk2bFnU0ERVaEUl81atXZ8GCBYwYMYJFi56kYcOGnH/+QsqVKxd1NBEVWhEpGsqUKUPv3r2jjiGyF50MJSIiko9UaEVERPKRCq2IiEg+UqEVERHJRyq0IiIi+UiFVkREJB+p0IqIiOQjFVoREZF8pEIrIiKSj3RlKCm0duzYwZtvvsncuXOpUaMGvXv35ogjjog6lojIflGhlULpt99+489//jMpKSmsXPkqW7Zs4ZFHjmPMmDG0atUq6ngiInHT0LEUSoMGDaJevXpMnz6dWrVq0bBhQ55++mmuvPJK3V9URBKKerRSKI0bN46MjCm0a2fMmLGr9Ry+/voIli1bRr169aKMJyISN/VopVAqXbo0mZkZe7S5O5mZTnJyckSpRET2n3q0UihdeOGFvPnmeUyYMIGOHVMA6NbtP8Boatacse8Hi4gUIiq0UihdffXVzJ49m/r165ORMYUtWzazatWTTJw4MepoIiL7RYVWCqWSJUvy2muv8dlnnzFnznRq1qxJx45LKFlSL1kRSSz61JJCrWnTpjRt2jTqGCIiByyuk6HMrKOZfWVmS83s5mymJ5vZ6+H0uWaWmtdBRUREElGuhdbMSgCDgU5AI6CXmTXKMltf4Bd3rw88CjyY10FFREQSUTw92hbAUnf/zt13AMOBblnm6QYMCX8fCXQwM8u7mCIiIokpnkJbE1gR8/fKsC3bedw9HdgEVM66IDO70szmm9n8devWHVhiERGRBFKgF6xw92fdvbm7N69atWpBrlpERCQS8RTaVUCtmL+PDNuyncfMSgKVgPV5EVBERCSRxVNo5wENzKyumZUGLgDGZplnLNAn/L07MNV15XcREZHcv0fr7ulm9nfgXaAE8KK7LzKzu4H57j4WeAEYamZLgQ0ExVhERKTYs6g6nma2Dvg+jlmrAD/nc5xEoW3xB22LP2hb/KE4bIs67q6TXBJIZIU2XmY2392bR52jMNC2+IO2xR+0Lf6gbSGFkW6TJyIiko9UaEVERPJRIhTaZ6MOUIhoW/xB2+IP2hZ/0LaQQqfQH6MVERFJZInQoxUREUlYKrQiIiL5qNAUWt3zNhDHdrjRzBab2edmNsXM6kSRsyDkti1i5jvPzNzMiuzXOuLZFmbWM3xtLDKzYQWdsaDE8R6pbWbTzOzT8H3SOYqcIru5e+Q/BFec+haoB5QGPgMaZZnnb8DT4e8XAK9HnTui7dAOKBv+/v+K4naId1uE81UAZgJzgOZR547wddEA+BQ4NPz78KhzR7gtngX+X/h7I2B51Ln1U7x/CkuPVve8DeS6Hdx9mrtvCf+cQ3CTh6IontcEwD3Ag8C2ggxXwOLZFlcAg939FwB3X1vAGQtKPNvCgYrh75WA1QWYT2QvhaXQ5tk9bxNcPNshVl9gQr4mik6u28LMTgRqufs7BRksAvG8Lo4CjjKzWWY2x8w6Fli6ghXPthgAXGxmK4HxwLUFE00ke7neVEAKJzO7GGgOpEWdJQpmlgQMAi6NOEphUZJg+LgtwSjHTDNr4u4bI00VjV7Ay+7+f2bWiuCGJ43dPTPqYFI8FZYere55G4hnO2BmfwZuA7q6+/YCylbQctsWFYDGwHQzWw6cDIwtoidExfO6WAmMdfed7r4M+Jqg8BY18WyLvsAbAO4+G0ghuNmASCQKS6HVPW8DuW4HMzsBeIagyBbV43CQy7Zw903uXsXdU909leB4dVd3nx9N3HwVz/tjNEFvFjOrQjCU/F1Bhiwg8WyLH4AOAGZ2DEGhXVegKUViFIpCGx5z3XXP2yXAGx7e89bMuoazvQBUDu95eyOQ49c9ElWc2+FhoDwwwswWmFnWD5kiIc5tUSzEuS3eBdab2WJgGvBPdy9qIz7xboubgCvM7DPgf8ClRXCnXBKILsEoIiKSjwpFj1ZERKSoUqEVERHJRyq0IiIi+UiFVkREJB+p0IqIiOQjFVoREZF8pEIrIiKSj/4/7NqVdyA99N0AAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Display\n", + "fig, ax = plt.subplots(1,1)\n", + "plt.title('Original process and its DPP(\"{}\") thinned version'.format(kernel_name))\n", + "# Poisson point process\n", + "ax.scatter(poiss_pts[:, 0], poiss_pts[:, 1], \n", + " label='Poisson',\n", + " marker='o', edgecolor='k', facecolor='none')\n", + "# Thinned (using DPP) Poisson point process\n", + "ax.scatter(poiss_pts[dpp_sample, 0], poiss_pts[dpp_sample, 1],\n", + " label='DPP-Poisson',\n", + " marker='+', facecolor='blue')\n", + "plt.legend(bbox_to_anchor=(1.0, 1.0), frameon=False)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Display underlying DPP kernel using DPPy" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAASwAAAETCAYAAACMUTsNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJztnXm8XdP5/99PJhnJ4IaMxBQhKiRVhBA0JabW0Ea1voZvUVpJ0aJfKqFV+tXW8FV+KhpaMYbS1pCoFCENmZABJYZMkpuQSebk+f2x1ol9T86z77n3nNx7zvW8X6/zys3+7LX22uvus+7eaz/r84iq4jiOUw40qu8GOI7j5IsPWI7jlA0+YDmOUzb4gOU4TtngA5bjOGWDD1iO45QNPmA5jlM2+IDlOE7Z4ANWHSMivxaRYfXdDgARmSkiRxaprg9F5Jj482sism8++9Y3pdQWp3qa1HcDGgIisiNQCXRV1fkp+1UAZwF71FXb0lBVc1ApkJuB64BTt1H9zpcUv8MqDn2AyrTBKnI28LSqrtnWDRKR+vxj9BQwUER23pYHqedzdOoBH7CKQx/gjTz2Ow54MbkhPpL8VETeFJHPRWSkiOwkIs+IyEoReV5E2sV9rxSR9+P2WSLyrRx1XSEibwKfi0gTETlQRKbFMo+KyMMi8svE/sdklb88tmV53Ld5PsdOoqprgSnAN6rrEBHpJSIfiMgZ8f+dRWSMiFTG7ZfkcY5p7U6tzykvfMAqDn2A6Xnstx/wTo7tpwJfB/YCTgSeAX4OVBB+R5kv2fvA4cAOwAjgLyLSKauuM4Djgbax7BPAKKA98CBgDjSRbwPHAj2ArxDuCvM9dpLZwP5pBxKRA4HngB+r6oMi0gj4G2Hw7wIcDQwTkeyBb8s5qupGq901qM8pE3zAKg75DlhtgZU5tt+uqoviI+XLwCRVnRbvVJ4ADgBQ1UdVdYGqblbVh4H/AAdl1XWbqs6Nj50HE+Ypb1PVDar6OPBaNW28LR7jU8KXvU8Njp1kZTxfi8MJj45nqerf47avAhWqep2qrlfVOcAfgSEp55jW7nzrc8oEnwMokPjo0ROYlsfunwFtcmxflPh5TY7/t47HOgu4FNg1aq2BHbPqmpv4uTMwX6t6CM0lnU8SP6+OdeR77CRtgGUp+oXAi6r6r8S2XYDOIpIs15gwiCfJdQ652p1vfU6Z4HdYhdMb2EDuR71s3iQ89tUYEdmFcHfwI6CDqrYFZgCStWtycFoIdBGR5D7dtuGxk/QifV7vQqC7iPw+sW0u8IGqtk182qjq4Kyy+Zq45VufUyb4gFU4fYBZQFMRaR4/2xn7Pg0cUcvjtCJ8USsBROQcwmCZxkRgE/CjODl9MumPcUU5drzr7AuMS6lzJWHOaYCI3Bi3vQasjJPqLUSksYj0FpGv1qLN26I+p57xAatw+hC+nGsSn1nGvvcDg0WkRU0PoqqzgN8SBqFFhAn8V6opsx44BTiP8Hj2PeDvwLptfOwTgX+p6oJq6l1GeNlwnIhcr6qbgBMIffoBsAS4hzDRX2OKXZ9T/4hbJNctInIDsFhVb6mn408C7lLVP23jY5ynqjO21TGcLyc+YDVwROQIwvzaEuBM4C5gN1VdWK8Nc5xa4G8JGz49gUcI81BzgNN8sHLKFb/DchynbPBJd8dxygYfsBzHKRt8wHIcp2zwActxnLLBByzHccoGH7AcxykbfMByHKd8UNV6+RAWvr4DvAdcmaXdCywGZmRt7waMJ6zVmwkMTWjNCYtd34jaiBzHbEywgfl71vYPgbcInlaTs7S2wGPA2wRTukPi9p5x/8xnBTAsUe4nsR0zCMZ5zRPa0Lh9JjAs5XyPBZYDG4FPcvTRmqjNyNFHn0VtYY4++jRqi3P00Z8I7hMrcvTRp1Fbk6OPPoh1rs3RR0tjuU05+ijTzmVZfXQtsCrWtzDrd/09goXMesLaxqT2cNy+NnmNxH55Ffg8u87YL9MTx1ucVWemT9cQFm0ntbmJcmuytH1jXevi5/eJfpkZ61ob++XxRLkRiXYuBy5LaJfFc8+cw4iEZn6nGsqnvgarxgQHy92AZoRBZp+EPgA4kK2/wJ2AA+PPbYB3M+UIViet489NgUnAwVnlLwVGk3vA2tFo633Af8efmxFcLnOdzyfALvH/XQhf4Bbx/48AZ8efexMGq5aElQbPExw0q5xvoo++Q3BYWJOjj86LF+eM7D6K+mHxi1Klj6J2ULzws/vo/wiOp7kGrBON38t9wE1Rm5ndR/F4/QiDVnYfHRPLLc/qo7eBQ2MfjY/77hP7JdOWZoQ/NB8mzvGbwHdjH2+5RmK/DIrHakMwIPwwagLsHrWmwOvAx4k6OxEWf4+OfZO87uYCRxnX5CPAdfHn9vF3tU/W76kxYdCdE9vSJR47M+g/Tri29uGLa6ci9ss/CQPtwVTznWoon/p6JDwIeE9V52hwFHgIODkjqupLhL/mVVDVhao6Nf68knDH0yX+X1V1Vdy1afxsCeMXka4EW9178m2kiOxA+LKNjMdYr8FhIJujgfdV9aPEtiZAi5gooSWQcS7oRXAUXa3B3vdFwl/w7PPN9NHDhL/Sy9m6j/5J+OtMYvtCVZ0a9XmEAatKH0VtBeGLmt1HvQiuErmYmN3ORB9dGTXN7qN4vF7A+hx9NDWem2T10Suq+mrsoxcIg2uX2C/vqOrf4rUzOqGhqn8l3ElVuUZiv4yNfbOScJe+IGqqqu/Ha6sp4cv/QabO+P8+hGtnE4nrLv7/zezjxX75KuFOEQ1OqDMT7cxcy0cTBrK3EnUK8E68dpoRBrMufHHtVMZ+eYUweCnVfKcaCvU1YHWhqmvkPL74ZeWFiOxKsA6elNjWWESmE77g41R1UqLILcDPgM05qlNgrIhMEZHzE9t7EDyg/hQTOdwjIq1ylB9CeOwLlQWr45sJfykXAstVdWyUZwCHi0gHEWkJDCa3qV52H22ghn0EdCU87uTqoynAKqOPcq3XUmAswX64XWL7lj4C/kFw+MzVRycSBqZQWdU+eg3YnNJH3wK6x/PI7pe1hP5LnkfmXHcl6xpJbO+XLJd17UyK55Upl7x2WmTVmbx2rkpo2dfOQ4Q7qux2DiH42h9AGIyyr52NhN/jpES/VIjIG8D/EKYKcvVLjb9T5UBZTrqLSGtgDGE+ZEVmu6puUtU+hF/wQSLSO+5/AmG+ZopR5WGqeiAhq83FIjIgbm9CuMjuVNUDCPMKV2a1pRlwEvBoYls7wl+3HgSr3lYi8r3YxtmEx6exwLOEW/oqd0nFIPbRnYQLOlcfHUK4A6xpH50NtM/VR4Q72M3k7qNjSAxYWX30NaCR0UdjCQPLv5LnkTjHS4CXszXCtb3VNRLLPEEY6IZmtES/9CTMkf1WVVdk9UsLgvdZss5Mv5wKXAPcHbVkvxwODAT+ndWWzLVzaqbOrH7ZEziSML+1ItEvz8a+HAV0zPwOvwzU14A1n6p3FV3jtmoRkaaEC/EBDUkVtiI+kownTEIC9AdOEpEPCbfKR4nIXxL7z4//LiZczBlXznnAvMRdyGOEizDJccBUVU36sB9DsOatVNUNhHmIQxPHG6mqfVV1AGHS+d0cp5HdR02peR/9lfDol4sVhAE4Vx/dDrTO1UeECfSV2H20gtx9NIOqA/OWPiLcRawgq48IczOfE+ZjnonSfKBb4hxnsbVHexNC31W5RmKZxwkDzx3Z10/U7wUmEPo7u1/+Sni8PyXZL7Hc/4vtWJvsF8Ij7xjCfFbLrHaeEP+9L9GWYwiPo8sILxCeJjwWbumXxLWziPAH71gK+E6VE/U1YL0O7CkiPeJfmSGEDCqpRG/ykcBsVf1dllYhIm3jzy0ITpZvA6jqVaraVVV3jcd6QVW/F/dtJSJtMj8TJmZnxHKfAHNFpGc8zNFs7SZ6BonHwcjHwMEi0jK2+WjC3EamrR3jv90JF//otD4ifHl2qGkfxX9z9hGwHWECPlcf/ZjwuLhVHxG+7K2x+6g1ufsou+1b+ij+vxVb99FIwpe+K1/00euEO4+HCW/Euibrjud/E7AueY0k+qUr8GyWlumXkYQ/HlX6hTAQvER4HB6bo19GEibx22T3C+HOezZhoJ+VOKYQJvKnZl3LHxMG6lGx3IasfuklIm3jtXMq0CG2tVbfqbJD62m2nzB38y7hzcb/ZGkPEp7fNxAu2PPi9sMIcwZv8kU4weCofYUQsvAm4aL5hXHcI0m8JSS8VXmDL8IhstvSB5gc6/0r0C6htSJciDvkOM4IwoU0A/gzsF1Ce5lw8b5BGMys8x1MuPPYSLg7SWoPEt4cavx8SnhrmOmjZbE+JZj3nZfoo4xWpc5EvUsJj3bzYrlMH32Wq1zso0zowmbCX/bzEn20jvCmK/v8RsTz25A8XtTeiG1fQ5iUTv6uL4/authvSW1cok/WEwaAwYl+ydS5hnDtDY798m7U1uaoM1P2fcKj2PRYbrfYNqvc2YnjLSdMrGe0Y6I2g62v5XsTdX4a+yKjTUm0/wMS1zkp36mG8nE/LMdxyoaynHR3HOfLiQ9YjuOUDT5gOY5TNviA5ThO2eADluM42xQJ2dBfE5E3RGSmiIyI23uIyCQReU9EHo7hGKnU+4CVtRTGNddcK7JWAqwjLBDfnxACc6yIHEyIl/u9qu5BCJk5r9qa6juugiw7F9dcc624Wil9CNH+UwnLsZYATeL2Q4Dnqitf73dYjuM0fLIWl48jBLcu0+A6AXku1i44cFREjgVuJVhw3KOqN6bt31JE2yb+v5ovFlh13rvqvpXLoCLuvLZV1cfbzyo30a6icayj6hKtFZXr2L5iOwDav1XVDaZyE8RiLN5vxyraqsq1tK5oDkDHRUuqlvscKqIHwfydOlXRVleupmVF9jKxrbUu86smXK5cDZlia7psZ57fHHpU0TZVLqNx7Jg2rKqiran8nBaxoV2ylpItrVQ6VAgAb2zav2pDly6BDqE/GjWuamihlUuRig7h2CuaVi23vBJ2qACg7+qq66aTfbZi59ZZxTawQ0WoqwVrstq5mQ4V4W/pPLpW0dZWrqJ5Rairx9sfVdEqN0CskiV7t6+iraxcR5t4TXy0qmp/sqwS2oZz6Nl6dpa0kbYVIUH6OxuyLtClldAhnvvmqVXb8ilUZJqQ5Q9S+RlURL+LFS3sftn+/aq/28p1EE+BDVkmR58SDLcyvAlLVLWCWrKHiK7Oc9+FYYXI2sSmu1X17lz7xiVQTxAWio/S8DiIiHQDnlHV1IXcBaWqF5HGwB2EdXvzgNdF5ClVzV5LtoW2wAWGdu0o+1hvf62TqU2hn6mducsYU7t18immNvR3OfsbgCsvtc4gnRuvHG5qb964i6mdadpTwYCt1v1+wa83XWVqO3/+gqm1bG1fqkuft/8ITp4ipjbuqj6mth9vmdpl/MLUHjjEnvL448TjTe38l+z+HDnAvpYOW/CKqU1e19zUzOXnwLj97X75+qkTTG1RzmX/X7AzfJS+RzprgIvz3PdqWKuqdsclUNVlIjKe8AjYVkSaxLusvBZrF/pI+KUwDXOcLxvCFy6Y1X2qrSu3McFsgqPKaXG3/wKerK6ugu6wyG0a9rUcDT4fOB+C5YDjOKWNUPjgkKATcF98ImsEPKKqfxeRWcBDIvJLwqL8kWmVUNw22cTn2bsBOov4amvHKXEyd1jFQFXfJDiqZm+fwxe+anlR6ID1pTANc5wvG40IxmelRqEDVtJkbj7BNOy7aQU6721Pro842C537Sn2HOLoMWfaBf/PloZ9fIepnXXpn03txl8Ptyvta0u0t6WvnJzLdDRwzpOjTO3SUX+wK73Wlt7+aG9Te5nDTa3joEWmxihb6oflvAynYM8gj3/lOFP7ZKI9wXAzl5vahJSJ9WsZYWofde5uamdhv6Q5JCWmc4KeY2qdxyw0tX0fe9/UADg9Xa6OIj8SFo2C2qSqG0XkPYLb4iZCSqOZRWmZ4zj1RjEfCYtJMQJHf0l4Dv2Pqv6qCPU5jlPPZO6w8vnUJQUfT1VfkpA2yXGcBkKp3mHVyQCZDGvovnNdHNFxnEL4Ug9YybCGfr08rMFxSh2hYb4ldBynAdIg3xLWhrWtmpnrAtNCF0akrJ36Fk/Yov0mmk7d55pau4VrTI2jbOmVr9lxDf1fsF/t356SQe4mrjC1C864y9QmnrPZ1I55abGpnbnWXn/51KBBpsbntnQ1vzS1tiwztWn9e5laM9abWtr6xGW0NbXdec/Uus+1+6xFN3v95Q/7mxIXvXSuqf14wO2m9vZp9trTQEFLCUv2kbCgt4Qi0k1EFhHCGvYVkWUiUr0Jl+M4JU2pviUsNKxhI3CcqjYFtid43UwsuFWO49QrxVz8XEwKDRxdSMh2i6quFJHZhAXRpr2M4zilT0NdmrOFGIt1ADAph7YlrKFz98bFOqTjONuIUp10L4pFsoi0BsYAw1R1K7syVb1bVfupar+Mi6bjOKVLg3wkBBCRpoTB6gFVrcYH0XGccqBU3xIWapEsBNOt2ar6u3zKrKalaWmc5rqQFrrwlMw2tf2W221Z8MZutpgy9MrJKbGvtqMveqJtIdz7xkNMbbHYN8KtH9hkas/oQFPbiftM7V7sV+0nn/6cqeld9vn94dxL7XIX2OXkH3Zf6w/tco/d8D27zqUpdXZOacvpKeUW2OVunDDM1NjHlg5+fLqpXbL3b+yCAPysGr16GuIj4UDg+8AFIrJGRD4RkcFFaJfjOPWIAE2b5PepSwo93Higjaquio+GEwgJPBzHKWMaNYIW21W/HxCCm+qIQsMaFLbkmcrMwflaQccpc0SgSQk+ExZj0r0xMAXYA7hDVVPDGjp0L8XoDsdxkmQeCUuNgsMaVHWTqvYh+LkfJCJbJUJMhjVkEpw6jlPCCCE1cj6fOqRoqepVdRlhTuvYYtXpOE49UaKLCQsNa6gANsRsrpkEiTellWn/1jI7G3NKwog014W00IVfpSRCHKh2hqEBp79masNHXGlqCw60M1RvsE0CGDgtZQnmfDvj+EWd7WiSY1/4l6ntfpSdxOCxLbktc3CNLZHyfvj6aXZSiBSDBE697gFT+zQlUOjTDV1t0c73Ae1saa+KN23xGVu68olbTO2dWT3tgimGDLcdkR62YPs85EkRQ91jGvr7gZ0Ic9x3q+qtIrI/cBfQGvgQODNX4HmSQpuUM0FigXU6jlPfCFC82ZuNwGWqOlVE2gBTRGQccA9wuaq+KCLnAj8l/U9iYY+Eqvqmqh5AWEO4gRomRXQcp0Qp4iOhqi5U1anx55WENPVdgL2Al+Ju44BTq6urWHNYQ2MjHMdpCGyjOawsk4SZwMlROp2qSZlzUvCAJSJdgeMJt3fWPueLyGQRmVxpryRxHKeUyP8t4Y6Z73f85Mwcm8Mk4VzgIhGZArSBFAvZSDGm1W4hLFxqY+1QJQnFdp6EwnFKnppNui9RVTudNrlNElT1bWBQ1Pci3PikUqhF8gnAYlW1zcodxyk/MpPu+Xyqq8owSRCRjvHfRsDVhDeG6XWF1TW1Q0R+TVj8vBFoTrBJflxVzaXy3ftV6GWTT8mpDfv4DvNYaQkj0lwXXtrffg8wXuzQhWv/Y0rIgJQ+W2hLD+q3TO0Hq/9oaqta2f0C15pKh03zTW1p2y52lXZTYIjtEnCN2vEl13e2k4LrphSHhMWfmNrmpXaSS7FzVyCLb7Pb0uESu9xS+/z0CDt+RkavtRvzPfsbr7fZ/XJS70fsOoG/ybenVHfXk0a/NqKT7XwqVZAXST2WiBwGvAy8BWQyo/wc2BO4OP7/ceAqrWZAKnQt4VXAVbFRRxJeUdq+Ho7jlAdFjMNS1QmxxlzcWpO6irGW8ENgJcECesdC63McpwQoUY/kYjVpoKouKVJdjuOUAiXoZl4nY2jSraFd99Z1cUjHcQqhRO+wihE4qsBYEZlixV8k3RpaVzQvwiEdx9mmNKJobwmLSTHG0MNUdX58RTlORN5W1ZeqLeU4TulSondYBYU1bFWZyHBglarebO3Tr5vo5J/k1j671Db3a7dwjX3gtOiNB1O0p21pxJ629nLK+u5ZKVkFFuxjh1+sm2wfr/kv7d9Rrxum2W1ZeKCp7dNpqql9h4dNbfjdN5qa/st+Db/H6Bmm9t7vtrJQ20KbCytNbWUv28WC39qSrE5JJvGRfQ4th9ru36t/0d7UPr/JfpBp//lnprZuvxSrETsXCACyX3qoQXX06yA6udowznisPxd2rJpQaOBoq7j6GhFpRYhata9Mx3HKgxI18Cv0pm8n4Mnod9MCWB4/juOUMyX6SFiovcwcYCohYHQ7gk2yuzY4TrnTQB1HdwAGAGcDqOp6cqy4ToY1dE9xdHQcp0QoroFf0Sg0rKEHUAn8SUSmicg9cS6rCsmwhoqtVMdxSo4SvcMqdMBqAhwI3BmdRz8HbMNzx3HKgxIdsAo93DxgXiIX4WNUM2DN36kTV156QU7txl8PtwseZUtysv2aOi1hxPDO9iv6o/Ufpna4nGBqz/c3JR6a9U1TGzL2r6Y28YYDTG2TTLcP+KgtzdreDnkYMegKU9M59mv/N0fvZWovcoSpyXspYQY97OMduugFU5t4ykC7zi4p7hBNU9pyT0q5TinlbrXLtRi22tR6zn/D1K5Iz/UCjK5Gz4MSXJpT6KT7J8BcEcmk/jgamFVwqxzHqV8a4h1WHKg6AtOiSVcT4BfFaJjjOPVII4LDXYlRqB/WO0BP2JKyfj5FuRd1HKfeKcFHwmLe0B0NvK+qH2ULybCG7bunLDdwHKc0aIiBo1kMwVi5lwxraFnRsoiHdBxnm1Cic1hFGbBEpBlwEqnvpRzHKSsa4FrCDMcBU1V1UUG1pJjev/K1FDFlmn/BgZ1sMSVhRJrrQlrowohXbG09+5laz0HvmNrBf7FDF1akTIx+cpr9+D0dO1TiD1xkateebr9On4y9YP/ckSlTmykJPyYttrWldLDFf9sSD6VoV9vSx7/raIt2l8GZtnRSs7+Z2lPrTzS1ts1sl4eiUKKPhMVq0hmkG7k4jlNOZAz8SoxiJKG4AjgN+IqIDALOUdWUvEaO45Q8JXqHVagfVhfgQqCNqu5DeKIdUoyGOY5TzxRp0l1EuonIeBGZJSIzRWRo3N5HRP4tItNjins7iWiiSYXSBGghIhuAlsCCHA32sAbHKSeKe4e1EbhMVadGw88pIjIO+A0wQlWfEZHB8f9HplVU6NKc+cDNwMeEKezlqjo2x34e1uA45UQRHUdVdaGqTo0/ryR45nUhJLDZPu62AzludrIpdGlOO+Bkgs3MMuBREfmeqv6lkHodx6lnttEclojsChwATAKGAc+JyM2Em6dDqytfaJOOAT5Q1crYmMfjQc0Bq8v8hdx45fDcou3jT/8Xppianmivht/Q1q5zoP7L1IbsY7snpLkupIUuNJPrTe2AG0wJOdV2AvjNmktM7afDbje1V2+xr40/8gNT27WvbSh7j/QytUv0N6ZWeV4bU7sOu9w7397f1GYu2N3UemK7IKwal3LLsA6GbndbTum8P9zBPQt/ZJc1EqXcP/Z89hmROyFI52YLmfVD21WDr9lSwdTMwG9HEUmmUblbVe/eqkqR1sAYYJiqrhCRXwI/UdUxIvJtYCRhTDEpdMD6GDhYRFoCawjLc1LyvzhO+WINVkCtBivAHKyA+husoKZ3WEuqy5ojIk0Jg9UDqvp43PxfwND486PAPdUdqNA5rEnAXOAzwoC1D7DVyOo4TplRxKU50cllJDBbVX+XkBbAFqO0o0gNIQ4UOofVm5B4oh3By/1ZoBvwXiH1Oo5TAhRv2U1/4PvAWyJbHCd/DvwAuFVEmgBriZEEaRT6SNgLmKSqqwFE5EXgFKg68VAlCcX22VU4jlNyFHHSXVUnxBpzkbLmbmsKXfw8AzhcRDrEeazBhDusKlRJQuFRDY5T+mQM/PL51CGFGvjNFpGbgLGEBBTTgU3FaJjjOPVMCRr4FWwvo6ojVbWvqg4gTL6/W3izHMepV0rUD0tU7RifvCoQ6aiqi0WkO+FO62BVXWbtv2+/5vrg5F1yal852R7rbn/KbkNvPcTUBn4w0dTa7FRpakuoMLXtJthtmTbIjkU64Nd2DNOIn9t13qP2y5OPP93T1P6Z4r5yTEp/psXfjB+U0tdn2X19yf12PNVbKbFrd3Khqa3Ejt+6ml+Z2i+4ztTSrHUemHGeqZ3Ze6Rd7jC7nFxnf//GHnW4qe1Eiu8OsL+8O6W6UIM0+u0nOvmJ/PaVPSnoWDUhrzssEblXRBaLyIys7ccCH4rIekLk6sVpg5XjOGVECRr45ftIOAo4NrkhJp24A+gNtAYWk2qJ5zhO2VCij4R5HU5VX4prgJIcBLynqnMAROQhwrpCz0voOOVOiab5KmTSvQshyj3DvLhtK0Tk/Oh3M/mzSn+J6DjlgDbO71OX1MkNXVwIeTeESfe6OKbjOLVHBTaVoONoIU2aT9Ug0a5xm+M45U4DHLBeB/YUkR6EgWoI8N3qCs2hB2dyf07tnCdHmeVu4gpTWywpT7bz7fAEWt1oSs2vsm8EJ95gZ5xJy3CTZhPT9Sp7+eV/ix26sHqtfe6/0OdN7euzBpoaJ9jSRXN+Z2qnj7KzvC05p6upyai5prb39K3y8n5RbmbKzfqZd5rSfmrbA42Wnqb2itphKR/J3qa2u75vak/yDVMbdNPLpjbsCvvaDVxVjZ6OCmxsnO+M0eaCjlUT8g1reJAwmb6viGwQkfNUdSPwI0J0+2pgZ1Wdue2a6jhOXaEibGrSJK9PXZLvW8IzRGQAsAq4X1VHxu1Pi8iJme3brpmO49QlirC+cbM8916/TduSJO/h0QhtMLc7jlO+KMLGElxMWCf3c0l7mabdd66LQzqOUyCbSjAxYcGLn/MhaS/TuCLFZN1xnJJAETbROK9PXVJ6Q6jjOPVOZsAqNep8wGrDKgaQ+3XtpaP+YJa74Azbyb/1A3b0/EWd7dfwf+BaU+t1wzRT2yR26MKKlOUMaRluLv/UznCTFrpwc3P7lfKrs4+yG5MS7fGdOaNM7Y7rLzO1wdc8bVdqmzXAUjvk4fk+drG99E1Te/c9KaeHAAAblklEQVTmH5razZfbWZZ++y07VOLDZ+xy8hO73HUP2SmR9hgyw9Q6/NQOa+zJO6ZWDBRhHflOutcdeT8Sisj7wPuE0IZ5InJeTEG9iGAev6+ILBMR20vDcZyyINxhNcnrU5fU5Gjn8EVYQ28AEekEHJdMQQ3YpkiO45QNZf1ImCt8QVUXEi1lVHWliGRSULtjg+OUMQ1+DisrBXW2tiWsoXX3dsU6pOM42wiFhhuHlZ2COltPujV07NfN3Rocp+SRos1PiUg3wkqYnQhj4d2qequIPAxkFm+2BZapasprliIMWEYKasdxyhhFWF+8t4QbgcuSc90iMk5Vv5PZQUR+CyyvrqJCMz9bKahNujCfX28yVpLbUQZMPMd+ff+M2s4Dx77wL1N7eNN3TG3WwgPtxtimBHxy2g6m9tNhdujC87fadaa5LqSFLoyw82FwrR2ZwcOcbWrjrjnM1I4/5p+mduvzdlLfE5+yO7QPre228HVTu3Cq7dbwEgeZ2tnY4TPYOUs4/zj7F7jGPgXmDtoqjecWXm//VVPruTo9OZUd1JEfxZzDqm6uO44j3yakq0+lJkko1hLCF3rGsIYLgJmEFNQXisgnIjJdRAbX5qQcxykdMmsJ8/nUBGOu+3BgkWpKeqhITZJQHAq8o6pNVbUrYU7qIFUVoA3wEXChqqZEDzqOUy7UIA5rx4wFevzkvJ1Omes+A3gwnzbVOgmFhoSGq+J/m8aPT6g7TgOgho+ES6rLS2jNdYtIE+AUoG8+Byp0DqsxIVh0D+AOVd0qpCHutyWsoVv3Qo7oOE5dUMylOdXMdR8DvK2q8/KpqyC3BlXdFF9DdgUOEpHexn5b3Bo6VNjrsRzHKQ2KvDSnP2Gu+6g4z52c6x5Cno+DUKQ4LFVdJiLjCclW7dWcjuOUBUV+SziBkJo1l3Z2TeqSMBWVx45hDuvviXWEFcCGOFi1AMYCN6nq39PqaXTAAdr8xRdyam9vbxv5d39psantNMBOVLA7dgKAidvb4RC9VthuDbPG2iEPzw460tRe5VBTu+5v9op+2d3+HWnaHesCWxqREp53rR2dAHYeBu67+NumNnS9/dp/+fG2qeOccZ1Nbbc5KUmahtv9ov1TXBdSElt8cLtdrsf37XKV97cxtYqKlabGc7bU8cCPbRFYLLtMqW5eKY3d+rXT6ycfnde+35MxBR2rJuR1hxWdGnYFGonIPELE1OvAfXEea0/g4+oGK8dxyoNyt0jeyqkhcoCIXAr0A7YvduMcx6kftIhLc4pJXpPuqvoS8Gn2dhHpChwP3FPkdjmOU49klubk86lLCh1CbwF+RggcNUmGNUg3213ScZzSoFTtZWod1iAiJwCLVXVKdfsmwxrosGNtD+k4Th2yLZbmFEohd1j9gZNiPEVzYHsR+Yuqfq84TXMcp74o1TmsWrdIVa8CrgIQkSOBy/MZrBo13kzL1qtzai9zuFnuzLVjTO1ezjW1xzjN1Cb+0Q5r+A4Pm9qIQVeY2h+4yNT+yA9Mje1siRNsKS1hRJrrQlrowoiUt9lL1c4mcRqPmdryK1PyUf7Slk7Bdizqu9urprbo/p1M7WJ+a2pdec/Umt9mO3GQYp4wia/Z4hI7kUbvA9fZh/u0p6kVg1J9JKx1WIOqjhSRD4GVQAvAn/UcpwFRtgMWdlgDwEBVXVLcZjmOU59sphHrUm/764dauzU4jtOwKcU7rEJT1SswVkSmWB44EMIaMl45Wrm0wEM6jrOtaaip6g9T1fki0hEYJyJvxyDTKiSTUDTpt797ZjlOGdDg7rBUdX78dzHwBKSYZTuOUzZsK4vkQqn1HZaItAIaRVP5VsAg4Lrqym1a0ZSlz3fJqXUctMgs99SgQaZ28ukpy9qvSWnMEPsV/fDlN5qazrFX7V97+k2mtmtf2+qgzSA7NOOiOXZ+jzuuv8zU0hJGfH32BFNLC13oID8ztV/pk6b27d/dZ2qPiP13btqth5iaLEm5Wb/eDhfYHbvP5qWEguz8jZSkLs+tMaWWmjuMB4DjvmJKM+wIGbjaDnkoBmFpTplOuhtuDeOBJ2POsRaEFD3VpulxHKf0Kes4LIywBhGZCtyqqveISDOg5TZoo+M4dUxZ28vkCmsQkR2AARDuoVV1PbC+uM1zHKe+KMWlOYVMuvcgpJb8k4hME5F74lzWViTDGlieko3ScZySoFTDGgoZsJoABwJ3quoBwOfAlbl2rOLWsENFAYd0HKcuKNUBq5B7vnnAvERqr8cwBizHccqLYqb5KiaFuDV8IiJzRaSnqr4DHA3Mqq5c39VTmDzFCAsYlVLw85S23JWSiGGwLf1Cf25q133XrvPN0XuZ2mRsL/57pJepDfy+KXH6qEdNbfA1dqLt44+xLRnue95OGJHmupAWunConGxq2+kwU7tY/2Bqd3Ceqen19u/ocT3O1EbxiKm9xhGm9r/82NTmsLup9Th5oan1efrfpjbtDjukg7G2BEaKmhpQqvYyeT0SxrCG94F9RWSeiJwnIj2BjsA0EVkD/A85bJQdxylPivVIKCLdRGS8iMwSkZkiMjSh/VhE3o7b7QDASKFuDT3jQRsD84HRedbnOE4JU+Q4rI3AZao6VUTaAFNEZBywE3AysL+qrotL/FIpllvD0cD7qmonCHQcp2woZhyWqi4EFsafV4rIbKAL8APgRlVdFzU7+WikULeGDKnpppNhDZUpc1GO45QORUxVv4V443MAMAnYCzhcRCaJyIsi8tXqyhc8qxYj3E8i2iXnIunW0K+ruFuD45Q4m2lUkxReO4rI5MT/747f+SqISGtgDDBMVVeISBOgPXAw8FXgERHZTVPS0RfjNcBxwFRVtVcuO45TdtTgkXBJdanqRaQpYbB6QFUzRv3zgMfjAPWaiGwmWK2b0eXFGLDOIOVxMJsVO7dm3FV9cmr9sDOGXZ2SqeAP515qatdPu9zWOv/K1EYv+K6pvZjy6vvckfZ7h0tSXBCeMhVYck5KLseU9yq3Pm96KnLt+hGmlpYwIs11IS10YZ3cYmpHXGXfdOtZ9gv6b11jX3Z/3WeIqZ096y5TO+i7b5la79Gvm9rP7r3N1OY82dnUpk2xQxcOvfgFU3u18ihTKwbFDGsQEQFGArNVNWk98ldgIDBeRPYCmgGpduuFJqG4AjgN+IqIDALOUdW1NT0hx3FKiyK/JewPfB94S0Smx20/B+4F7hWRGYR1yP+V9jgIBYQ1iEgX4EKgjaquEZFHCJPvo2p4Mo7jlCDFGrBUdQJ2LGuN8pgWGtbQBGghIhsI1jILanJwx3FKk80NcGnOfBG5GfgYWAOMVdWcCwZigorzATp2Lz0XQ8dxsinjpTm5EJF2hCjVHkBnoJWI5Ly9S7o17FDRtLaHdBynjihVt4ZCAkePAT5Q1UpV3QA8DhxanGY5jlPflOKAVcg938fAwSLSkvBIeDQwOb0ItGAN+5H71fEpPJ5zO0BblpmaXpCyNv09W7p6029t0c77gLyX8iLjP7ZUeV4bUzsd25FBRs21K11qhzyc+JRd5/Lj7dCFlAiS1IQRaa4LaaELw39t//7k2ZSQhyEpv/f2Z9h1ygemtrqVXWfL3ey2PPWBXe4H5/7d1P4583hT0z/bdT5wy6mmFhhTjZ5OWVskp4Q1zAU+IyRUfYs4T+U4TnlTqvYyhYQ19Aa6Au0IMRTPAt1IvadxHKccCGm+Su8tYV5zWDGbc7bXVS9gkqquVtWNwIvAKUVun+M49YAibNrcOK9PXVLIPd8M4Fci0oEwhzUYYw4rGdbQpXuhXoiO42xzFDZuLNM5rFyo6mwRuYlg1vo5MB3YZOy7xa1h/35N3K3BcUocVWHTxvKdw8qJqo4kLGpERG4grL52HKfMCQNWA7rDAhCRjqq6WES6E+avDq6uzDy6chm/yKmNf8VOHDCtv53AQf5h37Sdet0DpjZm8TGm1vpC+5elPezH2kkpnonXpVgr3MmFprb3dNvI9fncxhcA9KG1qa0cN9TU0sJLpt1quwukJoxIcV1IC10YPs0u1/bfdnKH5b+0wzb0hpS2PJoSRjEmpdxlKeWspCvAfme9Ztd55AZT+w22Q0lRUMpzwBKRboS5qfZAExFZRvBnHgm8ICJ7xl3vV1U7WMpxnLJBtRHr15beMrp83hJuBI5T1abA9sBiYGJMPNGC8LawDXCQiOyzzVrqOE7docDGxvl96pBq77BSDOR3AN5T1TkAIvIQYW1htbkJHccpcVTqfDDKhxqtJcwykO8CJNeMzIvbcpXbkoRibeWq2rXUcZy6Q4GNkt+nDsl7wMo2kK/JQZJuDc0r7Ilgx3FKiI15fuqQfNcS5jKQn09YipOha9zmOE65sxkoQbPzfN4SWgbyrwN7ikgPwkA1BLAzN0R6vP0RDxyS+/X3JxN3MMs1Y72p6Q/t29JPU+y3Hl1qa2JHUXDoIjs5wFI6mNo7397f1F5/pLepyUz7lfle+qapjePrprbbHPtvS9/dXrXbsiTl9f31tUsYkea6kBa68JPtOpnatXNMCXkw5Rx+nhK6MDelXAe73JK+9lNFT941tSt3sZ0xzrzJlAD4WbpcPQrYURX1Rj6PhBkD+f8WkbUiskZE7orrB39EiHBfDeysqjO3YVsdx6krlLBuJZ9PHZLPW8IJItIZ6KSqU0WkDTBFRPZR1adF5ESik8O2bqzjOHVIHc9P5UO+SSis0IZZKQkqHMcpV5SSHLBqbJGcFdqQb5ktYQ2VJfhc7DhOFpkBqwhvCUWkm4iMF5FZIjJTRIbG7cNFZL6ITI+fwdXVVaO1hLUNbUi6NfRrLe7W4DilTnHfEm4kLOdLTimNi9rvVfXmfCvKe8AyQhscx2moFOmRMGVKqcbkG4dlhTbUmCV7t+ePE3Mb79/M5WY5K3EFwGM32MljP91gJ2lov1OKG84dtjTxlIG2+G9bmrlgd1O7ml/ZBc+805TevfmHpnbhVLscw+3X8Ivu38kud70dRvG42m4bf91niF1nSsKINNeFtNCFEbvZGsw2lVfSiu2XotlRKQzEDoOZceBXTe2xCfZ1fWb/kSmNAa6wnTPyYhuFNWRNKfUHfiQiZxEMFi5T1c/Syuc7h/UtQmjDhTGsYZ6IDI7PposIuWL2FZFlIlJgTzmOU+/ULKxhx8wcdfzkTEaTY0rpTmB3oA/hDiwljVUg30fCiUDf5DMo8CFfODkkt0/Ms07HcUqVmr0lXKKq/dJ2yDWlpKqLEvofATsfWqSgsAZVnZVrO+7Y4DjlTREn3a0pJRHpFMcWCE9xM6qrq8aOo1ZYQ1q4QzIJRfvurWp6SMdx6oPixWFlVsu8JSLT47afA2eISB/C/dyHwAXVVVSUsIbqwh2SYQ279uvgYQ2OU+oUMXBUVScAud7yPF3TugoOa/BwB8dpgJRopHtBYQ21CXf4aFUPzn8p97LDCQPsebtltLXbtzTlps1eDA+Lb7e11T82Je2SYlr2kC315A1Tu5dzTW0//aap3Xy53ZaXsFf7P93/VFO7OOVlze5cZmqjeMTUzp51l6mJfGBqqQkjUlwX0kIXhmM7eR/Gcrst3WrXlgmj7ev6sN1zpvIMdY6x62zWq7q47dIMayiUasMaYhKKKXwR1jAvhtGfDMxMbP8k3/B6x3FKnBJ1a8g3CcV/q6oAFQQrme8CTwEHxe1tgI+AC1W1xs+ljuOUGEp4S5jPpw6pdRKKGNKQMWhvGj8+oe44DYESncMqJAkFItI4vqZcDIxT1ZwODkm3BpZVFtZix3G2PZk5rHw+dUhBSShUdZOq9iH4uR8kIjlXVCWTUNC2ohjtdhxnW1LGc1jVhi7EjM/jgWOL2zzHceqNcsyakxLSUAFsUNVlItIC+DpQjTU+9Gw9m5FG+MK1jDDL7c57pqadU8IM2qU0xs4XAR9dYkrSNGWq7mpbWjXOTkx5/nZ/NLXR0tPUfvstuy1nY4cSHDHT9l/smtLX8zjb1F7jCFM76Lu228bqVinhAo/WLmFEmutCWujCcOxEKLImpS272W0ZzBi7MW1sSXe16xzR/gq7IDA8Vc2DEp3DyicOK+PUsE5ELgCWEpbZzAPuiynr9wQ+VtVqFy86jlMGlGiar3weCTNODc2BjoSwhg9V9U1VPQAYBTxBsJhxHKchUESL5GJS7YClqgtVdWr8eSUhhLgLgIh0BY4H7tmWjXQcpx4owQGrpoufd6WqI8MthJyNKU/iVd0adurerKZtdBynrinXpTkZssMaROQEYLGqTqmubDKsoW1FjR1tHMepa0o0rCHfxc+5whr6AyfFtYPNge1F5C+qahtRO45THpTopHutwxpU9SrgqrjPkcDl+QxW72zYm8MW5H7p/FHn7ma57nMX22083X7dvFeFnTTh3aW7mlqLofafDr3Hft388e86mtpQbjO1B2bYq+tfUdt54MNnUkI6UhYVXHu7nbyi+W32q/2dv2GHBPwvtsNF79Gvm1rL3VLCBcakhDzMTQkvSUkYkeq6kBK6MHx8SrlD7HILU5xGOk0+xa7zm3adHfnY1ALVRhilU6KPhLUOa4hp6j8EVgItgB23WSsdx6l76vhxLx9qHdaQ0Aeq6h6qav8ZcRynvCjRsIZauzXgiSYcp+FSopHuBbk1EE5rrIhMsXKRxXJfuDUsdbcGxyl5StStoSae7rkSTRymqvNFpCMwTkTeVtWXsssmk1DI/n3dM8txSh0F1tV3I7amILcGVZ0f/11MWJ5jm4g7jlM+lOscVopbQyugUZzXagUMAq6rrr6+m6cyeV3znNpZ4SYsJy26rTY1XZDyav+ZlMbY5gLwC1uSTik3iRfZkl5jt/PM3iNN7SPZ227LT+y2nH/craZ29/eHmlpq4o7n1pjSHHY3tZ/da4d0PPVBSrjAZSkhDx1Sfu85ndlinSkJI9JcF9JCF4ZPtMtdxZ12Y2xDDfS6FDeKk/raBYHDUtU8KNGwhnzusDJhDReKyJqYhGIwsBPwqogsAz4F9oAU3w7HccqHIka6i0g3ERkvIrNEZKaIDM3SLxMRFZFqQ6MKcWuYA0wlBIxuR3AdtSMcHccpH4r7SLgRuExV9wEOBi4WkX1gS1auQVBtJCxQgFuDiOwADCA8LqKq66PzqOM4DYEiDVhpji/A7wkGCnm9jCvErWE3wuKPP4nI/oTchUNV9fMc5ba4NXTvXJMjOo5TL9RsLeGOIpLMCHt3jAzYiuQYEnObzlfVN8JUefXUOqxBRJoABwI/VtVJInIrcCVwTXbZZFhDv/3Ewxocp9SpWeDoElW101tHkmNIrP3nhMfBvCkkrGEeMC+R2usxwgDmOE65U+SwhhxjyO5AD+CNuCa5KzBVRHZOq6cQt4ZPRGSuiPRU1XeAo8lnuc5mYEVu6RA7WJ4f9rervHHCMFO78olbTE1G2/e8qzq0NDW9NeX29UxbSnuF/cA/bbeG3fV9U7vuoRtMbU1r+3i/WmVny5jE10ytpdrhJT1OXmhqc5605wJ+cK6dCkCn2H29pK99ggN5wdQmjLZvBtISRqS5LqSFLuwqtjPGcWq7e8x81g4T6f9MtTZ0hVHEsIZcY4iqvkV4iZfZ50Ogn6ouSaur1mENItIzHnCaiKwB/ocQ3uA4TkOgeAZ+/QljyFEiMj1+BtemSfnMYWXCGqaKSBvC5PqH8a6qJ4QM0MB8YHRtGuE4TglSpNlmVZ0ApM6qq+qu+dRVUBKKBEcD76vqR/kc1HEcpzYU6taQYQjwYEq5LW4NlZ/VtImO4ziBWiehSGxvBpwEPGqVTSahqEjLxOw4TolQmv4yhSShyHAcMFVVFxW7cY7j1Bel6eBX67CGBGeQ8jiYzYoWrRm3f5+c2gQ9xyx30Uvn2pXuY0vvzOppi0dtZ0rtH7OfXVsMs1/tn9Tsb6Y2ub/9On32EQeY2pN8w9T2GDLD1OYO6mZq6yu2NzWW2Ik7OO4rptTn6X+b2rQph5jaP2ceb2r7nfWaqfVMsZWYceBXTe2w3SebWlqGzbSEEWkhK2mhC1+TU03tR2pbjTx+XEpbALBdNfKjNO0aCnFrQESuAE4DfioiD4pIbt8Yx3HKjM2EQS+fT91Ra7cGEekCXAi0iauwGxMm3x3HKXvKdA4rJQnF8li+hYhsAFoCC7ZhWx3HqVPKcA4rSTKsIS6AvpngY7MGGKuqY41yW9waOna3540cxykVyncOC8jp1tAOOJmwgLEz0EpEcmZ+ToY17FDRtBjtdhxnm1Kapu6FuDUcA3ygqpWqugF4HDh02zTTcZy6pUznsFLCGj4GDhaRloRHwqOBlPfFge3fX8XXT52QU+s8xl7t/+MBt5vawY9Ptw+4iy396R/2Knr2s6We898wtafWn2hqy4Z2MrVxd9ppAwbd9LKpdfjpfFN7vb39an//594xtd4H2vmdZqQk2Zh2hx26cOjFtnuC/jkl8cOR9hfiyl3sJE2PTch5sx/qHJOShGLXlLZ8M6VcSsKINNeFtNCFI+U4U2v3Z1MqEpm3hKVFPnNYmbCGdSJyAbAUOF9VnxaRucBnhOH4LeI8leM45U5pBo4WEtbQm2C61Q5oTXhraEcqOo5TRpTmI2Ehbg29CG8LV6vqRuBFoLrwW8dxyoIynnTPkOXWMAM4XEQ6xHmswRh3WFXcGkow/bXjONmU5h1WrZNQACtE5CZgLPA5MB3Df7BKEop2noTCcUofpRQn3QsJa0BVR6pqX1UdQJh8T0ty7jhO2VCmd1hpbg0i0lFVF4tId8L81cHV1bdhGSzKNqiJ7PuYnWzh7dPs+IRL9v6Nqd12xM9M7aTej5jaU89929Su4CZTa9ssxaHQzu3ATiw2tWFX3GhqPbHDE3qutv9+dDzQTrT77qcpDhdXpzzT51znEHi18ihTe+AW27HgN1xqamfavwbO7D/S1Jr1MrKgACPaX2FqHVOSE79yUl9TS0sYkea6kBa6MOL7tlYcyvct4UBCWMMF0a3hk4SB/Asisg74D/CyZ352nIZCad5h5TNgjSc4MrQAtgc+Aj6NiSdaEN4WtgEOEpEUZyrHccqH0nxLmI9bgwKr4n+bxo8CBwHvqeocABF5iLC2sPrchI7jlDhlvPhZRBqLyHRgMTAuZnvuAsxN7DaPrbPpZMpvCWvwxIWOUw6Ur4EfqrpJVfsQItsPilHueZN0a2hfm1Y6jlPHlOkjYRJVXSYi44FjgVeoGijalZBM1XGcsqc0HwklTFGl7CBSAWyIg1ULwsvrm4BnCXFXRxMGqteB76rqzGrqqyRM3DuOs+3YRVUraltYRJ4Fdsxz9yWqemxtj1UT8rnD6gTcF98KNgIeUdW/A4jIj4DnCH7u91Y3WAEU0omO49QNdTUA1ZRq77Acx3FKhRotfnYcx6lPfMByHKds8AHLcZyywQcsx3HKBh+wHMcpG3zAchynbPABy3GcssEHLMdxyob/D6C6KXziIAUvAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "DPP.plot_kernel()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "anaconda-cloud": {}, + "kernelspec": { + "display_name": "Python [conda root]", + "language": "python", + "name": "conda-root-py" + }, + "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 +}