{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# First steps with the nullspace optimizer\n",
    "\n",
    "This notebook proposes you to familiar you with the nullspace optimization algorithm, which in the case of equality constrained optimization, reduces to the flow proposed by Yamashita in \n",
    "\n",
    "*H. Yamashita, A differential equation approach to nonlinear programming. Math. Program.18(1980) 155–168.*\n",
    "\n",
    "The flow reads \n",
    "\n",
    "\\begin{equation}\n",
    "\\newcommand{\\DD}{\\mathrm{D}}\n",
    "\\renewcommand{\\dagger}{T}\n",
    "\\renewcommand{\\x}{x}\n",
    "\\renewcommand{\\I}{I}\n",
    "\\renewcommand{\\g}{g}\n",
    "    \\label{eqn:flow}\n",
    "    \\left\\{\\begin{aligned}\n",
    "            \\dot\\x&={ -\\alpha_J(\\I-\\DD\\g^\\dagger(\\DD\\g\\DD\\g^ \\dagger)^{-1}\\DD\\g(x))\\nabla\n",
    "            J(x) }{ -\\alpha_C\\DD\\g^ \\dagger(\\DD\\g\\DD\\g^ \\dagger)^{-1}\\g(x) }\\\\\n",
    "    \\x(0)&=x_0\n",
    "        \\end{aligned}\\right.\n",
    "\\end{equation}\n",
    "\n",
    "where $J$ is the cost function and $g$ the equality constraint:\n",
    "\\begin{equation}\n",
    "    \\label{eqn:nlsvg}\n",
    "\\begin{aligned}\n",
    "\t\\min_{\\x\\in \\mathbb{R}^n}&  \\quad J(\\x)\\\\\n",
    "\t\\textrm{s.t.} & \\left.\\begin{aligned}\n",
    " \\g(\\x)&=0\n",
    "\t\t\\end{aligned}\\right.\n",
    "\\end{aligned}\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. A first optimization program"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write an optimization program to solve the constrained minimization problem on the hyperbola:\n",
    "            $$\n",
    "            \\begin{aligned}\n",
    "                \\min_{(x_1,x_2)\\in\\mathbb{R}^{2}} & \\qquad x_1+x_2\\\\\n",
    "                s.t. & \\left.\\begin{aligned}\n",
    "  x_1x_2 &= 1. \n",
    "                \\end{aligned}\\right.\n",
    "            \\end{aligned}\n",
    "$$\n",
    "\n",
    "In order to solve the optimization problem, we use the function `nlspace_solve` which solves the above optimization program and whose prototype is \n",
    "```python\n",
    "from nullspace_optimizer import nlspace_solve\n",
    "\n",
    "# Define problem\n",
    "\n",
    "results = nlspace_solve(problem: Optimizable, params=None, results=None)\n",
    "```\n",
    "The input variables are \n",
    "-  `problem` : an `Optimizable` object described below. This variable contains all the information about the optimization problem to solve (objective and constraint functions, derivatives...)\n",
    "- `params`  : (optional) a dictionary containing algorithm parameters.\n",
    "\n",
    "- `results` : (optional) a previous output of the `nlspace_solve` function.         The optimization will then keep going from the last input of           the dictionary `results['x'][-1]`. This is useful when one needs to restart an optimization after an interruption.\n",
    "\n",
    "The optimization routine `nlspace_solve` returns the dictionary `opt_results` which contains various information about the optimization path, including the values of the optimization variables `results['x']`.\n",
    "    \n",
    "    \n",
    " In our particular optimization test case in $\\mathbb R^n$ with $n=2$, we use the  `EuclideanOptimizable` class which inherits `Optimizable` which simplifies the definition of the optimization program (the inner product is specified by default). \n",
    "We allow the user to specify the initialization in the constructor `__init__`.\n",
    "\n",
    "Fill in the definition below with the right values for solving the above optimization program.\n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nullspace_optimizer import nlspace_solve, EuclideanOptimizable\n",
    "import numpy as np\n",
    "\n",
    "class problem1(EuclideanOptimizable):\n",
    "    def __init__(self,x0):\n",
    "        super().__init__(2)\n",
    "        self.xinit = x0\n",
    "        self.nconstraints = #to fill\n",
    "        self.nineqconstraints = # to fill\n",
    "\n",
    "    def x0(self):\n",
    "        return self.xinit\n",
    "\n",
    "    def J(self, x):\n",
    "        return #to fill\n",
    "\n",
    "    def dJ(self, x):\n",
    "        return #to fill\n",
    "\n",
    "    def G(self, x):\n",
    "        return #to fill\n",
    "\n",
    "    def dG(self, x):\n",
    "        return #to fill"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then write a routine to solve the problem several time with several initializations. Adapt the code below to specify the entries (0.1,0.1), (4.0,0.25), (4,1). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_problems():\n",
    "    xinits = #to fill \n",
    "    # Write xinits\n",
    "    problems = [problem1(x0=x0) for x0 in xinits]\n",
    "    params = {'dt': 0.1,  'alphaJ':2, 'alphaC':1, 'debug': -1}\n",
    "    return [nlspace_solve(pb, params) for pb in problems]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function nlspace_solve calls the null space optimizer on this problem. Get the result by calling the function run_problems:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = run_problems()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variable `results` contains various data about the optimization trajectories, e.g.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have access to the objective function and constraint histories, and the value of the Lagrange multiplier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i=1\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(results[i]['it'],results[i]['J'])\n",
    "plt.title('J')\n",
    "plt.figure()\n",
    "plt.plot(results[i]['it'],[ g[0] for g in results[i]['G']])\n",
    "plt.title('G')\n",
    "plt.figure()\n",
    "plt.plot(results[i]['it'][:-1],[mu[0] for mu in results[i]['muls']])\n",
    "plt.title('$\\lambda$');\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comment on the numerical values of the constraint G. Observe that the objective function keeps decreasing while maintaining the constraint satisfied.\n",
    "\n",
    "Plot the constraint $x_1x_2=1$ and the optimization trajectories:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nullspace_optimizer.examples.draw as draw\n",
    "t = np.linspace(1/3,5,100)\n",
    "plt.plot(t,1/t,color='red',linewidth=0.5,label=\"y=1/x\")\n",
    "plt.plot(-t,-1/t,color='red',linewidth=0.5)\n",
    "for i, r in enumerate(results):\n",
    "    draw.drawData(r, f'x{i}', f'C{i}', x0=True, xfinal=True, initlabel=None)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Further works\n",
    "\n",
    "1. Adapt the above codes to display other trajectories.\n",
    "2. Try with other initializations, change the parameter alphaC and alphaJ\n",
    "3. Comment on the trajectories\n",
    "4. Try another equality constrained optimization program below"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Another problem\n",
    "\n",
    "Do the same to solve \n",
    "\n",
    "$$\n",
    "        \\begin{aligned}\n",
    "                \\max_{(x_1,x_2)\\in\\mathbb{R}^{2}} & \\qquad x_2\\\\\n",
    "                s.t. & \\left\\{\\begin{aligned}\n",
    "                    (x_1-0.5)^{2}+x_2^{2} &= 2\\\\\n",
    "                    (x_1+0.5)^{2}+x_2^{2} &= 2.\n",
    "                \\end{aligned}\\right.\n",
    "            \\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class problem2(EuclideanOptimizable):\n",
    "    def __init__(self,x0):\n",
    "        super().__init__(2)\n",
    "        self.xinit = x0\n",
    "        self.nconstraints = 2\n",
    "        self.nineqconstraints = 0\n",
    "\n",
    "    def x0(self):\n",
    "        return self.xinit\n",
    "\n",
    "    def J(self, x):\n",
    "        return -x[1]\n",
    "\n",
    "    def dJ(self, x):\n",
    "        return [0, -1]\n",
    "\n",
    "    def G(self, x):\n",
    "        return [(x[0]-0.5)**2+x[1]**2-2,(x[0]+0.5)**2+x[1]**2-2]\n",
    "\n",
    "    def dG(self, x):\n",
    "        return [[2*(x[0]-0.5),2*x[1]],[2*(x[0]+0.5),2*x[1]]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def run_problems():\n",
    "    xinits = ([0,0], [1.0, 0],  [3, -1], [-1.5,-0.5])\n",
    "    # Write xinits\n",
    "    problems = [problem2(x0=x0) for x0 in xinits]\n",
    "    params = {'dt': 0.05,  'alphaJ':2, 'alphaC':1, 'debug': -1}\n",
    "    return [nlspace_solve(pb, params) for pb in problems]\n",
    "\n",
    "results=run_problems()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t=np.linspace(0,2*np.pi,100)\n",
    "plt.plot(0.5+np.sqrt(2)*np.cos(t),np.sqrt(2)*np.sin(t),color='red',linewidth='0.5')\n",
    "plt.plot(-0.5+np.sqrt(2)*np.cos(t),np.sqrt(2)*np.sin(t),color='red',linewidth='0.5')\n",
    "plt.axis('equal')\n",
    "for i, r in enumerate(results):\n",
    "    draw.drawData(r, f'x{i}', f'C{i}', x0=True, xfinal=True, initlabel=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is the behavior of the trajectories `x2` and `x3` surprising ?"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}