{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solve a box constrained optimization problem\n",
    "\n",
    "## Notebook to demonstrate how to solve an optimization problem subject to box constraints via oppy\n",
    "\n",
    "We consider the minimization problem of the Mishra's Bird function with box constraints\n",
    "\\begin{align*}\n",
    " \\min_{x,y}f(x,y) = \\min_{x,y}\\sin(y)e^{[(1-\\cos(x))^2]}+\\cos(x)e^{[(1-\\sin(y))^2]}+(x-y)^2,\n",
    "\\end{align*}\n",
    "\n",
    "\\begin{align*}\n",
    "\\mathrm{s.t.}\\quad -10\\leq x\\leq 0 \\quad \\mathrm{and} \\quad -6.5\\leq y\\leq 0.\n",
    "\\end{align*}\n",
    " \n",
    "This problem is known as a test function for optimization (see [Wikipedia](https://en.wikipedia.org/wiki/Test_functions_for_optimization)). The global minimum is located at\n",
    "\n",
    "\\begin{align*}\n",
    "f(-3.1302468, -1.5821422)=-106.7645367\n",
    "\\end{align*}\n",
    "\n",
    "Firstly, we solve the problem with the projected gradient method with default linesearch i.e. with the Armijo Goldstein method. Therefore, the gradient of the Mishra's Bird function is needed in all cases. The Mishra's Bird and it's gradient can be imported from the $\\small{\\textbf{costfunctions.py}}$ file from the $\\small{\\textbf{tests}}$ folder in the oppy package. \n",
    "\n",
    "First we import all dependencies needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from oppy.options import Options\n",
    "from oppy.conOpt import projected_gradient\n",
    "\n",
    "from oppy.tests.costfunctions import mishrabird, gradient_mishrabird\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initial point, maximum number of iterations, lower and upper bound for the optimization methods are defined in the following: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initial point\n",
    "x0 = np.array([-4, -2])\n",
    "\n",
    "# Maximum number of iterations:\n",
    "nmax = 1000\n",
    "\n",
    "# lower boundaries\n",
    "low = np.array([-10, -6.5])\n",
    "# upper boundaries\n",
    "up = np.array([0, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can execute the optimization algorithm and extract the residual history of each result. Furthermore we turn on the display option which shows more information in each iteration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "result_proj_grad = projected_gradient(mishrabird, gradient_mishrabird, x0, low, up,\n",
    "                                         Options(disp=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The algorithm stops at a total number of iteration of 14. Another approach for faster convergence is to use another method. We try the L-BFGS-B method and thus import the algorithm. Since we use the same input as before, we can execute the algorithm right away."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from oppy.conOpt import l_bfgs_b\n",
    "\n",
    "result_l_bfgs_b = l_bfgs_b(mishrabird, gradient_mishrabird, x0, low, up,\n",
    "                                         Options(disp=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The total number of iterations is 8. For less iterations we try the projected Newton method. Therefore, we define the action of the Hessian of the Mishra's Bird which is required for the projected Newton method and import the projected Newton and the hessian of the Mishra's Bird function. Then we execute the projected Newton algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from oppy.conOpt import projected_newton\n",
    "\n",
    "from oppy.tests.costfunctions import hessian_mishrabird\n",
    "\n",
    "def action_mishrabird(u, x, Act, Inact):\n",
    "    \"\"\"Compute the action of the Hessian.\"\"\"\n",
    "    ret = Act*u + Inact*((hessian_mishrabird(x)@(Inact*u)))\n",
    "    return ret\n",
    "\n",
    "result_proj_newton = projected_newton(mishrabird, gradient_mishrabird, action_mishrabird, x0, low, up,\n",
    "                                         Options(disp=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The projected Newton stops after 6 iterations. This result is better than the other methods. This can be compared in a residual semi y-axis log plot of the norm of the gradient with respect to the iteration number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "residual_proj_grad = result_proj_grad.res\n",
    "residual_l_bfgs_b = result_l_bfgs_b.res\n",
    "residual_proj_newton = result_proj_newton.res\n",
    "\n",
    "\n",
    "# Plot norm of the gradient in semilogy plot\n",
    "plt.semilogy(residual_proj_grad, label='proj_grad')\n",
    "plt.semilogy(residual_l_bfgs_b, label='l_bfgs_b')\n",
    "plt.semilogy(residual_proj_newton, label='proj_newton')\n",
    "plt.title('Norm of the gradient, different schemes')\n",
    "plt.grid(True)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this plot we notice that the results with the Armijo Goldstein and the projected Armijo linesearch method coincide in each optimization algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(result_proj_grad.x)\n",
    "print(result_proj_newton.x)\n",
    "print(result_l_bfgs_b.x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results show that the solutions of the optimization methods with linesearch methods are close to the analytic global minima $(-3.1302468, -1.5821422)$, which shows a succesfull execution of the algorithms."
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}