et02_tu.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Matching models with transferable utility</center>\n",
    "### <center>Alfred Galichon (NYU & Sciences Po)</center>\n",
    "## <center>'math+econ+code' masterclass series</center>\n",
    "#### <center>With python code examples</center>\n",
    "© 2018–2023 by Alfred Galichon. Past and present support from NSF grant DMS-1716489, ERC grant CoG-866274 are acknowledged, as well as inputs from contributors listed [here](http://www.math-econ-code.org/team).\n",
    "\n",
    "**If you reuse material from this masterclass, please cite as:**<br>\n",
    "Alfred Galichon, 'math+econ+code' masterclass series. https://www.math-econ-code.org/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References\n",
    "\n",
    "## Textbooks\n",
    "* [OTME] Alfred Galichon (2016). *Optimal transport methods in economics*. Princeton.\n",
    "* Pierre-André Chiappori (2017). *Matching with Transfers: The Economics of Love and Marriage*. Princeton.\n",
    "\n",
    "## Papers\n",
    "\n",
    "* Llyod Shapley and Martin Shubik (1971). \"The assignment game I: The core.\" *International Journal of Game Theory*  1, pp. 111–130.\n",
    "* Becker (1973). \"A theory of marriage part 1.\" *Journal of Political Economy*.\n",
    "* Eugene Choo and Aloysius Siow (2006). \"Who marries whom and why.\" *Journal of Political Economy* 114, no. 1, pp. 175-201.\n",
    "* Alfred Galichon and Bernard Salanié (2011-2021). \"Cupid's invisible hands: Social Surplus and Identification in Matching Models\". *Review of Economic Studies*, forthcoming.\n",
    "* Arnaud Dupuy, Alfred Galichon, Sonia Jaffe, and Scott Kominers (2020). \"Taxation in matching markets\". *International Economic Review*, forthcoming. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading the libraries\n",
    "\n",
    "First, let's load the libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.spatial.distance import cdist\n",
    "from time import time\n",
    "\n",
    "# !pip install gurobipy\n",
    "import gurobipy as grb\n",
    "import scipy.sparse as spr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Generating worker and firm data\n",
    "\n",
    "We will consider an academic job market, where candidate professors match with colleges--or rather, academic positions, as each position can only be filled by one applicant. \n",
    "\n",
    "A professor $x$ has characteristics $\\xi(x) \\in \\mathbb{R}^d$ where $\\xi$ includes measures of skills, such as general skills<br>\n",
    "* $\\xi_0$=teaching ability,<br>\n",
    "* $\\xi_1$=research ability, and<br>\n",
    "* $\\xi_2$=collegiality<br>\n",
    "as well as specific skills, like<br>\n",
    "* $\\xi_3,\\xi_4,\\xi_5$=specialization in micro, macro, and econometrics<br>\n",
    "and in addition, let's assume $\\xi$ also encodes geographic coordinates of the applicant's residential address<br>\n",
    "* $\\xi_6,\\xi_7$=residential address, (horizontal,vertical) geographic coordinate<br>\n",
    "\n",
    "\n",
    "College offer positions $y$ with characteristics $\\zeta(y) \\in \\mathbb{R}^d$, where $\\zeta_k$=how much one unit of characteristics $k$ is valued (in dollar terms) for this position, for $k\\in\\{0,...,7\\}$ and $\\zeta_6,\\zeta_7$ = college location's geographic coordinates."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preferences\n",
    "\n",
    "We assume that the amenity that worker $x$ associates with working for position $y$ is equal to minus the Euclidian distance between her location and the position's location, that is\n",
    "$$\\alpha_{xy} = - \\sqrt{(\\xi_{x6} - \\zeta_{y6})^2+(\\xi_{x7} - \\zeta_{y7})^2}.$$\n",
    "\n",
    "We assume that the output $\\gamma_{xy}$ that worker $x$ produces if she takes position $y$ is a CES function where specific/general skills are more or less substituable \n",
    "$$\\gamma_{xy}=\\left( \\left(\\sum_{k\\in\\{0,...,2\\} }(\\xi_{xk} \\zeta_{yk})^{r_g} \\right)^{r/{r_g}}+\n",
    "\\left(\\sum_{k\\in\\{3,...,5\\} }(\\xi_{xk} \\zeta_{yk})^{r_s} \\right)^{r/{r_s}} \\right)^{1/r}.$$\n",
    "\n",
    "The *total output* is defined as $\\Phi _{xy}=\\alpha _{xy}+\\gamma _{xy}$. It measures the complemenarity between workers and firms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generating demand and supply\n",
    "\n",
    "We shall take 50 applicants, 30 positions, $r_g = .8$, $r_s = .6$ and $r = .7$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "nbx, nby = 50,30 #1200, 1000\n",
    "\n",
    "d, rg,rs,r = 8,.8, .6, .7\n",
    "np.random.seed(777)\n",
    "ξ_x_k = np.random.rand(nbx,d)\n",
    "ζ_y_k = np.random.rand(nby,d)\n",
    "α_x_y = np.zeros((nbx,nby))\n",
    "γ_x_y = np.zeros((nbx,nby))\n",
    "for x in range(nbx):\n",
    "    for y in range(nby):\n",
    "        α_x_y[x,y] = - np.linalg.norm(ξ_x_k[x,6:7]-ζ_y_k[y,6:7]) / 50\n",
    "        γ_x_y[x,y] = (np.sum( (ξ_x_k[x,0:3]*ζ_y_k[y,0:3])**rg )**(r/rg)+(  np.sum(  ξ_x_k[x,3:6]*ζ_y_k[y,3:6] )**rs )**(r/rs))**r\n",
    "        \n",
    "Φ_x_y = α_x_y + γ_x_y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The class `TU_model`\n",
    "\n",
    "We encode the problem into:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class TU_model:\n",
    "    def __init__(self,Φ_x_y,n_x=np.array([]),m_y=np.array([]),T=1):\n",
    "        if n_x.size == 0:\n",
    "            n_x = np.ones(Φ_x_y.shape[0])\n",
    "        if m_y.size == 0:\n",
    "            m_y = np.ones(Φ_x_y.shape[1])    \n",
    "        self.n_x = n_x\n",
    "        self.m_y = m_y\n",
    "        self.nbx = n_x.size\n",
    "        self.nby = m_y.size\n",
    "        self.Φ_x_y = Φ_x_y\n",
    "        self.T = T\n",
    "        self.K_x_y = np.exp(Φ_x_y / (2*T) )\n",
    "        self.eq_u_x = np.array([])\n",
    "        self.eq_v_y = np.array([])        \n",
    "        self.eq_μ_x0 = np.array([])\n",
    "        self.eq_μ_0y = np.array([])\n",
    "        self.eq_μ_x_y = np.array([])\n",
    "        self.eq_deltamarg_x = np.array([])\n",
    "        self.eq_deltamarg_y = np.array([])\n",
    "        self.eq_deltainvmarg_x = np.array([])\n",
    "        self.eq_deltainvmarg_y = np.array([])\n",
    "        self.comp_code = -1\n",
    "        self.comp_nbsteps = -1\n",
    "        self.comp_time = -1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create a `mkt` object with the data generated above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "mkt = TU_model(Φ_x_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wage determination\n",
    "\n",
    "Let $w_{xy}$ be the wage (endogenous, determined at equilibrium) that college $y$ would be willing to pay to employee $x$, to be determined. Start by assuming that there are no taxes, so if $x$ and $y$ match, then their payoff are respectively:\n",
    "$U_{xy} = \\alpha_{xy} + w_{xy}\\\\\n",
    "V_{xy}= \\gamma_{xy} - w_{xy} $\n",
    "while we assume that these payoffs are zero if $x$ and $y$ do not match."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matching patterns\n",
    "\n",
    "Assume that there are $n_{x}$ employees of type $x$, and $m_{y}$ colleges of type $y$.\n",
    "The number of matched $xy$ pairs, denoted  $\\mu_{xy}$, is determined at equilibrium. It should satisfy the *populations constraints*:\n",
    "\n",
    "$\\left\\{\n",
    "\\begin{array}{l}\n",
    "\\sum_{y}\\mu _{xy}\\leq n_{x} \\\\\n",
    "\\sum_{x}\\mu _{xy}\\leq m_{y}.%\n",
    "\\end{array}\n",
    "\\right. $\n",
    "\n",
    "In the sequel it will be convenient to denote $\\mathcal{X}_0 = \\mathcal{X} \\cup \\{ 0 \\}$  and $\\mathcal{Y}_0 = \\mathcal{Y} \\cup \\{ 0 \\}$, and introduce\n",
    "\n",
    "$\\left\\{\n",
    "\\begin{array}{l}\n",
    "\\mu_{x0} = n_x - \\sum_{y}\\mu _{xy}\\leq n_{x} \\\\\n",
    "\\mu_{0y} = m_y - \\sum_{x}\\mu _{xy}\\leq m_{y}.\n",
    "\\end{array}\n",
    "\\right. $\n",
    "\n",
    "the number of unassigned individuals of each type.\n",
    "\n",
    "Next, we shall be looking for equilibrium conditions on $\\mu_{xy}$ and $w_{xy}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Becker model\n",
    "\n",
    "Introduce $u_x$ and $v_y$ the indirect utilities of $x$ and $y$, respectively:\n",
    "\n",
    "$u_x = \\max_{y}\\left\\{ \\alpha _{xy}+w_{xy}, 0 \\right\\}$\n",
    "\n",
    "$v_y = \\max_{x}\\left\\{ \\gamma _{xy}-w_{xy},0\\right\\}$\n",
    "\n",
    "We have\n",
    "\n",
    "$u_x + v_y \\geq \\Phi_{xy}$, with equality if $\\mu_{xy}>0$, and\n",
    "\n",
    "$u_x \\geq 0$ with equality if $\\mu_{x0} >0$, and\n",
    "\n",
    "$v_y \\geq 0$ with equality if $\\mu_{0y} >0$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equilibrium in the Becker model\n",
    "\n",
    "To recap, $\\left( \\mu ,u,v\\right) $ is an equilibrium matching iff\n",
    "\n",
    "(1) quantities $\\mu $ satisfy the populations constraints:\n",
    "\n",
    "$\\left\\{\n",
    "\\begin{array}{l}\n",
    "\\sum_{y}\\mu _{xy} + \\mu_{x0} = n_{x} \\\\\n",
    "\\sum_{x}\\mu _{xy} + \\mu_{0y} = m_{y}\n",
    "\\end{array}\n",
    "\\right. $\n",
    "\n",
    "(2) utilities $(u,v)$ satisfy pairwise stability:\n",
    "\n",
    "$\\left\\{\n",
    "\\begin{array}{l}\n",
    "u_{x}+v_{y}\\geq \\Phi _{xy}, \\forall x,y \\\\\n",
    "u_{x}\\geq 0, \\forall x \\\\\n",
    "v_{y}\\geq 0, \\forall y\n",
    "\\end{array}\n",
    "\\right. $\n",
    "\n",
    "(3) the complementarity conditions hold:\n",
    "\n",
    "$\\left\\{\n",
    "\\begin{array}{l}\n",
    "\\mu _{xy}>0\\implies u_{x}+v_{y}=\\Phi _{xy} \\\\\n",
    "\\mu _{x0}>0\\implies u_{x}=0 \\\\\n",
    "\\mu _{0y}>0\\implies v_{y}=0.\n",
    "\\end{array}\n",
    "\\right. $\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear programming formulation\n",
    "\n",
    "Noting that the above conditions are complementary slackness conditions in linear programming, one has:\n",
    "\n",
    "**Theorem (Becker-Shapley-Shubik)**. In the equilibrium problem above, $\\mu$ and $(u,v)$ are the resepctive optimal solutions associated with the following primal and dual linear programming\n",
    "problems:<br>\n",
    "$\n",
    "\\begin{array}{l}\n",
    "\\max_{\\mu \\geq 0} &&\\sum_{xy}\\mu _{xy}\\Phi _{xy} \\\\\n",
    "&&\\sum_{y}\\mu _{xy}\\leq n_{x}~\\left[ u_{x}\\geq 0\\right] \\\\\n",
    "&&\\sum_{x}\\mu _{xy}\\leq m_{y}~\\left[ v_{y}\\geq 0\\right]\n",
    "\\end{array}\n",
    "$\n",
    "\n",
    "and<br>\n",
    "$\\begin{array}{l}\n",
    "\\min_{u\\geq 0,v\\geq 0} &&\\sum n_{x}u_{x}+\\sum m_{y}v_{y} \\\\\n",
    "s.t.~ &&u_{x}+v_{y}\\geq \\Phi _{xy}~\\left[ \\mu _{xy}\\geq 0\\right]\n",
    "\\end{array}\n",
    "$\n",
    "<br>\n",
    "<br>\n",
    "Interpretation as a *welfare theorem*: coincidence of the central planner's solution and the decentralized solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computation\n",
    "\n",
    "Setting $z=vec\\left(  \\mu\\right)$, the Linear Programming problem then becomes\n",
    "\n",
    "$\n",
    "\\begin{array}{l}\n",
    "&  \\max_{z\\geq0}vec\\left(  \\Phi\\right)  ^{\\top}z\\\\\n",
    "s.t.~  &  \\left( I_{\\mathcal{X}} \\otimes 1_{\\mathcal{Y} }^{\\top}\\right)  z\\leq n\n",
    "\\nonumber\\\\\n",
    "&  \\left( 1_{\\mathcal{X}}^{\\top} \\otimes I_{\\mathcal{Y}}\\right)  z \\leq m\n",
    "\\nonumber\n",
    "\\end{array}\n",
    "$\n",
    "\n",
    "which is ready to be passed on to a linear programming solver. See chapter 3 of [OTME]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solveBecker(self):\n",
    "    A1 = spr.kron(spr.identity(self.nbx),np.array(np.repeat(1,self.nby)))\n",
    "    A2 = spr.kron(np.array(np.repeat(1,self.nbx)),spr.identity(self.nby))\n",
    "    \n",
    "    A = spr.vstack([A1, A2])\n",
    "    d = np.concatenate((self.n_x,self.m_y), axis = None)\n",
    "    obj = Φ_x_y.flatten()\n",
    "    m=grb.Model()\n",
    "    m.Params.OutputFlag = 0\n",
    "    x = m.addMVar(shape=len(obj), name=\"x\")\n",
    "    m.setObjective(obj @ x, grb.GRB.MAXIMIZE)\n",
    "    m.addConstr(A @ x <= d)\n",
    "    m.optimize()                \n",
    "    \n",
    "    if m.status == grb.GRB.Status.OPTIMAL:\n",
    "        self.μ_x_y = np.array(m.getAttr('x')).reshape(self.nbx,self.nby)\n",
    "        self.μ_x0 = self.n_x - self.μ_x_y.sum(axis = 1)\n",
    "        self.μ_0y = self.m_y - self.μ_x_y.sum(axis = 0)\n",
    "        self.eq_u_x = np.array(m.getAttr('pi'))[0:self.nbx]\n",
    "        self.eq_v_y = np.array(m.getAttr('pi'))[self.nbx:(self.nbx+self.nby)]\n",
    "        self.comp_code = 0\n",
    "    else:\n",
    "        self.comp_code = -1\n",
    "    return self.comp_code\n",
    "\n",
    "TU_model.solveBecker = solveBecker "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute the model and access to $(u,v)$ in the following manner:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "--------------------------------------------\n",
      "Warning: your license will expire in 3 days\n",
      "--------------------------------------------\n",
      "\n",
      "Academic license - for non-commercial use only - expires 2021-08-15\n",
      "Using license file C:\\Users\\alfre\\gurobi.lic\n",
      " u[0:6]= [0.         0.         0.         0.19230066 0.10627652 0.05214333] \n",
      " v[0:6]= [1.48349568 1.28994117 1.35458853 1.60641385 1.48682969 1.6658138 ]\n"
     ]
    }
   ],
   "source": [
    "mkt.solveBecker()\n",
    "print(' u[0:6]=' , mkt.eq_u_x[0:6],'\\n v[0:6]=', mkt.eq_v_y[0:6])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Choo-Siow model\n",
    "\n",
    "\n",
    "Introduce logit random utilities in the problems of workers and firms. We have:\n",
    "\n",
    "\n",
    "$u_x = \\mathbb{E} \\left[ \\max_{y}\\left\\{ \\alpha _{xy}+w_{xy} + T \\varepsilon_y, T \\varepsilon_0 \\right\\} \\right] $ \n",
    "\n",
    "$v_y = \\mathbb{E} \\left[ \\max_{x}\\left\\{ \\gamma _{xy}-w_{xy} + T \\eta_y,T \\eta_0\\right\\} \\right]$\n",
    "\n",
    "where $\\epsilon$ and $\\eta$ are random vectors of i.i.d. Gumbel distributions. We have\n",
    "\n",
    "$u_{x}=T\\log \\left( 1+\\sum_{y}\\exp \\left( \\frac{\\alpha _{xy}+w_{xy}}{T}%\n",
    "\\right) \\right)  $\n",
    "\n",
    "$v_{y}=T\\log \\left( 1+\\sum_{x}\\exp \\left( \\frac{\\gamma _{xy}-w_{xy}}{T}%\n",
    "\\right) \\right). $\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Choice probabilities\n",
    "\n",
    "We have\n",
    "\n",
    "$\\frac {\\mu_{xy}} {n_x} = \\Pr \\left( y \\text{ chosen by } x \\right) =\\frac{\\exp \\left( \\frac{\\alpha\n",
    "_{xy}+w_{xy}}{T}\\right) }{1+\\sum_{y}\\exp \\left( \\frac{\\alpha _{xy}+w_{xy}}{T}\n",
    "\\right) } = \\exp \\left( \\frac{\\alpha _{xy}+w_{xy}-u_{x}}{T}\\right)$\n",
    "\n",
    "$ \\frac {\\mu_{x0}} {n_x} = \\Pr \\left( 0 \\text{ chosen by } x \\right) = \\frac{1}{1+\\sum_{y}\\exp \\left(\n",
    "\\frac{\\alpha _{xy}+w_{xy}}{T}\\right) } \n",
    "=\\exp \\left( -\\frac{u_{x}} {T} \\right) $\n",
    "\n",
    "and \n",
    "\n",
    "$\\frac {\\mu_{xy}} {m_y} = \\Pr \\left( x \\text{ chosen by } y \\right) =\\frac{\\exp \\left( \\frac{\\gamma\n",
    "_{xy} - w_{xy}}{T}\\right) }{1+\\sum_{x}\\exp \\left( \\frac{\\gamma _{xy} - w_{xy}}{T}\n",
    "\\right) } = \\exp \\left( \\frac{\\gamma _{xy} - w_{xy}-v_{y}}{T}\\right)$\n",
    "\n",
    "$ \\frac {\\mu_{0y}} {m_y} = \\Pr \\left( 0 \\text{ chosen by } y \\right) = \\frac{1}{1+\\sum_{x}\\exp \\left(\n",
    "\\frac{\\gamma _{xy} - w_{xy}}{T}\\right) } \n",
    "=\\exp \\left( -\\frac{v_{y}} {T} \\right) $\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matching functions\n",
    "\n",
    "Write \n",
    "\n",
    "$ \\frac {\\mu_{xy}^2 } {n_x m_y} = \\Pr \\left( y \\text{ chosen by } x \\right) . \\Pr \\left( x \\text{ chosen by } y \\right) \\\\\n",
    "=\\exp \\left( \\frac{\\alpha _{xy} + w_{xy}-u_{x}}{T}\\right) . \\exp \\left( \\frac{\\gamma _{xy} - w_{xy}-v_{y}}{T}\\right) =  \\exp \\left( \\frac{\\Phi _{xy} - u_{x}-v_{y}}{T}\\right)  $\n",
    "\n",
    "and as a result, $\\mu_{xy} = M_{xy}(\\mu_{x0}\\mu_{0y})$, where the functions<br>\n",
    "$M_{xy}(\\mu_{x0}\\mu_{0y}):= \\sqrt{\\mu_{x0} \\mu_{0y}} \\exp \\left( \\frac{\\Phi _{xy}}{2T}\\right) $\n",
    "\n",
    "relate the number of matched pairs of type $xy$ to number of unmatched agents of respective types $x$ and $y$. They are called *matching functions*.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Introduce<br>\n",
    "$ a_{x} = u_{x}-T\\ln n_{x}\\text{ and }b_{y} = v_{y}-T\\ln m_{y}$<br>\n",
    "so that we have\n",
    "\n",
    "$\\mu_{xy} = \\exp(\\frac {\\Phi_{xy} - a_x - b_y} {2T} )$<br>\n",
    "$\\mu_{x0} = \\exp(\\frac {- a_x} {T}),~\\mu_{0y} = \\exp(\\frac {- b_y} {T})$\n",
    "\n",
    "\n",
    "**Theorem (Choo-Siow)**. At equilibrium, the quantities $a_x$ and $b_y$ defined above solve the following system of equations<br>\n",
    "$\\left\\{\n",
    "\\begin{array}\n",
    "[c]{l}%\n",
    "n_x = \\exp(\\frac {- a_x} {T}) + \\sum_{y\\in \\mathcal{Y}}\\exp(\\frac {\\Phi_{xy} - a_x - b_y} {2T} )\\\\\n",
    "m_y = \\exp(\\frac {- b_y} {T}) + \\sum_{x\\in \\mathcal{X}}\\exp(\\frac {\\Phi_{xy} - a_x - b_y} {2T} )\n",
    "\\end{array}\n",
    "\\right.$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create a method `M_x_y` which we append to the `TU_model` class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def M_x_y(self,μ_x0,μ_0y):\n",
    "    μ_x_y = np.zeros((self.nbx,self.nby))\n",
    "    for x in range(self.nbx):\n",
    "        for y in range(self.nby):\n",
    "            μ_x_y[x,y] = self.K_x_y[x,y] * np.sqrt(μ_x0[x] * μ_0y[y])\n",
    "    return(μ_x_y)\n",
    "####################\n",
    "TU_model.M_x_y = M_x_y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Choo-Siow model as an optimization problem\n",
    "\n",
    "One can reformulate the  previous equations are the first order conditions of the following optimization problem:<br>\n",
    "<br>\n",
    "**Theorem (Galichon-Salanié)**. Choo and Siow's equations are the first order conditions associated to the following convex optimization problem:<br>\n",
    "$W(\\Phi) = \\min_{a,b} f(a,b) $<br>\n",
    "where<br>\n",
    "$f(a,b):= \\left\\{ \\sum_x a_x n_x + \\sum_y b_y m_y + 2T \\sum_{xy} e^{ \\frac {\\Phi_{xy} - a_x - b_y} {2T} }+ T \\sum_x e^{\\frac {-a_x} {T}} + T \\sum_y e^{\\frac {-b_y} {T}} \\right\\}$\n",
    "<br>\n",
    "whose dual program is:<br>\n",
    "$\\begin{array}{l}W(\\Phi) = \n",
    "\\max_\\mu &&\\sum_{x,y} \\mu_{xy} \\Phi_{xy}- 2T \\sum_{x,y} \\mu_{xy} \\log \\mu_{xy} - T \\sum_x \\mu_{x0} \\log \\mu_{x0} - T \\sum_y \\mu_{0y} \\log \\mu_{0y}\\\\\n",
    "s.t.~ &&\\sum_{y} \\mu_{xy} + \\mu_{x0}  = n_x \\\\\n",
    "&&\\sum_{x} \\mu_{xy} + \\mu_{0y} = m_y. \\\\\n",
    "\\end{array}\n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remark\n",
    "\n",
    "Note that, setting $A_x =\\exp(-a_x / (2T)) $ and $B_y =\\exp(-b_y / (2T)) $, the problem rewrites equivalently as<br>\n",
    "$\\min_{A,B} F(A,B)$<br>\n",
    "where<br>\n",
    "$F(A,B):= \\left\\{ - \\sum_x n_x \\log A_x  - \\sum_y m_y \\log B_y + 2 \\sum_{xy} K_{xy} A_x B_y + \\sum_x A^2_x + \\sum_y B^2_y \\right\\}$<br>\n",
    "where $K_{xy} = T \\exp( \\Phi_{xy} / (2T) ).$<br>\n",
    "One has:<br>\n",
    "$\\frac{ \\partial F} {\\partial A_x}  = 2 A_x - \\frac {n_x} {A_x} + \\sum_y K_{xy} B_y$<br>\n",
    "$\\frac{ \\partial F} {\\partial B_y}  = 2 B_y - \\frac {m_y} {B_y} + \\sum_x K_{xy} A_x,$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize\n",
    "from time import time\n",
    "\n",
    "def f(self,ab_z):\n",
    "    a_x = ab_z[0:self.nbx]\n",
    "    b_y = ab_z[self.nbx:(self.nbx+self.nby)]\n",
    "    μ_x0 = np.exp( - a_x /  self.T )\n",
    "    μ_0y = np.exp( - b_y / self.T )\n",
    "    μ_x_y = self.M_x_y(μ_x0 , μ_0y)  \n",
    "    return( np.sum(self.n_x * a_x ) + np.sum(self.m_y * b_y ) + 2* np.sum( μ_x_y) + np.sum(μ_x0 ) + np.sum(μ_0y ) )\n",
    "\n",
    "TU_model.f = f\n",
    "\n",
    "def gradf_xy(self,ab_z):\n",
    "    a_x = ab_z[0:self.nbx]\n",
    "    b_y = ab_z[self.nbx:(self.nbx+self.nby)]\n",
    "    μ_x0 = np.exp( - a_x /  self.T )\n",
    "    μ_0y = np.exp( - b_y / self.T )\n",
    "    μ_x_y = self.M_x_y(μ_x0 , μ_0y)  \n",
    "    dGda_x = self.n_x - μ_x0 - np.sum( μ_x_y , 1)\n",
    "    dGdb_y = self.m_y - μ_0y - np.sum( μ_x_y , 0)\n",
    "    return(np.append(dGda_x,dGdb_y))\n",
    "\n",
    "TU_model.gradf_xy = gradf_xy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now build a method to solve for equilibrium via the optimization problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solveChooSiowOptim(self):\n",
    "    start_time = time()\n",
    "    uv0_z = np.zeros(self.nbx+self.nby)\n",
    "    res = minimize(self.f, uv0_z,  method='BFGS', jac = self.gradf_xy)\n",
    "    if res.success == False:\n",
    "        print('Failed to converge')\n",
    "    else:\n",
    "        comp_time = time() - start_time\n",
    "        a_x = res.x[0:self.nbx]\n",
    "        b_y = res.x[self.nbx:(self.nbx+self.nby)]\n",
    "        μ_x0 = np.exp( - a_x /  self.T )\n",
    "        μ_0y = np.exp( - b_y / self.T )\n",
    "        μ_x_y = self.M_x_y(μ_x0,μ_0y) \n",
    "        self.eq_u_x = a_x + np.log(self.n_x) \n",
    "        self.eq_v_y = b_y + np.log(self.m_y)\n",
    "        self.eq_μ_x0 = np.exp(- a_x / self.T)\n",
    "        self.eq_μ_0y = np.exp(- b_y / self.T)\n",
    "        self.eq_μ_x_y = μ_x_y       \n",
    "        self.eq_deltamarg_x = self.n_x - μ_x0 - np.sum( μ_x_y , 1)\n",
    "        self.eq_deltamarg_y = self.m_y - μ_0y - np.sum( μ_x_y , 0)\n",
    "        self.comp_time = comp_time\n",
    "    return(res.status)\n",
    "\n",
    "TU_model.solveChooSiowOptim = solveChooSiowOptim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run it on our example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " u[0:6]= [0.90023101 0.84255063 0.870189   0.99747671 0.96463365 0.9613982 ] \n",
      " v[0:6]= [8.19888742 8.03456204 8.13410101 8.36745667 8.25594198 8.38963211]\n"
     ]
    }
   ],
   "source": [
    "mkt.solveChooSiowOptim()\n",
    "print(' u[0:6]=' , mkt.eq_u_x[0:6],'\\n v[0:6]=', mkt.eq_v_y[0:6])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Choo-Siow model as an equilibrium problem with GS\n",
    "\n",
    "Can we reformulate as a competitive equilibrium problem with Gross Substitutes? Let's explore the model we are trying to solve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.19274941, 0.87381862, 0.31693207, 0.15259627],\n",
       "       [0.02816756, 0.62162143, 0.50866174, 0.15624843],\n",
       "       [0.00241425, 0.21030363, 0.34583779, 0.38794661]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.rand(3,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[3.3 0.  0.  ... 0.  0.  0. ]\n",
      " [0.  4.  0.  ... 0.  0.  0. ]\n",
      " [0.  0.  3.3 ... 0.  0.  0. ]\n",
      " ...\n",
      " [0.  0.  0.  ... 0.  0.  0. ]\n",
      " [0.  0.  0.  ... 0.  0.  0. ]\n",
      " [0.  0.  0.  ... 0.  0.  0. ]]\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(777)\n",
    "\n",
    "import numdifftools as nd\n",
    "toymodel = TU_model(np.array( np.random.rand(3,4)))\n",
    "theFun = lambda arg: (toymodel.f(arg) )\n",
    "theHessian = nd.Hessian(theFun)(np.zeros(nbx+nby)).round(1)\n",
    "H_11 = theHessian[0:nbx,0:nbx]\n",
    "H_12 = theHessian[0:nbx,nbx:(nbx+nby)]\n",
    "H_21 = theHessian[nbx:(nbx+nby),0:nbx]\n",
    "H_22 = theHessian[nbx:(nbx+nby),nbx:(nbx+nby)]\n",
    "print(theHessian)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the matrix is of the form\n",
    "\n",
    "$\\begin{pmatrix}\n",
    "H_{11} & H_{12} \\\\\n",
    "H_{21} & H_{22}\n",
    "\\end{pmatrix}$\n",
    "\n",
    "where $H_{11}$ and $H_{22}$ are diagonal matrices with positive diagonal terms, and $H_{12}$ and $H_{21}$ are matrices with nonnegative entries.\n",
    "\n",
    "Further, note that $H_{12}$ and $H_{21}$ are transpose to each other.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 0. 0. ... 0. 0. 0.]\n",
      " [0. 0. 0. ... 0. 0. 0.]\n",
      " [0. 0. 0. ... 0. 0. 0.]\n",
      " ...\n",
      " [0. 0. 0. ... 0. 0. 0.]\n",
      " [0. 0. 0. ... 0. 0. 0.]\n",
      " [0. 0. 0. ... 0. 0. 0.]]\n",
      "Norm of H_21^T - H_12 = 0.0\n"
     ]
    }
   ],
   "source": [
    "print(H_12)\n",
    "print('Norm of H_21^T - H_12 =', np.max(np.abs(H_21.T - H_12)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The change-of-sign trick\n",
    "\n",
    "Define $p\\in \\mathbb{R}^{\\mathcal{X}\\cup\\mathcal{Y}}$ by<br>\n",
    "$p_x = a_x, x\\in\\mathcal{X}$<br>\n",
    "$p_y = -b_y, y\\in\\mathcal{Y}$<br>\n",
    "and define $q \\in \\mathbb{R}^{\\mathcal{X}\\cup\\mathcal{Y}}$ by<br>\n",
    "$q_x = -n_x, x\\in\\mathcal{X}$<br>\n",
    "$q_y = b_y, y\\in\\mathcal{Y}.$<br>\n",
    "\n",
    "Define<br>\n",
    "$\\left\\{\n",
    "\\begin{array}\n",
    "[c]{l}%\n",
    "Q_x(p) = - \\exp(\\frac {- p_x} {T}) - \\sum_{y\\in \\mathcal{Y}}\\exp(\\frac {\\Phi_{xy} - p_x + p_y} {2T} )\\\\\n",
    "Q_y(p) = \\exp(\\frac {p_y} {T}) + \\sum_{x\\in \\mathcal{X}}\\exp(\\frac {\\Phi_{xy} - p_x +p_y} {2T} )\n",
    "\\end{array}\n",
    "\\right. $\n",
    "\n",
    "so that the equilibrium problem reformulates as<br>\n",
    "$Q(p)=q.$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In terms of Jacobian of the system, this amounts to \n",
    "\n",
    "$\\begin{pmatrix}\n",
    "H_{11} & - H_{12} \\\\\n",
    "- H_{21} & H_{22}\n",
    "\\end{pmatrix}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coordinate update algorithm\n",
    "\n",
    "Let's initialize $p_{x}^{0}=-T\\ln n_{x}$ and $p_{y}^{0}=-\\infty $.\n",
    "\n",
    "Then update $p_{y}^{1}$ so that $Q_{y}\\left( \\left( p_{x}^{0}\\right)\n",
    "_{x},p_{y}^{1}\\right) =q_{y}$($=m_{y})$.\n",
    "\n",
    "we solve $\\sum \\exp \\left( \\Phi _{xy}-p_{x}^{0}+p_{y}^{1}\\right) +\\exp\n",
    "\\left( p_{y}^{1}\\right) =m_{y}$\n",
    "\n",
    "Because the value of $p_{y}^{1}$ is finite, we have $p_{y}^{0}=-\\infty \\leq p_{y}^{1}$\n",
    "\n",
    "Then update $p_{x}^{1}$ so that $Q_{x}\\left( p_{x}^{1},\\left(\n",
    "p_{y}^{1}\\right) _{y}\\right) =q_{x}$($=-n_{x}$). We have\n",
    "\n",
    "$\\sum_{y}\\exp \\left( \\Phi _{xy}-p_{x}^{1}+p_{y}^{1}\\right) +\\exp \\left( -%\n",
    "\\frac{p_{x}^{1}}{T}\\right) =n_{x}\\geq \\exp \\left( -\\frac{p_{x}^{1}}{T}%\n",
    "\\right) $\n",
    "\n",
    "thus $\\exp \\left( -\\frac{p_{x}^{0}}{T}\\right) =n_{x}\\geq \\exp \\left( -\\frac{%\n",
    "p_{x}^{1}}{T}\\right) $\n",
    "\n",
    "thus $p_{x}^{0}\\leq p_{x}^{1}$\n",
    "\n",
    "We have $m_{y}\\geq \\exp \\left( p_{y}/T\\right) $, thus $p_{y}\\leq T\\ln m_{y}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set $p_{y}=-b_{y}$ for $y\\in Y$\n",
    "\n",
    "$\\left\\{\n",
    "\\begin{array}{l}\n",
    "Q_{x}\\left( p\\right) = -\\sum_{y}\\exp \\left( \\frac{\\Phi _{xy}-p_{x}+p_{y}}{2T}\\right) -\\exp \\left( -\\frac{p_{x}}{T}\\right) \\\\\n",
    "Q_{y}\\left( p\\right) = \\sum_{x}\\exp \\left( \\frac{\\Phi _{xy}-p_{x}+p_{y}}{2T}%\n",
    "\\right) +\\exp \\left( \\frac{p_{y}}{T}\\right)\n",
    "\\end{array}\n",
    "\\right.$\n",
    "\n",
    "We have\n",
    "\n",
    "$DQ=\n",
    "\\begin{pmatrix}\n",
    "\\frac{\\partial e_{x}}{\\partial p_{x}} & \\frac{\\partial e_{x}}{\\partial p_{y}}\n",
    "\\\\\n",
    "\\left( \\frac{\\partial e_{x}}{\\partial p_{y}}\\right) ^{\\top } & \\frac{%\n",
    "\\partial e_{y}}{\\partial p_{y}}%\n",
    "\\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "diag\\left( \\frac{1}{2T}\\exp \\left( \\frac{\\Phi _{xy}-p_{x}+p_{y}}{2T}\\right) +\n",
    "\\frac{1}{T}\\exp \\left( -\\frac{p_{x}}{T}\\right) \\right) &  \\\\\n",
    "-\\frac{1}{2T}\\exp \\left( \\frac{\\Phi _{xy}-p_{x}+p_{y}}{2T}\\right) &\n",
    "diag\\left( \\frac{1}{2T}\\exp \\left( \\frac{\\Phi _{xy}-p_{x}+p_{y}}{2T}\\right) +\n",
    "\\frac{1}{T}\\exp \\left( \\frac{p_{y}}{T}\\right) \\right)\n",
    "\\end{pmatrix}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gauss-Seidel algorithm\n",
    "\n",
    "The algorithm runs as follows:\n",
    "\n",
    "* Solve for $a_{x}$ in the 1st set of equations\n",
    "$\\sum_{y}\\exp \\left( \\frac{\\Phi _{xy}-a_{x}-b_{y}}{2T}\\right) + \\exp \\left( -\\frac{a_{x}}{T}\\right) = n_{x}\n",
    "$\n",
    "\n",
    "* Solve for $b_{y}$ in the 2nd set of equations\n",
    "$\\sum_{x}\\exp \\left( \\frac{\\Phi _{xy}-a_{x}-b_{y}}{2T}\\right) +\\exp \\left( -\\frac{b_{y}}{T}\\right) =m_{y}\n",
    "$\n",
    "\n",
    "These equations can be solved explicitly. Introduce $K_{xy}=\\exp \\left( \\frac{\\Phi _{xy}}{2T}\\right)$ and new unknowns<br>\n",
    "$A_{x}=\\exp \\left( -\\frac{a_{x}}{2T}\\right),~B_{y}=\\exp \\left( -\\frac{b_{y}}{2T}\\right)$\n",
    "\n",
    "we have\n",
    "\n",
    "$\\sum_{y}K_{xy}A_{x}B_{y}+A_{x}^{2} = n_{x} \\\\\n",
    "\\sum_{x}K_{xy}A_{x}B_{y}+B_{y}^{2} = m_{y}\n",
    "$\n",
    "\n",
    "thus\n",
    "\n",
    "$\n",
    "A_{x} = \\sqrt{n_{x}+\\left( \\frac{1}{2}\\sum_{y}K_{xy}B_{y}\\right) ^{2}}-\n",
    "\\frac{1}{2}\\sum_{y}K_{xy}B_{y} \\\\\n",
    "B_{y} = \\sqrt{m_{y}+\\left( \\frac{1}{2}\\sum_{x}K_{xy}A_{y}\\right) ^{2}}-\n",
    "\\frac{1}{2}\\sum_{x}K_{xy}A_{y}\n",
    "$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create a method implementing the Gauss-Seidel algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solveChooSiowIPFP(self, maxit = 1000,maxvaltol=1e-5,maxsteptol=1e-9,output=0):\n",
    "    start_time = time()\n",
    "    code = 0\n",
    "    deltainvmarg_x = np.zeros(self.nbx)\n",
    "    deltainvmarg_y = np.zeros(self.nby)\n",
    "    A_x = np.ones(self.nbx)\n",
    "    B_y = np.ones(self.nby)\n",
    "    for i in range(maxit):\n",
    "        \n",
    "        KThalfA_y = (self.K_x_y.T @ A_x) / 2\n",
    "        newB_y = np.sqrt(self.m_y + KThalfA_y * KThalfA_y ) - KThalfA_y\n",
    "        deltainvmarg_y = newB_y - B_y\n",
    "        B_y = newB_y\n",
    "        \n",
    "        KBhalf_x = (self.K_x_y @ B_y) / 2\n",
    "        newA_x = np.sqrt(self.n_x + KBhalf_x * KBhalf_x ) - KBhalf_x\n",
    "        deltainvmarg_x = newA_x - A_x\n",
    "        A_x = newA_x\n",
    "        \n",
    "        deltamarg_y = (self.K_x_y.T @ A_x) * B_y  + B_y * B_y  - self.m_y\n",
    "        steptol = np.max(np.abs( np.append(deltainvmarg_x,deltainvmarg_y) ) ) \n",
    "        valtol = np.max(np.abs(deltamarg_y) )\n",
    "        \n",
    "        if output > 1 :\n",
    "            print(\"μ_x0=\",μ_x0)\n",
    "        if valtol < maxvaltol :\n",
    "            code = 0\n",
    "            break\n",
    "        if steptol < maxsteptol :\n",
    "            code = 1\n",
    "            break\n",
    "        code = 2    \n",
    "    comp_time = time() - start_time\n",
    "    if output > 0 :\n",
    "        print( 'IPFP converged in', i, 'iterations and ',comp_time,' seconds.')\n",
    "        print('Max absolute margin discrepancy=', valtol)            \n",
    "        print('Max absolute update discrepancy =',steptol)\n",
    "        print('Code =',code)\n",
    "\n",
    "    μ_x0 = A_x * A_x\n",
    "    μ_0y = B_y * B_y\n",
    "    μ_x_y = self.M_x_y(μ_x0,μ_0y)\n",
    "    self.eq_μ_x0 = μ_x0\n",
    "    self.eq_μ_0y = μ_0y\n",
    "    self.eq_μ_x_y = μ_x_y\n",
    "    self.eq_u_x = - self.T * np.log(μ_x0 / self.n_x)\n",
    "    self.eq_v_y = - self.T * np.log(μ_0y / self.m_y)\n",
    "    self.eq_deltamarg_x = self.n_x - μ_x0 - np.sum( μ_x_y , 1)\n",
    "    self.eq_deltamarg_y = self.m_y - μ_0y - np.sum( μ_x_y , 0)\n",
    "    self.eq_deltainvmarg_x = deltainvmarg_x\n",
    "    self.eq_deltainvmarg_y = deltainvmarg_y\n",
    "    self.comp_code = code\n",
    "    self.comp_nbsteps = i\n",
    "    self.comp_time = comp_time\n",
    "    return code\n",
    "\n",
    "TU_model.solveChooSiowIPFP = solveChooSiowIPFP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run it using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IPFP converged in 12 iterations and  0.0  seconds.\n",
      "Max absolute margin discrepancy= 7.89435480441636e-06\n",
      "Max absolute update discrepancy = 5.22727313956306e-06\n",
      "Code = 0\n",
      "u[0:6]= [0.90021933 0.84253957 0.87017764 0.99746391 0.96462125 0.96138585] \n",
      "v[0:6]= [8.19891503 8.03458912 8.13412851 8.3674861  8.25597123 8.38965999]\n"
     ]
    }
   ],
   "source": [
    "mkt.solveChooSiowIPFP(output=1)\n",
    "print('u[0:6]=' , mkt.eq_u_x[0:6],'\\nv[0:6]=', mkt.eq_v_y[0:6])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Matching with linear taxes\n",
    "\n",
    "We now consider a version of the Becker model with a \"flat tax\": assume that the gross wage $w_{xy}$ is taxed at a constant rate $\\tau \\in [0,1)$:\n",
    "$U_{xy} = \\alpha_{xy} + (1-\\tau) w_{xy}\\\\\n",
    "V_{xy} = \\gamma_{xy} - w_{xy}.$<br>\n",
    "\n",
    "Dupuy, Galichon, Jaffe and Kominers (2020) have shown the following result:\n",
    "\n",
    "**Theorem (DGJK)**. The equilibrium matching in the problem above with a flat tax $\\tau$ above is the optimal matching in a TU matching problem with surplus function<br>\n",
    "$\\Phi^\\tau_{xy} =  \\alpha_{xy}   +  (1-\\tau )  \\gamma_{xy}.$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Proof\n",
    "\n",
    "At equilibrium there is a wage $w_{xy}$ such that:\n",
    "\n",
    "$x$ gets indidirect utility\n",
    "$u_{x}=\\max_{y}\\left\\{ \\alpha _{xy}+\\left( 1-\\tau \\right) w_{xy},0\\right\\} $, and\n",
    "\n",
    "$y$ gets indirect utility\n",
    "$v_{y}=\\max_{x}\\left\\{ \\gamma _{xy}-w_{xy},0\\right\\} $\n",
    "\n",
    "\n",
    "Equilibrium consists thus of looking for matching patterns $\\mu _{xy}$, wages $w_{xy}$, and indirect utilities $(u_x,v_y)$ such that:\n",
    "\n",
    "(i) populations constraints are met:<br>\n",
    "$\\left\\{\n",
    "\\begin{array}{l}\n",
    "\\sum_{y}\\mu _{xy}+\\mu _{x0}=n_{x} \\\\\n",
    "\\sum_{x}\\mu _{xy}+\\mu _{0y}=m_{y}\n",
    "\\end{array}\n",
    "\\right. $\n",
    "\n",
    "(ii) stability conditions hold:<br>\n",
    "$u_{x}\\geq \\alpha _{xy}+\\left( 1-\\tau \\right) w_{xy}\\text{ and }\n",
    "v_{y}\\geq \\gamma _{xy} - w_{xy} \\\\\n",
    "u_{x}\\geq 0,v_{y}\\geq 0$ \n",
    "\n",
    "(iii) complementary conditions hold:<br>\n",
    "$\\mu_{xy}>0 \\implies u_{x} = \\alpha _{xy}+\\left( 1-\\tau \\right) w_{xy}\\text{ and }\n",
    "v_{y} = \\gamma _{xy} - w_{xy} \\\\\n",
    "\\mu_{x0}>0 \\implies u_{x} =0 \\\\\n",
    "\\mu_{0y}>0 \\implies v_{y} = 0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Proof (ctd)\n",
    "\n",
    "Denote $\\tilde{v}_{y}=\\left( 1-\\tau \\right) v_{y}$ and $\\tilde{\\gamma}%\n",
    "_{xy}=\\left( 1-\\tau \\right) \\gamma _{xy}$ and $\\tilde{w}_{xy}=\\left( 1-\\tau\n",
    "\\right) w_{xy}$ the indirect utility of the firm and the output measured in\n",
    "post-tax dollars, then we have\n",
    "\n",
    "$\\sum_{y}\\mu _{xy}+\\mu _{x0}=n_{x}$\n",
    "\n",
    "$\\sum_{x}\\mu _{xy}+\\mu _{0y}=m_{y}$\n",
    "\n",
    "$u_{x}\\geq 0,\\tilde{v}_{y}\\geq 0$\n",
    "\n",
    "$u_{x}\\geq \\alpha _{xy}+\\tilde{w}_{xy}$ with equality if $\\mu _{xy}>0$\n",
    "\n",
    "$\\tilde{v}_{y}\\geq \\tilde{\\gamma}_{xy}-\\tilde{w}_{xy}$ with equality if $\\mu\n",
    "_{xy}>0$.\n",
    "\n",
    "Therefore, $\\mu$ is the optimal matching associated with surplus function<br>\n",
    "$\\tilde{\\Phi}_{xy} = \\alpha_{xy}+\\left( 1-\\tau \\right) \\gamma _{xy}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Embedding in a Choo-Siow model\n",
    "\n",
    "Consider the logit random utility version of the previous problem:\n",
    "\n",
    "\n",
    "$u_x = \\mathbb{E} \\left[ \\max_{y}\\left\\{ \\alpha _{xy}+(1 - \\tau) w_{xy} + T \\varepsilon_y, T \\varepsilon_0 \\right\\} \\right] $ \n",
    "\n",
    "$v_y = \\mathbb{E} \\left[ \\max_{x}\\left\\{ \\gamma _{xy}-w_{xy} + T \\eta_y,T \\eta_0\\right\\} \\right]$\n",
    "\n",
    "where $\\epsilon$ and $\\eta$ are random vectors of i.i.d. Gumbel distributions. \n",
    "\n",
    "Defining as before\n",
    "$\\tilde{v}_{y}=\\left( 1-\\tau \\right) v_{y}$ and $\\tilde{\\gamma}_{xy}=\\left( 1-\\tau \\right) \\gamma _{xy}$ and $\\tilde{w}_{xy}=\\left( 1-\\tau\n",
    "\\right) w_{xy}$\\, we have<br>\n",
    "$u_{x}=T\\log \\left( 1+\\sum_{y}\\exp \\left( \\frac{\\alpha _{xy}+\\tilde{w}_{xy}}{T}\\right) \\right)  \\\\\n",
    "\\tilde{v}_{y}=(1-\\tau) T\\log \\left( 1+\\sum_{x}\\exp \\left( \\frac{\\tilde{\\gamma} _{xy}-\\tilde{w}_{xy}}{(1-\\tau)T}\n",
    "\\right) \\right). $\n",
    "\n",
    "As a result, the model is a Choo-Siow model with surplus $\\alpha_{xy} + (1-\\tau) \\gamma_{xy}$ and scaling parameter $1-\\tau / 2$\n",
    "\n",
    "Define:<br>\n",
    "$\\begin{array}{l}W\\left( \\theta ,\\lambda \\right)  = \n",
    "\\max_\\mu &&\\sum_{x,y} \\mu_{xy} \\left( \\alpha \\theta +\\gamma\n",
    "\\lambda \\right)- (\\theta+\\lambda)T \\sum_{x,y} \\mu_{xy} \\log \\mu_{xy} - \\theta T \\sum_x \\mu_{x0} \\log \\mu_{x0} -  \\lambda T \\sum_y \\mu_{0y} \\log \\mu_{0y}\\\\\n",
    "s.t.~ &&\\sum_{y} \\mu_{xy} + \\mu_{x0}  = n_x \\\\\n",
    "&&\\sum_{x} \\mu_{xy} + \\mu_{0y} = m_y \\\\\n",
    "\\end{array}$<br>\n",
    "so that the equilibrium matching model with taxes solves for an equilibrium with $\\theta = 1$ and $\\lambda = 1-\\tau$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparative statics\n",
    "\n",
    "**Proposition (DGJK)**. $W$ is a convex and positive homogenous function, and we have:<br>\n",
    "$\\left\\{\n",
    "\\begin{array}{l}\n",
    "\\frac{\\partial W}{\\partial \\theta } (\\theta,\\lambda)=\\sum_{xy} \\mu_{xy} \\alpha_{xy} -T \\sum_x \\mu_{x0} \\log \\mu_{x0} \n",
    "=:A\\\\\n",
    "\\frac{\\partial W}{\\partial \\lambda } (\\theta,\\lambda)=\\sum_{xy} \\mu_{xy} \\gamma_{xy} -T \\sum_y \\mu_{0y} \\log \\mu_{0y} \n",
    "=:\\Gamma\n",
    "\\end{array}\n",
    "\\right.\n",
    "$\n",
    "\n",
    "therefore $\\frac{\\partial \\Gamma }{\\partial \\lambda }\\geq 0$\n",
    "\n",
    "Similarly, $W=\\theta A+\\lambda \\Gamma $ and $\\Gamma =\\theta \\frac{\\partial A%\n",
    "}{\\partial \\lambda }+\\Gamma +\\lambda \\frac{\\partial \\Gamma }{\\partial\n",
    "\\lambda }$\n",
    "\n",
    "hence $\\frac{\\partial \\Gamma }{\\partial \\lambda }=-\\frac{\\theta }{\\lambda }%\n",
    "\\frac{\\partial A}{\\partial \\lambda }$\n",
    "\n",
    "so with $\\theta =1$, one has\n",
    "\n",
    "$\\frac{\\partial A}{\\partial \\lambda }=-\\lambda \\frac{\\partial \\Gamma }{%\n",
    "\\partial \\lambda }\\leq 0$\n",
    "\n",
    "and\n",
    "\n",
    "$\\frac{\\partial \\left( A+\\Gamma \\right) }{\\partial \\lambda }=\\left(\n",
    "1-\\lambda \\right) \\frac{\\partial \\Gamma }{\\partial \\lambda }$\n",
    "\n",
    "As a result, the sum of amenities $A$ increases with the tax wage, but the total welfare $A+\\Gamma$ decreases, so the increase in $A$ is not enough to compensate the decrease in $\\Gamma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEDCAYAAAAiKuN6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAApr0lEQVR4nO3dd3xV9f3H8deHMMJGhohh76Esw1DqoE5wUKu02p8iOHDR4a9YwWpddVNRW39YrLutgyGioqg4qJugEBIwEIYQiCyRkZCQ8fn9cS/1GiO5gdzcm3vfz8cjD3LO+d57P+cQ8uZzprk7IiIi4aoV7QJERKRmUXCIiEilKDhERKRSFBwiIlIpCg4REakUBYeIiFSKggMws9FmlmlmpWaWeoBxT5jZFjPLKDP/BTNbEvxaZ2ZLKvi8/mb2cfAz083sl1W0KiIiEZdwwWFmJ5nZU2VmZwA/BxZW8PKngDPKznT3X7p7f3fvD8wCZlfwPvnAGHfvE3y/B82sWYXFi4jEgIQLjvK4+wp3zwpj3ELgmx9bbmYG/AJ4LjidZGb3m9miYGdxZfB9Vrr7quD3m4AtQKsqWBURkYirHe0C4szxwOb9oQBcBux090FmVg/40MzedPe1+19gZoOBusDq6i9XRKTyEiY4zOxToB7QCGgechziBnefX0UfcyHBbiPoNKCvmZ0fnG4KdAPWBmtqAzwLXOLupVVUg4hIRCVMcLj7EAgc4wDGuvvYqnx/M6tN4DjJMaGzgV+XF0xm1gR4DbjJ3T+pylpERCJJxziqzinAl+6eEzJvPnC1mdUBMLPuZtbQzOoCLwHPuPuMKNQqInLQFByAmZ1rZjnAscBrZjY/OP9IM5sXMu454GOgh5nlmNllIW9zAd/fTQXwD2A58HnwFN6/E+jyfgGcAIwNOY23f4RWT0SkSpluqy4iIpWhjkNERColIQ6Ot2zZ0jt27BjtMkREapTFixdvc/cfXGOWEMHRsWNH0tLSol2GiEiNYmZflTc/oruqzOwMM8sys2wzm1TOcjOzh4PL081sYMiy35pZRvB+Tr8Lmd/czN4ys1XBPw+L5DqIiMj3RSw4zCwJeAQYAfQGLjSz3mWGjSBwQVw3YDwwLfjao4ArgMFAP+AsM+sWfM0kYIG7dwMWBKdFRKSaRLLjGAxku/sad98HPA+MKjNmFIFrGTx4EVyz4NXUvYBP3D3f3YuB94FzQ17zdPD7p4GfRXAdRESkjEge40gBNoRM5wBDwhiTQuButXeaWQtgLzAS2H+QorW75wK4e66ZHV7eh5vZeAJdDO3bt//B8qKiInJycigoKKjkasWn5ORk2rZtS506daJdiojEuEgGh5Uzr+xFI+WOcfcVZnYv8BawB1gKFFfmw919OjAdIDU19QcXq+Tk5NC4cWM6duxI4Ka2icvd2b59Ozk5OXTq1Cna5YhIjIvkrqocoF3IdFtgU7hj3P1xdx/o7icQuJX5/jvObg7uztp/k8AtB1NcQUEBLVq0SPjQADAzWrRooe5LRMISyeBYBHQzs07BezNdAMwtM2YuMCZ4dtVQArcgzwXYvwvKzNoTuHngcyGvuST4/SXAywdboELjO9oWIhKuiAVH8KD2BAI3+lsBvOjumWZ2lZldFRw2D1gDZAOPAdeEvMUsM1sOvAJc6+47gvPvAU41s1XAqcFpEREJsW1PIbe9ksmugqIqf++IXgDo7vMIhEPovEdDvnfg2h957fE/Mn87cHIVlhkV27dv5+STA6vx9ddfk5SURKtWgQs0P/vsM+rWrfvfsQ8++CDjx4+nQYMGB3zPk046iSlTppCa+qOPTReROLevuJRnPl7HQ2+vYm9RCcO6tOSU3q2r9DMS4srxWNSiRQuWLFkCwK233kqjRo2YOHFiuWMffPBBLrroogqDQ0QS23tZW7j91eWs2ZrHid1bcfNZvel6eKMq/xzd5DCGLFiwgAEDBnD00Udz6aWXUlhYyMMPP8ymTZsYPnw4w4cPB+Dqq68mNTWVPn36cMstt1T4vvPmzaNfv3707NmTuXPnMmTIEIqLK3WSmojEsHXb8rjsqUWMfXIRpaXO45ek8tS4QREJDVDHAcBtr2SyfNOuKn3P3kc24Zaz+4Q9vqCggLFjx7JgwQK6d+/OmDFjmDZtGr/73e944IEHePfdd2nZsiUAd955J82bN6ekpISTTz6Z9PR0+vbtW+775ufnc/nll7No0SIAhgwZwvHHH0/t2vqrF6np9hQW89d3VvHEB2upm1SLSSN6Mm5YR+rVToro56rjiBElJSV06tSJ7t27A3DJJZewcOHCcse++OKLDBw4kAEDBpCZmcny5ct/9H2zsrLo3LkzKSkppKSk0LJlS84+++yIrIOIVI/SUmfm4hyGT3mPv7+/hnP6pfDuxJO46sQuEQ8NUMcBUKnOIFIaNmwY1ri1a9cyZcoUFi1axGGHHcbYsWMrvP6iXr16//2+cePGDBgw4JBqFZHoWbLhW26dm8mSDd/Sr10zpl98DAPaV++9XtVxxIiCggLWrVtHdnY2AM8++ywnnngiEPhlv3v3bgB27dpFw4YNadq0KZs3b+b1118/4Pt2796dNWvWkJeXx8aNG/nss8/IyMiI7MqISJXbsruAiTOW8rNHPmTjt3uZMrofL119XLWHBqjjiBnJyck8+eSTjB49muLiYgYNGsRVVwUudxk/fjwjRoygTZs2vPvuuwwYMIA+ffrQuXNnhg0bdsD3bdiwIXfddRdnnXUW+fn5vPDCC9x777107NiRQYMGVceqicgh2FdcypMfruWv72RTWFzClSd2ZsLwrjROjt595RLimeOpqale9kFOK1asoFevXlGqKDZpm4jElne+3Mwdr65g7bY8ftrzcG46sxedW0XmTKnymNlid//BhWHqOEREYszqrXu449XlvJe1lc4tG/LkuEEM71HujcCjQsEhIhIjdhcU8dd3snnyw7XUq53EH0f24pLjOlK3dmwdjk7o4HB33dwvKBF2WYrEqtJSZ+bnOdz3Rhbb8woZfUxbrj+9J60a16v4xVGQsMGRnJzM9u3bdWt1vnseR3JycrRLEUk4n6/fwW1zM1mas5OB7ZvxxNhU+rZtFu2yDihhg6Nt27bk5OSwdevWaJcSE/Y/AVBEqsfmXQXc+/qXzP5iI4c3rsfUX/ZjVL8UatWK/f/IJmxw1KlTR0+7E5FqV1hcwuMfrOWRd7IpKnGuPqkL1w7vSqN6NefXcc2pVESkBnN3FqzYwh2vLeer7fmc0qs1N53Zi44tw7trRCxRcIiIRFj2lj3c/upyFq7cSpdWDXnm0sGc0L1VtMs6aAoOEZEI2VVQxENvr+Lpj9ZRv24SN5/VmzHHdqBOUmydXltZCg4RkSpWUurMSNvA/fOz+CZ/HxcMasfvT+tBy0axeXptZSk4RESq0NIN33LTnAyWbdxJaofDeOrswRzdtmm0y6pSCg4RkSqwc28Rf3kzi2c/+YpWjerx0AX9OaffkXF5nZiCQ0TkELg7c5du4o5XV/BNXiGXHNuR/z2tO02iePfaSFNwiIgcpLXb8rh5TgYfZG+jb9umPDl2UNztliqPgkNEpJIKikp49P3V/N97q6mXVIvbR/Xhf4Z0IKkGXPVdFRQcIiKV8MGqbdz8cgZrt+Vxdr8jufnMXhzeJLHu86bgEBEJw5bdBdz52gpeXrKJDi0a1PiL+A6FgkNE5ABKSp1/f7ae+974ksKiUn57cjeuPqkLyXWSol1a1Cg4RER+RMbGnfxxTgZLN3zLsK4tuGPUUdX66NZYpeAQESljd0ERD7y1kqc/WkfzhnXj+pqMg6HgEBEJcndez/ia217JZMvuQi4a0oGJp/egaf34vSbjYCg4RESA9dvzufnlDN5fuZXebZrw94tT6d+uWbTLikkRvUWjmZ1hZllmlm1mk8pZbmb2cHB5upkNDFl2nZllmlmGmT1nZsnB+bea2UYzWxL8GhnJdRCR+FZYXMLf3lnFqVPfJ23dN/zprN7MnTBMoXEAEes4zCwJeAQ4FcgBFpnZXHdfHjJsBNAt+DUEmAYMMbMU4DdAb3ffa2YvAhcATwVfN9Xdp0SqdhFJDB+v3s5Nc5axemseI48+gj+d1YcjmibWNRkHI5K7qgYD2e6+BsDMngdGAaHBMQp4xt0d+MTMmplZm5Da6ptZEdAA2BTBWkUkgWzbU8hd81Yw+/ONtGtenyfHDWJ4j8OjXVaNEcngSAE2hEznEOgqKhqT4u5pZjYFWA/sBd509zdDxk0wszFAGvB7d99R9sPNbDwwHqB9+/aHui4iEgdKS50X0jZwz+tfkr+vmAnDu3Lt8K7Ur5u412QcjEge4yjvvDUPZ4yZHUagG+kEHAk0NLOLgsunAV2A/kAu8JfyPtzdp7t7qruntmqVmFd3ish3VuTu4vxHP2Ly7GX0PKIxr//2eCae3kOhcRAi2XHkAO1Cptvyw91NPzbmFGCtu28FMLPZwHHAP9198/7BZvYY8GrVly4i8SKvsJgH317JEx+uo2n9OvxldD9+PjBF12QcgkgGxyKgm5l1AjYSOLj9qzJj5hLY7fQ8gd1YO90918zWA0PNrAGBXVUnE9gthZm1cffc4OvPBTIiuA4iUkO5O28u38ytczPJ3VnAhYPbccMZPWnWoG60S6vxIhYc7l5sZhOA+UAS8IS7Z5rZVcHljwLzgJFANpAPjAsu+9TMZgKfA8XAF8D04FvfZ2b9Cez2WgdcGal1EJGaKWdHPrfOzeTtFVvoeURj/varARzToXm0y4obFjihKb6lpqZ6WlpatMsQkQgrKinl8Q/W8tDbqzCD607pzthhHamTFNFL1uKWmS1299Sy83XluIjEhS/W72DSrGVkbd7Nab1bc8s5fUhpVj/aZcUlBYeI1Gh5hcXcPz+Lpz9exxFNknlsTCqn9m4d7bLimoJDRGqsd7O2cNNLGWzauZcxQwM3JGycrBsSRpqCQ0RqnO17Crnj1eXMWbKJroc3YuZVx+rgdzVScIhIjeHuzFmykdtfWc6ewmJ+e3I3rhnehXq1dRFfdVJwiEiNsOGbfP44J4OFK7cyoH0z7j2vL91bN452WQlJwSEiMa2k1Hnqo3VMmZ9FLYPbzunDRUM7kFRLV35Hi4JDRGLWitxdTJqVztKcnQzv0Yo/n3u0TrGNAQoOEYk5BUUl/O2dbB59fzVN69fh4QsHcHbfNrq/VIxQcIhITPl0zXYmz17Gmm15nDewLTed2YvDGur+UrFEwSEiMWFXQRH3vP4l//50Pe2a1+fZywZzfDc9EiEWKThEJOrmZ37Nn17OYOvuQq44vhPXndqdBnX16ylW6W9GRKJmy64CbpmbyesZX9OrTRMeG5NK37bNol2WVEDBISLVzt15YdEG7py3gsLiUv5wRg+uOL6z7mJbQyg4RKRard2Wx+TZ6Xyy5huGdGrO3T8/ms6tGkW7LKkEBYeIVIuiklIe+88aHnp7FXVr1+Kenx/NL1LbUUsX8tU4Cg4RibhlOTu5YVY6y3N3MeKoI7jtnD4c3iQ52mXJQVJwiEjE7N1XwtS3V/KP/6yhZaN6PHrRMZxx1BHRLksOkYJDRCLig1XbmPxSOhu+2cuFg9szaURPmtbXszLigYJDRKrUjrx93DlvBTMX59C5ZUNeGD+UIZ1bRLssqUIKDhGpEu7OK+m53P5KJt/mFzFheFcm/LQryXX0rIx4o+AQkUO26du93DwngwVfbqFf26Y8e9kQerVpEu2yJEIUHCJy0Nyd5z7bwF3zVlBS6tx0Zi/GDeukZ2XEOQWHiByU9dvzmTQ7nY9Wb+e4Li2497y+tGveINplSTVQcIhIpZSWOk9/vI773sgiqZZx98+P5oJB7fSsjASi4BCRsK3euocbZqaT9tUOTurRirvOPZoj9US+hKPgEJEKFZeU8vgHa3ngrZUk10nigV/049wBKeoyEpSCQ0QOKOvr3fxh5lKW5uzktN6t+fPPjtLtQhKcgkNEylVUUsq091bz13dW0Ti5Dn/71QDOPFrP/RYFh4iUI2PjTq6fmc6K3F2c0+9Ibjm7Ny0a1Yt2WRIjIvrUFDM7w8yyzCzbzCaVs9zM7OHg8nQzGxiy7DozyzSzDDN7zsySg/Obm9lbZrYq+OdhkVwHkURSWFzClPlZjHrkQ7btKWT6xcfw8IUDFBryPRELDjNLAh4BRgC9gQvNrHeZYSOAbsGv8cC04GtTgN8Aqe5+FJAEXBB8zSRggbt3AxYEp0XkEH2xfgdnPvwBf3s3m3MHpPD2dSdyWh/dyVZ+KJK7qgYD2e6+BsDMngdGActDxowCnnF3Bz4xs2Zm1iaktvpmVgQ0ADaFvOak4PdPA+8BN0RwPUTi2t59JTzwVhaPf7CW1k2SeWrcIE7qcXi0y5IYFsngSAE2hEznAEPCGJPi7mlmNgVYD+wF3nT3N4NjWrt7LoC755pZuT/hZjaeQBdD+/btD3VdROLSp2u2c8OsdNZtz+dXQ9ozeURPGifr1udyYJE8xlHeqRcezpjgcYtRQCfgSKChmV1UmQ939+nunuruqa1atarMS0XiXl5hMX96OYNfTv+EEnf+ffkQ7jr3aIWGhCWSHUcO0C5kui3f7W6qaMwpwFp33wpgZrOB44B/ApvNrE2w22gDbIlQ/SJx6YNV27hhVjqbdu5l7HEd+cMZPWhQVydYSvgi2XEsArqZWSczq0vg4PbcMmPmAmOCZ1cNBXYGd0OtB4aaWQMLnDR+MrAi5DWXBL+/BHg5gusgEjd2FRQxaVY6Fz3+KfVq1+LFK4/l1nP6KDSk0iL2E+PuxWY2AZhP4KyoJ9w908yuCi5/FJgHjASygXxgXHDZp2Y2E/gcKAa+AKYH3/oe4EUzu4xAwIyO1DqIxIt3vtzMjbMz2LK7gCtP7Mx1p3TXA5bkoFnghKb4lpqa6mlpadEuQ6Ta7cjbx+2vLuelLzbSvXUj7j+/H/3aNYt2WVJDmNlid08tO189qkiceiMjl5vmZPJt/j5+c3I3rh3ehXq11WXIoVNwiMSZbXsKueXlTF5blkufI5vw9KWD6HNk02iXJXFEwSESJ9yduUs3cevcTPIKS7j+9B6MP6EzdZIiemchSUAKDpE4sHlXAX98aRlvr9hC/3bNuP/8vnRr3TjaZUmcUnCI1GDuzozFOdzx6nL2FZdy05m9GDesE0m1dOtziZywgsPMhrn7hxXNE5Hqs+nbvUyevYz3V25lcMfm3Ht+Xzq1bBjtsiQBhNtx/BUYGMY8EYkwd+fFtA38+dUVFJc6t57dmzHHdqSWugypJgcMDjM7lsCtPlqZ2f+GLGpC4KI+EalGoV3GkE7Nue/8vnRooS5DqldFHUddoFFwXOiRtl3A+ZEqSkS+r2yXcds5fbh4aAd1GRIVBwwOd38feN/MnnL3r6qpJhEJ8b1jGZ2ac7+6DImycI9xPGVmP7g3ibv/tIrrEZEgHcuQWBVucEwM+T4ZOI/AzQdFJAJyd+5l0ix1GRKbwgoOd19cZtaHZvZ+BOoRSWjuzoy0wHUZ6jIkVoV7HUfzkMlawDGAnmIvUoXUZUhNEe6uqsUEHvtqBHZRrQUui1RRIolEXYbUNOHuquoU6UJEEpG6DKmJwt1VlQxcA/yEQOfxATDN3QsiWJtI3Aq9x1RxiboMqVnC3VX1DLCbwG1GAC4EnkWPbRWptNydgesy3stSlyE1U7jB0cPd+4VMv2tmSyNRkEi8Upch8SLc4PjCzIa6+ycAZjYE0J1xRcKkLkPiSbjBMQQYY2brg9PtgRVmtgxwd+8bkepEajh1GRKPwg2OMyJahUgc+npnAZNnp/OuugyJM+EGx5/d/eLQGWb2bNl5IvL9LqOopJRbzu7NJeoyJI6EGxx9QifMrDaBq8dFJMT3uoyOgedldNRT+STOVPQgp8nAjUB9M9tF4MpxgH3A9AjXJlJjqMuQRFLR8zjuBu42s7vdfXI11SRSo6jLkEQT7q6q183shLIz3X1hFdcjUmO4OzMX53C7ugxJMOEGx/Uh3ycDgwnc+FAPcpKEpC5DElm4Nzk8O3TazNoB90WkIpEY5u7M/nwjt76SqS5DEla4HUdZOcBRVVmISKzbsquAG19axtsrtpDa4TCmjO6nLkMSUrh3x/0rgbviQuBBTgOACu9VZWZnAA8BScA/3P2eMsstuHwkkA+MdffPzawH8ELI0M7An9z9QTO7FbgC2BpcdqO7zwtnPUQOhrszd+km/vRyJgVFJdx0Zi/GDetEkroMSVDhdhzLCfzyd2An8Jy7H/BeVWaWBDwCnEqgQ1lkZnPdfXnIsBFAt+DXEGAaMMTds4D+Ie+zEXgp5HVT3X1KmLWLHLStuwu5ac4y5mduZkD7ZkwZ3Y8urRpFuyyRqKroOo7awF3ApcB6AtdxtAOeMLPP3L3oAC8fDGS7+5rgez0PjCIQQvuNAp5xdwc+MbNmZtbG3XNDxpwMrHb3ryq5biKH5NX0Tdw8J4O8fSVMHtGTy4/vrC5DhMBupwO5H2gOdHL3ge4+gMBuo2ZARf/jTwE2hEznBOdVdswFwHNl5k0ws3Qze8LMDivvw81svJmlmVna1q1byxsiUq7tewq59l+fM+HfX9C+eQNe+/VPuPLELgoNkaCKguMs4Ap3371/hrvvAq4mcFziQMr7V+aVGWNmdYFzgBkhy6cBXQjsysoF/lLeh7v7dHdPdffUVq1aVVCqSMAbGbmcNnUhby7/mutP78Gsq4+jW+vG0S5LJKZUdIzDg7uRys4sMbMfzC8jh8Burf3aApsqOWYE8Lm7bw757P9+b2aPAa9WUIdIhXbk7eOWuZnMXbqJo1Ka8K/RQ+h5RJNolyUSkyrqOJab2ZiyM83sIuDLCl67COhmZp2CncMFwNwyY+YSeM6HmdlQYGeZ4xsXUmY3lZm1CZk8F8iooA6RA3pr+WZOnbqQecty+d9Tu/PSNcMUGiIHUFHHcS0w28wuJXCluAODgPoEfmn/KHcvNrMJwHwCZ2Q94e6ZZnZVcPmjwDwCu7yyCZyOO27/682sAYEzsq4s89b3mVn/YC3rylkuEpad+UXc9koms7/YSK82TXj60kH0ObJptMsSiXlWzp6oHw4y+ymBW6sbkOnuCyJdWFVKTU31tLS0aJchMeSdLzczefYytu3Zx7XDuzJheFfq1q6oARdJLGa22N1Ty84P95Yj7wDvVHlVItVsV0ERd7yynBmLc+jRujH/GDOIo9uqyxCpjIO95YhIjfP+yq1MmpXO5l0FXDu8C785uRv1aidFuyyRGkfBIXFvd0ERd81bwXOfbaDr4Y2Yfc0w+rdrFu2yRGosBYfEtQ+zt/GHmenk7tzLlSd25rpTupNcR12GyKFQcEhcyiss5u7XV/DPT9bTuWVDZlx1HMd0KPcmAyJSSQoOiTsfr97OH2YtJWfHXi7/SScmnt5DXYZIFVJwSNzI31fMfW9k8dRH6+jQogEvXnksgzo2j3ZZInFHwSFxYdG6b5g4Yylfbc9n7HEd+cMZPWhQVz/eIpGgf1lSoxUUlXD//Cye+HAtbQ+rz3NXDOXYLi2iXZZIXFNwSI21+KsdXD9jKWu25XHx0A5MGtGThvX0Iy0SafpXJjVOQVEJU99ayWP/WUObpvX51+VDGNa1ZbTLEkkYCg6pUZZs+JaJM5aSvWUPFw5uz40je9I4uU60yxJJKAoOqREKi0t46O1VPPr+alo3SebpSwdzYnc9oEskGhQcEvMyNu7k9y8uJWvzbkYf05abzupN0/rqMkSiRcEhMWtfcSl/ezebR97NpkXDujwxNpWf9mwd7bJEEp6CQ2LS8k27mDhjKctzd3HugBRuPbsPTRuoyxCJBQoOiSlFJaU8+t5qHn5nFU3r1+HvFx/D6X2OiHZZIhJCwSExY+Xm3fz+xaUs27iTs/sdyW3n9KF5w7rRLktEylBwSNQVl5Ty2H/WMvWtlTRKrs3//c9ARh7dJtpliciPUHBIVGVv2cPEGUtZsuFbzuhzBH8+9yhaNqoX7bJE5AAUHBIVJaXOEx+s5f43s2hQN4mHLxzA2X3bYGbRLk1EKqDgkGq3dlse189YStpXOzilV2vu+vlRHN44OdpliUiYFBxSbUpLnac/Xse9b3xJ3aRaPPCLfpw7IEVdhkgNo+CQarF+ez7Xz1zKp2u/4aQerbjn5305oqm6DJGaSMEhEVVa6vzrs/XcPW8FSWbcd15fRqe2VZchUoMpOCRicnbkc8OsdD7M3s7x3Vpyz3l9SWlWP9plicghUnBIlXN3Xli0gT+/tgJ3565zj+bCwe3UZYjECQWHVKncnXu5YdYyFq7cyrGdW3Df+X1p17xBtMsSkSqk4JAq4e7MXJzD7a8up7jEuX1UHy4a0oFatdRliMQbBYccsi27Cpg8exkLvtzC4I7NuX90Xzq0aBjtskQkQmpF8s3N7AwzyzKzbDObVM5yM7OHg8vTzWxgcH4PM1sS8rXLzH4XXNbczN4ys1XBPw+L5DrIj3N35nyxkVOnLuSD7G3cfFZvnh8/VKEhEuciFhxmlgQ8AowAegMXmlnvMsNGAN2CX+OBaQDunuXu/d29P3AMkA+8FHzNJGCBu3cDFgSnpZpt3V3IVf9czO9eWEKXVg15/bfHc9lPOmnXlEgCiOSuqsFAtruvATCz54FRwPKQMaOAZ9zdgU/MrJmZtXH33JAxJwOr3f2rkNecFPz+aeA94IaIrYX8wKvpm7h5TgZ5+0q4cWRPLvtJZ5IUGCIJI5LBkQJsCJnOAYaEMSYFCA2OC4DnQqZb7w8Wd881s8PL+3AzG0+gi6F9+/YHU7+U8U3ePm6ek8Fry3Lp17YpU0b3o1vrxtEuS0SqWSSDo7z/gnplxphZXeAcYHJlP9zdpwPTAVJTU8t+rlTSGxm53DQng517i7j+9B5ceUJnaidF9BCZiMSoSAZHDtAuZLotsKmSY0YAn7v75pB5m/fvzjKzNsCWKqxZytiRt49b5mYyd+km+hzZhH9ePoSeRzSJdlkiEkWR/C/jIqCbmXUKdg4XAHPLjJkLjAmeXTUU2Fnm+MaFfH831f7XXBL8/hLg5aovXQDezPyaU6cuZN6yXP731O7MuXaYQkNEItdxuHuxmU0A5gNJwBPunmlmVwWXPwrMA0YC2QTOnBq3//Vm1gA4FbiyzFvfA7xoZpcB64HRkVqHRPVt/j5ue2U5L32xkV5tmvDMpYPpfaQCQ0QCLHBCU3xLTU31tLS0aJdRIyxYsZnJs5fxTd4+rh3elWuHd6VubR3LEElEZrbY3VPLzteV4wLAzr1F3P7KcmZ9nkPPIxrzxNhBHJXSNNpliUgMUnAI72ZtYfKsZWzdU8ivf9qVX/+0m7oMEflRCo4EtqugiD+/upwX03Lodngjpo85hr5tm0W7LBGJcQqOBLVw5VZumJXO5l0FXHNSF357Sjfq1U6KdlkiUgMoOBLM7oIi7pq3guc+20CXVg2Zfc0w+rdrFu2yRKQGUXAkkA9WbeOGWenk7tzLlSd05rpTu5NcR12GiFSOgiMB7Cks5u55K/jXp+vp3LIhM646jmM66G70InJwFBxx7qPV2/jDzHQ2fruXK47vxO9P66EuQ0QOiYIjTuXvK+be17/k6Y+/omOLBsy48lhSOzaPdlkiEgcUHHHo0zXbuX5mOht25DNuWEf+cHpP6tdVlyEiVUPBEUfy9xVz3xtZPPXROto3b8DzVwxlSOcW0S5LROKMgiNOLFr3DdfPWMq67flccmwHbhjRkwZ19dcrIlVPv1lquIKiEqbMz+LxD9eS0qw+/75iCMd1aRntskQkjik4arDFX+3g+hlLWbMtj4uGtmfyiF40rKe/UhGJLP2WqYEKikqY+tZKHvvPGto0rc+/Lh/CsK7qMkSkeig4apgv1u9g4oylrN6ax4WD23PjyJ40Tq4T7bJEJIEoOGqIwuISHnx7FX9/fzVHNEnmmUsHc0L3VtEuS0QSkIKjBkjP+Zbfv7iUVVv28MvUdvzxrF40UZchIlGi4IhhhcUlPLxgFY++v4aWjery5LhBDO9xeLTLEpEEp+CIUek53zJxxlJWbt7D+ce05eazetO0vroMEYk+BUeMKSwu4aG3V/H3hWto1ageT44dxPCe6jJEJHYoOGLI0g2BLmPVlj2MPqYtN6nLEJEYpOCIAaFnTB3eOFnHMkQkpik4oiy0y/hFalv+eKa6DBGJbQqOKCkoKuGhBeoyRKTmUXBEwZIN33J9SJdx01m9dV2GiNQYCo5qVFAUOJYxfeFqWjdJ5qlxgzhJXYaI1DAKjmqyJHgsI1tXf4tIDafgiDB1GSISbxQcEfTF+h1cPzNdXYaIxJVakXxzMzvDzLLMLNvMJpWz3Mzs4eDydDMbGLKsmZnNNLMvzWyFmR0bnH+rmW00syXBr5GRXIeDUVBUwt2vr+C8aR+RV1jM05cO5t7z+yo0RCQuRKzjMLMk4BHgVCAHWGRmc919eciwEUC34NcQYFrwT4CHgDfc/Xwzqws0CHndVHefEqnaD0Xo8zIuGNSOG89UlyEi8SWSu6oGA9nuvgbAzJ4HRgGhwTEKeMbdHfgk2GW0AfKAE4CxAO6+D9gXwVoPWUFRCVPfXsljC9fQukkyT186mBP1vAwRiUORDI4UYEPIdA7fdRMHGpMCFANbgSfNrB+wGPitu+cFx00wszFAGvB7d99R9sPNbDwwHqB9+/aHvjYHoC5DRBJJJI9xWDnzPMwxtYGBwDR3H0CgA9l/jGQa0AXoD+QCfynvw919urununtqq1aR+Z9/QVEJd88LHMvYu6+EZy4dzD3n6ViGiMS3SHYcOUC7kOm2wKYwxziQ4+6fBufPJBgc7r55/2Azewx4tWrLDs/n63dw/X+f/d2OySPVZYhIYohkcCwCuplZJ2AjcAHwqzJj5hLY7fQ8gd1YO909F8DMNphZD3fPAk4meGzEzNrsHwOcC2REcB1+oKCohKlvreSx/6zRs79FJCFFLDjcvdjMJgDzgSTgCXfPNLOrgssfBeYBI4FsIB8YF/IWvwb+FTyjak3IsvvMrD+BrmQdcGWk1qGsz4PHMtYEu4wbR/aisboMEUkwFjihKb6lpqZ6WlraQb++oKiEB95ayT+CXcY95/VVlyEicc/MFrt7atn5unK8Aou/2sH1M/d3Ge25cWRPdRkiktAUHAfw1wWrmPr2Sto0rc+zlw3m+G7qMkREFBwH0KFlQ345SF2GiEgoBccBnNPvSM7pd2S0yxARiSkRvcmhiIjEHwWHiIhUioJDREQqRcEhIiKVouAQEZFKUXCIiEilKDhERKRSFBwiIlIpCXGTQzPbCnx1kC9vCWyrwnJqOm2P72hbfJ+2x/fFw/bo4O4/uNdSQgTHoTCztPLuDpmotD2+o23xfdoe3xfP20O7qkREpFIUHCIiUikKjopNj3YBMUbb4zvaFt+n7fF9cbs9dIxDREQqRR2HiIhUioJDREQqRcERZGZnmFmWmWWb2aRylpuZPRxcnm5mA6NRZ3UIY1v0NLOPzazQzCZGo8bqFMb2+J/gz0S6mX1kZv2iUWd1CWN7jApuiyVmlmZmP4lGndWhom0RMm6QmZWY2fnVWV/EuHvCfwFJwGqgM1AXWAr0LjNmJPA6YMBQ4NNo1x3FbXE4MAi4E5gY7ZpjYHscBxwW/H5EvP5sVGJ7NOK746d9gS+jXXe0tkXIuHeAecD50a67Kr7UcQQMBrLdfY277wOeB0aVGTMKeMYDPgGamVmb6i60GlS4Ldx9i7svAoqiUWA1C2d7fOTuO4KTnwBtq7nG6hTO9tjjwd+YQEMgXs/ACef3BsCvgVnAluosLpIUHAEpwIaQ6ZzgvMqOiQeJsp7hquz2uIxAZxqvwtoeZnaumX0JvAZcWk21VbcKt4WZpQDnAo9WY10Rp+AIsHLmlf1fUjhj4kGirGe4wt4eZjacQHDcENGKoius7eHuL7l7T+BnwB2RLipKwtkWDwI3uHtJ5MupPrWjXUCMyAHahUy3BTYdxJh4kCjrGa6wtoeZ9QX+AYxw9+3VVFs0VOrnw90XmlkXM2vp7jX9hn9lhbMtUoHnzQwCNz0caWbF7j6nWiqMEHUcAYuAbmbWyczqAhcAc8uMmQuMCZ5dNRTY6e651V1oNQhnWySSCreHmbUHZgMXu/vKKNRYncLZHl0t+JsyePZhXSAew7TCbeHundy9o7t3BGYC19T00AB1HAC4e7GZTQDmEzgD4gl3zzSzq4LLHyVwRsRIIBvIB8ZFq95ICmdbmNkRQBrQBCg1s98ROJtkV7TqjpQwfzb+BLQA/i/4+7LY4/SuqGFuj/MI/CerCNgL/DLkYHncCHNbxCXdckRERCpFu6pERKRSFBwiIlIpCg4REakUBYeIiFSKgkNERCpFwSEiIpWi4BARkUr5f6B3su4IG4aJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEDCAYAAAAvNJM9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAApjElEQVR4nO3dd3xUdb7G8c83HZJQA9KLEEBQaaFIBxtYQCwoFsSGIIi97N7r6q53d+9eWUWKoKioWFjsqCi6SG8SFJBOaBJpASmhpP/uH5ldY3aAATI5yeR5v17zInPOb2aeHGOenG7OOURERAoL8zqAiIiUTCoIERHxSwUhIiJ+qSBERMQvFYSIiPilghAREb/KVEGY2Q1mtsbM8sws6STjtpnZj2a2wsySC0xvaWaLffM+M7MKp/i8Vr7xa8xslZndWJTfj4hIMIVsQZhZDzN7o9Dk1cC1wLwA3qKnc66Vc65gkbwKPOmcuwD4GHjsFO9xDBjknGsB9AZGm1mlQPKLiHgtZAvCH+fcOufchrN4i6b8Wi7fANcBmFm4mT1nZst8awr3+j5vo3Nuk+/rncBeoNpZfL6ISLEpUwVxGhzwtZktN7MhBaavBvr6vr4BqOv7+i7gkHOuHdAOuMfMGhZ8QzNrD0QBm4OaXESkiER4HaComdlSIBqIA6qY2QrfrCecczMDfJvOzrmdZlYd+MbM1jvn5gF3AmPM7A/AdCDLN/4y4EIzu973vCKQCGz1ZaoJTAFud87lnd13KCJSPEKuIJxzHSB/HwQw2Dk3+AzeY6fv371m9jHQHpjnnFtPfhlgZk2AK30vMeB+fwXk25H9BfDfzrklp5tFRMQr2sRUiJnFmln8v74mvxBW+55X9/0bBvw3MNH3spnAMDOL9M1v4nufKPJ3Zr/lnHu/eL8TEZGzU6YKwsz6m1kqcBHwhZnN9E2vZWYzfMPOARaY2UrgO+AL59xXvnkDzWwjsB7YCUz2TX8VWAt8b2argZfJXzsbAHQDBvsOmV1hZq2C/o2KiBQB0+W+RUTEnzK1BiEiIoELqZ3UCQkJrkGDBl7HEBEpNZYvX77POef3/KyQKogGDRqQnJx86oEiIgKAmW0/0TxtYhIREb9UECIi4pcKQkRE/AqpfRD+ZGdnk5qaSkZGhtdRPBcTE0OdOnWIjIz0OoqIlAIhXxCpqanEx8fToEEDzMzrOJ5xzrF//35SU1Np2LDhqV8gImVeyG9iysjIoGrVqmW6HADMjKpVq2pNSkQCFvIFAZT5cvgXLQcROR1loiBOZc/hDI5kZKPLjoiI/KrMF0RuXh6/HM1iy76jbE47yqHjRVcU+/fvp1WrVrRq1YoaNWpQu3btfz/Pysr6zdjRo0dz7NixU75njx49dDKgiBSLkN9JfSrhYWE0PSeeA8eySDuSyfb9R4mOCKdafDSVykcSdhabZapWrcqKFSsAeOaZZ4iLi+PRRx/1O3b06NHceuutlC9f/ow/T0SkKJX5NQiAsDCjalw0Tc+Jp16V8phB6oFjbNidzr4jmeTmFd2mp1mzZtG6dWsuuOAC7rzzTjIzMxkzZgw7d+6kZ8+e9OzZE4Bhw4aRlJREixYtePrpp0/6nk899RQDBgz49/N58+ZxxRVXFFlmESmbytQaxB8/W8PanYcDGpub58jOzSM3z2EGkeFhRISHUXh9onmtCjx9dYuA3jMjI4PBgwcza9YsmjRpwqBBg5gwYQIPPvggzz//PLNnzyYhIQGAP//5z1SpUoXc3FwuvvhiVq1axYUXXuj3fZ999lmuvvpqli5dSocOHXjttde44447AsokInIiWoM4gfAwIyYynHJR4YSZkZWTx7GsHLJy8jjTXRS5ubk0bNiQJk2aAHD77bczb948v2OnTZtGmzZtaN26NWvWrGHt2rUnfe9bbrmFadOmcfjwYb799lv69u17ZiFFRHzK1BpEoH/p+3M8O5e09EwOHcsGg8rlIqkWH010ZHjA7xEbGxvQuK1btzJq1CiWLVtG5cqVGTx48CnPX7jqqqt45plnSExM5JprriE6OjrgXCIi/mgNIkDlIsOpV6U8TWvEUaV8FAePZ7NhTzrb9x/lWFZOQO+RkZHBtm3bSElJAWDKlCl0794dgPj4eNLT0wE4fPgwsbGxVKxYkT179vDll1+e8r3j4uJo3rw5f/jDH7jzzjvP8LsUEflVmVqDKApREeHUrlyO6hWi2X8kk/1Hsjh0PJu46Aiqx0cTGx1xwhPSYmJimDx5MjfccAM5OTm0a9eOoUOHAjBkyBD69OlDzZo1mT17Nq1bt6ZFixace+65dO7cOaBs1113HZs2baJ169ZF9v2KSNkVUvekTkpKcoXPEVi3bh3nnXde0D4zNy+P/Uez2JeeRU5eHuWjIqgWH02FmBMXRbDMmTOHUaNG8fnnn59wTLCXh4iULma23DmX5G+e1iDOUnhYGNXjY0iIjQ7KuRQiIl5RQRSRf51LUSU2ikPHs9mbnknqgWPsORxGtfhoKpePIjwsuEXRo0cPevToEdTPEJGyo0wUhHOu2Db3mBmVykdRsVwk6Rk5pKVnsvPgcfYezqBqXDRVY6OICPfm2IBQ2pwoIsEX8gURExPD/v37i/2S32ZGhXKRVCgXydHM/KLYcziDtPRMqsZGkRAXTWRE8RXFv+4HERMTU2yfKSKlW8gXRJ06dUhNTSUtLc3rKJCbx6GMHHZn5YJB+ahw4qMjim2N4l93lBMRCUTIF0RkZGSJu4PaT/uP8cr8zUxLTiU7N48+59dgWPfGXFCnotfRRET+LeQPcy3J0tIzeX3hVt5evJ30zBy6JiZwf69E2jes4nU0ESkjTnaYqwqiBDickc3bS7bz+oKt7DuSRZfGCTx0aSJt66soRCS4VBClxPGsXN5esp2Jczez/2gWXRMTeOjSJrSpV9nraCISok5WEEHdO2pmvc1sg5mlmNmTfuY3M7PFZpZpZo+ezmtDUbmocO7pdi7zn+jJk32asfrnQ1z70iIGT/6OFTsOeh1PRMqYoK1BmFk4sBG4FEgFlgEDnXNrC4ypDtQHrgEOOOdGBfpaf0r7GkRhRzNzeHPxNl6Zt4WDx7Lp1aw6D13SRDuzRaTIeLUG0R5Icc5tcc5lAVOBfgUHOOf2OueWAdmn+9qyIDY6gvt6NGb+4z159LImLN9+gKvHLeDuN5ex+udDXscTkRAXzIKoDewo8DzVN61IX2tmQ8ws2cySS8S5DkEQHxPJiF6JzH+iJw9f2oTvtv7CVWMXMOSt5IDvkCcicrqCWRD+TlsOdHtWwK91zr3inEtyziVVq1Yt4HClUYWYSEZenMj8J3rx4CWJLN6ynyvGzGfolOWs362iEJGiFcyCSAXqFnheB9hZDK8NeRXLRfLgJU1Y8HgvRvZqzIKUffQePZ/h73zPxj3pXscTkRARzIJYBiSaWUMziwJuAqYXw2vLjIrlI3n4sqYseKInI3o2Zs6GvVw+eh4j3v2elL0qChE5O0E9D8LMrgBGA+HA6865P5vZUADn3EQzqwEkAxWAPOAI0Nw5d9jfa0/1eaF2FNPpOnA0i0nzt/DGom0cz86lb8tajLw4kUbV4ryOJiIllE6UK2P2H8nklflbeGvRdjJzcunXqjYjL06kYUKs19FEpIRRQZRR+45k8sq8Lby1eBtZOXn0b12HkRc3pn5VFYWI5FNBlHFp6Zm8PHczU5ZsJyfPcW3r2tzfK5F6Vct7HU1EPKaCEAD2Hs5gwtzNvLP0J/LyHNe1qcOIXo2pW0VFIVJWqSDkN/YczmDCnM28+11+UdyQVJcRvRpTu1I5r6OJSDFTQYhfuw9l8NKcFKZ+twOHY0BSXYb3bEwtFYVImaGCkJPaefA442enMC15B4ZxU/v8NYrq8bp/tUioU0FIQFIPHGP87M28n7yDqIgw7u56LkO6nUtcdMjfmVakzFJByGnZtu8oz329gS9W7aJqbBQjL05kYPt6REUE9fYhIuIBz24YJKVTg4RYxt/chk+HdybxnDienr6GS1+Yy+erdhJKf1CIyMmpIOSEWtatxHv3dGTy4HbERIQz4t0f6Dd+IYs27/M6mogUAxWEnJSZ0bNZdWY80JVRN7RkX3omN09ayuDJ37Fuly4xLhLKtA9CTktGdi5vLtrG+NkppGfmcG3rOjx8WROdQyFSSmkntRS5g8eymDBnM5MXbQNgcKcGDO/RmIrlI70NJiKnRQUhQfPzweM8//VGPvohlfjoCIb3bMztnRoQExnudTQRCYCOYpKgqV2pHH8f0JIZI7vSpn5l/vrlenqNmsMHy1PJzQudPz5EyiIVhBSJ82pW4I072vPuPR2oFh/No++v5Mox85m9fq8OjRUppVQQUqQ6NUrgk+GdGXdza45n53LHG8sYOGkJK3cc9DqaiJwmFYQUOTPjqgtr8c1D3flj3xZs2nOEfuMXMvyd79m276jX8UQkQNpJLUGXnpHNpHlbmDR/K9m5edzcoR4jL04kIS7a62giZZ6OYpISYW96Bi/+cxNTl+0gJiKMId0acXfXhsTqYoAinlFBSImyOe0Iz321ga/W7CYhLpoHLknkpnZ1iQzXFk+R4qbDXKVEaVQtjom3teXDYZ1omFCepz5ZzeUvzOPLH3fpiCeREkQFIZ5pW78y0+69iFcHJREeZgx753v6v7SIpVv2ex1NRFBBiMfMjEuan8OXD3Tlb9ddwK5Dx7nxlSXc9cYyNu5J9zqeSJmmgpASISI8jBvb1WPOoz157PKmfLf1F3qPnsfTn67m4LEsr+OJlEkqCClRykWFM7xnY+Y+3pNbOtRnypLt9Bw1hylLtuvSHSLFTAUhJVKV2CieveZ8vhjZlSbnxPPUJ6u5auwC7Z8QKUYqCCnRzqtZgalDOjLu5tYcOpbFja8sYcS737Pz4HGvo4mEPBWElHj/unTHrEd68MDFiXyzdg+9/j6HMbM2kZGd63U8kZClgpBSo1xUOA9d2oR/Ptydnk2r8/w3G7nk+bl8tVrnT4gEgwpCSp26Vcoz4da2vHt3B2KjIhj69vfc+tpSHRYrUsSCWhBm1tvMNphZipk96We+mdkY3/xVZtamwLyHzGyNma02s/fMLCaYWaX06dQ4gS9GduGPfVuw+ufD9HlxPs9MX8OhY9leRxMJCUErCDMLB8YDfYDmwEAza15oWB8g0fcYAkzwvbY2MBJIcs6dD4QDNwUrq5ReEeFh3N6pAbMf7cFN7ery5uJt9Bg1m3eX/qTDYkXOUjDXINoDKc65Lc65LGAq0K/QmH7AWy7fEqCSmdX0zYsAyplZBFAe2BnErFLKVYmN4s/9L+Dz+7uQWD2e33/8I33HLWDZtl+8jiZSagWzIGoDOwo8T/VNO+UY59zPwCjgJ2AXcMg597W/DzGzIWaWbGbJaWlpRRZeSqcWtSryj3s7MnZga345msUNExcz8r0f2HVIh8WKnK5gFoT5mVZ4nd/vGDOrTP7aRUOgFhBrZrf6+xDn3CvOuSTnXFK1atXOKrCEBjPj6pa1mPVId+7v1Ziv1uym16i5jPtWh8WKnI5gFkQqULfA8zr852aiE425BNjqnEtzzmUDHwGdgphVQlD5qAgeuawpsx7uTrcmCYz6eiOXvjCXmWt267BYkQAEsyCWAYlm1tDMosjfyTy90JjpwCDf0Uwdyd+UtIv8TUsdzay8mRlwMbAuiFklhNWtUp6Xb0vi7bs6EBMRzr1TlnPba9+xSYfFipxU0ArCOZcDjABmkv/LfZpzbo2ZDTWzob5hM4AtQAowCbjP99qlwAfA98CPvpyvBCurlA1dEhOY8UBXnr66OStTD9L7xfn86bO1HDquw2JF/NEtR6VM2n8kk1Ffb2Tqsp+oUj6KRy9vyoCkuoSH+dstJhK6dMtRkUKqxkXz12sv4LMRXWiYEMvvPvqRfuMXkKzDYkX+TQUhZdr5tSvy/tCLePGmVuxLz+L6iYt5cOoP7D6U4XU0Ec+pIKTMMzP6tarNrEe6M7xnI2b8uJtef5/D+NkpOixWyjQVhIhPbHQEj13ejG8e7kbnxgk8N3MDl70wj1nr9ngdTcQTKgiRQupXjWXSoCTeurM9keHGXW8mM+zt5ew5rM1OUraoIEROoFuTanz5QDceu7wps9bv5ZK/z2XKku3k6SKAUkaoIEROIioijOE9GzPzwW5cUKciT32ymusnLmLDbp1kJ6FPBSESgIYJsbxzdwf+fkNLtu47ypVj5vPczPXaiS0hTQUhEiAz47q2dZj1SA/6tqrF+Nmb6T16HgtT9nkdTSQoVBAip6lKbBTPD2jF23d1wAG3vLqUR6at5JejWV5HEylSKgiRM9QlMYGZD3ZjeM9GfLriZy7++xw++j5VV4qVkKGCEDkLMZHhPHZ5Mz4f2YUGCbE8PG0lt762lG37jnodTeSsqSBEikCzGhX4cGgnnu3XglU7DnH56HmMn51Cdm6e19FEzpgKQqSIhIUZt13UgG8e7k7PptV5buYGrhqzgOXbD3gdTeSMqCBEiliNijFMvK0tkwYlcTgjm+snLuKpT1ZzOEP3nZDSRQUhEiSXNj+Hbx7uzu0XNeDtpdu59Pm5fLV6l3ZiS6mhghAJorjoCJ7p24KP7+tMldhohr79Pfe8tZydB497HU3klAIqCDPrHMg0EfGvVd1KTB/Rmd/1acaClDQufX4ukxduJVfXdZISLNA1iLEBThORE4gMD+Pe7o345qHutG1QhT9+tpZrX1rI2p2HvY4m4lfEyWaa2UVAJ6CamT1cYFYFIDyYwURCVd0q5XnzjnZMX7mTP322lqvHLeDuLg158JImlIvS/1ZScpxqDSIKiCO/SOILPA4D1wc3mkjoKngXu+vb1OHleVu4bPRc5m5M8zqayL9ZIEdUmFl959z2YshzVpKSklxycrLXMURO25It+/n9xz+yJe0o/VrV4qmrmpMQF+11LCkDzGy5cy7J77wAC2I28B8DnXO9zj5e0VFBSGmWmZPLS7M3M2HOZspFhfP7K5oxIKkuZuZ1NAlhRVEQbQs8jQGuA3Kcc48XTcSioYKQUJCyN53ff7Sa77b9QoeGVfjLtRfQqFqc17EkRJ11QZzgTec657qfVbIipoKQUJGX55iWvIO/zFhHRnYew3s2ZmiPc4mO0E5sKVonK4hAz4OoUuCRYGaXAzWKNKWI/FtYmHFT+3r885HuXH5+DV7450b6jl3I6p8PeR1NypBAz4NYDiT7/l0MPALcFaxQIpKvenwMYwe25vXBSRw4lsU14xcydtYmcnSVWCkGZ7yJqSTSJiYJZQePZfGHT9cwfeVOWtapyN8HtKJxde2bkLNTFJuYYszsYTP7yMw+NLOHzCymaGOKyMlUKh/FmIGtGX9zG3765RhXjpnPawu2kqfLdUiQBLqJ6S2gBfmX1xgHnAdMCVYoETmxKy+sycyHutGlcQLPfr6WgZOWsOOXY17HkhAU6GGuK51zLU81zWvaxCRliXOO95NT+dPna3HO8dRVzbmxnc6bkNNz1puYgB/MrGOBN+wALAzgg3ub2QYzSzGzJ/3MNzMb45u/yszaFJhXycw+MLP1ZrbOd10oEfExMwa0q8tXD3blwjqVePKjH7nrzWT2Hs7wOpqEiEALogOwyMy2mdk28o9k6m5mP5rZKn8vMLNwYDzQB2gODDSz5oWG9QESfY8hwIQC814EvnLONQNaAusCzCpSptSpXJ537u7A01c3Z2HKPi59YR7TV+70OpaEgJNezbWA3mfw3u2BFOfcFgAzmwr0A9YWGNMPeMvlb+da4ltrqAkcBboBgwGcc1lA1hlkECkTwsKMOzo3pFuTajwybSUj3/uBmWt28z/9zqdybJTX8aSUCnQN4n+cc9sLPgpOO8FragM7CjxP9U0LZMy5QBow2cx+MLNXzSzW34eY2RAzSzaz5LQ0XQlTyrZG1eL4YOhFPHZ5U75es5vLRs/j2/V7vI4lpVSgBdGi4BMziwDanmDsv4f5mVZ4j/iJxkQAbYAJzrnW5K9R/Mc+DADn3CvOuSTnXFK1atVOEUkk9EWEhzG8Z2M+Hd6FqrFR3PlGMk98sIr0jGyvo0kpc9KCMLPfmVk6cKGZHTazdN/zPcCnp3jvVKBuged1gMIbRk80JhVIdc4t9U3/gPzCEJEANa9VgU9HdOa+Ho14f/kOeo+ez+LN+72OJaXISQvCOfdX51w88JxzroJzLt73qOqc+90p3nsZkGhmDc0sCrgJmF5ozHRgkO9opo7AIefcLufcbmCHmTX1jbuY3+67EJEAREeE83jvZrw/tBNREWEMnLSEP362hozsXK+jSSkQ6E7qL82sW+GJzrl5J3qBcy7HzEYAM8m/Penrzrk1ZjbUN38iMAO4AkgBjgF3FHiL+4F3fOWypdA8ETkNbetX5ouRXfjbl+uZvHAbczem8fyAVrSqW8nraFKCBXqi3GcFnsaQf4TSct0wSKT0WZiyj8feX8me9Ezu69GI+3slEhUR6O5ICTVnfaKcc+7qAo9LgfPJ3w8hIqVM58YJfPVQN/q3rs3Yb1O4ZvxC1u8+7HUsKYHO9M+GVPJLQkRKoQoxkYy6oSWTBiWxNz2DvmMXMmHOZnJ14T8pIKB9EGY2ll8PUQ0DWgMrgxVKRIrHpc3PoW397vz3Jz/yt6/W8891e/j7DS1pkOD3tCMpYwJdg1gLbAQ2AEuAx51ztwYtlYgUmyqxUYy/uQ0v3tSKTXvS6fPifN5avE2XEZeTr0H4Toj7C3An8BP5J7bVBV43s++cczrzRiQEmBn9WtWmQ8OqPPHhKv7w6Rq+XrOH/7v+QmpVKud1PPHIqdYgngOqAA2dc218ZzWfC1QCRgU5m4gUsxoVY3jjjnb8pf8FfP/TAS4fPY8Pl6cSSneelMCdqiCuAu5xzqX/a4Jz7jAwjPzzF0QkxJgZN3eox1cPdOO8GhV45P2V3DtlOfuOZHodTYrZqQrCOT9/OjjncvnP6yqJSAipV7U87w3pyH9dcR5zNqZx2Qvz+Gr1Lq9jSTE6VUGsNbNBhSea2a3A+uBEEpGSIjzMuKfbuXxxfxdqVyrH0Le/5/EPVnIsK8fraFIMTnomtZnVBj4CjgPLyV9raAeUA/o7534ujpCB0pnUIsGTnZvHmFmbGDc7hcbV4hh3cxua1oj3OpacpTM+k9o597NzrgPwJ2Ab+Ucy/ck5176klYOIBFdkeBiPXNaUt+/qwIFj2fQdt4Cp3/2kHdghLKBrMZUWWoMQKR5p6Zk8PG0F8zft4+qWtfhL//OJj4n0OpacgbO+FpOISEHV4qN58472PHZ5U2b8uIurxi5g9c+HvI4lRUwFISJnJCzMGN6zMVOHdCQrJ49rX1rE5IVbtckphKggROSstGtQhRkju9I1MYE/fraWe6cs5+CxLK9jSRFQQYjIWascG8Wrtyfx1FXNmb1hL1eOWcDy7Qe8jiVnSQUhIkXCzLirS0M+GNqJsDAY8PJiJs7drIv+lWIqCBEpUi3rVuKLkV3p3aIG//vlega/sUyX6SilVBAiUuQqxEQy7ubW/Ln/+SzZsp8rXpzP4s37vY4lp0kFISJBYWbc0qE+n9zXmbiYCG55dQmj/7lRd60rRVQQIhJUzWtV4LMRXbimdW1G/3MTt7y6hD2HM7yOJQFQQYhI0MVGR/D8gFaMuqElK3ccos+L85mzYa/XseQUVBAiUmyub1uHz+7vTPX4aAZPXsb/frme7Nw8r2PJCaggRKRYNa4ezyfDO3Nzh3pMnLuZG19eTOqBY17HEj9UECJS7GIiw/lL/wsYd3NrNu05whUvzmfmmt1ex5JCVBAi4pmrLqzF5yO7UL9qLPdOWc4z09eQmZPrdSzxUUGIiKfqV43lg2EXcWfnhryxaBvXTVjEtn1HvY4lqCBEpASIjgjnD1c3Z9KgJHb8cpyrxi7g0xW6J5nXVBAiUmJc2vwcZjzQlaY14nlg6gqe/HAVx7O0yckrKggRKVFqVyrH1CEdua9HI6Yu20G/8QvYtCfd61hlkgpCREqcyPAwHu/djLfubM8vR7O4etwCpi3boZsRFbOgFoSZ9TazDWaWYmZP+plvZjbGN3+VmbUpND/czH4ws8+DmVNESqZuTaoxY2RX2tSrzOMfruKhf6zgSGaO17HKjKAVhJmFA+OBPkBzYKCZNS80rA+Q6HsMASYUmv8AsC5YGUWk5KteIYYpd3Xg4UubMH3lTq4eu4A1O3X/6+IQzDWI9kCKc26Lcy4LmAr0KzSmH/CWy7cEqGRmNQHMrA5wJfBqEDOKSCkQHmaMvDiR9+7pyLGsHPqPX8S0ZTu8jhXyglkQtYGC/wVTfdMCHTMaeBw46YVazGyImSWbWXJaWtpZBRaRkq3DuVWZMbIr7RtW4fEPV/H0p6t1LacgCmZBmJ9phfcw+R1jZlcBe51zy0/1Ic65V5xzSc65pGrVqp1JThEpRarGRfPGHe24p2tD3ly8nVtfXcp+3bEuKIJZEKlA3QLP6wA7AxzTGehrZtvI3zTVy8zeDl5UESlNIsLD+K8rm/PCjS1ZseMgfcctZPXP2i9R1IJZEMuARDNraGZRwE3A9EJjpgODfEczdQQOOed2Oed+55yr45xr4Hvdt865W4OYVURKof6t6/DB0E7kOcf1Exfp7OsiFrSCcM7lACOAmeQfiTTNObfGzIaa2VDfsBnAFiAFmATcF6w8IhKaLqhTkekjunBh7Uo8MHUFf52xTrc1LSIWSieeJCUlueTkZK9jiIgHsnLyePbztUxZsp2uiQmMHdiaSuWjvI5V4pnZcudckr95OpNaREJCVEQYz15zPv977QUs2bKffuMXsmG3LtFxNlQQIhJSbmpfj6lDOnIsK5f+Ly3kq9W6EdGZUkGISMhpW78Kn43oQuI58Qx9eznPf7ORPO2XOG0qCBEJSTUqxvCPIR25vm0dxszaxJApy0nPyPY6VqmighCRkBUTGc5z11/I01c3Z/aGvfR/aRFb0o54HavUUEGISEgzM+7o3JApd7Vn/5FM+o1fyOz1e72OVSqoIESkTOjUKIHpI7pQt3J57nxzGS/NSdH9JU5BBSEiZUbdKuX5cFgnrrygJv/31QZGvPcDx7J0f4kTUUGISJlSLiqcsQNb80TvZsz4cRfXTVjMjl+OeR2rRFJBiEiZY2YM69GIyYPbkXrgGH3HLWBRyj6vY5U4KggRKbN6NK3O9BFdqBoXzW2vf8frC7Zqv0QBKggRKdMaJsTy8X2d6NWsOn/6fC2Pvr+KjOxcr2OVCCoIESnz4mMiefnWtjx4SSIffp/KjS8vZteh417H8pwKQkQECAszHrykCS/f1paUvUe4euxCkrf94nUsT6kgREQKuLxFDT4e3pm46HAGTlrCu0t/8jqSZ1QQIiKFNDknnk+Hd+GiRgn8/uMf+a+PfyQrJ8/rWMVOBSEi4kfF8pFMHtyOe7ufyztLf+KWV5eQlp7pdaxipYIQETmB8DDjd33O48WbWvHjz4foO24Bq1IPeh2r2KggREROoV+r2nwwtBNhZlw/cTEffZ/qdaRioYIQEQnA+bUrMn1EZ9rUq8TD01by7OdryckN7f0SKggRkQBVjYtmyl0dGNypAa8t2Mrtk7/j0PHQvQmRCkJE5DREhofxTN8W/N/1F/Ld1l+4YeIidh4MzZPqVBAiImdgQFJd3rijPbsOZtD/pYWs3XnY60hFTgUhInKGOjdOYNrQizCMAS8vZsGm0LoirApCROQsnFezAh8P70TtSuUYPPm7kDrCSQUhInKWalYsx/vDLqJdgyo8PG0l42eHxu1MVRAiIkWgQkwkb97Znmta1eK5mRv4r09Wl/rDYCO8DiAiEiqiIsJ4fkAralYqx4Q5m9lzKIOxN7emfFTp/FWrNQgRkSIUFmY80bsZz15zPrM37OWmV0rvNZxUECIiQXBbx/q8fFsSG/ekc92ERWxJO+J1pNOmghARCZJLm5/De/d05EhmDtdNWMTy7Qe8jnRagloQZtbbzDaYWYqZPelnvpnZGN/8VWbWxje9rpnNNrN1ZrbGzB4IZk4RkWBpXa8yHw3rRIVykdw8aQkz1+z2OlLAglYQZhYOjAf6AM2BgWbWvNCwPkCi7zEEmOCbngM84pw7D+gIDPfzWhGRUqFBQiwfDevEeTUrMPTt5by5aJvXkQISzDWI9kCKc26Lcy4LmAr0KzSmH/CWy7cEqGRmNZ1zu5xz3wM459KBdUDtIGYVEQmqqnHRvHdPRy5udg5PT1/DX2esIy+vZJ8rEcyCqA3sKPA8lf/8JX/KMWbWAGgNLPX3IWY2xMySzSw5LS3tbDOLiARNuahwXr6tLbd2rMfL87bwwD9WkJmT63WsEwpmQZifaYXr8qRjzCwO+BB40Dnn90pYzrlXnHNJzrmkatWqnXFYEZHiEB5mPNvvfJ7o3YzPVu5k0Gsl95LhwSyIVKBuged1gJ2BjjGzSPLL4R3n3EdBzCkiUqzMjGE9GvHiTa34/qcDXD9hET+XwEuGB7MglgGJZtbQzKKAm4DphcZMBwb5jmbqCBxyzu0yMwNeA9Y5554PYkYREc/0a1WbN+9sz+5DGfQfv5A1Ow95Hek3glYQzrkcYAQwk/ydzNOcc2vMbKiZDfUNmwFsAVKAScB9vumdgduAXma2wve4IlhZRUS80qlRAu8Pu4jwMOPGl5cwf1PJ2ZdqoXDFwX9JSkpyycnJXscQETltuw4d547Jy0jZe4S/XXch17WtUyyfa2bLnXNJ/ubpTGoRkRKgZsVyTBt6ER3OrcIj769k3LebPL9kuApCRKSEqBATyeTB7enfujajvt7I7z/+0dNLhpfOa9CKiISo/EuGt6RWpRjGz97M7kMZjLu5DbHRxf/rWmsQIiIljJnx2OXN+J9rzmfuxjQGTvLmkuEqCBGREurWjvV55bYkNu05wrUTFhb7JcNVECIiJdglzc/hvSEdOZaZ67tk+C/F9tkqCBGREq5V3Up8dF8nKpaL5OZJS/lqdfFcMlwFISJSCtSvGsuHwzrRvFYFhr2znDcWbg36Z6ogRERKiapx0bx7d0cuOe8cnvlsLX8J8iXDVRAiIqVIuahwJt7alkEX1eeVeVsYOfWHoF0yXOdBiIiUMuFhxh/7tqB2pXL89cv17E3P5I072lE+qmh/pasgRERKITPj3u6NqFExhoUp+ygXGV7kn6GCEBEpxfq1qk2/VsG5I7P2QYiIiF8qCBER8UsFISIifqkgRETELxWEiIj4pYIQERG/VBAiIuKXCkJERPwyr2+KXZTMLA3YfoYvTwD2FWGc0kzL4re0PH5Ly+NXobAs6jvnqvmbEVIFcTbMLNk5l+R1jpJAy+K3tDx+S8vjV6G+LLSJSURE/FJBiIiIXyqIX73idYASRMvit7Q8fkvL41chvSy0D0JERPzSGoSIiPilghAREb/KVEGYWW8z22BmKWb2pJ/5ZmZjfPNXmVkbL3IWlwCWRzMzW2xmmWb2qBcZi1MAy+MW38/FKjNbZGYtvchZHAJYFv18y2GFmSWbWRcvchaXUy2PAuPamVmumV1fnPmCxjlXJh5AOLAZOBeIAlYCzQuNuQL4EjCgI7DU69weL4/qQDvgz8CjXmcuAcujE1DZ93WfUP35CHBZxPHrPswLgfVe5/ZyeRQY9y0wA7je69xF8ShLaxDtgRTn3BbnXBYwFehXaEw/4C2XbwlQycxqFnfQYnLK5eGc2+ucWwZkexGwmAWyPBY55w74ni4B6hRzxuISyLI44ny/FYFYIJSPdgnkdwfA/cCHwN7iDBdMZakgagM7CjxP9U073TGhoix9r4E43eVxF/lrm6EooGVhZv3NbD3wBXBnMWXzwimXh5nVBvoDE4sxV9CVpYIwP9MK/9UTyJhQUZa+10AEvDzMrCf5BfFEUBN5J6Bl4Zz72DnXDLgGeDbYoTwUyPIYDTzhnMsNfpziE+F1gGKUCtQt8LwOsPMMxoSKsvS9BiKg5WFmFwKvAn2cc/uLKVtxO62fDefcPDNrZGYJzrnSfuE6fwJZHknAVDOD/Av4XWFmOc65T4olYZCUpTWIZUCimTU0syjgJmB6oTHTgUG+o5k6Aoecc7uKO2gxCWR5lCWnXB5mVg/4CLjNObfRg4zFJZBl0dh8vw19R/tFAaFamKdcHs65hs65Bs65BsAHwH2lvRygDK1BOOdyzGwEMJP8ow1ed86tMbOhvvkTyT/64AogBTgG3OFV3mALZHmYWQ0gGagA5JnZg+QfvXHYq9zBEuDPxx+AqsBLvt+NOS4Er+QZ4LK4jvw/prKB48CNBXZah5QAl0dI0qU2RETEr7K0iUlERE6DCkJERPxSQYiIiF8qCBER8UsFISIifqkgRETELxWEiIj49f9cvAFD1vYA1wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEDCAYAAAAvNJM9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAttklEQVR4nO3dd3hUZfr/8fc9KYTQq2AACRKkN0MHEbEAotgFK6wsgl3XXfW7RXdX97euDREFFZHVdRexoNgLNpAapEgRjYASQKrUEEKS+/dHxjVmAxkgk0kmn9d1zSXnnOfM3HPMlU/Oec55HnN3RERECgtEugARESmbFBAiIlIkBYSIiBRJASEiIkVSQIiISJEUECIiUqQKFRBm9oCZfWVmy8xsupnVLKLNSWa2pMBrt5ndUmD7jWa22sxWmNk/guvqmNnHZrbXzMaHWMsLwfdZbmaTzSyupL6niEhJiNqAMLNTzWxKodUfAG3dvT3wNXBX4f3cfbW7d3T3jsDJQCYwPfie/YAhQHt3bwM8GNwtC/gjcPsRlPgC0BJoB1QGRh7BviIiYRe1AVEUd3/f3XOCi/OARsXs0h/41t2/Cy6PAf7u7geC77cl+N997j6b/KD4BTM708zmmtkXZvaSmVUN7vO2BwELQqhFRKRUVaiAKORXwDvFtBkK/KfAcgugj5nNN7NPzazL4XY2s7rAH4DT3b0zkAbcVqhNHHAl8O4R1i8iElaxkS6gpJnZfKASUBWobWZLgpvucPf3gm1+D+SQf5nnUO8TD5zLLy9DxQK1gO5AF2CamTXzQ49X0h1oDXxuZgDxwNxCbZ4APnP3WaF+RxGR0hB1AeHu3SC/DwIY7u7DC243s6uBwUD/w/xiBxgIfOHumwusywBe/emykJnlAXWBrYd4DwM+cPdhRW40uxuoB1xbzNcSESl1FeoSk5kNAO4AznX3zGKaD+OXl5cAXgNOC75XC/LPCLYd5j3mAb3MrHlwn8TgfpjZSOAsYJi75x3hVxERCTuL1tFcizqDMLN08i8/bQ+umufuo83seGCSuw8KtksE1gPN3H1Xgf3jgclARyAbuN3dPwpuWwdUJz80dgJnuvtKMzsNuD/4uQB/cPcZZpYDfAfsCa5/1d3/UoKHQETkmERtQIiIyLGpUJeYREQkdFHVSV23bl1v2rRppMsQESk3Fi1atM3d6xW1LaoComnTpqSlpUW6DBGRcsPMvjvUNl1iEhGRIikgRESkSAoIEREpUlT1QRTl4MGDZGRkkJX1P+PoyRFISEigUaNGxMVpVHKRiiLqAyIjI4Nq1arRtGlTguMhyRFyd7Zv305GRgbJycmRLkdESknUX2LKysqiTp06CodjYGbUqVNHZ2EiFUzUBwSgcCgBOoYiFU+FCIjibN6dxe6sg2jYERGRn1X4gMjNc7bvzWbdtn2s/mEPW3ZncTC3ZAZX3b59Ox07dqRjx440aNCApKSk/y5nZ2f/ou3YsWPJzCxugFk49dRT9TCgiJSKqO+kLk5MwGjZsBq79x9kx75sftidxebdB6heOZbaVeKpWin2qC+v1KlThyVLlgBwzz33ULVqVW6/vehpq8eOHcsVV1xBYmLi0X6VYt1zzz00bdqU4cOHh+0zRCR6VPgzCICAGTUT42lWryotjqtGnarx7D2Qw9pt+1i9eQ9b9pTcWcXMmTPp1KkT7dq141e/+hUHDhxg3LhxbNy4kX79+tGvXz8AxowZQ2pqKm3atOHuu+8u9n3ffvttOnToQMuWLZkxYwbdunUjJyen2P1GjBjBbbf9PAvqc889x3XXXXf0X1BEokaFOoP48xsrWLlxd8jtc/KcnNw8cvMcDGIDRmwgQEzg5zOK1sdX5+5z2oT0fllZWQwfPpyZM2fSokULrrrqKiZMmMAtt9zCww8/zMcff0zdunUBuO+++6hduza5ubn079+fZcuW0b59+yLfNzMzk5EjR7Jw4UIAunXrRp8+fYiNLf5/7+TJk+nUqRObNm2iYcOGPPPMMzz44IMhfR8RiW46gziM2ICREBdD5fgY4gIBcvOcrIO57M/O5WBuHkfapZ2bm0tycjItWrQA4Oqrr+azzz4rsu20adPo3LkznTp1YsWKFaxcufKQ77t69WqaNWtGUlISSUlJ1K1bl3POOQeAL7/88r/9HhMnTuRPf/rTf5e3b9+OmTF06FBefvll0tPT2bFjB126dDnCbyYi0ahCnUGE+pf+oeTlOTuDfRWZ2TmYGTUrx7HvQA6J8THF9lVUqVIlpM9Zu3YtDz74IAsXLqRWrVoMHz682GcQKlWq9N9/V6tWjU6dOgHQrl27X/SDFNUHcfHFFzNixAg2btyo/gkR+S+dQRyBQMCoXSWe5vWrklK/GrUT49m9/yDfbt3LN5v3sm3vAXIO01eRlZXFunXrSE9PB+D555+nb9++QP4v9T178mcf3b17N1WqVKFGjRps3ryZd95557B1tWjRgjVr1rBv3z42bNjAggULWL58ecjf68QTTyQzM5MpU6Zw5ZVXhryfiES3CnUGUZIqx8eQFF+ZBjUS2LU/mx37stm4cz8/7MqiRuU4aleJ/5+zioSEBJ599lkuvvhicnJy6NKlC6NHjwZg1KhRDBw4kIYNG/Lxxx/TqVMn2rRpQ7NmzejVq9dha6lSpQp/+9vfGDx4MJmZmbz44ovcf//9NG3aNOTLReeffz7z58+nfv36R39QRCSqRNWc1KmpqV74GYFVq1bRqlWrUvn8/dk5bN+Xzc7Mg+S5kxAXQ50q8dRMjCMmULZP1qZMmUJaWhrjx48/ZJvSPJYiUjrMbJG7pxa1Lay/tcxsgJmtNrN0M7uziO1mZuOC25eZWefg+pPMbEmB124zuyWctZaEyvGxNKqVSKuG1UmqWRkDNuzcz6pNe8j4MZPM7OJvOxURKSvCdonJzGKAx4EzgAxgoZnNcPeCt+MMBFKCr27ABKCbu68GOhZ4nw3A9HDVWtJiAkadqpWoXSWe/Qdz2bE3/6xix75sKsfFULtqPDUrx//idtlIGz58uDqoReQXwtkH0RVId/c1AGY2FRgCFAyIIcBznn+da56Z1TSzhu6+qUCb/sC37n7IeVOL4+4RGWzOzEiMjyWxdiwN8/L4MRgSG37czw87s6iZGEftKpWoHB9T6rUdqWi6FCkioQlnQCQB6wssZ5B/llBcmySgYEAMBf5zqA8xs1HAKIAmTZr8z/aEhAS2b98e8SG/YwIB6latRJ0q8WRm57JjXzY/Zh5k+75sEuPzh/WoWTmOQBk6q/jJT/NBJCQkRLoUESlF4QyIon7TFf4z9LBtzCweOBe461Af4u5PAU9Bfid14e2NGjUiIyODrVu3hlJzqbI850B2Ljuyc1ib6wQMqlSKpUqlWGLLWFD8NKOciFQc4QyIDKBxgeVGwMYjbDMQ+MLdNx9tEXFxcWV+FjR3Z8HaHTw/7zveWZ4BwNntGjKyTzLtG9WMbHEiUmGFMyAWAilmlkx+J/NQ4LJCbWYANwT7J7oBuwr1PwzjMJeXooWZ0a1ZHbo1q0PGj5n8c846pi5Yz4ylG+maXJuRvZM5vdVxZfLyk4hEr7A+B2Fmg4CxQAww2d3vM7PRAO4+0fI7BcYDA4BMYIS7pwX3TSS/f6KZu+8K5fOKeg6ivNqTdZAXF67n2c/XsWHnfprWSeSa3slceHIjEuP1fKOIlIzDPQcR9Q/KlXc5uXm8u+IHnp61lqXrd1IzMY7LuzXh6h5NqV9dncYicmwUEFHA3Vn03Y88PWsN76/cTGzAOLdDEiP7JNOqYfVIlyci5dThAkLXKsoJMyO1aW1Sm9bmu+37mDx7LdPSMnjliwx6N6/LyD7J9G1RL6K38opIdNEZRDm2K/MgLyz4jn/OWcfm3QdIqV+VkX2SGdIxiYS4sv/wnYhEni4xRbnsnDzeXLaRp2etZdWm3dStGs+V3ZtyRfcm1Klaqfg3EJEKSwFRQbg7c7/dztOz1vDx6q1Uig1wQedGXNM7meb1q0a6PBEpg9QHUUGYGT2b16Vn87qkb9nDM7PX8eoXGfxnwfec1rI+I/sk06NZZIccEZHyQ2cQUW773gP8a973PDd3Hdv3ZdO6YXVG9klmcPvjiY8t23NUiEj46RKTkHUwl9cWb2DS7LWkb9nLcdUrMbxnMpd1bUKNxLhIlyciEaKAkP/Ky3M+/WYrz8xay+z0bSTGx3BJamNG9GrKCXWqRLo8ESllCggp0qpNu5k0ay0zlm4gJ885q3UDRvZJ5uQTaqmfQqSCUEDIYW3encVzc9fxr3nfs2v/QTo2rsmNpzXntJb1FRQiUU4BISHJzM7hlUUZPDVrDet37KddUg1u7p9C/1YKCpFopYCQI3IwN4/pizcw/qN0vt+RSduk6tx0WgpntD5OQSESZRQQclQO5ubx2uINjP84ne+2Z9K6YXVuPj2FMxUUIlFDASHHJCc3j9eXbOSxj75h3fZMWjWszs39m3Nm6waaxEiknFNASInIyc1jxtKNjP8onTXb9tGyQTVu7p/CWW0UFCLllQJCSlRunvPG0o2M++gb1mzdx0nHVeOm/ikMbKugEClvFBASFrl5zpvLNjJu5jd8u3UfLY6ryk39UxjUtqGCQqScOFxAhHUwHjMbYGarzSzdzO4sYruZ2bjg9mVm1rnAtppm9rKZfWVmq8ysRzhrlSMXEzCGdEzi/Vv7Mm5YJ/Icbvj3Ys4a+xkzlm4kNy96/vgQqYjCFhBmFgM8DgwEWgPDzKx1oWYDgZTgaxQwocC2R4F33b0l0AFYFa5a5djEBIxzOxzPe7ecwmPDOgFw03/yg+L1JRsUFCLlVDjPILoC6e6+xt2zganAkEJthgDPeb55QE0za2hm1YFTgGcA3D3b3XeGsVYpATEB45xgUIy/rBMBg5unLuHMRz7ltcUKCpHyJpwBkQSsL7CcEVwXSptmwFbgWTNbbGaTzKzIkeTMbJSZpZlZ2tatW0uuejlqgYAxuP3xvHvzKTxxeWdiAwFueXEJZzzyKdMXZ5CTmxfpEkUkBOEMiKJ6KQv/CXmoNrFAZ2CCu3cC9gH/04cB4O5PuXuqu6fWq1fvWOqVEhYIGIPaNeSdm/sw4fLOxMcEuPXFpZzxyGe8+oWCQqSsC2dAZACNCyw3AjaG2CYDyHD3+cH1L5MfGFIOBQLGwHYNefumPky84mQS4mK4bdpSTn/4U15epKAQKavCGRALgRQzSzazeGAoMKNQmxnAVcG7mboDu9x9k7v/AKw3s5OC7foDK8NYq5SCQMAY0LYBb93YmyevPJnE+Fhuf2kp/R/+lJfS1isoRMqYsD4HYWaDgLFADDDZ3e8zs9EA7j7R8gf0GQ8MADKBEe6eFty3IzAJiAfWBLf9eLjP03MQ5Yu78+GqLYz98GtWbNxNk9qJ3NCvOed3TiIuRtOhipQGPSgnZZq7M3PVFsbO/JrlG3bTuHZlbujXnAs6N1JQiISZAkLKBXfno6+28OjMb1iWsYvGtSvz27Nack77hho9ViRMIvYktciRMDP6tzqO16/vxeThqVSJj+Wm/yzmvCfmsGDtjkiXJ1LhKCCkzDEzTmt5HG/d1IcHLmrP5l1ZXPLkXK59Po01W/dGujyRCkOXmKTM25+dyzOz1zDhk285kJPH5d2acFP/FOpUrRTp0kTKPfVBSFTYuucAYz/8mqkL15MYF8N1/ZozoldTEuJiIl2aSLmlPgiJCvWqVeK+89vx3i196NasNve/+xX9H8of5ylP4zyJlDgFhJQ7zetXY9LVXfj3r7tRq0oct7y4hHMfn82cb7dFujSRqKKAkHKr54l1mXF9bx65tAM79mZz2dPzGfnPhaRv2RPp0kSiggJCyrVAwDi/UyM+uv1U7hjQkvlrdnDW2Fn8fvqXbN1zINLliZRr6qSWqLJ97wEe+yidf837jkqxAcaceiLX9G5G5Xh1ZIsURZ3UUmHUqVqJe85tw/u3nkLvlLo8+P7X9HvwE15KW68Ji0SOkAJColKzelV58spUXhrdg+NqJPDbl5cx+LHZzP5GHdkioVJASFTr0rQ2r13Xk8eGdWJP1kGueGY+V09ewOof1JEtUhwFhEQ9s/y5smf+pi+/H9SKxd//yMBHP+POV5axZXdWpMsTKbPUSS0Vzs7MbB77KJ3n5q4jNhDg2r7N+HWfZlSpFBvp0kRKnTqpRQqomRjPHwe35sPb+nJaq/qM/fAbTn3wE6Yu+F4d2SIFKCCkwjqhThUev6wzr4zpSZPaidz56pcMenQWH6/eQjSdWYscLQWEVHgnn1CLl0f3YMLlncnKyWXEswu58pkFrNi4K9KliURUWAPCzAaY2WozSzezO4vYbmY2Lrh9mZl1LrBtnZl9aWZLzEwdCxJWZsbAdg354Na+3H1Oa5Zv3MXgx2Zz+0tL2bRrf6TLE4mIsAWEmcUAjwMDgdbAMDNrXajZQCAl+BoFTCi0vZ+7dzxUB4pISYuPDTCiVzKf/rYfo/o0Y8aSjfR/6FMmzVpDTm5epMsTKVXhPIPoCqS7+xp3zwamAkMKtRkCPOf55gE1zaxhGGsSCUmNynHcNagVM3/Tlx7N6nDvW6s4Z/znfPH9j5EuTaTUhDMgkoD1BZYzgutCbePA+2a2yMxGHepDzGyUmaWZWdrWrVtLoGyRnzWuncikq1OZeMXJ7MzM5sIJc/i/6V+yK/NgpEsTCbtwBoQVsa7wrSGHa9PL3TuTfxnqejM7pagPcfen3D3V3VPr1at39NWKHIKZMaBtAz64rS/X9ErmxYXr6f/wJ0xfnKG7nSSqhTMgMoDGBZYbARtDbePuP/13CzCd/EtWIhFTtVIsfxjcmhk39CKpViK3vriUyyfN59uteyNdmkhYhDMgFgIpZpZsZvHAUGBGoTYzgKuCdzN1B3a5+yYzq2Jm1QDMrApwJrA8jLWKhKzN8TV4dUxP7j2vLV9u2MXAsbN4+IOvyTqYG+nSREpU2MYWcPccM7sBeA+IASa7+wozGx3cPhF4GxgEpAOZwIjg7scB083spxr/7e7vhqtWkSMVEzCu6H4CZ7Y5jr+9tYpxM79hxpIN/PW8tvRJ0aVOiQ4ai0mkBHyevo0/vLactdv2cW6H4/nD4FbUr5YQ6bJEiqWxmETCrFfzurxzcx9uOT2Fd5f/QP+HPuX5ues0tpOUawoIkRKSEBfDLae34N1b+tC+UQ3++PoKLpgwh+UbNGSHlE8KCJES1qxeVf51TTceHdqRDT9mcu742fzljZXsPZAT6dJEjogCQiQMzIwhHZOYedupDOvahGfnrOX0hz7lnS836dkJKTcUECJhVCMxjvvOb8crY3pSq0o8Y174gl9NWcj6HZmRLk2kWAoIkVLQuUkt3rihF384uxXz1+7gjEc+5YlP0snO0QCAUnYpIERKSWxMgJF9mvHhbX05tUV9/vHuas4eN4sFa3dEujSRIikgRErZ8TUrM/HKk3nm6lQys3O55Mm5/PalpezYlx3p0kR+QQEhEiH9Wx3HB7edwui+JzJ98Qb6P/QJ09LWqxNbygwFhEgEJcbHcufAlrx5U29OrFeV3728jEufnMfXm/dEujQRBYRIWdCyQXWmXduD+y9sx9db9jDo0Vn8492v2J+tAQAlchQQImVEIGBc2qUJM2/ry3mdknjik28545FP+firLZEuTSqokALCzHqFsk5Ejl2dqpV48OIOTB3VnYS4GEZMWciYfy3ih11ZkS5NKphQzyAeC3GdiJSQ7s3q8PZNffjtWSfx0VdbOPORT3l9yYZIlyUVyGHngzCzHkBPoJ6Z3VZgU3Xy53gQkTCKjw1wfb/mDGrXkNumLeHmqUv4YOVm7j2vLTUT4yNdnkS54s4g4oGq5AdJtQKv3cBF4S1NRH6SXLcKL13bg9vPbMG7y3/grLGf8enXWyNdlkS5kCYMMrMT3P27UqjnmGjCIKkIlm/Yxa0vLuGbLXu5svsJ3DWoJYnxYZscUqLc4SYMCvWnaoqZ/U+SuPtpx1SZiByxtkk1eOPG3jzw3mqemb2W2enbePiSDnRqUivSpUmUCbWT+nbgt8HXH4ElQLF/qpvZADNbbWbpZnZnEdvNzMYFty8zs86FtseY2WIzezPEOkUqhIS4GP44uDX//nU3DhzM5aKJc3n4/dUczNXgf1JyQgoId19U4PW5u98GdDvcPmYWAzwODARaA8PMrHWhZgOBlOBrFDCh0PabgVWh1ChSEfU8sS7v3noKQzoez7iP0rngiTmkb9FT2FIyQn0OonaBV10zOwtoUMxuXYF0d1/j7tnAVGBIoTZDgOc83zygppk1DH5mI+BsYNKRfCGRiqZ6QhwPX9KRCZd3JuPHTM4eN5tnP19LnubDlmMUah/EIsABA3KAtcA1xeyTBKwvsJzB/551FNUmCdgEjAV+R/5dU4dkZqPIP/ugSZMmxZQkEr0GtmvIyU1rccfLy/jzGyv5cNVmHrioA8fXrBzp0qScCvUSU7K7Nwv+N8Xdz3T32cXsZkW9VShtzGwwsMXdF4VQ21PunuruqfXq1SuuuUhUq18tgcnDu/C389ux+PudnDX2M6YvztAIsXJUQr3ElGBmt5nZq2b2ipndamYJxeyWATQusNwI2Bhim17AuWa2jvxLU6eZ2b9CqVWkojMzLuvWhHdu7kOL46px64tLueHfi/lR803IEQr1LqbngDbkD68xHmgFPF/MPguBFDNLNrN4YCgwo1CbGcBVwbuZugO73H2Tu9/l7o3cvWlwv4/c/YoQaxUR4IQ6VZh2bQ9+N+Ak3l+Z/3Ddx6s18J+ELtQ+iJPcvUOB5Y/NbOnhdnD3HDO7AXiP/GE5Jrv7CjMbHdw+EXgbGASkA5nAiCP9AiJyaDEB47pTm9O3RT1ufXEJI55dyOXdmvD7s1vp4TopVqhPUk8BJgbvNMLMugFXu/t14S3vyOhJapFDyzqYy0Pvr2bS7LWcUDuRhy7pyMkn6OG6iu5wT1KHeompGzDHzNYF+wXmAn3N7EszW1ZCdYpIGCXExfD7s1vzn19352Cuc/HEOTz43mqyc/RwnRQt1HPMAWGtQkRKTfdmdXj3lj78+Y2VjP84nY9Xb2HspR1JOe6wd5RLBRTqGcS97v5dwVfBdeEsUERKXrWEOB68uANPXnkym3ZlcfZjs5k0a40erpNfCDUg2hRcMLNY4OSSL0dEStNZbRrw3i2ncEpKXe59axWXT5rPhp37I12WlBGHDQgzu8vM9gDtzWy3me0JLm8GXi+VCkUkrOpVq8TTV6Vy/4XtWJaxkwGPfMYri/RwnRQTEO7+/9y9GvCAu1d392rBVx13v6uUahSRMDMzLu3ShHduPoWWDavxm5eWMuZfX7BDD9dVaKHe5npKUevd/bMSr+gY6DZXkWOXm+c8PWsND7//NdUrx/GPi9pxWsvjIl2WhMnhbnMNNSDeKLCYQP5IrYvK2oRBCgiRkrNq025ufXEJX/2wh2FdG/OHs1tTpZIeros2x/wchLufU+B1BtCW/H4IEYlSrRpW5/UbenFt32ZMXbiegY/OIm3djkiXJaUo1LuYCssgPyREJIpVio3hroGteHFUD/LcueTJuYz98GvdDltBhHS+aGaP8fNQ3QGgE3DYsZhEJHp0Ta7Nu7ecwp9eW87YD79hyfqdjL20IzUT4yNdmoRRqGcQK4GvgdXAPOB3Gl1VpGKpWimWhy7pwL3ntWVO+nbOHjebLzN2RbosCaPinoOINbN/AH8FfgWMJH+mtyFmFhf+8kSkLDEzruh+AtNG98DduXDiHKYu+D7SZUmYFHcG8QBQG0h2987u3gloBtQEHgxzbSJSRnVsXJM3b+pDt+Ta3Pnql/zu5aVkHcyNdFlSwooLiMHAr919z08r3H03MIb8eRxEpIKqXSWeKSO6ctNpzZmWlsGFE+bw/fbMSJclJai4gHAv4kEJd8/lf+eXFpEKJiZg3HbmSUwensr6HZkMfmwWM1fpDvhoUVxArDSzqwqvNLMrgK/CU5KIlDentTyOt27qQ+PaiVzzzzQeen81uboVttw77JPUZpYEvArsBxaRf9bQBagMnO/uG0qjyFDpSWqRyMo6mMufXl/OtLQM+qTU5dGhnahdRbfClmVH/SS1u29w927AX4B1wPfAX9y9ayjhYGYDzGy1maWb2Z1FbDczGxfcvszMOgfXJ5jZAjNbamYrzOzPIXxPEYmwhLgY/nFRB+6/sB3z1+5g8LhZLFm/M9JlyVEKdaiNj9z9MXcf5+4zQ9nHzGKAx4GBQGtgmJm1LtRsIJASfI0CJgTXHwBOc/cOQEdggJl1D+VzRSTyLu3ShFdG9yQQMC6eOIfn532n4cPLoaMdaiMUXYF0d1/j7tnAVGBIoTZDgOc83zygppk1DC7vDbaJC7700yVSjrRrVIM3b+xN7+Z1+eNry/nNtKXsz9atsOVJOAMiCVhfYDkjuC6kNmYWY2ZLgC3AB+4+v6gPMbNRZpZmZmlbt24tqdpFpATUTIznmau7cNsZLZi+ZAPnP/E5a7fti3RZEqJwBoQVsa7wWcAh27h7rrt3BBoBXc2syMEB3f0pd09199R69eodS70iEgaBgHFT/xSmjOjKD7uzOPex2by/4odIlyUhCGdAZACNCyw3AjYeaRt33wl8Agwo8QpFpNT0bVGPN2/sTXK9Kox6fhF/f+crcnLzIl2WHEY4A2IhkGJmyWYWDwwFZhRqMwO4Kng3U3dgl7tvMrN6ZlYTwMwqA6ej5y5Eyr1GtRJ5aXQPLuvWhImffsuVzyxg294DkS5LDiFsAeHuOcANwHvAKmCau68ws9FmNjrY7G1gDZAOPA1cF1zfEPjYzJaRHzQfuPub4apVREpPpdgY/nZ+Ox64qD1ffP8jg8fNZtF3P0a6LClCSFOOlhd6UE6kfFmxcRdj/vUFG3fu5w9nt+Lqnk0xK6prUsLlmKccFREJhzbH1+CNG3tz6kn1uOeNldw8dQmZ2TmRLkuCFBAiElE1Ksfx1JWp/Pask3hz2UbOe/xzvt26t/gdJewUECIScYGAcX2/5jz3q25s25vNkPGf886XmyJdVoWngBCRMqN3Sl3evLE3zetXZcwLX3DfWyt1K2wEKSBEpEw5vmZlpl3bg6t6nMDTs9Zy2aT5bNmTFemyKiQFhIiUOfGxAf4ypC1jL+3IsoydnD1uNgvW7oh0WRWOAkJEyqzzOiXx2vW9qFoplmFPz2PSrDUaFbYUKSBEpExr2aA6r9/Qi9Nb1efet1Zxw78Xs/eAboUtDQoIESnzqifEMfGKk7lrYEveWb6JIeNna1TYUqCAEJFywcy4tu+JvDCyOz9mHuSCJz5n0XfqlwgnBYSIlCs9TqzDq2N6UjMxnmFPz+dtPS8RNgoIESl3mtatwitjetIuqQbX//sLdV6HiQJCRMql2lXieWFkNwa2bcC9b63inhkryM1TSJQkBYSIlFsJcTGMH9aZUac0459zv+Pa5xdpsL8SpIAQkXItEDD+b1Ar/jKkDR99tZlhT81j6x5NQlQSFBAiEhWu6tGUJ69MZfXmPVwwQSPClgQFhIhEjTNaH8eLo3qwPzuXC56Yo+E5jpECQkSiSofGNZl+XS/qVI3niknzeWPpxkiXVG6FNSDMbICZrTazdDO7s4jtZmbjgtuXmVnn4PrGZvaxma0ysxVmdnM46xSR6NK4diKvjulJxyY1ufE/i5nwybe6DfYohC0gzCwGeBwYCLQGhplZ60LNBgIpwdcoYEJwfQ7wG3dvBXQHri9iXxGRQ6qZGM/z13TlnA7Hc/+7X/GH15ZrbokjFBvG9+4KpLv7GgAzmwoMAVYWaDMEeM7zo32emdU0s4buvgnYBODue8xsFZBUaF8RkcOqFBvDo5d2pFGtykz45Fs27crisWGdqFIpnL/6okc4LzElAesLLGcE1x1RGzNrCnQC5hf1IWY2yszSzCxt69atx1qziESZQMC4Y0BL7ju/LZ+s3sKlT81ly25NQBSKcAaEFbGu8EXAw7Yxs6rAK8At7r67qA9x96fcPdXdU+vVq3fUxYpIdLu82wk8c3UX1mzdx/lPzOGbzXsiXVKZF86AyAAaF1huBBS+neCQbcwsjvxweMHdXw1jnSJSQfRrWZ9p1/YgOzePCybMYe632yNdUpkWzoBYCKSYWbKZxQNDgRmF2swArgrezdQd2OXum8zMgGeAVe7+cBhrFJEKpm1SDaZf15MG1RO4avJ8Xlu8IdIllVlhCwh3zwFuAN4DVgHT3H2FmY02s9HBZm8Da4B04GnguuD6XsCVwGlmtiT4GhSuWkWkYmlUK5GXx/Qk9YTa3PLiEsZ/9I1ugy2CRdNBSU1N9bS0tEiXISLlRHZOHne8sozpizcwtEtj/npeW+JiKtbzw2a2yN1Ti9qme71EpMKKjw3w8CUdaFSrMo99lM7GXVk8cXlnquo2WEBDbYhIBWdm/ObMk7j/wnZ8nr6NSybOZbNugwUUECIiAFzapQmTh3fhu+37OO/xz/nqhyLvrK9QFBAiIkF9W9Rj2uge5Llz8YS5fJ6+LdIlRZQCQkSkgDbH12D6db1IqlWZqycv4OVFGZEuKWIUECIihRxfszLTRvege7M63P7SUh79sGLeBquAEBEpQvWEOJ4d0YWLTm7EIx9+zW9fXkZ2TsUaDVb3comIHEJcTIAHLmpP41qJPPLh1/ywK4snruhM9YS4SJdWKnQGISJyGGbGzaen8ODFHZi3ZjuXTJzLpl37I11WqVBAiIiE4KKTGzFlRFc2/Lif8x7/nJUbo/82WAWEiEiIeqfU5aUxPQiYccmTc/n06+ieg0YBISJyBFo2qM7063rRuHYi10xZyNtfbop0SWGjgBAROUINaiQw7drudGpSkxv/s5g3lhae6iY6KCBERI5CtYQ4pozoyskn1OLmqYt5fUn0zSuhgBAROUpVKsUyZUQXuiXX4dYXl/BKlD11rYAQETkGifGxTB7ehZ4n1uX2l5cyLW19pEsqMQoIEZFjVDk+hklXp9K7eV1+9/Iypi74PtIllYiwBoSZDTCz1WaWbmZ3FrHdzGxccPsyM+tcYNtkM9tiZsvDWaOISElIiIvh6atSOfWketz56pe8MP+7SJd0zMIWEGYWAzwODARaA8PMrHWhZgOBlOBrFDChwLYpwIBw1SciUtIS4mJ48sqT6d+yPr+fvpzn5q6LdEnHJJxnEF2BdHdf4+7ZwFRgSKE2Q4DnPN88oKaZNQRw98+AHWGsT0SkxFWKjWHCFSdzRuvj+NPrK5g8e22kSzpq4QyIJKBgb01GcN2RthERKVfiYwM8fllnBrRpwF/eXMmkWWsiXdJRCWdAWBHrCg+oHkqbw3+I2SgzSzOztK1bo/uxdxEpP+JjAzx2WSfObteQe99axcRPv410SUcsnMN9ZwCNCyw3Ago/bhhKm8Ny96eApwBSU1Mr3oweIlJmxcUEeHRoRwIB4+/vfEVunnN9v+aRLitk4QyIhUCKmSUDG4ChwGWF2swAbjCzqUA3YJe7R+/AJiJS4cTGBHjkkg7EBowH3ltNbp5zU/+USJcVkrAFhLvnmNkNwHtADDDZ3VeY2ejg9onA28AgIB3IBEb8tL+Z/Qc4FahrZhnA3e7+TLjqFREJl9iYAA9e3IGAGQ9/8DU5ec6tp6dgVtRV9rIjrDPKufvb5IdAwXUTC/zbgesPse+wcNYmIlKaYgLGAxe1JyYA42Z+Q16e85szW5TpkNCUoyIipSQQMP5+QXtiAsb4j9PJyXPuGHBSmQ0JBYSISCkKBIz7zmtHTMCY+Om35Obl8X+DWpXJkFBAiIiUskDA+OuQtsQGAjw9ay05ec6fBrcucyGhgBARiQAz4+5zWhMwY/Lna8nLc+45t02ZCgkFhIhIhJgZfxzcitgY46nP1pDrzl/ObUsgUDZCQgEhIhJBZsZdA1sSsJ/6JJz7zmtXJkJCASEiEmFmxh0DTiL2p7ubcp2/X5h/t1MkKSBERMoAM+M3Z7YgJmA8OvMbct154KIOEQ0JBYSISBlhZtx6Rn5IPPzB1+TlOQ9e3IHYmMhM/qmAEBEpY27qn5L/5PV7q8l18sdyikBIKCBERMqg6/s1JyY4CmxenjN2aEfiSjkkFBAiImXU6L4nEhsw7n1rFTl5eTw2rDPxsaUXEpG5sCUiIiEZ2acZd5/TmvdWbOa6F77gQE5uqX22AkJEpIwb0SuZvwxpw4erNjPmX1+QdbB0QkIBISJSDlzVoyn3nd+Wj77awrXPLyqVkFBAiIiUE5d3O4G/X9COz77Zyq+fSwt7SCggRETKkaFdm3D/he2Znb6Na/65kP3Z4QsJBYSISDlzSWpjHrq4A3O/3c6IKQvIzM4Jy+eENSDMbICZrTazdDO7s4jtZmbjgtuXmVnnUPcVEanILujciEcu7ciCtTsYPnlhWEIibM9BmFkM8DhwBpABLDSzGe6+skCzgUBK8NUNmAB0C3FfEZEKbUjHJAJmzP5mG5ViY0r8/cP5oFxXIN3d1wCY2VRgCFDwl/wQ4Dl3d2CemdU0s4ZA0xD2FRGp8M7pcDzndDg+LO8dzktMScD6AssZwXWhtAllXwDMbJSZpZlZ2tatW4+5aBERyRfOgChqjFoPsU0o++avdH/K3VPdPbVevXpHWKKIiBxKOC8xZQCNCyw3AjaG2CY+hH1FRCSMwnkGsRBIMbNkM4sHhgIzCrWZAVwVvJupO7DL3TeFuK+IiIRR2M4g3D3HzG4A3gNigMnuvsLMRge3TwTeBgYB6UAmMOJw+4arVhER+V+WfwNRdEhNTfW0tLRIlyEiUm6Y2SJ3Ty1qm56kFhGRIikgRESkSFF1icnMtgLfHeXudYFtJVhOeaZj8Us6Hr+k4/GzaDgWJ7h7kc8IRFVAHAszSzvUdbiKRsfil3Q8fknH42fRfix0iUlERIqkgBARkSIpIH72VKQLKEN0LH5Jx+OXdDx+FtXHQn0QIiJSJJ1BiIhIkRQQIiJSpAoVEMcyBWo0CuF4tDSzuWZ2wMxuj0SNpSmE43F58OdimZnNMbMOkaizNIRwLIYEj8OS4HwsvSNRZ2kJdQpkM+tiZrlmdlFp1hc27l4hXuQP+vct0Iz84cSXAq0LtRkEvEP+fBTdgfmRrjvCx6M+0AW4D7g90jWXgePRE6gV/PfAaP35CPFYVOXnPsz2wFeRrjuSx6NAu4/IH4T0okjXXRKvinQG8d8pUN09G/hpGtOC/jsFqrvPA36aAjUaFXs83H2Luy8EDkaiwFIWyvGY4+4/BhfnkT9PSTQK5Vjs9eBvRaAKh5jQK0qE8rsD4EbgFWBLaRYXThUpII5lCtRoVJG+ayiO9HhcQ/7ZZjQK6ViY2flm9hXwFvCrUqotEoo9HmaWBJwPTCzFusKuIgXEsUyBGo0q0ncNRcjHw8z6kR8Qd4S1osgJ6Vi4+3R3bwmcB/w13EVFUCjHYyxwh7vnhr+c0hPOKUfLmmOZAjUaVaTvGoqQjoeZtQcmAQPdfXsp1Vbajuhnw90/M7MTzayuu5f3geuKEsrxSAWmmhnkD+A3yMxy3P21UqkwTCrSGcSxTIEajTSt6y8VezzMrAnwKnClu38dgRpLSyjHorkFfxsG7/aLB6I1MIs9Hu6e7O5N3b0p8DJwXXkPB6hAZxB+DFOgRqNQjoeZNQDSgOpAnpndQv7dG7sjVXe4hPjz8SegDvBE8HdjjkfhSJ4hHosLyf9j6iCwH7i0QKd1VAnxeEQlDbUhIiJFqkiXmERE5AgoIEREpEgKCBERKZICQkREiqSAEBGRIikgRESkSAoIEREp0v8H9tX/O3yNjdEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "τ_k = np.arange(0,.5,.05)\n",
    "α_k = np.zeros(τ_k.size)\n",
    "γ_k = np.zeros(τ_k.size)\n",
    "Φ_k = np.zeros(τ_k.size)\n",
    "for k in range(τ_k.size ):\n",
    "    # Φτ_x_y = α_x_y + (1-τ_k[k])* γ_x_y\n",
    "    lintax_mkt = TU_model(α_x_y + (1-τ_k[k])* γ_x_y, T = 1-τ_k[k] / 2)\n",
    "    lintax_mkt.solveChooSiowIPFP()\n",
    "    α_k[k] = np.sum(lintax_mkt.eq_μ_x_y * (α_x_y - np.log(lintax_mkt.eq_μ_x_y)) )\n",
    "    γ_k[k] = np.sum(lintax_mkt.eq_μ_x_y * (γ_x_y - np.log(lintax_mkt.eq_μ_x_y)) )\n",
    "    Φ_k[k] = np.sum(lintax_mkt.eq_μ_x_y * (α_x_y+γ_x_y - 2 * np.log(lintax_mkt.eq_μ_x_y)))\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(τ_k,α_k)\n",
    "plt.legend(['Total α'], loc='upper left')\n",
    "plt.ylabel('Output')\n",
    "plt.show()\n",
    "\n",
    "plt.plot(τ_k,γ_k)\n",
    "plt.legend(['Total γ'], loc='upper left')\n",
    "plt.ylabel('Output')\n",
    "plt.show()\n",
    "\n",
    "plt.plot(τ_k,Φ_k)\n",
    "plt.legend(['Total α+γ'], loc='upper left')\n",
    "plt.ylabel('Output')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}