{
 "metadata": {
  "name": "",
  "signature": "sha256:44a320fd2431c789211a6ca4eb84e9c76e71fe2594baf5b13bfab0b49b75d2a7"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": [],
     "source": [
      "> This is one of the 100 recipes of the [IPython Cookbook](http://ipython-books.github.io/), the definitive guide to high-performance scientific computing and data science in Python.\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# 9.2. Minimizing a mathematical function"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "1. We import the libraries."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import scipy as sp\n",
      "import scipy.optimize as opt\n",
      "import matplotlib.pyplot as plt\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "2. First, let's define a simple mathematical function (the opposite of the **cardinal sine**). This function has many local minima but a single global minimum. (http://en.wikipedia.org/wiki/Sinc_function)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f = lambda _: 1-np.sin(_)/_"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "3. Let's plot this function on the interval $[-20, 20]$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.linspace(-20., 20., 1000)\n",
      "y = f(x)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.figure(figsize=(5,5));\n",
      "plt.plot(x, y);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "4. The `scipy.optimize` module comes with many function minimization routines. The `minimize` function offers a unified interface to many algorithms. The **Broyden\u2013Fletcher\u2013Goldfarb\u2013Shanno (BFGS) algorithm** (default in `minimize`) gives good results in general. The `minimize` function requires an initial point as argument. For scalar univariate functions, we can also use `minimize_scalar`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x0 = 3\n",
      "xmin = opt.minimize(f, x0).x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Starting from $x_0=3$, the algorithm was able to find the actual global minimum, as shown on the following figure."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.figure(figsize=(5,5));\n",
      "plt.plot(x, y);\n",
      "plt.scatter(x0, f(x0), marker='o', s=300);\n",
      "plt.scatter(xmin, f(xmin), marker='v', s=300);\n",
      "plt.xlim(-20, 20);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "5. Now, if we start from an initial point that is further away from the actual global minimum, the algorithm converges towards a *local* minimum only."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x0 = 10\n",
      "xmin = opt.minimize(f, x0).x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.figure(figsize=(5,5));\n",
      "plt.plot(x, y);\n",
      "plt.scatter(x0, f(x0), marker='o', s=300);\n",
      "plt.scatter(xmin, f(xmin), marker='v', s=300);\n",
      "plt.xlim(-20, 20);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "6. Like most function minimization algorithms, the BFGS algorithm is efficient at finding *local* minima, but not necessarily *global* minima, especially on complicated or noisy objective functions. A general strategy to overcome this problem consists in combining such algorithms with an exploratory grid search on the initial points. Another option is to use a different class of algorithms based on heuristics and stochastic methods. A popular example is the **simulated annealing method**."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "xmin = opt.minimize(f, x0, method='Anneal').x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.figure(figsize=(5,5));\n",
      "plt.plot(x, y);\n",
      "plt.scatter(x0, f(x0), marker='o', s=300);\n",
      "plt.scatter(xmin, f(xmin), marker='v', s=300);\n",
      "plt.xlim(-20, 20);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This time, the algorithm was able to find the global minimum."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "7. Now, let's define a new function, in two dimensions this time. This function is called the **L\u00e9vi function**. It is defined by\n",
      "\n",
      "$$f(x,y) = \\sin^{2}\\left(3\\pi x\\right)+\\left(x-1\\right)^{2}\\left(1+\\sin^{2}\\left(3\\pi y\\right)\\right)+\\left(y-1\\right)^{2}\\left(1+\\sin^{2}\\left(2\\pi y\\right)\\right)$$\n",
      "\n",
      "This function is very irregular and may be difficult to minimize in general. It is one of the many **test functions for optimization** that researchers have developed to study and benchmark optimization algorithms. (http://en.wikipedia.org/wiki/Test_functions_for_optimization)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def g(X):\n",
      "    # X is a 2*N matrix, each column contains\n",
      "    # x and y coordinates.\n",
      "    x, y = X\n",
      "    return np.sin(3*np.pi*x)**2+(x-1)**2*(1+np.sin(3*np.pi*y)**2)+(y-1)**2*(1+np.sin(2*np.pi*y)**2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "8. Let's display this function with `imshow`, on the square $[-10,10]^2$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 200\n",
      "k = 10\n",
      "X, Y = np.mgrid[-k:k:n*1j,-k:k:n*1j]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Z = g(np.vstack((X.ravel(), Y.ravel()))).reshape(n,n)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.figure(figsize=(5, 5));\n",
      "# We use a logarithmic scale for the color here.\n",
      "plt.imshow(np.log(Z), cmap=plt.cm.hot_r);\n",
      "plt.xticks([]); plt.yticks([]);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "9. The BFGS algorithm also works in multiple dimensions."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x0, y0 = opt.minimize(g, (8, 3)).x"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.figure(figsize=(5, 5));\n",
      "plt.imshow(np.log(Z), cmap=plt.cm.hot_r,\n",
      "           extent=(-k, k, -k, k), origin=0);\n",
      "plt.scatter(x0, y0, s=100);\n",
      "plt.xticks([]); plt.yticks([]);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "> You'll find all the explanations, figures, references, and much more in the book (to be released later this summer).\n",
      "\n",
      "> [IPython Cookbook](http://ipython-books.github.io/), by [Cyrille Rossant](http://cyrille.rossant.net), Packt Publishing, 2014 (500 pages)."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}