Assignment2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Assignemnt #2 - find genes with differential gene expression patterns. \n",
    "\n",
    "We want to find the genes that have statistically significant patterns of their gene expression levels. \n",
    "\n",
    "### Part 1 - all vs all \n",
    "For all gene pairs perform a test of significance. \n",
    "\n",
    "How do you vizualize the results? \n",
    "\n",
    "Do you need to account for multiple testuing (yes, you do...)\n",
    "\n",
    "### Part 2 - spatially differentual gene expression\n",
    "\n",
    "Define a measure that captures the idea of spatially varying gene expression. Think about it, there are many ways to do that. \n",
    "\n",
    "Hint: If expression pattern is uniform, what do you expect from two cells that are neighbors compared to cells that are not? \n",
    "\n",
    "For both parts, run first pass without taking into account cell types, and then chosing type level \"class\" repeat the analysis subseting on cells from the same type. Adjust code as needed. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Useful tool - KNN distances\n",
    "\n",
    "to find if expression patterns are local you need to know which cells are the neighbors of a given cell. \n",
    "We do that with sklean package (k nearest neighbors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.neighbors import NearestNeighbors\n",
    "np.random.seed(0)\n",
    "\n",
    "mu, sigma = 10, 5\n",
    "xy = np.random.normal(mu, sigma, (100, 2))\n",
    "\n",
    "\n",
    "k = 5  # or any other number you prefer\n",
    "nbrs = NearestNeighbors(n_neighbors=k, algorithm='ball_tree').fit(xy)\n",
    "distances, indices = nbrs.kneighbors(xy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.        , 0.10031233, 0.17354638, 0.33683095, 0.42193949],\n",
       "       [0.        , 0.22299133, 0.26091184, 0.35778274, 0.61538469],\n",
       "       [0.        , 1.76534673, 2.27033892, 2.37712608, 2.43643506],\n",
       "       ...,\n",
       "       [0.        , 0.02609462, 0.34825978, 0.38385796, 0.43303073],\n",
       "       [0.        , 0.18071092, 0.69765965, 0.80616595, 0.93103888],\n",
       "       [0.        , 0.20569891, 0.30451583, 0.32036085, 0.51223411]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  0, 788, 226, 965, 336],\n",
       "       [  1, 959, 614, 445, 444],\n",
       "       [  2, 204, 494, 412, 746],\n",
       "       ...,\n",
       "       [997,  38, 953, 549, 214],\n",
       "       [998, 403, 282, 831, 185],\n",
       "       [999, 773, 754, 938, 118]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA3/0lEQVR4nO3df3xU1Z3/8fckQII0mX4DwiQaaJbC1hhXUItgQSnfhYW2WRS/+/VH7dLtLq2I7FLWhapVoLpQ3JXt7lJpy6O1tYi6P6rI1kVpqSArLFikCrT9IsaCmGxKsDMYTSLJ/f4xTsgkM8n8uHPvufe+no9HHg8yuSRnJjdzPuecz/mckGVZlgAAABxS5HYDAABAsBB8AAAARxF8AAAARxF8AAAARxF8AAAARxF8AAAARxF8AAAARxF8AAAARw1yuwG9dXV16a233lJZWZlCoZDbzQEAABmwLEtnzpxRVVWVior6n9swLvh46623VF1d7XYzAABADk6cOKELL7yw32uMCz7KysokxRtfXl7ucmsAAEAmYrGYqquru/vx/hgXfCSWWsrLywk+AADwmExSJkg4BQAAjiL4AAAAjiL4AAAAjiL4AAAAjiL4AAAAjiL4AAAAjiL4AAAAjiL4AAAAjjKuyBgAoHA6uyztazit5jNtGllWqkk1FSou4hwtOIvgAwACYtuhRq3aekSN0bbuxyrDpVpRX6vZdZUutgxBw7ILAATAtkONWrjpQFLgIUlN0TYt3HRA2w41utQyBBHBBwD4XGeXpVVbj8hK8bXEY6u2HlFnV6orAPsRfACAz+1rON1nxqMnS1JjtE37Gk471ygEGsEHAPhc85n0gUcu1wH5IvgAAJ8bWVZq63VAvgg+AMDnJtVUqDJcqnQbakOK73qZVFPhZLMQYAQfAOBzxUUhraivlaQ+AUji8xX1tdT7gGMIPgAgAGbXVWrDLZcpEk5eWomES7Xhlsuo8wFHUWQMAAJidl2lZtZGqHAK1xF8AECAFBeFNGXscLebgYBj2QUAADiK4AMAADiK4AMAADiK4AMAADiK4AMAADiK4AMAADiK4AMAADgqq+BjzZo1+vjHP66ysjKNHDlS1157rX79618nXWNZllauXKmqqioNHTpU06dP1+HDh21tNAAA8K6sgo+dO3dq0aJF2rt3r7Zv366zZ89q1qxZam1t7b7mgQce0Lp167R+/Xrt379fkUhEM2fO1JkzZ2xvPAAA8J6QZVlWrv/5t7/9rUaOHKmdO3fq6quvlmVZqqqq0pIlS7R8+XJJUnt7u0aNGqW1a9fqS1/60oDfMxaLKRwOKxqNqry8PNemAQAAB2XTf+eV8xGNRiVJFRXxY5gbGhrU1NSkWbNmdV9TUlKia665Ri+++GI+Pwou6OyytOdYi7YcPKk9x1rU2ZVznAoAQLecz3axLEtLly7V1KlTVVdXJ0lqamqSJI0aNSrp2lGjRuk3v/lNyu/T3t6u9vb27s9jsViuTYKNth1q1KqtR9QYbet+rDJcqhX1tZ46/bKzy+IQLQAwTM7Bx+23365XXnlFu3fv7vO1UCj5zd2yrD6PJaxZs0arVq3KtRkogG2HGrVw0wH1nudoirZp4aYDnjl+2y8BVC4IugCYLKfgY/HixXr66ae1a9cuXXjhhd2PRyIRSfEZkMrKc2/uzc3NfWZDEu68804tXbq0+/NYLKbq6upcmgUbdHZZWrX1SJ/AQ5IsSSFJq7Ye0czaiNGdmZsBlNsdf5CDLgDekFXwYVmWFi9erCeffFLPP/+8ampqkr5eU1OjSCSi7du3a+LEiZKkjo4O7dy5U2vXrk35PUtKSlRSUpJj82G3fQ2nkzqt3ixJjdE27Ws4beyx3G4GUG53/H6ZtQJ6czuoh72yCj4WLVqkzZs3a8uWLSorK+vO8QiHwxo6dKhCoZCWLFmi1atXa9y4cRo3bpxWr16t8847TzfffHNBngDs1XwmfeCRy3VucCuAcrvj98usFdCb20E97JfVbpcNGzYoGo1q+vTpqqys7P544oknuq9ZtmyZlixZottuu01XXHGFTp48qeeee05lZWW2Nx72G1lWaut1bnAjgBqo45fiHX8hdwxlE3QBXpEI6nvf24mgftuhRpdahnxkvewykFAopJUrV2rlypW5tgkumlRTocpwqZqibSk70pCkSDg+5WkqNwIoE5ar/DBrBfTEbJ5/cbYLkhQXhbSivlZS/A+7p8TnK+prjf5DTwRQ6VoYUnzK1s4AyoSOP9Ng6uj/nKFuCzyB2Tz/IvhAH7PrKrXhlssUCSd3ZpFwqScSFt0IoExYrhoo6EpY/7NjumnjXk1du4MpaxjNhKAehZFznQ/42+y6Ss2sjXg2uzwRQPVOUosUKEnNhOWqRNC1cNMBhaSU7eiJHTAwnQlBPQojr7NdCoGzXWAnJ7fnJRLjpOSOP/HTnOrkU+0MSCcRFO1ePsMzgSWCo7PL0tS1OwYM6rl/zeDY2S6A6YqLQpoydrjmTrhAU8YOL+gblCnLVbPrKrV7+Qw9tmCybv/kR/u9ljVzmMwPOWhIjWUXwEamLFclgi6vrZlTSAq9Ob2ECmcQfAA2S3T8JvDSmjmFpJCOKUE97EPwAfiYCYmwmXC7OizMZ1JQj/yR8wH4mBfWzE2oDgvAWQQfgM+ZkgibDoWkgOBh2QUIAJPWzHsnlTZF38vo/5mSFAsgfwQfQECYsGaeKqm0YtiQjP6vCUmxAOxB8AHAEemSSt9u7ej3/5mSFAvAPuR8ACi4TJJKUzElKRaAvQg+ABTcQEmlCRXDBid9bkpSLAB7sewCoOAyTRa95zMXK1Je6npSLIDCIvgAUHCZJotGyktdT4oFUHgsuwAouESl1XRzGCHFS6mTVAoEA8EHgILzQqVVAM4h+ADgCNMrrQJwDjkfABxjUqVVAO4h+ADgKBMqrQJwF8suAADAUcx8AC7qfcgaSxAAgoDgA3BJqkPWKsOlWlFfS/IlAF9j2QVwQeKQtd4lx5uibVq46YC2HWp0qWUAUHgEH4DDMjlkbdXWI+rs6u/INQDwLoIPwGEDHbJmSWqMtmlfw2nnGgUADiLnA3BYpoesZXodAPeRPJ4dgg/AYZkespbpdUAQmdTZkzyePYIPwGGJQ9aaom0p8z5Cipcc55C1gZnUAcE5JnX2ieTx3n/LieRxjg5IjeADcFjikLWFmw4oJCW9aXHIWuZM6oDgHJM6+4GSx0OKJ4/PrI3w99wLCaeACzhkLT9sVQ4m03aKkTyeO2Y+AJdwyFpuGG0GVzadvRPnB5E8njuCD8BFHLKWPdM6IDjHtM6e5PHcsewCwFNM64DgHNM6+0TyeLr5tZDieUgkj/dF8AHAU0zrgOAc0zr7RPJ44mf3botE8ng6BB8APMW0DgjOMbGzJ3k8NyHLsow6QCIWiykcDisajaq8vNzt5gAwUGK3i5R6qzJv+v5m4jZras5k138TfADwpHw6IDoK7wvK79BLz5PgA0Ag5PLGbOKoGUjFa/cqwQcApJCuOibLNTCNF+/VbPpvEk4BBIJp1TGBdIJwrxJ8AAgESmHDK4JwrxJ8AAgEipPBK4JwrxJ8AAgEipPBK4JwrxJ8AAgEipPBK4JwrxJ8AAHU2WVpz7EWbTl4UnuOtXg6cS1TJlbHBFIJwr3KVlsgYLxWO8BuQX/+8A6v3avU+QCQkhdrBxSCl6pGegmvq/289Jpm038PcqhNAFw2UO2AkOK1A2bWRox9c7NLcVFIU8YOd7sZvuK1UbpX+PVeJecDCIgg1A6AOxIzar3vr6ZomxZuOqBthxpdahlMRfABBEQQagfAeUGoxgn7EXwAARGE2gFwHjNqyAXBBxAQQagdAOcxo4ZcEHwAefBSvYwg1A6A85hRQy7Y7QLkyIvZ/bPrKrXhlsv6tDtieLthrsSMWlO0LWXeR0jx+4sZNfREnQ8gB16vl+Gl2gEwX+LvQVLS34RX/h5gj2z6b5ZdgCz5Ibs/UTtg7oQLNGXscAIP5CUxoxYJJy+tRMKlBB5IiWUXIEvZZPf7sTgQkMrsukrNrI0wo4aMEHwg8LJdgiC7H0jNr9U4YT+CDwRaLkmjZPcDQH7I+UBg5VoSmnoZAJAfgg8EUj5Jo9TLAID8EHwgkPItCU12PwDkLuvgY9euXaqvr1dVVZVCoZCeeuqppK9//vOfVygUSvqYPHmyXe0FbGFH0ujsukrtXj5Djy2YrH+8cYIeWzBZu5fPIPAAgAFknXDa2tqqSy+9VH/2Z3+m66+/PuU1s2fP1sMPP9z9+ZAhQ3JvIZBGPoWy7EoaJbvfTBRRA8yWdfAxZ84czZkzp99rSkpKFIlEcm4UMJB8S5tTEtq/vFj2HgiaguR8PP/88xo5cqTGjx+vBQsWqLm5Oe217e3tisViSR9Af3LdpdITSaP+ZMe9AfiZKYdh2h58zJkzR48++qh27NihBx98UPv379eMGTPU3t6e8vo1a9YoHA53f1RXV9vdJPiInaXNSRr1Fz+UvQcKaduhRk1du0M3bdyrv3r8oG7auFdT1+5wJSjP62C5UCikJ598Utdee23aaxobGzVmzBg9/vjjmjdvXp+vt7e3JwUmsVhM1dXVHCyHlPYca9FNG/cOeN1jCyZnnItBfoA/FOLesBP3GdzkxGGY2RwsV/AKp5WVlRozZoyOHj2a8uslJSUqKSkpdDPgE7nsUhnoTZ+kUX8wuew9eSj+4rVAcqBZwZDis4IzayOOPY+CBx8tLS06ceKEKiv5A0P+st2lwpt+cJha9j7diDORh8ISn7d48T3FxMMws875eOedd3Tw4EEdPHhQktTQ0KCDBw/q+PHjeuedd3THHXdoz549euONN/T888+rvr5eI0aM0HXXXWd32xFA2ZQ2J/kwWEwse08eir949T3FxFnBrIOPl156SRMnTtTEiRMlSUuXLtXEiRN17733qri4WK+++qrmzp2r8ePHa/78+Ro/frz27NmjsrIy2xuP4Ml0l4ok3vQDxsQdTPlW0oU5vBxImjgrmPWyy/Tp09Vfjuqzzz6bV4OAgSR2qfSe+oz0mPrcc6zFuGlGFF4m94aTTBxxIjcmLl1kysS6RgXP+QAKYXZdpWbWRpKSvi4f87/089+8rS0HT+ro/5zJ6Pvwpu8/qe4NtxICTRxxIjdeDiQTs4ILNx1QSEoKQNyaFST4gGf13KWy7VCjrvm7n/U7MkmFN31/MmUHk10jTq/trvAjrweSps0KEnzA89LtJugP5dPhBDtGnF7cXeFHJi5dZMukWcGClFcHnNJfElg6lE/PjSllmb0mn0q6Xt1d4UcmJjTnIjErOHfCBZoydrhr7c2rwmkhZFMhDci0qmVPjBqzx+g7f9kunXR2WZq6dkfapcTESHv38hnGd3h+wt9CekZVOAUKKdPkrts/+VGNG/Uh1stzQJEse2Sbh+Ll3RV+ZtLShZcRfMDT3jjVmtF1n/joCN6gc2BiWeag8PLuCr8zJaHZy8j5gGdtO9Sof/hJ6jODEtyoauknFMlyj9d3VwD9IfiAJyVG5JnwQhKYqRh9u8fEcvGAXQg+4EkDjcgTlvzhePIR8sDo2z1+2V0BpELwAU/KdKT9kRHnFbgl/ha00bdp24nz2aYLmIyEU3gSI3JnmFiWuVBM3ULJ7gr4ETMf8KRCjshNG/263Z4gjL5NL+ZlSmEowC7MfMCTCjUiN230a0p7/Dz6Zjsx4LzAzHy4PXqE/ewekZs2+jWtPX4dfbOdGHBeIGY+TBk9wn52jchNG/2a1h7T5XPqK9uJAef5PvigNLT/2VFt0LRS1qa1x2T5Di5IXgac5+tll4FGj1J89MgSDEwb/ZrWHlPZsTQVtO3EgAl8HXywlotMmTb6Na09JrJrcEExL8B5vg4+GD0iU6aNfk1rj4nsHFwEYTsxYBJf53wwekSmTCumZVp7TGT34MLP24kB0/h65oPRI7Jh2ujXtPaYphCDC79uJwZM4+uZD0aPyJZpo1/T2mOSxOCiKdqWMu8jpHigxuACME/IsiyjtnrEYjGFw2FFo1GVl5fb8j2p8wGT5FOTAskSu12k1IMLZogA52TTfwci+JB4w4cZCITtx2sKmIHgAzBQuoJ3jNLzx+ACcF82/bevcz4AU1AuvbDsqHILwDm+3u0CmIKCdwBwDjMfyAnT3Nmh4B0AnEPwgayR4Jc9Ct4BwDksuyAr+Rzk1dllac+xFm05eFJ7jrUE6kA/Ct4BwDnMfCBj+SRNBn22hIJ36bGEBwQPwQcylk3SZM+dB+m2mCZmS4KyxTRRLr13EBYJQBCWLsAIelAqEXwhmAg+kLFckibZYposiOXS0wUYf3xppb6zqyHnoNQPnTbBl1n8cE95BcEHMpZL0mSusyVeku0bVpBqUqSb9WqMtunbuxpS/p9MglI/dNrMCJrFD/eUl5BwiozlkjTp9y2m2w41auraHbpp41791eMHddPGvZq6dke/ibdB0d+s10D6q3uSLum5MdqmWzcd0DOvmP/aDzQjKMWDryAlZbspn0R65IbgAxlLJE1K6hOApEua9PMWU96w+jfQrFcmegelmQQ0tz92QM+88lZeP7fQKDpnDgJBdxB8ICuJpMlIODlYiIRLU04T+3WLKW9YA7NjNqt3UJpJQNNlSbdtftno4M/vM4JeQiDoDnI+kLVskib9usU0CLks+cpnNiukeEDbOyjNpjMuRCKzXQmJfp4R9BoCQXcQfCAn2SRN+nGLKW9YA0vMejVF2/pdJskmKM2mM7Y7+LMzIXGg1yZd8JUNJ3Zu+GF3CIGgOwg+4Ai/bTHlDWtgmcx6ffHqGj39i8Z+g9KeHdyIYSWKlJeqKeZs8Gf3zpRCzwg6sXPDL7tDnAgE0VfIsiyjFqVjsZjC4bCi0ajKy8vdbg6QUmeXpalrdwz4hrV7+QzPBlh2GaiT6m/0nOr/fvi8wfrdu+9n9LMfWzA575mPxO863TJbPr/rQnTg6QKlRMvs2MLrxM9wUuL5SKkDQa89H7dk038TfAA54g0rc7lMz/fXwVmSQiEp3buXncHfnmMtumnj3gGvyzXQsXPpopCBkpM/ww1+mclxUzb9N8su8BxT1pn9mMtSKNkWVsukMm64dJB+997ZPl+3O5G50Pk9dhadcyIR2q/J1n5bGjYdwQc8xbTRCW9YhZFJB/e7987qy384To/vP1HQ4M9L+T1OJEL7Odk6SNWH3UbwAc8wtRw1b1j2y7Tj+siIYdq9fEZBgz8vJSQ6ESh5KRiDuSgyBk+gqFewZNPBJYK/uRMu0JSxw22fdcqlsq9bnCjq59fCgXAWwQc8gSqEwWJaB5dtZV+3OBEoeSkYg7lYdoEn+HmdGX2ZWBnXK/k9TiRCk2yNfBF8wBNYZw4eEzs4r+T3OBEoeSUYg5kIPuAJXkr6g33o4HLnRKDklWAM5iH4gCeYOA3vN6bUT+mNDg7wH4IPeIaJ0/B+YVr9FAD+Rnl1eI6pI3Sv8ts5HX7APQ4vorw6fI1pePtkUsZ81dYjmlkbofNzSKpZqIphQ3TthCrNrI0QiMAXqPMBBBj1U8ySmIXq/Ts53dqh7/3XG7pp415NXbtD2w41utRCwB4EH0CAUT/FHP3NQvXU+MFxAgQg8DKCDyDAqJ9ijoFmoXrL5jiBzi5Le461aMvBk9pzrIVjCOA6cj6AAKN+ijmymV3K5th6djLBRMx8IJDsGgl6fUTZ3zkdUryT+1RdvMiX156b1+QyuzRQwJIuh6SJpRu4jJkPBI5dI0G/jCjT1U8pCkldlvTd/3pD3/2vNzz53LxkoFmoVPoLWNjJBJMx84FAsWskWOgRpdMzKrPrKrV7+Qw9tmCyvvCJj0iKBx49MVourIFmoXrK5FRfdjLBZAQfCIyBRoJSZkl8dn2fdLYdatTUtTt008a9+qvHD9qyvTKTYKa4KKRJNRX6z0NNKb+HHc8N/UvMQkXC6Wc0Mj1OgJ1MMBnLLgiMbEaC/SXx2fV9UklXbTQx65BLtdFslocK+dyQmZ6H6f3kSJOePHhSp1vf7/56pscJsJMJJst65mPXrl2qr69XVVWVQqGQnnrqqaSvW5allStXqqqqSkOHDtX06dN1+PBhu9oL5MyukWChRpSFmFHJdnmI0bIZElV876m/WPvvnqnHFkzWP944QY8tmKzdy2dkFIAmckjSzY1ksnQDFErWwUdra6suvfRSrV+/PuXXH3jgAa1bt07r16/X/v37FYlENHPmTJ05cybvxgL5sGskWKgRpd1r9LkEM4yWzZMIROZOuEBTxg7PODm0vxwSToKG27IOPubMmaP7779f8+bN6/M1y7L0jW98Q3fffbfmzZunuro6/eAHP9C7776rzZs329JgIFd2jQQLNaK0e9Yhl2CG0bK/pMshiYRLOTAQrrI14bShoUFNTU2aNWtW92MlJSW65ppr9OKLL6b8P+3t7YrFYkkfQCHYNRIs1IjS7lmHXIIZRsv+03MnU7ZLN0Ch2Bp8NDXFs+RHjRqV9PioUaO6v9bbmjVrFA6Huz+qq6vtbBKQxK6RYCFGlHbPOuQazDBa9p9cl26AQinIbpdQKPnGtiyrz2MJd955p5YuXdr9eSwWIwBBQfXcTdB8pk0jy0pzOqbcru+TkJh1WLjpgEJSUq5GLrMO+ZROt/u5AUBPtgYfkUhEUnwGpLLy3Oioubm5z2xIQklJiUpKSuxsBjCgxEjQlO+TkK7aaKbbK3u3LZ9gxu7nBgAJtgYfNTU1ikQi2r59uyZOnChJ6ujo0M6dO7V27Vo7fxTgW3bOOtgZzACAXbIOPt555x299tpr3Z83NDTo4MGDqqio0OjRo7VkyRKtXr1a48aN07hx47R69Wqdd955uvnmm21tOOBnds46sIQCwDRZBx8vvfSSPvnJT3Z/nsjXmD9/vr7//e9r2bJleu+993Tbbbfp7bff1pVXXqnnnntOZWVl9rUaQFZYQgFgkpBlWUYd0hCLxRQOhxWNRlVeXu52cwBHdXZZzFAA8KRs+m/OdgEMkc0ZLEA2Os526Yd73tBvTr+rMRXn6XNTPqIhgzhXFO5h5gMwQLoD5RJzHtTXMJ+ps1ZrnjmijS80qOeRQEUhacG0Gt35qVr3GgbfYeYD8JCBzmAJKX4Gy8zaiBGdmVNM7cxTMXXWas0zR/TtXQ19Hu+y1P14rgGIl34/MA/BB+AyjrHvy9TOPJV0s1aJk4PdmrXqONuljS/0DTx62vhCg/561seyXoLx0u8HZmLRD3np7LK051iLthw8qT3HWrI67h1xHGOfLNGZ9w7IEp35tkONLrWsr96zVsVdnao5fVLSuaJudz95SE8eeNPxv48f7nlDA/24Lit+XTbc/P3wfuMfzHwgZ4x+7MEx9uekWoIqa2/VmZJhRi5Bdc9aWZZmvvbfWv789zX0/XbNWPAttQ8ukSWppbVDX/6XX0hy9u/jN6fftfU6yd0lQt5v/IWZD+TES6NT03GM/TlJS1CWpZsObtN/PfRnqmuKFzbsuQRlguYzbZrw1q/1xOavaOOP7tdHT7+poWfbNf7U8ZTXO/n3MabiPFuvk7JbIrQT7zf+Q/CBrA00+pHiox+mRDPDMfbn9F5amv76SyrveFff3PJ1lbW3pr3OFceO6RN3L9JTP/xrXfnmYbUNGqKHJv8fXfOljXq1clzK/+Lk38fnpnxEA90yRaH4dZlyY4mQ9xt/IvhA1twa/fgZx9jHJS0thUL6m08t0ZvlIzXmd01as2299EFlAFeXoE6dkpYskS66SCOe2aKuUEj/WveH+uSCb+uBaz6vMyXD+v3vTv19DBlUpAXTavq9ZsG0mqySTd1YIuT9xp/I+UDWSJAsDM5gObcE1RRtkyUpVvohLf7jZfqXzcv1mV+9oD2jL9GO6fPcWYJ67z3pn/5JWr1aisXij/3RH2nPgju0bH971t/Oib+PxDZau+p89P799BZSPGC28/fD+40/MfOBrJEgWTiJM1jmTrhAU8YOD1TgIaVegnr5go/pgavnS5Lu/elG/f14y9nXpbNTeuQR6fd/X/rKV+KBx4QJ0nPPSdu26RPX/2HKWauBOPX3ceenavWr++bonk9fpD+dMkb3fPoi/eq+OTnV93BjiZD3G3+iwimy1tllaeraHQOOfnYvnxG4ztNrTC0U1XtnQ8jq0qYtf6tP/Pq/pfHjpZdekpw4rPK556Rly6RfxHerqLpa+tu/lT77WakoeeyWeC2bYm267z8O63Tr+ym/pR/+PpzcecL7jXdk038TfCAniexzSUlvCJQD9w7Tty72CYzKulR82UTp5Ml45//DH0qhAnU2v/hFPOh47rn45+GwdNdd0uLF0tChA/73IPx9OBm4BuH19AOCDzjC9M4L6Xn2LJndu6Xp0+NLId/9rvSFL9j7/U+ckO65J77MYlnS4MHSokXSV78qDc+uuix/H/bi9TQfwQccY+q0PdJLTGOn20Fg/DT2mjXxWYihQ6X9+6WLL87/e0aj8e/7j/8otX3wutx4Y3yJ5fd+L+dvy9+HvXg9zUbwAXhczzfZER8qkSzpVGu7LW+4e4616KaNewe87rEFk808S6arS5ozJ74kUlsr7dsnDet/e2taHR3Shg3SffdJLS3xx66+Wvq7v5MmTbKvzUAAcKot4GGpppd7yneq2fNbF4uK4vkeEyZIR47E8zC+973svodlSf/6r9Kdd0qvvx5/7KKLpLVrpc98pnC5JAAksdUWMEq6MtI95VtS2hdbF0eOlDZvjgciDz8cD0Yy9cIL0uTJ0g03xAOPSET69relV16R6usJPAAHEHwAhuivjHRP+ZaU9s1ZMtOnS/feG//3woXqfPWQDj+6RS+t+aYOP7pFne+fTb7+V7+S5s6NL6sklmpWrZKOHpW++EVpEBPBgFMIPgBDDFRGuqd8Skr76iyZr35V+uQnpdZWdU2YoItvuVZX3HW7Lr7lWp0aUamXv/FdqalJWrhQqquTnn5aKi6Wbr1Veu21ePDyoQ+5/SyAwCH4AAyRS45FrnkZvjlLprhYr17zaVmSBnd1Jn3p/NgpXfrlv1DnmI9I3/pWfHvu3LnSoUPxJNNIxJUmAyDhFAHhhS16ueRYpPo/mT5XP5wl0/n+WY1c93VZ6juL0z2y6miXNWmSQn//99K0ac42EEBKBB/wPa8UJxro0K6e0h3gle1zTZwl41W/+pcf6+LYqQGvO7L4Tl1M4AEYg2UXeEZnl6U9x1q05eBJ7TnWklGyZbrdI/nuGCmE/nIxekqXl+Gl52qX946/mdl1J04WuCUAssHMBzwhl9mL/naPJKbpV209opm1EWOWGhK5GP3V+YikeN5efK52GDr6QluvA+AMgg8YL905JIkRfboEyYF2j/TcMWLS0kPvXIxMKpx69bnm62MTxulsUbEG9Uo2TbAk/bZsuD72fz/tbMMA9IvgA0bLZ0Tv5Uqe2eZiePm55uypp1T8p38qdXXKUvx+6LmOnLg/yj5UquLo76QRI9xoJYAUyPmA0bIZ0ffmi0qeGQrSc9XZs/Gy6NddJ505I02bpkOrHtRvy5ODi1MfqtD75WENbTwpzZol/e537rQXQB8EHzBaPiN631TyzEBgnutvfyvNni19/evxz7/8ZemnP9Ul9y7ViFONOrzpKb20er0Ob3pKFaf/R4P37pHOP196+eX4YXRnzrjbfgCSCD5guHxG9L6q5DmAQDzXffukyy+XfvrTeGn0xx+X1q2TBg+WJBUPHqSLPztXV9y5SBd/dq6KBw+KHxb3k59IFRXS3r3xQ+PefdflJwKA4ANGy3dE75tKnhnw7XO1rPjBb9OmSSdOSOPHS//93/GD4TLxB38gPfecVF4u7doVr3La5qPcF8CDQpZlZX8yVQHFYjGFw2FFo1GVl5e73RwYILHbRVJS4mkiIMmkY/VChVO7+Oq5vveedNtt0ve/H//8uuvi/87lveHFF+O5H62t0qc/Lf3oR9KQIXa2Fgi0bPpvgg94gleqlPqq43dbQ4N0/fXxfI2iImnNGulv/ia/I++ffz6e+9HWFv/ejz/OabaATQg+4Eumd+xeCZA84T//U/rsZ6W3345vkX3iCWnGDHu+93PPSfX1UkeHdPPN0iOPxE+6BZAXgg/AYekKoWWzNARJXV3S/fdLK1fGcz0mTZL+7d+k6mp7f87WrdK8efFtu1/4grRxY3x2BUDOsum/+WsD8jRQITQpXggtk7NoAu3tt+MzEitWxAOPW2+NJ4jaHXhI8Z+zeXM84Pje96TFi+M/E4AjCD6APOVTCA0fOHhQuuIK6ZlnpNJS6eGHpQ0bpJKSwv3MP/kT6Qc/iOeQPPRQPJ+EAARwBJlWQJ4CWdrcTo88In3pS/Ek0Joa6d//XZo40Zmffcst8R01X/yi9OCD0nnnSV/7mjM/20am50MBvRF8AHkKVGnzHKXsHN/viFco3bAhftGcOdKmTfGCYE5asCAe+PzlX0r33ScNHRov3+4RJDrDiwg+gDwlCqE1RdtS5n2EFC/05fnS5jlK1Tn+gc7okf/8O334lQPxZY8VK6R77nEv6XPx4vgMyPLl0l13xQOQJUvcaUsWcj3xGXAbOR9AnjIpbX7Ppy/SvobT2nLwpPYcawlM8mmic+wZeEz5zSv63j/fqg+/ckDvl4Wl//iPePDh9m6TZcviu2yk+IzMt76lzi5Le461GPl7I9EZXsbMB2CDRGnz3iP8SLhUf3xppe778S8DNy3ep3O0LH1p379r2c5HVGx16fDI39O9f7pS/zJ7joypsnHvvfEZkLVrpYULtfonr+u7H72m+8sm/d6ySXSeMna4cw0DMkDwAdhkdl2lZtZGknIb3m7t0KLNwZwW79k5Du1o07ofr9Oc//eiJOnf6v63vjprodqKS83qHEMhac0avfHmKX3k0e/qrh89qObPSFtr4wGISb83Ep3hZSy7ADYqLgppytjhmjvhAk2qqdB9Pw7utHjPTq9j0GCF295RR9Eg3T3rNt3xqSVqG1za5zoTdFrSTbU3avOls1Vsdekf/uNB/dEHQZNJvzcSneFlBB9AgZha/8OpPIaenV5nUbEW//Ey/d/PrtWjEz+VdD6LaZ3jvobTaoy16+4/uk3/XjdDg6wu/fOWB1TX9Jokc+q25HviM+Amll2AAjFxWtzJbZm9dwG1DPuwWoZ9uPvrpu4CSvw+rFCRls35K5WcfV/vDSrRkZE1Ka9zSyLR+dYPTnzuzZK0or6Weh8wEjMfQIGYNi2eaueJdC6PYduhRlt/Xia7gEzsHHvP2Cz5zF9r2af+Ul1FxWmvA5Adgg+gQEyaFndrW2ZiF1AknNxRR8KlRiRtptL793a2eJCs0Lm3SlOWMxK/03RCMiM3BUiFZRegQBIj/4WbDigkJXX8To/83dyWmWoXkMnlv036vfWHrbbwMmY+gAIyZeTvdv5Jz11AU8YOd73jHogpv7f+uP07BfLBzAdQYOlG/pK051iLI7MBpuWfeIHpMzb8TuFlBB+AAxIj/wSnDwPj/Jnc9P69mYTfKbyMZRfAYU7vOpG8u/Ok0Ew+u2Ug/E7hZSHLsoz6a4vFYgqHw4pGoyovL3e7OYCtOrssTV27I22iYGK0unv5jIJ0Ghy/fo5fXgu/PA94Xzb9N8EH4KA9x1p008a9A1732ILJBZvu7+yyjM1jcEq6o+gTr4IpSaWZ4ncKE2TTf5PzATjIhB0KJucxOGGgmieJ+hgzayOe6cCD/juF95DzATiIHQruM/XMHSBICD4AB5lU9TSoTJh9AoKO4ANwEDsU3MfsE+A+gg/AYV6onulnzD4B7iPhFHCB6dUz/cwrZ7cAfsZWWwCBRH0MwF5stQWAATD7hFxRVyV/BB8AAov6GMgWM2b2sD3hdOXKlQqFQkkfkUjE7h8DAICj3DiXya8KMvNx8cUX6yc/+Un358XFxYX4MQAAOMKPlXHdVJDgY9CgQcx2AAB8I5vKuCzlDawgdT6OHj2qqqoq1dTU6MYbb9Trr7+e9tr29nbFYrGkDwAATEJlXHvZHnxceeWVeuSRR/Tss89q48aNampq0lVXXaWWlpaU169Zs0bhcLj7o7q62u4mAQCQFyrj2qvgdT5aW1s1duxYLVu2TEuXLu3z9fb2drW3t3d/HovFVF1dTZ0PADAEW0vjr8HUtTvUFG1LmfcRUrxK8e7lMwL32iQYVedj2LBhuuSSS3T06NGUXy8pKVFJSUmhmwEAWaHDjWNraRyVce1V8OCjvb1dv/zlLzVt2rRC/ygAsMUzrzTqq1sO6XRrR/djQexwE1tLe4/0E1tLg3YWUeJcpt7BWCSA90a+bF92ueOOO1RfX6/Ro0erublZ999/v3bu3KlXX31VY8aMGfD/U14dgJvWPHNE397VkPJrISkwHW5imSHdDo8gLzMwK5aaq8sub775pm666SadOnVK559/viZPnqy9e/dmFHgAgJueeeWttIGHFJ9qD0otB7aWpkdl3PzZHnw8/vjjdn9LACi4zi5LX91yaMDrgtLhsrUUhVSQOh8A4DX7Gk7rdOv7GV0bhA6XraUoJA6WA+Aak9bOswkogtDhTqqpUGW4dMCtpZNqKpxuGnyA4APwGZM69P6YtoUz04Bi+LAhgehw2VqKQiL4AHzEtA49HRO3cCZG+v0lWUrSfXPrAtPhsrUUhVLwCqfZYqstkJt0HXqimzRli6jJWzjTvYYJX7q6Rnd+qtbRNpnAK7NpcFc2/TcJp4APDHTctxTfItrZ5f5YI5stnE5LjPQrw8lLMBXDBuuhmycGMvCQzm0tnTvhAk0ZO5zAA3lj2QXwAS/VZDB9C+fsukrNrI0w0u+F2Q/YieAD8AHTO/SevLCFkyJSybySSwTvYNkF8AEvdOgJicTOdGPmkOIdWxB2lHhBIg+m98xaIjl426FGl1oGLyP4AHzASx16YgunpD7tZQunWbyUSwRvIfgAfMBrHXoisTPSK7EzEi41ZlcOzE4OhreR8wH4hNdqMpDYaT4v5RLBWwg+AB/xWodOYqfZvJRLBG8h+AB8hg4dduF8FxQKOR8AgJS8lksE7yD4AACkRXJwap1dlvYca9GWgye151gLO36yxLILAKBfXsslKjSKruWPg+UASKJ8NpAJrxzg6IZs+m9mPgAwkgMyMFDRtZDiRddm1kYI3AdAzgcQcJTPBjJD0TX7EHwAHwhiAhnls4HMUXTNPiy7AAruskM2IzlqhyDoKLpmH2Y+EHhBXnZgJAdkzksHOJqO4AOBFvRlB0ZyQOZMLrrmtWVjll0QaEFfdqB8NpAdEw9w9OKyMcEHAi3oyw6JkdzCTQcUkpICELdHcoCpTCq6lq7uSGLZ2NS6Iyy7INBYdqB8NpCLxAGOcydcoCljh7u21OLVZWNmPhBoLDvEmTSSA5AZLy8bE3wg0Fh2OCcxkgPgDV5eNmbZBYHX37LDN2+eqPDQIZ7JIIczvLazAP7k5WVjZj4ApV52eLu1Q/f92FsZ5Cg8L+4sgD95edmYmQ/gAz0TyKLvdWjR5mAWHkN6QS5IB/OYXHdkIAQfQC9eziBH4XBfwERe3a3GsgvQi5czyFE43BcwlRd3qxF8AL14OYMchWPHfdHZZXmqg4B3eG23GsEH0IuXM8hROPneFySqAueQ8wH0wsmVSCWf+4JEVSAZwQfQi5czyFE4ud4XJKoCfRF8ACl4NYMchZXLfZFNoioQFOR8AGl4MYMchZftfUECM9AXwQfQD69lkMMZ2dwXJDADfbHsAgAFRAIz0BfBBwAUEAnMQF8EH0CBcQIqSGAGkpHzARQQhaWQQAIzcE7IsiyjhmGxWEzhcFjRaFTl5eVuNwfIWaKwVO8/sERXw4gXgJ9k03+z7AIUAIWlACA9gg+gACgsBQDpEXwABUBhKQBIj+ADKAAKSwFAegQfQAFQWAoA0iP4AAqAwlIAkB7BB1AgFJYCgNQoMgYUEIWlAKAvgg+gwDgZFwCSsewCAAAcRfABAAAcRfABAAAcRfABAAAcRfABAAAcRfABAAAcRfABAAAcRZ0PwAadXRaFxAAgQwQfQJ62HWrUqq1H1Bht636sMlyqFfW1lFAHgBQKtuzy0EMPqaamRqWlpbr88sv1wgsvFOpHAa7ZdqhRCzcdSAo8JKkp2qaFmw5o26FGl1oGAOYqSPDxxBNPaMmSJbr77rv18ssva9q0aZozZ46OHz9eiB8HuKKzy9KqrUdkpfha4rFVW4+osyvVFQAQXAUJPtatW6c///M/11/8xV/ooosu0je+8Q1VV1drw4YNhfhxgCv2NZzuM+PRkyWpMdqmfQ2nnWsUAHiA7cFHR0eHfv7zn2vWrFlJj8+aNUsvvvii3T8OcE3zmfSBRy7XAUBQ2J5weurUKXV2dmrUqFFJj48aNUpNTU19rm9vb1d7e3v357FYzO4mAQUxsqzU1usAICgKlnAaCiVvM7Qsq89jkrRmzRqFw+Huj+rq6kI1CbDVpJoKVYZLlW5DbUjxXS+TaiqcbBYAGM/24GPEiBEqLi7uM8vR3NzcZzZEku68805Fo9HujxMnTtjdJKAgiotCWlFfK0l9ApDE5yvqa6n3AQC92B58DBkyRJdffrm2b9+e9Pj27dt11VVX9bm+pKRE5eXlSR+AV8yuq9SGWy5TJJy8tBIJl2rDLZdR5wMAUihIkbGlS5fqc5/7nK644gpNmTJF3/nOd3T8+HHdeuuthfhxgKtm11VqZm2ECqcAkKGCBB833HCDWlpa9LWvfU2NjY2qq6vTM888ozFjxhTixwGuKy4KacrY4W43AwA8IWRZllEVkGKxmMLhsKLRKEswAAB4RDb9N6faAgAARxF8AAAARxF8AAAARxF8AAAARxF8AAAARxVkqy0Ae3R2WdQPAeA7BB+AobYdatSqrUfUGD13Km5luFQr6mupnArA01h2AQy07VCjFm46kBR4SFJTtE0LNx3QtkONLrUMAPJH8AEYprPL0qqtR5Sq+l/isVVbj6izy6j6gACQMYIPwDD7Gk73mfHoyZLUGG3TvobTzjUKAGxE8AEYpvlM+sAjl+sAwDQEH4BhRpaV2nodAJiG4AMwzKSaClWGS5VuQ21I8V0vk2oqnGwWANiG4AMwTHFRSCvqayWpTwCS+HxFfS31PgB4FsEHYKDZdZXacMtlioSTl1Yi4VJtuOUy6nwA8DSKjAGGml1XqZm1ESqcAvAdgg/AYMVFIU0ZO9ztZgCArVh2AQAAjiL4AAAAjiL4AAAAjiL4AAAAjiL4AAAAjiL4AAAAjiL4AAAAjiL4AAAAjiL4AAAAjjKuwqllWZKkWCzmcksAAECmEv12oh/vj3HBx5kzZyRJ1dXVLrcEAABk68yZMwqHw/1eE7IyCVEc1NXVpbfeektlZWUKhThAK51YLKbq6mqdOHFC5eXlbjfH13itncNr7Sxeb+cE4bW2LEtnzpxRVVWVior6z+owbuajqKhIF154odvN8Izy8nLf3sim4bV2Dq+1s3i9neP313qgGY8EEk4BAICjCD4AAICjCD48qqSkRCtWrFBJSYnbTfE9Xmvn8Fo7i9fbObzWyYxLOAUAAP7GzAcAAHAUwQcAAHAUwQcAAHAUwQcAAHAUwYcHPfTQQ6qpqVFpaakuv/xyvfDCC243yZdWrlypUCiU9BGJRNxuli/s2rVL9fX1qqqqUigU0lNPPZX0dcuytHLlSlVVVWno0KGaPn26Dh8+7E5jPW6g1/rzn/98n/t88uTJ7jTW49asWaOPf/zjKisr08iRI3Xttdfq17/+ddI13NtxBB8e88QTT2jJkiW6++679fLLL2vatGmaM2eOjh8/7nbTfOniiy9WY2Nj98err77qdpN8obW1VZdeeqnWr1+f8usPPPCA1q1bp/Xr12v//v2KRCKaOXNm99lPyNxAr7UkzZ49O+k+f+aZZxxsoX/s3LlTixYt0t69e7V9+3adPXtWs2bNUmtra/c13NsfsOApkyZNsm699dakxz72sY9ZX/nKV1xqkX+tWLHCuvTSS91uhu9Jsp588snuz7u6uqxIJGJ9/etf736sra3NCofD1re+9S0XWugfvV9ry7Ks+fPnW3PnznWlPX7X3NxsSbJ27txpWRb3dk/MfHhIR0eHfv7zn2vWrFlJj8+aNUsvvviiS63yt6NHj6qqqko1NTW68cYb9frrr7vdJN9raGhQU1NT0n1eUlKia665hvu8QJ5//nmNHDlS48eP14IFC9Tc3Ox2k3whGo1KkioqKiRxb/dE8OEhp06dUmdnp0aNGpX0+KhRo9TU1ORSq/zryiuv1COPPKJnn31WGzduVFNTk6666iq1tLS43TRfS9zL3OfOmDNnjh599FHt2LFDDz74oPbv368ZM2aovb3d7aZ5mmVZWrp0qaZOnaq6ujpJ3Ns9GXeqLQYWCoWSPrcsq89jyN+cOXO6/33JJZdoypQpGjt2rH7wgx9o6dKlLrYsGLjPnXHDDTd0/7uurk5XXHGFxowZox//+MeaN2+eiy3ztttvv12vvPKKdu/e3edr3NvMfHjKiBEjVFxc3CdCbm5u7hNJw37Dhg3TJZdcoqNHj7rdFF9L7CjiPndHZWWlxowZw32eh8WLF+vpp5/Wz372M1144YXdj3Nvn0Pw4SFDhgzR5Zdfru3btyc9vn37dl111VUutSo42tvb9ctf/lKVlZVuN8XXampqFIlEku7zjo4O7dy5k/vcAS0tLTpx4gT3eQ4sy9Ltt9+uH/3oR9qxY4dqamqSvs69fQ7LLh6zdOlSfe5zn9MVV1yhKVOm6Dvf+Y6OHz+uW2+91e2m+c4dd9yh+vp6jR49Ws3Nzbr//vsVi8U0f/58t5vmee+8845ee+217s8bGhp08OBBVVRUaPTo0VqyZIlWr16tcePGady4cVq9erXOO+883XzzzS622pv6e60rKiq0cuVKXX/99aqsrNQbb7yhu+66SyNGjNB1113nYqu9adGiRdq8ebO2bNmisrKy7hmOcDisoUOHKhQKcW8nuLrXBjn55je/aY0ZM8YaMmSIddlll3Vv44K9brjhBquystIaPHiwVVVVZc2bN886fPiw283yhZ/97GeWpD4f8+fPtywrviVxxYoVViQSsUpKSqyrr77aevXVV91ttEf191q/++671qxZs6zzzz/fGjx4sDV69Ghr/vz51vHjx91utielep0lWQ8//HD3NdzbcSHLsiznQx4AABBU5HwAAABHEXwAAABHEXwAAABHEXwAAABHEXwAAABHEXwAAABHEXwAAABHEXwAAABHEXwAAABHEXwAAABHEXwAAABHEXwAAABH/X9yILegpLP5HAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "focal_points=32\n",
    "\n",
    "plt.scatter(xy[:, 0], xy[:, 1])\n",
    "plt.scatter(xy[focal_points, 0], xy[focal_points, 1], color='red')\n",
    "\n",
    "for i in indices[focal_points]:\n",
    "    plt.plot([xy[focal_points, 0], xy[i, 0]], [xy[focal_points, 1], xy[i, 1]], 'r-')\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "physci200",
   "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}