06-SFE-Estimation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1362f56e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import copy\n",
    "import os\n",
    "import random\n",
    "import math as math\n",
    "import matplotlib.pylab as plt #Typically matplotlib.pylab is imported as plt (shortened name)\n",
    "import pdb\n",
    "from scipy.optimize import minimize\n",
    "import shutil\n",
    "import csv\n",
    "from MyLatexFunctions_SocialLearning import * \n",
    "from itertools import permutations \n",
    "from itertools import combinations\n",
    "import csv\n",
    "\n",
    "\n",
    "import networkx as nx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a841e70d",
   "metadata": {},
   "outputs": [],
   "source": [
    "datafolder = \"..\\\\RawData\\\\Session_Formal\"\n",
    "datafolder2 = \"..\\\\ImportedData\\\\\"\n",
    "outputfolder = \"..\\\\Output_Figures\\\\\" \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "662ad055",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAD3CAYAAAC+eIeLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+VklEQVR4nO3deXzM194H8M+syWQTSwQRJEhMiKwVVbHcRpespKidcrlapUoXfWi1dV20F/dWF0p7q1cUibQldBEk0XpoJbLJEGsWsolEIsns5/lDM48xoZjlN5P5vl+vvvoSyZnvSHycOXPO9/AYYwyEEEIsgs91AYQQYk8odAkhxIIodAkhxIIodAkhxIIodAkhxIIodAkhxIIodAkhxIIodO3Q/PnzsWrVKq7LMJCRkYGePXtyXYZJ9OnTB+np6VyXQawQhW4706dPH0gkEri6usLd3R3Dhg3D5s2bodVqdZ+zefNmvP322yZ/7FmzZmHFihUmH7cVj8dDYGCg3nNZsWIFZs2a9UBfP2rUKGzbts1M1Rnv448/Rnh4OBwcHB74ORHbQ6HbDu3fvx+NjY0oKSnBsmXLsG7dOsyZM+eBvlatVpu5OuNcu3YNu3bt4rqMezLmz69Hjx5YsWIFZs+ebcKKiLWh0G3HOnTogPj4eOzevRvbt29HYWEhAP0ZaetL+nXr1qFbt2544YUXoNVqsXbtWvTt2xedO3fGxIkTcePGDd24v/zyC4YNGwZ3d3d4e3vjq6++wueff46kpCR88MEHcHFxQVxcHIDbIfncc8/Bw8MDPj4++Oijj3TjtLS0YNasWejYsSMCAgLw+++//+lzeuONN7By5cp7htuJEyd0tQUFBSEjIwMAsHz5chw7dgwvv/wyXFxc8PLLL2PlypVYuHAhAEClUsHZ2RlvvPGGrjZHR0fU1dUBAPbt24eBAwfC3d0do0aNgkwm0z1mnz59sG7dOgwePBjOzs4GtZ09exY+Pj5/+o9FYmIixo4di86dO//pnwOxYYy0K71792aHDh0y+Li3tzf79NNPGWOMzZw5ky1fvpwxxtjRo0eZQCBgb7zxBpPL5ay5uZlt3LiRRUREsLKyMiaXy9m8efPYpEmTGGOMlZSUMBcXF7Zz506mVCrZ9evX2enTpw3GZYwxjUbDQkND2XvvvccUCgW7ePEi8/HxYT/++CNjjLE333yTDR8+nNXW1rLS0lI2cOBA5uXldc/nBoAVFxez0NBQtnXrVsYYY8uXL2czZ85kjDFWXl7OOnXqxA4cOMA0Gg37+eefWadOnVh1dTVjjLGRI0fqvo4xxg4fPswGDRrEGGPs119/Zb6+vmzIkCG63xs8eDBjjLFz584xJycn9vPPPzOlUsnWrVvH+vbtyxQKhe7PPCgoiJWWlrLm5ma970N2djbz9vZm+/fvf6Dv393PibQ/NNO1Ez169NCbrd6Jz+fjvffeg4ODAyQSCbZs2YLVq1ejZ8+ecHBwwLvvvouUlBSo1WokJSUhKioKkydPhkgkQufOnREcHNzmuL///jtqamrwzjvvQCwWw9fXF3PnztXN+Pbs2YPly5ejU6dO8Pb2xqJFi/70efB4PKxatQrvv/8+FAqF3u/t2LED0dHRiI6OBp/Px5gxYxAeHo6DBw+2Odbjjz+O8+fPo7a2FllZWZgzZw6uXr2KW7duITMzEyNHjgQA7N69GzExMRgzZgxEIhFee+01tLS04Pjx47qxFi1aBG9vb0gkEt3Hjh07hvj4eGzfvh2xsbF/+tyIfRByXQCxjKtXr6JTp05t/p6HhwccHR11vy4pKcG4cePA5///v8kCgQBVVVUoKytD3759H+gxS0pKcO3aNbi7u+s+ptFoEBkZCeD20oO3t7fu93r37v1A40ZHR6NXr174/PPPDR4vOTkZ+/fv131MpVJh9OjRbY4jkUgQHh6OzMxMZGVlYfny5cjNzcWvv/6KzMxM3dLDtWvX9Grj8/nw9vbG1atXdR+783m02rx5M0aOHHnPxyf2iWa6duD333/H1atXMXz48DZ/n8fj6f3a29sbP/zwA+rr63X/yeVyeHl5wdvbGxcvXnzgcXx8fPTGaWxs1M08u3fvjrKyMt3nl5aWPvBz+vvf/47Vq1ejublZ7/GmT5+u93hNTU1YtmxZm/UBwMiRI3HkyBGcPn0ajz32GEaOHImffvoJv/32G0aMGAHg9quEkpIS3dcwxlBWVgYvL697PnfgduiWlpbi1VdffeDnRdo/Ct12rKGhAWlpaZg0aRKmTZuGwMDAB/q6+fPnY/ny5bqgqampwffffw8AmDp1KtLT07Fnzx6o1WrU1tYiNzcXAODp6YlLly7pxhkyZAjc3Nywbt06tLS0QKPRoLCwUPeG2cSJE7FmzRrU1dWhvLwcmzZteuDnNmrUKAQGBmL79u26j02bNg379+/HTz/9BI1GA7lcjoyMDJSXl7dZH3A7dL/++msEBARALBbrtpX5+PjAw8NDV+eBAwdw+PBhqFQqrF+/Hg4ODhg2bNh9a3R1dcWPP/6IrKwsXfDfj1qthlwuh0aj0dVv7btJyCPgelGZmFbv3r2Zo6Mjc3FxYW5ubmzo0KHs448/Zmq1Wvc5d7+RdvebVxqNhq1fv575+fkxFxcX5uvry9566y3d72dlZbEhQ4YwV1dX1rNnT/bVV18xxhgrLi5mQUFBrEOHDiwhIYExxtjVq1fZpEmTmKenJ3N3d2cRERG6N/qamprY9OnTWYcOHZhUKmUffPDBn76Rdv78ed2vT5w4wQDovel04sQJNmLECNaxY0fWpUsXFh0dzUpKShhjjB0/fpz179+fubu7s4ULFzLGGGtsbGRCoZC9++67jDHGtFot8/DwYPPnz9d77NTUVCaVSpmbmxsbMWIEKyws1Pszv/vNyzs/VltbywYPHsxWrFhxz+fGGGMrV65kAPT+W7ly5X2/htgeHmN0cwQhhFgKLS8QQogFUegSYiGlpaVwcXFp87+HeROR2DZaXiCEEAuifbqE3EddkxLJ2WVIl1VBVtkIpVoLsZAPxhgYA3i829vFWj8u7eaKKKknJoR5o6OzmOvyiRWimS4hbcgtq8fWY5eQLqsCjwfIVdo//6I/OIr4YAyIknpibqQvgr3dzVcosTkUuoTcob5ZibdSC5BRXAOFWgOtEX87+DzAQSjAKD8PrEkMhLsTzXwJhS4hOoeKqrA0ORdylRZKzYPPbP+MWMCHo4iPDRODESX1NNm4xDZR6BK7xxjD6oMyJJ0sRYtKY7bHkYgEmBrRC8ujpW0eGyb2gUKX2DXGGJalFmBf3jWzBm4riUiA+KAeWJsYSMFrp2ifLrFrfz8os1jgAkCLSoN9edew+qDszz+ZtEsUusRuHSqqwk4zLym0pUWlQdLJUqTLqiz6uMQ6UOgSu1TfrMTS5FyLB26rFpUGS5NzUd+s5OTxCXcodIldeiu14KH23ppDi1KLt74t4LQGYnkUusTu5JbVI6O4xqTbwh6FUqNFxrka5JXVc1oHsSwKXWJ3th67BIWam2WFuynUGmw9dunPP5G0GxS6xK7UNSmRLqsy6qSZKWkZcEhWhbomWtu1FxS6xK4kZ5fB2rbH8nhASk4512UQC6HQJXYlXVbF+Rtod5OrtLR9zI5Qa0diV2SVjUZ9PVOrUPvzp5BfyYVWfgtC9+7oOHIGJH3DjaurosGorye2g2a6xG5UNcihVBs3y2VaDYSuXdBtylp4v7ob7iOmoeb7dVDXGzdTlau1qG6QGzUGsQ0UusRuyCoaIBYa9yPPFzvCPXIqhO6e4PH4cOo3BMIOnlBUXjBqXAch3+hZOLENFLrEbjQpNTB1fydNUx1UN65C7NHLqHEYY2hSqE1UFbFmFLrEbqjUWpgyc5lGjev7/gmXwCch6uxt3FgMRi99ENtAoUvshkjIN9l2Mca0uJ62HhAI0WnMfKPH4/Fg9NIHsQ30XSbtUn5+Pg4fPgy5/P/fnHIWC0zSw5YxhtqDH0HTVA+Pcf8DnsD4TUA8Hg/ODrSZyB7Qd5m0S59//jm2bNkCPp+P0NBQhIeHw6NXP8iV/Ywe+8ZPn0BVWwbPSX8HX+RggmoBhVoLaTdXk4xFrBuFLmlXtFotfv/9d1y+fBlarRZqtRonTpzAiRMn0LVrV3T86xdQ4dFnu+qb1biV+yMgEKF803Tdxzs9swAuA0c/8riOQj66ujk+8tcT20GhS2yWVqvFqVOn8O233+KXX37B2bNnUVtbCwBwd3fXfZ6DgwM++OADLFy4EM9//r/47UrdIz+msENX9F6WZmzpBqTd3Uw+JrFOFLrEJmi1WmRnZ+sF7PXr1wEAXbp0gb+/P+bOnYtx48YhLCwMjDE4OzvD2dkZBw4cwNChQwEAUVJP5F+9aVVHgR1FfLol2I7QxZTE6jDGkJOTg9TUVBw7dkw3g2WMoXPnzvD390dkZCTGjh2Lxx57DHx+2+8Hp6SkYMSIEejatavuY3VNSgxdexgKK9qe5SDk48SyJ9HRWcx1KcQCKHQJ53JycvDtt98iKytLN4NljKFTp07w9/fH8OHDMXbsWERERNwzYB/Ggp05+KGwwiraO/J5QPSg7vh4SijXpRALoeUFYlG5ublITU1FVlYWZDKZQcDOnDkTY8eOxdChQ00SsG2ZG+mLI2erObsf7U4OQgHmRvpyXQaxIApdYjb5+flITU1FZmYmZDIZampqdAHr5+eHGTNmICEhAcOGDTNbwLYl2Nsdo/w8cPhsNadX9ogFfIzy90CQtztnNRDLo+UFYhKFhYVISUnRzWCrq6vBGEPHjh3h5+eHJ554AmPHjrV4wN5LfbMSIz48igY5d/0OOkiEyHr9L+ggEXFWA7E8Cl3y0AoLC3Uz2KKiIr2A7d+/P5544gkkJCRg+PDhVhGw93KoqAqLdp3mZJlBIhJg0+QQ2rVghyh0yX2dOXMGe/fuRVZWFs6cOYOamhpotVq4u7vrAjY+Ph4jRoyw6oAFALlcjpMnT+LIkSPYs2cPiouL8cY3J/BdYa1Fg1ciEkBcehLFu9ciJiYGsbGxGDlyJPr06WOSY8rEulHoEh2ZTIa9e/ciMzMTZ86cQXV1tS5g+/XrhyeeeAJxcXEYNWqU1Qfs3Wpra+Hl5QWRSISWlhZoNBr07t0bly9fxrLUAuzLu2aR4JWIBIgP7oERDmWIi4vT7SdWKBRYsGAB/vWvf5m9BsItqwjdqgY5ZBUNaFJqoFJrIRLy4SwWIKC7Gx2NNJNz584hJSUFGRkZuoDVaDS6Gezjjz+O+Ph4jB492uYC9l5eeuklbNu2DSqVCk5OTti2bRsmT54MxhhWH5Qh6WSpWYNXIhJgakQvLI+WAgB8fHxQUlICAHB2dsbp06fRv39/sz0+sQ6c7F6oa1IiObsM6bIqyCoboVRrIRbywRgDY7fb3PF4PN3Hpd1cESX1xIQwb9pA/gjOnz+vC9jCwkJUVVXpArZfv34YP3484uLiMHr0aAiF7XNDi0KhQENDAxwdb/8jLhaLMX78eAC3f9ZWxAQgwqczlibnQq7SmnRXg1jAh0TMx/oJwXpruG+++SaWLFkCjUYDkUiEa9euUejaAYvOdHPL6rH12CWky6rA4+GhjmI6ivhg7PYxzrmRvgimbTZtunjxIpKTk3H06FGcOXMGlZWV0Gg06NChA/r164fHH38ccXFx+Mtf/tJuA/Zu169fR2JiIrp06YLPPvsMI0eOxJQpU/DOO+8YfG59sxJvpRYgo7gGCrXGqAMUfN7tfbij/D2wZlwg3J30JwxNTU3o2bMnNmzYAG9vb0yZMgUffPABZs2a9egPSqyeRULXLD/Ifh5Yk2j4g2xPLl68iJSUFBw9elQ3g1Wr1bqAjYiIQHx8PJ588km7Cdi7nTt3DjExMRg/fjz+8Y9/gM/nQ6VSQSAQ3HfZJO+PCcIhIyYIY/6YINxvH65SqYRYfPtnWCaTITY2FpMmTcKqVavazbIO0Wf20D1UVGW2l2yOIj42TAy2i203ly9f1s1gCwsLUVlZqQvYvn37IiIiAnFxcRgzZozdBuzdjh49ikmTJmHNmjWYPXv2I41R16RESk757aWwigbI1Vo4tLEUplBr4SjkQ9rdDVFST4wP7flIS2E1NTUYN24cunfvjq+//hoSieSR6ibWy2yhy8WbE+1lu83ly5d1M9iCggJdwLq5uekF7FNPPUUBew9ffvkl3nrrLezatQujRz96n9u7VTfIIatsRJNCrXvPwdlBCGk3V5O96SuXyzFnzhxcuHAB33//Pbp162aScYl1MEvoMsYsvw0nqAfWJgbaXPCWlJQgJSUFR44cQUFBASoqKnQB6+vri6FDhyImJgZPPfWU7mUouTetVov/+Z//QUpKCg4cOAB/f3+uS3okjDGsWrUKX375Jfbv34/AwECuSyImYpbQXXWgCDvNPMO9W+uMd0VMgMUe82GVlpbqAjY/Px+VlZVQqVRwdXXVzWCjo6PxzDPPUMA+gubmZkyfPh01NTVITU1Fly5duC7JaDt37sTixYuxfft2PPvss1yXQ0zA5KFLRytvKy8vR0pKCg4fPoz8/HxUVFToAtbX1xdDhgxBTEwMnn32WQpYE6ioqEB8fDykUim2bt0KBwfT3F1mDX799VeMHz8ey5cvx8svv8x1OcRIJg1da2kikvnaaIvuaigvL8fevXt1AXvt2jWoVCq4uLjoAjY6OhrPPvusbp8oMZ28vDzExcVh3rx5WL58uc0tMT2IS5cuITY2FlFRUdiwYQOt5dswk4buizuyraJd3pPSrvhsaphZxr927Rr27t2L9PR05OXloaKiAkqlEi4uLvDx8dEFbHR0NAWsBRw4cACzZs3Cxx9/jOeff57rcsyqvr4eEyZMgFgsxq5du+DqSrcH2yKThW5uWT0mbz1hFY2hJSIBds0danSf0oqKCt0SQV5eHq5du6YXsOHh4YiJiUFMTAwFrIUxxrBp0yasXbsWqampujvQ2juVSoWXX34Z//u//4u0tDT06tWL65LIQzJZ6Fr7FSgKhQIrVqzA888/j/DwcIOvqays1JvBXr16FUqlEs7OzvDx8cFjjz2G6OhoxMTE0N5JjqnVarzyyivIzMxEWloa+vTpw3VJFsUYw8aNG7F+/Xp89913eOyxx7guiTwEkywM1TUpkS6rsorABQAtAw7JqlDXpERHZzHKy8vx7LPPoqioCDweD97e3rqAzc3N1QvYPn36YOTIkXj22WcRGxsLJycnrp8OuUNDQwOef/55aLVa/Prrr+jQoQPXJVkcj8fDkiVL0K9fP0RHR2Pz5s147rnnuC6LPCCTzHQ/z7qIDenFVnet9dIx/uhUfRrTpk2DXC4HcPsHtrWdXu/evREeHo7o6GjExcVRwFq5kpISxMbGIjIyEh999BG9mYTbl3omJCRgwYIFePPNN9vlm4jtjUlCd+KW4/jtSp1RY1zf/0/Ir+RBq5JD4NwRbkOfg2vQ00aNyb9+EZe3vaL3MQcHB9TU1NCbEDbmxIkTSExMxJtvvolFixZRuNzh6tWriIuLQ3BwMDZv3kxbEK2cSTpqyCobjR7DbegEeL34JXotSUbX8W+jPuu/UFReMGpMSff+KC0tRVJSEsaPHw93d3coFArdrJfYht27dyMuLg5btmzBK6+8QoF7Fy8vL2RlZaG2thZPPfUUbty4wXVJ5D6MDt2qBjmUauOXFcQevcETtl7QxwMPPKjrKowaU8UAhw4emDJlCpKTk3Hjxg2Ulpa2i5NK9oAxhtWrV+P111/HoUOHEBcXx3VJVsvFxQWpqakIDw/H0KFDcf78ea5LIvdg9KKYrKIBYiEfChMEb+1Pn6Kp4DCYWgGxZ19I+hruMngYDkI+ZJWNukYkrW+iEeunUCgwb948nDlzBidOnECPHj24LsnqCQQC/POf/4Sfnx+GDx+OPXv2YOTIkVyXRe5idOg2KTUw1fmKzk+/hE5j/gbF1bOQlxaAJzDuamrGGJoU3J2OI4+mtrYW48aNQ5cuXZCZmQlnZ2euS7Ip8+bNg4+PDyZMmIAPP/wQM2fO5LokcgejlxdUai1M2b2BxxfA0XsgNI3X0Xj6oFFjMQaTLH0Qyzl37hyGDh2Kxx9/HCkpKRS4j2jMmDHIyMjAe++9hxUrVkCrpb8H1sLo0BUJ+TDL+xpardFrujweIBbqP0WFQmHUmMR8MjIyMGLECCxbtgzr1q2jmxOMFBAQgJMnT+qaube0tHBdEoEJQtdZLDD63WRNUz2aijKhVbaAaTVouZSNJlkmHHsHGTUuj8eDmM9w5MgRLF26FP3794dEIkFNTY1R4xLT+/LLL/H888/jm2++wZw5c7gup93w8PDA4cOHIRKJMHr0aFRVVXFdkt0zek1X2t3N+JfwPB4aT/+A2p8+BZgWwg5d0fHJuXDyM+48vUKlwTNDA8Ga6/+4XoXB0dERnTt3Nq5eYjJ3Nh3PzMzEgAEDuC6p3XF0dMSOHTvw/vvvIyIigpqic8wkhyMC3/sJjRy2c7wXN0ch5nYqxuuvv65bVhAIBOjatSvCw8MRFham+z9diWJ5rU3Hq6ur8e2339JWPgtISkrCq6++iq+//hrPPPMM1+XYJZMsmkm7WefpLml3NyxcuBDZ2dnw8vICn8/Hyy+/jOPHj2PGjBmQy+X497//jYCAAHh5eSEhIQHvv/8+Dhw4QC/DzKyiogIjR46Ek5MT0tPTKXAtZOrUqfj222/xwgsv4JNPPuG6HLvU7nsvzI30BXC7Ucr8+fMxf/58jBgxQu9zGWO4cuUKsrOzcerUKWRnZyM7OxtOTk66mXDrf56e3N9KYevy8vIQHx+PuXPnttum49bu0qVLurv3NmzYAIFAwHVJdsMkoVvXpMTQtYdNckDCVByEfJxY9uQjXYMN3A7iy5cvGwSxi4uL3rJEWFgYunbtauLq268DBw7ghRdewKZNm9p903FrV19fj/Hjx8PR0RHffPMN9SOxELvpp2sKjDFcunRJL4hzcnLg6upqEMQeHh4mfez2YNOmTVizZg327t2Lxx9/nOtyCG43RV+wYAFOnjyJ/fv3U1N0C6CbI4zEGMPFixd1M+FTp04hJycHHTp0MAhie123VKvVWLx4MY4ePYq0tDT4+PhwXRK5A2MMGzZswIYNG6gpugXQHWlmoNVq2wzijh076gK4NYzb+/a1O5uO79mzxy6bjtuK7777DnPnzqWm6GbWLm8Dznr9L+ggMa5vg6lptVpcuHBBL4hPnz6NTp06GQRxp06duC7XJFqbjg8fPhybNm2ipuM2gJqim59JQxcADhVVYdGu05wsM0hEAmyaHIIoqW3sMGgN4tb14dYg7ty5s94+4tDQUJsL4pMnT2LcuHF44403qAeujaGm6OZl8tAFgL8fKELSyVKLBq9EJMDUiF5YERNgscc0B61Wi/Pnz+vtmDh9+jS6dOliEMQdO3bkutw27dmzBwsWLMCXX35JPXBt1K1btzB16lQ0NDRg7969NvePvjUzS+gyxrAstQD78q5ZJHglIgHig3tg7bjAdjmj0mq1KC4uNgjiO0/WhYWFcR7EjDGsWbMGmzdvxr59+xAcHMxZLcR4Go0Gb7zxBtLS0pCWlob+/ftzXVK7YJbQBf7o+n9QZvYZb+sMd3m0tF0G7r1oNBoUFxfrbV/Lzc2Fp6en3vpwaGgo3N3dzV6PQqHA3/72NxQWFmLfvn3UdLwd2bJlC1auXIk9e/YYHCwiD89sodvqUFEVlibnQq7SmnRXg1jAh0TMx/oJwTazhmtuGo0G586d0wvivLw8dOvWzSCITbmLoLa2FomJiejUqRN27NhBPXDboUOHDmHq1Kn45z//iRkzZnBdjk0ze+gCt3c1vJVagIziGijUGqMOUPB5gINQgFH+HlgzLhDuTrTIfz8ajQZnz57V2zWRl5eHHj16GASxm5vbQ49fXFyMmJgYJCYmYs2aNdQDtx0rKipCbGwspkyZgvfff5++14/IIqHbKq+sHluPXcIhWRV4PDxUrwZHER+MAWOknpgb6WuRgw/tlVqtNgji/Px8eHl56QVxSEjIfYM4IyMDzz//PFavXo2//vWvFnwGhCvV1dUYO3YsvL298dVXX0EikXBdks2xaOi2qmtSIiWnHOmyKsgqGiBXa+Eg5P/R8/b2jQ88Hg8KtRaOQj6k3d0QJfXE+NCej9xLgdyfWq2GTCYzCGJvb2+9fcQhISFwdXXFf/7zH7z55pv45ptv8OSTT3JdPrEguVyO2bNn49KlS/j++++pCdRD4iR071bdIIesshFNCjWUai3EQj6cHYSQdnPV3eRLLK81iO/cR5yfnw+JRAKlUokFCxYgJiYGISEhcHFx4bpcYkGMMbz33nv46quvkJaWhkGDBnFdks2witAltqG5uRnTpk3DlStXMHPmTN3uiYKCAvTu3VtvRhwcHExBbAeSkpKwePFi/Pe//6Wm6A+IQpc8kIqKCiQkJMDf3x/btm2Dg4OD7vdUKhWKior09hEXFhaid+/eevuIQ0JCaGdDO/TLL79g/PjxeOedd/DSSy9xXY7Vo9Alfyo/Px9xcXH461//ihUrVjzQfmiVSoUzZ87obV8rLCyEj4+PXhAHBwdTELcDFy9eRGxsLJ5++mmsX7+emqLfB4Uuua+DBw9i1qxZ+OijjzBp0iSjxlIqlQZBfObMGfj6+urtmggODoaTk5OJngGxlLq6OowfPx4SiYSaot8HhS65p02bNuEf//gHUlNTzdZ0XKlUorCwUG/XRFFREfr27asXxEFBQRTENkClUuGll17Cb7/9hrS0NHh7e3NdktWh0CUGuG46rlAoDIJYJpOhX79+BkFM+0StD2MM69evx8aNG6kpehsodIme1qbjGo0GycnJVtN0XKFQoKCgQC+Iz549i/79++sF8eDBgymIrURrU/QtW7YgMTGR63KsBoUu0bmz6fhHH30Ekci6GsHfrTWI79w1cfbsWfj5+eltXxs8eDAcHWm/NxdycnIQHx+PhQsX4o033rCrplT3QqFLALSfpuNyudwgiM+dOwd/f3+9IA4MDKQgtpDy8nLExcUhNDQUn332md03RafQJUhOTsZLL73UbpuOy+Vy5Ofn6wVxcXExBgwYYBDEd+4/JqZz69YtTJkyBbdu3cLevXuttgG/JVDo2jHGGNauXYvPPvvM7pqOt7S0ID8/X2/72vnz5zFgwAC9GzoGDRpEQWwiGo0Gr7/+Og4cOIADBw6gX79+XJfECQpdO6VUKjFv3jwUFBRg//791HQct4M4Ly9PL4gvXLgAqVSqd6AjMDDQ7l8iG2Pz5s149913kZycjMjISK7LsTgKXTtETccfXHNzs0EQX7x4EQEBAXpBPGjQIArih/Dzzz9j2rRpdtkUnULXzrQ2HR83bhzWrl1LjagfQXNzM3Jzc/W2r12+fLnNILb2HSBcOnPmDGJjYzF16lS7aopOoWtHqOm4+TQ1NRkE8ZUrVzBw4EC9IB44cCAF8R2qq6uRkJCAXr16mbQpelWDHLKKBjQpNVCptRAJ+XAWCxDQ3Y3zdrEUunbiP//5D5YtW4adO3dS03ELaQ3iO3dNXLlyBYMGDdI70BEQEGDXQdzS0oLZs2fj8uXLj9wUva5JieTsstsXI1Q26vpy330xQuvHpd1cESX1xIQwb4tfjECh285ptVqsWLECe/bsQVpaGgYMGMB1SXbt1q1bBkFcUlKCwMBAgyAWCoVcl2sxjDG8++67+Prrr7F///4Hboqe+8cVYOlGXAEW9ccVYMEWugKMQrcda25uxsyZM1FRUYHvvvsOXbp04bok0obGxkaDIC4tLcXgwYP1glgqlbb7IN6xYweWLFmC//73v3j66afv+XlmuezWzwNrEs1/2S2FbjtVWVmJ+Ph4+Pn54YsvvqC9pjamsbERp0+f1ts1UV5ergvi1jAeMGBAuwvi1qboK1euxIsvvmjw+4eKqrA0ORdylRZKzYPPbP+MWMCHo4iPDRODESU1371vFLrtUH5+PuLj4zF79my8/fbbNnukl+hraGgwCOKrV68iKCjIIIhtvYn4xYsXERMTg2eeeUbXFJ0xhtUHZUg6WYoWlcZsjy0RCTA1oheWR0vN8neHQredOXjwIGbOnImPPvoIkydP5rocYmY3b940COKKigqDIPb397e5IG5tiu7k5ISkpCSsPnQF+/KumTVwW0lEAsQH9cDaxECTBy+FbjvS2nR87969GDZsGNflEI7cvHkTOTk5etvXKisrERQUpHfE2c/Pz+qDWKVSYcGCBbjSMRRlDn0sEritWme8K2ICTDouhW47oFar8eqrr+LIkSOcNB0n1q++vh6nT5/WzYZPnTqFqqoqBAcHGwSxtR1SOFRUhUW7Tls0cFtJRAJsmhxi0jVeCl0b19DQgEmTJkGtVltV03Fi/err65GTk6MXxDU1NXpBHBYWxmkQ1zcrMeLDo2iQqzl5fADoIBEi87XRJtvVQKFrw0pLSxEbG4thw4Zh06ZNdr3BnphGXV2dXhBnZ2ejpqYGISEhekHcv39/iwTxizuycfhstUl3KTwssYCPJ6Vd8dnUMJOMR6Fro3777TeMGzcOr732GhYvXkw7FIjZ3LhxQ7dG3BrGtbW1BkHcr18/kwZxblk9Jm89wcmywt0kIgF2zR2KIBMcoKDQtUEpKSl48cUX8cUXXyA+Pp7rcogdqq2tNQjiuro6hISE6B3o6Nu37yMH8YKdOfihsMKogw+mwucB0YO64+MpoUaPRaFrQ1qbjn/66afYt28fQkJCuC6JEJ3a2lq9HRPZ2dmor69HaGioQRD/2SuzuiYlhq49DIWau2WFuzkI+Tix7EmjezVQ6NqI1qbj+fn52L9/P7y8vLguiZA/df36dV0Qt4ZxQ0ODQRD7+vrqBfHnWRexIb34oXopmJujiI+lY/wxN9LXqHEodG1AbW0tnnvuObi7uyMpKYmajhObVlNTYxDEjY2NurXhsLAw7L7eA9llNx/5MRqy96Op4DCUNVfgLB2JLrGvmqT2CJ9O2D3vcaPGoNC1cufPn0dMTAwSEhKwdu1aq9/MTsijqK6u1u2ayMnJgWzgHLQYsUus+dxxgMdDy+UcMJXSZKHr5ihE/sp7N+J5EBS6ViwzMxPPP/88Vq1ahblz53JdDiEWUdUgx4gPj5pkPbcu67/QNFw3WeiKhXz88vpooxqhW9fRE6Lz1VdfYeLEiUhKSqLAJXZFVtEAsdA6o8lByIesstGoMdpXT7h2QKvV4u2338auXbuQkZEBqVTKdUmEWFSTUgNrfQHOGEOTwrjTcRS6VqSlpQUzZsxARUUFTpw4AQ8PD65LIsTiVGotrDRzwRigNHLZwzrn8HaosrISo0aNglgsRnp6OgUusVsiIR/WesCSx4PRSx8UulagoKAAQ4cORUxMDHbs2AFHR25vKyWES85igdHH2plWA6ZWAloNwLRgaiWY1vjjxDweD84Oxi0Q0PICx3744QfMmDED//73vzFlyhSuyyGEE4wxXL58GdnZ2fglpxBy/hCjxrv56y7c/PUb3a+bzhxFhycmwz1yqlHjKtRaSLu5GjUGhS6HPv74Y6xevRrff/89NR0ndoMxhitXrugdjsjJyYFEItEdjnDQ8KFSPvraqXvkVKMDti2OQr5R28UACl1OqNVqLFmyBOnp6Th+/Dg1HSftFmMMJSUlBj0ZHB0ddceAFy9ejLCwMHTr1k33dYVbjuO3K3UcVt42aXc3o8eg0LWwxsZGTJo0CSqVCsePH4e7uzvXJRFiEowxlJaW6oVrdnY2xGKxLmAXLVqEsLAwdO/e/b5jRUk9kX/1ptX1XjDFDRJ0Is2CqOk4aS8YYygrKzMIWKFQqNfIJiwsDD169Hjo8dtzlzGa6VrI77//jrFjx2Lp0qV49dVXqek4sRmMMZSXl+sF7KlTpyAQCHQB+9JLLyE8PPyRArYtHZ3FiJJ6WlU/3TFST6MDF6DQtQhqOk5sBWMMV69eNQhYHo+nm73Onz9fF7DmnDzMjfTFkbPVVnFzhINQYHRLx1a0vGBGdzYd//777xEaanzXeUJMhTGGa9euGQQsY0xveSA8PBxeXl6cvDqjO9LIA1Mqlfjb3/6GvLw8ajpOrEJbAavRaPSuYA8LC0PPnj2tZvnLWm4Dznr9L+ggMc17MLS8YAY3btxAYmIi3N3dkZWVBRcXF65LInamoqJC7w2uU6dOQaVS6YJ1zpw5+PTTT+Ht7W01AdsWJyEQrpYhC75Qw/K9pCUiAdZPCDZZ4AI00zU5ajpOLK2ystIgYJVKpd7yQFhYGHr16mXVAQvc3sP+22+/4ejRo9i3bx9+//138Pl8LNnxC/YX1Vt0fVciEmBqRC+siAkw6bgUuibU2nT8/fffx7x587guh7RDlZWVBgcN5HK5QcD27t3b6gO2Ldu3b8esWbMgEomgUqkAAOvWrcPrr7+OZakF2Jd3zSLBKxEJEB/cA2vHBZr8z5FC10S2b9+O119/HTt37kRUVBTX5ZB2oKqqyiBgm5ubDfbB9unTxyYDti0qlQqhoaEoLCwEADg7O6O6uhpOTk5gjGH1QRmSTpaaNXhbZ7jLo6Vm+XOlNV0jabVavPPOO/jmm2+QmZlJTcfJI6murjYI2Fu3bukCdtq0adi4cSN8fHzaTcC25YsvvsCZM2fg5OQEhUKBOXPmwMnJCcDtDl8rYgIQ4dMZS5NzIVdpTbqrQSzgQyLmY/2EYJOcPLsXmukaoaWlBTNnzsTVq1fx3XffUQ9c8kDuvA23NWAbGxsRGhqqN4O9+1ry9m7BggX47LPP8N5772H69Ol49tln8cMPP6BPnz4Gn1vfrMRbqQXIKK6BQq0x6gAFn3d7H+4ofw+sGRcIdyfjD0DcD4XuI6qsrERCQgL69euHL774gnrgkjZdv37dIGBv3rxpELB9+/a1q4C9k1arxZNPPoljx45hz549SExMfOCvzSurx9Zjl3BIVgUeDw/Vq8FRxAdjt0+azY30RZC3+yNU//AodB9BQUEB4uLi8MILL+Cdd96x278sRF9tba1BL4K6ujq9gA0PD4evry/4fLo/AADq6+sRHByM69ev4/jx4xg8ePAjjVPXpERKTjnSZVWQVTRArtbCQcgHYwyM3b7xgcfjQaHWwlHIh7S7G6Kknhgf2tMkR3sfBoXuQ/rxxx8xY8YM/Otf/6Km43bsxo0bBgF748YNhIaG6u0k6Nu3LwXsPchkMkRERKBDhw7Iy8tDp06dTDZ2dYMcsspGNCnUUKq1EAv5cHYQQtrN1eh+uMai0H0In3zyCf7+978jJSUFTzzxBNflEAu5ceMGcnJy9E5y1dbWGgRsv379KGAfUFpaGsaNG4eIiAhkZGRAKLSf9/QpdB+ARqPBkiVLcOjQIaSlpcHX1zSNL4j1qaurMwjY69evIyQkRC9g+/fvTwH7iD744AMsW7YMc+bMwdatW7kux+IodP9EY2MjJk+eDIVCgeTkZGo63o60Buydb3RVV1frArZ1HdbPz48C1kSmT5+OnTt3YsOGDXjllVe4LocTFLr3UVpairi4OAwdOhQff/wxNR23YfX19QYBW1VVheDgYIOApaPbpqdUKjFs2DDk5+fjwIEDGDNmDNclcYZC9x6o6bjtunnzpkHAVlZWIigoSC9g/f39KWAtoLq6GkFBQWhpaUF2djb69u3LdUmcotBtw969ezF//nxs27YNCQkJXJdD7qOhoUEXsK3rsNeuXTMI2AEDBlDAcuDUqVMYMWIEvLy8cPr0aeq4BwpdPYwxrFu3Dp988gk1HbdCDQ0NOH36tF7AXr16FYMHD9Y7aDBgwAC7ejfcWn3zzTeYPn06oqKicPDgQVoX/wOF7h+USiXmz5+P3NxcajpuBRobGw0CtqyszCBgpVIpBawVWr58OdasWYNXX30V69ev57ocq0Khi9v7MJ977jm4ubkhKSmJXgJZ2K1bt3D69Gm9gwalpaUIDAzUC9iAgAAKWCun1Woxbtw4pKWlYdu2bXjhhRe4Lsnq2H3onj9/HrGxsYiLi8O6deto3c/Mbt26hdzcXL2ALSkpwaBBg/SOykqlUtotYmOam5sRHh6Oy5cv48iRI3j88ce5Lskq2XXoZmVlYeLEidR03EyampoMAvbKlSsYNGiQ3kGDgIAAClgbV1JSgtDQUAgEAuTk5KBnz55cl2S17DZ0v/76a7z22mtISkqy6z2DptIasHdu07p8+bJBwA4cOJACtp3JzMzEU089hQEDBuDkyZPUce9P2F3otjYd37lzJ9LS0hAQYNr7j+xBc3OzQcBeunQJAwcO1NumNXDgQIjFlu3gRCxry5YteOmll5CYmIjk5GSuy7EJdhW6rU3Hy8vL8d1336Fr165cl2T1mpubkZeXp7eL4OLFiwgICNAL2EGDBlHA2pmFCxfik08+wcqVK7Fy5Uquy7EZdhO6VVVVSEhIgK+vL7788kt6CdSGlpYWg4C9cOECpFKpQcA6ODhwXS7hyJ1Nx3fv3o3nnnuO65JsilWEblWDHLKKBjQpNVCptRAJ+XAWCxDQ3c0kvS8LCwsRGxtLTcfv0NLSgvz8fL2APX/+PAYMGKAXsIGBgRSwRKe+vh4hISGoqakxqum4PeMkdOualEjOLrvd5b2yUddk+O4u760fl3ZzRZTUExPCvB+6y3tr0/GNGzdi6tSpZnpG1k0ulyM/P19vF0FxcTH8/f319sEGBgbSKwByT3c2Hc/NzUXnzp25LskmWTR0c/+4zyjdiPuMov64zyj4Ae4z+vTTT7Fq1Sq7ajoul8tRUFCgF7Dnzp2Dn5+f3j5YCljyMA4ePIiEhAQMGTIEmZmZdEjFCBYJXbPc3OnngTWJbd/cqdFosHTpUvz00084cOBAu206rlAoDAL27Nmz8PPz09umNXjwYApY8sham47Pnj0b27Zt47ocm2f20D1UVGW2O+odRXxsmKh/R31r03G5XI6UlJR203S8NWDv3KZ19uxZ9O/f3yBgJRIJ1+WSdmLGjBnYsWMHNm7caLdNx03NbKHLGMPqgzIknSxFi0pjjocAAEhEAkyN6IXl0VKUlZUhLi4OERER+OSTT2x2E75SqTQIWJlMhn79+ukFbFBQEAUsMQuVSoUnnnhC1wDq6aef5rqkdsMsocsYw7LUAuzLu2bWwG0lEQkworcEae9MwZIlS7BkyRKb2aGgVCpRWFioF7BFRUXo27ev3i6CoKAgODk5cV0usQN3Nh0/deoU+vXrx3VJ7YpZQnfVgSLsNPMM924OAh4iuwPbFkRb7DEflkql0gVs6zrsmTNn4OvraxCwzs7OXJdL7FB2djYiIyPRo0cP5ObmUsc9MzB56B4qqsKiXactGritJCIBNk0O0Vvj5YpKpcKZM2cMArZPnz56ARscHEwBS6zCrl27MG3aNGo6bmYmDd36ZiVGfHgUDXK1qYZ8aB0kQmS+NrrNXQ3molKpUFRUpLeLoLCwEL1799bbBxscHEwzB2KV3n77baxevRqLFy/Ghg0buC6nXTNp6L64IxuHz1abdJfCwxIL+HhS2hWfTQ0zy/hqtdogYAsKCtCrVy+9gA0JCaGAJVbvzqbjW7duxezZs7kuqd0z2Q7n3LJ6ZBTXcBq4AKDUaJFxrgZ5ZfUIeoADFPejVqshk8n0AjY/Px/e3t66gJ00aRKCg4Ph6upqmidAiIU0Nzfjsccew8WLF5GVlWU3B4i4ZrLQ3XrsEhRqy6/jtkWh1mDrsUv4eMr/Xyyp0WjwwQcfIDY2FoGBgQZf0xqwd+4iaA3Y1tnrxIkTERISQgFLbF5JSQnCwsLA5/Nx4cIFajpuQSYJ3bomJdJlVUadNDMlLQMOyapQ16RER2cxampqkJCQgJMnT6KxsREBAQFtBqyXl5cuYMePH4+QkBC4ublx/XQIMamsrCyMGTMG/v7++O233+i0ooWZZE3386yL2JBe/FC9FMzNUcTH0jH+GCSsQkxMDBoaGqDRaODq6gqtVosePXro7SIICQlBhw4duC6bELP6/PPP8eKLL1LTcQ6ZZKabLqsyWeCqblzFtS9ehvOAJ9Al7rVHHkeu0uJAbinmLRqt93GNRoPy8nJ07NjR2FIJsSnUdNw6mCR0ZZWNphgGAHDj581w6N7fJGNduqFAUVERjhw5gtTUVBw/fhzNzc1Qq7nb0kaIpWm1WkRFRSErKwt79uzB+PHjuS7JrhkdulUNcijVppnlNhVlgu/oDFHnAVDXVxg9nlytRWcvHyxYsAALFiyASqXC2bNn0aVLFxNUS4j1a2hoQHBwMKqrq5GdnY2goCCuS7J7Rh85kVU0QCw0/uSKVtGM+mNJ6PiXOUaP1cpByNebhYtEIgQGBtpMXwZCjHH27Fl4e3tDqVTiypUrFLhWwui0bFJqYIrzFfVZ/4VL0FMQunkYPVYrxhiaFLSUQOzPwYMHERgYiIEDB+LKlSv06s6KGB26KrUWxmausuoS5CV5cHsswdhy9DAGky19EGIrPvzwQ8TGxmL69Ok4fvw43fJgZYz+boiEfBj7al1eWgD1zSqUf/oCAIAp5QDTouL6K+j+wr8feVweDyZZ+iDEVrQ2Hd+wYQMWL17MdTmkDUaHrrNYYPQaqUvw03CWjtD9uuG3VKhvVqHT0wuMGpfH48HZgf6VJ+2fWq3GsGHDkJubix9++IGajlsxoxNJ2t3N6JfwfJEjIPr/UzE8kSN4QjEETsYdVlCotZB2oyO7pH2rrq5GcHAwmpubcebMGfTvb5otl8Q8jA5dTzdHiIV8KEy4duoeaZqr0h2FfHR1oyOOpP3KycnB8OHD0aNHD5w7d476gtgAkyx4WutsUtqd+iaQ9mv37t0YMmQIIiMjUVxcTIFrI0wSulFSTziKrOsNK0cR3ypukCDEHN5++21MnjwZCxcuxE8//US3PNgQk3ynJoR5G71tzNQYA8aHUrs60r5otVqMHTsW//jHP7B161Zs3LiR65LIQzLJW/sdncWIknrih8IKq2jvyOcBY6Se6OhsuSt7CDG3lpYWhIeHU9NxG2ey1yRzI33hIBSYajijOAgFmBvpy3UZhJhMaWkpvL29UV1djQsXLlDg2jCThW6wtztG+XlALOB2bUks4GOUv4fRV/UQYi2ysrLg5+eHHj16oKysjG55sHEmTcg1iYGcv6EmEfOxNnEwpzUQYiqff/45Ro8ejdjYWOTn59MtD+2ASRPS3UmM9ROCIRFxs8wgEQmwfkIwOkhEnDw+Iaa0aNEizJ8/H2+//TZSUlK4LoeYiMnPyI4J8MTUiF5IOlmKFpXlLqqUiASYGtGLtokRm6fVajFmzBhkZmZS0/F2yCyNCZZHS9EoV2Nf3jWLBK9EJEB8cA8sj5aa/bGIfatqkENW0YAmpQYqtRYiIR/OYgECuruZ5PTjnU3HT506heDgYOOLJlbFLKHL4/GwNjEQro5Cs894W2e4y6Ol1JycmFxdkxLJ2WVIl1VBVtkIpVoLsZAPxhgYu93Jjsfj6T4u7eaKKKknJoR5P/SWxXPnzmHIkCFwdXWlHrjtmEluA76fQ0VVWJqcC7lKC6XGdP0ZxAI+JGI+1k8IpiUFYnK5ZfXYeuwS0mVV4PHwUBevOor4YOz2Sc25kb4IvsdOGpVKhYiICKxfvx4tLS1ISEhAeHg4jh07Rj1w2zGzhy4A1Dcr8VZqATKKa6BQa4w6QMHn3d6HO8rfA2vGBcLdiQ5AENMxy8+qnwfWJBr+rO7evRszZ84EYwwqlQqzZs3Cl19+aeQzINbOIqHbKu+P2cMhI2YPY/6YPdA+XGJq5nxV5ijiY8NE/VdlwcHByMvLAwC4ubnh4sWLtKRgBywauq3qmpRIySm/vU5W0QC5WguHNtbJFGotHIV8SLu7IUrqifGhPeloLzE5xhhWH5RZ9P2H7OxsDBkyBIwx8Hg8CAQCzJ07F59++qnZHp9YB05C927VDXLIKhvRpFDr3pBwdhBC2s2V+uESs2KMYVlqgUV32jw70ANfL3gatbW16NevH+Lj4xEVFYXIyEi4uLiYvQbCLasIXUK4supAEXZaeE+5mM/QtbEYqSumoWvXrhZ7XGIdqAknsVuHiqosHrgAoNTyUOsegPxamu/YIwpdYpfqm5VYmpxr8cBt1aLSYGlyLuqblZw8PuEOhS6xS2+lFjzU7hlzaFFq8da3BZzWQCyPQpfYndyyemQU15h0W9ijUGq0yDhXg7yyek7rIJZFoUvsztZjl6BQc7OscDeFWoOtxy5xXQaxIApdYlfqmpRIl1VZxbVSAKBlwCFZFeqaaG3XXlDoEruSnF0Ga+uLxOMBKTnlXJdBLIS6ahC7ki6rMvoNtMqkZVBcOwce/3azfoFrZ3jN2/LI48lVWqTLquhePztBoUvsiqyy0STjdHpqPlyDnjbJWAAgq2gw2VjEutHyArEbVQ1yKNXc7li4F7lai+oGOddlEAugmS6xG7KKBoiFfChMELz1GdtRn7Edok5ecB8xHY69jbsM1UHIh6yykXqN2AEKXWI3mpQamKLVSMfRL0DU2Rs8gQhNsixU712F7i98BFHH7o88JmMMTQq10bUR60fLC8RuqNRamKK9k0MPf/AdnMATiuAS+CQcvKRouXjKqDEZg9UufRDTotAldkMk5JtnuxiPB8C4NOfxALGQ/jraA/ouE7vhLBYYfXmpVn4LLZeywdRKMK0Gt84chaKsEBKfUKPG5fF4cHag1T57QN9lYjek3d2MfgnPtBrUZ+2A6kY5wOND1LknPBJXQNS5p1HjKtRaSLu5GjUGsQ0UusRueLo5Gr17QeDUAd1nbTRhVbc5Cvm0c8FO0PICsSvWOpuUdnfjugRiIRS6xK5EST3hKLKuH3tHEV/vlmDSvlnXTx8hZjYhzNsk28ZMiTFgfKhxa8LEdlDoErvS0VmMKKkn+FbSaYzPA8ZIPdHRWcx1KcRCKHSJ3Zkb6QsHoYDrMgAADkIBdRezMxS6xO4Ee7tjlJ8HxAJuf/zFAj5G+XsgyNud0zqIZVHoEru0JjGQ8zfUJGI+1iYa1yiH2B4KXWKX3J3EWD8hGBIRN8sMEpEA6ycEo4NExMnjE+5Q6BK7NSbAE1Mjelk8eCUiAaZG9KJtYnaKQpfYteXRUsQH9bBY8EpEAsQH98DyaKlFHo9YHx4zRYNRQmwYYwyrD8qQdLIULSrzXc3eOsNdHi01uvEOsV0UuoT84VBRFZYm50Ku0kKpMV1vW7GAD4mYj/UTgmlJgVDoEnKn+mYl3kotQEZxDRRqDbRG/O3g827vwx3l74E14wLh7kQHIAiFLiFtyiurx9Zjl3BIVgUeDw91bbujiA/Gbp80mxvpS/twiR4KXULuo65JiZSccqTLqiCraIBcrYWDkA/GGBi7feMDj8eDQq2Fo5APaXc3REk9MT60Jx3tJW2i0CWEEAuiLWOEEGJBFLqEEGJBFLqEEGJBFLqEEGJBFLqEEGJB/werA5L3CtNhPwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAD3CAYAAAC+eIeLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9NUlEQVR4nO3de3zMV/4/8Ndcc5OISwQRl7hOiEwuK2rrtmK76xKbUJeiqNXVuvRCWzb0Zi3aL7bVomh3uz/UrbZ1iXalmsTWUnIPQxS5iGQSkUgkmevn/P7QmTUmlMxn5jOX9/Px8IdJ5sx7iJfPnM857yNijDEQQghxCLHQBRBCiCeh0CWEEAei0CWEEAei0CWEEAei0CWEEAei0CWEEAei0CWEEAei0PVA8+fPx6pVq4Quw0paWhq6dOkidBm86N69O1JTU4UugzghCl030717d/j4+MDf3x+BgYEYMmQItm7dCo7jzN+zdetWrFy5kvfXnj17NlasWMH7uCYikQgREREW72XFihWYPXv2Iz1/xIgR2LFjh52qs41Wq8XcuXPRrVs3+Pv7IyoqCseOHRO6LGIHFLpu6PDhw6ivr0dxcTGWLVuGdevWYe7cuY/0XIPBYOfqbHPjxg3s2bNH6DIeqKV/fgaDAaGhoUhPT8ft27exatUqTJ48GUVFRfwWSARHoevGWrdujYSEBOzduxeff/45CgoKAFhekZo+0q9btw4dO3bEnDlzwHEc1q5di549e6Jdu3aYPHkybt26ZR73P//5D4YMGYLAwECEhobiH//4B7Zt24Zdu3bhvffeQ6tWrTB+/HgAd0Ny4sSJCAoKQo8ePfDhhx+ax2lqasLs2bPRpk0bhIeH4+zZs7/4nl5//XW89dZbDwy306dPm2uLjIxEWloaACA5ORknT57EwoUL0apVKyxcuBBvvfUWFi1aBADQ6/Xw8/PD66+/bq7N29sbNTU1AIBDhw6hf//+CAwMxIgRI6BSqcyv2b17d6xbtw4DBw6En5+fVW0XL15Ejx49HvqfhZ+fH95++210794dYrEY48aNQ48ePZCZmfmLfybExTDiVrp168aOHz9u9XhoaCjbvHkzY4yxWbNmseTkZMYYY99//z2TSCTs9ddfZxqNhjU2NrKNGzeyuLg4VlpayjQaDXv++efZ1KlTGWOMFRcXs1atWrHdu3cznU7Hbt68ybKzs63GZYwxo9HIoqOj2TvvvMO0Wi27cuUK69GjB/vmm28YY4y98cYb7Mknn2TV1dWspKSE9e/fn4WEhDzwvQFghYWFLDo6mm3fvp0xxlhycjKbNWsWY4yx69evs7Zt27KjR48yo9HI/v3vf7O2bduyyspKxhhjw4cPNz+PMca+++47NmDAAMYYYz/88AMLCwtjgwYNMn9t4MCBjDHGLl26xHx9fdm///1vptPp2Lp161jPnj2ZVqs1/5lHRkaykpIS1tjYaPH3kJmZyUJDQ9nhw4cf6e/PpKKignl5eTGVSvVYzyPOj650PUTnzp0trlbvJRaL8c4778DLyws+Pj745JNPsHr1anTp0gVeXl54++23ceDAARgMBuzatQvx8fGYNm0aZDIZ2rVrB6VS2ey4Z8+eRVVVFd58803I5XKEhYVh3rx55iu+ffv2ITk5GW3btkVoaCgWL178i+9DJBJh1apVePfdd6HVai2+tnPnTowZMwZjxoyBWCzG6NGjERsbi5SUlGbHeuKJJ3D58mVUV1cjIyMDc+fORVlZGe7cuYP09HQMHz4cALB3716MHTsWo0ePhkwmw9KlS9HU1IRTp06Zx1q8eDFCQ0Ph4+NjfuzkyZNISEjA559/jnHjxv3iezPR6/WYPn06Zs2ahX79+j3y84hrkApdAHGMsrIytG3bttmvBQUFwdvb2/z74uJiJCYmQiz+3//JEokEarUapaWl6Nmz5yO9ZnFxMW7cuIHAwEDzY0ajEUOHDgVwd+ohNDTU/LVu3bo90rhjxoxB165dsW3bNqvX279/Pw4fPmx+TK/XY+TIkc2O4+Pjg9jYWKSnpyMjIwPJycnIycnBDz/8gPT0dPPUw40bNyxqE4vFCA0NRVlZmfmxe9+HydatWzF8+PAHvn5zOI7DzJkzIZfL8dFHHz3y84jroCtdD3D27FmUlZXhySefbPbrIpHI4vehoaE4duwYamtrzb80Gg1CQkIQGhqKK1euPPI4PXr0sBinvr7efOXZqVMnlJaWmr+/pKTkkd/TX/7yF6xevRqNjY0Wrzdz5kyL12toaMCyZcuarQ8Ahg8fjhMnTiA7Oxu/+tWvMHz4cHz77bf48ccfMWzYMAB3PyUUFxebn8MYQ2lpKUJCQh743oG7oVtSUoJXXnnlkd4TYwxz586FWq3Gl19+CZlM9mh/GMSlUOi6sbq6Ohw5cgRTp07FjBkzEBER8UjPmz9/PpKTk81BU1VVha+//hoAMH36dKSmpmLfvn0wGAyorq5GTk4OACA4OBhXr141jzNo0CAEBARg3bp1aGpqgtFoREFBgfmG2eTJk7FmzRrU1NTg+vXr2LRp0yO/txEjRiAiIgKff/65+bEZM2bg8OHD+Pbbb2E0GqHRaJCWlobr1683Wx9wN3T/+c9/Ijw8HHK53LysrEePHggKCjLXefToUXz33XfQ6/VYv349vLy8MGTIkIfW6O/vj2+++QYZGRnm4H+YF154ASqVCocPH7aYpiBuRuhJZcKvbt26MW9vb9aqVSsWEBDABg8ezD766CNmMBjM33P/jbT7b14ZjUa2fv161qdPH9aqVSsWFhbGli9fbv56RkYGGzRoEPP392ddunRh//jHPxhjjBUWFrLIyEjWunVrNmHCBMYYY2VlZWzq1KksODiYBQYGsri4OPONvoaGBjZz5kzWunVrplAo2HvvvfeLN9IuX75s/v3p06cZAPONNNNjw4YNY23atGHt27dnY8aMYcXFxYwxxk6dOsV69+7NAgMD2aJFixhjjNXX1zOpVMrefvttxhhjHMexoKAgNn/+fIvXPnjwIFMoFCwgIIANGzaMFRQUWPyZ33/z8t7Hqqur2cCBA9mKFSse+N6KiooYAObl5cX8/PzMv3bu3PnA5xDXJGKMTo4ghBBHoekFQghxIApdQhykpKQErVq1avbX49xEJK6NphcIIcSBaJ0uIQ9R06DD/sxSpKrUUFXUQ2fgIJeKwRgDY4BIdHe5mOlxRUd/xCuC8XRMKNr4yYUunzghutIlpBk5pbXYfvIqUlVqiESARs/98pN+5i0TgzEgXhGMeUPDoAwNtF+hxOVQ6BJyj9pGHZYfzEdaYRW0BiM4G/51iEWAl1SCEX2CsCYpAoG+dOVLKHQJMTt+QY0l+3Og0XPQGR/9yvaXyCVieMvE2DBZiXhFMG/jEtdEoUs8HmMMq1NU2HWmBE16o91ex0cmwfS4rkgeo2h22zDxDBS6xKMxxrDsYD4O5d6wa+Ca+MgkSIjsjLVJERS8HorW6RKP9pcUlcMCFwCa9EYcyr2B1SmqX/5m4pYodInHOn5Bjd12nlJoTpPeiF1nSpCqUjv0dYlzoNAlHqm2UYcl+3McHrgmTXojluzPQW2jTpDXJ8Kh0CUeafnB/Mdae2sPTToOy/+VL2gNxPEodInHySmtRVphFa/LwlpCZ+SQdqkKuaW1gtZBHItCl3ic7SevQmsQZlrhflqDEdtPXv3lbyRug0KXeJSaBh1SVWqbdprxiWPAcZUaNQ00t+spKHSJR9mfWQpnWx4rEgEHsq4LXQZxEApd4lFSVWrBb6DdT6PnaPmYB6HWjsSjqCrqbXo+M+hR/e/N0BTlgNPcgTSwE9oMfxY+PWNtq6u8zqbnE9dBV7rEY6jrNNAZbLvKZZwRUv/26PjMWoS+sheBw2ag6ut1MNTadqWqMXCorNPYNAZxDRS6xGOoyusgl9r2Iy+WeyNw6HRIA4MhEonh22sQpK2Doa34yaZxvaRim6/CiWug0CUeo0FnBN/9nYwNNdDfKoM8qKtN4zDG0KA18FQVcWYUusRj6A0c+MxcZjTg5qH/Q6uIUZC1C7VtLAabpz6Ia6DQJR5DJhXztlyMMQ43j6wHJFK0HT3f5vFEItg89UFcA/0tE7eUl5eH7777DhrN/25O+cklvPSwZYyhOuVDGBtqEZT4Z4gkti8CEolE8POixUSegEKXuKVt27bhd7/7HVq3bo0nnngCK1asAHerlJeP8Le+/Rj66lJ0mPQmxDIvHqoFtAYOio7+vIxFnBv910rcCmMM169fh1gsvns0uk6H06dP4/Tp08jLy4M8ZgG0NgSv4XYl7uR8A0hkuL5ppvnxtr9bgFb9R7Z4XG+pGB0CvFv8fOI6KHSJyzIFbGZmJjIzM3Hu3DlkZmZCJBKhb9++5qkEb29vrFu3DosWLcKUbf/Fj0U1LX5NaesO6LbsCF9vwUzRKYD3MYlzotAlLoExhrKyMquABYCYmBjExsZi/vz5iImJQUhICDiOg5+fH9q2bYujR49i8ODBAIB4RTDyym471VZgb5mYTgn2IHQwJXE6jDHcuHHDIlwzMzPBcZw5YGNiYhATE4MuXbo88ObYgQMHMGzYMHTo0MH8WE2DDoPXfmfTFAPfvKRinF42Cm385EKXQhyAQpcIrrmANRgMVgEbGhrKy+qDBbuzcKyg3CnaO4pFwJgBnfDRM9FCl0IchEKXOFR5eblFwJ47dw56vd4iYGNjY3kL2ObklNZi2vbTgp2Pdi8fmQR75g1GZGig0KUQB6HQJXZTUVFhFbBardbi6jU2NhZdu3a1W8A+yAs7M/HdxUpBj+yRS8QYpeiALdNjBKuBOB6FLuGFWq22mB44d+4cmpqaLK5eY2Ji0K1bN4cHbHNqG3UY9v73qNMI1++gtY8UGa/9Bq19ZILVQByPQpc8NrVabbWKoLGx0eLqNSYmBt27d3eKgH2Q4xfUWLwnW5BpBik4bJn5K4wO7+jw1ybCotAlD1VZWWkVsHfu3LEK2B49ejh1wAKARqPBmTNncOLECezbtw+FhYV4/YvT+Kqg2qHB6y0Vo+bsIVSnbkdCQgISEhIwfPhwp/9PivCDQpeYVVVVWQVsfX09oqOjLQI2LCzM5cKhuroaISEhkMlkaGpqgtFoRLdu3XDt2jUsO5iPQ7k3HBK8PjIJEpSdMYj9hEmTJgIAfH19odPpsGDBAvztb3+zew1EWE4Ruuo6DVTldWjQGaE3cJBJxfCTSxDeKYC2RtrJzZs3rZZp3b59G9HR0RY3unr27OlyAfsgL774Inbs2AG9Xg9fX1/s2LED06ZNA2MMq1NU2HWmxK7B6yOTYHpcVySPUQAAunfvjpKSkrtf8/FBbm4uevfubbfXJ85BkNCtadBhf2YpUlVqqCrqoTNwkEvFYIyBsbtt7kQikflxRUd/xCuC8XRMKC0gb4Hq6mqrgK2pqWk2YMVi9+yBpNVqMXfuXBw6dAgajQZ+fn6orKyETPa/m1jHL6ixZH8ONHqO11UNcokYPnIx1j+ttNh5tmXLFrz66qvQ6/XgOA4bN27ESy+9xNvrEufk0NDNKa3F9pNXkapSQyTCY23F9JaJwdjdbZzzhoZBSesam3Xr1i2rZVo1NTWIioqyWEngzgF7v5s3byIpKQnt27fHli1bMHz4cDzzzDN48803rb63tlGH5QfzkVZYBa3BaNMGCrEI8JJKMKJvENYkRiDQ1/KCoaGhAV26dMGGDRtQVVWFZcuW4bnnnsOOHTta/qLE6TkkdO3yg9wnCGuSrH+QPYkpYO+dh62urjbPwZoCtlevXh4TsPe7dOkSxo4di0mTJuGvf/0rxGIx9Ho9JBLJQ/9Mcn++QDhuwwXC6J8vEB628UGn00Euv/szfOTIESQmJmLQoEFIT0+HVEqtUdyR3UPXnh/ZvGVibJis9IhmITU1NVYBe/PmTURFRVkEbO/evT02YO/3/fffY+rUqVizZg2ee+65Fo1R06DDgazrd6fCyuugMXDwamYqTGvg4C0VQ9EpAPGKYEyK7tKiqTCVSoW4uDgEBAQgLy8Pbdu2bVHdxHnZLXSFuDnhLjd8ampqkJWVZTFNUFlZaQ5Y0zRBnz59KGAf4LPPPsPy5cuxZ88ejBzZ8j6396us00BVUY8GrcF8z8HPSwpFR3/ebvrW1tZCqVTi5s2bOHXqFAYOHMjLuMQ52CV0GWOOX4YT2RlrkyJcLnhra2utAlatVkOpVFoFrEQiEbpcp8dxHP785z/jwIEDOHr0KPr27St0SS3CcRxGjRqFkydPYu/evZg4caLQJRGe2CV0Vx29gN12vsK9n+mKd8XYcIe95uO6ffu2VcBWVFQgMjLSYhVB3759KWBboLGxETNnzkRVVRUOHjyI9u3bC12SzRYsWIAtW7bgnXfewcqVK4Uuh/CA99AVcmulj0yCTdOinGKOt66uDllZWRbLtG7cuGEVsP369aOA5UF5eTkSEhKgUCiwfft2eHnxc3aZM9iyZQsWLlyIiRMnYt++fUKXQ2zEa+g6SxOR9KUjHbqqoa6uDtnZ2RYBW1ZWhoEDB1oFLN2R5l9ubi7Gjx+P559/HsnJyS43xfQo0tLS8NRTT0GhUODMmTNu9Z+Kp+E1dD2hXV59fb1VwJaWliIyMtJiFQEFrGMcPXoUs2fPxkcffYQpU6YIXY5dFRcXIzo6GhKJBDk5OejcubPQJZEW4C103bExtClg712mVVpaioEDB1oErEKhoIB1MMYYNm3ahLVr1+LgwYPmM9DcnambW1FREU6cOIEnnnhC6JLIY+ItdF39CJQ7d+5YBWxJSQkiIiIsVhGEh4dTwArMYDDgpZdeQnp6Oo4cOYLu3bsLXZJDcRyHxMREHDlyBJ999hlmzZoldEnkMfCSHjUNOqSq1E4RuADAMeC4So2aBl2zC9Tv3LmDnJwci1UExcXFGDBgAGJiYjBy5EgsXboU4eHhFnvzifDq6uowZcoUcByHH374Aa1btxa6JIcTi8X4+uuvkZycjDlz5qCgoADvv/++0GWRR8TLle62jCvYkFrodMdaLxndF89EB1sFbFFREfr3729xk6t///4UsE6uuLgY48aNw9ChQ/Hhhx/SJw4Au3btwqxZsxAfH4+UlBTaLOMCeAndyZ+cwo9FNTaNcfPw/0FTlAtOr4HErw0CBk+Ef+RTNo0pqb6K8p1vNBuwpv3uxDWcPn0aSUlJeOONN7B48WK3XKHQUufOncOwYcPQpUsXZGVloVWrVkKXRB6Cl0sFVUW9zWMEDH4a7X7/EkRSGfTVpajYvRzy4J7w6tirxWN6d+yF2tpaClgXt3fvXixcuBCfffYZxo8fL3Q5Tic2NhZXr15FVFQUunTpgszMTPTs2VPossgD2PxZRF2ngc5g+7SCPKgbRFLTx3sRRBDBUFNu05h6BtRqnGfKgzwexhhWr16N1157DcePH6fAfYiOHTuiuLgYvXr1gkKhwPHjx4UuiTyAzVe6qvI6yKViaHkI3upvN6Mh/zswgxby4J7w6Rlr03heUjFUFfV0+oQL0mq1eP7553H+/HmcPn2a1qQ+ArlcjnPnzmHGjBn43e9+h40bN2Lx4sVCl0XuY3PoNuiM4Gt/RbunXkTb0X+CtuwiNCX5EElsu7HFGEODVrjdcaRlqqurkZiYiPbt2yM9PR1+fn5Cl+RSdu7ciQEDBuDll19GQUEBtm3bJnRJ5B42Ty/oDRz47N4gEkvgHdofxvqbqM9OsWksxsDL1AdxnEuXLmHw4MF44okncODAAQrcFlq2bBm++uor/P3vf8eTTz4Jg4EuPpyFzaErk4phlxvJHGfznK5IBMillm9Rq9XaNCaxn7S0NAwbNgzLli3DunXraPmTjRISEpCbm4u8vDz06NEDt27dErokAh5C108usXn5jrGhFg0X0sHpmsA4I5quZqJBlQ7vbpE2jSsSiSAXM5w4cQJLlixB79694ePjg6qqKpvGJfz77LPPMGXKFHzxxReYO3eu0OW4jfDwcJSUlEAikaBr164oKCgQuiSPZ/OcrqJTgO0f4UUi1GcfQ/W3mwHGQdq6A9qMmgffPrbtp9fqjfjd4Aiwxtqfj1dh8Pb2Rrt27Wyrl/Dm3qbj6enp6Nevn9AluZ3AwEBcvXoVo0aNglKpxL59+5CUlCR0WR6Ll80REe98i3oB2zk+SIC3FPPaFuK1114zTytIJBJ06NDB4mTcmJgYdOzYUeBqPY+p6XhlZSX+9a9/uUXTcWf34osvYuvWrXj33XexYsUKocvxSLxMmik6+vMxDO8UnQKwaNEiZGZmIiQkBGKxGAsXLsSpU6fw7LPPQqPR4IMPPkB4eDhCQkIwYcIEvPvuuzh69CjUarXQ5bu18vJyDB8+HL6+vkhNTaXAdZDNmzfj448/xltvveX2rTCdldv3Xpg3NAzA3UYp8+fPx/z58zFs2DCL72WMoaioyKI/Q2ZmJnx9fS22D8fExCA4WPhTKVxdbm4uEhISMG/ePLdtOu7sTE3Rw8PDcebMGdq16UC8hG5Ngw6D137HywYJvnhJxTi9bFSLjsEG7gbxtWvXrIK4VatWFtMSMTEx6NChA8/Vu6+jR49izpw52LRpE11pCay4uBhRUVGQyWTIzs6mDSgOQv10HwNjDFevXrUI4qysLPj7+1sFcVBQEK+v7Q42bdqENWvW4Msvv6Tm206CmqI7Hp0cYSPGGK5cuWLR/DwrKwutW7e2CmJPnbc0GAx4+eWX8f333+PIkSPo0aOH0CWRe3AchwkTJiAlJYWaojsAnZFmBxzHNRvEbdq0sTjmJyYmxu2Xr93bdHzfvn0e2XTcVSxfvhzr1q3DkiVLqCm6HbnlacAZr/0GrX2cqyE5x3H46aefLII4Ozsbbdu2tQritm3bCl0uL0xNx5988kls2rSJmo67AFNT9NGjR+Po0aO0K9AOeA1dADh+QY3Fe7IFmWbwkUmwaVoU4hWuscLAFMSm+WFTELdr185iHXF0dLTLBfGZM2eQmJiI119/HS+99BKtUHAh1BTdvngPXQD4y9EL2HWmxKHB6yOTYHpcV6wYG+6w17QHjuNw+fJlixUT2dnZaN++vVUQt2nTRuhym7Vv3z4sWLCAmo67sIqKCkRFRUGj0SAzMxNhYWFCl+Q27BK6jDEsO5iPQ7k3HBK8PjIJEpSdsTYxwi2vqDiOQ2FhoVUQ37uzLiYmRvAgZoxhzZo12Lp1Kw4dOgSlUilYLcR2Op0OQ4YMQV5eHlJSUhAfHy90SW7BLqEL/Nz1P0Vl9yte0xVu8hiFWwbugxiNRhQWFlosX8vJyUFwcLDF/HB0dDQCAwPtXo9Wq8Wf/vQnFBQU4NChQ7Tm041Mnz4de/bswd/+9jcsWrRI6HJcnt1C1+T4BTWW7M+BRs/xuqpBLhHDRy7G+qeVLjOHa29GoxGXLl2yCOLc3Fx07NjRKoj5XEVQXV2NpKQktG3bFjt37qQeuG5o7dq1+POf/4w//vGP1BTdRnYPXeDuqoblB/ORVlgFrcFo0wYKsQjwkkowom8Q1iRGINCXti8+jNFoxMWLFy1WTeTm5qJz585WQRwQEPDY4xcWFmLs2LFISkrCmjVr6G63Gzt06BAmTpyIuLg4pKWl0WqUFnJI6JrkltZi+8mrOK5SQyTCY/Vq8JaJwRgwWhGMeUPDHLLxwV0ZDAarIM7Ly0NISIhFEEdFRT00iNPS0jBlyhSsXr0af/zjHx34DohQzp8/j8GDB6NNmzbIyclxuVU1zsChoWtS06DDgazrSFWpoSqvg8bAwUsq/rnn7d0TH0QiEbQGDt5SMRSdAhCvCMak6C4t7qVAHs5gMEClUlkFcWhoqMU64qioKPj7++Pvf/873njjDXzxxRcYNWqU0OUTB6qtrYVSqcTNmzdx+vRpDBgwQOiSXIogoXu/yjoNVBX1aNAaoDNwkEvF8POSQtHRn07yFZApiO9dR5yXlwcfHx/odDosWLAAY8eORVRUFK3l9DAcx+E3v/kN/vOf/2D//v1ITEwUuiSX4RShS1xDY2MjZsyYgaKiIsyaNcu8eiI/Px/dunWzuCJWKpUUxB7A1BR91apVSE5OFrocl0Az4eSRlJeXY8KECejbty/++9//wsvLy/w1vV6PCxcumK+Iv/jiCxQUFKBbt24W64ijoqJoZYOb2bx5MwYMGIBFixYhPz8fe/bsEbokp0dXuuQX5eXlYfz48fjjH/+IFStWPNJ6aL1ej/Pnz1ssXysoKECPHj0sglipVFIQu4ETJ07g97//PTVFfwQUuuShUlJSMHv2bHz44YeYOnWqTWPpdDqrID5//jzCwsIsVk0olUr4+vry9A6Io1y7dg0xMTHUFP0XUOiSB9q0aRP++te/4uDBg3Zrbq3T6VBQUGCxauLChQvo2bOnRRBHRkZSELuAxsZGREdHo7i4GGlpaYiLixO6JKdDoUusCN10XKvVWgWxSqVCr169rILYx8fHobWRX0ZN0R+OQpdYMDUdNxqN2L9/v9M0HddqtcjPz7cI4osXL6J3794WQTxw4EAKYidhaoq+dOlSvPfee0KX4zQodInZvU3HP/zwQ8hkztUI/n6mIL63+9rFixfRp08fi+VrAwcOhLc3rfcWgqkp+m9/+1scOXKEtomDQpf8zF2ajms0GqsgvnTpEvr27WsRxBERERTEDnL27FkMHz6cmqL/jEKXYP/+/XjxxRfdtum4RqNBXl6eRRAXFhaiX79+VkF87/pjwp+KigoolUpotVpkZWV59OGkFLoejDGGtWvXYsuWLR7XdLypqQl5eXkWy9cuX76Mfv36WZzQMWDAAApinuh0OjzxxBPIz8/HsWPHPLZnB4Wuh9LpdHj++eeRn5+Pw4cP05pK3A3i3NxciyD+6aefoFAoLDZ0RERE0OJ/GzzzzDPYu3cvPvjgAyxcuFDochyOQtcDUdPxR9fY2GgVxFeuXEF4eLhFEA8YMICC+DGsWbMGycnJmDdvHj755BOhy3EoCl0PY2o6npiYiLVr19Ld5BZobGxETk6OxfK1a9euNRvEzr4CREie2hSdQteDUNNx+2loaLAK4qKiIvTv398iiPv3709BfA97NUVX12mgKq9Dg84IvYGDTCqGn1yC8E4BgreLpdD1EH//+9+xbNky7N6922NvYDiaKYjvXTVRVFSEAQMGWGzoCA8P9+ggrqmpgVKpxK1bt/Df//63RU3Raxp02J9ZevdghIp6c1/u+w9GMD2u6OiPeEUwno4JdfjBCBS6bo7jOKxYsQL79u3DkSNH0K9fP6FL8mh37tyxCuLi4mJERERYBbGnfNwG7v6cjhw5EqdOncK+ffswcuRIzJ07F5988gnat2//wOfl/HwEWKoNR4DF/3wEmNJBR4BR6LqxxsZGzJo1C+Xl5fjqq68e+sNLhFNfX28VxCUlJRg4cKBFECsUCrcP4hdeeAFbt25FaGgobty4gbfffhsrVqyw+j67HHbbJwhrkux/2C2FrpuqqKhAQkIC+vTpg08//ZTWmrqY+vp6ZGdnW6yauH79ujmITWHcr18/twpixhji4uJw9uxZAEC7du1QUVFh8R6PX1Bjyf4caPQcdMZHv7L9JXKJGN4yMTZMViJeEczbuPej0HVDeXl5SEhIwHPPPYeVK1e67JZeYqmurs4qiMvKyhAZGWkVxBKJROhyW2Tv3r0WfZvFYjF2796NKVOmgDGG1Skq7DpTgia90W41+MgkmB7XFcljFHb5t0Oh62ZSUlIwa9YsfPjhh5g2bZrQ5RA7u337tlUQl5eXWwVx3759XSKIb9++jS+//BIpKSk4ceIEampq4OXlhcrKSqw+XoRDuTfsGrgmPjIJEiI7Y21SBO/BS6HrRkxNx7/88ksMGTJE6HKIQG7fvo2srCyL5WsVFRWIjIy02OLcp08fpw/iixcv4qWXXgJTJqLUq7tDAtfEdMW7Ymw4r+NS6LoBg8GAV155BSdOnBCk6ThxfrW1tcjOzjZfDZ87dw5qtRpKpdIqiJ1tw8zxC2os3pPt0MA18ZFJsGlaFK9zvBS6Lq6urg5Tp06FwWBwqqbjxPnV1tYiKyvLIoirqqosgjgmJkbQIK5t1GHY+9+jTmMQ5PUBoLWPFOlLR/K2qoFC14WVlJRg3LhxGDJkCDZt2uTRC+wJP2pqaiyCODMzE1VVVYiKirII4t69ezskiF/YmYnvLlbyukrhccklYoxSdMCW6TG8jEeh66J+/PFHJCYmYunSpXj55ZdphQKxm1u3bpnniE1hXF1dbRXEvXr14jWIc0prMW37aUGmFe7nI5Ngz7zBiORhAwWFrgs6cOAAXnjhBXz66adISEgQuhzigaqrq62CuKamBlFRURYbOnr27NniIF6wOwvHCspt2vjAF7EIGDOgEz56JtrmsSh0XYip6fjmzZtx6NAhREVFCV0SIWbV1dUWKyYyMzNRW1uL6OhoqyD+pU9mNQ06DF77HbQG4aYV7uclFeP0slE292qg0HURpqbjeXl5OHz4MEJCQoQuiZBfdPPmTXMQm8K4rq7OKojDwsIsgnhbxhVsSC18rF4K9uYtE2PJ6L6YNzTMpnEodF1AdXU1Jk6ciMDAQOzatYuajhOXVlVVZRXE9fX15rnhmJgY7L3ZGZmlt1v8GnWZh9GQ/x10VUXwUwxH+3Gv8FJ7XI+22Pv8EzaNQaHr5C5fvoyxY8diwoQJWLt2rdMvZiekJSorK82rJrKysqDqPxdNNqwSa7x0ChCJ0HQtC0yv4y10A7ylyHvrKZvGoNB1Yunp6ZgyZQpWrVqFefPmCV0OIQ6hrtNg2Pvf8zKfW5Px/2Csu8lb6MqlYvzntZE2NUJ3rq0nxOwf//gHJk+ejF27dlHgEo+iKq+DXOqc0eQlFUNVUW/TGO7TE85NcByHlStXYs+ePUhLS4NCoRC6JEIcqkFnhLN+AGeMoUFr2+44Cl0n0tTUhGeffRbl5eU4ffo0goKChC6JEIfTGzg4aeaCMUBn47SHc17De6CKigqMGDECcrkcqampFLjEY8mkYjjrBkuRCDZPfVDoOoH8/HwMHjwYY8eOxc6dO+HtLexppYQIyU8usXlbO+OMYAYdwBkBxoEZdGCc7duJRSIR/LxsmyCg6QWBHTt2DM8++yw++OADPPPMM0KXQ4ggGGO4du0aMjMz8Z+sAmjEg2wa7/YPe3D7hy/Mv284/z1a/3oaAodOt2lcrYGDoqO/TWNQ6Aroo48+wurVq/H1119T03HiMRhjKCoqstgckZWVBR8fH/PmCC+jGHpdy+dOA4dOtzlgm+MtFdu0XAyg0BWEwWDAq6++itTUVJw6dYqajhO3xRhDcXGxVU8Gb29v8zbgl19+GTExMejYsaP5eQWfnMKPRTUCVt48RacAm8eg0HWw+vp6TJ06FXq9HqdOnUJgYKDQJRHCC8YYSkpKLMI1MzMTcrncHLCLFy9GTEwMOnXq9NCx4hXByCu77XS9F/g4QYJ2pDkQNR0n7oIxhtLSUquAlUqlFo1sYmJi0Llz58ce3527jNGVroOcPXsWf/jDH7BkyRK88sor1HScuAzGGK5fv24RsOfOnYNEIjEH7IsvvojY2NgWBWxz2vjJEa8Idqp+uqMVwTYHLkCh6xDUdJy4CsYYysrKrAJWJBKZr17nz59vDlh7XjzMGxqGExcrneLkCC+pxOaWjiY0vWBH9zYd//rrrxEdbXvXeUL4whjDjRs3rAKWMWYxPRAbG4uQkBBBPp3RGWnkkel0OvzpT39Cbm4uNR0nTqG5gDUajRZHsMfExKBLly5OM/3lLKcBZ7z2G7T24eceDE0v2MGtW7eQlJSEwMBAZGRkoFWrVkKXRDxMeXm5xQ2uc+fOQa/Xm4N17ty52Lx5M0JDQ50mYJvjKwViDSpkIAwGOL6XtI9MgvVPK3kLXICudHlHTceJo1VUVFgFrE6ns5geiImJQdeuXZ06YIG7a9h//PFHfP/99zh06BDOnj0LsViMV3f+B4cv1Dp0ftdHJsH0uK5YMTac13EpdHlkajr+7rvv4vnnnxe6HOKGKioqrDYaaDQaq4Dt1q2b0wdscz7//HPMnj0bMpkMer0eALBu3Tq89tprWHYwH4dybzgkeH1kEiQoO2NtYgTvf44Uujz5/PPP8dprr2H37t2Ij48XuhziBtRqtVXANjY2Wq2D7d69u0sGbHP0ej2io6NRUFAAAPDz80NlZSV8fX3BGMPqFBV2nSmxa/CarnCTxyjs8udKc7o24jgOb775Jr744gukp6dT03HSIpWVlVYBe+fOHXPAzpgxAxs3bkSPHj3cJmCb8+mnn+L8+fPw9fWFVqvF3Llz4evrC+Buh68VY8MR16MdluzPgUbP8bqqQS4Rw0cuxvqnlbzsPHsQutK1QVNTE2bNmoWysjJ89dVX1AOXPJJ7T8M1BWx9fT2io6MtrmDvP5bc3S1YsABbtmzBO++8g5kzZ+L3v/89jh07hu7du1t9b22jDssP5iOtsApag9GmDRRi0d11uCP6BmFNYgQCfW3fAPEwFLotVFFRgQkTJqBXr1749NNPqQcuadbNmzetAvb27dtWAduzZ0+PCth7cRyHUaNG4eTJk9i3bx+SkpIe+bm5pbXYfvIqjqvUEInwWL0avGViMHZ3p9m8oWGIDA1sQfWPj0K3BfLz8zF+/HjMmTMHb775psf+YyGWqqurrXoR1NTUWARsbGwswsLCIBbT+QEAUFtbC6VSiZs3b+LUqVMYOHBgi8apadDhQNZ1pKrUUJXXQWPg4CUVgzEGxu6e+CASiaA1cPCWiqHoFIB4RTAmRXfhZWvv46DQfUzffPMNnn32Wfztb3+jpuMe7NatW1YBe+vWLURHR1usJOjZsycF7AOoVCrExcWhdevWyM3NRdu2bXkbu7JOA1VFPRq0BugMHORSMfy8pFB09Le5H66tKHQfw8cff4y//OUvOHDgAH79618LXQ5xkFu3biErK8tiJ1d1dbVVwPbq1YsC9hEdOXIEiYmJiIuLQ1paGqRSz7mnT6H7CIxGI1599VUcP34cR44cQVgYP40viPOpqamxCtibN28iKirKImB79+5NAdtC7733HpYtW4a5c+di+/btQpfjcBS6v6C+vh7Tpk2DVqvF/v37qem4GzEF7L03uiorK80Ba5qH7dOnDwUsT2bOnIndu3djw4YNeOmll4QuRxAUug9RUlKC8ePHY/Dgwfjoo4+o6bgLq62ttQpYtVoNpVJpFbC0dZt/Op0OQ4YMQV5eHo4ePYrRo0cLXZJgKHQfgJqOu67bt29bBWxFRQUiIyMtArZv374UsA5QWVmJyMhINDU1ITMzEz179hS6JEFR6Dbjyy+/xPz587Fjxw5MmDBB6HLIQ9TV1ZkD1jQPe+PGDauA7devHwWsAM6dO4dhw4YhJCQE2dnZ1HEPFLoWGGNYt24dPv74Y2o67oTq6uqQnZ1tEbBlZWUYOHCgxUaDfv36edTdcGf1xRdfYObMmYiPj0dKSgrNi/+MQvdnOp0O8+fPR05ODjUddwL19fVWAVtaWmoVsAqFggLWCSUnJ2PNmjV45ZVXsH79eqHLcSoUuri7DnPixIkICAjArl276COQg925cwfZ2dkWGw1KSkoQERFhEbDh4eEUsE6O4zgkJibiyJEj2LFjB+bMmSN0SU7H40P38uXLGDduHMaPH49169bRvJ+d3blzBzk5ORYBW1xcjAEDBlhslVUoFLRaxMU0NjYiNjYW165dw4kTJ/DEE08IXZJT8ujQzcjIwOTJk6npuJ00NDRYBWxRUREGDBhgsdEgPDycAtbFFRcXIzo6GhKJBFlZWejSpYvQJTktjw3df/7zn1i6dCl27drl0WsG+WIK2HuXaV27ds0qYPv3708B62bS09Px29/+Fv369cOZM2eo494v8LjQNTUd3717N44cOYLwcH7PP/IEjY2NVgF79epV9O/f32KZVv/+/SGXO7aDE3GsTz75BC+++CKSkpKwf/9+octxCR4Vuqam49evX8dXX32FDh06CF2S02tsbERubq7FKoIrV64gPDzcImAHDBhAAethFi1ahI8//hhvvfUW3nrrLaHLcRkeE7pqtRoTJkxAWFgYPvvsM/oI1IympiargP3pp5+gUCisAtbLy0vocolA7m06vnfvXkycOFHoklyKU4Suuk4DVXkdGnRG6A0cZFIx/OQShHcK4KX3ZUFBAcaNG0dNx+/R1NSEvLw8i4C9fPky+vXrZxGwERERFLDErLa2FlFRUaiqqrKp6bgnEyR0axp02J9ZerfLe0W9ucnw/V3eTY8rOvojXhGMp2NCH7vLu6np+MaNGzF9+nQ7vSPnptFokJeXZ7GKoLCwEH379rVYBxsREUGfAMgD3dt0PCcnB+3atRO6JJfk0NDN+fk8o1QbzjOK//k8I+UjnGe0efNmrFq1yqOajms0GuTn51sE7KVLl9CnTx+LdbAUsORxpKSkYMKECRg0aBDS09Npk4oNHBK6djm5s08Q1iQ1f3Kn0WjEkiVL8O233+Lo0aNu23Rcq9VaBezFixfRp08fi2VaAwcOpIAlLWZqOv7cc89hx44dQpfj8uweuscvqO12Rr23TIwNky3PqDc1HddoNDhw4IDbNB03Bey9y7QuXryI3r17WwWsj4+P0OUSN/Hss89i586d2Lhxo8c2Heeb3UKXMYbVKSrsOlOCJr3RHi8BAPCRSTA9riuSxyhQWlqK8ePHIy4uDh9//LHLLsLX6XRWAatSqdCrVy+LgI2MjKSAJXah1+vx61//2twA6qmnnhK6JLdhl9BljGHZwXwcyr1h18A18ZFJMKybD468+QxeffVVvPrqqy6zQkGn06GgoMAiYC9cuICePXtarCKIjIyEr6+v0OUSD3Bv0/Fz586hV69eQpfkVuwSuquOXsBuO1/h3s9LIsLQTsCOBWMc9pqPS6/XmwPWNA97/vx5hIWFWQWsn5+f0OUSD5SZmYmhQ4eic+fOyMnJoY57dsB76B6/oMbiPdkODVwTH5kEm6ZFWczxCkWv1+P8+fNWAdu9e3eLgFUqlRSwxCns2bMHM2bMoKbjdsZr6NY26jDs/e9RpzHwNeRja+0jRfrSkc2uarAXvV6PCxcuWKwiKCgoQLdu3SzWwSqVSrpyIE5p5cqVWL16NV5++WVs2LBB6HLcGq+h+8LOTHx3sZLXVQqPSy4RY5SiA7ZMj7HL+AaDwSpg8/Pz0bVrV4uAjYqKooAlTu/epuPbt2/Hc889J3RJbo+3Fc45pbVIK6wSNHABQGfkkHapCrmltYh8hA0UD2MwGKBSqSwCNi8vD6GhoeaAnTp1KpRKJfz9/fl5A4Q4SGNjI371q1/hypUryMjI8JgNRELjLXS3n7wKrcHx87jN0RqM2H7yKj565n8HSxqNRrz33nsYN24cIiIirJ5jCth7VxGYAtZ09Tp58mRERUVRwBKXV1xcjJiYGIjFYvz000/UdNyBeAndmgYdUlVqm3aa8YljwHGVGjUNOrTxk6OqqgoTJkzAmTNnUF9fj/Dw8GYDNiQkxBywkyZNQlRUFAICAoR+O4TwKiMjA6NHj0bfvn3x448/0m5FB+NlTndbxhVsSC18rF4K9uYtE2PJ6L4YIFVj7NixqKurg9FohL+/PziOQ+fOnS1WEURFRaF169ZCl02IXW3btg0vvPACNR0XEC9XuqkqNW+Bq79VhhufLoRfv1+j/filLR5Ho+dwNKcEzy8eafG40WjE9evX0aZNG1tLJcSlUNNx58BL6Koq6vkYBgBw699b4dWpNy9jXb2lxYULF3DixAkcPHgQp06dQmNjIwwG4Za0EeJoHMchPj4eGRkZ2LdvHyZNmiR0SR7N5tBV12mgM/BzldtwIR1ibz/I2vWDobbc5vE0Bg7tQnpgwYIFWLBgAfR6PS5evIj27dvzUC0hzq+urg5KpRKVlZXIzMxEZGSk0CV5PJu3nKjK6yCX2r5zhdM2ovbkLrT5zVybxzLxkootrsJlMhkiIiJcpi8DIba4ePEiQkNDodPpUFRURIHrJGxOywadEXzsr6jN+H9oFflbSAOCbB7LhDGGBi1NJRDPk5KSgoiICPTv3x9FRUX06c6J2By6egMHWzNXp74KTXEuAn41wdZyLDAG3qY+CHEV77//PsaNG4eZM2fi1KlTdMqDk7H5b0MmFcPWT+uaknwYbqtxffMcAADTaQDGofzmS+g054MWjysSgZepD0Jchanp+IYNG/Dyyy8LXQ5phs2h6yeX2DxH2kr5FPwUw8y/r/vxIAy31Wj71AKbxhWJRPDzov/lifszGAwYMmQIcnJycOzYMWo67sRsTiRFpwCbP8KLZd6A7H+7YkQyb4ikckh8bdusoDVwUHSkLbvEvVVWVkKpVKKxsRHnz59H7978LLkk9mFz6AYHeEMuFUPL49xp4FB+jkr3lorRIYC2OBL3lZWVhSeffBKdO3fGpUuXqC+IC+BlwtNZryYVnahvAnFfe/fuxaBBgzB06FAUFhZS4LoIXkI3XhEMb5lz3bDylomd4gQJQuxh5cqVmDZtGhYtWoRvv/2WTnlwIbz8TT0dE2rzsjG+MQZMiqZ2dcS9cByHP/zhD/jrX/+K7du3Y+PGjUKXRB4TL7f22/jJEa8IxrGCcqdo7ygWAaMVwWjj57gjewixt6amJsTGxlLTcRfH22eSeUPD4CWV8DWcTbykEswbGiZ0GYTwpqSkBKGhoaisrMRPP/1EgevCeAtdZWggRvQJglwi7NySXCLGiL5BNh/VQ4izyMjIQJ8+fdC5c2eUlpbSKQ8ujteEXJMUIfgNNR+5GGuTBgpaAyF82bZtG0aOHIlx48YhLy+PTnlwA7wmZKCvHOufVsJHJsw0g49MgvVPK9HaRybI6xPCp8WLF2P+/PlYuXIlDhw4IHQ5hCe875EdHR6M6XFdsetMCZr0jjuo0kcmwfS4rrRMjLg8juMwevRopKenU9NxN2SXxgTJYxSo1xhwKPeGQ4LXRyZBgrIzksco7P5axLOp6zRQldehQWeE3sBBJhXDTy5BeKcAXnY/3tt0/Ny5c1AqlbYXTZyKXUJXJBJhbVIE/L2ldr/iNV3hJo9RUHNywruaBh32Z5YiVaWGqqIeOgMHuVQMxhgYu9vJTiQSmR9XdPRHvCIYT8eEPvaSxUuXLmHQoEHw9/enHrhujJfTgB/m+AU1luzPgUbPQWfkrz+DXCKGj1yM9U8raUqB8C6ntBbbT15FqkoNkQiPdfCqt0wMxu7u1Jw3NAzKB6yk0ev1iIuLw/r169HU1IQJEyYgNjYWJ0+epB64bszuoQsAtY06LD+Yj7TCKmgNRps2UIhFd9fhjugbhDWJEQj0pQ0QhD92+VntE4Q1SdY/q3v37sWsWbPAGINer8fs2bPx2Wef2fgOiLNzSOia5P589XDchquH0T9fPdA6XMI3e34q85aJsWGy5acypVKJ3NxcAEBAQACuXLlCUwoewKGha1LToMOBrOt358nK66AxcPBqZp5Ma+DgLRVD0SkA8YpgTIruQlt7Ce8YY1idonLo/YfMzEwMGjQIjDGIRCJIJBLMmzcPmzdvttvrE+cgSOjer7JOA1VFPRq0BvMNCT8vKRQd/akfLrErxhiWHcx36Eqb3/cPwj8XPIXq6mr06tULCQkJiI+Px9ChQ9GqVSu710CE5RShS4hQVh29gN0OXlMuFzN0qC/EwRUz0KFDB4e9LnEO1ISTeKzjF9QOD1wA0HEiVAeGI6+arnc8EYUu8Ui1jTos2Z/j8MA1adIbsWR/DmobdYK8PhEOhS7xSMsP5j/W6hl7aNJxWP6vfEFrII5HoUs8Tk5pLdIKq3hdFtYSOiOHtEtVyC2tFbQO4lgUusTjbD95FVqDMNMK99MajNh+8qrQZRAHotAlHqWmQYdUldopjpUCAI4Bx1Vq1DTQ3K6noNAlHmV/ZimcrS+SSAQcyLoudBnEQairBvEoqSq1zTfQKnYtg/bGJYjEd5v1S/zbIeT5T1o8nkbPIVWlpnP9PASFLvEoqop6XsZp+9v58I98ipexAEBVXsfbWMS50fQC8RjqOg10BmFXLDyIxsChsk4jdBnEAehKl3gMVXkd5FIxtDwEb23a56hN+xyytiEIHDYT3t1sOwzVSyqGqqKeeo14AApd4jEadEbw0Wqkzcg5kLULhUgiQ4MqA5VfrkKnOR9C1qZTi8dkjKFBa7C5NuL8aHqBeAy9gQMf7Z28OveF2MsXIqkMrSJGwStEgaYr52wakzE47dQH4ReFLvEYMqnYPsvFRCIAtqW5SATIpfTP0RPQ3zLxGH5yic2Hl3KaO2i6mglm0IFxRtw5/z20pQXw6RFt07gikQh+XjTb5wnob5l4DEWnAJs/wjPOiNqMndDfug6IxJC164KgpBWQteti07haAwdFR3+bxiCugUKXeIzgAG+bVy9IfFuj0+yNPFZ1l7dUTCsXPARNLxCP4qxXk4pOAUKXQByEQpd4lHhFMLxlzvVj7y0TW5wSTNybc/30EWJnT8eE8rJsjE+MAZOibZsTJq6DQpd4lDZ+csQrgiF2kk5jYhEwWhGMNn5yoUshDkKhSzzOvKFh8JJKhC4DAOAllVB3MQ9DoUs8jjI0ECP6BEEuEfbHXy4RY0TfIESGBgpaB3EsCl3ikdYkRQh+Q81HLsbaJNsa5RDXQ6FLPFKgrxzrn1bCRybMNIOPTIL1TyvR2kcmyOsT4VDoEo81OjwY0+O6Ojx4fWQSTI/rSsvEPBSFLvFoyWMUSIjs7LDg9ZFJkKDsjOQxCoe8HnE+IsZHg1FCXBhjDKtTVNh1pgRNevsdzW66wk0eo7C58Q5xXRS6hPzs+AU1luzPgUbPQWfkr7etXCKGj1yM9U8raUqBUOgScq/aRh2WH8xHWmEVtAYjOBv+dYhFd9fhjugbhDWJEQj0pQ0QhEKXkGblltZi+8mrOK5SQyTCYx3b7i0Tg7G7O83mDQ2jdbjEAoUuIQ9R06DDgazrSFWpoSqvg8bAwUsqBmMMjN098UEkEkFr4OAtFUPRKQDximBMiu5CW3tJsyh0CSHEgWjJGCGEOBCFLiGEOBCFLiGEOBCFLiGEOBCFLiGEOND/B2ypyi1ow/pLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "network = [\"c1\",\"c2\"]\n",
    "Alldata =  pd.DataFrame()\n",
    "for i in range(1,3):\n",
    "    data = pd.read_csv(\"..\\\\ImportedData\\\\network_c\"+str(i)+\".csv\", sep=\",\", header=0)\n",
    "    #data.columns = [ \"from\", \"to\" ]\n",
    "    data[\"network\"] = i\n",
    "    Alldata = Alldata.append(data)\n",
    "\n",
    "    \n",
    "for i in range(1,3):\n",
    "    # Build your graph. Note that we use the DiGraph function to create the graph!\n",
    "    G=nx.from_pandas_edgelist(Alldata[Alldata.network==i], 'from', 'to', create_using=nx.DiGraph() )\n",
    "\n",
    "    # Make the graph\n",
    "    f1 = plt.figure()\n",
    "    plt.title(\"Directed Network_\"+str(i))\n",
    "    nx.draw(G, pos=nx.circular_layout(G), with_labels=True, node_size=1500, alpha=1, arrows=True) \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71685cb3",
   "metadata": {},
   "source": [
    "# Content\n",
    "\n",
    "- [Prepare data](#Prepare-data)\n",
    "\n",
    "- [Define rules](#Define-rules)\n",
    "- [Set up the loglikelihood function](#Set-up-the-loglikelihood-function)\n",
    "- [Run likelihood maximization](#Run-likelihood-maximization)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "788ad15b",
   "metadata": {},
   "source": [
    "## Background \n",
    " -- *From Romero and Rosokha ([2018](https://www.sciencedirect.com/science/article/pii/S0014292118300394),[2019a](https://www.aeaweb.org/articles?id=10.1257/mic.20160220),[2019b](https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3290732)) and Rosokha and Wei ([2020](https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3526505))*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "339bf944",
   "metadata": {},
   "source": [
    "Strategy frequency estimation method (SFEM) is a finite-mixture estimation aproach to estimate the proportion of strategies in experimental data proposed by Dal Bo and Frechette (2011). The method works by first specifying the set of $K$ strategies considered by the modeler. Then, for each subject $n\\in \\{1,...,N\\}$, and each strategy $k\\in \\{1,...,K\\}$, the method prescribes to compare subject $n$'s actual play with how strategy $k$ would have played in her place. Let $C(k,n)$ denote the number of periods in which subject $n$'s play correctly matches the play of strategy $k$. Then, let $C$ denote a $K \\times N$ matrix of the number of correct matches for all combinations of subjects and strategies. Similarly, let $E$ denote a $K \\times N$ matrix of the number of mismatches when comparing subjects' play with what the strategies would do in their place. Then, define a Hadamard-product $P$:\n",
    "\n",
    "\\begin{equation}\n",
    "    P = \\beta^{C}\\circ(1-\\beta)^{E},\n",
    "\\end{equation}\n",
    "\n",
    "where $\\beta$ is the probability that a subject plays according to a strategy and $(1-\\beta)$ the probability that the subjects deviates from that strategy. Thus,\n",
    "each entry $P(k,n)$ is the likelihood that the observed choices by subject $n$ were generated by strategy $k$. Then, using the matrix dot product, we define the log-likelihood function $\\mathcal{L}$:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\mathcal{L}(\\beta,\\phi)= ln \\big( \\phi' \\cdot P \\big) \\cdot  \\mathbf{1} .\n",
    "\\end{equation}\n",
    "\n",
    "For this example, the set of strategies is the 20 strategies commonly used in the literature on repeated Prisoner's Dilemma (Fudenberg, Rand, and Dreber 2012). For details, see Appendix H of Romero and Rosokha (2018)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a33d67c",
   "metadata": {},
   "source": [
    "# Prepare data\n",
    "\n",
    "- [Return](#Content)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "id": "a7ed30e6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['subject_id', 'subject_idnum', 'pcode', 'session', 'network_type',\n",
       "       'waiting_cost', 'round_id', 'pid_in_session', 'ppayoff', 'pid_ingroup',\n",
       "       'p_payoff', 'p_secspent', 'p_currentball', 'p_node', 'p_period',\n",
       "       'p_lastinforset', 'p_currentinforset', 'p_history', 'p_time_on_page',\n",
       "       'p_guessedurn', 'g_id', 'g_period', 'g_trueurn', 'g_round', 'g_match',\n",
       "       'prnd1playerid_in_group', 'prnd1playerpayoff', 'prnd1playerpayoff_list',\n",
       "       'prnd1groupid_in_subsession', 'prnd1subsessionround_number',\n",
       "       'id_in_group', 'payoff', 'date', 'app', 'prnd1playerpayoff_total',\n",
       "       'pay_sig1', 'pay_sig2', 'pay_sig3', 'pay_sig4', 'bid_a', 'bid_b',\n",
       "       'bid_c', 'bid_d', 'bid_e', 'treatment', 'treatment_degree1',\n",
       "       'prediction', 'prediction_acc', 'prediction_welfare', 'in_degree',\n",
       "       'out_degree', 'in_degree_m', 'max_turn', 'wait2long', 'min_turn',\n",
       "       'wait2short', 'rational', 'Know_network', 'Know_prob', 'in_degree_1',\n",
       "       'Settle1', 'Settle_sum', 'within_group_round', 'Settle_latter2',\n",
       "       'BR_long', 'A_period', 'B_period', 'C_period', 'D_period', 'E_period',\n",
       "       'direct_neighbor_period', 'level_2_neighbor', 'level_3_neighbor',\n",
       "       'level_1_neighbor_wait', 'level_1_neighbor_wait2',\n",
       "       'level_2_neighbor_wait', 'inter_level2s', 'optimal_turn_empirical'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 157,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data  = pd.read_csv(datafolder2+\"DataAnalysis_v4_EmpiricalOptimalDecisionTurn.csv\")\n",
    "data.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a530cf0a",
   "metadata": {},
   "source": [
    "# Define rules\n",
    "\n",
    "- [Return](#Content)\n",
    "\n",
    "\n",
    "## optimal_turn_empirical \n",
    "\n",
    "Defined in appendix\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "0f94ea87",
   "metadata": {},
   "outputs": [],
   "source": [
    "def EmpiricalOptimal(row):\n",
    "    fit = 0\n",
    "    if row['network_type']=='a':\n",
    "        if row['p_node'] == 'A' or row['p_node']=='E':\n",
    "            fit = row['p_period']==0\n",
    "        elif row['p_node'] == 'D':\n",
    "            fit = row['p_period']==0\n",
    "            if row['A_period']==0 and row['waiting_cost']=='L':\n",
    "                fit = row['p_period']==1\n",
    "            elif row['A_period']>=1 and row['waiting_cost']=='L':\n",
    "                fit = row['p_period']==2\n",
    "        elif row['p_node'] == 'C':\n",
    "            if (row['A_period']==0 or row['B_period']==0) and row['waiting_cost']=='L':\n",
    "                fit = row['p_period']==1\n",
    "            elif row['waiting_cost']=='H':\n",
    "                fit = row['p_period']==1\n",
    "            elif (row['A_period']>=1 and row['B_period']>=1) and row['waiting_cost']=='L':\n",
    "                fit = row['p_period']==2\n",
    "        elif row['p_node'] == 'B':\n",
    "            if (row['A_period']==0 and row['waiting_cost']=='L') or row['waiting_cost']=='H' :\n",
    "                fit = row['p_period']==1\n",
    "            elif row['A_period']>=1  and row['waiting_cost']=='L':\n",
    "                fit = row['p_period']==2\n",
    "                \n",
    "    if row['network_type']=='b':\n",
    "        if row['p_node'] == 'C' or  row['p_node'] == 'D' or row['p_node']=='E':\n",
    "            if row['waiting_cost']=='H' or ( row['waiting_cost']=='L' and row['A_period']==0):\n",
    "                fit = row['p_period']==0\n",
    "            elif  row['waiting_cost']=='L' and row['A_period']==0:\n",
    "                fit = row['p_period']==1\n",
    "            elif row['waiting_cost']=='L' and (row['A_period']>=1 or row['B_period']==0):\n",
    "                fit = row['p_period']==2\n",
    "            elif row['waiting_cost']=='L' and (row['A_period']>=2 and row['B_period']>=1):\n",
    "                fit = row['p_period']==3\n",
    "        elif row['p_node'] == 'A' :\n",
    "            if row['waiting_cost']=='H' or ( row['waiting_cost']=='L' and row['B_period']==0):\n",
    "                fit = row['p_period']==1\n",
    "            elif row['waiting_cost']=='L' and row['B_period']>=1 :\n",
    "                fit = row['p_period']==2\n",
    "        elif row['p_node'] == 'B' :\n",
    "            if row['waiting_cost']=='H' or ( row['waiting_cost']=='L' and row['A_period']==0):\n",
    "                fit = row['p_period']==1\n",
    "            elif row['waiting_cost']=='L' and row['A_period']>=1 :\n",
    "                fit = row['p_period']==2    \n",
    "    return fit "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "id": "e82219c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def MoreBetter(row):\n",
    "    fit = 0\n",
    "    if row['network_type']=='a':\n",
    "        if  row['p_node'] == 'D':\n",
    "            if row['A_period']==0:\n",
    "                fit = row['p_period']==1\n",
    "            elif row['A_period']>=1 :\n",
    "                fit = row['p_period']==2\n",
    "        elif row['p_node'] == 'A' or row['p_node']=='E' :\n",
    "            fit = row['p_period']==1\n",
    "        elif row['p_node'] == 'C':\n",
    "            if row['A_period']==0 and row['B_period']==0:\n",
    "                fit = row['p_period']==1\n",
    "            else:\n",
    "                fit = row['p_period']==2\n",
    "        elif row['p_node'] == 'B':\n",
    "            if row['A_period']==0  :\n",
    "                fit = row['p_period']==1\n",
    "            else:\n",
    "                fit = row['p_period']==2\n",
    "                \n",
    "    if row['network_type']=='b':\n",
    "        if row['p_node'] == 'C' or  row['p_node'] == 'D' or row['p_node']=='E':\n",
    "            if row['A_period']==0:\n",
    "                fit = row['p_period']==1\n",
    "            elif row['A_period']>=1 and row['B_period']==0:\n",
    "                fit = row['p_period']==2\n",
    "            elif row['A_period']>=2 or row['B_period']>=1:\n",
    "                fit = row['p_period']==3\n",
    "        elif row['p_node'] == 'A' :\n",
    "            if  row['B_period']==0:\n",
    "                fit = row['p_period']==1\n",
    "            else:\n",
    "                fit = row['p_period']==2\n",
    "        elif row['p_node'] == 'B' :\n",
    "            if row['A_period']==0:\n",
    "                fit = row['p_period']==1\n",
    "            else :\n",
    "                fit = row['p_period']==2    \n",
    "    return fit "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "id": "aab4ea95",
   "metadata": {},
   "outputs": [],
   "source": [
    "def MoreBetterNaive(row):\n",
    "    fit = 0\n",
    "    if row['network_type']=='a':\n",
    "        if row['p_node'] == 'D':\n",
    "            fit = row['p_period']>=2\n",
    "        elif row['p_node'] == 'A' or row['p_node']=='E' :\n",
    "            fit = row['p_period']>=1\n",
    "        elif row['p_node'] == 'C' or row['p_node'] == 'B':\n",
    "                fit = row['p_period']>=2\n",
    "                \n",
    "    if row['network_type']=='b':\n",
    "        if row['p_node'] == 'C' or  row['p_node'] == 'D' or row['p_node']=='E':\n",
    "            fit = row['p_period']>=3         \n",
    "        elif row['p_node'] == 'A' or row['p_node'] == 'B' :\n",
    "            fit = row['p_period']>=2\n",
    "        \n",
    "    return fit "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "id": "e3dcd2ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "data['Empirical_Optimal']= data.apply(EmpiricalOptimal, axis=1)\n",
    "data['More_Better']= data.apply(MoreBetter, axis=1)\n",
    "data['More_Better_Naive']= data.apply(MoreBetterNaive, axis=1)\n",
    "data['Lone_Wolf'] = data['p_period']==0\n",
    "data['Always_1Turn'] = data['p_period']==1\n",
    "\n",
    "\n",
    "strats = ['Empirical_Optimal','More_Better','More_Better_Naive','Lone_Wolf','Always_1Turn',]\n",
    "n_strats=len(strats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "id": "0d7acbef",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Empirical_Optimal</th>\n",
       "      <th>More_Better</th>\n",
       "      <th>Lone_Wolf</th>\n",
       "      <th>Always_1Turn</th>\n",
       "      <th>network_type</th>\n",
       "      <th>waiting_cost</th>\n",
       "      <th>A_period</th>\n",
       "      <th>B_period</th>\n",
       "      <th>p_period</th>\n",
       "      <th>p_node</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>True</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>b</td>\n",
       "      <td>L</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>a</td>\n",
       "      <td>L</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>b</td>\n",
       "      <td>L</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>a</td>\n",
       "      <td>L</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>a</td>\n",
       "      <td>L</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>A</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3995</th>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>b</td>\n",
       "      <td>L</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>E</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3996</th>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>b</td>\n",
       "      <td>L</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>E</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3997</th>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>a</td>\n",
       "      <td>L</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>E</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3998</th>\n",
       "      <td>False</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>b</td>\n",
       "      <td>L</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>E</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3999</th>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>False</td>\n",
       "      <td>True</td>\n",
       "      <td>a</td>\n",
       "      <td>L</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>E</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>4000 rows × 10 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      Empirical_Optimal  More_Better  Lone_Wolf  Always_1Turn network_type  \\\n",
       "0                  True         True      False         False            b   \n",
       "1                 False         True      False          True            a   \n",
       "2                 False        False      False          True            b   \n",
       "3                 False         True      False          True            a   \n",
       "4                 False         True      False          True            a   \n",
       "...                 ...          ...        ...           ...          ...   \n",
       "3995              False         True      False          True            b   \n",
       "3996              False        False       True         False            b   \n",
       "3997              False         True      False          True            a   \n",
       "3998              False        False       True         False            b   \n",
       "3999              False         True      False          True            a   \n",
       "\n",
       "     waiting_cost  A_period  B_period  p_period p_node  \n",
       "0               L         2         2         2      A  \n",
       "1               L         1         1         1      A  \n",
       "2               L         1         1         1      A  \n",
       "3               L         1         1         1      A  \n",
       "4               L         1         1         1      A  \n",
       "...           ...       ...       ...       ...    ...  \n",
       "3995            L         0         1         1      E  \n",
       "3996            L         2         2         0      E  \n",
       "3997            L         2         1         1      E  \n",
       "3998            L         1         0         0      E  \n",
       "3999            L         1         2         1      E  \n",
       "\n",
       "[4000 rows x 10 columns]"
      ]
     },
     "execution_count": 170,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data[['Empirical_Optimal','More_Better','Lone_Wolf','Always_1Turn','network_type','waiting_cost', 'A_period','B_period', 'p_period', 'p_node']]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d248fab",
   "metadata": {},
   "source": [
    "# Set up the loglikelihood function\n",
    "\n",
    "- [Return](#Content)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "id": "66198829",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Likelhood function takes as an input a vector of proportions of strategies and returns the likelihood value\n",
    "#Note cMat and eMat are global matrices that are updated externally for each treatment.\n",
    "def objective(x,args):\n",
    "    \n",
    "    C = args[0]\n",
    "    E = args[1]\n",
    "    \n",
    "    bc=np.power(x[0],C) #beta to the power of C\n",
    "    be=np.power(1-x[0],E) #beta to the power of E\n",
    "    prodBce = np.multiply(bc,be) #Hadamard product\n",
    "    \n",
    "    #maximum is taken so that there is no log(0) warning/error\n",
    "    res = np.log(np.maximum(np.dot(x[1:],prodBce),np.nextafter(0,1))).sum() \n",
    "    \n",
    "    return -res\n",
    "\n",
    "def constraint1(x):\n",
    "    \n",
    "    return x[1:].sum()-1\n",
    "\n",
    "#Set up the boundaries and constraints\n",
    "b0 = (np.nextafter(0.5,1),1-np.nextafter(0,1))\n",
    "b1 = (np.nextafter(0,1),1-np.nextafter(0,1))\n",
    "bnds = tuple([b0]+[b1]*n_strats) #Beta is at least .5\n",
    "con1 = {'type': 'eq', 'fun': constraint1} \n",
    "cons = ([con1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5fa7958b",
   "metadata": {},
   "source": [
    "# Run likelihood maximization\n",
    "\n",
    "- [Return](#Content)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "id": "c35d8ac8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# data=data[data['round_id']>=15]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "id": "502a6c6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "Count = data.groupby(['subject_id', 'network_type','waiting_cost'])[strats].sum().reset_index()\n",
    "Number = data.groupby(['subject_id', 'network_type','waiting_cost'])['round_id'].count().reset_index()\n",
    "\n",
    "# CC = Count[(Count['network_type']=='b') & (Count['waiting_cost']=='L' )]\n",
    "# NN = Number[(Number['network_type']=='b') & (Number['waiting_cost']=='L' )]\n",
    "CC = Count\n",
    "NN = Number\n",
    "C = CC[strats].to_numpy()\n",
    "N = NN['round_id'].to_numpy()\n",
    "E = np.tile(N.transpose(), (1, n_strats)).reshape((len(N),n_strats))-C\n",
    "C = C.transpose()\n",
    "E = E.transpose()\n",
    "# Count = Count.merge(N)\n",
    "# Count['']\n",
    "# pd.merge([Count,N], keys=['subject_id'])\n",
    "# Count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "id": "726d410c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\optimize.py:282: RuntimeWarning: Values in x were outside bounds during a minimize step, clipping to bounds\n",
      "  warnings.warn(\"Values in x were outside bounds during a \"\n",
      "C:\\Users\\xinxi\\AppData\\Local\\Temp/ipykernel_18340/3692682034.py:9: RuntimeWarning: divide by zero encountered in power\n",
      "  be=np.power(1-x[0],E) #beta to the power of E\n",
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\_numdiff.py:557: RuntimeWarning: invalid value encountered in subtract\n",
      "  df = fun(x) - f0\n"
     ]
    }
   ],
   "source": [
    "#Some random starting point\n",
    "x0 = np.zeros(n_strats+1)\n",
    "x0[0] = .5+.5*np.random.random()\n",
    "temp = np.random.random(n_strats)\n",
    "x0[1:]=temp/temp.sum()\n",
    "\n",
    "bestX=x0\n",
    "bestObjective=objective(x0,[C,E])\n",
    "\n",
    "for k in range(150): #Do many times so that there is low chance of getting stuck in local optimum\n",
    "\n",
    "    x0 = np.zeros(n_strats+1)\n",
    "    x0[0] = .5+.5*np.random.random()\n",
    "    temp = np.random.random(n_strats)\n",
    "    x0[1:]=temp/temp.sum()\n",
    "\n",
    "    #Notice that we are minimizing the negative\n",
    "    solution = minimize(objective,x0,method='SLSQP',bounds=bnds,constraints=cons,args=([C,E]))\n",
    "    x = solution.x\n",
    "    obj = solution.fun\n",
    "    check = solution.success #solution.success checks whether optimizer exited succesfully\n",
    "\n",
    "    if bestObjective>obj and check: \n",
    "        bestObjective=obj\n",
    "        bestX=x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "id": "ea135c08",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                   Estimates\n",
      "beta                  0.7231\n",
      "Always_1Turn          0.2788\n",
      "Empirical_Optimal     0.2632\n",
      "More_Better           0.1803\n",
      "Lone_Wolf             0.1730\n",
      "More_Better_Naive     0.1048\n",
      "LL                -2536.7000\n"
     ]
    }
   ],
   "source": [
    "results=pd.DataFrame(bestX.round(4).tolist()+[np.round(-bestObjective,4)],index=['beta']+strats+['LL'])\n",
    "results=results.rename(columns={0:'Estimates'})\n",
    "results=results.sort_values(by=['Estimates'],ascending=False)\n",
    "print(results)\n",
    "# results.to_csv(\"output\\\\06-sfem_results.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7148022",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77ccb22e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}