et05_matching-with-contracts.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Matching with contracts</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": [
    "\n",
    "# References\n",
    "\n",
    "\n",
    "## Papers\n",
    "\n",
    "\n",
    "* [A00] Hiroyuki Adachi (2000). \"On a characterization of stable matchings.\" *Economics Letters* 68 pp. 43–49.\n",
    "\n",
    "* KELSO-CRAWFORD\n",
    "\n",
    "* HATFIELD-MIGROM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Setting up the model\n",
    "\n",
    "## Populations\n",
    "\n",
    "XXX"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Generic algorithm\n",
    "\n",
    "Let's first load the libraries we will need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create the `Contract_market` class to encompass that information and store it in a convenient way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Contract_market:\n",
    "    def __init__(self,α_x_y_w,γ_x_y_w,n_x = None, m_y = None):\n",
    "        self.nbx,self.nby,self.nbw = α_x_y_w.shape[0],α_x_y_w.shape[1],α_x_y_w.shape[2]\n",
    "        self.α_x_y_w = α_x_y_w\n",
    "        self.γ_x_y_w = γ_x_y_w\n",
    "        if n_x is None:\n",
    "            n_x = np.ones(self.nbx)\n",
    "        if m_y is None:\n",
    "            m_y = np.ones(self.nby) \n",
    "        self.n_x,self.m_y = n_x, m_y\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CONVENTION: $\\alpha_{00w} = \\gamma_{00w} = +\\infty$ so that these are always available.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([20, 10, 10, 10])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tst = np.ones(3)\n",
    "np.where(np.append(np.NaN,tst)==1,10,20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def μ_from_u(self, u_x):\n",
    "    return np.where( (self.α_x_y_w[1:,:,:] == u_x.reshape((-1,1,1))) , 1 , 0 )\n",
    "\n",
    "Contract_market.μ_from_u = μ_from_u\n",
    "\n",
    "def μ_from_v(self, v_y):\n",
    "    return np.where( self.γ_x_y_y[:,1:,:] == v_y.reshape((1,-1,1)) , 1 , 0 )\n",
    "\n",
    "Contract_market.μ_from_v = μ_from_v\n",
    "\n",
    "def u_from_μ(self,μ_x_y_w):\n",
    "    return np.sum(μ_x_y_w * self.α_x_y_w[1:,:,:], axis =(1,2))\n",
    "\n",
    "Contract_market.u_from_μ = u_from_μ  \n",
    "\n",
    "def v_from_μ(self,μ_x_y_w):\n",
    "    return np.sum(μ_x_y_w * self.γ_x_y_w[:,1:,:], axis =(0,2))\n",
    "\n",
    "Contract_market.v_from_μ = v_from_μ  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "def uv_is_stable(self,u_x , v_y, output = 0):\n",
    "    blocking = ( self.α_x_y_w > u_x.reshape((-1,1,1)) ) & ( self.γ_x_y_w > v_y.reshape((1,-1,1)))\n",
    "    #print (blocking)\n",
    "    if np.any (blocking ):\n",
    "        if (output>0):\n",
    "            print(\"A blocking pair is (x\" + str(np.nonzero(blocking)[0][0]) + \",y\" + str(np.nonzero(blocking)[1][0]) +\",w\" + str(np.nonzero(blocking)[2][0]) +\").\" )\n",
    "        return False\n",
    "    else:\n",
    "        return True\n",
    "        \n",
    "Contract_market.uv_is_stable = uv_is_stable\n",
    "\n",
    "def μ_is_stable(self,μ_x_y = None, output=0 ):\n",
    "    if (μ_x_y is None):\n",
    "        μ_x_y = self.eq_μ_x_y\n",
    "    return(self.uv_is_stable(self.u_from_μ(μ_x_y),self.v_from_μ(μ_x_y),output))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "nbx,nby,nbw = 4,3,2\n",
    "np.random.seed(seed=1000)\n",
    "α_x_y = np.random.rand(nbx,nby)\n",
    "γ_x_y = np.random.rand(nbx,nby)\n",
    "w_z = np.linspace(0,1, num = nbw)\n",
    "w_z[np.newaxis,np.newaxis,:].shape\n",
    "\n",
    "tst_mkt = Contract_market(α_x_y[:,:,np.newaxis]+w_z[np.newaxis,np.newaxis,:],γ_x_y[:,:,np.newaxis] - w_z[np.newaxis,np.newaxis,:])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tst_mkt.uv_is_stable(np.zeros(4),np.zeros(3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "self = tst_mkt\n",
    "μA_x_y_w = np.ones((self.nbx, self.nby,self.nbw), dtype = 'int64')\n",
    "μP_x_y_w = np.zeros((self.nbx, self.nby,self.nbw), dtype = 'int64')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "μA_x_y_w = np.zeros((self.nbx, self.nby,self.nbw), dtype = 'int64')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'props_x' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-8-58b7bb488028>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[1;33m[\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mp\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mp\u001b[0m \u001b[1;32min\u001b[0m \u001b[0menumerate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mprops_x\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'props_x' is not defined"
     ]
    }
   ],
   "source": [
    "[(x,p) for x,p in enumerate(props_x)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'props_x' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-9-4fcddc9115e7>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mα_x_y_w\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnbx\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnbx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mprops_x\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'props_x' is not defined"
     ]
    }
   ],
   "source": [
    "self.α_x_y_w.reshape((self.nbx,-1))[range(self.nbx),props_x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "μP_x_y_w.reshape((self.nbx,-1))[range(self.nbx),props_x]=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'u_x' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-10-14c4e862a5d0>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mμP_x_y_w\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mμ_from_u\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mu_x\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'u_x' is not defined"
     ]
    }
   ],
   "source": [
    "μP_x_y_w = self.μ_from_u(u_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[0, 0],\n",
       "        [0, 0],\n",
       "        [0, 1]],\n",
       "\n",
       "       [[0, 0],\n",
       "        [0, 1],\n",
       "        [0, 0]],\n",
       "\n",
       "       [[0, 0],\n",
       "        [0, 1],\n",
       "        [0, 0]],\n",
       "\n",
       "       [[0, 1],\n",
       "        [0, 0],\n",
       "        [0, 0]]])"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "μP_x_y_w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[0, 0],\n",
       "        [0, 0],\n",
       "        [0, 0]],\n",
       "\n",
       "       [[1, 1],\n",
       "        [1, 1],\n",
       "        [1, 1]],\n",
       "\n",
       "       [[0, 0],\n",
       "        [0, 0],\n",
       "        [0, 0]],\n",
       "\n",
       "       [[0, 0],\n",
       "        [0, 0],\n",
       "        [0, 0]]], dtype=int64)"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "μA_x_y_w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "μA_x_y_w[1,:,:]=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "masked_array(\n",
       "  data=[[[0.6535895854646095, 1.6535895854646094],\n",
       "         [0.11500694312440574, 1.1150069431244058],\n",
       "         [0.9502828643490245, 1.9502828643490244]],\n",
       "\n",
       "        [[--, --],\n",
       "         [--, --],\n",
       "         [--, --]],\n",
       "\n",
       "        [[0.040709624769089126, 1.040709624769089],\n",
       "         [0.3971944613457058, 1.397194461345706],\n",
       "         [0.23313219734837998, 1.23313219734838]],\n",
       "\n",
       "        [[0.8417407242530616, 1.8417407242530617],\n",
       "         [0.20708234438674689, 1.2070823443867469],\n",
       "         [0.7424695335824086, 1.7424695335824087]]],\n",
       "  mask=[[[False, False],\n",
       "         [False, False],\n",
       "         [False, False]],\n",
       "\n",
       "        [[ True,  True],\n",
       "         [ True,  True],\n",
       "         [ True,  True]],\n",
       "\n",
       "        [[False, False],\n",
       "         [False, False],\n",
       "         [False, False]],\n",
       "\n",
       "        [[False, False],\n",
       "         [False, False],\n",
       "         [False, False]]],\n",
       "  fill_value=-10.0)"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tst = np.ma.masked_array(self.α_x_y_w , μA_x_y_w  != 0,fill_value=-10)\n",
    "tst"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.95028286,       -inf, 1.39719446, 1.84174072])"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "themax = tst.max(axis=(1,2),fill_value = -10)\n",
    "#themax.set_fill_value(-np.inf)\n",
    "np.ma.filled(themax, fill_value=-np.inf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.95028286e+00, 1.00000000e+20, 1.39719446e+00, 1.84174072e+00])"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "themax.data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "u_x = np.ma.masked_array(self.α_x_y_w , μA_x_y_w  == 0).max(axis = (1,2),fill_value=-10).data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.e+20, 1.e+20, 1.e+20, 1.e+20])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.ma.masked_array(self.α_x_y_w , μA_x_y_w  == 0).max(axis = (1,2)).data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-inf"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.ma.maximum_fill_value(u_x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.e+20, 1.e+20, 1.e+20, 1.e+20])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u_x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "μA_x_y_w = np.ones((self.nbx, self.nby,self.nbw), dtype = 'int64')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "μA_x_y_w[0,:,:]=0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "max() got an unexpected keyword argument 'max_filler'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-89-2b9085c28a15>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mma\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmasked_array\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mα_x_y_w\u001b[0m \u001b[1;33m,\u001b[0m \u001b[0mμA_x_y_w\u001b[0m  \u001b[1;33m==\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmax\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0maxis\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mmax_filler\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1000\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m: max() got an unexpected keyword argument 'max_filler'"
     ]
    }
   ],
   "source": [
    "np.ma.masked_array(self.α_x_y_w , μA_x_y_w  == 0).max(axis = (1,2),max_filler=-1000).data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1e+20"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.ma.masked_array(self.α_x_y_w , μA_x_y_w  == 0).get_fill_value()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[0, 0],\n",
       "        [0, 0],\n",
       "        [0, 0]],\n",
       "\n",
       "       [[1, 1],\n",
       "        [1, 1],\n",
       "        [1, 1]],\n",
       "\n",
       "       [[1, 1],\n",
       "        [1, 1],\n",
       "        [1, 1]],\n",
       "\n",
       "       [[1, 1],\n",
       "        [1, 1],\n",
       "        [1, 1]]], dtype=int64)"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "μA_x_y_w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[0.65358959, 1.65358959],\n",
       "        [0.11500694, 1.11500694],\n",
       "        [0.95028286, 1.95028286]],\n",
       "\n",
       "       [[0.4821914 , 1.4821914 ],\n",
       "        [0.87247454, 1.87247454],\n",
       "        [0.21233268, 1.21233268]],\n",
       "\n",
       "       [[0.04070962, 1.04070962],\n",
       "        [0.39719446, 1.39719446],\n",
       "        [0.2331322 , 1.2331322 ]],\n",
       "\n",
       "       [[0.84174072, 1.84174072],\n",
       "        [0.20708234, 1.20708234],\n",
       "        [0.74246953, 1.74246953]]])"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "self.α_x_y_w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Kelso-Crawford algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solveKelsoCrawford(self):\n",
    "    self.comp_nbsteps = 0\n",
    "    #tracemax = self.nbx*self.nby\n",
    "    #if trace:\n",
    "    #    self.traceu0_x_t = np.zeros((tracemax,self.nbx))    \n",
    "    μA_x_y_w = np.ones((self.nbx, self.nby,self.nbw), dtype = 'int64') # initially all offers are non rejected\n",
    "    while True :\n",
    "        u_x = np.maximum(np.ma.filled(np.ma.masked_array(self.α_x_y_w , μA_x_y_w  == 0).max(axis = (1,2)).filled_value=-np.inf),np.zeros(self.nbx))\n",
    "        μP_x_y_w = self.μ_from_u(u_x)\n",
    "        # μP_x_y_w = np.zeros((self.nbx, self.nby,self.nbw), dtype = 'int64')\n",
    "        # props_x = np.ma.masked_array(self.α_x_y_w, μA_x_y_w[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",
    "        \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",
    "\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.traceu0_x_t[self.comp_nbsteps,x] = np.sum(self.αo_x_y [x,:] * μP_x_y[x,:])\n",
    "        self.comp_nbsteps +=1\n",
    "        self.eq_μ_x_y = μE_x_y[0:-1,0:-1]\n",
    "        if np.max(np.abs(rej_x_y)) == 0: \n",
    "            if trace:\n",
    "                self.traceu0_x_t = self.traceu0_x_t[0:min(self.comp_nbsteps,tracemax),:]\n",
    "            break # if all offers have been accepted, then algorithm stops\n",
    "    return (0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function detects stable matchings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_stable(self, μ_x_y_w = None, output=0 ):\n",
    "    if (μ_x_y_w is None):\n",
    "        μ_x_y_w = self.eq_μ_x_y_w\n",
    "    μext_x_0 = 1-np.sum(μ_x_y_w,axis = (1,2))\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",
    "    # end for x and y\n",
    "    if output > 0 :\n",
    "        print('The matching is stable.')\n",
    "    return (True)\n",
    "\n",
    "NTU_market.is_stable = is_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_stable(μ_x_y,output=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tst = np.array([[0,1,2],[3,4,5]])\n",
    "tst[1,range(2)]=-1\n",
    "tst"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tst = np.array([[0,1,2],[3,4,5]])\n",
    "tst[1,0:2]=-1\n",
    "tst"
   ]
  },
  {
   "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_stable(μ_x_y,output=1)"
   ]
  },
  {
   "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.traceu0_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",
    "        \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",
    "\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.traceu0_x_t[self.comp_nbsteps,x] = np.sum(self.αo_x_y [x,:] * μP_x_y[x,:])\n",
    "        self.comp_nbsteps +=1\n",
    "        self.eq_μ_x_y = μE_x_y[0:-1,0:-1]\n",
    "        if np.max(np.abs(rej_x_y)) == 0: \n",
    "            if trace:\n",
    "                self.traceu0_x_t = self.traceu0_x_t[0:min(self.comp_nbsteps,tracemax),:]\n",
    "            break # if all offers have been accepted, then algorithm stops\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_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_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].\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 `u_from_v` and `v_from_u` 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 u_from_v(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",
    "\n",
    "NTU_market.u_from_v = u_from_v\n",
    "\n",
    "def  v_from_u(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.v_from_u = v_from_u\n",
    "\n",
    "def μ_from_u(self, uo_x):\n",
    "    return np.where( (self.αo_x_y == uo_x.reshape((-1,1))) , 1 , 0 )[:,0:-1]\n",
    "\n",
    "NTU_market.μ_from_u = μ_from_u\n",
    "\n",
    "def μ_from_v(self, vo_y):\n",
    "    return np.where( self.γo_x_y == vo_y , 1 , 0 )[0:-1,:]\n",
    "\n",
    "NTU_market.μ_from_v = μ_from_v\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.traceu0_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.u_from_v(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.traceu0_x_t = self.traceu0_x_t[0:min(self.comp_nbsteps,tracemax),:]\n",
    "            break\n",
    "        uo_x = uonew_x      \n",
    "        vo_y = self.v_from_u(uo_x) # each y proposes to favorite x among those willing to consider them\n",
    "        if output >= 2:\n",
    "            μP_x_y = self.μ_from_u(uo_x)\n",
    "            μE_x_y = self.μ_from_v(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.traceu0_x_t[self.comp_nbsteps,:] = uo_x\n",
    "        self.comp_nbsteps += 1\n",
    "    self.eq_μ_x_y = self.μ_from_v(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_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_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_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:\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  v_from_u_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.v_from_u_bis = v_from_u_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.traceu0_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.u_from_v(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.traceu0_x_t = self.traceu0_x_t[0:min(self.comp_nbsteps,tracemax),:]\n",
    "            break\n",
    "        uo_x = uonew_x      \n",
    "        vo_y = self.v_from_u_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.traceu0_x_t[self.comp_nbsteps,:] = uo_x\n",
    "        self.comp_nbsteps += 1\n",
    "    self.eq_μ_x_y = self.μ_from_v(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.traceu0_x_t\n",
    "    mkt_tst.solveAdachiBis(trace=True)\n",
    "    t2 = mkt_tst.traceu0_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",
    "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.u_from_v(vo_y)\n",
    "    vonew_y = self.v_from_u(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",
    "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.u_from_v(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.v_from_u(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.traceu0_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.traceu0_x_t = self.traceu0_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.traceu0_x_t[self.comp_nbsteps,:] = p_z[0:self.nbx]\n",
    "        self.comp_nbsteps += 1\n",
    "    self.eq_μ_x_y = self.μ_from_v(- 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.traceu0_x_t\n",
    "\n",
    "mkt_tst.solveCU(trace=True)\n",
    "t2 = mkt_tst.traceu0_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",
    "Gale and Shapley interpets 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.u_from_v(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.traceu0_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.traceu0_x_t = self.traceu0_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.traceu0_x_t[self.comp_nbsteps,:] = p_z[0:self.nbx]\n",
    "        self.comp_nbsteps += 1\n",
    "    self.eq_μ_x_y = self.μ_from_v(- 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.traceu0_x_t\n",
    "#print(mkt_tst.comp)\n",
    "\n",
    "mkt_tst.solveDampedCU(trace=True)\n",
    "t2 = mkt_tst.traceu0_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",
    "Cf ROTH-SOTOMAYOR P. 62"
   ]
  },
  {
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def next_us(self,uo_x,down = True):\n",
    "    μ_x_y = self.μ_from_u(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.transpose() == uo_x + incr).transpose() , 1 , 0 )[:,0:-1]\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": "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": "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",
    "    # H = nx.DiGraph()\n",
    "    us_list = [uo_x.tolist()]\n",
    "    i = 0\n",
    "    # H.add_nodes_from([i])\n",
    "    # H.nodes[i]['pos'] = (0,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",
    "        # H.add_nodes_from(inds_toadd)\n",
    "        # for node in inds_toadd:\n",
    "            # H.nodes[node]['pos'] = (0,H.nodes[i]['pos'][1]+1)\n",
    "        #H.add_edges_from([(i,j) for j in inds_alreadyhere+inds_toadd])\n",
    "        us_list = us_list + us_toadd \n",
    "        i += 1\n",
    "    return([np.array(u,dtype = 'int64') for u in us_list])\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",
    "the_u_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for u_x in the_u_s:\n",
    "    μ_x_y = rs_ex_2_17.μ_from_u(u_x)\n",
    "    print(rs_ex_2_17.is_stable(μ_x_y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Aggregate stable matchings\n",
    "\n",
    "The reference for this section is [GH19].\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def proposeNoHet(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",
    "        prefs_x_y = self.prefslistα_x_y\n",
    "    else:\n",
    "        n_x = self.m_y\n",
    "        nbx = self.nby\n",
    "        nby = self.nbx\n",
    "        prefs_x_y = self.prefslistγ_x_y.transpose()    \n",
    "    μ_x_y = np.zeros((nbx,nby))\n",
    "    for x in range(nbx):\n",
    "        nxres = n_x[x]\n",
    "        for yind in prefs_x_y[x,]:\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",
    "    return(μ_x_y)\n",
    "\n",
    "NTU_market.proposeNoHet = proposeNoHet\n",
    "\n",
    "def propose(self,axis,μbar_x_y,heterogeneity = 'none'):\n",
    "    if heterogeneity == 'none':\n",
    "        μ_x_y = self.proposeNoHet(axis,μbar_x_y)\n",
    "    else:\n",
    "        raise Exception(\"Heterogeneity \" + heterogeneity + \" is not supported.\")\n",
    "    return (μ_x_y)\n",
    "\n",
    "NTU_market.propose = propose\n",
    "\n",
    "def solveDARUM(self,output=0,trace=False,tol = 1e-5):\n",
    "    self.comp_nbsteps = 0\n",
    "    tracemax = self.nbx*self.nby\n",
    "    if trace:\n",
    "        self.traceu0_x_t = np.zeros((self.nbx,tracemax))    \n",
    "    μA_x_y = np.zeros((self.nbx, self.nby)) # initially all offers are non rejected\n",
    "    for x in range(self.nbx):\n",
    "        for y in range(self.nby):\n",
    "            μA_x_y[x,y] = max(self.n_x[x],self.m_y[y])\n",
    "    while True :\n",
    "        μP_x_y = self.propose(0 , μA_x_y,'none') # the x's pick their preferred offers among those not rejected\n",
    "        μE_x_y = self.propose(1, μP_x_y.transpose(),'none').transpose() # the y's pick their preferred offers among those made\n",
    "        rej_x_y = μP_x_y - μE_x_y # compute rejected offers\n",
    "        if np.max(np.abs(rej_x_y)) < tol: \n",
    "            if trace:\n",
    "                self.traceu0_x_t = self.traceu0_x_t[:,0:min(self.comp_nbsteps,tracemax)]\n",
    "            break # if all offers have been accepted (within tolerange), then algorithm stops\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(\"mu_P=, mu_E=\" )\n",
    "            print(μP_x_y)\n",
    "            print(μE_x_y)        \n",
    "        if trace and self.comp_nbsteps<tracemax:\n",
    "            for x in range(self.nbx):\n",
    "                self.traceu0_x_t[x,self.comp_nbsteps] = np.sum(self.αo_x_y [x,:] * μP_x_y[x,:])\n",
    "\n",
    "        self.comp_nbsteps +=1\n",
    "    self.eq_μ_x_y = μE_x_y\n",
    "    return (0)\n",
    "\n",
    "NTU_market.solveDARUM = solveDARUM "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}