D4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Day 4: matching models with nontransferable utility</center>\n",
    "### <center>Alfred Galichon (NYU+ScPo)</center>\n",
    "## <center>'math+econ+code' masterclass on equilibrium transport and matching models in economics</center>\n",
    "<center>© 2020-2022 by Alfred Galichon.  Support from  NSF DMS-1716489 and ERC CoG-866274 EQUIPRICE grants is acknowledged.</center>\n",
    "\n",
    "#### <center>with Python code</center>\n",
    "\n",
    "**If you reuse material from this masterclass, please cite as:**<br>\n",
    "Alfred Galichon, 'math+econ+code' masterclass on equilibrium transport and matching models in economics, June 2022. https://github.com/math-econ-code/mec_equil"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# References\n",
    "\n",
    "## Textbooks\n",
    "\n",
    "* [RS90] Alvin Roth and Marilda Sotomayor (1990). *Two-sided matching*. Econometric Society Monographs, Cambridge University Press.\n",
    "\n",
    "## Papers\n",
    "\n",
    "* [GS62] David Gale and Lloyd Shapley (1962). \"College Admissions and the Stability of Marriage.\" *American Mathematical Monthly* 69 (1), pp. 9–14.\n",
    "\n",
    "* [A00] Hiroyuki Adachi (2000). \"On a characterization of stable matchings.\" *Economics Letters* 68 pp. 43–49.\n",
    "\n",
    "* [GL21] Alfred Galichon, Flavien Léger (2021). \"Matching algorithms as coordinate update\". In progress. \n",
    "\n",
    "* [NY09] Muriel Niederle and Leeat Yariv (2009). \"Decentralized Matching with Aligned Preferences.\" *NBER Working Paper* No. w14840.\n",
    "\n",
    "* [GGH21] Alfred Galichon, Octavia Ghelfi and Marc Henry (2021). \"Stable and extremely unequal.\" Working paper.\n",
    "\n",
    "* [GH19] Alfred Galichon, and Yu-Wei Hsieh.\"Aggregate stable matching with money burning.\" SSRN id=2887732. \n",
    "\n",
    "* [GHS21] Alfred Galichon, Yu-Wei Hsieh and Maxime Sylvestre. \"Monotone comparative statics for submodular functions, with an application to aggregated deferred acceptance.\" In progress. \n",
    "\n",
    "* [GKW19] Alfred Galichon, Scott Kominers, and Simon Weber (2019). \"Costly Concessions: An Empirical Framework for Matching with Imperfectly Transferable Utility.\" *Journal of Political Economy* 127 no. 6, pp. 2875-2925.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Setting up the model\n",
    "\n",
    "## Populations\n",
    "\n",
    "* As yesterday we consider a population of workers $x\\in\\mathcal{X}$ and firms $y\\in\\mathcal{Y}$, and assume that $n_x$ is the number of workers of type $x$, while $m_y$ is the number of firms of type $y$.\n",
    "* We shall sometimes assume that there is one individual per type, that is $n_x = 1$ and $m_y = 1$, in which case we shall talk about *invidividual matching,* and sometimes allow for more than one agent per type, a situation we shall refer to as *aggregate matching*. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preferences\n",
    "\n",
    "* Assume the wages are set. We describe the preferences as follows:<br>\n",
    "    $\\alpha_{xy}$ = valuation of an $xy$ match by $x$<br>\n",
    "    $\\gamma_{xy}$ = valuation of an $xy$ match by $y$<br>\n",
    "    If $x$ and $y$ remain unmatched, they get respective utilities $\\alpha_{x0}$ and $\\gamma_{0y}$.<br>\n",
    "* As before we denote $\\mathcal{X}_0 = \\mathcal{X} \\cup \\{0\\}$ and $\\mathcal{Y}_0 = \\mathcal{Y} \\cup \\{0\\}$.\n",
    "* Assume strict preferences:<br>\n",
    "    $\\alpha_{xy} \\neq \\alpha_{xy^\\prime}$ for $y\\neq y^\\prime$<br>\n",
    "    $\\gamma_{xy} \\neq \\gamma_{x^\\prime y}$ for $x \\neq x^\\prime$.\n",
    "\n",
    "* The preferences are usually described in an ordinal way, by specifying a order relation for each agent on the other side of the market."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Representing the model\n",
    "\n",
    "Let's first load the libraries we will need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import networkx as nx\n",
    "import scipy.optimize as opt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create the `NTU_market` class to encompass that information and store it in a convenient way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class NTU_market:\n",
    "    def __init__(self,α_x_y,γ_x_y,n_x = np.array([]), m_y = np.array([])):\n",
    "        nbx,nby = α_x_y.shape\n",
    "        self.α_x_y = np.hstack((α_x_y,np.zeros((nbx,1) )))\n",
    "        self.γ_x_y = np.vstack((γ_x_y,np.zeros( (1,nby) )))\n",
    "        if n_x.size == 0:\n",
    "            n_x = np.ones(nbx)\n",
    "        if m_y.size ==  0:\n",
    "            m_y = np.ones(nby) \n",
    "        self.n_x,self.m_y = n_x, m_y\n",
    "        self.largex,self.largey = nby+1, nbx+1\n",
    "        self.smallx,self.smally = -1, -1\n",
    "        self.nbx,self.nby = nbx, nby\n",
    "        self.αo_x_y = np.zeros((nbx,nby+1), dtype = 'int64') \n",
    "        self.γo_x_y = np.zeros((nbx+1,nby), dtype = 'int64') \n",
    "        self.prefslistα_x_y = np.zeros((nbx,nby+1), dtype = 'int64') \n",
    "        self.prefslistγ_x_y = np.zeros((nbx+1,nby), dtype = 'int64') \n",
    "        self.traceuo_x_t = np.array([])\n",
    "        self.traceu_x_t = np.array([])\n",
    "        for x in range(nbx):\n",
    "            thelistx = (- self.α_x_y )[x,:].argsort()\n",
    "            self.αo_x_y[x, thelistx] =  nby  - np.arange(nby+1)\n",
    "            self.prefslistα_x_y[x,:] = (thelistx) % (self.nby+1) \n",
    "        for y in range(nby):\n",
    "            thelisty = ( - self.γ_x_y)[:,y].argsort()\n",
    "            self.γo_x_y[thelisty, y] = nbx  - np.arange(nbx+1)\n",
    "            self.prefslistγ_x_y[:,y] = (thelisty) % (self.nbx+1) \n",
    "        self.comp_nbsteps = -1\n",
    "        self.comp_time = -1\n",
    "        self.eq_μ_x_y = np.array([])   \n",
    "        self.eq_u_x = np.array([])\n",
    "        self.eq_v_y = np.array([])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following `print_prefs` method prints the ordinal preferences of the agents:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def list_prefs(x,xletter,yletter,αo_y):\n",
    "    nby = len(αo_y) - 1\n",
    "    xsprefs = xletter + str(x) +' : ' \n",
    "    for y in range(nby):\n",
    "        ind = np.where(αo_y== (nby - y ))[0][0]\n",
    "        if (ind == nby):\n",
    "            break        \n",
    "        xsprefs = xsprefs + yletter + str( ind ) + ' > '\n",
    "    return(xsprefs[:-3])\n",
    "\n",
    "def print_prefs(self,xs=[],ys=[]):\n",
    "    if xs == [] and ys==[] :\n",
    "        xs = range(self.nbx)\n",
    "        ys = range(self.nby)\n",
    "    if xs != [] :\n",
    "        for x in xs :\n",
    "            print(list_prefs(x,'x','y',self.αo_x_y[x,:]))\n",
    "        \n",
    "    if xs != [] and ys != [] :\n",
    "        print('===')\n",
    "        \n",
    "    if ys != [] :\n",
    "        for y in ys :\n",
    "            print(list_prefs(y,'y','x',self.γo_x_y[:,y]))\n",
    "\n",
    "NTU_market.print_prefs = print_prefs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We simulate such a market and print the corresponding preferences:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(seed=1000)\n",
    "running_mkt = NTU_market(np.random.rand(3,2)-0.1,np.random.rand(3,2)-0.2)\n",
    "running_mkt.print_prefs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following example is taken from [RS90], example 2.17"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# example 2.17 in [RS90]\n",
    "α_x_y = np.array([[0,1,2,3],[1,0,3,2],[2,3,0,1],[3,2,1,0]])\n",
    "γ_x_y = np.array([[3,2,1,0],[2,3,0,1],[1,0,3,2],[0,1,2,3]])\n",
    "rs_ex_2_17 =  NTU_market(α_x_y,γ_x_y)\n",
    "rs_ex_2_17.print_prefs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stable matchings\n",
    "\n",
    "We shall see two notions of stable matchings with non-transferable utility (NTU):<br> \n",
    "* First, the historical notion (due to [GS62]) of *Gale-Shapley NTU stable matchings*. \n",
    "* Next, a related notion of stability that will fit into the framework we saw yesteray, which we shall call *aggregate NTU stable matchings*, following [GH19] and [GHS21]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stability in the Gale-Shapley sense\n",
    "\n",
    "**Assume that there is one individual per type:**<br>\n",
    "$n_x = 1$ and $m_y = 1$ for all $x\\in\\mathcal{X}$ and $y\\in\\mathcal{Y}$.<br>\n",
    "\n",
    "As before, $\\mu_{xy}\\in\\{0,1\\}$ will be a dummy variable equal to 1 iff $x$ is matched with $y$,<br>\n",
    "    and $u_x$ and $v_y$ are the payoffs of worker $x$ and firm $y$ at equilibrium.\n",
    "\n",
    "**Definition.** $(\\mu,u,v)$ is a stable outcome in the Gale-Shapley sense if:\n",
    "\n",
    "(i) Population constraints are satisfied<br>\n",
    "$\\left\\{\n",
    "\\begin{array}[l]\n",
    "~\\sum_{y} \\mu_{xy} + \\mu_{x0} = n_{x}=1 \\\\\n",
    "\\sum_{x}\\mu_{xy} + \\mu_{0y} = m_{y}=1\n",
    "\\end{array}\n",
    "\\right.$<br>\n",
    "    \n",
    "(ii) There is no blocking pair and all individuals are rational:<br>\n",
    "$ \\max \\{ u_x - \\alpha_{xy},v_y - \\gamma_{xy} \\} \\geq 0$, and<br>\n",
    "$u_x \\geq 0$ and $v_y \\geq 0$\n",
    "\n",
    " \n",
    "(iii) Strong complementarity:<br>\n",
    "$\\mu_{xy}>0$ implies $u_x = \\alpha_{xy}$ and $v_y = \\gamma_{xy}$<br>\n",
    "$\\mu_{x0}>0$ implies $u_x = 0$<br>\n",
    "$\\mu_{0y}>0$ implies $v_y = 0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remark\n",
    "\n",
    "* In the Gale-Shapley definition:<br>\n",
    "$\\mu_{xy}>0$ implies $u_x = \\alpha_{xy}$ and $v_y = \\gamma_{xy}$<br>\n",
    "\n",
    "* This is stronger than assuming<br>\n",
    "$\\mu_{xy}>0$ implies $\\max \\{ u_x - \\alpha_{xy},v_y - \\gamma_{xy} \\} \\geq 0$<br>\n",
    "Indeed, in the latter case, $(u_x,v_y)$ needs to be on the NTU frontier but doesn't have to be the efficient point $(\\alpha_{xy},\\gamma_{xy})$ as required in [GS62].\n",
    "\n",
    "* [GH19] and [GHS21] impose the latter condition as the definition of stability, which is weaker than Gale-Shapley. We shall see later why it may make sense not to impose the outcome to be Pareto efficient. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function detects stable matchings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_GS_stable(self, μ_x_y = None, output=0 ):\n",
    "    if (min(self.n_x)<1) or (max(self.n_x)>1) or (min(self.m_y)<1) or (max(self.m_y)>1):\n",
    "        if output > 0 :\n",
    "            print('n_x or m_y do not only contain ones.')\n",
    "        return(False)\n",
    "    if (μ_x_y is None):\n",
    "        μ_x_y = self.eq_μ_x_y\n",
    "    μext_x_y0 = np.hstack([μ_x_y, 1-np.sum(μ_x_y,axis = 1).reshape(-1,1) ])\n",
    "    μext_x0_y = np.vstack([μ_x_y, 1-np.sum(μ_x_y,axis = 0).reshape(1,-1) ])\n",
    "    if (np.logical_and(μext_x_y0 != 0  , μext_x_y0 != 1 )).any() or (np.logical_and(μext_x0_y != 0  , μext_x0_y != 1 )).any():\n",
    "        if output > 0 :\n",
    "            print('The μ is not feasible.')\n",
    "        return(False)\n",
    "    uo_x = np.sum(μext_x_y0 * self.αo_x_y, axis = 1)\n",
    "    vo_y = np.sum(μext_x0_y * self.γo_x_y, axis = 0)\n",
    "    for x in range(self.nbx):\n",
    "        for y in range(self.nby):\n",
    "            if (self.αo_x_y[x,y] > uo_x[x]) and (self.γo_x_y[x,y] > vo_y[y]):\n",
    "                if output > 0 :\n",
    "                    print('The matching is not stable. One blocking pair is: x'+str(x)+',y'+str(y)+'.')\n",
    "                return(False)\n",
    "\n",
    "    for x in range(self.nbx):\n",
    "        if (self.αo_x_y[x,self.nby] > uo_x[x]):\n",
    "            if output > 0 :\n",
    "                print('The matching is not stable. One blocking pair is: x'+str(x)+',y0.')\n",
    "            return(False)\n",
    "\n",
    "    for y in range(self.nby):\n",
    "        if (self.γo_x_y[self.nbx,y] > vo_y[y]):\n",
    "            if output > 0 :\n",
    "                print('The matching is not stable. One blocking pair is: x0,y'+str(y)+'.')\n",
    "            return(False)\n",
    "            \n",
    "    if output > 0 :\n",
    "        print('The matching is stable.')\n",
    "    return (True)\n",
    "\n",
    "NTU_market.is_GS_stable = is_GS_stable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's go back to example 2.17 in [RS90], and let's check the stability of the following matching: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "μ_x_y = np.array([[0,0,0,1],[1,0,0,0],[0,1,0,0],[0,0,1,0]])\n",
    "_ = rs_ex_2_17.is_GS_stable(μ_x_y,output=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try again with the following alternative matching:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "μ_x_y = np.array([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]])\n",
    "_ = rs_ex_2_17.is_GS_stable(μ_x_y,output=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some useful functions\n",
    "\n",
    "Let's construct some useful functions that will be useful later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def μ_from_uo(self, uo_x):\n",
    "    return np.where( (self.αo_x_y == uo_x.reshape((-1,1))) , 1 , 0 )[:,0:-1]\n",
    "NTU_market.μ_from_uo = μ_from_uo\n",
    "\n",
    "def μ_from_vo(self, vo_y):\n",
    "    return np.where( self.γo_x_y == vo_y , 1 , 0 )[0:-1,:]\n",
    "NTU_market.μ_from_vo = μ_from_vo\n",
    "\n",
    "def uo_from_μ(self, μ_x_y = None):\n",
    "    if (μ_x_y is None):\n",
    "        μ_x_y = self.eq_μ_x_y\n",
    "    return np.sum(np.hstack([μ_x_y, (self.n_x-np.sum(μ_x_y,axis = 1)).reshape(-1,1) ]) * self.αo_x_y, axis = 1)\n",
    "NTU_market.uo_from_μ = uo_from_μ\n",
    "\n",
    "def vo_from_μ(self, μ_x_y = None):\n",
    "    if (μ_x_y is None):\n",
    "        μ_x_y = self.eq_μ_x_y\n",
    "    return np.sum(np.vstack([μ_x_y, (self.m_y-np.sum(μ_x_y,axis = 0)).reshape(1,-1) ]) * self.γo_x_y, axis = 0)\n",
    "NTU_market.vo_from_μ = vo_from_μ\n",
    "\n",
    "def u_from_μ(self, μ_x_y = None):\n",
    "    if (μ_x_y is None):\n",
    "        μ_x_y = self.eq_μ_x_y\n",
    "    return np.sum(np.hstack([μ_x_y, (self.n_x-np.sum(μ_x_y,axis = 1)).reshape(-1,1) ]) * self.α_x_y, axis = 1)\n",
    "NTU_market.u_from_μ = u_from_μ\n",
    "\n",
    "def v_from_μ(self, μ_x_y = None):\n",
    "    if (μ_x_y is None):\n",
    "        μ_x_y = self.eq_μ_x_y\n",
    "    return np.sum(np.vstack([μ_x_y, (self.m_y-np.sum(μ_x_y,axis = 0)).reshape(1,-1) ]) * self.γ_x_y, axis = 0)\n",
    "NTU_market.v_from_μ = v_from_μ\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Deferred acceptance: Gale and Shapley's algorithm\n",
    "\n",
    "\n",
    "The reference for this section is [GS62].\n",
    "\n",
    "**Principle:** Workers make offers to firms that have not rejected them yet. Define:\n",
    "\n",
    "* $A^{t} \\subseteq \\mathcal{X} \\times \\mathcal{Y} =$set of available firms to worker at\n",
    "time $t$\n",
    "\n",
    "* $P^{t} \\subseteq \\mathcal{X} \\times \\mathcal{Y} =$set of proposals made by workers at\n",
    "time $t$ \n",
    "\n",
    "* $E^{t}\\subseteq \\mathcal{X} \\times \\mathcal{Y} =$set of proposals kept by firms at the\n",
    "end of round $t$ \n",
    "\n",
    "For  $B \\subseteq \\mathcal{X} \\times \\mathcal{Y}$, we introduce the following notations:\n",
    "\n",
    "* for $x \\in \\mathcal{X}$, denote $B(x)$ the set of $y \\in \\mathcal{Y}$ such that $xy \\in B$.\n",
    "\n",
    "* for $y \\in \\mathcal{Y}$, denote $B(y)$ the set of $x \\in \\mathcal{X}$ such that $xy \\in B$.\n",
    "\n",
    "\n",
    "**Deferred acceptance algorithm of Gale and Shapley.** \n",
    "\n",
    "At time $t=0$, all firms are avaiable to anyone $A^{0} =\\mathcal{X} \\times \\mathcal{Y}$.\n",
    "\n",
    "Iterate over $t$:\n",
    "\n",
    "$\n",
    "\\left\\{\n",
    "\\begin{array}{l}\n",
    "P^{t}=\\left\\{ xy\\in \\mathcal{X}\\times \\mathcal{Y}:y\\in \\arg \\max_{y\\in\n",
    "A^{t}\\left( x\\right) \\cup \\left\\{ 0\\right\\} }\\left\\{ \\alpha _{xy}\\right\\}\n",
    "\\right\\} \\\\\n",
    "E^{t}=\\left\\{ xy\\in \\mathcal{X}\\times \\mathcal{Y}:x\\in \\arg \\max_{x\\in\n",
    "P^{t}\\left( y\\right) \\cup \\left\\{ 0\\right\\} }\\left\\{ \\gamma _{xy}\\right\\}\n",
    "\\right\\} \\\\\n",
    "A^{t+1} =A^{t} \\backslash \\left\\{ P^{t} \\backslash  E^{t} \\right\\} \n",
    "\\end{array}\n",
    "\\right.\n",
    "$\n",
    "\n",
    "\n",
    "\n",
    "Repeat until $P^{t} =E^{t} $, i.e. no offer is rejected."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's implement as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solveGaleShapley(self,output=0, trace=False):\n",
    "    if (output>=2):\n",
    "        print(\"Offers made and kept are denoted +1; offers made and rejected are denoted -1.\")\n",
    "    self.comp_nbsteps = 0\n",
    "    tracemax = self.nbx*self.nby\n",
    "    if trace:\n",
    "        self.traceuo_x_t = np.zeros((tracemax,self.nbx))    \n",
    "    μA_x_y = np.ones((self.nbx+1, self.nby+1), dtype = 'int64') # initially all offers are non rejected\n",
    "    while True :\n",
    "        μP_x_y = np.zeros((self.nbx+1, self.nby+1), dtype = 'int64')\n",
    "        props_x = np.ma.masked_array(self.αo_x_y, μA_x_y[0:-1,:] ==0).argmax(axis = 1) # x's makes an offer to their favorite y\n",
    "        μP_x_y[range(self.nbx),props_x] = 1 \n",
    "        μP_x_y[self.nbx,0:self.nby] = 1\n",
    "        \n",
    "        μE_x_y = np.zeros((self.nbx+1, self.nby+1), dtype = 'int64') # y's retains their favorite offer among those made:\n",
    "        rets_y = np.ma.masked_array(self.γo_x_y, μP_x_y[:,0:-1] ==0).argmax(axis = 0)\n",
    "        μE_x_y[rets_y,range(self.nby)] = 1\n",
    "\n",
    "        rej_x_y = μP_x_y - μE_x_y # compute rejected offers\n",
    "        rej_x_y[self.nbx,:] = 0\n",
    "        rej_x_y[:,self.nby] = 0\n",
    "\n",
    "        μA_x_y = μA_x_y - rej_x_y  # offers from x that have been rejected are no longer available to x\n",
    "        if output >= 2:\n",
    "            print(\"Round \"+str(self.comp_nbsteps)+\":\\n\" , ( (2*μE_x_y-1) * μP_x_y)[0:-1,0:-1])\n",
    "        if trace and self.comp_nbsteps<tracemax:\n",
    "            for x in range(self.nbx):\n",
    "                self.traceuo_x_t[self.comp_nbsteps,x] = np.sum(self.αo_x_y [x,:] * μP_x_y[x,:])\n",
    "        self.comp_nbsteps +=1\n",
    "        if np.max(np.abs(rej_x_y)) == 0: \n",
    "            if trace:\n",
    "                self.traceuo_x_t = self.traceuo_x_t[0:min(self.comp_nbsteps,tracemax),:]\n",
    "            break # if all offers have been accepted, then algorithm stops\n",
    "    self.eq_μ_x_y = μE_x_y[0:-1,0:-1]\n",
    "    return (0)\n",
    "\n",
    "NTU_market.solveGaleShapley = solveGaleShapley \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test on our example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "running_mkt.solveGaleShapley(output=2,trace=True)\n",
    "running_mkt.is_GS_stable()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We investigate a bit more systematically:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(seed=77)\n",
    "nbx,nby=50,35\n",
    "for i in range(10):\n",
    "    tstmkt = NTU_market(np.random.rand(nbx,nby)-0.3,np.random.rand(nbx,nby)-0.2)\n",
    "    tstmkt.solveGaleShapley()\n",
    "    print('Stable = '+str(tstmkt.is_GS_stable())+'; nb of steps = '+str(tstmkt.comp_nbsteps))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Deferred acceptance revisited: Adachi's algorithm\n",
    "\n",
    "The reference for this section is [A00] and [GL21].\n",
    "\n",
    "**Adachi's algorithm.**\n",
    "\n",
    "Set initially $v^0_y$: say $v^0_y < \\min_{x\\in \\mathcal{X}_0} \\left\\{ \\gamma_{xy} \\right\\}$.\n",
    "\n",
    "Iterate over t:\n",
    "\n",
    "$\\left\\{ \n",
    "\\begin{array}{l}\n",
    "u_{x}^{t+1}=\\max \\left\\{ \\max_{y\\in \\mathcal{Y}}\\left\\{ \\alpha _{xy}:\\gamma\n",
    "_{xy}\\geq v_{y}^{t}\\right\\} ,\\alpha _{x0}\\right\\}  \\\\ \n",
    "v_{y}^{t+1}=\\max \\left\\{ \\max_{x\\in \\mathcal{X}}\\left\\{ \\gamma _{xy}:\\alpha\n",
    "_{xy}\\geq u_{x}^{t+1}\\right\\} ,\\gamma _{0y}\\right\\}\n",
    "\\end{array}\n",
    "\\right. $\n",
    "\n",
    "until $u^{t+1} = u^{t}$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now implement Adachi's algorithm. It will be based on a pair of functions `uo_from_vo` and `vo_from_uo` which implement the maximization problem above, namely:<br>\n",
    "$\\max \\left\\{ \\max_{y\\in \\mathcal{Y}}\\left\\{ \\alpha _{xy}:\\gamma_{xy}\\geq v_{y}\\right\\} ,\\alpha _{x0}\\right\\}$, and<br>\n",
    "$\\max \\left\\{ \\max_{x\\in \\mathcal{X}}\\left\\{ \\gamma _{xy}:\\alpha\n",
    "_{xy}\\geq u_{x}\\right\\} ,\\gamma _{0y}\\right\\}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def uo_from_vo(self,vo_y):\n",
    "    excluded = np.hstack([(self.γo_x_y < vo_y)[0:-1,:],np.zeros((self.nbx,1))])\n",
    "    return(np.ma.masked_array(self.αo_x_y,  excluded ).max(axis = 1).data)\n",
    "NTU_market.uo_from_vo = uo_from_vo\n",
    "\n",
    "def  vo_from_uo(self,uo_x):\n",
    "    excluded = np.vstack([(self.αo_x_y < uo_x.reshape((-1,1)))[:,0:-1],np.zeros((1,self.nby))])\n",
    "    return(np.ma.masked_array(self.γo_x_y,  excluded ).max(axis = 0).data)\n",
    "NTU_market.vo_from_uo = vo_from_uo\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adachi's algorithm iterates the loop above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solveAdachi(self,output=0,trace=False, startu_x = None, startv_y = None ):\n",
    "    self.comp_nbsteps = 0\n",
    "    tracemax = self.nbx*self.nby\n",
    "    if trace:\n",
    "        self.traceuo_x_t = np.zeros((tracemax,self.nbx))    \n",
    "    if startu_x is None:\n",
    "        uo_x = self.largex * np.ones(self.nbx, dtype = 'int64') # x's utilities are highest\n",
    "    else:\n",
    "        uo_x = startu_x\n",
    "    if startv_y is None:\n",
    "        vo_y = self.smally * np.ones(self.nby, dtype = 'int64') # y's utilities are lowest\n",
    "    else:\n",
    "        vo_y = startv_y\n",
    "    while True :\n",
    "        uonew_x = self.uo_from_vo(vo_y) # each x proposes to favorite y among those willing to consider them\n",
    "        if (uonew_x == uo_x).all() :\n",
    "            if trace:\n",
    "                self.traceuo_x_t = self.traceuo_x_t[0:min(self.comp_nbsteps,tracemax),:]\n",
    "            break\n",
    "        uo_x = uonew_x      \n",
    "        vo_y = self.vo_from_uo(uo_x) # each y proposes to favorite x among those willing to consider them\n",
    "        if output >= 2:\n",
    "            μP_x_y = self.μ_from_uo(uo_x)\n",
    "            μE_x_y = self.μ_from_vo(vo_y)\n",
    "            print(\"Round \"+str(self.comp_nbsteps)+\":\\n\",\"μ_P=\\n\" ,μP_x_y,\"\\n μ_E=\\n\", μE_x_y)\n",
    "        if trace and self.comp_nbsteps<tracemax:\n",
    "            self.traceuo_x_t[self.comp_nbsteps,:] = uo_x\n",
    "        self.comp_nbsteps += 1\n",
    "    self.eq_μ_x_y = self.μ_from_vo(vo_y)\n",
    "    return(0)\n",
    "\n",
    "NTU_market.solveAdachi = solveAdachi "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We try it on `running_mkt`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "running_mkt.solveAdachi()\n",
    "running_mkt.is_GS_stable()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we investigate a bit more Adachi's algorithm similarly to what we did above with Gale and Shapley's algorithm, and in comparison with the latter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Times: Gale-Shapley; Adachi\")\n",
    "\n",
    "np.random.seed(seed=77)\n",
    "nbx,nby=50,35\n",
    "for i in range(10):\n",
    "    tstmkt = NTU_market(np.random.rand(nbx,nby)-0.3,np.random.rand(nbx,nby)-0.2)\n",
    "    res=tstmkt.solveGaleShapley()\n",
    "    if not tstmkt.is_GS_stable() :\n",
    "        raise Exception('Output of Gale Shapley is not stable')\n",
    "    i1 = tstmkt.comp_nbsteps\n",
    "    t1 = tstmkt.comp_time\n",
    "    mu1 = tstmkt.eq_μ_x_y\n",
    "    res=tstmkt.solveAdachi()\n",
    "    i2 = tstmkt.comp_nbsteps\n",
    "    t2 = tstmkt.comp_time\n",
    "    mu2 = tstmkt.eq_μ_x_y\n",
    "    if not tstmkt.is_GS_stable() :\n",
    "        raise Exception('Output of Adachi is not stable')\n",
    "    if  np.max(np.abs(mu1 - mu2)) != 0:\n",
    "        raise Exception(seed,': results differ.')\n",
    "    print(\"GS steps = \"+ str(i1)+\"\\t; A steps = \"+str(i2))\n",
    "print('All matching outcomes from Gale-Shapley and Adachi are stable and coincide.')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Testing an alternative representation of Adachi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, we show experimentally that Adachi's algorithm can equivalently be expressed by the following formulation, see [GL21]:\n",
    "\n",
    "**Adachi Bis algorithm.**\n",
    "\n",
    "Set initially $v^0_y$: say $v^0_y < \\min_{x\\in \\mathcal{X}_0} \\left\\{ \\gamma_{xy} \\right\\}$.\n",
    "\n",
    "Iterate over t:\n",
    "\n",
    "$\\left\\{ \n",
    "\\begin{array}{l}\n",
    "u_{x}^{t+1}=\\max \\left\\{ \\max_{y\\in \\mathcal{Y}}\\left\\{ \\alpha _{xy}:\\gamma\n",
    "_{xy}\\geq v_{y}^{t}\\right\\} ,\\alpha _{x0}\\right\\}  \\\\ \n",
    "v_{y}^{t+1}=\\max \\left\\{ \\max_{x\\in \\mathcal{X}}\\left\\{ \\gamma _{xy}:\\alpha\n",
    "_{xy}= u_{x}^{t+1}\\right\\} ,\\gamma _{0y}\\right\\}\n",
    "\\end{array}\n",
    "\\right. $\n",
    "\n",
    "until $u^{t+1} = u^{t}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def  vo_from_uo_bis(self,uo_x):\n",
    "    excluded = np.vstack([(self.αo_x_y != uo_x.reshape((-1,1)))[:,0:-1],np.zeros((1,self.nby))])\n",
    "    return(np.ma.masked_array(self.γo_x_y,  excluded ).max(axis = 0).data)\n",
    "\n",
    "NTU_market.vo_from_uo_bis = vo_from_uo_bis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solveAdachiBis(self,output=0,trace=False):\n",
    "    self.comp_nbsteps = 0\n",
    "    tracemax = self.nbx*self.nby\n",
    "    if trace:\n",
    "        self.traceuo_x_t = np.zeros((tracemax,self.nbx))    \n",
    "    uo_x = self.largex * np.ones(self.nbx, dtype = 'int64') # x's utilities are highest\n",
    "    vo_y = self.smally * np.ones(self.nby, dtype = 'int64') # y's utilities are lowest\n",
    "    while True :\n",
    "        uonew_x = self.uo_from_vo(vo_y) # each x proposes to favorite y among those willing to consider them\n",
    "        if (uonew_x == uo_x).all() :\n",
    "            if trace:\n",
    "                self.traceuo_x_t = self.traceuo_x_t[0:min(self.comp_nbsteps,tracemax),:]\n",
    "            break\n",
    "        uo_x = uonew_x      \n",
    "        vo_y = self.vo_from_uo_bis(uo_x) # each y proposes to favorite x among those willing to consider them\n",
    "        if trace and self.comp_nbsteps<tracemax:\n",
    "            self.traceuo_x_t[self.comp_nbsteps,:] = uo_x\n",
    "        self.comp_nbsteps += 1\n",
    "    self.eq_μ_x_y = self.μ_from_vo(vo_y)\n",
    "    return(0)\n",
    "\n",
    "NTU_market.solveAdachiBis = solveAdachiBis "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(seed=77)\n",
    "nbx,nby=50,35\n",
    "for i in range(10):\n",
    "    mkt_tst = NTU_market(np.random.rand(4,3)-0.2,np.random.rand(4,3)-0.3)\n",
    "    mkt_tst.solveAdachi(trace =True)\n",
    "    t1  = mkt_tst.traceuo_x_t\n",
    "    mkt_tst.solveAdachiBis(trace=True)\n",
    "    t2 = mkt_tst.traceuo_x_t\n",
    "    print(\"Discrepancy between Adachi and Adachi Bis=\", np.max(np.abs(t2-t1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reformulating Adachi and Gale-Shapley coordinate updates algorithm\n",
    "\n",
    "[GL21] reformulate Adachi's and Gale-Shapley algorithms as a variants of the (blockwise) Gauss-Seidel coordinate update algorithm.\n",
    "\n",
    "As usual, change the sign of $u_{x}$ take $p=\\left( u,-v\\right) $ \n",
    "and define<br>\n",
    "$\\left\\{ \n",
    "\\begin{array}{l}\n",
    "Q_{x}\\left( p\\right) =p_{x}-\\max_{y\\in \\mathcal{Y}}\\left\\{ \\alpha\n",
    "_{xy}:p_{y}\\geq -\\gamma _{xy},0\\right\\} ,x\\in \\mathcal{X} \\\\ \n",
    "Q_{y}\\left( p\\right) =p_{y}-\\min_{x\\in \\mathcal{X}}\\left\\{ -\\gamma\n",
    "_{xy}:\\alpha _{xy}\\geq p_{x},0\\right\\} ,y\\in \\mathcal{Y}%\n",
    "\\end{array}%\n",
    "\\right.$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Q_z(self,p_z):\n",
    "    uo_x = p_z[0:self.nbx]\n",
    "    vo_y = - p_z[self.nbx:(self.nbx+self.nby)]\n",
    "    uonew_x = self.uo_from_vo(vo_y)\n",
    "    vonew_y = self.vo_from_uo(uo_x)\n",
    "    return(np.append(uo_x - uonew_x,vonew_y - vo_y))\n",
    "\n",
    "NTU_market.Q_z = Q_z    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adachi as Gauss-Seidel\n",
    "\n",
    "[GL21] reformulate Adachi's algorithm as a (blockwise) Gauss-Seidel algorithm.\n",
    "\n",
    "\n",
    "Start by $p_{x}^{0}=\\max_{y\\in \\mathcal{Y}_{0}}\\left\\{ \\alpha _{xy}\\right\\} $\n",
    "and $p_{y}^{0}=\\max_{x\\in \\mathcal{X}_{0}}\\left\\{ -\\gamma _{xy}\\right\\} $,\n",
    "which are such that $e\\left( p\\right) \\geq 0$.\n",
    "\n",
    "Adachi reinteprets as a blockwise Gauss-Seidel algorithm<br>\n",
    "$\\left\\{ \n",
    "\\begin{array}{l}\n",
    "p_{x}^{t+1}:Q_{x}\\left( p_{x}^{t+1},\\left( p_{y}^{t}\\right) _{y}\\right) =0,\n",
    "\\\\ \n",
    "p_{y}^{t+1}:Q_{y}\\left( p_{y}^{t+1},\\left( p_{x}^{t+1}\\right) _{x}\\right) =0.%\n",
    "\\end{array}%\n",
    "\\right.$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cux_z(self,p_z):\n",
    "    uo_x = p_z[0:self.nbx]\n",
    "    vo_y = - p_z[self.nbx:(self.nbx+self.nby)]\n",
    "    uonew_x = self.uo_from_vo(vo_y)\n",
    "    return (np.append(uonew_x , - vo_y))\n",
    "        \n",
    "NTU_market.cux_z = cux_z\n",
    "\n",
    "\n",
    "def cuy_z(self,p_z):\n",
    "    uo_x = p_z[0:self.nbx]\n",
    "    vo_y = - p_z[self.nbx:(self.nbx+self.nby)]\n",
    "    excluded = np.vstack([(self.αo_x_y < np.repeat([uo_x],self.nby+1,axis = 1).reshape((self.nbx,-1)))[:,0:-1],np.repeat(False,self.nby).reshape((-1,self.nby))])\n",
    "    vonew_y = self.vo_from_uo(uo_x)\n",
    "    return( np.append(uo_x, - vonew_y ) )\n",
    "\n",
    "NTU_market.cuy_z = cuy_z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solveCU(self,output=0,trace=False):\n",
    "    self.comp_nbsteps = 0\n",
    "    tracemax = self.nbx*self.nby\n",
    "    if trace:\n",
    "        self.traceuo_x_t = np.zeros((tracemax,self.nbx))    \n",
    "    p_z = np.append(self.largex * np.ones(self.nbx), - self.smally * np.ones(self.nby) )\n",
    "    while True :\n",
    "        pnew_z = self.cux_z(p_z) \n",
    "        pnew_z = self.cuy_z(pnew_z)\n",
    "        if (pnew_z == p_z).all() :\n",
    "            if trace:\n",
    "                self.traceuo_x_t = self.traceuo_x_t[0:min(self.comp_nbsteps,tracemax),:]\n",
    "            break\n",
    "        p_z = pnew_z      \n",
    "        if trace and self.comp_nbsteps<tracemax:\n",
    "            self.traceuo_x_t[self.comp_nbsteps,:] = p_z[0:self.nbx]\n",
    "        self.comp_nbsteps += 1\n",
    "    self.eq_μ_x_y = self.μ_from_vo(- p_z[self.nbx:(self.nbx+self.nby)])\n",
    "    return(0)\n",
    "\n",
    "NTU_market.solveCU = solveCU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We test numerically that Adachi and Gauss-Seidel coincide:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1000)\n",
    "\n",
    "mkt_tst = NTU_market(np.random.rand(80,120)-0.2,np.random.rand(80,120)-0.3)\n",
    "\n",
    "mkt_tst.solveAdachi(trace =True)\n",
    "t1  = mkt_tst.traceuo_x_t\n",
    "\n",
    "mkt_tst.solveCU(trace=True)\n",
    "t2 = mkt_tst.traceuo_x_t\n",
    "\n",
    "print(\"Discrepancy between Adachi and Gauss-Seidel=\", np.max(t2-t1))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gale-Shapley as damped Gauss-Seidel\n",
    "\n",
    "[GL21] reformulate Gale and Shapley's algorithm as a damped blockwise Gauss-Seidel algorithm.\n",
    "\n",
    "Gale and Shapley's proposal step reformulates as<br>\n",
    "$\\left\\{ \n",
    "\\begin{array}{l}\n",
    "p_{x}^{t+1}=decr_{x}\\left( p_{x}^{t}\\right) \\text{ if }Q_{x}\\left(\n",
    "p_{x}^{t+1},\\left( p_{y}^{t}\\right) _{y}\\right) >0\\text{, } \\\\ \n",
    "p_{x}^{t+1}=incr_{x}\\left( p_{x}^{t}\\right) \\text{ if }Q_{x}\\left(\n",
    "p_{x}^{t+1},\\left( p_{y}^{t}\\right) _{y}\\right) <0 \\\\ \n",
    "p_{x}^{t+1}=p_{x}^{t}\\text{ if }Q_{x}\\left( p_{x}^{t+1},\\left(\n",
    "p_{y}^{t}\\right) _{y}\\right) =0, \\\\ \n",
    "p_{y}^{t+1}:Q_{y}\\left( p_{y}^{t+1},\\left( p_{x}^{t+1}\\right) _{x}\\right) =0%\n",
    "\\end{array}%\n",
    "\\right.$<br> \n",
    "where<br>\n",
    "$decr_{x}\\left( p\\right) =\\max_{y\\in \\mathcal{Y}_{0}}\\left\\{ \\alpha\n",
    "_{xy}:\\alpha _{xy}<p\\right\\} $ is the next value below $p$ (in terms of the $\\alpha _{xy}$'s), and<br> \n",
    "$incr_{x}\\left(p\\right) =\\min_{x\\in \\mathcal{X}_{0}}\\left\\{ \\alpha _{xy}:\\alpha\n",
    "_{xy}>p\\right\\} $ is the next value above $p$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def damped_cux_z(self,p_z):\n",
    "    uo_x = p_z[0:self.nbx]\n",
    "    vo_y = - p_z[self.nbx:(self.nbx+self.nby)]\n",
    "    uonew_x = self.uo_from_vo(vo_y)\n",
    "    return (np.append(uo_x - np.where(uo_x > uonew_x, 1,0)+ np.where(uo_x < uonew_x, 1,0) , - vo_y))\n",
    "        \n",
    "NTU_market.damped_cux_z = damped_cux_z\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solveDampedCU(self,output=0,trace=False):\n",
    "    self.comp_nbsteps = 0\n",
    "    tracemax = self.nbx*self.nby\n",
    "    if trace:\n",
    "        self.traceuo_x_t = np.zeros((tracemax,self.nbx))    \n",
    "    p_z = np.append(self.largex * np.ones(self.nbx), - self.smally * np.ones(self.nby) )\n",
    "    while True :\n",
    "        pnew_z = self.damped_cux_z(p_z)   # each x proposes to favorite y among those willing to consider them:\n",
    "        pnew_z = self.cuy_z(pnew_z)\n",
    "        if (pnew_z == p_z).all() :\n",
    "            if trace:\n",
    "                self.traceuo_x_t = self.traceuo_x_t[0:min(self.comp_nbsteps,tracemax),:]\n",
    "            break\n",
    "        p_z = pnew_z      \n",
    "        # each y proposes to favorite x among those willing to consider them:\n",
    "        if trace and self.comp_nbsteps<tracemax:\n",
    "            self.traceuo_x_t[self.comp_nbsteps,:] = p_z[0:self.nbx]\n",
    "        self.comp_nbsteps += 1\n",
    "    self.eq_μ_x_y = self.μ_from_vo(- p_z[self.nbx:(self.nbx+self.nby)])\n",
    "    return(0)\n",
    "\n",
    "NTU_market.solveDampedCU = solveDampedCU     \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We test numerically that Gale and Shapley and damped Gauss-Seidel coincide:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1000)\n",
    "\n",
    "mkt_tst = NTU_market(np.random.rand(80,120)-0.2,np.random.rand(80,120)-0.3)\n",
    "\n",
    "mkt_tst.solveGaleShapley(trace =True)\n",
    "t1  = mkt_tst.traceuo_x_t\n",
    "\n",
    "mkt_tst.solveDampedCU(trace=True)\n",
    "t2 = mkt_tst.traceuo_x_t\n",
    "\n",
    "print(\"Discrepancy between Gale and Shapley and coordinate update=\", np.max(t2-t1))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing running times\n",
    "\n",
    "We build a wrapper that keeps track of the times.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from time import time\n",
    "\n",
    "def solveDeferredAcceptance(self, algorithm = 'Adachi', output=0, trace=False):\n",
    "    start_time = time()\n",
    "    if algorithm == 'GS' :\n",
    "        self.solveGaleShapley(output,trace)\n",
    "    elif algorithm == 'Adachi' :\n",
    "        self.solveAdachi(output,trace)\n",
    "    elif algorithm == 'DARUM' :\n",
    "        self.solveDARUM(output,trace)   \n",
    "    elif algorithm == 'CU' :\n",
    "        self.solveCU(output,trace)\n",
    "    else:\n",
    "        raise Exception(\"Algorithm \" + algorithm + \" is not implemented.\")\n",
    "    self.comp_time =  time() - start_time\n",
    "    if output >= 1:\n",
    "        print(\"Converged in \",self.comp_nbsteps,\" steps and \", self.comp_time, \" seconds.\")\n",
    "    if output == 1:\n",
    "        print(\"mu_x_y=\",self.eq_μ_x_y) \n",
    "    return (0)\n",
    "\n",
    "    \n",
    "NTU_market.solveDeferredAcceptance = solveDeferredAcceptance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Enumeration of stable matchings\n",
    "\n",
    "The following algorithm, described in [RS] p. 62, enumerates all stable matchings."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we implement a convenient method that returns a list of matched pairs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def matched_pairs(self, μ_x_y = None):\n",
    "    if μ_x_y is None:\n",
    "        μ_x_y = self.eq_μ_x_y\n",
    "    nzx,nzy = np.nonzero(μ_x_y)\n",
    "    return [(nzx[i],nzy[i]) for i in range(len(nzx) )]\n",
    "\n",
    "NTU_market.matched_pairs = matched_pairs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following method returns the matchings (as designed by the payoff vectors) that are immediately above (or below) a given matching."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def next_us(self,uo_x,down = True):\n",
    "    μ_x_y = self.μ_from_uo(uo_x)\n",
    "    vo_y = self.vo_from_uo(uo_x)\n",
    "    if down:\n",
    "        incr= -1\n",
    "    else:\n",
    "        incr = 1\n",
    "    g = nx.DiGraph()\n",
    "    xs = list(range(self.nbx))\n",
    "    ys = list(range(self.nbx,self.nby))\n",
    "    g.add_nodes_from(xs+ys)\n",
    "    g.add_edges_from([(x,self.nbx+y) for x,y in self.matched_pairs(μ_x_y)])\n",
    "    #μprime_x_y = np.where( (self.αo_x_y.T == uo_x + incr).T , 1 , 0 )[:,0:-1]\n",
    "    μprime_x_y = np.where( (self.αo_x_y[:,0:-1].T == uo_x + incr).T &  (self.γo_x_y[0:-1,:]*incr <= vo_y * incr), 1 , 0 ) \n",
    "    compatible_pairs = self.matched_pairs(μprime_x_y)\n",
    "    g.add_edges_from([(self.nbx+j,i) for (i,j) in compatible_pairs])\n",
    "    cycles = list(nx.simple_cycles(g))\n",
    "    newus = []\n",
    "    for c in cycles:\n",
    "        xs_to_update = [x for x in c if (x <self.nbx)]\n",
    "        nextu_x = uo_x.copy()\n",
    "        nextu_x[xs_to_update] += incr\n",
    "        newus.append(nextu_x)\n",
    "    return(newus)\n",
    "\n",
    "NTU_market.next_us = next_us"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's run an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rs_ex_2_17.solveGaleShapley()\n",
    "uo_x = np.sum(rs_ex_2_17.αo_x_y [:,0:-1] * rs_ex_2_17.eq_μ_x_y,axis = 1)\n",
    "print(rs_ex_2_17.next_us(uo_x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following method lists all stable matchings, starting with the lattice upper bound for side $\\mathcal{X}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def enumerate_us(self):\n",
    "    self.solveGaleShapley()\n",
    "    uo_x = np.sum(self.αo_x_y [:,0:-1] * self.eq_μ_x_y,axis = 1)\n",
    "    us_list = [uo_x.tolist()]\n",
    "    i = 0\n",
    "    while (i < len(us_list)):\n",
    "        inds_alreadyhere = []\n",
    "        us_toadd = [arr.tolist() for arr in self.next_us(np.array(us_list[i],dtype='int64'))]\n",
    "        for u_x in us_toadd:\n",
    "            if u_x in us_list:\n",
    "                inds_alreadyhere.append(us_list.index(u_x))\n",
    "                us_toadd.remove(u_x)\n",
    "        \n",
    "        inds_toadd = list(range(len(us_list),len(us_list)+len(us_toadd) ))\n",
    "        us_list = us_list + us_toadd \n",
    "        i += 1\n",
    "    return(np.array([u for u in us_list],dtype = 'int64') )\n",
    "\n",
    "NTU_market.enumerate_us = enumerate_us"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "the_u_s = rs_ex_2_17.enumerate_us()\n",
    "for u_x in the_u_s:\n",
    "    print('u_x = '+str(u_x)+ ' ; stable = '+str(rs_ex_2_17.is_GS_stable(μ_x_y)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aligned preferences\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following [NY09], we define aligned preferences as the situation when $\\alpha_{xy}=\\gamma_{xy}=\\varphi_{xy}$. \n",
    "\n",
    "We setup a derived class for this called `Aligned_NTU_market`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Aligned_NTU_market(NTU_market):\n",
    "    def __init__(self,ϕ_x_y):\n",
    "        NTU_market.__init__(self,α_x_y = ϕ_x_y, γ_x_y = ϕ_x_y)\n",
    "        self.ϕ_x_y = ϕ_x_y\n",
    "        self.n,_ =  ϕ_x_y.shape\n",
    "        ϕo_x_y = np.zeros(self.n* self.n)\n",
    "        ϕo_x_y[ϕ_x_y.flatten().argsort()] = 1+np.arange(self.n*self.n)\n",
    "        self.ϕo_x_y=ϕo_x_y.reshape((self.n,self.n))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `MaxMaxLex` algorithm was proposed in [GGH21] to look for a stable matching in this type of markets. It looks for the pair $xy$ that maximizes $\\varphi_{xy}$, matches $x$ and $y$, and removes them from the available individuals. We implement into: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solveMaxMaxLex(self):\n",
    "    theϕo_x_y = self.ϕo_x_y.copy()\n",
    "    self.eq_μ_x_y = np.zeros((self.n,self.n), dtype = 'int64')\n",
    "    for k in range(self.n):\n",
    "        x,y = np.unravel_index(theϕo_x_y.argmax(),(self.n,self.n))\n",
    "        self.eq_μ_x_y[x,y]=1\n",
    "        theϕo_x_y[x,:]=0\n",
    "        theϕo_x_y[:,y]=0\n",
    "    return 0\n",
    "    \n",
    "Aligned_NTU_market.solveMaxMaxLex = solveMaxMaxLex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compare the `MaxMaxLex` algorithm with Adachi:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(seed=77)\n",
    "aligned_prefs_ex = Aligned_NTU_market(np.random.rand(20,20))\n",
    "\n",
    "aligned_prefs_ex.solveMaxMaxLex()  \n",
    "print('MaxMaxLex:',aligned_prefs_ex.uo_from_μ(),'; stable = '+str(aligned_prefs_ex.is_GS_stable()))\n",
    "aligned_prefs_ex.solveAdachi()\n",
    "print('Adachi   :',aligned_prefs_ex.uo_from_μ(),'; stable = '+str(aligned_prefs_ex.is_GS_stable()))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Aggregate stable matchings\n",
    "\n",
    "The reference for this section are [GH19] and [GHS21].\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aggregate stability\n",
    "\n",
    "**We now allow for the possibility that there are more than one individual per type:**<br>\n",
    "$n_x$ and $m_y$ don't need to be equal to one for all $x\\in\\mathcal{X}$ and $y\\in\\mathcal{Y}$.<br>\n",
    "\n",
    "As in the previous lecture, $\\mu_{xy}$ is the mass of $xy$ pairs at equilibrium,<br>\n",
    "and $u_x$ and $v_y$ are the payoffs of worker $x$ and firm $y$ at equilibrium.\n",
    "\n",
    "**Definition.** $(\\mu,u,v)$ is an aggregate stable outcome if:\n",
    "\n",
    "(i) Population constraints are satisfied<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.$<br>\n",
    "    \n",
    "(ii) There is no blocking pair and all individuals are rational:<br>\n",
    "$ \\max \\{ u_x - \\alpha_{xy},v_y - \\gamma_{xy} \\} \\geq 0$, and<br>\n",
    "$u_x \\geq 0$ and $v_y \\geq 0$\n",
    "\n",
    " \n",
    "(iii) Strong complementarity:<br>\n",
    "$\\mu_{xy}>0$ implies $ \\max \\{ u_x - \\alpha_{xy},v_y - \\gamma_{xy} \\} = 0$<br>\n",
    "$\\mu_{x0}>0$ implies $u_x = 0$<br>\n",
    "$\\mu_{0y}>0$ implies $v_y = 0$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We implement the *distance-to-frontier function* in this NTU matching setting as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def DTF(self ,u_x,v_y):\n",
    "    return np.maximum(u_x.reshape((-1,1)) - self.α_x_y[:,0:-1], v_y - self.γ_x_y[0:-1,:])\n",
    "\n",
    "NTU_market.DTF = DTF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function detects aggregate stable matchings:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_aggregate_stable(self, μ_x_y =None ,u_x=None,v_y =None, output=0 ):\n",
    "    if μ_x_y is None:\n",
    "        μ_x_y = self.eq_μ_x_y\n",
    "    if u_x is None:\n",
    "        u_x = self.eq_u_x\n",
    "    if v_y is None:\n",
    "        v_y = self.eq_v_y\n",
    "    if (np.min(μ_x_y) < 0) or (np.min(self.n_x-np.sum(μ_x_y,axis = 1))<0) or (np.min(self.m_y-np.sum(μ_x_y,axis = 0))<0) :\n",
    "        if output > 0 :\n",
    "            print('The μ is not feasible.')\n",
    "        return(False)\n",
    "    if (np.min(u_x)<0) or  (np.min(v_y)<0):\n",
    "        if output > 0 :\n",
    "            print('u_x < 0 for some x or v_y < 0 for some y.')\n",
    "        return(False)\n",
    "    D_x_y = self.DTF(u_x,v_y)\n",
    "    if np.min(D_x_y)<0 :\n",
    "        if output > 0 :\n",
    "            print('There is a blocking pair.')\n",
    "        return(False)\n",
    "    if np.sum(μ_x_y * D_x_y)>0 :\n",
    "        if output > 0 :\n",
    "            print('Complementary slackness does not hold.')\n",
    "        return(False)\n",
    "    if output > 0 :\n",
    "        print('The matching is stable.')\n",
    "    return (True)\n",
    "\n",
    "NTU_market.is_aggregate_stable = is_aggregate_stable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A first example\n",
    "\n",
    "Consider the following example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_ex=NTU_market(np.array([[.2,.6],[.3,-.1],[.7,.5]]),np.array([[.7,.3],[.3,.6],[0.2,-.2]]) )\n",
    "print('α_x_y=\\n',  new_ex.α_x_y[:,0:-1],'\\nγ_x_y=\\n',new_ex.γ_x_y[0:-1,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve this using Gale and Shapley into:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_ex.solveGaleShapley()\n",
    "theu,thev = new_ex.u_from_μ(),new_ex.v_from_μ()\n",
    "print('\\nu_x=',theu,'and v_y=',thev,'\\n\\nμ_x_y=\\n',new_ex.eq_μ_x_y,'\\n\\n D_x_y=\\n',new_ex.DTF(theu,thev))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is interesting to note that while u_x= [0.6 0.3 0. ] and v_y= [0.3 0.3] are stable, one can decrease utility of some agents and remain stable (in the aggregate sense). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theu,thev = new_ex.u_from_μ(),new_ex.v_from_μ()\n",
    "print(\"Is stable with full u and v :\",new_ex.is_aggregate_stable(new_ex.eq_μ_x_y,theu,thev))\n",
    "theu[0] = theu[0] - 0.3\n",
    "theu[1] = theu[1] - 0.2\n",
    "print(\"Is stable with altered u and v :\",new_ex.is_aggregate_stable(new_ex.eq_μ_x_y,theu,thev))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A second example: one car, two passengers\n",
    "\n",
    "In the next example we have one car and two passengers of the same type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "one_car_ex = NTU_market(np.array([[1]]),np.array([[1]]),n_x = np.array([1]), m_y = np.array([2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sensible solution consists in matching the car to one of the passengers, and leave the other passenger unmatched. However, it is not a stable solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "one_car_ex.is_GS_stable(np.array([[1]]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, the above solution is aggregate stable if $u=1$ for the car and $v=0$ for *both* passengers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "one_car_ex.is_aggregate_stable(np.array([[1]]),np.array([1]),np.array([0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Deferred acceptance with aggregate stable matchings\n",
    "\n",
    "The DARUM algorithm is an aggregate version of deferred acceptance. It is based on the specification of *aggregate choice functions* $C_{\\mathcal{X}}$ and $C_{\\mathcal{Y}}$ on each sides of the market:\n",
    "\n",
    "* $\\mu \\in C_{\\mathcal{X}}\\left( \\overline{\\mu }\\right) $ means that $\\mu $ is preferred (under $x$'s preferences) over all $\\tilde{\\mu%\n",
    "}$ such that<br>\n",
    "    $\\tilde{\\mu}\\leq \\overline{\\mu }$ and $\\sum_{y}\\tilde{\\mu}_{xy}\\leq\n",
    "n_{x}\\forall x\\in \\mathcal{X}$.\n",
    "\n",
    "* $\\mu \\in C_{\\mathcal{Y}}\\left( \\overline{\\mu }\\right) $ means that $\\mu $ is preferred (under $y$'s preferences) over all $\\tilde{\\mu%\n",
    "}$ such that<br>\n",
    "    $\\tilde{\\mu}\\leq \\overline{\\mu }$ and $\\sum_{x}\\tilde{\\mu}_{xy}\\leq m_{y}\\forall y\\in \\mathcal{Y}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The DARUM algorithm is then:\n",
    "\n",
    "**Initialization**: $\\mu _{xy}^{A,0}=\\min \\left( n_{x},m_{y}\\right) $\n",
    "\n",
    "**Repeat**:<br>\n",
    "* Proposal phase: $\\mu ^{P,t}\\in C_{\\mathcal{X}}\\left( \\mu ^{A,t}\\right) $<br>\n",
    "* Disposal phase: $\\mu ^{D,t}\\in C_{\\mathcal{Y}}\\left( \\mu ^{P,t}\\right) $<br>\n",
    "* Adjustment phase: $\\mu ^{A,t+1}=\\mu ^{A,t}-\\left( \\mu ^{P,t}-\\mu ^{D,t}\\right)$<br>\n",
    "\n",
    "**Until** $|  \\mu ^{P,t}-\\mu ^{D,t} |$  becomes below tolerance.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We implement the DARUM algorithm generically as follows -- we'll have to specify our aggregate choice functions later:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solveDARUM(self,het1 = 'none',het2 = 'none',output=0,trace=False,tol = 1e-5):\n",
    "    self.comp_nbsteps = 0\n",
    "    tracemax = 100*self.nbx*self.nby\n",
    "    if (output>=2):\n",
    "        print(\"Offers made and kept are denoted +1; offers made and rejected are denoted -1.\")\n",
    "    if trace:\n",
    "        self.traceu_x_t = np.zeros((tracemax,self.nbx))\n",
    "        self.traceuo_x_t = np.zeros((tracemax,self.nbx))    \n",
    "\n",
    "    μA_x_y = np.array([[min(self.n_x[x],self.m_y[y]) for y in range(self.nby)] for x in range(self.nbx)]) \n",
    "    # initially all offers are non rejected\n",
    "    while True :\n",
    "        μP_x_y,self.eq_u_x = self.aggregateChoice(0 , μA_x_y,het1) # the x's pick their preferred offers among those not rejected\n",
    "        μE_y_x,self.eq_v_y = self.aggregateChoice(1, μP_x_y[:,0:-1].T,het2)\n",
    "        μE_x_y = μE_y_x.T # the y's pick their preferred offers among those made\n",
    "        rej_x_y = μP_x_y[:,0:-1] - μE_x_y[0:-1,:] # compute rejected offers\n",
    "        μA_x_y = μA_x_y - rej_x_y  # offers from x that have been rejected are no longer available to x\n",
    "        if output >= 2:\n",
    "            print(\"Round \"+str(self.comp_nbsteps)+\":\\n\" , ( (2*μE_x_y[0:-1,:]-1) * μP_x_y[:,0:-1]))            \n",
    "        if trace and self.comp_nbsteps < tracemax:\n",
    "                self.traceu_x_t[self.comp_nbsteps,:] = self.eq_u_x\n",
    "                for x in range(self.nbx):\n",
    "                    self.traceuo_x_t[self.comp_nbsteps,x] = np.sum(self.αo_x_y [x,:] * μP_x_y[x,:])\n",
    "        self.comp_nbsteps +=1\n",
    "        if np.max(np.abs(rej_x_y)) < tol: \n",
    "            if trace:\n",
    "                self.traceu_x_t = self.traceu_x_t[0:min(self.comp_nbsteps,tracemax),:]\n",
    "                self.traceuo_x_t = self.traceuo_x_t[0:min(self.comp_nbsteps,tracemax),:]\n",
    "            break # if all offers have been accepted (within tolerange), then algorithm stops\n",
    "    self.eq_μ_x_y = μE_x_y[0:-1,:]\n",
    "    return (0)\n",
    "\n",
    "NTU_market.solveDARUM = solveDARUM "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `aggregateChoice` above has not been defined yet; it will be a wrapper function to various aggregate choice functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def aggregateChoice(self,axis,μbar_x_y,heterogeneity = 'none'):\n",
    "    if heterogeneity == 'none':\n",
    "        μ_x_y,u_x = aggregateChoice_noHet(self,axis,μbar_x_y)\n",
    "    elif heterogeneity == 'logit':\n",
    "        μ_x_y,u_x = aggregateChoice_logit(self,axis,μbar_x_y)\n",
    "    else:\n",
    "        raise Exception(\"Heterogeneity \" + heterogeneity + \" is not supported.\")\n",
    "    return (μ_x_y,u_x)\n",
    "\n",
    "NTU_market.aggregateChoice = aggregateChoice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We shall implement the case with no heterogeneity (`heterogeneity = 'none'`) first, and the logit case later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Case with no heterogeneity\n",
    "\n",
    "We have in that case:<br>\n",
    "$\\begin{array}{l}\n",
    "C_{\\mathcal{X}}\\left( \\overline{\\mu }\\right) =\\arg \\max_{\\mu \\geq 0}\n",
    "&\\left\\{ \\sum_{xy}\\mu _{xy}\\alpha _{xy}\\right\\}  \\\\\n",
    "s.t.~ &\\sum_{y\\in \\mathcal{Y}}\\mu _{xy}\\leq n_{x} \\\\\n",
    "~ &\\mu _{xy}\\leq \\overline{\\mu }_{xy}\n",
    "\\end{array}$<br>\n",
    "and<br>\n",
    "$\\begin{array}{l}\n",
    "C_{\\mathcal{y}}\\left( \\overline{\\mu }\\right) =\\arg \\max_{\\mu \\geq 0}\n",
    "&\\left\\{ \\sum_{xy}\\mu _{xy}\\gamma _{xy}\\right\\}  \\\\\n",
    "s.t.~ &\\sum_{x\\in \\mathcal{X}}\\mu _{xy}\\leq m_{y} \\\\\n",
    "~ &\\mu _{xy}\\leq \\overline{\\mu }_{xy}\n",
    "\\end{array}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We build this into:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def aggregateChoice_noHet(self,axis,μbar_x_y):\n",
    "    if axis == 0 : # if proposing side = x\n",
    "        n_x = self.n_x\n",
    "        nbx, nby = self.nbx, self.nby\n",
    "        prefs_x_y,α_x_y = self.prefslistα_x_y, self.α_x_y\n",
    "    else:\n",
    "        n_x = self.m_y\n",
    "        nbx, nby = self.nby, self.nbx\n",
    "        prefs_x_y, α_x_y = self.prefslistγ_x_y.T, self.γ_x_y.T\n",
    "    μ_x_y = np.zeros((nbx,nby+1))\n",
    "    u_x = np.zeros(nbx)\n",
    "    for x in range(nbx):\n",
    "        nxres = n_x[x]\n",
    "        for yind in prefs_x_y[x,]:\n",
    "            if yind == nby:\n",
    "                μ_x_y[x,yind] = nxres\n",
    "                break\n",
    "            if μbar_x_y[x,yind] > 0:\n",
    "                μ_x_y[x,yind] = min(nxres , μbar_x_y[x,yind])\n",
    "                nxres -= μ_x_y[x,yind]\n",
    "            if nxres == 0:\n",
    "                break\n",
    "        u_x[x]=α_x_y[x,yind]\n",
    "    return(μ_x_y,u_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We verify that in the case when there is one agent per type, DARUM coincides with Gale and Shapley at every step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(seed=1000)\n",
    "mkt_tst = NTU_market(np.random.rand(2,3)-0.2,np.random.rand(2,3)-0.3)\n",
    "mkt_tst.solveDARUM(trace=True,output=0)\n",
    "Darum_is_GS_stable = mkt_tst.is_GS_stable()\n",
    "trace_Darum = mkt_tst.traceuo_x_t\n",
    "mkt_tst.solveGaleShapley(trace=True,output=0)\n",
    "GaleShapley_is_GS_stable = mkt_tst.is_GS_stable()\n",
    "trace_GaleShapley = mkt_tst.traceuo_x_t\n",
    "\n",
    "print(\"Both DARUM and Gale-Shapley outputs are stable=\", (Darum_is_GS_stable & GaleShapley_is_GS_stable))\n",
    "print(\"Discrepancy between DARUM and Gale-Shapley =\", np.max(np.abs(trace_GaleShapley-trace_Darum)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "one_car_ex.solveDARUM()\n",
    "one_car_ex.eq_μ_x_y,one_car_ex.eq_u_x,one_car_ex.eq_v_y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aggregate stable matching with heterogeneity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logit case\n",
    "\n",
    "We next implement the logit case, in which:<br>\n",
    "$\\begin{array}{l}\n",
    "C_{\\mathcal{X}}\\left( \\overline{\\mu }\\right) =\\arg \\max_{\\mu \\geq 0}\n",
    "&\\left\\{ \\sum_{xy}\\mu _{xy}\\alpha _{xy} - \\sum_{xy}\\mu _{xy} \\log \\mu _{xy} \\right\\}  \\\\\n",
    "s.t.~ &\\sum_{y\\in \\mathcal{Y}}\\mu _{xy}\\leq n_{x} \\\\\n",
    "~ &\\mu _{xy}\\leq \\overline{\\mu }_{xy}\n",
    "\\end{array}$\n",
    "\n",
    "The solution to the above problem has\n",
    "\n",
    "$\\mu _{xy}=\\min \\left\\{ \\mu _{x0}e^{\\alpha _{xy}},\\overline{\\mu }%\n",
    "_{xy}\\right\\} $\n",
    "\n",
    "where $\\mu _{x0}$ solves\n",
    "\n",
    "$\\mu _{x0}+\\sum_{y\\in \\mathcal{Y}}\\min \\left\\{ \\mu _{x0}e^{\\alpha _{xy}},%\n",
    "\\overline{\\mu }_{xy}\\right\\} =n_{x}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We build this into:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def aggregateChoice_logit(self,axis,μbar_x_y):\n",
    "    if axis == 0 : # if proposing side = x\n",
    "        n_x = self.n_x\n",
    "        nbx = self.nbx\n",
    "        nby = self.nby\n",
    "        α_x_y = self.α_x_y\n",
    "    else:\n",
    "        n_x = self.m_y\n",
    "        nbx = self.nby\n",
    "        nby = self.nbx\n",
    "        α_x_y = self.γ_x_y.T    \n",
    "    μ_x_y = np.zeros((nbx,nby+1))\n",
    "    u_x = np.zeros(nbx)\n",
    "    for x in range(nbx):\n",
    "        thesolμ = opt.brentq(lambda theμ : (theμ+ np.minimum(theμ*np.exp(α_x_y[x,0:-1]),μbar_x_y[x,:]).sum() - n_x[x]) ,0,n_x[x])\n",
    "        μ_x_y[x,0:-1] = np.minimum(thesolμ*np.exp(α_x_y[x,0:-1]),μbar_x_y[x,:])\n",
    "        μ_x_y[x,nby] = thesolμ\n",
    "        u_x[x] = - np.log(thesolμ / n_x[x])\n",
    "    return(μ_x_y,u_x)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run it on the one car example, for which we have:<br>\n",
    "$\\left\\{ \\begin{array}{l}\n",
    "\\mu_{x0}+\\min\\{\\mu_{x0},\\mu_{0y} \\} e = 1 \\\\\n",
    "\\mu_{0y}+\\min\\{\\mu_{x0},\\mu_{0y} \\} e = 2\n",
    "\\end{array}\\right.\n",
    "$<br>\n",
    "thus $\\mu_{0y}=\\mu_{x0}+1$,  so<br>\n",
    "$\\mu_{x0} = (1+e)^{-1}$, $\\mu_{0y} = (2+e)/(1+e)$<br>\n",
    "and $u_x = \\log(1+e)$ and $v_y =\\log(1+e) - \\log(1+e/2)$.    \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "one_car_ex.solveDARUM(het1='logit',het2='logit')\n",
    "print('Solution from Darum: μ=',one_car_ex.eq_μ_x_y,', u=', one_car_ex.eq_u_x,', v=',one_car_ex.eq_v_y)\n",
    "print('Predicted solution: μ=',np.exp(1)/(1+np.exp(1)),', u=', np.log(1+np.exp(1)),', v=',np.log(1+np.exp(1))-np.log(1+np.exp(1)/2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution using IPFP\n",
    "\n",
    "As shown in the previous lecture, one can solve the problem using IPFP = Gauss-Seidel algorithm, as done in [GKW19].\n",
    "\n",
    "This consists of iterately:<br>\n",
    "a. solving in $\\mu_{x0}$ the equation<br>\n",
    "$\\mu_{x0}+\\sum_y \\min \\left\\{ \\mu _{x0}e^{\\alpha _{xy}}, \\mu _{0y}e^{\\gamma _{xy}} \\right\\}=n_x,$<br>\n",
    "and:<br>\n",
    "b. solving in $\\mu_{0y}$ the equation<br>\n",
    "$\\mu_{0y}+\\sum_x \\min \\left\\{ \\mu _{x0}e^{\\alpha _{xy}}, \\mu _{0y}e^{\\gamma _{xy}} \\right\\}=m_y.$<br>\n",
    "until the update become below tolerance.\n",
    "\n",
    "\n",
    "The solution is given by<br>\n",
    "$\\mu_{xy} = \\min \\left\\{ \\mu _{x0}e^{\\alpha _{xy}}, \\mu _{0y}e^{\\gamma _{xy}} \\right\\}.$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We build this into:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solveIPFP(self,output=0,trace=False,tol = 1e-5):\n",
    "    self.comp_nbsteps = 0\n",
    "    tracemax = 100*self.nbx*self.nby\n",
    "    if trace:\n",
    "        self.traceu_x_t = np.zeros((tracemax,self.nbx))\n",
    "    μ_x0 = np.zeros(self.nbx)\n",
    "    μ_0y = self.m_y.astype(np.float)\n",
    "    while True :\n",
    "        for x in range(self.nbx):\n",
    "            μ_x0[x] = opt.brentq (lambda theμ : (theμ+ np.minimum(theμ*np.exp(self.α_x_y[x,0:-1]),μ_0y*np.exp(self.γ_x_y[x,:])).sum() - self.n_x[x]) ,0,1.1*self.n_x[x])\n",
    "        μP_x_y = np.minimum(μ_x0.reshape((-1,1))*np.exp(self.α_x_y[:,0:-1]),μ_0y*np.exp(self.γ_x_y[0:-1,:]))        \n",
    "        self.eq_u_x = -np.log(μ_x0 / self.n_x)\n",
    "        for y in range(self.nby):\n",
    "            μ_0y[y] = opt.brentq(lambda theμ : (theμ+ np.minimum(μ_x0*np.exp(self.α_x_y[:,y]),theμ*np.exp(self.γ_x_y[0:-1,y])).sum() - self.m_y[y]) ,0,1.1*self.m_y[y])\n",
    "        μE_x_y = np.minimum(μ_x0.reshape((-1,1))*np.exp(self.α_x_y[:,0:-1]),μ_0y*np.exp(self.γ_x_y[0:-1,:]))        \n",
    "        self.eq_v_y = -np.log(μ_0y / self.m_y)\n",
    "        rej_x_y = μP_x_y - μE_x_y # compute rejected offers\n",
    "        if output >= 2:\n",
    "            print('μP_x_y=\\n',μP_x_y)            \n",
    "            print('μE_x_y=\\n',μE_x_y)            \n",
    "        if trace and self.comp_nbsteps < tracemax:\n",
    "            self.traceu_x_t[self.comp_nbsteps,:] = self.eq_u_x\n",
    "        self.comp_nbsteps +=1\n",
    "        if np.max(np.abs(rej_x_y)) < tol: \n",
    "            if trace:\n",
    "                self.traceu_x_t = self.traceu_x_t[0:min(self.comp_nbsteps,tracemax),:]\n",
    "            break # if all offers have been accepted (within tolerance), then algorithm stops\n",
    "    self.eq_μ_x_y = μE_x_y\n",
    "    return (0)\n",
    "\n",
    "NTU_market.solveIPFP = solveIPFP "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We run the IPFP on the one-car example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "one_car_ex.solveIPFP()\n",
    "print('Solution from IPFP: μ=',one_car_ex.eq_μ_x_y,', u=', one_car_ex.eq_u_x,', v=',one_car_ex.eq_v_y)\n",
    "print('Predicted solution: μ=',np.exp(1)/(1+np.exp(1)),', u=', np.log(1+np.exp(1)),', v=',np.log(1+np.exp(1))-np.log(1+np.exp(1)/2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the `running_mkt` example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "running_mkt.solveDARUM(het1='logit',het2='logit',trace=True)\n",
    "traceD = running_mkt.traceu_x_t\n",
    "running_mkt.eq_μ_x_y,running_mkt.eq_u_x,running_mkt.eq_v_y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "running_mkt.solveIPFP()\n",
    "traceI = running_mkt.traceu_x_t\n",
    "running_mkt.eq_μ_x_y,running_mkt.eq_u_x,running_mkt.eq_v_y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interestingly, in this example, the steps taken by both algorithms are exactly the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.min(np.abs(traceD-traceI))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, this is far from being the case in general. This is not surprising given that IPFP is a Gauss-Seidel algorithm, while DARUM is a deferred acceptance algorithm. See for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(77)\n",
    "\n",
    "nbx,nby=3,2\n",
    "mkt_tst = NTU_market(np.random.rand(nbx,nby)-0.2,np.random.rand(nbx,nby)-0.1)\n",
    "\n",
    "mkt_tst.solveIPFP(trace =True,output=0)\n",
    "t1  = mkt_tst.traceu_x_t\n",
    "μ1 = mkt_tst.eq_μ_x_y\n",
    "u1 = mkt_tst.eq_u_x\n",
    "n1 = mkt_tst.comp_nbsteps\n",
    "\n",
    "mkt_tst.solveDARUM(het1='logit',het2='logit',trace=True)\n",
    "t2 = mkt_tst.traceu_x_t\n",
    "μ2 = mkt_tst.eq_μ_x_y\n",
    "u2 = mkt_tst.eq_u_x\n",
    "n2 = mkt_tst.comp_nbsteps\n",
    "\n",
    "print('Discrepancy betweeo outcomes of IPFP and DARUM = '+str(np.sum(np.abs(μ1.flatten()-μ2.flatten()))))\n",
    "\n",
    "print('IPFP took '+str(n1)+' steps, while DARUM tool '+str(22)+' steps.')\n",
    "\n",
    "print('trace IPFP\\n',t1)\n",
    "print('trace Darum\\n',t2)"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}