BA_parallel.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parallel case"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns\n",
    "import _lib.pr_func as pr\n",
    "from _lib.utility import *\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "sns.set()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dimensions\n",
    "N,M,K = 20,5,20\n",
    "pr.set_dims([('w',N),('x',M),('a',K)])\n",
    "\n",
    "# environment\n",
    "U = pr.func(vars=['w','a'],val=gauss_utility(N,K,sigma=2.0))\n",
    "pw = pr.func(vars=['w'],val='unif').normalize()\n",
    "\n",
    "# Lagrange multipliers\n",
    "beta1  = 100.0                               # first step, p(x) -> p(x|w)\n",
    "beta2 = 10.0*pr.func(vars=['x'],val='unif') # second step, p(a|x) -> p(a|x,w);\n",
    "# different beta2s for each agent labelled by 'x' can be set by \n",
    "# beta2 = pr.func(vars=['x'], val=np.array([...])), or beta2.val = np.array([...]) (must have a length of M) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Blahut-Arimoto iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initializations\n",
    "px = pr.func(vars=['x'], val='unif').normalize()\n",
    "pagx = pr.func(vars=['a','x'], val='unif').normalize(['a'])\n",
    "F = pr.func('f(w,x)', val='rnd')\n",
    "pagxw_temp = 0\n",
    "\n",
    "# BA iteration\n",
    "for i in range(0,10000):\n",
    "    pxgw = (px*pr.exp(beta1*F)).normalize(['x'])\n",
    "    px = pr.sum(pw*pxgw,['x'])\n",
    "    pagxw = (pagx*pr.exp(beta2*U)).normalize(['a'])\n",
    "    pagx = pr.sum(['w'],pw*pxgw*pagxw)/px\n",
    "    F = pr.sum(['a'],pagxw*(U-pr.log(pagxw/pagx)/beta2))\n",
    "    \n",
    "    # stopping condition\n",
    "    if np.linalg.norm(pagxw.val-pagxw_temp)<1e-10: break\n",
    "    pagxw_temp = pagxw.val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+IAAAKbCAYAAABxQDD2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAWJQAAFiUBSVIk8AAAIABJREFUeJzs3XuU1XW9+P/XXBhAYDCUNEdwMJlBIRRZmeaxtJVpmkEs\nwQEzlDTIHKWTl256VGp5upkK4a2j0yAeRPSQdjmxKkE92QHNNKGZKYRFTEGkkQwNc2Hv3x8t5sSP\n4bJhz3vj18djLf/53Pbr41KY53w++/Mpymaz2QAAAACSKC70AAAAAPBWIsQBAAAgISEOAAAACQlx\nAAAASEiIAwAAQEJCHAAAABIS4gAAAJCQEAcAAICEhDgAAAAkJMQBAAAgISEOAAAACQlxAAAASEiI\nAwAAQEKl+TjIpk2bYvbs2bFs2bJ47bXXYuDAgXHaaafFNddcE0OGDNlp28WLF0ddXV2sXbs2ysvL\n48Mf/nBcffXV0a9fv336rNbW1rj33nvjhz/8YWzcuDGOPvrouPjii2PKlClRVFSUj9MBAACAHlOU\nzWazB3KATZs2xcSJE+NPf/pTnH766VFdXR1r1qyJpUuXxsCBA+ORRx6JysrKiIi499574/bbb4/q\n6up43/veF01NTbFs2bIYM2ZM1NfXR1lZ2R4/a/v27fHpT386li1bFu9///ujqqoqnn766WhsbIxp\n06bFDTfccCCnAgAAAD3ugEP8pptuikceeSQ+//nPx2WXXda1/Pvf/35cf/31cdZZZ8U999wTzc3N\n8aEPfSje9a53xbx586JXr14REXHnnXfG3Llz48Ybb4yPf/zje/ysJ598Mq699tqdorujoyMuv/zy\n+N///d/4/ve/H9XV1QdyOgAAANCjDvg74j/96U9j0KBBMXXq1J2Wjxs3LoYOHRrPPvtsZDKZWLhw\nYXR2dsb06dO7IjwiYsaMGdG/f/949NFH9/pZ8+fPj9LS0pgxY0bXsl69esXMmTMjm83GokWLDvR0\nAAAAoEcd0HfEt2/fHtOnT4/S0tIoLt616cvKyqKjoyM6OztjxYoVERFxyimn7LRN796946STTopn\nn302tmzZEgMGDOj2s9rb2+M3v/lNjBgxIgYOHLjTutGjR0ffvn27PgMAAAAOVgcU4iUlJbtcCd9h\n9erV8eqrr8bQoUOjrKws1q1bF4cffni3D2WrqKiIiIg1a9bE6NGjuz1ec3NzdHZ2xtChQ7ud48gj\nj4y1a9fu/8kAAABAAj3y+rJMJhOzZs2KTCYTkyZNioiIzZs37/Zq947lLS0tuz3m5s2bd9q2u2O0\ntrZGZ2fngYwOAAAAPSrvIZ7NZuOmm26K5557LkaNGtV1xbyzs3O3T0XfsbytrW23x90R2AdyDAAA\nACi0vLxHfIfOzs648cYb4/HHH48hQ4bE3LlzuwK5T58+0dHR0e1+7e3tERHRt2/f3R67d+/eERF7\nPEZRUdEejwEAALw1ZDKZQo+Qk3nz5hV6hJy9GV8fvWHDhkKPEBF5DPHW1ta45pprYtmyZVFZWRkP\nPvhgHHHEEV3ry8vLY8uWLd3uu2P57m47j4iuB7Tt7vb1LVu2xCGHHNLtQ+MAAADgYJGXav3b3/4W\nU6dOjWXLlsUJJ5wQDz/8cBx11FE7bVNZWRmvvfZabNu2bZf9m5ubo7i4OI455pjdfkZFRUX06tUr\n1q9fv8u67du3x4YNG2LYsGEHfjIAAADQgw44xNva2mL69Onx0ksvxSmnnBLz5s2Lww47bJftxo4d\nG5lMJp5//vld9v/1r38dxx13XPTv33+3n1NaWhonnnhirFq1aper4i+//HK0trbGmDFjDvR0AAAA\noEcdcIjffvvt8eKLL8aYMWPi/vvv321Mf+QjH4mSkpKYM2dO13fCIyLuueeeaGlpiYsuumivnzV+\n/Phob2+P2bNndy3r6OiIO++8MyIiJk6ceIBnAwAAAD3rgL4jvmnTppg/f35ERBx77LFx//33d7vd\npz71qXjnO98Z06ZNi/vvvz/Gjx8fZ511Vvz+97+PpUuXxsknn9z1mrMd6urqYsuWLTF16tQoLy+P\niIgJEybEY489FnV1ddHU1BQjR46MZ555JhoaGmLatGlRXV19IKcDAAAAPa4om81m93fnn/70p/GZ\nz3xmr9utWLEiysvLI5vNxsMPPxwPP/xwrFu3LgYPHhxnn312XHXVVbs8qO0DH/hANDc3x89+9rM4\n+uiju5a3tLTE7Nmz48c//nFs3rw5hg4dGpMnT47Jkyd7UBsAABARnpqegqem778DCnEAAICDkRDv\neUJ8/7mEDAAAAAkJcQAAAEhIiAMAAEBCQhwAAAASEuIAAACQkBAHAACAhIQ4AAAAJCTEAQAAICEh\nDgAAAAkJcQAAAEhIiAMAAEBCQhwAAAASEuIAAACQkBAHAACAhIQ4AAAAJCTEAQAAICEhDgAAAAkJ\ncQAAAEhIiAMAAEBCQhwAAAASEuIAAACQkBAHAACAhIqy2Wy20EMAAAAHr0wmU+gRctbc3FzoEXIy\ndOjQQo+Qs/Xr1xd6hJxVVFQUeoSIcEUcAAAAkhLiAAAAkJAQBwAAgISEOAAAACQkxAEAACAhIQ4A\nAAAJCXEAAABISIgDAABAQkIcAAAAEhLiAAAAkJAQBwAAgISEOAAAACQkxAEAACAhIQ4AAAAJCXEA\nAABISIgDAABAQkIcAAAAEhLiAAAAkJAQBwAAgISEOAAAACQkxAEAACAhIQ4AAAAJCXEAAABISIgD\nAABAQkIcAAAAEhLiAAAAkJAQBwAAgISEOAAAACQkxAEAACCh0kIPAAAAbzWZTKbQI+Rk5cqVhR4h\nZ7NmzSr0CDl5+umnCz1Czo466qhCj/Cm5Yo4AAAAJCTEAQAAICEhDgAAAAkJcQAAAEhIiAMAAEBC\nQhwAAAASyvvryzZu3BjnnXde1NbWxqWXXtq1/AMf+EA0Nzfvcd/bbrstJkyYsMdtrrvuunjiiSe6\nXXfFFVfEtddem/PMAAAAkEpeQ3zr1q1RW1sbLS0tu6z7xCc+EVu2bNll+bZt2+KBBx6I3r17x7ve\n9a69fkZDQ0McfvjhUVNTs8u6sWPH7t/gAAAAkEjeQry5uTlqa2tj5cqV3a7/56vj/+wrX/lKZDKZ\n+NKXvhTDhw/f42d0dHTEmjVr4swzz4za2toDHRkAAACSy8t3xOvq6uKCCy6IhoaGOPXUU/d5v+ef\nfz4eeuihOP3002PixIl73X716tXR0dER1dXVBzIuAAAAFExeQry+vj4qKirioYceinHjxu3zfl/7\n2teipKQkvvzlL+/T9o2NjRERQhwAAIA3rbzcmn7LLbfEe9/73igpKYm1a9fu0z5LliyJl19+OSZN\nmhTHHnvsPu2zI8TXrFkTNTU10djYGH369IkzzzwzZs6cGUccccT+ngIAAAAkkZcr4meccUaUlJTk\ntM+DDz4YxcXF8clPfnKf99kR4nPnzo2jjz46LrrooqisrIzHH388Jk6cGBs2bMhpBgAAAEgt768v\n2xerVq2KX/3qV/GhD30oKisr93m/Pn36RGVlZcyZM2enB7vdfffdcccdd8RXvvKVmDNnTg9MDAAA\nAPlRkBBfvHhxRERMmjQpp/2+853vdLt8+vTpsWjRonjqqadi69at0a9fvwOeEQAAAHpCXm5Nz9VT\nTz0Vhx56aJx22ml5OV5xcXGMGDEiOjs73Z4OAADAQS15iL/66quxbt26+MAHPhClpft+Qb61tTV+\n/etfR0NDQ7frt23bFhERvXv3zsucAAAA0BOS35r+0ksvRUTE2LFjc9rvL3/5S1x00UVRVVUVTz75\n5E7rWltbY9WqVTFo0KCoqKjI26wAAACQb8mviK9atSoiIkaOHJnTfkOGDImRI0dGU1NTPPHEE13L\ns9lsfOtb34rXX389Jk+eHEVFRXmdFwAAAPIp+RXxP/zhDxER8fa3v32P29XV1cWWLVti6tSpUV5e\nHhERt956a1xyySVx/fXXx5IlS6KioiKef/75eOWVV+Ld7353zJgxo8fnBwAAgAOR/Ir45s2bIyJi\nwIABe9yuvr4+5syZE2+88UbXslGjRsWiRYvinHPOiRUrVsT8+fOjpaUlrr766njggQeirKysR2cH\nAACAA1WUzWazhR4CAADeSjKZTKFHyMnKlSsLPULOZs2aVegRclJbW1voEXL2L//yL4UeIWcHy1eZ\nC/L6MgAAAHirEuIAAACQkBAHAACAhIQ4AAAAJCTEAQAAICEhDgAAAAkJcQAAAEjIe8QBAHhTe7O9\nkzsi4ne/+12hR8jJOeecU+gRcrZ06dJCj5CToUOHFnqEnBUXu667v/ybAwAAgISEOAAAACQkxAEA\nACAhIQ4AAAAJCXEAAABISIgDAABAQkIcAAAAEhLiAAAAkJAQBwAAgISEOAAAACQkxAEAACAhIQ4A\nAAAJCXEAAABISIgDAABAQkIcAAAAEhLiAAAAkJAQBwAAgISEOAAAACQkxAEAACAhIQ4AAAAJCXEA\nAABISIgDAABAQkIcAAAAEhLiAAAAkJAQBwAAgISEOAAAACQkxAEAACAhIQ4AAAAJCXEAAABIqLTQ\nAwAAcPDIZDKFHiFnjY2NhR4hZ1OnTi30CDn5wQ9+UOgRcnbMMccUeoScFBUVFXoEEnJFHAAAABIS\n4gAAAJCQEAcAAICEhDgAAAAkJMQBAAAgISEOAAAACQlxAAAASEiIAwAAQEJCHAAAABIS4gAAAJCQ\nEAcAAICEhDgAAAAkJMQBAAAgISEOAAAACQlxAAAASEiIAwAAQEJCHAAAABIS4gAAAJCQEAcAAICE\nhDgAAAAkJMQBAAAgobyH+MaNG2Ps2LFRV1e3y7pHH300qquru/1n0qRJ+3T81tbWuOOOO+Lss8+O\n0aNHx3nnnRfz58+PbDab5zMBAACA/CvN58G2bt0atbW10dLS0u36xsbGiIi44ooronfv3jutO/LI\nI/d6/O3bt8c111wTy5Yti/e///1xzjnnxNNPPx233nprrF+/Pm644YYDPwkAAADoQXkL8ebm5qit\nrY2VK1fudpvGxsY49NBD49prr92vz/jRj34Uy5Yti2nTpnVF9zXXXBOXX355PPjggzF+/Piorq7e\nr2MDAABACnm5Nb2uri4uuOCCaGhoiFNPPXW32zU1NUVVVdV+f878+fOjtLQ0ZsyY0bWsV69eMXPm\nzMhms7Fo0aL9PjYAAACkkJcQr6+vj4qKinjooYdi3Lhx3W6zYcOG2Lx5835fsW5vb4/f/OY3MWLE\niBg4cOBO60aPHh19+/aNFStW7NexAQAAIJW83Jp+yy23xHvf+94oKSmJtWvXdrvNju+Hd3R0xJVX\nXhkvvvhibNu2LU4++eS45pprYvTo0Xv8jObm5ujs7IyhQ4fusq6kpCSOPPLI3X42AAAAHCzyckX8\njDPOiJKSkj1usyPEFyxYEG1tbTFhwoQ4/fTT47nnnospU6bEM888s8f9N2/eHBERAwYM6Hb9gAED\norW1NTo7O/fjDAAAACCNvD41fU8ymUxUVFTEzJkz46Mf/WjX8uXLl8ell14aX/jCF+JnP/vZLk9T\n32FHYJeVlXW7fsfytra2KC1NdloAAACQk2TFOmPGjJ0esrbDKaecEhdccEEsXrw4li9fHmeccUa3\n++8I9I6Ojm7Xt7e3R1FRUfTt2zd/QwMAHKBMJlPoEXLS1NRU6BFyNmnSpEKPkLM320OGhw8fXugR\nclZUVFToEWC38nJr+oE64YQTIiJi/fr1u91mxwPadveO8i1btsQhhxwSxcUHxSkBAABAt5JV68qV\nK3f7VPO2traIiN3elh4RUVFREb169eo21rdv3x4bNmyIYcOG5WdYAAAA6CHJQvwzn/lMfOITn4jX\nX399l3UvvPBCRESMGjVqt/uXlpbGiSeeGKtWrdrlqvjLL78cra2tMWbMmPwODQAAAHmWLMTPPffc\nyGQy8e1vfzuy2WzX8h//+MexdOnSePe73x1VVVV7PMb48eOjvb09Zs+e3bWso6Mj7rzzzoiImDhx\nYs8MDwAAAHmS7GFtV155ZTz99NOxcOHCaGxsjLFjx8aaNWti6dKlMXjw4Ljtttt22r6uri62bNkS\nU6dOjfLy8oiImDBhQjz22GNRV1cXTU1NMXLkyHjmmWeioaEhpk2bFtXV1alOBwAAAPZLsivi5eXl\nsWDBgpg6dWps2rQp5s2bF6+88kpceOGF8fjjj8eQIUN22r6+vj7mzJkTb7zxRteykpKS+O53vxuX\nXnpprF69Ourr62P79u1x0003xXXXXZfqVAAAAGC/FWX/+T5xAADyyuvLet6b8euJXl/W87xNiYOZ\n/zoBAAAgISEOAAAACQlxAAAASEiIAwAAQEJCHAAAABIS4gAAAJCQEAcAAICEhDgAAAAkJMQBAAAg\nISEOAAAACQlxAAAASEiIAwAAQEJCHAAAABIS4gAAAJCQEAcAAICEhDgAAAAkJMQBAAAgISEOAAAA\nCQlxAAAASEiIAwAAQEJCHAAAABIS4gAAAJBQaaEHAADYV5lMptAj5KypqanQI+Skpqam0CPkbOHC\nhYUeIWfDhw8v9Ag5KS52/Q7yyf9RAAAAkJAQBwAAgISEOAAAACQkxAEAACAhIQ4AAAAJCXEAAABI\nSIgDAABAQkIcAAAAEhLiAAAAkJAQBwAAgISEOAAAACQkxAEAACAhIQ4AAAAJCXEAAABISIgDAABA\nQkIcAAAAEhLiAAAAkJAQBwAAgISEOAAAACQkxAEAACAhIQ4AAAAJCXEAAABISIgDAABAQkIcAAAA\nEhLiAAAAkJAQBwAAgISEOAAAACQkxAEAACCh0kIPAAAURiaTKfQIOWtqair0CDmrqakp9Ag5WbBg\nQaFHyFlVVVWhR8hZcbHrYfBW5k8AAAAASEiIAwAAQEJCHAAAABIS4gAAAJCQEAcAAICEhDgAAAAk\nJMQBAAAgobyH+MaNG2Ps2LFRV1e3y7qWlpb4+te/HmeffXaMGjUq3vOe98SVV14Zv/3tb/f5+Ndd\nd11UV1d3+883v/nNPJ4JAAAA5F9pPg+2devWqK2tjZaWll3Wtba2xsUXXxwNDQ0xZsyY+OAHPxgb\nNmyIJUuWxLPPPhsPPvhgjB07dq+f0dDQEIcffnjU1NTssm5f9gcAAIBCyluINzc3R21tbaxcubLb\n9fPmzYuGhoa45JJL4stf/nLX8uXLl8ell14aN998czz55JN7/IyOjo5Ys2ZNnHnmmVFbW5uv0QEA\nACCZvNyaXldXFxdccEE0NDTEqaee2u02S5YsiaKiopg5c+ZOy0855ZQ45ZRToqmpKTZu3LjHz1m9\nenV0dHREdXV1PsYGAACA5PJyRby+vj4qKirilltuibVr18Yvf/nLXbapqamJ1157Lfr377/LurKy\nsoj4x63te9LY2BgRIcQBAAB408pLiN9yyy3x3ve+N0pKSmLt2rXdbnPhhRd2u/z111+P559/Pg45\n5JA4+uij9/g5O0J8zZo1UVNTE42NjdGnT58488wzY+bMmXHEEUcc0HkAAABAT8vLrelnnHFGlJSU\n7Ne+3/jGN2Lr1q0xbty4rivju7MjxOfOnRtHH310XHTRRVFZWRmPP/54TJw4MTZs2LBfMwAAAEAq\neX1qeq7mzp0bjz/+eFRUVMRnP/vZvW7fp0+fqKysjDlz5sTw4cO7lt99991xxx13xFe+8pWYM2dO\nT44MAAAAB6RgIX7nnXfG3Llz49BDD4177703Bg4cuNd9vvOd73S7fPr06bFo0aJ46qmnYuvWrdGv\nX798jwsAAAB5kZdb03Oxffv2+NKXvhRz586Nww47LL73ve/tdHV7fxQXF8eIESOis7PT7ekAAAAc\n1JJeEW9vb4+rr746nnrqqaioqIgHHnggKisr92nf1tbWroezjRgxYpf127Zti4iI3r1753NkAAAA\nyKtkIZ7NZuNzn/tcPPXUUzF8+PD4j//4j5yecv6Xv/wlLrrooqiqqoonn3xyp3Wtra2xatWqGDRo\nUFRUVOR7dAAAAMibZLemz5s3L5YsWRLHHHNM1NfX5/yqsSFDhsTIkSOjqakpnnjiia7l2Ww2vvWt\nb8Xrr78ekydPjqKionyPDgAAAHmT5Ip4e3t7zJ07NyIiqqurY/78+d1uV1NTE4MHD46IiLq6utiy\nZUtMnTo1ysvLIyLi1ltvjUsuuSSuv/76WLJkSVRUVMTzzz8fr7zySrz73e+OGTNmpDgdAAAA2G9J\nQnz16tXx17/+NSIilixZEkuWLOl2uw9+8INdIV5fXx/Nzc3xsY99rCvER40aFYsWLYq77rorfvnL\nX8bSpUujoqIirr766rjiiiv2+h5yAAAAKLSibDabLfQQAEB6mUym0CPkrKmpqdAj5KympqbQI+Rk\nwYIFhR4hZ1VVVYUeIWfFxclfXgQcRPwJAAAAAAkJcQAAAEhIiAMAAEBCQhwAAAASEuIAAACQkBAH\nAACAhJK8RxwA3grebK8D8yqwNN5srwPzKjCAnudPLQAAAEhIiAMAAEBCQhwAAAASEuIAAACQkBAH\nAACAhIQ4AAAAJCTEAQAAICEhDgAAAAkJcQAAAEhIiAMAAEBCQhwAAAASEuIAAACQkBAHAACAhIQ4\nAAAAJCTEAQAAICEhDgAAAAkJcQAAAEhIiAMAAEBCQhwAAAASEuIAAACQkBAHAACAhIQ4AAAAJCTE\nAQAAICEhDgAAAAkJcQAAAEhIiAMAAEBCQhwAAAASEuIAAACQkBAHAACAhEoLPQAAdCeTyRR6hJw1\nNTUVeoSc1NTUFHqEnC1YsKDQI+Ssqqqq0CPkpLjYdRqAnuZPWgAAAEhIiAMAAEBCQhwAAAASEuIA\nAACQkBAHAACAhIQ4AAAAJCTEAQAAICEhDgAAAAkJcQAAAEhIiAMAAEBCQhwAAAASEuIAAACQkBAH\nAACAhIQ4AAAAJCTEAQAAICEhDgAAAAkJcQAAAEhIiAMAAEBCQhwAAAASEuIAAACQkBAHAACAhPIe\n4hs3boyxY8dGXV1dt+sXL14c48ePj5NOOine9773xW233RZbt27d5+O3trbGHXfcEWeffXaMHj06\nzjvvvJg/f35ks9k8nQEAAAD0nLyG+NatW6O2tjZaWlq6XX/vvffGDTfcEJlMJj7+8Y/HiBEjoq6u\nLj75yU9Ge3v7Xo+/ffv2uOaaa+Luu++OYcOGxSc+8YkoLS2NW2+9Nb7+9a/n81QAAACgR+QtxJub\nm+OSSy6Jl156abfr77rrrhgzZkw89thjce2118Z9990XV155Zbz44ouxcOHCvX7Gj370o1i2bFlM\nmzYt7rvvvrj22mvjsccei1NPPTUefPDBaGxszNfpAAAAQI/IS4jX1dXFBRdcEA0NDXHqqad2u83C\nhQujs7Mzpk+fHr169epaPmPGjOjfv388+uije/2c+fPnR2lpacyYMaNrWa9evWLmzJmRzWZj0aJF\nB34yAAAA0IPyEuL19fVRUVERDz30UIwbN67bbVasWBEREaeccspOy3v37h0nnXRSNDQ0xJYtW3b7\nGe3t7fGb3/wmRowYEQMHDtxp3ejRo6Nv375dnwEAAAAHq7yE+C233BKLFy+Ok08+ebfbrFu3Lg4/\n/PDo16/fLusqKioiImLNmjW73b+5uTk6Oztj6NChu6wrKSmJI488MtauXZv78AAAAJBQXkL8jDPO\niJKSkj1us3nz5hgwYEC363Ys391D3nbs/8/bdneM1tbW6Ozs3JeRAQAAoCCSvUe8s7MzysrKul23\nY3lbW9se9//nbffnGAAAAFBopak+qE+fPtHR0dHtuh2vLuvbt+9u9+/du3dExB6PUVRUtMdjALxV\nZTKZQo+Qs6ampkKPkLOamppCj5CTBQsWFHqEnFVVVRV6hJwVFye77gHAm0SyvxnKy8t3+zC2Hct3\nd9t5RHQ9oG13t69v2bIlDjnkEH/ZAQAAcFBLVq2VlZXx2muvxbZt23ZZ19zcHMXFxXHMMcfsdv+K\nioro1atXrF+/fpd127dvjw0bNsSwYcPyOjMAAADkW7IQHzt2bGQymXj++ed3Wt7W1ha//vWv47jj\njov+/fvvdv/S0tI48cQTY9WqVbtcFX/55ZejtbU1xowZ0yOzAwAAQL4kC/GPfOQjUVJSEnPmzOn6\nTnhExD333BMtLS1x0UUX7fUY48ePj/b29pg9e3bXso6OjrjzzjsjImLixIn5HxwAAADyKNnD2t75\nznfGtGnT4v7774/x48fHWWedFb///e9j6dKlcfLJJ8ekSZN22r6uri62bNkSU6dOjfLy8oiImDBh\nQjz22GNRV1cXTU1NMXLkyHjmmWeioaEhpk2bFtXV1alOBwAAAPZLshCPiPjc5z4X73jHO+Lhhx+O\n+vr6GDx4cFx66aVx1VVX7fJasvr6+mhubo6PfexjXSFeUlIS3/3ud2P27Nnx4x//OF544YUYOnRo\n3HTTTTF58uSUpwIAAAD7pSibzWYLPQQAPcvry9Lw+rKe5/VlAPy/wN8MAAAAkJAQBwAAgISEOAAA\nACQkxAEAACAhIQ4AAAAJCXEAAABISIgDAABAQkIcAAAAEhLiAAAAkJAQBwAAgISEOAAAACQkxAEA\nACAhIQ4AAAAJCXEAAABISIgDAABAQkIcAAAAEhLiAAAAkJAQBwAAgISEOAAAACQkxAEAACAhIQ4A\nAAAJCXEAAABIqLTQAwC8GWUymUKPkJOmpqZCj5CzmpqaQo+QswULFhR6hJxUVVUVeoScFRe7hgDA\nm5+/zQAAACAhIQ4AAAAJCXEAAABISIgDAABAQkIcAAAAEhLiAAAAkJAQBwAAgISEOAAAACQkxAEA\nACAhIQ4AAAAJCXEAAABISIgDAABAQkIcAAAAEhLiAAAAkJAQBwAAgISEOAAAACQkxAEAACAhIQ4A\nAAAJCXEAAABISIgDAABAQkIcAAAAEhLiAAAAkJAQBwAAgISEOAAAACQkxAEAACAhIQ4AAAAJCXEA\nAABISIgDAABAQkIcAAAAEiot9AAAmUym0CPkrKmpqdAj5KSmpqbQI+RswYIFhR4hZ1VVVYUeISfF\nxX4fDwDphSQxAAAgAElEQVSF4G9gAAAASEiIAwAAQEJCHAAAABIS4gAAAJCQEAcAAICEhDgAAAAk\nlPT1ZdXV1Xvdpr6+Pt7znvfscZspU6bECy+80O26m2++OSZPnrxf8wEAAEBPSxriV111VbfLX3vt\ntfjP//zPOOyww+LYY4/d63EaGxtj2LBhcf755++ybtSoUQc8JwAAAPSUpCFeW1vb7fJPf/rTUVRU\nFN/4xjdi8ODBezzG+vXro6WlJS688MLdHg8AAAAOVgX/jvgTTzwRP//5z2PixIlx+umn73X7xsbG\niNi329wBAADgYFPQEG9ra4tvf/vbMWDAgPjXf/3XfdpHiAMAAPBmVtAQf/jhh+OPf/xjXH755fG2\nt71tn/ZpbGyMoqKieOGFF+JjH/tYnHTSSfG+970vvvrVr8aWLVt6eGIAAAA4MAUL8e3bt0d9fX30\n69cvpkyZss/7NTY2RjabjbvuuitOOOGEmDhxYgwaNCjq6+tjypQp0dLS0oNTAwAAwIFJ+rC2f/bz\nn/88/vjHP8Zll10W5eXl+7RPJpOJ8vLyOP744+Pee++NI444omv5zTffHI888kjMnj07vvCFL/Tk\n6AAAALDfCnZFfPHixRERMWnSpH3ep7i4OBYuXBiLFy/uivAdy2+44Ybo27dv/PCHP8z7rAAAAJAv\nBQnxtra2+MUvfhFVVVX79N7wfdGvX7+orKyMTZs2RVtbW16OCQAAAPlWkBBfvnx5/P3vf49zzjkn\np/3eeOON+NWvfhVr1qzpdv22bduiuLg4SksLdsc9AAAA7FFBQvyll16KiIixY8fmtN/KlStj8uTJ\n8bWvfW2XdX/+859j/fr1cfzxx0dJSUle5gQAAIB8K0iIr1q1KiIiRo4cmdN+Y8eOjcGDB8fTTz8d\nK1as6Fre3t4es2bNio6Ojrj44ovzOisAAADkU0Hu4f7DH/4Qffr02evT0mfPnh0REbW1tRERUVZW\nFrNmzYqrrroqLrvssjj33HPj0EMPjV/84hexevXqOP/882PChAk9Pj8AAADsr4KE+F//+tcYMGDA\nXrebM2dORPxfiEdEnHXWWTF//vyYO3duLF26NNra2mLYsGFx4403xpQpU6KoqKjH5gYAAIADVZAQ\nf/bZZ/dpu8bGxm6Xn3TSSXHfffflcyQAAABIomDvEQcAAIC3IiEOAAAACQlxAAAASEiIAwAAQEJC\nHAAAABIS4gAAAJCQEAcAAICECvIecaDnZDKZQo+Qs8bGxkKPkLPJkycXeoScLFiwoNAj5KyqqqrQ\nI+SsuNjvtwGAvfMTAwAAACQkxAEAACAhIQ4AAAAJCXEAAABISIgDAABAQkIcAAAAEhLiAAAAkJAQ\nBwAAgISEOAAAACQkxAEAACAhIQ4AAAAJCXEAAABISIgDAABAQkIcAAAAEhLiAAAAkJAQBwAAgISE\nOAAAACQkxAEAACAhIQ4AAAAJCXEAAABISIgDAABAQkIcAAAAEhLiAAAAkJAQBwAAgISEOAAAACQk\nxAEAACAhIQ4AAAAJCXEAAABISIgDAABAQkXZbDZb6CHgYJbJZAo9Qk5+97vfFXqEnF144YWFHiFn\njz76aKFHyElVVVWhR8hZcbHfFQMA/2/yUw4AAAAkJMQBAAAgISEOAAAACQlxAAAASEiIAwAAQEJC\nHAAAABIS4gAAAJCQEAcAAICEhDgAAAAkJMQBAAAgISEOAAAACQlxAAAASEiIAwAAQEJCHAAAABIS\n4gAAAJCQEAcAAICEhDgAAAAkJMQBAAAgISEOAAAACQlxAAAASEiIAwAAQEKlqT/wjjvuiLvvvrvb\ndeedd158+9vf3uP+mzdvjrvuuiuWLl0ar732Wrzzne+Myy+/PM4777yeGBcAAADyKnmINzQ0RFlZ\nWXzqU5/aZd3w4cP3uO/f//73mDZtWvz2t7+Nc889N97xjnfEkiVL4rOf/Wy8/vrr8fGPf7ynxgYA\nAIC8SB7iTU1Ncdxxx0VtbW3O+9bX18fKlSvjpptuiosvvjgiIq688sqoqamJb37zm/HhD384Djvs\nsHyPDAAAAHmT9DviLS0t0dzcHNXV1fu1/8MPPxyHH3541NTUdC3r379/zJgxI1pbW+PJJ5/M16gA\nAADQI5KGeENDQ0TEfoX4unXrYuPGjTF27NgoKSnZad173vOeiIhYsWLFgQ8JAAAAPSjpremNjY0R\nEfH666/HZZddFq+88kpERJx22mkxc+bMOPbYY3e777p16yIiYujQobusGzx4cPTu3TvWrl2b/6EB\nAAAgj5JeEd8R4g888ED0798/Jk6cGKNHj46f/OQnMWnSpPjtb3+72303b94cERHl5eXdru/fv39s\n2bIl/0MDAABAHiW9Il5SUhIVFRVx2223dd1OHhHxxBNPxHXXXRdf/OIX47/+67+63bezszMiIsrK\nyrpdX1ZWFq2trfkfGgAAAPIoaYj/27/9W7fLP/rRj8bChQtjxYoV8eqrr3Z7i3rv3r0jIqK9vb3b\nY7S3t8chhxySv2HpEdlsttAj5GzVqlWFHiEn06ZNK/QIOVu4cGGhR8hZVVVVoUfISXFx0hugAADY\ng4PmJ7MTTjghIiLWr1/f7fqBAwdGxD+evN6dlpaW6N+/f88MBwAAAHmS7Ip4Z2dnrFq1KrLZbJx4\n4om7rN+2bVtE/N+V7/+/ysrKiOg+1P/85z9HW1tbDBs2LH8DAwAAQA9IFuKZTCamTJkShxxySDz3\n3HM7vYIsm83Giy++GKWlpXH88cd3u/9RRx0VRx11VLzwwguRyWR2us1y+fLlERExZsyYnj0JAAAA\nOEDJbk0vKyuLs846K/72t7/Ffffdt9O6Bx54IJqamuIjH/nIbp+KHvGP75Jv2LAhHnrooa5lLS0t\ncc8990SfPn1i3LhxPTY/AAAA5EPSh7XdcMMN8eKLL8Ydd9wRy5cvjxEjRsQrr7wSy5cvj+OOOy4+\n//nPd207e/bsiIiora3tWnbFFVfEf//3f8dXv/rVWLFiRQwZMiSWLFkSf/jDH+LGG2+MQYMGpTwd\nAAAAyFlRNvFjrDdu3Bh33nlnPP3007F58+Z4+9vfHuecc05ceeWVMWDAgK7tqqurI+L/3j2+w1/+\n8pe4/fbb46mnnorW1tY49thj45Of/GScf/75KU+D/fRmfGr6ypUrCz1CTt6MT03/3ve+V+gRcrbj\nz6g3C09NBwA4eCQPcd7a3oz/uQnxnifEe54QBwA4ePjJDAAAABIS4gAAAJCQEAcAAICEhDgAAAAk\nJMQBAAAgISEOAAAACQlxAAAASEiIAwAAQEJCHAAAABIS4gAAAJCQEAcAAICEhDgAAAAkJMQBAAAg\nISEOAAAACQlxAAAASEiIAwAAQEJCHAAAABIS4gAAAJCQEAcAAICEhDgAAAAkJMQBAAAgISEOAAAA\nCRVls9lsoYdg/2UymUKPkJN169YVeoScnXnmmYUeISc/+clPCj1CzoYPH17oEXJWXOz3mAAA7B8/\nSQIAAEBCQhwAAAASEuIAAACQkBAHAACAhIQ4AAAAJCTEAQAAICEhDgAAAAkJcQAAAEhIiAMAAEBC\nQhwAAAASEuIAAACQkBAHAACAhIQ4AAAAJCTEAQAAICEhDgAAAAkJcQAAAEhIiAMAAEBCQhwAAAAS\nEuIAAACQkBAHAACAhIQ4AAAAJCTEAQAAICEhDgAAAAkJcQAAAEhIiAMAAEBCQhwAAAASEuIAAACQ\nkBAHAACAhEoLPcDBJJvNFnqEnP3P//xPoUfIyezZsws9Qs6efPLJQo+Qk+HDhxd6hJwVF/udIAAA\nbx1++gUAAICEhDgAAAAkJMQBAAAgISEOAAAACQlxAAAASEiIAwAAQEJCHAAAABJK/h7xTZs2xezZ\ns2PZsmXx2muvxcCBA+O0006La665JoYMGbLX/adMmRIvvPBCt+tuvvnmmDx5cr5HBgAAgLxJGuKb\nNm2KiRMnxp/+9Kc4/fTT47zzzos1a9bED37wg3jmmWfikUceicrKyj0eo7GxMYYNGxbnn3/+LutG\njRrVQ5MDAABAfiQN8dmzZ8ef/vSn+PznPx+XXXZZ1/Lvf//7cf3118e///u/xz333LPb/devXx8t\nLS1x4YUXRm1tbYqRAQAAIK+Sfkf8pz/9aQwaNCimTp260/Jx48bF0KFD49lnn41MJrPb/RsbGyMi\norq6ukfnBAAAgJ6S7Ir49u3bY/r06VFaWhrFxbv2f1lZWXR0dERnZ2eUlZV1ewwhDgAAwJtdshAv\nKSnZ5Ur4DqtXr45XX301hg4dutsIj/hHiBcVFcULL7wQX/7yl2PNmjVRXl4e55xzTlx99dUxYMCA\nnhofAAAA8qLgry/LZDIxa9asyGQyMWnSpD1u29jYGNlsNu6666444YQTYuLEiTFo0KCor6+PKVOm\nREtLS6KpAQAAYP8kf33ZP8tms3HTTTfFc889F6NGjdrtFfOIfwR7eXl5HH/88XHvvffGEUcc0bX8\n5ptvjkceeSRmz54dX/jCF1KNDwAAADkr2BXxzs7O+OIXvxiPPvpoDBkyJObOnbvH29KLi4tj4cKF\nsXjx4q4I37H8hhtuiL59+8YPf/jDFKMDAADAfitIiLe2tsaVV14Zjz/+eFRWVkZ9ff1OcZ2rfv36\nRWVlZWzatCna2tryOCkAAADkV/Jb0//2t7/FFVdcES+99FKccMIJ8d3vfjcOO+ywve73xhtvxO9/\n//t429veFsOGDdtl/bZt26K4uDhKSwt6tz0AAADsUdIr4m1tbTF9+vR46aWX4pRTTol58+btU4RH\nRKxcuTImT54cX/va13ZZ9+c//znWr18fxx9/fJSUlOR7bAAAAMibpCF+++23x4svvhhjxoyJ+++/\nP/r377/P+44dOzYGDx4cTz/9dKxYsaJreXt7e8yaNSs6Ojri4osv7omxAQAAIG+S3ce9adOmmD9/\nfkREHHvssXH//fd3u92nPvWp6N27d8yePTsiImprayMioqysLGbNmhVXXXVVXHbZZXHuuefGoYce\nGr/4xS9i9erVcf7558eECRPSnAwAAADsp2Qh/tJLL0VHR0dERDz22GO73W7q1KnRu3fvmDNnTkT8\nX4hHRJx11lkxf/78mDt3bixdujTa2tpi2LBhceONN8aUKVOiqKioZ08CAAAADlCyEP/gBz8YjY2N\n+7z97rY96aST4r777svXWAAAAJBUwd4jDgAAAG9FQhwAAAASEuIAAACQkBAHAACAhIQ4AAAAJCTE\nAQAAIKGibDabLfQQB4vm5uZCj5Czo48+utAj5GTdunWFHiFnFRUVhR4hJ8XFfr8GAAAHMz+xAwAA\nQEJCHAAAABIS4gAAAJCQEAcAAICEhDgAAAAkJMQBAAAgISEOAAAACQlxAAAASEiIAwAAQEJCHAAA\nABIS4gAAAJCQEAcAAICEhDgAAAAkJMQBAAAgISEOAAAACQlxAAAASEiIAwAAQEJCHAAAABIS4gAA\nAJCQEAcAAICEhDgAAAAkJMQBAAAgISEOAAAACf1/7d17UFTnGcfx34ICKiKixJRVgyZZvFVqGLXG\ndqq2DQYNWlu8oPValVFRMTXGNKaptrWZtiho8DajFBGrpuqYxjQk8T5JqyBqvABTxUFMQJQhiiIX\nOf3DgYawu1yEI+r38+f7nmf3Pfvy8O5z9lwoxAEAAAAAMBGFOAAAAAAAJqIQBwAAAADARBTiAAAA\nAACYiEIcAAAAAAATUYgDAAAAAGCiFg97AM1JUFDQwx5CvcXHxz/sIdSL1Wp92EOoNxcXjlcBAAAA\naDxUGAAAAAAAmIhCHAAAAAAAE1GIAwAAAABgIgpxAAAAAABMRCEOAAAAAICJKMQBAAAAADARhTgA\nAAAAACaiEAcAAAAAwEQU4gAAAAAAmIhCHAAAAAAAE1GIAwAAAABgIgpxAAAAAABMRCEOAAAAAICJ\nKMQBAAAAADARhTgAAAAAACaiEAcAAAAAwEQU4gAAAAAAmIhCHAAAAAAAE1GIAwAAAABgIgpxAAAA\nAABMRCEOAAAAAICJTC/Ey8vLFR8fr5CQEPXt21c//vGP9e6776qsrKxO8YWFhVq+fLmGDRumwMBA\njRkzRvv372/iUQMAAAAA0DhML8SXL1+ulStXytvbW5MnT1anTp0UGxurV199tdbYO3fuaPr06dq+\nfbsCAwM1ceJE3bx5U1FRUUpMTDRh9AAAAAAAPJgWZr7ZyZMntWPHDgUHBysmJkYWi0WGYej111/X\n3r17dfDgQQ0dOtRhfEJCgs6dO6e33npLEydOlCTNmTNH48eP11/+8he9/PLL6tChg1m7AwAAAABA\nvZn6i/i2bdskSfPmzZPFYpEkWSwWLVq0SBaLRbt27XIan5SUpI4dO2r8+PFVbZ6enoqIiFBxcbHe\nf//9phs8AAAAAACNwNRCPCUlRe3bt5fNZqvW3qlTJ/n7++vEiRMOY7Ozs5WXl6egoCC5urpW6xs4\ncKAkOY0HAAAAAKA5MK0QLy0tVW5urrp27Wq332q16ubNmyooKLDbn52dLUl24319feXu7q7Lly83\n2ngBAAAAAGgKphXihYWFkqS2bdva7a9sv3XrltN4Ly8vu/2enp4OYwEAAAAAaC5MK8TLy8slSW5u\nbnb7K9tLSkoaHO8oFgAAAACA5sK0u6Z7eHhIksPnhZeWlkqSWrVqZbff3d292nb24lu3bv1AY8zN\nzX2geAAAAAAAamPaL+Kenp5ycXFRUVGR3f7K08odnbrerl07SXIYX1RUJE9Pz0YYKQAAAAAATce0\nQtzNzU1+fn7Kycmx25+TkyMfHx95e3vb7ff396/a7tuuXbumkpISdevWrdHGCwAAAABAUzD18WVB\nQUHKz89XVlZWtfa8vDxdvnxZgYGBDmP9/Pzk5+en1NRUVVRUVOs7fvy4JKlfv36NP2gAAAAAABqR\nqYX46NGjJUmrVq2qKqYNw1B0dLQkady4cU7jQ0NDlZubq8TExKq2oqIirV+/Xh4eHho1alQTjRwA\nAAAAgMZhMQzDMPMNo6KitH//fvXt21cDBw5UWlqaUlJSFBwcrJiYGFksFknSmjVrJEmRkZFVsUVF\nRfr5z3+uy5cv66WXXlKXLl2UnJysK1euaNmyZZo0aZKZuwIAAAAAQL2ZXoiXlZVp48aN2rNnj/Ly\n8uTn56fQ0FDNnDmz2qPJAgICJEkZGRnV4q9fv67o6GgdPHhQxcXF6t69u2bMmKERI0aYuRsAAAAA\nADSI6YU4AAAAAABPMlOvEQcAAAAA4ElHIQ4AAAAAgIkoxAEAAAAAMBGFOAAAAAAAJqIQBwAAAADA\nRBTiAAAAAACYqMXDHkBTKy8vV2Jionbu3KmcnBz5+vpqzJgxmjVrllq2bFlrfGFhoWJjY3Xo0CHd\nuHFDzz77rH71q18pJCTEhNEjPz9fa9as0eHDh3Xjxg21a9dOgwYN0oIFC9SlS5da48PDw5Wammq3\n7+2339aECRMae8j4htWrV2vdunV2+0JCQrRq1Sqn8eSf+QICAmrdJiEhQQMHDnS6Dblnrry8PIWE\nhCgyMlJTp06t0b93717Fx8fr8uXL8vLy0ssvv6z58+erTZs2dXr94uJibdiwQR988IHy8vLUuXNn\nTZw4UeHh4bJYLI28N08eZ/NXVFSkuLg4ffzxx/rqq6/Upk0bBQUFKTIyUj179qzT6y9evFj79u2z\n2zdz5kz9+te/ftBdeKI5m79du3bpzTfftBsXGBionTt31vr65F/TcTR3w4YN09WrV53Grly5UmPG\njHG6DbnX+OpTGzT3te+xL8SXL1+uHTt2KCgoSMOGDdPJkycVGxurjIwMxcbGOo29c+eOpk+frgsX\nLmj48OH6zne+o+TkZEVFRamgoECTJk0yaS+eTPn5+QoLC9NXX32lwYMHKyQkRFlZWfrnP/+po0eP\naseOHfL393f6GhkZGerWrZtGjBhRo69Pnz5NNHJUSk9Pl5ubm2bNmlWj7/nnn3caS/49HPPmzbPb\nfuPGDW3fvl0dOnRQ9+7da30dcs88t2/fVmRkpIqKiuz2b9iwQdHR0QoICNCkSZOUmZmp+Ph4nT59\nWgkJCXJzc3P6+vfu3dOCBQt0+PBh/ehHP1JwcLCOHDmi5cuXKycnR0uWLGmK3XpiOJu/4uJiTZw4\nUenp6erXr59+8pOfKDc3V8nJyTp27Ji2bNmioKCgWt8jPT1dHTt21Pjx42v01SUejtWWfxkZGZLu\nF13u7u7V+p5++ulaX5/8azrO5m7y5Mm6detWjfa7d+9q8+bNcnd313e/+91a34Pca1z1qQ0eibXP\neIylpqYaNpvNiIyMNCoqKgzDMIyKigrjtddeM2w2m3HgwAGn8evWrTNsNpuRmJhY1Xbr1i1jxIgR\nRmBgoHH9+vUmHf+TbtmyZYbNZjM2b95crX3v3r2GzWYzZs+e7TT+ypUrhs1mM/74xz825TDhxNCh\nQ43Ro0c3KJb8a14iIiKMgIAA49ixY7VuS+6ZJycnx/jZz35m2Gw2w2azGVu2bKnR36tXL2PcuHFG\naWlpVfvq1asNm81mbN26tdb32Ldvn2Gz2Yw//elPVW2lpaXG5MmTjYCAACM9Pb3R9udJU9v8bdiw\nwbDZbMaKFSuqtf/nP/8xevbsaYwcObLW9ygtLTV69+5tzJ07tzGHDqP2+TMMw5g0aZIxYMCABr8H\n+dc06jJ39qxYscKw2WzGzp07a92W3Gt8da0NHpW177G+Rnzbtm2S7v/CU3n6gMVi0aJFi2SxWLRr\n1y6n8UlJSTWOYnl6eioiIkLFxcV6//33m27w0CeffCIfHx9NmTKlWvuoUaPUtWtXHTt2TBUVFQ7j\nK49C1+VUWzS+oqIiXb16tcGfP/nXfOzbt08HDhxQWFiYBg8eXOv25J454uPj9corryg9PV3f//73\n7W6zc+dOlZeXa/bs2dUux4qIiJCnp2et66B0fy1t0aKFIiIiqtpatmyphQsXyjAMvffeew++M0+g\nusxfcnKyLBaLFi5cWK19wIABGjBggDIzM5WXl+f0fS5evKiysjLysZHVZf4kKTMzUzabrcHvQ/41\nvrrO3belpKQoMTFRgwcPVlhYWK3bk3uNr661waOy9j3WhXhKSorat29f4x9gp06d5O/vrxMnTjiM\nzc7OVl5enoKCguTq6lqtr/LaSGfxeDD37t3T7NmzNW/ePLm41PwzdXNzU1lZmcrLyx2+BsXAw5We\nni6pYZ8/+dd8lJSUaNWqVWrbtq0WLVpUpxhyzxwJCQmyWq1KTEzUqFGj7G5TmScDBgyo1u7u7q7v\nfe97Sk9Pt3v6ZaXS0lJ98cUX6tGjh9q1a1etr2/fvmrVqhW52EB1mb/x48crKipKnp6eNfoqT6u8\nffu20/chH5tGXeYvNzdXhYWFDf7syb+mUZe5s+edd96Rq6urw2v+v43ca1z1qQ0elbXvsb1GvLS0\nVLm5uQoMDLTbb7ValZWVpYKCAvn4+NToz87OliR17dq1Rp+vr6/c3d11+fLlRh0z/s/V1bXG0a5K\nFy9e1KVLl9S1a1en13dkZGTIYrEoNTVVb775prKysuTl5aXg4GDNnz9fbdu2barhQ/9fgAoKCjRt\n2jSdPXtWkjRo0CAtXLjQ6XXG5F/zkZSUpC+//FJRUVFq3759nWLIPXP87ne/04svvihXV1eH+ZCd\nna2OHTvavTGN1WqVJGVlZalv3752469evary8nK7uejq6qqnn36aXGyguszfL37xC7vtBQUFSklJ\nUevWrdW5c2en71P5vzgrK0vjx49XRkaGPDw8NGTIEC1cuFCdOnV6oP14UtVl/io/+7KyMs2ZM0dp\naWm6e/euXnjhBS1YsMBh3lUi/5pGXebu25KTk3XmzBmNHTu2TvdJkci9xlaf2uBRWfse21/ECwsL\nJcnhF77KdkdHQyrjvby87PZ7eno6PZKCplFRUaEVK1aooqJCY8eOdbptRkaGDMNQbGysevXqpbCw\nMPn4+CghIUHh4eEOb6yCxlG5AG3evFmenp4KCwtT37599dFHH2ns2LG6cOGCw1jyr3m4d++eEhIS\n1KZNG4WHh9c5jtwzxw9/+MMaZ4x8W2FhYa3roLP5qMtaWlxc7PTsJNhXl/lz5M9//rNu376tUaNG\n1XrDocr/xXFxcercubPGjRsnf39/7d69W2FhYcrNzW3QGJ50dZm/ys/+73//u0pKSjRmzBgNHjxY\nn3/+ucLDw3X06FGn8eRf02hI7m3ZskUuLi6aMWNGnWPIPXPYqw0elbXvsf1FvPKDcbRAVbaXlJQ0\nOL64uPhBh4l6MAxDb731lj7//HP16dPH4VEx6X5Senl5qWfPntqwYUPVUceKigq9/fbb2rFjh9as\nWaOlS5eaNfwnjqurq6xWq1auXFntUVf79u3T4sWL9cYbb2jPnj12Y8m/5uHAgQP68ssvNW3aNIcH\nRb6N3GteysvLG7wOVsZ/c1tnr9GixWP7laJZiYuL0+7du2W1WhUVFVXr9h4eHvL399fatWurPa1i\n3bp1Wr16tX7/+99r7dq1TTnkJ1ZFRYWsVqsWLlyo0NDQqvbjx49r6tSpWrp0qT799NMad1OvRP41\nD+fPn9fJkyf10ksv1fq0nm8i95qeo9rgUVn7HttfxD08PCTdPx3IntLSUklSq1at7PZX/lOs3M5e\nfOvWrR90mKij8vJyvfHGG9q1a5e6dOmiuLg4p78CuLi4aOfOndq7d2+1U39cXFy0ZMkStWrVSh98\n8IEZQ39i/fa3v9WBAwdqPG86NDRU/fv31/nz53Xp0iW7seRf87B3715JqvXsk28i95oXDw+PBq+D\n0v9z0dlrWCwWp6+BxhMTE6OYmBh5e3trw4YNNa5dtOfdd9/VRx99VOORkbNnz1bnzp118ODBWq8z\nR8NERETowIED1Ypw6f51q6+88ory8/N1/Phxh/HkX/PQkLVQIveamrPa4FFZ+x7bQtzT01MuLi4O\nT3oH2YUAAAeySURBVDuoPK3V0SkHlYubo/iioiK7N09B4ysuLtacOXO0e/du+fv7KyEh4YGuq2nT\npo38/f2Vn5/v9GgYmk6vXr0kSTk5OXb7yb+Hr6SkRJ999plsNludr4erDblnPi8vL4eXcdS2Dkq1\n5+KtW7fUunVruzfOQeO5d++efvOb3yguLk4dOnTQ3/72txpf7uvLxcVFPXr0UHl5OafIPgS1rYMS\n+ddcHDx4UN7e3ho0aFCjvB659+Bqqw0elbXvsc1cNzc3+fn5OfwHl5OTIx8fH3l7e9vtrzz1xF78\ntWvXVFJSom7dujXaeGHf119/rSlTpujw4cPq1auXkpKS5OfnV2vczZs3dfLkSWVlZdntv3v3rlxc\nXDiVq4mUl5frzJkzOn36tN3+u3fvSpLD0/HIv4fv+PHjunPnjoKDg+sVR+41L/7+/rpx40ZVzn3T\n1atX5eLiomeeecZhvNVqVcuWLe3m4r1795Sbm0suNrHS0lLNnTtX7733nqxWq5KSktSjR486xRYX\nF+vUqVNVT7H4ttr+F+PBnDt3zuGdlSsPRjr77Mm/h+/SpUvKzs7WsGHD6rVukXtNpy61waOy9j22\nhbgkBQUFKT8/v8YXwry8PF2+fNnhHdUlyc/PT35+fkpNTa3xrOrK04j69evX+INGlZKSEs2ePVun\nT5/WgAEDtHXrVnXo0KFOsefOndOECRP0zjvv1Oi7du2acnJy1LNnzwbfKAfOVVRUKDw8XDNnztS9\ne/eq9RmGobS0NLVo0UI9e/a0G0/+PXyVB1GCgoLqFUfuNS9BQUGqqKhQSkpKtfaSkhKdOnVKzz33\nnNOzS1q0aKHAwECdP3++xi8DZ86cUXFxMbnYhAzD0KuvvqqDBw/q+eef1/bt2+t1jer169c1btw4\nLV68uEZfcXGxzp8/Lx8fn6q7CKNxzZ07V5MnT1ZBQUGNvtTUVElSnz59HMaTfw9fQ9dCcq9p1LU2\neFTWvse6EB89erQkadWqVVVf5g3DUHR0tCRp3LhxTuNDQ0OVm5urxMTEqraioiKtX79eHh4e9Xr2\nIOovOjpaaWlp6tevnzZt2lSvU5GDgoLk6+urI0eOVDsaXVpaqhUrVqisrEwTJ05simFD989IGTp0\nqL7++mtt3LixWt/mzZuVmZmpkSNHOr0BGPn3cJ0/f16S1Lt373rFkXvNy8iRI+Xq6qq1a9dWu+fC\n+vXrVVRUVOs6KN1fS0tLS7VmzZqqtrKyMsXExEiSwsLCGn/gkCRt3bpVycnJeuaZZxp0WVaXLl3U\nu3dvZWZmat++fVXthmHor3/9qwoKCjRhwgRZLJbGHjokDR8+XBUVFVq1apUMw6hq//DDD3Xo0CH1\n799fNpvN6WuQfw9XQ9dCcq9p1LU2eFTWvsf63MAXX3xRISEh2r9/v8aNG6eBAwcqLS1NKSkpCg4O\n1pAhQ6q2rfyQIyMjq9pmzpypf/3rX/rDH/6gEydOqEuXLkpOTtaVK1e0bNkyu88fR+PIz8/Xtm3b\nJEndu3fXpk2b7G43a9Ysubu715g/Nzc3rVixQvPmzdO0adM0fPhweXt767PPPtPFixc1YsQIjRkz\nxpydeUItWbJEaWlpWr16tY4fP64ePXro7NmzOn78uJ577jm9/vrrVduSf83PlStX5OHhUevd0sm9\n5u3ZZ5/V9OnTtWnTJo0ePVpDhw7Vf//7Xx06dEgvvPBCjZsPxcfH69atW5oyZUrV3I8ZM0b/+Mc/\nFB8fr8zMTPXu3VtHjx5Venq6pk+froCAgIexa4+90tJSxcXFSZICAgKq1sRvGz9+vHx9fSXZn7/l\ny5frl7/8pV577TUlJyfLarUqJSVFZ8+eVf/+/RUREWHODj2B5syZoyNHjmjnzp3KyMhQUFCQsrKy\ndOjQIfn6+mrlypXVtif/mp8rV65Ikp566imn25F7Ta8+tcGjsvZZjG8eonsMlZWVaePGjdqzZ4/y\n8vLk5+en0NBQzZw5s9pdtys/zMpn/lW6fv26oqOjdfDgQRUXF6t79+6aMWOGRowYYep+PGk++eQT\nzZ07t9btTpw4IS8vL4fzd+rUKcXFxenkyZNV1xWPHTtW4eHh3NzEBHl5eYqJidGRI0dUWFiop556\nSsHBwZozZ061m2SQf83PD37wA0nSsWPHnG5H7jUPu3fv1tKlS7V06VJNnTq1Wp9hGEpKSlJSUpKy\ns7Pl6+urn/70p5o3b16Nm9UMGzZMV69e1aeffqrOnTtXtRcVFWnNmjX68MMPVVhYqK5du2rChAma\nMGEC89kI7M3fhQsXqs7sc2bv3r1Vl/k4mr+LFy8qNjZW//73v3X79m1ZrVa734XQMM7y7+bNm1q7\ndq0+/vhj5efny9vbW0OGDNH8+fNrFHfkn/mczZ10/0BXWlqavvjiC6e5Qu41vfrWBo/C2vfYF+IA\nAAAAADQnHEYDAAAAAMBEFOIAAAAAAJiIQhwAAAAAABNRiAMAAAAAYCIKcQAAAAAATEQhDgAAAACA\niSjEAQAAAAAwEYU4AAAAAAAmohAHAAAAAMBEFOIAAAAAAJiIQhwAAAAAABNRiAMAAAAAYCIKcQAA\nAAAATEQhDgAAAACAiSjEAQAAAAAwEYU4AAAAAAAmohAHAAAAAMBE/wMJ8hj3jU03sgAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x114a10410>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 333,
       "width": 497
      }
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DKL1 = 0.362, DKL2 = 2.504, EU = 0.947, FE = 0.770\n"
     ]
    }
   ],
   "source": [
    "# calc and plot p(a|w)\n",
    "pagw = pr.sum(['x'],pxgw*pagxw)\n",
    "plt.pcolor(pagw.val,cmap=\"Greys\")\n",
    "plt.show()\n",
    "\n",
    "# calc and print stats\n",
    "DKL1 = pr.sum(pw*pxgw*pr.log(pxgw/px)).val/np.log(2)             # DKL of the coordinating agent \n",
    "DKL2_x = pr.sum(pw*pxgw*pagxw*pr.log(pagxw/pagx)/px, ['x'])      # DKL2 as a function of x\n",
    "DKL2 = pr.sum(px*DKL2_x).val/np.log(2)                           # DKL2 averaged over all agents in the second step\n",
    "EU = pr.sum(pw*pagw*U).val                                       # expected utility\n",
    "FE = (EU - DKL1/beta1 - pr.sum(px*(DKL2_x/beta2))).val           # free energy\n",
    "print \"DKL1 = {:.3f}, DKL2 = {:.3f}, EU = {:.3f}, FE = {:.3f}\".format(DKL1,DKL2,EU,FE)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}