BA_serial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Serial 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": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dimensions\n",
    "N,M,K = 20,20,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 = 20.0    # first step, p(x) -> p(x|w)\n",
    "beta2 = 10.0    # second step, p(a) -> p(a|x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Blahut-Arimoto iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initializations\n",
    "px = pr.func(vars=['x'], val='unif').normalize()\n",
    "pa = pr.func(vars=['a'], val='unif').normalize()\n",
    "F = pr.func('f(w,x)', val='rnd')\n",
    "pxgw_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",
    "    pwgx = (pxgw*pw)/px\n",
    "    pagx = (pa*pr.exp(beta2*pr.sum(['w'],pwgx*U))).normalize(['a'])\n",
    "    pa = pr.sum(pxgw*pw*pagx,['a'])\n",
    "    F = pr.sum(['a'],pagx*(U-pr.log(pagx/pa)/beta2))\n",
    "    \n",
    "    # stopping condition\n",
    "    if np.linalg.norm(pxgw.val-pxgw_temp)<1e-15: break\n",
    "    pxgw_temp = pxgw.val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+IAAAKbCAYAAABxQDD2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAWJQAAFiUBSVIk8AAAIABJREFUeJzs3XmY13W9///HMAO4sJhLWiMEloJKKHJlLsdSL02PZpoJ\ngksoaZKJ0smlTY9L5/K0uUG4dZRQLNdjesrrcFmCemVHNJcUGa4QLoPCSA/K0MDMMJ/fHx34xo8B\nWT7zGidvt+vin/f6fF920dx5v+f9rqlUKpUAAAAARXTr7AEAAADgvUSIAwAAQEFCHAAAAAoS4gAA\nAFCQEAcAAICChDgAAAAUJMQBAACgICEOAAAABQlxAAAAKEiIAwAAQEFCHAAAAAoS4gAAAFCQEAcA\nAICC6qpxkCVLlmTixImZOXNm3njjjfTt2zcHHnhgLrjggvTr12+tbR988MFMmTIlCxYsSJ8+ffLP\n//zPOf/887Pttttu1Lmamppy88035+c//3lef/317Lrrrjn11FNzyimnpKamphqXAwAAAB2mplKp\nVLbkAEuWLMmIESPypz/9KQcffHAGDRqU+fPnZ8aMGenbt2/uvvvuDBgwIEly880355prrsmgQYPy\niU98InPnzs3MmTMzbNiwTJ06NT169NjguVatWpUvfelLmTlzZj75yU9mjz32yOOPP56GhoaMHTs2\nl1xyyZZcCgAAAHS4LQ7xyy67LHfffXe+9rWv5cwzz1yz/Gc/+1kuvvjiHHbYYbnpppuyaNGifOpT\nn8pHP/rR3HHHHenevXuS5Prrr8/kyZNz6aWX5rTTTtvguR5++OFceOGFa0V3S0tLzjrrrPzP//xP\nfvazn2XQoEFbcjkAAADQobb4d8QfffTRbL/99hkzZsxay48//vj0798/Tz75ZNra2nLPPfektbU1\n55xzzpoIT5Jx48alV69euffee9/xXNOmTUtdXV3GjRu3Zln37t0zYcKEVCqV3HfffVt6OQAAANCh\ntuh3xFetWpVzzjkndXV16dZt3abv0aNHWlpa0tramlmzZiVJ9t9//7W26dmzZ/bdd988+eSTWbZs\nWXr37t3uuZqbm/O73/0ugwcPTt++fddaN3To0Gy99dZrzgEAAADvVlsU4rW1tevcCV9t3rx5efXV\nV9O/f//06NEjr732Wnbcccd2X8pWX1+fJJk/f36GDh3a7vEWLVqU1tbW9O/fv905dtlllyxYsGDz\nLwYAAAAK6JDPl7W1teWqq65KW1tbRo4cmSRZunTpeu92r17e2Ni43mMuXbp0rW3bO0ZTU1NaW1u3\nZHQAAADoUFUP8UqlkssuuyxPPfVUhgwZsuaOeWtr63rfir56+cqVK9d73NWBvSXHAAAAgM5Wle+I\nr9ba2ppLL700DzzwQPr165fJkyevCeStttoqLS0t7e7X3NycJNl6663Xe+yePXsmyQaPUVNTs8Fj\nAAAAvBvV1NR09gibbAs/wPWeVrUQb2pqygUXXJCZM2dmwIABuf3227PzzjuvWd+nT58sW7as3X1X\nL1/fY+dJ1rygbX2Pry9btizbbLNNuy+NAwAAgHeLqlTrW2+9lTFjxmTmzJnZa6+9ctddd+WDH/zg\nWtsMGDAgb7zxRlasWLHO/osWLUq3bt3yoQ99aL3nqK+vT/fu3bNw4cJ11q1atSqLFy/OwIEDt/xi\nAAAAoANtcYivXLky55xzTl544YXsv//+ueOOO7LDDjuss93w4cPT1taWZ555Zp39n3/++XzkIx9J\nr1691nueurq67LPPPpk9e/Y6d8VffPHFNDU1ZdiwYVt6OQAAANChtjjEr7nmmjz33HMZNmxYbr31\n1vXG9Kc//enU1tZm0qRJa34nPEluuummNDY25uSTT37Hc51wwglpbm7OxIkT1yxraWnJ9ddfnyQZ\nMWLEFl4NAAAAdKyayhb8hv2SJUty2GGHpaWlJZ/73OfygQ98oN3tvvjFL6Znz575/ve/n1tvvTUf\n/vCHc9hhh+X3v/99ZsyYkf322y8//vGP13oj+pQpU7Js2bKMGTMmffr0SfK3R9BPPfXUPPfcczno\noIOy995754knnsicOXMyduzYXHLJJZt7KQAAAJ3Gy9reW7YoxB999NF8+ctffsftZs2alT59+qRS\nqeSuu+7KXXfdlddeey077bRTjjzyyJx33nnrvKjt8MMPz6JFi/LLX/4yu+6665rljY2NmThxYh55\n5JEsXbo0/fv3z+jRozN69GgvagMAALokIf7eskUhDgAAwJYT4u8tbiEDAABAQUIcAAAAChLiAAAA\nUJAQBwAAgIKEOAAAABQkxAEAAKAgIQ4AAAAFCXEAAAAoSIgDAABAQUIcAAAAChLiAAAAUJAQBwAA\ngIKEOAAAABQkxAEAAKAgIQ4AAAAFCXEAAAAoSIgDAABAQUIcAAAAChLiAAAAUJAQBwAAgIKEOAAA\nABQkxAEAAKAgIQ4AAAAFCXEAAAAoSIgDAABAQUIcAAAAChLiAAAAUJAQBwAAgIKEOAAAABQkxAEA\nAKAgIQ4AAAAFCXEAAAAoSIgDAABAQUIcAAAAChLiAAAAUJAQBwAAgIKEOAAAABQkxAEAAKAgIQ4A\nAAAFCXEAAAAoSIgDAABAQUIcAAAAChLiAAAAUJAQBwAAgIKEOAAAABQkxAEAAKAgIQ4AAAAFCXEA\nAAAoSIgDAABAQUIcAAAAChLiAAAAUJAQBwAAgILqOnsAAACA97qamprOHoGC3BEHAACAgoQ4AAAA\nFCTEAQAAoCAhDgAAAAUJcQAAAChIiAMAAEBBVf982euvv55jjjkm48ePzxlnnLFm+eGHH55FixZt\ncN+rr746J5544ga3ueiii/LQQw+1u+7ss8/OhRdeuMkzAwAAQClVDfHly5dn/PjxaWxsXGfd5z//\n+Sxbtmyd5StWrMhtt92Wnj175qMf/eg7nmPOnDnZcccdM2rUqHXWDR8+fPMGBwAAgEKqFuKLFi3K\n+PHj8/LLL7e7/u/vjv+9b3/722lra8s3v/nN7L777hs8R0tLS+bPn59DDz0048eP39KRAQAAoLiq\n/I74lClTctxxx2XOnDk54IADNnq/Z555JnfeeWcOPvjgjBgx4h23nzdvXlpaWjJo0KAtGRcAAAA6\nTVVCfOrUqamvr8+dd96Z448/fqP3+853vpPa2tp861vf2qjtGxoakkSIAwAA0GVV5dH0K664Igcd\ndFBqa2uzYMGCjdpn+vTpefHFFzNy5MjstttuG7XP6hCfP39+Ro0alYaGhmy11VY59NBDM2HChOy8\n886bewkAAABQRFXuiB9yyCGpra3dpH1uv/32dOvWLV/4whc2ep/VIT558uTsuuuuOfnkkzNgwIA8\n8MADGTFiRBYvXrxJMwAAAEBpVf982caYPXt2fvvb3+ZTn/pUBgwYsNH7bbXVVhkwYEAmTZq01ovd\nbrzxxlx33XX59re/nUmTJnXAxAAAAFAdnRLiDz74YJJk5MiRm7TfD3/4w3aXn3POObnvvvvy2GOP\nZfny5dl22223eEYAAADoCFV5NH1TPfbYY9luu+1y4IEHVuV43bp1y+DBg9Pa2urxdAAAAN7Viof4\nq6++mtdeey2HH3546uo2/oZ8U1NTnn/++cyZM6fd9StWrEiS9OzZsypzAgAAQEco/mj6Cy+8kCQZ\nPnz4Ju33l7/8JSeffHL22GOPPPzww2uta2pqyuzZs7P99tunvr6+arMCAABAtRW/Iz579uwkyd57\n771J+/Xr1y9777135s6dm4ceemjN8kqlkh/84Ad58803M3r06NTU1FR1XgAAAKim4nfE//CHPyRJ\n3v/+929wuylTpmTZsmUZM2ZM+vTpkyS58sorc/rpp+fiiy/O9OnTU19fn2eeeSYvvfRSPvaxj2Xc\nuHEdPj8AAABsieJ3xJcuXZok6d279wa3mzp1aiZNmpS33357zbIhQ4bkvvvuy1FHHZVZs2Zl2rRp\naWxszPnnn5/bbrstPXr06NDZAQAAYEvVVCqVSmcPAQAA8F7WrVunfNBqi7S1tXX2CF1W1/uvDQAA\nAF2YEAcAAICChDgAAAAUJMQBAACgICEOAAAABQlxAAAAKEiIAwAAQEG+Iw4AAIV1tW9G+140VFfX\n+hsAAAAAujghDgAAAAUJcQAAAChIiAMAAEBBQhwAAAAKEuIAAABQkBAHAACAgoQ4AAAAFCTEAQAA\noCAhDgAAAAUJcQAAAChIiAMAAEBBQhwAAAAKEuIAAABQkBAHAACAgoQ4AAAAFCTEAQAAoCAhDgAA\nAAUJcQAAAChIiAMAAEBBQhwAAAAKEuIAAABQkBAHAACAgoQ4AAAAFCTEAQAAoCAhDgAAAAUJcQAA\nAChIiAMAAEBBQhwAAAAKquvsAQAA4L2mUql09ghAJ3JHHAAAAAoS4gAAAFCQEAcAAICChDgAAAAU\nJMQBAACgICEOAAAABQlxAAAAKEiIAwAAQEFCHAAAAAoS4gAAAFCQEAcAAICChDgAAAAUJMQBAACg\nICEOAAAABQlxAAAAKEiIAwAAQEFCHAAAAAoS4gAAAFCQEAcAAICChDgAAAAUJMQBAACgoKqH+Ouv\nv57hw4dnypQp66y79957M2jQoHb/jBw5cqOO39TUlOuuuy5HHnlkhg4dmmOOOSbTpk1LpVKp8pUA\nAABA9dVV82DLly/P+PHj09jY2O76hoaGJMnZZ5+dnj17rrVul112ecfjr1q1KhdccEFmzpyZT37y\nkznqqKPy+OOP58orr8zChQtzySWXbPlFAAAAQAeqWogvWrQo48ePz8svv7zebRoaGrLddtvlwgsv\n3Kxz/OIXv8jMmTMzduzYNdF9wQUX5Kyzzsrtt9+eE044IYMGDdqsYwMAAEAJVXk0fcqUKTnuuOMy\nZ86cHHDAAevdbu7cudljjz02+zzTpk1LXV1dxo0bt2ZZ9+7dM2HChFQqldx3332bfWwAAAAooSoh\nPnXq1NTX1+fOO+/M8ccf3+42ixcvztKlSzf7jnVzc3N+97vfZfDgwenbt+9a64YOHZqtt946s2bN\n2qxjAwAAQClVeTT9iiuuyEEHHZTa2tosWLCg3W1W/354S0tLzj333Dz33HNZsWJF9ttvv1xwwQUZ\nOnToBs+xaNGitLa2pn///uusq62tzS677LLecwMAAMC7RVXuiB9yyCGpra3d4DarQ/ynP/1pVq5c\nmRNPPDEHH3xwnnrqqZxyyil54oknNrj/0qVLkyS9e/dud33v3r3T1NSU1tbWzbgCAAAAKKOqb03f\nkLa2ttTX12fChAn5zGc+s2b5008/nTPOOCNf//rX88tf/nKdt6mvtjqwe/To0e761ctXrlyZurpi\nlwUAAACbpOrfEV+fcePG5Ve/+tVaEZ4k+++/f4477rgsWbIkTz/99Hr3Xx3oLS0t7a5vbm5OTU1N\ntt566+oNDQAAHaCmpqZL/QGqq1iIb8hee+2VJFm4cOF6t1n9grb1faN82bJl2WabbdKt27vikgAA\nAKBdxar15ZdfXu9bzVeuXJkk630sPUnq6+vTvXv3dmN91apVWbx4cQYOHFidYQEAAKCDFAvxL3/5\ny/n85z+fN998c511zz77bJJkyJAh692/rq4u++yzT2bPnr3OXfEXX3wxTU1NGTZsWHWHBgAAgCor\nFuJHH3102tracu2116ZSqaxZ/sgjj2TGjBn52Mc+lj322GODxzjhhBPS3NyciRMnrlnW0tKS66+/\nPkkyYsSIjhkeAAAAqqTY68XPPffcPP7447nnnnvS0NCQ4cOHZ/78+ZkxY0Z22mmnXH311WttP2XK\nlCxbtixjxoxJnz59kiQnnnhi7r///kyZMiVz587N3nvvnSeeeCJz5szJ2LFjM2jQoFKXAwAAAJul\n2B3xPn365Kc//WnGjBmTJUuW5I477shLL72Uk046KQ888ED69eu31vZTp07NpEmT8vbbb69ZVltb\nmx/96Ec544wzMm/evEydOjWrVq3KZZddlosuuqjUpQAAAMBmq6n8/XPiAABAh+tqX/ppa2vr7BHg\nH0rX+hsAAAAAujghDgAAAAUJcQAAAChIiAMAAEBBQhwAAAAKEuIAAABQkBAHAACAgoQ4AAAAFCTE\nAQAAoCAhDgAAAAUJcQAAAChIiAMAAEBBQhwAAAAKEuIAAABQkBAHAACAgoQ4AAAAFCTEAQAAoCAh\nDgAAAAUJcQAAAChIiAMAAEBBQhwAAAAKEuIAAABQUE2lUql09hAAAP+ounXrWvc92traOnsEgH94\nXev/GQAAAKCLE+IAAABQkBAHAACAgoQ4AAAAFCTEAQAAoCAhDgAAAAUJcQAAAChIiAMAAEBBQhwA\nAAAKEuIAAABQkBAHAACAgoQ4AAAAFCTEAQAAoCAhDgAAAAUJcQAAAChIiAMAAEBBQhwAAAAKEuIA\nAABQkBAHAACAgoQ4AAAAFCTEAQAAoCAhDgAAAAUJcQAAAChIiAMAAEBBQhwAAAAKEuIAAABQkBAH\nAACAgoQ4AAAAFFTX2QMAAPwjq1QqnT0CAO8y7ogDAABAQUIcAAAAChLiAAAAUJAQBwAAgIKEOAAA\nABQkxAEAAKAgIQ4AAAAFVT3EX3/99QwfPjxTpkxZZ11jY2O++93v5sgjj8yQIUPy8Y9/POeee25e\neeWVjT7+RRddlEGDBrX75/vf/34VrwQAAACqr66aB1u+fHnGjx+fxsbGddY1NTXl1FNPzZw5czJs\n2LAcccQRWbx4caZPn54nn3wyt99+e4YPH/6O55gzZ0523HHHjBo1ap11G7M/AAAAdKaqhfiiRYsy\nfvz4vPzyy+2uv+OOOzJnzpycfvrp+da3vrVm+dNPP50zzjgjl19+eR5++OENnqOlpSXz58/PoYce\nmvHjx1drdAAAACimKo+mT5kyJccdd1zmzJmTAw44oN1tpk+fnpqamkyYMGGt5fvvv3/233//zJ07\nN6+//voGzzNv3ry0tLRk0KBB1RgbAAAAiqvKHfGpU6emvr4+V1xxRRYsWJDf/OY362wzatSovPHG\nG+nVq9c663r06JHkb4+2b0hDQ0OSCHEAAAC6rKqE+BVXXJGDDjootbW1WbBgQbvbnHTSSe0uf/PN\nN/PMM89km222ya677rrB86wO8fnz52fUqFFpaGjIVlttlUMPPTQTJkzIzjvvvEXXAQAAAB2tKo+m\nH3LIIamtrd2sfb/3ve9l+fLlOf7449fcGV+f1SE+efLk7Lrrrjn55JMzYMCAPPDAAxkxYkQWL168\nWTMAAABAKVV9a/qmmjx5ch544IHU19fnK1/5yjtuv9VWW2XAgAGZNGlSdt999zXLb7zxxlx33XX5\n9re/nUmTJnXkyAAAALBFOi3Er7/++kyePDnbbbddbr755vTt2/cd9/nhD3/Y7vJzzjkn9913Xx57\n7LEsX7482267bbXHBQAAgKqoyqPpm2LVqlX55je/mcmTJ2eHHXbIj3/847Xubm+Obt26ZfDgwWlt\nbfV4OgAAAO9qRe+INzc35/zzz89jjz2W+vr63HbbbRkwYMBG7dvU1LTm5WyDBw9eZ/2KFSuSJD17\n9qzmyAAAAFBVxUK8Uqnkq1/9ah577LHsvvvu+Y//+I9Nesv5X/7yl5x88snZY4898vDDD6+1rqmp\nKbNnz87222+f+vr6ao8OAAAAVVPs0fQ77rgj06dPz4c+9KFMnTp1kz811q9fv+y9996ZO3duHnro\noTXLK5VKfvCDH+TNN9/M6NGjU1NTU+3RAQAAoGqK3BFvbm7O5MmTkySDBg3KtGnT2t1u1KhR2Wmn\nnZIkU6ZMybJlyzJmzJj06dMnSXLllVfm9NNPz8UXX5zp06envr4+zzzzTF566aV87GMfy7hx40pc\nDgAAAGy2IiE+b968/O///m+SZPr06Zk+fXq72x1xxBFrQnzq1KlZtGhRPvvZz64J8SFDhuS+++7L\nDTfckN/85jeZMWNG6uvrc/755+fss89+x++QAwAAQGerqVQqlc4eAgDgH1VX+7U5PxoCdLziny8D\nAACA9zIhDgAAAAUJcQAAAChIiAMAAEBBQhwAAAAKEuIAAABQUJHviAMAvFd1tc+XAdDx3BEHAACA\ngoQ4AAAAFCTEAQAAoCAhDgAAAAUJcQAAAChIiAMAAEBBQhwAAAAKEuIAAABQkBAHAACAgoQ4AAAA\nFCTEAQAAoCAhDgAAAAUJcQAAAChIiAMAAEBBQhwAAAAKEuIAAABQkBAHAACAgoQ4AAAAFCTEAQAA\noCAhDgAAAAUJcQAAAChIiAMAAEBBQhwAAAAKEuIAAABQkBAHAACAgoQ4AAAAFCTEAQAAoCAhDgAA\nAAUJcQAAACioplKpVDp7CAD4R9CtW9f69+22trbOHgEA3pO61k8MAAAA0MUJcQAAAChIiAMAAEBB\nQhwAAAAKEuIAAABQkBAHAACAgoQ4AAAAFCTEAQAAoCAhDgAAAAUJcQAAAChIiAMAAEBBQhwAAAAK\nEuIAAABQkBAHAACAgoQ4AAAAFCTEAQAAoCAhDgAAAAUJcQAAAChIiAMAAEBBQhwAAAAKEuIAAABQ\nUNVD/PXXX8/w4cMzZcqUdtc/+OCDOeGEE7LvvvvmE5/4RK6++uosX758o4/f1NSU6667LkceeWSG\nDh2aY445JtOmTUulUqnSFQAAAEDHqWqIL1++POPHj09jY2O762+++eZccsklaWtry2mnnZbBgwdn\nypQp+cIXvpDm5uZ3PP6qVatywQUX5MYbb8zAgQPz+c9/PnV1dbnyyivz3e9+t5qXAgAAAB2iaiG+\naNGinH766XnhhRfWu/6GG27IsGHDcv/99+fCCy/MLbfcknPPPTfPPfdc7rnnnnc8xy9+8YvMnDkz\nY8eOzS233JILL7ww999/fw444IDcfvvtaWhoqNblAAAAQIeoSohPmTIlxx13XObMmZMDDjig3W3u\nueeetLa25pxzzkn37t3XLB83blx69eqVe++99x3PM23atNTV1WXcuHFrlnXv3j0TJkxIpVLJfffd\nt+UXAwAAAB2oKiE+derU1NfX584778zxxx/f7jazZs1Kkuy///5rLe/Zs2f23XffzJkzJ8uWLVvv\nOZqbm/O73/0ugwcPTt++fddaN3To0Gy99dZrzgEAAADvVlUJ8SuuuCIPPvhg9ttvv/Vu89prr2XH\nHXfMtttuu866+vr6JMn8+fPXu/+iRYvS2tqa/v37r7OutrY2u+yySxYsWLDpwwMAAEBBVQnxQw45\nJLW1tRvcZunSpendu3e761YvX99L3lbv//fbtneMpqamtLa2bszIAAAA0CmKfUe8tbU1PXr0aHfd\n6uUrV67c4P5/v+3mHAMAAAA6W12pE2211VZpaWlpd93qT5dtvfXW692/Z8+eSbLBY9TU1GzwGADQ\nkSqVSmePAAB0AcXuiPfp02e9L2NbvXx9j50nWfOCtvU9vr5s2bJss8026dat2CUBAADAJitWrQMG\nDMgbb7yRFStWrLNu0aJF6datWz70oQ+td//6+vp07949CxcuXGfdqlWrsnjx4gwcOLCqMwMAAEC1\nFQvx4cOHp62tLc8888xay1euXJnnn38+H/nIR9KrV6/17l9XV5d99tkns2fPXueu+IsvvpimpqYM\nGzasQ2YHAACAaikW4p/+9KdTW1ubSZMmrfmd8CS56aab0tjYmJNPPvkdj3HCCSekubk5EydOXLOs\npaUl119/fZJkxIgR1R8cAAAAqqjYy9o+/OEPZ+zYsbn11ltzwgkn5LDDDsvvf//7zJgxI/vtt19G\njhy51vZTpkzJsmXLMmbMmPTp0ydJcuKJJ+b+++/PlClTMnfu3Oy999554oknMmfOnIwdOzaDBg0q\ndTkAAACwWYqFeJJ89atfzQc+8IHcddddmTp1anbaaaecccYZOe+889b5LNnUqVOzaNGifPazn10T\n4rW1tfnRj36UiRMn5pFHHsmzzz6b/v3757LLLsvo0aNLXgoAAABslpqKb60AQFXU1NR09gibxI8A\nANA5fOsLAAAAChLiAAAAUJAQBwAAgIKEOAAAABQkxAEAAKAgIQ4AAAAFCXEAAAAoSIgDAABAQUIc\nAAAAChLiAAAAUJAQBwAAgIKEOAAAABQkxAEAAKAgIQ4AAAAFCXEAAAAoSIgDAABAQUIcAAAAChLi\nAAAAUJAQBwAAgIKEOAAAABQkxAEAAKAgIQ4AAAAF1XX2AADwj6KmpqazRwAAugB3xAEAAKAgIQ4A\nAAAFCXEAAAAoSIgDAABAQUIcAAAAChLiAAAAUJAQBwAAgIKEOAAAABQkxAEAAKAgIQ4AAAAFCXEA\nAAAoSIgDAABAQUIcAAAAChLiAAAAUJAQBwAAgIKEOAAAABQkxAEAAKAgIQ4AAAAFCXEAAAAoSIgD\nAABAQUIcAAAAChLiAAAAUJAQBwAAgIKEOAAAABQkxAEAAKAgIQ4AAAAFCXEAAAAoSIgDAABAQUIc\nAAAACqqpVCqVzh4CoKvp1q1r/TtmW1tbZ48AAMD/6Vo/SQIAAEAXJ8QBAACgICEOAAAABQlxAAAA\nKEiIAwAAQEFCHAAAAAqqK3myQYMGveM2U6dOzcc//vENbnPKKafk2WefbXfd5ZdfntGjR2/WfAAA\nANDRiob4eeed1+7yN954Iz/5yU+yww47ZLfddnvH4zQ0NGTgwIE59thj11k3ZMiQLZ4TAAAAOkrR\nEB8/fny7y7/0pS+lpqYm3/ve97LTTjtt8BgLFy5MY2NjTjrppPUeDwAAAN6tOv13xB966KH86le/\nyogRI3LwwQe/4/YNDQ1JNu4xdwAAAHi36dQQX7lyZa699tr07t07//Iv/7JR+whxAAAAurJODfG7\n7rorf/zjH3PWWWflfe9730bt09DQkJqamjz77LP57Gc/m3333Tef+MQn8m//9m9ZtmxZB08MAAAA\nW6amUqlUOuPEq1atyhFHHJG33norM2bMSJ8+fTZqv6OPPjrz589P7969c9RRR2WbbbbJrFmz8sor\nr2SPPfbIT37yk/Tq1auDpwfe67p16/Tf7NkkbW1tnT0CAAD/p+jL2v7er371q/zxj3/MmWeeudER\n3tbWlj6J9IwjAAAgAElEQVR9+mTPPffMzTffnJ133nnN8ssvvzx33313Jk6cmK9//esdOToAAABs\ntk67I/7lL385jz76aB555JGN+mTZO1m+fHkOPvjg9OrVK08++WQVJgRYP3fEAQDYXJ3yk+TKlSvz\n61//OnvssUdVIjxJtt122wwYMCBLlizJypUrq3JMAAAAqLZOCfGnn346f/3rX3PUUUdt0n5vv/12\nfvvb32b+/Pntrl+xYkW6deuWurpOe+IeAAAANqhTQvyFF15IkgwfPnyT9nv55ZczevTofOc731ln\n3Z///OcsXLgwe+65Z2pra6syJwAAAFRbp4T47NmzkyR77733Ju03fPjw7LTTTnn88ccza9asNcub\nm5tz1VVXpaWlJaeeempVZwUAAIBq6pSXtR133HF57bXX1twZX5+JEycmScaPH79m2WOPPZbzzjsv\nNTU1Ofroo7Pddtvl17/+debNm5djjz02P/jBD1JTU9Oh8wN4WRsAAJurU0L8n/7pn5LkHd9uPmjQ\noCRJQ0PDWsuff/75TJ48Ob/97W+zcuXKDBw4MCNHjswpp5zS5X44BrqmrvZ3jRAHAHj36LTPlwF0\nZUIcAIDN1bV+kgQAAIAuTogDAABAQUIcAAAAChLiAAAAUJAQBwAAgIKEOAAAABQkxAEAAKCgus4e\nAKArqlQqnT0CAABdlDviAAAAUJAQBwAAgIKEOAAAABQkxAEAAKAgIQ4AAAAFCXEAAAAoSIgDAABA\nQUIcAAAAChLiAAAAUJAQBwAAgIKEOAAAABQkxAEAAKAgIQ4AAAAFCXEAAAAoSIgDAABAQUIcAAAA\nChLiAAAAUJAQBwAAgIKEOAAAABQkxAEAAKAgIQ4AAAAFCXEAAAAoSIgDAABAQUIcAAAAChLiAAAA\nUJAQBwAAgIKEOAAAABQkxAEAAKAgIQ4AAAAF1XX2AABdUU1NTWePAABAF+WOOAAAABQkxAEAAKAg\nIQ4AAAAFCXEAAAAoSIgDAABAQUIcAAAAChLiAAAAUJAQBwAAgIKEOAAAABQkxAEAAKAgIQ4AAAAF\nCXEAAAAoSIgDAABAQUIcAAAAChLiAAAAUJAQBwAAgIKEOAAAABQkxAEAAKAgIQ4AAAAFCXEAAAAo\nSIgDAABAQXWlT3jdddflxhtvbHfdMccck2uvvXaD+y9dujQ33HBDZsyYkTfeeCMf/vCHc9ZZZ+WY\nY47piHEBAACgqoqH+Jw5c9KjR4988YtfXGfd7rvvvsF9//rXv2bs2LF55ZVXcvTRR+cDH/hApk+f\nnq985St58803c9ppp3XU2AAAAFAVNZVKpVLyhIcffnj69u2b//zP/9zkfW+66aZce+21ueyyy3Lq\nqacmSRobGzNq1KgsXLgwv/zlL7PDDjtUe2SAdXTr1rV+s6etra2zRwAA4P8U/UmysbExixYtyqBB\ngzZr/7vuuis77rhjRo0atWZZr169Mm7cuDQ1NeXhhx+u1qgAAADQIYqG+Jw5c5Jks0L8tddey+uv\nv57hw4entrZ2rXUf//jHkySzZs3a8iEBAACgAxX9HfGGhoYkyZtvvpkzzzwzL730UpLkwAMPzIQJ\nE7Lbbrutd9/XXnstSdK/f/911u20007p2bNnFixYUP2hAQAAoIqK3hFfHeK33XZbevXqlREjRmTo\n0KH57//+74wcOTKvvPLKevddunRpkqRPnz7tru/Vq1eWLVtW/aEBAACgioreEa+trU19fX2uvvrq\nNY+TJ8lDDz2Uiy66KN/4xjfW+xK31tbWJEmPHj3aXd+jR480NTVVf2gAAACooqIh/q//+q/tLv/M\nZz6Te+65J7Nmzcqrr77a7iPqPXv2TJI0Nze3e4zm5uZss8021Ru2i+hqb0Luam+aTtLl/nf117/+\ntbNHeE9YtWpVZ48AAEAX9a6por322itJsnDhwnbX9+3bN8nf3rzensbGxvTq1atjhgMAAIAqKXZH\nvLW1NbNnz06lUsk+++yzzvoVK1Yk+X93vv//BgwYkKT9UP/zn/+clStXZuDAgdUbGAAAADpAsRBv\na2vLKaeckm222SZPPfXUWp8gq1Qqee6551JXV5c999yz3f0/+MEP5oMf/GCeffbZtLW1rfWI89NP\nP50kGTZsWMdeBAAAAGyhYo+m9+jRI4cddljeeuut3HLLLWutu+222zJ37tx8+tOfXu9b0ZO//S75\n4sWLc+edd65Z1tjYmJtuuilbbbVVjj/++A6bHwAAAKqhplKpVEqdbOHChRk1alSWLFmSgw46KIMH\nD85LL72Up59+Oh/5yEdy55135n3ve1+SZOLEiUmS8ePHr9m/sbExn/vc57JgwYJ86lOfSr9+/TJ9\n+vT84Q9/yKWXXprTTjut1KW8a3hZW8fzsjbaU/Cvzqqoqanp7BEAAPg/RUM8SV5//fVcf/31efzx\nx7N06dK8//3vz1FHHZVzzz03vXv3XrPdoEGDkvy/b4+v9pe//CXXXHNNHnvssTQ1NWW33XbLF77w\nhRx77LElL+NdQ4h3PCFOe4Q4AACbq3iIU11CvOMJcdrT1f7qFOIAAO8eXa+KAAAAoAsT4gAAAFCQ\nEAcAAICChDgAAAAUJMQBAACgICEOAAAABQlxAAAAKEiIAwAAQEFCHAAAAAoS4gAAAFCQEAcAAICC\nhDgAAAAUJMQBAACgICEOAAAABQlxAAAAKEiIAwAAQEFCHAAAAAoS4gAAAFCQEAcAAICChDgAAAAU\nJMQBAACgICEOAAAABdV19gDwbtfW1tbZI/AuVFNT09kjAADQRbkjDgAAAAUJcQAAAChIiAMAAEBB\nQhwAAAAKEuIAAABQkBAHAACAgoQ4AAAAFCTEAQAAoCAhDgAAAAUJcQAAAChIiAMAAEBBQhwAAAAK\nEuIAAABQkBAHAACAgoQ4AAAAFCTEAQAAoCAhDgAAAAUJcQAAAChIiAMAAEBBQhwAAAAKEuIAAABQ\nkBAHAACAgoQ4AAAAFCTEAQAAoCAhDgAAAAUJcQAAAChIiAMAAEBBQhwAAAAKquvsAdgykyZN6uwR\nNsn555/f2SNssh/+8IedPQIAAPAPxB1xAAAAKEiIAwAAQEFCHAAAAAoS4gAAAFCQEAcAAICChDgA\nAAAUJMQBAACgoOLfEV+yZEkmTpyYmTNn5o033kjfvn1z4IEH5oILLki/fv3ecf9TTjklzz77bLvr\nLr/88owePbraIwMAAEDVFA3xJUuWZMSIEfnTn/6Ugw8+OMccc0zmz5+f//qv/8oTTzyRu+++OwMG\nDNjgMRoaGjJw4MAce+yx66wbMmRIB00OAAAA1VE0xCdOnJg//elP+drXvpYzzzxzzfKf/exnufji\ni/Pv//7vuemmm9a7/8KFC9PY2JiTTjop48ePLzEyAAAAVFXR3xF/9NFHs/3222fMmDFrLT/++OPT\nv3//PPnkk2lra1vv/g0NDUmSQYMGdeicAAAA0FGK3RFftWpVzjnnnNTV1aVbt3X7v0ePHmlpaUlr\na2t69OjR7jGEOAAAAF1dsRCvra1d5074avPmzcurr76a/v37rzfCk7+FeE1NTZ599tl861vfyvz5\n89OnT58cddRROf/889O7d++OGh8AAACqotM/X9bW1parrroqbW1tGTly5Aa3bWhoSKVSyQ033JC9\n9torI0aMyPbbb5+pU6fmlFNOSWNjY6GpAQAAYPMU/3zZ36tUKrnsssvy1FNPZciQIeu9Y578Ldj7\n9OmTPffcMzfffHN23nnnNcsvv/zy3H333Zk4cWK+/vWvlxofAAAANlmn3RFvbW3NN77xjdx7773p\n169fJk+evMHH0rt165Z77rknDz744JoIX738kksuydZbb52f//znJUYHAACAzdYpId7U1JRzzz03\nDzzwQAYMGJCpU6euFdebatttt82AAQOyZMmSrFy5soqTAgAAQHUVfzT9rbfeytlnn50XXnghe+21\nV370ox9lhx12eMf93n777fz+97/P+973vgwcOHCd9StWrEi3bt1SV9epT9sDAADABhW9I75y5cqc\nc845eeGFF7L//vvnjjvu2KgIT5KXX345o0ePzne+85111v35z3/OwoULs+eee6a2trbaYwMAAEDV\nFA3xa665Js8991yGDRuWW2+9Nb169drofYcPH56ddtopjz/+eGbNmrVmeXNzc6666qq0tLTk1FNP\n7YixAQAAoGqKPce9ZMmSTJs2LUmy22675dZbb213uy9+8Yvp2bNnJk6cmCQZP358kqRHjx656qqr\nct555+XMM8/M0Ucfne222y6//vWvM2/evBx77LE58cQTy1wMAAAAbKZiIf7CCy+kpaUlSXL//fev\nd7sxY8akZ8+emTRpUpL/F+JJcthhh2XatGmZPHlyZsyYkZUrV2bgwIG59NJLc8opp6SmpqZjLwIA\nAAC2ULEQP+KII9LQ0LDR269v23333Te33HJLtcYCAACAojrtO+IAAADwXiTEAQAAoCAhDgAAAAUJ\ncQAAAChIiAMAAEBBQhwAAAAKqqlUKpXOHoLN16NHj84eYZM0Nzd39gibrK2trbNH2CTduvn3NQAA\neDfzEzsAAAAUJMQBAACgICEOAAAABQlxAAAAKEiIAwAAQEFCHAAAAAoS4gAAAFCQEAcAAICChDgA\nAAAUJMQBAACgICEOAAAABQlxAAAAKEiIAwAAQEFCHAAAAAoS4gAAAFCQEAcAAICChDgAAAAUJMQB\nAACgICEOAAAABQlxAAAAKEiIAwAAQEFCHAAAAAoS4gAAAFCQEAcAAICChDgAAAAUJMQBAACgICEO\nAAAABQlxAAAAKEiIAwAAQEFCHAAAAAoS4gAAAFCQEAcAAICChDgAAAAUJMQBAACgICEOAAAABQlx\nAAAAKEiIAwAAQEFCHAAAAAoS4gAAAFCQEAcAAICChDgAAAAUJMQBAACgICEOAAAABQlxAAAAKEiI\nAwDw/7V3r0FRnWccwP8LCqgrIkpMWSBokkXUQs2OWGM7VdoGCwYtLXKzXoswCgqmxpjGNIW2NtMp\nCBgEnVGKiBVTZUhjGhLB2yQtclFjEJgqDJcERBkiqyvLsqcfnKVB9gbuHm7/38f3Pc/ue/bl4d3n\n7LkQEZGIWIgTERERERERiYiFOBEREREREZGIWIgTERERERERiYiFOBEREREREZGIWIgTERERERER\niYiFOBEREREREZGIRC/ENRoNcnJyEBgYCB8fH/z4xz/Ge++9h56eHrPiOzs7kZSUBH9/f/j6+iIk\nJARnz5618qiJiIiIiIiILEP0QjwpKQn79u2Dk5MT1q1bh1mzZiE9PR2vvfaaydiHDx9i06ZNOHHi\nBHx9fREVFYX79+8jMTEReXl5IoyeiIiIiIiI6OlIBEEQxHqzyspKREREICAgAGlpaZBIJBAEAW+8\n8QYKCwuRlZWF5cuXG4zPyspCamoq3n77bURFRQEAlEolwsPD0dzcjHPnzmHGjBli7c6IYGdnN9xD\nGBS1Wj3cQxg0rVY73EMYFBsbXnFCRERERDSSifqN/fjx4wCAuLg4SCQSAIBEIsHOnTshkUhw6tQp\no/H5+fmYOXMmwsPD+9qkUiliY2OhUqnwwQcfWG/wRERERERERBYgaiFeXl6O6dOnQy6X92ufNWsW\nPD09ceXKFYOxjY2NaGtrg0KhgK2tbb++xYsXA4DReCIiIiIiIqKRQLRCXK1Wo7W1FR4eHnr7ZTIZ\n7t+/j46ODr39jY2NAKA33sXFBfb29mhoaLDYeImIiIiIiIisQbRCvLOzEwAwdepUvf269q6uLqPx\njo6OevulUqnBWCIiIiIiIqKRQrRCXKPRADB8czFde3d395DjDcUSERERERERjRQTxHojBwcHADD4\nvHDd3bQnTZqkt9/e3r7fdvriJ0+e/LTDHHVG413IRxvehZyIiIiIiCxJtApDKpXCxsYGSqVSb7/u\ntHJDp65PmzYNAAzGK5VKSKVSC4yUiIiIiIiIyHpEK8Tt7Ozg6uqK5uZmvf3Nzc1wdnaGk5OT3n5P\nT8++7Z50584ddHd3Y/bs2RYbLxEREREREZE1iHrOrUKhQHt7O+rr6/u1t7W1oaGhAb6+vgZjXV1d\n4erqioqKCmi12n59ZWVlAICFCxdaftBEREREREREFiRqIb569WoAQGpqal8xLQgCUlJSAABhYWFG\n44ODg9Ha2oq8vLy+NqVSiaysLDg4OGDVqlVWGjkRERERERGRZUgEQRDEfMPExEScPXsWPj4+WLx4\nMaqqqlBeXo6AgACkpaVBIpEAADIyMgAA8fHxfbFKpRK/+MUv0NDQgFdeeQXu7u4oLi5GU1MT9u7d\ni7Vr14q5K0RERERERESDJnoh3tPTg0OHDuHMmTNoa2uDq6srgoODER0d3e/RZF5eXgCA2trafvF3\n795FSkoKSktLoVKpMGfOHGzevBlBQUFi7gYRERERERHRkIheiBMRERERERGNZ3xAMhEREREREZGI\nWIgTERERERERiYiFOBEREREREZGIWIgTERERERERiYiFOBEREREREZGIWIgTERERERERiWjCcA/A\n2jQaDfLy8lBQUIDm5ma4uLggJCQEW7ZswcSJE03Gd3Z2Ij09HefPn8e9e/fw/PPP49e//jUCAwNF\nGD21t7cjIyMDFy5cwL179zBt2jQsWbIEO3bsgLu7u8n4yMhIVFRU6O175513EBERYekh07fs378f\nBw8e1NsXGBiI1NRUo/HMP/F5eXmZ3CY3NxeLFy82ug1zT1xtbW0IDAxEfHw8NmzYMKC/sLAQOTk5\naGhogKOjI372s59h+/btmDJlilmvr1KpkJ2djQ8//BBtbW1wc3NDVFQUIiMjIZFILLw344+x+VMq\nlcjMzMQnn3yCr7/+GlOmTIFCoUB8fDy8vb3Nev1du3ahqKhIb190dDR+85vfPO0ujGvG5u/UqVN4\n66239Mb5+vqioKDA5Osz/6zH0Nz5+/ujpaXFaOy+ffsQEhJidBvmnuUNpjYY6WvfmC/Ek5KScPLk\nSSgUCvj7+6OyshLp6emora1Fenq60diHDx9i06ZNuHnzJlasWIHvfOc7KC4uRmJiIjo6OrB27VqR\n9mJ8am9vR2hoKL7++mssXboUgYGBqK+vxz//+U9cunQJJ0+ehKenp9HXqK2txezZsxEUFDSgb8GC\nBVYaOenU1NTAzs4OW7ZsGdD34osvGo1l/g2PuLg4ve337t3DiRMnMGPGDMyZM8fk6zD3xPPgwQPE\nx8dDqVTq7c/OzkZKSgq8vLywdu1a1NXVIScnB9euXUNubi7s7OyMvn5vby927NiBCxcu4Ec/+hEC\nAgJw8eJFJCUlobm5Gbt377bGbo0bxuZPpVIhKioKNTU1WLhwIX7yk5+gtbUVxcXFuHz5Mo4ePQqF\nQmHyPWpqajBz5kyEh4cP6DMnngwzlX+1tbUAHhdd9vb2/fqeffZZk6/P/LMeY3O3bt06dHV1DWh/\n9OgRjhw5Ant7e3z3u981+R7MPcsaTG0wKtY+YQyrqKgQ5HK5EB8fL2i1WkEQBEGr1Qqvv/66IJfL\nhZKSEqPxBw8eFORyuZCXl9fX1tXVJQQFBQm+vr7C3bt3rTr+8W7v3r2CXC4Xjhw50q+9sLBQkMvl\nQkxMjNH4pqYmQS6XC3/605+sOUwyYvny5cLq1auHFMv8G1liY2MFLy8v4fLlyya3Ze6Jp7m5Wfj5\nz38uyOVyQS6XC0ePHh3QP2/ePCEsLExQq9V97fv37xfkcrlw7Ngxk+9RVFQkyOVy4c9//nNfm1qt\nFtatWyd4eXkJNTU1Ftuf8cbU/GVnZwtyuVxITk7u1/6f//xH8Pb2FlauXGnyPdRqtTB//nxh27Zt\nlhw6CabnTxAEYe3atYKfn9+Q34P5Zx3mzJ0+ycnJglwuFwoKCkxuy9yzPHNrg9Gy9o3pa8SPHz8O\n4PEvPLrTByQSCXbu3AmJRIJTp04Zjc/Pzx9wFEsqlSI2NhYqlQoffPCB9QZP+PTTT+Hs7Iz169f3\na1+1ahU8PDxw+fJlaLVag/G6o9DmnGpLlqdUKtHS0jLkz5/5N3IUFRWhpKQEoaGhWLp0qcntmXvi\nyMnJwauvvoqamhp8//vf17tNQUEBNBoNYmJi+l2OFRsbC6lUanIdBB6vpRMmTEBsbGxf28SJE5GQ\nkABBEPD+++8//c6MQ+bMX3FxMSQSCRISEvq1+/n5wc/PD3V1dWhrazP6Prdu3UJPTw/z0cLMmT8A\nqKurg1wuH/L7MP8sz9y5e1J5eTny8vKwdOlShIaGmtyeuWd55tYGo2XtG9OFeHl5OaZPnz7gH+Cs\nWbPg6emJK1euGIxtbGxEW1sbFAoFbG1t+/Xpro00Fk9Pp7e3FzExMYiLi4ONzcA/Uzs7O/T09ECj\n0Rh8DRYDw6umpgbA0D5/5t/I0d3djdTUVEydOhU7d+40K4a5J47c3FzIZDLk5eVh1apVerfR5Ymf\nn1+/dnt7e3zve99DTU2N3tMvddRqNb744gvMnTsX06ZN69fn4+ODSZMmMReHyJz5Cw8PR2JiIqRS\n6YA+3WmVDx48MPo+zEfrMGf+Wltb0dnZOeTPnvlnHebMnT7vvvsubG1tDV7z/yTmnmUNpjYYLWvf\nmL1GXK1Wo7W1Fb6+vnr7ZTIZ6uvr0dHRAWdn5wH9jY2NAAAPD48BfS4uLrC3t0dDQ4NFx0z/Z2tr\nO+Bol86tW7dw+/ZteHh4GL2+o7a2FhKJBBUVFXjrrbdQX18PR0dHBAQEYPv27Zg6daq1hk/4/wLU\n0dGBjRs34saNGwCAJUuWICEhweh1xsy/kSM/Px9fffUVEhMTMX36dLNimHvi+P3vf4+XX34Ztra2\nBvOhsbERM2fO1HtjGplMBgCor6+Hj4+P3viWlhZoNBq9uWhra4tnn32WuThE5szfL3/5S73tHR0d\nKC8vx+TJk+Hm5mb0fXT/i+vr6xEeHo7a2lo4ODhg2bJlSEhIwKxZs55qP8Yrc+ZP99n39PRg69at\nqKqqwqNHj/DSSy9hx44dBvNOh/lnHebM3ZOKi4tx/fp1rFmzxqz7pADMPUsbTG0wWta+MfuLeGdn\nJwAY/MKnazd0NEQX7+joqLdfKpUaPZJC1qHVapGcnAytVos1a9YY3ba2thaCICA9PR3z5s1DaGgo\nnJ2dkZubi8jISIM3ViHL0C1AR44cgVQqRWhoKHx8fPDxxx9jzZo1uHnzpsFY5t/I0Nvbi9zcXEyZ\nMgWRkZFmxzH3xPHDH/5wwBkjT+rs7DS5DhqbD3PWUpVKZfTsJNLPnPkz5C9/+QsePHiAVatWmbzh\nkO5/cWZmJtzc3BAWFgZPT0+cPn0aoaGhaG1tHdIYxjtz5k/32f/9739Hd3c3QkJCsHTpUnz++eeI\njIzEpUuXjMYz/6xjKLl39OhR2NjYYPPmzWbHMPfEoa82GC1r35j9RVz3wRhaoHTt3d3dQ45XqVRP\nO0waBEEQ8Pbbb+Pzzz/HggULDB4VAx4npaOjI7y9vZGdnd131FGr1eKdd97ByZMnkZGRgT179og1\n/HHH1tYWMpkM+/bt6/eoq6KiIuzatQtvvvkmzpw5ozeW+TcylJSU4KuvvsLGjRsNHhR5EnNvZNFo\nNENeB3Xx397W2GtMmDBmv1KMKJmZmTh9+jRkMhkSExNNbu/g4ABPT08cOHCg39MqDh48iP379+MP\nf/gDDhw4YM0hj1tarRYymQwJCQkIDg7uay8rK8OGDRuwZ88enDt3bsDd1HWYfyNDdXU1Kisr8cor\nr5h8Ws+3Mfesz1BtMFrWvjH7i7iDgwOAx6cD6aNWqwEAkyZN0tuv+6eo205f/OTJk592mGQmjUaD\nN998E6dOnYK7uzsyMzON/gpgY2ODgoICFBYW9jv1x8bGBrt378akSZPw4YcfijH0cet3v/sdSkpK\nBjxvOjg4GIsWLUJ1dTVu376tN5b5NzIUFhYCgMmzT76NuTeyODg4DHkdBP6fi8ZeQyKRGH0Nspy0\ntDSkpaXByckJ2dnZA65d1Oe9997Dxx9/POCRkTExMXBzc0NpaanJ68xpaGJjY1FSUtKvCAceX7f6\n6quvor29HWVlZQbjmX8jw1DWQoC5Z23GaoPRsvaN2UJcKpXCxsbG4GkHutNaDZ1yoFvcDMUrlUq9\nN08hy1OpVNi6dStOnz4NT09P5ObmPtV1NVOmTIGnpyfa29uNHg0j65k3bx4AoLm5WW8/82/4dXd3\n47PPPoNcLjf7ejhTmHvic3R0NHgZh6l1EDCdi11dXZg8ebLeG+eQ5fT29uK3v/0tMjMzMWPGDPzt\nb38b8OV+sGxsbDB37lxoNBqeIjsMTK2DAPNvpCgtLYWTkxOWLFlikddj7j09U7XBaFn7xmzm2tnZ\nwdXV1eA/uObmZjg7O8PJyUlvv+7UE33xd+7cQXd3N2bPnm2x8ZJ+33zzDdavX48LFy5g3rx5yM/P\nh6urq8m4+/fvo7KyEvX19Xr7Hz16BBsbG57KZSUajQbXr1/HtWvX9PY/evQIAAyejsf8G35lZWV4\n+PAhAgICBhXH3BtZPD09ce/evb6c+7aWlhbY2NjgueeeMxgvk8kwceJEvbnY29uL1tZW5qKVqdVq\nbNu2De+//z5kMhny8/Mxd+5cs2JVKhWuXr3a9xSLJ5n6X0xP58svvzR4Z2XdwUhjnz3zb/jdvn0b\njY2N8Pf3H9S6xdyzHnNqg9Gy9o3ZQhwAFAoF2tvbB3whbGtrQ0NDg8E7qgOAq6srXF1dUVFRMeBZ\n1brTiBYuXGj5QVOf7u5uxMTE4Nq1a/Dz88OxY8cwY8YMs2K//PJLRERE4N133x3Qd+fOHTQ3N8Pb\n23vIN8oh47RaLSIjIxEdHY3e3t5+fYIgoKqqChMmTIC3t7feeObf8NMdRFEoFIOKY+6NLAqFAlqt\nFuXl5f3au7u7cfXqVbzwwgtGzy6ZMGECfH19UV1dPeCXgevXr0OlUjEXrUgQBLz22msoLS3Fiy++\niM6c/5QAAAUsSURBVBMnTgzqGtW7d+8iLCwMu3btGtCnUqlQXV0NZ2fnvrsIk2Vt27YN69atQ0dH\nx4C+iooKAMCCBQsMxjP/ht9Q10LmnnWYWxuMlrVvTBfiq1evBgCkpqb2fZkXBAEpKSkAgLCwMKPx\nwcHBaG1tRV5eXl+bUqlEVlYWHBwcBvXsQRq8lJQUVFVVYeHChTh8+PCgTkVWKBRwcXHBxYsX+x2N\nVqvVSE5ORk9PD6KioqwxbMLjM1KWL1+Ob775BocOHerXd+TIEdTV1WHlypVGbwDG/Bte1dXVAID5\n8+cPKo65N7KsXLkStra2OHDgQL97LmRlZUGpVJpcB4HHa6larUZGRkZfW09PD9LS0gAAoaGhlh84\nAQCOHTuG4uJiPPfcc0O6LMvd3R3z589HXV0dioqK+toFQcBf//pXdHR0ICIiAhKJxNJDJwArVqyA\nVqtFamoqBEHoa//oo49w/vx5LFq0CHK53OhrMP+G11DXQuaedZhbG4yWtW9Mnxv48ssvIzAwEGfP\nnkVYWBgWL16MqqoqlJeXIyAgAMuWLevbVvchx8fH97VFR0fjX//6F/74xz/iypUrcHd3R3FxMZqa\nmrB37169zx8ny2hvb8fx48cBAHPmzMHhw4f1brdlyxbY29sPmD87OzskJycjLi4OGzduxIoVK+Dk\n5ITPPvsMt27dQlBQEEJCQsTZmXFq9+7dqKqqwv79+1FWVoa5c+fixo0bKCsrwwsvvIA33nijb1vm\n38jT1NQEBwcHk3dLZ+6NbM8//zw2bdqEw4cPY/Xq1Vi+fDn++9//4vz583jppZcG3HwoJycHXV1d\nWL9+fd/ch4SE4B//+AdycnJQV1eH+fPn49KlS6ipqcGmTZvg5eU1HLs25qnVamRmZgIAvLy8+tbE\nJ4WHh8PFxQWA/vlLSkrCr371K7z++usoLi6GTCZDeXk5bty4gUWLFiE2NlacHRqHtm7diosXL6Kg\noAC1tbVQKBSor6/H+fPn4eLign379vXbnvk38jQ1NQEAnnnmGaPbMfesbzC1wWhZ+yTCtw/RjUE9\nPT04dOgQzpw5g7a2Nri6uiI4OBjR0dH97rqt+zB1z/zTuXv3LlJSUlBaWgqVSoU5c+Zg8+bNCAoK\nEnU/xptPP/0U27ZtM7ndlStX4OjoaHD+rl69iszMTFRWVvZdV7xmzRpERkby5iYiaGtrQ1paGi5e\nvIjOzk4888wzCAgIwNatW/vdJIP5N/L84Ac/AABcvnzZ6HbMvZHh9OnT2LNnD/bs2YMNGzb06xME\nAfn5+cjPz0djYyNcXFzw05/+FHFxcQNuVuPv74+WlhacO3cObm5ufe1KpRIZGRn46KOP0NnZCQ8P\nD0RERCAiIoLzaQH65u/mzZt9Z/YZU1hY2HeZj6H5u3XrFtLT0/Hvf/8bDx48gEwm0/tdiIbGWP7d\nv38fBw4cwCeffIL29nY4OTlh2bJl2L59+4DijvknPmNzBzw+0FVVVYUvvvjCaK4w96xvsLXBaFj7\nxnwhTkRERERERDSS8DAaERERERERkYhYiBMRERERERGJiIU4ERERERERkYhYiBMRERERERGJiIU4\nERERERERkYhYiBMRERERERGJiIU4ERERERERkYhYiBMRERERERGJiIU4ERERERERkYhYiBMRERER\nERGJiIU4ERERERERkYhYiBMRERERERGJiIU4ERERERERkYhYiBMRERERERGJiIU4ERERERERkYhY\niBMRERERERGJiIU4ERERERERkYj+B/odx+8LZD4ZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108882fd0>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 333,
       "width": 497
      }
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DKL1 = 2.732, DKL2 = 2.707, EU = 0.918, FE = 0.510\n"
     ]
    }
   ],
   "source": [
    "# calc and plot p(a|w)\n",
    "pagw = pr.sum(['x'],pxgw*pagx)\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 first step\n",
    "DKL2 = pr.sum(pw*pxgw*pagx*pr.log(pagx/pa)).val/np.log(2)   # DKL of the second step\n",
    "EU = pr.sum(pw*pagw*U).val                                  # expected utility\n",
    "FE = EU - DKL1/beta1 - DKL2/beta2                           # free energy\n",
    "print \"DKL1 = {:.3f}, DKL2 = {:.3f}, EU = {:.3f}, FE = {:.3f}\".format(DKL1,DKL2,EU,FE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}