et06_matching-as-lcp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Matching as a LCP</center>\n",
    "### <center>Alfred Galichon (NYU & Sciences Po)</center>\n",
    "## <center>'math+econ+code' masterclass</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Lemke-Howson algorithm\n",
    "\n",
    "### References\n",
    "\n",
    "XXX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!pip install gurobipy\n",
    "import numpy as np\n",
    "import scipy.sparse as spr\n",
    "import gurobipy as grb\n",
    "import sympy\n",
    "import math\n",
    "import lemkelcp as lcp\n",
    "from sympy.solvers import solve\n",
    "from sympy import *\n",
    "\n",
    "def round_expr(expr, num_digits):\n",
    "    return expr.xreplace({n : round(n, num_digits) for n in expr.atoms(Number)})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setting up and solving the LCP\n",
    "\n",
    "We shall solve the problem as a LCP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class LCP: # z >=0, w = M z + q >= 0, z.w=0\n",
    "    def __init__(self,M_k_k,q_k):\n",
    "        self.M_k_k = M_k_k\n",
    "        self.nbk,_ = M_k_k.shape\n",
    "        self.q_k = q_k\n",
    "    \n",
    "    def qp_solve(self,silent = True):\n",
    "        qp = grb.Model()\n",
    "        if silent:\n",
    "            qp.Params.OutputFlag = 0\n",
    "        qp.Params.NonConvex = 2\n",
    "        zqp_k = qp.addMVar(shape = self.nbk)\n",
    "        wqp_k = qp.addMVar(shape = self.nbk)\n",
    "        qp.addConstr(self.M_k_k @ zqp_k - wqp_k == - self.q_k)\n",
    "        qp.setObjective(zqp_k @ wqp_k, sense = grb.GRB.MINIMIZE)\n",
    "        qp.optimize()\n",
    "        z_k = np.array(qp.getAttr('x'))[:self.nbk]\n",
    "        return(z_k)         "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Set parameter Username\n",
      "Academic license - for non-commercial use only - expires 2023-12-23\n",
      "[0. 1. 0.]\n",
      "[4. 0. 0.]\n",
      "0.0\n"
     ]
    }
   ],
   "source": [
    "kostrava_ex = LCP(np.array([[ 0,  3,  1],[-3,  2, -2],[-1, -2,  2]]) , np.array([1,-2,2]))\n",
    "zsol = kostrava_ex.qp_solve()\n",
    "print(zsol)\n",
    "wsol = kostrava_ex.M_k_k @ zsol + kostrava_ex.q_k\n",
    "print(wsol)\n",
    "print(zsol @ wsol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "murty_ex_2_10 = LCP(np.array([[ 1,0, 0],[2, 1,0],[2,2,1]]) , np.array([-8,-12,-14]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving LCP using lemkelcp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#def NashEquilibrium_lcplemke_solve(self, silent = True):\n",
    "#    sol_z = self.lcp_pb.solve_lemkelcp()\n",
    "#    x_i = sol_z[:self.nbi]\n",
    "#    y_j = sol_z[self.nbi:]\n",
    "#    val1 = 1 / y_j.sum()\n",
    "#    val2 = 1 /  x_i.sum()\n",
    "#    p_i = x_i * val2\n",
    "#    q_j = y_j * val1\n",
    "#    sol_dict = {'val1':val1, 'val2':val2, 'p_i':p_i,'q_j':q_j}\n",
    "#    return(sol_dict)\n",
    "\n",
    "#NashEquilibrium.lcplemke_solve = NashEquilibrium_lcplemke_solve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lemke algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tableaux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Tableau():\n",
    "    def __init__(self, names_basic,names_nonbasic, M_i_j,c_i,obj_j = None): #z_B = c - M @ z_N\n",
    "        self.names_initial_basic = names_basic\n",
    "        self.names_initial_nonbasic = names_nonbasic\n",
    "        self.nonbasic = list(symbols(names_nonbasic))\n",
    "        if obj_j is None:\n",
    "            self.base = dict()\n",
    "        else:\n",
    "            self.base = { Symbol('obj') : obj_j @ self.nonbasic }\n",
    "        basic = list(symbols(names_basic))\n",
    "        self.base.update( { basic[i]: c_i[i]  - (M_i_j @ self.nonbasic)[i]  for i in range(len(c_i))} ) \n",
    "        \n",
    "        \n",
    "    def variables(self):\n",
    "        return( list(self.base.keys())[1:] +self.nonbasic) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Tableau_determine_departing(self,entering_var):\n",
    "    runmin = float('inf')\n",
    "    departing_var = None \n",
    "    for var in self.base.keys() - {Symbol('obj')}:\n",
    "        the_expr_list = solve(self.base[var] - var,entering_var)\n",
    "        if the_expr_list: # if one can invert the previous expression\n",
    "            the_expr = the_expr_list[0] # express entering variable as a function of the other ones:\n",
    "            val_entering_var = the_expr.subs([ (variable,0) for variable in [var]+self.nonbasic])\n",
    "            if (val_entering_var > 0) & (val_entering_var < runmin) :\n",
    "                runmin,departing_var = val_entering_var, var\n",
    "    return departing_var # if no variable is found, None retunrned\n",
    "\n",
    "Tableau.determine_departing = Tableau_determine_departing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Tableau_determine_departing_complementary(self,entering_var):\n",
    "    runmax = - float('inf')\n",
    "    departing_var = None \n",
    "    for var in self.base.keys() - {Symbol('obj')}:\n",
    "        the_expr_list = solve(self.base[var] - var,entering_var)\n",
    "        if the_expr_list: # if one can invert the previous expression\n",
    "            the_expr = the_expr_list[0] # express entering variable as a function of the other ones:\n",
    "            val_entering_var = the_expr.subs([ (variable,0) for variable in [var]+self.nonbasic])\n",
    "            if (val_entering_var > 0) & (val_entering_var > runmax) :\n",
    "                runmax,departing_var = val_entering_var, var\n",
    "    return departing_var # if no variable is found, None retunrned\n",
    "\n",
    "Tableau.determine_departing_complementary = Tableau_determine_departing_complementary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Tableau_pivot(self,entering_var,departing_var, verbose = 0):\n",
    "    expr_entering = solve(self.base[departing_var] - departing_var,entering_var)[0]\n",
    "    for var in self.base:\n",
    "        self.base[var] = self.base[var].subs([(entering_var, expr_entering)])  \n",
    "    self.base[entering_var] = expr_entering\n",
    "    del self.base[departing_var]\n",
    "    self.nonbasic.remove(entering_var)\n",
    "    self.nonbasic.append(departing_var)\n",
    "    if verbose >0:\n",
    "        print('Entering = ' + str( entering_var)+'; departing = '+ str( departing_var))\n",
    "    if verbose >1:\n",
    "        print(str( entering_var)+' = '+str(round_expr(expr_entering,2)))\n",
    "    return expr_entering\n",
    "\n",
    "Tableau.pivot = Tableau_pivot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Tableau_evaluate(self,thevar):\n",
    "    if thevar in set(self.nonbasic):\n",
    "        return 0.0\n",
    "    else:\n",
    "        return np.float(self.base[thevar].evalf(subs = {variable:0.0 for variable in self.nonbasic} ))\n",
    "\n",
    "Tableau.evaluate=Tableau_evaluate\n",
    "\n",
    "def Tableau_print_solution(self,title=None):\n",
    "    if not (title is None):\n",
    "        print(title)\n",
    "    for var in self.base:\n",
    "        print(str(var)+'='+str(self.base[var].subs([ (variable,0) for variable in self.nonbasic])))\n",
    "\n",
    "Tableau.print_solution = Tableau_print_solution\n",
    "\n",
    "def Tableau_output_solution(self):\n",
    "    xB = np.zeros(len(self.names_initial_basic))\n",
    "    xN = np.zeros(len(self.names_initial_nonbasic))\n",
    "    for i,name in enumerate(self.names_initial_basic):\n",
    "        xB[i]=self.evaluate(Symbol(name))\n",
    "    for i,name in enumerate(self.names_initial_nonbasic):\n",
    "        xN[i]=self.evaluate(Symbol(name))\n",
    "    return xB,xN\n",
    "    \n",
    "Tableau.output_solution = Tableau_output_solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def LCP_lemke_howson_solve(self,verbose=2):\n",
    "    zks = ['z_' + str(k+1) for k in range(self.nbk)]+['z_0']\n",
    "    wks = ['w_' + str(k+1) for k in range(self.nbk)]\n",
    "    tab = Tableau(wks, zks, - np.block([self.M_k_k,np.ones((self.nbk,1))]), self.q_k)\n",
    "    keys = zks[:-1]+wks\n",
    "    labels = wks+zks[:-1]\n",
    "    complements = {Symbol(keys[t]): Symbol(labels[t]) for t in range(len(keys))}\n",
    "    entering_var = Symbol('z_0')\n",
    "    departing_var=tab.determine_departing_complementary(entering_var)\n",
    "    while True:\n",
    "        tab.pivot(entering_var,departing_var,verbose=verbose)\n",
    "        if departing_var == Symbol('z_0'):\n",
    "            break\n",
    "        else:\n",
    "            entering_var = complements[departing_var]\n",
    "        departing_var = tab.determine_departing(entering_var)\n",
    "    _,z_k = tab.output_solution()\n",
    "    return z_k[:-1]\n",
    "\n",
    "LCP.lemke_howson_solve = LCP_lemke_howson_solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Entering = z_0; departing = w_3\n",
      "z_0 = w_3 - 2.0*z_1 - 2.0*z_2 - z_3 + 14.0\n",
      "Entering = z_3; departing = w_2\n",
      "z_3 = -w_2 + w_3 - z_2 + 2.0\n",
      "Entering = z_2; departing = z_3\n",
      "z_2 = -w_2 + w_3 - z_3 + 2.0\n",
      "Entering = w_3; departing = w_1\n",
      "w_3 = -w_1 + 2.0*w_2 - z_1 + z_3 + 2.0\n",
      "Entering = z_1; departing = w_3\n",
      "z_1 = -w_1 + 2.0*w_2 - w_3 + z_3 + 2.0\n",
      "Entering = z_3; departing = z_2\n",
      "z_3 = -w_2 + w_3 - z_2 + 2.0\n",
      "Entering = w_2; departing = z_3\n",
      "w_2 = w_3 - z_2 - z_3 + 2.0\n",
      "Entering = w_3; departing = z_0\n",
      "w_3 = 2.0*w_1 - z_0 + 2.0*z_2 + z_3 + 2.0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([8., 0., 0.])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "murty_ex_2_10.lemke_howson_solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "murty_ex_2_9 = LCP(np.array([[ 1,0, 3],[-1, 2,5],[2,1,2]]) , np.array([-3,-2,-1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Entering = z_0; departing = w_1\n",
      "z_0 = w_1 - z_1 - 3.0*z_3 + 3.0\n",
      "Entering = z_1; departing = w_2\n",
      "z_1 = 0.5*w_1 - 0.5*w_2 + z_2 + z_3 + 0.5\n",
      "Entering = z_2; departing = z_0\n",
      "z_2 = 0.5*w_1 + 0.5*w_2 - z_0 - 4.0*z_3 + 2.5\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([3. , 2.5, 0. ])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "murty_ex_2_9.lemke_howson_solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0. , 0. , 7.5])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(murty_ex_2_9.M_k_k @  np.array([3. , 2.5, 0. ]) +murty_ex_2_9.q_k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3.0000000e+00, 2.5000000e+00, 2.2014867e-10])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "murty_ex_2_9.qp_solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Matching model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Matching_market():\n",
    "    def __init__(self, n_x=None,m_y=None):\n",
    "        if n_x is None:\n",
    "            self.n_x = np.ones(self.nbx)\n",
    "        else:\n",
    "            self.n_x = n_x\n",
    "        if m_y is None:\n",
    "            self.m_y = np.ones(self.nby)\n",
    "        else:\n",
    "            self.m_y = m_y\n",
    "        self.nba = self.nbx * self.nby + self.nbx + self.nby  \n",
    "        self.MOM_z_xy = spr.bmat([[spr.kron(spr.identity(self.nbx), np.ones((1,self.nby)))],\n",
    "                                 [spr.kron(np.ones((1,self.nbx)),spr.identity(self.nby))]]) \n",
    "        \n",
    "        \n",
    "    def lcp_solve(self,silent=True,precision = 16):\n",
    "        qp = grb.Model()\n",
    "        if silent:\n",
    "            qp.Params.OutputFlag = 0\n",
    "        qp.Params.NonConvex = 2\n",
    "        self.qp_setup(qp)\n",
    "        qp.optimize()\n",
    "        sol = qp.getAttr('x')\n",
    "        ustar_x = np.array(sol[:self.nbx])\n",
    "        vstar_y = np.array(sol[self.nbx:(self.nbx+self.nby)])\n",
    "        mustar_x_y =np.array(sol[(self.nbx+self.nby):(self.nbx+self.nby+self.nbx*self.nby)]).reshape((self.nbx,self.nby)) \n",
    "        return {'scalprod':qp.getAttr('objval'),\n",
    "                'u_x': ustar_x.round(precision),\n",
    "                'v_y': vstar_y.round(precision),\n",
    "                'mu_x_y':mustar_x_y.round(precision)}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TU markets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "class TU_market(Matching_market):\n",
    "    def __init__(self,Φ_x_y,n_x=None,m_y = None):\n",
    "            self.nbx,self.nby = Φ_x_y.shape\n",
    "            self.Φ_x_y = Φ_x_y\n",
    "            Matching_market.__init__(self,n_x,m_y)\n",
    "        \n",
    "    def qp_setup(self,qp):\n",
    "        nbx,nby = self.nbx,self.nby\n",
    "\n",
    "        p_z = qp.addMVar(shape = nbx+nby)\n",
    "        mu_xy = qp.addMVar(shape = nbx*nby)\n",
    "\n",
    "        dualp_z = qp.addMVar(shape = nbx+nby)\n",
    "        dualmu_xy = qp.addMVar(shape = nbx*nby)\n",
    "\n",
    "        qp.addConstr(dualp_z == np.concatenate([self.n_x,self.m_y]) - self.MOM_z_xy @ mu_xy )\n",
    "        qp.addConstr(dualmu_xy == self.MOM_z_xy.T @ p_z -self.Φ_x_y.flatten() )\n",
    "\n",
    "        qp.setObjective(p_z @ dualp_z + mu_xy @ dualmu_xy   , sense = grb.GRB.MINIMIZE)\n",
    "        return(qp)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'scalprod': 2.343516718615645e-08,\n",
       " 'u_x': array([0.25053, 0.     , 0.09625]),\n",
       " 'v_y': array([0.80375, 0.64947]),\n",
       " 'mu_x_y': array([[0., 1.],\n",
       "        [0., 0.],\n",
       "        [1., 0.]])}"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "TU_ex=TU_market(np.array([[.2,.6],[.3,-.1],[.7,.5]])+np.array([[.7,.3],[.3,.6],[0.2,-.2]]))\n",
    "TU_ex.lcp_solve(precision=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### LTU markets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "class LTU_market(Matching_market):\n",
    "    def __init__(self,ζ_x_y,Φ_x_y,n_x=None,m_y = None):\n",
    "            self.nbx,self.nby = Φ_x_y.shape\n",
    "            self.Φ_x_y = Φ_x_y\n",
    "            self.ζ_x_y = ζ_x_y\n",
    "            Matching_market.__init__(self,n_x,m_y)\n",
    "\n",
    "    def qp_setup(self,qp):\n",
    "        nbx,nby = self.nbx,self.nby\n",
    "\n",
    "        p_z = qp.addMVar(shape = nbx+nby)\n",
    "        mu_xy = qp.addMVar(shape = nbx*nby)\n",
    "\n",
    "        dualp_z = qp.addMVar(shape = nbx+nby)\n",
    "        dualmu_xy = qp.addMVar(shape = nbx*nby)\n",
    "\n",
    "        zetaMXT = spr.diags(self.ζ_x_y.flatten()) @ spr.kron(spr.identity(self.nbx), np.ones((self.nby,1)))\n",
    "        oneminuszetaMYT = spr.diags(1 - self.ζ_x_y.flatten()) @ spr.kron(np.ones((self.nbx,1)),spr.identity(self.nby))\n",
    "        \n",
    "        qp.addConstr(dualp_z == np.concatenate([self.n_x,self.m_y]) - self.MOM_z_xy @ mu_xy )\n",
    "        qp.addConstr(dualmu_xy == spr.hstack([zetaMXT,oneminuszetaMYT]) @ p_z -self.Φ_x_y.flatten()/2 )\n",
    "\n",
    "        qp.setObjective(p_z @ dualp_z + mu_xy @ dualmu_xy , sense = grb.GRB.MINIMIZE)\n",
    "        return(qp)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'scalprod': 2.425545405086363e-08,\n",
       " 'u_x': array([0.10396, 0.     , 0.02131]),\n",
       " 'v_y': array([0.42869, 0.34604]),\n",
       " 'mu_x_y': array([[0., 1.],\n",
       "        [0., 0.],\n",
       "        [1., 0.]])}"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LTU_ex=LTU_market(np.array([[.5,.5],[.5,.5],[.5,.5]]),np.array([[0.9, 0.9],[0.6, 0.5],[0.9, 0.3]])/2)\n",
    "LTU_ex.lcp_solve(precision=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NTU markets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "class NTU_market(Matching_market):\n",
    "    def __init__(self,α_x_y,γ_x_y,n_x=None,m_y = None):\n",
    "        self.nbx,self.nby = α_x_y.shape\n",
    "        self.α_x_y = α_x_y\n",
    "        self.γ_x_y = γ_x_y\n",
    "        Matching_market.__init__(self,n_x,m_y)\n",
    "\n",
    "    \n",
    "    def qp_setup(self,qp):\n",
    "        nbk = 2\n",
    "        nbx,nby = self.nbx,self.nby\n",
    "        \n",
    "        p_z = qp.addMVar(shape = nbx+nby)\n",
    "        mu_xy = qp.addMVar(shape = nbx*nby)\n",
    "        pi_kxy = qp.addMVar(shape = (nbk*nbx*nby))\n",
    "        D_xy = qp.addMVar(shape = nbx*nby, lb = -grb.GRB.INFINITY)\n",
    "\n",
    "        dualp_z = qp.addMVar(shape = nbx+nby)\n",
    "        dualmu_xy = qp.addMVar(shape = nbx*nby)\n",
    "        dualpi_kxy = qp.addMVar(shape = (nbk*nbx*nby))\n",
    "\n",
    "        theDks = [spr.hstack([spr.kron(spr.identity(nbx), np.ones((nby,1))),spr.csr_matrix((nbx*nby,nby))]),\n",
    "                spr.hstack([spr.csr_matrix((nbx*nby,nbx)), spr.kron(np.ones((nbx,1)),spr.identity(nby))])]\n",
    "        the_cks = [self.α_x_y.flatten(),self.γ_x_y.flatten()]\n",
    "\n",
    "        qp.addConstr(dualp_z == np.concatenate([self.n_x,self.m_y]) - self.MOM_z_xy @ mu_xy )\n",
    "        qp.addConstr(dualmu_xy == D_xy)\n",
    "        for k in range(nbk):\n",
    "            qp.addConstr(dualpi_kxy[(k*nbx*nby):((k+1)*nbx*nby)] == D_xy - theDks[k] @ p_z + the_cks[k]  ) \n",
    "        qp.addConstr(spr.kron(np.ones((1,nbk)), spr.identity(nbx*nby)) @ pi_kxy == np.ones((1,nbx*nby)))\n",
    "\n",
    "        qp.setObjective(p_z @ dualp_z + mu_xy @ dualmu_xy + dualpi_kxy @ pi_kxy  , sense = grb.GRB.MINIMIZE)\n",
    "        return(qp)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'scalprod': 5.398852163033417e-08,\n",
       " 'u_x': array([0.26858, 0.05265, 0.     ]),\n",
       " 'v_y': array([0.3, 0.3]),\n",
       " 'mu_x_y': array([[0., 1.],\n",
       "        [1., 0.],\n",
       "        [0., 0.]])}"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "NTU_ex=NTU_market(np.array([[.2,.6],[.3,-.1],[.7,.5]]),\n",
    "                  np.array([[.7,.3],[.3,.6],[0.2,-.2]]))\n",
    "NTU_ex.lcp_solve(precision=5)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ETU markets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ETU_market(Matching_market):\n",
    "    def __init__(self,α_x_y,γ_x_y,τ,n_x=None,m_y = None):\n",
    "        self.nbx,self.nby = α_x_y.shape\n",
    "        self.α_x_y = α_x_y\n",
    "        self.γ_x_y = γ_x_y\n",
    "        self.τ = τ\n",
    "        Matching_market.__init__(self,n_x,m_y)\n",
    "    \n",
    "    def qp_setup(self,qp):\n",
    "        nbx,nby = self.nbx,self.nby\n",
    "        ζ_x_y = np.exp(- self.α_x_y / self.τ) / (np.exp(- self.α_x_y / self.τ) + np.exp(- self.γ_x_y / self.τ))\n",
    "        Φ_x_y = 4 / (np.exp(- self.α_x_y / self.τ) + np.exp(- self.γ_x_y / self.τ)) - 2\n",
    "        p_z = qp.addMVar(shape = nbx+nby)\n",
    "        mu_xy = qp.addMVar(shape = nbx*nby)\n",
    "\n",
    "        dualp_z = qp.addMVar(shape = nbx+nby)\n",
    "        dualmu_xy = qp.addMVar(shape = nbx*nby)\n",
    "\n",
    "        zetaMXT = spr.diags(ζ_x_y.flatten()) @ spr.kron(spr.identity(self.nbx), np.ones((self.nby,1)))\n",
    "        oneminuszetaMYT = spr.diags(1 - ζ_x_y.flatten()) @ spr.kron(np.ones((self.nbx,1)),spr.identity(self.nby))\n",
    "        \n",
    "        qp.addConstr(dualp_z == np.concatenate([self.n_x,self.m_y]) - self.MOM_z_xy @ mu_xy )\n",
    "        qp.addConstr(dualmu_xy == spr.hstack([zetaMXT,oneminuszetaMYT]) @ p_z - Φ_x_y.flatten()/2 )\n",
    "\n",
    "        qp.setObjective(p_z @ dualp_z + mu_xy @ dualmu_xy , sense = grb.GRB.MINIMIZE)\n",
    "        return(qp)\n",
    "\n",
    "        \n",
    "    def lcp_solve(self,silent=True,precision = 16):\n",
    "        interim_sol = Matching_market.lcp_solve(self,silent,16)\n",
    "        return {'scalprod':interim_sol['scalprod'],\n",
    "                'u_x': (self.τ*np.log(interim_sol['u_x'])).round(precision),\n",
    "                'v_y': (self.τ*np.log(interim_sol['v_y'])).round(precision),\n",
    "                'mu_x_y':(interim_sol['mu_x_y']).round(precision)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'scalprod': 1.9848958712543625e-08,\n",
       " 'u_x': array([  2.60582,   1.4281 , -20.22011]),\n",
       " 'v_y': array([3.52637, 3.64928]),\n",
       " 'mu_x_y': array([[0.49988, 0.50012],\n",
       "        [0.50012, 0.49988],\n",
       "        [0.     , 0.     ]])}"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ETU_ex=ETU_market(np.array([[2,6],[3,1],[7,5]]),\n",
    "                  np.array([[7,3],[3,6],[2,2]]),1)\n",
    "ETU_ex.lcp_solve(precision=5)\n"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}