{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Interior Point\n",
    "For this training course we'll be using the optimization routines in SciPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.optimize as opt\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need to define a callback function so that we can plot the optimization steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def callback(x):\n",
    "    global xprev\n",
    "    plt.plot([xprev[0],x[0]],[xprev[1],x[1]],'b.-')\n",
    "    xprev = x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.1 Logarithmic Barrier\n",
    "\n",
    "Recall that the logarithmic barrier function, for a single inequality constraint $c$, is given by\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "\\Phi_{\\mu}(x) = f(x) - \\mu \\log c(x)\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "where $\\mu >0$ is a barrier parameter. The gradient for the barrier function is\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "\\nabla \\Phi_{\\mu}(x) = \\nabla f(x) - \\mu\\frac{1}{c(x)}\\nabla c(x)\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "and the Hessian is given by\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "\\nabla^2 \\Phi_{\\mu}(x) = \\nabla^2 f(x) - \\mu\\frac{1}{c(x)}\\nabla^2 c(x)  +\\mu\\frac{1}{c(x)^2}\\nabla c(x)\\nabla c(x)^T\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "#### Coding Task:\n",
    "\n",
    "For different fixed values of $\\mu$, apply your favourite SciPy optimization routine to minimize $\\Phi_{\\mu}(x)$ for the following objectives:\n",
    "\n",
    "1. $f(x) = \\frac{1}{2}(ax_1^2 + x_2^2) \\qquad\\qquad\\qquad$ with $\\; c(x) = x_1 + x_2^2 - 1$ \n",
    "\n",
    "2. $f(x) = (a - x_1)^2 + b(x_2 - x_1^2)^2 \\qquad$ with $\\; c(x) = x_1 + x_2^2 - 1$ \n",
    "\n",
    "starting at $(x_1,x_2) = (1,1)$ for suitable choices of parameters $a,b$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "a = 1\n",
    "mu = \n",
    "\n",
    "# Initial guess\n",
    "x0 = np.array(...)\n",
    "\n",
    "# Objective function, gradient and Hessian\n",
    "obj = lambda x: \n",
    "gobj = lambda x: \n",
    "hobj = lambda x: \n",
    "\n",
    "# Constraint function, gradient and Hessian\n",
    "con = lambda x:\n",
    "gcon = lambda x: \n",
    "hcon = lambda x:\n",
    "\n",
    "# Barrier function, gradient and Hessian\n",
    "fun = lambda x: \n",
    "jac = lambda x:\n",
    "hess = lambda x: \n",
    "\n",
    "# Plot function contours\n",
    "plt.figure()\n",
    "X = np.linspace(...)\n",
    "Y = np.linspace(...)\n",
    "Z = np.meshgrid(X,Y)\n",
    "plt.contour(...)\n",
    "plt.colorbar()\n",
    "\n",
    "# Plot constraint\n",
    "plt.contour(X,Y,con(Z),[0],colors='r')\n",
    "\n",
    "# Call SciPy's Newton-CG\n",
    "xprev = x0 # for plotting\n",
    "res = opt.minimize(..., method=..., callback=callback)\n",
    "\n",
    "# Print results and show plot\n",
    "print(res)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "a = 1\n",
    "b = 100\n",
    "mu = \n",
    "\n",
    "# Initial guess\n",
    "x0 = np.array(...)\n",
    "\n",
    "# Objective function, gradient and Hessian\n",
    "obj = lambda x: \n",
    "gobj = lambda x: \n",
    "hobj = lambda x: \n",
    "\n",
    "# Constraint function, gradient and Hessian\n",
    "con = lambda x: \n",
    "gcon = lambda x: \n",
    "hcon = lambda x: \n",
    "\n",
    "# Barrier function, gradient and Hessian\n",
    "fun = lambda x: \n",
    "jac = lambda x: \n",
    "hess = lambda x: \n",
    "\n",
    "# Plot function contours\n",
    "plt.figure()\n",
    "X = np.linspace(...)\n",
    "Y = np.linspace(...)\n",
    "Z = np.meshgrid(X,Y)\n",
    "plt.contour(...)\n",
    "plt.colorbar()\n",
    "\n",
    "# Plot constraint\n",
    "plt.contour(X,Y,con(Z),[0],colors='r')\n",
    "\n",
    "# Call SciPy's Newton-CG\n",
    "xprev = x0 # for plotting\n",
    "res = opt.minimize(..., method=..., callback=callback)\n",
    "\n",
    "# Print results and show plot\n",
    "print(res)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.2 Central Paths\n",
    "\n",
    "Recall that the perturbed first-order optimality conditions for a single inequality constrained problem are given by\n",
    "\n",
    "$\n",
    "\\begin{align}\n",
    "\\nabla f(x) - \\nabla c(x) y &= 0 \\\\\n",
    "c(x)y &= \\mu\n",
    "\\end{align}\n",
    "$\n",
    "\n",
    "for $c(x) > 0$ and $y > 0$.\n",
    "\n",
    "Interior point methods track solutions $(x^*,y^*)$ to this nonlinear system of equations whilst maintaining feasibility (i.e. $c(x^*), y^* > 0$) as $\\mu$ shrinks to zero. These correspond to solutions in the interior of the feasible domain (hence the name) which follow a \"central path\" through the feasible domain.\n",
    "\n",
    "#### Coding Task:\n",
    "\n",
    "For progressively smaller values of $\\mu$, use Scipy's root finding to find and plot solutions to the above nonlinear system for the following objectives:\n",
    "\n",
    "1. $f(x) = \\frac{1}{2}(ax_1^2 + x_2^2) \\qquad\\qquad\\qquad$ with $\\; c(x) = x_1 + x_2^2 - 1$ \n",
    "\n",
    "2. $f(x) = (a - x_1)^2 + b(x_2 - x_1^2)^2 \\qquad$ with $\\; c(x) = x_1 + x_2^2 - 1$ \n",
    "\n",
    "starting at $(x_1,x_2,y) = (1,1,1)$ for suitable choices of parameters $a,b$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "a = 1\n",
    "mus = np.arange(...)\n",
    "\n",
    "# Initial guess\n",
    "x0 = np.array(...)\n",
    "y0 = \n",
    "\n",
    "# Objective function, gradient and Hessian\n",
    "obj = lambda x: \n",
    "gobj = lambda x: \n",
    "hobj = lambda x: \n",
    "\n",
    "# Constraint function, gradient and Hessian\n",
    "con = lambda x: \n",
    "gcon = lambda x: \n",
    "hcon = lambda x: \n",
    "\n",
    "# Nonlinear system as a vector function and its Jacobian\n",
    "sys = lambda x,y,mu: \n",
    "Jsys = lambda x,y,mu: \n",
    "\n",
    "# Plot function contours\n",
    "plt.figure()\n",
    "X = np.linspace(...)\n",
    "Y = np.linspace(...)\n",
    "Z = np.meshgrid(X,Y)\n",
    "plt.contour(...)\n",
    "plt.colorbar()\n",
    "\n",
    "# Plot constraint\n",
    "plt.contour(X,Y,con(Z),[0],colors='r')\n",
    "\n",
    "# Call SciPy's root for each value of mu\n",
    "v0 = np.hstack((x0,y0)) # starting points \n",
    "xprev = None # for plotting\n",
    "for mu in mus:\n",
    "    \n",
    "    # nonlinear system and its Jacobian\n",
    "    fun = lambda v: sys(...)\n",
    "    jac = lambda v: Jsys(...)\n",
    "    \n",
    "    # find root of nonlinear system \n",
    "    res = opt.root(...)\n",
    "    \n",
    "    # plot root\n",
    "    if xprev is not None:\n",
    "        plt.plot([xprev[0],res.x[0]],[xprev[1],res.x[1]],'b-')\n",
    "    xprev = res.x # for plotting\n",
    "    \n",
    "# Show plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "a = 1\n",
    "b = 100\n",
    "mus = np.arange(...)\n",
    "\n",
    "# Initial guess\n",
    "x0 = np.array(...)\n",
    "y0 = \n",
    "\n",
    "# Objective function, gradient and Hessian\n",
    "obj = lambda x: \n",
    "gobj = lambda x: \n",
    "hobj = lambda x: \n",
    "\n",
    "# Constraint function, gradient and Hessian\n",
    "con = lambda x: \n",
    "gcon = lambda x: \n",
    "hcon = lambda x: \n",
    "\n",
    "# Nonlinear system as a vector function and its Jacobian\n",
    "sys = lambda x,y,mu: \n",
    "Jsys = lambda x,y,mu:\n",
    "\n",
    "# Plot function contours\n",
    "plt.figure()\n",
    "X = np.linspace(...)\n",
    "Y = np.linspace(...)\n",
    "Z = np.meshgrid(X,Y)\n",
    "plt.contour(...)\n",
    "plt.colorbar()\n",
    "\n",
    "# Plot constraint\n",
    "plt.contour(X,Y,con(Z),[0],colors='r')\n",
    "\n",
    "# Call SciPy's root for each value of mu\n",
    "v0 = np.hstack((x0,y0)) # starting points \n",
    "xprev = None # for plotting\n",
    "for mu in mus:\n",
    "    \n",
    "    # nonlinear system and its Jacobian\n",
    "    fun = lambda v: sys(...)\n",
    "    jac = lambda v: Jsys(...)\n",
    "    \n",
    "    # find root of nonlinear system \n",
    "    res = opt.root(...)\n",
    "    \n",
    "    # plot root\n",
    "    if xprev is not None:\n",
    "        plt.plot([xprev[0],res.x[0]],[xprev[1],res.x[1]],'b-')\n",
    "    xprev = res.x # for plotting\n",
    "    \n",
    "# Show plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}