{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Module biogeme.algorithms "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examples of use of each function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This webpage is for programmers who need examples of use of the functions of the class. The examples are designed to illustrate the syntax. They do not correspond to any meaningful model. For examples of models, visit  [biogeme.epfl.ch](http://biogeme.epfl.ch)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Important note: the functions in the module ``algorithms`` belonged to the module ``optimization`` in version 3.2.6 of Biogeme. Since version 3.2.7, the generic optimization algorithms are contained in the module ``algorithms`` while the module ``optimization`` contains function that can be directly used for the estimation by Biogeme.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2021-08-04 16:35:52.221727\n"
     ]
    }
   ],
   "source": [
    "import datetime\n",
    "print(datetime.datetime.now())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "biogeme 3.2.8 [2021-08-04]\n",
      "Version entirely written in Python\n",
      "Home page: http://biogeme.epfl.ch\n",
      "Submit questions to https://groups.google.com/d/forum/biogeme\n",
      "Michel Bierlaire, Transport and Mobility Laboratory, Ecole Polytechnique Fédérale de Lausanne (EPFL)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import biogeme.version as ver\n",
    "print(ver.getText())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import biogeme.algorithms as algo\n",
    "import biogeme.biogeme as bio\n",
    "import biogeme.database as db\n",
    "import biogeme.models as models\n",
    "from biogeme.expressions import Beta, Variable\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Defne the verbosity of Biogeme"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import biogeme.messaging as msg\n",
    "logger = msg.bioMessage()\n",
    "logger.setSilent()\n",
    "#logger.setDetailed()\n",
    "#logger.setDebug()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modified Cholesky factorization by Schnabel and Eskow (1999)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example by Eskow and Schnabel, 1991"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.3571, -0.103 ,  0.0274, -0.0459],\n",
       "       [-0.103 ,  0.2525,  0.0736, -0.3845],\n",
       "       [ 0.0274,  0.0736,  0.234 , -0.2878],\n",
       "       [-0.0459, -0.3845, -0.2878,  0.5549]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[ 0.3571, -0.1030,  0.0274, -0.0459],\n",
    "              [-0.1030,  0.2525,  0.0736, -0.3845],\n",
    "              [ 0.0274,  0.0736,  0.2340, -0.2878],\n",
    "              [-0.0459, -0.3845, -0.2878,  0.5549]])\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.7449161 ,  0.        ,  0.        ,  0.        ],\n",
       "       [-0.06161768,  0.59439319,  0.        ,  0.        ],\n",
       "       [-0.38635223,  0.00604629,  0.46941286,  0.        ],\n",
       "       [-0.51616551, -0.22679419, -0.26511936,  0.00152451]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L, E, P = algo.schnabelEskow(A)\n",
    "L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The factor $L$ is such that $A + E = PLL^TP^T$. Therefore, the expression below should be the null matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-5.55111512e-17,  0.00000000e+00,  0.00000000e+00,\n",
       "         0.00000000e+00],\n",
       "       [ 0.00000000e+00, -5.55111512e-17,  1.38777878e-17,\n",
       "         0.00000000e+00],\n",
       "       [ 0.00000000e+00,  1.38777878e-17,  0.00000000e+00,\n",
       "         0.00000000e+00],\n",
       "       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
       "         0.00000000e+00]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P @ L @ L.T @ P.T - E - A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example by Schnabel and Eskow (1999)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1890.3, -1705.6,  -315.8,  3000.3],\n",
       "       [-1705.6,  1538.3,   284.9, -2706.6],\n",
       "       [ -315.8,   284.9,    52.5,  -501.2],\n",
       "       [ 3000.3, -2706.6,  -501.2,  4760.8]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A2 = np.array([[ 1890.3, -1705.6, -315.8,  3000.3],\n",
    "               [-1705.6,  1538.3,  284.9, -2706.6],\n",
    "               [ -315.8,   284.9,   52.5,  -501.2],\n",
    "               [ 3000.3, -2706.6, -501.2,  4760.8]])\n",
    "A2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 6.89985507e+01,  0.00000000e+00,  0.00000000e+00,\n",
       "         0.00000000e+00],\n",
       "       [-7.26392069e+00,  3.19413321e-01,  0.00000000e+00,\n",
       "         0.00000000e+00],\n",
       "       [-3.92269109e+01, -1.28891153e-01,  4.44731720e-01,\n",
       "         0.00000000e+00],\n",
       "       [ 4.34835220e+01,  1.90522168e-01,  3.34584739e-01,\n",
       "         1.71484739e-03]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L, E, P = algo.schnabelEskow(A2)\n",
    "L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The factor $L$ is such that $A + E = PLL^TP^T$. Therefore, the expression below should be the null matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-2.27373675e-13,  0.00000000e+00,  0.00000000e+00,\n",
       "         0.00000000e+00],\n",
       "       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
       "         0.00000000e+00],\n",
       "       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
       "         0.00000000e+00],\n",
       "       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
       "        -9.09494702e-13]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P @ L @ L.T @ P.T - E - A2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rosenbrock problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "class rosenbrock(algo.functionToMinimize):\n",
    "    def __init__(self):\n",
    "        self.x = None\n",
    "        \n",
    "    def setVariables(self, x):\n",
    "        self.x = x\n",
    "\n",
    "    def f(self, batch=None):\n",
    "        if batch is not None:\n",
    "            raise excep.biogemeError('This function is not '\n",
    "                                     'data driven.')\n",
    "        n = len(self.x)\n",
    "        f = sum(100.0 * (self.x[i + 1]-self.x[i]**2)**2\n",
    "                + (1.0 - self.x[i])**2 for i in range(n - 1))\n",
    "        return f\n",
    "\n",
    "    def g(self):\n",
    "        n = len(self.x)\n",
    "        g = np.zeros(n)\n",
    "        for i in range(n - 1):\n",
    "            g[i] = g[i] - 400 * self.x[i] * \\\n",
    "                (self.x[i + 1] -self.x[i]**2) - \\\n",
    "                2 * (1 - self.x[i])\n",
    "            g[i + 1] = g[i + 1] + \\\n",
    "                200 * (self.x[i + 1] - self.x[i]**2)\n",
    "        return g\n",
    "    \n",
    "    def h(self):\n",
    "        n = len(self.x)\n",
    "        H = np.zeros((n, n))\n",
    "        for i in range(n - 1):\n",
    "            H[i][i] = H[i][i] - 400 * self.x[i + 1] \\\n",
    "                + 1200 * self.x[i]**2 + 2\n",
    "            H[i+1][i] = H[i+1][i] - 400 * self.x[i]\n",
    "            H[i][i+1] = H[i][i+1] - 400 * self.x[i]\n",
    "            H[i+1][i+1] = H[i+1][i+1] + 200\n",
    "        return H \n",
    "\n",
    "    def f_g(self, batch=None):\n",
    "        if batch is not None:\n",
    "            raise excep.biogemeError('This function is not '\n",
    "                                     'data driven.')\n",
    "        return self.f(), self.g()\n",
    "\n",
    "    def f_g_h(self, batch=None):\n",
    "        if batch is not None:\n",
    "            raise excep.biogemeError('This function is not '\n",
    "                                     'data driven.')\n",
    "        return self.f(), self.g(), self.h()\n",
    "    \n",
    "    def f_g_bhhh(self, batch=None):\n",
    "        raise excep.biogemeError('This function is not data driven.')\n",
    "\n",
    "theFunction = rosenbrock()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 5.8 from Bierlaire (2015)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is an infinite number of local optima. For each integer $k$,\n",
    "\n",
    "\\\\[x^* = ((-1)^{k+1}, k \\pi) \\\\]\n",
    "\n",
    "is a local optimum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "class example58(algo.functionToMinimize):\n",
    "    def __init__(self):\n",
    "        self.x = None\n",
    "        \n",
    "    def setVariables(self, x):\n",
    "        self.x = x\n",
    "\n",
    "    def f(self, batch=None):\n",
    "        if batch is not None:\n",
    "            raise excep.biogemeError('This function is not '\n",
    "                                     'data driven.')\n",
    "        n = len(self.x)\n",
    "        f = 0.5 * self.x[0] * self.x[0] + \\\n",
    "            self.x[0] * np.cos(self.x[1])\n",
    "        return f\n",
    "\n",
    "    def g(self):\n",
    "        n = len(self.x)\n",
    "        g = np.array([self.x[0] + np.cos(self.x[1]),\n",
    "                      -self.x[0]*np.sin(self.x[1])])\n",
    "        return g\n",
    "\n",
    "    def h(self):\n",
    "        n = len(self.x)\n",
    "        H = np.array([[1, \n",
    "                       -np.sin(self.x[1])],\n",
    "                      [-np.sin(self.x[1]),\n",
    "                       -self.x[0] * np.cos(self.x[1])]])\n",
    "        return H \n",
    "\n",
    "    def f_g(self, batch=None):\n",
    "        if batch is not None:\n",
    "            raise excep.biogemeError('This function is not '\n",
    "                                     'data driven.')\n",
    "        return self.f(), self.g()\n",
    "\n",
    "    def f_g_h(self, batch=None):\n",
    "        if batch is not None:\n",
    "            raise excep.biogemeError('This function is not '\n",
    "                                     'data driven.')\n",
    "        return self.f(), self.g(), self.h()\n",
    "    \n",
    "    def f_g_bhhh(self, batch=None):\n",
    "        raise excep.biogemeError('This function is not data driven.')\n",
    "\n",
    "ex58 = example58()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Line search algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha=0.0009765625 nfev=12\n"
     ]
    }
   ],
   "source": [
    "x = np.array([-1.5, 1.5])\n",
    "theFunction.setVariables(x)\n",
    "f, g = theFunction.f_g()\n",
    "alpha, nfev = algo.lineSearch(theFunction, x, f, g, -g)\n",
    "print(f\"alpha={alpha} nfev={nfev}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Newton with linesearch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rosenbrock"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 1.])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = np.array([-1.5, 1.5])\n",
    "xstar, messages = algo.newtonLineSearch(theFunction, x0)\n",
    "xstar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Algorithm:\tUnconstrained Newton with line search\n",
      "Relative gradient:\t5.571454408920588e-11\n",
      "Number of iterations:\t23\n",
      "Number of function evaluations:\t80\n",
      "Number of gradient evaluations:\t80\n",
      "Number of hessian evaluations:\t24\n",
      "Cause of termination:\tRelative gradient = 5.6e-11 <= 6.1e-06\n"
     ]
    }
   ],
   "source": [
    "for k, v in messages.items():\n",
    "    print(f'{k}:\\t{v}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 5.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.        , 3.14159265])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = np.array([1, 1])\n",
    "xstar, messages = algo.newtonLineSearch(ex58, x0)\n",
    "xstar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Algorithm:\tUnconstrained Newton with line search\n",
      "Relative gradient:\t2.202540372794277e-10\n",
      "Number of iterations:\t5\n",
      "Number of function evaluations:\t28\n",
      "Number of gradient evaluations:\t28\n",
      "Number of hessian evaluations:\t6\n",
      "Cause of termination:\tRelative gradient = 2.2e-10 <= 6.1e-06\n"
     ]
    }
   ],
   "source": [
    "for k, v in messages.items():\n",
    "    print(f'{k}:\\t{v}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Newton with trust region"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rosenbrock"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.99999995, 0.99999989])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = np.array([-1.5, 1.5])\n",
    "xstar, messages = algo.newtonTrustRegion(theFunction, x0)\n",
    "xstar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Algorithm:\tUnconstrained Newton with trust region\n",
      "Relative gradient:\t9.752977092603033e-10\n",
      "Cause of termination:\tRelative gradient = 9.8e-10 <= 6.1e-06\n",
      "Number of iterations:\t28\n",
      "Number of function evaluations:\t49\n",
      "Number of gradient evaluations:\t22\n",
      "Number of hessian evaluations:\t22\n"
     ]
    }
   ],
   "source": [
    "for k, v in messages.items():\n",
    "    print(f'{k}:\\t{v}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 5.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.00000000e+00, -1.56439954e-09])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = np.array([1.0, 1.0])\n",
    "xstar, messages = algo.newtonTrustRegion(ex58, x0)\n",
    "xstar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Algorithm:\tUnconstrained Newton with trust region\n",
      "Relative gradient:\t1.5037932097369974e-09\n",
      "Cause of termination:\tRelative gradient = 1.5e-09 <= 6.1e-06\n",
      "Number of iterations:\t5\n",
      "Number of function evaluations:\t10\n",
      "Number of gradient evaluations:\t6\n",
      "Number of hessian evaluations:\t6\n"
     ]
    }
   ],
   "source": [
    "for k, v in messages.items():\n",
    "    print(f'{k}:\\t{v}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BFGS and linesearch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rosenbrock"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.99999897, 0.99999797])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = np.array([-1.5, 1.5])\n",
    "xstar, messages = algo.bfgsLineSearch(theFunction, x0, maxiter=10000)\n",
    "xstar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Algorithm:\tInverse BFGS with line search\n",
      "Relative gradient:\t1.531551364614023e-07\n",
      "Cause of termination:\tRelative gradient = 1.5e-07 <= 6.1e-06\n",
      "Number of iterations:\t32\n",
      "Number of function evaluations:\t152\n",
      "Number of gradient evaluations:\t33\n"
     ]
    }
   ],
   "source": [
    "for k, v in messages.items():\n",
    "    print(f'{k}:\\t{v}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 5.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.00000142e+00,  3.57636862e-06])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = np.array([1, 1])\n",
    "xstar, messages = algo.bfgsLineSearch(ex58, x0, maxiter=10000)\n",
    "xstar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Algorithm:\tInverse BFGS with line search\n",
      "Relative gradient:\t3.4378215567736776e-06\n",
      "Cause of termination:\tRelative gradient = 3.4e-06 <= 6.1e-06\n",
      "Number of iterations:\t10\n",
      "Number of function evaluations:\t48\n",
      "Number of gradient evaluations:\t11\n"
     ]
    }
   ],
   "source": [
    "for k, v in messages.items():\n",
    "    print(f'{k}:\\t{v}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BFGS and trust region"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rosenbrock"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.00000052, 1.00000078])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = np.array([-1.5, 1.5])\n",
    "xstar, messages = algo.bfgsTrustRegion(theFunction, x0, maxiter=10000)\n",
    "xstar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Algorithm:\tBFGS with trust region\n",
      "Relative gradient:\t1.7562432530828692e-06\n",
      "Cause of termination:\tRelative gradient = 1.8e-06 <= 6.1e-06\n",
      "Number of iterations:\t50\n",
      "Number of function evaluations:\t88\n",
      "Number of gradient evaluations:\t38\n"
     ]
    }
   ],
   "source": [
    "for k, v in messages.items():\n",
    "    print(f'{k}:\\t{v}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 5.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-9.99999972e-01,  1.58776353e-08])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = np.array([1, 1])\n",
    "xstar, messages = algo.bfgsTrustRegion(ex58, x0, maxiter=10000)\n",
    "xstar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Algorithm:\tBFGS with trust region\n",
      "Relative gradient:\t2.7200302000474536e-08\n",
      "Cause of termination:\tRelative gradient = 2.7e-08 <= 6.1e-06\n",
      "Number of iterations:\t17\n",
      "Number of function evaluations:\t30\n",
      "Number of gradient evaluations:\t13\n"
     ]
    }
   ],
   "source": [
    "for k, v in messages.items():\n",
    "    print(f'{k}:\\t{v}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Newton for simple bounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Rosenbrock"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bound constraints do not exclude the unconstrained optimum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.99999961, 0.99999918])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = np.array([-1.5, 1.5])\n",
    "theBounds = algo.bioBounds([(-1000, 1000), (-1000, 1000)])\n",
    "xstar, messages = algo.simpleBoundsNewtonAlgorithm(theFunction, \n",
    "                                                  theBounds,\n",
    "                                                  x0)\n",
    "xstar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Algorithm:\tNewton with trust region for simple bound constraints\n",
      "Proportion analytical hessian:\t100.0%\n",
      "Relative projected gradient:\t2.0646130807833174e-07\n",
      "Relative change:\t0.00038043149628330664\n",
      "Number of iterations:\t32\n",
      "Number of function evaluations:\t77\n",
      "Number of gradient evaluations:\t23\n",
      "Number of hessian evaluations:\t23\n",
      "Cause of termination:\tRelative gradient = 2.1e-07 <= 6.1e-06\n"
     ]
    }
   ],
   "source": [
    "for k, v in messages.items():\n",
    "    print(f'{k}:\\t{v}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the bound constraints do exclude the unconstrained optimum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.99497475,  1.        ])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = np.array([-1.5, 1.5])\n",
    "theBounds = algo.bioBounds([(-1000, 0), (1, 1000)])\n",
    "xstar, messages = algo.simpleBoundsNewtonAlgorithm(theFunction, \n",
    "                                                  theBounds,\n",
    "                                                  x0)\n",
    "xstar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Algorithm:\tNewton with trust region for simple bound constraints\n",
      "Proportion analytical hessian:\t100.0%\n",
      "Relative projected gradient:\t7.676708776216401e-09\n",
      "Relative change:\t2.0046022513153794e-05\n",
      "Number of iterations:\t7\n",
      "Number of function evaluations:\t20\n",
      "Number of gradient evaluations:\t7\n",
      "Number of hessian evaluations:\t7\n",
      "Cause of termination:\tRelative gradient = 7.7e-09 <= 6.1e-06\n"
     ]
    }
   ],
   "source": [
    "for k, v in messages.items():\n",
    "    print(f'{k}:\\t{v}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 5.8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bound constraints do not exclude the unconstrained optimum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.        , 3.14159265])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = np.array([1, 1])\n",
    "theBounds = algo.bioBounds([(-1000, 1000), (-1000, 1000)])\n",
    "xstar, messages = algo.simpleBoundsNewtonAlgorithm(ex58, \n",
    "                                                  theBounds,\n",
    "                                                  x0)\n",
    "xstar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Algorithm:\tNewton with trust region for simple bound constraints\n",
      "Proportion analytical hessian:\t100.0%\n",
      "Relative projected gradient:\t2.389957069442179e-11\n",
      "Relative change:\t2.245052812572204e-06\n",
      "Number of iterations:\t6\n",
      "Number of function evaluations:\t19\n",
      "Number of gradient evaluations:\t7\n",
      "Number of hessian evaluations:\t7\n",
      "Cause of termination:\tRelative change = 2.25e-06 <= 1e-05\n"
     ]
    }
   ],
   "source": [
    "for k, v in messages.items():\n",
    "    print(f'{k}:\\t{v}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bound constraints do exclude the unconstrained optimum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.54030231,  1.        ])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0 = np.array([-1, 0])\n",
    "theBounds = algo.bioBounds([(-1000, 0), (1, 1000)])\n",
    "xstar, messages = algo.simpleBoundsNewtonAlgorithm(ex58, \n",
    "                                                  theBounds,\n",
    "                                                  x0)\n",
    "xstar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Algorithm:\tNewton with trust region for simple bound constraints\n",
      "Proportion analytical hessian:\t100.0%\n",
      "Relative projected gradient:\t8.881784197001252e-16\n",
      "Relative change:\t0.4596976941318611\n",
      "Number of iterations:\t1\n",
      "Number of function evaluations:\t4\n",
      "Number of gradient evaluations:\t2\n",
      "Number of hessian evaluations:\t2\n",
      "Cause of termination:\tRelative gradient = 8.9e-16 <= 6.1e-06\n"
     ]
    }
   ],
   "source": [
    "for k, v in messages.items():\n",
    "    print(f'{k}:\\t{v}')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}