figure3-mmd-stats.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7e78ff44",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.spatial.distance import cdist, pdist\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "import os\n",
    "from matplotlib.pyplot import *\n",
    "\n",
    "def mmd(X, Y, sigma_sq=None):\n",
    "    \"\"\"\n",
    "    Implements Gretton's test for equality of\n",
    "    distributions in high-dimensional settings\n",
    "    using concentration bounds on the maximum\n",
    "    mean discrepancy (MMD). This function uses\n",
    "    the unbiased estimator of the MMD (see\n",
    "    Lemma 6, Gretton et al., 2012) and upper\n",
    "    bounds the p-value using a Hoeffding\n",
    "    large-deviation bound (see Theorem 10,\n",
    "    Gretton et al., 2012).\n",
    "\n",
    "    The test considers two sets of observed\n",
    "    datapoints, X and Y, which are assumed to\n",
    "    be drawn i.i.d. from underlying probability\n",
    "    distributions P and Q. The null hypothesis\n",
    "    is that P = Q.\n",
    "\n",
    "    Note that this function assumes that the number\n",
    "    of samples from each distribution are equal.\n",
    "\n",
    "    Reference\n",
    "    ---------\n",
    "    Gretton et al. (2012). A Kernel Two-Sample Test.\n",
    "    Journal of Machine Learning Research 13: 723-773.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    X : ndarray (num_samples x num_features)\n",
    "        First set of observed samples, assumed to be\n",
    "        drawn from some unknown distribution P.\n",
    "\n",
    "    Y : ndarray (num_samples x num_features)\n",
    "        Second set of observed samples, assumed to be\n",
    "        drawn from some unknown distribution Q.\n",
    "    \n",
    "    sigma_sq: float, optional\n",
    "        Kernel bandwidth parameter. If not provided,\n",
    "        it is set to the median distance between samples.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    mmd : float\n",
    "        MMD Distance\n",
    "    \"\"\"\n",
    "\n",
    "    assert X.shape == Y.shape\n",
    "    m = X.shape[0]\n",
    "\n",
    "    # Compute pairwise distances\n",
    "    xd = pdist(X, metric=\"euclidean\")\n",
    "    yd = pdist(Y, metric=\"euclidean\")\n",
    "    xyd = cdist(X, Y, metric=\"euclidean\").ravel()\n",
    "    \n",
    "    if sigma_sq is None:\n",
    "        # Set kernel bandwidth (Gretton et al. suggest to use\n",
    "        # the median distance).\n",
    "        sigma_sq = np.median(\n",
    "            np.concatenate((xd, yd, xyd))\n",
    "        ) ** 2\n",
    "\n",
    "    # Compute unbiased MMD distance.\n",
    "    kxx = np.mean(np.exp(-(xd**2) / (2 * sigma_sq)))\n",
    "    kyy = np.mean(np.exp(-(yd**2) / (2 * sigma_sq)))\n",
    "    kxy = np.mean(np.exp(-(xyd**2) / (2 * sigma_sq)))\n",
    "    mmd_obs = kxx + kyy - 2 * kxy\n",
    "\n",
    "    return mmd_obs\n",
    "\n",
    "def get_mmd_sigma_sq(X, Y):\n",
    "    assert X.shape == Y.shape\n",
    "    m = X.shape[0]\n",
    "    # Compute pairwise distances\n",
    "    xd = pdist(X, metric=\"euclidean\")\n",
    "    yd = pdist(Y, metric=\"euclidean\")\n",
    "    xyd = cdist(X, Y, metric=\"euclidean\").ravel()\n",
    "\n",
    "    # Set kernel bandwidth (Gretton et al. suggest to use\n",
    "    # the median distance).\n",
    "    sigma_sq = np.median(\n",
    "        np.concatenate((xd, yd, xyd))\n",
    "    ) ** 2\n",
    "    \n",
    "    return sigma_sq\n",
    "\n",
    "def permutation_test_mmd(X, Y, n_samples=1000, n_iters=1000):\n",
    "    \n",
    "    mmds = []\n",
    "    for i in range(n_iters):\n",
    "        #generate random indicies\n",
    "        x_rand_idx = np.random.choice(np.arange(len(X)), size=n_samples, replace=False)\n",
    "        y_rand_idx = np.random.choice(np.arange(len(Y)), size=n_samples, replace=False)\n",
    "        \n",
    "        #randomly sample data\n",
    "        x_rand = X[x_rand_idx]\n",
    "        y_rand = Y[y_rand_idx]    \n",
    "        \n",
    "        working_mmd = mmd(x_rand, y_rand, sigma_sq=13.65)\n",
    "        mmds.append(working_mmd)\n",
    "    \n",
    "    return np.array(mmds)\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ecd0e88d",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style('white')\n",
    "sns.set_style('ticks')\n",
    "sns.set_context('talk')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "832ed295",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_feather('vocalization_df.feather')\n",
    "\n",
    "df_c2 = df[df.cohort=='c2']\n",
    "df_c4 = df[df.cohort=='c4']\n",
    "df_c5 = df[df.cohort=='c5']\n",
    "\n",
    "c2_latents = np.vstack(df_c2.latent_means.values)\n",
    "c4_latents = np.vstack(df_c4.latent_means.values)\n",
    "c5_latents = np.vstack(df_c5.latent_means.values)\n",
    "\n",
    "all_latents = np.concatenate((c2_latents, c4_latents, c5_latents))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "986e43a8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e55f881b",
   "metadata": {},
   "outputs": [],
   "source": [
    "shuffle_mmd = permutation_test_mmd(all_latents, all_latents, n_iters = 1000)\n",
    "c2_c4_mmd = permutation_test_mmd(c2_latents, c4_latents, n_iters = 1000)\n",
    "c2_c5_mmd = permutation_test_mmd(c2_latents, c5_latents, n_iters = 1000)\n",
    "c4_c5_mmd = permutation_test_mmd(c4_latents, c5_latents, n_iters = 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "16494a6e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgoAAAEZCAYAAAD2aw39AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABS00lEQVR4nO3deVhUZfsH8O+AAgqoDCogJbkOMSC4LygIooJKpli4545LljtIaS6IorhCppiaGv4ylVzDBUVLLRXN0gS0RNEETHFhWIZtfn/wMjkMwzrDgHw/1zVXzjnPOec+z/vquec5zyKQyWQyEBERERVDR9sBEBERUfXFRIGIiIhUYqJAREREKjFRICIiIpWYKBAREZFKTBSKcHFxgYuLi7bDICIiqhbqaDuA6iYtLU3bIRAREVUbbFEgIiIilZgoEBERkUpMFIiIiEglJgpERESkEhMFIiIiUomJAhEREanE4ZFEVKPIZDI8ffoUWVlZyM/P13Y4RDWWjo4ODAwM0LhxYwgEApXlmCgQUY0hk8nwzz//IC0tDfr6+tDV1dV2SEQ1Vk5ODiQSCaRSKSwtLVUmC0wUKmHp0T9x+/Er2DRrgC88xdoOh+iN9/TpU6SlpcHMzAxCoVDb4RDVeKmpqUhJScHTp0/RpEmTYsswUaiE249f4XJCqrbDIKo1srKyoK+vzySBSE2EQiFevHiBrKwslWXYmZGIaoz8/Hy+biBSM11d3RL7+zBRICIiIpWYKBAREZFK7KNARDVeYcdibWPHZnoTsUWBiGq8249f4XaSdhOF20mvKp2sxMfHY/bs2XB0dIStrS169uyJWbNmIS4uTl5mzJgxGDduXCWjLbBnzx707NkT7dq1Q1hYGO7evYuhQ4fC1tYWnp6eCAkJgY2NjVquRTUXWxSI6I1gY9EA+3y6a+363lt/qdTxcXFxGDFiBDp06IBFixZBKBQiOTkZu3fvxocffojdu3fDwcFBPcECyMjIwMqVK+Hs7IwJEybg7bffRlBQEB49eoQvv/wSpqamiI6OVtv1qOZiokBEVA3s2rULpqamCAsLUxjZ0adPH3h4eGDz5s0ICwtT2/XS0tKQl5cHNzc3dO7cGQDw4sULtG3bFs7OzgDARIEA8NUDEVG18OzZM8hkMqVhaoaGhvD394eHh4d8m0wmw9atW9G7d2+0a9cOw4cPx61bt+T7/fz80LdvX4XzPHr0CCKRCIcPH0ZERAScnJwAAP7+/hCJRBCJRLh06RKuXr0KkUiEiIiIYuM8ffo0hg4dCjs7O/Ts2RNBQUHIzs5WVzVQNcREgYioGnBycsKjR48wfPhwhIeH4++//5bvc3d3x5AhQ+Tfr1y5grNnz2LRokVYvXo1UlJSMG3aNOTl5ZXpWr1798ZXX30FAJg2bRr27duHffv2wc7ODjY2Nti3bx969+6tdNzRo0fx8ccfo02bNvjyyy8xdepU7Nu3D3Pnzq3czVO1xlcPRETVwKhRo/Dvv/9i586dWLZsGYCCWfN69uyJMWPGoF27dvKyBgYG2LZtGxo0aACg4DXC559/joSEBLRu3brUawmFQnknxebNm8v7PhgZGSEvL6/YvhAymQzBwcFwcXFBUFCQfLu5uTlmzJiBa9euoWPHjhW9farG2KJARFQNCAQCzJ49GxcuXMC6deswbNgwGBoa4siRI/jwww8RHh4uL9u2bVt5kgAAb731FgDg1SvNjfy4d+8ekpOT4erqitzcXPmnV69eqFu3Li5duqSxa5N2sUWBiKgaadCgAQYOHIiBAwcCAG7fvo0FCxYgKCgIgwYNAgDUq1dP4RgdnYLffJpcdvvFixcAgEWLFmHRokVK+588eaKxa5N2MVEgItKy5ORkDBs2DJ9++ik++OADhX02NjaYNWsWZsyYgUePHpXpfAKBQKm/QkZGRqViNDY2BgAsXLiw2FcMJiYmlTo/VV989UBEpGVNmjSBrq4u9u7dC6lUqrT/3r17MDAwQPPmzct0PkNDQ6SmpiqMRrh27VqlYmzVqhWEQiH++ecf2NnZyT8mJiYIDg5W6HxJbxa2KBARaZmuri4WL16MmTNnwsvLC6NGjUKrVq2QmZmJixcvIjw8HHPmzJH/qi+Ni4sL9uzZA39/fwwbNgx37tzBzp07K7Xypq6uLmbNmoWlS5dCR0cHTk5OePnyJTZt2oS0tDTO4PgGY6JARG+E20mvKj07YmWvb2PRoPSCKvTp0wfff/89tm/fji1btuDZs2fQ19eHjY0NNmzYoDQvQkkcHR3h6+uLPXv24OTJkxCLxQgNDcXw4cMrHB8AeHt7w8jICF9//TX27t0LIyMjdO7cGXPmzEGTJk0qdW6qvgQymUym7SCqk06dOgEAYmJiSi3rvfUXXE5IRdcWQq1OHUtUWzx48AAAYGVlpbCdi0IRVZyqv1eFtNqikJ+fj3379mHv3r149OgRTE1N0adPH8ycORNGRkYAgAsXLmD9+vX466+/YGpqitGjR2PChAkK57l58yZWr16NW7duwdDQEEOHDsXMmTNRt25dbdwWEVUxPpyJNEericLXX3+NDRs2YOLEiejevTsSEhKwadMm/PXXX9i+fTuuX7+OqVOnwsPDA59++imuXbuG1atXQyaTYeLEiQAKMqFx48ahffv22LBhA/7++2+sX78eEokEixcv1ubtERER1XhaSxRkMhm+/vpreHt7y6f/7NGjB0xMTDB79mzExsZi06ZNsLGxwZo1awAUTHGam5uLLVu2YMyYMdDT00NYWBiMjY2xefNm6OnpwdnZGQYGBggICICPjw/MzMy0dYtEREQ1ntaGR6anp+O9996TTyBSqGXLlgCAu3fvIiYmBv369VPY379/f7x69QrXr18HAFy8eBEuLi7Q09OTl3F3d0deXh4uXLig4bsgIiJ6s2mtRcHIyAiff/650vaoqCgABZOM5OTkoEWLFgr7CztbJCQkwN7eHklJSUplhEIhjIyMkJCQoHT+ws6KqqSlpZV5CBIREdGbrlpNuPT7778jLCwMbm5uSEtLAwB5p8ZChoaGAACJRKKyTGE5iUSi4YiJiIjebNVmHoVr165h6tSpeOuttxAQECBvDRAIBMWW19HRQeHIzuLKyGQy+fznrytt2GNpLQ5ERES1SbVoUfjxxx8xfvx4WFhY4JtvvoGJiYm8+b9oq0Dhd2NjY3lLQnEtBxkZGXyFQEREVElaTxR27tyJOXPmwMHBAeHh4WjatCmAgjXSdXV1kZiYqFC+8HuLFi1gaGgIMzMz+WQRhZ49ewaJRKLUd4GIiIjKR6uvHvbv349Vq1ZhwIABCAoKUhi5oK+vj06dOuHUqVP46KOP5K8XTp48CWNjY9ja2gIomKo0OjoaCxYskB9/8uRJ6OrqokuXLlV/U0RU9SL9gOSb2o4CMLcDPFZpOwoitdJai8KzZ8+wYsUKWFpaYtSoUbh9+zZu3Lgh/6SmpmLatGm4fv06Zs+ejfPnz2PDhg3Yvn07fHx85OuxT5o0Cf/++y+mTJmC6Oho7Ny5EytXrsSHH36IZs2aaev2iKgqJd/UfqJQiRj8/PwgEolUfi5duqTmYJW5urris88+AwBcvnwZIpGoTFPZl0V6ejpcXV1x+PBhtZxPHUqqbycnJ41fv2gd+/n5lWs9j+JkZ2djy5YtcHd3h4ODA9zd3bF582aFVUQrQmstCj///DMyMzPxzz//YNSoUUr7V69ejcGDByMkJASbNm3CjBkzYGZmhgULFihM4dyqVSvs2LEDq1evxieffAITExOMHz8eM2fOrMrbISJtM7cDxh/X3vV3DqzU4ebm5ti4cWOx+1q3bl2pc5dFaGioRvp1SSQSTJ8+Hf/884/az11Z3t7eGDp0qNL211u3NUUsFmPfvn1q/d82ICAAR44cwfTp02FnZ4ebN2/iyy+/RHJyMpYtW1bh82otUXj//ffx/vvvl1qub9++pWZZnTp1wvfff6+myIiIqp6enh4cHBy0dn1NLBN94cIFLF++HC9evFD7udXB3Nxca3VuZGSk1mu/fPkS33//PebNm4dJkyYBALp3L1iscO3atZg7dy4aNmxYoXNrvTMjERGVTUZGBtasWYN+/frB1tYWHTp0wMSJExEXFycv4+fnBx8fH4SHh8PV1RX29vaYOHEi/v33Xxw4cABubm5o3749xo0bh0ePHsmPe/3Vw+vu3r0LkUiEgwcPKmy/d+8eRCIRzp07pzLeSZMmwdbWFtu2bSvT/S1cuBC9evVCfn6+wnZ/f3+4urpCJpMhNTUVc+fOhaOjI9q1a4fBgwfj0KFDZTp/RcTGxmLGjBno1q0bxGIxnJycsGLFCkilUnkZkUiEffv2Yd68eWjfvj26deuG0NBQSCQSLFy4EB07doSjoyPWrFkjH9Zf0uudoKAgODg4ID09XWF7cHAwevbsiby8PKVjJBIJhg8fDldXV4XthbMdP3z4sMJ1wESBiKiayM3NVfoUPlgAYMGCBTh06BB8fHywY8cOLFy4EPHx8Zg3b55CuatXr+LgwYNYtGgRFi9ejCtXrmDMmDHYs2cP/Pz88Nlnn+H3339HQEBAqTG1adMGdnZ2Sv0LDh06hCZNmqBXr14qjz1y5AjWrl0LoVBYpvsfPHgwnjx5ovDwzM7ORlRUFDw9PSEQCDB//nz8/fffWLp0KcLCwmBjYwNfX19cvny5TNd4XX5+frF1XiglJQWjRo2CVCpFUFAQtm3bhgEDBmD37t3YvXu3wrmCgoJgYmKCzZs3w8XFBSEhIRg2bBjq1auHjRs3wsXFBV9//TVOnTpValxeXl7IzMzE6dOnFWI9evQoPD09oaurq3SMpaUllixZIk8MCp05cwZ169ZVuYR0WVSbCZeIiGqzxMREiMXKy2UvWbIEI0aMgFQqRWZmJhYtWgR3d3cAQJcuXSCRSLBq1So8f/5c/kBOT0/Hxo0b8fbbbwMATp8+jejoaERFRcm3xcbG4tixY2WKzcvLC0uXLkVSUhIsLCyQn5+PI0eOqHxoFWrbtm256qBr166wsLDA8ePH5aPWLly4gJcvX+K9994DAFy5cgUzZsyAm5ubvA4aNWqEunXrlutaABASEoKQkBCl7X/88Qf09fURHx8PGxsbbNy4UT4rcI8ePXDx4kVcvXoVkydPlh8jFovlLTLW1taIiIiAqampfBVjR0dHREZG4saNG+jfv3+JcbVu3Rr29vY4fPiw/BX9r7/+iuTk5GL7VKhy+vRp/PDDDxg7dmyl+p8wUSAiqgbMzc0RGhqqtN3S0hJAwZDx7du3Ayj4pZuQkID79+8jOjoaAJCTkyM/xtTUVJ4QFH4XCoUK2xo1aiSfBr80gwYNwqpVq3D06FFMmTIFly9fRlJSUrkeWmUhEAjg6emJAwcOYNGiRahTpw6OHz8OsViMVq1aAShIJkJCQnD79m306tULzs7O8PX1rdD1RowYAS8vL6XthZ0ZnZyc4OTkhJycHPz111948OAB7ty5g9TUVDRu3FjhmHbt2sn/bGJiAl1dXdjb2yvcW8OGDfHq1asyxebl5YUlS5YgJSUFZmZm+OGHH2Bra4s2bdqU6fhTp05h7ty56Nixo3yF5opiokBEVA3o6enBzs6uxDI///wzAgMDce/ePRgaGsLa2hr169cHAIVXD4W/fl9XOKS8IoyNjeHm5oYjR45gypQpOHToEOzs7Mr80CqPwYMHIywsDL/88gs6d+6Ms2fPYtasWfL969evx5YtWxAZGYmTJ09CR0cHPXr0wLJly+RJVVk1bdq0xDrPz8/HunXrEB4ejoyMDFhYWKBdu3bQ19dXqG9A/XU+cOBABAYG4tixYxg+fDiioqIwb968Mh37zTffICgoCF26dMGXX34JfX39CscBsI8CEVGNkJiYiBkzZsDGxgZRUVG4fv069u7dCxcXlyq5/rBhw3D37l3ExsYiKioKQ4YM0ch1WrduDbFYjBMnTiA6OhpSqRQDB/439NTY2Bjz58/H2bNnERkZiTlz5uD69etYvny52mMJCwvDN998g0WLFiEmJgbnzp3Dpk2bytznojKMjIzQr18/REZG4vz588jNzcWgQYNKPS4wMBArV67EgAEDsG3btmIXTSwvJgpERDXArVu3IJVKMXXqVIVXCD///DMAKI0UULdu3brB0tISgYGByM7OLtNDq6Lee+89nD9/HpGRkXB0dJQ38ycnJ8PZ2RknTpwAUNCjf/LkyejRoweSkpLUHse1a9cgEokwdOhQ+Tv+lJQU3LlzR+P1DRS8frh16xb27t0LV1fXUoc3btiwAbt27cL48eMRHBystvkg+OqBiKgGEIvFqFOnDtasWYNx48ZBKpUiIiJCPjwxMzNTo9cXCAQYMmQIQkND4e7uXuEx+WUxaNAgrF69GlFRUQgKCpJvNzc3h6WlJQICAiCRSNC8eXPcunUL58+fx/Tp0+Xlbty4AaFQiObNm1cqjnbt2mHz5s3Ytm0b7O3t8eDBA2zduhXZ2dkar2+goD+GpaUlrl69irCwsBLLxsfHY+vWrbC1tYW7uzt+//13hf2tW7eucOsCEwUiejMk36z07IiVvr55yX0MKsPKygpr165FaGgopk6dioYNG8LBwQF79uzBmDFjEBMTI+/wpym9e/dGaGio2jsxFtW4cWM4OjoiJiZGPrqhUEhICIKDg7Fx40Y8f/4cFhYWmDlzpsIIBG9vbwwZMgSrVlVu3Q0fHx88f/4cu3btQlpaGiwsLDB48GAIBAKEhYVBIpGopWlfFYFAAGdnZ5w6dQo9e/Yssezp06eRn5+PW7duwdvbW2l/eHg4OnXqVLE4ZEV7ZNRyhRVZljnOvbf+gssJqejaQoh9Pt01HRpRrVe4UqzSmHAuClUlQkNDsX//fkRHR0NHh2+uNS0/Px/u7u7w8PDA7NmzNXYdlX+v/octCkRU873BD+fqICIiAnfu3EF4eDjmzp3LJEHD0tLSsGvXLvz2229ISUnByJEjtRoPEwUiIipRXFwcvv/+e3h4eGD06NHaDueNZ2BggO+++w4ymQwrV66EmZmZVuNhokBERCXy9/eHv7+/tsOoNerWrYsLFy5oOww5th8RERGRSkwUiIiISCUmCkRERKQSEwUiIiJSiYkCERERqcREgYiIiFRiokBEREQqcR4FIqrxgq4EIS41TtthwFpoDd8uvtoOg0it2KJARDVeXGoc4lPjtRpDfGp8hZMVPz8/iEQilZ9Lly6pOVplrq6u+OyzzwAAly9fhkgkKtOaN2WRnp4OV1dXHD58WC3nU4eS6tvJyUnj1y9ax35+fujbt2+lzpmdnY0NGzbA1dUV9vb28Pb2VsvETWxRIKI3gkgowk73nVq7/vgT4yt1vLm5OTZu3FjsvtatW1fq3GURGhoKY2NjtZ9XIpFg+vTp+Oeff9R+7sry9vYudiVMPT09jV9bLBZj3759av3fdvny5Th27Bjmzp2Lli1b4tChQ/Dx8cGePXvQoUOHCp+XiQIRUTWgp6cHBwcHrV3fxsZG7ee8cOECli9fjhcvXqj93Opgbm6utTo3MjJS67WfPXuGgwcPYt68efL1OLp164Zr167h//7v/yqVKPDVAxFRDZGRkYE1a9agX79+sLW1RYcOHTBx4kTExf33ysPPzw8+Pj4IDw+XN0FPnDgR//77Lw4cOAA3Nze0b98e48aNw6NHj+THvf7q4XV3796FSCTCwYMHFbbfu3cPIpEI586dUxnvpEmTYGtri23btpXp/hYuXIhevXohPz9fYbu/vz9cXV0hk8mQmpqKuXPnwtHREe3atcPgwYNx6NChMp2/ImJjYzFjxgx069YNYrEYTk5OWLFiBaRSqbyMSCTCvn37MG/ePLRv3x7dunVDaGgoJBIJFi5ciI4dO8LR0RFr1qyBTCYDUPLrnaCgIDg4OCA9PV1he3BwMHr27Im8vDylYxo0aID9+/fDy8tLvk1HRwd16tRRiLUimCgQEVUTubm5Sp/CBwsALFiwQN6cvGPHDixcuBDx8fGYN2+eQrmrV6/i4MGDWLRoERYvXowrV65gzJgx2LNnD/z8/PDZZ5/h999/R0BAQKkxtWnTBnZ2dkr9Cw4dOoQmTZqgV69eKo89cuQI1q5dC6FQWKb7Hzx4MJ48eaLw8MzOzkZUVBQ8PT0hEAgwf/58/P3331i6dCnCwsJgY2MDX19fXL58uUzXeF1+fn6xdV4oJSUFo0aNglQqRVBQELZt24YBAwZg9+7d2L17t8K5goKCYGJigs2bN8PFxQUhISEYNmwY6tWrh40bN8LFxQVff/01Tp06VWpcXl5eyMzMxOnTpxViPXr0KDw9PaGrq6t0TN26dSEWi9GwYUPk5+cjKSkJK1asQGJiIj744INy183r+OqBiKgaSExMhFgsVtq+ZMkSjBgxAlKpFJmZmVi0aBHc3d0BAF26dIFEIsGqVavw/Plz+QM5PT0dGzduxNtvvw0AOH36NKKjoxEVFSXfFhsbi2PHjpUpNi8vLyxduhRJSUmwsLBAfn4+jhw5ovKhVaht27blqoOuXbvCwsICx48fR5cuXQAUvL54+fIl3nvvPQDAlStXMGPGDLi5ucnroFGjRqhbt265rgUAISEhCAkJUdr+xx9/QF9fH/Hx8bCxscHGjRthaGgIAOjRowcuXryIq1evYvLkyfJjxGKxvEXG2toaERERMDU1xeLFiwEAjo6OiIyMxI0bN9C/f/8S42rdujXs7e1x+PBhvP/++wCAX3/9FcnJycX2qSjqm2++QVBQEICCfhg9evQovTJKwESBiKgaMDc3R2hoqNJ2S0tLAIC+vj62b98OoOCXbkJCAu7fv4/o6GgAQE5OjvwYU1NTeUJQ+F0oFCpsa9SoEdLS0soU26BBg7Bq1SocPXoUU6ZMweXLl5GUlFSmh1Z5CAQCeHp64sCBA1i0aBHq1KmD48ePQywWo1WrVgAKkomQkBDcvn0bvXr1grOzM3x9KzYkdcSIEQpN9YUKOzM6OTnByckJOTk5+Ouvv/DgwQPcuXMHqampaNy4scIx7dq1k//ZxMQEurq6sLe3V7i3hg0b4tWrV2WKzcvLC0uWLEFKSgrMzMzwww8/wNbWFm3atCn1WCcnJ9ja2uL69ev48ssvkZOTg5UrV5bpusVhokBEVA3o6enBzs6uxDI///wzAgMDce/ePRgaGsLa2hr169cHAIVXD4W/fl9Xr169CsdmbGwMNzc3HDlyBFOmTMGhQ4dgZ2dXpodWeQ0ePBhhYWH45Zdf0LlzZ5w9exazZs2S71+/fj22bNmCyMhInDx5Ejo6OujRoweWLVsmT6rKqmnTpiXWeX5+PtatW4fw8HBkZGTAwsIC7dq1g76+vkJ9A+qv84EDByIwMBDHjh3D8OHDERUVhXnz5pXp2MKRFF26dEFubi5CQ0Mxa9YsmJmZVSgW9lEgIqoBEhMTMWPGDNjY2CAqKgrXr1/H3r174eLiUiXXHzZsGO7evYvY2FhERUVhyJAhGrlO69atIRaLceLECURHR0MqlWLgwIHy/cbGxpg/fz7Onj2LyMhIzJkzB9evX8fy5cvVHktYWBi++eYbLFq0CDExMTh37hw2bdpU5j4XlWFkZIR+/fohMjIS58+fR25uLgYNGqSyfEpKCg4ePKjUAVIsFkMmk+HJkycVjoWJAhFRDXDr1i1IpVJMnTpV4RXCzz//DABKIwXUrVu3brC0tERgYCCys7NLfGhV1nvvvYfz588jMjISjo6O8mb+5ORkODs748SJEwCAli1bYvLkyejRoweSkpLUHse1a9cgEokwdOhQ+RwTKSkpuHPnjsbrGyh4/XDr1i3s3bsXrq6uaNiwocqy//77L/z9/eV1U+jixYvQ09NDixYtKhwHXz0Q0RshPjW+0pMeVfb6IqFIY+cXi8WoU6cO1qxZg3HjxkEqlSIiIkI+PDEzM1Nj1wYK3rEPGTIEoaGhcHd3L/GhVVmDBg3C6tWrERUVJe+UBxT047C0tERAQAAkEgmaN2+OW7du4fz585g+fbq83I0bNyAUCtG8efNKxdGuXTts3rwZ27Ztg729PR48eICtW7ciOztb4/UNFPTHsLS0xNWrVxEWFlZiWbFYDBcXF6xatQpZWVlo0aIFzp07h/DwcHzyyScwMjKqcBxMFIioxrMWWms7BIiEIo3GYWVlhbVr1yI0NBRTp05Fw4YN4eDggD179mDMmDGIiYmRd/jTlN69eyM0NFTtnRiLaty4MRwdHRETEyMf3VAoJCQEwcHB2LhxI54/fw4LCwvMnDlTYQSCt7c3hgwZglWrVlUqDh8fHzx//hy7du1CWloaLCwsMHjwYAgEAoSFhUEikVTqAVwagUAAZ2dnnDp1Cj179iy17Lp167B582Zs374dT548wTvvvINly5ZVenikQFa0R4aWxMbGYtiwYThz5gzMzc3l2/v27YvExESl8r/88ov8PdHNmzexevVq3Lp1C4aGhhg6dChmzpxZoeEynTp1AoAyzXHuvfUXXE5IRdcWQuzz6V7uaxFR+Tx48ABAwUOTql5oaCj279+P6Oho6OjwzbWm5efnw93dHR4eHpg9e7bGrlPa36tq0aJw7949+Pj4KEx0ARSMBX748CHmzp0rH1NbqEGDBgAKbnDcuHFo3749NmzYgL///hvr16+HRCKRj18lIqKKi4iIwJ07dxAeHo65c+cySdCwtLQ07Nq1C7/99htSUlIwcuRIrcaj1UQhNzcX+/btw9q1a4v99R8fHw+ZTIY+ffqobFILCwuDsbExNm/eDD09PTg7O8PAwAABAQHw8fGp8HAQIiIqEBcXh++//x4eHh7ydQRIcwwMDPDdd99BJpNh5cqVWn+OaTUtvHbtGoKDgzFhwoRix4fGxsZCX18f77zzjspzXLx4ES4uLgqrfbm7uyMvL08ty2sSEdV2/v7+uHHjBlavXo06dapFQ/QbrW7durhw4QIuXryIAQMGaDsc7bYotGrVClFRUTA1NUVERITS/vj4eDRq1Ahz5szBxYsXkZeXh969e8Pf3x9NmjRBZmYmkpKSlIZ9CIVCGBkZISEhQemchX0QVElLS9PIUqtEREQ1kVZbFBo3bgxTU1OV++Pi4vD06VO0adMGW7ZswcKFC3H16lWMHTsWWVlZ8ulHi+t1amhoCIlEorHYiYiIaoNq3Yb0+eefQyaTyefL7tSpE1q1aoWRI0fiyJEjcHZ2BlAwLKQomUxWbIeb0kYzlNbiQEREVJtU666r7dq1U1hUAwA6duwIY2NjxMXFyVsSims5yMjI4CsEIiKiSqq2iUJGRgYOHjyIuLg4he0ymQw5OTkwMTGBoaEhzMzM5GNACz179gwSiaRSU1YSERFRBROF16eufP78OcLDw/F///d/ePHihbrigr6+PoKCgpSWXT1z5gyysrLk8yo4OjoiOjoa2dnZ8jInT56Erq6u0twLREREVD7l6qPw6tUrzJ49G69evcL+/fshkUjg5eWFpKQkyGQybN68GXv37lVYsKSidHV1MW3aNKxatQoBAQFwdXXFnTt3EBISgj59+qBr164AgEmTJuH48eOYMmUKPvroI9y/fx/r1q3Dhx9+iGbNmlU6DiKq/pIDAyGNjSu9oIbpv2sNc39/bYdBpFblalHYsGEDLl++jF69egEADhw4gMePH2P+/PnYvXs3dHR0sGHDBrUFN378eAQEBODy5cuYNm0aduzYgeHDh2PdunXyMq1atcKOHTuQkZGBTz75BDt37sT48ePx2WefqS0OKoNIP2DnQMVPpJ+2o6JaQhobh6w47SYKWXFxFU5W/Pz8IBKJVH4uXbqk5miVubq6yv/dvHz5MkQiUZmmsi/JL7/8gtGjR6Nz585wdHTEzJkz8fDhQ3WEW2kl1beTk5PGr1+0jv38/NC3b99KnTM7OxsbNmyAq6sr7O3t4e3trZb5hMrVonD27FmMHj0an3zyCQDI50CYMGECAGDUqFHYuXNnhQIZOnRosQuNfPDBB6UuaNGpUyd8//33FbouqUnyzYKPud1/34mqkIG1Naz27Nba9R+MGVup483NzbFx48Zi97Vu3bpS5y6L0NBQtXYAv379OiZOnIg+ffogODgYGRkZ2Lx5M0aMGIGjR4/CxMREbdeqKG9v72KfO69P4KcpYrEY+/btU+v/tsuXL8exY8cwd+5ctGzZEocOHYKPjw/27NmDDh06VPi85UoUnj17hjZt2gAomJjoxo0bCrNGmZiYVMnSm1RNmdsB448X/HnnQMV9kX6KyYO5HeBRuZXdiN4kenp6cHBw0Nr1bWxs1Hq+r7/+Gq1atcLGjRvlQ9U7dOiA3r174/Dhwxg3bpxar1cR5ubmWqtzIyMjtV772bNnOHjwIObNmyefZrtbt264du0a/u///q9SiUK5Xj2YmZnJm42ioqLkMyUWun79OiwsLCocDL3BClsciv65GEFXgjD+xPgSP0FXglQeT/SmysjIwJo1a9CvXz/Y2tqiQ4cOmDhxosLoMD8/P/j4+CA8PFzeBD1x4kT8+++/OHDgANzc3NC+fXuMGzcOjx49kh/3+quH1929excikQgHDx5U2H7v3j2IRCKcO3eu2FjbtWuHjz76SGE+GzMzMxgbG6t8/bBw4UL06tUL+fn5Ctv9/f3h6uoKmUyG1NRUzJ07F46OjmjXrh0GDx6MQ4cOlVZ1FRYbG4sZM2agW7duEIvFcHJywooVKyCVSuVlRCIR9u3bh3nz5qF9+/bo1q0bQkNDIZFIsHDhQnTs2BGOjo5Ys2YNChdsLun1TlBQEBwcHJCenq6wPTg4GD179kReXp7SMQ0aNMD+/fvh5eUl36ajo4M6deooxFoR5UoUXFxcsGvXLgQEBGD16tVo2LAhXF1dkZKSgoCAABw+fBgDBw4s/URUOxW2OBS+nlAhLjUO8anxKvfHp8YjLlX7HdeI1C03N1fpU/hgAYAFCxbIm5N37NiBhQsXIj4+HvPmzVMod/XqVRw8eBCLFi3C4sWLceXKFYwZMwZ79uyBn58fPvvsM/z+++8ICAgoNaY2bdrAzs4Ohw8fVth+6NAhNGnSRN5nraipU6di2LBhCtuuXLmCly9fqmxuHzx4MJ48eaLw8MzOzkZUVBQ8PT0hEAgwf/58/P3331i6dCnCwsJgY2MDX19fXL58udR7KSo/P7/YOi+UkpKCUaNGQSqVIigoCNu2bcOAAQOwe/du7N6t+JorKCgIJiYm2Lx5M1xcXBASEoJhw4ahXr162LhxI1xcXPD111/j1KlTpcbl5eWFzMxMnD59WiHWo0ePwtPTE7q6ukrH1K1bF2KxGA0bNkR+fj6SkpKwYsUKJCYmlvr6vjTlevUwf/58ZGZm4sCBAzAzM8OSJUtgYGAgX37U09MTkydPrlRARAAgEoqw0734/i7jT4yv4miINC8xMRFisVhp+5IlSzBixAhIpVJkZmZi0aJFcHd3BwB06dIFEokEq1atwvPnzyEUCgEA6enp2Lhxo3wE2unTpxEdHY2oqCj5ttjYWBw7dqxMsXl5eWHp0qVISkqChYUF8vPzceTIEZUPreKkpqZi0aJFMDc3x+DBg4st07VrV1hYWOD48ePy4e0XLlzAy5cv8d577wEoSDZmzJgBNzc3eR00atSo2BWISxMSEoKQkBCl7X/88Qf09fURHx8PGxsbbNy4EYaGhgCAHj164OLFi7h69arC804sFstbZKytrREREQFTU1MsXrwYQMFQ/sjISNy4cQP9+/cvMa7WrVvD3t4ehw8fxvvvvw8A+PXXX5GcnFxsn4qivvnmGwQFFbS6ent7o0ePHqVXRgnKlSjo6ekhICBAKQu1trbGTz/9hCZNmlQqGCKi2src3Fxp3hgAsLS0BFAwt8z27dsBFPzSTUhIwP379xEdHQ0AyMnJkR9jamqqMEzd1NQUQqFQYVujRo3k6+WUZtCgQVi1ahWOHj2KKVOm4PLly0hKSirTQwsAnjx5gokTJ+LJkyf45ptvUL9+/WLLCQQCeHp64sCBA1i0aBHq1KmD48ePQywWo1WrVgAKkomQkBDcvn0bvXr1grOzM3x9fcsUR1EjRoxQaKovVNiZ0cnJCU5OTsjJycFff/2FBw8e4M6dO0hNTUXjxo0VjmnXrp38zyYmJtDV1VWYWVggEKBhw4Z49epVmWLz8vLCkiVLkJKSAjMzM/zwww+wtbWV9xMsiZOTE2xtbXH9+nV8+eWXyMnJwcqVK8t03eKU69XD2LFj8csvvyht19PTQ5MmTXD27Fm+eiAiqgA9PT3Y2dkpfQpbCQDg559/hoeHB5ycnDB9+nQcOXJE/lB7/dVD4a/f19WrV6/CsRkbG8PNzQ1HjhwBUPDawc7OrkwPrfj4eHh7eyMlJQVff/210rT8RQ0ePBipqan45ZdfkJWVhbNnzyq0QKxfvx7jxo3DrVu38Pnnn8PZ2RkTJ07EP//8U+77atq0abF1Xrh+UH5+PoKDg9GlSxcMHDgQy5cvR2xsLPT19RXqG1B/nQ8cOBB6eno4duwY0tPTERUVVebErHXr1ujSpQumTp0KHx8f/PDDD0hJSalwLCW2KGRmZuL58+fy71euXEHfvn1hZWWlVDY/Px8//fSTQucYIk2JT41XeAVhLbSGb5eK/aogqgkSExMxY8YM9O3bF2FhYfLWgfDwcPz8888av/6wYcMwbtw4xMbGIioqCnPmzCn1mCtXrmDatGkwNjZGeHh4mRKL1q1bQywW48SJE5BIJJBKpQo/QI2NjTF//nzMnz8f9+7dw5kzZ7B582YsX74cW7ZsqdQ9FhUWFoZvvvkGy5YtQ9++feXDR4v2vdAEIyMj9OvXD5GRkbCwsEBubi4GDRqksnxKSgouXLgAd3d3haRFLBZDJpPhyZMnMDMzq1AspSYK77//vrx5SiAQIDAwEIGBgcWWl8lkcHR0rFAgRGVlLbRW+F5Sx0eiN8WtW7cglUoxdepUhVcIhUlC0ZEC6tatWzdYWloiMDAQ2dnZJT60ACAuLg4+Pj6wtLTE9u3by/WQeu+99/D1118jLS0Njo6O8mb+5ORkeHt7Y+HChXB3d0fLli3RsmVL3LhxQyM/Uq9duwaRSKTwSz4lJQV37typkmGVXl5eGDduHPbu3QtXV1c0bNhQZdl///0X/v+bFfT11ykXL16Enp5epdY+KjFREAqFWLNmDW7evAmZTIYvv/wSffv2hUgkUiqro6MDoVDIVw+kcUVbDti5kYCCmRErO+lRZa9vYG1desEKEovFqFOnDtasWYNx48ZBKpUiIiJCPjxR03PYCAQCDBkyBKGhoXB3dy/xoQUAn3/+OXJycvDxxx8jKSkJSUlJ8n1F+1AUNWjQIKxevRpRUVHyTnlAQT8OS0tLBAQEQCKRoHnz5rh16xbOnz+P6dOny8vduHEDQqEQzZs3r8QdF/Q72Lx5M7Zt2wZ7e3s8ePAAW7duRXZ2dpXMGdS1a1dYWlri6tWrCAsLK7GsWCyGi4sLVq1ahaysLLRo0QLnzp1DeHg4PvnkE/lqyxVRamdGZ2dnODs7AwAeP36M4cOHl/qOiYioKum/q7kHdFkZWFtrNA4rKyusXbsWoaGhmDp1Kho2bAgHBwfs2bMHY8aMQUxMjLzDn6b07t0boaGhpb4rf/z4MW7eLJgr5dNPP1XaP2zYMKxYsULl8Y0bN4ajoyNiYmLkoxsKhYSEIDg4GBs3bsTz589hYWGBmTNnKoxA8Pb2xpAhQ7BqVeUmdfPx8cHz58+xa9cupKWlwcLCAoMHD4ZAIEBYWBgkEkmlHsClEQgEcHZ2xqlTp9CzZ89Sy65btw6bN2/G9u3b8eTJE7zzzjtYtmxZpYdHCmRFe2TUcp06dQKAMs1x7r31F1xOSEXXFkLs8+mu6dCqt8KZGIvOzFjc96L7iihsIVA1PLKy5anmKlxSvrh+UqR5oaGh2L9/P6KjoxUmUiLNyM/Ph7u7Ozw8PDB79myNXae0v1flGh4JAD/99BOOHj2Kp0+fFjs7lEAgwK5du8p7WiIiqqYiIiLk8+XMnTuXSYKGpaWlYdeuXfjtt9+QkpKCkSNHajWeciUK4eHh8jkUTE1Nq2ThDCIi0q64uDh8//338PDwkK8jQJpjYGCA7777DjKZDCtXrqzwaAV1KVeisHv3blhbW2Pbtm1Kk00QEdGbyd/fX96jnjSvbt26alkeWl3KlSgkJSXB39+fSQKVTfLN//ojvL4E9WuCrgQprdsQnxoPkVB5ZA0REVW9cr1oat68OZ4+faqpWOhNYm6nmBgU/f4/xS0AJRKKlOZKICIi7ShXi8KUKVOwYsUK9O/fv0wzbFEt5lH2YUklLQBF9DodHR2FNQ2IqPLy8vJKXFSrXInCtWvXYGhoiMGDB6NFixYQCoXyObELcdQDEWmKgYEBJBIJUlNTFdZAIKKKSU1NhVQqlU9PXZxyJQqFU4Wam5sjMzOzQotwEBFVVOPGjSGVSpGSkoIXL16UeYljIlKWl5cnTxJK6ntYrkTh7NmzlQ6MiKiiBAIBLC0t8fTpU2RlZWl8fQOiN1ndunXlSULRtwOvK/eES0RE2iQQCNCkSRNth0FUa5QrURg7tmwLruzevbtCwRDRf5IDAyGNjSu1nP671jDnGHci0pByJQrFLeOZn5+P58+fQyqVwtLSkqMhiNREGhtX6oqEWXGlJxJERJWhlj4KeXl5OHPmDD7//HNMnDhRLYHVJJcTUrH06J/4wlOs7VCohiqu9aAwSbDao7qFTpvLKhNR7aCWlT10dXXRr18/fPDBBwgODlbHKWuc249faTsEqsEKWw9ep+lli4mIykKtnRnfeecdfPvtt+o8ZbXXtQXHcpN6lNZ6QESkDWpbKzQ7OxtHjhyBqampuk5JREREWqaWUQ/Z2dlISEjAq1evMHPmTLUERkRERNpX6VEPQEEfhZYtW2LQoEEYOXKkWgIjIiIi7ePMjERERKRShToz5uXl4datW/jnn3+gp6cHCwsLiMUcGkhERPSmKXeiEB0djaVLlyIlJQUymQxAwZSqTZs2xRdffAFXV1e1B0lERETaUa5EISYmBjNnzoSpqSlmz56NVq1aQSaT4d69e9i7dy8++eQT7N69Gx06dNBUvERERFSFyjU8MiQkBJaWljh27BimTJmCPn36wM3NDVOmTMGxY8dgaWmJr776qkKBxMbGQiwWIzk5WWH7hQsX4OXlBXt7e7i6umLHjh1Kx968eRNjxoxB+/bt0bNnT6xbtw45OTkVioOIiIj+U65E4Y8//sAHH3wAY2NjpX1GRkYYNmwYfv/993IHce/ePfj4+CA3N1dh+/Xr1zF16lS0bNkSISEh8PT0xOrVq7F9+3Z5mQcPHmDcuHHQ19fHhg0bMGHCBOzcuRMrV64sdxxUhZJvAjsHFvw3+SYQ6aftiIiIqBhqnZlRIBCU65d8bm4u9u3bh7Vr16Ju3bpK+zdt2gQbGxusWbMGAODk5ITc3Fxs2bIFY8aMgZ6eHsLCwmBsbIzNmzdDT08Pzs7OMDAwQEBAAHx8fGBmZqa2+yM1MbdT/J6dXpAsEBFRtVOuFgV7e3scOHAAGRkZSvskEgn2798POzu7Yo4s3rVr1xAcHIwJEyZg3rx5CvukUiliYmLQr18/he39+/fHq1evcP36dQDAxYsX4eLiAj09PXkZd3d35OXl4cKFC+W5PaoqHquA8ccLPuZ2gJ6htiMiIiIVytWi8PHHH2Ps2LEYNGgQRo8ejXfeeQcA5J0ZU1JSsHTp0jKfr1WrVoiKioKpqSkiIiIU9j18+BA5OTlo0aKFwnYrKysAQEJCAuzt7ZGUlKRURigUwsjICAkJCUrX7NSpU4kxpaWlFftqhYiIqDYqV6LQqVMnhISEYNmyZVi9ejUEAoF8iGSTJk2wfv16dOvWrczna9y4scp9aWlpAAr6PrzO0LDg16dEIlFZprCcRCIpcyxUAZF+/70ySL6p/EqBiIhqvHL3UejTpw969+6NP//8Uz6ls5mZGezt7VGnjvq6PLw+R0NxdHR0Siwjk8mgo6P8ZiUmJqbE65bW4qDK5YRULD36J77wrEUTTxV2RDS3++9DRERvlDL1Ufj222/h6ekpH5Wgq6uLdu3aYcCAAThz5gxmzpyp9uWlC5v/i7YKFH43NjaWtyQU13KQkZFRZa8QbJo1AADcfvyqSq5XrZjb/dffwGOVtqMhIiI1KzFRkMlkWLBgAQICAvDkyRM8fvxYqcxbb70FHR0dBAUFYc6cOWoLrHnz5tDV1UViYqLC9sLvLVq0gKGhIczMzPDgwQOFMs+ePYNEIlHqu6ApX3iK0bWFsEquRUREVJVKTBT279+PI0eOYOTIkfjpp5/QvHlzpTKzZ8/GmTNnMHjwYERGRuLQoUNqCUxfXx+dOnXCqVOn5K8YAODkyZMwNjaGra0tAMDR0RHR0dHIzs5WKKOrq4suXbqoJRaq/uJT4zH+xHj5J+hKkLZDIiJ6I5TYqWD//v3o3LkzFi9eXOJJ9PX1ERgYiPj4eHz33Xd4//331RLctGnTMH78eMyePRtDhgzBb7/9hu3bt2Pu3LmoV68eAGDSpEk4fvw4pkyZgo8++gj379/HunXr8OGHH6JZs2ZqiYOqN2uhtcL3+NR4LUVSNsmBgZDGxilsy4qLg4G1tYojiIi0p8RE4a+//sKnn35aphPp6Oigf//+2Lp1q1oCA4Du3bsjJCQEmzZtwowZM2BmZoYFCxZgwoQJ8jKtWrXCjh07sHr1anzyyScwMTHB+PHjMXPmTLXFQdWbbxdfhe/jT4zXUiTFK5oYZFy9CgCo37mzfJuBtTX0361YopAVF4cHY8bKv+u/aw1zf/8KRktEpKjEREFXV1dhIqPSmJiYFDvSoCyGDh2KoUOHKm3v27cv+vbtW+KxnTp1wvfff1+h6xJpmjQ2TqHFoH7nzmp7mBdNLrLi4lSUJCKqmBITBSsrK9y6davMJ7t58yab+4mKYWBtDas9u9V+3qLJxustC0RE6lDiz/+BAwfi6NGjuHv3bqknunv3Lo4ePQonJye1BUdERETaVWKi4O3tjWbNmmHMmDE4cuQI8vLylMrk5+fj2LFjGD9+PAwNDfHRRx9pLFgiIiKqWiW+ejA0NMRXX32F6dOnw9fXF0uXLoVYLEaTJk2Qn5+PZ8+e4c8//0RGRgYsLCzw5ZdfomnTplUVOxEREWlYqXMut2zZEkeOHEF4eDiOHz+O69evy2dorFu3LhwcHNCvXz94e3uXq+MjERERVX9lWpxBT08P48ePx/jxBcPOUlNToauri4YNG2o0OHrzBF0JQlzqfz3z41PjIdJiPEREVLIKreIkFHK6YqqYuNS4guRAWJAeiIQiWCfHFywutXPgfwXN7bh2BBFRNaC+5R6JykgkFGGn+87/NkT6AbKb/31Pvql8EBERaQUTBdK+oi0Hr7csEBGRVjFRoDdS4SJRhayF1kpTPRMRUemYKNAbp6YtEqVuRdd+ALj+AxFVHBMFeuNU90WiNKm4haW4/gMRVQYTBaI3SHGtBlz/gYgqg4kClV2kn+KIhOSbBcMYiYjojcVEgcou+aZicmBux0ShiOTAQEhjFZv6X19imoiopmGiUEFLj/6Jywmp6Nqilk0+ZW4HjD+u7SiqLWlsnFJiYGBtXWzfASKimoCJQgXdfvwKAGDTrIGWI6HqxsDaGlZ7dms7DCIitShxmWkqWdcWQnzhKdZ2GERERBrDRIGIiIhUYqJAREREKjFRICIiIpWYKBAREZFKTBSIiIhIJQ6PpOop+eZ/y02b2ykvRU1ERFWCiQJVP6/P9vj6lNFERFTlmChQ9fN660FhqwIREWkF+ygQERGRSkwUiIiISCUmCkRERKQSEwUiIiJSiYkCERERqcREgYiIiFRiokBEREQqVft5FHJzc9GhQwdIpVKF7fXr18dvv/0GALhw4QLWr1+Pv/76C6amphg9ejQmTJigjXCJiIjeKNU+UUhISIBUKkVQUBDeeecd+XYdnYLGkOvXr2Pq1Knw8PDAp59+imvXrmH16tWQyWSYOHGilqImIiJ6M1T7RCEuLg46Ojro378/6tWrp7R/06ZNsLGxwZo1awAATk5OyM3NxZYtWzBmzBjo6elVdcj0mqArQYhLjZN/j0+Nh0go0mJERERUHtW+j0JsbCyaN29ebJIglUoRExODfv36KWzv378/Xr16hevXr1dVmKRCXGoc4lPj5d9FQhGshdZajKh2yoqLw4MxYxU+yYGB2g6LiGqAat+iEB8fDz09PUycOBHXr19HnTp14OHhgQULFiA5ORk5OTlo0aKFwjFWVlYACl5bdOvWTWFfp06dSrxeWloajI2N1XsTtZxIKMJO951ajSE+NR7jT4yXf7cWWsO3i68WI6o6+u8qJ2ZZcXHFlCQiUlbtE4W4uDhIJBJ88MEHmDp1Km7duoWQkBAkJCRgzpw5AAAjIyOFYwwNDQEAEomkyuOl6qdoC8brLRy1gbm/v9K2B2PGaiESIqqJqn2isH79ejRs2BAiUcF77c6dO8PU1BTz58/HxYsXAQACgaDYYws7PL4uJiamxOuV1uJANU/RloPXWxaIiKhk1T5R6NKli9K23r17K3wv2nJQ+J2vEEjTkgMDIY39rxk/Ky4OBtbsg0FEb45q3Znx2bNn2L9/Px4+fKiwPSsrCwBgamoKXV1dJCYmKuwv/F607wKRuklj4xTe9xtYWxfbJ4CIqKaq1i0KAoEAixcvxtixY7Fw4UL59h9//BG6urro0aMHOnXqhFOnTuGjjz6Sv4I4efIkjI2NYWtrq63QqRYxsLaG1Z7d2g6DiEgjqnWiIBQKMWrUKOzZswdGRkbo1KkTrl27hi1btmDUqFGwsrLCtGnTMH78eMyePRtDhgzBb7/9hu3bt2Pu3LnFDqkkIiKisqvWiQIA+Pr6wszMDAcPHkRYWBjMzMzwySefYNKkSQCA7t27IyQkBJs2bcKMGTNgZmaGBQsWcApndYj0A5Jv/vc9+SZgbqe9eIiIqMpV+0Shbt26mDx5MiZPnqyyTN++fdG3b98qjKqWSL6pmByY2zFRICKqZap9okBaZm4HjD+u7SiIiEhLqvWoByIiItIuJgpERESkEl89ENVShQtFFdJ/17rY6Z6JqHZjokBUCxWdFIqLRBGRKkwUqPpLvgnsHPjfd3M7wGOV9uJ5AxRtOeAiUUSkChMFqt6KDsd8fV4HIiLSOCYKVL0VbTl4vWWBiIg0jqMeiIiISCUmCkRERKQSEwUiIiJSiYkCERERqcTOjKRWQVeCEJf635j8+NR4iIQiLUZERESVwUSB1CouNU4hORAJRbAWWpdyVNWLT43H+BPj5d+thdbw7eKrxYiIiKonJgqkdiKhCDvdd2o7DJWKJi7xqfFaioSIqPpjokC1TtGWg9dbFoiISBE7MxIREZFKTBTU6HJCKpYe/VPbYRAREakNEwU1sWnWAABw+/ErLUdCRESkPuyjoCZfeIqZJFQVriZJRFRlmChQzcLVJDUmKy5Oablp/XetlZakJqLahYkCKYr0++/hm3xT+cGsbVxNUiP031We6yIrLq6YkkRU2zBRqO1eTwwA4MGFgv9a9SxIEqpbokAaUVyrQdHWBSKqnZgo1HbJNxVbDgoTBL7zV5IcGAhprOKv7Ky4OBhYV7+ZJ4mI1IWJAhUkBuOPazsKrSo6pbPXjy9h99xIoUzG1asAgPqdO8u3GVhbF9tsT0T0pmCiQLVe0SmdY1JiMDAuF6lPdPHEsj4AoH7d+mjeuTM79xFRrcNEgWq9olM6B10JgnHdg3hiCXw3U/y/Ra5aVOv1KzSl6EgIJkpEtQ8nXFIzzs5Y8/l28YW10BrWQmvsdN9Za5fJ1n/XWqH/RVZcnFIfDSJ687FFQY1smjXA5YRUTrxEb4SiLQccBUFUOzFRUKPaODtj0JUgxKX+9yuzoJm+in+Bvz5TI0dsEBGpFRMFqpS41DiF5EAkFCl1DtSo1+d54CyNRERqx0SBKk0kFGmvo9/rrQc7B1Z4HYiicyRwfgQiogJMFGqbojMxVsdpmiuqEutASGPjFJIDzo9QPK4HQVT7vDGJwrFjx/DVV1/h4cOHsLS0hI+PD95//31th1X9FJ2J8U2aprmS60AYWFvDas9uNQb0ZuF6EES10xuRKERGRmLevHkYO3YsevXqhaioKPj6+sLAwADu7u7aDq/64UyMVAFcD4KodnojEoV169bBw8MD/v/7h6xXr154+fIlNm7cyEShkoqOaihKK6McyqNonwUAyVFPIX2SrbAt60k2DN4yUXmaolM8WwutlSZqqq34OoLozVbjE4WHDx8iMTERc+bMUdjev39/REZG4uHDh3j77be1FB2w9Oif8iGTNs0a4AtPcdUHUYmlo4uOaiiqykc5lIeK+5Q+yS5IDJrqybcZNMqBfqOcYssXN8VzTEpMiQlUced4ExOL4l5HZFy9ioyrVxU6hzJxIKq5BDKZTKbtICrj/PnzmDJlCg4fPgzr13qp3759G0OGDMG2bdvg5OQk396pU6cSz5eWlgYAMDY2LrFcenYuAMBQr47S9rz84qtUV0egVL5KZKcD+bmAzv+uLdAF6hqU6dCM3AwAQP069TUVXZXLzyi4J536r91TdnrBf/UMSz1emidFniyvzNfLy8+Dro7uG1WHJZFlZUGWn//f97w8CHR1FetbBWNjY0RHR2syPCIqpxrfolD4YDcyUlzpz9Cw4B98iURS7nMKBIJSyxQ+8IsmFlpJBEpThoefKhV9uJU14dKGYh9Y5agjfV39Cl+7OteLuggMDFD63yBFtaFeiGqqavhUK5/CBpGiD/fC7To6istZxMTEqPX6hS0U6j5vTcd6KR7rpXiF9cLWBKLqp8YvClX4C6Roy0F6errCfiIiIiq/Gp8otGjRAgCQmJiosP3BgwcK+4mIiKj8anyiYGVlhbfeegsnTpxQ2H7q1Cm88847aNasmZYiIyIiqvlqfB8FAJgxYwYWLlyIhg0bonfv3jh79iwiIyOxfv16bYdGRERUo70RicLQoUORnZ2NHTt2YP/+/Xj77bcRFBSEAQMGaDs0IiKiGu2NSBQAYPjw4Rg+fLi2wyAiInqj1Pg+CkRERKQ5NX5mRiIiItIctigQERGRSkwUiIiISCUmCkRERKQSE4Uijh07hoEDB6Jdu3bw8PDAoUOHSiyfnp6OpUuXwtHREe3bt8fkyZNx//59hTK5ubnYsGEDnJ2dYW9vj5EjR+KPP/7Q3E1ogCbqRSKRICgoCG5ubnBwcICnpyf27t2LmtRtRhP18jqJRAIXFxd89tln6g1cwzRVL9999x08PDxgZ2eH/v37Y/fu3Zq5ASL6j4zkfvzxR5lIJJKtWLFC9tNPP8kWL14sa9u2rSwyMlLlMZMnT5Z169ZNFhERITt58qTM09NT1qtXL9mrV6/kZZYsWSKzt7eX7dmzR3bmzBnZ6NGjZe3bt5clJiZWxW1VmqbqZdKkSbIuXbrIvv32W9mlS5dkwcHBMmtra9mWLVuq4rYqTVP18jp/f39Z27ZtZf7+/pq6DbXTVL3s2LFDZm1tLVu/fr3s0qVLsrVr18ratm0r27t3b1XcFlGtxUThNW5ubrJZs2YpbPv0009l7u7uxZa/evWqrG3btrLz58/Ltz179kzm4OAg27p1q0wmk8kePnwoe/fddxX+MZNKpbLevXvLFi9erIG7UD9N1Mvt27dlbdu2lf34448Kxy5evFjWsWNHNd+BZmiiXl537tw5Wfv27WUdO3asUYmCJupFIpHIHBwcZOvXr1c4ds6cObIZM2ao9waISAFfPfzPw4cPkZiYiH79+ils79+/P+7du4eHDx8qHXPx4kUYGhrC0dFRvk0oFKJz58746aefAAC//vor8vLy0L9/f3kZPT099O7dW16mOtNUvchkMnh7e6N79+4Kx7Zs2RJpaWl4/vy5Bu5GfTRVL4VevnyJzz//HPPnz0eDBg00cxMaoKl6uXDhAjIyMjBy5EiFY9euXYvQ0FAN3AkRFWKi8D/37t0DoLzapJWVFQAgISGh2GOsrKygq6ursL158+by8vfu3UPDhg0hFAqVzvv48WNkZWWp7R40QVP1YmNjg2XLlqFRo0YKZaKiotCkSROl7dWNpuql0PLly9GqVasaN9uopuolPj4ejRo1QlJSEoYPHw5bW1s4OzuzjwJRFWCi8D9paWkAACMjI4XthoaGAAo6lRUlkUiUyhceU1i+pDJAQSeu6kxT9VKcXbt24cqVK5g8eTIEAkFlwtY4TdbL6dOncebMGQQEBFT7eihKU/WSmpqKnJwcTJs2Df369cO2bdvQp08frFixAhEREeq+DSJ6zRuz1kNlyf7X077oP8yF23V0lHMqWQm98wvLqyqj6nrVjabqpahvv/0WK1euhIeHB8aOHVvRcKuMpuolNTUVX3zxBRYsWIC33npLXeFWGU3VS05ODtLT0zFnzhyMHj0aANC9e3c8fvwYISEhGDp0qFriJyJlbFH4H2NjYwDKv3gKf/EX7n+dkZFRsS0C6enp8l9IJZUp3F+daapeCuXn5yMoKAjLly/HwIEDERwcXO2TJ0Bz9bJkyRK0atUKw4YNQ25uLnJzcwEUPEwL/1ydaapeClsknJ2dFcr06tULjx8/lrdkEJH6MVH4n8J3qomJiQrbHzx4oLC/6DEPHz5U+kX04MEDefmWLVvixYsXePnypVKZt956C3p6emq7B03QVL0ABb8SZ82ahR07dmDChAkIDg5GnTo1o5FLU/Vy8uRJXLlyBba2thCLxRCLxfjnn39w8OBBiMViPHr0SBO3ozaaqpfCPg7Z2dkKZXJycgBU/5Y5opqMicL/WFlZ4a233sKJEycUtp86dQrvvPMOmjVrpnRMz5498erVK1y6dEm+LTU1FTExMejRowcAyP978uRJeZns7GycP39evq8601S9AIC/vz9OnTqFhQsXwtfXt0b9Y6+pejlw4IDSp0mTJujTpw8OHDiApk2bavbGKklT9dKrVy8AwPHjxxWOjY6OhkgkqvYtc0Q1Wc34+VZFZsyYgYULF6Jhw4bo3bs3zp49i8jISKxfvx5AwT9eiYmJaN26NYyMjNC5c2d06dIFc+bMwbx589CoUSOEhITA2NgYI0aMAABYWlpiyJAhCAgIQEZGBqysrLBz5068fPkSkyZN0ubtlpkm6uXcuXM4cuQIXF1d4eDggBs3bihc08bGptq3tmiiXuzs7JSuo6enBxMTk2L3VUeaqJfmzZtjxIgR2Lp1K+rUqQMHBwccP34cv/76KzZv3qzN2yV681XttA3V3//93//J+vbtK7O1tZV5eHjIfvjhB/m+gwcPytq2bSv79ddf5dtevHgh8/Pzk3Xq1EnWoUMH2eTJk2V///23wjmlUqlsxYoVsu7du8vs7e1lI0eOlN24caOqbkkt1F0vfn5+srZt26r8JCUlVeXtVZgm/v9SlIuLS42acEkm00y95OXlybZs2SJzdXWV2drayjw9PWUnT56sqlsiqrUEMlkNmlifiIiIqhT7KBAREZFKTBSIiIhIJSYKREREpBITBSIiIlKJiQIRERGpxESBiIiIVGKiUEv4+flBJBLh3XffRWpqqspy7733HkQiEfz8/NR67OufwiWC58+fj7t375Yr/vKcZ8yYMXB1dS3T+V+XnZ2NlJSUch9HRPQm4syMtUx+fj6io6Ph5eWltO/hw4eIj4/XyLELFy6EiYkJACAzMxMPHjxAREQETp48iW3btqFr165lir8855k6dSoyMzPLdN5C//zzDyZMmAAfHx+uSEhEBCYKtc5bb72FM2fOFPuwj4qKglAoVNlqUJlj3dzclJZNHjt2LLy8vDBr1ixERUXJVwgsSXnO4+joWOr5inr06BHu379f7uOIiN5UfPVQy/Tp0weXLl1CVlaW0r7Tp0+X2FRfmWOLY2FhAV9fX6SmpuLgwYPlOlYT5yEiImVMFGoZNzc3ZGZmKqzUBwDPnj3Db7/9hn79+mnkWFXc3d2hp6eHn3/+udzHlnaeon0UsrOzsWLFCvTp00fev2Hp0qXyJcAjIiIwduxYAAWvOEQikfzYP//8EzNnzkSPHj0gFovRvXt3zJ07F8nJyfIyISEhsLOzw/379+Hj44P27dujc+fO8PX1xfPnzxXilUgkCAwMRO/evWFvbw9PT0/s379foczLly+xfPly9OrVC7a2tvDw8MCuXbuUlmMmItIkvnqoZTp27AgTExOcOXNG4SF65swZ1KtXD927d9fIsaro6+ujefPmiIuLK/ex5T3PsmXLcOzYMYwdOxZvv/027t69i/DwcDx48AA7duxA586dMXXqVGzZsgXe3t7o2LEjACA+Ph4jR46ElZUVpkyZgnr16uH69es4fPgwnjx5gj179sivkZ+fj7Fjx6JTp07w9fXFzZs3ceDAAWRlZWHjxo0AChKWUaNG4e7du/jwww9hbW2N8+fP4/PPP0dmZibGjh2LjIwMjB49GklJSRg5ciTMzc3x66+/IjAwEPfv38cXX3xRqfoiIiorJgq1jK6uLlxcXBAdHY38/Hzo6BQ0Kp0+fRq9e/cucWnnyhxbkgYNGiAxMbFCx5bnPEePHoWXlxfmzJkj31a/fn38/PPPSE9Px9tvv40ePXpgy5YtcHBwwODBgwEAe/fuhUAgwO7du9GoUSMAgLe3N3JycnD8+HG8ePFCvj03NxcDBgyQj/wYPnw4UlJSEBUVhczMTNSrVw8HDhxAXFwcgoOD4enpKT/f6NGjERYWhlGjRmH79u1ISEjAwYMH5S0bI0eOxLp167B161Z4e3vD2tq60nVGRFQavnqohfr06YNnz57hxo0bAAqawX/55Re4ublp9FhVcnNzIRAIKnx8Wc9jbm6OH3/8EREREXj16hUAYNasWTh48GCJHSmXLFmCs2fPypMBoOC+9fX1AQAZGRkK5T08PBS+v/vuu8jNzcWLFy8AAOfOnYNQKMSgQYPkZQQCAVavXo3w8HDo6Ojg1KlTaNu2LZo0aYLU1FT5p7Ceo6OjS68QIiI1YItCLdSzZ0/Uq1cPZ8+eRYcOHXD+/Hno6OjA2dlZo8eq8uLFCwiFwgofX9bzLFmyBLNmzcLChQuxaNEiODg4oG/fvvDy8oKxsbHK4wQCAZ4/f46tW7ciPj4eiYmJePz4sbyvQH5+vkL5ojEUtrTk5eUBKBiC2bx5c6WkxtLSUv7nxMREZGVlqXydk5SUpDJeIiJ1YqJQCxkYGKBHjx44c+YM5s2bh9OnT6NHjx5lGp5YmWOLI5FI8PDhQ/Tu3btCx5fnPN27d0d0dLT8c/HiRaxcuRLffPMNIiIiVCYZ586dw/Tp09G0aVN069YNTk5OsLW1xYULF7B161al8qW1juTl5ZWpTMeOHfHxxx8Xu79p06YlHk9EpC589VBLubm54d69e7hz5w5++ukn9O3bt0qOLerEiROQyWTo06dPhc9RlvNkZ2fj999/R1paGgYOHIjg4GBcvHgRCxYsQFJSEo4fP67y3MuXL4eVlRV+/PFHrFq1ChMmTECXLl2URjKUVbNmzfDw4UOl7efPn4evry+ePn0KS0tLpKeno0ePHgofsViMV69eoV69ehW6NhFReTFRqKVcXFygq6uLoKAgZGVllWsOhMoc+7onT55g06ZNMDMzk3fq09R5Xrx4AW9vb4UWAB0dHdjZ2cn/DBR02AQUXye8ePECzZo1Q/369eXbkpKScOrUKQD/vVIoKycnJzx9+hSnT59W2L5r1y6cO3cOJiYmcHV1RVxcHM6dO6dQ5quvvsKnn35a5qmviYgqi68eaikTExN07NgRFy5cQNeuXeXTImvq2KioKHk5qVSKe/fu4dChQ5BKpdi2bRsMDAzKdO2Knqdp06bw9PTE3r17kZmZifbt2+PFixf49ttv0bhxY3kHxMJzHzlyBDKZDEOGDIGTkxN+/PFHLF68GHZ2dnj06BG+//57+fTQ6enpZYq90PDhw3Hw4EHMnj0bo0aNQosWLXDu3DlcvHgRgYGB0NXVhY+PD06dOoWPP/4Yw4cPR5s2bXDt2jUcPnwYTk5OcHJyKtc1iYgqiolCLdanTx9cuXKlQhMllffYlStXyv9saGgICwsLuLq6YvLkyWjRokWZr1uZ8yxfvhxvv/02jh8/juPHj8vnfpg9e7a8f0KrVq0wZswYRERE4ObNm+jatSuWLFmC+vXr4+zZszh8+DDMzc3x/vvvo2/fvhgxYgR+/fVX2NjYlPkeDAwMsGfPHmzYsAHHjx9HWloaWrVqhQ0bNsgTlkaNGmHfvn3YtGkTTpw4gX379qFZs2aYPn06pkyZIm8BISLSNIGM07wRERGRCvxZQkRERCoxUSAiIiKVmCgQERGRSkwUiIiISCUmCkRERKQSEwUiIiJSiYkCERERqcREgYiIiFRiokBEREQqMVEgIiIilf4f7fNzWmg/sIUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure(figsize=(4,4))\n",
    "\n",
    "_bins = 20\n",
    "hist(shuffle_mmd, histtype='step', bins=_bins, label='Shuffle')\n",
    "\n",
    "hist(c2_c4_mmd, histtype='step', bins=_bins, label='Family 1 vs. Family 2')\n",
    "hist(c2_c5_mmd, histtype='step', bins=_bins, label='Family 1 vs. Family 3')\n",
    "hist(c4_c5_mmd, histtype='step', bins=_bins, label='Family 2 vs. Family 3')\n",
    "legend(bbox_to_anchor=(1,1))\n",
    "xlabel('MMD Distance')\n",
    "ylabel('Counts')\n",
    "sns.despine();\n",
    "# savefig('mmd_distance.svg', dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6b23f109",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ttest_indResult(statistic=-212.62560321226508, pvalue=0.0)\n",
      "Ttest_indResult(statistic=-251.6347167517595, pvalue=0.0)\n",
      "Ttest_indResult(statistic=-321.3126580974118, pvalue=0.0)\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import ttest_ind\n",
    "print(ttest_ind(shuffle_mmd, c2_c4_mmd, equal_var=False))\n",
    "print(ttest_ind(shuffle_mmd, c2_c5_mmd, equal_var=False))\n",
    "print(ttest_ind(shuffle_mmd, c4_c5_mmd, equal_var=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd36b825",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env: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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}