ch06-code-listing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 6: Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Robert Johansson\n",
    "\n",
    "Source code listings for [Numerical Python - A Practical Techniques Approach for Industry](http://www.apress.com/9781484205549) (ISBN 978-1-484205-54-9).\n",
    "\n",
    "The source code listings can be downloaded from http://www.apress.com/9781484205549"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "matplotlib.rcParams['mathtext.fontset'] = 'stix'\n",
    "matplotlib.rcParams['font.family'] = 'serif'\n",
    "matplotlib.rcParams['font.sans-serif'] = 'stix'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import sympy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sympy.init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy import optimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import cvxopt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import division"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Univariate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "r, h = sympy.symbols(\"r, h\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Area = 2 * sympy.pi * r**2 + 2 * sympy.pi * r * h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Volume = sympy.pi * r**2 * h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "h_r = sympy.solve(Volume - 1)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Area_r = Area.subs(h_r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rsol = sympy.solve(Area_r.diff(r))[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rsol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "_.evalf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# verify that the second derivative is positive, so that rsol is a minimum\n",
    "Area_r.diff(r, 2).subs(r, rsol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Area_r.subs(r, rsol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "_.evalf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def f(r):\n",
    "    return 2 * np.pi * r**2 + 2 / r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "r_min = optimize.brent(f, brack=(0.1, 4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "r_min"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f(r_min)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "optimize.minimize_scalar(f, bracket=(0.1, 4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "r = np.linspace(0, 2, 100)[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 4))\n",
    "\n",
    "ax.plot(r, f(r), lw=2, color='b')\n",
    "ax.plot(r_min, f(r_min), 'r*', markersize=15)\n",
    "ax.set_title(r\"$f(r) = 2\\pi r^2+2/r$\", fontsize=18)\n",
    "ax.set_xlabel(r\"$r$\", fontsize=18)\n",
    "ax.set_xticks([0, 0.5, 1, 1.5, 2])\n",
    "ax.set_ylim(0, 30)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig('ch6-univariate-optimization-example.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two-dimensional"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x1, x2 = sympy.symbols(\"x_1, x_2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f_sym = (x1-1)**4 + 5 * (x2-1)**2 - 2*x1*x2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fprime_sym = [f_sym.diff(x_) for x_ in (x1, x2)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Gradient\n",
    "sympy.Matrix(fprime_sym)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fhess_sym = [[f_sym.diff(x1_, x2_) for x1_ in (x1, x2)] for x2_ in (x1, x2)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Hessian\n",
    "sympy.Matrix(fhess_sym)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f_lmbda = sympy.lambdify((x1, x2), f_sym, 'numpy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fprime_lmbda = sympy.lambdify((x1, x2), fprime_sym, 'numpy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fhess_lmbda = sympy.lambdify((x1, x2), fhess_sym, 'numpy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def func_XY_X_Y(f):\n",
    "    \"\"\"\n",
    "    Wrapper for f(X) -> f(X[0], X[1])\n",
    "    \"\"\"\n",
    "    return lambda X: np.array(f(X[0], X[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f = func_XY_X_Y(f_lmbda)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fprime = func_XY_X_Y(fprime_lmbda)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fhess = func_XY_X_Y(fhess_lmbda)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X_opt = optimize.fmin_ncg(f, (0, 0), fprime=fprime, fhess=fhess)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X_opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 4))\n",
    "x_ = y_ = np.linspace(-1, 4, 100)\n",
    "X, Y = np.meshgrid(x_, y_)\n",
    "c = ax.contour(X, Y, f_lmbda(X, Y), 50)\n",
    "ax.plot(X_opt[0], X_opt[1], 'r*', markersize=15)\n",
    "ax.set_xlabel(r\"$x_1$\", fontsize=18)\n",
    "ax.set_ylabel(r\"$x_2$\", fontsize=18)\n",
    "plt.colorbar(c, ax=ax)\n",
    "fig.tight_layout()\n",
    "fig.savefig('ch6-examaple-two-dim.pdf');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Brute force search for initial point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def f(X):\n",
    "    x, y = X\n",
    "    return (4 * np.sin(np.pi * x) + 6 * np.sin(np.pi * y)) + (x - 1)**2 + (y - 1)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_start = optimize.brute(f, (slice(-3, 5, 0.5), slice(-3, 5, 0.5)), finish=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_start"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f(x_start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_opt = optimize.fmin_bfgs(f, x_start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f(x_opt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def func_X_Y_to_XY(f, X, Y):\n",
    "    s = np.shape(X)\n",
    "    return f(np.vstack([X.ravel(), Y.ravel()])).reshape(*s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 4))\n",
    "x_ = y_ = np.linspace(-3, 5, 100)\n",
    "X, Y = np.meshgrid(x_, y_)\n",
    "c = ax.contour(X, Y, func_X_Y_to_XY(f, X, Y), 25)\n",
    "ax.plot(x_opt[0], x_opt[1], 'r*', markersize=15)\n",
    "ax.set_xlabel(r\"$x_1$\", fontsize=18)\n",
    "ax.set_ylabel(r\"$x_2$\", fontsize=18)\n",
    "plt.colorbar(c, ax=ax)\n",
    "fig.tight_layout()\n",
    "fig.savefig('ch6-example-2d-many-minima.pdf');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nonlinear least square"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def f(x, beta0, beta1, beta2):\n",
    "    return beta0 + beta1 * np.exp(-beta2 * x**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "beta = (0.25, 0.75, 0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xdata = np.linspace(0, 5, 50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y = f(xdata, *beta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ydata = y + 0.05 * np.random.randn(len(xdata))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def g(beta):\n",
    "    return ydata - f(xdata, *beta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "beta_start = (1, 1, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "beta_opt, beta_cov = optimize.leastsq(g, beta_start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "beta_opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.scatter(xdata, ydata, label=\"samples\")\n",
    "ax.plot(xdata, y, 'r', lw=2, label=\"true model\")\n",
    "ax.plot(xdata, f(xdata, *beta_opt), 'b', lw=2, label=\"fitted model\")\n",
    "ax.set_xlim(0, 5)\n",
    "ax.set_xlabel(r\"$x$\", fontsize=18)\n",
    "ax.set_ylabel(r\"$f(x, \\beta)$\", fontsize=18)\n",
    "ax.legend()\n",
    "fig.tight_layout()\n",
    "fig.savefig('ch6-nonlinear-least-square.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "beta_opt, beta_cov = optimize.curve_fit(f, xdata, ydata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "beta_opt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constrained optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def f(X):\n",
    "    x, y = X\n",
    "    return (x-1)**2 + (y-1)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_opt = optimize.minimize(f, [0, 0], method='BFGS').x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bnd_x1, bnd_x2 = (2, 3), (0, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x_cons_opt = optimize.minimize(f, [0, 0], method='L-BFGS-B', bounds=[bnd_x1, bnd_x2]).x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 4))\n",
    "x_ = y_ = np.linspace(-1, 3, 100)\n",
    "X, Y = np.meshgrid(x_, y_)\n",
    "c = ax.contour(X, Y, func_X_Y_to_XY(f, X, Y), 50)\n",
    "ax.plot(x_opt[0], x_opt[1], 'b*', markersize=15)\n",
    "ax.plot(x_cons_opt[0], x_cons_opt[1], 'r*', markersize=15)\n",
    "bound_rect = plt.Rectangle((bnd_x1[0], bnd_x2[0]), \n",
    "                           bnd_x1[1] - bnd_x1[0], bnd_x2[1] - bnd_x2[0],\n",
    "                           facecolor=\"grey\")\n",
    "ax.add_patch(bound_rect)\n",
    "ax.set_xlabel(r\"$x_1$\", fontsize=18)\n",
    "ax.set_ylabel(r\"$x_2$\", fontsize=18)\n",
    "plt.colorbar(c, ax=ax)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig('ch6-example-constraint-bound.pdf');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lagrange multiplier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = x1, x2, x3, l = sympy.symbols(\"x_1, x_2, x_3, lambda\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f = x1 * x2 * x3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "g = 2 * (x1 * x2 + x2 * x3 + x3 * x1) - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "L = f + l * g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "grad_L = [sympy.diff(L, x_) for x_ in x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sols = sympy.solve(grad_L)\n",
    "sols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "g.subs(sols[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f.subs(sols[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def f(X):\n",
    "    return -X[0] * X[1] * X[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def g(X):\n",
    "    return 2 * (X[0]*X[1] + X[1] * X[2] + X[2] * X[0]) - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "constraints = [dict(type='eq', fun=g)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "result = optimize.minimize(f, [0.5, 1, 1.5], method='SLSQP', constraints=constraints)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "result.x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inequality constraints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def f(X):\n",
    "    return (X[0] - 1)**2 + (X[1] - 1)**2\n",
    "\n",
    "def g(X):\n",
    "    return X[1] - 1.75 - (X[0] - 0.75)**4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%time x_opt = optimize.minimize(f, (0, 0), method='BFGS').x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "constraints = [dict(type='ineq', fun=g)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%time x_cons_opt = optimize.minimize(f, (0, 0), method='SLSQP', constraints=constraints).x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%time x_cons_opt = optimize.minimize(f, (0, 0), method='COBYLA', constraints=constraints).x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 4))\n",
    "x_ = y_ = np.linspace(-1, 3, 100)\n",
    "X, Y = np.meshgrid(x_, y_)\n",
    "c = ax.contour(X, Y, func_X_Y_to_XY(f, X, Y), 50)\n",
    "ax.plot(x_opt[0], x_opt[1], 'b*', markersize=15)\n",
    "\n",
    "ax.plot(x_, 1.75 + (x_-0.75)**4, 'k-', markersize=15)\n",
    "ax.fill_between(x_, 1.75 + (x_-0.75)**4, 3, color=\"grey\")\n",
    "ax.plot(x_cons_opt[0], x_cons_opt[1], 'r*', markersize=15)\n",
    "\n",
    "ax.set_ylim(-1, 3)\n",
    "ax.set_xlabel(r\"$x_0$\", fontsize=18)\n",
    "ax.set_ylabel(r\"$x_1$\", fontsize=18)\n",
    "plt.colorbar(c, ax=ax)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig('ch6-example-constraint-inequality.pdf');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear programming"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "c = np.array([-1.0, 2.0, -3.0])\n",
    "\n",
    "A = np.array([[ 1.0, 1.0, 0.0],\n",
    "              [-1.0, 3.0, 0.0],\n",
    "              [ 0.0, -1.0, 1.0]])\n",
    "\n",
    "b = np.array([1.0, 2.0, 3.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A_ = cvxopt.matrix(A)\n",
    "b_ = cvxopt.matrix(b)\n",
    "c_ = cvxopt.matrix(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sol = cvxopt.solvers.lp(c_, A_, b_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = np.array(sol['x'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sol['primal objective']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quandratic problem with cvxopt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quadratic problem formulation:\n",
    "    \n",
    "$\\min \\frac{1}{2}x^TPx + q^T x$\n",
    "\n",
    "$G x \\leq h$ \n",
    "\n",
    "For example, let's solve the problem\n",
    "\n",
    "min $f(x_1, x_2) = (x_1 - 1)^2 + (x_2 - 1)^2 =$\n",
    "\n",
    "$x_1^2 -2x_1 + 1 + x_2^2 - 2x_2 + 1 = $\n",
    "\n",
    "$x_1^2 + x_2^2 - 2x_1 - 2x_2 + 2 =$\n",
    "\n",
    "$= \\frac{1}{2} x^T P x - q^T x + 2$\n",
    "\n",
    "and\n",
    "\n",
    "$\\frac{3}{4} x_1 + x_2 \\geq 3$, $x_1 \\geq 0$\n",
    "\n",
    "where\n",
    "\n",
    "$P = 2 [[1, 0], [0, 1]]$ and $q = [-2, -2]$\n",
    "\n",
    "and\n",
    "\n",
    "$G = [[-3/4, -1], [-1, 0]]$ and $h = [-3, 0]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cvxopt import matrix, solvers\n",
    "\n",
    "P = 2 * np.array([[1.0, 0.0],\n",
    "                  [0.0, 1.0]])\n",
    "q = np.array([-2.0, -2.0])\n",
    "\n",
    "G = np.array([[-0.75, -1.0],\n",
    "              [-1.0,  0.0]])\n",
    "\n",
    "h = np.array([-3.0, 0.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_P = cvxopt.matrix(P)\n",
    "_q = cvxopt.matrix(q)\n",
    "_G = cvxopt.matrix(G)\n",
    "_h = cvxopt.matrix(h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time sol = solvers.qp(_P, _q, _G, _h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = sol['x']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol['primal objective'] + 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 4))\n",
    "x_ = y_ = np.linspace(-1, 3, 100)\n",
    "X, Y = np.meshgrid(x_, y_)\n",
    "c = ax.contour(X, Y, func_X_Y_to_XY(f, X, Y), 50)\n",
    "\n",
    "\n",
    "\n",
    "y_ = (h[0] - G[0, 0] * x_)/G[0, 1]\n",
    "angle = -np.arctan((y_[mask][0] - y_[mask][-1]) / (x_[mask][-1]- x_[mask][0])) * 180 / np.pi\n",
    "mask = y_ < 3\n",
    "ax.plot(x_, y_, 'k')\n",
    "ax.add_patch(plt.Rectangle((0, 3), 4, 3, angle=angle, facecolor=\"grey\"))\n",
    "\n",
    "ax.plot(x_opt[0], x_opt[1], 'b*', markersize=15)\n",
    "ax.plot(x[0], x[1], 'r*', markersize=15)\n",
    "\n",
    "ax.set_ylim(-1, 3)\n",
    "ax.set_xlim(-1, 3)\n",
    "ax.set_xlabel(r\"$x_1$\", fontsize=18)\n",
    "ax.set_ylabel(r\"$x_2$\", fontsize=18)\n",
    "plt.colorbar(c, ax=ax)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig('ch6-example-quadratic-problem-constraint-inequality.pdf');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Arbitrary function callback API"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cvxopt import modeling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1 = modeling.variable(1, \"x1\")\n",
    "x2 = modeling.variable(1, \"x2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1 * x2 + x2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(modeling)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(solvers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(solvers.cp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CVXPY \n",
    "\n",
    "cvxpy.org"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Versions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%reload_ext version_information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%version_information numpy, scipy, cvxopt, sympy, matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.6 (NPM ed2)",
   "language": "python",
   "name": "py3.6-npm-e2"
  },
  "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}