mmca_sociopatterns_interactive.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "07e9d4e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "from itertools import combinations\n",
    "import json\n",
    "import numpy as np\n",
    "import random as rd\n",
    "from collections import defaultdict\n",
    "\n",
    "%matplotlib inline \n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "20f95073",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_simplicial_complex_from_cliques(cliques):\n",
    "    \n",
    "    G = nx.Graph()\n",
    "    triangles_list = set() #will contain list of triangles (2-simplices)\n",
    "    \n",
    "    for c in cliques:\n",
    "        d = len(c)\n",
    "        \n",
    "        if d==2:\n",
    "            i, j = c\n",
    "            G.add_edge(i, j)\n",
    "        \n",
    "        elif d==3:\n",
    "            #adding the triangle as a sorted tuple (so that we don't get both ijk and jik for example)\n",
    "            triangles_list.add(tuple(sorted(c)))\n",
    "            #adding the single links\n",
    "            for i, j in combinations(c, 2):\n",
    "                G.add_edge(i, j)\n",
    "            \n",
    "        else: #d>3, but I only consider up to dimension 3\n",
    "            #adding the triangles\n",
    "            for i, j, k in combinations(c, 3):\n",
    "                triangles_list.add(tuple(sorted([i,j,k])))\n",
    "\n",
    "            #adding the single links\n",
    "            for i, j in combinations(c, 2):\n",
    "                G.add_edge(i, j)\n",
    "                \n",
    "    if nx.is_connected(G)==False:\n",
    "        print('not connected')\n",
    "                \n",
    "    #Creating a dictionary of neighbors\n",
    "    node_neighbors_dict = {}\n",
    "    for n in G.nodes():\n",
    "        node_neighbors_dict[n] = G[n].keys()\n",
    "        \n",
    "    #converting the triangle set of tuples into a triangle list of lists\n",
    "    triangles_list = [list(tri) for tri in triangles_list]\n",
    "\n",
    "    return node_neighbors_dict, triangles_list\n",
    "\n",
    "\n",
    "def import_sociopattern_simcomp_SCM(dataset_dir, dataset, n_minutes, thr):\n",
    "    filename = dataset_dir+'random_'+str(n_minutes)+'_'+str(thr)+'min_cliques_'+dataset+'.json'\n",
    "    SCM_cliques_list = json.load(open(filename,'r'))\n",
    "    \n",
    "    #considering one realization of the SCM\n",
    "    realization_number = rd.choice(range(len(SCM_cliques_list)))\n",
    "    cliques = SCM_cliques_list[realization_number] \n",
    "    node_neighbors_dict, triangles_list = create_simplicial_complex_from_cliques(cliques)\n",
    "    \n",
    "    N = len(node_neighbors_dict.keys())\n",
    "    avg_k1 = 1.*sum([len(v) for v in node_neighbors_dict.values()])/N\n",
    "    avg_k2 = 3.*len(triangles_list)/N \n",
    "    #ass = nx.degree_assortativity_coefficient(facet_list_to_graph(cliques))\n",
    "\n",
    "    #print N, avg_k1, avg_k2, \"Assortativity %.2f\"%ass\n",
    "    \n",
    "    return node_neighbors_dict, triangles_list, avg_k1, avg_k2\n",
    "\n",
    "def get_tri_neighbors_dict(triangles_list):\n",
    "    tri_neighbors_dict = defaultdict(list)\n",
    "    for i, j, k in triangles_list:\n",
    "        tri_neighbors_dict[i].append((j,k))\n",
    "        tri_neighbors_dict[j].append((i,k))\n",
    "        tri_neighbors_dict[k].append((i,j))\n",
    "    return tri_neighbors_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "246c053b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def markovChain(beta_A, beta_D_A, mu_A, beta_B, beta_D_B, mu_B, epsilon_A, epsilon_B, node_neighbors_dict, tri_neighbors_dict, NSteps, i0_A, i0_B, i0_AB):\n",
    "    N = len(node_neighbors_dict)\n",
    "    \n",
    "    #initial seeds\n",
    "    p_A = np.zeros(N)\n",
    "    p_A[rd.sample(range(N),int(N*i0_A))] = 1\n",
    "    p_A_new = np.copy(p_A)\n",
    "    \n",
    "    p_B = np.zeros(N)\n",
    "    p_B[rd.sample(range(N),int(N*i0_B))] = 1\n",
    "    p_B_new = np.copy(p_B)\n",
    "    \n",
    "    p_AB = np.zeros(N)\n",
    "    p_AB[rd.sample(range(N),int(N*i0_AB))] = 1\n",
    "    p_AB_new = np.copy(p_AB)\n",
    "    \n",
    "    q_BA = 1\n",
    "    q_AB = 1\n",
    "    q_S = 1\n",
    "    q_A = 1\n",
    "    q_B = 1\n",
    "    \n",
    "    pATime = [np.mean(p_A)]\n",
    "    pBTime = [np.mean(p_B)]\n",
    "    pABTime = [np.mean(p_AB)]\n",
    "    \n",
    "    for k in range(0,NSteps):\n",
    "        for i in range(0,N):\n",
    "            #updating the infections - simple contagion\n",
    "            for j in node_neighbors_dict[i]:\n",
    "                q_BA *= 1-epsilon_B*beta_A*(p_A[j]+p_AB[j])\n",
    "                q_AB *= 1-epsilon_A*beta_B*(p_B[j]+p_AB[j])\n",
    "                q_S *= 1-beta_A*(p_A[j]+p_AB[j])-beta_B*(p_B[j]+p_AB[j])+beta_A*beta_B*p_AB[j]*p_AB[j]\n",
    "                q_A *= 1-beta_A*(p_A[j]+p_AB[j])\n",
    "                q_B *= 1-beta_B*(p_B[j]+p_AB[j])\n",
    "                \n",
    "            #updating the infections - triangle contagion\n",
    "            for j, k in tri_neighbors_dict[i]:\n",
    "                q_BA *= 1-epsilon_B*beta_D_A*(p_A[j]+p_AB[j])*(p_A[k]+p_AB[k])\n",
    "                q_AB *= 1-epsilon_A*beta_D_B*(p_B[j]+p_AB[j])*(p_B[k]+p_AB[k])\n",
    "                q_S *= 1-beta_D_A*(p_A[j]+p_AB[j])*(p_A[k]+p_AB[k])-beta_D_A*(p_B[j]+p_AB[j])*(p_B[k]+p_AB[k])+beta_D_A*beta_D_B*p_AB[j]*p_AB[j]*p_AB[k]*p_AB[k]\n",
    "                q_A *= 1-beta_D_A*(p_A[j]+p_AB[j])*(p_A[k]+p_AB[k])\n",
    "                q_B *= 1-beta_D_B*(p_B[j]+p_AB[j])*(p_B[k]+p_AB[k])\n",
    "            \n",
    "            #calculating f_A and f_B\n",
    "            if ((1-q_A)*(1-0.5*(1-q_B))) == 0.0:\n",
    "                f_A = 0.\n",
    "            else:\n",
    "                f_A=(1-q_A)*(1-0.5*(1-q_B))/((1-q_A)*(1-0.5*(1-q_B))+(1-q_B)*(1-0.5*(1-q_A)))\n",
    "            if ((1-q_B)*(1-0.5*(1-q_A))) == 0.0:\n",
    "                f_B = 0.\n",
    "            else:            \n",
    "                f_B=(1-q_B)*(1-0.5*(1-q_A))/((1-q_A)*(1-0.5*(1-q_B))+(1-q_B)*(1-0.5*(1-q_A)))\n",
    "            \n",
    "            #Updating the vectors\n",
    "            p_AB_new[i] = p_B[i]*(1-mu_B)*(1-q_BA)+p_A[i]*(1-mu_A)*(1-q_AB)+p_AB[i]*(1-mu_A)*(1-mu_B)\n",
    "            p_A_new[i] = p_B[i]*mu_B*(1-q_BA)+p_A[i]*(1-mu_A)*q_AB+p_AB[i]*mu_B*(1-mu_A)+(1-p_A[i]-p_B[i]-p_AB[i])*f_A*(1-q_S)\n",
    "            p_B_new[i] = p_A[i]*mu_B*(1-q_AB)+p_B[i]*(1-mu_B)*q_BA+p_AB[i]*mu_A*(1-mu_B)+(1-p_A[i]-p_B[i]-p_AB[i])*f_B*(1-q_S)\n",
    "            \n",
    "            #resetting the parameters\n",
    "            q_BA = 1\n",
    "            q_AB = 1\n",
    "            q_S = 1\n",
    "            q_A = 1\n",
    "            q_B = 1\n",
    "        \n",
    "        p_A = np.copy(p_A_new)\n",
    "        p_B = np.copy(p_B_new)\n",
    "        p_AB = np.copy(p_AB_new)\n",
    "        \n",
    "        pATime.append(np.mean(p_A))\n",
    "        pBTime.append(np.mean(p_B))\n",
    "        pABTime.append(np.mean(p_AB))\n",
    "    \n",
    "    return pATime, pBTime, pABTime    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "316f8420",
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_dir = '.../sociopatterns_thr_data_random/'\n",
    "dataset = 'Thiers13'\n",
    "n_minutes = 5 #aggregation time\n",
    "thr = 0.80 #fraction of removed cliques (0.80: retaining the 20% most weighted)\n",
    "\n",
    "node_neighbors_dict, triangles_list, avg_k1, avg_k2 = import_sociopattern_simcomp_SCM(dataset_dir, dataset, n_minutes, thr)\n",
    "tri_neighbors_dict = get_tri_neighbors_dict(triangles_list)\n",
    "\n",
    "mu_A = 0.05 #rec. prob. A\n",
    "mu_B = 0.05 #rec. prob. A\n",
    "\n",
    "lambda1_A = 1.2 #rescaled 1-simplex inf. prob. A\n",
    "lambdaD_A = 0.8 #desired rescaled 2-simplex inf. prob. A\n",
    "\n",
    "lambda1_B = 1.2 #rescaled 1-simplex inf. prob. B\n",
    "lambdaD_B = 0.8 #desired rescaled 2-simplex inf. prob. B\n",
    "\n",
    "i0_A = 0.4 #initial conditions (% of infectious nodes at t=0) A\n",
    "i0_B = 0.4 #initial conditions (% of infectious nodes at t=0) B\n",
    "i0_AB = 0.2 #initial conditions (% of infectious nodes at t=0) AB\n",
    "\n",
    "beta_A = 1.*(mu_A/avg_k1)*lambda1_A #1-simplex inf. prob. A\n",
    "beta_D_A = 1.*(mu_A/avg_k2)*lambdaD_A #2-simplex inf. prob. A\n",
    "\n",
    "beta_B = 1.*(mu_B/avg_k1)*lambda1_B #1-simplex inf. prob. A\n",
    "beta_D_B = 1.*(mu_B/avg_k2)*lambdaD_B #2-simplex inf. prob. A\n",
    "\n",
    "epsilon_A = 1.5\n",
    "epsilon_B = 2.\n",
    "\n",
    "NSteps = 300"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e832e2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "pt_A, pt_B, pt_AB = markovChain(beta_A, beta_D_A, mu_A, beta_B, beta_D_B, mu_B, epsilon_A, epsilon_B, node_neighbors_dict, tri_neighbors_dict, NSteps, i0_A, i0_B, i0_AB)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "17826b85",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'density of infected')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7CElEQVR4nO3dd3xUZfb48c/JpEMIJPQSSujSq4DriiIiCoplwd5Z17auu34Fd9ey/nQtu7quqyL2CrooCEpVRNRFINJ7SQKElgKkkD5zfn/MgCFOkglkMinn7eu+Zm557pzrJXPmPve5zyOqijHGGFNaUKADMMYYUzNZgjDGGOOVJQhjjDFeWYIwxhjjlSUIY4wxXgUHOoCq1LRpU+3QoUOgwzDGmFrjp59+SlfVZt7W1akE0aFDBxISEgIdhjHG1BoisqesdVbFZIwxxitLEMYYY7yyBGGMMcarOnUPwpuioiJSUlLIz88PdCinJTw8nLZt2xISEhLoUIwx9UydTxApKSlERUXRoUMHRCTQ4VSKqpKRkUFKSgodO3YMdDjGmHqmzlcx5efnExsbW+uSA4CIEBsbW2uvfowxtVudTxBArUwOJ9Tm2I0xtVudr2IyxtReLnVR7CqmyFV08rXIWeK9Zyp2FeNS188TrlPnPZOq4lSne73r1O1OrvO8P/mfZ0gE5eehEU6sP/G+3PXl7MPb+pPbKKfso2SZ0iJDIrm1161n8H/aO0sQ1WT27NlcccUVbN26le7duwc6HGOqlEtdHC86TnZhNlmFWae8nnifW5RLfnE++c588orzyCvOc8+XWJZfnE+hs/DkF79TnYE+tBpH+GWtQmxErCWI2mzGjBmcc845zJw5k8ceeyzQ4RjjE1UlqzCLlOwU9ufsJy0vjfS89JNTRl6Ge8rPqPDLPDI4kvDgcCKCI4gIjiDcEU54cDiNwxufnI8IjiDUEUpIUAjBQcGEOEIICQr5eb7U+xOvjiAHQRKEQ9yvJyZBfrlM3MtOvAYRRFBQEEGcuu7EF/Ep75FTqn1PvPe2vuQXeXn787wpcx+BrGa2BFENcnJy+OGHH/jmm28YP368JQhT4+QV55GYmcjOozvZfWw3+7L3sT9nPynZKeQU5ZyybbAEExMRQ9OIpjSLbEaP2B7EhsfSOKwxUaFRNAptRKOwRkSFRp2cGgQ3wBHkCNDRmdNVrxLE4/M2s+VAVpXus2frRjw67qxyt5kzZw5jxoyha9euxMTEsGbNGgYMGFClcRjjq9yiXDZnbGZD2gY2pW9i57Gd7M3ae7KeOzQolLZRbWnTsA39m/enTcM2tI1qS9uGbWke2ZzosGiCpF60b6n36lWCCJQZM2Zw//33AzBp0iRmzJhhCcJUm7ziPNYcXsOKAytYdWgVO47uOFkdFBcVR7eYblzS8RI6N+lM58adiYuKs1/7BqhnCaKiX/r+kJGRwdKlS9m0aRMigtPpRER49tlnrQmr8ZvU3FSW7FnCN/u+Yc3hNRS5iggJCqF/8/7c1vs2+jbrS++mvWkS3iTQoZoarF4liECYNWsWN954I6+99trJZb/+9a/5/vvv+dWvfhXAyExdk5GXwcLkhSxOXsza1LUoSnx0PNd0v4bhrYczoMUAIoIjAh2mqUUsQfjZjBkzmDJlyinLrrzySj766CNLEOaMudTFyoMrmbVjFkv3LaXYVUyXJl24q99djG4/mk6NOwU6RFOLWYLws2XLlv1i2X333Vf9gZg6pchZxJdJX/LOpnfYnbmb6LBorul+DVd0voLOTToHOjxTR1iCMKYWKXIW8d8d/+XNjW+SmpdK1yZdeeqcpxjdYTRhjrBAh2fqGEsQxtQCLnWxIGkB/1n7H1JyUhjYYiBPjHiCYa2HWWMH4zeWIIyp4bYf2c7fVvyNDekb6NakG6+OepURrUdYYjB+59enXURkjIhsF5FdIjLFy/rLRGSDiKwTkQQROcfXssbUdblFufwz4Z9M/GIiKTkpPHnOk3wy7hPOaXOOJQdTLfx2BSEiDuBl4EIgBVgtInNVdUuJzb4G5qqqikgf4BOgu49ljamzNqdv5qHvHmJP1h6u7HIlfxj4B6LDogMdlqln/FnFNATYpaqJACIyE7gMOPklr6olO3lpACf7wq2wrDF1kUtdvLP5HV5a8xKxEbG8ddFbDG45ONBhmXrKnwmiDbCvxHwKMLT0RiIyAfg70By4pDJlPeUnA5MB4uLizjhof3A4HPTu3RtVxeFw8J///Ifhw4cHOixTw2QXZvPQ8of4bv93XNj+Qh4d9qhdNZiA8meC8FZJ+ovRLlR1NjBbRM4FngBG+VrWU346MB1g0KBB3kfTCLCIiAjWrVsHwKJFi5g6dSrffvttYIMyNcq+rH3cs/Qe9mbt5S9D/8Jvuv3G7jOYgPNngkgB2pWYbwscKGtjVV0uIvEi0rSyZWuTrKwsmjSx/m/MzxIOJXD/svsBmD56ulUpmRrDnwliNdBFRDoC+4FJwLUlNxCRzsBuz03qAUAokAEcq6jsaVkwBQ5tPOPdnKJlb7j46XI3ycvLo1+/fuTn53Pw4EGWLl1atTGYWmt5ynIeWPYArRu25uXzX6Zdo3YVFzKmmvgtQahqsYjcAywCHMBbqrpZRO70rJ8GXAncKCJFQB4wUd2Drnot669Y/a1kFdOKFSu48cYbT/buauqvJXuW8H/L/48ujbsw/cLpNA5vHOiQjDmFXx+UU9X5wPxSy6aVeP8M8IyvZc9YBb/0q8OwYcNIT08nLS2N5s2bBzocEyBfJn7Jw98/TJ+mfXhl1CtEhUYFOiRjfsGGhapm27Ztw+l0EhsbG+hQTIAsT1nOn7//MwNbDOS1C1+z5GBqLOtqoxqcuAcB7kHg3333XRwOG7GrPlqXuo4/Lvsj3WK68dL5LxEZEhnokIwpkyWIauB0OgMdgqkBdh7dyV1f30WLBi14ddSrNAhpEOiQjCmXVTEZUw2O5B/h7q/vJtwRzmsXvkZMeEygQzKmQnYFYYyfFbuKefDbB8nIy+C9se/RpmGbQIdkjE8sQRjjZ/9M+CerDq3iyXOe5KzYswIdjjE+syomY/xo3u55fLD1A67rcR3j48cHOhxjKsUShDF+sidrD0/8+AQDmg/gj4P+GOhwjKk0SxDG+EGRq4ip300lOCiYZ859hpCgkECHZEylWYKoJrNnz0ZE2LZtGwDJyclERETQr18/+vbty/Dhw9m+fXuAozRV5fUNr7MxfSOPDHuElg1aBjocY06LJYhqMmPGDM455xxmzpx5cll8fDzr1q1j/fr13HTTTTz11FMBjNBUlfVp65m+YTrjOo1jTIcxgQ7HmNNmCaIa5OTk8MMPP/Dmm2+ekiBKsm7A64YCZwF//v7PtIhswdShUwMdjjFnpF41c31m1TNsO7KtSvfZPaY7Dw15qNxt5syZw5gxY+jatSsxMTGsWbOGmJgYdu/eTb9+/cjOziY3N5eVK1dWaWym+r2+4XX2ZO1h+oXTrY8lU+vZFUQ1mDFjBpMmTQJg0qRJzJgxA/i5imn37t3861//YvLkyYEM05yhxMxE3tz0Jpd0uoRhrYcFOhxjzli9uoKo6Je+P2RkZLB06dKT4z84nU5EhLvuuuuU7caPH88tt9xS7fGZqqGqPLHiCSKDI3lw0IOBDseYKmFXEH42a9YsbrzxRvbs2UNycjL79u2jY8eOpKSknLLd999/T3x8fICiNGfq892fk3A4gQcGPkBshHXlbuqGenUFEQgzZsxgypQppyy78soreeqpp07eg1BVQkNDeeONNwIUpTkT2YXZPJ/wPP2b92dClwmBDseYKmMJws+WLVv2i2X33Xcf9913X/UHY/zijY1vcKzgGFOGTCFI7KLc1B32r9mYM3Ag5wAfbPmAcfHj6BnbM9DhGFOlLEEYcwZeXPMiIsK9/e8NdCjGVLl6kSBUNdAhnLbaHHtdtyl9E/OT5nNjzxutOw1TJ9X5BBEeHk5GRkat/KJVVTIyMggPDw90KKYUVeWfCf8kJjyG23rfFuhwjPELv96kFpExwIuAA3hDVZ8utf464MTDCTnA71R1vWddMpANOIFiVR10OjG0bduWlJQU0tLSTu8gAiw8PJy2bdsGOgxTyspDK0k4nMDDQx+2saVNneW3BCEiDuBl4EIgBVgtInNVdUuJzZKAX6vqURG5GJgODC2xfqSqpp9JHCEhIXTs2PFMdmHMKVSVV9e9SvPI5lzZ5cpAh2OM3/izimkIsEtVE1W1EJgJXFZyA1X9n6oe9cz+CNhPZVPj/XjwR9akruGO3ncQ6ggNdDjG+I0/E0QbYF+J+RTPsrLcBiwoMa/AYhH5SUTK7KRIRCaLSIKIJNTWaiRTe6gqr653Xz1c0eWKQIdjjF/5M0GIl2Ve7xSLyEjcCaJkZ0kjVHUAcDFwt4ic662sqk5X1UGqOqhZs2ZnGrMx5frx4I+sTV1rVw+mXvBngkgB2pWYbwscKL2RiPQB3gAuU9WME8tV9YDnNRWYjbvKypiAsasHU9/4M0GsBrqISEcRCQUmAXNLbiAiccBnwA2quqPE8gYiEnXiPTAa2OTHWI2p0Lq0daxNXcutvW61qwdTL/itFZOqFovIPcAi3M1c31LVzSJyp2f9NOARIBZ4RUTg5+asLYDZnmXBwEequtBfsRrji7c3vU10WDQTOluHfKZ+8OtzEKo6H5hfatm0Eu9vB273Ui4R6OvP2IypjKTMJJbtW8bkPpOJDIkMdDjGVIs6/yS1MVXh/S3vExIUwqTukwIdijHVxhKEMRXIyMvg812fMy5+HE0jmgY6HGOqjSUIYyowc/tMCl2F3HTWTYEOxZhqVeY9CBF5iTKeWwBQVRvxxtR5+cX5zNw2k/PanUfHaOuyxdQv5V1BJAA/AeHAAGCnZ+qHuwM9Y+q8hckLOVZwjBt63BDoUIypdmVeQajquwAicjPuTvOKPPPTgMXVEp0xATZz20zio+MZ3HJwoEMxptr5cg+iNRBVYr6hZ5kxddrGtI1sztjMxO4T8TyTY0y94stzEE8Da0XkG8/8r4HH/BaRMTXEzO0ziQyOZFyncYEOxZiAqDBBqOrbIrKAn8dpmKKqh/wbljGBdTT/KAuTFjKhywQahjYMdDjGBESFVUzivrYeBfRV1c+BUBGxjvNMnTZ712wKXYVM6mYPxpn6y5d7EK8Aw4BrPPPZuEeKM6ZOcrqcfLL9Ewa1GETnJp0DHY4xAeNLghiqqncD+QCeEeCsK0tTZ608uJL9OfuZ2G1ioEMxJqB8uUld5BlfWgFEpBng8mtUxgTQpzs/pXFYY86PO9/nMi6XkpVfxNHcIo4cL+To8UJyCorJL3K6p2IXeYVO8oudFDsVlyqq4FL1TO7xJlwuUNzzNYXWpFjKfna3+tWgUKLCg3n8sl5Vvl9fEsS/cQ/Y01xEngSuAv5a5ZEYUwMczT/K0n1LmdRtktcxH/IKnWw5mMnGlEyS0o+z90gue4/ksu9oHoXFFf9uCg0OIiRICAoSgkQIEggSQUQQ4ed5sKa1ZahJ/1tqSiwxkf6p1PGlFdOHIvITcAHuYUQvV9WtfonGmAD7IvELil3FJ0eMKyh28lPyUb7dmcZ3O9LZdijr5K/7qLBg4mIj6doiigt6tKBFo3BiGoTQJDKUJpGhRIUHEx7iIDzEQUSIg7DgIIKCasg3ijE+qDBBiMj7qnoDsM3LMmPqDFXls52f0atpL7KzmzL1m418sf4A2QXFhDiEQe1juOf8LvRuE03vNtG0aBRmv/JNneZLFdNZJWc89yMG+iccYwJnbep6dh3bRfTxa7jiu/8RHhLE2F6tuLh3K4bFx9IwzK/jaxlT45TXm+tU4GEgQkSyTiwGCoHp1RCbMdVCVZm34SCP/e8VNDyEyMKBPH1FDy7p04qo8JBAh2dMwJTXWd/fgb+LyN9VdWo1xmRMtdlxOJupn23kp72HadRtDUObnc8bN4+2qiNj8O05iFUiEn1iRkQai8jl/gvJGP9zupRXlu3i0n9/T1L6ca4/PwuVAu4adI0lB2M8fEkQj6pq5okZVT0GPOq3iIzxs4ycAm5+exXPLtzOqJ7NWfyHc0lnBW0atqF/8/6BDs+YGsOXBOFtG5/u1onIGBHZLiK7RGSKl/XXicgGz/Q/Eenra1ljTsfmA5lc+tL3rEw6wjNX9uaV6waiQVmsOLiCsR3HEiQ2Cq8xJ/jy15AgIs+LSLyIdBKRF3CPNFcuT2unl4GLgZ7ANSLSs9RmScCvVbUP8ASem98+ljWmUpbvSOM301YA8NnvhjNxcBwA85Pm41IX4+KtW29jSvIlQdyLu+XSx8AnQB5wtw/lhgC7VDVRVQuBmcBlJTdQ1f95+nYC+BFo62tZYypjwcaD3PrOatrFRDL7rhH0anPythpfJH5Br9heNua0MaX48iT1cWCKiDRU1ZxK7LsNsK/EfAo/jynhzW3AgtMsa0yZFmw8yD0z1tKvXWPevmUwjUo0Xd15dCfbjmxjyhCrxTSmNF/GgxguIluALZ75viLyig/79tYUxGv3ViIyEneCeOg0yk4WkQQRSUhLS/MhLFOffL318Mnk8O6tQ05JDuC+enCIgzEdxgQoQmNqLl+qmF4ALgIyAFR1PXCuD+VSgHYl5tsCB0pvJCJ9gDeAy1Q1ozJlPfFMV9VBqjqoWbNmPoRl6os1e49y90drOKt1I969dcgvnoR2qYsvE79kRJsRxEbEBihKY2oun5psqOq+UoucPhRbDXQRkY4iEgpMAuaW3EBE4oDPgBtUdUdlyhpTnuT049z2zmpaNArnrZsHe+0mY/Wh1RzOPWxjThtTBl+aq+4TkeGAer6s7wMq7M1VVYtF5B5gEeAA3lLVzSJyp2f9NOARIBZ4xfNwUrHnasBr2dM4PlMPHS8oZvL7CSjw3q1DaNowzOt2XyR+QYOQBpzX7rxqjc+Y2sKXBHEn8CLuG8cpwGJ8a8WEqs4H5pdaNq3E+9uB230ta0xFVJX/+3QDu1JzeO/WobSPbeB1u/zifJbsWcKF7S8kPDi8mqM0pnYor7O+Z1T1IWCkql5XjTEZc9qmL0/kyw0HmXJxd87p0rTM7ZanLOd40XEu6XRJNUZnTO1S3j2IsSISAlhHfaZW+DExg2cWbuOS3q347bmdyt12YfJCYsNjGdxicDVFZ0ztU14V00IgHWjg6e5bcDc1FUBVtVE1xGeMTzLzinjg43W0j23As1f1KbfDveNFx1mespwJnSfgCHJUY5TG1C5lXkGo6oOqGg18qaqNVDWq5Gs1xmhMhR75fBOHswt4YWI/GlQwsM+yfcsocBYwpqM9+2BMeSps5qqq1sWFqdE+X7efz9cd4PcXdKFfu8YVbr8weSHNI5tbz63GVMCXJ6mvEJGdIpIpIlkikl1ihDljAurAsTz+MmcTA+Iac9d58RVun1WYxQ/7f+CiDhdZz63GVMCXZq7PAuNUtcJnH4ypTqrKX+ZsotipvDCxH8GOir/wv9n7DUWuIutawxgf+PIT6rAlB1MTfbHhIEu3pfLH0V3LfN6htAXJC2jTsA29m/b2c3TG1H6+XEEkiMjHwByg4MRCVf3MX0EZU5FjuYU8Pm8zfdpGc8sI37rpPpZ/jJUHVnLDWTfYsKLG+MCXBNEIyAVGl1imuPtQMiYgnpq/laO5Rbx361AcQb592X+19yuKtdiql4zxkS/jQdxSHYEY46sVuzP4JCGF350XT8/Wvre4Xpi8kLioOHrE9PBjdMbUHeV1tfF/qvqsiLyEl7EYVPU+v0ZmjBfFThePzt1E2yYR/P6CLj6XS89LZ/Wh1dze+3arXjLGR+VdQZy4MZ1QHYEY44sPftzDjsM5vHbDQMJDfH8K+qs9X+FSl1UvGVMJZSYIVZ3neX23+sIxpmwZOQU8v2QHv+rSlNE9W1Sq7IKkBcRHx9Olie9XHcbUd77cpDamRvjH4u3kFjp5dFzPSlUTHT5+mLWpa/ldv9/5MboazFkERXlQXADF+SVeS04Fp27jKv55chaBywmuohLLSqx3edY7S6xXl2fSEu9dgJZaXmr9KducWF/GPk7QkjXg6sNyylhe2f3oaWxb3vIz0KAZ3LOqavZVgiUIUytsTMlk5up93DaiI52bR1Wq7OI9i1G0dlYvqUJ+JuQdgVzPlHcEcjMg7xgUZENhNhTkQGGO+7X0Mmdh1cUTFAKOEAgKPnVylJoPcgACElRiKjWPQFAQSEj525yyTE5dflKJ91W2nDKWe9m+qj7zdIU1PPN9eFHheBAicrWq/tcvn26MD1SVR+duIrZBKPeNqnwV0cLkhXRr0o2O0b49L1Etigsh5xBkHYTsEtOJ+ZzUn5OCljXCr0BoQwiLcn9BhDZ0v0Z2OHU+tCGEREBwOASHQXCE59UzH1Jq3hEGjtBSX/oh7i9zU6+UdwUxVkT+gns8CEsQJmA+X3eANXuP8exVfWgUHlKpsgdyDrAhbQO/H/B7P0VXBlXIOQxHk93TkaSf3x9NdieH0hyhENXKPTXrBpGxEBnjfo2IKfG+ift9WLR9aRu/svEgTI2WX+Tk2YXb6NWmEVcNaFvp8ouSFwFwUYeLqjo0N1XI3Adp2yF1K6Rt80zb3dU7JwlEt4UmHaDLKIiOg0atfk4IjVq7v/itCa6pQcprxfQg8KCIfG5dfptAefd/yRzIzOcfv+lLkI9PTJe0IGkBvWJ70S6q3ZkH4yyG9B1wYK17OrjOnRRKJoKGLd2//vtfD7GdoUlHd1Jo3M5dfWNMLeLLk9SXiUgL4MTYjCtVNc2/YRkDR48X8p9vdjGyWzOGx5c9vnRZ9mbtZeuRrfxp0J9OL4CcNNi7wj3tXwOHNkBRrntdaENo1dedCJp190zd3FU/xtQRFSYIEbka+AewDHf10ksi8qCqzvJzbKae+883uzheUMyUi0+va4yFyQsBGN1+dAVbehxPh91LIfl7d1JI3+FeHhwOrfrBgJugdX/3FNvZ6v9NnedLM9e/AINVNRVARJoBXwEVJggRGQO8CDiAN1T16VLruwNvAwOAP6vqP0qsSwayASdQrKqDfDkgUzfsO5LLeyuSuWpgW7q1rFyz1hMWJS+iX7N+tGrYyvsGLpe7qmjnYti1xH2VgLpv/sadDf2uhfYj3MkhOPS0j8WY2sqXBBF0Ijl4ZODbSHQO4GXgQiAFWC0ic1V1S4nNjgD3AZeXsZuRqpruQ4ymjnlu0XYcQcIDF3Y7rfKJmYnsOLqDhwY/dOoKlwv2/QgbZ8HWuXA8DRBoOxhGPgydR7mrjoJ878bDmLrKlwSxUEQWATM88xOB+T6UGwLsUtVEABGZCVwGnEwQnsSTKiKXVCpqU6dtSDnG3PUHuGdkZ1pGh5/WPhYlLUIQRnfwVC+l7YA178KmzyD7gPtZgG5joNslEH8+NIitwiMwpm7w5Sb1gyJyBXAO7nsQ01V1tg/7bgPsKzGfAgytRGwKLBYRBV5T1eneNhKRycBkgLi4uErs3tREqspT87cS0yCU3/6602nvZ1HyIgY070/z5BWw+g1IWu5+2KvLhdDrCeg6xm9PnxpTV/jU1YZn9LjKDhDkrU1iZToeGaGqB0SkObBERLap6nIvsU0HpgMMGjSoijo2MYHyzfZUfkw8wuPjzyKqkg/FnbDz8Hp2Z+7m4RwnrJzjfubggkeh/w3QsFnVBmxMHebPvphSgJKNz9sCB3wtrKoHPK+pIjIbd5XVLxKEqTuKnS7+Pn8bHZs24Nqhp3E1WJADK6excOPrBDUI4cKGHWDU8+6rBrunYEyl+TNBrAa6iEhHYD8wCbjWl4Ii0gD3zfFsz/vRwN/8FqmpET5dk8LO1BxevW4AIY5KNCF1OWHdh7D0/6E5h1ncsTODm3Si6XjrIcaYM+FLa6RLRaTSDb5VtRi4B1iEe/ChT1R1s4jcKSJ3evbdUkRSgAeAv4hIiog0AloA34vIemAV8KWqLqxsDKb2yC0s5vklO+gf15gxvVr6XnDvj/DauTD3Xmjcnu2/eYtkChnd7Wr/BWtMPeHLFcQk4EUR+RR4W1W3VlTgBFWdT6kWT6o6rcT7Q7irnkrLAvr6+jmm9nvzuyQOZxXwn2sH+DbWQ34WfP24+wZ0dDu46m04awIL17yIQxxc2P5C/wdtTB3nSyum6z2/6q8B3va0KnobmKGq2f4O0NR96TkFTPt2N6N7tmBwBx+6qti9FObc7e4W++y7YOSfIawhqsqi5EUMbTWUJuFN/B+4MXWcT1VHqpoFfArMBFoBE4A1InKvH2Mz9cS/v95JfrGLhy7uXv6GxYWw5BF4fwKEN4Lbv4Ixfz/ZXHVLxhZSclL813OrMfWML30xjQduAeKB94EhnpZFkbjvLbzk3xBNXZaYlsNHK/cyaXA74puV81zC0WT47y1wYA0MuhVGPwmhkadssjB5IcESzAVxF/g3aGPqCV/uQVwFvFD6GQRVzRWRW/0Tlqkvnl24nbDgIO4f1bXsjRK/hf/e5B6H+DfvQ8/xv9jkRPXSsNbDiA6L9mPExtQfvlQxHSydHETkGQBV/dovUZl64ac9R1i4+RC//XU8zaK8jJWgCj9Oc1cpNWwBd3zjNTkAbEjfwMHjB616yZgq5EuC8NYc5OKqDsTUL6rKk19upXlUGLf/ystY0S4nfPkALHzI3S3G7V9BbHyZ+1uYtJCQoBDOjzvfj1EbU7+UWcUkIr8D7gLiRWRDiVVRwA/+DszUbYs2H2LN3mM8fUVvIkNL/TMsyoNZt8H2L2HE7+GCx8ode8GlLhbvWcyINiOICj29rsGNMb9U3j2Ij4AFwN+BKSWWZ6vqEb9GZeq0IqeLZxZup0vzhlw1sNRjMLlHYMYk2LcKLn4Ohk6ucH/rUteRmpvKAwMf8FPExtRP5SUIVdVkEbm79AoRibEkYU7XjFV7SUo/zps3DSK4ZJcax9PhvcvcI7ld/Q6cdblP+5ufNJ9wRzjntTvPH+EaU29VdAVxKfAT7l5YSz7eqsDp98Vs6q3s/CJe/GonZ3eK4fzuzX9ekZMK7453N2e99hOIH+nT/oqcRSxKXsTIdiNpENLAP0EbU0+VmSBU9VLPq5c7iMacnte+TSTjeCFvj+3xc5ca2Yfg3XGQmQLX/Rc6/srn/f1w4AeOFRzjkk425pQxVc2XzvpGeHpURUSuF5HnRcRG5jGVdigznze+T2R839b0advYvTDrALw91v16/aeVSg4AXyZ+SeOwxgxvM7zqAzamnvOlmeurQK6I9AX+D9iD+4lqYyrln4u343LBgxd5xpnOSXNfOeSkwvWfQfvKfckfLzrOsn3LuKjDRYQEnd7gQsaYsvmSIIpVVXGPJ/2iqr6Iu6mrMT7bmJLJrDUp3DyiA+1iIiHvqPsBuKwDcP0siKvMaLRuX+/9mnxnvlUvGeMnvnS1kS0iU4HrgXNFxAHYzzXjM1Xl8XmbiYkM5Z7zO7tHfvvwakjfDtd+DHFnn9Z+v9j9BW0atqFfs35VG7AxBvDtCmIiUADc5hm/oQ3wnF+jMnXKvA0HSdhzlD9d1I1GDifMvAb2r4Gr3oL403vyOT0vnZWHVjK241jfxo8wxlSaL+NBHAKeLzG/F3jPn0GZuiOv0MnT87fSs1UjfjOgNcy6CZKWw4Tp0GPcae93QdICXOri0k6XVmG0xpiSfGnFdIWI7BSRTBHJEpFsEcmqjuBM7Td9eSIHMvN59NIeOBY+CNu+gIufhb4Tz2i/XyZ+SY+YHnRqbI/jGOMvvlQxPQuMV9VoVW2kqlGq2sjfgZna78CxPF79dhdje7dk6P63IeEtGHE/DP3tGe03KTOJzRmb7ea0MX7mS4I4XJlxqI054ekF23ApPNF+Ayz9f9BnEox67Iz3O2/3PIIkiDEdxpx5kMaYMvnSiilBRD4G5uC+WQ2Aqn7mr6BM7ffDrnTmrj/A8/1Tif36j9BpJIx/Cc7whrLT5WTu7rkMaz2MFg1aVFG0xhhvfEkQjYBcYHSJZQpYgjBeFRQ7+evnmxgdvZ8Jux6FFmfBxPchOPSM973y4EoO5x7mT4P/VAWRGmPK40srpltOd+ciMgZ4EXAAb6jq06XWdwfeBgYAf1bVf/ha1tRcry9PpDg9kf80+jvSoClcNwvCqubZytm7ZhMdFs357WxgIGP8zZdWTF1F5GsR2eSZ7yMif/GhnAN4Gffocz2Ba0SkZ6nNjgD3Af84jbKmBtp3JJcPl/7ErIbPERqEuwuNqKqpCsosyGTp3qWM7TiWUMeZX40YY8rny03q14GpQBGAqm4AJvlQbgiwS1UTVbUQmIm7u46TVDVVVVef2HdlypqaR1V5ck4CrzmepZkedXfb3bRLle1/QdICCl2FXN758irbpzGmbL4kiEhVXVVqWbEP5doA+0rMp3iW+cLnsiIyWUQSRCQhLS3Nx90bf1iyaT9XJz1CL0lCrnoL2g2u0v3P2TWHbk260SOmR5Xu1xjjnS8JIl1E4nHfmEZErgIO+lDOW3MV9TEun8uq6nRVHaSqg5o1a+bj7k1Vy8orJHf2/VzgWIuO/Qd0H1ul+99xdAebMzZzeefLrWsNY6qJL62Y7gamA91FZD+QhLvjvoqkAO1KzLcFDvgY15mUNQHw4ztTudy1hMN976bFkNuqfP+f7fyM4KBgezjOmGrkSyumRGCUZ9CgIFXN9nHfq4EuItIR2I/7vsW11VDWVLPtC6cx+vAbbGx6Mb0vf7LK959XnMfcXXO5sP2FNAlvUuX7N8Z4V2aCEJEHylgOgKo+7239CapaLCL3AItwN1V9S1U3i8idnvXTRKQlkID7WQuXiNwP9FTVLG9lK3twxv/yti4h/seHSXD0pdftb5/xg3DeLExaSHZRNhO7nVn/TcaYyinvCuJEw/VuwGBgrmd+HLDcl52r6nxgfqll00q8P4S7+sinsqaGObieoP/ewE5XGxzXvU94eIRfPuaT7Z/QuXFnBjQf4Jf9G2O8KzNBqOrjACKyGBhwompJRB4D/lst0Zma69heCt+7knRnJAv6/ocHurT3y8dsztjMpoxNTB0y1W5OG1PNfLlJHQcUlpgvBDr4JRpTO+Qewfn+FRTk5fJw5N+ZNu4cv33UJ9s/ISI4gnHxpz92hDHm9PiSIN4HVonIbNxNTScA7/o1KlNzFebCzGvRI8ncUTiFP908johQh18+KqswiwVJCxjbcSxRoTYMujHVzZdWTE+KyALgV55Ft6jqWv+GZWokZxH892Z074/8vvBehowcz6AOMX77uHm755FXnMfV3a7222cYY8rmyxUEqroGWOPnWExN5nLB53fDzkU8wR0caDuGF8/v7LePc7qcfLDlA/o268tZsWf57XOMMWXz5UlqU9+pwqKpsOFjZkTdzMeuUfxrYj+CHf775/PNvm9IyUnhprNu8ttnGGPKZwnCVGz5c7ByGgmtrmFq2oX87bJetI9t4NePfG/Le7Rp2Ma69TYmgCxBmPKteh2+eZKDHSZwddIlTBocx5UDvT66UmU2pG1gbeparu9xPY4g/9wAN8ZUzKd7EKaeWvM+zP8TeZ0u4tKkifRo1ZDHxvv/fsB7W94jKiSKCV0m+P2zjDFlsysI4926j2DuvTg7nc8NWb+j0BXEK9cNIDzEv7/o9+fsZ8meJVzV9SoahPi3GssYUz5LEOaXNnwCc+5CO/2ah0OnkpCSy3NX96VDU/9/Yb+/5X0E4doe1jejMYFmCcKcatOnMPu30OEc3m73FB+vS+MPo7oypldLv390Wm4as3bM4tJOl9Kygf8/zxhTPksQ5meb58Cnd0C7s/l28H/4f4uTGdu7Jff68XmHkt7Z/A7FrmIm95lcLZ9njCmfJQjjtnEWzLoV2g5ix6g3ueeT7XRv2Yh/XN2XoCD/d5KXnpfOJ9s/4ZJOlxDXKM7vn2eMqZglCAM/vQuf3g5xZ3Pw0g+48YOtNAgL5o2bBhEZWj0N3d7b/B6FrkLu6H1HtXyeMaZiliDquxWvwLz7oPMoMq+cwU0fbeV4QTHv3DqY1o39M75DaUfyjzBz+0wu7ngxHaI7VMtnGmMqZgmivlKFb591d6HRYzz5V73P5BlbSEo/zms3DqR7y0bVFso7m94hvzjf7j0YU8PYg3L1kcsJ8x+EhDeh77UUXPIvfvvBelYlH+HFSf0ZHt+02kLZn7OfD7d+yLj4cXSK7lRtn2uMqZgliPqmMBc+vQ22z4cR91M08q/c89E6vt2RxjNX9mZ839bVGs6/1/wbEeHe/vdW6+caYypmCaI+yUmDGRPhwFoY+w+KB97GHz5ex5Ith3l8/FlMHFy9rYc2pW9iftJ87uh9hz33YEwNZAmivkjbDh9NhOxDMPEDCjqP4fcfrWXh5kNMvbg7Nw3vUK3hqCrPrX6OmPAYbut9W7V+tjHGN369SS0iY0Rku4jsEpEpXtaLiPzbs36DiAwosS5ZRDaKyDoRSfBnnHXe9gXw+gVQmAM3zSM/fgyT3/uJhZsP8cilPfntr+OrPaSl+5ayJnUNd/e72/pcMqaG8tsVhIg4gJeBC4EUYLWIzFXVLSU2uxjo4pmGAq96Xk8Yqarp/oqxzlOF7/4JS/8ftOoLkz4kK6wFk99excqkIzx9RW8mDan+h9LyivN4bvVzxEfHc0WXK6r9840xvvFnFdMQYJeqJgKIyEzgMqBkgrgMeE9VFfhRRBqLSCtVPejHuOqHgmz4/B7YMgd6Xw3jX+LAcbjl1RXsTsvhXxP7cVm/NgEJbfqG6ezP2c9bF71FcJDVchpTU/mziqkNsK/EfIpnma/bKLBYRH4SkTIbyIvIZBFJEJGEtLS0Kgi7Dji0EaafB1vnwqjH4YrX2ZJWxIRXfuDAsTzevXVIwJLD7mO7eWfzO4yPH8/gloMDEoMxxjf+/PnmrQMfrcQ2I1T1gIg0B5aIyDZVXf6LjVWnA9MBBg0aVHr/9YsqJLwFC6dCRBO4cS50/BVLthzmDx+vIyo8mP/+bli1PgRXktPl5PEVjxMZHMkfB/0xIDEYY3znzwSRArQrMd8WOODrNqp64jVVRGbjrrL6RYIwHnlH4Ys/wObZEH8BTHgNV2RTXlyygxe/3kmfttFMv2EQLaPDAxbiB1s/YG3qWp465yliwmMCFocxxjf+rGJaDXQRkY4iEgpMAuaW2mYucKOnNdPZQKaqHhSRBiISBSAiDYDRwCY/xlq77fwKXhkGW+bCBY/CdbPIDGrM5Pd/4sWvd3LVwLZ88tthAU0OiZmJ/HvNvxnZbiSXdro0YHEYY3zntysIVS0WkXuARYADeEtVN4vInZ7104D5wFhgF5AL3OIp3gKYLSInYvxIVRf6K9ZaqyAbFv0Z1rwLzbrDpI+gzQASko/w+5nrOJyVz2PjenLT8A54/l8GRKGzkKnfTSUiJIJHhj0S0FiMMb7zaxMSVZ2POwmUXDatxHsF7vZSLhHo68/Yar2dX8GXf4Bj+2DE7+G8h3E6wnh16U5e+GonbRpHMOt3w+nXrnGgI+WFn15gS8YW/jXyXzSNqL5+nowxZ8baGNY2mfth4RR3C6XYLnDrIogbyo7D2Tw46yfW7zvGuL6teXJCLxqFhwQ6WpbtW8YHWz/g2u7XckHcBYEOxxhTCZYgagtnEax8DZb9HVzFcP5fYPh9FEkIry3dyb+/3kWDMAcvTurH+L6ta0Q1TlJmEg9/9zA9YnpYqyVjaiFLEDWdKmz7Er56FDJ2QZeLYOyz0KQDa/Ye5S+zN7HlYBaX9mnFY+PPomnDsEBHDEB2YTb3Lb2P4KBg/jXyX4Q6QgMdkjGmkixB1GQpCbD4r7D3f9C0K1zzMXS9iENZBTzz8Tpmr91P86gwpl0/kDG9ak5vqEWuIh5c/iAp2SlMHz2d1g2rtwtxY0zVsARREx1c7x7tbdsX0KAZXPoC9L+RPKfwxtJdvLJsN05V7hnZmd+dF0+DsJpzGlWVx/73GD/s/4HHhj1mT0sbU4vVnG8W4x6nYdkzsGMBhEfDeVNh2N3kB0Xywf/2MO3bRNJzCri4V0seHtuDdjGRgY74FKrKC2teYO7uudzV7y6u7HploEMyxpwBSxCB5nLBrq/gx1cg8RsIbwwj/wxDf0u+oyEzV+3llWW7Sc0uYHh8LK9eP4DBHWreU8iqyktrX+LtTW8zsdtE7uxzZ6BDMsacIUsQgVKQA+tnwMpp7pvPUa3cT0EPvp3UwlDe+3YPH67cw9HcIoZ0jOHf1/Tn7E6xgY7aK1XlxTUv8uamN7myy5U8PPThGtGKyhhzZixBVLdDm2DtB7D+I8jPhNYD4Mo30R7jWXfgOO/PTWTe+gMUu5RRPVpw+zkdGdIxpsZ+4Ra5injsf48xd/dcrup6FX89+68EiV/HoTLGVBNLENUh7xhsmgVr3oeD6yAoBHpcCkN/R3qTvsxee4BP/r2Cnak5RIY6uG5oe24e3oEOTWv2SGs5hTk8sOwBVhxcwV197+LOvnfW2ERmjKk8SxD+UpADOxfDls9hx0IozocWvWDMM2R1uZzFycXMX3qQ5TuWUuxSBsQ15ukrenNJn1ZE1YAnoCuSlJnEH7/9I4nHEvnb8L8xocuEQIdkjKliliCqUn4m7FjkTgq7vnInhQbNof8NpHe5mq+OtmTx1lS+m7eGIqfSpnEEt53TkasGtqVLi6hAR++zLxK/4G8r/ka4I5xXRr3C8NbDAx2SMcYPLEGcqbyjsH2BOynsXgrOQohqhbP/jWyLOZ95R9uxbMcRtn2XDqTTpnEENw/vwCV9WtO3bXStqpLJKsziudXPMWfXHAY0H8Cz5z5LiwYtAh2WMcZPLEGcjuMZ7ofYts6FxGXgKsbVqC0Hu97Ad8HDmJvehjU/ZpJf5CI4aC+DOjRh6sXdOa9bc7q2aFirkgK4Wykt3rOYp1c9zZH8I9zR+w7u6neXjSdtTB1nf+G+yj0CW+fB5s8g6TtQJ7kN2rG+xSQ+zRvEnNTmFKeCCPRs5WTS4DjO7hTLiM6xteKeQll2Hd3FC2teYHnKcnrE9ODlC16mZ2zPQIdljKkGliDKk58J275EN30KicsQVzHpoW1ZEjqBD7IGsDm/PaHHHPRtF83kc2MY3DGGge2b1Ihuts/UwZyDvLzuZebunkuDkAb8adCfuK7HdXbVYEw9Yn/tpanCvlXkr3yTkK1zcLgKOEQzPi++mHnOs9mrnRnYIYaxZ8fwSPsm9G3XmPAQR6CjrjLbj2zn3c3vsiBpAUESxI09b+T23rfTOLxxoEMzxlQzSxAnuJykr/oYvnuepsd3UqzhfOo8h8WhFxAVfzZDOsXyjw4xdG0RhSOodt1DqEh+cT5f7/2a2Ttns/LQSiKCI5jYfSI39byJVg1bBTo8Y0yA1PsEUVCQz/Nzfs+5Sd8xIn8/O11t+CT6PqT3VYzo2Z5rWkcTVMcSAkCRs4hVh1axZM8SFiUvIqcohzYN23D/gPu5qutVRIdFBzpEY0yA1fsEkZN/lHnHv+ejVg76BI3kjnP+yl0dBwY6rCqnqqRkp7Dq0Cp+PPgj3+//npyiHCKCIxgVN4rLO1/OoJaDrJsMY8xJoqqBjqHKDBo0SBMSEipdbvvO1SzKWsGH2z4ktziXXrG9uDT+Uka2G1lrB7vJL85n+9HtbMnYwqb0Taw+tJqDxw8CEBsey7ltz+WCuAs4u/XZhDlqxih0xpjqJyI/qeogr+ssQfwsqzCLebvnMWvHLHYd2wVA95juDGs1jD7N+tCnWR+aRzavqnCrRGZBJnuz9pKcleyeMpNJzEwkKTMJpzoBiAmPYWCLgQxuOZihLYfSMbpjrXsWwxjjHwFLECIyBngRcABvqOrTpdaLZ/1YIBe4WVXX+FLWmzNNECUlZyazbN8yvtn3DRvTN1LkKgKgSVgT4hrF0b5Re+Ki4mjVsBUx4TEnp0ahjQgPDj+tqpoiVxEFxQXkFeeRWZBJVmEWmQWZZBZmklmQybGCY6TmpnI49zCHjx8mNTeV3OLck+Ud4qBNwzZ0iO5A95ju9IztyVmxZ9EisoUlBGOMVwFJECLiAHYAFwIpwGrgGlXdUmKbscC9uBPEUOBFVR3qS1lvqjJBlFToLGTbkW1sSNvArmO72Ju9lz1Ze0jNTS2zTJgjjPDgcMId4YQEeX8uwqlOCpwF5BfnU+AsOPmLvywOcdAsshnNI5vTIrLFyal9o/a0j25Pu4btCHHU/mcwjDHVp7wE4c+b1EOAXaqa6AliJnAZUPJL/jLgPXVnqR9FpLGItAI6+FC22oQ6Qk9WMZWUV5xHWm4aR/KPkJGfwZH8I2QXZruvApx55Bfnk1+cT6GrEOGXv+CDJIgwR9jJZBLmCCPcEU5EcATRYdE0CmtEdFg00aHu9w1DGtpNZGNMtfFngmgD7Csxn4L7KqGibdr4WBYAEZkMTAaIi4s7s4grKSI4grhGccQ1qt7PNcaY6uDPn6PeKr1L12eVtY0vZd0LVaer6iBVHdSsWbNKhmiMMaYs/ryCSAHalZhvCxzwcZtQH8oaY4zxI39eQawGuohIRxEJBSYBc0ttMxe4UdzOBjJV9aCPZY0xxviR364gVLVYRO4BFuFuqvqWqm4WkTs966cB83G3YNqFu5nrLeWV9VesxhhjfskelDPGmHqsvGau1mbSGGOMV5YgjDHGeGUJwhhjjFd16h6EiKQBe06zeFMgvQrDCaS6cix15TjAjqUmqivHAWd2LO1V1etDZHUqQZwJEUko60ZNbVNXjqWuHAfYsdREdeU4wH/HYlVMxhhjvLIEYYwxxitLED+bHugAqlBdOZa6chxgx1IT1ZXjAD8di92DMMYY45VdQRhjjPHKEoQxxhiv6n2CEJExIrJdRHaJyJRAx1NZIpIsIhtFZJ2IJHiWxYjIEhHZ6XltEug4vRGRt0QkVUQ2lVhWZuwiMtVznraLyEWBidq7Mo7lMRHZ7zk36zxD7J5YVyOPRUTaicg3IrJVRDaLyO89y2vdeSnnWGrVeRGRcBFZJSLrPcfxuGe5/8+JqtbbCXdPsbuBTrjHoFgP9Ax0XJU8hmSgaallzwJTPO+nAM8EOs4yYj8XGABsqih2oKfn/IQBHT3nzRHoY6jgWB4D/uRl2xp7LEArYIDnfRTuseF71sbzUs6x1KrzgnsAtYae9yHASuDs6jgn9f0K4uS42apaCJwY+7q2uwx41/P+XeDywIVSNlVdDhwptbis2C8DZqpqgaom4e4ifkh1xOmLMo6lLDX2WFT1oKqu8bzPBrbiHgK41p2Xco6lLDXyWNQtxzMb4pmUajgn9T1BlDUmdm2iwGIR+ckzPjdAC3UPvITntXnAoqu8smKvrefqHhHZ4KmCOlEFUCuORUQ6AP1x/2Kt1eel1LFALTsvIuIQkXVAKrBEVavlnNT3BOHz2Nc12AhVHQBcDNwtIucGOiA/qY3n6lUgHugHHAT+6Vle449FRBoCnwL3q2pWeZt6WVbTj6XWnRdVdapqP9zDLw8RkV7lbF5lx1HfE4Qv42bXaKp6wPOaCszGfSl5WERaAXheUwMXYaWVFXutO1eqetjzh+0CXufny/wafSwiEoL7C/VDVf3Ms7hWnhdvx1JbzwuAqh4DlgFjqIZzUt8TRK0e+1pEGohI1In3wGhgE+5juMmz2U3A54GJ8LSUFftcYJKIhIlIR6ALsCoA8fnsxB+vxwTc5wZq8LGIiABvAltV9fkSq2rdeSnrWGrbeRGRZiLS2PM+AhgFbKM6zkmg79AHesI9JvYO3Hf6/xzoeCoZeyfcrRXWA5tPxA/EAl8DOz2vMYGOtYz4Z+C+xC/C/avntvJiB/7sOU/bgYsDHb8Px/I+sBHY4PmjbVXTjwU4B3d1xAZgnWcaWxvPSznHUqvOC9AHWOuJdxPwiGe538+JdbVhjDHGq/pexWSMMaYMliCMMcZ4ZQnCGGOMV5YgjDHGeGUJwhhjjFeWIIw5TSLSWETu8rxvLSKzAh2TMVXJmrkac5o8/ft8oarldXtgTK0VHOgAjKnFngbiPZ2o7QR6qGovEbkZd8+aDqAX7r5+QoEbgAJgrKoeEZF44GWgGZAL3KGq26r7IIwpi1UxGXP6pgC71d2J2oOl1vUCrsXdz8+TQK6q9gdWADd6tpkO3KuqA4E/Aa9UR9DG+MquIIzxj2/UPQZBtohkAvM8yzcCfTw9jA4H/uvuMghwD/BiTI1hCcIY/ygo8d5VYt6F++8uCDjmufowpkayKiZjTl827qEsK03d4xIkicjV4O55VET6VmVwxpwpSxDGnCZVzQB+EJFNwHOnsYvrgNtE5ERvvHVhuFtTh1gzV2OMMV7ZFYQxxhivLEEYY4zxyhKEMcYYryxBGGOM8coShDHGGK8sQRhjjPHKEoQxxhiv/j8jnn1tCYvHlQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(pt_A, label='A')\n",
    "plt.plot(pt_B, label='B')\n",
    "plt.plot(pt_AB, label='AB')\n",
    "plt.legend()\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('density of infected')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "7e668f53",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'density of infected')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAxMklEQVR4nO3deXgUVdb48e/JDiRhCRCWAAEE2QMYQUVRVBRUBlEZwN3R19dxG2feccTX+anjjDNuo477NvoiKi6AiMriMiogIgQEZJUlLAESSEIWyNrp8/ujOhJilk5Ip9Pp83mefqrr1u3KKYr0ya2qe6+oKsYYY4JXiL8DMMYY41+WCIwxJshZIjDGmCBnicAYY4KcJQJjjAlyYf4OoK7at2+viYmJ/g7DGGMCyurVqzNVtUNV2wIuESQmJpKSkuLvMIwxJqCIyO7qttmlIWOMCXKWCIwxJshZIjDGmCAXcPcIqlJaWkpaWhpFRUX+DqVRRUVFkZCQQHh4uL9DMcYEsGaRCNLS0oiJiSExMRER8Xc4jUJVycrKIi0tjZ49e/o7HGNMAGsWl4aKioqIi4sLmiQAICLExcUFXSvIGNPwmkUiAIIqCZQLxmM2xjS8ZnFpyBjTdLnK3BS73LjciqvMTZlbcbn156WrzH3cepnbTWnZ8etlbudyqFudpQJuVVSPLRXF7abCtvLtnm3ln9WqP1uufGT+igP0Vxytv6q6HFe36u112d9xu61QOTmxHaP7Vtkn7IRYImgAOTk5vPPOO9x666011tu1axfLly/nyiuv9Gq/P/zwA8OHD2fRokVceOGFDRGqMdUqcyu5haVkHy0ht7CE/CIXR4pdHC12kV/k4mhxGUdLyt87r2KXm2JXmbMsPfa+qNRT5nK++M2JKW/833J2b0sETVVOTg4vvPCCV4ngnXfe8ToRzJo1izPPPJNZs2ZZIjD1VlRaxoHcIg7kFpKeW8SB3CLSc4vIyCvicEEJ2UdLOFxQSk5BCbV9Z0eFhxAdGU50ZCgtI8JoERFKZFgI0ZFhRIaFEhkeQmRYiPM+LMSzHkpEWAjhoSGEhQihIUJ4qBAacmw9LEQIC/3leqhnPVQEEQj5xdK5RCr8cpuIONuRY/XKt5fXDwHh2D7KlX/xViyt7kpsxfLy+seXVawr1ZRXXaexWCJoANOnT2fHjh0MHTqUsWPH8thjj/GnP/2JhQsXIiL8+c9/ZsqUKUyfPp3NmzczdOhQrrvuOn7/+99Xu09VZfbs2Xz++eecddZZFBUVERUV1YhHZQJJmVtJzTzKtox8UrOOsivzKLsyC9iZeZTMI8W/qN+mZTjxMVG0axVBv06xtG0VTruWEbRrFUHbVhG0aRlBdGQYMVFhtIoMIzoijFaRoYSFNpvbiqaCZpcI/vLxRjbtz2vQfQ7oEssDEwZWu/2RRx5hw4YNrF27FoA5c+awdu1a1q1bR2ZmJqeeeiqjR4/mkUce4YknnuCTTz6p9Wd+++239OzZk969e3POOeewYMECLrvssoY6JBPAXGVuNh/IZ11aDpsO5LFpfx5b0vMoKnX/XKd9dCS92rfi3H4d6N6uJZ1bt6Bz6yg6t2lBp9goWkSE+vEITFPT7BJBU7Bs2TKmTZtGaGgo8fHxnH322axatYrY2Fiv9zFr1iymTp0KwNSpU5k5c6YlgiBVWuZmze7DrNiZzapd2azZc5iCkjIAYqPCGNilNVeN7MGAzrH0jY8hsX1LYqKsk6HxXrNLBDX95d5YtKpHCeqgrKyMOXPmMH/+fB5++OGfO4/l5+cTExPTQFGapuxQfjFfbz3I11sPsWTbIfKLXIhAv06xTD4lgeTEdgzr3oaubVrYY8TmhDW7ROAPMTEx5Ofn/7w+evRoXn75Za677jqys7NZsmQJjz/+OPv27TuuXnW++OILkpKSWLx48c9l1113HfPmzeOaa67xyTEY/8srKmXxhnTmr9vPt9szcSvEx0Zy8eDOnHNyR07vHUfrFvaXvml4lggaQFxcHKNGjWLQoEGMHz+exx57jO+++46kpCREhMcee4xOnToRFxdHWFgYSUlJXH/99UyZMoWbbrqJBQsWHLe/WbNmMWnSpOPKLr/8cl588UVLBM2MqvJ9ajYzV+zm800ZlLjcdG/XktvGnMS4QZ0Y0DnW/uI3PicnehmjsSUnJ2vliWk2b95M//79/RSRfwXzsQeyYlcZc1bvY8byXWzNyKd1i3AuHdqFicO6MqxbG/vyNw1ORFaranJV26xFYEwjKiwpY9bKPby8ZAcZecUM6BzLY5cPYUJSF3uSx/iNJQJjGkGZW5m9ei9PfPYTh/KLGdGzHf+cPJRRJwXXYImmabJEYIyPLd+eyUOfbGJLej6n9GjLc9OGMbJXnL/DMuZnlgiM8ZGcghIe+mQTc9fsI6FtC567chgXD+5sLQDT5Pg0EYjIOOBfQCjwmqo+Umn7OcBHQKqnaK6qPuTLmIxpDIs2HODP8zZyuKCE28ecxO3nnkRUuN0DME2TzxKBiIQCzwNjgTRglYjMV9VNlaouVdVLfBWHMY2pqLSMv3y8iVkr9zCwSywzfnMqA7u09ndYxtTIlyNIjQC2q+pOVS0B3gUm+vDn+U356KO1KR991BuJiYkMHjyYoUOHMnjwYD766KMTDdP4WGrmUSa9sJxZK/fw23N6M++2UZYETEDwZSLoCuytsJ7mKavsdBFZJyILRcT/40PUgy8SAcBXX33F2rVrmT17NnfeeeeJhGh87KutB5nw7DIO5BbyxvWncs+4foTbSJ0mQPjyf2pVd8Qq915bA/RQ1STgWWBelTsSuVlEUkQk5dChQw0bZQOoOAz13Xffjapy9913M2jQIAYPHsx77733c72lS5cydOhQnnrqKa/3n5eXR9u2bX0VvjlBb63YzU0zUugR15IFd57FmH4d/R2SMXXiy5vFaUC3CusJwP6KFVQ1r8L7BSLygoi0V9XMSvVeAV4Bp2dxjT914XRI//EEQ6+k02AY/0i1m30xDDXAmDFjUFV27tzJ+++/3xBHYhqQ2608umgLLy/Zybn9OvLstGG0irQH8Uwlx81LedyEldWXHVdeoUxCIbTh/4/58n/tKqCPiPQE9gFTgeOm5hKRTkCGqqqIjMBpoWT5MKZG0RDDUINzaah9+/bs2LGD8847j3POOYfo6GgfRW3qwu1W7pv3I7NW7uWa03rwwIQBwTVpi6sYivOhOM9ZlhaBq8gpdxU6y1LP0lVU4VUMZaXgLgW3C9xlnnXXL19lpc52t6tCfTegoG7Pq8L7n8v1WHmNdbXqcqjbl3N1dX1h1F0w9i8NvlufJQJVdYnI7cBinMdHX1fVjSJyi2f7S8AVwG9FxAUUAlP1RAc/quEv98bS0OM39e7dm/j4eDZt2sSIESMadN+m7txu5d65P/Jeyl5uH3MS/3NB38DtG+Aqhvx0OJoJBZmVlllQkO35ws/1LD2vspK6/ywJhbAoCA2HkLBjy5BQCCl/H+b8xRtS4RUWWWE91JnXUUIAz1JCqinj2Puq6v5cLpU+X34uj5s/suKBVF9W77rVzG1Z+fMJp1b1L3vCfNqOVdUFwIJKZS9VeP8c8JwvY2gMDT0MdWUHDx4kNTWVHj16NGTYph7cbmX63PW8n5LGneeexO/HNvEk4CqG7FTI2uYs8/ZBbtqx5dFq7rmFRkKr9tCyHUS2htgEiIyp9Ir1LKMhvIXzJV/xFV5p3QeXNEzDsDPTABp6GOpyY8aMITQ0lNLSUh555BHi4+Mb+chMZQ8v2OwkgfP68Iexff0dzjFuN2TvhPR1cGAdZGxyvvxz9ngueXhEREPrBIjtCp2GOO9jOkOrDp4v/jhnGRFd/WztptmxYagDXDAfe2N7belO/vbpZm4Ylcj9lwzwb0ug8DDsWQG7v4W01c4DEiWe1mZoBLQ/Gdr3gbiTji3b9YSoNvYFH6RsGGpjTtAn6/fzt083M35QJ/58sR+SQFkp7PkOfloMO7+BjA2AOl/6nZMgaaqz7JwEHfpBWETjxmcCmiUCY2qxZs9h/vDeOpJ7tOWpKUMJDWmkJFBSAFsXwJZPYfuXzg3b0Ajofhqccy8kjoKupzjX5405Ac0mEahq075p5wOBdlkvEGUeKebWt9bQqXUUr16b7PuB49xlsGsprH8fNs13Lve06ggDJkDf8dDrHOfmrDENqFkkgqioKLKysoiLC55JPlSVrKwsoqKi/B1Ks+Uqc3PHOz9wuKCEubeeQdtWPrzcUpANa96EVa9B7l6IiIGBE2HIVOgxCkKCqI+CaXTNIhEkJCSQlpZGUxx+wpeioqJISEjwdxjN1hOf/cR3O7N4YnKS7waPO7wLlj0N6951OmIlnuV0GDr5IrvkYxpNs0gE4eHh9OzZ099hmGbky80ZvPTNDq4c2Z0rTvFBss1OhaVPOAlAQpybvSNvgfiAHHfRBLhmkQiMaUhZR4q5Z856+nWK4YEJAxp254U5sORx+P5lJwEk3whn3gWxXRr25xhTB5YIjKlAVbnvww3kFbqYeeNIIsMa6Oaw2w2r34D//M3pAzD8GjjnfyG2c8Ps35gTYInAmArmrtnHoo3pTB/fj/6d6zZIYLUyt8FHt8PeFdDjTBj3D+g8pGH2bUwDsERgjMf+nEIenL+RUxPb8l9n9TrxHbrLYPmz8NXfnRu/l74ISdOsZ69pciwRGONx/0cbcbmVf05ugE5jRw7CnJsg9Rvodwlc/CTE2FhRpmmyRGAM8NnGdL7YnMG94/vRPa7lie1s5zcw97+gKBd+9RwMu9paAaZJs0Rggt7RYhcPzt/IyfEx/ObME3gMWdV5Gmjxvc4gb9fMg/gGfurIGB+wRGCC3jNfbmN/bhGzpw2r/4TzZS5YNB1WvepcCpr0sg0FYQKGJQIT1Lam5/PvZalMPbUbyYnt6reT4iPwwfWw/XM44w44/yEbEsIEFEsEJqj97dNNREeFcc+4fvXbQVEuvP1rSFsJlzwNyTc0aHzGNAZLBCZoffPTIZZuy+T+SwbUb0C5gmx46zJnUpgr3oCBlzZ4jMY0BksEJiiVuZW/f7qZHnEtufq0eswFXZANMyZA5k8w5W04eVzDB2lMI7ELmSYozVmdxtaMfO4Z14+IsDr+GhQfgbevcHoMX/meJQET8KxFYIJOQYmLJz7byvDubRg/qFPdPuwqhveuhv1rYcpb0Ptcn8RoTGOyFoEJOq8vS+VgfjH3Xdy/bhMZuctg7s2w8yuY+Bz0u8h3QRrTiCwRmKCSV1TKK0t2cn7/eE7pUcfHRb94EDbNgwv+BkOv9EV4xviFJQITVN5Ytou8Ihd3nd+nbh9cOwuWPwOn3uT0FTCmGbFEYIJGbmEpry3bydgB8QzqWoepJ/eugo/vdKaRHPeI7wI0xk8sEZig8fqyVPLr2hrI2w/vXQWxXeHXb0JouO8CNMZP7KkhExRyC0p5fVkqFw6M934i+jIXzL7ReVz02o+gZT2HoDCmibNEYILC69+mkl/s4nfn9fX+Q988AnuWw6RXoGN/3wVnjJ/59NKQiIwTka0isl1EptdQ71QRKRORK3wZjwlOBSUuZny3i/P7xzOgi5fTT+74CpY8AUOvhqQpvg3QGD/zWSIQkVDgeWA8MACYJiK/GJzdU+9RYLGvYjHB7f1Ve8kpKOW353g5/eSRg05/gfZ94aLHfBucMU2AL1sEI4DtqrpTVUuAd4GJVdS7A5gDHPRhLCZIlZa5eXVpKsk92nrXb0AVPv6dM6ro5P+DiFY+j9EYf/NlIugK7K2wnuYp+5mIdAUmAS/VtCMRuVlEUkQk5dChQw0eqGm+Fvx4gH05hdxydm/vPrD+Pdi6AM77fza7mAka1d4sFpFnAa1uu6reWcu+q+q7X3l/TwP3qGpZTV39VfUV4BWA5OTkamMypiJV5aVvdtKnYzTn9utY+wfy9sOCP0G30+C0W30foDFNRE0tghRgNRAFDAe2eV5DgTIv9p0GdKuwngDsr1QnGXhXRHYBVwAviMilXuzbmFot2ZbJ5gN53Dy6FyEhtYwppArz7wR3KVz6AoSENk6QxjQB1bYIVHUGgIhcD4xR1VLP+kvAZ17sexXQR0R6AvuAqcBxA7So6s8zhYvI/wGfqOq8Oh2BMdV4belO4mMjmTi0a+2V181yppoc/zjEeXkZyZhmwpt7BF2AmArr0Z6yGqmqC7gd52mgzcD7qrpRRG4RkVvqE6wx3tp+MJ+l2zK55rQetc83cDQLFt/nXBI69abGCdCYJsSbDmWPAD+IyFee9bOBB73ZuaouABZUKqvyxrCqXu/NPo3xxozlu4kIC2HaiO61V/7ifijOg0uesknnTVCqNRGo6hsishAY6Smarqrpvg3LmPrLLSxlzpo0JiZ1IS46subKu5fDD2/BqLvsKSETtGr980ecx3nOB5JU9SMgQkRG+DwyY+rpg5S9FJSUcf2oxJorukrg47ugdXc4+0+NEZoxTZI37eAXgNOBaZ71fJwew8Y0OWVuZcZ3uxjRs13tg8uteAEyt8LFT1jHMRPUvEkEI1X1NqAIQFUPAxE+jcqYevpycwZ7swu54YzEmivmZzhjCfUdB30vbJTYjGmqvEkEpZ7xgBRARDoAbp9GZUw9zVyxm86toxg7IL7miv/5K7iK4IKHGycwY5owbxLBM8CHQEcReRhYBvzDp1EZUw+7s46ydFsm00Z0Jyy0hv/aB9Y5N4hH/je0P6nxAjSmifLmqaG3RWQ1cB7OsBGXqupmn0dmTB29s3IPoSHClFO7VV9JFRZOdyaZGX134wVnTBNWayIQkZmqeg2wpYoyY5qEEpeb2SlpnNevI/GxUdVX3PSRM9nMJU9BizaNFp8xTZk3l4YGVlzx3C84xTfhGFM/izemk3W0hCtH1tCBrKwUvnwIOvSH4dc1XnDGNHHVJgIRuVdE8oEhIpLneeXjzBvwUaNFaIwX3v5+NwltWzC6T4fqK/0wE7J3wPkP2KByxlRQbSJQ1X+oagzwuKrGel4xqhqnqvc2YozG1GjHoSOs2JnNtBHdqx9ltKQAvn7UGU+o77jGDdCYJs6bS0MrReTnnjki0saGijZNyazv9xAWIkxOTqi+0vcvwZF0OP9BqGHuC2OCkTeJ4AFVzS1fUdUc4AGfRWRMHRSVljF7TRoXDIynY0w1N4kLsmHZ005LoMfpjRqfMYHAm0RQVR1vRi01xucWbjhATkEpV43sUX2lb592Rhc97/5Gi8uYQOJNIkgRkSdFpLeI9BKRp3BmLjPG72at3EtiXEtO7xVXdYXcffD9yzBkCsQPrLqOMUHOm0RwB1ACvAe8DxQCt/kyKGO8sTvrKCtTs5mc3K36m8RLHgd3GYz538YNzpgA4k3P4qPAdBGJVtUjjRCTMV6ZszoNEZg0rJqpKHP2OENJnHIdtK3h0pExQc6b+QjOEJFNwCbPepKIvODzyIypgdutzFmzjzNPak+XNi2qrrT0SecJoTN/37jBGRNgvLk09BRwIZAFoKrrgNG+DMqY2izfkcW+nEImJ1czrlB5a2D4tdC6hsdKjTFeJQJUdW+lojIfxGKM12av3ktMVBgXVDfctLUGjPGaN4lgr4icAaiIRIjIHwEbfdT4TV5RKQs3pPOrpC5EhVcxVIS1BoypE28SwS04Twl1BdKAodhTQ8aPPll3gGKXu/rLQtYaMKZOqn1qSEQeVdV7gDGqelUjxmRMjWav3kufjtEkJVQxJ3HFJ4WsNWCMV2pqEVwkIuGADTBnmoztB4+wZk8Ok5MTkKrGDFr2lLUGjKmjmvoRLAIygVYikoczO5mWL1U1thHiM+Y4s1enERoiXFpV34H8dKc1MPQqaw0YUwc1DUN9t6q2Bj6tMAT1z8tGjNEYAFxlbuauSeOcvh2qHmDuu+fA7YJRv2v84IwJYLXeLFbViY0RiDG1Wbo9k4P5xVUPN12QDateh0FXQLuejR+cMQHMm57Fl4nINhHJLZ+lzHOpyJhGNTsljXatIji3XxV9B1a+AqVH7d6AMfXgzXDSjwETVNX6Dhi/ySko4fNNGVx1Wnciwir9/VKcDytehJMvhvgB/gnQmADmTT+CjPomAREZJyJbRWS7iEyvYvtEEVkvImtFJEVEzqzPzzHN30dr91NS5uaKU6q4LLT6/6AoB876Q2OHZUyz4E2LIEVE3gPmAcXlhao6t6YPiUgo8DwwFqcj2ioRma+qmypU+xKYr6oqIkNwhrnuV7dDMMFg9uo0BnSOZWCXSn0HSotg+XPQ82xISPZPcMYEOG8SQSxQAFxQoUyBGhMBMALYrqo7AUTkXWAinlFMASoNa93Ks19jjrMlPY8f9+XywIQqLvuse8eZi/iyVxo/MGOaCW/mI7ihnvvuClQcrC4NGFm5kohMAv4BdAQurmpHInIzcDNA9+7d6xmOCVQfpKQRHipMHFqp70CZy5mLuGsy9LQBcY2pr5qGmPiTqj4mIs9SxV/qqnpnLfuuasqoqvbzIfChiIwG/gqcX0WdV4BXAJKTk63VEERKy9zM+2Ef5/WLp12riOM3bpwLObth3D+c3sTGmHqpqUVQfoM4pZ77TgMqjgqWAOyvrrKqLvHMi9xeVTPr+TNNM/PVloNkHS35Zd8Bt9sZTqJDf+g73j/BGdNMVJsIVPVjz3JGPfe9CugjIj2BfcBU4MqKFUTkJGCH52bxcCACzwQ4xgB8sDqNDjGRnN23w/EbfloEBzfBZa9CiFfTahhjquHNzeJ6UVWXiNwOLAZCgddVdaOI3OLZ/hJwOXCtiJQChcAUVbVLPwaAzCPFfLXlIDee2ZOw0Apf9qqw9Alo0wMGXua/AI1pJnyWCABUdQGwoFLZSxXePwo86ssYTOCa98M+XG79Zd+B1CWwbzVc8hSE+vS/sDFBodo2tYg86llObrxwjHGoKh+kpJHUrQ194mOO37j0nxDdCZKurPrDxpg6sfkITJO0YV8eWzPymVy5NZC2GlK/gdNvg/AqRiA1xtSZzUdgmqQPVu8lIiyECUldjt+w7EmIagPJ9e3eYoypzOYjME1OUWkZH63dz4UDO9G6RfixDQe3wJZPYMTNEBlT/Q6MMXXiTc/iiSISD5zqKfpeVQ/5NiwTzD7flEFuYSlTKk9O/+3TEN4SRt7il7iMaa68mY9gMrASmAz8GlgpIlf4OjATvN5P2UvXNi04o3fcscLDu2H9+3DK9dAqrtrPGmPqzptn7/4MnKqqBwFEpAPwBTDbl4GZ4LQvp5Bl2zO549w+hIRUGDZi+bMgIXD67f4LzphmypsumSHlScAjy8vPGVNnc1enocrxTwsdOQg/zISkKdC6iknrjTEnxJsWwSIRWQzM8qxPoVInMWMagtutfLA6jTN6x9GtXctjG1a8CK5iGHWX32Izpjnz5mbx3SJyGXAmzqOjr3hGDDWmQX2fms2e7AL+MLbvscKiXFj1GgyYCO37+C84Y5oxr/rne2Yjq20iGmNOyAcpe4mJCmPcoE7HCle9BsV5Ng2lMT5k1/pNk5BXVMqCDQf4VVIXosJDncLSQueyUO/zoHOSfwM0phmzRGCahE/WHaCo1M2vK/Yd+OEtOHrIWgPG+Jg3/QguERFLGMan3k/ZS9/4aIYkeCanLyuFb5+BhBHQY5R/gzOmmfPmC34qsE1EHhOR/r4OyASfnzLyWbs3h18nd0PKp5zcMAdy98BZ/2PTUBrjY7UmAlW9GhgG7ADeEJHvRORmEbHBXkyD+CBlL2EhwqXDPH0E3GWw9EnoOBD6Xujf4IwJAl5d8lHVPGAO8C7QGZgErBGRO3wYmwkCJS43c9fs49x+HWkfHekUbvoIMrfCaGsNGNMYvLlH8CsR+RD4DxAOjFDV8UAS8Ecfx2eaucUb08k6WsK0kd2dArcbvnkM2p8MAy71a2zGBAtv+hFcATylqksqFqpqgYj8xjdhmWDxzvd7SGjbgtF9PJPTb54PhzbD5f+GkFD/BmdMkPDm0tCBykmgfBpLVf3SJ1GZoLDz0BG+25nFtBHdCQ0RpzWw5HGI6wMDJ/k7PGOChjeJYGwVZeMbOhATfGat3ENYiDA52TPA3NZPIWMDjL7bWgPGNKJqLw2JyG+BW4HeIrK+wqYY4FtfB2aat6LSMj5YncbYAfF0jIkCVfjmUWjXCwZd7u/wjAkqNd0jeAdYCPwDmF6hPF9Vs30alWn2Fm1IJ6eglCvLbxJvXQjpP8KlL0KoV0NgGWMaSE2/caqqu0TktsobRKSdJQNzIt75fg894loyqnf7Y62Btokw+Nf+Ds2YoFNbi+ASYDWgOENQl1Oglw/jMs3Ytox8Vu7KZvr4fs4sZD8thgNr4VfPWWvAGD+o9rdOVS/xLHs2XjgmGLyzcg/hocIVpyQ4rYGvH4E23SFpqr9DMyYoedOhbJSItPK8v1pEnhSR7r4PzTRHBSUu5qxO48KBnZyexFsXwP41zpNCoeH+Ds+YoOTN46MvAgUikgT8CdgNzPRpVKbZ+vCHfeQVubjujERnTKH//A3iToKkK/0dmjFBy5tE4FJVBSYC/1LVf+E8QlorERknIltFZLuITK9i+1Uist7zWu5JNqaZUlVmLN/FwC6xJPdo64wwenATjLnP7g0Y40feJIJ8EbkXuBr4VERCccYcqpGn3vM4nc8GANNEZEClaqnA2ao6BPgr8EpdgjeB5bsdWfyUcYTrz0hE3C746mHoNNjGFDLGz7xJBFOAYuBGVU0HugKPe/G5EcB2Vd2pqiU4I5dOrFhBVZer6mHP6gogwevITcB5Y/ku2rWKYEJSF/hhJhzeBefeDyE275Ex/uTNfATpqvqkqi71rO9R1Te92HdXYG+F9TRPWXVuxOnA9gue+Q9SRCTl0KFDXvxo09TszS7gy80ZTBvRjShKnBFGu42EPlWNYGKMaUzePDV0mYhsE5FcEckTkXwRyfNi31UNJK/V/IwxOIngnqq2q+orqpqsqskdOnTw4kebpmbmit2ICFef1gNWvQb5B+C8+22+AWOaAG/u0D0GTFDVzXXcdxpQYSZyEoD9lSuJyBDgNWC8qmbV8WeYAHC02MW7K/cwbmAnOkcUO7OP9RoDiWf6OzRjDN7dI8ioRxIAWAX0EZGeIhKBM/fx/IoVPP0R5gLXqOpP9fgZJgC8t2oveUUubjqrJyx9AgoPw9iH/B2WMcbDmxZBioi8B8zDuWkMgKrOrelDquoSkduBxUAo8LqqbhSRWzzbXwLuB+KAFzyTlrtUNbk+B2KaptIyN/9elsqIxHYMi8mF71+GoVdC5yH+Ds0Y4+FNIogFCoALKpQpzl/yNVLVBcCCSmUvVXh/E3CTV5GagLTgxwPsyynkL78aCF/cAyFhcO6f/R2WMaaCWhOBqt7QGIGY5kdVefmbnfTu0IpzW+2CjR/C2fdAbBd/h2aMqcCbp4b6isiXIrLBsz5EROxPOlOrb7dnselAHjef1ZOQz++D6Hg4405/h2WMqcSbm8WvAvcCpQCquh7nxq8xNXp5yQ46xERyWdQqSFvlXBKKjPZ3WMaYSrxJBC1VdWWlMpcvgjHNx/q0HJZuy+Sm0zoR/p+HIH4QDL3K32EZY6rgzc3iTBHpjaczmIhcARzwaVQm4D3z5XZatwjnevc8yNkN1863CemNaaK8SQS34QwG109E9uEMFHe1T6MyAW3Dvly+2JzBg6NaErniGWcy+l5n+zssY0w1vHlqaCdwvmdymhBVzfd9WCaQPfef7cREhnLV4eecyWYueNjfIRljalBtIhCRP1RTDoCqPumjmEwA25Kex6KN6fxr6D7Ct3zhJIHYzv4OyxhTg5paBOWTz5wMnMqx4SEmAEt8GZQJXM/9ZztxES4m7H8GOg6Akf/t75CMMbWoafL6vwCIyGfA8PJLQiLyIPBBo0RnAsqW9Dw+/fEAb/f6ipB9e+GyBTYPsTEBwJvHR7sDJRXWS4BEn0RjAtrji7YyJDKd09PfhiFTIHGUv0MyxnjBm6eGZgIrReRDnEdIJwEzfBqVCTirdmXz1ZZ0lnd4A3FF2w1iYwKIN08NPSwiC4GzPEU3qOoPvg3LBBJV5dGFW7i95Zd0yv8RLnsVom0CIWMChTctAlR1DbDGx7GYAPXl5oNk7NnCnS1nQZ8LYPBkf4dkjKkDmzXcnJAyt/L4oi083fINQkPD4ZKnbPpJYwKMVy0CY6rz7qo9jMiayynh62Hck9A6wd8hGWPqyBKBqbfcglJmL/oP74bPQk86H0n+jb9DMsbUg10aMvX2r8838ZeyZwiNbIFMfN4uCRkToKxFYOpla3o+bVOeZkjoTpj4JsR08ndIxph6shaBqTNV5b0573Jr6IcUD5wCAyb6OyRjzAmwRGDqbPGqjfzXwb9zpGU3Iic84e9wjDEnyC4NmTo5fKSI2AW3ESf5hF49D6Ji/R2SMeYEWYvA1MnyN/8fZ7CWrLMeIrRLkr/DMcY0AEsExmvrl85nXMarbI67gM7n/tbf4RhjGoglAuOVwoM76P7lrewN6UrPG161R0WNaUYsEZjalRwl5/XJiJaRM3EGUdFt/B2RMaYBWSIwNVMl/c0b6Vi4k0X9Hmbo0GR/R2SMaWCWCEyNji7+K53SFjKj5XVcOvk6f4djjPEBSwSmWrrmTVqt+Cez3ecw+vq/EhkW6u+QjDE+4NNEICLjRGSriGwXkelVbO8nIt+JSLGI/NGXsZg62v4F+vFdLCkbTOGFT3BSvPUXMKa58lkiEJFQ4HlgPDAAmCYiAypVywbuBKx7alOyfy1l717LlrIE5vX5O1efcZK/IzLG+JAvWwQjgO2qulNVS4B3geMGpVHVg6q6Cij1YRymLg5uxv3mJDJcLflLzAM8NOUMxB4VNaZZ82Ui6ArsrbCe5imrMxG5WURSRCTl0KFDDRKcqUL2TvTNS8kpgd+4/8zD111AdKSNQmJMc+fLRFDVn5Fanx2p6iuqmqyqyR062KToPpGzF53xK44WFjKlcDp3XHEhJ3WM8XdUxphG4MtEkAZ0q7CeAOz34c8z9ZWdCm9cRMmRw0wtuJsJ55/LxUM6+zsqY0wj8WUiWAX0EZGeIhIBTAXm+/DnmfrI3A7/dzHFBblcXjCdfsNHc8e5dnPYmGDiswvAquoSkduBxUAo8LqqbhSRWzzbXxKRTkAKEAu4ReQuYICq5vkqLlNBxiaYeSklpaVcXnAvsb2G8fdJg+3msDFBxqd3AlV1AbCgUtlLFd6n41wyMo1t1zKYdSUlIZFMKvhfyjqczItXn0JEmPUxNCbY2G99MNr4IcycRGGLDlx09H6K2/blrZtG0rpFuL8jM8b4gT0bGExUYfmz8Pn9HO04nAsO3kZYTDvev2kk7aMj/R2dMcZPLBEEi9Ii+Ph3sP5dsrqP54LdV9GiZSvevmkk8bFR/o7OGONHlgiCQd5+eO9q2LeaHQN/x/h1p9GjXStm3jiSTq0tCRgT7CwRNHc/fQbzboHSIr495Wmu/S6egV1imXHDCNq2ivB3dMaYJsBuFjdXrhL47M/wzmQ0phMv9Xudq77tyOm94nj7ppGWBIwxP7MWQXOUnQpzboJ9KZQOu4E7s69g4apcrjmtB/dPGEB4qOV/Y8wxlgiaE7cbVr4MXz4EIeGkjX2Ja7/rzO7sPB6aOJBrT0/0d4TGmCbIEkFzkbkNProd9q5A+1zAxwl3c/fCLGJbuHjrxpGc3jvO3xEaY5ooSwSBrrQIvnsWljwBYVEcueh57vmpP58uTOfMk9rz1JShdIixPgLGmOpZIghUqrDlE1h8H+TshgET+arXH7l7UQa5hRncfeHJ3HJ2b0JDbNwgY0zNLBEEovQNsPh/IfUb6NCfrMtnc//6OD6dvZeBXWJ566aR9OtkcwwbY7xjiSCQZG6Dr/8BG+ZCVCxl4x7l9aIxPPVBKmXuDP5nbF9uOae3PRVkjKkTSwSBIHuncw9g3SwIa4Ge+QcWtb6CR7/OYFfWds7v35EHJgykW7uW/o7UGBOALBE0ZWkpsPwZ2PwxhISjI3/L8s7X8PDXmWw6kEq/TjG8cf2pjOnX0d+RGmMCmCWCpqasFLYuhBUvwJ7vIKo17jN+xzdtJvHsqiOs+Xon3dq14OkpQ/lVUhdC7GawMeYEWSJoKrJTYc2b8MNbcPQgtO5O6di/M0fH8PKKg6RmptGtXQv+eukgpiR3swlkjDENxhKBPxUehk3zYcNsSF0CEgJ9LmRPz1/zekZv5n6eTl5RKkMSWvPclcMYN7ATYXYj2BjTwCwRNLaCbNj2uTNL2PYvwF0K7Xpx5PS7WRwxlhkbS1m/PpeI0P2MG9SJK0d2Z2TPdjaPsDHGZywR+JoqHNoCPy2CnxbD3u9B3RDThYJhN/JV+Ghm7m7L918fRjWTfp1ieGDCAC4d2tVGCDXGNApLBA1NFbJ2wO5vYfdyZ5m7FwB3pyHsH3wbX5YN47197dn07REATupYyh3n9uGSIZ3pGx/jz+iNMUHIEsGJOpoJB9Z5Xmthzwo4kgGAu2UHstqfwtq4q5h3dBBf7gujaJebsBAhOTGCe8b147z+He3L3xjjV5YIvKHqfLlnbT/2ytwO6T9CXtrP1QpbJbC75XBWt+jPoiO9WJrdFrKFEIEBXWKZemo7Tu8dxxm944iJCvfjARljzDFBkwiO5OeQd2g/rSKEVuEhhIUAbheUHIXiPCjOd15Fuc6Xfn465KejeQfQ/P2ElBb8vC9XSARZEV3ZGdKXtRFjWV6QwDpXN/KKoiELEuNaMiAhlv9JjmVo9zYM696W6Mig+ac2xgSYoPl2+mnZXIZ//3uv6hYTzkFtR7q2IUPbk6F92K0dSdXOpGpn9mscrcoi6B7Xkh5xLRnYrhXj27Wkb3w0/TrH2pe+MSagBM03Vqf+o1hZ/A8KSt0UupTCUqXQpRTQggJpSYG04Ki0xB0RDZGtiQoPIyo8hBbhobRtFcEZ0RFMiI6kfXQk7VpFEBsVZo90GmOahaBJBF0ST6ZL4sn+DsMYY5oc66ZqjDFBzqeJQETGichWEdkuItOr2C4i8oxn+3oRGe7LeIwxxvySzxKBiIQCzwPjgQHANBEZUKnaeKCP53Uz8KKv4jHGGFM1X7YIRgDbVXWnqpYA7wITK9WZCLypjhVAGxHp7MOYjDHGVOLLRNAV2FthPc1TVtc6xhhjfMiXiaCqZyu1HnUQkZtFJEVEUg4dOtQgwRljjHH4MhGkAd0qrCcA++tRB1V9RVWTVTW5Q4cODR6oMcYEM18mglVAHxHpKSIRwFRgfqU684FrPU8PnQbkquoBH8ZkjDGmEp91KFNVl4jcDiwGQoHXVXWjiNzi2f4SsAC4CNgOFAA31Lbf1atXZ4rI7nqG1R7IrOdnmxo7lqapuRxLczkOsGMp16O6DaL6i0vyzZaIpKhqsr/jaAh2LE1TczmW5nIcYMfiDetZbIwxQc4SgTHGBLlgSwSv+DuABmTH0jQ1l2NpLscBdiy1Cqp7BMYYY34p2FoExhhjKrFEYIwxQS5oEkFtQ2I3dSKyS0R+FJG1IpLiKWsnIp+LyDbPsq2/46xMRF4XkYMisqFCWbVxi8i9nnO0VUQu9E/UVavmWB4UkX2e87JWRC6qsK0pH0s3EflKRDaLyEYR+Z2nPKDOTQ3HEXDnRUSiRGSliKzzHMtfPOW+Pyeq2uxfOB3adgC9gAhgHTDA33HV8Rh2Ae0rlT0GTPe8nw486u84q4h7NDAc2FBb3DjDla8DIoGennMW6u9jqOVYHgT+WEXdpn4snYHhnvcxwE+emAPq3NRwHAF3XnDGXov2vA8HvgdOa4xzEiwtAm+GxA5EE4EZnvczgEv9F0rVVHUJkF2puLq4JwLvqmqxqqbi9Dgf0RhxeqOaY6lOUz+WA6q6xvM+H9iMM/JvQJ2bGo6jOk3yOADUccSzGu55KY1wToIlETSH4a4V+ExEVovIzZ6yePWMzeRZdvRbdHVTXdyBep5u98yw93qFZnvAHIuIJALDcP4CDdhzU+k4IADPi4iEisha4CDwuao2yjkJlkTg1XDXTdwoVR2OM6vbbSIy2t8B+UAgnqcXgd7AUOAA8E9PeUAci4hEA3OAu1Q1r6aqVZQ1meOp4jgC8ryoapmqDsUZiXmEiAyqoXqDHUuwJAKvhrtuylR1v2d5EPgQpwmYUT6jm2d50H8R1kl1cQfceVLVDM8vrxt4lWNN8yZ/LCISjvPl+baqzvUUB9y5qeo4Avm8AKhqDvA1MI5GOCfBkgi8GRK7yRKRViISU/4euADYgHMM13mqXQd85J8I66y6uOcDU0UkUkR64sxlvdIP8XlNjp9adRLOeYEmfiwiIsC/gc2q+mSFTQF1bqo7jkA8LyLSQUTaeN63AM4HttAY58Tfd8ob8Y78RThPFOwA7vN3PHWMvRfO0wHrgI3l8QNxwJfANs+ynb9jrSL2WThN81Kcv2BurClu4D7POdoKjPd3/F4cy0zgR2C95xezc4Acy5k4lxHWA2s9r4sC7dzUcBwBd16AIcAPnpg3APd7yn1+TmyICWOMCXLBcmnIGGNMNSwRGGNMkLNEYIwxQc4SgTHGBDlLBMYYE+QsERhTCxFpIyK3et53EZHZ/o7JmIZkj48aUwvPGDafqGpN3f2NCVhh/g7AmADwCNDbMxjYNqC/qg4SketxRoIMBQbhjGcTAVwDFAMXqWq2iPQGngc6AAXAf6nqlsY+CGOqY5eGjKnddGCHOoOB3V1p2yDgSpyxbB4GClR1GPAdcK2nzivAHap6CvBH4IXGCNoYb1mLwJgT85U64+Dni0gu8LGn/EdgiGdUzDOAD5xhcQBnIhFjmgxLBMacmOIK790V1t04v18hQI6nNWFMk2SXhoypXT7ONIh1ps7Y+KkiMhmc0TJFJKkhgzPmRFkiMKYWqpoFfOuZtP7xeuziKuBGESkfPbY5TJNqmhF7fNQYY4KctQiMMSbIWSIwxpggZ4nAGGOCnCUCY4wJcpYIjDEmyFkiMMaYIGeJwBhjgtz/BwzyzjIt0WnjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "tot_A=[]\n",
    "tot_B=[]\n",
    "for i in range(len(pt_A)):\n",
    "    tot_A.append(pt_A[i]+pt_AB[i])\n",
    "    tot_B.append(pt_B[i]+pt_AB[i])\n",
    "\n",
    "plt.plot(tot_A, label='tot. A')\n",
    "plt.plot(tot_B, label='tot. B')\n",
    "plt.legend()\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('density of infected')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58a46578",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf33db58",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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": 5
}