estimate.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Derivation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Discretely many importance weights and rewards, maximum likelihood of sample $\\{ (w_i, r_i) \\}$ from $h$ is $$\n",
    "\\begin{alignat}{2}\n",
    "&\\!\\max_{Q \\succeq 0} &\\qquad& \\sum_n \\log(Q_{w_n, r_n}),\\label{eq:mle}\\\\\n",
    "&\\text{subject to} &  & \\vec{w}^\\top Q \\vec{1} = 1, \\tag{$\\beta$} \\label{eq:mlesumw} \\\\\n",
    "&                  &  & \\vec{1}^\\top Q \\vec{1} = 1. \\tag{$\\gamma$} \\label{eq:mlesum}\n",
    "\\end{alignat}\n",
    "$$\n",
    "Estimate is $\\hat V(\\pi) = \\vec{w}^\\top \\hat{Q} \\vec{r}$. \n",
    "\n",
    "Dual (ignoring constants) is $$\n",
    "\\begin{aligned}\n",
    "\\sup_{\\beta,\\gamma}& -\\beta - \\gamma + \\sum_{n} \\log\\left(w_n \\beta + \\gamma\\right)\\; \\text{ s.t. } \\; \\forall w,r: w \\beta + \\gamma \\geq 0.\n",
    "\\end{aligned}\n",
    "$$ One dual variable can be eliminated by summing the KKT stationarity conditions and leveraging complementary slackness.  Introducing $\\phi \\succeq 0$ as the (matrix of) dual variables associated with $Q \\succeq 0$: $$\n",
    "\\begin{aligned}\n",
    "\\frac{c_{w_i,r_j}}{q_{w_i,r_j}} &= \\phi_{w_i,r_j} + w_i \\beta + \\gamma \\implies n = 0 + \\beta + \\gamma, \\\\\n",
    "\\end{aligned}\n",
    "$$ resulting in the 1-D dual $$\n",
    "\\begin{aligned}\n",
    "\\sup_{\\beta} & \\sum_{n} \\log\\left((w_n - 1) \\beta + n\\right) \\; \\text{ s.t. } \\;\\forall w,r: (w - 1) \\beta + n \\geq 0.\n",
    "\\end{aligned}\n",
    "$$  This can be solved by 1-D bracketed search on the gradient followed by recovery of the primal values.\n",
    "\n",
    "Primary recovery begins with the primal-dual relationship for observed $(w, r)$ pairs: $$\n",
    "\\hat Q_{w,r} = \\sum_n \\frac{\\mathbb{1}_{w=w_n,r=r_n}}{\\beta^* (w_n - 1) + N}.\n",
    "$$  The MLE will sometimes put mass on unobserved importance weights, in which case the distribution over rewards for that importance weight is not determined.  The unobserved mass can be determined by solving the linear feasibility problem $$\n",
    "\\begin{alignat}{2}\n",
    "& &  & w_{\\min} \\hat{q}_{\\min} + w_{\\max} \\hat{q}_{\\max} = 1 - \\sum_n \\frac{w_n}{\\beta^* (w_n - 1) + N}, \\notag \\\\\n",
    "&                  &  & \\hat{q}_{\\min} + \\hat{q}_{\\max} = 1 - \\sum_n \\frac{1}{\\beta^* (w_n - 1) + N}, \\notag \\\\\n",
    "& & & {\\hat{q}_{\\min} \\geq 0, \\hat{q}_{\\max} \\geq 0},\\notag\n",
    "\\end{alignat}\n",
    "$$\n",
    "where $\\hat{q}_{\\min}$ and $\\hat{q}_{\\max}$ are associated with\n",
    "$w_{\\min}$ and $w_{\\max}$ respectively.  For robustness we convert this into a non-negative least squares problem $$\n",
    "\\begin{alignat}{2}\n",
    "&\\!\\min_{q_{\\min} \\geq 0, q_{\\max} \\geq 0} &\\qquad& \\left\\| \\left(\\begin{array}{cc} 1 & 1 \\\\ w_{\\min} & w_{\\max} \\end{array} \\right) \\left(\\begin{array}{c} q_{\\min} \\\\ q_{\\max} \\end{array}\\right) - \\left(\\begin{array}{c} 1 - \\sum_n \\frac{1}{\\beta^* (w_n - 1) + N} \\\\ 1 - \\sum_n \\frac{w_n}{\\beta^* (w_n - 1) + N} \\end{array} \\right) \\right\\|^2. \\notag\n",
    "\\end{alignat}\n",
    "$$\n",
    "When $q_{\\min} + q_{\\max} > 0$, the MLE is actually an interval; the center of this interval is found using $1/2 (r_{\\min} + r_{\\max})$ as the reward for unobserved importance weights.\n",
    "\n",
    "**Using a baseline:** When using a baseline, pass in shifted rewards and then add the correction to the result.  Given reward predictor $\\hat r: \\mathcal{X} \\times A \\to [r_{\\min}, r_{\\max}]$, construct data for the MLE $$\n",
    "\\begin{aligned}\n",
    "(w_n, \\tilde r_n) &\\leftarrow \\left(\\frac{\\pi(a_n|x_n)}{h(a_n|x_n)}, r_n - \\hat\n",
    "r(x_n, a_n) \\right),\n",
    "\\end{aligned}\n",
    "$$ apply the MLE on this data (with modified $\\tilde r_{\\min}$ and $\\tilde r_{\\max}$), and then adjust the result via $$\n",
    "\\begin{aligned}\n",
    "\\hat V^{\\text{(rpmle)}} &= \\hat V^{\\text{(mle)}} + \\sum_n \\sum_a \\pi(a_n|x_n) \\hat r(x_n, a_n).\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "**With censorship**: Suppose some $r_j = \\varnothing$ implying the reward was exogenously censored, and suppose we want to estimate $$\n",
    "\\begin{aligned}\n",
    "\\mathbb{E}\\left[r | r \\neq \\varnothing\\right] = \\frac{\\mathbb{E}\\left[r 1_{r \\neq \\varnothing}\\right]}{\\mathbb{E}\\left[1_{r \\neq \\varnothing}\\right]}.\n",
    "\\end{aligned}\n",
    "$$ One possible estimate is $$\n",
    "\\begin{aligned}\n",
    "\\hat V(\\pi) = \\frac{w^\\top Q (r 1_{r \\neq \\varnothing})}{w^\\top Q 1_{r \\neq \\varnothing}}\n",
    "\\end{aligned}\n",
    "$$ which is straightforward when there is no mass assigned to unobserved importance weights.  When there is mass assigned to unobserved importance weights, the MLE is again an interval and we can choose the center point of the interval as the estimate.\n",
    "\n",
    "In python we represent censored rewards with `None`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Limiting Behaviour"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $w_{\\max} \\to \\infty$ but $\\sum_n w_n < 1$, then $$\n",
    "\\begin{aligned}\n",
    "\\left. \\frac{\\partial}{\\partial \\beta} \\left( \\sum_{n} \\log\\left((w_n - 1) \\beta + N\\right) \\right) \\right|_{\\beta=0} &= \\sum_n \\frac{w_n - 1}{N} < 0,\n",
    "\\end{aligned}\n",
    "$$ therefore $\\beta^* < 0$, but $\\beta^* (w_{\\max} - 1) +n \\geq 0$ therefore $\\beta^* \\uparrow 0$.  Thus $$\n",
    "\\begin{aligned}\n",
    "\\hat V(\\pi) &\\to \\frac{1}{N} \\sum_n w_n r_n + \\left(1 - \\frac{1}{N} \\sum_n w_n\\right) \\rho\n",
    "\\end{aligned}\n",
    "$$ where $\\rho \\in [ r_{\\min}, r_{\\max} ]$.\n",
    "\n",
    "If $w_{\\min} = 0$ but $\\sum_n w_n > 1$, then additional mass does not affect the value estimate which is therefore a point estimate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Deviations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "\\hat V(\\pi) &= \\rho + \\frac{1}{N} \\sum_n \\frac{w_n (r_n - \\rho)}{1 + \\beta^* (w_n - 1)}, \\\\\n",
    "\\beta^* &= \\arg \\sup_{\\beta} \\left. \\left\\{ \\sum_n \\log\\left(1 + \\beta (w_n - 1) \\right) \\right| \\forall w: 1 + \\beta (w - 1) \\geq 0 \\right\\}.\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Lemma 1\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\left| \\beta^* \\right| \\sum_n \\frac{(w_n - 1)^2}{1 + \\beta^* (w_n - 1)} &\\leq \\left| \\sum_n (w_n - 1) \\right|\n",
    "\\end{aligned}\n",
    "$$ with equality for the unconstrained maximizer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "#### Proof"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "For the unconstrained maximizer, $$\n",
    "\\begin{aligned}\n",
    "0 &= \\sum_n \\frac{w_n - 1}{1 + \\beta^* (w_n - 1)} \\\\\n",
    "&= \\sum_n (w_n - 1) \\left( 1 - \\frac{\\beta^* (w_n - 1)}{1 + \\beta^* (w_n - 1)} \\right), \\\\\n",
    "\\beta^* \\sum_n \\frac{(w_n - 1)^2}{1 + \\beta^* (w_n - 1)} &= \\sum_n (w_n - 1), \\\\\n",
    "\\left| \\beta^* \\right| \\sum_n \\frac{(w_n - 1)^2}{1 + \\beta^* (w_n - 1)} &= \\left| \\sum_n (w_n - 1) \\right|.\n",
    "\\end{aligned}\n",
    "$$ For the constrained maximizer, first note the sign of $\\beta^*$ is the sign of $\\sum_n (w_n - 1)$, $$\n",
    "\\begin{aligned}\n",
    "\\left. \\sum_n \\frac{w_n - 1}{1 + \\beta (w_n - 1)} \\right|_{\\beta=0} &= \\sum_n (w_n - 1).\n",
    "\\end{aligned}\n",
    "$$\n",
    "If the constrained maximizer is positive than $$\n",
    "\\begin{aligned}\n",
    "0 &< \\left. \\frac{\\partial}{\\partial \\beta} \\sum_n \\log\\left(1 + \\beta (w_n - 1)\\right)\\right|_{\\beta=\\beta^*} \\\\\n",
    "&= \\sum_n \\frac{w_n - 1}{1 + \\beta^* (w_n - 1) } \\\\\n",
    "&= \\sum_n (w_n - 1) \\left( 1 - \\frac{\\beta^* (w_n - 1)}{1 + \\beta^* (w_n - 1)} \\right), \\\\\n",
    "\\beta^* \\sum_n \\frac{(w_n - 1)^2}{1 + \\beta^* (w_n - 1)} &< \\sum_n (w_n - 1), \\\\\n",
    "\\left| \\beta^* \\right| \\sum_n \\frac{(w_n - 1)^2}{1 + \\beta^* (w_n - 1)} &< \\left| \\sum_n (w_n - 1) \\right|.\n",
    "\\end{aligned}\n",
    "$$ If the constrained maximizer is negative than $$\n",
    "\\begin{aligned}\n",
    "0 &> \\left. \\frac{\\partial}{\\partial \\beta} \\sum_n \\log\\left(1 + \\beta (w_n - 1)\\right)\\right|_{\\beta=\\beta^*} \\\\\n",
    "&= \\sum_n \\frac{w_n - 1}{1 + \\beta^* (w_n - 1)} \\\\\n",
    "&= \\sum_n (w_n - 1) \\left( 1 - \\frac{\\beta^* (w_n - 1)}{1 + \\beta^* (w_n - 1)} \\right), \\\\\n",
    "\\beta^* \\sum_n \\frac{(w_n - 1)^2}{1 + \\beta^* (w_n - 1)} &> \\sum_n (w_n - 1), \\\\\n",
    "\\left| \\beta^* \\right| \\sum_n \\frac{(w_n - 1)^2}{1 + \\beta^* (w_n - 1)} &< \\left| \\sum_n (w_n - 1) \\right|.\n",
    "\\end{aligned}\n",
    "$$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lemma 2 (Freedman)\n",
    "\n",
    "For an arbitrary sequence of logging policies, $$\n",
    "\\begin{aligned}\n",
    "\\mathbb{E}\\left[\\frac{1}{N} \\left| \\sum_n (w_n - 1) \\right|\\right] &\\leq 5 \\sqrt{\\frac{w_{\\max}}{N}} + 4 \\frac{w_{\\max}}{N}.\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Proof\n",
    "\n",
    "$\\sum_n (w_n - 1)$ is a martingale with bounded predictable quadratic variation $\\mathbb{E}_{n-1}[(w_n - 1)^2] \\leq w_{\\max}$.  Hence Freedman's inequality applies."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "(* can't be done apparently *)\n",
    "a + 2 Integrate[Exp[-n t^2 / (2 s^2 + 2 c t / 3)], { t, a, Infinity }]\n",
    "$Aborted\n",
    "\n",
    "(* upper bound the bound by making the denominator bigger *)\n",
    "Reduce[2 s^2 + 2 c t / 3 <= 2 c t && t > 0 && s > 0, t, Reals] /. s > 0 -> True /. c > 0 -> True\n",
    "\n",
    "        2\n",
    "     3 s\n",
    "t >= ----\n",
    "     2 c\n",
    "\n",
    "2 Integrate[Exp[-n t / (2 c)], { t, 3 s^2 / (2 c), Infinity }] /. Re[n / c] > 0 -> True\n",
    "a + 2 Integrate[Exp[-n t^2 / (3 s^2)], { t, a, 3 s^2 / (2 c) }] \n",
    "Solve[D[%, a] == 0, a, Reals] /. n > 0 -> True // FullSimplify // PowerExpand // FullSimplify\n",
    "%%% + (%% /. %[[2]]) // FullSimplify // PowerExpand // FullSimplify\n",
    "% /. (Erf[(Sqrt[3] Sqrt[n] s)/(2 c)] - Erf[Sqrt[Log[2]]]) -> 1 // Collect[#, n]&\n",
    "(* in our case, s^2 <= c because of Radon-Nikodym *)\n",
    "%  // FullSimplify // List@@#& // { (#[[1]] /. s -> 0), (#[[2]] /. s -> Sqrt[c] // N)}& // Plus @@ #&\n",
    "\n",
    "4 c   4.51201 Sqrt[c]\n",
    "--- + ---------------\n",
    " n        Sqrt[n]\n",
    "\n",
    "(* can also not substitute s -> Sqrt[c] to get a bound based upon the quadratic variation *)\n",
    "\n",
    "2 Integrate[Exp[-n t / (2 c)], { t, 3 s^2 / (2 c), Infinity }] /. Re[n / c] > 0 -> True\n",
    "a + 2 Integrate[Exp[-n t^2 / (3 s^2)], { t, a, 3 s^2 / (2 c) }] \n",
    "Solve[D[%, a] == 0, a, Reals] /. n > 0 -> True // FullSimplify // PowerExpand // FullSimplify\n",
    "%%% + (%% /. %[[2]]) // FullSimplify // PowerExpand // FullSimplify\n",
    "% /. (Erf[(Sqrt[3] Sqrt[n] s)/(2 c)] - Erf[Sqrt[Log[2]]]) -> 1 // Collect[#, n]&\n",
    "%  // FullSimplify // List@@#& // { #[[1]], #[[2]] // N}& // Plus @@ #&\n",
    "\n",
    "       4 c           4.51201 s\n",
    "------------------ + ---------\n",
    "       2      2       Sqrt[n]\n",
    " (3 n s )/(4 c )\n",
    "E                n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bias Finite Sample Result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider estimator of the value difference with logging policy $$\n",
    "\\begin{aligned}\n",
    "\\hat{R}(\\pi) &= \\hat{V}(\\pi) - \\hat{V}(h) = \\rho + \\frac{1}{N} \\sum_n \\frac{w_n}{1 + \\beta^* (w_n - 1)} (r_n - \\rho) - \\left(\\rho + \\frac{1}{N} \\sum_n \\frac{1}{1 + \\beta^* (w_n - 1)} (r_n - \\rho)\\right) \\\\\n",
    "&= \\frac{1}{N} \\sum_n \\frac{(w_n - 1)}{1 + \\beta^* (w_n - 1)} (r_n - \\rho).\n",
    "\\end{aligned}\n",
    "$$ We can determine the bias of $\\hat{R}(\\pi)$ by taking differences with the doubly-robust estimator which is unbiased. $$\n",
    "\\begin{aligned}\n",
    "\\Delta &= \\frac{1}{N} \\sum_n \\frac{(w_n - 1)}{1 + \\beta^* (w_n - 1)} (r_n - \\rho) - \\frac{1}{N} \\sum_n (w_n - 1) (r_n - \\rho) \\\\\n",
    "&= \\frac{1}{N} \\sum_n (w_n - 1) (r_n - \\rho)  \\left(1 - \\frac{1}{1 + \\beta^* (w_n - 1)}\\right) \\\\\n",
    "&= \\frac{1}{N} \\sum_n (r_n - \\rho) \\frac{\\beta^* (w_n - 1)^2}{1 + \\beta^* (w_n - 1)}, \\\\\n",
    "|\\Delta| &\\leq |r_n - \\rho|_{\\infty} \\frac{1}{N} |\\beta^*| \\sum_n \\frac{(w_n - 1)^2}{1 + \\beta^* (w_n - 1)} \\\\\n",
    "&\\leq |r_n - \\rho|_{\\infty} \\frac{1}{N} \\left| \\sum_n (w_n - 1) \\right|. \\\\\n",
    "\\text{Bias} &= \\mathbb{E}\\left[ \\Delta \\right] \\leq \\mathbb{E}\\left[ \\left| \\Delta \\right| \\right] \\\\\n",
    "&\\leq |r_n - \\rho|_{\\infty} \\mathbb{E}\\left[ \\frac{1}{N} \\left| \\sum_n (w_n - 1) \\right| \\right] \\\\\n",
    "&\\leq |r_n - \\rho|_{\\infty} \\left( 5 \\sqrt{\\frac{w_{\\max}}{N}} + 4 \\frac{w_{\\max}}{N} \\right).\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tests"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Censorship changes results\n",
    "\n",
    "We learned this the hard way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "code_folding": [],
    "hidden": true,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.17414127154917453,\n",
      " {'betastar': -421.93139841688657,\n",
      "  'num': 159912,\n",
      "  'qfunc': <function estimate.<locals>.<lambda> at 0x7f70b41fbef0>,\n",
      "  'vmax': 0.276316821372124,\n",
      "  'vmin': 0.07196572172622505})\n",
      "(0.1522250873888096,\n",
      " {'betastar': -708.0158311345647,\n",
      "  'num': 268338,\n",
      "  'qfunc': <function estimate.<locals>.<lambda> at 0x7f7034069950>,\n",
      "  'vmax': 0.2276442757816804,\n",
      "  'vmin': 0.07680589899593883})\n"
     ]
    }
   ],
   "source": [
    "data, wmin, wmax, censored = None, None, None, None\n",
    "for data, wmin, wmax, censored in [\n",
    "    # some data where exogenous censorship is discarded\n",
    "   ([ (c, w, r) for c, w, r in [ \n",
    "      (86, -59.999996, 0.0), (44884, -1.0497237, 0.0), (16331, -1.0447762, 0.0), (31257, -1.0344828, 0.0), \n",
    "      (15868, -1.0, 0.0), (41332, 0.0, 0.0), (1958, 0.0, 1.0), (17763, 1.0, 0.0), (1339, 1.0, 1.0), \n",
    "      (30726, 1.0344828, 0.0), (3867, 1.0344828, 1.0), (2034, 1.0447762, 1.0), (16728, 1.0447762, 0.0), \n",
    "      (40629, 1.0497237, 0.0), (3445, 1.0497237, 1.0), (85, 59.999996, 0.0), (6, 59.999996, 1.0), \n",
    "     ]\n",
    "     if w >= 0\n",
    "    ], 0, 380, False),\n",
    "    # same data where exogenous censorship is modeled\n",
    "   ([ (c, -w if w < 0 else w, None if w < 0 else r) for c, w, r in [ \n",
    "      (86, -59.999996, 0.0), (44884, -1.0497237, 0.0), (16331, -1.0447762, 0.0), (31257, -1.0344828, 0.0), \n",
    "      (15868, -1.0, 0.0), (41332, 0.0, 0.0), (1958, 0.0, 1.0), (17763, 1.0, 0.0), (1339, 1.0, 1.0), \n",
    "      (30726, 1.0344828, 0.0), (3867, 1.0344828, 1.0), (2034, 1.0447762, 1.0), (16728, 1.0447762, 0.0), \n",
    "      (40629, 1.0497237, 0.0), (3445, 1.0497237, 1.0), (85, 59.999996, 0.0), (6, 59.999996, 1.0), \n",
    "     ]], 0, 380, True),\n",
    "]:\n",
    "    import MLE.MLE\n",
    "\n",
    "    from pprint import pformat\n",
    "    print(pformat(MLE.MLE.estimate(datagen=lambda: data, \n",
    "                                   wmin=wmin, wmax=wmax, rmin=0, rmax=1, raiseonerr=True, censored=censored)))\n",
    "  \n",
    "del data, wmin, wmax, censored"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison with CVX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "code_folding": [
     2
    ]
   },
   "outputs": [],
   "source": [
    "# CVXPY (primal) implementation\n",
    "\n",
    "class MLETest:\n",
    "    @staticmethod\n",
    "    def cvxestimate(data, wmin, wmax, rmin, rmax):\n",
    "        import cvxpy as cp\n",
    "        import numpy as np\n",
    "        from collections import defaultdict\n",
    "        \n",
    "        cdict = defaultdict(int)\n",
    "        n = 0\n",
    "        for (ci, wi, ri) in data:\n",
    "            assert ci >= 0\n",
    "            assert wi >= wmin and wi <= wmax\n",
    "            assert ri >= rmin and ri <= rmax\n",
    "            if ci > 0:\n",
    "                cdict[(wi, ri)] += ci\n",
    "            n += ci\n",
    "        assert n >= 1\n",
    "        cdict[(wmin, rmin)] += 0\n",
    "        cdict[(wmin, rmax)] += 0\n",
    "        cdict[(wmax, rmin)] += 0\n",
    "        cdict[(wmax, rmax)] += 0\n",
    "        cdict.default_factory = None\n",
    "        \n",
    "        wvec = np.array(list(set(w for (w, _), _ in cdict.items())))\n",
    "        wmaxvec = np.max(wvec)\n",
    "        rvec = np.array(list(set(r for (_, r), _ in cdict.items())))\n",
    "        C = np.array([ [ cdict.get((w, r), 0)/n for r in rvec ] for w in wvec ])\n",
    "        Q = cp.Variable((len(wvec), len(rvec)))\n",
    "            \n",
    "        prob = cp.Problem(cp.Maximize(cp.sum(cp.multiply(C, cp.log(Q)))), [\n",
    "                                cp.sum(cp.matmul((wvec/wmaxvec).T, Q)) == 1/wmaxvec,\n",
    "                                cp.sum(Q) == 1\n",
    "                          ])\n",
    "        prob.solve(solver='ECOS')\n",
    "            \n",
    "        vhat = 0\n",
    "        for i, wi in enumerate(wvec):\n",
    "            for j, rj in enumerate(rvec):\n",
    "                if cdict.get((wi, rj), 0) > 0:\n",
    "                    vhat += wi * Q.value[i, j] * rj\n",
    "                else:\n",
    "                    vhat += wi * Q.value[i, j] * 0.5 * (rmax - rmin)\n",
    " \n",
    "        from scipy.special import xlogy\n",
    "    \n",
    "        return vhat, { \n",
    "            'qstar': { (wvec[i], rvec[j]): Q.value[i, j] for i in range(len(wvec)) for j in range(len(rvec)) },\n",
    "            'likelihood': np.sum(xlogy(C, Q.value)),\n",
    "            'sumofone': np.sum(Q.value),\n",
    "            'sumofw': np.sum(wvec.dot(Q.value)),\n",
    "        }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "code_folding": [
     0
    ],
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 14/14 [02:04<00:00,  9.07s/it]\n"
     ]
    }
   ],
   "source": [
    "def testestimate():\n",
    "    from importlib import reload\n",
    "    from math import ceil\n",
    "    import environments.ControlledRangeVariance\n",
    "    import MLE.MLE\n",
    "    import numpy as np\n",
    "    from tqdm import tqdm as tqdm\n",
    "    import sys\n",
    "\n",
    "    reload(environments.ControlledRangeVariance)\n",
    "    reload(MLE.MLE)\n",
    "\n",
    "    wsupport = [ 0, 2, 20 ]\n",
    "    wmax = wsupport[-1]\n",
    "    env = environments.ControlledRangeVariance.ControlledRangeVariance(seed=45, wsupport=wsupport, expwsq=5)\n",
    "\n",
    "    for ndata in tqdm(map(ceil, np.logspace(1, 7, 14)), file=sys.stderr, total=14):\n",
    "        for i in range(1001):\n",
    "            (truevalue, data) = env.sample(ndata)\n",
    "            vhat, qstar = MLE.MLE.estimate(datagen = lambda: data, wmin=0, wmax=wmax, rmin=0, rmax=1, raiseonerr=True)\n",
    "            cvxvhat, cvxqstar = MLETest.cvxestimate(data, wmin=0, wmax=wmax, rmin=0, rmax=1)\n",
    " \n",
    "            from pprint import pformat\n",
    "            assert np.allclose(vhat, cvxvhat, atol=1e-4), pformat(\n",
    "            {\n",
    "                'data': [(c, w, r) for c, w, r in data if c > 0],\n",
    "                'vhat': vhat,\n",
    "                'cvxvhat': cvxvhat,\n",
    "                'qstar': qstar,\n",
    "                'cvxqstar': cvxqstar,\n",
    "            })\n",
    "                                    \n",
    "testestimate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 14/14 [02:04<00:00,  9.15s/it]\n"
     ]
    }
   ],
   "source": [
    "def megatestestimate():\n",
    "    from importlib import reload\n",
    "    from math import ceil\n",
    "    import environments.ControlledRangeVariance\n",
    "    import MLE.MLE\n",
    "    import numpy as np\n",
    "    from tqdm import tqdm as tqdm\n",
    "    import sys\n",
    "\n",
    "    reload(environments.ControlledRangeVariance)\n",
    "    reload(MLE.MLE)\n",
    "    \n",
    "    def getenv():\n",
    "        import numpy\n",
    "        wsupport = numpy.geomspace(0.5, 1000, 10)\n",
    "        env = environments.ControlledRangeVariance.ControlledRangeVariance(seed=45, wsupport=wsupport, expwsq=100)\n",
    "        return env, env.getpw(), env.range(), env.expectedwsq()\n",
    "\n",
    "    env = getenv()[0]\n",
    "    wmin, wmax = env.range()\n",
    "    \n",
    "    for ndata in tqdm(map(ceil, np.logspace(1, 7, 14)), file=sys.stderr, total=14):\n",
    "        for i in range(1001):\n",
    "            (truevalue, data) = env.sample(ndata)\n",
    "            vhat, qstar = MLE.MLE.estimate(lambda: data, wmin=wmin, wmax=wmax, rmin=0, rmax=1, raiseonerr=True)\n",
    "            try:\n",
    "                cvxvhat, cvxqstar = MLETest.cvxestimate(data, wmin=wmin, wmax=wmax, rmin=0, rmax=1)\n",
    "            except:\n",
    "                continue\n",
    "                \n",
    "            from pprint import pformat\n",
    "            assert np.allclose(vhat, cvxvhat, atol=1e-4) or not np.isfinite(cvxqstar['likelihood']), pformat(\n",
    "            {\n",
    "                'data': [(c, w, r) for c, w, r in data if c > 0],\n",
    "                'vhat': vhat,\n",
    "                'cvxvhat': cvxvhat,\n",
    "                'qstar': qstar,\n",
    "                'cvxqstar': cvxqstar,\n",
    "            })\n",
    "                                    \n",
    "megatestestimate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison with baselines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "code_folding": [
     0,
     43,
     50
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "****** Constant 0.5 ******\n",
      "****** ClippedDR ******\n",
      "****** SNIPS ******\n",
      "****** EMP ******\n",
      "****** MLE ******\n"
     ]
    }
   ],
   "source": [
    "def produceresults(env, method, maxexp=5, numpts=20, ndataperpt=10000):\n",
    "    from math import ceil\n",
    "    import numpy as np\n",
    "    \n",
    "    wmin, wmax = env.range()\n",
    "\n",
    "    for ndata in map(ceil, np.logspace(1, maxexp, numpts)):\n",
    "        estimates=[]\n",
    "        for i in range(1, ndataperpt+1):\n",
    "            (truevalue, data) = env.sample(ndata)\n",
    "            try:\n",
    "                estimate = None\n",
    "                estimate = method(data=data, wmin=wmin, wmax=wmax)\n",
    "                assert np.isfinite(estimate)\n",
    "            except:\n",
    "                print('truevalue was {}'.format(truevalue))\n",
    "                print('data was {}'.format(data))\n",
    "                print('estimate was {}'.format(estimate))\n",
    "                raise\n",
    "            \n",
    "            essden = sum(c*w*w for (c, w, _) in data)\n",
    "            essnum = sum(c*w for (c, w, _) in data)\n",
    "            ess = 0 if essden == 0 else essnum*(essnum/essden)\n",
    "                                                \n",
    "            estimates.append(\n",
    "                ( truevalue,\n",
    "                  truevalue - estimate,\n",
    "                  (truevalue - estimate)**2,\n",
    "                 ess\n",
    "                )  \n",
    "            )\n",
    "            \n",
    "        yield (ndata,\n",
    "                { \n",
    "                    'bias': np.abs(np.mean([ x[1] for x in estimates])),\n",
    "                    'biasstd': np.std([ x[1] for x in estimates ], ddof=1) / np.sqrt(len(estimates) - 1),\n",
    "                    'mse': np.mean([ x[2] for x in estimates ]),\n",
    "                    'msestd': np.std( [ x[2] for x in estimates ], ddof=1) / np.sqrt(len(estimates) - 1),\n",
    "                    'ess': np.mean([ x[3] for x in estimates ]),\n",
    "                    'essstd': np.std([ x[3] for x in estimates ], ddof=1) / np.sqrt(len(estimates) - 1),\n",
    "                },\n",
    "              )\n",
    "        \n",
    "class ClippedDR:\n",
    "    @staticmethod\n",
    "    def estimate(data, baseline=0.5, **kwargs):\n",
    "        import numpy as np\n",
    "        n = sum(c for c, _, _ in data)\n",
    "        return baseline if n == 0 else np.clip(sum(c*w*(r-baseline)+c*baseline for c, w, r in data) / n, a_min=0, a_max=1)\n",
    "    \n",
    "class SNIPS:\n",
    "    @staticmethod\n",
    "    def estimate(data, **kwargs):\n",
    "        effn = sum(c*w for c, w, _ in data)\n",
    "        return 0.5 if effn == 0 else sum(c*w*r for c, w, r in data) / effn\n",
    "    \n",
    "class EMP:\n",
    "    @staticmethod\n",
    "    def estimate(data, **kwargs):\n",
    "        vhat, mledict = MLE.MLE.estimate(datagen=lambda: data, **kwargs)\n",
    "        qstar = mledict['qfunc']\n",
    "        q = [(qstar(c,w,r), w, r) for c, w, r in data]\n",
    "        effn = sum(qi for qi, wi, ri in q)\n",
    "        return 0.5 if effn == 0 else sum(qi*wi*ri for qi, wi, ri, in q) / effn\n",
    "\n",
    "from importlib import reload\n",
    "import environments.ControlledRangeVariance\n",
    "import MLE.MLE\n",
    "\n",
    "reload(environments.ControlledRangeVariance)\n",
    "reload(MLE.MLE)\n",
    "\n",
    "def getenv():\n",
    "    wsupport = [ 0, 2, 1000 ]\n",
    "    env = environments.ControlledRangeVariance.ControlledRangeVariance(seed=45, wsupport=wsupport, expwsq=100)\n",
    "    return env, env.getpw(), env.range(), env.expectedwsq()\n",
    "\n",
    "allres = []\n",
    "for (name, method) in [ ('Constant 0.5', lambda **kwargs: 0.5),\n",
    "                        ('ClippedDR', ClippedDR.estimate),\n",
    "                        ('SNIPS', SNIPS.estimate),\n",
    "                        ('EMP', EMP.estimate),\n",
    "                        ('MLE', lambda data, **kwargs: MLE.MLE.estimate(datagen=lambda: data, **kwargs)[0]),\n",
    "                      ]:\n",
    "    print('****** {} ******'.format(name))\n",
    "    res = []\n",
    "    for zzz in produceresults(getenv()[0], method, numpts=14):\n",
    "        res.append(zzz)\n",
    "#         print('{}'.format(zzz), flush=True)\n",
    "    wmax = getenv()[2][1]\n",
    "    allres.append((name, [(x[0] / wmax, x[1]) for x in res]))\n",
    "    del wmax\n",
    "import pickle\n",
    "pickle.dump( allres, open( \"epsilongreedy_estimate_allres.p\", \"wb\" ) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "code_folding": [
     2
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaIAAAEfCAYAAAADX8LEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3wU1d7H8c+ZLem9BwgJoAgKiICAcBUVEVBAEUQfFZUOInrVey0XEQQUFdu1ISAoClYUBQtSFRG9giJiQEJLCIT0nmy2neeP2YQQEklCkk0579cr7GZ2dve3S3a+e86cOSOklCiKoiiKu2juLkBRFEVp2VQQKYqiKG6lgkhRFEVxKxVEiqIoilupIFIURVHcSgWRoiiK4lYqiBRFURS3arFBJIR4XAhxQAjhFELc4O56FEVRWiqjuwtwow3ASmBZTe4UGhoqY2Nj66UgRVGU5mrXrl0ZUsqwym5rFEEkhGgNPAz0BLoBXkCclPJoJeu2AV4ErgEEsBG4X0qZVJPnlFL+5Hq8GtUaGxvLzp07a3QfRVGUlk4IkVjVbY2la64DcDOQDWyraiUhhDewGbgAuBO4AzgP2CKE8GmAOhVFUZQ61ihaRMD3UsoIACHEBGBQFetNBNoBHaWUB13r7wESgMnAC65lG4GLq3iMEVLK7XVYu6IoinIOGkUQSSmd1Vx1OPBTaQi57ntECLEdGIEriKSUA+u+SkVRFKU+NJauueq6ENhbyfI/gc719aRCiElCiJ1CiJ3p6en19TSKoigtUlMLomD0/UgVZQFBNXkgIcRsIUQy0BdYKoRIdg2aOIOUcrGUsqeUsmdYWKWDPhRFUZRaampBVGeklLOllK2llB5SylDX9WR316UoitLSNLUgyqbylk9VLaU6I4QYJoRYnJubW59PoyiK0uI0isEKNfAn+n6iijoD8fX5xFLKtcDanj17TqzP56kJp1NSYndSYnfolzYnEomH0YCHUcNs1PAwahgNTe37RtMnpcThlDhcl3anxOHQL52y/O/OsvXsDond4cTqcGJz3SYlmA36/6XZqGHQBCaDfmnUBEaDhlET+nJNw2BwLXctq+lxckrDkVIiJTilxFl2eeq6dOqXDtfy09Z1Vli33HWTQcNs0PRLo4bJIPS/H4PWaP8emloQfQEsFEK0k1IeBhBCxAL9gEfcWNdZZRdaWbvnBJkFVoptDixWBxa7A4tNDxGLzYnVFShWhxOrXd8Y6Zf6j93p2li5Nl7Oap7lXYC+4TIIjJqG0aBvtEwGfaNW9gdrFHgYDJhNeoB5GA14mjQ8jQa8zAYAbKX1OZzYHdK10dSvl106nWV12h2nNsRlG92KP1KC1AsV6B+U0s+L4NRBx0JA6cdICNeaosI6ZeuJU7ehP7xE/0fq18rePyn1K6duP3O5LF0Xeep6ufuf2iC4Lqv3X1PvBKBpAk2AJvRgMojSZfp7dtrr5/TXStn1ireVW1r+PSu3WAj9704Tp4KxNEANpX+Lrr9Lk0ErC9byf5f636Y4bcNqdq1TXxtVIcBml5TYHWWfRatD/6Jnc31J0P/+T30OTvt8ll4v9zdud0qcztNDxR1K/w6MmkAr94WldLtQ8f/AaDj13nuaDFzdKZzberet87oaTRAJIUa5rvZwXQ4RQqQD6VLK71zLlgDTgc+FEDPR/+TnAseANxuy3po6mWdh1ud/nrG8/IfTqJX+52uYXX8Ivp5GzAbttBaOp8lwKiRMBrxMBrzMGl4m/b/TYndQbNVDzmJ3UmI7FXilLSf9A+bAatc/VEVWO7YKIVL6AarsQ1O6gTMIgabhunRtbMSZGxwPo4aPZjgtDI2aHn4mTUMIfSNXuqF3VtjIA+U+wPqy8uuX3a7fjJNT3zhLU06UC63Trpe+JleACVE+9M68n1YhGEs/yAYN/dJQ+gHXMFX4AqC3aFwfdteXAUPZ++DaABs1/b0x6LVZXRu68l9K9N/1/59Ty/SNYOlGsfSLgd1Z/kuMa8PoWqfsdZZ/Hyp5jzjtPTr1vpV9YSj/3rm+AThKa7M7y9V56guK3Smx2B04SsrVVW6DXboRd8rqf+mqL5o49Vkt/fs2nLERd/2dexgwaq7ANJaGqsBsMGA26v/fBk3/O9JEacu13GfItbz0C0Tp85SGSPl1DAbXOkJ/TAGUuEKzfG+J1eZ0BasesKV/PxX/pmyu/5diq4N8p/20vx+nlAR5m+sliIR0VzRXIISoqpDvpJQDyq0Xw+lT/GxCn+LnaD3XNwwY1qFDh4kJCQk1vr/V7iQt36J3m7laHI25qVye3aH/UQNlrSlNa/x1K82HdHVn2su1Rqp93xo/F2XdWaUttabwOW3shBC7pJQ9K72tsQRRU9GzZ0+p5ppTFEWpmb8LIrUXW1EURXErFUSKoiiKW6kgqiZ1HJGiKEr9UEFUTVLKtVLKSQEBAe4uRVEUpVlRQaQoiqK4lQoiRVEUxa1UECmKoihupYKomtRgBUVRlPqhgqia1GAFRVGU+qGCSFEURXErFUSKoiiKW6kgUhRFUdxKBZGiKIriViqIqkmNmlMURakfKoiqSY2aUxRFqR8qiBRFURS3UkGkKIqiuJUKIkVRFMWtVBA1oO3Ht6NOza4oinI6FUQNZOfJnUzZOIXpm6erMFIURSlHBVE1nevw7R4RPbiyzZV8n/w992y6R4WRoiiKiwqiajrX4dtCCF6+8mV6RvRk2/FtKowURVFcVBA1ICEEy65dVhZGUzdOVWGkKEqLp4KogZUPo+0ntqswUhSlxVNB5AYVw2jKhikqjBRFabFUEDWkkvyyq6Vh1CuiFz+m/MjkDZNVGCmK0iKpIGooqfHwQifY8nTZIiEESwctpVdEL3ak7GDShkkqjBRFaXFUEDUUv0gIbAvfLYBVY8DpAEDTNJYOWsqlkZfyU8pPKowURWlxVBA1FO9g5ITNEHcFHPgGXusNBemAHkZLrllSFkYTv52owkhRlBZDBVEDsaWlcfiGm8hv+wD0mQaZB+HVnpC8Ezg9jH4++TMTvp2gwkhRlBZBBVE1nevMCs7CQhw5OSRPu4fMkx3hhjfAVgTLBsOud4DTw+h/J/+nwkhRlBZBBVE1nevMCh5xccR+/BHG8HDSFjxDyqf7kHevBw8/WHsfrP0nSHlGGI1fPx6n01nHr0ZRFKXxUEHUgMzR0bT7+is8LuxMzkcfk/ToSzgn/Qih58GuZbBsENiKTwujX1J/Yfy3KowURWm+VBA1MIOXF3Eff4zftddS9MsvHL7lbuw3fwkdh8Cx/8Erl0Dm4bIw6h3Zm52pOxn37TgVRoqiNEsqiNxAaBqtX36JkClTsCUlcXj4CCzdZ8Pl/4b8k7CoPxz4Fk3TWHzNYnpH9WZX6i7uXn+3CiNFUZodFURuFH7/fUQvWICjsJCjY24hn74w6m3ACe/fAt8/r4fRQD2Mfk37VYWRoijNjgoiNwsYMZy2K95BmEwkT7+XzJ+zYMIW8A6BzU/Ch3eggR5GkXoY3bX+LhVGiqI0GyqIGgHv7t2J+/xzfUTds89y4r/vIaf/ApFdYd8XsOgytOJsFl+zmD5Rffgt7Tfu/OZOHK7ZGRRFUZoyFUSNhLmVPqLO86ILyV39KYkT78V55wa4aBSk7YPXeqGl/M6bA9+kT1Qfdqfv5q5v7lJhpChKk6eCqBExeHkR+/HH+A2+luKdOzk8fDj2Ac/AwDlgyYVl16Lt+eC0MFItI0VRmjoVRI2MEILWL71E6LSp2JKTOXTd9RSHXAu3fgAGE3x+D9r6R8vC6Pf037nzmzuxO+3uLl1RFKVWVBA1UmEzZhD9zAJkURGJt9xKXqKAydvALxp+XoT2zjDevOJl+kb11cPoaxVGiqI0TSqIqulc55qrjYDhw4lxjag7PuM+Mj7dBPf8D9pcCok/oL3Rm0U9HqVvVF/2ZOxh7NdjsTlsDVafoihKXVBBVE3nOtdcbXl3707cF59jjIgg/fkXOPHEPOS4b+GSsZBzDG3x5SyKvZm+UX35I+MPxqwbQ2ZxZoPWqCiKci5UEDUB5uho2n31JZ4XXUTumjUk3nYbzmsXwtDnwG5B+2AMi7w60S+6Hwk5Cdz4+Y38kvKLu8tWFEWpFhVETYQ+ou4j/IYMpvjX3zh83fXY2t0Id3wGZl+0TbN5I8fKzeffTK41lykbp7BkzxJ3l60oinJWKoiaECEErV98kZBpU7EdP87hoddRXBAIU36A4HaIP1fz+G9fMqfLdAyagVd+e4UZm2dgsVvcXbpSG/YSsKn/O6X5E+rEazXTs2dPuXPnTneXQe7ataQ89h8QgqhnnyXg6ivg/TFweCuYfTnQ7x6mpW0mtSiVWP9YXr36Vdr6t3V32crZFKTDiV/h+K+Q9idoJugwEM4fDD4h7q5OUWpNCLFLStmz0ttUENVMYwkigKLff+fY+Ak4i4oInTGDsCmTYesz8MPz4LBSGDeAyUGe/J71J35mP2b1mcXguMHuLlspz+mEjAOu8NkFucnlbnMdqKwZQBggpjdccD2EtHdPrYpyDlQQ1aHGFEQAtpQUjv7fbdhTUvAfPpzoZxYgknfCx3dBXjJO3yievvAffJz+C5rQuOWCW3igxwMYNaO7S2+5rIVwYrcePid2g7Xg1G1OBxSm66cDKUzXl/lFQkAb8AwAISD0fD2Q2lyq/64oTYAKojrU2IIIwFFcTNLYsVj+2Itn167ELHsLg9EJqyfAwQ2gGfm64xXMth2jyF5Er4hePH/F8wR5Bbm79JYjN1nvbjvxK6QfAFluWqaK4SMd+iwavpEgJeSn6Ms8/CGgNfhHg2YEnzD9hIrtrwKTl/tem6JUgwqiOtQYgwhASsnxhx4i/8uvMISGELN0KZ4dO8JPb8CW+WAtIDmsA3eH+HKyJIton2ievfxZuoV3c3fpzZPDDmnxenfbid+gIPX02ysNHzP4RugtIK8gEK6xRA6bHkY5x8Car3fT+Ue7Wkn+YPKG9lfC+UPAN6zhX6uiVIMKojrUWIOoVOay5aS98AJC04icPZvAkTdCyh74ZBxkJmD18Gda7Hn8bE3Hx+TDjO4zuPWCWxGqi+fcFefooXPiV/09rzhasSbhA+ATCgYPyDuu/y6lPvlt7jFXK8kJHgEQ2Ea/v8Gsd9ddcD2Entdwr1tRqkEFUR1q7EEEULhzJ8nT7sGZl0fg6FFEPvkkwloAX9wH+75ASieLWndgsckKCIbEDWFmn5l4m7zdXXrTk5cCidv1bresw0CFz5PTDoUZfxM+wafv5zGYoE1vaDcAIi7Sb8s6Akd/0J+nOFtfz2HTAyo3Wd/HpBn1eQgD24CHH4R0cO1H6g2aOkpDcT8VRHWoKQQRgD0ri8Q778SacBCPzp1pu3wZBn9/2PUObJoDxVnsDAhnRmgA+c4SOod05rnLnyPGP8bdpTcNtmL44xM48I0eNuWVhU+Kfnm28AG9BdNuAMRcBuYqvhBICal79VA69rNeg5RgyYGcJL37TzrBM1APJN9I8IvQh363vwrMPvXxTihKtaggqkNNJYhA32904l//Jm/dOgxBQbRe/CbeXbpAajx8NgVO7iHXZObO1q05hI0wrzAe7/s4V7a50t2lN25HtsHuladaJ+AKn/Ldbs6zh49XEMRdrgeQf3TNarBb9f1PR3+AlN368zuseisp5xjYivRjkPxdrSTvUH0/Useh4Bt+ru+AotSYCqI61JSCqFTWypWkLngGgIjHHiP41lugJB++eRT++ASn3cITUdGs8TTgZfTijs53MKXbFEyayc2VNzLZibBruX7GXNBbIwVpkH+iQvi4WiKVhY9mhNY99fCJ7FY33WYl+ZD0Exz5HjIS9DqKs1z7klIBqYdegKuV1LYPdLwOwi849+dWlGpSQVSHmmIQARTv2cOxyVNwZGefOt4I4Pf3YdNcyD/BV/6BPB4SiEMI+rfqz5zL5hDipY7mx1oIez6ChA2nhl1b8vRRcZaccuFTOuCgkoEfwe308Gl7mb4Pp74UpMPRbXpLKe+4Pk1Q6b4kW5G+D8q/lR5KrS6BSyfpQ8IVpZ6pIKpACBEEvAucDxQDqcA0KeXBs923qQYRgCMnh8S7x1Gybx8e559PzDtvYwwKgtR98OU/kcf+R5LRwLjoSNI0aB/Qnrn959IltIu7S3cPKeHId/DbSijJ05c5bHqrIzdJD6DQ8/UNe2Xh4+EPcf+AuCsgyA3TK2Uddg1y+BGKsqAoU28lFaQBErxDILyzfkqRC0eCQR3krNQfFUQVCCECgZ5Syo2u32cAI6WUA85236YcRKDvN0qZ+Ti5n36K5u9PmzfewPuS7vo3/E1zYc8HlFjzmBEZyY+eZoI9gpje/V5GnT+qZQ3xzjoCO9/SQwf0UMo7Dhl/6WEUGAMh5+ktjPI0I0RdrO+Pibq4cWzcpYSTf5wa5GDJgdzjkJOo71cKbKvvq7rsXgjr6O5qlWaq0QeREKI18DDQE+gGeAFxUsqjlazbBngRuAYQwEbgfill0jk8f0/gEyll7NnWbepBVCr7k09IfXIu0ukk/F8PEXLnnfoGa+9q2DwPmX2EZYGB/DfQH7PRg+HtR3Bv93sJ9Ax0d+n1q6QA9nwACRspG4ptyXV1w+XqXW/hnfTWTnmBMXrXW2x/fSqexspuheM79VBK3qm/rtxjYPTUW0fdb4OLb1MzNSh1rikE0QDgQ2AXYAAGUUkQCSG8gd+BEmAm+pZiHuANdJVSFtby+d8DMqWU951t3eYSRADF+/ZxbMJEHJmZ+A2+llauA2FJ2wdfP4Iz6Ud+NQnujQinSNPoGHwBYzqOYXDcYHxMzWwosJRwcJMeQiX5+jKH1dUNd0zvhgvrqB+rU9YyFHrX2/mDm+ZEpIWZ8NsK+OsbSP1Tn7XBJ1zfj3XZDGjdw90VKs1IUwgiTUrpdF2fACyh8iC6D3gB6Fi6P0cIEQckAP+WUr7gWrYRuLiKpxshpdxe7jGfAK4FBkopi85Wa3MKIgBHfj5J4ydg2bMHc7t2xLzzNqawMP3b/3fP4vz9AzJLspgaGcFfZhN+Bi86h3ZhVMfRXBlzJR4GD3e/hHOXcRB2LoOsQ/rvUuo79zMO6MOiA2P0A0TLd8MFxUKvCc1jBoOUPfDLUn3kXWYCIPTX1WW0/hq9mnkrWGkQjT6IyjtLEG0CPKWU/Sos/w5ASnlFDZ9rJjAMGCSlzK3OfZpbEIG+3+jk3HnkvP8+mq8vrV59Bd/evU911X2/EHv6ftb4+fBicBB5miDU5EfX8B6M7DiKfq36Nc3ZvC15sHuVfg6n8t1wqfFQUtoN1/n0UW4mb+h2C5w3qHnNfO2ww74v4Lf34OQefTi6h78+1PyyGfo+L0U5B80piE4Cn0spJ1dY/jowWkpZ7RkfXS2hodQghKB5BlGp3LXrSJk5E2m3E3b//YROnKDfkBoPm56k+Pj/kEVZLA4MYEWAH3YhaG0O4uLo3txw/mh6RfZqGgManE59VvI9H+pDs8HVDXdAbwkZPFzdcFGnd8O1uwIu/r/GvQ/oXBWkw87l8NeXeheto0QfzNBpuD6YwS/C3RUqTVRzCiIr8IKU8pEKy+cBj0gpq/W1XAhxIbAXOASUngzGXuWbJMQkYBJATExMj8TExOq/oCam5NAhku4ehz0tDd+rrqLVf19GMxr1lsKvK5Dx68jN3E+JNZ+FwUF84+uNCUGsRwjdW/djRMeb6RrW1d0vo2rpf+ndcNlH9d+l1PcBZSTo3XBBbSG4w+mj3YJioed4CDvfHRW7x/Fd8PObkLhDH6pu9ICILnrrqNMwNX+dUmMqiOpQc24RlXIUFXFs0mSKd+7EFBND2xXvYIqM1G/MSIBdb2NJ+omCzP2kOkt4KiSYPZ4eeKPR3jOMnm2uYFinWzgvqBHtPynO0aflObKNsm644hx91FhJnj4LQninM7vhuo7Ru+Eq2fBKKUkt0k/v4GX0wtvojanicO6mzG6F+DWw6229u64kXz8H0nnXwuUPQXCcuytUmpDmFESpwJq66JqrrZYQRKVSn3mWrLffRvP2ptXLL+Hbv79+g5T6TNC/rSQt/U8cGQf40wjPhgSRYjQSjIF2nmH0irmSYRfdQRu/Nu55AVLqk4Ge+E3foNqK9eV2q348UN5xVzfcBfqsCKeNhrtcH8pcrhsu25JNYl4iR3KPkJiXSFJ+ElaH9bSnNGpGvI3eeJu8y8Kp9LqX0eu05aW/l173Mno1zq7N/JP6YIb4LyDTdcx36PnQZ5q+v8xodm99SpPQnIJoM2CWUvavsHwr+mup0WCFGtY1DBjWoUOHiQkJCfX1NI1O3rffcuLhh5ElVkKmTCF8xr2nbrRb9dNK/LmG/TkJ+GQeYbOXmTeCAigWgihhItYzlD5tB3J913GEedfz94TibH0EXGaC3nLLOqRPcVNKSr2bKSNBPzdQUKw+7Fqr2A03DktQWxLzEjmad5TEvEQS8xLJLan2rsRaEUKcEVI+Jh98Tb74mf3KLstfb9BTdyT9DD+9ps/UUJiutx5jL4crH4PIixquDqVJak5BdD+wEDhfSnnYtSwWffj2I1LK5+u7vpbUIipVkphI0l13Y09Jwbt3b1q/uQiDp+epFYqyYPcqCg5v5fe8Q4TnHOdDP28+9vPFiKC1MBNjDqZf3DUM6T6ZAI862Nlvt0L2ET1UMhP0ACrKqHr94mxXN1y+PrVNWCfw8AXAISUnjEaOtr2Uo37BJOYnkVqYiqx4biHAIR0UWAvIt+aTZ82jwKbvYtSEhiY0DMKAQRjKrmtCw6AZqr5NGNC0ym+rTuvIIAz4mn3xNfnia/bF3+xfdr2y8PI0ep71Mf+WzaKf/mLnMv2UFI4SfXh7j3Fw6cSy91RRKmoSQSSEGOW6ejUwBZgGpAPpUsrvXOv4oB/QWsypA1rnAn7oB7QWVHzcutYSgwjAWVJC8vR7Kdy2DWNYGG2WLtFPRV5e5iHY9TaH035nf+5hIvPSWRzox3ZvL/ycEKN5EmXyZ0C7wVzTfSreNdlo5aWcaulkHtKnp6l4HqDypNTPkFpSAAUpkHdC3+Ee1okMnxASsXIUK4nSyjFvf2z+UfppE057CEmxvZg8a15Z8BTaCssCysPggZ/JD01oOKQDp3SedumQDpxO/XploXY2GnqIGYQBo2Y8dakZMIqzXJZbXyt31leTZioLrhCvEDqHdKZraFd8zTUMkNxk2PE67F+rd38aPKBVT7hqJsReVuPXqjR/TSWIqirku/JzwAkhYjh9ip9N6FP8HK3vGqHlBlGpjMVLSP/vywjNQMR//kPQmJvPXOnodpy/rWRbwWGyc47ibcnl9UB/DpnNRDihleZFuMmXgXFDuLLHNMwVT9hmLXQFzsFTl9a/+Y7hsOm3l+TjKMmHknyEtQDNFVROBEcCo9gcFM4h4aAA1wzaJm/9XD2ujbDVYT0tdPKt+Thcs20bhKGsZeFv9sfP7Fejg3mllJWHlfPMZWWXrhCzSzsOpwO7067/Xu6yOgEnEJUGmUkzEegRSIhXCJ1COtEtrBsXh11cs2mcjv4AP7ykz2FXkqcPZugyBi5/ELyDq/84SrPXJIKosWup+4gqU7hrF8n3TMeZk4Pf0CG0WrhQnxqoPLsV9q8l98/VrLWk4JWbTKHdwtJAP3I0jVipEaF5EW704drYwbT2iaIk+zAl2UexFmVQggMrkhLXj1U6sUkH3tZi/KyFBFqLCbYWE1piIcBhK3vaIk0jxeTBCbMnJ8wenDB7cNzDE4tmOFWbMODwi6TA5EmeLb8seEocp/YnlXZplYaOt9G7UQ4kcErnacFUet3h1APstN+d9lOhJu1YHVZsTv298zf7E+IZQqhXKBcEX8DF4RfTPbx79fbrWYv004n8b4l+LBboIxBHvAHRjXgov9KgVBDVoZbeIirlyMkhacIELHv/xBwbS8zyZZiios5csSgLfv+AA4e/ZY09nbjcVA5j5wN/XwxAe6eBYIMXhnIbeSElwXYb0dYSoq0W12UJ4bYSSuPEDpw0nx44J8ye5BiMZ8x4IKWkCEk+DvKMJvKFpMB2alpCD4NHWeCU7mMxlA+uCjyNnsT4xRAXEEcbvzZ4GDxw4kRKiVM6T/uRSH05FZa71i293SEdZzwGcEawOJyOsm6/iiFTfrndaS9rWdmctrKWV/lWlJSSAlsBmcWZZFoyy/Z3eRo8CfEKIcQzhI7BHbkk4hIuDruY1n5nOW9R9lG9dbR/nT6YwTsUbn0f2lz69/dTWgQVRHVIBdEpUkpS5z9F9sqVCC8vop97Dv+rr6p85cxD2Hct57u0XWyxZdI5P5OfjbDZx5sgh5M4h8TL6cDT6cTs1LvDHIBTCIo0jWJNo0gzUKxpWIRGiSaQcMYPrs1s+R8HsrQzrqyLrXzwmA1VDz8WQtDKtxVt/dsS6x9LXEAcEd4RjbJ1VF1O6cThdFDsKCYhO4H9WfuJz4wnrSitLJRyLDk4cWIQBoI8gwjxDOG8wPO4NOpSuoV3I84/rvL3QEp9yqQt8yH5F32apDErIbbfmesqLYoKojqkguhMeRs2cOLhR5DFxQSNHUvEIw9XvaFO3EH2r8v5tOAQB215tCvI4jsPA8dMRgwSpACJwCkETgQOIZBCIKDCj/74lS0XgBAaQjOAMGAwmPH1CsXPw/+sXWyBHoHEBsQS5x9HW/+2xPjH/G1QNScpBSnEZ8WzP2s/f2X+RVpxGlmWLDKLM7E69eOlSrvw2vq3pV90P7pHdKdDYIczW5Al+fDpJPjra/D0h9EroP2Ahn9RSqNRr0EkhOgAhKGfRuHAOT1YI6b2Ef092/HjJN49DltSEp5duxLz1lIMflWcEttuhf3r2Ld3FR/ZUimyWzBKSW4l3WqVEpp+WgaDWT+Y0mDWR22VXTcBZ38cs8FMrH+s3toJiCXWP7ZuhpY3AzaHjYM5B9mXtY/4jHgSchKq7MKL8oni8taX0zOiJxeEXICpdPShvQQ+m6wfTGz2g1HL4byBbnxVijvVeRAJIYzAY8A9QKhr8TtSynGu229z3TZJSrm3VlU3UqpFVDWnw8GJBx8k/5v1GIKCaL3oDby7dav6DkVZ2HavZNPhrzggLdiRejea0LAbPbAbPXEYPYgWXB4AACAASURBVHAYTdgNZuyaEYdmwKkZajwYWghBlE8Usf6xZaET5RPVpLvYGlKOJYd9WfvYl7WP3Wm7Sc5PJsuSRXZJNk55qgsvwjuCftH96NeqH93DuyOcdvjiXn2CWZM3jFwCFwx198tR3KBOg8gVQl+hH+9jRz+YtDPwdrkgigUOA3OklHNqXXkjpILo7LJWrSLt6QVIKQl74AFCx919ljscgfwUfeivTyh4Bp61ZVR+R33FIc1ll64d/EZhpJVfq+Zx7qRGQEpJUn4S+zL3sSdjD7vTdpNRnEGmJbNsyiN/sz89Inowt99cAky+8OWD8OsK/ViuGxbBhSPc/CqUhlbXQXQ/+snpNgJ3SilThBBOygWRa70D6AejNqu9lCqIqqc4Pp5jkybjyMjAd8AAWr3yXzRTM5oQVClTbC/mQPYB4jPi2ZGyg8O5h0kvSqfIXkQr31bM6z+PnuE9YP1j8L/FetfpsP9C10qOQVOarboOol1ADHCelDLHtayyIPoc6C6ljKl15Y2QCqLqcxQVcWzyFIp/+QVjdDQxy5bhEdvW3WUp9SytKI296Xt57ffXOJB9AC+jF3d2vpPJXSdh3Po0bP+vPpv5dS9A99vdXa7SQP4uiGpzUpGOwM+lIfQ38tEHMTQLQohhQojFubn1O/Flc2Lw9ib23RWETJmM/eRJjtwwgtx169xdllLPwr3DuartVawYsoKRHUbikA7e3PMmEzdMIq3vFBjwiD7Me90/4Zdl7i5XaQRqE0QScFZjvWjAUovHb5SklGullJMCAtSoqpoKv/9+2ixZjDCZOPGvf5My83HUYQPNn4/Jhzn95jC/33xCvULZmbqTMevG8F1cLxg4BxDw9b/hp0XuLlVxs9oE0RGgmxCiyvsKIbyArsC+2hamNC++/frR7ptvMHdoT84nn3BkxA3YsrLcXZbSAAbHDWbl0JX0iuhFZnEm/9z6T57ScrANXgCaAb79D2x/2d1lKm5UmyD6AmgNPPg36/wbCAI+r01RSvNkCgmh3RdfEDDyRkoSEjh87WAKdvzk7rKUBhDlG8XiQYuZ2m0qJs3E+/vf57aT33Di2jn6sV+b5sB3z7q7TMVNahNELwAngQVCiFVCiJGu5aFCiCFCiGXALCAJeL2O6lSaCSEE0U89RdSCp5FWK8cmTCDtv6+4uyylARg1I1Mvnsqb17xJW/+27Mvax+iD7/J1/ylg9IStC2DTXHeXqbhBbQ9o7YLe2omFM44tFMAx4LrmdjArqFFzdankyBGSxo3HnpKCV/eLaf366xiDgtxdltIACqwFzN4xm02Jm5BIhgV3Y84fm9CshdB3OgxSgdTc1MsUP0IIT+BuYAjQDjCgB9DXwGIpZeHf3L3JUVP81A+n1crxGfdRsHUrmo8PkbNnEzDseneXpTSQzxI+4/ldz5Nbkkt7z3DePnKAQEs+XDoJhqquuuZETXpah1SLqH7kfPoZJ+fNQxYV4TtgANEvPI/B29vdZSkN4Hj+cR767iH2Zu7FR/PgldQMehXkQI87YZgaxNBc1PVxRIpS5wJH3kiHb9fj1a0bBVu3cujqgRT88IO7y1IaQCu/VqwYuoK7Ot+FTUjGh/nzakgozl1vw2dT3V2e0gBqHERCCLMQItzVNVd+ua8QYp4QYq0Q4hUhRJu6K1NpCYyhocR++AHh/3oIZ2EhxyZO4sQjj+C0291dmlLPTJqJB3s9yOtXv064dzhv+nsxLjqKzD8+gE/G6wfAKs1WbVpEjwMpQPfSBa5jir4HHgWuQ595e4cQIqQuilRalpDx44lbtxZzu3bkrvmcQwOvoXhvsxv3olSid1RvPh3+KZe3vpzfPMzc0DqaHQfXIj+8XYVRM1abILoaOC6l3FFu2Y3AxcBeYALwGfrMClPOuUKlRfKIiaHdurUEjx+PPSODo2NuIfWZZ9SMDC2Av4c/r139Gg9f+gg2Dx8mR4bzUuo2Ct8bqcKomapNEMUCf1VYNgJ9GPftUsplwGj0VtON51Sd0qIJIYj410PEfvgBxvBwspa/zeHrrqckMdHdpSkN4P86/R+rrvuA9v5xLAsMYGLRn+xZfhU4HGe/s9Kk1CaIgoHUCssuAxKllH8ASCmdwM/os3Q3C2rSU/fxuvBC2m/cQMBNI7EePcqRYcPJWPqWu8tSGkC7wHZ8NOJTRrUbziGzmUkijeXLeuO0W91dmlKHahNENqBs5k8hRDj6cUQVhzgVAb61L61xUZOeupdmMBA9fz4xS5eg+fqSvnAhR24egy0jw92lKfXMZDDxxD/m88wVC4l0wgvmEqYv78HJnCPuLk2pI7UJogNAv3Kj5m5C75arGERRQNo51KYoZ/C57DI6bNqI79VXYfnjDw4NupacTz9zd1lKAxjQbjDLbt7MDVbBDybJXauvZ9Nfq91dllIHahNEHwOBwPdCiBeAZwArsKZ0BSGEAbgEOFgXRSpKeZqXF21ee43o559HaBopjz1G0vgJOAoK3F2aUs+C/SJ48q6dPOEMolhIHv7xCZ7e8pAaxNLE1SaIXgS2AD2B+wEv4CEpZfnWzyD07rvvz7lCRalCwNAhtN+4Ae9evSjcvp2DVw8kf8tWd5el1DNhMnPT3d/zht8lXFhSwqqk9dy2eigZRaqbtqmqcRBJKUuAgcAVwM1ARynlaxVWswD/BN495woV5W8YAwNp++4KImY9jiwpIXnaNJIffAinzebu0pT6JASdR7/Hax3v4u7cPP4oTGbU6iFsS97m7sqUWlBzzdWQmmuu8bKlpHBs2jRK9u3HGB5Oq5dfwrt797PfUWnS5O4P2bzxAZ4IDqDIYOTWznfwUM+HEEK4uzSlHDXXnNIimKKiaPfZZ4ROn449O5vE224nZe48tf+gmRMXj+HqUR+yKj2fThYLK+JXcNtXt5FjyXF3aUo1nctpINqgd89FA55VrCallM3qxCKqRdQ0lCQkcOye6diSkjC1aUPUgqfx6dHD3WUp9Sn9AJZ3rudVcwnvBPgT6BHIwisW0juqt7srU6jj00AIIYzAq+hT+ZS2fSu2gaVrmZRSGmpWbuOkzkfU9EgpSZ03j+wPPgSnE+/evYl66inM0VHuLk2pLwXpyOVD2FyczKPhodg0E3dfdDczLpnh7spavLoOonnAY4Ad+ApIAKocNyulnFOjJ2jkVIuo6bEcPEjKY//BsmcPwmQi4MYbiHj0UTQvL3eXptQHmwXevYHkE//jnshIDpsMXBJ+Ca9d/Rq+5mZzjH2TU9dBlIg+zU8/KeWeOqivSVFB1HTlbdxE6vz52FNS0Pz8CL3nHoLvHKt2ajdHUsIn4yn5czVzw8P53MeTYM9gXr7yZS4Ov9jd1bVIdT1YIRz4riWGkNK0+Q+8mg6bNxH2z/tBStIWLODQoEHkb1Mn4Gt2hIDRy/Dodx9z0zOYn5lHniWHcevH8ebvb7q7OqWC2gRRElBS14UoSkMQQhA6eTLnfbeVgBEjsKWcJHniRI7edhslSUnuLk+pa9c8iRiykOEFhXxyIhV/NF7d/SoTv51Ika3I3dUpLrUJog+AK4QQqrNVabI0Hx+in1lA+6++xKtHD4p//Y3DQ6/j+L8fxlFY6O7ylLp06XgYs5L2TsHXRw7TyxjETyk/cf1n17M3Q51wsTGoTRA9hX4+oi+FEOfXcT2K0qDMMTHErnyPNkuWYIqOJu+LLzh4+RWkL1qkjj9qTjpeC3d/jZdHAG8l7GGaMZpMSyZjvx7LO3++4+7qWrxaHUckhPABdgCdgEQgGXBWsqqUUl59ThU2MmqwQvOWtWIF6a++hjMvD2NkJJEz/4PfwIHuLkupK7nHYdlgyE1iV1wf7jXlkW8roH+r/rw44EU8jVUdEqmcq7oeNRcKbAC6cubxQxU1m+OISqkgav6cFgupCxaQ++mnSKsNzy5diH5qPh7nnefu0pS6YMmHt4fCyT3kRnZhXGQoB3IPEekdyWsDX+P8INXRUx/qOoiWAuPQu+cWoZ/q4e+OI/quRk/QyKkgajmsJ1M5+dhjFO7YAZqG3zUDiZozB4M6OWLT53DA+2Pg4AZkYFvmdh/M6qPrMRlM/KvXvxjTcYy7K2x26jqIUtC74TpLKVvcebNVELU8Rbt2kTLrCayHDiG8vAi+4w5C75uBZmhWjf2Wad39sPNt8A5hw6DHeDx+CYW2QkaeN5I5lzWrY/Hdrq6PI/IDfmyJIaS0TN49etD+y3VEzp2L5ulJ5uLFHLxiALnr1rm7NOVcXf8SXD0LirO55qvHWXPxw7T1a8unCZ8ydeNUnM7Kdn0rda02QbQPPYwUpUUJGj2KDtu+J/iuu3Dm53PioX9xeMQNFPz4o7tLU87FPx6AG98Eh53I1ZNYHTuGbmHd+OH4D9z65a0U24vdXWGzV5sgeg0Y0NKGbgshhgkhFufmqoZgS6YZjUQ88jAdNm/C98or9Vm+x43n4KBryVq1Sn2Dbqq6joaxa8DogcfaGazw7cFVba4iPiuem764iaziLHdX2KzVdvj2AmAs8DiwXkqZXNeFNVZqH5FSnuXgQdIXPk/hjz8irVa0gAACbxhB6L33YvBVx3w3ORkHYfkQKEyDy+5jvr8HH/z1ASGeIbwz5B3a+rd1d4VNVl0PVnDUYHUppTTW6AkauaqCyOl0kpGRQU5ODg5HTd4ipbHz9PSkdevWmEymKtexZ2eT/sor5K37EmdeHsLDA98rLif8gQcwx8Y2XLHKuSvMgCVXQU4idL+Dpe0u5pXfXsHL6MWigYvUpKm1VNdBVKO+ByllszoLbFVBlJSUhBCCiIgITCaTmtG5mZBSkpmZSX5+PnFxcWdd3+lwkP3uu2S/9x625OOgaXh26ULYjHvx7devASpW6kRJASwdCOn7oONQvuh1C7N3zEEIwTP/eIaBbdVBzjVVp6PmpJRaTX7OvfymobCwkFatWmE2m1UINSNCCEJCQrBYLNVaXzMYCLnrLjps3EjrRYvwvOgiLH/8wbHxEzh4zSCy3lup9iM1BR6+MHkbtO4Ff33F8G2LeHXAyxiEgYe+e4j397/v7gqblRYTFA1B09Tb2RzV9ouF34AriPvoQ+LWrcV34EDsqamkzptHQt/LOPn00zjy8+u4UqVOGU0w7lvocA0k/shlXz7KuwOX4mvy5emfn+alX19yd4XNhtpyKko982zXjjavvkKH7T8QdOdYALLfWUFCv/4k33svJUePurdApWqaBrd/Al1vhpN/0PHDO/lk0DLCvcN564+3eHTbo+6usFlQQaQoDcTo50fko49y3o4fiZg1C2NkJPkbNnJ46HUcGTOGgh+2u7tEpSojl0DvqZB1mMi3h7PmykW0D2jPusPrmLB+Ag6nGqB0LlQQKUoD0zSN4P+7lQ7frqfN8mV4deuG5Y+9HJswgYPXXEPmu++p/UiN0ZAFcNVMKEjF961r+LjvU/QI78HPJ3/m5nU3U2hT57GqLRVELciqVavo2bMnvr6+REVFMWTIEH74of5Ok71161Zat25dZ4/39ttv079//79dp6SkhHHjxuHv709kZCQvvPDC3z6ewWDA19e37Gfr1q11Vm91+PbtS+z7q2j/zTf4DR6MPS2dtPnzSejTl5Pzn8Kend2g9Shncfm/YOjzYMnF9NYglneZwbWx13Ig+wAjPx9JWlGauytsklQQtRAvvPAC999/P4899hipqakkJSUxbdo0Pv/8c3eXVqdmz55NQkICiYmJbNmyhWeffZZvvvmmyvX79u1LQUFB2c+AAQMarthyzDFtaP3Si5y/40eCJ4xHGAxkv/suCZf140C//iTeMZa0F16k6NdfVWvJ3XqNg5veAkcJ4u2hLGxzPWM7jyWlMIWbvriJQ9mH3F1h0yOlVD81+OnRo4esTHx8fKXLG4OcnBzp4+MjP/rooyrXsVgs8r777pNRUVEyKipK3nfffdJisUgppdyyZYts1aqVXLhwoQwLC5ORkZFy2bJlZff98ssvZadOnaSvr6+Mjo6Wzz33nCwoKJCenp5SCCF9fHykj4+PPH78uPz5559lnz59ZEBAgIyMjJT33HOPLCkpKXssQL7xxhuyQ4cOMiAgQE6bNk06nU4ZHx8vPTw8pKZp0sfHRwYEBFT6OqKiouT69evLfp85c6YcM2ZMpesuX75c9uvXr1rvYUP//zqdTpn10cfy6J13yQP9/yHjO3WW8R0vkPEdL5D7unaTCdcOlsfuu19mvv++tKamNmhtisvBLVLOjZDyyVAp934q39n7juz2Tjd56XuXyl9SfnF3dY0OsFNWsV1tVrMeNCZz1v5J/Im8en2OztH+PDHswrOut2PHDiwWCzfeeGOV68yfP5+ffvqJ3bt3I4RgxIgRzJs3j7lz5wJw8uRJcnNzOX78OBs2bGDUqFHccMMNBAUFMX78eD766CP+8Y9/kJ2dzZEjR/Dx8eHrr7/m9ttvJzn51AxQKSkpvPjii/Ts2ZPk5GSGDBnC66+/zv3331+2zrp16/jll1/Iy8ujR48eDBs2jMGDB7No0SKWLl1aZXdidnY2KSkpdOvWrWxZt27dWLNmTZWv+7fffiM0NJTg4GDuuOMOHn30UYxG938shBAEjR5F0OhRADjy8ijYto3CH3dg2bsXW3Iy+UePkv/NN6TOnoMWEIA5Lg6vrl3x/Ud/vPv0QfubmSCUOtB+ANz9NawYBp+MZ+x1zxPa/yke3/44kzZMYn6/+QxpN8TdVTYJqmuuBcjMzCQ0NPRvN7ArV65k1qxZhIeHExYWxhNPPMG7775bdrvJZGLWrFmYTCaGDh2Kr68vf/31V9lt8fHx5OXlERQUxCWXXFLl8/To0YM+ffpgNBqJjY1l8uTJfPfd6edOfOSRRwgMDCQmJoYrr7yS3bt3V+t1FhTo52cMKHfiuoCAAPKrOF7n8ssvZ+/evaSlpbF69Wref/99nnvuuWo9V0Mz+PsTcN11RM+fR7vP19Bx107arVtL2IMP4HPFFWg+Plj27iV7xQqOTZzEX90v4eDVAzk2dRqZy5ZjTUx090tonlp1h0lbwSsQ1v2Tocf38cbANzBpJh754RFW/LnC3RU2Ce7/6tdMVael0lBCQkLIyMjAbrdXGUYnTpygbdtTEzq2bduWEydOnPYY5e/r7e1dtuFfvXo18+bN45FHHqFr164sWLCAvn37Vvo8Bw4c4IEHHmDnzp0UFRVht9vp0aPHaetERkZW+jxn4+uaZDQvLw9PT8+y635+lZ+1pF27dmXXu3TpwqxZs3juued49NGmcWyIR4cOeHToABMnAuAsKaFwxw4Kf/iB4j1/YE1MpGDLFgq2bCHt2WfRfH0xxcTgddGF+PTrh2///mg+Pm5+Fc1ASAeYsh2WXgWbnuTSoum8N/Q9xq0fx3M7nyO1KJV/9fqXu6ts1FQQtQB9+/bFw8ODNWvWMGrUqErXiY6OJjExkQsv1AM0KSmJ6Ojoaj1+r169+Pzzz7HZbLz66qvcfPPNHDt2rNIZCaZOnUr37t15//338fPz46WXXuKTTz6p1vOcbYaDoKAgoqKi+P3337nmmmsA+P3338teU3UeX9ZiNvrGQvPwwG/AAPzKDbiwnThB/patFP3vf1j278d68CAl8fHkfPQxaBrC0xPNxweDvz+GwECMoaGYIiMwRkVhahODObYt5jZt0Mxm972wpsA/CqbugCVXw45XOa8ok9XDPuGOr8eyIn4FJwtPsvCKhWr6ryq02CASQnwIdAIcgA14VEq5yb1V1Y+AgACefPJJ7rnnHoxGI4MGDcJkMrFx48aykWW33nor8+bNo1evXgghePLJJ7n99tvP+thWq5WPP/6Y66+/noCAAPz9/cumOoqIiCAzM5Pc3Nyy7rL8/Hz8/f3x9fVl//79vPHGG4SFhVXrdURERJCcnIzVasVcxYZx7NixzJs3j549e5KamsqSJUtYvnx5pet+/fXXXHLJJURERLB//37mzp3L6NGjq1VLU2GKjib4tv8j+Lb/A/RJWYt/203hd1sp3vsnjowMHHl52FNTsR49CpXNHC8EwsPDFVh+GAKDMISEYAwPxxQdhbl1Gz2wYmLQvL3rrHZpt+N0OMBiwWmzIS0WpM2G02JBGI2YWrXCUIfPd868AmHqdlg2CH5/n/CiLD4b9Sl3rB/Lt4nfctc3d7F00FJMBrXvrqIWG0TAZCllDoAQojuwSQgRKqVslmNjH3zwQSIjI5k3bx633XYbfn5+9OjRg//85z8AzJw5k7y8PLp27QrA6NGjmTlzZrUe+91332X69Ok4HA46duzIypUrAbjgggu49dZbadeuHQ6Hg/j4eBYuXMikSZN49tln6d69O2PGjGHz5s3Vep6rrrqKCy+8kMjISDRNIyMj44x15syZw9SpU2nbti1eXl48/PDDDB48GNBbeZ07dyY+Pp6YmBg2bdrEXXfdRUFBAREREdx+++089thj1aqlqdIMBnx69sCnZ48zbpNOJ/a0NKxHj2JNOobtxAlsJ1NwpGdgz8rCkZuLPTMLa/JxsNkqfXzh4YHm7Y3m54cwGpEOBzidSIcdHE6k0wkOh37p+jntupT6dSlP/fwNYTYjvLz0kAxwtepCQjCGhWOMjMTcKhpTmxhMrVthaIhuSJMnTNwCK0ZAwnq83xvJh2PXMGXzvfx88mdGrx3Ne0Pfw9eszlVVXq1OjFfnRQjRGngY6Al0A7yAOCnl0UrWbQO8CFwDCGAjcL+UMukcnn8A8Clw1iCq6jQQ+/bto1OnTrUtQWnk1P/v6ezZOViPHsGalIQt+Ti2kyexp6fhyMzCkZODs6AA6XAgNE2fr81gKLsuDAb999JLTUMYjWXLhNEARhPCoIHBiDAZEUbTqeVGAzic2DMzcWRn48jNxVlQgLOoCGm16kFWidLQMvj4oAX4YwgMwhgSrIdWVCSm6GjMbdroLa1zPamhlPDRWNj3BYRfCBM28fBPs/nqyFdEeEfw3tD3iPSJPPvjNCN/dxqIxtIi6gDcDOwCtgGDKltJCOENbAZKgDsBCcwDtgghukopazTHhhDiRWAEEADc1FxbQ4pS14xBgRiDuuPdvbu7SzmNw2LRW3LJydiOJWNPPYk9LQ17egb2nGycuXk4CwuxZ2YirX9V3uISAmNUFH7XDCTk7rsxRdYiMISAMe/C2vth13J4ow/PTNxCpE8ky/cuZ9QXo3jr2rfoGNzx3F90M9BYWkRaaQgIISYAS6ikRSSEuA94AegopTzoWhYHJAD/llK+4Fq2EajqNIojpJSnzS4phBgMzAX6SSmtf1erahG1TOr/t/mx5xdgTzmB9dgx7MddrbrUNGypqVj+/BNZXAyahrldO/yvu47g22/DUMUIzL+1aS5sex78ImDCFlae2MJzvzyH2WDm5Stfpm905SNMm5s6PUNrfTtLEG0CPKWU/Sos/w5ASnnFOTzvQWCMlHLX362ngqhlUv+/LYvTbid3zRpyPv6Ekvh4pM0GJhOenTsReOONBIwcWbORhDteg29ngmcAjFvPxsJjPLztYSSSOZfNYVj7YfX3YhqJOj1Dq5tdCOytZPmfQOfqPogQwsvVkir9vS8QAhyuYv1JQoidQoid6enpNSxZUZSmRjMaCRo1irgPP+C87T8Qdt8MzDExWP7Yy8nZczjQuw+J48aTt2FD9Yb8970HbngDSvJhyVUMNAaxdNBSzAYzM7fPZNkfy+r/RTViTS2IgoHKpiPOAoJq8DhewCohxF4hxG5gIfo+okqnOpZSLpZS9pRS9qzuUGNFUZoHg78/oVOn0v7LdbT/dj1Bt9+OITCQoh9/5Pi9MzjQpy/J991HUSU9JafpdgvcsgqcNnh7KN0Lclg1dBWBHoG8+OuLPP3z0w3zghqhphZEdUJKmSWl7CulvEhKebGUsp+UsnpjiBVFabHMrVsTOfM/nLdlM20//gj/665DCEH++m9JvP0OEi6/ghOPP47l4MHKH+D8a2HsWhAGWDmadsm7WT1sNa18W7Fq/yr+ueWfTfqg6tpqakGUTeUtn6paSnVGCDFMCLE4Nze3Pp9GUZQmwrtLF1o9v5DzdvxI69dexad/PxwFBeR+/AlHhg3n4KBrSV24EGtq6ul3jOkNEzaA2RtWTyA0fi2fDv+UTsGd2Ji0kTu/uRObo/LjtJqrphZEf6LvJ6qoMxBfn08spVwrpZxUfkJNRVEUIQR+V19NzNKl+mngn5iFZ5cu2FNSyFr6FoeuHsjhETeQsXQpjtIJeCMuhMnbwDsY1v0T7x2v88F1H3BZ1GX8lvYbN629ifySyifrbY6aWhB9AfQRQpTNVimEiAX6uW5TamH27Nll0/kkJSXh6+uLo7KpXupZbGwsGzdubPDnVZS6YvDwIPjWW4n76EM6fP8dodPvwRwTQ0lCAukLnyeh/z84/tBD+mwSQW1h2k8Q0Bo2z0XbMIs3B73J8PbDOZJ7hBu+uIGUwhR3v6QG0WiCSAgxSggxCiide2SIa1n5IdlLgKPA50KIEUKI4cDnwDHgzQYtuAmqzqnCY2JiKCgowGAwuKlK3ezZszGZTPj5+eHn58f555/P9OnTSUk59cHcunUrmqbh6+uLn58fHTt2rHJeOUVpaMagIMKmT6f9V1/S7qsvCbr1VgwBAeSt+5Kjt9yK02oFn1B9stTQ82HHK/DZVOb3m8eELhNIL0pn1Bej2J+1390vpd41miACPnb9THH9/rrr9zmlK7hmTrgKOAC8C6wEjgBXSSmrd66AWmrq+4ia4qnCx4wZQ35+PllZWXz22WecPHmSHj16nBZG0dHRFBQUkJeXx4svvsjEiRPLzpOkKI2FR2wskU/Mov3WLfhdNxTLnj0cHj4CR04OePrpp5GIvgR+XwXv38p93WfwWO/HKLQVcsdXd7D9xPazP0kT1miCSEopqvgZUGG9JCnlTVJKfymln5TyhsrmpKuH+prsPqLc3FxmzZrFa6+9lzYn3AAAGthJREFUxsiRI/Hx8cFkMjFs2LAzTgR39OhRhBDY7XYABgwYwKOPPsqll16Kv78/I0aMICsr67R1Fy9eTHR0NFFRUSxcuLDssZxO5/+3d+dxUdf5A8dfbxAREBUQL0BErbw17UDN2+xQ8cjNzPDK+qWWtR1am4lptukmpeaVFmoamrRpurlumbZq5WamebXWpuJ95hWiIJ/fHzMQIOgMMHwH5v18POYxzPfzPd7zZpg3n+/14Y033qBOnTqEhITw4IMPZi0LtpulRkZGEhISwsSJE/ON38fHh4YNG7J06VJCQ0OZMmXKNfOICPfffz/BwcH8+OOPhcqXUq7i5eVF+JQpBA8aRNr+/fyvewxXDh+GMmVh6Fqo3R72roaE+3jopj5MaTcFg+HJtU/y6S+l9+iDu9xrrvRZ/SIc2+HabVRrDPe9ccPZHBkq/HoWLlzImjVriIqKYsCAAYwcOZJFixZlta9bt46ff/6ZX3/9lY4dO9KsWTM6d+7M9OnTWb58OV999RWhoaGMHDmSESNGkJiYyO7duxk2bBifffYZd955Jy+99FKOIcXz4u3tTY8ePVizZs01bRkZGaxatYpTp05Rt27dAr1PpYpL1RdHU6ZKFU68+Sb7evai5vwE/Bo2hNjlkDQYdn0C77aj09C1zOsyj2FfDOOVr1/hxKUTDG081Orwi5zb9IiU6zgyVPj1xMbG0qhRIwICApgwYQIfffRRjpMZ4uLiCAgIoHHjxgwePJjExEQAZs+ezcSJEwkPD8fX15dx48aRlJREeno6SUlJdOvWjbZt2+Lr68uECROyxjG6nho1auToVR05coRKlSrh5+dHr169iI+P51Y3uxGnUnkJGTKYGpMmkZGayoGH+3Nxw0bbzVL/NB9uGwLHd8KsljQLrMWSbksI8g1i6tapTPw2/70HJZX2iFzFgZ5KcXFkqPDriYiIyPo5MjKStLS0HGMB5W7fscPWEzxw4AC9evXKUWC8vb05fvw4R44cybFcQEAAISEhN4zl8OHDBAcHZ72uUaMGhw4d4vLly7z44ot8+eWXPPPMM06/R6WsULF7N7xDgjk04kkODhtG9fHjqdS7F3R7C/xD4N9vwsxoIh9fz8cxHxO7OpYl/13CiZQTvN3h7VIz4qv2iBxUkk9WyD5UeEEcPHgw6+fk5GR8fHyoXLlyvu2ZQ4xHRESwevVqzp49m/VITU0lLCyM6tWr51guJSWF06dPXzeOjIwMVq5cSZs2ba5p8/X1ZdKkSezYsaPA71MpK5Rv1YrIxA/x8vfn6Msvc3L2bFtDxzFwz+tw8TjMbElIylmSYpJoENyALw9+SezqWK5cve5gASWGFiIHleSTFbIPFb58+XJSUlJIS0tj9erVjBo16obLL1q0iN27d5OSksLYsWPp06dPjtO7J0yYQEpKCrt27SIhIYG+ffsC8MQTT/Dyyy9z4MABAE6ePJl1ll6fPn1YtWoVGzdu5MqVK4wdO5aMfAY0S09PZ8+ePfTr149jx47x7LPP5jlf2bJlee655xg/frxT+VHKan716lF7+SeUCQ3l1NtTOfqq/WThlsOhx0y4fB7mtMP/xM8kdk3krrC72H5yOw98+gDnL5+3NvgioIXIQzz33HPEx8fz2muvERoaSkREBO+88w49e/a84bKxsbEMGjSIatWqkZqayrRp03K0t2vXjrp169KpUyeef/55unSxjWv49NNPExMTQ5cuXQgMDCQ6OprNmzcD0LBhQ2bMmMHDDz9M9erVCQoKIjw8PMd6ly5dSvny5alYsSIxMTGEhITw/fffZ/W48jJkyBCSk5NZuXKlsylSylI+NWpQe9VKykZFcTZxCQeHj7Ddd65ZP+i7CK5egffvwevAJmZ1nkWvur3Yf34/PVf05MjFI1aHXyhuNx6Ru/O08Yjat2/PI488wtCh156ps3//fqKiokhLSyvwiRAlRWn9/Sr3k3H5MskDBnJp+3bKNW1K5AcLbWMfHfgGFvW23b27TwLU786MH2Yw58c5BJYNZO7dc2lQ2eHRcIpdaRqPyDIl+RiRUqrk8PL1JXJJIuU7diR1+3b29ehB+rlzENkSHv0cyvjBRwNh6weMuHUEY6LHkJKWwsB/DmTT4ZJ54asWIgeV5GNESqmSRUSImDmDSn0f5Mq+/ezr1p0rhw5DtUb2m6UGwadPwaZpPHjLg8S3j8+68HXNvmuvs3N3WojUda1fvz7P3XJgu0mpMabU75ZTyirVX32Vyk8+SfqpU+zr3ZtLu3dDcC3b/ekqhMHnr8AXr9KhZgcS7knAt4wvozeM5h//+4fVoTtFC5FSSrmx0CdHUG3cODJSUuwXvm6A8lVg+DcQchNsjIdPR9I4tDGL71+Mfxl/Xt70Mit+cd/7SOamhUgppdxcUN8HCZs2FYzh4LDhnP3kEyhXAYZtgurNYOsCWPoIdSrWZnHXxQT4BDD267Ek7U2yOnSHaCFSSqkSoELHjtRcMB+vcuU4+vIYTs2eDWV84bEvIaot7FkJC2KICqxJYtdEAssGMuHbCSz5aYnVod+QFiIH6VlzSimr+TdrRq1lH+EdFMTJt6dydPx48PKGAZ9Cgx6w/9+QcB81A2qwtOtSKpatyOubX2fxnsVWh35dWogcpGfNKaXcgW9UFLU/XYFPeDhnP0zk4IgRZBgDDy6EZv3h4GZ4/x7C/KuytNtSgsoFMek/k5i/a77VoedLC5FSSpUwZUJCqL1iOb7163Nx7ZckP9yfjMuXoedMaD4QDm+BeZ2p7leZj7p9RIhfCPFb4pn34zyrQ8+TFiIPsnHjRlq1akXFihUJDg6mdevWfPfdd8yfPx8RYfLkyTnmDw8PZ/369YBt6O5HHnkkq01ECAgIoHz58oSFhfHss89mDQ2R33aUUkXHKyCAWknL8G/ZkkvbtrGvZy/Sz5+HmGlw+1A4+gPM7UhV3yA+6v4Rof6hTPthGrO2zbI69GtoIfIQ58+fp1u3bjz11FOcOXOGw4cPExcXh6+vLwDBwcFMnjyZCxcuOLzO7du3c/HiRdauXcuHH37I3Llzb7gdpVTR8fL2JjLhfSp068aVffvYF9OD9DNnoOsUuPMJOPYjzG1PqE8FkronUdW/KjO3z2T6D9OtDj0HLUQeYu/evQD069cPb29v/Pz86NKlC02aNAGgfv36tGzZkvj4eKfXXa9ePdq0acPOnTtvuB2lVNELe/NvBA8ZTPrx4xx4uD8ZaWlw3yRo+SQc3wVz2hHkXY6kmCRqBNTg3R/f5a3v37I67Cx6SbyLTPrPJH4685NLt1EvuB6j7xjt0Lw333wz3t7eDBw4kIceeojo6GiCgoJyzDNhwgQ6dOjAU089lWPwuRvZvXs3GzZsYOLEiQ5tRylV9KqOGgVXr3JmwUIODn2MmvMTkHsmgrcPbHwL5rSl4uNfsaz7Mvqu6sv7O98nPSOdF25/werQtUfkqJJ++naFChXYuHEjIsJjjz1GaGgoMTExHD9+PGueZs2acffddzNp0iSH1tm8eXOCgoLo3r07Q4cOZfDgwQ5tRynlGlVfeonyHdqTsnkzx+LG2SZ2HgdtX4BTe2FOGyqIN8tillEzsCYLdy/kr5v/amHEdsYYfTjxaNGihcnL7t2785zurvbs2WNatGhhHnroIZOQkGBat26dNT0wMNAcO3bMhIWFmXXr1hljjImLizP9+/fPWh4wP//8s1PbKclK2u9Xea6Mq1fN/2J6mN231DOnEhL+aFj3V2PiKhgztZkxqRfM71d+N93/3t00mt/IvPr1qy6PC9hi8vle1R6Rh6pXrx6DBg1i586d10zv3bs3EydOdOl2lFKuIV5eRC5eRJkqoZyY/DcurFtna2j/InQcC2d+hVmt8E9PY2m3pdSpVIdle5cxdtNYy2LWQuQhfvrpJ6ZMmcKhQ4cAOHjwIImJiURHR18zb1xcHAkJCZw9e9al21FKuYZ3+fLUXLQYLz8/Dj/zZ1LtJxHR9jm4ezycTYbZd+GXnsqSrku4qdJNfPLLJ7y04SVL4tVC5CECAwPZvHkzd955JwEBAURHR9OoUSOmTJlyzbxRUVHExsby+++/u3Q7SinX8a0ZQfismZirV0keMJD006dtDa2fhnsmwrmDMOsuyl25xJJuS6gfXJ9Vv67iha+K/+QFHSrcSZ42VLiy0d+vKqnOJn3M0VdewScigtqrVtqGHQfY/C78cxSUrw7/t4F0v0rEro5l5+mddKrZibfav4WIFFkcOlS4Ukp5qEp9HiDk0SGkJSeTPORRsjofdz4O98fDxWMwuzVlLv3GovsX0SS0CWuT1zLyy5EUV0dFC5FSSpVyVZ5/nvKdOnFpyxaOjnnlj4bbh0C3t+H3EzC7Nd6/n2ThvQtpXqU56w+tZ/gXw4ulGGkhUkopDxA+fRq+9etz7uOPOf3ee380tBgIMe/A7ydh9l14XzjG/Hvnc0e1O9h4ZCOPf/64y4uRFiKllPIA4uVF5KIPKFO1KiemxHNh7do/Gm/tDz1nQ8oZmNMGOXeYeV3m0ap6K749+i1D1gwhw2S4LDYtRA4q6XdWUEop74AAIhcvsp3W/exzpO7Z80dj077Qey5cOmsrRmcPMvvu2bQJa8OW41sY9M9BLitGWogcZHRgPKVUKVA2PJyIObMxGRkcGDSYtFOn/mhs/AD0SYDL5+HdtsiZ/czoNIMOER344cQPjPt6nEti0kKklFIexv+226g+/lUyzp/nQL9+tkH1MjXsAX9aAJcvwNx2yJlfmdZxGl2junJbtTzPvi40LURKKeWBKvXqRchjj5F28BDJg4fkPCGhfjfouxiupMC77eHUXt5o+wYxdWJcEosWIqWU8lBVnv0zgV3u5tLWrRx96S85G2+5F/otgfRUmNsRTuzJeyVFQAuRh6hVqxZly5blVPb9wcCtt96KiLB//34GDRrEmDFj8lw++9DgmY/cQ4srpUqesKlTKdewIeeWL+fUnHdzNt7UGfovg/TLMK8znP6fS2LQQuRBoqKiSExMzHq9Y8cOUlJSHF4+c2jwzMeoUaNcEaZSqhiJCDUXLqBMtWqcnDqV82v+lXOG2u2h/98hIx1WPu2SGLQQeZDY2FgWLlyY9XrBggUMGDDAwoiUUu4g67Ruf3+OvPACl3btyjlD7TYwaBX0muOS7etQ4S5y7PXXubzHtUOF+9avR7W//OXGM9pFR0fzwQcfsGfPHm6++WaWLFnCpk2b8t0dp5TyHGXDwoiYN5fk2AEkDx5C7X+swic09I8Zwm932ba1R+RhMntFn3/+OfXr1ycsLMzhZZs3b06lSpWyHmvWrHFhpEqp4ubfrBnVXnuNjAsXrj2t24W0R+QizvRUilNsbCxt27Zl3759Tu+W27p1K3Xr1nVRZEopd1CpRwxXDuzn9MxZHBg4iFqJHxbpcBB50R6Rh4mMjCQqKorPPvuM3r17Wx2OUsoNVRk5ksB77yV12zaOjBrt8u1pj8gDvffee/z2228EBASQnp6eo+3q1aukpqZmvfby8qJs5kBaSimPEfZWPPsPHeL8ypWUjYoidPgwl21Le0QOKk03Pa1Tpw633Zb3rTreeOMN/Pz8sh4dO3bMamvatGmO64ieeeaZ4gpZKVXMRISaHyykTPXqnJo+nXOrV7tuWzpUuHN0qHDPpL9f5anSjh3j1+4xmMuXqf3pCsrWqlWg9ehQ4UoppQrEp1o1IubNBeDo2DiXbEOPESmllLou/6ZNCXv7Lco1aOCS9WshUkopdUOB2Y4XFzXdNVeE9Hhb6aS/V6VcSwtREfHx8eHSpUtWh6FcIC0tjTJldOeBUq6ihaiIVKlShcOHD5OSkqL/QZciGRkZHD9+HB0iXinX0X/zikiFChUAOHLkCGlpaRZHo4pSQEAAlStXtjoMpUotLURFqEKFClkFSSmllGN015xSSilLaSFSSillKS1ESimlLKWFSCmllKX0pqdOEpGTwIFskyoC55x4XRk45YLQcm+nKJe70Tz5tec13V3ylde2imoZzZfzy1xvPs2Xc/MVJl+5pxVlviKNMaF5thhj9FGIB/Cuk6+3FEccRbncjebJrz2v6e6Sr4LmTPPlmmWuN5/mq/jylXtaceVLd80V3konXxdXHEW53I3mya89r+nukq+Cbkvz5Zplrjef5su5+QqTr9zTiiVfumuumInIFpPPmBzqWpov52i+nKP5co6r8qU9ouL3rtUBlDCaL+dovpyj+XKOS/KlPSKllFKW0h6RUkopS2khUkopZSktRG5MRIJEZJWI7BWR7SLyLxGpa3Vc7kxEXrHnK0NEelodjzsRkToistGenx9ERA/S34B+nhxXmO8rLUTuzQBvG2NuNsY0BVYB8yyOyd19DtwL/NvqQNzQbGCBMeZmYBSwWETE4pjcnX6eHFfg7ystRE4QkXARmS4i34hIiogYEamVz7wRIpIkIudE5LyI/F1EajqzPWPMWWPMF9kmfQ3kuT13VNz5AjDGfGuM+bWwsbuDosyfiIQC0cB8AGPM54AALVz+RopRUX/mStPnKS9Fma/CfF9pIXJOXeBB4DdgQ34ziYg/8CVQDxgIxAI3AetEJKAQ238GWFGI5Yub1fkq6YoyfzWBo8aY7KM27rdPL030M+ccV+bL8e8rV93eojQ+AK9sPw/F1hWtlcd8TwNXgbrZpkUB6cCz2aZ9ge2+TXk9WudaZxy2/zD8rc5DCcnXeqCn1Tlwl/xh6/nszbXcv4DeVr9Pd81Zafs8FXO+nPq+0h6RE4wxGQ7OGgN8a4z5Jduy+4BNQI9s0zobYyrn89iUOZ+IjAHuB+4zxqQUzbtxPavyVVoUcf6Sgeoi4pNtuVr26aVGUX/mSjtX5Ksg31daiFyjIbAzj+m7gAbOrEhE4oDuQBdjTEHurl0SFFm+PNQN82eMOQn8BxgEICJ3YztG9H3xhOh29DPnHIfyVdDvKy1ErhGMbZ9rbmeAIEdXIiINgXFACPCViGwTkS1FEqF7KZJ8AYjIOBE5BLQE5onIIREJL4IY3Zmj+XsCGCwie4G/Af2NfT+KB3IoZx76ecrLDfNVmO+rMkUUpHIBY8wubP+1KgcZY8Zh+2NQuRhjfgZaWR1HSaKfJ8cV5vtKe0Su8Rt5/yef338Vnk7zVTiaP+dpzpzj0nxpIXKNXdj2qebWANhdzLGUBJqvwtH8OU9z5hyX5ksLkWt8CkSLSO3MCfaLxFrb21ROmq/C0fw5T3PmHJfmS4eBcJKI9LH/2Anbwd/hwEngpDHmK/s8AcB24BIwBtu5+ROAQKCJMeZiccdtFc1X4Wj+nKc5c45b5MvqC6pK2sP+C8jrsT7XfDWBj4HzwAVgOXlcKFbaH5ovzZ/mzL0f7pAv7REppZSylB4jUkopZSktREoppSylhUgppZSltBAppZSylBYipZRSltJCpJRSylJaiJRSSllKC5FSSilLaSFSSillKS1ESllMRKaKiBGRdlbHopQV9BY/SllMRJKBckA1Y0yG1fEoVdy0R6SUhUTkdiACWKFFSHkqLURKWau3/fkTS6NQykJaiJRygP0YjrH/3FdEvhGRiyJyQUTWishdBVx1L2y31f/CwTiq2mPZm0fbE5lxikjdXG2N7dM3FzBOpVxGC5FSThCR8cCHwBXgH8AhoCOwVkRaOrmuBsAtwGfGmCsOLvab/Tkw17q8geezTQrOtdyf7c+TnYlRqeKghUgp54wA7jDGtDPG9AUaAnOBssB4J9eVuVvu744uYC9YKeQqRMADQB3ga/vrrEIkIlWAh4Ff0F2Ayg1pIVLKOXHGmO8zX9hPMHjF/rKNiPg4sa7eQCqw2skYzgABIpL973c0cAKYZn+dvUc0HPAF4vWECOWOtBAp5ZxVuScYY45j22XmC4Q4shIRqQXcCnxujLnoZAw5ds+JyN1Ac2xF6Ji9Ldje5gsMA04C853cjlLFQguRUs5Jzmf6eftzOQfX4/RuuWxyHycaDVwEZgLn7NOC7M/9gSrAO8aYSwXYllIup4VIKScU4a6t3kA6sLIAy56xPweKyG1AJ2CuMeY34Ky9LXPX3DPYjinNyL0SETkqImNFZIyIHLSfAThXRLxFpJWIfCUiv4vIdhFpkmvZB0TkCxE5LiKXRGSniPwpW3sL+1l6A7NNqygiO0Rkg4g4WrCVB9BCpFQxE5GqQEvg38aY0wVYRWaPqAK23lAa8JZ9WmaPKFhEOgONgfdzb8d+AkM14FGghv15OjDU/jwLmA30w9a7is8VQxMgCXgE6A6sBz4UkVsA7MfRlgMv2wubj31+H6CHMSa1AO9blVJlrA5AKQ/UE9s/gQXZLQd/FKLm2HpWi40xB+3TzgMGW4/oz8BVri0iAE3tzx8YY8bYf/6XiDyJrbA0sfewsPe6Hs2+sDEmLvNn+6nj6+3z3AH8194UB2zDVsw6Yite0caYMyiVjRYipYpfb2zFYnkBl8/8Iv8LIGS7NsgYc1VELmIrUtWBZcaYfXmsowlwGZiUOUFEymA74WJmZhGyC+SPXX6ZJ0D8HzAYqI2tZ5YpJVssP4rIMmCOfVL7fGJRHk4LkVLFSEQqAR2A/xhjDhdwNZlFIhzbxbA7c7Wfs7dB/hewNgW+M8ZcyDatAbbrodbmmrcJsBNARARYga3QvQ1sBU4D7e3b+inXsr8A/sBkY8x3N3pjyjPpMSKlild3bMdJCnNhafbeSl6FJvM40brs1zzl0hTbbrPc064CO/KYvt3+c0vgHqC/MeZ1Y8w/7QWmDrZrorIKkYj0w3YM6zsgVkT8rvuulMfSQqSUA4wxYoyR67TXss+z/warKsxp25nbSsyMxxjzVR7tjextHfNa3n7iQD2uLUTNgP9mP81bRMKwXRuVOW+E/Tl7wWkADAJ2GmOu2qe1w3bd0migj30dw518q8pDaCFSqnh9A7xgjPnZwhjqY9sFl1ch+iHXtMyTGjJ7RFuxnwAhIp3sJzeswNYb2gYgIvWxHf+aZ4yZYoxJBhKA0SISUNRvRpV8WoiUKkbGmMnGmDctDqMptmuYch9bym933anM41n2AvoYtl10K4Cu2O5zVw7YJiLVsN2yaCMwMtt6XgcqAk8X6TtRpYKO0KqUUspS2iNSSillKS1ESimlLKWFSCmllKW0ECmllLKUFiKllFKW0kKklFLKUlqIlFJKWUoLkVJKKUv9P26ojxaFaHH2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "class FlassPlot:\n",
    "    @staticmethod\n",
    "    def pic(x, y, label):\n",
    "        import matplotlib.pyplot as plt\n",
    "        \n",
    "        plt.loglog(x, y, label=label)\n",
    "        plt.legend()\n",
    "        \n",
    "    @staticmethod\n",
    "    def forpaper():\n",
    "        import matplotlib.pyplot as plt\n",
    "\n",
    "        LEGEND_SIZE = 12\n",
    "        SMALL_SIZE = 16\n",
    "        MEDIUM_SIZE = 22\n",
    "        BIGGER_SIZE = 24\n",
    "\n",
    "        plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "        plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title\n",
    "        plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "        plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "        plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "        plt.rc('legend', fontsize=LEGEND_SIZE)    # legend fontsize\n",
    "        plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "        \n",
    "    @staticmethod\n",
    "    def axeslabel(xlabel, ylabel):\n",
    "        import matplotlib.pyplot as plt\n",
    "        \n",
    "        plt.xlabel(xlabel)\n",
    "        plt.ylabel(ylabel)\n",
    "        \n",
    "    @staticmethod\n",
    "    def title(title):\n",
    "        import matplotlib.pyplot as plt\n",
    "        \n",
    "        plt.title(title)\n",
    "        \n",
    "    @staticmethod\n",
    "    def savefig(filename):\n",
    "        import matplotlib.pyplot as plt\n",
    "        \n",
    "        plt.savefig(filename, bbox_inches='tight')\n",
    "        \n",
    "    @staticmethod\n",
    "    def plt():\n",
    "        import matplotlib.pyplot as plt\n",
    "        \n",
    "        return plt\n",
    "  \n",
    "import pickle\n",
    "allres = pickle.load(open( \"epsilongreedy_estimate_allres.p\", \"rb\" ) )\n",
    "\n",
    "renameit = { 'MLE': 'EL' }\n",
    "FlassPlot.forpaper()\n",
    "for name, res in allres:\n",
    "    x = [ x[0] for x in res ]\n",
    "    y = [ x[1]['mse'] for x in res ]\n",
    "    ylo = [ x[1]['mse'] - 1.96 * x[1]['msestd'] for x in res ]\n",
    "    yhi = [ x[1]['mse'] + 1.96 * x[1]['msestd'] for x in res ]\n",
    "    FlassPlot.plt().loglog([ x[0] for x in res ], [ x[1]['mse'] for x in res ], label=renameit.get(name, name))\n",
    "    FlassPlot.plt().fill_between(x, ylo, yhi, alpha=0.7)\n",
    "FlassPlot.plt().legend()\n",
    "\n",
    "FlassPlot.axeslabel('n / $w_{max}$', 'mse')\n",
    "FlassPlot.plt().savefig(\"epsilongreedy_mse.pdf\", bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}