{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear complementarity problem methods\n",
    "\n",
    "**Randall Romero Aguilar, PhD**\n",
    "\n",
    "This demo is based on the original Matlab demo accompanying the  <a href=\"https://mitpress.mit.edu/books/applied-computational-economics-and-finance\">Computational Economics and Finance</a> 2001 textbook by Mario Miranda and Paul Fackler.\n",
    "\n",
    "Original (Matlab) CompEcon file: **demslv07.m**\n",
    "\n",
    "Running this file requires the Python version of CompEcon. This can be installed with pip by running\n",
    "\n",
    "    !pip install compecon --upgrade\n",
    "\n",
    "<i>Last updated: 2022-Sept-04</i>\n",
    "<hr>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from compecon import LCP, tic, toc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate problem test data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n = 8\n",
    "z = np.random.randn(n, 2) - 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Boundaries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "a = np.min(z, 1)\n",
    "b = np.max(z, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Objective function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "q = np.random.randn(n)\n",
    "M = np.random.randn(n, n)\n",
    "M = - np.dot(M.T, M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define the problem by creating an LCP instance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "L = LCP(M, q, a, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set 100 random initial points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nrep = 100\n",
    "x0 = np.random.randn(nrep, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve by applying Newton method to semi-smooth formulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t0 = tic()\n",
    "it1 = 0\n",
    "L.opts.transform = 'ssmooth'\n",
    "for k in range(nrep):\n",
    "    L.newton(x0[k])\n",
    "    it1 += L.it\n",
    "t1 = toc(t0)\n",
    "n1 = L.fnorm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solve by applying Newton method to minmax formulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "t0 = tic()\n",
    "it2 = 0\n",
    "L.opts.transform = 'minmax'\n",
    "for k in range(nrep):\n",
    "    L.newton(x0[k])\n",
    "    it2 += L.it\n",
    "t2 = toc(t0)\n",
    "n2 = L.fnorm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Print results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hundredths of seconds required to solve randomly generated linear \n",
      " complementarity problem on R^8 using Newton and Lemke methods\n",
      "\n",
      "Algorithm           Time      Norm   Iterations  Iters/second\n",
      "------------------------------------------------------------\n",
      "Newton semismooth   0.78     1e-08       1447    1858.1\n",
      "Newton minmax       0.13     1e-16        365    2808.7\n"
     ]
    }
   ],
   "source": [
    "print('Hundredths of seconds required to solve randomly generated linear \\n',\n",
    "      'complementarity problem on R^8 using Newton and Lemke methods')\n",
    "print('\\nAlgorithm           Time      Norm   Iterations  Iters/second\\n' + '-' * 60)\n",
    "print('Newton semismooth {:6.2f}  {:8.0e}   {:8d}  {:8.1f}'.format(t1, n1, it1, it1/t1))\n",
    "print('Newton minmax     {:6.2f}  {:8.0e}   {:8d}  {:8.1f}'.format(t2, n2, it2, it2/t2))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.7 ('base')",
   "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.9.7"
  },
  "vscode": {
   "interpreter": {
    "hash": "ad2bdc8ecc057115af97d19610ffacc2b4e99fae6737bb82f5d7fb13d2f2c186"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}