notebooks/02.07-Fancy-Indexing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fancy Indexing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The previous chapters discussed how to access and modify portions of arrays using simple indices (e.g., `arr[0]`), slices (e.g., `arr[:5]`), and Boolean masks (e.g., `arr[arr > 0]`).\n",
    "In this chapter, we'll look at another style of array indexing, known as *fancy* or *vectorized* indexing, in which we pass arrays of indices in place of single scalars.\n",
    "This allows us to very quickly access and modify complicated subsets of an array's values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exploring Fancy Indexing\n",
    "\n",
    "Fancy indexing is conceptually simple: it means passing an array of indices to access multiple array elements at once.\n",
    "For example, consider the following array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[90 40  9 30 80 67 39 15 33 79]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "rng = np.random.default_rng(seed=1701)\n",
    "\n",
    "x = rng.integers(100, size=10)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we want to access three different elements. We could do it like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[30, 15, 9]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[x[3], x[7], x[2]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we can pass a single list or array of indices to obtain the same result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([30, 15, 80])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ind = [3, 7, 4]\n",
    "x[ind]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When using arrays of indices, the shape of the result reflects the shape of the *index arrays* rather than the shape of the *array being indexed*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[30, 15],\n",
       "       [80, 67]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ind = np.array([[3, 7],\n",
    "                [4, 5]])\n",
    "x[ind]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fancy indexing also works in multiple dimensions. Consider the following array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2,  3],\n",
       "       [ 4,  5,  6,  7],\n",
       "       [ 8,  9, 10, 11]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = np.arange(12).reshape((3, 4))\n",
    "X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Like with standard indexing, the first index refers to the row, and the second to the column:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2,  5, 11])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "row = np.array([0, 1, 2])\n",
    "col = np.array([2, 1, 3])\n",
    "X[row, col]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the first value in the result is `X[0, 2]`, the second is `X[1, 1]`, and the third is `X[2, 3]`.\n",
    "The pairing of indices in fancy indexing follows all the broadcasting rules that were mentioned in [Computation on Arrays: Broadcasting](02.05-Computation-on-arrays-broadcasting.ipynb).\n",
    "So, for example, if we combine a column vector and a row vector within the indices, we get a two-dimensional result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2,  1,  3],\n",
       "       [ 6,  5,  7],\n",
       "       [10,  9, 11]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X[row[:, np.newaxis], col]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, each row value is matched with each column vector, exactly as we saw in broadcasting of arithmetic operations.\n",
    "For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0, 0, 0],\n",
       "       [2, 1, 3],\n",
       "       [4, 2, 6]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "row[:, np.newaxis] * col"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is always important to remember with fancy indexing that the return value reflects the *broadcasted shape of the indices*, rather than the shape of the array being indexed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combined Indexing\n",
    "\n",
    "For even more powerful operations, fancy indexing can be combined with the other indexing schemes we've seen. For example, given the array `X`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0  1  2  3]\n",
      " [ 4  5  6  7]\n",
      " [ 8  9 10 11]]\n"
     ]
    }
   ],
   "source": [
    "print(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can combine fancy and simple indices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([10,  8,  9])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X[2, [2, 0, 1]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also combine fancy indexing with slicing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 6,  4,  5],\n",
       "       [10,  8,  9]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X[1:, [2, 0, 1]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can combine fancy indexing with masking:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  2],\n",
       "       [ 4,  6],\n",
       "       [ 8, 10]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mask = np.array([True, False, True, False])\n",
    "X[row[:, np.newaxis], mask]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of these indexing options combined lead to a very flexible set of operations for efficiently accessing and modifying array values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Selecting Random Points\n",
    "\n",
    "One common use of fancy indexing is the selection of subsets of rows from a matrix.\n",
    "For example, we might have an $N$ by $D$ matrix representing $N$ points in $D$ dimensions, such as the following points drawn from a two-dimensional normal distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100, 2)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mean = [0, 0]\n",
    "cov = [[1, 2],\n",
    "       [2, 5]]\n",
    "X = rng.multivariate_normal(mean, cov, 100)\n",
    "X.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the plotting tools we will discuss in [Introduction to Matplotlib](04.00-Introduction-To-Matplotlib.ipynb), we can visualize these points as a scatter plot (see the following figure):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWwAAAD0CAYAAAC/3RwjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcC0lEQVR4nO3de2xU150H8O9g4xcma6d11XikpHVVwOoGYc8qSosQwQlxlJRdKyHgklqgaFWFooYEagIIWEgTHkIoFc3ySKNmrYQ0Fg11iNiWJTVqiiPYdGJHgAZXiSuiDDQxBS9+4tfdP+iYmfG9d+7j3Me58/38BeM7d86JyW/O/Z3fOSekKIoCIiLyvSleN4CIiIxhwCYikgQDNhGRJBiwiYgkwYBNRCQJBmwiIknkOnXjaDTq1K2JiAItEomovu5YwNb7UJnFYjFUVlZ63QxHsG/yCnL/sq1veoNdpkSIiCTBgE1EJAkGbCIiSTBgExFJggGbiEgSjlaJEBFli5b2OHYf78SlnkGUlxSisXYm6qrCQj+DAZuIyKaW9jg2HDmLwZExAEC8ZxAbjpwFAKFBmykRIiKbdh/vnAjWCYMjY9h9vFPo5zBgExHZdKln0NTrVjFgExHZVF5SaOp1qywH7IMHD2Lp0qV49NFHcfjwYZFtIiKSSmPtTBROzUl5LQRgwawyoZ9jKWCfOXMG7e3t+PWvf43XX38df/vb34Q2iohIJnVVYTwWCSOU9JoC4O1oHC3tcWGfYylgnzp1CjNmzMCqVavw1FNP4b777hPWICIiGZ280I30E81FTzxaKuu7du0aLl26hAMHDuDzzz/HypUr8fvf/x6hUCjlulgsJqSRfjI0NBTIfgHsm8yC3D9Z+qY38ajVfrN9sxSwS0pKUFFRgby8PFRUVCA/Px9Xr17FV77ylZTrgrglYrZt9RgUQe4bEOz+ydK38pLLiKsE7fKSQs32u7K9aiQSwZ/+9CcoioIvvvgCg4ODKCkpsXIrIqJAUJt4LJyag8bamcI+w9IIe8GCBfjwww+xePFiKIqCLVu2ICcnJ/MbiYgCKrGi0cnl6ZaXpq9bt05YI4iIgqCuKix8/5BkXDhDRCQJBmwiIkkwYBMRSYIBm4hIEgzYRESSYMAmIpIEAzYRkSQYsImIJMGATUQkCQZsIiJJMGATEUmCAZuISBIM2EREkmDAJiKSBAM2EZEkGLCJiCTBgE1EJAkGbCIiSVg+IoyISCYt7XFHz1t0AwM2Efme3WDb0h7HhiNnMTgyBgCI9wxiw5GzACBV0GbAJiJfa+3qxcunLxoKtlqBfffxzon3JwyOjGH38U4GbCIiUZo+umYo2OqNoi/1DKreW+t1v+KkIxH5Wnf/qOrr6cFWbxRdXlKoeg+t1/2KAZuIhGlpj2PuzlZ8c/0xzN3Zipb2uO17lk1TTwSkB1u9UXRj7UwUTs1Jeb1wag4aa2fabp+bGLCJSIhESiLeMwgFt1ISdoP28upSQ8FWbxRdVxXGjkfvRrikECEA4ZJC7Hj0bqny14DNgP33v/8d8+fPx6effiqqPUQkKb2UhB01FdMNBdtMo+i6qjDa1tfgrzsfQdv6GumCNWBj0nFkZARbtmxBQUGByPYQkaScnNirqwpnDLCJn8tea63HcsDetWsX6uvr8corr4hsDxFJqrykEHGV4OzmxJ6RwC4zSwH7yJEjuP322zFv3jzdgB2LxSw3zK+GhoYC2S+AffO71q5eNH10Dd39oyiblovl1aWoqZgOwB/9W3Z3MfZ+MIQbY8rEa/k5ISy7u9hW2/zQN6eY7VtIURQl82WpnnjiCYRCIYRCIcRiMXzjG9/A/v37UVZWNnFNNBpFJBIxe2vfi8ViqKys9LoZjmDf/Cu9xhi4mZ9N5HL90j8nln/7pW9OUOubXuy0NMI+dOjQxJ8bGhqwdevWlGBNRNqsBDVZVuoFPSXhNa50JHKR1T0tgrJSj+yxXYf9+uuv41vf+paIthAFntXSt6Cs1CN7uHCGyEVWR8pBWalH9jBgE7nI6kg5KCv1yB7msIlc1Fg7U7XaI3mkrDUpyQk9YsAmclGm1XhB2WifnMGATeQyvZGyLOV75A3msIl8hOV7pIcBm8hHWL5HehiwiXyE5XukhzlsIh8RuUWoE/t6kLcYsIl8RkT5HqtNgokBm8hlbox8WW0STAzYRC5ya+TLapNg4qQjkYucOvcwHatNgokBm8hFWiPceM+g7dPFk4msNmlpj2PuzlZ8c/0xzN3ZKrSdZA5TIkQu0jr3EIDQ1IioahNOXvoLAzaRi9Q2f0oYHBnD1qPnhU1Iiqg24eSlvzBgU9bxsj458TnPNHeo/rxncAQ9gyMA/DGa5eSlvzCHTVkl8Ygf7xmEgltB0c28bF1VGGGDk39OTEiawclLf2HApqziVpVGJmqTglq8HM0umFWGUNprXCrvHaZEKKv45RE/keLY9u55XBsY0b3Wq9FsS3scb0fjUJJeCwF4LMKDFLzCETZlFT894tdVhVGUpz9m8nI0q/Y0ogA4eaHbk/YQAzZlGS92w9OrY9Yb2Xt9bqNfnkboFgZsyipuH2abaZJTb2Q/MDyKZ5s7HF+sovWF4qenEbqJOWzKOm4eZpupjnnBrDK8cfoz1fcmcttOlvfpLYwxcmAwuYsjbCIH6aUVEpN6RjhVyZLpC8XNpxHKzNIIe2RkBBs3bkQ8Hsfw8DBWrlyJ+++/X3TbiKSntRS9vKRQNVjqSew34uaufm4+jVBmlkbYR48eRUlJCd588028+uqr+NnPfia6XUS+ZHYjJL1JTiuTd6IX+TBPLRdLAfuhhx7C6tWrAQCKoiAnx9gCAMo+QdrpzcoqSb20gpWgKDo1wjMk5RJSFEXJfJm6vr4+rFy5EkuWLMGiRYtSfhaNRlFUVGS7gX4zNDSEgoICr5vhCNF9a+3qxd4PruDG2K1/Yvk5ITz9va+ipmK6sM8xQkTflv/mM3zZPzrp9a9Ny0XT4jtN30/tv48RIQD/vbwi5TU7/Wvt6kXTR9fQ3T+Ksmm5WF5d6vrvR0+2/T83MDCASCSier3lKpHLly9j1apVWLZs2aRgnVBZWWn19r4Vi8UC2S9AfN/+/Z3WScHoxpiCN8/2YdUj9wj7HCNE9K27v0vj9VFL966sBMLlqRtR6VWNJJSXFE76vOT+md3cqrISWPWI6ea7Jtv+n4tGo5rXWwrYV65cwZNPPoktW7bgu9/9rpVbUBYI2sILvQlEs9KD6ktL50wE1ZMXujX3zM6UruD+1cFmKYd94MABXL9+Hfv27UNDQwMaGhowNDQkum0kuaBNaInK92bKhWttDGVkHw+/bG5FzrA0wt60aRM2bdokui0UMEFbeCHqFJdMtc9aG0MpAN6OxvEvd92u+ZlBe6qhVFzpSI4RFeD8RERdspGgWlcVxu7jnZN28st02ovItA35DwM2OYoLLyYzGlStjJaD9lRDqbg0nchlRnPhVuYAuJw82DjCJt/z8gxGJxhNFVkdLfOpJrgYsMnXglqmZiSoBnEOgOxhwCZfy1RRIZJfRvJ+aQf5DwM2+ZpbZWp+Gcn7pR3kT5x0JF9za/GN6AUnVje94sIX0sOATb7m1m5yIkfyVnb1y/R58Z5B6Xc7JPsYsMnX3CpTEzmStzNK1vs8M4Gfgok5bPKUkQk2N8rURC44sTNaV2tHMqcmXEkODNjkmUwTbG5WSySX0MV7BpETCqWMis18rp3l4entUMN9QbIXUyLkGb3UgVoe+NnmDmxqOetYe+qqwhM587F/nOthNg3R0h5H/43JhxyYGa3XVYXRtr4G4YDtdkj2MWCTZ/RSB2rBXAFw6PRnlnK4rV29hqo2tr173nL+OfEl0zOYumFTadFUS3l3Ht9F6RiwyTN6E31awVwBTJe4tbTHsfeDKymj9WeaOzBn2/+kBO6W9vik3fESjKQhtE5BL8rLtZTK4b4glI45bPKM3kSfyBzu7uOdqucm9gyOpOTM9b4IjKQhnFjkw31BKBlH2OQZvRFkY+1MhDTeZzaHqxcwk9MdmbYtzURkaWCQTpsncTjCJk9pjSDrqsL488WrOHT6MySPja3kcLWqNhISgVrrupLCqYZGuaJKA7k8nbRwhE2+9ULd3Xhp6RzbOdzG2pnIz9Ear98aAWtN8m391+8AyDzqTX5iAJBSGmhmhMzl6aSFI2zyNRE53LqqMOKX4ng12jNpUjF5BKy3nanRUW/iz3ZGyDyXkbQwYJP0jCywqamYjlWP3JPxWq0vCDPbvNrdEpbnMpIWBmySmtl8r9URe6ZRb/IXweR6FP17pOO5jKSFOWySmlv5Xr0KkPRVmWbvkU6teuaxyM2yQ1aNZDeOsElqbuV7M9WMa23WlH6tUclPAqwaoQSOsElqWqPWKaGQ0FGoXs243peDiBWKrBqhBMsj7PHxcWzduhWdnZ3Iy8vDCy+8gLvuuktk24gy0tqOdExRhI9CtfLfWpOE4ZJCtK2vsf25rBqhBMsj7Pfeew/Dw8Nobm7G2rVrsXPnTpHtIjIkMfLNCU2us3ZrFOr0Jk1uHZNG/mc5YEejUcybNw8AMGfOHJw7d05Yo4jMLM2uqwpjXFGf7nNjFOr0Jk3ctY8SLKdE+vr6UFxcPPH3nJwcjI6OIjf31i1jsZi91vnQ0NBQIPsF+KdvrV292PvBlYkNm+I9g3juNx8jfimOmorpqu8pm5aLL/sn70NdNi0X/3nsf/Ff0au4MtCFsmm5WF5dqnkfq2YWAK/+2x1Jr1xHLHbd9n1bu3rR9NE1DI6MYUoIGFeAr/2jDzMLbn2GX353TmDfbrEcsIuLi9Hf3z/x9/Hx8ZRgDQCVlZVWb+9bsVgskP0CvOtb+mKW/hujk3bXuzGm4M2zfVj1yD2q99j4/dtUqzgevLscL5+OT7z+Zf8oXj59FeFy/++C19Iex8unL060fVy52aeN3//nSW3nv0s5qfUtGo1qXm85JVJdXY33338fANDR0YEZM2ZYvRVlMbWTZdIPAEjQS29opSVOXuiWtsKC1SGUzvIIe+HChWhra0N9fT0URcH27dtFtos8ZPQsRRFnLhqpYU7INMmmVsXxTHOH6rUyVFiwOoTSWQ7YU6ZMwfPPPy+yLeQDRhdpiFrMYTT4TJ0S0p1kU/vyAG7WQatNR8pQYcE9RSgdF85QCqOP4aIe1w0HH+3dUVXTKhuOnMXWo+dVg3UIxg4k8BqrQygdAzalMPoYLupxXS0oqcXmkTFF88tA68tDKxeuwJ0l3XZPjeGZjpSOe4lQCqOP4VYe1/Vy3smva50OI+qMx3BaG83k4s3k90WkjHimIyXjCJtSGH0MN/u4rpW2aGmPo64qjLb1NfjrzkfQtr5GddUiAM3XzeR009uo1y4zfUjHCg9yAkfYlELv1BUr1yWY2dR/TGPVotbrWvuJpAurtFHUwQSJn2d6SmCFB9nBgE2TGH0MN/O4bibnHdbZTEmrHcDNgKkVKEOA6kZMZtqldW1ipJ2c/pC5OoX8iykRcoWZDYysVEck0ipaQb1smvrYxEy7tK5NHLabTMHkyVNWeJBdDNiUwm5lgxYzQdhOdYTW5yyvLp34e3If+2+MYmraiepq7Wppj6P/xuS9Sgqn5mimapR/tJ0VHiQKUyI0obWrN2XvCpEnm5jNeVutjtD6nJkFNzdJSq/eSJT+JTZWUstzp78nobRoKv5j0Xc0UzGi9sMmSmDAzgJGS9ESu8IlGxwZw9aj520vQQfcK1FT+5zErnZaS+ETGyup9U3rPUV5uRPX8tBccgMDdsCZqQfuVtmeFLg5Ck2MRNXeL2JPEbfoVWloVYdkmpg0+/RAZBUDdsCZKVvT2lM6XfL7ZTsgVq/kDlAPzkYWCXGBC7mBk44BZ6ZsbXl16aQJu0z3tbJAxKmJTSPUJiWTiapaIXICA3bAmSlbq6mYPqk6o7Roqu77ze4pYma1oBMSFSglher96r8xOqkt3NOD/IIpkYBTWwWoNzpMf7RXq5BIfr/ZPUUypWjcyIcn+tjSHse2d8/j2sCtTaJ6BkdUUzpMeZAfcIQdcHZHh5nebzZdoDcid3v0XVcVRlHe5DEL9/wgv+IIOwvYHR2mV0EkglnyfY2OivVG5GYmSK1QG73zVBeSCQM2ZZSpEsTMF4JeiuZZB4/z0loU9E+FU1X3zeaeH+RHTIlIxosKi23vnhe2VaheisXMBKlZWouCQiGwAoSkwRG2RLyoeW5pj6dMyiWzOvLVGpGbnSA1Q3NR0MAIXlo6h4teSAoM2BJxOser9ZlaRKcNnFwxqLUoqLykkBUgJA0GbIlkqrBwItDpjaKdSBs4FTyXV5fi5dNXud8HSY0BWyJaFRYlRVMdS5XoLeVOrhYRTfQXUE3FdITLw0x9kNQYsCWileNVFDiWKtE7fsupHLpTuXqmPkh2rBKRiFaFxf+plKUBYsrhkj9TjROLTHiALZE6SyPs3t5eNDY2oq+vDyMjI1i/fj2qqqpEt41UqI0StTbQFzUpmPjMb64/pnpOoehFJlzMQqTO0gj7tddew7333os33ngDO3bswPPPPy+6XWSCW7vJOVkn7cXnEMnGUsBesWIF6uvrAQBjY2PIz88X2igyx63d5Nz6YuB2pkTqMqZEDh8+jKamppTXtm/fjtmzZ6O7uxuNjY3YuHGjYw0kY9yYUHPrZBWe4EKkLqQoGkc+Z9DZ2Yk1a9Zg3bp1mD9//qSfR6NRFBUV2W6g3wwNDaGgoMDrZjjCat9au3rR9NE1dPePomxaLpZXl6KmYroDLbQuyL83INj9y7a+DQwMIBKJqF5vadLxk08+werVq/Hzn/8cs2bN0ryusrLSyu19SaZzC62KxWKmf2ct7fGUTZW+7B/Fy6evIlzurxI6K32TSZD7l219i0ajmtdbCth79uzB8PAwXnzxRQBAcXEx9u/fb+VWUpDt3EI3ebFcnihbWQrYQQ7OahiUtLEEj8g9XDhjAIOSNpbgEbmHAdsABiVtLMEjcg8DtgEMStp4ojiRe7j5kwGsC9bHTZWI3MGAbVAiKHlVYpQNZYVEpI8BWwIsKyQigAFbCiLKCjlCJ5IfA7YE7JYVcoROFAysEpGA3bJCHghAFAwM2BKwW1bIhT9EwcCALQG7tc5c+EMUDMxhS8JOrbPW4b1c+EMkFwbsLMCFP0TBwICdJbgakUh+DNgeYm00EZnBgO0RrdroP1+8ipMXuhnEiWgSBmyPaNVGHzr9GRKHbHKBCxElkyJgBzF1oFUDnX4iMk+2IaIE39dhJ1IH8Z5BKLg16mxpj3vdNFvM1EBzgQsRARIE7KAuq1ZbvRjSuJYLXIgIkCAlEtRl1Wq10QtmleHtaFx1gUsQ00JEZI7vA3Z5SSHiKsE5CKPO9KB98kI3HouEJ1WJADC82x4DO1Fw+T5gB3lZtVpp39vR+KR9QububDW0H7Zb26jyS4HIG77PYQf5kFej+XmjaSE38v1BnQQmkoHvR9hAcJdVq6V6gMmB2GhayI18v4jTb4jIGlsj7E8//RSRSAQ3btwQ1Z6s0dIeN1wVYnQ/bDe2UQ3qJDCRDCwH7L6+PuzatQt5eXki25M1dh/vnLRIBrhZ2pceiI2mhewedGAE99Ym8o6llIiiKNi8eTPWrFmDH//4x6LblBX0VjqqpRaMpIXc2EY1yJPARH6XMWAfPnwYTU1NKa+Vl5fj4YcfxqxZsxxrWNBp5aXDNkeqTuf7ubc2kXdCiqKoPZnrWrhwIb7+9a8DADo6OjB79mwcOnQo5ZpoNIqioiIxrfSRoaEhFBQU2L5Pa1cv9n5wBTfGbv3nz88J4envfRU1FdNt398KUX3zoyD3DQh2/7KtbwMDA4hEIqrXW0qJnDhxYuLPNTU1+NWvfqV6XWVlpZXb+1osFhPSr8pKIFzur3pmUX3zoyD3DQh2/7Ktb9FoVPN6Kcr6grpQI6jlikTkDNsBu7W1VUQ7NLm1eo+IyO98v9IxqLv1ERGZ5fuAzYUaREQ3+T6HHdTd+oKalyci5/h+hO3G6j23cQMlIrLC9wE7iLv1MS9PRFb4PiUCBK/8jXl5IrLC9yPsIOIGSkRkBQO2B4KYlyci50mREgkabqBERFYwYHskaHl5InIeUyJERJJgwCYikgQDNhGRJBiwiYgkwYBNRCSJrK4S4QZMRCQTqQO2nYDLgxGISDbSpkTs7njHDZiISDbSBmy7AZcbMBGRbHyVEjGT4rAbcIN6MAIRBZdvRthmUxx2d7zjBkxEJBvfBGyzKQ67ATeIByMQUbD5JiViNsUhYsc7bsBERDLxTcC2klNmwCWibOKblAhzykRE+iyNsMfGxrBjxw6cO3cOw8PD+MlPfoIFCxbYagg39Sci0mcpYL/zzjsYHR3FW2+9hS+++AK/+93vhDSGKQ4iIm2WAvapU6fw7W9/Gz/60Y+gKAo2b94sul1ERJQmY8A+fPgwmpqaUl4rLS1Ffn4+Dh48iA8//BAbNmzAoUOHHGskEREBIUVRFLNvevbZZ/HQQw+htrYWADB37ly0tbWlXBONRlFUVCSmlT4yNDSEgoICr5vhCPZNXkHuX7b1bWBgAJFIRPV6SymRSCSCP/7xj6itrcWFCxdwxx13qF5XWVlp5fa+FovFAtkvgH2TWZD7l219i0ajmtdbKutbsmQJFEXBkiVLsHnzZmzbts3KbYiIyARLKREj9L4liIhIm1ZKxLGATUREYvlmpSMREeljwCYikgQDtkm9vb146qmn8MMf/hBLly5Fe3u7100S7sSJE1i7dq3XzRBifHwcW7ZswdKlS9HQ0ICLFy963SThPv74YzQ0NHjdDKFGRkbQ2NiIZcuWYfHixfjDH/7gdZOEGhsbw4YNG1BfX48f/OAH+Mtf/mLofb7ZrU8Wr732Gu69916sWLECXV1dWLt2LX7729963SxhXnjhBZw6dSowZVTvvfcehoeH0dzcjI6ODuzcuRP79+/3ulnC/PKXv8TRo0dRWBisk5KOHj2KkpIS7N69Gz09Pairq8P999/vdbOEOXnyJADgrbfewpkzZ/DSSy8Z+nfJgG3SihUrkJeXB+Dmt2R+fr7HLRKruroaDzzwAJqbm71uihDRaBTz5s0DAMyZMwfnzp3zuEVi3XnnnfjFL36BdevWed0UoZIX5imKgpycnAzvkMsDDzyA++67DwBw6dIl3HbbbYbex4CtQ21Z/vbt2zF79mx0d3ejsbERGzdu9Kh19mj17eGHH8aZM2c8apV4fX19KC4unvh7Tk4ORkdHkZsbjH/6tbW1+Pzzz71uhnDTpk0DcPP39/TTT+OZZ57xtkEOyM3NxXPPPYcTJ05g7969xt7jcJuk9vjjj+Pxxx+f9HpnZyfWrFmDdevW4Z577vGgZfZp9S1oiouL0d/fP/H38fHxwATroLt8+TJWrVqFZcuWYdGiRV43xxG7du3CT3/6UyxZsgTHjh3LuJ0HJx1N+uSTT7B69Wrs2bMH8+fP97o5lEF1dTXef/99AEBHRwdmzJjhcYvIiCtXruDJJ59EY2MjFi9e7HVzhGtpacHBgwcBAIWFhQiFQpgyJXM45lDDpD179mB4eBgvvvgigJsjuCBNYgXNwoUL0dbWhvr6eiiKgu3bt3vdJDLgwIEDuH79Ovbt24d9+/YBuDnBGpRNoB588EFs2LABTzzxBEZHR7Fx40ZDfeNKRyIiSTAlQkQkCQZsIiJJMGATEUmCAZuISBIM2EREkmDAJiKSBAM2EZEkGLCJiCTx//EHjvtouHWEAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('seaborn-whitegrid')\n",
    "\n",
    "plt.scatter(X[:, 0], X[:, 1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use fancy indexing to select 20 random points. We'll do this by first choosing 20 random indices with no repeats, and using these indices to select a portion of the original array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([82, 84, 10, 55, 14, 33,  4, 16, 34, 92, 99, 64,  8, 76, 68, 18, 59,\n",
       "       80, 87, 90])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "indices = np.random.choice(X.shape[0], 20, replace=False)\n",
    "indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(20, 2)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "selection = X[indices]  # fancy indexing here\n",
    "selection.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now to see which points were selected, let's overplot large circles at the locations of the selected points (see the following figure):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWwAAAD0CAYAAAC/3RwjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABGDElEQVR4nO3deUBUVfvA8e8MMCzDMrIJiuKGSlqalktqLrlluVZuhbnnkrmlvqlZLpVmvf7KUsvMUsmFUnNJy7XU1NfItYjcRUUQAYGBgWFmfn8YEyM7DMvg8/mnunPn3nNGe+bMc895jsJkMpkQQghR4SnLuwFCCCEKRwK2EELYCAnYQghhIyRgCyGEjZCALYQQNkICthBC2Aj70rpweHh4aV1aCCEqtebNm+d6vNQCdn43tWUREREEBweXdzNKhfTNdlXm/j1ofctvsCspESGEsBESsIUQwkZIwBZCCBshAVsIIWyEBGwhhLARErCFEKKEkpOTWfjhRzRr/SS1GjSiYZPm9B/0EkeOHMGaBVFLdVqfEEJUZkajkblz5/Lxx0up17Ql3foPJaB6dZJTdZz5/QQhQ17G3c2Vzz//nBYtWpT4fhKwhRCiGIxGIyEhIVy/fp2Pw/bg6lkVNycH8+sNHmnOSyPHkhxxmGeeeYZNmzbRsWPHEt1TUiJCCFEMCxYsICoqih9//BE7N2/UjpbjX7WjPYlpmQwYMIBNmzYxcOBArl27VqJ7SsAWQogi0mq1fPTRR6xZswYnJyc81Sq06ZmW56Rn4qlWAdCxY0deeuklli1bVqL7Fjtgf/bZZwwYMIB+/foRFhZWokYIIYQt2bBhA23atKFWrVoANAnQkKzL5Fq8lt+uxPPjH7c4fjkef3cn83vGjBnDl19+iU6nK/Z9ixWwjx8/zsmTJ1m/fj1r167l1q1bxW6AEELYmq1btxISEmL+b3+NM00CPPg7Jpk72nS8XR1oUNWN09fvEp2YBkBQUBD16tXjyJEjxb5vsR46Hj58mPr16zN+/HhSUlKYPn16sRsghBC25s6dO1SrVs3iWHSSjpa1vSwePCbr9Jy+noi/xhmAatWqER8fX+z7FitgJyQkcPPmTVasWMH169cZO3Ysu3fvRqFQWJwXERFR7IZVVDqdrlL2C6Rvtqwy968i9s1gMBAZGYmnp6f52B8XkqjirCQxWxw0mUxcSjMSaJcIQFxcHLGxseb+FLVvxQrYGo2GOnXqoFKpqFOnDo6OjsTHx+Pl5WVxXmUsifiglXqsLCpz36By968i9u3xxx8nKiqK4cOHm49dNUSTlmHIMcKuprIjONif9PR0/vrrL7p27UpQUBBQRuVVmzdvzqFDhzCZTMTExJCWloZGoynOpYQQwuaMHj2alStXkpn578yQrAePyTo9RpOJZJ2eZF0mTQI0AGzevJlGjRqZg3VxFCtgd+zYkeDgYJ5//nnGjh3LnDlzsLOzK3YjhBDCljzyyCPUrl2bL7/80nzMX+NMl4eq4qyyIy4lHWeVHV0eqoq/xhmdTsfixYsZP358ie5b7JWO8qBRCPEg++yzz+jYsSN+fn706tULuBe0sx4wZklLS2PgwIE0aNCAPn36lOiesnBGCCGK4aGHHmLHjh2MGTOG0aNHc+rUKYvXdToda9asoVWrVri4uPDVV1+hVJYs5EotESGEKKbHH3+cU6dOsXLlSnr16kWVKlUICAggPT2d06dP89hjj7FgwQKeeeaZEgdrkIAthBAl4uvry6xZs5gxYwa//fYbcXFxODo6EhQUZF4JaS0SsIUQwgrs7e1p1apVqd5DcthCCGEjJGALIYSNkIAthBA2QgK2EELYCAnYQghhIyRgCyGEjZCALYQQNkICthBC2AgJ2EIIYSMkYAshhI2QgC2EEDZCArYQQtgICdhCCGEjJGALIYSNkPKqQogSS05OJjQ0lIMHD5KUlISrqyutW7dm6NChVKlSpbybB0B0YhqnrycSr83AU62iSYAmx3ZeFZ2MsIUQxZaens7UqVMJDAxkz549PPvss4wbN46+ffsSHh5OnTp1eOWVV0hOTi7RfaIT09h9Lppvjl9l97loohPTivz+PX/GkJZhwNvVkbQMA3v+jCnydcqbjLCFEMWSmppKjx498Pb25syZMwQEBFi8PmjQIGJjY3njjTdo3749e/fuxdPTs8j3ua3Vc+LPGNyc7PF2dUSbnsmeP2PMO5Jnl9co+vT1RNyc7HFzcgAw//P09USbGmXLCFsIUSzDhg2jRo0abNq0KUewzuLr68sXX3xBhw4deO655zCZTEW+z/m4DHOwVSoUuDk54OZkz+nriRbn5TeKjtdmoHa0HJ+qHe2J12YUuT3lSQK2EKLIzp07x6FDh/jiiy8K3FxWoVDwwQcfcPv2bQ4ePFjkeyXqDIUKttlH0fcHdk+1Cm16psX52vRMPNWqIrenPEnAFkIU2fLlyxk9ejSOjo4Wx/PKNSuVSsaNG8eyZcuKfC+Nk12hgm1+o+gmARqSdZkk6/QYTSaSdXqSdZk0CdAUuT3lSQK2EKLINm7cyLBhwyyOFfRg76WXXmLHjh2kp6cX6V5B3qpCBdv8RtH+Gme6PFQVZ5UdcSnpOKvscs2BV3Qleuh4584d+vXrx5dffkndunWt1SYhKi2TycTBgwc5ceIEqampeHh40LlzZx5++OHyblqhGQwGEhISqFGjhsXxgh7subu74+bmRkJCAn5+foW+n4/agaB6VTl9PZG4lHQ81Spa1fHKEWybBGjY82cMcG9krU3PJFmXSas6XgD4a5xtLkDfr9gjbL1ez5w5c3BycrJme4SolAwGA5988gkPPfQQEyZMICYmBpPJxMWLF+nevTtPPvkk33//fXk3s1CUSiUKhQKDwWBxvDAP9vR6PSpV0fPG/hpnujf2Z3DLQLo39s818FaWUXR+ij3CXrRoEQMHDuTzzz+3ZnuEqHTS09MZMGAA8fHxfPbZZ7Rr1w6FQmF+fcmSJWzbto2pU6fy22+/MW/ePIvXKxqFQkHdunX5/fffadmypfl4Vkoia2QNlrnmixcvYmdnh4eHR6m1rTKMovOjMBVjns3mzZu5desW48aNIyQkhLfffjtHSiQ8PBwXFxerNbSi0Ol0lfZXhfTN+kwmE9OnTycjI4PFixfnO7qMj49n2LBh9O3bl6FDh+Z4/bZWz/m4DBJ1BjROdgR5q/BR3wuOZd2/VatWcfHiRd59912L9h2LSkXtoMTZQUGa3oRWb6RVDRd81A4sXrwYgGnTphXpXg/a38vU1FSaN2+e6/nFCtgvvvgiCoUChUJBREQEtWrVYvny5fj4+JjPCQ8Pz/OmtiwiIoLg4ODybkapkL5Z36+//sqQIUM4e/Yszs4Fj/yuXr1K06ZNuXz5MhqNxnw864Gem5O9RX426yd/WfcvLi6OoKAgfv/9d2rXrm3RztwWrsTGxtKoUSOOHTtW5OddD9rfy/xiZ7FSIqGhoeZ/zxphZw/WQoh7Pv30U8aPH28RrPOraREYGMjTTz/N119/zcSJE83vqWgr9by9vZk/fz49evRg3759VKtWDcg9JREfH0/Pnj0ZM2aMTE4oIZnWJ0QpSUlJYfv27RbpjcLUtBgzZgyrV6+2uFZFXKn36quv8vLLL9O6dWs+//xzUlJSLF7X6XSsXbuWVq1a0bZtW+bNm1dOLa08SlxLZO3atdZohxCVzq1bt/D29raoVleYkXKjRo24evWqxbUKeqBXXv7zn//QqlUrPv74Y9544w06depElSpVSE5OZt++fTRr1owlS5bwzDPPlGs7Kwsp/iREKTEajTmWbcdrM/B2tVwdqHa0Jy7l38UkSqUSo9FocU5Bc4zLU4cOHejQoQNRUVEcOnSIpKQk3NzcmD9/PvXq1Svv5lUqErCFKCVVq1YlNjYWrVaLWq0GCjdSvnDhgjknnCVrjnFBi0fKU40aNRg8eHB5N6NSk4AtRCnx8PCgY8eOfPPNN4waNQoo3Ej5o0+X06JLb745ftXioWRln2MsCiYPHYUoRePHj+eTTz4xrwosaDXeuYtRbNn8He2f7W/ThfZF6ZARthClqHPnzvj6+jJhwgQ+/fRTFApFniPl1NRUBg3oT7d+L1EzoDpQ/tP3RMUiI2whSpFSqeTbb7/l5MmTDBw4kAsXLuQ4x2Qycfz4cTp27IinXwCjp86yeL28p++JikNG2EKUMg8PD/bv38/cuXNp3bo1jz32GN26dUOtVhMfH8+mTZtITExk8uTJ1G3fjzS9ETc7O/P7K8L0PVExyAhbiDLg7OzMwoULuXbtGoMGDeLSpUscO3aMmJgY3nnnHc6fP8+rr75K0xpVKkWhfVE6ZIQtRBlydnZmyJAhDBkyJNfXrTl9L78l8MI2ScAWooKxxvS97MWiCtppXNgOCdhClLGyGPlWtGJRwjokhy1EGSpM8SdrqIjFokTJScAWogxlH/kqFQrcnBxwc7Ln9PVEq94nvw1phe2SlIgQZSh78ad4bTpX4lJJ0ukxYbJqasSaxaLk4WXFISNsIcpQ1sg3XpvOqai7pGcaUdkrUNkprZoasdaGtGWVwhGFIyNsIcpQ1sj3clwKTvZKUJjQZZhoWkODg52Cg5ExVFE7WmU0a43ZJvLwsmKREbZ44EQnprH7XDTfHL/K7nPRZTpazBr5pmcayTAYcLRX0rSGBk+1ivRMA0cvxVeo0aw8vKxYJGCLB0pF+Invr3Hmyfo+NA/0pFlNT/ODwMhbyXipVaX+QLIo5OFlxSIBWzxQymqWRkGaBGjMS9DjUtI5cuE24dcSyDQYLUav5T2a9Xd34vjleH784xa/XYnnWrxWlsqXIwnY4oFSUX7iZ6VGdHoDRy7EAfBIdQ+MwKmoRHN7ynM0e292yF0aVHXD29WBO9p0/o5JpkmAh+Svy4k8dBQPlIq0ma2/xpkqahWdGvri5uRgnjmiUJi4HJeMg517ue7bmP3XSA1PFwCSdXqik3Q0KZcWCQnY4oFSkvnJly9fZu/evSQmJuLi4kLz5s1p2bIlCoUi3/flN485a1521pxsnT4Trc7A7WQdTWpUKdd9GwuzYbAoW5ISEQ+U4sxPPnDgAM8++yyPP/44hw8f5tatW5w9e5aQkBCaNWvGqlWrcuxynqWgh5yeahXXE1LNc7L9PZzxdlfhaG9HgjaDA5GxpT6TJa9ZM/LAseKREbZ44BRlfvJ///tflixZwrx589i0aRMuLi7m14xGI3v37uWtt95i9+7dhIaGolJZBrOC5jH7uzux9uhV9AYjnmoV9nYKYu7qSM0wkqK7RbfGfuYgXxqV9vKr6mfN1ZLCOmSELUQeVq1axbJlyzh69CjDhg2zCNZwb/uvrl27cvDgQTIyMhgxYgQmk8ninPwecmY91NM4O+Dp4sDtZB0R0cmYAD93R9INRs5cT0JvMBZ6JktMTAzr169nxYoVrFmzhnPnzuV7fn6zZqy1WlJYT7FG2Hq9npkzZ3Ljxg0yMjIYO3YsTz31lLXbJkS50Wq1TJ8+ncOHDxMQEJDvuY6Ojqxfv54mTZrw66+/0qZNG/Nr+T3kzAqWgV5q0jONKJQKXB0NxKWko3FxQOPsgMFgZF9ELD5ujvnWGzlz5gzz5s1j9+7ddOrUCS8vL9LS0pg5cya1a9dm/PjxDBgwIEe+vaA8tTVWSwrrKVbA3rZtGxqNhsWLF5OYmEifPn0kYItK5ZtvvqFdu3YEBwdbHM/rAaKLiwvjx49n2bJlFgE7v7TCgchYvF0dqeXtwqmou9xNy8TD0Y5Mo5H4lAxU9kr+jklBqQQPZzvUjg65pkZWrlzJG2+8wZtvvsny5cvRaDTm1/R6Pdu3b2f+/Pn88MMPrFq1CgeHf788KtKsGVGwYqVEunfvzsSJE4F7Oz7bZdswVIjsynMZeEmsWrWKMWPGWBwr6AHiyy+/zM6dO7l79675PfmlFbKCpafakaY1PHBztOdOqh6NswOZJhNJaXpUdgrsFQou3U7F38M5R2pkw4YNzJ8/n3Xr1jFx4kSLYA3g4OBAv379OHLkCDExMYwbN84ibZN9AY/sIVnxKUz3J92KICUlhbFjx9K/f3969uxp8Vp4eHiOnF9loNPpcHJyKu9mlApr9+22Vs+xqFTUDkqcHRSk6U1o9UZa1XDBR+1Q8AWsqKh9a9euHZs3b8bHx8d87NerWnSZRlxU/45zUjOMONkreSJQDUCPHj345JNPqFOnToH3uP/ziUnW8+ftDFz++d0beScdgxH83BzwdFbi4WRPQx9HEtKM9GjgTnp6Op06dWLlypXUqVOnwP5ptVr69OnD4sWLadq0qUU7zsdlkKgzoHGyI8hbVeZ/Pvl50P6fS01NpXnz5rmeX+xZItHR0YwfP57BgwfnCNZZ7v85WRlERERUyn6B9ft29Vw0QYEGi5/byTo9qSo7goP9rXafwihq3xQKBQ0bNsTX19d87GTSVWq4OqLMlgc2mkzEpaQTHBwIgIuLC4GBgYW6VzAQVO/fFMuj1VX0aOHExt+iUCoApwyquKjwdnXEZDKRpNNTxacK1f75/NasWUPLli157rnnLPqX37zvyZMns2vXLgYNGmTRjicL/cmUvQft/7nw8PA8zy9WwI6Li2P48OHMmTOH1q1bF+cS4gFgywsvqlatyqVLlywCdkH53vT0dG7evGnxnrzcH1Q7NvA1B9XoJB1pGQb0BhOnohJJ02diMoK9UmExre6LL77g9ddfz3Hd/DbfHTZsGHXq1OHu3bt4eHiU+HMSZatYOewVK1aQlJTEsmXLCAkJISQkBJ1OZ+22CRtnywsvBg0axKpVqyyOFZTv/e6773jssccKDNgF5cKz7uNgp6CWlzOXbmv57Vo8GZkGizoeFy5cyPHTuaDiVp6envj7+3Pjxg0rfEqirBVrhD179mxmz55t7baISsaWF16MGDGChg0bsnDhQry87rU36wHi6euJxKWk46lWmZeOm0wmli5dyvTp0wu8doGLaf65z8HIWM7eSKKWlwvdGlXF0d6O09fv4uvuhL/GmczMTOztLf8XLsyvGnt7e/R6fYk+H1E+ZKWjKDX5BbiKzs/Pj1GjRvHCCy+wc+dOnJ3vtTm3eckmk4lZs2ZhMpnyfJ6TXUFBNS0tDVNqAvaZWjrU98bDxfLcrMDu6+vL1atXqVq1qvm1gtI2er2+0GkbUfHISkdRqvw1znRv7M/gloF0b+xvE8E6y8KFC6lWrRqdO3fmzJkzuZ4THR3NK6+8wo4dO9i2bVuOEW9ucksVJWl1RB7fT7du3fD09KR58+YMfbo1o3u0YM2ni4mLiQYsS8H279+f1atXW1ynoLTNtm3baNSoEf7+ZfvQV1iHBGwh8mBnZ8eaNWvo168fPXr0oG3btixdupT169fzxRdfMGDAABo1aoS9vT2HDh0q9Kj1/qB66fIVJg/uzs51KwgJCSEhIYHo6Gi++/Uv3vxkHYnxdxjdpyPb1q+2GC2PHDmSDRs2kJCQYL52fvO+TSYTn3zyCePGjSuVz0uUPkmJiAovv2lqpU2pVDJ16lQmTpzI9u3b+fHHH83lVZ988klWrlyJu7t7ka6ZPVV0/vI15r3yPGPHvcrbs2ZYnNckQENsUj2GTZvPcy+P4Y1XBhGfnMbit/8DQLVq1QgJCWHw4MEsXLjQ4vq5fT6LFi0iPj6efv36FeOTEBWBBGxRoRU0Ta2s2Nvb07dvX/r27WuV62UF1e6vD2P8K6OYc1+wzjonK7CnVfHjvZUbmT6kJzFD++L/z8KX//73v/Tv359XXnmFtWvXEhQUlOM6CQkJvPPOO3z//fccOHAgR0VBYTskYItyo9fr2bZtG+fOnUOn0+Hp6UnPnj1p2LCh+ZyCZlRYU1mP5CMjIzl58iTff/99IdsRyPUJE/j0009ZuXIlcO+LJCwsjIkTJ/LEE0/QrFkzBgwYgI+PD1qtlj179vDdd9/x7LPPcvToUby9vUutP6L0ScAWZU6r1fL++++zcuVK6tevz5NPPombmxvXrl2jQ4cONGrUiClTppCSksL6XYewM+lxdfOgZfvO1G/UpFQW35THSH7FihWMGDECR8d/Z4EU1I6s6YaLFy821w2xs7Nj/PjxfPDBB4SFhbF7927u3r2Li4sLjz76KJGRkRYzSYTtkoAtytSdO3d4+umnqV27Nnv27KFRo0YWr8+fP59hw4bRq1cv6tatS+tufVHYq0hLSmDepBFU8fKhZ8gYOj1d8PS5orD2SL4wo/WjR4/y4YcfFqkdfn5+NGrUiLNnz9KuXTuL9zo5OZkXsonKSWaJiDKTkZFB7969adeuHRs2bMgRrOPj4+nWrRsqlYqffvqJjIwMWjzckE4vjGDg+P+wetcx+g1/ja8/epc96z6xatusuZt6QSsZs6SkpODm5pZvO+K1Gfwdk8Tuc9Hmaodubm4kJycXuV3C9knAFmUmLCwMe3t7Fi9enKOQfmZmJn379qV169Zs2LCBp556is2bN/Pu3Nl0CPLEWWVHQlom7Tp3Y++Bn9m1bQuffvqp1dpmzWX0BS0Pz+Lu7k5iouWx7O2I12ZwKiqRJF0mfu7O5sAfG3enyDNTROUgAVuUmWXLljFp0iSUyn//2mXVy566eCWxiVqmzXnHHMybNWtGvXr1OH7wR4vFNw8HBbJ161beeustUlNTrdI2a9aFLuxovWPHjmzevDnPdlyOS0ahMGEyKajjo8bNyYG0hFtcuHCRRx99tMjtErZPArYoE+fPn+fy5cs8++yz5mPZUwc/bw2lx8Bh7Pvrtjl1EJ2YRqse/Xn3/5bn2Pygfv36tG7dmvXr11ulfdkXnFyITeavW0lo0zM5fT2xyJsuFHa0Pnr0aNauXYtWq821HdF3dbg5OtC0hgee6nsPJg9+v54nuvVGrVYXs6fClknAFmUiKiqKBg0aWCzdzkodpCXEcuV8BJ2f7mVOHWQFc2ffQG7dvMG+iFg+2nee09f+XdU3evRovv76a6u10V/jTJMADa6ODjT0c6eur2ue+ee8RCemkaBN50BkLEcu3CYuJT3P0XpgYCCdOnVi5syZOdrRvbE/3Rv708DPzRysr1z4ix++DeX5l4Zbpb/C9kjAFmXCZDLlugGs2tGeuNhb+FWviYNKZU4dnL6eSKbRSFSCDqPJhK+bI3ZKWH/imjl4BgcHF7pM6G2tvlBblR2MjOVyXAq/X0vgVFQCeoOp0DuWZ33JODnY80TdexX+jlyIQ6c35Dk9cOXKlezfv59p06ZhNBotXsueHvn7zzPMGDWQQRNm80y7xwrVZ1H5SMAWZaJ69eqcP38eg8FgPpZf6iBem0FsUjopMVFovKuiUCjQOKswGE2FCp7ZRSemcSwqlbQMA0qFguOX7vDergjWH79iEbijE9P49eIdFApwd3IgPdPIqahE0jONhZotkv1ho7erE23q+dCpoS9V1Ko8pwZqNBp+/vlnTp48SVBQEIsXLyY6Ohqj0UgVJwWOsX/wwfTRzBjRnzEz5vHW5FdsqoCWsC4J2KJMNGzYED8/P3788UfzsawRpLOHN7duXCM+WWtOHXiqVdzRZnByTxjNn+oNgC7TgLeryhw8IyIiCAgIKPDep68nonZQojeYOHP9LkqlAm+1I3/HpFikO05fT8TbVYUCJQqFAmcHe5xVSiJvJRVqtkhxpwZ6enqyZ88e1q9fzx9//EFwcDAODg54eHiwcO5shrzQm+gbUbwzZZQE6wecLJwRZWbcuHF89NFHPP300ygUimy1MhypXqcBv//yIxNHDTUHpbC9R7lxKZJhbbqSps8kLcNIgK+rOXh+/vnnDBkypMD7xmszcHZQcOVOCs4qJc4O9uY9ErPSHf4aZ+K1GTTwc+PM9SQAnByUmIxwR5tRqNkiBdWizo9CoaBFixa0aNECgOt3Ujh7M4mEVD2eahXJmXa4FngVUdnJCFuUmUGDBhEXF8e8efMwmUzAvw/Y5r0xlZ83r6Gq+70HbMr0JPYtnUG7geOITzeislMQ5OuKnVJBkwANf//9N0ePHrXYTDYvnmoVaXoTybpMnOztANDpjbg5OliMgD3VKhzt7WhawwNHeyVJOj0mTLSu41moka21pgZGJ6axPzIOnd6It6sjF8+fZ+a7S3hz3rssX76c33//vUjXE5WHBGxRZpycnNixYwdhYWGMHj2aq1evml/r3bs3jo6OTJo0iR07dtC6dWtGDB3C6oUzeSq4KjU81fhrnOjyUFXsMpLp3bs38+bNw8XFpcD7NgnQoNUbsVcq/hmpG0jTG6jl7WIxAv53L0UlTWtqaFazCrW9XenQ4F4djqw543k9uLTW1MCsXPif//uZN0YNYPao57gWeZo/L17j1KlT9O3bl1atWhEaGmr+4hMPBkmJiDLl7+/P4cOHmTt3Ls2aNaNNmza0a9cOJycnmjVrxqeffsq6det4++23ee211+6955/RrcFg4IcffmDy5MmEhIRYFOK/evUqmzdvJjY2FgcHB+rVq8dzzz2HWq3GX+NMqxou3DC48evFO3i7qngkwB0HO6XFHpP5bWlW2OJQWf8em5SOv4ezeS/LohSSupOSzg9ff8z+Hd8xdMIM2nV9FnsHFXEp6QxuGWj+HObMmcPevXv54osvrPbnIyo2CdiizGk0GpYsWcKCBQvYtGkTf/zxB2lpaXh6erJ3715++ukn3nvvPbZu3UrHjh1xdnYmNjaWjRs3Uq1aNd5//31zEf7jx4/z5tvzOHb0KI937E5AjZr4qI2EhYWZA/usWbPwUTvwZHAgHRr4mosyOavscuwxmVfx/6IUhyppIamD363m8J4f+PibnVTx8gEgWac3/xKws7OjZ8+edOrUiV69ejFp0iTZReYBIQFblBu1Ws2wYcNyHG/Xrh1vvvkmW7du5cyZMyQkJFClShW2bNlCs2bNzOdt3LiR8a9OoOfwSXw2ewleGnfzzuxTZ1Ql424sH3zwAa1bt2bZsmUEBwfnGZALUtDGudmr8/1x8y4PV/fADYdcz81PXFwcmz7/iLdX78RercFoMuW527xarWbz5s0EBwfTpUsXgoODi9wvYVskYIsKSaVS0b9/f/r375/r6/v37+e1115j/mcb8KtVP9fRbPfGgSxdupRPPvmEUaNGcfz4cfz8/IrVnvxmgNyfLlHZKTlxJYGWtT3NqxQLO1tk9erV9O3bh4GdmlmkZup4qzl9PZEDkbEW5Vo9PDwYPXo0GzZsoFevXsXqm7Ad8tBR2ByTycS0adNYsWIFHtXr5jn3OTIykkmTJvHWW28RFRVFQEAAPj4+vP7661y8eLFI98xvBsj91fmC/T0wmSAiOqnIs0W+/PJLXnnlFa5FnmH94v8wuVcLejevxRONazFjRH8iju4hWauzmD8+atQoduzYQUZG0UvBCtsiAVvYnBMnTpCQkEDv3r1zlCP9/Vo8u89E8dHc6bR78klcXFwIDw8nLCwMf39/Dh06hFKppFWrVkycONFi5WV+8tuN/P4FM55qFS1qe5KeacxxbkEuXbrE66+/zuDBg3n44Yf57bff2HQkgmVbDvHsCy/x/bpVjO/TlhuRp8wrPqtXr45KpeLOnTtF/zCFTSl2SsRoNPL2228TGRmJSqViwYIFBAYGWrNtQuTq66+/ZtSoUSiVSpoEaNjzZwyJqXr+jklGgZFdS2fjYEpnQeh+ej5WB3+NM2lpaVSrVo2oqCjef/99Zs2aRb9+/Rg+fDhfffVVjjonuckr/51busTRXsmT9X3o3ti/0P26desWGRkZdO7cmbfeestchvbQzav4+/lSvUcfOvbow/8O7WPR6yN5df4ndG88ELj3IDIzMzO/y4tKoNgj7L1795KRkcHGjRuZOnUqCxcutGa7RAWn0+k4fPgwO3bsYP/+/cTExJTZvaOionjooYeAf0e+t5LSyDSaOPvTRkhNYOHyNfh4eljUHQkODiYqKgoADw8Ptm/fTkRERImnxVlrwczgwYNxc3Nj0KBBFjXD76+50qLdU0x551M+fXMCd+7cISkpieTkZLy8vHK7rKhEij3CDg8PN+8p17RpU86dO2e1RomK68qVKyxfvpzVq1cTGBiIj48PqampnD59mm7dujF+/Pgcew0WR357It5f+c9f40ygl5pHaziwbso3zFq8ApXKEXuTyWJmxv2jaBcXFxYtWsSECRMYOXJkoUbZuclv/nZhnTx5kgsXLjBq1ChWrVrF4sWLza9l/YoAzPO6az3Sks5duvLVV1/h7OxMu3btCrWISNi2YgfslJQUXF3/rW6Q9ZMse73jiIiIkrWuAtLpdJWyX1Bw337++WdmzpxJ7969WbNmjUUKLCkpiW3btjFw4EB69OjB5MmTix0Ab2v1HItKRe2gxNlBwcXbJs6cv0qrGi74qB1wcXHhl19+ISgoyPyetAQt23YdxtFZjUsVb6KuR5GaYcTJXskvv93lj+gUDh4LR1OvGb/8dgYf9b30RdWqVUlLS2Pt2rU8/vjjxWpvlkA7CHQHSCUxOpHE6MK/97333qNv37507dqVQYMGMWDAANRqNbe1es7HZRB1Nx1thgm1SkEND0eCvFU816cn06dPx8HBgRkzZjywfy9tWVH7VuyA7erqarFThtFotAjWQKWcFxoREVEp+wX59+3AgQPMmTOHH374gVatWuV6TsuWLZk8eTI9evQgNDSUd955p1D3vX80nWBIp7qXA+lJdzCmG/Dz8kHp7E6qyo7gYH9ee+01hg0bxvvvv29OHWj809gS+iWtn+pBQEANtOmZ2OkyaRLgwenrd7l+7Q+S42/zROdnuZKhIKjevw8B+/fvz7Vr1wpVSKq0/Pbbb+zatYsGDRrw3HPPMX/+fJZ+sZYTCfF4+tpTo4a9eT521gNMY7PGjBw5kiZNmtC2bdsH8u+lrcutb+Hh4XmeX+wcdrNmzfjll18AOHXqFPXr1y/upUQFl5mZycsvv0xoaGiewTqLt7c3O3fu5KuvvuLUqVMFXjv7NmFeahW//+8Yb08Zy6geLXlrwlDmTx7Fy0+35p0JL/HTrh1kZmbyxBNP4OzszO7du83X8dc44+NkRKPRWMzMiE7S4eZkz+EdG3nmhZeo4uaSY0MCDw8PkpKSivvxWMXdu3fx9PQEMG8u3KfXs9yNvpLrZr7Xrl0jJCQEo9HIwoULi/1rRtiWYo+wu3TpwpEjRxg4cCAmk4l3333Xmu0S5ej+EW/0mUMEBATQpUuXfM/LyjP7+PgwduxYli9fzmeffZbvvbLmMDsqjSx+YxIRZ37nka796Tv+TTo3rQtARrqOn3Z+zw9rl3N861d8//33LFq0iGHDhnHw4EEaNGgAgJ93FdRqGNzy31TN1pM32Pf9ek7+dpwOQ6cTr81A4+Jgkdu+e/cubm5u1vr4ikWtVpOSkoKPjw8qlYrvvvuOQeP/w9tjBlCzThBtuzyDm0cVUrXJHD6wh8vnwhkyZAh+fn74+fnJDJEHRLFH2Eqlknnz5rFhwwY2btxI3bp1rdkuUU6yj3i9XR1JyzDwf0uXM/jlkQWel30xx8iRI9m0aRMpKSn53i9em4GTvYL3po0jLVXL51v28/KoccRl2HPkwm0ORsZyIiqZuq26s+/gIdq2bUvnzp154oknWLhwIe3btyc0NJQrMYk413iIr9Z/a66kd+bvK3z90TvsCV1O/+kfYK/24FRUIjcS0syrDk0mE9u3b6dNmzal84EW0uOPP27xi8He3p4RE6ayfMdRerwQwqXIPzl6YDfnToXTpmNXrl27xpgxY9BqtTKd9gEiS9OFhdwKF928ch63Wo0LPC/ruL/G2Tzyu3btmnkKXm481Sq2bVhLfFwsi1d/i0rlSGx8Ki7mhSgmQAFKBUqlkoULF3L79m1mzpzJJ598QmBgIG++PY/xr02iXbdexN64ytJ3ZnEnLo6zx3+mefvuDHlvLfYqFU72duj0BiJjkuje+N4S9UOHDmE0GunYsaN1P8giGjduHJMnT2bMmDHm9Ma92SHpPNbpGdo/3dsih61WO7NixQpGjBiBSlXwkndROchKR2Eht22uMjPS0Rpy30A3u/u3w3JyckKn0+V7v0eqe7Bj42qeHzUJewcVyTo9kTFJ1PZ2wdnB3rxVl1plx+nriSgUCubNm0doaCjJycl06NCBWZ+EMv/zMLy8vKlZJ4gTB3+k8aOPMWbZTmYv/Ih2TYJwUCru7TDj6EBNTxfzYpoZM2bw2muvlUoOuKD62dl16tQJg8HAunXrzMfyW135119/sXbtWl555RWrt1tUXDLCFsC9WT6//PILYT/s487tWFQO9vgHBNK932BcPapgTIm3OL+g7bCMRiO3bt0yP0iD3HPeVyJO4YCBlm3bm+cwa5wdiE3KwEVlj7uTAzq9kfMxKej095aRBwQE0KlTJ0JDQxkzZsy9rb0aNiQ4OJhBoyfy5oTh/PrzQRoEtOHoxTiC/d0J9nWiRoAvyTo9Krt79bOHDh1KYGAgPfuHsPtcdK5zvu+X3/zw+88rTP3sLAqFgo0bN/LUU0/h5OTECy+8AOS+ujIiIoLu3bvzwQcfSDrkASMjbME333xDUFAQH3/8MU+1e4JOzw3lqb4v4uHpxZvjh5CSnMS+Tass3lPQ6r59+/bh6+trDih55bz3Hz5G965d6PFIdQa3DKR7Y39MgFIJziq7eyNslR1KJdxN05vv37VrV/NWWdlXAt7VZdJx7Luo3DzYN3cQu9YsZW/4XySm6YmIiGDJnCk81zKI2rVrc+zYMb7fto0XBr3IyRPH8VKrcuTisysob5/d/QWhss/wyEvjxo358ccfmTJlCj179mTXrl0YjUbz62fPnmXs2LG0bduW+fPnM3To0ML88YpKREbYD7iFCxfy+eefs27dOjQaDQ899NC/o8jWneg/YgLhuzbw9qwZbNu2zVzCs6DVfcuWLWPcuHHmVENeOe9LN+PwVqst2qRxVpGUlkma/t4ejLpMA0bjveNZ1Gq1eR1A9pWAl25rsVfZ02PcPLwyYvjpu7V8Nbkvmfp0TAYDKpUjgwcPYsaMGTRs2JCwI3/yw3cb+OTtKdSsXY+Zi1fg5qQq8sYEWf/MGnlfup1CPV/LmSeFqYndtGlTIiMjWb9+PbNmzWLQoEF4eXmh1Wqxt7dn9OjRnD17lmrVqhX0RysqIQnYD7CNGzfy2Wef8euvv+Lv729ecXX/z/Bnmk7j15/30b9/f37//XeLOh65/bxft24dJ0+eZO3ateZjeW0AgKOauNuWpU5r+6hxcrDjdorOvLN5gMYFf42T+Zy4uDg8PDzM7cj68riVlIa/hxO1vd3wVPvyyMOLUJoy+evPP/jlx+0EBARY5Kv19moGjRjHgJdH89G86cwcM5j3Vm4gXpuzil9efbgQm0xsUrpF+uNafCpODnbU9Pz3y6iwNbFdXFwYMWIEw4cP586dOyQmJuLi4oKvr2+OxWniwSIpkQeUyWRi7ty5fPnll/j7F1xRbsuWLfj5+dGzZ08SEhJyPSczM5OlS5cybdo0du7caVG64P4CRnAvgLXt0Int27eTlvZvWqFJgAY7pYIgXzfaBfkQ5Otm3i09y6ZNm+jWrZv5v7N2X+/e2J/6Vd3NgfHAzi1EnDnJ6+98TI0aNXI8XMxql72DA5PnfkgVL29Wffx+roE1rz7cTdPnSH80qOrO3zEpJSoIpVAo8Pb2pl69elSrVk2CtZAR9oPq559/RqFQ0KFDB4vjeT1Uc3Z25sCBAzRs2JBatWrx/PPPM3DgQLy9vUlNTWXfvn18/vnn1KlTh8OHD+eYl5+Vtrh++TzH9+4gNjYGg0lBm0cf4uGHHyYsLMy8NLygdMvJkye5fv06zzzzTI5+3V8o6du1K+kzYjIP19Dk2kcFJuK0egI0984f9OobTB/Skw/fW5Djc0nQpnP0UjxeahUN/NxxtL+3ia/GWZVjxkz1Ks7o9AbzDI/iFIQS4n4SsB9QmzZtYujQoRYjzttaPSfymdlQu3ZtunfvztNPP83du3dZsGCB+ef6o48+yq5du3j44Ydzvd/fp//HZ7PnEPHXXzzRrS8NGj2Mv5uKy3//yYkTJzh27Bh169Y1L2DJK92SlpbG+PHjmTRpUq4jzuzB/sRv4dy9E8vrIweRHHuvrGr22RtKBUTeSuFafCoxd9OornGhdp06tGzRgl9+3Eadfx7qZX/PE3W9iLyVzJELcTxR18t8r9xmzNT2URepHrYQBZGA/QDIbdR8+/btHKPr83EZePrmfKh2MDKGKmpH4rUZmNRe3E5M4c3/zGDGjBmFun9oaChTpkxhyZIlPP/88zkWenzwwQcMHDiQ9u3bs3HjRp577rlcr3Pnzh2ef/556taty6RJk/K8X1awv/jzFfr2epYAL1ciYu+9lvXgUG8wcuZ6Es4OdgR6uWAy3RthNwnQ8Fzf3vz666/mWRjZHza64YB3PSeSdXqcVXbmL5X7y5/mtmmuECUlOexKLq+paJkmJXq93uLcRJ0hx0/79EwDRy/Fm9+fnqHnwh2dxVS2/BaI7N+/n6lTp3LgwAEGDx6c66o8Dw8Pdu3axciRI3nhhRfo3r07P/30E7du3eL27dscP36cUaNGUa9ePR577DG++uoriwL/eUlJSclRIyRrwc+VuFScHexwVtnh7GBPptFknnbn7u5OcnJyjvdkl32RUH4LXISwJhlhV3J5TUVz8anOiRMnePHFF83napzscvy0j7yVjJdahZuTAyaTicsRZ2jV6WnztLeCFojMmjWLTz/9NN/l6VlWrFiBk5MTZ8+eZfbs2Vy9ehWDwYCfnx+DBw/mr7/+omrVqoXuu7u7OxcuXLA4lvXgMDldj/s//dRlGnBzsjdPu0tISMDd3T3He/JaJAR5p3CEsCYJ2JVcXlPRWnR7jnkjevHOO++g/mcedJC3iiu6TPM52vRM4lIyaFPPG4DIsydJTkqkVdv25tFlfnOToy9FcPPmTfr06WNx//xWC06ePJnmzZtz9epVc7uK68knn2T+/PkWu4lnPZS0VypIyzCgUEJahpEGVd3NQXjVli2MHDkyx3uyfy6S8hDlQVIilVxeU9Ea1KtDmzZtWLFihfm4j9ohx0/75tVdOPxDGB/Nnc78yaPw9Pbll30/4uFkB+SfLggNDWXYsGHY2dmZXytotWBgYCAtWrRg586dJe57o0aNqF+/Plu3bjUfy0pf1K/qSpw2HaPRRC0vZyKi73IgMpZz5/7g7LlzFnl0SXmIikICdiWX3xLyDz74gMWLF7N582bz+Vnzmfs87MPerz5kSr+2HPppJ9euXkahVNDkiY6EfbGUV559giVLlqBxts/1C8FTreLmzZsW23hBwUu2oxPTcKjiz/aj5wosmFQYkyZNYvbs2SQmJlr0cVDLWrzxdDD1q7px9sa9zQtaBLoT+tF82vZ8kTuplgtnsj6XrOXzEqxFeZCAXcnlNzqsX78+O3fuZMKECbz22mtcuXIFgMTERDp27MiNGzdYvuxTfN2d0GmTeeuLrQwZO5lDR35l2/dbCQsL44v5U7mrTc/1C0GpVGIymSzak9+IPGv0rTcYcXPKv65HYfXp04c+ffowcuRIoqMtN1n01zhTRa2iU0NfmldXs+KtiahUDgwe/Vq+NT+EKC8SsB8A+Y0OmzdvzokTJ3B1dSUkJIR27doRHBxMQkICZ86c4e2336bXM0/zw497qV2jGvHajHsPHOsEs2/fPu7G3+bYxo9z/UKoWbMm586ds2hLXikaT7XKPPq+eelvqlYLKFTBpMJYtGgRnTt3Jjj4IXq8EMK7a3ay8/R1bsRr+ev8RTYuf5+Qbi1wUbvy1kdf4q52sigTK0RFIQ8dBdWqVePdd99lwIABbN++nStXrjB79mzq1q3LE088wa27ujxngnz77bfUrVuXOTNnUK2xZanPl19+mQ4dOjB37lwcHe89+MzvAd6ByFhSY69x7dJ5Hm/XyXxOQQWTCqJQKHguZAR1n3qRY7vD+HjmWG7fugmAq7uGdt378OFXW6hZ5176JlmnL1TNDyHKmgRsG1PYeszFoVKpOHHiBHPmzLHYQfxgZCyX41LMc5VrebmaR77dG/szcOBAVq5cyVtvvWVxvYYNG9KoUSNCQ0MZPnw4kP+yc0+1iq+//pzuzw1GpboX4AtbMKkg5+MyqBFQjYdefZ3hr76O0WgkKTUdvUmB3nCvX0aTSWaAiApNUiI2pCj1mIsjJSWF/fv3M3jwYIt7/nrxDgoFuDs5kJ5p5FRUIumZRnPaYOTIkYSGhuZ6zUWLFjFjxgyOHDliPpZXiub0T2GcPPoLXfsPK3bBpLzcvyhIqVTirnbCBDIDRNgMGWHbkIL2USypO3fu4OPjYzH/+fT1RLxdVShQmrfrgkwibyXR8p9RaL169XI80Mvy2GOPsW7dOvr06cOMGTMYMWIEVapUsTjn8uXL/Pe//2Xnzp1s2baDZEdvqxdMym1RUNboXRa9CFshAduG5LUIJi4l3SqpEqVSabHDSdY9G/i5ceb6valvTg5KTEa4o80wj3wNBkO+S8W7devGgQMHWLhwIXXq1OGZZ56hdu3aGAwGTp06xf/+9z9efvlljh07hq+vb5HaXFi5LQqS1IewNRKwbUheS6QVUKT9A/Pi5eVFfHw8d+7cwcvLy3zPtAwDTWt4cCUulei7qSRo9Xi4OJhnb1w4d45atWrle+3GjRuzbt06YmNj2bJlC7du3cLJyYmXXnqJb7/9FhcXF4vzrZ2r91E7EFQv75KtQtgCCdg2JK8ZFg52WCVV4uLiQt++fVm9ejWvv/66xT3dnOyp6eVMTLKOKmpHWtT2NOfQt378qfmhYkF8fX0L3Om7qBvYFpakPoStk4eONiSvRTAmFPlWkyuKsWPHsmzZMlJTU3Pc8+yNu7g72dOqjhfero64OTmQFh/Nnh938/LLL1ulj1C8DWyFeBAUa4SdnJzMtGnTSElJQa/X85///IdHH33U2m0TuchtlFiYanKF1bJlS9q2bcugQYPYtGkTjo6O5ntm5dCV/2x6kHDnNu9NHkafEa+h0WhK1K/s8svVC/EgK9YIe/Xq1bRq1Yp169bx3nvvMW/ePGu3SxRBfvVCikqhUPDFF1/g4OBAx44d2bdvn3l5edYXgz4jg4O7tjLxxWdp2fFpBg8fY9X+5LcaUogHWbFG2EOHDjUXojcYDOZVbKJ8FLQHYlGpVCo2bdrEypUrmTRpEnq9/t6Gt/aOnL14g5OH9xJYN4hhU96mQYsONK1RpcBrFoWUMxUidwUG7LCwML7++muLY++++y6PPPIIt2/fZtq0acycObPUGigKx9oP1JRKJa+88gqjR4/myJEjnDhxAq1WS6cn/Bg4/BVcqwZafaVlFmt/AQlRWShM95dTK6TIyEimTJnC9OnTad++fY7Xw8PDc0zVqgx0Oh1OTk7l3YxSUdy+3dbqOR+XQaLOgMbJjiBvFT5qh4LfWIYq858bVO7+PWh9S01NpXnz5rmeX6yUyIULF5g4cSL/93//R8OGDfM8Lzg4uDiXr5Cy5gX/cfUajer5lcrIsrxFREQU+c8sOjGNE3/G4OlrT41/0hdXdJkE1atYy7uL0zdbUpn796D1LTw8PM/zixWwP/zwQzIyMnjnnXcAcHV1Zfny5cW5lE3IPi+4irPSPP9Yak6U/nJ5IcS/ihWwK3Nwzk32oJT4z7zgrOMPelCSKXhClB1ZOFMI+e2S8qCTKXhClB0J2IUgQSlv1pwDLoTInwTsQsgelEwSlCzIjuJClB0p/lQI2ecFX0ozUk1lJ/OCs5GiSkKUDQnYhZQVlALtEgkO9i/z+5fm1mBCCNsgKREbUNpbgwkhbIOMsG2ANeY6ywhdCNsnI2wbUNJphTJCF6JykIBtA0o6rVA2BBCicpCAbQNKOtdZFv4IUTlIwLYBJZ3rLAt/hKgc5KGjjSjJXGfZEECIyqHCB+yrV6+yYcMGbt68iVKppFatWgwePBgfH5/ybprNkA0BhKgcKmxKJDw8nF69etGsWTOuXbtGrVq1CAgI4OTJk9SvX5+XXnqJCxculHczbYa/xpnujf0Z3DKQ7o39JVgLYYMq5Ah7x44dDB8+nPnz57Nhw4YcO9fEx8ezYsUK2rZty44dO3jsscfKqaUlI3OjhRBFUeEC9okTJxg+fDg7duygRYsWuZ7j6enJzJkzadSoET179uTo0aPUqlWrbBtaQtk3RfB2dUSbnsmeP2NoEuBBdJJOgrgQIocKlxKZPXs2ixYtyjNYZ9e7d29GjBjB+++/XwYts67c5kYbjCbWn4iSBS5CiFxVqID9999/c/LkSQYPHmxxPDoxjd3novnm+FV2n4u2CGDjxo1j/fr1JCUllXVzSyS3udGxyWkYjCZZ4CKEyFWFCtihoaEMGTIER8d/t5wqaFl1tWrV6NChA1u2bCmvZhdLbnOj41Iy8LpvbrQscBFCZKlQOewbN27QsmVLi2OFKXzUoEEDbt68WbaNLQGTyYT2yhlWffMdGanJuLg4U71OA5wbPomvv7vFubLARQiRpUIFbLgXzLIrzCav97+nIluzZg2LFi0C4Nm+z6Ozr09CUgpnjv3MuU8W0rJLb4a8Og1fby+LBS4yo0QIUaECdo0aNfjzzz8tjmWlDrJG1pBz1BkREUH//v3LrJ3FYTKZmD59Ojt37mTZsmW0b98ehUKRLRAPw5B8m61fLuXNkc8z/aO11K5Z3bwaMbcZJbktT5fALkTlVaFy2CEhIaxbt460tH8fKhZU+CgqKorDhw/Tp0+f8ml0IS1ZsoSffvqJI0eO0KFDB3Owzp6fd/X0o9f4t+nTtw9fvjWWjkGe+GucC11tr6zKqOb3EFgIUXoqVMCuU6cOLVu2ZM2aNeZjBRU+Wrp0KS+99BKurq7l1ewCpaSksGDBArZu3UqVKlXMx3MLxO7ODnQNmYCbmxvffvstUPhqe2VRRlVqawtRfipUSgRgwYIFdO3alUaNGtG2bVsg78JHGzdu5JtvvuHYsWNl3cwiCQ0NpX379tSuXdvi+OXbWpJ0GaSkG3BzsqeWlysaFwfiUtJ57bXXeP/993nxxRcLlRaCwuX7S8oau98IIYqnRCPsixcv0rx5c9LTrRcQHn30Ub755hv69u3LkiVLuHv3bo5zbt26xaxZs5gyZQo7d+4kICDAavcvDaGhoYwcOdLiWHRiGlfjU0nSZeLu5EB6ppFTUYncSEjDU63i2Wef5dKlS1y6dKnQ9bDLooyq1NYWovwUO2CnpKSwaNEiVCrrTznr0qUL+/fv5/jx49SqVYthw4axYMEC5s2bR//+/QkODiYuLo6jR4/SpEkTq9/f2m7evElQUJDFsdPXE6lf1RWTSYFOb8TJ3g6FwkRkTBJNAjTY29tTu3ZtoqOjC10Pu6QbHRSG1NYWovwUKyViMpl48803mTJlCuPGjbN2mwB4+OGH2bBhA7du3SIsLIzo6GiUSiVdunTh888/R6PRlMp9S4NSqcx1umJAFRdcHe25EpdKkk6Pm6MD7s725kBsMplQKu99pxamHnZZlFGV2tpClJ8CA3ZYWBhff/21xbFq1arRo0cPGjZsWGoNy+Ln58eECRNK/T6lqWbNmpw7d44GDRqYj2WNVD3Vjniq7+Wdk3V6nFV2AKSnp3PhwgVq1KhRpHuVZKODwl5famsLUT4UpmKsOunSpQt+fn4AnDp1ikceeYTQ0FCLc8LDw3OURa0MdDodTk5ORXrP9u3b+f777/niiy/Mx25r9RyLSkXtoMTZQUGa3oRWb6RVDRd81A65vqe0FadvtqIy9w0qd/8etL6lpqbSvHnz3N9gKqGOHTuadDpdjuO//fZbSS9dIf35559Ffo9OpzP5+vqaTp06ZXH8ZkKqadfZm6bQY1dMu87eNN1MSDWZTCaTXq83PfbYY6YtW7ZYo8mFVpy+2YrK3DeTqXL370HrW36xs8JN68uNra/ec3R05IMPPqBPnz78/PPP1KxZE8g9fWEwGBg7diyenp707NmzPJorhKigSrxwZv/+/RbV9aytsizUCAkJYdKkSbRp04Z169ah0+lynHPixAl69erFhQsXCAsLw87OrhxaKoSoqCr8CLsyLdSYOHEijRo1YvHixUyZMoV+/frh6+tLWloaBw8eJC4ujrFjxzJx4sRS/RIUQtimCh+wy2L1Xlnq3LkznTt35vz58/zwww8kJCTg7e3N3Llz6datm4yqhRB5qvABu7DLsm2Nq08ADZ7qb5GXl2AthMhPhSr+lJuyWL1X1ipLXl4IUbYqfMAu7LJsW1IWVfWEEJVPhU+JQOmv3itrlS0vL4QoGxV+hF0ZSQElIURxSMAuB5UxLy+EKH0SsMtBZczLCyFKn03ksCujypaXF0KUPhlhCyGEjZCALYQQNkICthBC2AgJ2EIIYSMkYAshhI14oGeJ2PrGCEKIB4tNB+ySBNysAkxuTvZ4uzqiTc9kz58xMh9aCFFh2WxKpKQV76QAkxDC1thswC5pwI3XZqB2tPyBoXa0J16bUQqtFUKIkqtQKZGipDhKWvGusm6MIISovCrMCLuoKY6SVryTAkxCCFtTYQJ2UVMcJQ24UoBJCGFrKkxKpKgpjqyAe/p6InEp6XiqVbSq41WkgCsFmIQQtqTCBOzi5JQl4AohHiQVJiUiOWUhhMhfsQK2wWBgwYIFDBw4kH79+nHgwIESN0RyykIIkb9ipUS+//57MjMz2bBhAzExMezatcsqjZEUhxBC5K1YAfvw4cMEBQUxevRoTCYTb775prXbJYQQ4j4FBuywsDC+/vpri2NVqlTB0dGRzz77jBMnTvDGG28QGhpaao0UQggBCpPJZCrqmyZPnkz37t3p1q0bAG3atOHIkSMW54SHh+Pi4mKdVlYgOp0OJyen8m5GqZC+2a7K3L8HrW+pqak0b9481/OLlRJp3rw5P//8M926deOvv/7C398/1/OCg4OLc/kKLSIiolL2C6Rvtqwy9+9B61t4eHie5xdrlkj//v0xmUz079+fN998k7lz5xbnMkIIIYqgWCmRwsjvW0IIIUTe8kqJlFrAFkIIYV0VZqWjEEKI/EnAFkIIGyEBu4iSk5MZM2YML730EgMGDODkyZPl3SSr27NnD1OnTi3vZliF0Whkzpw5DBgwgJCQEK5evVreTbK606dPExISUt7NsCq9Xs+0adMYPHgwzz//PPv27SvvJlmVwWDgjTfeYODAgQwaNIi///67UO+rMNX6bMXq1atp1aoVQ4cO5dKlS0ydOpUtW7aUd7OsZsGCBRw+fLjSTKPau3cvGRkZbNy4kVOnTrFw4UKWL19e3s2ympUrV7Jt2zacnStXSYdt27ah0WhYvHgxiYmJ9OnTh6eeeqq8m2U1WfWXNmzYwPHjx1myZEmh/l5KwC6ioUOHolLdK/lqMBhwdHQs4B22pVmzZnTu3JmNGzeWd1OsIjw8nHbt2gHQtGlTzp07V84tsq6aNWuydOlSpk+fXt5NsarsC/NMJhN2dnbl3CLr6ty5Mx06dADg5s2buLu7F+p9ErDzkduy/HfffZdHHnmE27dvM23aNGbOnFlOrSuZvPrWo0cPjh8/Xk6tsr6UlBRcXV3N/21nZ0dmZib29pXjr363bt24fv16eTfD6tRqNXDvz++1115j0qRJ5dugUmBvb8+MGTPYs2cPH3/8ceHeU8ptsmkvvPACL7zwQo7jkZGRTJkyhenTp9OiRYtyaFnJ5dW3ysbV1RWtVmv+b6PRWGmCdWUXHR3N+PHjGTx4MD179izv5pSKRYsW8frrr9O/f3927txZYDkPeehYRBcuXGDixIl8+OGHtG/fvrybIwrQrFkzfvnlFwBOnTpF/fr1y7lFojDi4uIYPnw406ZN4/nnny/v5ljd1q1b+eyzzwBwdnZGoVCgVBYcjmWoUUQffvghGRkZvPPOO8C9EVxleohV2XTp0oUjR44wcOBATCYT7777bnk3SRTCihUrSEpKYtmyZSxbtgy494C1shSB6tq1K2+88QYvvvgimZmZzJw5s1B9k5WOQghhIyQlIoQQNkICthBC2AgJ2EIIYSMkYAshhI2QgC2EEDZCArYQQtgICdhCCGEjJGALIYSN+H9gFra5Dx1ZjgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X[:, 0], X[:, 1], alpha=0.3)\n",
    "plt.scatter(selection[:, 0], selection[:, 1],\n",
    "            facecolor='none', edgecolor='black', s=200);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This sort of strategy is often used to quickly partition datasets, as is often needed in train/test splitting for validation of statistical models (see [Hyperparameters and Model Validation](05.03-Hyperparameters-and-Model-Validation.ipynb)), and in sampling approaches to answering statistical questions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modifying Values with Fancy Indexing\n",
    "\n",
    "Just as fancy indexing can be used to access parts of an array, it can also be used to modify parts of an array.\n",
    "For example, imagine we have an array of indices and we'd like to set the corresponding items in an array to some value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0 99 99  3 99  5  6  7 99  9]\n"
     ]
    }
   ],
   "source": [
    "x = np.arange(10)\n",
    "i = np.array([2, 1, 8, 4])\n",
    "x[i] = 99\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use any assignment-type operator for this. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0 89 89  3 89  5  6  7 89  9]\n"
     ]
    }
   ],
   "source": [
    "x[i] -= 10\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice, though, that repeated indices with these operations can cause some potentially unexpected results. Consider the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[6. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n"
     ]
    }
   ],
   "source": [
    "x = np.zeros(10)\n",
    "x[[0, 0]] = [4, 6]\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where did the 4 go? This operation first assigns `x[0] = 4`, followed by `x[0] = 6`.\n",
    "The result, of course, is that `x[0]` contains the value 6.\n",
    "\n",
    "Fair enough, but consider this operation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([6., 0., 1., 1., 1., 0., 0., 0., 0., 0.])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i = [2, 3, 3, 4, 4, 4]\n",
    "x[i] += 1\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You might expect that `x[3]` would contain the value 2 and `x[4]` would contain the value 3, as this is how many times each index is repeated. Why is this not the case?\n",
    "Conceptually, this is because `x[i] += 1` is meant as a shorthand of `x[i] = x[i] + 1`. `x[i] + 1` is evaluated, and then the result is assigned to the indices in `x`.\n",
    "With this in mind, it is not the augmentation that happens multiple times, but the assignment, which leads to the rather nonintuitive results.\n",
    "\n",
    "So what if you want the other behavior where the operation is repeated? For this, you can use the `at` method of ufuncs and do the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 0. 1. 2. 3. 0. 0. 0. 0. 0.]\n"
     ]
    }
   ],
   "source": [
    "x = np.zeros(10)\n",
    "np.add.at(x, i, 1)\n",
    "print(x)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `at` method does an in-place application of the given operator at the specified indices (here, `i`) with the specified value (here, 1).\n",
    "Another method that is similar in spirit is the `reduceat` method of ufuncs, which you can read about in the [NumPy documentation](https://numpy.org/doc/stable/reference/ufuncs.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Binning Data\n",
    "\n",
    "You could use these ideas to efficiently do custom binned computations on data.\n",
    "For example, imagine we have 100 values and would like to quickly find where they fall within an array of bins.\n",
    "We could compute this using `ufunc.at` like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(seed=1701)\n",
    "x = rng.normal(size=100)\n",
    "\n",
    "# compute a histogram by hand\n",
    "bins = np.linspace(-5, 5, 20)\n",
    "counts = np.zeros_like(bins)\n",
    "\n",
    "# find the appropriate bin for each x\n",
    "i = np.searchsorted(bins, x)\n",
    "\n",
    "# add 1 to each of these bins\n",
    "np.add.at(counts, i, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The counts now reflect the number of points within each bin—in other words, a histogram (see the following figure):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWsAAAD0CAYAAABdAQdaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAARNklEQVR4nO3df0zV9aPH8dcB/AXsXNa0RXpBK0vQGTec2B/iWhHMzWm7KIrSFNbMsatMQ5QhWhTKZdXKfVHzj9Yyr1q3kj+aW7E2NnFs90x0wsnKzG7ImvYdO8CJH8rn/uH98pXfcM75cHhzno+/4sP5fHi9+Rxfvvv4Pp+Pw7IsSwCASS0s2AEAAKOjrAHAAJQ1ABiAsgYAA1DWAGAAyhoADBBh14FdLpddhwaAKS05OXnQNtvKergfOJm53W4lJCQEO8aEYsyhgTGbY7iJLpdBAMAAlDUAGICyBgADUNYAYADKGgAMQFkDgAEoawAwgK3rrAHTnK7/Tecbmv06xtqkucpOiQtQIuABZtbAQ843NKupxePz/k0tHr/LHhgKM2tggMRYp85uf96nfbNOXApwGuABZtYAYADKGgAMQFkDgAEoawAwAGUNAAagrAHAAJQ1ABhg1HXW9+/fV0lJiW7evCmHw6E333xTM2bM0L59++RwOLRw4UIdPHhQYWH0PgDYZdSy/v777yVJZ86cUX19vd5//31ZlqWCggKlpKSotLRUNTU1SktLsz0sAISqUafDL730ksrKyiRJt2/fltPpVGNjo5YvXy5JSk1NVV1dnb0pASDEjenj5hERESoqKtK3336rDz/8UBcvXpTD4ZAkRUVFqa2tbcj93G534JJOgM7OTuMy+4sx9+f1eiX5/t71d3+7cJ7NN+Z7g1RUVOiNN97Qhg0b1NXV1be9o6NDTqdzyH1Me7KwqU9D9gdj7i+ytlWS7+9df/e3C+fZHD4/3fzrr7/WiRMnJEmzZs2Sw+HQkiVLVF9fL0mqra3VsmXLAhgVADDQqDPrl19+Wfv379fmzZt17949FRcX68knn9SBAwf03nvv6YknnlB6evpEZAWAkDVqWUdGRuqDDz4YtP3UqVO2BAIADMbiaAAwAGUNAAagrAHAAJQ1ABiAsgYAA1DWAGAAyhoADEBZA4ABKGsAMABlDQAGoKwBwACUNQAYgLIGAANQ1gBgAMoaAAxAWQOAAShrADAAZQ0ABqCsAcAAlDUAGICyBgADUNYAYADKGgAMEDHSN3t6elRcXKzm5mZ1d3drx44dio2N1fbt2zV//nxJ0qZNm7R69eqJyAoAIWvEsq6urlZMTIwqKyvV2tqqdevWKT8/X9u2bVNubu5EZQSAkDdiWWdkZCg9PV2SZFmWwsPDde3aNd28eVM1NTWKj49XcXGxoqOjJyQsAIQqh2VZ1mgvam9v144dO7RhwwZ1d3frmWee0ZIlS3Ts2DF5PB4VFRUN2sflcikyMtKW0Hbp7OzUzJkzgx1jQjHm/vZeuC1J+s+Mx306tr/724XzbA6v16vk5ORB20ecWUtSS0uL8vPzlZ2drTVr1sjj8cjpdEqS0tLSVFZWNuy+CQkJfkSeeG6327jM/mLM/UXWtkry/b3r7/524Tybw+VyDbl9xNUgd+/eVW5urgoLC5WZmSlJysvL09WrVyVJly5d0uLFiwMcFQAw0Igz6+PHj8vj8aiqqkpVVVWSpH379qm8vFzTpk3T7NmzR5xZAwACY8SyLikpUUlJyaDtZ86csS0QAGAwPhQDAAagrAHAAJQ1ABiAsgYAA1DWAGAAyhoADEBZA4ABKGsAMABlDQAGoKwBwACUNQAYgLIGAANQ1gBgAMoaAAxAWQOAAShrADAAZQ0ABqCsAcAAlDUAGICyBgADUNYAYADKGgAMQFkDgAEiRvpmT0+PiouL1dzcrO7ubu3YsUNPPfWU9u3bJ4fDoYULF+rgwYMKC6PzAcBOI5Z1dXW1YmJiVFlZqdbWVq1bt06LFi1SQUGBUlJSVFpaqpqaGqWlpU1UXgAISSNOiTMyMrRr1y5JkmVZCg8PV2Njo5YvXy5JSk1NVV1dnf0pASDEjTizjoqKkiS1t7dr586dKigoUEVFhRwOR9/329raht3f7XYHMKr9Ojs7jcvsL8bcn9frleT7e9ff/e3CeTbfiGUtSS0tLcrPz1d2drbWrFmjysrKvu91dHTI6XQOu29CQkJgUk4Qt9ttXGZ/Meb+ImtbJfn+3vV3f7twns3hcrmG3D7iZZC7d+8qNzdXhYWFyszMlCQlJiaqvr5eklRbW6tly5YFOCoAYKARy/r48ePyeDyqqqpSTk6OcnJyVFBQoKNHjyorK0s9PT1KT0+fqKwAELJGvAxSUlKikpKSQdtPnTplWyAAwGAskAYAA1DWAGAAyhoADEBZA4ABKGsAMABlDQAGoKwBwACUNQAYgLIGAANQ1gBgAMoaAAxAWQOAAShrADAAZQ0ABqCsAcAAlDUAGICyBgADUNYAYADKGgAMQFkDgAEoawAwAGUNAAagrAHAAGMq6ytXrignJ0eS1NTUpJUrVyonJ0c5OTn65ptvbA0IAJAiRnvByZMnVV1drVmzZkmSGhsbtW3bNuXm5toeDgDwwKhlHRcXp6NHj2rv3r2SpGvXrunmzZuqqalRfHy8iouLFR0dbXtQmOF0/W8639Ds1zHWJs1VdkpcgBIBU8OoZZ2enq7ff/+97+ulS5dq/fr1WrJkiY4dO6a//e1vKioqGnJft9sduKQToLOz07jM/gr0mP+r7rZ++Xu3nnhkuk/7//L3bnm9Xv2bsyNgmQYaacxer1eS7+9df/e3C+9t841a1gOlpaXJ6XT2/XdZWdmwr01ISPA9WRC43W7jMvsr0GOOrG3VkshInd3+vE/7Z524JMne985IY46sbfXr5/u7v114b5vD5XINuX3cq0Hy8vJ09epVSdKlS5e0ePFi/5IBAEY17pn1oUOHVFZWpmnTpmn27NkjzqwBAIExprKeN2+ezp07J0lavHixzpw5Y2soAEB/fCgGAAxAWQOAAShrADAAZQ0ABqCsAcAAlDUAGGDc66wBuzW1ePo+yegL7i2CqYiyxqSyNmmuX/s3tXgkibLGlENZY1LJTonzq2j9mZEDkxnXrAHAAJQ1ABiAsgYAA3DNGlPOaKtJvF5v332nh9o3MdZpUzLAd5Q1phR/V5Mkxjr9PgZgB8oaU8pYVpOY+gQRhDauWQOAAShrADAAZQ0ABqCsAcAAlDUAGICyBgADUNYAYIAxlfWVK1eUk5MjSbp165Y2bdqk7OxsHTx4UL29vbYGBACMoaxPnjypkpISdXV1SZIOHz6sgoICnT59WpZlqaamxvaQABDqRi3ruLg4HT16tO/rxsZGLV++XJKUmpqquro6+9IBACSNoazT09MVEfHPT6VbliWHwyFJioqKUltbm33pAACSfLg3SFjYP/u9o6NDTufwdyhzu92+pQqSzs5O4zL7K9Bj9nq9kib3ubfzPE/W8fPeNt+4yzoxMVH19fVKSUlRbW2tVqxYMexrTbtZTije4CfQY/7HrUcn8+/RzvM8WcfPe9scLpdryO3jXrpXVFSko0ePKisrSz09PUpPT/c7HABgZGOaWc+bN0/nzp2TJC1YsECnTp2yNRQAoD8+FAMABqCsAcAAlDUAGICyBgADUNYAYADKGgAMQFkDgAHG/QlGAPY5Xf+bzjc0+3WMtUlzlZ0SF6BEmCyYWQOTyPmGZjW1eHzev6nF43fZY3JiZg1MMomxTp3d/rxP+2aduBTgNJgsmFkDgAEoawAwAGUNAAbgmjX68Xc1QlOLR4mxwz+QIhQ0tXh8vnbM7w/DYWaNfvxdjZAY69TapLkBTGSWtUlz/SrbUP/9YXjMrDGIP6sRQl12ShxrnGELZtYAYADKGgAMQFkDgAEoawAwAGUNAAagrAHAAJQ1ABjA53XWr7zyiqKjoyVJ8+bN0+HDhwMWCgDQn09l3dXVJcuy9OmnnwY6DwBgCD5dBvnhhx/0119/KTc3V6+++qoaGhoCHAsA8DCfZtYzZ85UXl6e1q9fr19//VWvvfaaLly4oIgIPr0OAHbwqV0XLFig+Ph4ORwOLViwQDExMbpz545iY2P7vc7tdgck5ETp7Ow0LrO/Bo7Z6/VKMu/cjcdUPs/Dnb+pPObhTLUx+1TWX3zxhX788UcdOnRIf/zxh9rb2zVnzpxBr0tISPA74ERyu93GZfbXwDFH1rZKMu/cjcdUPs/Dnb+pPObhmDpml8s15HafyjozM1P79+/Xpk2b5HA4VF5eziUQALCRTw07ffp0vfvuu4HOAgAYBh+KAQADUNYAYADKGgAMQFkDgAEoawAwAGUNAAagrAHAAJQ1ABiAsgYAA1DWAGAAyhoADEBZA4ABKGsAMAD3NQWmmKYWj7JOXOq3zev19t3r2m5rk+YqOyVuQn5WKKGsgSlkbdLcoP78phaPJFHWNqCsgSkkOyVuyKKcqKemDJzRI3C4Zg0ABqCsAcAAlDUAGIBr1lPM6frfdL6hecyvH7hKoKnFo8RYpw3JECqGWo0yHqwmGRoz6ynmfENz37/I+yIx1hn0FQUw19qkuX79Zd/U4hnXZCOUMLOeghJjnTq7/fkxvXaiVgkgNAy3GmWsWE0yPGbWAGAAn2bWvb29OnTokK5fv67p06fr7bffVnx8fKCzAQD+n08z6++++07d3d06e/as9uzZoyNHjgQ6FwDgIT7NrF0ul1auXClJSkpK0rVr1wIW6L9dv+vc//xvwI43HhN5/wS7sJoDpvN3Nck/BOvP84Zl/6p/T54X8OP6VNbt7e2Kjo7u+zo8PFz37t1TRET/w7nd7nEf+3ZLm7xery+x/Nbb2xu0nx0o8/8lQssfCxvz776zs9On82Qyxjx5LX8sTF5vRED+HAbrz/Ptlttyu9sCflyfyjo6OlodHR19X/f29g4qakk+rTJISJD+w5dQARCKKyMYc2gwZcwJCdKeAB3LlDEP5HK5htzu0zXr5557TrW1tZKkhoYGPf30074nAwCMyqeZdVpami5evKiNGzfKsiyVl5cHOhcA4CE+lXVYWJjeeuutQGcBAAyDD8UAgAEoawAwAGUNAAagrAHAAJQ1ABjAYVmWZceBh1vYDQAYWXJy8qBttpU1ACBwuAwCAAagrAHAAJT1EG7cuKHk5GR1dXUFO4rt2tra9Prrr2vLli3KysrS5cuXgx3JNr29vSotLVVWVpZycnJ069atYEeyXU9PjwoLC5Wdna3MzEzV1NQEO9KE+PPPP7Vq1SrduHEj2FEChmcwDtDe3q6KigpNnz492FEmxMcff6wVK1Zo69at+uWXX7Rnzx599dVXwY5li4cfmtHQ0KAjR47o2LFjwY5lq+rqasXExKiyslKtra1at26dXnzxxWDHslVPT49KS0s1c+bMYEcJKGbWD7EsSwcOHNDu3bs1a9asYMeZEFu3btXGjRslSffv39eMGTOCnMg+dj40Y7LKyMjQrl27JD14f4eHhwc5kf0qKiq0ceNGPfroo8GOElAhO7P+/PPP9cknn/Tb9vjjj2v16tVatGhRkFLZa6gxl5eXa+nSpbpz544KCwtVXFwcpHT2G+tDM6aSqKgoSQ/GvnPnThUUFAQ3kM2+/PJLPfLII1q5cqU++uijYMcJKJbuPSQtLU2PPfaYpAf36V66dKk+++yzIKey3/Xr17V7927t3btXq1atCnYc2xw+fFjPPvusVq9eLUlKTU3tuy/7VNbS0qL8/Py+69ZT2ebNm+VwOORwOOR2uzV//nwdO3ZMc+bMCXY0/1kY0gsvvGB1dnYGO4btfvrpJys9Pd1yu93BjmK7CxcuWEVFRZZlWdbly5etvLy8ICey3507d6yMjAyrrq4u2FEm3JYtW6yff/452DECZur+/x/G5N1331V3d7feeecdSQ8e2TZV/9EtFB+acfz4cXk8HlVVVamqqkqSdPLkySn3j2+hgMsgAGAAVoMAgAEoawAwAGUNAAagrAHAAJQ1ABiAsgYAA1DWAGAAyhoADPB/B94FbuHEJ84AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the results\n",
    "plt.plot(bins, counts, drawstyle='steps');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, it would be inconvenient to have to do this each time you want to plot a histogram.\n",
    "This is why Matplotlib provides the `plt.hist` routine, which does the same in a single line:\n",
    "\n",
    "```python\n",
    "plt.hist(x, bins, histtype='step');\n",
    "```\n",
    "\n",
    "This function will create a nearly identical plot to the one just shown.\n",
    "To compute the binning, Matplotlib uses the `np.histogram` function, which does a very similar computation to what we did before. Let's compare the two here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NumPy histogram (100 points):\n",
      "33.8 µs ± 311 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
      "Custom histogram (100 points):\n",
      "17.6 µs ± 113 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
     ]
    }
   ],
   "source": [
    "print(f\"NumPy histogram ({len(x)} points):\")\n",
    "%timeit counts, edges = np.histogram(x, bins)\n",
    "\n",
    "print(f\"Custom histogram ({len(x)} points):\")\n",
    "%timeit np.add.at(counts, np.searchsorted(bins, x), 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our own one-line algorithm is twice as fast as the optimized algorithm in NumPy! How can this be? If you dig into the `np.histogram` source code (you can do this in IPython by typing `np.histogram??`), you'll see that it's quite a bit more involved than the simple search-and-count that we've done; this is because NumPy's algorithm is more flexible, and particularly is designed for better performance when the number of data points becomes large:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NumPy histogram (1000000 points):\n",
      "84.4 ms ± 2.82 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "Custom histogram (1000000 points):\n",
      "128 ms ± 2.04 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "x = rng.normal(size=1000000)\n",
    "print(f\"NumPy histogram ({len(x)} points):\")\n",
    "%timeit counts, edges = np.histogram(x, bins)\n",
    "\n",
    "print(f\"Custom histogram ({len(x)} points):\")\n",
    "%timeit np.add.at(counts, np.searchsorted(bins, x), 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What this comparison shows is that algorithmic efficiency is almost never a simple question. An algorithm efficient for large datasets will not always be the best choice for small datasets, and vice versa (see [Big-O Notation](02.08-Sorting.ipynb#Big-O-Notation)).\n",
    "But the advantage of coding this algorithm yourself is that with an understanding of these basic methods, the sky is the limit: you're no longer constrained to built-in routines, but can create your own approaches to exploring the data.\n",
    "Key to efficiently using Python in data-intensive applications is not only knowing about general convenience routines like `np.histogram` and when they're appropriate, but also knowing how to make use of lower-level functionality when you need more pointed behavior."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "jupytext": {
   "formats": "ipynb,md"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}