ComputingQREForDirectedNetworkb.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "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": "markdown",
   "metadata": {},
   "source": [
    "# This file is used to compute the QRE in network (b). \n",
    "## Note： this is preliminary because only assign prob of deciding at turn 0 and for A to decide at turn 1, in total, eq only contains 6 unknown probs\n",
    "## Check the later version!!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "outputfolder = \"C:\\\\Users\\\\xinxi\\\\Anaconda3\\\\Xinxin's Files\\\\Networks\\\\4O3I_Distance\\\\\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAD3CAYAAAC+eIeLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deViUVf8/8Pes7IsLoiIuoOCgCAqJG2pPWE9shqm55FZpZmalpfZVs/LxcenRvqVZaps9UipIpmiZmIrpzw0FQUZRSUB2EQSB2c/vD2K+jqOmzj1zzzCf13V1dTkyZz4j+Paec5/zOQLGGAMhhBCLEPJdACGE2BMKXUIIsSAKXUIIsSAKXUIIsSAKXUIIsSAKXUIIsSAKXUIIsSAKXTs0c+ZMLFu2jO8yjBw+fBidOnXiuwxOdO3aFWlpaXyXQawQhW4L07VrVzg5OcHNzQ2enp4YNGgQvvzyS+h0Ov3XfPnll1iyZAnnrz116lQsXryY83GbCQQCBAcHG7yXxYsXY+rUqQ/1/OHDh+Orr74yU3WmW79+PcLDw+Hg4PDQ74nYHgrdFmjPnj2oq6tDQUEBFi5ciFWrVuHll19+qOdqNBozV2eakpISbNu2je8y7suUP7+OHTti8eLFeOmllzisiFgbCt0WzMPDA/Hx8di+fTu2bNmCnJwcAIZXpM0f6VetWoX27dtj2rRp0Ol0WLlyJfz9/dGmTRuMHTsWN2/e1I/7xx9/YNCgQfD09ISvry++++47bNq0CYmJiVi9ejVcXV0RFxcHoCkkn3/+eXh5eaFbt2747LPP9OM0NjZi6tSpaNWqFYKCgnD69Om/fU/z58/H0qVL7xtuJ06c0NcWEhKCw4cPAwAWLVqEo0ePYvbs2XB1dcXs2bOxdOlSvPHGGwAAtVoNFxcXzJ8/X1+bo6MjqqurAQC7d+9Gr1694OnpieHDh0Mul+tfs2vXrli1ahX69OkDFxcXo9ouXryIbt26/e0/FqNGjcJzzz2HNm3a/O2fA7FhjLQoXbp0YQcOHDB63NfXl23YsIExxtiUKVPYokWLGGOMHTp0iIlEIjZ//nymUChYQ0MD++STT1hERAQrKipiCoWCzZgxg40bN44xxlhBQQFzdXVlP/zwA1OpVOzGjRvs3LlzRuMyxphWq2X9+vVjH374IVMqlezq1ausW7du7Ndff2WMMbZgwQI2ZMgQVlVVxQoLC1mvXr2Yj4/Pfd8bAJaXl8f69evHNm/ezBhjbNGiRWzKlCmMMcauX7/OWrduzfbu3cu0Wi377bffWOvWrVlFRQVjjLFhw4bpn8cYYwcPHmS9e/dmjDF27Ngx5ufnx/r376//vT59+jDGGLt06RJzdnZmv/32G1OpVGzVqlXM39+fKZVK/Z95SEgIKywsZA0NDQbfh4yMDObr68v27NnzcN/Au94TaXnoStdOdOzY0eBq9U5CoRAffvghHBwc4OTkhI0bN2L58uXo1KkTHBwc8MEHHyA5ORkajQaJiYmIiorC+PHjIZFI0KZNG4SGht5z3NOnT6OyshLvv/8+pFIp/Pz8MH36dP0V344dO7Bo0SK0bt0avr6+mDNnzt++D4FAgGXLluGjjz6CUqk0+L2tW7ciOjoa0dHREAqFGDFiBMLDw7Fv3757jjVw4EBcvnwZVVVVSE9Px8svv4zi4mLcvn0bR44cwbBhwwAA27dvR0xMDEaMGAGJRIJ33nkHjY2NOH78uH6sOXPmwNfXF05OTvrHjh49ivj4eGzZsgWxsbF/+96IfRDzXQCxjOLiYrRu3fqev+fl5QVHR0f9rwsKCpCQkACh8P/+TRaJRCgvL0dRURH8/f0f6jULCgpQUlICT09P/WNarRaRkZEAmqYefH199b/XpUuXhxo3OjoanTt3xqZNm4xeLykpCXv27NE/plar8eSTT95zHCcnJ4SHh+PIkSNIT0/HokWLkJmZiWPHjuHIkSP6qYeSkhKD2oRCIXx9fVFcXKx/7M730ezLL7/EsGHD7vv6xD7Rla4dOH36NIqLizFkyJB7/r5AIDD4ta+vL3755RfU1NTo/1MoFPDx8YGvry+uXr360ON069bNYJy6ujr9lWeHDh1QVFSk//rCwsKHfk//+te/sHz5cjQ0NBi83qRJkwxer76+HgsXLrxnfQAwbNgw/P777zh37hyeeOIJDBs2DPv378epU6cwdOhQAE2fEgoKCvTPYYyhqKgIPj4+933vQFPoFhYW4u23337o90VaPgrdFqy2thapqakYN24cXnzxRQQHBz/U82bOnIlFixbpg6ayshI///wzAGDixIlIS0vDjh07oNFoUFVVhczMTACAt7c38vPz9eP0798f7u7uWLVqFRobG6HVapGTk6O/YTZ27FisWLEC1dXVuH79OtatW/fQ72348OEIDg7Gli1b9I+9+OKL2LNnD/bv3w+tVguFQoHDhw/j+vXr96wPaArd77//HkFBQZBKpfplZd26dYOXl5e+zr179+LgwYNQq9VYs2YNHBwcMGjQoAfW6Obmhl9//RXp6en64H8QjUYDhUIBrVarr9/aV5OQx8D3pDLhVpcuXZijoyNzdXVl7u7ubMCAAWz9+vVMo9Hov+buG2l337zSarVszZo1LCAggLm6ujI/Pz/23nvv6X8/PT2d9e/fn7m5ubFOnTqx7777jjHGWF5eHgsJCWEeHh5s5MiRjDHGiouL2bhx45i3tzfz9PRkERER+ht99fX1bNKkSczDw4PJZDK2evXqv72RdvnyZf2vT5w4wQAY3HQ6ceIEGzp0KGvVqhVr27Yti46OZgUFBYwxxo4fP8569OjBPD092RtvvMEYY6yuro6JxWL2wQcfMMYY0+l0zMvLi82cOdPgtVNSUphMJmPu7u5s6NChLCcnx+DP/O6bl3c+VlVVxfr06cMWL1583/fGGGNLly5lAAz+W7p06QOfQ2yPgDE6OYIQQiyFphcIIcSCKHQJsZDCwkK4urre879HuYlIbBtNLxBCiAXROl1CHqC6XoWkjCKkycshL6uDSqODVCwEYwyMAQJB03Kx5sdl7d0QJfPGmDBftHKR8l0+sUJ0pUvIPWQW1WDz0XykycshEAAKte7vn/QXR4kQjAFRMm9Mj/RDqK/n3z+J2A0KXULuUNOgwnsp2TicVwmlRgudCX87hALAQSzC8AAvrBgVDE9nuvIlFLqE6B3ILce8pEwo1DqotA9/Zft3pCIhHCVCrB0biiiZN2fjEttEoUvsHmMMy/fJkXiyEI1qrdlex0kiwsSIzlgULbvntmFiHyh0iV1jjGFhSjZ2Z5WYNXCbOUlEiA/piJWjgil47RSt0yV27V/75BYLXABoVGuxO6sEy/fJ//6LSYtEoUvs1oHccvxg5imFe2lUa5F4shBp8nKLvi6xDhS6xC7VNKgwLynT4oHbrFGtxbykTNQ0qHh5fcIfCl1il95LyX6ktbfm0KjS4b2fsnmtgVgehS6xO5lFNTicV8npsrDHodLqcPhSJbKKanitg1gWhS6xO5uP5kOp4Wda4W5KjRabj+b//ReSFoNCl9iV6noV0uTlJu0045KOAQfk5aiup7lde0GhS+xKUkYRrG15rEAAJJ+9zncZxEIodIldSZOX834D7W4KtY6Wj9kRau1I7Iq8rM6k5zONGlW/bYDiWiZ0itsQe3ZAq2GT4eQfblpdpbUmPZ/YDrrSJXajvFYBlca0q1ym00Ls1hbtJ6yE79vb4Tn0RVT+vAqaGtOuVBUaHSpqFSaNQWwDhS6xG/LSWkjFpv3IC6WO8IycCLGnNwQCIZy794fYwxvKsismjesgFpp8FU5sA4UusRv1Ki247u+kra+G+mYxpF6dTRqHMYZ6pYajqog1o9AldkOt0YHLzGVaDW7s/g9cg5+CpI2vaWMxmDz1QWwDhS6xGxKxkLPlYozpcCN1DSASo/WImSaPJxDA5KkPYhvou0xapPPnz+PgwYNQKP7v5pSLVMRJD1vGGKr2fQZtfQ28Ev4HApHpi4AEAgFcHGgxkT2g7zJpkTZt2oSNGzdCKBSiX79+CAsLg1fnHqhX+MHUa42b+z+HuqoI3uP+BaHEgZN6lRodZO3dOBmLWDcKXdKi6HQ6nDp1Cn/++Sd0Oh00Gg1OnDiBEydOoF27dnB/aRN0JoSu5lYFbmf+CogkuL5ukv7x1v98Ha69nnzscR3FQrRzd3zs5xPbQaFLbJZOp8Pp06exa9cu/PHHH7h48SKqqqogEAjg4eGh/zoHBwesXr0ab7zxBl7Y9P9w6lr1Y7+m2KMduixM5aJ8A7IO7pyPSawThS6xCTqdDhkZGfjpp59w9OhRXLp0CTdu3AAAtG3bFj179sT06dORkJCAsLAwMMbg4uICFxcX7N27FwMGDAAARMm8cb74llVtBXaUCOmUYDtCB1MSq6PT6XD27Fl9wDZfwTLG0KZNG/Ts2RNDhgzBc889hyeeeAJC4b2nC5KTkzF06FC0a9dO/1h1vQoDVh6E0oqWZzmIhTix8Cm0cpHyXQqxAApdwrvmgE1PT8fFixdx48YNfcAGBgZi8ODBSEhIQP/+/e8bsI/i9R/O4pecUqto7ygUANG9O2D9hH58l0IshKYXiEVlZmYiJSUF6enpkMvl+oBt3bo1AgMDMWXKFCQkJCAiIoKTgL2X6ZF++P1iBW/no93JQSzC9Eg/vssgFkShS8zm/Pnz2Llzpz5gKysr9QEbEBCAyZMn47nnnsPAgQPNFrD3EurrieEBXjh4sYLXI3ukIiGGB3ohxNeTtxqI5dH0AuFETk4OkpOT9QFbUVEBxhhatWqFgIAADB48GM899xwGDRpk0YC9n5oGFYZ+fAi1Cv76HXg4iZH+7j/g4SThrQZieRS65JHl5OTor2Bzc3MNArZHjx4YPHgwRo4ciSFDhlhFwN7PgdxyzNl2jpdpBieJCOvG96VVC3aIQpc80IULF/QBe+HCBVRWVkKn08HT0xMBAQEYNGgQ4uPjMXToUKsOWABQKBQ4efIkfv/9d+zYsQN5eXmY/+MJ7MqpsmjwOklEkBaeRN72lYiJiUFsbCyGDRuGrl27crJNmVg3Cl2iJ5fLsXPnThw5cgQXLlxARUWFPmB79OiBQYMGIS4uDsOHD7f6gL1bVVUVfHx8IJFI0NjYCK1Wiy5duuDPP//EwpRs7M4qsUjwOklEiA/tiKEORYiLi9OvJ1YqlZg1axY+/fRTs9dA+GUVoVteq4C8tBb1Ki3UGh0kYiFcpCIEdXCnrZFmcunSJSQnJ+Pw4cP6gNVqtfqAHThwIOLj4/Hkk0/aXMDez6xZs/DVV19BrVbD2dkZX331FcaPHw/GGJbvkyPxZKFZg9dJIsLEiM5YFC0DAHTr1g0FBQUAABcXF5w7dw49evQw2+sT68DL6oXqehWSMoqQJi+HvKwOKo0OUrEQjDEw1tTmTiAQ6B+XtXdDlMwbY8J8aQH5Y7h8+bI+YHNyclBeXq4P2O7du2P06NGIj4/H8OHDIRa3zAUtSqUStbW1cHRs+kdcKpVi9OjRAJp+1hbHBCGiWxvMS8qEQq3jdFWDVCSEk1SINWNCDeZwFyxYgLlz50Kr1UIikaCkpIRC1w5Y9Eo3s6gGm4/mI01eDoEAj7QV01EiBGNN2zinR/ohlJbZ3NPVq1eRlJSkD9iysjJotVp4eHige/fuGDhwIOLi4vCPf/yjxQbs3W7cuIFRo0ahbdu2+OKLLzBs2DBMmDAB77//vtHX1jSo8F5KNg7nVUKp0Zq0gUIoaFqHOzzQCysSguHpbHjBUF9fj06dOmHNmjXo3LkzJkyYgNWrV2Pq1KmP/6LE6lkkdM3ygxzghRWjjH+Q7Ul+fj6Sk5Px+++/669gNRqNPmAjIiIQHx+Pp556ym4C9m6XLl1CTEwMRo8ejX//+98QCoVQq9UQiUQPnDbJ+usC4YAJFwgj/rpAeNA6XJVKBam06WdYLpcjNjYW48aNw7Jly1rMtA4xZPbQPZBbbraPbI4SIdaODbWLZTd//vknkpKScOjQIf0VrEajgbu7O7p3744BAwYgNjYWI0aMsNuAvduhQ4cwbtw4rFixAi+99NJjjVFdr0Ly2etNU2GltVBodHC4x1SYUqODo1gIWQd3RMm8Mbpfp8eaCqusrERCQgI6dOiA77//Hk5OTo9VN7FeZgtdPm5OtJTlNteuXdNfwWZnZxsErL+/PyIiIhAXF4enn36aAvY+vvnmG7z33nvYtm0bnnzy8fvc3q2iVgF5WR3qlRr9PQcXBzFk7d04u+mrUCjw8ssv48qVK/j555/Rvn17TsYl1sEsocsYs/wynJCOWDkq2OaCt6CgwCBgS0tL9QHr5+eHiIgIxMbG4umnn9Z/DCX3p9Pp8D//8z9ITk7G3r17ERgYyHdJj4UxhmXLluGbb77Bnj17EBwczHdJhCNmCd1le3Pxg5mvcO/WfMW7OCbIYq/5qAoLC40CVq1Ww83NTX8FGx0djX/+858UsI+hoaEBkyZNQmVlJVJSUtC2bVu+SzLZDz/8gLfeegtbtmzBs88+y3c5hAOchy5trWxy/fp1JCcn4+DBgzh//rxBwPr5+aF///6IiYnBs88+SwHLgdLSUsTHx0Mmk2Hz5s1wcODm7DJrcOzYMYwePRqLFi3C7Nmz+S6HmIjT0LWWJiJH3nnSoqsarl+/jp07d+oDtqSkBGq1Gq6urvqAjY6OxrPPPqtfJ0q4k5WVhbi4OMyYMQOLFi2yuSmmh5Gfn4/Y2FhERUVh7dq1NJdvwzgN3de2ZlhFu7ynZO3wxcQws4xfUlKCnTt3Ii0tDVlZWSgtLYVKpYKrqyu6deumD9jo6GgKWAvYu3cvpk6divXr1+OFF17guxyzqqmpwZgxYyCVSrFt2za4udHpwbaIs9DNLKrB+M0nrKIxtJNEhG3TB5jcp7SsrAzJycn6gC0pKTEI2PDwcMTExCAmJoYC1sIYY1i3bh1WrlyJlJQU/RloLZ1arcYbb7yB48ePIzU1FZ07d+a7JPKIOAtdaz8CRalUYvHixXjhhRcQHh5u9JyysjKDK9ji4mKoVCq4uLgYBGx0dDScnZ0t+XbIXTQaDd58800cOXIEqamp6Nq1K98lWRRjDJ988gnWrFmDXbt24YknnuC7JPIIOJkYqq5XIU1ebhWBCwA6BhyQl6O6XoVWLlJcv34dzz77LHJzcyEQCNC5c2eDK9jr16/rA7Zr164YOnQooqOjERsbSwFrZWpra/HCCy9Ap9Ph2LFjBket2wuBQIC5c+eie/fuiI6Oxpdffonnn3+e77LIQ+LkSndT+lWsTcuzumOt540IROuKc3jxxRehUCgANP3ANrfT69KlC8LDwxEdHY24uDgKWCtXUFCA2NhYREZG4rPPPqObSWg61HPkyJF4/fXXsWDBghZ5E7Gl4SR0x248jlPXqk0a48ae/0BxLQs6tQIil1ZwH/A83EKeMWlMQeUVXPv6LYPHHBwcUFlZSTchbMyJEycwatQoLFiwAHPmzKFwuUNxcTHi4uIQGhqKL7/8kpYgWjlOOmrIy+pMHsN9wBj4vPYNOs9NQrvRS1CT/l8oy66YNKZzxx4oLCxEYmIiRo8eDQ8PDyiVSv1VL7EN27dvR1xcHDZu3Ig333yTAvcuPj4+SE9PR1VVFZ5++mncvHmT75LIA5gcuuW1Cqg0pk8rSL26QCBuPqBPAAEE0FSXmjSmmgng4OGFCRMmICkpCdXV1SgsLGwRO5XsAWMMy5cvx7vvvosDBw4gLi6O75KslqurK1JSUhAeHo4BAwbg8uXLfJdE7sPkSTF5aS2kYiGUHARv1f4NqM8+CKZRQurtDyd/41UGj8JBLIS8rE7fiEQgEMDX19fkOon5KZVKzJgxAxcuXMCJEyfQsWNHvkuyeiKRCP/5z38QEBCAyMhIbN++HcOGDeO7LHIXk0O3XqUFV/sr2jwzC61HvApl8UUoCrMhEJl2NDVjDPVK/nbHkcdTVVWFhIQEtG3bFkeOHIGLiwvfJdmUGTNmwM/PD2PGjMHHH3+MKVOm8F0SuYPJ0wtqjQ5cdm8QCEVw9O0Fbd0N1J3bZ9JYjIGTqQ9iOZcuXcKAAQMwcOBAJCcnU+A+pqioKBw5cgQfffQRFi9eDJ2O/h5YC5NDVyIWwiz3NXQ6k+d0BQJAKjZ8i0ql0qQxifkcPnwYQ4cOxcKFC7Fq1So6OcFEMpkMJ06c0Ddzb2xs5LskAg5C10UqMvlusra+BvW5R6BTNYLptGjMz0C9/Agcu4SYNK5AIIBUyPD7779j3rx56N69O5ycnHDjxg2TxiXc+/bbb/HCCy/gxx9/xMsvv8x3OS2Gl5cXDh48CIlEguHDh6OsrIzvkuyeyXO6sg7upn+EFwhQd+4XVO3fADAdxB7t0Oqp6XAOMG0/vVKtxT8HBIM11Px1vAqDo6Mj2rRpY1q9hDN3Nh0/cuQIevbsyXdJLY6joyO2bt2Kjz76CAMGDKCm6DzjZHNE8If7UcdjO8f7cXcU45VWlzB//nz9tIJIJIK3tzfCwsIQFhaG8PBwhIWF0ZEoPGhuOl5RUYGffvqJlvJZQGJiIt5++218//33+Oc//8l3OXaJk0kzWXvr3N0l6+COOXPmICMjAz4+PhAKhZg9ezaOHTuGyZMnQ6FQ4NNPP0VQUBB8fHwwcuRIfPTRR9i3bx/Ky8v5Lr9FKy0txbBhw+Ds7Iy0tDQKXAuZOHEifvrpJ0ybNg2ff/453+XYpRbfe2F6pB+ApkYpM2fOxMyZMzF06FCDr2WM4dq1a8jIyMCZM2eQkZGBjIwMODs766+Em//frl07Pt5Oi5KVlYX4+HhMnz69xTYdt3bNTdFHjBiBtWvXQiQS8V2S3eAkdKvrVRiw8iAnGyS44iAW4sTCpx7rGGygKYj//PNPoyB2dXXVB3DzfxTED2/v3r2YNm0a1q1b1+Kbjlu75qboDg4O+PHHH6kfiYXYTT9dLjDGkJ+fbxDEZ8+ehZubm1EQe3l5cfraLcG6deuwYsUK7Ny5EwMHDuS7HIKmpuizZ8/GiRMnsGfPHmqKbgF0coSJGGO4evWqURB7eHgYTEuEhYXZ7bylRqPBW2+9hUOHDiE1NRXdunXjuyRyB2qKbll0RpoZ6HQ6fRA3h/HZs2fRqlUroyBu6cvX7mw6vmPHDrtsOm4rfv75Z7zyyivUFN3MWuRpwOnv/gMeTqb1beCaTqfDlStXDIL43LlzaN26tVEQt27dmu9yOdHcdHzIkCFYt24dNR23AdQU3fw4DV0AOJBbjjnbzvEyzeAkEWHd+L6Iknlb/LUfh06nw+XLl42CuG3btgbriPv162dzQXzy5EkkJCRg/vz51APXxlBTdPPiPHQB4F97c5F4stCiweskEWFiRGcsjgmy2GuaQ3MQN88PNwdxu3btjIK4VatWfJd7Tzt27MDrr7+Ob775hnrg2qjbt29j4sSJqK2txc6dO23uH31rZpbQZYxhYUo2dmeVWCR4nSQixId2xMqE4BZ5RaXVag2COCMjQx/Ed64j7tevHzw9zX/z8H4YY1ixYgW+/PJL7N69G6GhobzVQkyn1WqxYMEC7NmzB6mpqejRowffJbUIZgld4K+u//vkZr/ibb7CXRQta5GBez9arRZ5eXkGQZyZmQlvb2+jILbEzSulUolXX30VOTk52L17NzUdb0E2bdqE999/Hzt27DDaWEQendlCt9mB3HLMS8qEQq3jdFWDVCSEk1SINWNCbWYO19y0Wi0uXbpksHwtMzMTHTp0MLhRx3UQV1VVYdSoUWjdujW2bt1KPXBboAMHDmDixIn4z3/+g8mTJ/Ndjk0ze+gCTasa3kvJxuG8Sig1WpM2UAgFgINYhOGBXliREAxPZ5rkfxCtVouLFy8aBHFWVhY6duxoFMTu7u6PPH5eXh5iYmIwatQorFixgnrgtmC5ubmIjY3FhAkT8NFHH9H3+jFZJHSbZRXVYPPRfByQl0MgwCP1anCUCMEYMELmjemRfhbZ+NBSaTQaoyA+f/48fHx8DJav9e3b94FBfPjwYbzwwgtYvnw5XnnlFQu+A8KXiooKPPfcc/D19cV3330HJycnvkuyORYN3WbV9Sokn72ONHk55KW1UGh0cBAL/+p523Tig0AggFKjg6NYCFkHd0TJvDG6X6fH7qVAHkyj0UAulxssXzt//jx8fX2NgtjNzQ3ffvstFixYgB9//BFPPfUU3+UTC1IoFHjppZeQn5+Pn3/+Gd7eNL33KHgJ3btV1CogL6tDvVIDlUYHqVgIFwcxZO3d9Cf5EsvTaDTIzc01CmInJyeoVCrMnj0b0dHR6Nu3L1xdXfkul1gQYwwffvghvvvuO6SmpqJ37958l2QzrCJ0iW1oaGjAiy++iGvXrmHKlCn61RM5OTno0qWLwTri0NBQCmI7QE3RHx2FLnkopaWlGDlyJAIDA/HVV1/BwcFB/3tqtRq5ubkGy9dycnLQtWtXoyCmlQ0tz7FjxzB69GgsWbIEs2bN4rscq0ehS/7W+fPnERcXh1deeQWLFy9+qPXQarUaFy5cMApiPz8/oyB2dna2wLsg5pSfn4+YmBg888wzWLNmDTVFfwAKXfJA+/btw9SpU/HZZ59h3LhxJo2lUqmMgvjChQvw8/Mz2NAREhJCQWyDampqMHr0aDg6OlJT9Aeg0CX3tW7dOvz73/9GSkqK2ZqOq1Qq5OTkGCxfy83Nhb+/v8E6Ygpi26BWq/H666/j5MmTSE1Nha+vL98lWR0KXWKE76bjSqXSKIjlcjm6d+9uFMS0TtT6MMawdu1afPLJJ9i1axfCw8P5LsmqUOgSA81Nx7VaLZKSkqym6bhSqUR2drZBEF+8eBE9evQwWEccEhICR0daZmgNdu3ahenTp2Pjxo0YNWoU3+VYDQpdondn0/HPPvsMEol1NYK/m0KhMAriS5cuISAgwCCI+/TpQ0HMk+am6LNnz8b8+fPtqinV/VDoEgAtp+m4QqHA+fPnDTZ05OXlITAw0CiI71z2RsynuSl637598cUXX9h9U3QKXYKkpCTMmjWrxTYdb2xsNAriy5cvo2fPngbL14KDgymIzaS5KXpdXR127txptQ34LYFC1x9vFGEAAB5+SURBVI4xxrBy5Up88cUXdtd0vDmI71y+dvnyZchkMoMg7t27NwUxR7RaLebPn4/U1FTs3bsX3bt357skXlDo2imVSoUZM2YgOzsbe/bsoabjaArirKwsgyC+cuUKgoKCjILY3j8im2Ljxo1YunQpkpKSEBkZyXc5Fkeha4eo6fjDa2hoMAriq1evIigoyGD5GgXxo2luir5mzRpMmjSJ73IsikLXzjQ3HU9ISMDKlSupEfVjqK+vR1ZWlsGqifz8fPTq1csoiK19BQifmpuiT5w4ER9++KHd/CxS6NqR5qbj//rXvzB9+nS+y2lR6uvrkZmZaRDE165dMwriXr16URDfobkpeufOnfHtt99yttmlvFYBeWkt6lVaqDU6SMRCuEhFCOrgznu7WApdO/Htt99i4cKF+OGHH6jpuIXcvn3bKIgLCgrQu3dvg+VrQUFBdh3ECoUC06ZNw7Vr17Br167HaopeXa9CUkZR08EIZXX6vtx3H4zQ/LisvRuiZN4YE+Zr8YMRKHRbOJ1Oh8WLF2PHjh1ITU1Fz549+S7Jrt2+fRvnzp0zWL5WWFiI4OBgoyAWi8V8l2sxjDF88MEH+P7775GamopevXo91PMy/zoCLM2EI8Ci/joCLNRCR4BR6LZgDQ0NmDJlCkpLS7Fr1y60bduW75LIPdTV1RkFcVFREfr06WOwakImk7X4IG5uiv7f//4XzzzzzH2/ziyH3QZ4YcUo8x92S6HbQpWVlSE+Ph4BAQH4+uuvaa2pjWkO4uZpiTNnzqC4uNgoiHv27NnigviPP/7AmDFj8P777+O1114z+v0DueWYl5QJhVoHlfbhr2z/jlQkhKNEiLVjQxElM9+5bxS6LVB2djbi4uLw0ksvYcmSJTa7pZcYqq2tNQrikpIShISEGAWxrTcRv3r1KmJjYw2aojPGsHyfHIknC9Go1prttZ0kIkyM6IxF0TKz/N2h0G1hmpuOf/rppxg/fjzf5RAzu3XrlkEQZ2Rk6IP4zlUTthjE1dXVGD16NJydnZGYmIjlB65hd1aJWQO3mZNEhPiQjlg5Kpjz4KXQbUGam47v3LkTgwYN4rscwpOamhqjIC4rKzMK4sDAQKsPYrVajVmzZiG9ti0EAZFQaCwXV81XvItjgjgdl0K3BdBoNHj77bdx8OBB7N271+JNx4n1q6mpwdmzZw2Wr5WXlyM0NNQgiAMCAqwuiH/LLcOs/56GBpbfPOEkEWHd+L6czvFS6Nq42tpajBs3DhqNBjt27ICnp2WWvRDbV11dbRTElZWVCA0NNVi+FhAQwNtusZoGFYZ+fAi1Cg0vrw8AHk5iHHnnSc5WNVDo2rDCwkLExsZi0KBBWLdunV0vsCfcuHnzplEQ37hxA3379jUI4h49elgkiF/bmoGDFys4XaXwqKQiIZ6StcMXE8M4GY9C10adOnUKCQkJeOedd/DWW2/RCgViNlVVVfogbg7jqqoq9OvXzyCIu3fvzmkQZxbVYPzmExa5cfZ3nCQibJs+ACEcbKCg0LVBycnJeO211/D1118jPj6e73KIHaqqqtKHcHMQV1dX64O4OYz9/f0fO4hf/+EsfskpNWnjA1eEAiC6dwesn9DP5LEodG1Ic9PxDRs2YPfu3ejbty/fJRGid+PGDZw9e9ZgHfGtW7fuGcR/98msul6FASsPQqnhb1rhbg5iIU4sfMrkXg0UujZCpVLh1VdfRVZWFvbs2QMfHx++SyLkb1VWVhoEcUZGhj6I71w1cXcQb0q/irVpeY/US8HcHCVCzBsRiOmRfiaNQ6FrA6qqqvD888/D09MTiYmJ1HSc2LTKykqDG3UZGRmoq6szCOL/lrZFdlnjY79GbcYe1GcfhKryGlxkw9A29m1Oao/o1hrbZww0aQwKXSt3+fJlxMTEYOTIkVi5cqXVraEkhAsVFRUG88Pn/CdCIHV+7PEaLh0HBAI0/nkWTK3iLHTdHcU4v/T+jXgeBoWuFTty5AjGjh2LZcuWYcaMGXyXQ4hFlNcqMPTjQ5zM51an/xfa2hucha5ULMQf7z5pUiN0+zgfwwZ99913GDNmDBITEylwiV2Rl9ZCKrbOaHIQCyEvqzNpjJbVE64F0Ol0WLJkCbZt24YjR45AJpPxXRIhFlWv0sJaP4AzxlCvNG13HIWuFWlsbMTkyZNRWlqKEydOwMvLi++SCLE4tUYHK81cMAaoTJz2sM5reDtUVlaG4cOHQyqVIi0tjQKX2C2JWAhr3WApEMDkqQ8KXSuQnZ2NAQMGICYmBlu3boWjI7+nlRLCJxepyORt7UynBdOoAJ0WYDowjQpMZ/p2YoFAABcH0yYIaHqBZ7/88gumTJmC//3f/8WECRP4LocQXjDGcO3aNWRkZODomfO4zcIA4ePH061j23Dr2I/6X9dfOASPwePhGTnRpDqVGh1k7d1MGoNCl0fr16/H8uXL8dNPP2Hw4MF8l0OIRTDGUFBQYLA5IiMjA46OjvrNEY4qMRpNuDD1jJxocsDei6NYaNJyMYBClxcajQZz585FWloajh07Bj8/07YVEmKtGGMoLCw06MeQkZEBBwcHfcC++eabCAsLQ/v27fXPy9l4HKeuVfNY+b3JOribPAaFroXV1dVh3LhxUKlUOH78ODUdJy1Gc8DeubMsIyMDEolEH7Bz5sxBWFgYOnTo8MCxomTeOF98y+p6L3BxggTtSLOg5qbjAwcOxPr166npOLFZjDEUFRUZBaxIJNIHbPP/O3bs+Mjjt+QuY3SlayGnT5/Gc889h3nz5uHtt9+mpuPEZjDGcP36daOAFQqF+nCdNWuWPmC5+Nlu5SJFlMzbqvrpjpB5mxy4AIWuRezcuRMzZ86kpuPE6jHGUFxcbNQFDIA+YGfOnImwsDD4+PiY9eJheqQffr9YYRUnRziIRSa3dGxGoWtGjDGsWrUKn3/+Ofbv349+/UzvOk8IVxhjKCkpMQjYM2fOgDGmD9gZM2YgLCwMnTp1svins1BfTwwP8LKKM9KGB3pxclQPQHO6ZqNSqTBz5kxkZmZS03FiFe4VsFqt1mD+NTw8nJeAvR9rOQ04/d1/wMOJm3swdKVrBjdv3sTzzz8Pd3d3pKenw9XVle+SiJ0pLS01Cli1Wq0P15deegkbNmyAr6+v1QTsvTiLgTD1Bfyu7gqBxMHir+8kEWHNmFDOAhegK13OXb58GbGxsYiLi8OqVauo6Tgxu7KyMoP51zNnzkCpVBqtIujcubNVByzQtIb91KlTOHToEHbv3o3Tp09DKBRi7tY/sCe3xqLzu04SESZGdMbimCBOx6XQ5RA1HSfmVl5ebhSwjY2NRgHbpUsXqw/Ye9myZQumTp0KiUQCtVoNAFi1ahXeffddLEzJxu6sEosEr5NEhPjQjliZEMz5nyOFLke2bNmCd999F4mJiRgxYgTf5ZAWoLy83GiZVkNDg8HJumFhYejatatNBuy9qNVq9OvXDzk5OQAAFxcXVFRUwNnZGYwxLN8nR+LJQrMGb/MV7qJomVn+XGlO10Q6nQ7vv/8+fvzxRxw+fBhBQdx+FCH24e4zwjIyMnD79m19uL744ov45JNP0K1btxYTsPfyzTff4MKFC3B2doZSqcQrr7wCZ+ems9IEAgEWxwQholsbzEvKhEKt43RVg1QkhJNUiDVjQjnZeXY/dKVrgsbGRkyZMgXFxcXYtWsX9cAlD6X5NNw7A/bu03DDwsLg5+fXogP2brNnz8aGDRvwwQcfYNKkSYiOjsYvv/yCrl27Gn1tTYMK76Vk43BeJZQarUkbKISCpnW4wwO9sCIhGJ7Opm+AeBAK3cdUXl6OkSNHwt/fH19//TX1wCX3dOPGDaONBrdu3TIKWH9/f7sK2DvpdDo89dRTOHr0KLZv347nn3/+oZ+bVVSDzUfzcUBeDoEAj9SrwVEiBGNNO82mR/pxtg7371DoPoacnBzExsZi2rRpeP/99+32LwsxVFVVZRSw1dXVBgEbHh4OPz8/CIV0fgAA1NTUoG/fvqisrMTx48fRp0+fxxqnul6F5LPXkSYvh7y0FgqNDg5iIRhjYKzpxAeBQAClRgdHsRCyDu6IknljdL9OnGztfRQUuo/o119/xeTJk/HJJ59g4kTu+3US23Dz5k2jdbDNAXvnjS5/f38K2PuQy+WIiIiAh4cHsrKy0Lp1a87GrqhVQF5Wh3qlBiqNDlKxEC4OYsjau5ncD9dUFLqPYMOGDVi2bBmSk5Op6bgduXnzJs6ePWsQsFVVVUYB2717dwrYh5SamoqEhARERETg8OHDEIvt554+he5D0Gq1mDt3Ln777Tfs3buXmo63YNXV1QanGZw5cwY3btxA3759DZZp9ejRgwL2Ma1evRoLFy7Eyy+/jM2bN/NdjsVR6P6Nuro6jB8/HgqFAsnJydR0vAWprq7G2bNnDaYJKisrERoaahCwAQEBFLAcmTx5MhITE7F27Vq8+eabfJfDCwrdBygsLERcXBwiIiLw+eefU9NxG1ZTU2MUsOXl5QgNDTVYRRAQEEBbt81ApVJh8ODByMrKwt69e+16AxGF7n00Nx2fO3cu5s6dSysUbMitW7eMArasrAwhISEGARsYGEgBawEVFRUICQlBY2MjMjIy4O/vz3dJvKLQvYeUlBS8+uqr+OqrrzBy5Ei+yyEPUFtbaxSwJSUlRgHbs2dPClgenDlzBkOHDoWPjw/OnTtHHfdAoWuAMYbVq1dj/fr1+Pnnn6npuJWpra3FuXPnDNbBFhcXo0+fPgbrYClgrcOPP/6ISZMmISoqCvv27aN58b9Q6P5FpVLhtddew7lz56jpuBWoq6szCtiioiKEhIQYLNPq2bOnXS03shWLFi3CihUr8Pbbb2PNmjV8l2NVKHTRdBf7+eefh5ubGxITE+kjkIXdvn3bIGDPnDmDoqIi9OnTxyBgZTIZBayV0+l0SEhIQGpqKr766itMmzaN75Ksjt2H7pUrVxATE4PY2FisXr2aPpaaWXPA3rkOtrCwEMHBwUYBS6tFbEtDQwPCw8Px559/4vfff8fAgQP5Lskq2XXopqenY+zYsfjwww/x6quv8l1Oi1NfX28UsAUFBejdu7fBOtigoCAKWBtXUFCAfv36QSQS4ezZs+jUqRPfJVktuw3d77//Hu+88w41HedIfX09MjMzDVYRXLt2Db169TII2F69elHAtjDp6ekYMWIEevbsiZMnT1LHvb9hd6Hb3HT8hx9+QGpqKjUdfwwNDQ1GAZufn49evXoZLNPq1asXpFLLdnAilrVx40bMmjULo0aNQlJSEt/l2AS7Ct3GxkZMnToVRUVF2LVrF9q1a8d3SVavoaEBWVlZBgF79epVBAUFGQRs7969KWDtzBtvvIHPP/8cS5cuxdKlS/kux2bYTeg2Nx338/PDN998Qx+B7qGxsdEoYK9cuQKZTGYUsA4Olj8Om1iHO5uOb9u2DaNHj+a7JJtiFaFbXquAvLQW9Sot1BodJGIhXKQiBHVw56T3ZU5ODuLi4jBlyhQsXbqUtvSiKWDPnz9vsA728uXL6Nmzp8FGAwpYcieumo7bM15Ct7pehaSMoqYu72V1+ibDd3d5b35c1t4NUTJvjAnzfeQu7/v378ekSZPsuum4QqEwCti8vDz07NnTYJlWcHAwBSy5rzubjmdmZqJNmzZ8l2STLBq6mX+dZ5RmwnlGUX+dZxT6EOcZffHFF/joo4+QlJSEIUOGmFK6zVAoFMjOzjbYaJCXl4fAwECjgKUpFvKw9u3bh5EjR+KJJ55Aeno6bVIxgUVC1ywndwZ4YcWoe5/cqdVqMW/ePOzfvx+pqakttquRUqnE+fPnDdbBXrp0CQEBAQYB26dPHwpY8tg+/vhjLFiwANOmTcPXX3/Ndzk2z+yheyC33Gxn1DtKhFg71vCM+rq6OkyYMAENDQ1ITk5Gq1atOHtNPimVSmRnZxsE7MWLF9GjRw+DdbB9+vSBk5MT3+WSFmLy5MnYunUrPvnkE7ttOs41s4UuYwzL98mReLIQjWqtOV4CAOAkEWFiRGcsipbh+vXriI2NRf/+/bFhwwabXYSvUqmMAlYul6N79+4GARsSEkIBS8xCrVZj8ODByMzMxJ49e/DMM8/wXVKLYZbQZYxhYUo2dmeVmDVwmzlJRBjoI8Fvy6bgrbfewrx582xmhYJKpUJOTo7BMq3c3Fz4+/sbBayzszPf5RI7cGfT8dOnT6NHjx58l9SimCV0l+3NxQ9mvsK9G1Mr8KSvBN+9GW+x13xUKpUKFy5cMArYbt26GayDDQ0NpYAlvMjIyEBkZCQ6duyIzMxM6rhnBpyH7oHccszZds6igdvMSSLCuvF9DeZ4+aJWq3HhwgWDZVo5OTn3DFgXFxe+yyWEmo5bCKehW9OgwtCPD6FWoeFqyEfm4STGkXeevOeqBnNRq9XIzc01CtguXboYbDSggCXWasmSJVi+fDk1HbcATkP3ta0ZOHixgtNVCo9KKhLiKVk7fDExzCzjazQao4DNzs5Gly5dDJZphYaG0kczYvV0Oh1GjRqFPXv2YPPmzXjppZf4LqnF42yFc2ZRDQ7nVfIauACg0upw+FIlsopqEPIQGygeRKPRQC6XG2w0yM7ORufOnfUBO27cOISGhsLNzY2jd0CIZTQ0NOCJJ55Afn4+0tPTMXjwYL5Lsguche7mo/lQaiw/j3svSo0Wm4/mY/2E/ztYUqvVYtWqVYiLi0NwcLDRc5oD9s5lWtnZ2ejUqZP+6nXs2LHo27cvBSyxeQUFBQgLC4NQKMTly5ep6bgFcRK61fUqpMnLTdppxiUdAw7Iy1Fdr0IrFykqKysRHx+PU6dO4fbt2wgKCjIK2PPnz8PHx0cfsKNHj0bfvn3h7u7O99shhFPNTccDAwNx6tQp2q1oYZzM6W5Kv4q1aXmP1EvB3BwlQswbEYje4nLExMSgtrYWWq0Wbm5u0Ol06Nixo8E62L59+8LDw4Pvsgkxq02bNuG1115DQkICkpOT+S7HLnFypZsmL+cscNU3i1Hy9Wy49ByMtnHvPPY4CrUOezMLMWPOkwaPa7VaFBcXw9PTtPleQmwNNR23DpyErrysjothAAA3f/sSDh242QGTf1OJCxcu4NChQ0hJScHx48fR0NAAjYa/JW2EWJpOp0NUVBTS09OxY8cOajrOM5NDt7xWAZWGm6vc+twjEDq6QNKmJzQ1pSaPp9Do0LaTH15/PQivv/461Go1Ll68SH1Aid2ora1FaGgoKioqkJGRgZCQEL5LsnsmbzmRl9ZCKjZ954pO2YCao4lo9Y+XTR6rmYNYaHAVLpFIEBwcbDN9GQgxxaVLl+Dr6wuVSoVr165R4FoJk9OyXqUFF/sratL/C9eQpyF29zJ5rGaMMdQraSqB2J99+/ahd+/e6NWrF65du4a2bdvyXRL5i8mhq9boYGrmqsrzoSjIgvsTI00txwBj4GzqgxBb8fHHHyM2NhaTJk3C8ePH6ZQHK2Pyd0MiFsLUT+uKwmxobpXj+oZpAACmUgBMh9Ibb6LDtE8fe1yBAJxMfRBiK5qbjq9duxZvvfUW3+WQezA5dF2kIpPnSF1Dn4GLbKj+17WnUqC5VY7Wz7xu0rgCgQAuDvSvPGn5NBoNBg0ahMzMTPzyyy/UdNyKmZxIsg7uJn+EF0ocAcn/7YoRSBwhEEshcjZts4JSo4OsPW3ZJS1bRUUFQkND0dDQgAsXLlDTcStncuh6uztCKhZCyeHcqWckN0elO4qFaOdOWxxJy3X27FkMGTIEHTt2xKVLl6gviA3gZMLTWq8mZR2obwJpubZv347+/fsjMjISeXl5FLg2gpPQjZJ5w1FiXTesHCVCqzhBghBzWLJkCcaPH485c+Zg//79dMqDDeHkOzUmzNfkZWNcYwwY3Y/a1ZGWRafTISEhAf/+97+xefNmrF27lu+SyCPi5NZ+KxcpomTe+CWn1CraOwoFwAiZN1q5WO7IHkLMrbGxEeHh4bh69So1HbdhnH0mmR7pBwexiKvhTOIgFmF6pB/fZRDCmcLCQvj6+qKiogJXrlyhwLVhnIVuqK8nhgd4QSrid25JKhJieKCXyUf1EGIt0tPTERAQgI4dO6KoqIhOebBxnCbkilHBvN9Qc5IKsXJUH15rIIQrmzZtwpNPPonY2FicP3+eTnloAThNSE9nKdaMCYWThJ9pBieJCGvGhMLDScLL6xPCpTlz5mDmzJlYsmQJnfLQgnC+R3ZEkDcmRnRG4slCNKotd1Clk0SEiRGdaZkYsXk6nQ4jRozAkSNHqOl4C2SWxgSLomWoU2iwO6vEIsHrJBEhPrQjFkXLzP5axL6V1yogL61FvUoLtUYHiVgIF6kIQR3cOdn9eGfT8TNnziA0NJSDqok1MUvoCgQCrBwVDDdHsdmveJuvcBdFy6g5OeFcdb0KSRlFSJOXQ15WB5VGB6lYCMYYGGvqZCcQCPSPy9q7IUrmjTFhvo+8ZPHSpUvo378/3NzcqAduC8bJacAPciC3HPOSMqFQ66DSctefQSoSwkkqxJoxoTSlQDiXWVSDzUfzkSYvh0CARzp41VEiBGNNOzWnR/oh9D4radRqNSIiIvDxxx9DqVRi5MiRCA8Px9GjR6kHbgtm9tAFgJoGFd5LycbhvEooNVqTNlAIBU3rcIcHemFFQjA8nWkDBOGOWX5WA7ywYpTxz2pSUhImTZoExhjUajWmTp2Kb775xsR3QKydRUK3WdZfVw8HTLh6GPHX1QOtwyVcM+enMkeJEGvHGn4q69u3LzIzMwEAbm5uyM/PpykFO2DR0G1WXa9C8tnrTfNkpbVQaHRwuMc8mVKjg6NYCFkHd0TJvDG6Xyfa2ks4xxjD8n1yi95/yMjIQP/+/cEYg0AggEgkwiuvvIIvvvjCbK9PrAMvoXu3iloF5GV1qFdq9DckXBzEkLV3o364xKwYY1iYkm3RlTbP9vLC968/g6qqKvj7+2PkyJGIiopCZGQkXF1dzV4D4ZdVhC4hfFm2Nxc/WHhNOTQqtL99Bbs/nIJ27dpZ7nWJVaAmnMRuHcgtt3zgAoBYiltewThfRdc79ohCl9ilmgYV5iVlWj5w/9Ko1mJeUiZqGlS8vD7hD4UusUvvpWQ/0uoZc2hU6fDeT9m81kAsj0KX2J3MohoczqvkdFnY41BpdTh8qRJZRTW81kEsi0KX2J3NR/Oh1PAzrXA3pUaLzUfz+S6DWBCFLrEr1fUqpMnLreJYKQDQMeCAvBzV9TS3ay8odIldScoogrX1RRIIgOSz1/kug1gIddUgdiVNXm7yDbSyxIVQllyCQNjUrF/k1gY+MzY+9ngKtQ5p8nI6189OUOgSuyIvq+NknNZPz4RbyDOcjAUA8tJazsYi1o2mF4jdKK9VQKXhd8XC/Sg0OlTUKvgug1gAXekSuyEvrYVULISSg+CtObwFNYe3QNLaB55DJ8Gxi2mHoTqIhZCX1VGvETtAoUvsRr1KCy5ajbR6chokbXwhEElQL09Hxc5l6DDtM0hadXjsMRljqFdqTK6NWD+aXiB2Q63RgYv2Tg4dAyF0cIZALIFr8FNw8JGh8eoZk8ZkDFY79UG4RaFL7IZELDTPcjGBAIBpaS4QAFIx/XW0B/RdJnbDRSoy+fBSneI2GvMzwDQqMJ0Wty8cgrIoB07d+pk0rkAggIsDzfbZA/ouE7sh6+Bu8kd4ptOiJn0r1DevAwIhJG06wWvUYkjadDJpXKVGB1l7N5PGILaBQpfYDW93R5NXL4icPdBh6iccVtXEUSyklQt2gqYXiF2x1qtJWQd3vksgFkKhS+xKlMwbjhLr+rF3lAgNTgkmLZt1/fQRYmZjwnw5WTbGJcaA0f1MmxMmtoNCl9iVVi5SRMm8IbSSTmNCATBC5o1WLlK+SyEWQqFL7M70SD84iEV8lwEAcBCLqLuYnaHQJXYn1NcTwwO8IBXx++MvFQkxPNALIb6evNZBLItCl9ilFaOCeb+h5iQVYuUo0xrlENtDoUvskqezFGvGhMJJws80g5NEhDVjQuHhJOHl9Ql/KHSJ3RoR5I2JEZ0tHrxOEhEmRnSmZWJ2ikKX2LVF0TLEh3S0WPA6SUSID+2IRdEyi7wesT4CxkWDUUJsGGMMy/fJkXiyEI1q8x3N3nyFuyhaZnLjHWK7KHQJ+cuB3HLMS8qEQq2DSstdb1upSAgnqRBrxoTSlAKh0CXkTjUNKryXko3DeZVQarTQmfC3QyhoWoc7PNALKxKC4elMGyAIhS4h95RVVIPNR/NxQF4OgQCPdGy7o0QIxpp2mk2P9KN1uMQAhS4hD1Bdr0Ly2etIk5dDXloLhUYHB7EQjDEw1nTig0AggFKjg6NYCFkHd0TJvDG6Xyfa2kvuiUKXEEIsiJaMEUKIBVHoEkKIBVHoEkKIBVHoEkKIBVHoEkKIBf1/GkiLr+X05ecAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAD3CAYAAAC+eIeLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXyM5/438M+s2SRiiSCJJSHJJCKTpYLW0tN42sf6C6VUFVWqqrroQlFtHQftQ09Lceh6ftQSdVpLVAVJnDqW7AlDVJpFJJNFIpFk9uv5QzPHmFAy98w9k/m+X6+++hLmur9DfNxz3df1vQSMMQZCCCE2IeS7AEIIcSYUuoQQYkMUuoQQYkMUuoQQYkMUuoQQYkMUuoQQYkMUuoQQYkMUuk5o/vz5WLVqFd9lmElJSYG/vz/fZXCiT58+SE5O5rsMYocodNuZPn36wM3NDZ6envD29sbQoUOxdetWGAwG46/ZunUrVqxYwfm1Z82aheXLl3M+bguBQICIiAiT97J8+XLMmjXrgV4/cuRIfPnll1aqzjJqtRpz5sxB79694enpiaioKBw5coTvsogVUOi2QwcPHkRDQwOKi4uxZMkSrFu3DnPmzHmg1+p0OitXZ5nr169j9+7dfJdxT239/dPpdAgICEBqaipu3ryJVatWYcqUKSgqKuK2QMI7Ct12rGPHjhg/fjz27NmD7777Dvn5+QBM70hbPtKvW7cO3bt3x+zZs2EwGLB27VoEBQWhS5cumDJlCm7cuGEc99///jeGDh0Kb29vBAQE4Ntvv8W2bduwc+dOfPzxx+jQoQPGjRsH4HZITpo0CT4+Pujbty8+//xz4zjNzc2YNWsWOnXqhLCwMJw/f/5P39M777yDlStX3jPczpw5Y6wtMjISKSkpAIBly5bh1KlTWLhwITp06ICFCxdi5cqVePXVVwEAWq0WHh4eeOedd4y1ubq6ora2FgBw4MABhIeHw9vbGyNHjoRCoTBes0+fPli3bh0GDhwIDw8Ps9ouXbqEvn373vcfCw8PD3zwwQfo06cPhEIhxo4di759+yIjI+NPf0+Ig2GkXenduzc7duyY2dcDAgLY5s2bGWOMzZw5ky1btowxxtjJkyeZSCRi77zzDlOpVKypqYl9+umnLC4ujpWWljKVSsXmzZvHpk6dyhhjrLi4mHXo0IF9//33TKPRsOrqapaVlWU2LmOM6fV6Fh0dzT788EOmVqvZ1atXWd++fdnPP//MGGPs3XffZY899hirqalhJSUlLDw8nPn5+d3zvQFgBQUFLDo6mm3fvp0xxtiyZcvYzJkzGWOMXbt2jXXu3JkdPnyY6fV69ssvv7DOnTuzyspKxhhjI0aMML6OMcaOHz/OBgwYwBhj7Ndff2WBgYFs0KBBxp8bOHAgY4yxy5cvM3d3d/bLL78wjUbD1q1bx4KCgpharTb+nkdGRrKSkhLW1NRk8ueQkZHBAgIC2MGDBx/sD/APFRUVzMXFhSkUiod6HbF/dKfrJHr27Glyt3onoVCIDz/8EC4uLnBzc8M//vEPrF69Gv7+/nBxccEHH3yAffv2QafTYefOnYiPj8e0adMgkUjQpUsXyOXyVsc9f/48qqqq8P7770MqlSIwMBBz58413vHt3bsXy5YtQ+fOnREQEIBFixb96fsQCARYtWoVPvroI6jVapOf27FjB0aPHo3Ro0dDKBRi1KhRiI2NRVJSUqtjDRkyBFeuXEFNTQ3S0tIwZ84clJWV4datW0hNTcWIESMAAHv27MGYMWMwatQoSCQSvPXWW2hubsbp06eNYy1atAgBAQFwc3Mzfu3UqVMYP348vvvuO4wdO/ZP31sLrVaL6dOnY+bMmQgNDX3g1xHHIOa7AGIbZWVl6Ny5c6s/5+PjA1dXV+OPi4uLkZCQAKHwv/8mi0QiKJVKlJaWIigo6IGuWVxcjOvXr8Pb29v4Nb1ej2HDhgG4PfUQEBBg/LnevXs/0LijR49Gr169sG3bNrPrJSYm4uDBg8avabVaPP74462O4+bmhtjYWKSmpiItLQ3Lli1DdnY2fv31V6SmphqnHq5fv25Sm1AoREBAAMrKyoxfu/N9tNi6dStGjBhxz+u3xmAwYMaMGZBKpdi0adMDv444DrrTdQLnz59HWVkZHnvssVZ/XiAQmPw4ICAAR44cQV1dnfE/lUoFPz8/BAQE4OrVqw88Tt++fU3GaWhoMN559ujRA6WlpcZfX1JS8sDv6a9//StWr16NpqYmk+vNmDHD5HqNjY1YsmRJq/UBwIgRI3DixAlkZWXhkUcewYgRI3D06FGcO3cOw4cPB3D7U0JxcbHxNYwxlJaWws/P757vHbgduiUlJXjjjTce6D0xxjBnzhwolUr88MMPkEgkD/abQRwKhW47Vl9fj0OHDmHq1Kl47rnnEBER8UCvmz9/PpYtW2YMmqqqKvz0008AgOnTpyM5ORl79+6FTqdDTU0NsrOzAQC+vr4oLCw0jjNo0CB4eXlh3bp1aG5uhl6vR35+vvGB2ZQpU7BmzRrU1tbi2rVr2Lhx4wO/t5EjRyIiIgLfffed8WvPPfccDh48iKNHj0Kv10OlUiElJQXXrl1rtT7gduj+85//RFhYGKRSqXFZWd++feHj42Os8/Dhwzh+/Di0Wi3Wr18PFxcXDB069L41enp64ueff0ZaWpox+O/n5ZdfhkKhwMGDB02mKUg7w/ekMuFW7969maurK+vQoQPz8vJigwcPZps2bWI6nc74a+5+kHb3wyu9Xs/Wr1/PgoODWYcOHVhgYCBbunSp8efT0tLYoEGDmKenJ/P392fffvstY4yxgoICFhkZyTp27MgmTJjAGGOsrKyMTZ06lfn6+jJvb28WFxdnfNDX2NjIZsyYwTp27MhkMhn7+OOP//RB2pUrV4w/PnPmDANgfJDW8rXhw4ezTp06sa5du7LRo0ez4uJixhhjp0+fZv3792fe3t7s1VdfZYwx1tDQwMRiMfvggw8YY4wZDAbm4+PD5s+fb3Lt/fv3M5lMxry8vNjw4cNZfn6+ye/53Q8v7/xaTU0NGzhwIFu+fPk931tRUREDwFxcXJiHh4fxvx07dtzzNcQxCRijkyMIIcRWaHqBEEJsiEKXEBspKSlBhw4dWv3vYR4iEsdG0wuEEGJDtE6XkPuobdQgMaMUyQolFBUN0OgMkIqFYIyBMUAguL1crOXrsu6eiJf5YnJMADp5SPkun9ghutMlpBXZpXXYfqoQyQolBAJApTX8+Yv+4CoRgjEgXuaLucMCIQ/w/vMXEadBoUvIHeqaNFi6Pw8pBVVQ6/QwWPC3QygAXMQijAz2wZqJEfB2pztfQqFLiNGxi0osTsyGSmuARv/gd7Z/RioSwlUixIYpcsTLfDkblzgmCl3i9BhjWJ2kwM6zJWjW6q12HTeJCNPjemHZaFmr24aJc6DQJU6NMYYl+/NwIOe6VQO3hZtEhPGRPbF2YgQFr5OidbrEqf01SWGzwAWAZq0eB3KuY3WS4s9/MWmXKHSJ0zp2UYnvrTyl0JpmrR47z5YgWaG06XWJfaDQJU6prkmDxYnZNg/cFs1aPRYnZqOuScPL9Ql/KHSJU1q6P++h1t5aQ7PGgKX/yuO1BmJ7FLrE6WSX1iGloIrTZWFtodEbkHK5CjmldbzWQWyLQpc4ne2nCqHW8TOtcDe1To/tpwr//BeSdoNClziV2kYNkhVKi3aaccnAgGMKJWobaW7XWVDoEqeSmFEKe1seKxAA+zKv8V0GsREKXeJUkhVK3h+g3U2lNdDyMSdCrR2JU1FUNFj0eqbTouaXzVAVZcOgugWxdw90GvE83IJiLaurvN6i1xPHQXe6xGko61XQ6Cy7y2UGPcSeXdH92bUIeGMPvIc/h6qf1kFXZ9mdqkpnQGW9yqIxiGOg0CVOQ1FeD6nYsm95odQV3sOmQ+ztC4FACPd+gyDu6At1xW8WjesiFlp8F04cA4UucRqNGj247u+kb6yF9kYZpD69LBqHMYZGtY6jqog9o9AlTkOrM4DLzGV6HaoP/D90iHgCki4Blo3FYPHUB3EMFLrEaUjEQs6WizFmQPWh9YBIjM6j5ls8nkAAi6c+iGOgP2XSLuXm5uL48eNQqf77cMpDKuKkhy1jDDVJn0PfWAefhPcgEFm+CEggEMDDhRYTOQMKXdIubdu2DU899RQ6duyIIUOGYPny5TDcKOXkI/yNo19AW1OKbk+/D6HEhYNqAbXOAFl3T07GIvaN/mkl7QpjDNeuXYNQKLx9NLpGgzNnzuDMmTPIzc2FNOYVqC0IXt3NStzK/hkQSXBt4wzj1zs/9Qo6hD/e5nFdxUJ083Jt8+uJ46DQJQ6rJWAzMjKQkZGB9PR0ZGRkQCAQICQkxDiV4OrqinXr1uHVV1/FM9v+g3NFtW2+prhjN/Recoirt2Ak6+HF+ZjEPlHoEofAGENZWZlZwAJATEwMYmNjMX/+fMTExMDPzw8GgwEeHh7o3LkzDh8+jMGDBwMA4mW+yC27aVdbgV0lQjol2InQwZTE7jDGcP36dbOANRgMxoCNiYlBTEwM/P397/lwbN++fRg+fDi6detm/FptowaD1x63aIqBay5iIc4seQKdPKR8l0JsgEKX8K61gNXpdGYBGxAQwMnqg1e+z8SR/HK7aO8oFACjB/TApmej+S6F2AiFLrGp8vJyk3DNyMiARqMxCdjY2FjOArY12aV1mLb9DG/no93JTSLC7rmDERngzXcpxEYodInVVFRUmAWsSqUyuXuNjY1Fr169rBaw9/Lyjgwcv1TJ65E9UpEQT8i6Ycv0GN5qILZHoUs4oVQqTQI2PT0dzc3NZgHbu3dvmwdsa+qaNBj+yUnUq/jrd9DRTYy0t/+Cjm4S3mogtkehSx5aS8DeOQ/b1NRkEq4xMTHo06ePXQTsvRy7qMSi3Vm8TDOIoceWGYMwKqy7za9N+EWhS+6rsrLSLGBv3bplFrB9+/a164AFAJVKhbNnz+LEiRPYu3cvCgoK8M6uM/gxv8amwesqFqL2/AHUJG/H+PHjMX78eIwYMcLu/5Ei3KDQJUZVVVVmqwgaGhoQHR1tErCBgYEOFw41NTXw8/ODRCJBc3Mz9Ho9evfujd9//x1L9ufhQM51mwSvm0SE8fKeGMR+w9NPTwIAuLu7Q6PRYMGCBfjss8+sXgPhl12ErrJeBUV5PRo1emh1BkjEQnhIRQjr4UVbI62kurraLGBv3ryJ6Ohok3nYoKAghwvYe1mwYAG+/PJLaLVauLu748svv8S0adPAGMPqJAV2ni2xavC6SUSYHtcLy0bLAAB9+vRBSUnJ7Z9zc0NOTg769+9vtesT+8BL6NY2apCYUYpkhRKKigZodAZIxUIwxsDY7TZ3AoHA+HVZd0/Ey3wxOSaAFpC3QU1NjVnA1tbWthqwQmH77IGkVqsxZ84cHDhwACqVCh4eHqisrIRE8t+HWMcuKrE4MRsqrYHTVQ1SkRBuUiHWT5ab7DzbsmUL3nzzTWi1WhgMBnz66ad47bXXOLsusU82Dd3s0jpsP1WIZIUSAgEeaiumq0QIxm5v45w7LBByWtfYqhs3bpgF7I0bNxAVFWUSsP369Wu3AXu36upqTJw4EV27dsWWLVswYsQIPPvss3j//ffNfm1dkwZL9+chpaAKap3eog0UQgHgIhZhZIgP1iREwNvd9IahsbER/v7+WL9+Paqrq7FkyRK88MIL+PLLL9t+UWL3bBK6VvlGDvbBmonm38jOpLa21mwdbHV1tUnAxsbGOlXA3u3y5csYM2YMnn76afztb3+DUCiEVquFSCS67+9Jzh83CMcsuEEY9ccNwv02Pmg0Gkilt7+HDx06hISEBAwaNAipqakQi6k1Sntk9dC15kc2V4kQG6bInaJZSG1tLTIzM00CtqqqCnK53CRg+/fv77QBe7eTJ09i6tSpWLNmDV544YU2jVHbqMG+zGu3p8LK66HSGeDSylSYWmeAq1gIWQ8vxMt88XS0f5umwhQKBeLi4uDl5YXc3Fx07ty5TXUT+2W10OXj4UR7eeBTV1dnErDp6emorKxEVFSUyVKt4OBgCth7+Prrr7F06VLs3r0bjz/e9j63d6usV0FR0YBGtc74zMHDRQxZd0/OHvrW1dVBLpejuroap0+fxsCBAzkZl9gHq4QuY8z2y3Aie2LtxAiHC96WgL1zHlapVEIul5sFrEgk4rtcu2cwGPDee+9h3759OHz4MEJCQvguqU0MBgPi4+ORlpaGPXv2YNKkSXyXRDhildBddfgivrfyHe7dWu54l48Js9k1H9bNmzfNAraiogKRkZEm62BDQkIoYNugqakJM2bMQFVVFfbv34+uXbvyXZLFXnnlFWzZsgUffvghVqxYwXc5hAOchy6fWyvdJCJsnBZlF3O89fX1xoBtmSa4fv26WcCGhoZSwHKgvLwc48ePh0wmw/bt2+Hiws3ZZfZgy5YtWLhwISZNmoS9e/fyXQ6xEKehay9NRFLfetymqxrq6+uRlZVlErBlZWUYOHCgyTKt0NBQeiJtBTk5ORg3bhzmzZuHZcuWOdwU04NISUnBk08+idDQUJw9exaurrRpyFFxGrrO0C6voaHBLGBLS0vNAlYmk1HA2sDhw4cxa9YsbNq0Cc888wzf5VhVcXExoqOjIRKJkJ2djZ49e/JdEmkDzkK3PTaGvnXrFrKyskyWaZWUlCAiIsIkYMPCwihgbYwxho0bN2Lt2rXYv3+/8Qy09q6pqQmxsbH4/fffceLECQwZMoTvkshD4ix0Hf0IlFu3biE7O9skYIuLixEREWGyikAmk5lsHSW2p9Pp8NprryE1NRWHDh1Cnz59+C7JpgwGAxISEnDo0CF8/fXXmDlzJt8lkYfAye1ZbaMGyQqlXQQuABgYcEyhRG2jptUF6o2NjWYBW1RUhAEDBiAmJgYjR47EW2+9hbCwMApYO1NfX49nnnkGBoMBv/76Kzp27Mh3STYnFArx008/YdmyZZg9ezby8/PxySef8F0WeUCc3OluS7uKDckFdnes9eJRIZge090kYNPT0/H7778bA7blDjY8PJwC1s4VFxdj7NixGDZsGD7//HOa0gGwc+dOzJw5E/Hx8UhKSqLNMg6Ak9Cd8o/TOFdUa9EY1Qf/H1RFOTBoVRB5dILX4EnwjHzSojGF1VdRsXMJwsPDTZZphYeHG/e7E8dw5swZTJw4Ee+++y4WLVrULlcotFV6ejqGDx8Of39/ZGZmokOHDnyXRO6Dk1sFRUWDxWN4DZ6MLv/3NQjEEmhrSlHx/VJIfYPg0r1fm8d069EPdXV1FLAObs+ePVi4cCG+/vprjBs3ju9y7E5sbCwKCwsRFRUFf39/ZGRkICgoiO+yyD1Y/FlEWa+CRmf5tILUpzcE4paP9wIIIICuttyiMbVMgDqV/Ux5kIfDGMPq1avx9ttv49ixYxS499G9e3cUFxejX79+kMlkOHbsGN8lkXuw+E5XUV4PqVgINQfBW3N0MxrzjoPp1JD6BsEtKNai8VzEQigqGuj0CQekVqsxb948XLhwAWfOnKE1qQ9AKpUiPT0dzz33HJ566ils2LCBmqLbIYtDt1GjB1f7K7o8uQCdR70EddklqEryIBBZ9mCLMYZGNX+740jb1NTUICEhAV27dkVqaio8PDz4Lsmh7NixAxEREXjjjTeQn5+P7du3810SuYPF0wtanQFcdm8QCEVwDQiHvqEaDVlJFo3FGDiZ+iC2c/nyZQwePBhDhgzBvn37KHDb6N1338VPP/2Eb7/9Fo899hh0Orr5sBcWh65ELIRVHiQbDBbP6QoEgFRs+hbVarVFYxLrSUlJwfDhw7FkyRKsW7eOlj9ZaNy4ccjNzUVubi769u2LGzdu8F0SAQeh6yEVWbx8R99Yh8aLqTBomsEMejQXZqBRkQrX3pEWjSsQCCAVMpw4cQKLFy9Gv3794ObmhurqaovGJdz75ptv8Mwzz2DXrl2YM2cO3+W0GzKZDCUlJRCJROjVqxfy8vL4LsnpWTynK+vhZflHeIEADVlHUHN0M8AMEHfshk5PzIV7sGX76dVaPZ4aHAHWVPfH8SoMrq6u6NKli2X1Es7c2XQ8NTUVoaGhfJfU7nh7e6OwsBBPPPEEoqKisHfvXkycOJHvspwWJ5sjIj48igYe2znei5erGC92uox33nnHOK0gEong6+trshstJiYG3bt357la59PSdLyyshL/+te/2kXTcXu3YMECbN26FR999BGWL1/OdzlOiZNJM1l3Ty6G4ZyshxcWLVqEjIwM+Pn5QSgUYuHChfj111/x/PPPQ6VS4bPPPkNYWBj8/PwwYcIEfPTRR0hKSoJSqeS7/HatvLwcI0aMgLu7O5KTkylwbWTz5s344osvsHLlSkyZMoXvcpxSu++9MHdYIIDbjVLmz5+P+fPnY/jw4Sa/ljGGoqIis+PM3d3dTU7ajYmJQbdu3fh4O+1KTk4Oxo8fj7lz57bbpuP2rqUpukwmw7lz52jXpg1xErq1jRoMXnuckw0SXHERC3FmyRNtOgYbuB3Ev//+u1kQd+jQwaSXLgXxwzl8+DBmz56NjRs3tvum4/aupSm6WCxGVlYWbUCxEeqn+xAYYygsLDQJ4szMTHh6epoFsY+PD6fXbg82btyINWvW4IcffqDm23aiqakJMTExKCoqoqboNkInR1iIMYarV6+aBXHHjh1NpiViYmKcdt5Sp9Ph9ddfx8mTJ3Ho0CH07duX75LIHagpum3RGWlWYDAYjEHcEsaZmZno1KmTWRC39+VrdzYd37t3r1M2HXcU7733HtauXYvFixdTU3QrapenAae9/Rd0dLOvhuQGgwG//fabSRBnZWWhc+fOZkHcuXNnvsvlREvT8cceewwbN26kpuMOoKUp+qhRo3D48GHaFWgFnIYuABy7qMSi3Vm8TDO4SUTYOC0K8TJfm1+7LQwGA65cuWIWxF27djVZRxwdHe1wQXz27FkkJCTgnXfewWuvvUYrFBwINUW3Ls5DFwD+evgidp4tsWnwuklEmB7XC8vHhNnsmtbQEsR3Hi+UlZWFbt26mQVxp06d+C63VXv37sUrr7xCTccdWEVFBaKioqBSqZCenk5N0TlkldBljGHJ/jwcyLluk+B1k4gwXt4TaxMi2uUdlV6vNwnijIwMYxDfuY44Ojoa3t7Wf3h4L4wxrFmzBlu3bsWBAwcgl8t5q4VYTqPRYOjQocjNzcXhw4cxatQovktqF6wSusAfXf+TFFa/4225w102WtYuA/de9Ho9CgoKTII4Ozsbvr6+ZkFsi4dXarUaL730EvLz83HgwAFa89mOPPfcc9i1axf+/ve/49VXX+W7HIdntdBtceyiEosTs6HSGjhd1SAVCeEmFWL9ZLnDzOFam16vx+XLl02Wr2VnZ6NHjx4mD+q4DuKamhpMnDgRnTt3xo4dO6gHbju0du1avPfee3jxxRexbds2vstxaFYPXeD2qoal+/OQUlAFtU5v0QYKoQBwEYswMsQHaxIi4O1O2xfvR6/X49KlSyZBnJOTg549e5oFsZeX10OPX1BQgDFjxmDixIlYs2YNPe1uxw4cOIBJkyYhLi4OKSkptBqljWwSui1ySuuw/VQhjimUEAjwUL0aXCVCMAaMkvli7rBAm2x8aK90Op1ZEOfm5sLPz89k+VpUVNR9gzglJQXPPPMMVq9ejRdffNGG74Dw5cKFCxg8eDA6deqE7Oxsh1tVYw9sGrotahs12Jd5DckKJRTl9VDpDHARC//oeXv7xAeBQAC1zgBXsRCyHl6Il/ni6Wj/NvdSIPen0+mgUChMlq/l5uYiICDALIg9PT3xzTff4N1338WuXbvwxBNP8F0+saG6ujrI5XJUV1fjzJkzGDBgAN8lORReQvdulfUqKCoa0KjWQaMzQCoWwsNFDFl3TzrJl0c6nQ4XL140C2I3NzdoNBosXLgQo0ePRlRUFK3ldDIGgwF/+ctf8O9//xuJiYlISEjguySHYRehSxxDU1MTnnvuORQVFWHmzJnG1RP5+fno3bu3yTpiuVxOQewEqCn6w6OZcPJAysvLMWHCBISEhOA///kPXFxcjD+n1Wpx8eJF4/zwrl27kJ+fjz59+pgFMa1saF82b96MAQMG4NVXX0V+fj52797Nd0l2j+50yZ/Kzc3FuHHj8OKLL2L58uUPtB5aq9XiwoULJuuI8/PzERgYaBbE7u7uNngXxJpOnjyJp556CmFhYTh79iw1Rb8PCl1yX0lJSZg1axY+//xzTJ061aKxNBqNWRBfuHABgYGBJhs6IiMjKYgdUHFxMaKioiCRSKgp+n1Q6JJ72rhxI/72t79h//79VmturdFokJ+fb7J87eLFiwgKCjJZR0xB7BjubIqekpKCuLg4vkuyOxS6xAzfTcfVarVZECsUCvTr188siN3c3GxaG/lzBoMBEyZMQFJSEr755hs8//zzfJdkVyh0iYmWpuN6vR6JiYl203RcrVYjLy/PJIgvXbqE/v37m6wjjoyMhKsrLTO0B0uXLsW6devw1ltv4eOPP+a7HLtBoUuM7mw6/vnnn0Misa9G8HdTqVRmQXz58mUEBwebBPHAgQMpiHlCTdHNUegSAO2n6bhKpUJubq7Jho6CggKEhISYBfGdy96I9Zw/fx4jRoygpuh/oNAlSExMxIIFC9pt0/Hm5mazIL5y5QpCQ0NNlq9FRERQEFtJRUUF5HI51Go1MjMznfpwUgpdJ8YYw9q1a7FlyxanazreEsR3Ll+7cuUKZDKZSRAPGDCAgpgjGo0GQ4YMQV5eHo4cOeK0PTsodJ2URqPBvHnzkJeXh4MHD9KaStwO4pycHJMg/u233xAWFmYWxLT4v+2mT5+O3bt347PPPsPChQv5LsfmKHSdEDUdf3BNTU1mQXz16lWEhYWZLF+jIH44LU3R582bh61bt/Jdjk1R6DqZlqbjCQkJWLt2LT1NboPGxkbk5OSYrJooLCxEeHi4WRDb+woQPrU0RR88eDBOnjzpNE3RKXSdSEvT8b/+9a+YO3cu3+W0K42NjcjOzjYJ4qKiIrMgDg8PpyC+g7WaoivrVVCU16NRo4dWZ0zIHcQAAB2ISURBVIBELISHVISwHl68t4ul0HUS33zzDZYsWYLvv//eaR9g2NqtW7fMgri4uBgDBgwwWb4WFhbm1EFcV1eHyMhI3LhxA//5z3/a1BS9tlGDxIzS2wcjVDQY+3LffTBCy9dl3T0RL/PF5JgAmx+MQKHbzhkMBixfvhx79+7FoUOHEBoayndJTu3WrVvIysoyWb5WUlKCiIgIsyB2lo/bwO3v08cffxynT59GYmIiRowYgTlz5mDbtm3o2rXrPV+X/ccRYMkWHAEW/8cRYHIbHQFGoduONTU1YebMmSgvL8ePP/54329ewp+GhgazIC4tLcXAgQNNVk3IZLJ2H8Qvv/wytm7dCn9/f5SXl+ODDz5otTm6VQ67DfbBmonWP+yWQredqqiowPjx4xEcHIyvvvqK1po6mJYgbpmWSE9PR1lZmVkQh4aGtqsgZowhLi4O58+fBwB06dIFFRUVJu/x2EUlFidmQ6U1QKN/8DvbPyMVCeEqEWLDFDniZb6cjXs3Ct12KC8vD+PGjcMLL7yAFStWOOyWXmKqvr7eLIivX7+OyMhIsyAWiUR8l9sme/bsMenbLBQK8f333+OZZ54BYwyrkxTYebYEzVq91Wpwk4gwPa4Xlo2WWeXvDoVuO9PSdPyzzz7DtGnT+C6HWNnNmzdNgjgjI8MYxHeumnCUIL558yZ++OEHJCUl4cSJE6itrYWLiwsqKyux+lgRDuRct2rgtnCTiDA+sifWTozgPHgpdNuRlqbjP/zwA4YOHcp3OYQndXV1ZkFcUVFhFsQhISF2H8SXLl3CokWLUNwpBoLgYVDpbBdXLXe8y8eEcTouhW47oNPp8MYbb+D48eM4fPiwUzcTIa2rq6tDZmamyfI1pVIJuVxuEsTBwcF2F8S/XKzAgv89Dx1sv5HHTSLCxmlRnM7xUug6uPr6ekydOhU6nQ579+6Ft7dtlr0Qx1dbW2sWxFVVVZDL5SbL14KDg3nbuVjXpMHwT06iXqXj5foA0NFNjNS3HudsVQOFrgMrKSnB2LFjMXToUGzcuNGpF9gTbty4ccMsiKurqxEVFWUSxP3797dJEL+8IwPHL1VyukrhYUlFQjwh64Yt02M4GY9C10GdO3cOCQkJeOutt/D666/TCgViNTU1NcYgbgnjmpoaREdHmwRxv379OA3i7NI6TNt+xiYPzv6Mm0SE3XMHI5KDDRQUug5o3759ePnll/HVV19h/PjxfJdDnFBNTY0xhFuCuLa21hjELWEcFBTU5iB+5ftMHMkvt2jjA1eEAmD0gB7Y9Gy0xWNR6DqQlqbjmzdvxoEDBxAVFcV3SYQYVVdXIzMz02Qd8c2bN1sN4j/7ZFbbqMHgtceh1vE3rXA3F7EQZ5Y8YXGvBgpdB6HRaPDSSy8hJycHBw8ehJ+fH98lEfKnqqqqTII4IyPDGMR3rpq4O4i3pV3FhuSCh+qlYG2uEiEWjwrB3GGBFo1DoesAampqMGnSJHh7e2Pnzp3UdJw4tKqqKpMHdRkZGWhoaDAJ4v8t74q8iuY2X6M+4yAa845DU1UED9kIdB37Bie1x/XtjD3zhlg0BoWunbty5QrGjBmDCRMmYO3atXa3hpIQLlRWVprMD2cFTYdA6t7m8ZounwYEAjT/ngmm1XAWul6uYuSufNKiMSh07VhqaiqmTJmCVatWYd68eXyXQ4hNKOtVGP7JSU7mc2vT/hf6+mrOQlcqFuLfbz9uUSN0OqvFTn377beYPHkydu7cSYFLnIqivB5SsX1Gk4tYCEVFg0VjtJ+ecO2EwWDAihUrsHv3bqSmpkImk/FdEiE21ajRw14/gDPG0Ki2bHccha4daW5uxvPPP4/y8nKcOXMGPj4+fJdEiM1pdQbYaeaCMUBj4bSHfd7DO6GKigqMHDkSUqkUycnJFLjEaUnEQtjrBkuBABZPfVDo2oG8vDwMHjwYY8aMwY4dO+Dqyu9ppYTwyUMqsnhbOzPowXQawKAHmAFMpwEzWL6dWCAQwMPFsgkCml7g2ZEjRzBz5kz8/e9/x7PPPst3OYTwgjGGoqIiZGRk4FR6Lm6xGEDY9ni6+etu3Px1l/HHjRdOouOj0+A9bLpFdap1Bsi6e1o0BoUujzZt2oTVq1fjX//6Fx599FG+yyHEJhhjKC4uNtkckZGRAVdXV+PmCFeNGM0W3Jh6D5tuccC2xlUstGi5GEChywudToc333wTycnJ+PXXXxEYaNm2QkLsFWMMJSUlJv0YMjIy4OLiYgzY1157DTExMejevbvxdfn/OI1zRbU8Vt46WQ8vi8eg0LWxhoYGTJ06FRqNBqdPn6am46TdaAnYO3eWZWRkQCKRGAN20aJFiImJQY8ePe47VrzMF7llN+2u9wIXJ0jQjjQbamk6PmTIEGzatImajhOHxRhDaWmpWcCKRCJjwLb8v2fPng89fnvuMkZ3ujZy/vx5/M///A8WL16MN954g5qOE4fBGMO1a9fMAlYoFBrDdcGCBcaA5eJ7u5OHFPEyX7vqpztK5mtx4AIUujbxww8/YP78+dR0nNg9xhjKysrMuoABMAbs/PnzERMTAz8/P6vePMwdFogTlyrt4uQIF7HI4paOLSh0rYgxhnXr1uGLL77A0aNHER1tedd5QrjCGMP169dNAjY9PR2MMWPAzps3DzExMfD397f5pzN5gDdGBvvYxRlpI0N8ODmqB6A5XavRaDSYP38+srOzqek4sQutBaxerzeZf42NjeUlYO/FXk4DTnv7L+joxs0zGLrTtYIbN25g0qRJ8PLyQlpaGjp06MB3ScTJlJeXmwWsVqs1husLL7yAzZs3IyAgwG4CtjXuYiBGewEntH0gkLjY/PpuEhHWT5ZzFrgA3ely7sqVKxg7dizGjRuHdevWUdNxYnUVFRUm86/p6elQq9Vmqwh69epl1wEL3F7Dfu7cOZw8eRIHDhzA+fPnIRQK8eaOf+PgxTqbzu+6SUSYHtcLy8eEcTouhS6HqOk4sTalUmkWsM3NzWYB27t3b7sP2NZ89913mDVrFiQSCbRaLQBg3bp1ePvtt7Fkfx4O5Fy3SfC6SUQYL++JtQkRnP8+Uuhy5LvvvsPbb7+NnTt3YtSoUXyXQ9oBpVJptkyrqanJ5GTdmJgY9OnTxyEDtjVarRbR0dHIz88HAHh4eKCyshLu7u5gjGF1kgI7z5ZYNXhb7nCXjZZZ5feV5nQtZDAY8P7772PXrl1ISUlBWBi3H0WIc7j7jLCMjAzcunXLGK7PPfccPv30U/Tt27fdBGxrvv76a1y4cAHu7u5Qq9V48cUX4e5++6w0gUCA5WPCENe3CxYnZkOlNXC6qkEqEsJNKsT6yXJOdp7dC93pWqC5uRkzZ85EWVkZfvzxR+qBSx5Iy2m4dwbs3afhxsTEIDAwsF0H7N0WLlyIzZs344MPPsCMGTMwevRoHDlyBH369DH7tXVNGizdn4eUgiqodXqLNlAIBbfX4Y4M8cGahAh4u1u+AeJ+KHTbSKlUYsKECQgKCsJXX31FPXBJq6qrq802Gty8edMsYIOCgpwqYO9kMBjwxBNP4NSpU9izZw8mTZr0wK/NKa3D9lOFOKZQQiDAQ/VqcJUIwdjtnWZzhwVytg73z1DotkF+fj7Gjh2L2bNn4/3333favyzEVE1NjVnA1tbWmgRsbGwsAgMDIRTS+QEAUFdXh6ioKFRVVeH06dMYOHBgm8apbdRgX+Y1JCuUUJTXQ6UzwEUsBGMMjN0+8UEgEECtM8BVLISshxfiZb54Otqfk629D4NC9yH9/PPPeP755/Hpp59i+nTu+3USx3Djxg2zdbAtAXvng66goCAK2HtQKBSIi4tDx44dkZOTg86dO3M2dmW9CoqKBjSqddDoDJCKhfBwEUPW3dPifriWotB9CJs3b8aqVauwb98+ajruRG7cuIHMzEyTgK2pqTEL2H79+lHAPqBDhw4hISEBcXFxSElJgVjsPM/0KXQfgF6vx5tvvolffvkFhw8fpqbj7Vhtba3JaQbp6emorq5GVFSUyTKt/v37U8C20ccff4wlS5Zgzpw52L59O9/l2ByF7p9oaGjAtGnToFKpsG/fPmo63o7U1tYiMzPTZJqgqqoKcrncJGCDg4MpYDny/PPPY+fOndiwYQNee+01vsvhBYXufZSUlGDcuHGIi4vDF198QU3HHVhdXZ1ZwCqVSsjlcpNVBMHBwbR12wo0Gg0effRR5OTk4PDhw069gYhC9x5amo6/+eabePPNN2mFggO5efOmWcBWVFQgMjLSJGBDQkIoYG2gsrISkZGRaG5uRkZGBoKCgvguiVcUuq3Yv38/XnrpJXz55ZeYMGEC3+WQ+6ivrzcL2OvXr5sFbGhoKAUsD9LT0zF8+HD4+fkhKyuLOu6BQtcEYwwff/wxNm3ahJ9++omajtuZ+vp6ZGVlmayDLSsrw8CBA03WwVLA2oddu3ZhxowZiI+PR1JSEs2L/4FC9w8ajQYvv/wysrKyqOm4HWhoaDAL2NLSUkRGRpos0woNDXWq5UaOYtmyZVizZg3eeOMNrF+/nu9y7AqFLm4/xZ40aRI8PT2xc+dO+ghkY7du3TIJ2PT0dJSWlmLgwIEmASuTyShg7ZzBYEBCQgIOHTqEL7/8ErNnz+a7JLvj9KH722+/YcyYMRg7diw+/vhj+lhqZS0Be+c62JKSEkRERJgFLK0WcSxNTU2IjY3F77//jhMnTmDIkCF8l2SXnDp009LSMGXKFHz44Yd46aWX+C6n3WlsbDQL2OLiYgwYMMBkHWxYWBgFrIMrLi5GdHQ0RCIRMjMz4e/vz3dJdstpQ/ef//wn3nrrLWo6zpHGxkZkZ2ebrCIoKipCeHi4ScCGh4dTwLYzaWlpGDVqFEJDQ3H27FnquPcnnC50W5qOf//99zh06BA1HW+DpqYms4AtLCxEeHi4yTKt8PBwSKW27eBEbOsf//gHFixYgIkTJyIxMZHvchyCU4Vuc3MzZs2ahdLSUvz444/o1q0b3yXZvaamJuTk5JgE7NWrVxEWFmYSsAMGDKCAdTKvvvoqvvjiC6xcuRIrV67kuxyH4TSh29J0PDAwEF9//TV9BGpFc3OzWcD+9ttvkMlkZgHr4mL747CJfbiz6fju3bvx9NNP812SQ7GL0FXWq6Aor0ejRg+tzgCJWAgPqQhhPbw46X2Zn5+PcePGYebMmVi5ciVt6cXtgM3NzTVZB3vlyhWEhoaabDSggCV34qrpuDPjJXRrGzVIzCi93eW9osHYZPjuLu8tX5d190S8zBeTYwIeusv70aNHMWPGDKduOq5SqcwCtqCgAKGhoSbLtCIiIihgyT3d2XQ8OzsbXbp04bskh2TT0M3+4zyjZAvOM4r/4zwj+QOcZ7RlyxZ89NFHSExMxGOPPWZJ6Q5DpVIhLy/PZKNBQUEBQkJCzAKWpljIg0pKSsKECRPwyCOPIC0tjTapWMAmoWuVkzuDfbBmYusnd+r1eixevBhHjx7FoUOH2m1XI7VajdzcXJN1sJcvX0ZwcLBJwA4cOJAClrTZJ598gnfffRezZ8/GV199xXc5Ds/qoXvsotJqZ9S7SoTYMMX0jPqGhgY8++yzaGpqwr59+9CpUyfOrskntVqNvLw8k4C9dOkS+vfvb7IOduDAgXBzc+O7XNJOPP/889ixYwc+/fRTp206zjWrhS5jDKuTFNh5tgTNWr01LgEAcJOIMD2uF5aNluHatWsYO3YsBg0ahM2bNzvsInyNRmMWsAqFAv369TMJ2MjISApYYhVarRaPPvoosrOzcfDgQTz55JN8l9RuWCV0GWNYsj8PB3KuWzVwW7hJRBjiJ8Evq2bi9ddfx+LFix1mhYJGo0F+fr7JMq2LFy8iKCjILGDd3d35Lpc4gTubjp8/fx79+/fnu6R2xSqhu+rwRXxv5TvcuzGtCo8HSPDta+Ntds2HpdFocOHCBbOA7du3r8k6WLlcTgFLeJGRkYFhw4ahZ8+eyM7Opo57VsB56B67qMSi3Vk2DdwWbhIRNk6LMpnj5YtWq8WFCxdMlmnl5+e3GrAeHh58l0sINR23EU5Dt65Jg+GfnES9SsfVkA+to5sYqW893uqqBmvRarW4ePGiWcD27t3bZKMBBSyxVytWrMDq1aup6bgNcBq6L+/IwPFLlZyuUnhYUpEQT8i6Ycv0GKuMr9PpzAI2Ly8PvXv3NlmmJZfL6aMZsXsGgwETJ07EwYMHsX37drzwwgt8l9TucbbCObu0DikFVbwGLgBo9AakXK5CTmkdIh9gA8X96HQ6KBQKk40GeXl56NWrlzFgp06dCrlcDk9PT47eASG20dTUhEceeQSFhYVIS0vDo48+yndJToGz0N1+qhBqne3ncVuj1umx/VQhNj3734Ml9Xo91q1bh3HjxiEiIsLsNS0Be+cyrby8PPj7+xvvXqdMmYKoqCgKWOLwiouLERMTA6FQiCtXrlDTcRviJHRrGzVIVigt2mnGJQMDjimUqG3UoJOHFFVVVRg/fjzOnTuHW7duISwszCxgc3Nz4efnZwzYp59+GlFRUfDy8uL77RDCqZam4yEhITh37hztVrQxTuZ0t6VdxYbkgofqpWBtrhIhFo8KwQCxEmPGjEF9fT30ej08PT1hMBjQs2dPk3WwUVFR6NixI99lE2JV27Ztw8svv4yEhATs27eP73KcEid3uskKJWeBq71RhutfLYRH6KPoOu6tNo+j0hpwOLsE8xY9bvJ1vV6PsrIyeHtbNt9LiKOhpuP2gZPQVVQ0cDEMAODGL1vh0oObHTCFN9S4cOECTp48if379+P06dNoamqCTsffkjZCbM1gMCA+Ph5paWnYu3cvNR3nmcWhq6xXQaPj5i638WIqhK4ekHQJha6u3OLxVDoDuvoH4pVXwvDKK69Aq9Xi0qVL1AeUOI36+nrI5XJUVlYiIyMDkZGRfJfk9CzecqIor4dUbPnOFYO6CXWndqLTX+ZYPFYLF7HQ5C5cIpEgIiLCYfoyEGKJy5cvIyAgABqNBkVFRRS4dsLitGzU6MHF/oq6tP9Fh8j/A7GXj8VjtWCMoVFNUwnE+SQlJWHAgAEIDw9HUVERunbtyndJ5A8Wh65WZ4ClmatRFkJVnAOvRyZYWo4JxsDZ1AchjuKTTz7B2LFjMWPGDJw+fZpOebAzFv9pSMRCWPppXVWSB91NJa5tng0AYBoVwAwor34NPWZ/1uZxBQJwMvVBiKNoaTq+YcMGvP7663yXQ1phceh6SEUWz5F2kD8JD9lw44/rz+2H7qYSnZ98xaJxBQIBPFzoX3nS/ul0OgwdOhTZ2dk4cuQINR23YxYnkqyHl8Uf4YUSV0Dy310xAokrBGIpRO6WbVZQ6wyQdactu6R9q6yshFwuR1NTEy5cuEBNx+2cxaHr6+UKqVgINYdzp97DuDkq3VUsRDcv2uJI2q/MzEw89thj6NmzJy5fvkx9QRwAJxOe9no3KetBfRNI+7Vnzx4MGjQIw4YNQ0FBAQWug+AkdONlvnCV2NcDK1eJ0C5OkCDEGlasWIFp06Zh0aJFOHr0KJ3y4EA4+ZOaHBNg8bIxrjEGPB1N7epI+2IwGJCQkIC//e1v2L59OzZs2MB3SeQhcfJov5OHFPEyXxzJL7eL9o5CATBK5otOHrY7socQa2tubkZsbCyuXr1KTccdGGefSeYOC4SLWMTVcBZxEYswd1gg32UQwpmSkhIEBASgsrISv/32GwWuA+MsdOUB3hgZ7AOpiN+5JalIiJEhPhYf1UOIvUhLS0NwcDB69uyJ0tJSOuXBwXGakGsmRvD+QM1NKsTaiQN5rYEQrmzbtg2PP/44xo4di9zcXDrloR3gNCG93aVYP1kONwk/0wxuEhHWT5ajo5uEl+sTwqVFixZh/vz5WLFiBZ3y0I5wvkd2VJgvpsf1ws6zJWjW2u6gSjeJCNPjetEyMeLwDAYDRo0ahdTUVGo63g5ZpTHBstEyNKh0OJBz3SbB6yYRYby8J5aNlln9WsS5KetVUJTXo1Gjh1ZngEQshIdUhLAeXpzsfryz6Xh6ejrkcjkHVRN7YpXQFQgEWDsxAp6uYqvf8bbc4S4bLaPm5IRztY0aJGaUIlmhhKKiARqdAVKxEIwxMHa7k51AIDB+XdbdE/EyX0yOCXjoJYuXL1/GoEGD4OnpST1w2zFOTgO+n2MXlVicmA2V1gCNnrv+DFKREG5SIdZPltOUAuFcdmkdtp8qRLJCCYEAD3XwqqtECMZu79ScOywQ8nuspNFqtYiLi8Mnn3wCtVqNCRMmIDY2FqdOnaIeuO2Y1UMXAOqaNFi6Pw8pBVVQ6/QWbaAQCm6vwx0Z4oM1CRHwdqcNEIQ7VvleDfbBmonm36uJiYmYMWMGGGPQarWYNWsWvv76awvfAbF3NgndFjl/3D0cs+DuYdQfdw+0DpdwzZqfylwlQmyYYvqpLCoqCtnZ2QAAT09PFBYW0pSCE7Bp6LaobdRgX+a12/Nk5fVQ6QxwaWWeTK0zwFUshKyHF+Jlvng62p+29hLOMcawOklh0+cPGRkZGDRoEBhjEAgEEIlEePHFF7FlyxarXZ/YB15C926V9SooKhrQqNYZH0h4uIgh6+5J/XCJVTHGsGR/nk1X2vzfcB/885UnUVNTg6CgIEyYMAHx8fEYNmwYOnToYPUaCL/sInQJ4cuqwxfxvY3XlEOnQfdbv+HAhzPRrVs3212X2AVqwkmc1rGLStsHLgCIpbjpE4HcGrrfcUYUusQp1TVpsDgx2/aB+4dmrR6LE7NR16Th5fqEPxS6xCkt3Z/3UKtnrKFZY8DSf+XxWgOxPQpd4nSyS+uQUlDF6bKwttDoDUi5XIWc0jpe6yC2RaFLnM72U4VQ6/iZVribWqfH9lOFfJdBbIhClziV2kYNkhVKuzhWCgAMDDimUKK2keZ2nQWFLnEqiRmlsLe+SAIBsC/zGt9lEBuhrhrEqSQrlBY/QKvYuQTq65chEN5u1i/y7AK/ef9o83gqrQHJCiWd6+ckKHSJU1FUNHAyTuf/Mx+ekU9yMhYAKMrrORuL2DeaXiBOQ1mvgkbH74qFe1HpDKisV/FdBrEButMlTkNRXg+pWAg1B8Fbl/Id6lK+g6SzH7yHz4Brb8sOQ3URC6GoaKBeI06AQpc4jUaNHly0Gun0+GxIugRAIJKgUZGGyh9WocfszyHp1KPNYzLG0KjWWVwbsX80vUCchlZnABftnVx6hkDo4g6BWIIOEU/AxU+G5qvpFo3JGOx26oNwi0KXOA2JWGid5WICAQDL0lwgAKRi+uvoDOhPmTgND6nI4sNLDapbaC7MANNpwAx63LpwEurSfLj1jbZoXIFAAA8Xmu1zBvSnTJyGrIeXxR/hmUGPurQd0N64BgiEkHTxh8/E5ZB08bdoXLXOAFl3T4vGII6BQpc4DV8vV4tXL4jcO6LHrE85rOo2V7GQVi44CZpeIE7FXu8mZT28+C6B2AiFLnEq8TJfuErs69veVSI0OSWYtG/29d1HiJVNjgngZNkYlxgDno62bE6YOA4KXeJUOnlIES/zhdBOOo0JBcAomS86eUj5LoXYCIUucTpzhwXCRSziuwwAgItYRN3FnAyFLnE68gBvjAz2gVTE77e/VCTEyBAfRAZ481oHsS0KXeKU1kyM4P2BmptUiLUTLWuUQxwPhS5xSt7uUqyfLIebhJ9pBjeJCOsny9HRTcLL9Ql/KHSJ0xoV5ovpcb1sHrxuEhGmx/WiZWJOikKXOLVlo2UYH9nTZsHrJhFhvLwnlo2W2eR6xP4IGBcNRglxYIwxrE5SYOfZEjRrrXc0e8sd7rLRMosb7xDHRaFLyB+OXVRicWI2VFoDNHruettKRUK4SYVYP1lOUwqEQpeQO9U1abB0fx5SCqqg1ulhsOBvh1Bwex3uyBAfrEmIgLc7bYAgFLqEtCqntA7bTxXimEIJgQAPdWy7q0QIxm7vNJs7LJDW4RITFLqE3Edtowb7Mq8hWaGEorweKp0BLmIhGGNg7PaJDwKBAGqdAa5iIWQ9vBAv88XT0f60tZe0ikKXEEJsiJaMEUKIDVHoEkKIDVHoEkKIDVHoEkKIDVHoEkKIDf1/GPrB+X0GBu8AAAAASUVORK5CYII=\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(outputfolder+\"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",
   "metadata": {},
   "source": [
    "# Part 1: Preliminaries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CRRA utility function\n",
    "def CRRA(x,gamma): #can use offset by xMin in case of negative values of x\n",
    "    if gamma == 1:\n",
    "        return np.log(x)\n",
    "    else:\n",
    "        return ((x)**(1-gamma) - 1)/(1-gamma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the payoff, assuming that the utility depends solely on the payoff\n",
    "# Estimate the prob of decising at turn 0, each subject has 11 strategies, at each node (A->E), wait or decide at turn 0\n",
    "# The difference (compared to network (a)) here is that when a subject is placed at A, his strategies are \"decide at 0, at 1, or later\"\n",
    "# pOther = 【P_A0, P_A1, P_B0, P_C0, P_D0, P_E0】\n",
    "# Calculate expected utility of action a given other players' mixing prob. (Turocy,for network)\n",
    "\n",
    "def payoff(pOther, cost):\n",
    "    '''This is for network (a)'''\n",
    "    diff_pay = 100\n",
    "    q = 0.7\n",
    "    if cost=='h':\n",
    "        c = 8\n",
    "    elif cost == 'l':\n",
    "        c = 1\n",
    "    # We know there are 5 nodes, pay matrix is 5X3 \n",
    "    pay =  np.empty((5,6),dtype=float)\n",
    "    pay[:,0] = ProbCorrect_Bayesian(1,q)*diff_pay\n",
    "    \n",
    "    # Note: since here the payoff is about turn 0,1,2, need to do the max in the \n",
    "    # indegree = 2, node A\n",
    "    pay_1 = ProbCorrect_Bayesian(3,q)*diff_pay - c\n",
    "    pay_2 = pOther[2]*ProbCorrect_Bayesian(3,q)*diff_pay + (1- pOther[2])*ProbCorrect_Bayesian(5,q)*diff_pay - 2*c\n",
    "    pay[0,1] = pay_1\n",
    "    pay[0,2] = pay_2\n",
    "    j = 2\n",
    "    i = 1\n",
    "    while j < 5:        \n",
    "        j = j+ 1\n",
    "        pay[0,j] = pay_2 - i*c\n",
    "        i += 1\n",
    "    \n",
    "    # indegree = 3, node B\n",
    "    pay_1 = ProbCorrect_Bayesian(4,q)*diff_pay - c\n",
    "    pay_2 = pOther[0]*ProbCorrect_Bayesian(4,q)*diff_pay + (1- pOther[0])*ProbCorrect_Bayesian(5,q)*diff_pay - 2*c\n",
    "    pay[1,1] = pay_1\n",
    "    pay[1,2] = pay_2\n",
    "    j = 2\n",
    "    i = 1\n",
    "    while j < 5:        \n",
    "        j = j+ 1\n",
    "        pay[1,j] = pay_2 - i*c\n",
    "        i += 1\n",
    "    \n",
    "    \n",
    "    # indegree = 1, node C or D, (position 2,3)\n",
    "    pay_1 = ProbCorrect_Bayesian(2,q)*diff_pay - c\n",
    "    pay_2 = pOther[0]*ProbCorrect_Bayesian(2,q)*diff_pay + (1-pOther[0])*ProbCorrect_Bayesian(4,q)*diff_pay - 2*c\n",
    "    pay_3 = pOther[0]*ProbCorrect_Bayesian(2,q)*diff_pay + (pOther[1] + (1-pOther[0]-pOther[1])*pOther[2])*ProbCorrect_Bayesian(4,q)*diff_pay + (1-pOther[0]-pOther[1])*(1-pOther[2])*ProbCorrect_Bayesian(5,q)*diff_pay - 3*c\n",
    "       \n",
    "    pay[2,1] = pay_1\n",
    "    pay[3,1] = pay_1\n",
    "    pay[2,2] = pay_2\n",
    "    pay[3,2] = pay_2\n",
    "    pay[2,3] = pay_3\n",
    "    pay[3,3] = pay_3\n",
    "    j = 3\n",
    "    i = 1\n",
    "    while j < 5:        \n",
    "        j = j+ 1\n",
    "        pay[2,j] = pay_3 - i*c\n",
    "        pay[3,j] = pay_3 - i*c\n",
    "        i += 1\n",
    "    \n",
    "    \n",
    "    # indegree = 1, node E, (position 4)\n",
    "    pay_2 = pOther[0]*ProbCorrect_Bayesian(2,q)*diff_pay + (1- pOther[0])*ProbCorrect_Bayesian(3,q)*diff_pay - 2*c\n",
    "    pay_3 = pOther[0]*ProbCorrect_Bayesian(2,q)*diff_pay + (pOther[1] + (1-pOther[0]-pOther[1])*pOther[2])*ProbCorrect_Bayesian(3,q)*diff_pay + (1-pOther[0]-pOther[1])*(1-pOther[2])*ProbCorrect_Bayesian(5,q)*diff_pay - 3*c\n",
    "\n",
    "    pay[4,1] = pay_1\n",
    "    pay[4,2] = pay_2\n",
    "    pay[4,3] = pay_3\n",
    "    j = 3\n",
    "    i = 1\n",
    "    while j < 5:        \n",
    "        j = j+ 1\n",
    "        pay[4,j] = pay_3 - i*c\n",
    "        \n",
    "        i += 1\n",
    "\n",
    "    return pay\n",
    "    \n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[70.   , 77.4  , 81.692, 80.692, 79.692, 78.692],\n",
       "       [70.   , 77.4  , 81.692, 80.692, 79.692, 78.692],\n",
       "       [70.   , 69.   , 76.4  , 80.692, 79.692, 78.692],\n",
       "       [70.   , 69.   , 76.4  , 80.692, 79.692, 78.692],\n",
       "       [70.   , 69.   , 76.4  , 80.692, 79.692, 78.692]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "payoff([0,0,0,0,0,0], 'l')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[70.   , 77.4  , 76.4  , 75.4  , 74.4  , 73.4  ],\n",
       "       [70.   , 77.4  , 79.046, 78.046, 77.046, 76.046],\n",
       "       [70.   , 69.   , 72.2  , 71.2  , 70.2  , 69.2  ],\n",
       "       [70.   , 69.   , 72.2  , 71.2  , 70.2  , 69.2  ],\n",
       "       [70.   , 69.   , 72.2  , 71.2  , 70.2  , 69.2  ]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "payoff([0.5,0.3,1,1,1,1], 'l')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Compute QRE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assumptions:\n",
    "    * Five agents\n",
    "    * Five nodes, the prob of deciding at turn 0\n",
    " #   <font color='red'> *  Next: prob of deciding at turn 1, given turn 0's outcome! </font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-35-f1cce35e0ea9>:23: RuntimeWarning: invalid value encountered in log\n",
      "  H= p_0s_t - np.log(1 - np.exp(p_0s_t))  - l*(payoff(p_real, cost)[:,0]- payoff(p_real, cost)[:,1:].max(1))\n",
      "<ipython-input-35-f1cce35e0ea9>:25: RuntimeWarning: invalid value encountered in log\n",
      "  p_A2 = np.log(1 - np.exp(p[0]) - np.exp(p[1]))\n",
      "<ipython-input-35-f1cce35e0ea9>:60: UserWarning: You have mixed positional and keyword arguments, some input may be discarded.\n",
      "  lgd = ax.legend(handles, labels = ['A','A_1','B', 'C','D','E'], loc = 0)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEYCAYAAABMVQ1yAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZhcdZX/8feppbs66c6eANlDCEsCJmJAUBCUZeLKICoRJ4jKYMbBfQMfnxF/zjDoz+GnM+AwCIrjwjIgiAwakE1WTViUsIewpMnWWTudXqvq/P64t5NKp7q7ulNdlXv783qefrqq7nZuVfepb33vt77H3B0REYm+RLUDEBGR8lBCFxGJCSV0EZGYUEIXEYkJJXQRkZhQQhcRiQkldBGRmFBCjxgzczM7pNpxVJKZXWdm/1yB41xiZr8o074OMLM/mtkOM/u3EtY/z8we2sdjTjezFjNL9rPex8zsrn05luyflNB7Ef6DPW1mrWa23sx+ZGajC5ZfYmZd4T/QNjN7xMyOL1h+spnlw+WFP8cXP2K8lSNh9bHfXPjcNpvZU2b2vnIfZxAuADYBo9z9y5U4oLu/7u717p7rZ71fuvvplYipkJndb2bn97F8ZthgSVU4ri+G/+PbzewnZlZbyeOXkxJ6EWb2ZeC7wFeB0cBxwEzgLjNLF6x6o7vXAxOA+4D/6bGrteE/WOHPo0N/BsPOo+HrMAa4FrjJzMb1XKnCiWIG8Kzrq9gVY4EB5TQz+xvgIuAUgv/xg4Fvlz+6ylBC78HMRhG8oJ9199+7e5e7vwp8BJgFnNNzG3fPAr8EppjZxEEc8xNm9tuC+6vM7KaC+2vMbEHBJqea2UtmttXMrjQzK1j3k2b2XLhsmZnNKFjmZra0t217xJQ0s2+Y2ctht8HjZjYtXPY2M1setmiWm9nbCrY7z8xWh9u8En68PwK4Cji++xPNQJ+jUrh7HvgJUAccHH6KutnMfmFmzcB5ZjbZzG43sy3h8/z3PXaTMbMbw/ifMLP5vR2vt+fBzK4DPg58LTzfU4tsOz6Mo9nM/gzM7rH8cDO7O4zzBTP7SMGyOjP7NzN7LTz2Q+Fje7Rwi70WBY8/VLC/vl7P+83sO2b2cLifu8xsQi/Px1gzu8PMmsK/rzvMbGq47F+AE4ErwufkiiK7+GP4e1u4zvHWoxusyDneb2b/YmYPA60Er3vJf+fh63Stuz/j7luB7wDn9bLu/s/d9VPwAywCskCqyLKfAb8Mb18C/CK8XQNcRvAROxU+djLQWOIxDwa2EbzBHgS8BrxRsGwrkAjvO3AHQWt0OtAELAqX/S2wCjgCSAHfBB4pOE6v2xaJ6avA08BhgAHzgfHAuDCeJeExPhreHw+MBJqBw8J9HATMC2+fBzzU4xgXhedd9KdgveuAf+4lzl37DeP5PLCD4JPVJUBX+LwkCBL9A8CPgAywIHwOTil4TbuADwFp4CvAK0C6yHF7fR76izlcfgNwU/icHQm8UXAeI4E1wCfCfR9N8LfV/VxeCdwPTAGSwNuAWoIWpofblPRalHAe9wMvA4eGz9/9wGW9nNN44CxgBNBA8In1toLl9wPn9/Gc7Iq/4LFLCP/Piq0T7vN1YF4Yf5qB/Z3/BTi74P6EcPvx1c5Fg/lRC31vE4BNHrS6e1oHFLbAPxK2NtuAvwc+1GO7yRb0rxf+jOy5U3dfTZCEFgAnAcuAN8zs8PD+gx60Prtd5u7b3P11gq6e7tb7p4F/dffnwjguBRYUttL72Lan84FvuvsLHviLu28G3gu85O4/d/esu18PPA+8P9wuDxxpZnXuvs7dn+ll/7j7Ze4+pref3rYr4rjwdVhPkJDOdPft4bJH3f228PmbAJwAfN3d2939KeAagmTW7XF3v9ndu4DLCRL/cUWO2d/z0CsLLlqeBfyTu+9095UEjYVu7wNedfefhvt+ArgF+JAFXQqfBD7v7m+4e87dH3H3jiKHKuW1KOU8furuL7p7G8GbUNG/GXff7O63uHuru+8A/oXg73eoXedBCzsbvm5Q+t95PbC94H737YYhinVIKaHvbRMwwYr3tx5E8G7f7aYw8RwArATe0mP9tUUS1c5ejvsAQav+HeHt+wn+GU4K7xdaX3C7leCPEoJ+2x92v3kAWwha11NK2LanaQQts54mE3yCKPQaMCU8t7OBpcA6M/vf8E1pqD0WPrcT3P04d/9DwbI1BbcnA1vCZNPtNfZ8fnatH74JNIbb9dTr81BCvBMJWpOFsRXuawbw1sKGAPAx4ECCN6UMxV+bXQbwWpRyHiX9zZjZCDP7r7ArqJmgC2WM9TPqpgzWFHms1L/zFmBUwf3u2zuKrLvfU0Lf26NAB/DBwgfDlvW72Tu54u6bCFrHl5jZQYM8bndCPzG8/QC9J/TerAE+3eMNpM7dHxlEPGvo0a8bWkuQcApNJ+gywN2XuftpBG9+zwM/DtfZ6+KgBX30PUcB7foZRMzFFB53LTDOzApbX7tiD00riC8BTA2366nP56EfTQTdetMKHptecHsN8ECP17He3f+BoMHRTvHXZg99vBblOo+evkzQRfdWdx9F0DiBoFEBRf4GeoZc5LGdBF043Q4scbtSPUPQndhtPrAh/DQaOUroPYQf1b8N/IeZLTKztJnNJOgP3ERw8bPYds8TdJV8bZCHfgB4J1Dn7o3AgwT9+eOBJ0vcx1XAxWY2D8DMRpvZhwcZzzXAd8xsjgXeZGbjgTuBQ83sHDNLmdnZwFzgDgvGXn8gfPPrIGj9dA+h2wBMNbOa7gO4+6W+9yigXT+DjLtX7r4GeAT4VzPLmNmbgE+x52v6FjP7YPgJ7QvheTxWZHe9Pg8lxJEDfk3QABhhZnMJLs51uyPc95Lw7y9tZseY2RG++8Lv5RZc4E2GFw/3GGrXz2tRlvMoooGg+3GbBaOMvtVj+QaCa0K9aSLoJipc5yngHRaMsR8NXDyIuPry38CnzGyumY0luO50XZmPUTFK6EW4+/eAbwDfJ/jo9QpBK+HUPrpMAP4vcIGZTQrvTy7S8jyrl2O+SPBP92B4vxlYDTzs/YwrLtjHrQTDLW8IP/KuJPhUMRiXE/SX3kVwce1agjebzQR9vF8GNhO8gb0v/JSSCB9fS9DdcxLwmXB/9xK0htab2aZBxlQOHyW4sLYWuBX4lrvfXbD8NwRdFd0XCj9Y0C+7Sz/PQykuJOgGWE+QQH5asO8dwOnA4jDO9QSva3fS/grBBevlBM/zd9n7f7mv16Kc51HoBwQXTjcRvAn+vsfyHxJcB9hqZv9eJJZWgn73h8OupuPC1+ZG4K/A4wzujaZX7v574HsE/eyvhT8934giw9w1TLY/ZvZJglb728OLLCIi+x0l9BKZ2RKgy91vqHYsIiLFKKGLiMSE+tBFRGKiopPgFJowYYLPnDmzWocXEYmkxx9/fJO7F51ipGoJfebMmaxYsaJahxcRiSQz6/lFsF3U5SIiEhNK6CIiMaGELiISE1XrQxcRKZeuri4aGxtpb2+vdihlk8lkmDp1Kul0uv+VQ0roIhJ5jY2NNDQ0MHPmTHqvZREd7s7mzZtpbGxk1qxZJW9XUpdLOEnVCxZUeLmoyPKTLah28lT4808DiF1EZJ+0t7czfvz4WCRzADNj/PjxA/7E0W8LPZzL+ErgNIK5oZeb2e3u/myPVR909/2hOK+IDENxSebdBnM+pbTQjwVWuftqd+8kKJ11xoCPVCYP3ngjV33qs7zy9F+qFYKIyH6plIQ+hT0rgjRSvCrL8Wb2FzP7Xfd83D2Z2QVmtsLMVjQ1NRVbpV+NzzzPzpZXeG1lr5XNRESq4tZbb8XMeP7556ty/FISerF2f88ZvZ4AZrj7fOA/gNuK7cjdr3b3he6+cOLEot9c7VfNyKB4SevWISkcLyIyaNdffz0nnHACN9xQnUlZS0nojexZKmuvklzu3uzuLeHtO4G0mU0oW5QF6kYFhWzat/dVZ0JEpLJaWlp4+OGHufbaa6uW0EsZtrgcmGNmswjqDC4GzilcwcwOJKjD52Z2LMEbxZDU5BsxdjQAHa2tQ7F7EYm4b//2GZ5d21zWfc6dPIpvvb9oT/Iut912G4sWLeLQQw9l3LhxPPHEExx99NFljaM//bbQ3T1LUC5rGfAcQaX7Z8xsqZktDVf7ELDSzP4C/Duw2IdoovX6CeMA6GqNzxcIRCT6rr/+ehYvXgzA4sWLuf766yseQ0lfLAq7Ue7s8dhVBbevAK4ob2jFjTnwAACyHXuVeRQR6bclPRQ2b97Mvffey8qVKzEzcrkcZsb3vve9ig6njNxcLmMOOgiAXKcSuojsH26++WbOPfdcXnvtNV599VXWrFnDrFmzeOihhyoaR/QS+tgJQIJcNlvtUEREgKC75cwzz9zjsbPOOotf/epXFY0jcnO5pBIpIE0+m6t2KCIiANx///17Pfa5z32u4nFEroUOYJYmn89XOwwRkf1KNBM6SfJ5tdBFRApFM6FbEncldBGRQtFM6CRxV5eLiEihaCZ0Mxy10EVECkU0oSdwNGxRRKRQNBN6Mom7vlgkIlIoogk9AXSRy6nbRUT2H6XOh75o0SLGjBnD+95X3iJvkUzoiVQCcDq2lndGNRGRfVHqfOhf/epX+fnPf17240fum6IAiZog7JYNTYyYMLbK0YjIfuV3F8H6p8u7zwOPgndf1ucq3fOh33fffXzgAx/gkksu6XXdU045pei3S/dVJFvoydo0ANs3bqpyJCIigWLzoVdaJFvoqUwNAM2bBleXVERirJ+W9FC5/vrr+cIXvgDsng+90gUuopnQR9YC0LJla5UjERHRfOj7pGbUSAB2NuuiqIhUn+ZD3weZ0WGh6B0qFC0i1TfQ+dBPPPFEPvzhD3PPPfcwdepUli1bVpY4ItnlUheObOloa6tyJCIiA58P/cEHHxySOCLZQm84YDwAXe0dVY5ERGT/EckWeoMKRYvIfu7pp59myZIlezxWW1vLn/70pyE7ZiQT+uhRY4EkWdUVFZH91FFHHcVTTz1V0WNGssulLlUHpMnnlNBFRLpFNqGbpTQ5l4hIgUgm9EwqE9QVdbXQRUS6RTKhJywRlKHLqwydiOwfkskkCxYsYP78+Rx99NE88sgjFY8hkhdFAbAEeZWhE5H9RF1d3a6LoMuWLePiiy/mgQceqGgMkWyhQ1iGzpXQRWT/09zczNixlZ/aO8ItdMPROHQR2dN3//xdnt/Sd8WggTp83OF8/div97lOW1sbCxYsoL29nXXr1nHvvfeWNYZSRDahWyKBaxy6iOwnCrtcHn30Uc4999xdsy9WSmQTOkmju65oMpmsdjQisp/oryVdCccffzybNm2iqamJSZMmVey4JfWhm9kiM3vBzFaZ2UV9rHeMmeXM7EPlC7GXY4V1Rds2bxvqQ4mIDMjzzz9PLpdj/PjxFT1uvy10M0sCVwKnAY3AcjO73d2fLbLed4HyzAPZX1zpoFXesmET9ZMq+6SJiPTU3YcO4O787Gc/q3jvQSldLscCq9x9NYCZ3QCcATzbY73PArcAx5Q1wl5YbfBE7WzaUonDiYj0aX/45nopXS5TgDUF9xvDx3YxsynAmcBVfe3IzC4wsxVmtqKpad/qgSYywXtR6xZ1uYiIQGkJvdglWu9x/wfA172fgeHufrW7L3T3hRMnTiw1xqKSI4JC0W3bVYZORARK63JpBKYV3J8KrO2xzkLghnB4zgTgPWaWdffbyhJlEen6OgB2bt8+VIcQEYmUUhL6cmCOmc0C3gAWA+cUruDus7pvm9l1wB1DmcwB0mNGANC6o2UoDyMiEhn9JnR3z5rZhQSjV5LAT9z9GTNbGi7vs9+83J58fSs/f/Q1FoQt9PZWFYoWEYESv1jk7ncCd/Z4rGgid/fz9j2s3m1o7uDXT77B3Hc0ANDR1j6UhxMRiYzITc41Khzdkh01GoDOzs5qhiMiAsD69etZvHgxs2fPZu7cubznPe/hxRdfrGgMkUvoDZk0AJ6qA5Jku5TQRaS63J0zzzyTk08+mZdffplnn32WSy+9lA0bNlQ0jsjN5VIfttBzuRrM0uQ0QZeIVNl9991HOp1m6dKlux7r/tZoJUUuoTeECb2zKwWkyOWV0EVkt/WXXkrHc+WdPrf2iMM58Bvf6HX5ypUrectb3lLWYw5GBLtcwoTemcJIkc9X/+u2IiL7g8i10GtTSWpSCTo6U2RIKKGLyB76akkPlXnz5nHzzTdX/Lg9Ra6FDsFIl9aOVFhXVF0uIlJd73rXu+jo6ODHP/7xrseWL1+umqKlqK9N0dpJUIZOdUVFpMrMjFtvvZW7776b2bNnM2/ePC655BImT55c0Tgi1+UCwdDFHe1dYV1RtdBFpPomT57MTTfdVNUYItlCb8ik2NGeDVvoSugiIhDphN4VzCwT1hUVERnuIprQ0+xoz+JJA5w2VS0SEYlmQq+vTdHSnoV0UHujZcOmKkckIlJ9kUzoozIpWjqzsKuu6NYqRyQiUn2RTOgNmTTu4DVB+DtVV1REJKoJPRht6WELvW37jmqGIyJCMplkwYIFzJs3j/nz53P55ZeTz+crGkNkx6ED5EakSQLtKkMnIlVWV1fHU089BcDGjRs555xz2L59O9/+9rcrFkMkW+i7ptAdWQNAe4vK0InI/mPSpElcffXVXHHFFbh7xY4b0RZ6WLVoRIYaoL21rboBich+48GbXmTTmvJ+ap8wrZ4TP3LogLY5+OCDyefzbNy4kQMOOKCs8fQmki307jJ0HQ1BoeiOttZqhiMiUlQlW+cQ2RZ60IfeXj8SgM5OFYoWkcBAW9JDZfXq1SSTSSZNmlSxY0ayhb6ry8XqgBSdqisqIvuRpqYmli5dyoUXXoiZVey4kWyh16WTJBNGLpfCLEVWdUVFpMra2tpYsGABXV1dpFIplixZwpe+9KWKxhDJhG5m1Nem6OxKA2lyua5qhyQiw9z+MElgJLtcIOh26exKYSTJqQydiEiUE3qa9s4kRoK85kQXEYlyQk/R3l1XVGXoREQinNBrU7S2J1W1SEQkFN2EnkmxsyOBg+qKiogQ1YTe9GJQtagtgSfAXaNcRERKSuhmtsjMXjCzVWZ2UZHlZ5jZX83sKTNbYWYnlD/U0FO/giuPYVb+VVra83gCoItcl1rpIlI93dPndv9cdtllFY+h33HoZpYErgROAxqB5WZ2u7s/W7DaPcDt7u5m9ibgJuDwoQiYQxdBegTHb7yRXH4x+WTwrtTatIWGyZX7iq2ISKHC6XOrpZQW+rHAKndf7e6dwA3AGYUruHuL756FZiQwdDPSjBgHC87h0A2/YwLbyaeCQ6muqIgMd6V8U3QKsKbgfiPw1p4rmdmZwL8Ck4D3FtuRmV0AXAAwffr0gca621v/geTya1iSupvWsFD0zk1bBr8/EYmN+667mo2vrS7rPifNOJh3nndBn+t0f/W/28UXX8zZZ59d1jj6U0pCLzazzF4tcHe/FbjVzN4BfAc4tcg6VwNXAyxcuHDwrfgJh9A0+V383Rt386OaU6kBdm7ZPujdiYjsq/2hy6WUhN4ITCu4PxVY29vK7v5HM5ttZhPcfcj6QTYf9fccvvZeSAYzLbZtax6qQ4lIhPTXko6zUvrQlwNzzGyWmdUAi4HbC1cws0MsnCPSzI4GaoDN5Q52DzPfzjP5GSRsG6AydCIi/bbQ3T1rZhcCy4Ak8BN3f8bMlobLrwLOAs41sy6gDTjbh7hUR30mzb9l38Ok5MPAaNpbVLVIRKqnZx/6okWLKj50saTpc939TuDOHo9dVXD7u8B3yxta3xoyae7IH8/HR/wRgE6VoRORKtL0ufugvjZFFym2jwnGnne0qA9dRIa3yCb0ZCIoctE4/hAAOncMbZe9iMj+LrIJHYIJujoSo4AUXZ3t0LKx2iGJSJUM8WW7ihvM+UQ6odfXpshmazBLk3Vg+bXVDklEqiCTybB58+bYJHV3Z/PmzWQymQFtF8maot0aMinaupJAimwiBcuvgRO+COmBPQkiEm1Tp06lsbGRpqamaodSNplMhqlTpw5om4gn9DTbOsO6ogCtm+Dpm+Doc6sdmohUUDqdZtasWdUOo+oi3eXSkEnR0ZnGSJAjDwccBY/+CGLysUtEZCAintDTtHUkw7qiWTj+H6HpOXj53mqHJiJScRFP6Cla2xNhXdEcHHkW1B8Ij15Z7dBERCou2gm9NkVnVzqsK9oFqRqYvxhW3wfZzmqHJyJSUdFO6JkUnq/Zs67ohEPB89DcWN3gREQqLOIJPQ1eQ97yQJZcZxeMnREs3PpaVWMTEam0iCf0oIWeS4Zl6DZugTFhJaRtr1cxMhGRyot0Qq/PpMBTZJPBLGctGzZBw2RIpGCbWugiMrxEOqGPyqQBIxsWit65aQskUzBqilroIjLsRDqhN2SCL7p21QYJvW1rWFd07Az1oYvIsBPxhJ4GoLMmuN/aXVd0zHS10EVk2Il0Qq+vDVro3Qm9fUdYV3TMTGhZD11t1QlMRKQKIp3Qa1IJalMJWjPBaXS0dif0cKTLdo1FF5HhI9IJHYJul9a67oTeHjyosegiMgxFevpcgFGZFNvDFnpne5jQd41FV0IXkeEjBi30FDvC0S6dnR3Bg/UHQrJGCV1EhpUYJPQ0WTJAmq6ucEKuRAJGT9NIFxEZViKf0OtrU3RlU5ilyGa7di/QWHQRGWYin9AbMim6utJAilyuIKFrLLqIDDMxSOhpOrqCuqLZfHb3gjEzghqjHS3VC05EpIJikNBTdHalMMIydN12jUVfU53AREQqLBYJ3fM1u+uKdhs7M/itfnQRGSZikdDJ14BZ8Ra6hi6KyDARg4SeDsrQGTgFCX3kREjV6cKoiAwbMUjoQQvdraCuKIBZ0Erf+mrVYhMRqaSSErqZLTKzF8xslZldVGT5x8zsr+HPI2Y2v/yhFtfdQs8ngrqiXR0duxeOnaEWuogMG/0mdDNLAlcC7wbmAh81s7k9VnsFOMnd3wR8B7i63IH2piGTwr2GbCKsWrR+8+6FY6arD11Eho1SWujHAqvcfbW7dwI3AGcUruDuj7j71vDuY8DU8obZu4baoMulK6wrunNTYUKfAe3boW1bpcIREamaUhL6FKBwMHdj+FhvPgX8rtgCM7vAzFaY2YqmpqbSo+xD0OVSS1cquCC6s2nr7oW7Rrqo20VE4q+UhG5FHvOiK5q9kyChf73Ycne/2t0XuvvCiRMnlh5lHzLpBClq6EgHCb21u64o7J4XXQldRIaBUuZDbwSmFdyfCqztuZKZvQm4Bni3u2/uuXyomBkNmVo62oMul7btzbsXjulO6OpHF5H4K6WFvhyYY2azzKwGWAzcXriCmU0Hfg0scfcXyx9m3xoyadprgg8Nbc07dy+oGws1DWqhi8iw0G8L3d2zZnYhsAxIAj9x92fMbGm4/Crgn4DxwI/MDCDr7guHLuw91demaM3kAejYWZDQd41FVwtdROKvpBJ07n4ncGePx64quH0+cH55QytdULUouL2rDF03zYsuIsNE5L8pCkGXy9a64HZHe8eeC7vHonvR67giIrERi4Q+KpNiW21wKl2dPRP6DOhsgbatRbYUEYmPWCT0+kyKnNeyR13Rbt1j0TWni4jEXCwSekMmRTabwiy9Z11R0Fh0ERk2Sroour/rnqDLcLK5Hgld86KLyDARk4QeVi0iR66wrihAZjRkxqiFLiKxF5Mul7CF3rMMXTeNRReRYSAeCT2ccRGseELXvOgiMgzEI6GHc6LvVVe025gwoWssuojEWEwSenp3GTp6SejZNthZnil7RUT2RzFJ6Ck8X0vefM+6ot12jUVXP7qIxFdsEjr5NPkEQJau1iLzuYCGLopIrMUioY+sSYHXkEsEc6K3bOjRtTI6nM5dCV1EYiwWCT2RMOqSI+hKBVPo7ty4Zc8VauthxASNdBGRWItFQgcYka6js7uu6JYiRaE1Fl1EYi42Cb2+ZuSuhN5WLKFrLLqIxFxsEnpDzQjaaoMRLq3NO/ZeYcx02L4G8vkKRyYiUhnxSeiZkbTWBFPntje37L3CmBmQ64Qd6yocmYhIZcQmoY+uHUFLJkzora17rzBG0+iKSLzFJqGPqquhuTboTulsa9t7BY1FF5GYi01Cb8ik2FYXjEPvbO/ce4VdY9HVQheReIpNQh+VSbO91gCjs2ddUYB0BuoP1NBFEYmt2CT0+toUbkFd0Wy2SAsdwqGLSugiEk+xSejBfC41mKX2rivabcx0JXQRia0YJfQ07jUYKbK5IlPoQjDSZfsb0NtyEZEIi1FCT4WFopN71xXtNmY6eA6a36hscCIiFRCrhE6+BixBrljVIigYuqiRLiISP/FJ6LVhoWh6KRQNuwtdqB9dRGIoPgk9bKG7Gd5bQh81FSyhFrqIxFJsEnp9WCjazXF6GeWSqoGGyRqLLiKxVFJCN7NFZvaCma0ys4uKLD/czB41sw4z+0r5w+xfOpkgRS35hOPeSa6zl6Q+dgZsXlXZ4EREKqDfhG5mSeBK4N3AXOCjZja3x2pbgM8B3y97hAOQSdbRkukA8rx832PFV5p+HKx9Etq2VjQ2EZGhVkoL/VhglbuvdvdO4AbgjMIV3H2juy+H3vo6KiOTquPVSUH5udWPPVl8pTl/EwxdfPm+CkYmIjL0SknoU4A1Bfcbw8cGzMwuMLMVZraiqamp/w0GqC5dx19mrAOSbGhcU3ylqQuhbiy8dHfZjy8iUk2lJHQr8pgP5mDufrW7L3T3hRMnThzMLvpUnx5BV62TSoymeefm4islkjD7FFh1t6oXiUislJLQG4FpBfenAmuHJpx9M7JmBAB1taPozG0jl8sVX3HO6bCzCdb10i0jIhJBpST05cAcM5tlZjXAYuD2oQ1rcBrSIwEYMW400MnrDz1efMVDTgVM3S4iEiv9JnQPvqVzIbAMeA64yd2fMbOlZrYUwMwONLNG4EvAN82s0cxGDWXgxYzKBAm97tCgO2fVQ38uvuLI8UFf+kt3VSo0EZEhlyplJXe/E7izx2NXFdxeT9AVU1Wja4OEnl04Ee5LsP71Pr4ROud0uO9SaGmC+vL354uIVFpsvqqVKoMAAAskSURBVCkKMLauHoDmdDa4MNqypfeV55wGOLx8T2WCExEZYrFM6NvbdzKydiwdfV0YPXA+jJwELy6rYIQiIkMnZgk9g3uC7e0tTJh0AO7trP3zX4uvnEgE3S4v36OCFyISC7FK6A2ZGsjXsKOzlelvmgfAS3/8U+8bzDkN2rdD4/IKRSgiMnRiltCDqkU7u1o5fNFJgLHu1T5mVpz9TrCkRruISCzELqGTr6Gls5URE8aSTIxmW/Om3jfIjIbpx2s8uojEQswSelAoui3bCsDImjF05Lb1vdGc02DD00HxaBGRCItZQg+6XNqybQCMn3AA7m2se/LZ3jc69G+C36vUSheRaItVQq9NJTCvoSMXJPRp8w4H4IV7H+l9o4mHw+hp6nYRkciLVUI3M1KWoTPfDhBeGIW1L7/S10ZBt8vq+yHbUYEoRUSGRqwSOkDaMnTlg8TcMHkSSevnwigE49E7W+D1RysQoYjI0IhdQq9JZsh6+677dTVjaM/2c2F01jsgWQsvaviiiERX7BJ6JjGSLG20Z4OkPm7cRNx3svGZF3vfqGYkzDxB49FFJNJil9Anpo4Ey/GH1/8AwLQjDgPghXv66U6Zczpsfgm2rB7qEEVEhkTsEvqUzJFYdjy3vnQrAIedfgIAb7y4qu8N55wW/H7pD0MZnojIkIldQh9dV4s3H8Of1/+ZNc1rGDtrGglrYOu2fopSj58N42bDS5p9UUSiKXYJffbEkbRsfjNGgltXBa30unQJF0Yh6HZ55UHobB3iKEVEyi92Cf2DR09lVHoCozmK36z6Ddl8lrFjJ5L3Fra+3EcFI4BDT4dcB7z6YGWCFREpo9gl9JG1KT7+tpmsazyKjW0beWTtI0w99BAAnrvrj31vPOPtUFMPf/ovyOcrEK2ISPnELqEDnPe2maTa51HDKG558RYOO+1EABqf7+fCaKoWTr0kKHrx4PeHPE4RkXKKZUIfN7KGxcfMonXLm3mg8QGY3oBZPVu3bOx/42POhzedHRSQ1vwuIhIhsUzoAJ86YRZd2xeS8xy/ffm31KXG0NpVwoVRM3jfD+CAeXDL+bD11SGPVUSkHGKb0KeNG8H7j1iAt83kf164hTGjx5P3ZravWdf/xjUj4OyfgzvcuAS62oY+YBGRfRTbhA7w6ZMOpn3rQta0vEZm6lgAXvh9PxdGu407GD54Naz/K/zvl4PkLiKyH4t1Qj/8wFG8/aBTIF/L09PXA/D6sy+UvoPDFsFJX4enfgmP/3SIohQRKY9YJ3SAC0+eS+f2+dyaux+zEWzZvGFgOzjp63DIqXDn16BxxdAEKSJSBrFP6MfMHMfsuneR9Q7SiYbSLowWSiThgz+GUQfBTedCSz9TCIiIVEnsEzrAF084hVz7AWTTRi7fTMv6ASblEePg7F9A62a4+ROQyw5NoCIi+2BYJPRTjjiAsbkT2FC/FXCe//0gvtp/0Hx47+XBtABXLIQ//xg6d5Y9VhGRwRoWCT2RMD79lg/z7JSgFN2alc8Nbkdv/ljQUh8xHu78Clw+F+75P7BjfRmjFREZnGGR0AE+uvAI1oyeilmG9evX0NG8Y3A7OuL9cP4f4JPLgipHD14O/+9IuO0zsOGZ8gYtIjIA5iWMrzazRcAPgSRwjbtf1mO5hcvfA7QC57n7E33tc+HChb5iRWVHjXzrrl8z9rrfkss1AWlG1hzIzNmH8dYlH2Ts7OmD2+nml+Gx/wyGNna1wswTYeoxMOHQ8OcQyIwu63mIyPBlZo+7+8Kiy/pL6GaWBF4ETgMageXAR9392YJ13gN8liChvxX4obu/ta/9ViOht3R08s5rTuddz03koC2jyOebcW8DEqQTE5g4aSrHf+j9zDzxmIHvvHVLMFb9LzfClpchX3DhtP6AMLnPgTEzIDMKakcFMzvWNkBt+LumAdJ1kExDIg2JYfMBSkRK1FdCT5Ww/bHAKndfHe7sBuAM4NmCdc4A/tuDd4fHzGyMmR3k7iV8z75y6mtr+Nzx3+HnDf/DvV3r6Mht5O0vNXDI+nHkutpZu/4JbrniCbgiQdAblQAMs+A3JDAsvN2bQ8KfYjaEPyIynNXWjeQzP/1x2fdbSkKfAqwpuN9I0Arvb50pwB4J3cwuAC4AmD59kF0c+2jJ0Sey5OgTd93fsrOVJ9auZmXTanY89hQNK9eR7Mrv+qq/dX+C2fVJpp8uqpJmCCh1GoEwhhLXFpFoSGfSQ7LfUhJ6sXzSMyOVsg7ufjVwNQRdLiUce8iNGzmCU+ccyalzjoS3faDa4YiIDFopnbSNwLSC+1OBtYNYR0REhlApCX05MMfMZplZDbAYuL3HOrcD51rgOGD7/tZ/LiISd/12ubh71swuBJYRDFv8ibs/Y2ZLw+VXAXcSjHBZRTBs8RNDF7KIiBRTSh867n4nQdIufOyqgtsO/GN5QxMRkYHQQGcRkZhQQhcRiQkldBGRmFBCFxGJiZIm5xqSA5s1Aa8NcvMJwKYyhhMlw/Xch+t5w/A99+F63tD3uc9w94nFFlQtoe8LM1vR2+Q0cTdcz324njcM33MfrucNgz93dbmIiMSEErqISExENaFfXe0Aqmi4nvtwPW8Yvuc+XM8bBnnukexDFxGRvUW1hS4iIj0ooYuIxETkErqZLTKzF8xslZldVO14KsHMppnZfWb2nJk9Y2afr3ZMlWZmSTN70szuqHYslRKWcrzZzJ4PX/vjqx1TpZjZF8O/9ZVmdr2ZZaod01Aws5+Y2UYzW1nw2Dgzu9vMXgp/jy11f5FK6GHB6iuBdwNzgY+a2dzqRlURWeDL7n4EcBzwj8PkvAt9Hniu2kFU2A+B37v74cB8hsn5m9kU4HPAQnc/kmDa7sXVjWrIXAcs6vHYRcA97j4HuCe8X5JIJXQKCla7eyfQXbA61tx9nbs/Ed7eQfCPPaW6UVWOmU0F3gtcU+1YKsXMRgHvAK4FcPdOd99W3agqKgXUmVkKGEFMK6C5+x+BLT0ePgP4WXj7Z8Dflrq/qCX03opRDxtmNhN4M/Cn6kZSUT8Avgbkqx1IBR0MNAE/DbuarjGzkdUOqhLc/Q3g+8DrBIXmt7v7XdWNqqIO6K74Fv6eVOqGUUvoJRWjjiszqwduAb7g7s3VjqcSzOx9wEZ3f7zasVRYCjga+E93fzOwkwF89I6ysM/4DGAWMBkYaWZ/V92ooiFqCX3YFqM2szRBMv+lu/+62vFU0NuBD5jZqwRdbO8ys19UN6SKaAQa3b37k9jNBAl+ODgVeMXdm9y9C/g18LYqx1RJG8zsIIDw98ZSN4xaQi+lYHXsmJkR9KU+5+6XVzueSnL3i919qrvPJHi973X32LfW3H09sMbMDgsfOgV4toohVdLrwHFmNiL82z+FYXJBOHQ78PHw9seB35S6YUk1RfcXvRWsrnJYlfB2YAnwtJk9FT72jbDWq8TXZ4Ffho2X1QyT4uvu/iczuxl4gmCE15PEdBoAM7seOBmYYGaNwLeAy4CbzOxTBG9uHy55f/rqv4hIPESty0VERHqhhC4iEhNK6CIiMaGELiISE0roIiIxoYQuIhITSugiIjGhhC5SwMyOMrPXzOwfqh2LyEApoYsUcPenCaYYOLfasYgMlBK6yN42AvOqHYTIQCmhi+ztMqDWzGZUOxCRgVBCFylgZouAkcD/ola6RIwSukgoLET8PeAzwNPAkdWNSGRglNBFdvsm8N/u/ipK6BJBSugiQFhI4jSC+qWghC4RpPnQRURiQi10EZGYUEIXEYkJJXQRkZhQQhcRiQkldBGRmFBCFxGJCSV0EZGY+P/7LwAV7aJtdAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "#Define the function whose 0 we want to find  (Turocy)\n",
    "def F(p, l = 1,cost='h'): ## probability vector and value of lambda\n",
    "    # The part where it is about prob of deciding at turn 0\n",
    "  \n",
    "    p_0s = np.delete(p, 1)\n",
    "    p_0s_t =  np.transpose(p_0s)\n",
    "    p_real = np.exp(p)\n",
    "#     p_1s_t = np.log(1 - np.exp(p_0s_t))  \n",
    "#     p_1s_t[0] = p[1]\n",
    "\n",
    "#     pay_turn1 = payoff(p, cost)[:,1:].max(1)\n",
    "#     pay_turn1[0] = payoff(p, cost)[0,1]\n",
    "    \n",
    "#     H= p_0s_t - p_1s_t - l*(payoff(p, cost)[:,0]- pay_turn1)\n",
    "    \n",
    "#     p_A2 = np.log(1 - np.exp(p[0]) - np.exp(p[1]))\n",
    "   \n",
    "#     A_1 = p[1] - p_A2 - l*(payoff(p,cost)[0,1] - np.max(payoff(p,cost)[0,2:]))\n",
    "\n",
    "\n",
    "    # This version: the probability is the difference between prob of deciding at turn 0 vs later\n",
    "    \n",
    "    H= p_0s_t - np.log(1 - np.exp(p_0s_t))  - l*(payoff(p_real, cost)[:,0]- payoff(p_real, cost)[:,1:].max(1))\n",
    "    # But additionally, we need subjects to expect the prob of A deciding at turn 1 or later\n",
    "    p_A2 = np.log(1 - np.exp(p[0]) - np.exp(p[1]))\n",
    "    A_12 =  p[1] - p_A2 - l*(payoff(p_real,cost)[0,1] - np.max(payoff(p_real,cost)[0,2:]))\n",
    "    \n",
    "    H_A = np.append(H, A_12)\n",
    "    A_02 = p[0] - p_A2 - l*(payoff(p_real,cost)[0,0] - np.max(payoff(p_real,cost)[0,2:]))\n",
    "    H_A = np.append(H_A, A_02)\n",
    "    A_01 = p[0] - p[1] - l*(payoff(p_real,cost)[0,0] - payoff(p_real,cost)[0,1])\n",
    "    H_A = np.append(H_A, A_01)\n",
    "    \n",
    "    res =  np.dot(np.transpose(H_A),H_A)\n",
    "#     res =  np.dot(np.transpose(H),H) \n",
    "    return res\n",
    "\n",
    "cons = ({'type': 'ineq', 'fun' : lambda p: 1- np.exp(p[0]) - np.exp(p[1])-0.00000001 } ,\n",
    "       {'type': 'ineq', 'fun' : lambda p: 1- np.exp(p[2])-0.00000001 },\n",
    "        {'type': 'ineq', 'fun' : lambda p: 1- np.exp(p[3])-0.00000001 },\n",
    "        {'type': 'ineq', 'fun' : lambda p: 1- np.exp(p[4]) -0.00000001},\n",
    "        {'type': 'ineq', 'fun' : lambda p: 1- np.exp(p[5])-0.00000001},\n",
    "       )\n",
    "\n",
    "\n",
    "cost='l'\n",
    "lSeq = np.arange(0,10,.25)\n",
    "QRE = np.zeros((len(lSeq),6))\n",
    "p0 = np.log([0.25,0.25,0.5,0.5,0.5,0.5])\n",
    "for (i,l) in enumerate(lSeq):\n",
    "    res = minimize(F, p0, args=(l,cost,), method='SLSQP',constraints=cons)\n",
    "    p0 = res.x # use the new value as the new guess for next lambda\n",
    "    QRE[i] =  np.exp(res.x)\n",
    "    \n",
    "#Plot\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(lSeq, QRE[:,0],lSeq, QRE[:,1],lSeq, QRE[:,2],lSeq, QRE[:,3],lSeq, QRE[:,4],lSeq, QRE[:,5])\n",
    "ax.set_xlabel('$\\lambda$')\n",
    "handles, labels = ax.get_legend_handles_labels()\n",
    "lgd = ax.legend(handles, labels = ['A','A_1','B', 'C','D','E'], loc = 0)\n",
    "ax.set_title(\"QRE when cost=\"+str(cost)+\" Prob of decision at turn 0\")\n",
    "plt.show()    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "     fun: nan\n",
       "     jac: array([nan, nan, nan, nan, nan, nan])\n",
       " message: 'Inequality constraints incompatible'\n",
       "    nfev: 7\n",
       "     nit: 1\n",
       "    njev: 1\n",
       "  status: 4\n",
       " success: False\n",
       "       x: array([-9.56420080e-01, -4.92303448e-01, -1.72831904e+00,  2.25477965e-04,\n",
       "        2.24761876e-04,  2.24739599e-04])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-35-f1cce35e0ea9>:23: RuntimeWarning: invalid value encountered in log\n",
      "  H= p_0s_t - np.log(1 - np.exp(p_0s_t))  - l*(payoff(p_real, cost)[:,0]- payoff(p_real, cost)[:,1:].max(1))\n",
      "<ipython-input-35-f1cce35e0ea9>:25: RuntimeWarning: invalid value encountered in log\n",
      "  p_A2 = np.log(1 - np.exp(p[0]) - np.exp(p[1]))\n",
      "<ipython-input-35-f1cce35e0ea9>:7: RuntimeWarning: overflow encountered in exp\n",
      "  p_real = np.exp(p)\n",
      "<ipython-input-35-f1cce35e0ea9>:23: RuntimeWarning: overflow encountered in exp\n",
      "  H= p_0s_t - np.log(1 - np.exp(p_0s_t))  - l*(payoff(p_real, cost)[:,0]- payoff(p_real, cost)[:,1:].max(1))\n",
      "<ipython-input-5-439751b583a0>:48: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  pay_3 = pOther[0]*ProbCorrect_Bayesian(2,q)*diff_pay + (pOther[1] + (1-pOther[0]-pOther[1])*pOther[2])*ProbCorrect_Bayesian(4,q)*diff_pay + (1-pOther[0]-pOther[1])*(1-pOther[2])*ProbCorrect_Bayesian(5,q)*diff_pay - 3*c\n",
      "<ipython-input-5-439751b583a0>:67: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  pay_3 = pOther[0]*ProbCorrect_Bayesian(2,q)*diff_pay + (pOther[1] + (1-pOther[0]-pOther[1])*pOther[2])*ProbCorrect_Bayesian(3,q)*diff_pay + (1-pOther[0]-pOther[1])*(1-pOther[2])*ProbCorrect_Bayesian(5,q)*diff_pay - 3*c\n",
      "<ipython-input-35-f1cce35e0ea9>:25: RuntimeWarning: overflow encountered in exp\n",
      "  p_A2 = np.log(1 - np.exp(p[0]) - np.exp(p[1]))\n",
      "<ipython-input-35-f1cce35e0ea9>:38: RuntimeWarning: overflow encountered in exp\n",
      "  cons = ({'type': 'ineq', 'fun' : lambda p: 1- np.exp(p[0]) - np.exp(p[1])-0.00000001 } ,\n",
      "<ipython-input-35-f1cce35e0ea9>:40: RuntimeWarning: overflow encountered in exp\n",
      "  {'type': 'ineq', 'fun' : lambda p: 1- np.exp(p[3])-0.00000001 },\n",
      "<ipython-input-35-f1cce35e0ea9>:41: RuntimeWarning: overflow encountered in exp\n",
      "  {'type': 'ineq', 'fun' : lambda p: 1- np.exp(p[4]) -0.00000001},\n",
      "<ipython-input-35-f1cce35e0ea9>:42: RuntimeWarning: overflow encountered in exp\n",
      "  {'type': 'ineq', 'fun' : lambda p: 1- np.exp(p[5])-0.00000001},\n",
      "<ipython-input-37-f1a579afe64a>:15: UserWarning: You have mixed positional and keyword arguments, some input may be discarded.\n",
      "  lgd = ax.legend(handles, labels = ['A','A_1','B', 'C','D','E'], loc = 0)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEYCAYAAABMVQ1yAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXhU5d3/8fc3k40AgZCAAmETExFUUHGpS3FpFVdq3bAtrk8tj1q72Cq21mp93Gq12p9trSvWVqhatdpaqfuGC2hRQZFNlgCRsBPIPvfvj3OCw5BlQmbmZGY+r+vKlZkzZ875nlk+c899ztzHnHOIiEjqywq6ABERiQ8FuohImlCgi4ikCQW6iEiaUKCLiKQJBbqISJpQoIuIpAkFeoowM2dmewZdRzKZ2VQz+78A1nuUmVXEaVlmZg+Z2QYzey+G+Yf6z3V2J9c7z8yOameewWZWbWahzqxLug4FehQzO9/MPjazbWZWaWZ/MLNeEbdfZ2YN/htho5nNNLOvRNx+lJmF/dsj/77S8hrTm/94vpmE9TQHYfPjvdTMpiR6vTE4Avg6UOqcOzhZK3XOjXLOvdrOPMudcz2cc01JKgvY/h76SzvzLDWzryWrJn+dx5rZfP+9/4qZDUnm+uNBgR7BzK4AbgV+CvQCDgWGAv8xs5yIWf/mnOsBlACvAI9HLWqV/0aJ/Hs78VsgQG//uTkHuNbMxkfP0NnWbwcNAZY657YmcZ0Zr6PPsZmVAE8CvwD6ALOBvyWgtMRyzunPG/6gEKgGzoqa3gNYA5znX78O+EvE7SMBB/T1rx8FVMS4zguAZyOuLwIei7i+AhjjX3bAZGAhsAH4PWAR814IfOrfNgMYEnFbm/eNqikE/AxYDGwB3gcG+bcdBswCNvn/D4u43/nAEv8+nwPfBvYGaoEm/7Hd2MHnZKpf67/85b4LDG9l3qH+dmZHTJsF/KT5OQGuAiqBR4A84E5glf93J5AX+Rz6j8NaYCnw7TbqHAA8A6z3n8Pv+tMvitr+61t5vH/jr2cJcGnkduA1LB4AVgMrgf8DQhH3/67/vG8BPgEO8KcvBb7mXz4YL6A2A18Ad7T0mLW2HRGv+8eAP/vrmgeMbeMxuQvv9bvZfw0d6U8fD9QDDf5j8mEL930ECAM1/jxX0sL7KmobrwOeAP7ir/N/gFeBG4C3/Jr/A5S0Uu/FwMyI69399Y8IOps69J4JuoCu8ue/0BqJCISI2x4G/hrxwvmLfzkXuMV/Mza/KXZ64bWxzj2AjXjflPoDy4CVEbdtALL86w74J9AbGAxUAeP9277hvwH3BrKBa6JenK3et4Wafgp8DOwFGDAaKMZrtWwAJvnrOMe/Xuy/+DcDe/nL6A+M8i+fD7wZtY4p/na3+Bcx31S8cDnYX+dfgemt1D3U385sv+7DgW3Asf5z0oj37SsP6Ab8CngH6Af0BWYCN0Q8h43AHf7844CtzdvXwrpfA/4A5ANj/Mf32Na2P+q+k4H5wCD/MX6FHUP2aeBP/mPcD3gP+J5/25l4IX+Qv8174n+Qs2PYvQ1M8i/3AA6Nfsxi2I7r8D6cTsT7ELoZeKeN7fqO/9rIBq7A+yDNj34PtXH/7fW39r5i50BvwHsvZPnP8at4DZPyiOu3tLK+u4A/Rk2bC5wedDZ15C/wArrKn/8CrGzltluA/0S8cOrxwqcJWAccFTHvUXiti+ig6t7KslcABwATgXv9N+wIvNb7MxHzOeCIiOuPAVP8y/8GLoq4LQsvzIa0d98W6vkMmNDC9EnAe1HT3sYLrO7+Np4OdIua53zaCLR2npOpwP0R108E5rcy71B/OzfifdB8Clwe8ZzUNweKP20xcGLE9ePxukaa52+MfM78x+wXLax3kP866Bkx7WZgaizbD7wMTI64fhxffjDtBtRFPqZ4H6Sv+JdnAD9oZblL+TLsXgeuJ6p1yo4fgu1tx3XAixG3jQRqOvBcbgBGRywrEYH+etTtrwLXRFy/BHi+lfU9QFTY47Xsz9+V125Qf+pD/9JaoKSVvrf+eK2VZo8553rjveHmAgdGzb/KOdc76q+1PtTX8F6sX/Uvv4rXIhznX49UGXF5G15rC7x+2rv8nbQb8Vq1BgyM4b7RBuGFXbQBeN8gIi0DBvrbdjZea3O1mf3LzEa0svyOirXuZiXOuSLn3N7Oud9FTK9yztVGXI/enmX+tGYbop6z6Nsjl7PeObclat6BLczbkgF4H+qR9202BMjBe0ybn9s/4bXUofXnKtpFeK3U+WY2y8xObqWO9rYj+rnIb62v2syuMLNPzWyTX3cvvH1OibSihWmxvn6q8bpdIxXiddWkDAX6l97Gaw19M3KimXUHTmDncMU5txb4HnCdmfXfxfU2B/qR/uXXaD3QW7MC72t45AdIN+fczF2oZwUwvIXpq/ACJtJgvK/8OOdmOOe+jvfhNx+4z5/HRS/IzH7WwlFA2/92oeZYRNcRvT2D/WnNivznvrXbI5fTx8x6Rs27Msa6VuMFc+R9m63Ae02WRDyvhc65URG3t/Rc7cA5t9A5dw7eB8GtwBNR2xaP7djOzI7E219xFlDkN3424TUyoIXXREtlR13fChRErCOE11XW1n06Yh5e92Lz8rvjPbbzOrHMpFOg+5xzm/C+lv4/MxtvZjlmNhTvCJa1eP23Ld1vPt5X3yt3cdWvAUfjfa2uAN7A688vBv4b4zLuAa42s1EAZtbLzM7cxXruB24wszL/GOr9zKwYeA4oN7NvmVm2mZ2N97X7n2a2m5md6r8J6vBaO82Hwn0BlJpZbvMKnHM3uZ2PAtr+t4t1d9Q04Boz6+sf4XAt3g61SNebWa4fUCez89FMOOdW4PW/32xm+Wa2H16LuMXXSwseAy43s1IzK8Lbv9C87NV4O/JuN7NCM8sys+FmNs6f5X7gJ2Z2oP9c7dnSoXZm9h0z6+uca+4KhC+fn3htR6SeeF1WVUC2mV3Ljq3fL4ChZtZW/nyBtx+p2QK8bwQn+UecXYO3fyNengL2MbPTzSwf7/Xwkf/+ThkK9AjOuV/jHdnwG748WqMAr5+urcPObgMuNrPmr8IDWmh5nt7KOhfgBeAb/vXNeEc7vOViPD7YOfcUXstrupltxusGOiGW+7bgDryQ+Q/ejs4H8D5s1uGF2hV4+w2uBE72v6Vk+dNX4XX3jMPrrwSvj3geUGlma3expkT4P7wjPz7C2wn8gT+tWSVev+8qvFCb3Mab+xy8/uhVeMHwS+fcCzHWcR9eg+BDv4Yno24/F2/n+yd+PU/gfQvCOfc4cCPwKN7r9Wm8HavRxgPz/G8/dwETo7qf4rEdkWbg7ddZgNdtU8uO3SHNH4zrzOyDVpZxM94H7kYz+4nf4LoE70NsJV6LPS4//gJwzlXh7QO6Ee9xPgRvv1ZKMb/zX1pgZhfitdoPd84tD7oeEZG2KNDbYWaTgAbn3PSgaxERaYsCXUQkTagPXUQkTSRzTIsdlJSUuKFDhwa1ehGRlPT++++vdc5FH7IJBBjoQ4cOZfbs2UGtXkQkJZlZ9A/8tlOXi4hImlCgi4ikCQW6iEiaUKCLiKQJBbqISJpoN9DN7EEzW2Nmc1u53czsd2a2yMw+MrMD4l+miIi0J5YW+lS8wX1acwJQ5v9dDPyx82WJiEhHtXscunPudX8Y2dZMAP7svDEE3jGz3mbW3x/6M2U1NTVRs24jm1asYtPKNWxdu57q9Rup3VJNOBxu9X7OOXCOcNiBC/tnEgnjwnj/OzXSgsM5cGFvuWH35TrCYUfnhoMWkWQZMnJvjvj+uXFfbjx+WDSQHYfGrPCn7RToZnYxXiuewYMHR98cqKb6Bl6/aypz//s2DeEtOFeHdyY5EZH4Cs9t4ogELDcegW4tTGuxqeicuxfvvJmMHTu2SzQnG7bV8tJtf2L+/Fk0hTeSZT0ozO9Pbm4eeXndyC8ooFvPHnTvXUj34iJ69i0mKzenzWVm52RjoRAWChHKziIrO4SFsskKZZGV1cn90NkhskLZhHKzycrJJpTj/c/KziYUCnVu2SKS0uIR6BXseAqtUlo+VVeXUrNxCy/c8nsWL/2QsNtCKKsXI0ccy9eu+l9yCvKDLk9EpMPiEejPAJeZ2XS8s3xs6sr959WVVTx/6x9Zvnouzm0jO6sPo/c/mXE/uohQOy1vEZGurN1AN7NpeCcxLjGzCuCXeGcixzl3D965Jk8EFuGdVfuCRBUbD4/85Bq2NawkN9SPAw8/hUMmf0tdFSKSFmI5yuWcdm53wKVxqyiB1i1cyraG1ZT0HMF59/8m6HJEROIqo34p+u4jTwJhRh87rt15RURSTUYF+tIln2HWnX3POjHoUkRE4i5jAn3LqjXUNFRS1H2g+sxFJC1lTKC/89DjQBMjDzk46FJERBIiYwJ9yWfzMMvngPNOC7oUEZGEyIhAr9m4heq61RTmDyAnLy/ockREEiIjAv29Bx8DGthrf43sKyLpKyMCfeFHc4BcDr7g9KBLERFJmLQP9IZttWyuXU2PvP7kFfYMuhwRkYRJ+0Cf/ecnca6WPffeJ+hSREQSKu0Dff57s4BsDrngjKBLERFJqLQO9Kb6BjZsW0lBzu702L1v0OWIiCRUWgf6h3/7J85tY9ieI4IuRUQk4dI60Oe+MRPI4pDzdHSLiKS/tA30pqYm1m2uID97N4qGDWr/DiIiKS5tA33+sy8TdlsYNGjPoEsREUmKtA30OTNeAeDQSRq7RUQyQ9oGetXG5eSG+tFvVHnQpYiIJEVaBvriF2fSFN7IgN2HBV2KiEjSpGWgz37meQAOPvuUgCsREUmetAz0yqplZGcVM+iQMUGXIiKSNGkX6CvenUNjeB27lwwJuhQRkaRKu0Cf/fi/ADjw1OMCrkREJLnSLtA3rF8L5LHn148IuhQRkaRKu0Cvb6glZPlBlyEiknRpF+gNTXWEsnTeUBHJPGkX6E3hOnJCCnQRyTxpF+hh6sjNUaCLSOZJq0Bv2FaLc7Xk5xcEXYqISNKlVaCvW7wMgIIe3QOuREQk+dIq0NcvWQFA9169A65ERCT50irQN66qBKBXv5KAKxERSb60CvTNa9YD0HtQ/4ArERFJvrQK9K2bNgFQPHxwwJWIiCRfWgV6zdZqwOg9ZGDQpYiIJF1aBXptXQ1m+YRyc4IuRUQk6dIq0Osb68jSOC4ikqHSKtAbm+rIzsoNugwRkUCkVaA3uTpyQmqhi0hmiinQzWy8mX1mZovMbEoLt/cys2fN7EMzm2dmF8S/1PaFXS15uQp0EclM7Qa6mYWA3wMnACOBc8xsZNRslwKfOOdGA0cBt5tZUvs+tm3YBDRoHBcRyVixtNAPBhY555Y45+qB6cCEqHkc0NPMDOgBrAca41ppO9Yt+ByAgp49k7laEZEuI5ZAHwisiLhe4U+LdDewN7AK+Bj4gXMuHL0gM7vYzGab2eyqqqpdLLll65euBKBnn6K4LldEJFXEEujWwjQXdf14YA4wABgD3G1mhTvdybl7nXNjnXNj+/bt2+Fi27Jp9RcA9Npd47iISGaKJdArgEER10vxWuKRLgCedJ5FwOfAiPiUGJst6zYAUDRYvxIVkcwUS6DPAsrMbJi/o3Mi8EzUPMuBYwHMbDdgL2BJPAttz9bNWwAoKRuazNWKiHQZ2e3N4JxrNLPLgBlACHjQOTfPzCb7t98D3ABMNbOP8bpornLOrU1g3Tup2bYVCFGwW3EyVysiXUBDQwMVFRXU1tYGXUrc5OfnU1paSk5O7EOZtBvoAM6554DnoqbdE3F5FXBczGtNgPp6fxyXUCjIMkQkABUVFfTs2ZOhQ4fiHWyX2pxzrFu3joqKCoYNGxbz/dLml6L1jXWETCeHFslEtbW1FBcXp0WYA5gZxcXFHf7GkTaB3hiuIztLgS6SqdIlzJvtyvakTaA3uTpys/WzfxHJXGkR6E1NTThXS16eAl1EgvPUU09hZsyfPz+Q9adFoFevWgM00a2ge9CliEgGmzZtGkcccQTTp08PZP1pEejrFi0DoHvhTj9OFRFJiurqat566y0eeOCBwAI9psMWu7oNy70frvYs7hNwJSIStOufnccnqzbHdZkjBxTyy1NGtTnP008/zfjx4ykvL6dPnz588MEHHHDAAXGtoz1p0ULf9IX3G6Ze/fsFXImIZKpp06YxceJEACZOnMi0adOSXkNatNCr13vjuPQZWhpwJSIStPZa0omwbt06Xn75ZebOnYuZ0dTUhJnx61//OqmHU6ZFC31btTeOS3F57L+oEhGJlyeeeIJzzz2XZcuWsXTpUlasWMGwYcN48803k1pHWgR6Tc02IJduvXVyCxFJvmnTpnHaaaftMO3000/n0UcfTWodadHlUl9fS5Z+9i8iAXn11Vd3mnb55ZcnvY60aKE3NGkcFxGRtAj0RldHTkiBLiKZLS0CPezqyM3Rz/5FJLOlfKA31Tf447h0C7oUEZFApXygr1+yAnAUdO8RdCkiIoFK/UBfvByA7r00jouIZLaUD/SNKysBKOync4mKSGZL+UDfXOWP4zJg94ArEZFMF+t46OPHj6d3796cfPLJcV1/ygd69cZNAJQMHxRwJSKS6WIdD/2nP/0pjzzySNzXn/K/FN1WXQ1A0Z5DAq5ERLqEf0+Byo/ju8zd94UTbmlzlubx0F955RVOPfVUrrvuulbnPfbYY1v8dWlnpXwLvba2BrN8cvL0wyIRCU5L46EnW8q30OsbNI6LiERopyWdKNOmTeOHP/wh8OV46Mk+wUXKB3pjWOO4iEiwNB56nDSGNY6LiARL46HHSdjVkadxXEQkQB0dD/3II4/kzDPP5KWXXqK0tJQZM2bEpY6U7nKp27wFqCO/W0HQpYhIBuvoeOhvvPFGQupI6Rb62oXez/4LemgcFxGRlG6hb1i6AoAevXsHXImIyI4+/vhjJk2atMO0vLw83n333YStM6UDfeOqNQAU7lYScCUiIjvad999mTNnTlLXmdJdLtXr1gNQNKh/wJWIiAQvtQN902YAivWzfxGR1A70mm3VQBaFgwcEXYqISOBSOtDr6rxxXEKhUNCliEiGC4VCjBkzhtGjR3PAAQcwc+bMpNeQ0jtF6xv1s38R6Rq6deu2fSfojBkzuPrqq3nttdeSWkNKt9AbmurIzlKgi0jXsnnzZoqKipK+3pRuoYddHTkh/ahIRL5063u3Mn9922cM6qgRfUZw1cFXtTlPTU0NY8aMoba2ltWrV/Pyyy/HtYZYpHig15KX2y3oMjqmsR5wu3Zft4v3E0l3zkE4/OXlRC6/Fd26dWOOPwb622+/zbnnnsvcjz5qebRFAyz+HSQxBbqZjQfuAkLA/c65nQYcNrOjgDuBHGCtc25cHOvcSXVlFdBIfkEKjOPiHCx8Ad68A5a/HXQ1Iunn+Megsg6AqwYdD4OOj/86Kj9s+3YX3j7PV4YVsHZNJVXzXqVfSZ+d5+3RDwoHxr3EdgPdzELA74GvAxXALDN7xjn3ScQ8vYE/AOOdc8vNrF/cK42ydpE3jkv3nj0TvapdF26CeU/Bm3fCFx9DYSkc+RPI7cyHUPLGVhZJGfm9oWfAPzA0217D/AULaQo7igePgJaOwsvtnpASYmmhHwwscs4tATCz6cAE4JOIeb4FPOmcWw7gnFsT70KjbVxWAUCPPsnf8dCuxjqY8yi8dRds+BxKymHCH2DfMyE7N+jqRNLPp59Cz90DLaGmppYxR44HwDnHw39+hFDv+LfC2xJLoA8EVkRcrwAOiZqnHMgxs1eBnsBdzrk/Ry/IzC4GLgYYPHjwrtS73cbKKgB67963U8vpsOo1sLUKmuqhqcH731j35eV1i+DdP0F1JQw4AI67AfY6CbJS+oAiEWlHU1NT0CXEFOgtfceP3uuQDRwIHAt0A942s3eccwt2uJNz9wL3AowdO7ZTey6q128AoPeQ0s4spmPmPAr/uNTrK2vLsHHwzT95/5N4+ikRyWyxBHoFMCjieimwqoV51jrntgJbzex1YDSwgATZumULACV7dq6lH7P37oPnfgJ7HAVjL4RQLoRyIJQXcTkX8guhd5JqEhGJEEugzwLKzGwYsBKYiNdnHukfwN1mlg3k4nXJ/DaehUar3bYNyKZHMrpc3rgDXroe9joRzngIdMo7EemC2g1051yjmV0GzMA7bPFB59w8M5vs336Pc+5TM3se+AgI4x3aODeRhdfV15BlCQ5W5+DlG+CN22GfM+C0e7yWuIhIFxTTcejOueeA56Km3RN1/TbgtviV1raGpjqyEjmOSzgMM66Gd++BA86Dk38LWRoETES6rpT9pWhjuI6cUIJ+JRpugmcuhzl/gUMvheNv1M5NEenyUvZYuiZXR252AlrojfXw94u8MB93lcJcRGJSWVnJxIkTGT58OCNHjuTEE09kwYKEHRfSopRsoTc1NeFcLXl5cW6hh8PwxAUw/5/w9Rvg8Mvju3wRSUvOOU477TTOO+88pk+fDsCcOXP44osvKC8vT1odKRnom5evAsJ0K4jzz2dn3uWF+fE3w1cuie+yRSRtvfLKK+Tk5DB58uTt08aMGZP0OlIy0NctWgZA98LC+C10xSx46QYY+Q049H/jt1wRSarKm26i7tP4Dp+bt/cIdv/Zz1q9fe7cuRx44IFxXeeuSMk+9A0rVgPQs6Q4Pgus2Qh/vxB6DYRT7lKfuYikpJRsoW/+Yi0AvQfEYVBH5+DZH8DmVXDB89Ctd+eXKSKBaaslnSijRo3iiSeeSPp6o6VkC71640YAioYOamfOGLw/FT55Go75BQw6qPPLE5GMc8wxx1BXV8d99923fdqsWbN0TtFYbNtSDUBJWSfHTPniE3h+Cgw/Bg7TES0ismvMjKeeeooXXniB4cOHM2rUKK677joGDBiQ1DpSssultnYbkEteYSdOblG/zTtEMa8QTvuThrcVkU4ZMGAAjz32WKA1pGSg1zXUdn4cl+enQNV8mPSUdzooEZEUl5LN0oamOkJZnfiV6Nwn4YOH4Ygfed0tIiJpICUDvcnVkbOrgb5hqXdUS+lBcPTP41qXiEiQUjLQw66O3F0dk/xfVwAGpz+goXBFJK2kXKA31NXhXC35+bswjsu6xbDoRTjsMigaEv/iREQClHKBvsH/2X9Bjx4dv/PsByErGw44N85ViYgEL+UCfe3iFQB079WrY3dsqIE5f4URJ0HP3RNQmYhkslAoxJgxYxg1ahSjR4/mjjvuIBxu54TycZZyhy2ur1gJQGHfDo7jMu9pqNkAYy9KQFUikum6devGnDlzAFizZg3f+ta32LRpE9dff33Saki5Fvri5Z8B8O/qV6jcWhn7HWc/AMVlMOyrCapMRMTTr18/7r33Xu6++26cc0lbb8q10E+95FKee/YR/pH3FC/845tcfcjVnLzHyVhbIySu/ggqZnnjnGskRZG09sZjC1i7ojquyywZ1IMjz+rYiSr22GMPwuEwa9asYbfddotrPa1JuRZ60eABfPvSq/jrWY+zZ9Ge/OzNn/HjV3/M+tr1rd9p9gOQnQ+jJyavUBHJeMlsnUMKttCbDSocxEPHP8TDnzzM3f+9mw/+8QHXfeU6jh589I4z1m6Gjx6HfU6Hgj7BFCsiSdPRlnSiLFmyhFAoRL9+yRtaJOVa6JFCWSEu3OdCpp88nX4F/bj8lcu55s1r2FK/5cuZPvobNGzVzlARSZqqqiomT57MZZdd1nZ3cJylbAs9UnlROY+e+Cj3fHQPD3z8ALMqZ/G7Y37HXkXl3rHn/UfDwAOCLlNE0lhNTQ1jxoyhoaGB7OxsJk2axI9//OOk1pAWgQ6QE8rh+/t/n3Gl4/jRqz9i0r8nceuICzh6zSdwyu+0M1REEqqpqSnoElK7y6Ul+/Xdj2knTWN4r+H8YO4feKC4L26f04MuS0Qk4dIu0AH6FfTjoa/+hvFba7izsBvXzLqZ+qb6oMsSEUmotAx0gPyPHufWNVVcWnYWzyx+hotmXMS6mnVBlyUikjDpGejhMLz/EDbkcCYf9gtuH3c789fP55x/ncNn6z8LujoRkYRIz0Bf/LJ3IouxFwJw3NDjmHrCVJpcE5P+PYmXlr8UbH0iIgmQnoE++wHo3hf2PnX7pFHFo7bvLP3hKz/kzvfvpDHcGGCRIiLxlX6BvnEFLHge9p8E2bk73NSvoB9TT5jKGeVn8MDcB/jeC99jbc3agAoVkXTSPHxu898tt9yS9BrS5jj0Zu7DaeAcn/Q/jfULq9hS28iW2ga21Day2b88rOQ8rj10X26ddRNnPXsWtx91O/v32z/o0kUkhUUOnxuUtAv0pbP+zZbwUE59ZAWwYqfbC3JDbKtvYrfCQs457DZeXHcbFz5/IT868EdMGjkpqT/TFRGJp7QK9HcXrmL0lo9ZXnIaDx43lp75OfTIy6Znfvb2y6EsY+bitdzxnwX8v+c3MKDoEvYo+we3zb6NOVVz+NVhv6JH7i6c3k5EuoRXpt7LmmVL4rrMfkP24OjzL25znuaf/je7+uqrOfvss+NaR3vSJtBrG5p45O9Pc4g1cOjRp5A3ovXxhw8bXsJXJhfz2oIqbv/PAt5/bwK7DdqNF5c9y4INC/jtUb+lrKgsidWLSKpTl0sc3fniQgZvmQM5kLfH4e3Ob2YctVc/xpX35YVPvuCOFwpZu6wfKwZN5+x/nsMvDv0539jzG+qCEUkx7bWk01laHOUyd+Um7ntjCROKlkLfvaF77OcbNTOOG7U7z11+JHdOOIPe66dQs6WUa2dey5WvXc22hm2JK1xEJI5SPtAbm8Jc9feP6FsQorx+Hgw5bJeWk5VlnDJ6ADMuP4nvDL2R+qqv8fzS55jw1Jks3LAwzlWLSLpp7kNv/psyZUrSa4ipy8XMxgN3ASHgfudciwdYmtlBwDvA2c65J+JWZRvue+Nz5q3azKMn5WIvVe9yoDfLzwnx8xNHMX7Utfzg6RGsaprKmc9M5KqDpjBx7zPUBSMiLUqJ4XPNLAT8HjgBGAmcY2YjW5nvVmBGvItszZKqau58cQHjR+3OYSF/jJZOBnqzA4cU8eIl/8OEvrdTt3UwN836FRc/f4W6YESky4qly+VgYJFzbkr6Mc8AAAx/SURBVIlzrh6YDkxoYb7vA38H1sSxvlaFw44pT35MXnYWv5owCpbNhKJhUDggbuvIzwlx06mH8dD4e+lWfSJvf/Eix04/jSc/msdbi9by/rL1zF25iUVrqqnYsI211XXUNgT/KS0imSmWLpeB7PgLnQrgkMgZzGwgcBpwDHBQawsys4uBiwEGDx7c0Vp3MG3Wct77fD2/Pn0/+vXIheUzYcRJnVpmaw7doy+vfvdGrnhmf97ceiNTnn+Y+vXjWpw3O8s4eb/+/M+Re7DPwF4JqUdEduacS6suUedch+8TS6C39AhFr+lO4CrnXFNbD6hz7l7gXoCxY8d2vFpf5aZabnluPofvWcyZY0thzadQswGGtH+44q4qyM3mj2dM5Oi/3cde+zVxfvmh1DY0UdcY9v43hKltbGJJ1VYen72Cp+es4vA9i7n4q8P5allJWr3QRLqa/Px81q1bR3FxcVq815xzrFu3jvz8/A7dL5ZArwAGRVwvBVZFzTMWmO4/kCXAiWbW6Jx7ukPVxMA5xzVPz6UhHObm0/bznrxlb3k3xqn/vC0jistZs20ph+7R+qGRP/p6OdPeW85Db33OeQ++x4jde/LdI/fglNEDyM1O+QOLRLqc0tJSKioqqKqqCrqUuMnPz6e0tLRD94kl0GcBZWY2DFgJTAS+FTmDc25Y82Uzmwr8MxFhDvDcx5W8+OkX/PzEvRlcXOBNXDYTCgdC7yGJWOUOyovKeWf1OzSEG8jJymlxnl7dcpg8bjgXHj6MZz9cxX1vLOGKxz/k1zPm840xA8nPCSWsvjRonIh0QuLeW/E0dkgfjigrifty2w1051yjmV2Gd/RKCHjQOTfPzCb7t98T96racPCwPlxy1HAuOHxoc4FeoA87MilpVlZURmO4kaWblrY7PEBudhanH1jKNw8YyOsL13Lv64v50+vxHWNCRFLP5HHDgwl0AOfcc8BzUdNaDHLn3PmdL6t1fXvmceX4EV9OWL8EqiuT0t0CXgsdYMGGBTGP92JmjCvvy7jyvoksbZd2oohI+kj9sVyWzfT+J3CHaKRhhcPIzsrukr8eTYedQSKy61J/D92ymVBQDCXlSVldTiiHYb2GsWDDgqSsT0QkVmkQ6G963S1JbJ2WF5WzcGPXa6GLSGZL7UDfuAI2Lk9ad0uzst5lVG6tZFPdpqSuV0SkLakd6Mvf9v4naYdos+Ydo12xH11EMldqB/qytyCvF+y2T1JXuz3Q1e0iIl1Iigf6TBh8KGQl98cE/Qr6UZhbqB2jItKlpG6gV1fB2gVJ724B7/DA8qJyBbqIdCmpG+jLk3v8ebSyojIWbVhE2IUDWb+ISLTUDfRlMyGnAPqPDmT15UXlbGvcxsrqlYGsX0QkWgoH+ltQehBk5wayeh3pIiJdTWoGes1GqJwbWHcLwJ699wRQP7qIdBmpGegr3gVcIDtEmxXkFDCo5yAFuoh0GakZ6EvfhKwcKB0baBllvcvU5SIiXUZqBvqymTDwQMjpFmgZ5X3KWb5lObWNtYHWISICqRjoddWweg4MDa7/vFl5UTlhF2bxpsVBlyIikoKBXjELwo2B9p83K+vtneBiwXr1o4tI8FIv0HN7wMgJUHpw0JUwqOcg8kP52jEqIl1C6p2xaNBBMOjPQVcBQCgrxPDewzVIl4h0CanXQu9iyovKdaSLiHQJCvROKi8qZ33tetbWrA26FBHJcAr0Tior8neMqh9dRAKmQO+k5kBXt4uIBE2B3kl98vtQ0q1ELXQRCZwCPQ60Y1REugIFehyUF5WzeONiGsONQZciIhlMgR4HZUVl1IfrWb55edCliEgGU6DHQfPJLtSPLiJBUqDHwR699iBkIQW6iARKgR4HuaFchhYO1Y5REQmUAj1OyovK1UIXkUAp0OOkrKiMVVtXsaV+S9CliEiGUqDHSfOO0UUbFwVciYhkKgV6nDQHuvrRRSQoCvQ42b377vTM6al+dBEJjAI9TsyMsqIyBbqIBEaBHkdlRWUs3LAQ51zQpYhIBlKgx1F5UTnVDdUs2bQk6FJEJAMp0OPoyIFH0iOnB1e/cTW1jbVBlyMiGSamQDez8Wb2mZktMrMpLdz+bTP7yP+baWaj419q19e/R39uOuImPl3/KTe+e6O6XkQkqdoNdDMLAb8HTgBGAueY2cio2T4Hxjnn9gNuAO6Nd6Gp4ujBR3Pxfhfz9KKneXzB40GXIyIZJJYW+sHAIufcEudcPTAdmBA5g3NupnNug3/1HaA0vmWmlktGX8LhAw7n5vdu5sOqD4MuR0QyRCyBPhBYEXG9wp/WmouAf3emqFQXygpx61dvZbeC3fjxqz9mbc3aoEsSkQwQS6BbC9Na7Bw2s6PxAv2qVm6/2Mxmm9nsqqqq2KtMQb3yevHbo37LprpNXPn6lTqbkYgkXCyBXgEMirheCqyKnsnM9gPuByY459a1tCDn3L3OubHOubF9+/bdlXpTyt7Fe/PLr/ySWZWzuPP9O4MuR0TSXCyBPgsoM7NhZpYLTASeiZzBzAYDTwKTnHP6qWSEU4afwsS9JvLwJw/z/NLngy5HRNJYu4HunGsELgNmAJ8Cjznn5pnZZDOb7M92LVAM/MHM5pjZ7IRVnIKuPOhKxvQdw7VvXcuiDRqNUUQSw4I6Vnrs2LFu9uzMyf0129Zw1rNnUdNYQ6+8Xp1aVpZlbf9v2E7XRaRrO63sNM4bdd4u3dfM3nfOjW3ptuxOVSUx61fQjz9+7Y9M/2w6TeGmXV6O8/dHh10YhyPswuAgTNi7LCJdXnG34oQsV4GeRHsX7831h10fdBkikqY0louISJpQoIuIpAkFuohImlCgi4ikCQW6iEiaUKCLiKQJBbqISJpQoIuIpInAfvpvZlXAsl28ewmQqYOMZ+q2Z+p2Q+Zue6ZuN7S97UOccy0OVxtYoHeGmc1ubSyDdJep256p2w2Zu+2Zut2w69uuLhcRkTShQBcRSROpGuj3Bl1AgDJ12zN1uyFztz1Ttxt2cdtTsg9dRER2lqotdBERiaJAFxFJEykX6GY23sw+M7NFZjYl6HqSwcwGmdkrZvapmc0zsx8EXVOymVnIzP5rZv8MupZkMbPeZvaEmc33n/uvBF1TspjZj/zX+lwzm2Zm+UHXlAhm9qCZrTGzuRHT+pjZC2a20P9fFOvyUirQzSwE/B44ARgJnGNmI4OtKikagSucc3sDhwKXZsh2R/oB3knKM8ldwPPOuRHAaDJk+81sIHA5MNY5tw8QAiYGW1XCTAXGR02bArzknCsDXvKvxySlAh04GFjknFvinKsHpgMTAq4p4Zxzq51zH/iXt+C9sQcGW1XymFkpcBJwf9C1JIuZFQJfBR4AcM7VO+c2BltVUmUD3cwsGygAVgVcT0I4514H1kdNngA87F9+GPhGrMtLtUAfCKyIuF5BBgUbgJkNBfYH3g22kqS6E7gSyKSzYO8BVAEP+V1N95tZ96CLSgbn3ErgN8ByYDWwyTn3n2CrSqrdnHOrwWvMAf1ivWOqBbq1MC1jjrs0sx7A34EfOuc2B11PMpjZycAa59z7QdeSZNnAAcAfnXP7A1vpwFfvVOb3GU8AhgEDgO5m9p1gq0oNqRboFcCgiOulpOlXsWhmloMX5n91zj0ZdD1JdDhwqpktxetiO8bM/hJsSUlRAVQ455q/iT2BF/CZ4GvA5865KudcA/AkcFjANSXTF2bWH8D/vybWO6ZaoM8CysxsmJnl4u0oeSbgmhLOzAyvL/VT59wdQdeTTM65q51zpc65oXjP98vOubRvrTnnKoEVZraXP+lY4JMAS0qm5cChZlbgv/aPJUN2CPueAc7zL58H/CPWO2YnpJwEcc41mtllwAy8Pd8POufmBVxWMhwOTAI+NrM5/rSfOeeeC7AmSbzvA3/1Gy9LgAsCricpnHPvmtkTwAd4R3j9lzQdBsDMpgFHASVmVgH8ErgFeMzMLsL7cDsz5uXpp/8iIukh1bpcRESkFQp0EZE0oUAXEUkTCnQRkTShQBcRSRMKdBGRNKFAFxFJEwp0kQhmtq+ZLTOz/w26FpGOUqCLRHDOfYw3xMC5Qdci0lEKdJGdrQFGBV2ESEcp0EV2dguQZ2ZDgi5EpCMU6CIRzGw80B34F2qlS4pRoIv4/BMR/xq4BPgY2CfYikQ6RoEu8qVrgD8755aiQJcUpEAXAfwTSXwd7/yloECXFKTx0EVE0oRa6CIiaUKBLiKSJhToIiJpQoEuIpImFOgiImlCgS4ikiYU6CIiaeL/A7OyFUf1i/GhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cost='h'\n",
    "lSeq = np.arange(0,10,.25)\n",
    "QRE = np.zeros((len(lSeq),6))\n",
    "p0 = np.log([0.25,0.25,0.5,0.5,0.5,0.5])\n",
    "for (i,l) in enumerate(lSeq):\n",
    "    res = minimize(F, p0, args=(l,cost,), method='SLSQP',constraints=cons)\n",
    "    p0 = res.x # use the new value as the new guess for next lambda\n",
    "    QRE[i] =  np.exp(res.x)\n",
    "    \n",
    "#Plot\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(lSeq, QRE[:,0],lSeq, QRE[:,1],lSeq, QRE[:,2],lSeq, QRE[:,3],lSeq, QRE[:,4],lSeq, QRE[:,5])\n",
    "ax.set_xlabel('$\\lambda$')\n",
    "handles, labels = ax.get_legend_handles_labels()\n",
    "lgd = ax.legend(handles, labels = ['A','A_1','B', 'C','D','E'], loc = 0)\n",
    "ax.set_title(\"QRE when cost=\"+str(cost)+\" Prob of decision at turn 0\")\n",
    "plt.show()    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "     fun: nan\n",
       "     jac: array([nan, nan, nan, nan, nan, nan])\n",
       " message: 'Inequality constraints incompatible'\n",
       "    nfev: 7\n",
       "     nit: 1\n",
       "    njev: 1\n",
       "  status: 4\n",
       " success: False\n",
       "       x: array([-0.97515601, -0.47620391, -0.97424505,  0.02006188,  0.01977326,\n",
       "        0.0198609 ])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.30675987, 0.12251035, 0.0412845 , 0.01373314, 0.00475128,\n",
       "       0.00173234, 0.00173234, 0.00173234, 0.00173234, 0.00173234,\n",
       "       0.00173234, 0.00173225, 0.00172967, 0.00172967, 0.00172967,\n",
       "       0.00172967, 0.00172967, 0.00172967, 0.00172967, 0.00172967,\n",
       "       0.00172967, 0.00172967, 0.00172967, 0.00172967, 0.00172967,\n",
       "       0.00172967, 0.00172967, 0.00172967, 0.00172967, 0.00172967,\n",
       "       0.00172967, 0.00172967, 0.00172967, 0.00172967, 0.00172967,\n",
       "       0.00172967, 0.00172967, 0.00172967, 0.00172967, 0.00172967])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1- QRE[:,0] - QRE[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "h\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([0.49999   , 0.47778543, 0.450728  , 0.42560805, 0.40133081,\n",
       "       0.37756595, 0.37756595, 0.37756595, 0.37756595, 0.37756595,\n",
       "       0.37756595, 0.37756027, 0.37747722, 0.37747722, 0.37747722,\n",
       "       0.37747722, 0.37747722, 0.37747722, 0.37747722, 0.37747722,\n",
       "       0.37747722, 0.37747722, 0.37747722, 0.37747722, 0.37747722,\n",
       "       0.37747722, 0.37747722, 0.37747722, 0.37747722, 0.37747722,\n",
       "       0.37747722, 0.37747722, 0.37747722, 0.37747722, 0.37747722,\n",
       "       0.37747722, 0.37747722, 0.37747722, 0.37747722, 0.37747722])"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(cost)\n",
    "QRE[:,2] \n",
    "# compare with the matrix one, not comparable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-14-685574865e82>:5: UserWarning: You have mixed positional and keyword arguments, some input may be discarded.\n",
      "  lgd = ax.legend(handles, labels = ['C'], loc = 0)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEYCAYAAABMVQ1yAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAfXElEQVR4nO3de5xcdX3/8dc7m82F3JAkQq4kYEASbso+QPSnpl7agCJqtQ3eUYtRUWu9FK1W1NYHVX9tbaVSVIp4gaLCT1QUfCgEb60kgrLLzYBANpuEJMgmu3E32d3P749zdnMymdmd3Ux2Mue8n4/HPHbOzDlnPmfm7Hu+850z56uIwMzMGt+EehdgZma14UA3M8sJB7qZWU440M3McsKBbmaWEw50M7OccKCbmeWEA/0wIykkPa3eddiBJD0i6UU1WtcrJG2U1CXpGVXMf7uktx7kY35Y0peqmO8Hkt54MI9l9VHYQJf0Jkn3SNotaYuk/5A0K3P/pZL2pv9wT0r6haSzM/evlDSQ3p+9nF3+EfMtfT5/NsZlJemfJO1IL5+WpArzLknf9LLP+UeHWfftknrS+bZLukHSvLHUWWOfBS6OiOkRcdd4PGBEfCoiRnxTiIhzIuIr41HToMzrOnGYeS6V9LVxrusoSTdK6pb0qKTXjOfjj1YhA13S+4B/Aj4AzAKeBSwBbpXUnJn1vyNiOjAHuA34ZsmqOtJ/yOzll4d+C3LnIuDlwGnAqcBLgbeNsMyRmef8kyPMe3H6Op4AHAn8S+kMwwXJIXIs0DbOj1loY3yNLwf2AEcDrwW+IGlFTQurpYgo1AWYCXQBf1Fy+3TgceCN6fSlwNcy9y8HApibTq8E2qt8zAuB72amNwDXZ6Y3Aqen1wNYA/wO+APJDqXMvG8G7kvvuwU4NnPfsMuW1NQEfBh4CNgFrAcWpfc9G7gT6Ez/Pjuz3JuAh9Nlfk+yk58E9AD96XP75Chfk18AF2Wm3wL8T4V5l6TbObHKdd8OvDUz/U6gNb3+CPC3wG+BXmAi8DKSoH0yXfakzLKPAB8C7k2f3/8CplR43AnAR4BH0/3qGpLGw+T0OQqgG3iowvIvBu5PX4PPA2tLtmO4/WAF8CPgCWAr8OHSfRqYAnwN2JFu653A0aXPWaXtKHkt3gg8BmwH/m6Y1+IlwF3ATpJ9/tLMfY+l6+pKL2eXLLuKJFj3pvf/JvOavCgzX3YbB+t7S7r+O0j235+RfEL6A8k+fE6Feqelj3lC5ravApfVMpNqeal7AeO+wcmO0UeZQAC+Any9zI4xCbgs3WEnpretpPpAPy79p5kAzEv/OTZl7vsDMCGdDuB7JC3JxcA2YFV638tJ3gxOIgmfjwC/yDxOxWXL1PQB4B7gREAkrePZwFFpPa9PH+OCdHp2uoPvBE5M1zEPWJFefxPws5LHuCTd7rKXzHydwFmZ6RZgV4W6B/9JNwHtJKE6Z5jn/nb2hdMc4CfAV9PpR4C7gUXAVJIWfDdJmDYDH0yf70mZ+VvT+Y8Cfg78Q4XHfXO67HEkjYUbBh8381o9rcKyc9Ln+VVpHe8l2WcHt6PifgDMADYD7yMJ7RmDzy3779NvA74LHEHy5n4GMLPMc1ZxOzKvxRfT5+80kjfGkyps10rgFJL/g1NJ3mxeXrKuim/UlDSyMq/JSIF+Dcm+O5VkP90L/FW63W8HOijT8AGeAfyx5Lb3k2mcHW6Xuhcw7hsMrwO2VLjvMuDWzI6xhyR8+klaMitLds4BDgyqaRXWvRF4JrAauBL4FfB0ktb7TZn5Avg/menrgUvS6z8A3pK5bwKwm7R1NtyyZep5ADi/zO2vB35Vctsv03+Eaek2/jkwtWSeN1ES6KN4TfqBp2eml6XbUu6fbDpJ4E8k+Rj8LeCWYdZ9e/ocPUnyJvB19n3KegR4c2bej7L/J6cJ6TIrM/Ovydx/LpVb2D8G3pGZPjENksEGwXCB/gYyn1BI3nDb2ReyFfcDkjfguyqs91L2hd2bST4ZnVrhOXvrSNvBvsBcmLn/V8DqKl/3fwX+Jb0+uK5DEejHleynGzLTR6TzHFPm8Z5LSVaQvBHcPpb9fDwuRexD3w7MqdCfNo+kVTvo+og4kiQ4WklaMVkdEXFkyaW7wuOuJXkTeF56/Xbg+ellbcm8WzLXd5OEGCT/sJ9Lv6R9kuQjtYAFVSxbahFJd0up+SSfILIeBRak2/aXJN06myV9X9LTK6x/NLpIusIGzQS6Iv0PyoqIrohYFxF9EbEVuBj4U0kzS+fNeHf62iyIiNdGRPY13pi5vt+2R8RAev+CCvM/mi5TTunz+Cj73oRGMj/7OOnzkH3c4faDSq9rqa+SdNVcJ6kj/SK6ucx81WxHVfucpLMk3SZpm6ROkv1oThW1HqyNJdND9UbE7vRquZpL90vS6V21K622ihjovyT5WPjK7I2SpgHncGC4EhHbST6iXnoQR0gMBvpz0+trqRzolWwE3lbyBjI1In4xhno2AseXub2DJDCyFpO0VImIWyLixSRvfveTfNyGpJWzn/QwudKjgIYumVnbSD6uDzqN6r8wHHzcskfFjGJ5KNn29EibRaTbnlqUub44Xaac0udxMUm3ydYqatqcfZxMHYOG2w8qva77iYi9EfHxiFhO8p3JS0k+GdRyO0p9A7iJ5LuaWcAV7HvdDth/ypVd5rZuklb2oGOqXK4aDwITJS3L3DaafXPcFS7QI6IT+Djw75JWSWqWtITkCJbtJB/Jyy13P0mL5oNjfOi1wJ+QdFW0Az8l6c+fTfJFUTWuAD40+C27pFmSXj3Ger4EfFLSsvSwwVMlzQZuBk6Q9BpJEyX9JckXwt+TdLSkl6Vvfr0kLZj+dH1bgYWSJg0+QCSHyZUeBTR0ydRyDfA3khZImk/S/3t1uaLTVt6Jkiak9f4byUfgzjE+D1nXAy+R9MK0tfq+dDuzb5jvlLRQ0lEkXyr/d4V1XQu8V9JSSdOBT5EcNdVXRR3fB1ZIemX6SfLd7B9Uw+0H3wOOkfTXkiZLmiHprNIHkPQnkk6R1ETSX7+Xfa9lrbaj1AzgiYjokXQmkD0EcBtJF+Zxwyy/FVgiKZtbdwOr0//jFpLvHWoi/UR6A/AJSdMkPQc4n+TTzWGpcIEOEBGfJvln/Cz7jtY4gqQvrlKXCcBngIskPTWdnl+m5fnnFR7zQZIA/Gk6vZPkaJGfR0S5f6Ry67iR5HDL6yTtJOkGOqeaZcv4Z5IAu5XkH/rLJG82O0haa+8j+d7gg8BL008pE9LbO0g+5j8feEe6vp+QtFy2SNo+ylr+k+QLunvSbfp+ehsAktokvTadPA74Icnr1koSuBeM8vHKiogHSL5j+XeSN/fzgPMiYk9mtm+QPGcPp5d/qLC6q0j+8e8g2b96gHdVWcd24NUk3+nsIPlO4eeZ+yvuBxGxi+RL3fNIuhZ+R9KQKHUMyfcPO0mOlllLctRLzbajjHeQhOMu4O9J9r/BbdoN/CPw87Qr6Vlllh88bHiHpF+n1z9K8onkDyQNtW+Msbbhap5KcoTPtcDbI+KwbaGrTDdl4Uh6M8nO8JyIeKze9ZiZjYUDPSXp9cDeiLiu3rWYmY2FA93MLCcK2YduZpZH433+iiFz5syJJUuW1Ovhzcwa0vr167dHxNxy99Ut0JcsWcK6devq9fBmZg1JUukP/4a4y8XMLCcc6GZmOeFANzPLCQe6mVlOONDNzHJixECXdJWkxyW1Vrhfkv5N0gZJv5X0zNqXaWZmI6mmhX41yVkBKzmH5ORBy0jGhvzCwZdlZmajNeJx6BFxR3p62UrOB65JT8L/P5KOlDQvIjbXqMa669nbz9adPWzp7GHrrl76BwYqzhsB/QPBQAT9AzAQg9eTS0RyW6TzDqSnXhgYiDGftLlWDuYsEIPVD65jaFU+tYTZAVqWHMXzTij726CDUosfFi1g/xFB2tPbDgh0SReRtOJZvHhxDR66tiKCb65v5+6NT7Kls4fNnT1s3dnDE917Rl7YKtJYh54wy6k1zz/+sA30cv+uZZtlEXElyXiatLS0HFZNt109e/nAN3/LD9u28JQjmpk3ayrzZ03hmYuP5JiZUzhm1hTmzZrKU2dOZlLT8D1VTRPEhAmiSWLCBJigwetigpJppX9h37QA1Tn9DubRB0uv9zaYFVUtAr2d/YfHWkjlYbkOSxse38XbvrqeR3bs5u/OPYm3PnepQ8nMGk4tAv0m4GJJ1wFnAZ2N1H/+g3s28/5v/oYpzU187S1ncfbxs+tdkpnZmIwY6JKuJRnceI6kduBjQDNARFxBMgblucAGkhG/LzxUxdZSX/8An7n1Af5z7cOcvuhIvvC6ZzJv1tR6l2VmNmbVHOUy7HiN6dEt76xZReNgR1cv77r2Ln7x0A5ec9ZiPnbeciZPbKp3WWZmB6Vup8+tl9ZNnVx0zTq2d+/h0686lb9oWTTyQmZmDaBwgf6J797L3oHg22uezSkLZ9W7HDOzminUuVz6B4LWjk5ecso8h7mZ5U6hAv3327vZvaefFfNn1rsUM7OaK1Sgt3V0AnDyArfOzSx/ChXorZs6mTRxAk976vR6l2JmVnMFC/SdnHTMDJpH+Om+mVkjKkyyRSRfiK5wd4uZ5VRhAn3jE39kV08fJ893oJtZPhUm0FuHvhD1ES5mlk/FCfRNnUycIE44eka9SzEzOySKE+gdO1l29AymNPucLWaWT4UI9IigbVMnJ/sHRWaWY4UI9C07e9jRvcc/KDKzXCtEoLdu2gngn/ybWa4VItDbOjqR4KR5DnQzy69CBHrrpp0cN2ca0yYX7mzBZlYghQj0to5O95+bWe7lPtC3d/WyubPHvxA1s9zLfaC3daRfiPoXomaWc7kP9NZNyU/+V7iFbmY5l/tAb+voZPFRRzBranO9SzEzO6RyH+itm3b6hFxmVgi5DvTO3Xt57Ind7m4xs0LIdaC3bfYYomZWHPkOdP/k38wKJNeB3trRybxZU5gzfXK9SzEzO+TyHeibOt1/bmaFkdtA7+7t4+Ht3T7CxcwKI7eBft/mnUTgn/ybWWHkNtCHfiHqFrqZFURVgS5plaQHJG2QdEmZ+58i6UZJv5X0K0kn177U0Wnr2MnsaZM4ZuaUepdiZjYuRgx0SU3A5cA5wHLgAknLS2b7MHB3RJwKvAH4XK0LHa3Wjp2sWDALSfUuxcxsXFTTQj8T2BARD0fEHuA64PySeZYDPwaIiPuBJZKOrmmlo9Czt5/fbd3lQaHNrFCqCfQFwMbMdHt6W9ZvgFcCSDoTOBZYWLoiSRdJWidp3bZt28ZWcRUe3LqLvoHwL0TNrFCqCfRyfRZRMn0Z8BRJdwPvAu4C+g5YKOLKiGiJiJa5c+eOuthqDQ4K7SNczKxIqhlksx1YlJleCHRkZ4iIncCFAEo6rX+fXuqitaOTGVMmsuioqfUqwcxs3FXTQr8TWCZpqaRJwGrgpuwMko5M7wN4K3BHGvJ10bapk5Pn+wtRMyuWEQM9IvqAi4FbgPuA6yOiTdIaSWvS2U4C2iTdT3I0zHsOVcEj2ds/wH1bdvkXomZWONV0uRARNwM3l9x2Reb6L4FltS1tbDY83sWevgF/IWpmhZO7X4p6DFEzK6rcBfpD27ppbhJL50yrdylmZuMqd4G+q2cvM6c00zTBX4iaWbHkLtC7evuYNrmqrwbMzHIld4He7UA3s4LKXaB39fYxw4FuZgWUu0Dv7u1n2uSmepdhZjbuchfo7kM3s6LKZaBPd6CbWQHlLtC7HehmVlC5CvT+gWD3nn53uZhZIeUq0Lv3JKdgdwvdzIooX4HemwS6W+hmVkS5DPTpUxzoZlY8uQr0XT2DXS4+Dt3MiidXgd7d2w/AtEluoZtZ8eQq0Lvc5WJmBZarQB/qQ/eXomZWQLkK9C4f5WJmBZbLQHcL3cyKKFeB3t3bx8QJYvLEXG2WmVlVcpV8g4NbSB5+zsyKJ1eBvssn5jKzAstVoCctdP+oyMyKKWeB3u8WupkVVq4CfZdHKzKzAstVoHtwCzMrMge6mVlO5CrQPUC0mRVZbgI9ItxCN7NCqyrQJa2S9ICkDZIuKXP/LEnflfQbSW2SLqx9qcP7495+BsLncTGz4hox0CU1AZcD5wDLgQskLS+Z7Z3AvRFxGrAS+L+SJtW41mH51LlmVnTVtNDPBDZExMMRsQe4Dji/ZJ4AZij5zf104Amgr6aVjmBwcAuPVmRmRVVNoC8ANmam29Pbsj4PnAR0APcA74mIgdIVSbpI0jpJ67Zt2zbGksvrSoef82hFZlZU1QR6uTNdRcn0nwF3A/OB04HPS5p5wEIRV0ZES0S0zJ07d9TFDsenzjWzoqsm0NuBRZnphSQt8awLgRsisQH4PfD02pRYnW73oZtZwVUT6HcCyyQtTb/oXA3cVDLPY8ALASQdDZwIPFzLQkfi0YrMrOhGTL+I6JN0MXAL0ARcFRFtktak918BfBK4WtI9JF00fxsR2w9h3Qdwl4uZFV1V6RcRNwM3l9x2ReZ6B/CntS1tdDxAtJkVXW5+Kdrd24cER0zyYYtmVky5CfRdvX1Mm+Th58ysuHIT6B6tyMyKLkeB7tGKzKzYchPoXT7TopkVXK4C3cegm1mR5SbQux3oZlZwuQn0rt4+ZjjQzazAchXobqGbWZHlJtDd5WJmRZeLQO/t62dvfzDDZ1o0swLLRaAPjlY0zT/7N7MCy0WgD41W5C4XMyuwfAS6z7RoZpaPQO/e49GKzMxyEegercjMLC+B3uMuFzOzXAR6t1voZmb5CHR/KWpmlrNA93HoZlZkuQj07t4+pjRPYGJTLjbHzGxMcpGAXR6tyMwsH4He7dGKzMzyEeg+da6ZmQPdzCw3chHo3R6tyMwsP4HuFrqZFV0uAt1dLmZmOQr06ZP9oyIzK7aGD/S+/gF69g4wfXJzvUsxM6urhg/0oeHn3EI3s4KrKtAlrZL0gKQNki4pc/8HJN2dXlol9Us6qvblHqhrj0/MZWYGVQS6pCbgcuAcYDlwgaTl2Xki4jMRcXpEnA58CFgbEU8cioJL+dS5ZmaJalroZwIbIuLhiNgDXAecP8z8FwDX1qK4agydOtfDz5lZwVUT6AuAjZnp9vS2A0g6AlgFfLvC/RdJWidp3bZt20Zba1kercjMLFFNoKvMbVFh3vOAn1fqbomIKyOiJSJa5s6dW22NwxrqcpnkQDezYqsm0NuBRZnphUBHhXlXM47dLbCvy2WGu1zMrOCqCfQ7gWWSlkqaRBLaN5XOJGkW8HzgO7UtcXj+UtTMLDFiCkZEn6SLgVuAJuCqiGiTtCa9/4p01lcAt0ZE9yGrtoyh4ed8HLqZFVxVzdqIuBm4ueS2K0qmrwaurlVh1erq7ae5SUye6EA3s2LLwS9FPVqRmRnkINB9pkUzs0QuAt0tdDOzHAS6B7cwM0vkItDdQjczy0Gg73Kgm5kBOQj0pMvFhyyameUg0Ps9WpGZGQ0e6AMDQfcejydqZgYNHui79/YT4fO4mJlBgwe6T8xlZrZPQwe6T51rZrZPYwd6jwe3MDMb1NCB7i4XM7N9GjrQhwaIdqCbmTV2oHfvSQPdfehmZo0d6EN96D4O3cyswQO9tx9wl4uZGTR4oHf39jFBMLXZLXQzs4YO9MHRiiTVuxQzs7pr+EB3d4uZWaKhA92jFZmZ7dPQge4WupnZPg50M7OcaOhA92hFZmb7NHig97sP3cws1dCB3tXbxwwHupkZ0MCBHhFDx6GbmVkDB3pv3wD9A+FANzNLNWyge7QiM7P9VRXoklZJekDSBkmXVJhnpaS7JbVJWlvbMg80NLiFRysyMwNgxDSU1ARcDrwYaAfulHRTRNybmedI4D+AVRHxmKSnHqqCB+3q8WhFZmZZ1bTQzwQ2RMTDEbEHuA44v2Se1wA3RMRjABHxeG3LPFC3RysyM9tPNYG+ANiYmW5Pb8s6AXiKpNslrZf0hnIrknSRpHWS1m3btm1sFac8WpGZ2f6qCfRy56aNkumJwBnAS4A/Az4q6YQDFoq4MiJaIqJl7ty5oy42a7DLZbp/KWpmBlTRh07SIl+UmV4IdJSZZ3tEdAPdku4ATgMerEmVZXSnoxW5D93MLFFNC/1OYJmkpZImAauBm0rm+Q7wXEkTJR0BnAXcV9tS9zd0lIsD3cwMqKKFHhF9ki4GbgGagKsiok3SmvT+KyLiPkk/BH4LDABfiojWQ1l4lw9bNDPbT1VpGBE3AzeX3HZFyfRngM/UrrThdfX2ccSkJpomePg5MzNo4F+KerQiM7P9NWyge3ALM7P9NWygdzvQzcz207CB3uXRiszM9tPAgd7vFrqZWUbDBrq7XMzM9tewge7RiszM9tfQge4WupnZPg0Z6Hv7B9jTN+AWuplZRkMGus+FbmZ2oIYM9H2nznWgm5kNashAHxzcwl0uZmb7NGagD5061z8sMjMb1JCB3pUObjHDw8+ZmQ1pzEDvcZeLmVmphgz0bg9uYWZ2gIYM9MHRitzlYma2T0MHurtczMz2achA7+7tY9LECTQ3NWT5ZmaHREMmos/jYmZ2oIYMdJ8618zsQA2Zij51rpmN1d69e2lvb6enp6fepQxrypQpLFy4kObm5qqXachUTLpc/CtRMxu99vZ2ZsyYwZIlS5BU73LKigh27NhBe3s7S5curXq5Bu1y6XcL3czGpKenh9mzZx+2YQ4gidmzZ4/6U0SDBrr70M1s7A7nMB80lhobMtB3OdDNzA7QkIHe7S9FzazBbdmyhdWrV3P88cezfPlyzj33XB588MGDWmfDBXr/QLB7T79b6GbWsCKCV7ziFaxcuZKHHnqIe++9l0996lNs3br1oNbbcKk4OLiFA93MDtbHv9vGvR07a7rO5fNn8rHzVgw7z2233UZzczNr1qwZuu30008/6MduuBZ6t8/jYmYNrrW1lTPOOKPm6224VPRoRWZWKyO1pBtNVS10SaskPSBpg6RLyty/UlKnpLvTy9/XvtSERysys0a3YsUK1q9fX/P1jhjokpqAy4FzgOXABZKWl5n1pxFxenr5RI3rHDI0WpEHtzCzBvWCF7yA3t5evvjFLw7dduedd7J27dqDWm81LfQzgQ0R8XBE7AGuA84/qEc9CD4Xupk1OknceOON/OhHP+L4449nxYoVXHrppcyfP/+g1ltNKi4ANmam24Gzysx3tqTfAB3A+yOirXQGSRcBFwEsXrx49NUCc2dM4txTjmHO9MljWt7M7HAwf/58rr/++pqus5pAL/f70yiZ/jVwbER0SToX+H/AsgMWirgSuBKgpaWldB1VOePYozjj2KPGsqiZWa5V0+XSDizKTC8kaYUPiYidEdGVXr8ZaJY0p2ZVmpnZiKoJ9DuBZZKWSpoErAZuys4g6RilZ5KRdGa63h21LtbMrBYixtRBMK7GUuOIXS4R0SfpYuAWoAm4KiLaJK1J778CeBXwdkl9wB+B1dEIz5iZFc6UKVPYsWPHYX0K3cHzoU+ZMmVUy6leudvS0hLr1q2ry2ObWXE1+ohFktZHREu5ZXzsn5kVSnNz86hGAWokDXcuFzMzK8+BbmaWEw50M7OcqNuXopK2AY+OcfE5wPYaltNIirrtRd1uKO62F3W7YfhtPzYi5pa7o26BfjAkrav0LW/eFXXbi7rdUNxtL+p2w9i33V0uZmY54UA3M8uJRg30K+tdQB0VdduLut1Q3G0v6nbDGLe9IfvQzczsQI3aQjczsxIOdDOznGi4QB9pwOo8krRI0m2S7pPUJuk99a5pvElqknSXpO/Vu5bxIulISd+SdH/62p9d75rGi6T3pvt6q6RrJY3utIMNQtJVkh6X1Jq57ShJP5L0u/TvU6pdX0MF+igGrM6bPuB9EXES8CzgnQXZ7qz3APfVu4hx9jnghxHxdOA0CrL9khYA7wZaIuJkktN2r65vVYfM1cCqktsuAX4cEcuAH6fTVWmoQOcwG7B6vETE5oj4dXp9F8k/9oL6VjV+JC0EXgJ8qd61jBdJM4HnAV8GiIg9EfFkfasaVxOBqZImAkdQMkpaXkTEHcATJTefD3wlvf4V4OXVrq/RAr3cgNWFCTYASUuAZwD/W99KxtW/Ah8EBupdyDg6DtgG/Ffa1fQlSdPqXdR4iIhNwGeBx4DNQGdE3FrfqsbV0RGxGZLGHPDUahdstECvZsDq3JI0Hfg28NcRsbPe9YwHSS8FHo+I9fWuZZxNBJ4JfCEingF0M4qP3o0s7TM+H1gKzAemSXpdfatqDI0W6CMOWJ1XkppJwvzrEXFDvesZR88BXibpEZIuthdI+lp9SxoX7UB7RAx+EvsWScAXwYuA30fEtojYC9wAPLvONY2nrZLmAaR/H692wUYL9BEHrM6jdADuLwP3RcQ/17ue8RQRH4qIhRGxhOT1/klE5L61FhFbgI2STkxveiFwbx1LGk+PAc+SdES677+QgnwhnLoJeGN6/Y3Ad6pdsKGGoKs0YHWdyxoPzwFeD9wj6e70tg9HxM11rMkOvXcBX08bLw8DF9a5nnEREf8r6VvAr0mO8LqLnJ4GQNK1wEpgjqR24GPAZcD1kt5C8ub26qrX55/+m5nlQ6N1uZiZWQUOdDOznHCgm5nlhAPdzCwnHOhmZjnhQDczywkHuplZTjjQzTIknSLpUUlvr3ctZqPlQDfLiIh7SE4x8IZ612I2Wg50swM9DqyodxFmo+VANzvQZcBkScfWuxCz0XCgm2VIWgVMA76PW+nWYBzoZql0IOJPA+8A7gFOrm9FZqPjQDfb5yPANRHxCA50a0AOdDMgHUjixSTjl4ID3RqQz4duZpYTbqGbmeWEA93MLCcc6GZmOeFANzPLCQe6mVlOONDNzHLCgW5mlhP/H7eMhcDuM/LnAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(lSeq, QRE[:,3])\n",
    "ax.set_xlabel('$\\lambda$')\n",
    "handles, labels = ax.get_legend_handles_labels()\n",
    "lgd = ax.legend(handles, labels = ['C'], loc = 0)\n",
    "ax.set_title(\"QRE when cost=\"+str(cost)+\" Prob of decision at turn 0\")\n",
    "plt.show()  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--------------------------------\n",
      "Difference: 1\n",
      "Guess: [0.5 0.5]\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'Utility' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-12-4da642b6bd15>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      9\u001b[0m     \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"Difference:\"\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mdiff\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     10\u001b[0m     \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"Guess:\"\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mp0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 11\u001b[1;33m     \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"Expected Utility:\"\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mUtility\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mp0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcost\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mgamma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     12\u001b[0m     \u001b[0mp1\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ml\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mUtility\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mp0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcost\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mgamma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mexp\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ml\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mUtility\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mp0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcost\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mgamma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     13\u001b[0m     \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"New Policy:\"\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mp1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'Utility' is not defined"
     ]
    }
   ],
   "source": [
    "#Define the function whose 0 we want to find (PolicyIteration)\n",
    "l=2\n",
    "cost=.5\n",
    "gamma=.5\n",
    "p0=np.array([.5,.5])\n",
    "diff=1\n",
    "while diff>.0001:\n",
    "    print(\"--------------------------------\")\n",
    "    print(\"Difference:\",diff)\n",
    "    print(\"Guess:\",p0)\n",
    "    print(\"Expected Utility:\",Utility(p0,cost,gamma))\n",
    "    p1=np.exp(l*Utility(p0,cost,gamma))/np.sum(np.exp(l*Utility(p0,cost,gamma)))\n",
    "    print(\"New Policy:\",p1)\n",
    "    diff=np.sum(np.abs(p0-p1))\n",
    "    p0=.9*p0+.1*p1\n",
    "print(\"QRE: \",p1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": 2
}