hse2021task6.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Platform market equilibrium with Newton\n",
    "\n",
    "In this exercise you will compare the performance of successive\n",
    "approximation solver with the Newton solver and devise a poly-algorithm.\n",
    "\n",
    "Recall the model of platform market equilibrium from episode/video 22.\n",
    "\n",
    "- $ m $ products, $ i=1,\\dots,m $  \n",
    "- unit mass of consumers of $ n $ types with different preferences over the product  \n",
    "- $ p_j $ are fractions of consumer types in the population, $ j=1,\\dots,n $  \n",
    "- utility function of consumers of type $ j \\in \\{1,\\dots,n\\} $ from product $ i \\in \\{1,\\dots,m\\} $ is given by  \n",
    "\n",
    "\n",
    "$$\n",
    "u_{ij} = c_{ij} + s_i,\n",
    "$$\n",
    "\n",
    "- $ c_{ij} $ are *valuations* of each product by each of the consumer type  \n",
    "- $ s_i $ are market shares of each product, increase utility when more people are using the same platform  \n",
    "\n",
    "\n",
    "Random utility framework\n",
    "\n",
    "- value of choice $ U_{ij} = u_{ij} + \\epsilon[i] $ is effected by the random factors  \n",
    "- one $ \\epsilon[i] $ for each product choice, forming the vector $ \\epsilon = (\\epsilon[1],\\dots,\\epsilon[m]) $  \n",
    "- elements of $ \\epsilon $ are extreme value type I (EV1) random components of the utility, iid between consumers and products $ \\Rightarrow $ leading to the logit choice probabilities  \n",
    "\n",
    "\n",
    "$$\n",
    "P_{ij} = \\frac{\\exp(u_{ij})}{\\sum_{k=1}^m \\exp(u_{kj})}\n",
    "= \\frac{\\exp(u_{ij}-\\alpha)}{\\sum_{k=1}^m \\exp(u_{kj}-\\alpha)}, \\; \\forall \\alpha\n",
    "$$\n",
    "\n",
    "- assume choice is made once, and no changes are made afterwards  \n",
    "- under the assumption of unit mass of consumers we have  \n",
    "\n",
    "\n",
    "$$\n",
    "s_i = \\sum_{j=1}^n p_j P_{ij} = P \\cdot p\n",
    "$$\n",
    "\n",
    "- combining the last three expressions  \n",
    "\n",
    "\n",
    "$$\n",
    "u_{ij} = c_{ij} + \\sum_{t=1}^n p_t \\frac{\\exp(u_{it})}{\\sum_{k=1}^m \\exp(u_{kt})}\n",
    "$$\n",
    "\n",
    "- $ mn $ fixed point equations  \n",
    "- fixed point in the ($ mn $)-dimensional space of $ u_{ij} $  \n",
    "\n",
    "\n",
    "Using the code from the lecture notebook in episode 22, write a class for the platform market\n",
    "equilibrium model that would contain the following solvers for the equilibrium:\n",
    "- successive approximations (SA) solver (`solve_sa()` made into the method of the class)\n",
    "- multivariate Newton solver (see example in episode 23) that would solve the same fixed point equation\n",
    "- poly-algorithmic solver that runs a pre-specified number of SA iterations and then switches to the Newton method\n",
    "\n",
    "Using the appropriate callback function, illustrate the runs of the three solvers (initialized at the same starting point) to show their convergence speed.\n",
    "\n",
    "Find the number of SA iterations in the poly-algorithm that makes it the best method of the three (for a given draw of the utilities and market shares)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hints:\n",
    "\n",
    "Before coding up the Newton solver, let’s make it clear what dimensions different objects have.\n",
    "\n",
    "- valuations and utilities are in $ m \\times n $ matrix that is occasionally reshaped into a vector (one-dimensional by default in Python, we don’t need to specify whether it is column or row vector unless absolutely necessary). With the Pythonian row-major reshaping, `np.reshape(m*n)` and `np.reshape((m,n))` convert the utilities between  \n",
    "\n",
    "\n",
    "$$\n",
    "\\text{Row-major: }\n",
    "u =\n",
    "\\begin{bmatrix}\n",
    "u_{11} \\dots u_{1n}\\\\\n",
    "\\vdots \\ddots \\vdots\\\\\n",
    "u_{m1} \\dots u_{mn}\\\\\n",
    "\\end{bmatrix}\n",
    "\\leftrightarrow\n",
    "\\big(\n",
    "\\underbrace{u_{11},\\dots,u_{1n}}_{\\text{first row}},\n",
    "\\underbrace{u_{21},\\dots,u_{2n}}_{\\text{second row}},\\dots,\n",
    "\\underbrace{u_{m1},\\dots,u_{mn}}_{m\\text{-th row}}\n",
    "\\big)\n",
    "$$\n",
    "\n",
    "It is however more convenient to convert between matrix and vector form using column-major scheme, so that the choice probabilities of each consumer type form a continuous sub-vector. In Python this is possible using `np.reshape(m*n,order='F')` and `np.reshape((m,n),order='F')`, in which case the conversion is\n",
    "\n",
    "$$\n",
    "\\text{Column-major: }\n",
    "u =\n",
    "\\begin{bmatrix}\n",
    "u_{11} \\dots u_{1n}\\\\\n",
    "\\vdots \\ddots \\vdots\\\\\n",
    "u_{m1} \\dots u_{mn}\\\\\n",
    "\\end{bmatrix}\n",
    "\\leftrightarrow\n",
    "\\big(\n",
    "\\underbrace{u_{11},\\dots,u_{m1}}_{\\text{first column}},\n",
    "\\underbrace{u_{12},\\dots,u_{m2}}_{\\text{second column}},\\dots,\n",
    "\\underbrace{u_{1n},\\dots,u_{mn}}_{n\\text{-th column}}\n",
    "\\big)\n",
    "$$\n",
    "\n",
    "With the `order='F'` option the reshape may have to copy the data in memory which is inefficient.  To avoid this, the whole problem could have been set up transposed with $ m $ consumer types and $ n $ products.\n",
    "Try the two conversion styles in a separate cell with some arbitrary numbers.\n",
    "\n",
    "- conditional choice probabilities computed by `model.ccp(u)` is a $ m \\times n $ matrix with $ n $ columns containing choice probabilities by different consumer types  \n",
    "\n",
    "\n",
    "$$\n",
    "{ccp}_{ij} = [\\text{prob of choosing product } i \\text{ by consumer type j}] =\n",
    "\\frac{\\exp(u_{ij})}{\\sum_{k=1}^m \\exp(u_{kj})}\n",
    "$$\n",
    "\n",
    "- market shares computed by `model.shares(ccp)` is a column-vector of $ m $ elements computed as linear combination of the columns of the ccp matrix using consumer type weights  \n",
    "- fixed point operator `model.F(u)` takes $ mn $ utilities and returns the updated $ mn $ vector of utilities: output is one-dimensional array which can be used as input to `model.F(u)` as well  \n",
    "\n",
    "\n",
    "The Newton solver requires an different operator to solve for zeros. Let us define $ G(u) $ to take an array of $ mn $ utilities, and return a $ mn $ array which $ (i,j) $ element is computed as\n",
    "\n",
    "$$\n",
    "g_{ij} = u_{ij} - c_{ij} - \\sum_{t=1}^n  {ccp}_{it} p_t,\n",
    "i\\in\\{1,\\dots,m\\}, j\\in\\{1,\\dots,n\\},\n",
    "$$\n",
    "\n",
    "and which has the same structure as the **column-major vector form** of the utilities above. (Note that the successive approximation solver both types of conversions work equivalently.)\n",
    "\n",
    "The Jacobian of a multivalued function $ G(u) $ of multiple variables is a matrix of first derivatives of all the values of the function (rows) with respect to all of its variables (columns).\n",
    "Jacobian of $ G(u) $ is given by a $ mn \\times mn $ block matrix\n",
    "\n",
    "$$\n",
    "dG(u) =\n",
    "\\begin{bmatrix}\n",
    "A_{11} \\dots A_{1n}\\\\\n",
    "\\vdots \\ddots \\vdots\\\\\n",
    "A_{n1} \\dots A_{nn}\n",
    "\\end{bmatrix}, \\text{ where }\n",
    "A_{jl} =\n",
    "\\begin{bmatrix}\n",
    "\\frac{\\partial g_{1j}}{\\partial u_{1l}} \\dots \\frac{\\partial g_{mj}}{\\partial u_{1l}}\\\\\n",
    "\\vdots \\ddots \\vdots\\\\\n",
    "\\frac{\\partial g_{1j}}{\\partial u_{ml}} \\dots \\frac{\\partial g_{mj}}{\\partial u_{ml}}\n",
    "\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "Note that blocks $ A_{jl} $ correspond to the consumer types and the cross-effects between their utilities, and thus collects derivatives of the choice probabilities of the same consumer type.  This is due to the column-major vector conversion, and will be convenient for computations.\n",
    "\n",
    "To compute $ \\frac{\\partial g_{ij}}{\\partial u_{kl}} $, first note\n",
    "that the derivatives of choice probabilities with respect to utilities take convenient form\n",
    "\n",
    "$$\n",
    "\\frac{\\partial{ccp}_{ij}}{\\partial u_{kl}} =\n",
    "\\begin{cases}\n",
    "{ccp}_{ij}(1 - {ccp}_{ij}),& \\text{ when } i=k \\text{ and } j=l, \\\\\n",
    "- {ccp}_{ij} {ccp}_{kj},& \\text{ when } i \\ne k \\text{ but } j=l, \\\\\n",
    "0 ,& \\text{ when } j \\ne l.\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Then differentiating the expression for $ g_{ij} $ above, after some algebra we have\n",
    "\n",
    "$$\n",
    "\\frac{\\partial g_{ij}}{\\partial u_{kl}} =\n",
    "\\begin{cases}\n",
    "1 - {ccp}_{il}(1 - {ccp}_{il}) p_l,& \\text{ when } i=k \\text{ and } j=l, \\\\\n",
    "- {ccp}_{il}(1 - {ccp}_{il}) p_l,& \\text{ when } i=k \\text{ but } j \\ne l, \\\\\n",
    "{ccp}_{il} {ccp}_{kl} p_l,& \\text{ when } i \\ne k.\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "This formula can be applied directly in a nested loop to fill all values of the Jacobian $ dG(u) $. Such algorithm is clear but not efficient, as using Numpy’s matrix operations produces much faster code. Yet, it may be worth to start with a clear implementation, and then after testing and making sure that the code works, a better implementation can be written with direct comparison to the already existing one.\n",
    "\n",
    "Note that in the formula for $ \\frac{\\partial g_{ij}}{\\partial u_{kl}} $ all three cases have the generic component $ {ccp}_{il} {ccp}_{kl} p_l $ for all $ (i,j,k,l) $, with an addition of term $ -{ccp}_{il} p_l $ for the cases when $ i=k $, and 1 for the cases when in addition $ j=l $.\n",
    "This implies that the Jacobian matrix $ dG(u) $ can be constructed in three steps:\n",
    "- first, fill the $ mn \\times mn $ matrix with the generic component $ {ccp}_{il} {ccp}_{kl} p_l $,\n",
    "- second, subtract the term $ {ccp}_{il} p_l $ to the main diagonal of all the blocks,\n",
    "- third, add 1 to the main diagonal of the whole matrix.\n",
    "\n",
    "In reality first step is probably going to be done block-by-block in a loop for each $ l $, and the second step of adding diagonal elements (`numpy.diag()`) can then be performed simultaneously. Indexes $ i $ and $ k $ are then the indexes within the block, and the common component $ {ccp}_{il} {ccp}_{kl} p_l $ is given by an outer product of the vectors of choice probabilities (`numpy.outer()`).\n",
    "Observe as well that the right-hand-side of the formula does not depend on $ j $, implying that the blocks can be copied (`numpy.tile()`) along the corresponding dimension.\n",
    "However, other implementations avoiding loops entirely and instead relying on more efficient Numpy functions are also possible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0  1  2]\n",
      " [ 3  4  5]\n",
      " [ 6  7  8]\n",
      " [ 9 10 11]]\n",
      "[ 0  3  6  9  1  4  7 10  2  5  8 11]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "# example of column-major reshaping\n",
    "a = np.arange(12).reshape(4,3)\n",
    "print(a)\n",
    "b = a.reshape(12,order='F')\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of platform products = 2\n",
      "Number of consumer types = 3\n",
      "Population composition = [0.09865415 0.53618974 0.36515611]\n",
      "Valuations:\n",
      "[[0.09831954 0.67436529 0.65062023]\n",
      " [0.41568115 0.75820323 0.37061744]]\n",
      "x0= [0.09831954 0.41568115 0.67436529 0.75820323 0.65062023 0.37061744]\n",
      "dG loops :\n",
      "[[ 0.9759472   0.0240528  -0.13381216  0.13381216 -0.08952285  0.08952285]\n",
      " [ 0.0240528   0.9759472   0.13381216 -0.13381216  0.08952285 -0.08952285]\n",
      " [-0.0240528   0.0240528   0.86618784  0.13381216 -0.08952285  0.08952285]\n",
      " [ 0.0240528  -0.0240528   0.13381216  0.86618784  0.08952285 -0.08952285]\n",
      " [-0.0240528   0.0240528  -0.13381216  0.13381216  0.91047715  0.08952285]\n",
      " [ 0.0240528  -0.0240528   0.13381216 -0.13381216  0.08952285  0.91047715]]\n",
      "dG matrix:\n",
      "[[ 0.9759472   0.0240528  -0.13381216  0.13381216 -0.08952285  0.08952285]\n",
      " [ 0.0240528   0.9759472   0.13381216 -0.13381216  0.08952285 -0.08952285]\n",
      " [-0.0240528   0.0240528   0.86618784  0.13381216 -0.08952285  0.08952285]\n",
      " [ 0.0240528  -0.0240528   0.13381216  0.86618784  0.08952285 -0.08952285]\n",
      " [-0.0240528   0.0240528  -0.13381216  0.13381216  0.91047715  0.08952285]\n",
      " [ 0.0240528  -0.0240528   0.13381216 -0.13381216  0.08952285  0.91047715]]\n",
      "Differences:\n",
      "[[False False False False False False]\n",
      " [False False False False False False]\n",
      " [False False False False False False]\n",
      " [False False False False False False]\n",
      " [False False False False False False]\n",
      " [False False False False False False]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def solve_sa(F,x0,tol=1e-6,maxiter=100,callback=None,raise_error=True):\n",
    "    '''Computes the solution of fixed point equation x = F(x)\n",
    "    with given initial value x0 and algorithm parameters\n",
    "    Method: successive approximations\n",
    "    '''\n",
    "    for i in range(maxiter):  # main loop\n",
    "        x1 = F(x0)  # update approximation\n",
    "        err = np.amax(np.abs(x0-x1))  # allow for x to be array\n",
    "        if callback != None: callback(iter=i,err=err,x=x1,x0=x0)\n",
    "        if err<tol:\n",
    "            break  # break out if converged\n",
    "        x0 = x1  # get ready to the next iteration\n",
    "    else:\n",
    "        if raise_error:\n",
    "            raise RuntimeError('Failed to converge in %d iterations'%maxiter)\n",
    "    return x1\n",
    "\n",
    "def mnewton(fun,grad,x0,tol=1e-6,maxiter=100,callback=None):\n",
    "    '''Multinomial Newton method for solving systems of equation F(x)=0,\n",
    "    x is vector of m elements, F(x) is a m-valued function.\n",
    "    Given tolerance and number of iterations are used.\n",
    "    Callback function is invoked at each iteration if given.\n",
    "    '''\n",
    "    # conversion to array function of array argument\n",
    "    npfun = lambda x: np.asarray(fun(x))\n",
    "    npgrad = lambda x: np.asarray(grad(x))\n",
    "    for i in range(maxiter):\n",
    "        x1 = x0 - np.linalg.inv(npgrad(x0)) @ npfun(x0)  # matrix version\n",
    "        err = np.amax(np.abs(x1-x0)) # max over vector elements\n",
    "        if callback != None: callback(iter=i,err=err,x0=x0,x1=x1,fun=fun)\n",
    "        if err<tol: break\n",
    "        x0 = x1\n",
    "    else:\n",
    "        raise RuntimeError('Failed to converge in %d iterations'%maxiter)\n",
    "    return x1\n",
    "\n",
    "class model:\n",
    "    '''Simple platform equilibrium model'''\n",
    "\n",
    "    def __init__(self,m=2,n=2,seed=None):\n",
    "        '''Define default model parameters'''\n",
    "        self.m,self.n = m,n  # number of products and consumer types\n",
    "        np.random.seed(seed)\n",
    "        self.c = np.random.uniform(size=(m,n))    # valuations (random uniform)\n",
    "        self.p = np.random.dirichlet(np.ones(n))  # population composition (from symmetric Dirichlet distribution)\n",
    "\n",
    "    def __repr__(self):\n",
    "        return 'Number of platform products = {:d}\\nNumber of consumer types = {:d}\\nPopulation composition = {}\\nValuations:\\n{}'.format(self.m,self.n,self.p,self.c)\n",
    "\n",
    "    def ccp(self,u):\n",
    "        '''Conditional choice probabilities\n",
    "        Input: m*n array of utilities, to be reshaped\n",
    "        Output: m by n matrix\n",
    "        '''\n",
    "        u = np.asarray(u).reshape((self.m,self.n),order='F')   # convert to matrix, column-major reshaping\n",
    "        u = u - np.amax(u,axis=0) # de-max by column (avoid exp of large numbers)\n",
    "        e = np.exp(u)\n",
    "        esum = e.sum(axis=0)  # sums of exps\n",
    "        return e/esum         # matrix of choice probabilities\n",
    "\n",
    "    def shares(self,ccp):\n",
    "        '''Market shares from choice probabilities\n",
    "        Input: m by n matrix of ccps\n",
    "        Output: market shares, m by 1 column vector\n",
    "        '''\n",
    "        out = ccp @ self.p  # one-dim vector\n",
    "        return out[:,np.newaxis]  # column vector\n",
    "\n",
    "    def F(self,u):\n",
    "        '''Fixed point equation u=F(u)\n",
    "        Input: m*n array of utilities\n",
    "        Output: m*n array of utilities\n",
    "        '''\n",
    "        ccp = self.ccp(u)     # matrix of choice probabilities\n",
    "        sh = self.shares(ccp) # market shares\n",
    "        u1 = self.c + sh     # updated utilities\n",
    "        return u1.reshape(self.m*self.n,order='F') # return one dimensional array\n",
    "\n",
    "    def G(self,u):\n",
    "        '''LHS of the equation is standard form u-F(u)=0.\n",
    "        Input: m*n array of utilities\n",
    "        Output: m*n array of equation residuals\n",
    "        '''\n",
    "        return u - self.F(u)\n",
    "\n",
    "    def dGa(self,u):\n",
    "        '''Jacobian of G(u) computed in a direct and inefficient but more clear way\n",
    "        Input: m*n array of utilities\n",
    "        Output: m*n by m*n matrix of first derivatives\n",
    "        '''\n",
    "        ccp = self.ccp(u)\n",
    "        out = np.ones((self.m*self.n,self.m*self.n))\n",
    "        for i in range(self.m):\n",
    "            for j in range(self.n):\n",
    "                for k in range(self.m):\n",
    "                    for l in range(self.n):\n",
    "                        # formula for derivatives of dg(i,j)/du(k,l)\n",
    "                        if i==k and j==l:\n",
    "                            val = 1 - ccp[i,j]*(1-ccp[i,j])*self.p[l]\n",
    "                        elif i==k and j!=l:\n",
    "                            val = - ccp[i,l]*(1-ccp[i,l])*self.p[l]\n",
    "                        elif i!=k:\n",
    "                            val = ccp[i,l]*ccp[k,l]*self.p[l]\n",
    "                        else:\n",
    "                            raise error\n",
    "                        # column-major indexing\n",
    "                        out[j*self.m+i,l*self.m+k] = val\n",
    "        return out\n",
    "\n",
    "    def dGb(self,u):\n",
    "        '''Jacobian of G(u)\n",
    "        Input: m*n array of utilities\n",
    "        Output: m*n by m*n matrix of first derivatives\n",
    "        '''\n",
    "        ccp = self.ccp(u)\n",
    "        out = np.empty((self.m*self.n,self.m*self.n))\n",
    "        # step 1 and 2\n",
    "        for l in range(self.n):\n",
    "            block = self.p[l]*np.outer(ccp[:,l],ccp[:,l])\n",
    "            block -= np.diag(self.p[l]*ccp[:,l]) # add diagonal elements\n",
    "            out[:,l*self.m:(l+1)*self.m] = np.tile(block,(self.n,1))\n",
    "        # step 3\n",
    "        out += np.eye(self.m*self.n)\n",
    "        return out\n",
    "\n",
    "\n",
    "def printiter(**kwargs):\n",
    "    print('iter %d, err = %1.3e'%(kwargs['iter'],kwargs['err']))\n",
    "\n",
    "md = model(m=2,n=3,seed=291)\n",
    "# md = model(m=3,n=2)\n",
    "print(md)\n",
    "\n",
    "x0 = md.c.reshape(md.m*md.n,order='F')\n",
    "print('x0=',x0)\n",
    "\n",
    "print('dG loops :',md.dGa(x0),sep='\\n')\n",
    "print('dG matrix:',md.dGb(x0),sep='\\n')\n",
    "print('Differences:',md.dGa(x0)-md.dGb(x0) > 1e-8,sep='\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test of Jacobian PASSED (diff = 3.0322590144e-07 > 1.0000000000e-05)\n",
      "Test of Jacobian PASSED (diff = 3.0322590145e-07 > 1.0000000000e-05)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# derivative check\n",
    "from scipy.optimize import check_grad\n",
    "def test_jackobian(fun,jac,x,epsilon=1e-8,test_tol=1e-5):\n",
    "    '''Test the matrix of analytic derivatives against the numerical ones'''\n",
    "    n = np.size(x[0]) # number of variables\n",
    "    data = np.empty(x.shape) # test results\n",
    "    for k,ix in enumerate(x): # over all points (which are in rows)\n",
    "        for i in range(n): # over all function outpus\n",
    "            func = lambda x: fun(x)[i]\n",
    "            grad = lambda x: jac(x)[i,:]\n",
    "            data[k,i] = check_grad(func,grad,ix,epsilon=epsilon)\n",
    "    res = np.max(np.absolute(data))\n",
    "    if res<test_tol:\n",
    "        print(f'Test of Jacobian PASSED (diff = {res:1.10e} > {test_tol:1.10e})')\n",
    "        return True\n",
    "    else:\n",
    "        print(f'Test of Jacobian FAILED (diff = {res:1.10e} > {test_tol:1.10e})')\n",
    "        return False\n",
    "\n",
    "points = -10 + 20*np.random.uniform(size=(md.m*md.n,md.m*md.n))\n",
    "test_jackobian(md.G,md.dGa,points)\n",
    "test_jackobian(md.G,md.dGb,points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SA:\n",
      "iter 0, err = 5.064e-01\n",
      "iter 1, err = 3.167e-03\n",
      "iter 2, err = 1.567e-03\n",
      "iter 3, err = 7.748e-04\n",
      "iter 4, err = 3.832e-04\n",
      "iter 5, err = 1.895e-04\n",
      "iter 6, err = 9.373e-05\n",
      "iter 7, err = 4.635e-05\n",
      "iter 8, err = 2.292e-05\n",
      "iter 9, err = 1.134e-05\n",
      "iter 10, err = 5.606e-06\n",
      "iter 11, err = 2.773e-06\n",
      "iter 12, err = 1.371e-06\n",
      "iter 13, err = 6.781e-07\n",
      "iter 14, err = 3.354e-07\n",
      "iter 15, err = 1.658e-07\n",
      "iter 16, err = 8.202e-08\n",
      "iter 17, err = 4.056e-08\n",
      "iter 18, err = 2.006e-08\n",
      "iter 19, err = 9.920e-09\n",
      "iter 20, err = 4.906e-09\n",
      "iter 21, err = 2.426e-09\n",
      "iter 22, err = 1.200e-09\n",
      "iter 23, err = 5.934e-10\n",
      "iter 24, err = 2.935e-10\n",
      "iter 25, err = 1.451e-10\n",
      "iter 26, err = 7.177e-11\n",
      "SA: Equilibrium found!\n",
      "Equilibrium choice probabilities:\n",
      "[[0.42750812 0.48537877 0.57574719]\n",
      " [0.57249188 0.51462123 0.42425281]]\n",
      "Equilibrium market shares:\n",
      "[0.51266817 0.48733183]\n"
     ]
    }
   ],
   "source": [
    "print('SA:')\n",
    "x = solve_sa(md.F,x0=x0,tol=1e-10,callback=printiter)\n",
    "print('SA: Equilibrium found!')\n",
    "ccp = md.ccp(x)\n",
    "shares = md.shares(ccp).squeeze()  # make one-dim array\n",
    "print('Equilibrium choice probabilities:',ccp,'Equilibrium market shares:',shares,sep='\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Newton:\n",
      "iter 0, err = 5.127e-01\n",
      "iter 1, err = 2.587e-06\n",
      "iter 2, err = 1.616e-13\n",
      "Newton: Equilibrium found!\n",
      "Equilibrium choice probabilities:\n",
      "[[0.42750812 0.48537877 0.57574719]\n",
      " [0.57249188 0.51462123 0.42425281]]\n",
      "Equilibrium market shares:\n",
      "[0.51266817 0.48733183]\n"
     ]
    }
   ],
   "source": [
    "print('Newton:')\n",
    "x = mnewton(md.G,md.dGb,x0=x0,tol=1e-10,callback=printiter)\n",
    "print('Newton: Equilibrium found!')\n",
    "ccp = md.ccp(x)\n",
    "shares = md.shares(ccp).squeeze()  # make one-dim array\n",
    "print('Equilibrium choice probabilities:',ccp,'Equilibrium market shares:',shares,sep='\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "First doing 3 iterations of successive approximation solver\n",
      "iter 0, err = 5.064e-01\n",
      "iter 1, err = 3.167e-03\n",
      "iter 2, err = 1.567e-03\n",
      "Afterwards switching to Newton solver\n",
      "iter 0, err = 1.533e-03\n",
      "iter 1, err = 5.441e-08\n",
      "iter 2, err = 2.220e-16\n",
      "Newton: Equilibrium found!\n",
      "Equilibrium choice probabilities:\n",
      "[[0.42750812 0.48537877 0.57574719]\n",
      " [0.57249188 0.51462123 0.42425281]]\n",
      "Equilibrium market shares:\n",
      "[0.51266817 0.48733183]\n"
     ]
    }
   ],
   "source": [
    "def solve_poly(funSA,funN,gradN,x0,sa_iter=5,tol=1e-6,maxiter=100,callback=None):\n",
    "    '''Polyalgorithm to solve by first doing sa_iter and then switching to Newton method'''\n",
    "    print(f'First doing {sa_iter} iterations of successive approximation solver')\n",
    "    x1 = solve_sa(funSA,x0=x0,maxiter=sa_iter,tol=tol,callback=callback,raise_error=False)\n",
    "    # TODO: would be good to check for convergence here and not call Newton if already converged\n",
    "    print(f'Afterwards switching to Newton solver')\n",
    "    mnewton(funN,gradN,x0=x1,tol=tol,callback=callback,maxiter=maxiter-sa_iter)\n",
    "\n",
    "number_sa_iter = 3 # change to see how performance is affected\n",
    "solve_poly(md.F,md.G,md.dGb,x0=x0,sa_iter=number_sa_iter,tol=1e-10,callback=printiter)\n",
    "print('Newton: Equilibrium found!')\n",
    "ccp = md.ccp(x)\n",
    "shares = md.shares(ccp).squeeze()  # make one-dim array\n",
    "print('Equilibrium choice probabilities:',ccp,'Equilibrium market shares:',shares,sep='\\n')"
   ]
  }
 ],
 "metadata": {
  "date": 1638093387.437004,
  "filename": "hse2021task6.rst",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "title": "Platform market equilibrium with Newton"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}