Fig1-scaling-multiple-3state-systems.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import msmtools\n",
    "from tqdm.notebook import tqdm\n",
    "\n",
    "import itertools\n",
    "from pyemma.util.statistics import confidence_interval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def transition_matrix_3states(p, n_subsystems):\n",
    "    \"\"\"\n",
    "    compiles a transition matrix of n independent 3 state systems.\n",
    "    :param p: float, 0 <= p <= 1, transition probability to switch state\n",
    "    :param n_subsystems: int, number of sub-systems\n",
    "    :return: np.ndarray; reversible transition probability matrix\n",
    "    \"\"\"\n",
    "    T_ = np.array([[(1-p), p/2, p/2],\n",
    "                  [p/2, (1-p), p/2],\n",
    "                  [p/2, p/2, (1-p)]])\n",
    "    \n",
    "    assert n_subsystems >= 2\n",
    "    \n",
    "    T = T_\n",
    "    for _ in range(n_subsystems-1):\n",
    "        T = np.kron(T_, T)\n",
    "       \n",
    "    assert msmtools.analysis.is_tmatrix(T)\n",
    "    assert msmtools.analysis.is_reversible(T)\n",
    "    return T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "trials = 30\n",
    "N = int(1e5)\n",
    "eps = 0.1\n",
    "n_systems = np.arange(2, 7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b9d86facf2204ca6a3b3042084c28fc1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/30 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sampled_independently = np.zeros((len(n_systems), trials), dtype=int)\n",
    "sampled_combinatorics = np.zeros((len(n_systems), trials), dtype=int)\n",
    "sampled_pairwise = np.zeros((len(n_systems), trials), dtype=int)\n",
    "\n",
    "for _trial in tqdm(range(trials)):\n",
    "    for num_subsys_n, num_subsys in enumerate(n_systems):\n",
    "        T = transition_matrix_3states(eps, num_subsys)\n",
    "        traj_comb = msmtools.generation.generate_traj(T, N)\n",
    "        \n",
    "        independent_trajs = list(np.unravel_index(traj_comb, tuple((3 for _ in range(num_subsys)))))\n",
    "        \n",
    "        # independent reversible sampling\n",
    "        for fin_n, fin in enumerate(range(5, N, 5)):\n",
    "            have_sampled_both_independently = True\n",
    "            for n, t in enumerate(independent_trajs):\n",
    "                c = msmtools.estimation.cmatrix(t[:fin], 1, sparse_return=False, nstates=3)\n",
    "\n",
    "                if 0 in c.sum(axis=1):\n",
    "                    have_sampled_both_independently = False\n",
    "                    break\n",
    "\n",
    "                T = msmtools.estimation.transition_matrix(c)\n",
    "                have_sampled_both_independently *= msmtools.analysis.is_connected(T)\n",
    "            if have_sampled_both_independently:\n",
    "                sampled_independently[num_subsys_n, _trial] = fin\n",
    "                break\n",
    "\n",
    "        # all combinatorial state reversible sampling\n",
    "        for fin_n, fin in enumerate(range(5, N, 5)):\n",
    "            have_sampled_combinatorics = True\n",
    "            c = msmtools.estimation.cmatrix(traj_comb[:fin], 1, sparse_return=False, nstates=3**num_subsys)\n",
    "            if 0 in c.sum(axis=1):\n",
    "                continue\n",
    "\n",
    "            T = msmtools.estimation.transition_matrix(c)\n",
    "            if msmtools.analysis.is_connected(T):\n",
    "                sampled_combinatorics[num_subsys_n, _trial] = fin\n",
    "                break\n",
    "                \n",
    "        # all pairwise reversible sampling\n",
    "        for fin_n, fin in enumerate(range(5, N, 5)):\n",
    "            have_sampled_pairwise = True\n",
    "            for traj_a, traj_b in itertools.combinations(independent_trajs, 2):\n",
    "                pairwise_traj_comb = np.ravel_multi_index(np.stack([traj_a, traj_b]), (3, 3))\n",
    "                \n",
    "                \n",
    "                c = msmtools.estimation.cmatrix(pairwise_traj_comb[:fin], 1, sparse_return=False, nstates=9)\n",
    "                if (0 in c.sum(axis=1)):\n",
    "                    have_sampled_pairwise = False\n",
    "                    break\n",
    "                T = msmtools.estimation.transition_matrix(c)\n",
    "\n",
    "                have_sampled_pairwise *= msmtools.analysis.is_connected(T)\n",
    "            if have_sampled_pairwise:\n",
    "                sampled_pairwise[num_subsys_n, _trial] = fin\n",
    "                break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.matplotlib.rcParams['axes.titlesize'] = 12\n",
    "plt.matplotlib.rcParams['axes.labelsize'] = 20\n",
    "plt.matplotlib.rcParams['legend.fontsize'] = 11\n",
    "plt.matplotlib.rcParams['xtick.labelsize'] = 11\n",
    "plt.matplotlib.rcParams['ytick.labelsize'] = 11\n",
    "plt.matplotlib.rcParams['xtick.minor.pad'] = 2\n",
    "plt.matplotlib.rcParams['xtick.major.pad'] = 3\n",
    "plt.matplotlib.rcParams['ytick.minor.pad'] = 2\n",
    "plt.matplotlib.rcParams['ytick.major.pad'] = 3\n",
    "plt.matplotlib.rcParams['axes.labelpad'] = 1\n",
    "plt.matplotlib.rcParams['axes.labelsize'] = 11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATMAAADCCAYAAAAl4SVWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABG/ElEQVR4nO2dd3xc5ZW/n3PvnaJiyd0ybsK4EeMm2+ASYsD09tsEwpLQkhDKkg2BTbLJ7kJC2oYUCJAGJIHQwiYhWYIJIWBjJ4DXFBuDCWDLFRtXuUiyypR7z++Pe2c0I42kkTQqI93n8xnp1ve+c2fud8457/ueV1QVHx8fn3zH6O0K+Pj4+OQCX8x8fHz6Bb6Y+fj49At8MfPx8ekX+GLm4+PTL/DFzMfHp19g9XYFOsPw4cO1vLy8t6vh4+PTzaxdu7ZKVUdkc2xeill5eTmvv/56b1fDx8enmxGRHdke67uZPj4+/YK8EjMRuUBE7q+uru7tqvj4+PQx8krMVHWZql5bWlra21Xx8fHpY+SVmPn4+Pi0hi9mPj4+/YK8bM308fHpv6gqGolg19RQEg6XZHueL2Y+Pj69TlLAamuxDx3CicZAhIBhZK1Rvpj5+Pj0Ck0CdhT70MGkgBnhEOagkHsMZJ1w0Y+ZdSOzZ8+moaGhw+eVl5fz9ttvd0ONOsb27dsZPnx4l8u4//77c1Sj7Mn2HooIR48e7YEa+YArYE5jI7EDVUQ2VdJYuZnYvn1gmpiDijGLixCrczbWgLfMNuysZ8U/jlLd4FBaYLB0ejEzxhXmpOz169fnpJx8JiFm1157bW9XxaeXSLfADuFEoykWWHHOrjOgLbMNO+tZ9kYN1Q0OANUNDsveqGHDzvqclJ/6q19eXs7XvvY1Fi5cSHl5OT/5yU+Sx7344ovMmDGDE088kRtvvJHUVOYbN27knHPOYf78+cyaNYsHH3wwrfzbbruNRYsWMXXqVP7whz8k973yyiuceuqpzJ07l7lz5/LnP/8ZaLK2/uu//os5c+YwdepUXnrppeR5P/3pT5k0aRInn3wyv/rVr9LezzPPPMPixYuZO3cuCxcuZM2aNQCsWrWK2bNnc9111zFz5kxmzZrFu+++C8DnPvc53nnnHWbPns3FF1/c4h79+te/5swzz+SSSy5h2rRpLF26lHfeeYfzzjuPKVOmcNlllyXvx759+/joRz/KzJkzmTFjBg8//HCX72ECx3G44YYbmDZtGrNmzWLx4sWtfq4+7dO6BWZkZYFFtm7lyBN/ZGzp4Ikdumi+vebOnavt8Zc3q/XBv1e1+frWk3v0tj+2fH3ryT1tnveXN6vbvb66T5PW1taqquqECRP0i1/8oqqqbtu2TYuKirS2tlYbGxv1mGOO0ZUrV6qq6m9/+1sFdMOGDRqLxbSiokLfffddVVWtqanRKVOmJNcB/cY3vqGqqu+9954OHTpU9+3bp4cPH9bZs2fr7t27VVV19+7dOmbMGD18+LBu27ZNAV22bJmqqj766KO6aNEiVVV98803dfTo0bp3715VVf2Xf/kXHTZsmKqqbt68WRcsWKDV1e57f/vtt3XcuHGqqrpy5Uq1LEvXrVunqqrf/va39ZOf/GRyX1uf14MPPqiDBw/WnTt3qqrqeeedp7NmzdIjR45oLBbTGTNm6PPPP6+qqpdcconecsstyfdUVlamGzZs6PI9rK2t1XXr1unkyZPVtm1VVT106FBWn7FPE47jqN3QoNH9B7ThvY1a99YGrdvwtjZUVmrjtm1Zv6pXrNDd3/qO7r7tmzpz9DGqWepCn3IzReQhIKaqn+2J69lOx7Z3lUsvvRRwrbQhQ4awa9cuotEohYWFnHLKKQBccsklSZds06ZNvPvuu8nzACKRCO+++y7Tpk0D4OqrrwZg6tSpVFRUsGbNGizLYtu2bZxzzjnJ80SEzZs3M3z4cIqLizn//PMBWLBgAV/84hcB18I677zzGDVqFADXXnstv/vd7wD461//ypYtW/jIRz6SLDMej7Nv377k9efMmZMsc9myZVnfl8WLFzN27FgA5syZQ3l5OYlRHrNmzWLz5s2cfvrpLF++nDvuuAOA0aNHc95557Fy5Uocx+nSPQSYOHEitm1z9dVXc9pppyXvj0/baMKFPOq5kJEIGAZGqPMuZMO6N8C2O3xezsVMRH4IXASUAzNU9W1v+xTgIWAYcBC4UlUrU877V+CvwGm5qMfZM9vvnnLXs/uTLmYqpQUGnzp5WC6qkUY4HE4um6ZJPB5Pc4eao6oMHz4869ibqiIiqCozZ87k73//e4tjtm/fTigUalGPxPltlX322WenuXYJ3n333YzvLVuan9tWWSKSdm7i/bZV72zuYWlpKf/4xz9YtWoVK1as4Ctf+Qrr1q2jrKws6/cxUGhbwAZ1uly7ro5I5Wacus6FebojZvYk8BGgeeqOe4GfquoU4KfAfYkdIjIXKARWd0N9WmXp9GICZvq2gOlu7ymmTZtGQ0NDUnieeOIJEgPpp06dSmFhIY888kjy+Pfee4+amprkeiL+U1lZyfr16znppJNYtGgRlZWVrFy5Mnnca6+91uZDD3DqqafyzDPPsH//foC0mNmZZ57Js88+yz/+8Y+0MtujpKSEXCUGOP3005Mto3v37uWZZ57h1FNP7fI9BDhw4AANDQ2cffbZ3H777ZSWlrJ169ac1Ls/oIkYWFUVkcpKGisr3RiYYWAOGoRZ1LlWSHUcojt3UbviBar/8L80vvkWmGb7J2Yg55aZqr4E6b+gIjISqADO8DY9DvxEREao6gHgPGAc8B2gQkROVtUXc1235iRaLburNTMbQqEQjz/+ODfccAMFBQWcdtppjB8/HgDLsli2bBk33XQTP/jBD7Btm1GjRiVdv8T5ixcvpqqqivvuu4+RI0cC8NRTT/HlL3+Zm266iWg0ysSJE9t1/WbOnMl//ud/snjxYsrKyjjvvPOS+yZPnsyjjz7K1VdfTUNDA9FolMWLFzN//vx2y5w6dSonnHAC06ZN44knnujsreKee+5JNjKoKrfffjvTp08H6NI9BNi5cyfXXHMN8XiceDzOOeecw4IFCzpd1/5ARgtMDK8VsvMWGIDjWWGRys049fVIQQHhE6YTmjyJ+IED1K1e02FXU9r7te4sIrIdOF9V3/Ysr4dVdXrK/neAy1V1Xcq2cuCWTDEzEbkWuBZg/Pjxc3fsyDpnW79FRKitraW4uOcsSZ/+TasCFgoigUDXynYcYrt3E9lUSWzXB6CKdcxowlMmExg3DjGaHMXI1q00rFvPmXfewZu7P5A2ik3SpxoAVHU7kDH4r6r3A/cDzJs3z5+G3ccnR6gqGo02DSVKEbCuWmAATl09kc2eFVZXh4TDSSustfJDEycSmjiRXV+/NWtfv6fEbCcwRkRMVbVFxASO8bb7dJLusqp9BgZO6ljIiNeRNUcC5lphe4hs2tRkhY0eTeG8uQTGjUWyiIup4yCQlVUGPSRmqrpfRNYDnwAe9f6/4cXLskZELgAumDRpUu4r6eMzAHBSx0KmCVhuQhVOfX1TLCxhhU3/EKHJkzFLshNJJxpFI1HEMGiIxbIeD5jzmJmI3AN8DCgDqoCDqjpdRKbhds0YAhzG7ZqxsTPXmDdvnvoTmvj4ZIeTOpQoEkkKWFdjYAlUtckK27nLs8LKCE+Z0iErTBsaUdvGLCrEHD4cs7gYw7LWquq8bOrRbQ0A3UGKZXZNZWVlu8f7+AxUulvAAJyGBs8Kq8Q5WoeEQ4QmTXJjYSXZpSFzolE0GkXEwBo2FHPwYIyUfoYi0j/FLIFvmfn4tKSlgIERCuVUwFSV+J49NG6qJPb+TtcKKxtFaMoUguPHZW2FOQ2N4DiYhYWYw4dhFhdnPLcjYtanWjN9fHw6RlLADh/CaWwSsFxmowDPCtu8hcimSpyjR5FQiPCHjndjYaUdsMK8WFhg2DDMwaVpVlhXyausGfk21dyTTz7J8ccfz5w5c9i4sfXwYPO8Yd2dY2v9+vUtOo3mC7/+9a8zZt9ozm233caXvvSlHqhRz+NEIsSqDtJYWelaSHv3AHjZKIpzHgurXfU3jvz+DzSsewOjuIiikz/M4I9fROG8ue0KmToOdl09dm0thmURHD+O8PHTCJSNyqmQQZ5ZZqq6DFg2b968a3JVZv2GDdSuWIlTXY1RWsqgpadSOGNGTsq+7777+OY3v8nHP/7xnJSXK9avX8/TTz/NJZdc0ttV8cmSdAus0YuB5d4CA88K27LVtcJqa5FQkPDx0whNmYyZxTSPiX5rGo25VtjwYW4sLGVMcHeQV5ZZrqnfsIHqZX/G8Sw9p7qa6mV/pn7Dhi6XffPNN/Piiy/yla98hVNPPbWF9dWZLK6rV6+moqKC2bNnM336dB5//HF2797N6NGjaWxsTB534YUX8pvf/Ib9+/dz+umnM2PGDGbMmMHNN9/MwYMH+drXvsby5cuZPXs2N954I9B+/rP/+I//YM6cOUybNo21a9dyzTXXMHPmTE466ST27t2bsb7l5eXccsstLFy4kPHjx/Ob3/yGu+66ixNPPJFJkybx4otNI9YefvhhZsyYwcyZM/noRz+aHB8ajUa57rrrmDp1Kqeddhqvvvpq2jW+//3vc+KJJ1JRUcEFF1zQal3yEXUc7NpaIlu20rhpU4oFNiinFhh4VtiePRz929858sQfaVi7DqOwwLPCLqZw/rx2hUxtG/toHc7RoxiBQJMVNmpUtwsZ5Jll1hGqn/0r8b372jwmumtXy/FfsRjVf1pGw9o3Wj3PKhtF6dlntVn2j370I9544w2+9KUvcf7557N9+/Zsq94q3/ve97j55pu54oorUFWqq6sZPHgwS5Ys4be//S1XXXUVO3bs4PXXX+eJJ57gpz/9KRMmTGD58uUAHD58mCFDhvDNb36Tp59+OjlO8siRI1x//fU888wzjB49mj179jB//vxk2umDBw/y4Q9/mO9+97v84Ac/YOnSpaxatYpf/OIX3HDDDfzkJz/h29/+dsY6RyIR/u///o/XXnuNU045he9///u8+uqr/O53v+M//uM/eOmll3j77bf56le/ytq1axk9ejS33norn//85/ntb3/Lfffdx7Zt23j77beJxWJ85CMfoby8HIBHH32UzZs3s2bNGgzD4Oc//zlf/OIXeeyxx7p8r3sTjcWIV1cTP1CFxuJIODcdWTPhNDY2xcI8Kyw0bSrhyZMxB3fECosipkVg5AjM0tIeEa/m5JWY5bzTbGsDWTuRS6knOPXUU/nud7/Ljh07OOOMMzjppJMAuPHGG7n55pu56qqr+PnPf85nPvMZgsEgCxYs4M477+TLX/4yS5Ys4ayzMgvw6tWr281/lhh0XlFRwdixY5k9ezYAc+fO5fnnn2+1zv/8z/+cPK++vj65PnfuXDZv3gzAypUrOffccxk9ejQA1113HbNmzUruu+qqqwgEAgQCAS6//PJkZtynnnqK119/nYqKCsDNr5avs92rKtrYSPzQIezDhwFBCsIYBbmNKyWuFd+7j8imSqLvvw+OgzVyJAWzZhIsn5Bdi6Rtuy2S6mAWFWONHo1RXJw2vrKnySsx60jMrD3LCWDfXfckXcxUjNJShn3qyk7VsTUsy8JxmnKnpbqF2XLTTTdxwQUXsHz5cj7/+c9z5pln8u1vf5tFixZh2zYvv/wyDz30UNIVW7hwIevXr+f555/nkUce4fbbb09LkZ2go/nPOpK7LHGs6T0gqeupedQy5SlL7GsNVeWWW27hM5/5TKvH9HUSrqRdVYVdX4+YJlJU1OJ+5AKnsbEpFlZTgwSDhKZOITxlMubgwe3XNTUWZpq9aoVlYkDHzAYtPRWaxx0CAXd7jikrKyMWiyWtkd/85jcdLmPTpk0cd9xxXHfddXzhC19Iix99/vOf59JLL2XhwoWMGzcOgG3btlFSUsKll17KnXfeydq1a3Ecp0WOsc7mP8sVS5cu5ZlnnknGu37xi19w+umnJ/c98sgjxONxGhoa0u7bhRdeyM9+9jMOHz4MuC7tm2++2SN17ioaixE7UEXjxo1E338fx7YxBw3CKCzMqZCpKrG9+zj69xfdFsnX12KEQxR9eDGDP34RRSfOb1fIkrGw2qMYgSCh8gmEp00lMHJknxEyyDPLLNckWi27qzUzFcuyuPvuuznjjDOYMGECp57accG85557WLlyJcFgkFAoxI9//OPkvksvvZTPfe5z3HDDDcltq1at4o477khahffeey+GYbB06VJ++MMfMmvWLJYsWcI999zTqfxnuWL69Ol897vf5YwzzkBEmDhxIvfd5+buvPbaa3nrrbeYPn06Y8eOZcmSJWzbtg2AK664gqqqKpYsWQI0TUqScFH7GqqKNjQQO3QI58gRmlzJgpxfy2mMEN2yhcbKSpzqGiQQIDRlMqEpk7GGDMmurn3YCstEXo0A8Icztc5LL73E9ddfz4YNG7rFRfHpPAlXMn7gANrQAJaFhMM5/5xUlfj+/W4sbPsONxY2YgShKZPdWFgWmWDVtt2uH45iFhdhDR+OUVTUa7GwfjsCoDv6mfUHrr76ap5//nkefvhhX8j6EBqLET9yxG2VtG0kFMTohlZJJxIhumUrjZsqcaqrO2+FRWKIZRIY0fetsEzklZj5ZKb5/JY+vUfSlTzoupIqYBQUYHQyr31b14kfOEBkYyXRHTvAtjFHDKdo0UKCx5Z33AobVIx1zDG9aoV1FV/MfHxygDoOdk2N60o2NrquZHERRhct5UT6aKeuDqOoiNCMExDHJrKpEvuIZ4VNOs61woYObb+eqmgkisa8WFieWmGZyCsx85Mz+vQ1nGgUO+FKOg5GDl3JyNataRN7OHV1NKx5BQBz+DDXCiufkNVIgLR+YYNKsMbktxWWibwSMz9m5tMXSLqSVVU41TWoCEZBOOeuZP3adRk7cEtBAaXnnZtdPSNRNOb1zh810rXCgsGc1rOvkFdi5uPTm6htN7VK5tCVTLtGLEb0/feJbNmK1mfOGK0NbWeSHghWWCb697vLQ15//XUuu+yyLpVx77338qMf/ShHNeoeTjnlFJ5++ul2jysvL0+OEe0tnGiU2P79NL630U0LDRiDBmEUFOSk9Vgdh+gHH3D07y9y+He/p+6l1Ti1R1t26PYwiopalqGK0+hm1tBolMCokYSnTiVUPgFz0KB+L2TgW2Zs2LCBFStWUF1dTWlpKUuXLmVGN3SazZZ58+a1OlA6Ho9jZdFKdf311+e6WgMOVcWpryd+8CB2dTUYBkY4nNW4xWzLtw8dIrJlK9Ft29HGRnd40cTjCB53LNaIEUS3bWs5Ga5pUlAxu6mc5lbY2DHuKIIBIF7NGXjvOIUNGzawbNmy5NCe6upqli1bxoYcpAACd3zhbbfdxqJFi5g6dSp/+MMfkvsuu+wy5s2bx4wZM/joRz+aHJKzatUq5s1z+wgm0u/cdtttfPjDH+aXv/wlY8aMSabHOffcc5MDwPfv38+YMWOA9MSEmdIGAdTU1PDZz36WE088kZkzZ/KFL3wBO0N8ZtWqVcyaNYtrrrmGGTNmUFFRwT/+8Q8uueQSPvShD3HWWWdRV1cHwNGjR/n0pz/NCSecwAknnMD3vve9ZDnvvPMOJ510EhUVFVx++eVpY1P37NnDxRdfzIknnsiMGTP47//+74z38xvf+AbTpk1j9uzZzJkzhyNHjnT8Q2kHtW3ihw8TqawkunUrTn09RnExZlFRToTMPnqUhrc2UP2nZdQ8/QyRjZuwRo6k+NQlDL7kYooWnkRg5EhEhNDEiRQtWpC0xIyiIooWLSB47LE4jY2eFRYjMGpUkxXWy4O9e5N+a5k9++yz7ea22rVrV4sHOBaL8ac//Ym1a9e2el5ZWRlnn312VvUwDIPVq1ezceNGFi1axMknn8zIkSO5++67k/nMbrnlFr73ve9x++23tzj/4MGDHH/88dx2222A29P/hRde4KKLLmL79u3u2LtYjBUrVnDaaae1OD9T2iCAf/u3f2PJkiX88pe/xHEcLrvsMh544AGuuaZl28o777zDQw89xC9+8Qs+97nPcdZZZ7FmzRrGjh3Lueeey+OPP85nP/tZvvWtb+E4Dhs2bKC2tpaFCxcyc+ZMzjnnHK644gpuvPFGrrrqKtasWcPixYuT5V955ZXceuutfOQjHyEajbJ06VLmz5/PGWeckTzm8OHD/PCHP2T//v0UFBRQW1tLQQ6HATmRiNsqWXUQ1EFCoZy1SjrRKNHtO4hu3UZ8n5uWyho5kvDCkwhOmNBmt4jEZLiQsMIa0Lo6jEGDsMYOG7BWWCbySsxy3TUjkyXS1vbOcPXVVwMwdepUKioqWLNmDRdeeCEPP/wwjz32GNFolLq6OqZMmZLx/HA4nJYR9rTTTmP58uWMGTOGk046CVXllVdeYfny5SxdurTF+a2lDXrqqad49dVXueOOOwCor69n7NixGeswderUZMqfiooKduzYkTw2NZXP8uXLufvuuxERSkpK+MQnPsHy5ctZvHgxb7/9NldccQUACxYsSLrydXV1rFq1igMHmqZQra2t5d13300Ts5KSEqZOncrll1/O2Wefzfnnn8+gLopN0pWsOohdU+3O4l1YkBNxUNsm9sFuolu3Ed25ExwHo6SEgjmzCR5bnnV+MrdFMoJG40jAIlBWhllS0m9bJLtCXolZR7pmZGM53XXXXWSaT6C0tJRPfepTnalimyRS3bz44ov8/Oc/Z/Xq1YwYMYLf/OY33H///RnPKWqWDmbp0qV861vfYuzYsSxduhRVZcWKFaxYsYKvf/3rLc5vLW2QqvLkk08y0fvVb4vmKX+arzd4rWttpfJpLVDuOA4iwmuvvUagjf5SpmmyZs0aXn75ZV544QXmzp3Ls88+y8yZM9utf3PUtps6uEaiELDcXFxdDOarKnZVlRsH274DjUTc6demTiE08VjMYcOyvkZT7/xELGyY2yLpD1drlQFtny5durTFAxQIBDJaOJ3lwQcfBKCyspL169dz0kknceTIEUpLSxk2bBiRSIQHHngg6/ImTJiAaZo89NBDLF26lNNPP51f//rXBAIBxo8f3+L41tIGXXjhhdx+++1JK7SqqiqZjaKznHHGGfzyl79EVamtreV//ud/OP300ykpKeGEE05Ipu959dVXk3HJQYMGcfLJJ6e52Dt37mwRIqitreXAgQMsWbKEb3zjG5xwwgkdbuV0IhFie/fS+N57xHZ94ObRH1TsBva7IBJ2TS0N69+k+n//RM0zzxLZvIXA6NEULz2VwR+/mKIT52MNH57VNZxoFLu2Fo1ECYwYkR4L84WsTfLKMss1CVenO1szQ6EQixcvpqqqivvuu4+RI0dyzjnn8OijjzJt2jTGjh3LvHnzWuS2b4ulS5fy0ksvJTOzFhQUcPLJJ2c8trW0QXfddRf//u//zqxZs9xgcyjEXXfdxbHHHtvp93rrrbfyr//6r8n7d8UVVyQt5IcffphPf/rT3HnnncydO5cFCxYkz3vssce4+eabk+cNGjSIBx54gLKysuQx1dXVXHTRRTQ0NOA4DhUVFXzsYx9rt06qilNX57qStbVuq2RB111JpzFCdPt2Nw7muchWWRkFM04gOGE80gE3MG027+IizFGjBkx3ilySVymAEuTLJMAiQm1tLcXFuZ9Bx6dt1LaJ19Rg79+PE40igQASCnXJulHbJrZrF5EtW4l9sNt1AQcPJnjcRDcOlqH/V5vlxWKoN9u4NXQY5pDBOZ9+Ld/ptymAfHzaw4lEsA8dJn7oIDiKhENdmgxEVYnv209061aiO953J+4oKCA8bRrB447FHDKkQwKZyPVPPI6EwwTGjMEsKclZ/7WBjC9m3Ug+Wr35SJMrWYVdW4sYJtJFV9KurnYD+Vu34dTVgWURHD+e0HHHYpWVdbhsjcfRhkZ3xvHBgwkMHerW0Y+D5QxfzHzyFo3HiVfXYFcdwInEkGCgS62STkMD0e3biWzZhn3wIIgQGD2agorZBMeN6/A8lcluFbE4EghgHTMaq6Qkp/Nd+jThi5lP3uE0NmIfPkL8YBUoSEEYs6Rz+bg0Hif6/k6iW7cR270bVDGHDqVw3lyCE4/tVH5+v1tF75BXYubnMxu4qOO4HVz3H8CuO+q6kp3s/a6OQ3zfPiJbtrnzRsZiGEWFhKd/iODEiVhDBneqjk406vYtMy036eHgwX7n1h7Eb8306bNoPI7T0EC8utod7O0oEgwgwWCnrJz44cNEt24jsnUbWl+PBAIEJowndNxErFGjOlVmi24Vw4b53SpyiN+a6ZO3ONEoTl0d9pEjOHV1KCCW1em+YU59PZGt24hu3ebOFC5CYMwxBOfPJTh2bFa58jPhd6voe7T7SYrIvwG/UtVqEXkEmA/cqKrPdXvtfPo9iSC5ffQo9uHDOI0RRIBgECnqXOLD1ASH8T3uSAJz+HAKT5xP8NjyTouO362ib5PNz9KnVPVOETkVGAl8BrgH8MXMp1Oo4+A0NODU1hI/fATiMVQMjFAQc1DnOhir4xDbs4folq3uwO64jVFcTHjmDEITJ2KWlnS+vn63irwgGzFLpJA4FXhMVVeLiB8Q8OkQyfjXkSPYNTVu/Ms03F75BZ23lDInOJxI8LiJWCNGdFpw/G4V+Uc2YtYgIv8FXAYsFvfb4TfR+LRLi/iXggQ6Hv9qMd3ah44HO05kyzYcLwtsYOwYQhMnEhg7pktuX1q3ihK/W0U+kZWbCdwAfFlV94rIcUDmvM4+A5pETMmuq8tZ/CvjdGuvuS3Z1siRhBecRLC87QSH2eB3q8h/2hUzVd0E3CQixSJSrKpbgO92f9V88oH0+NdhiMe7HP9Klq1K/WtrW51ureScs7pWvp+tol+RTWvmNOARYAagIrIBuEJVN3Z35Xz6JhqPux1Yq6szxL+6nspao1EiW7cR2VTpth5mOqad6dbaLN/vVtEvycbN/DXwY1xBAzd29hCwoLUTfPofTjSKXVeH08X4V1vEDx4ksnETkW3bIR7HHDoUCQbRaLTFsZmmW2uLpm4VNhIO+d0q+iHZiFlAVR9OWX9URG7KZSVE5EPAFwDTq9OnNR+HJvQjkvGvo178KxpBoEvxr4zXicWIbNtOZNMm7IOHwDQJHltOeMoUzOHDsppurc3y4/Gkded3q+jfZCNmb4rIh1X1JQARWQysae1gEfkhcBFQDsxQ1be97VNwLbphwEHgSlWtBFDVd4DrvON+DxQBRzv5nnw6STL+VVNL/MhhiNuoiBv/ynGCyfihQ0Q2VRLZug1iMczBg91OrcdNTAu8J2YmSm3NLKiYndye8X0ku1XEkEAQa7TfrWIgkI2YVQBXichmb30SsF5EXgVQ1RObHf8kcDfwYrPt9wI/VdVHReRy4D4gOTea1yn3s0AVUN/B9+HTSZLxryNe/AtFjET8K7cumMbjRLfvoHHTJuwDVWAYBMvLCU2d3GafsNTp1tosv0W3irF+t4oBRDZi9oWOFJhiwSW3ichIXFFMzB32OPATERmhqge881YCK0Xkx8BsYF1HruuTPcn41+EjOPV1oLgzFOVomrXm2EeO0LipkuiWrWg0ilFaQuH8ea4V1sUuFZDoVhFFTNPvVjGAyaZrxt8AUoWnE4wDPlBV2yvTFpHd3vYDInIKcDEgQABoMe2OiFwLXAtknIXIp3VaxL8ije6PjRf/6g7LRW2b6I4dRDZWEt+/37XCJownNGUK1qiRXZ/Wze9W4dOMbLpmnAT8DndaunEiMg+4VlWvzVUlVHUVsKqdY+4H7gc3BVCurt1faTP+laOZujNhV9cQ2bSJyJataCSCMWgQBXMrCE06LifdH/xuFT6tkY2beSdwDl6vf1V9XUQe6uB1dgJjRMT0rDITOMbbnjV+csa26cn4V9p1bZvozp1ENm4ivnefm2Zn/DjCU6ZgjS7LyeS6frYKn/bIRsyCqvpOsy9ky44/baCq+0VkPfAJ4FHv/xsddVs7MqP5QKGn41+p2LW1bovk5i1oYyNGcREFc2YTmjypU+mmm6O2jdPQiKi63SqG+d0qfFonGzGLiEgx7mOS6BOWuVu2u/8e4GNAGbBcRA6q6nTgeuAhEfkacBi4squVH4ikx78O4UQiCAKh7ot/pV3fcYjt3EXjpk3Ed+9xrbCxYwhNnUJg9OicCKhGoziNUcQyCYz0A/o+2ZGNmH0HN3fZMSLya+Bs4PLWDlbVG4EbM2x/Dzipc9V0GahupsZiOA0N2LW1rvvYQ/GvVOyjR4lUbiZSuRltaMAoLKRg9ixCkyZhFBV2ufw0V7KggOCEcX5A36dDZDUHgIhMBM7CbW18TlU3t3NKt9Lf5wBI9Jdy0+dUo9GIaxdbppv/vodiReo4xD74gMimSmK7PgBwrbApkwmMGZMbK8x3JX3aIKdzAIjIv6vq94GfZ9jmkwOSrmN9PU51DU69N/bRMJBQECPHve/bw6mrJ7J5M5FNlTj19e4M3jNnEJo8KWcjATQaxUn0DfNdSZ8ckI2beSnQXLgybet2+oubqaqu61hfj11Ti11bg6iiiDv7UA7HPnakTrHdu4ls3ORaYapYx4ym8MT5BMaNzY0VlnAlY3GkIExw3Fi3VdJ3JX1yQKtiJiJnAGfixspShau022vVCvncmplIG+3UHiVeU+OOG0QgYOY080RHcRoavFhYJc7ROiQcIjz9Q4QmT8YsyU08znclfXqCtiyzKO5gbwXqUrbvwU/O2C7qOF7W1Xrs6iNoY6PrOnp5v4xw14fxdLpuqsT37qVx4yZi7+90rbCyURRUVBAcPy5nMTnflfTpSVoVM28Y099E5A+JzBe9TV92M5NTpjU04FTXYB896rqOXqtjb7iOzXEaG4ls3uLGwmprkVCQ8PHTCE2ZjFmaG4O7eQdX35X06SmyiZmdKSI7+8K8mX3NzUx0mYjX1OLU1KC27bpOAQujqLBPuFGqSnzffiKbNhHd8T44DtbIkRTMnklwwoTcWWG+K+nTy/jzZnYA94FtwDlah11T7bpQApgmEg5h9CHrw4lEiG7ZSuOmTTjVNUggQGjqFEKTJ2MNGZyz6/iupE9fwZ83sw3UcVzXsa4ep6Yap74+rctEVyfsyDWqSvzAASIbK4nu2AG2jTliOEWLFxEsn4BY2Xzc2V3HdyV9+hp5NW9md8fMVNW1NBoaWnSZMHqpy0Q2ONEo0a3biGzchH3kCAQChCYdR2jKZKyhQ3N2Hd+V9OnL5NW8md0RM0t0mXCHCtW6XSYEsAK92mWiOc0nwg3PmY1VWkLjxk1Et2+HuI05bCiFCxcQOrY8pymik66kYfiupE+fJavhTH2NrgxnSub5qqvDrqlxE/wJiNmzQ4U6QvOJcNOwLELHlrtJD4cPy9k1m7uS1ogRvivp0+PkdDhTvpPsMuENFbLr6lzX0TBc17G4b7qOqTSseyPzRLjBIIMv+iiSQyvJdyV9ehLHcbBtO+0Vj8eJx+NEo1FCoVDWgel+KWZONIo2NnpdJqrBcVBAAoE+02WiPdRxiO/fT3TbDpy6zPO7aDSaMyHTWAynMeK7kj5dRlVbCJRt28RiMaLRKLFYLPmyM/1Ii2AYBoZhYJpm1l/CvBKz1hoAmrpMHMWursGJel0mLAsJh/PGNUq0Rka3bSe643131m7LBNPMaJl1dCLcTNfzWyV9skFVM1pRqcKUEKp4PJ7RYBARTNNMClU4HM6mO1PWcbBssmaco6p/ybbA7iS1ASBeV+d2laipdf+L12Ui2Pe6TLSFqmJXHSS6fTvR7Ttw6uvBMAiMHUuwfALBsWOI7tzZpYlwW1wz1ZUsLSUwfFiXXclE7NVxHFckO/hyHKfFq/l2VU17GEzTTL4S2xK/6qn/M23zcWnNzWtuQcVisRbnqmryvideoVCIcDjcK/e4rYHmY1V1F3Ar8Bdv2/25nMiks0Tq63nzb38DEcQ0MSzLFTLvi5t8gRsPSywbRnJZRJKxsuSXHlcQDQBvvxgGiY8lcUxqGZm2N//f/DgFnEOHie14n9iOHThHj7oCdswxFMydQ3DcuLTWyI5MhJsmEs3WnVgMu9GdDMQYOgSjpASxLBpjMbfFshURaU1gmi+n1iFBW1/q5selvppvS6zH4/FWBbO1ayX2pR6TEMDm//uDSGZy8+LxeAsLKhaL4ThOWr0T9yj1XliWRTAY7DPvrzXasswe8ua7HC4inwTW4g5l6nXUcQgXFhIIh5MPhAJ4DzApy473RQfXInF3pT94mrJMs2UFr68ZSYFKe3BUk6KVZg9721P/W/X1hPcfIFxVhdXQiIoQGVxKZMpkIsOHQSDgPiR79qSLrQiEwxiLFyLeezqsirNrF06q8Hh1TeJdV2Mx16oLuplpjYICqK11X6Qe3raQpL5Mr9W3+XH5QHMRTDzsfVkk23LzsolDQUs3r6CgoE+NWukqbQ00XyoiBbgiNgn4Z2CqiDwBLFfVe3uojm2SfNjclV6tSyaMo3VYe/dg7d2HWVeHAvbQoTQeeyyxUSPduStVCZMuoKSsJwTZpukBSrxf0zCwTLOFFYhtu7FDBBk8BKu0xJ2lqQ/eo56mu8S3O0WyrfeSKpqhUN8aVteTtOVmPgG8AERU9ZvetvXAt4ClPVK7PEXq6wns3Ye1Zy/m0aOugA0ZTOP4acRHjUSbzeKdav10mXgcJxpDRLAGD8YoLsbIYQdan9bpTpFMlO/TOm25md8ATgMmiEglsAV3oHkxcHcP1K0FidbM8WPG9Mbl20QaGrD27iOwdy9mjeu+2YNLaZw21RWw7pyo1huGhW0jwSCBkSPc0Qt9sAOwT8fxRSw72nIzNwAbRORSVV0oIpOB54FPA/cBJ/RQHVPrtAxYdsK0aX0iBZA0NnoCtg+zuhoAu6SEximTiZeNQnMwd2SbJF1JkMJCrNJS35X06Rds3bqVdevWUVpa2rKVqxWy6Wf2MICqVorIYVX9bKdr2A+QSARr336svXsxDx9BAHtQMZHJk4iVjUILuz7tWrv4rqRPP2br1q2sXr261YaM1mhXzFT15ymr/9TBevULJBr1BGwf5qFDroAVFRE9biKx0WVoFzuvZoXvSvr0MxzHIRqNEo1GiUQiRCIRotEor7zySoeFDDo4AkBVd3T4CvlKLIa1fz+BPZ6AqeIUFhKdeCzxsjKcnuqY67uSPn2chCilClJiOdO2xHI0Gs1pPfJqOFO3E4+7ArZ3H2bVQVfACgqIlk8gXjYKZ9Cgnuv+EYvhxOK+K0lT/KSuro6ioiIqKiqYmKHDsE/XsG07TWiyEaRIJJJxdEAqoVCIYDCYHB1QUlJCKBRKvhL7EsvPPfcc9fWZxyO3hS9mcRvrwAGsvXuxqg4ijoMTDhObMJ5Y2SickpKeEzDflWxB8/hJXV0dq1evBhgQgtYZIbdtu1UrqS0LKh6Pt1qmiKSJTjgcprS0tFVBSl3uqBcxd+7c7omZ9UtsG+tAFdbefVgHDrgCFgoSGzuGWFkZzuDSnu2Am3AlFaS4GKtkUAtXsj9ZJ631Zs/0evXVV1t8qW3b5pVXXkk+fKn3KVOfvfb2d/T4zi539LwPPviAt956K03IX375ZXbu3ElxcXGrYtWeKKUKTWFhIUOGDGlTkEKhEAFvdEpPkPher1u3rkPn5WVyxhOmTdPf3nsvgY703XIczKoqtzPr/gOIbeMEA8RHjSJeNgp7yJCeH0GQcCUNA7O0pFVXMlPrjmmaLFq0qNOC1pqYNN8ej8fbFZzEOe0dm1p2Pn7v+hKJ3v7NRae9dcuy8ireOn369K11dXXHZXNsXllmHe406ziYBw8R2LvXFbB4HA0EiJWVER/tCVhPD/1IuJJxGwmHXFeysDAt7U5qK080GuW1117LaJ2sWbOGqqqqNoWktawIXRWTxDi/1l6WZREKhdo8JpvXihUraGhoaHH9wsJCzj33XO+WZhhT28rA90z7O3p8a+fmuoyVK1fSGpdffnleiVJPkFdillWnWcfBPHwYa88+Avv3I7EYalnER44gVlaGPWxotwuYqhK3baKxGFEvW0E0EiHS2EgsbhMzhLgI0ZSBwqlN1G25CanEYjG2bNnSqhAEAoEui0lrr54a/zdv3ryMVuncuXMp6okuMb1IUVERdXV1Gbf7QtaSvBKzBIdqavjL3/7O7LFjGT/9Q6DqCtjefVj79mFEY6hpegI2Cnv48A4LmO04RBMpU+LxdGHy1mPNtyeWs7B8AoEAwWAw+SouLk5bT7gIwWCQ1atX09jY2KKMoqIiLr744g69r3wjNX7SH+KFHaGioiKjkFdUVPRirfoueSlmAHUor+zaiXX4MJPjcYxIFDUM4iNGECkbRXT4MOKqrsAcPdq2+GQQKNtx2ry+YRgELYtgIEAwECAUDDKosJBgIEDAsgiKEDBNguEw4cGDCQ8aRCgcJhgMEggEOmTZzJ8/f0B/qSdOnDggxKs5A1nIO0PeihlAHHi57ijbgiEiRQVERIhVHyZ68ACxdlw1wbOOLCv5v7S4OG09KUyeYCWXLSuZzyuN1FbJQcVYg1q2SnaGTF/qOXPm+F/qAcBAEPLUbMO2bRO3bZymH+6sH568FjNwBa066IpNUWtClEGYLNPMWdxBYzHUa5W0hgzBLC7u1Ozhjiq2o9g22KrE40rUVqJxxRg0hrlLvIYPz4PdfiCCZQqmAZYhmIYk1w1xkzoahvdf/OwLfR33gQYHRR2w1dum4LjZR0n+9RyHZO67ZBkp5TXbpmhyY/MoiOIlKW3tfE3flnqMNt+X4Vj3vdnYtoPj2DiODZra4CGYAYtAIIQVKCAQDGFZAcQKDxwxK0Y4d/Hinr+w1yqpto0RCmFlaJVseYr7ZY076oqVDVFbicUdYnGwHfcrJYg7pwFgChiGELRcQUrFUfeLFLchFne89UQJXpbclBS4hiFYBq74mQaWAZYpGYXPaH6xXiA1LXeiO0fznPN9VaAdR1F1f5RU3R8qR8Fx3M857rhxWfeHy93maCLJKCR0JyWhddNaK29ZWllObmvrVrUoOsNSs/Ml8fK2u5+Va1U5tu1+l1POMS2LUDhEKFhMwOsmYlkBTNPCbKXLiCptx3tSyGsxs4BZY8f27EVt2+1agdvBNVBS4k4e7H0QtqMpL4jGlZjtEI1D3E5pmsf9Ehi4wmEaEDA7ltzP8L5NXgJrb2tr2Us98cMV0EjMxknUI/FLKu5yom6mJ35WwuIzwZSE5SfuMZ4IdlRUUkUqVaxSM6yqarKLRzgcxjTNtDz2mXqIJ841Pcu7eYrqjpKwjJKC5LhWjO14lrTt/Th5n7ftNAlY8ockIU7S9OkY3j1O/HiYGX6s0ul90XYcxxUqx8Z2vM8LI2HTYZoWgWCQQGEhwWAIKxDE9EIypml1ewt43opZEdLUmtkDuK5kDMSA0sFoQRGOYdIQV6KRuCdaqaa6+y12v7CuMIQD0muZvV3hSa4lFLBVEuLnONBoO2i0SQwT1oNq8jnF9ATZFEVQDFFEHQQHcDBIWH5NqZ6DwSDhcDitZTe1n5rZTiigrZEEmWYXSljGrkBJ0p3DELxae9vdZdtxLaw0MdLU9+5N/iHu/RVP6APSe59zV2iyrLz/qkkJdT9jk0AwSKigkEAw5HX9sTzB6n6xao+8FLOhJSWcs+QjHRsBkC3eF941/R2cxijxuE3cDBIrHEI8EEbjgtQqKvE2XME8/Dan0EL8SPmyOw5qu5ZVQt3iCjFVHMT7YgcwrRBm0HUnMCzUdOcANQwLK2AipiAWbrcZExzLQDwrz1ZQB0xDM8447ziKreA4goOFLSaOCTaKLYojoAEH4iCOYjhKPO7WXx0btR1Ubex4HDsex7FjOHYc2/uPOu78CaKINFkfIu6MXYYY7v+USUj6OqpOU8yq2SgMVbeFPhAMUVBY4IlVimVldawFvjfISzHrKurFJ2zvgbBtJW67lpWdMEdiMUDRwiKk1G2VNAwhLH1y3pQukQjONk0pZ7sWi+cWJSwwwzSxAkFCnjVlBYLuL3OiI61pYhjtN6w0BbWVhqgXF9NmQWlpWjfEdcENQ7BtJea0dONUm+orJBpAXDfOTMQGgxbZfuUdx0EdBzsRA/LuSzwWw47HicejyenbbDueqHILDMNIil9iubuET1VxbLupzqmRec9yDASChMOeGxgMejErzw3M84QGfULMROTDuOm4Q8ARVf3XLhWoXkzDi2vEbTfg6v53WrTaJFxBsW2CdsydQXzYECgscpfzFLfJ2/HchoRVpWkiBSCGELACbkuSFcAKBglYAYxmIpWrX2YRwWwR62sdJyWIbpkQtLrfEjIMAwwDM4tHJO3HwG76b9tx4vFYUvDsWIx4NOLFmlqUgojh/hgYkhLnM1pex7OumvpCuj6wGK5YhcKuG+i67U1uYMbuRP2InIuZiPwQuAgoB2ao6tve9inAQ8Aw4CBwpapWAqjqS8BL3nF/EpFiVT3a1nVsRyGuXlBWiXliZSesq2SF3H+JX+mA2cy0UtcKE9txZ00aPATCBbQTje11kpPvprh+QFqgWSHZYhQKF2AFvF9iy0oKlWmYff4X2Z3Iud0wX6/hNjRY7u9eOynnVNX7zDxBSjaAxN3xtfFY0vqLRCM4jt0k3CIErCDBUJhAyHUDk5aVZWVlFfdnusMyexJ39qYXm22/F/ipqj4qIpfjTopyWuoBInIe8G77QgZVtXGsoJcUThKtQ00tQ+36grYDsSiCokWD0KJiCAT7TKgr8QuccBlSW+YVMMTACgQIhtxWo0AwlOYuJERqIH+5+yIigpiJH5D2k20mRE8Q//Nsh5yLmWdlNc/RNBKoAM7wNj0O/ERERqjqAe+YTwHlqvrVdq8BWAYEA51we+IxJBZ3g9Glg9FeciWbN3PjdUxUz/8zrYDbchRoaua2khZV77cc+fQMCXfTp316KmY2DvhAVW0AVbVFZLe3/YCInA98G3haRO4Fbk2IXAIRuRa4FmDUyFEdu7o6EI0hjutK6oihEAx3qyvZspnbi6h7FlbGZu7UYKz/Bfbx6RB9ogFAVZ8G2uz9qqr3A/cDTJsyre2UFAlsB2IRN36UcCWDwa5W16tPejO342ha/6tkM3dRAYFAiECzYKwvVj4+rbN1f4R1OxopHTk2p/Nm5oKdwBgRMT2rzASO8bZnTSI545jR7SRnjMXcRIymCaVD3bksO+hKtmzmTres3JzooWQzdyAYSu/t3MeD6j75QeKhros4FIUMKiaEmTgy1NvV6la27o+wenM9dtYDmVx6RMxUdb+IrAc+ATzq/X+juSuZRTnLgGXTpmRIztjclRzStivZkWZuV6wGVjN3X2IgPtDQ8qGuizis3uzOWtTe+3fH6KYMY0sblpWynjxGU/Z18fjE9VP2px1P5vITy7uPxDosZNA9XTPuAT4GlAHLReSgqk4HrgceEpGvAYeBK3NyQdv2OrgCRcVJVzI5SDlme31+nKZBvG5FvWbuEIFQ2G/m7qNs3tfImi0N6Q90ZT11EYexQ4NutofEAwUpy5m2a7NjUh6iVrYnj0/dlrY9ZZtXDinbneZlpl4n9fwMdT941Ca1lxG4kZOXNtWzbkdDipBkEJee+HCyJDHcyxC8URWp28QbCtY01K0zQgbd05p5I3Bjhu3vASd1pexUN1NsG62rIy6CU1iEEwqjpgHxKBJ3xc0KBJIuoJtaJGFVWXnXzN2XrRPHywASd7wOyl4n5Vhy2dufXG7qF5i6P5ayP9HJufnDDO5Qs3U7Glm3o2X23Z5EEiMOEmMzaRqnmdif+gA3P7ZpuWlspzt6QRCDjO8dXKEqKw20EITUMaKpA9mbi0XyPNLXM5WVtt7s+ER5acfT8nod5YnXqqmLdFzR+kQDQLYk3Mzjpx5/TVQEhg4lUFziWlZe/iPXFXTjVvkkVm3RFXcjgWrqKIgmcYklBcgVo5jdbMRE6vGZ9nvZIjpCYniSZUoyDZFlCAVBIWAaySwdlglv74q0Ws6SqUXNBCFb0WgSnsRoiOSDmXZMa2V27iHtKK091EUhgw9P6b/zH1RMCPfdmFmuscIFjJv/kc71M8sDHMdNyhi33f+vbWto8cHaDqzZUk/VUTspKrFm1k+atdQJwXFFpUlsLNPN/GGFDE9sXDEKpOy3zKaUQYFUsUoeQ8aB462x7UCs1Qe6fERuWqb7KpkeatNwt/dnEj/QHbW880rMEm7mhPLj+pzVlRgDGo277lLaK64tt9lKLNOxdvaWTsyGzfvcbLOBFDEJmkJh0EgTo0B7YpPc5u7vC8kZYeA+0JD+UPfF8EJ3MnFkiIkjQ3x9/66t2Z6Tl5MAHzNplt5893PMLS/o0ger6lozWYlPKyKVsKCytXwSSRgTr6DlikgwRYgClqQds2ZLPY2xlp9TUUi4eP7gTr//fKEvxwt9upfpxx2zta7mYP+bBDiV+qiyurKe2gabkaWBZkLkCk/SVYs3uWypYhRvmag0I4kYT8BKiI5QEDAYFHYzOKQKTyCDGAVMkts64mIlsB1txTop6HBZ+UjiV9rHpy3ySswSbubo42YCbqvW+p0R2NkySCziWTlJgYFwQBgUNtoRnnTxCZjuRCG9yUB2N3x8siWvxCzRmnnMpFlpnWbPmlHcwnXrb7MR+daJj0/b5JWYZaIoZFBW2n4qFR8fn/5NXvdtGCitWj4+Pu2TV2ImIheIyP0AhUFh0aRC3/Xy8fEB8kzMVHWZql47rNjinypKfSHz8fFJkldi5uPj49Maed8A4OPj0/skMoLQItNIemYRaJZFRMSdWDlBYiaexHIH8MXMx6cfkzmVUZPAACmph5oEBqUpZVbqpKTJgr3/SfHxZnc3JDlnqSmCeGNxDXFntDelaWZ7wzu+aaB/+mB/Q6Ch7kh1tu81r8QsdWymT/6Q+kClJvNLTx7oPRMZnhlJXUlFEuW3/BFX7+FqtjHt3ExlKy0n9spUfnKD96CnldVWvZLnpRTRrNjEfQAQT3ESAkOifilikrx26jZ1a9VcYAwRDNNN1WMYCWFxJ0o2cCeXNwyjmcCkZx8xmq13ZyYRTc6h2D55JWaJTrMzZs9tmWnWJ2dkJz7eQywZHvS0stwvu2m4E/8a3uB3d90dXWEabjcbwzDS8mg1z5tlSJOrknaNFguumGWqS8bzUlYyjVRu85reUvb16kjZ6Qc0zx3WWwLTV8krMfPJTFqKYtLFJ1WMUsUnITLp5bj/DU90WohPUoDSxadF0j9pnuSvNx6ogfMQ+7j4YtYLdNjy0QzuT0pq5EziYxlGMkaRH+Lj49M1fDHLgmzFR72AT2sxn45aPqZhtHCzMqdK9sXHx6dfilnm4HJT83ByW0ogNWPQN7G/gzEf3/Lx8el58krMkq2Zxx5HQ8yhMS5NItSsNcr0WnBML5NqQnwMT3ysFgHnloJj+OLj45M35GWm2blz5+nfV7/iTkphZLKEfPHx8ekPiMhaVZ2XzbF5ZZklEIGikD8Jr4+PTxP+2EwfH59+gS9mPj4+/QJfzHx8fPoFvpj5+Pj0C/KyNVNEaoGNvV2PXmI4UNXblegl/Pc+8JiqqoOyOTAvWzOBjdk21/Y3ROR1/70PPAbqexeR17M91nczfXx8+gW+mPn4+PQL8lXM7u/tCvQi/nsfmAzU9571+87LBgAfHx+f5uSrZebj4+OTRl6JmYgME5FnRGSjiLwlIn8UkRG9Xa+eQESeFJE3ReQNEXlRRGb3dp16GhH5uoioiJzQ23XpKURku4i8JyLrvddZvV2nnkJEwiLycxGpFJENiQnAWyPfumYo8H1VXQUgIj8Abgeu7s1K9RBXqWo1gIj8P+ABoKJ3q9RziEgFsAB4v7fr0gtcrKpv93YleoHvA43AFFVVERnV1sF5ZZmp6qGEkHmsASb0UnV6lISQeZQCWc9ak++ISAj4KXADmecc8elniEgxcCVwq3qBfVXd19Y5eSVmqYiIAfwL8FRv16WnEJFfisj7wHeAq3q7Pj3IN4FHVXVbb1ekl3jMC6v8TEQG93ZleojjgIPA10XkdRFZJSIfbuuEvBUz4MfAUeAnvV2RnkJVP6uq44H/BH7Q2/XpCURkITAf+Flv16WXOFlVZ+HeA2HgfN8tYCLwhjfy4SvAH0WkpLUT8rJrhoj8EJgJXKCqkd6uT28gIg3AWFU92Nt16U5E5KvAjUDU2zQW2Ad8WlWf67WK9QIiMgN4SlWP7e26dDciMhzYAwQTbqaIvANcqaoZhzjlnWUmIt8B5gL/NFCETESKRWRcyvoFwCHv1a9R1dtV9RhVLVfVcmAXcNZAEDIRKRKRUm9ZgEuB9b1aqR5CVauAlcAZACIyBRgJbG7tnLxqzRSR6bgu1iZgtTfRyDZV/WivVqz7KQJ+LyJFgI0rYhdoPprVPh1hFPAHETEBE3gHtxFkoHA98ICI3AHEgCtU9UhrB+elm+nj4+PTnLxzM318fHwy4YuZj49Pv8AXsz6CiMwUkT96yy+IyPhWjpsnIo91ovxyEekTmUpF5FMi8kQXyzhFRM7MVZ28MtXrrNltiMhtIhLszmsMVHwx6zucBrwgIgXASFXNOGxHVV9X1ct6tmp9klOAnIpZD/F1wBezbsAXs15GRG4QkTW4X/JrgHXACBFZIyLnZTj+lEQq4YS1JSLf8Qagb0ztJS0inxORzSLyIs3Gr4rIuSLysoisFZH/E5EFKeW/KSIPisg6EXlVRD6Uct5VIvKKd94LIjLV2/4pEXlORH4rIv/wyi7z9gVF5D6vfi8AJzary79711knIstSzrtNRB73kgu8JyJ/FpFCr7/V9cCV3uDrr2a4T9eKyLve/rdEZJq3Pc36ymCNfUlEVnt1vcg7plBEfi8i73j35nfe9mdE5OKUsj4mIs95y1+XpgHib4jIYBH5qXfoam/7YBEpEXdkx6tePe/2Wi8Rt9f7HeImFtgpIl8SkU949dsuIh9vq34DDlX1X738wu3ZvcVb/jpweRvHngK87i2X445VPN9bvwx42VueCewGRnnrPwOqvOXjgP8DSrz16cD7KeUrsMRbvyrleicDfwZC3vo5Kdf7FHAYGOet/wL4jrf8eeA5IAAUAq8DT3j7LsdNwGd46/8CPOYt3wZUAoO9e/QccE3Kvh+2cZ+qU+oSAgq9ZQWKU45LrnvLX/OWp+IOpxkJfBRYnnLOEO//2cDKlO0rgP8HDAFqgQJv+yDAauX6v8TtcgCucfF4yntcBfzW234MUJ9yT08EdnnLGes30F551c+sHzMJ2OItzwU68st6VFWf9pbXAHd4y6cAf9amwbn3A5d4y2fhCtrfxe2rB2BJU1aCzar6N2/5EeB+cYeRXADMAl7xzhPcBzfBy6q6M6UuZ3jLpwIPqWoMiInIo0DCgrwQmAes88q0cIUowV/V61skIq949c6GF4AHReRP3n3YmuV5vwJQ1Y0isg43U8ebwDTPslqFK+gAfwV+JCLH44rUcUDis9gIPCoizwJPq2ptK9e7EDhRRL7orRfidgxO8HtVdYDdInIQ+F9v+1pgjIiE26jfgMIXs15ERMqBJ3EFISQibwJTgAki8pqqfjaLYlJHQdg0faaS4VhS9j2rqldmqNPx7Zz3gKp+rZX9jZ2sy7dV9YEsyyxoo6xUPoY7nvE0YKWIXK+qf/HKMMDNl9VOGQKoqm717stSXGv0v0Vkhqo2egKS6Mh6n6raXtkLgMXe9deKyNmq+lYr1/inNsS2+ftvxK2UnRD/turXzvvrV/gxs15EVber6mxgOa7VdBXwF1WdlaWQtcVK4FwRGemtp8bMngPOFndEBQAiMj9l/yQROdlb/iSwQVVrgGW4caqx3jmmiMzNoi4rgCtExBK3geOTKfueAm4QkSFemSERmZVFmTW4qZBaICIWMFFVX1XV2733O8fbvQVX5BLvrTmf9sqYDMzGtULHAraqPgncDIwAhnrHPwT8E/DPuC4jIjIIGKGqf1PVrwNvA4mEkrXN6v0U8NWUONlwEenQ2Mt26jdg8C2zvsFC3ID2jbjC1mVU9S0R+W/gZRHZS4rroaqVInI58CtPXILAy8Br3iHrgU+IyF241sCV3nl/F5H/Ap7yHr4g8Htcl6ct7seN4f0D14X6G3CsV+Yj4g4q/ptnaRi48b032ynzf3EFcj3wP55oJTCBX4ubLscBdgKJRoKbgftEZA9NLmEqERF5GXfS3etUdb+InAPc7tXPBL6rqru9+td6rmSBqh7wyijFHYZU4L2fdcAfvX134LZaN+CGAm7CTUL4pogorqV9E9CRdEczWqvfQMIfzuSThoicghtYH3ATznYGzwp8CzcT8GvtHe/Tffhupo9PJxGRC3Hd1ud8Iet9fMvMx8enX+BbZj4+Pv0CX8x8fHz6Bb6Y+fj49At8MfPx8ekX+GLm4+PTL/DFzMfHp1/w/wHtmqPpPwwzWQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 324x201.6 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(4.5, 2.8))\n",
    "names = ['independent models', 'full system model', 'pairwise models']\n",
    "colors = ['#80aee7ff', '#e68082ff', 'grey']\n",
    "data = [sampled_independently, sampled_combinatorics, sampled_pairwise]\n",
    "for _name, _data, _col in zip(names, data, colors):\n",
    "    ax.plot(n_systems, _data.mean(axis=1), 'o-', label=_name, color=_col)\n",
    "    l, h = confidence_interval(_data.T)\n",
    "    ax.fill_between(n_systems, l, h, alpha=.25, color=_col)\n",
    "\n",
    "\n",
    "ax.set_ylabel('# steps')\n",
    "ax.set_xlabel('# independent subsystems')\n",
    "ax.set_xticks(n_systems)\n",
    "ax.set_xlim(n_systems[0] - .1, n_systems[-1] + .1)\n",
    "\n",
    "ax.semilogy()\n",
    "ax.legend(loc='upper left', frameon=False)\n",
    "ax.set_xlim(2, 6)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}