{
"metadata": {
"name": "",
"signature": "sha256:0c4193e9147ea2979632ab81eef54c9508b9b59b07a1bf09c81bd40e58019617"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "code",
"collapsed": false,
"input": [
"from __future__ import division\n",
"from numpy import *\n",
"from matplotlib import pyplot as plt\n",
"from matplotlib import animation\n",
"from mpl_toolkits.mplot3d import axes3d\n",
"%pylab inline\n",
"\n",
"#from http://jakevdp.github.io/blog/2013/05/12/embedding-matplotlib-animations/\n",
"from tempfile import NamedTemporaryFile\n",
"\n",
"VIDEO_TAG = \"\"\"\"\"\"\n",
"\n",
"def anim_to_html(anim):\n",
" if not hasattr(anim, '_encoded_video'):\n",
" with NamedTemporaryFile(suffix='.mp4') as f:\n",
" anim.save(f.name, fps=20,extra_args=['-vcodec', 'libx264', '-pix_fmt', 'yuv420p'])\n",
" video = open(f.name, \"rb\").read()\n",
" anim._encoded_video = video.encode(\"base64\")\n",
" return VIDEO_TAG.format(anim._encoded_video)\n",
"\n",
"from IPython.display import HTML\n",
"def display_animation(anim):\n",
" plt.close(anim._fig)\n",
" return HTML(anim_to_html(anim))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Populating the interactive namespace from numpy and matplotlib\n"
]
}
],
"prompt_number": 219
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def rosenbrock(v):\n",
" return (1.-v[...,0])**2 + 100.*(v[...,0]-v[...,1]**2)**2\n",
"def d_rosenbrock(v):\n",
" if len(shape(v))==1:\n",
" return array([\n",
" 2.*(-100*(v[1]-v[0]**2)*2*v[0] + v[0]-1),\n",
" 200.*(v[1]-v[0]**2)\n",
" ])\n",
" else:\n",
" return array([\n",
" 2.*(-100*(v[...,1]-v[...,0]**2)*2*v[...,0] + v[...,0]-1),\n",
" 200.*(v[...,1]-v[...,0]**2)\n",
" ])\n",
"def plot_ros_simplices(simplices,n_frames = 20, figsize=(8, 6), dpi=80):\n",
" n_simplices = shape(simplices)[0]\n",
" fig = plt.figure(figsize=figsize, dpi=dpi)\n",
" ax = fig.add_subplot(111, projection='3d')\n",
" ax.axis('off')\n",
" ax.set_xlim((-1, 1))\n",
" ax.set_ylim((-1, 1))\n",
" N = 100;\n",
" x = linspace(-1,1.5,N); y = linspace(-1.,1.5,N);\n",
" X,Y = meshgrid(x,y)\n",
" Z = rosenbrock(dstack((X,Y)))\n",
" simplex_plots = []\n",
" for s in simplices: simplex_plots += ax.plot([], [], [], marker='.',linestyle='-', c='r')\n",
" \n",
" def init():\n",
" ax.plot_wireframe(X, Y, Z, rstride=10, cstride=10)\n",
" for s in simplex_plots:\n",
" s.set_data([], [])\n",
" s.set_3d_properties([])\n",
" def animate(i):\n",
" ax.azim = 360/n_frames*i\n",
" ax.elev = 30\n",
" si = ceil(i*n_simplices/n_frames)\n",
" loop = array([0,1,2,0])\n",
" for s,sp in zip(simplices[:si],simplex_plots):\n",
" sp.set_data(s[loop,0],s[loop,1]); sp.set_3d_properties(rosenbrock(s[loop]))\n",
" return animation.FuncAnimation(fig, animate, init_func=init,\n",
" frames=n_frames, interval=20, blit=True)\n",
"display_animation(plot_ros_simplices(array([[[-1,-1],[-1,-1],[-1,-1]]])))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
""
],
"metadata": {},
"output_type": "pyout",
"prompt_number": 220,
"text": [
""
]
}
],
"prompt_number": 220
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"start = array([-1,-1]);\n",
"simplex = vstack((start+.1*eye(2),array(start)))\n",
"def nelder_mead(simplex,f,d):\n",
" simplices = simplex[None,...] #keep track\n",
" while(True):\n",
" simplex = simplices[-1] #our working simplex\n",
" simplex = simplex[argsort(-f(simplex))] #sort the vertices, descending by f value\n",
" xmean = sum(simplex[1:],axis=0)/d\n",
" if f(2*xmean-simplex[0]) < f(simplex[-1]): #reflection is best\n",
" if f(3*xmean-2*simplex[0]) < f(2*xmean-simplex[0]):\n",
" simplex[0] = 3*xmean-2*simplex[0] #grow works, keep it\n",
" else:\n",
" simplex[0] = 2*xmean-simplex[0] #grow doesn't work, keep original reflection\n",
" elif f(2*xmean-simplex[0]) < f(simplex[1]): #reflection helped a bit, keep and iterate\n",
" simplex[0] = 2*xmean-simplex[0]\n",
" else: #reflection doesn't help at all\n",
" if f(1.5*xmean+.5*simplex[0]) < f(simplex[1]): #try reflect and shrink\n",
" simplex[0] = 1.5*xmean+.5*simplex[0]\n",
" elif f(.5*xmean+.5*simplex[0]) < f(simplex[1]): #try just shrinking\n",
" simplex[0] = .5*xmean+.5*simplex[0]\n",
" else: simplex[:-1] = .5*(simplex[:-1]+simplex[-1]) #nothing is working, shrink all towards best\n",
" simplices = vstack((simplices,simplex[None,...]))\n",
" if abs( f(sum(simplices[-1],axis=0)/(d+1)) - f(sum(simplices[-2],axis=0)/(d+1)) ) < 1e-6:\n",
" break\n",
" return simplices\n",
"anim = plot_ros_simplices(nelder_mead(simplex,rosenbrock,2),n_frames=50,dpi=150,figsize=(10,8))\n",
"display_animation(anim)\n",
"#anim.save('nelder-mead.gif', writer='imagemagick', fps=10)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
""
],
"metadata": {},
"output_type": "pyout",
"prompt_number": 221,
"text": [
""
]
}
],
"prompt_number": 221
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"phi = .5*(1+sqrt(5))\n",
"def find_decrease(x0,d,f,alpha=1):\n",
" #from a starting point, direction, and initial step length, return a step length that gives descent\n",
" a = alpha\n",
" while(True):\n",
" x1 = x0+a*d\n",
" if f(x1) < f(x0): break\n",
" a /= phi\n",
" if a<1e-4:\n",
" a = -alpha #if we shrink all the way to the start, switch direction\n",
" return a\n",
"def bracket(x0,d,f,alpha=1):\n",
" #find a bracketing triple from starting point, direction, and intitial step length\n",
" alpha = find_decrease(x0,d,f,alpha=alpha)\n",
" while(True):\n",
" x1,x2 = x0+alpha*d,x0+alpha*(1+1/phi)*d\n",
" if (f(x1) < min(f(x0),f(x2))):\n",
" break\n",
" alpha *= phi\n",
" return x0,x1,x2 \n",
"#test\n",
"plt.figure()\n",
"b = array(bracket(0,1,cos))\n",
"t = linspace(0,b[2],100)\n",
"plt.plot(t,cos(t)); plt.plot(b,cos(b),linestyle='',marker='.',color='r')"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 431,
"text": [
"[]"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEACAYAAABbMHZzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtclVW+x/EPitWoBeYcMZHxAhiQCphKNx0w0cRLlpaa\nqaPOZJajdrWmOSedsjSnmsqmo11Mu6jljGlqdHkV4iU0FbWySU09ASp57ZhaIDznj3VyUkFwb/Ze\nz977+3699suQx72/LfG31/4961lPmOM4DiIiElJq2Q4gIiL+p+IvIhKCVPxFREKQir+ISAhS8RcR\nCUEq/iIiIcjr4j9ixAiioqJo06ZNpceMHTuW+Ph4kpOTyc/P9/YlRUTES14X/+HDh5OdnV3p95ct\nW8b27dvZtm0bM2fOZPTo0d6+pIiIeMnr4t+pUycaNGhQ6fcXL17MsGHDAEhLS+Pw4cMUFxd7+7Ii\nIuIFn/f8i4qKiImJOfl106ZNKSws9PXLiojIWfjlhO/pO0iEhYX542VFRKQS4b5+gejoaAoKCk5+\nXVhYSHR09BnHxYWF8Y2vw4iIBJnY2Fi2b99+zn/O5zP/Pn36MGfOHADy8vKIjIwkKirqjOO+AZz2\n7XEOHcJxnAofxcUOixY5PPCAw1VXOVx4oUPPng7PP++we3fFf8abx8MPP1zjz6lMoZ1LmZSpph/f\nfOPZtNnrmf+gQYNYvnw5+/fvJyYmhkmTJlFaWgrAqFGjyMrKYtmyZcTFxVGvXj1mzZpV+ZN9+CFE\nRlb67UaNoE8f8wA4eND8kSVL4KGH4PLLYdAguOkmuOgib//PRESCl9fFf+7cuVUeM3369Oo92VkK\nf0UuvhgGDDCP48dh6VJ44w24917ze6NHQ3LyOT2liEhICJorfH/1K+jfHxYuhC+/hCZNoGdPSE+H\n998Hx4O7FqSnp9d0TK8pU/W5MZcyVY8y+V6Y43hSFmteWFgYNR2ltBTmz4fHH4cLLjCtoRtuAC02\nEpFg4WntDOri/7Pycnj3XZg4EerUgSlToEsXn7yUiIhfqfhXQ3m5+STw5z9DXBz87W+QmOjTlxQR\n8SlPa2fQ9Pyro1Ytsxroq6+gRw/o3Bnuuw+OHLGdTETEv0Kq+P/svPNg/Hj44gvYv9/M/hcutJ1K\nRMR/QqrtU5mVK2HkSGjbFqZPhwquQRMRcSW1fbxwzTWwcaM5D9C2rTkvICISzDTzP81nn8HgwXDl\nlfDcc7pSWETcTTP/GtKhA+Tnw/nnQ0oKrFljO5GISM1T8a9AvXowcyY8+ST07g3PPuvZFcIiIm6l\ntk8Vduww20bExcFLL6kNJCLuoraPj7RsCatXQ4MGkJYG27bZTiQi4j0V/2q44AKYMcNcG3DNNfDB\nB7YTiYh4R22fc5Sba7aLnjABxo3TJnEiYpf29vGjXbugVy/IyDD7A9WubTuRiIQqFX8/O3wY+vWD\nCy+EN9+EunVtJxKRUKQTvn4WGQnvvQcREeYTwP79thOJiFSfir8XzjsPXn0VunaFTp2goMB2IhGR\n6vH6Hr6hLiwMJk+Ghg3NSqD334eEBNupRETOTsW/htx9t3kDyMgwN5Jv1852IhGRyqn416Bhw8wV\nwD16mNtGduxoO5GISMVU/GvYDTeY+wT36gXvvANXXWU7kYjImXTC1wd69YLXXoO+fc1FYSIibqPi\n7yPdu5v1//37m72BRETcRMXfh7p2/fcngLVrbacREfk3FX8f694dZs0y9wVYv952GhERQ8XfD3r2\nNLuC9uoFW7bYTiMiotU+ftO3Lxw5Yj4J5OZCixa2E4lIKFPx96MhQ+D77yEzE1asgEsusZ1IREKV\nir+fjRljdgT9+RNAZKTtRCISirSlswWOY24Es3kzZGebO4WJiHhC+/kHmPJyGDTI/Dpvnm4IIyKe\n0X7+AaZWLZgzBw4cMJ8CQuh9T0RcQMXfovPPh4ULTe//qadspxGRUKITvpZFRJgtoK+6Cpo1M9tB\niIj4moq/C8TEmC2gu3WD6Gi48krbiUQk2Knt4xIpKeaWkDfeCDt22E4jIsFOxd9FsrLgz382+wB9\n/73tNCISzLTU04XuvNPM/t99F8LVmBORs9BSzyDyzDNw4gTce6/tJCISrFT8XSg8HN5+21z9++KL\nttOISDBS28fFtm6FTp3gn/+Eq6+2nUZE3EhtnyDUqhXMng033wwFBbbTiEgwUfF3ueuug/Hj4YYb\n4Phx22lEJFio7RMAHAduvdVs/jZ7NoSF2U4kIm6htk8QCwszJ343b4bnn7edRkSCgdfFPzs7m4SE\nBOLj45k6deoZ38/JySEiIoLU1FRSU1N59NFHvX3JkFS3rjnx+8gjsHKl7TQiEui8uoSorKyMMWPG\n8NFHHxEdHU2HDh3o06cPiYmJpxz329/+lsWLF3sVVKBlS7MFxIAB8Nln0KSJ7UQiEqi8mvmvXbuW\nuLg4mjdvTp06dRg4cCCLFi064zj18mtOjx4wejTcdBOUltpOIyKByqviX1RURExMzMmvmzZtSlFR\n0SnHhIWFsXr1apKTk8nKymLLli3evKQAf/oTNGgAEybYTiIigcqrtk9YNZadtGvXjoKCAurWrct7\n771H37592bp1a4XHTpw48eR/p6enk56e7k28oPXzXcAuv9zcB0D3ABAJHTk5OeTk5Hj9PF4t9czL\ny2PixIlkZ2cD8Pjjj1OrVi0mnGVK2qJFC9avX8/FF198ahAt9Txn69aZNtCqVeaCMBEJPVaWerZv\n355t27axa9cuSkpKmD9/Pn369DnlmOLi4pPB1q5di+M4ZxR+8Uz79vDoo9CvHxw7ZjuNiAQSr9o+\n4eHhTJ8+ne7du1NWVsbIkSNJTExkxowZAIwaNYoFCxbwwgsvEB4eTt26dZk3b16NBBfjtttgxQoY\nOxZeesl2GhEJFLrCNwj88IP5FPDQQzBkiO00IuJPntZOFf8g8fnn0KUL5ObCaZdZiEgQ0/YOIa5N\nG5gyxaz/V/9fRKqimX8QcRzT9qlbF2bOtJ1GRPxBM38hLAxeeAE++QTmz7edRkTcTDP/ILRhg7kP\nQF6e2Q9IRIKXZv5yUrt2ZuXPwIFQUmI7jYi4kWb+Qcpx4PrrzcqfCnbaFpEgoaWecob9+yElxWwD\n3bWr7TQi4gtq+8gZfv1rc9vHYcNg3z7baUTETTTzDwETJsCXX8K77+r+vyLBRjN/qdQjj8B33+n+\nvyLyb5r5h4ht28ze/zk5cNllttOISE3RzF/OKj7ebP9wyy3w00+204iIbZr5hxDHMXf9at4cnnzS\ndhoRqQla6inVcuCAWf75yiuQmWk7jYh4S20fqZaGDc26/+HD4eBB22lExBbN/EPUXXdBUZHZAE7L\nP0UCl2b+ck4efxy2bIHXX7edRERs0Mw/hG3aZLZ9WLcOmjWznUZEPKGZv5yz5GS47z4YOhTKymyn\nERF/UvEPcffcY5aAPv207SQi4k9q+wg7d0LHjvDxx+ZewCISONT2EY+1aGH2/B8yRFf/ioQKzfwF\nMK2fvn0hKcmsBBKRwKArfMVr331nTgL/4x9mEzgRcT+1fcRrjRqZbZ+HDYOjR22nERFf0sxfzjBk\nCFx0kfb/FwkEavtIjTl82Kz60eZvIu6nto/UmMhIePllGDnSvBGISPDRzF8qdccdcOyY2QVURNxJ\nM3+pcU88AStWwOLFtpOISE3TzF/OKjcXBg6Ezz839wIQEXfRCV/xmbvvht27Yd4820lE5HRq+4jP\nTJ5stn9+6y3bSUSkpmjmL9WyZg1cf715E4iKsp1GRH6mto/43IMPwtdfm+0fdOtHEXdQ20d8buJE\n2LoV5s61nUREvKWZv5yT9eshKws2boRLLrGdRkQ08xe/uPxyuP12uO02sw20iAQmFX85Zw89BAUF\n8NprtpOIiKfU9hGPbNwI3bqZX5s0sZ1GJHSp7SN+lZICd96p9o9IoFLxF4/96U9QVASzZ9tOIiLn\nSm0f8cqmTWbP//x8iI62nUYk9KjtI1YkJ8OYMfCHP6j9IxJIVPzFaw8+CHv3at9/kUCito/UiM2b\n4dprYcMGiImxnUYkdFhr+2RnZ5OQkEB8fDxTp06t8JixY8cSHx9PcnIy+fn53r6kuFDbtjB2rNo/\nIoHCq+JfVlbGmDFjyM7OZsuWLcydO5evvvrqlGOWLVvG9u3b2bZtGzNnzmT06NFeBRb3euAB2LcP\nZs2ynUREquJV8V+7di1xcXE0b96cOnXqMHDgQBYtWnTKMYsXL2bYsGEApKWlcfjwYYqLi715WXGp\nOnVM4Z8wwVwBLCLu5VXxLyoqIuYXDd6mTZtSVFRU5TGFhYXevKy4mNo/IoEh3Js/HFbNTd1PPxlR\n2Z+bOHHiyf9OT08nPT3d02hi0QMPwBVXmE8BI0bYTiMSXHJycsjJyfH6ebwq/tHR0RT84vN9QUEB\nTZs2PesxhYWFRFdyNdAvi78Erjp1zLLPLl3MBWBa/SNSc06fGE+aNMmj5/Gq7dO+fXu2bdvGrl27\nKCkpYf78+fTp0+eUY/r06cOcOXMAyMvLIzIykijdBzDotWkD48ap/SPiVl4V//DwcKZPn0737t1J\nSkpiwIABJCYmMmPGDGbMmAFAVlYWLVu2JC4ujlGjRvH3v/+9RoKL+02YYFb/vPKK7SQicjpd5CU+\n9fnnpv2ji79EfEN7+4grqf0j4k4q/uJzav+IuI/aPuIXav+I+IbaPuJqbdrA+PFq/4i4hYq/+M2E\nCbB/P7z8su0kIqK2j/jVF19ARgasXw+/+Y3tNCKBT20fCQitW8Pdd8PIkWr/iNik4i9+d9998P33\nMHOm7SQioUttH7Fiyxbo3Bk++wxatLCdRiRwqe0jASUpCe6/37R/ysttpxEJPSr+Ys0998Dx46Dt\nnkT8T20fserrr+HqqyEvD+LibKcRCTxq+0hAuvRSeOghGD4cyspspxEJHSr+Yt24cRAWBs88YzuJ\nSOhQ20dc4ZtvIC0NVq6EhATbaUQCh9o+EtBiY+Evf4Fhw+DECdtpRIKfir+4xu23w0UXwRNP2E4i\nEvzU9hFXKSiAyy+HDz+E5GTbaUTcT20fCQoxMTBtGgwdCiUlttOIBC8Vf3GdoUOheXOYNMl2EpHg\npbaPuNLevZCSAu+8A1dcYTuNiHup7SNBpXFjeP558yng6FHbaUSCj2b+4mpDh5oVQNOn204i4k6e\n1k4Vf3G1w4ehbVt46SXo1s12GhH3UdtHglJkJLzyitn6+dAh22lEgodm/hIQxo2DffvgzTdtJxFx\nF838JahNmQL5+TB3ru0kIsFBM38JGOvXQ48e5teYGNtpRNxBM38Jepdfbto/w4fr1o8i3lLxl4Ay\nYQIcO6a9/0W8pbaPBJxvvjFX/X78MbRpYzuNiF1q+0jIiI012z4PHgw//mg7jUhg0sxfApLjQP/+\n0KwZPPWU7TQi9ugKXwk5Bw6YPf9nzYLMTNtpROxQ20dCTsOG8OqrZvXP/v2204gEFs38JeDddx9s\n3Wq2fw4Ls51GxL8085eQNXkyFBbCjBm2k4gEDs38JSh8/TVcfTXk5kJSku00Iv6jmb+EtEsvNfv/\nDByo5Z8i1aGZvwQNx4EBA6BRI938RUKHZv4S8sLCYOZMWLrUnPwVkcpp5i9B59NPoW9fWLdOu39K\n8NPMX+T/XXkljB9vtn84ccJ2GhF3UvGXoDRhAlxwAUyaZDuJiDup7SNBq7gY2rWD2bOha1fbaUR8\nQ20fkdNERcGcOTB0KOzdazuNCBy99TbKO6dDVhYcPmw1i2b+EvQefhhWroQPPoDatW2nkVB1/Dh8\n8et0Ohxbbn7jppvgrbe8fl6/z/wPHjxIZmYmrVq1olu3bhyu5F2sefPmtG3bltTUVDp27Ojpy4l4\n7L/+y9z28S9/sZ1EQtkf/wh1IuqaL9q3N+uSLfK4+E+ZMoXMzEy2bt3Ktddey5QpUyo8LiwsjJyc\nHPLz81m7dq3HQUU8Vbs2zJ0LL71kZv8i/jZ7tvn0GbvmTTPj//BDiIy0msnjtk9CQgLLly8nKiqK\nvXv3kp6ezr/+9a8zjmvRogXr1q2jYcOGZw+ito/4WE6O2f5h3Tpo2tR2GgkVX3wBGRm+u+2o39s+\nxcXFREVFARAVFUVxcXGlwbp27Ur79u158cUXPX05Ea+lp8O4cXDzzVBaajuNhIIjR8xEf9o0991v\nOvxs38zMzGRvBcskJk+efMrXYWFhhFWykfqqVau45JJL2LdvH5mZmSQkJNCpU6cKj504ceLJ/05P\nTyc9Pb2K+CLnZsIEWLUK7r8fnn7adhoJZo4DI0fCNdfA735Xc8+bk5NDTk6O18/jVdsnJyeHxo0b\ns2fPHjIyMips+/zSpEmTqF+/Pvfcc8+ZQdT2ET85dMicb5s82bSBRHzh6afhjTdMr/+CC3z3On5v\n+/Tp04fZs2cDMHv2bPr27XvGMceOHePIkSMAHD16lA8++IA2bvvsIyGnQQP4xz/M6osvv7SdRoLR\nihUwdSosWODbwu8Nj2f+Bw8e5Oabb+bbb7+lefPmvPXWW0RGRrJ7927+8Ic/sHTpUnbs2MGNN94I\nwIkTJxg8eDAPPvhgxUE08xc/mz3bzP4/+wwiImynkWCxezd06AAvvwzXXef71/O0duoiLwlpd9wB\nRUWwcCHU0vXu4qWffjILC3r1goce8s9rqviLeKCkBLp0gcxMcyWwiDdGjYJ9+0y7x1+TCU9r51lX\n+4gEu/POM/9QO3SAlBS4/nrbiSRQzZxpev1r1gTGp0jN/EUw/2B79TI3gE9MtJ1GAs3q1eYGQitX\nQqtW/n1t7eop4oW0NHjiCTPzP3jQdhoJJN9+C/37mwUE/i783tDMX+QX7r4bNm+G996DOnVspxG3\nO3rUXMR1661QweVLfqETviI14MQJ6N0b4uLguedspxE3cxyzVUjduvDqq1DJJgc+p7aPSA0ID4d5\n8+Cjj2DGDNtpxM0mTYLCQvNzYqvwe0OrfUROExEBixdDp04QG6tbQMqZ3njD9Pjz8tx7BW9V1PYR\nqURurjmRl5MDSUm204hbrFoFN9xgtmhu3dp2GrV9RGpc587w17+aJaDffWc7jbjBjh1mQvDaa+4o\n/N5Q8Rc5i6FDYfBgswT0+HHbacSmAwegRw9zW9Du3W2n8Z7aPiJVcBwYMsQs61uwQDeBD0XHj5tz\nP9dcY3brdBMt9RTxoZISM+tLTDRLQANxdYd4pqwMBgwwW4G8/rr7tm5Qz1/Eh847D/75T3MSeNo0\n22nEXxzHXLx14ADMmuW+wu8NLfUUqaaICHPl71VXQePG5nyABLcpU8yqnuXL4fzzbaepWSr+Iucg\nOhqysyEjw9wRrHdv24nEV156CV580WzW1qCB7TQ1Tz1/EQ+sXQs9e5rbQXbubDuN1LSFC+HOO82M\nPz7edpqzU89fxI86doS5c82a7w0bbKeRmvTRR+amLEuWuL/we0PFX8RDXbvCf/83ZGXpRvDBYuVK\nuOUWc3K/XTvbaXxLPX8RL9x4o1kD3q0bfPJJYO3nLqdat878fb7xhlnPH+xU/EW8NHgw/Pij+SSw\nfDm0aGE7kZyrTZvMNh4vv2zu5xwKVPxFasDIkeYTQJcu5hNA8+a2E0l1bdoE110H06eH1uotFX+R\nGjJmjPk1Pd28AegTgPv9XPife86cvA8lKv4iNWjMGLP1Q0aGuTioZUvbiaQy+fnmZH0oFn5Q8Rep\ncXfeabYBSE+HDz6AhATbieR0n34KffvC3/8O/frZTmOHir+ID4webe7tmpEBy5ZBaqrtRPKzTz4x\n996dM8ds1heqVPxFfGTYMKhf3+z9/s47Zk8gsWvJEhgxAt5+23wyC2W6yEvEh/r1MzPM6683hUfs\nmTULfv97ePddFX7Q3j4ifpGXZ3rMkyebZaHiP45jduecOdNsynfppbYT1SxPa6faPiJ+cMUV5gKw\nHj2gqAj+8z91Qxh/OHEC7rrLjP2qVdCkie1E7qGZv4gf7d1rdgNt3drMRINtj3g3+d//hYEDzRvA\nW29BZKTtRL6hXT1FAkDjxuZuYEeOmO0g9u+3nSg4ffut2Z/nN7+BpUuDt/B7Q8VfxM/q1TM3gu/U\nCdLS4IsvbCcKLrm5ps32u9/BCy9AnTq2E7mT2j4iFr3+uulJT59ubhIunnMceP55eOQRmD3bbNsQ\nCjytnSr+IpZt3Gi2Eu7bF6ZO1UzVE8eOwR13mC0bFi4MrW011PMXCVApKWYv+a5v38bnDdM5npEF\nhw/bjhUwvvwSOnSAsjJYvTq0Cr83VPxFXODii6FHy620O7KcX+W8R2HWbbYjuZ7jmP3309Ph3nvN\nxXT16tlOFTi0zl/EJcLq1QXgh8T29CyayRWj4K9/hQsvtBzMhfbtg9tvh6+/Nmv4k5JsJwo8mvmL\nuMWbb8JNN1F/9Yfkbo7kxAlIToacHNvB3GXxYjMuLVuadpkKv2d0wlfExZYsgVGjzB5Bjz4KF11k\nO5E9330Hd99t+vqzZ5ulsqITviJBqVcv+Pxz+OEHM8NdsMD0ukPJz7391q3hkkvMeKjwe08zf5EA\nkZtr7hPQrBk89VRo3CTms8/MdRAlJWY7jJQU24ncRzN/kSDXubNZx96li5n53nmnOfEZjAoLYehQ\nsxX2iBHmzlsq/DVLxV8kgJx3nlnW+NVXULs2JCbCxInBc1nA3r1mpt+2LcTEmNU8I0aY/1epWSr+\nIgHo17+GZ5+FNWvMJmZxcfDww3DggO1kniksNG9qSUmmx//ll+beB1rm6jsq/iIBLDYWXnnFvAkU\nFpo3gdtvN58MAsGGDXDrrWamf+IEbN4Mf/ubObErvuVx8X/77be57LLLqF27Nhs2bKj0uOzsbBIS\nEoiPj2fq1KmevpyInEVsrFkR869/mW2jMzLg2mvNxnHHjtlOd6rvv4cZM6BjR9PTb9sWduwwRb9p\nU9vpQofHxb9NmzYsXLiQzp07V3pMWVkZY8aMITs7my1btjB37ly+CpQpCZDjwqtrlKn63JjL15mi\nosw5gP/5H/MJ4M03IToahg83F0cdP+7/TGDuXzB/Ptx8s1mt9OGHJueuXXD//Wfutx+Kf3f+5nHx\nT0hIoFWrVmc9Zu3atcTFxdG8eXPq1KnDwIEDWbRokacv6Xdu/MtWpupzYy5/ZTr/fLjpJli2zPTP\nU1LMzLpxY7jhBnjuObNevrzcN5nKykxLZ9o0s7VydLS5MKtbN9i+3VyvkJVV+YncUP678xef7u1T\nVFRETEzMya+bNm3KmjVrfPmSInKaJk1g3DjzOHAA3n/fbBnx3HNw8KC5avjoUfMGER9vZuaNGlXv\nHsOOY1bo7NhhivqmTbB+vdmmOjraLEsdNQrmzdPdtNzmrMU/MzOTvXv3nvH7jz32GL17967yycN0\nh2oRV2nYEG65xTwA9uwxq2waNIBFi0wR37XLnCf4j/+AiAhTtH/1K/MpwXGgtBQOHTJvHAcOmDeP\nli3NIznZ3Jy+XTuzU6m4mOOl9PR0Z/369RV+79NPP3W6d+9+8uvHHnvMmTJlSoXHxsbGOoAeeuih\nhx7n8IiNjfWodtdI28ep5NLi9u3bs23bNnbt2kWTJk2YP38+c+fOrfDY7du310QUERGpBo9P+C5c\nuJCYmBjy8vLo2bMnPXr0AGD37t307NkTgPDwcKZPn0737t1JSkpiwIABJCYm1kxyERHxmGs2dhMR\nEf/x6xW+1bnga+zYscTHx5OcnEx+fr4rcuXk5BAREUFqaiqpqak8+uijPs0zYsQIoqKiaNOmTaXH\n+Hucqsrk7zECKCgoICMjg8suu4zWrVvz7LPPVnicv8eqOrn8PV4//vgjaWlppKSkkJSUxIMPPljh\ncf4cq+pksvFzBeYapdTU1EoXttioU2fL5NE4eXSmwAMnTpxwYmNjnZ07dzolJSVOcnKys2XLllOO\nWbp0qdOjRw/HcRwnLy/PSUtLc0WuTz75xOndu7fPs/wsNzfX2bBhg9O6desKv29jnKrK5O8xchzH\n2bNnj5Ofn+84juMcOXLEadWqlSt+pqqTy8Z4HT161HEcxyktLXXS0tKcFStWnPJ9G2NVVSYb4+Q4\njvPkk086t9xyS4WvbWOcqsrkyTj5beZfnQu+Fi9ezLBhwwBIS0vj8OHDFBcXW88FlZ/U9oVOnTrR\noEGDSr9vY5yqygT+HSOAxo0bk/L/+/zWr1+fxMREdu/efcoxNsaqOrnA/+NVt665R3BJSQllZWVc\nfNpaTBtjVVUm8P84FRYWsmzZMn7/+99X+No2xqmqTHDu4+S34l/RBV9FRUVVHlNYWGg9V1hYGKtX\nryY5OZmsrCy2bNni00xVsTFOVbE9Rrt27SI/P5+0tLRTft/2WFWWy8Z4lZeXk5KSQlRUFBkZGSSd\ndvNbG2NVVSYb43TXXXcxbdo0atWquDzaGKeqMnkyTn4r/tW94Ov0dy9fXyhWnedv164dBQUFbNq0\niT/+8Y/07dvXp5mqw9/jVBWbY/TDDz/Qv39/nnnmGerXr3/G922N1dly2RivWrVqsXHjRgoLC8nN\nza1wuwJ/j1VVmfw9TkuWLKFRo0akpqaedSbtz3GqTiZPxslvxT86OpqCgoKTXxcUFND0tC38Tj+m\nsLCQ6Oho67kuvPDCkx9Pe/ToQWlpKQcPHvRprrOxMU5VsTVGpaWl9OvXj1tvvbXCH3hbY1VVLps/\nUxEREfTs2ZN169ad8vs2f64qy+TvcVq9ejWLFy+mRYsWDBo0iI8//pihQ4eecoy/x6k6mTwaJ89P\nP5yb0tJSp2XLls7OnTudn376qcoTvp9++qlfTqRUJ9fevXud8vJyx3EcZ82aNU6zZs18nmvnzp3V\nOuHrr3GqKpONMSovL3eGDBnijB8/vtJjbIxVdXL5e7z27dvnHDp0yHEcxzl27JjTqVMn56OPPjrl\nGH+PVXUy2fi5+llOTo7Tq1evM37f1r+/s2XyZJx8urHbL/3ygq+ysjJGjhxJYmIiM2bMAGDUqFFk\nZWWxbNky4uLiqFevHrNmzXJFrgULFvDCCy8QHh5O3bp1mTdvnk8zDRo0iOXLl7N//35iYmKYNGkS\npaWlJ/MDwkF5AAAAjUlEQVTYGKeqMvl7jABWrVrF66+/Ttu2bUlNTQXMvlPffvvtyVw2xqo6ufw9\nXnv27GHYsGGUl5dTXl7OkCFDuPbaa63++6tOJhs/V7/0czvHdp2qKpMn46SLvEREQpBu4ygiEoJU\n/EVEQpCKv4hICFLxFxEJQSr+IiIhSMVfRCQEqfiLiIQgFX8RkRD0f0EODBB4qr4PAAAAAElFTkSu\nQmCC\n",
"text": [
""
]
}
],
"prompt_number": 431
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def linemin_old(x0,d,f,alpha=1,eps=1e-6,just_result=False):\n",
" #golden section line minimization, starting from bracketing triple\n",
" x0,x1,x2 = bracket(x0,d,f,alpha=alpha)\n",
" x = array([[x0,x1,x2]])\n",
" l = lambda x: sum(x**2,axis=0)\n",
" while(l(x2-x0) > eps):\n",
" if l(x1-x0) > l(x2-x1):\n",
" x3 = x0 + (x1-x0)/phi*d\n",
" if f(x3)\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 24\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 25\u001b[0m \u001b[0;31m#test\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 26\u001b[0;31m \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlinemin\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcos\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 27\u001b[0m \u001b[0mt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mpi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 28\u001b[0m \u001b[0mfig\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfigure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdpi\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m150\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mfigsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m8\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mNameError\u001b[0m: name 'linemin' is not defined"
]
}
],
"prompt_number": 5
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"sgn = lambda x: (1, -1)[x<0]\n",
"def mnbrak(ax,bx,f,glimit=1e2,eps=1e-6):\n",
" #given distinct initial points ax and bx, search for a bracketing triple ax,bx,cx\n",
" #from numerical recipes, 10.1\n",
" if f(bx)>f(ax):\n",
" tmp = ax; ax = bx; bx = tmp; #switch\n",
" cx = bx + phi*(bx-ax);\n",
" while True:\n",
" r=(bx-ax)*(f(bx)-f(cx)) #do parabolic fit\n",
" q=(bx-cx)*(f(bx)-f(ax))\n",
" denom = q-r if q-r!=0 else eps*sgn(q-r)\n",
" u=bx-((bx-cx)*q-(bx-ax)*r)/(2*denom)\n",
" ulim=bx+glimit*(cx-bx)\n",
" if (bx-u)*(u-cx) > 0:\n",
" if f(u) < f(cx):\n",
" return bx,u,cx\n",
" elif f(u) > f(bx):\n",
" return ax,bx,u\n",
" u = cx+phi*(cx-bx) #parabolic fit didn't help, use default magnification\n",
" elif (cx-u)*(u-ulim) > 0:\n",
" if f(u) < f(cx):\n",
" bx=cx\n",
" cx=u\n",
" u = cx+phi*(cx-bx)\n",
" elif (u-ulim)*(ulim-cx) > 0:\n",
" u = ulim\n",
" else:\n",
" u = cx+phi*(cx-bx)\n",
" ax=bx; bx=cx; cx=u;\n",
" if f(bx) < f(cx): break\n",
" return ax,bx,cx"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 222
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"phi = .5*(1+sqrt(5))\n",
"x = array([mnbrak(0,1,cos)])\n",
"t = linspace(0,2*pi,100)\n",
"fig = plt.figure()\n",
"#ax = plt.axes(xlim=(0, 2*pi), ylim=(-1.1, 1))\n",
"plt.plot(t,cos(t))\n",
"plt.plot(x,cos(x),linestyle='',marker='.',color='r')\n",
"print x"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"[[ 1. 2.61803399 5.23606798]]\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlY1WX+xvE3iGZqipViIxSGC2iOYji0WViapRNqaa5p\n7jWZbVNNNZU25ZItrpk6/Vwyt5xMM6R0DHVcc2kxV1QMqZjRCXMNgfP740knTfTAWZ7vOed+Xde5\nivzSuTP88OFZw1wulwsREQla4bYDiIiIb6nQi4gEORV6EZEgp0IvIhLkVOhFRIKcCr2ISJDzuND3\n7t2bqKgoGjZsWOwzgwYNok6dOjRq1IjNmzd7+pYiIlICHhf6Xr16kZ6eXuyvp6WlkZmZya5du5g0\naRIPPvigp28pIiIl4HGhb9asGVWrVi321xcuXEjPnj0BSE5OJi8vj9zcXE/fVkRE3OTzMfqcnBxi\nYmJOfxwdHc3+/ft9/bYiIvILv0zGnn3KQlhYmD/eVkREgAhfv0HNmjXJzs4+/fH+/fupWbPmb567\n8sraZGfv9nUcEZGgEhcXR2Zm5nmf8XlHn5qayvTp0wFYu3YtkZGRREVF/ea57OzduFyuM175+S4W\nLXJx770uqlRx0b+/i2+/df3mOSe8XnzxResZQjG78tt/BUv+zEwXPXqYWtOtm4tPPnFRUGA/34Ve\nu3dfuEH2uKPv0qULy5cv58CBA8TExDBkyBBOnjwJwIABA2jdujVpaWnUrl2bihUrMmXKFLf/3WXL\nQps25nXwILz+OjRuDF27wgsvQLVqnqYXkVB3+DD07w8ffAADB0JWFkRG2k7lXR4X+lmzZl3wmXHj\nxnn6Nlx2GQwdCo88Yv7asCGMGgWdOoGG/EWkpFwumDIFJkyAhx+GnTvh0kttp/INn4/Re1tUFIwe\nbbr63r1h9myYONH8c5tSUlLsBvBAIGcH5bctEPNnZ0PfvvCf/8CYMSn07Ws7kW+FuVwuR1w8EhYW\nRkmj/PwzDB4MM2bAvHmQnOybbCISPDIyoEsXeOghePppM0QcyNypnQFd6E9ZuNB8dx42DPr08XIw\nEQkKLheMGWPqxIwZ0KKF7UTeETKFHmD7dmjfHu68E157DcJ1XJuI/KKwEB58ED7/3Ey61qplO5H3\nhFShB/jxR7NCp149mDwZIgJuBkJEvC0/H7p3h//+Fz78ECpVsp3Iu9ypnUHV91atCkuWQE6OWY3z\n88+2E4mITUePQmoqnDwJixYFX5F3V1AVeoCKFeGjj8x4XMeO5n+wiISeEydMka9eHd5/H8qXt53I\nnqAauvm1/HwzZl+1KkyfrjF7kVBSUAAdOsBFF8HMmVCmjO1EvhNyQze/Vq6c+S6enW02Qzjj25mI\n+FpRkdljk58P774b3EXeXUFb6AEqVDDDOOvWwUsv2U4jIv7w5JPmGIN580zDJwG4M7akKleGjz82\nm6ni480krYgEp8mTzaTr2rWm0RMjaMfoz/bll9CypSn6TZv67G1ExJKMDNPIrVwJdevaTuM/IT1G\nf7ZGjcx3+/btQRdciQSXzExT5N97L7SKvLtCpqM/ZdgwWLAAVqzQ+J1IMDh2zAzNPvCAOb8m1ITc\nzlh3FBVBu3ZQuza88YbP305EfKxPH7NmfsaM0Dyy3J3aGfSTsWcLD4epU6FJE2jWzAzliEhgmjYN\nVq82Z9iEYpF3V8h19KesWwd33WVm56++2m9vKyJe8s03kJICn30G11xjO409mow9j+RkePZZcy51\nQYHtNCJSEidOQOfOMGJEaBd5d4VsRw9mvP6OO8wQzvPP+/WtRcQDTz4Je/aYTVGhPmSjyVg35ORA\nYiKkpUFSkt/fXkRKaPlyc5Xol1/C5ZfbTmOfhm7cULOmuXXmvvvMMi0Rca5Dh6BnT5g0SUW+JEK+\noz+la1dznOmoUdYiiMgF9O5t7nidONF2EufQ0E0JHDxoJnXmz4frrrMWQ0SKsWSJuRt6yxa45BLb\naZxDQzclcNllppvv29ccbyoiznH0KAwYAG+/rSJfGurof8XlgrZtzaTsCy9YjSIiv/LEE5Cba3a/\nypk0dFMK+/ebVTjLl0P9+rbTiMjnn5vNjV9/DdWq2U7jPBq6KYXoaHNJSf/+upVKxLaCAvNn8bXX\nVOQ9oUJ/Dv37w88/m2vIRMSet9+GyEjo1s12ksCmoZtifP65uUF+2zbzhSYi/pWba1bCZWRAgwa2\n0ziXxug9NGCAuUV+zBjbSURCz/33m+GakSNtJ3E2FXoPHTxoJmQ/+QQaN7adRiR0/Otf5tCybdu0\nnPJCNBnrocsug7/9DQYO1MSsiL8UFsLDD5sJWBV571Chv4A+fcxmjXnzbCcRCQ1Tp0KlSuYOWPEO\nDd244bPPzBkb27ZB+fK204gEr8OHoV49c69z06a20wQGDd14SfPm0KgRjB5tO4lIcBs+HFq2VJH3\nNnX0btq1C66/3lxfFhVlO41I8MnKgmuvha++MseHi3u06sbLHn8cjhwxZ2GLiHd16QLx8fDii7aT\nBBYVei/78UeoWxdWrjRfkCLiHRs2mAMFd+6EihVtpwksGqP3sqpV4amnzKXiIuIdLhc8/bTp5FXk\nfUOFvoQGDjTHI6xZYzuJSHBYssScGtu7t+0kwUuFvoQuvhiGDDEdiMNHmkQcr6jI/FkaOhQiImyn\nCV4q9KXQo4c5HiEtzXYSkcA2Zw6UKwd33207SXDTZGwpLVgAzz8PX3wB4fp2KVJiJ0+aRQ3vvAMp\nKbbTBC5NxvpQaqoZxnn/fdtJRALT1KkQF6ci7w/q6D2wZIk5fGnLFo0vipTEiRNmqfL770Nysu00\ngU0dvY+1aGF2yb73nu0kIoFl8mRz9LeKvH+oo/fQihWQ06Y/nRJ3El6pAsycqSupRM7j2DGoXdss\nZtA9D57zS0efnp5OfHw8derUYcSIEb/59YyMDKpUqUJiYiKJiYm8/PLLnr6lo9x8M9SP2En4yuWw\neLG5cFZEijV+PNx4o4q8P3k0slxYWMjAgQNZunQpNWvWpGnTpqSmppKQkHDGc7fccgsLFy70KKiT\nXZlQAdZA0bVJhOsgHJFiHTliLhRZtsx2ktDiUUe/fv16ateuTWxsLGXLlqVz584sWLDgN88F4pBM\nSVRNm8mKqI5M7bpEwzYi5/HWW+bYb1327V8eFfqcnBxiYmJOfxwdHU1OTs4Zz4SFhbF69WoaNWpE\n69at2bp1qydv6UyRkZT7cC5DRkeSn287jIgznezVn5v+msLknNaQl2c7TkjxaOgmLCzsgs80adKE\n7OxsKlSowOLFi2nXrh07d+4857ODBw8+/fcpKSmkBNAC2+uuM5s/pk2Dfv1spxFxnh9W7OSGk8vh\nX5i5rLlzbUcKSBkZGWRkZJToczxadbN27VoGDx5Meno6AMOGDSM8PJynn3662M+pVasWGzdu5NJL\nLz0zSICuuvm1Vauge3dz1GrZsrbTiDjHsWOwpmprbstfDElJZhOKhjm9wuerbpKSkti1axdZWVnk\n5+czZ84cUlNTz3gmNzf3dIj169fjcrl+U+SDxY03mp1+775rO4mIs0yaBFNazoSOHVXkLfBo6CYi\nIoJx48bRqlUrCgsL6dOnDwkJCUycOBGAAQMGMG/ePCZMmEBERAQVKlRg9uzZXgnuVC++CPffbw4+\n025ZEbMLduRIWLQoEhI1XGODNkz5wM03w4AB0K2b7SQi9k2YAB9/DIsW2U4SnHSVoCWffAJPPGEu\nOdbJlhLKTp40Z9rMnAnXX287TXDSWTeW3H47lC8PQbxHTMQts2ZBrVoq8rapo/eRDz6A4cNh3Tpw\nYxWqSNApKjIbo8aONQcAim+oo7eoXTs4fBiWLrWdRMSO+fOhcmW47TbbSUSF3kfCw+GZZ+CVV2wn\nEfE/l8vcA/vss/qJ1glU6H2oSxfIyjLDNyKhZOlSs6zyrrtsJxFQofepsmXh8cfhHKc3iwS1ESPg\nqae06swpNBnrY0ePmlUHK1dCvXq204j43saN0L49ZGZCuXK20wQ/TcY6QMWK8Kc/mTO4RULBiBHw\n2GMq8k6ijt4PDhwwm0a++QauuMJ2GhHfycw0a+b37oVKlWynCQ3q6B3i8svNqZajRtlOIuJbr71m\njv9QkXcWdfR+sm8fNGkCe/ZAlSq204h4X26uuZNhxw6oXt12mtChjt5BrroKWrWCv//ddhIR3xg/\nHjp1UpF3InX0frRxo9kxu2ePLiaR4HLsGMTGanWZDeroHebaa6FOHd2gJsFn6lS44QYVeadSR+9n\naWnw3HOwaZO2hktwKCw0BX7aNHPLmviXOnoHuuMOyM+Hf/7TdhIR71iwAKpVMx29OJMKvZ+Fh5tL\nSUaOtJ1ExDteew3+/Gf9hOpkKvQWdOtmbp/assV2EhHPrFljllW2a2c7iZyPCr0FF10EDz2kDVQS\n+N58Ex59FMqUsZ1EzkeTsZYcOGBW4GzfDlFRttOIlFxWFiQlmb9qJ6w9mox1sMsvh3vvhQkTbCcR\nKZ0xY6B3bxX5QKCO3qJt2yAlxRyPUL687TQi7jt0CK6+Gr74AmJibKcJberoHS4hwfzo+957tpOI\nlMw775gjPVTkA4M6esv++U8YNMiswNHyNAkEBQUQFwf/+IdpVMQudfQB4NZbISICliyxnUTEPfPn\nm05eRT5wqNBbFhZmlqdpqaUEilGjzA1SEjg0dOMAJ06Yk/8yMsx53iJOtX69WS2WmWl+EhX7NHQT\nIMqXN7fyjBljO4nI+Y0eDQ8/rCIfaNTRO8QPP5hVOLt3w6WX2k4j8ls5OdCwoblPITLSdho5RR19\nAKlRA1JTYfJk20lEzu2tt8w5TSrygUcdvYNs2gRt2+oGKnGe48fNdZirVpmjO8Q51NEHmCZNoFYt\ns3xNxEneew/+8AcV+UClQu8wjzxiJrxEnMLlMl+TjzxiO4mUlgq9w7RtC/v3w4YNtpOIGJ99BkVF\n0KKF7SRSWir0DhMRAQMHqqsX5xg92hzToSM6ApcmYx3oxx/NyYBbt8IVV9hOI6Fs92647jpzwmqF\nCrbTyLloMjZAVa0KnTvDxIm2k0ioGz/enDmvIh/Y1NE71Nat5sCzffvM1YMi/nb4sDmaY/NmuPJK\n22mkOOroA1j9+vD738PcubaTSKiaPh2aN1eRDwYq9A52aqmlftARfysqgrFjtaQyWKjQO9idd0Je\nHqxdazuJhJpPP4WLL4abbrKdRLxBhd7BwsPNSYE61VL8bcwYLakMJpqMdbhDh8yxCF9/DTVr2k4j\noWDnTtPJf/utLq0PBJqMDQJVqpgTA99+23YSCRXjxkG/firywUQdfQDYsQNuvtkstdQfPvGln34y\nSyq/+gqio22nEXeoow8S9eqZky3nzLGdRILd1KnmTBsV+eDicaFPT08nPj6eOnXqMGLEiHM+M2jQ\nIOrUqUOjRo3YvHmzp28ZkgYN0lJL8a1TSyoHDbKdRLzNo0JfWFjIwIEDSU9PZ+vWrcyaNYtt27ad\n8UxaWhqZmZns2rWLSZMm8eCDD3oUOFS1agVHjsDq1baTSLBKT4dLLoEbb7SdRLzNo0K/fv16ateu\nTWxsLGXLlqVz584sWLDgjGcWLlxIz549AUhOTiYvL4/c3FxP3jYkhYebUy211FJ85VQ3ryWVwcej\nQp+Tk0NMTMzpj6Ojo8nJybngM/v37/fkbUPW/ffDkiXmvHoRb9qxw1xl2bmz7STiCxGefHKYm9/6\nz54RLu7zBg8efPrvU1JSSElJKW20oFS58v+WWr78su00EkzGjYO+fbWqKxBkZGSQkZFRos/xqNDX\nrFmT7Ozs0x9nZ2cTfdZ0/dnP7N+/n5rF7Pz5daGXcxs40Cy1/Otf9YdSvOOnn8ydsF99ZTuJuOPs\nJnjIkCEX/ByPhm6SkpLYtWsXWVlZ5OfnM2fOHFJTU894JjU1lenTpwOwdu1aIiMjiYqK8uRtQ5qW\nWoq3aUll8POoo4+IiGDcuHG0atWKwsJC+vTpQ0JCAhN/uTFjwIABtG7dmrS0NGrXrk3FihWZMmWK\nV4KHskGD4LnnoEcPTZyJZ04tqZw61XYS8SXtjA1ARUUQHw9TpmgpnHgmLQ2ef95cRq+mITBpZ2yQ\n0lJL8RadUhka1NEHKJ1JIp7avh1uuUVnKAU6dfRBrHJl6N4dJkywnUQC1bhx0L+/inwoUEcfwHbu\nhGbN1JFJyemeg+Chjj7I1a0L114Ls2bZTiKBZsoUuP12FflQoY4+wKWnwzPPmO3rmlATdxQWmiZh\nxgy4/nrbacRT6uhDwO23w/HjsHKl7SQSKNLS4LLL4LrrbCcRf1GhD3C6QFxKSksqQ4+GboLA4cNm\nqeWmTXDVVbbTiJN984057iArCy66yHYa8QYN3YSISy6Bnj3hrbdsJxGnGzsWHnhART7UqKMPErt3\nQ3IyfPstVKhgO4040Y8/wtVXw7ZtUKOG7TTiLeroQ0hcnDn35t13bScRp5o8Gf74RxX5UKSOPogs\nW2YmZrds0USbnKmgwDQDH3xg9l5I8FBHH2KaN4cyZWDpUttJxGk+/BBiYlTkQ5UKfRAJC4NHHoFR\no2wnEacZNQoefdR2CrFFQzdB5vhxs9Ry5Uqz+1Fkwwa45x4zYR/h0VVD4kQauglBF18M/fppA5X8\nz+jR8NBDKvKhTB19EPruO7jmGtizByIjbacRm77/HurXN938pZfaTiO+oI4+RP3ud9C6Nfz977aT\niG1vvQWdO6vIhzp19EFK47Jy/Lg5EmPlSqhXz3Ya8RV19CEsKQmuvNKsm5bQ9N570LSpiryo0Ae1\nxx6DN9+0nUJscLnMksrHHrOdRJxAhT6ItW0Lubmwdq3tJOJvS5aYI6xvu812EnECFfogVqaMOXdc\nXX3oefNNs0FKR2EIaDI26P30k7kEevNmM2YvwW/bNnMcRlaWLo0PBZqMFSpXNmfVawNV6HjzTXPm\nvIq8nKKOPgTs2wdNmsDevabwS/D697/NKpsdO6B6ddtpxB/U0Qtg1lK3bAnvvGM7ifjaW29Bx44q\n8nImdfQh4vPPTQHIzNQGqmB16kC75cshPt52GvEXdfRyWtOmZjL2H/+wnUR85d13zf9nFXk5mwp9\nCHniCXj9dbOZRoJLURG88Yb5fyxyNhX6EPLHP0JeHvzrX7aTiLelpZlL4VNSbCcRJ1KhDyFlysDj\nj8PIkbaTiLeNHAlPPqkNUnJumowNMcePmw1Uy5aZc8ol8K1bB506aaI9VGkyVn7j4ovNbUOvv247\niXjLyJHmJzUVeSmOOvoQdPAg1KkDW7aYS0okcO3aBTfcYI47qFjRdhqxQR29nNNll0H37joWIRi8\n8YY57kBFXs5HHX2IysqCa6/VsQiBTMcdCKijl/OIjYVWrWDiRNtJpLTGjjWTsCryciHq6EPYl1+a\nS8T37IGLLrKdRkrip5/g6qvNipu4ONtpxCZ19HJejRqZ17RptpNISU2aZA6qU5EXd6ijD3ErV0Kv\nXrB9u5bnBYqffzbdfFqa+UYtoU0dvVxQs2ZQo4YOOwsk06b976cxEXeooxcWLYK//tVcN6gt9M5W\nWGhW2kyZYr5Ji6ijF7e0aWNOP0xPt51ELmTePIiKgptusp1EAokKvRAWBs8+Cy+/rCOMnayoCF55\nBZ57Tj95Scmo0Atgbp86cAAyMmwnkeJ89BGULQt33mk7iQSaUo/R//e//6VTp07s27eP2NhY5s6d\nS2Rk5G+ei42NpXLlypQpU4ayZcuyfv36cwfRGL1106aZ17JltpPI2Vwuc3vUc89B+/a204iT+HSM\nfvjw4bRs2ZKdO3dy2223MXz48GJDZGRksHnz5mKLvDhD167maIRVq2wnkbN98olZVtm2re0kEohK\nXegXLlxIz549AejZsycffvhhsc+qUw8MZcvCX/5ixurFOVwu+NvfTDcfrsFWKYVSf9nk5uYSFRUF\nQFRUFLm5ued8LiwsjBYtWpCUlMTkyZNL+3biJz17muOLN2ywnUROycgw8ycdO9pOIoHqvHshW7Zs\nyQ8//PCbf/7KK6+c8XFYWBhhxSwDWLVqFVdccQX/+c9/aNmyJfHx8TQrZgHw4MGDT/99SkoKKboA\n0+8uugiefhqGDDGTf2KXywWDB5tuvkwZ22nECTIyMsgo4aqJUk/GxsfHk5GRQY0aNfj+++9p3rw5\n27dvP+/nDBkyhEqVKvHEOa6q12Ssc5w4YS4m+eADMwEo9vzzn/Dgg7B1q46okHPz6WRsamoq0345\nDWvatGm0a9fuN88cO3aMw4cPA3D06FE+/fRTGjZsWNq3FD8pXx6eecZ0kmKPywUvvggvvKAiL57x\naHnlvffey7fffnvG8srvvvuOfv368fHHH7Nnzx7uvvtuAAoKCujWrRvPPPPMuYOoo3eUn382Xf37\n70Nysu00oWnJEnj4YfjmGw3bSPHcqZ0660aKNXEizJ+voxFscLngxhtNoe/SxXYacTKddSMeOXV8\n8erVtpOEnk8+gbw8uPde20kkGKjQS7HKlYPnnzcnW+qHLf9xuczv+eDBGrIR71Chl/Pq2RNycszq\nD/GPDz4wB5h16GA7iQQLjdHLBc2dC6+9Zu4n1amJvlVQAA0bwptvwh132E4jgUBj9OIVHTrAyZNm\nYlZ86913oVo1aNXKdhIJJuroxS2LF8MTT8DXX2vc2Fd+/hnq1oWZM82KGxF3qKMXr7njDrj8cpg+\n3XaS4PX222bYRkVevE0dvbhtzRqz3G/HDqhQwXaa4JKXZ+6CXbrUFHsRd6mjF6+6/nq47joYNcp2\nkuAzbBj88Y8q8uIb6uilRDIzTbHfuhWqV7edJjjs2wdNmsBXX0HNmrbTSKDREQjiE488YpYBjh9v\nO0lwuO8+qFULXnrJdhIJRCr04hMHDkB8vLlysF4922kC26ZN0KYN7NwJl1xiO40EIo3Ri09cfrm5\nnOQc1wpICbhc8Nhj5qgDFXnxJRV6KZVBg8zqm8WLbScJXHPnwqFD0Lev7SQS7DR0I6X28cfw+ONm\nE1W5crbTBJajRyEhAd57D4q5WVPELRq6EZ9q0wbi4mDsWNtJAs+IEXDDDSry4h/q6MUjO3bATTeZ\nrr5GDdtpAsPevZCUBF98ATExttNIoNOqG/GLp56C3Fz45QphuYD27c26+eeft51EgoEKvfjFkSPQ\noIEp9CkpttM428KF8OSTZnPURRfZTiPBQGP04heVKsGYMfDAA+YERjm3I0fMHbBvv60iL/6lQi9e\n0bat2UT16qu2kzjX4MFwyy3QvLntJBJqNHQjXpOdDYmJ5pTLOnVsp3GWL7+Eli1hyxadESTepaEb\n8auYGHjuOejf39x5KkZBgdkUNWyYirzYoUIvXjVokBmnnzDBdhLnGDkSqlaF3r1tJ5FQpaEb8bpT\na+vXrYOrr7adxq4tW8yY/MaNcOWVttNIMNLQjVhRrx785S+mgw3lIZyCAujVC4YOVZEXu1ToxSce\nfRTy80P7zPpXXzVDNjq0TGzT0I34zK5d5jyXzz6Da66xnca/1q83VwNu2KBuXnxLQzdiVZ06pqvt\n3BmOH7edxn9++gm6dDET0iry4gTq6MWnXC5T9C67LHSGce67Dy6+GCZNsp1EQoE7tTPCT1kkRIWF\nmS3/iYmwYIHZQRvMZswwwzUbNthOIvI/6ujFL9asMUV+9WqoXdt2Gt/4+mu49VZYuhQaNbKdRkKF\nxujFMa6/3pz10r69Odwr2Pz4o/lve/NNFXlxHnX04jcul1lbf+wYzJ5thnWCQWEh3HWXmXwePdp2\nGgk16ujFUcLCzEqU3bvNsQDB4sUXzR2wr71mO4nIuamjF7/LzoYV8f1pedVOqsdWgJkzITLSdqxS\n+b//g5dfNnMQUVG200go0g1T4lhHklKotHG5+aBjR5g7126gUkhPh/vvh+XLzbEPIjZo6EYcq1L1\nCgBsjkhi558Db8H5pk3Qowd88IGKvDifCr3YMXMmdOzI9jFLaNkxkt27bQdy3zffmOMN3n7bHPEg\n4nTaMCV2REbC3Ll0AX4KN+vPP/vM+ccab91qbooaORLuvtt2GhH3qNCLdQMGmKWXzZs7u9hv3Qot\nWpjze7p1s51GxH0q9OIIDzxgiv0tt8BHH0HjxrYTnWnNGtPBv/oqdO9uO41IyWiMXhzjwQfhjTfM\n0MjixbbT/M8//gGpqfD3v5sDy0QCjZZXiuOsXm265xdeMMXf1g5al8tsgho92vyUkZhoJ4fI+Wgd\nvQSszExT7BMSYOJE/++nOnDAXAP4ww9mCWVMjH/fX8RdWkcvAat2bXO5eLVqppNevdp/771smZkj\nqF8fVq1SkZfAV+pC//7779OgQQPKlCnDpk2bin0uPT2d+Ph46tSpw4gRI0r7dhKCLr4Yxo2DUaPg\nnnvM3av//rfv3i8nx0y09uhhjjYYMQLKlfPd+4n4S6kLfcOGDZk/fz4333xzsc8UFhYycOBA0tPT\n2bp1K7NmzWLbtm2lfUtHy8jIsB2h1JyevW1b2L4dqlSBBg3g9dfPPOrY0/yHDsHQoeZ44auuMu91\n++2eZS4Jp//+X4jyO1+pC318fDx169Y97zPr16+ndu3axMbGUrZsWTp37syCBQtK+5aOFshfLIGQ\nvUoVU+BXrDBLHWNj4cknYd++0uffvRseeQRq1YItW8xQ0SuvQKVKXo1+QYHw+38+yu98Pl1Hn5OT\nQ8yvBjijo6NZt26dL99SglxCAsybB3v3wtixZvy+fHnza61amY6/cuVzf25enino6emwaJEZqunT\nB776CqKj/fffIOJv5y30LVu25IcffvjNPx86dCh33XXXBf/lYcFys4Q4Tq1aZs39iBHQr585D/6h\nh2DHDrjkEvPrZcuaZ/PzYc8eOH4c4uPN7tbx4yE5GSK0ZVBCgctDKSkpro0bN57z19asWeNq1arV\n6Y+HDh3qGj58+DmfjYuLcwF66aWXXnqV4BUXF3fBOu2VfsZVzBrOpKQkdu3aRVZWFr/73e+YM2cO\ns2bNOuezmZmZ3ogiIiJnKfVk7Pz584mJiWHt2rW0adOGO++8E4DvvvuONm3aABAREcG4ceNo1aoV\n9evXp1MCJ8dbAAAEWklEQVSnTiQkJHgnuYiIuMUxO2NFRMQ3rO+MDeQNVb179yYqKoqGDRvajlIq\n2dnZNG/enAYNGnDNNdcwZswY25FK5MSJEyQnJ9O4cWPq16/PM888YztSiRUWFpKYmOjW4gYnio2N\n5fe//z2JiYn84Q9/sB2nRPLy8ujQoQMJCQnUr1+ftWvX2o7kth07dpCYmHj6VaVKlfP/+S3F/KvX\nFBQUuOLi4lx79+515efnuxo1auTaunWrzUglsmLFCtemTZtc11xzje0opfL999+7Nm/e7HK5XK7D\nhw+76tatG1C//y6Xy3X06FGXy+VynTx50pWcnOxauXKl5UQl8/rrr7u6du3quuuuu2xHKZXY2FjX\nwYMHbccolR49erjeeecdl8tlvn7y8vIsJyqdwsJCV40aNVzffvttsc9Y7egDfUNVs2bNqFq1qu0Y\npVajRg0a/3Lwe6VKlUhISOC7776znKpkKlQwd8/m5+dTWFjIpZdeajmR+/bv309aWhp9+/YN6AP9\nAjH7oUOHWLlyJb179wbMfGKVKlUspyqdpUuXEhcXd8aepbNZLfTn2lCVk5NjMVHoysrKYvPmzSQn\nJ9uOUiJFRUU0btyYqKgomjdvTv369W1Hcttjjz3GyJEjCQ+3PoJaamFhYbRo0YKkpCQmT55sO47b\n9u7dS7Vq1ejVqxdNmjShX79+HDt2zHasUpk9ezZdu3Y97zNWv8K0ocoZjhw5QocOHRg9ejSV/L3/\n30Ph4eF88cUX7N+/nxUrVgTMdvZFixZRvXp1EhMTA7IjPmXVqlVs3ryZxYsXM378eFauXGk7klsK\nCgrYtGkTf/rTn9i0aRMVK1Zk+PDhtmOVWH5+Ph999BEdO3Y873NWC33NmjXJzs4+/XF2djbR2ovu\nVydPnuSee+6he/futGvXznacUqtSpQpt2rRhw4YNtqO4ZfXq1SxcuJBatWrRpUsXli1bRo8ePWzH\nKrErrrgCgGrVqtG+fXvWr19vOZF7oqOjiY6OpmnTpgB06NDhvKfwOtXixYu59tprqVat2nmfs1ro\nf72hKj8/nzlz5pCammozUkhxuVz06dOH+vXr8+ijj9qOU2IHDhwgLy8PgOPHj7NkyRISA+QaqKFD\nh5Kdnc3evXuZPXs2t956K9OnT7cdq0SOHTvG4cOHATh69CiffvppwKxAq1GjBjExMezcuRMw49wN\nGjSwnKrkZs2aRZcuXS74nNWTPn69oaqwsJA+ffoE1IaqLl26sHz5cg4ePEhMTAwvvfQSvXr1sh3L\nbatWrWLGjBmnl8cBDBs2jDvuuMNyMvd8//339OzZk6KiIoqKirjvvvu47bbbbMcqlUAcxszNzaV9\n+/aAGQrp1q0bt/vzfGcPjR07lm7dupGfn09cXBxTpkyxHalEjh49ytKlS92aG9GGKRGRIBe40/0i\nIuIWFXoRkSCnQi8iEuRU6EVEgpwKvYhIkFOhFxEJcir0IiJBToVeRCTI/T9TqXjlkYjS/AAAAABJ\nRU5ErkJggg==\n",
"text": [
""
]
}
],
"prompt_number": 8
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def linemin(xi,d,f,alpha=1,eps=1e-3,just_result=False):\n",
" #golden section line minimization\n",
" fs = lambda a: f(xi+a*d) #scalar\n",
" x0,x1,x2 = mnbrak(0,alpha,fs)\n",
" #x0,x1,x2 = map(lambda a: x0+a*d,mnbrak(0,alpha,fs))\n",
" x = array([[x0,x1,x2]])\n",
" #l = lambda x: sum(x**2,axis=0)\n",
" i=0\n",
" while(x2-x0 > eps):\n",
" #print 'linemin',l(x2-x0)\n",
" i = i+1\n",
" if x1-x0 > x2-x1:\n",
" x3 = x0 + (x1-x0)/phi\n",
" if fs(x3)abs(bx-ax):\n",
" x1 = bx; x2 = bx+C*(bx-ax)\n",
" else:\n",
" x2 = bx; x1 = bx-C*(bx-ax)\n",
" f1 = fs(x1); f2= fs(x2)\n",
" while abs(x3-x0) > tol*(abs(x1)+abs(x2)) and abs(x1)+abs(x2) > tol:\n",
" if f2"
]
}
],
"prompt_number": 247
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def powell(x0,f,alpha=1.,eps=1e-4,max_iter=100):\n",
" #powell direction set method, start x0, function f\n",
" x = array([x0])\n",
" n_dim = shape(x)[-1]\n",
" d = eye(n_dim) #take coordinate axes as initial directions\n",
" #it is important to make sure these are descent directions, else linemin will fail\n",
" i=0\n",
" while(True):\n",
" df = [] #changes in objective along each direction\n",
" i = i+1\n",
" for di in d:\n",
" a_min = golden(x[-1],di,f,alpha=alpha) #minimize in this direction\n",
" xnew = x[-1]+a_min*di\n",
" x = vstack((x,array([xnew])))\n",
" df.append(f(x[-2])-f(x[-1])) #save change in objective\n",
" if max(df) < eps or i>max_iter: break\n",
" d[argmax(df)] = sum(d,axis=0)/n_dim #replace old direction with overall change\n",
" #d[argmax(df)] /= linalg.norm(d[argmax(df)])\n",
" return x"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 267
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def plot_ros_powell(pts,n_frames = 20, figsize=(8, 6), dpi=80):\n",
" n_pts = shape(pts)[0]\n",
" fig = plt.figure(figsize=figsize, dpi=dpi)\n",
" ax = fig.add_subplot(111, projection='3d')\n",
" ax.axis('off')\n",
" ax.set_xlim((-1, 1))\n",
" ax.set_ylim((-1, 1))\n",
" N = 100;\n",
" x = linspace(-1.5,2,N); y = linspace(-1.5,2,N);\n",
" X,Y = meshgrid(x,y)\n",
" Z = rosenbrock(dstack((X,Y)))\n",
" x_plot, = ax.plot([], [], [], marker='.',linestyle='-', c='r')\n",
" def init():\n",
" ax.plot_wireframe(X, Y, Z, rstride=10, cstride=10)\n",
" x_plot.set_data([], [])\n",
" x_plot.set_3d_properties([])\n",
" def animate(i):\n",
" ax.azim = 360/n_frames*i\n",
" ax.elev = 30\n",
" xi = ceil(i*n_pts/n_frames)\n",
" x_plot.set_data(pts[:xi,0],pts[:xi,1]);\n",
" x_plot.set_3d_properties(rosenbrock(pts[:xi]))\n",
" return animation.FuncAnimation(fig, animate, init_func=init,\n",
" frames=n_frames, interval=20, blit=True)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 268
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"x=powell([-1,-1],rosenbrock,alpha=.1,eps=1e-5,max_iter=100);\n",
"anim = plot_ros_powell(x,n_frames=50,dpi=150,figsize=(10,8))\n",
"#anim.save('direction-set.gif', writer='imagemagick', fps=10)\n",
"display_animation(anim)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
""
],
"metadata": {},
"output_type": "pyout",
"prompt_number": 270,
"text": [
""
]
}
],
"prompt_number": 270
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def mag(x):\n",
" return sqrt(sum(x**2,axis=0))\n",
"def conj_grad(f,df,x0,alpha=.1,eps=1e-6,max_iter=100):\n",
" dnew = -df(x0) #take starting direction as steepest descent\n",
" dnew /= mag(dnew)\n",
" x1 = x0 + golden(x0,dnew,f)*dnew\n",
" x = array([x0,x1])\n",
" i=0\n",
" def gamma(x):\n",
" #g= dot(df(x[-1]),df(x[-1])) / dot(df(x[-2]),df(x[-2])) #Fletcher\u2013Reeves\n",
" g= dot(df(x[-1]), df(x[-1])-df(x[-2]))/dot(df(x[-2]), df(x[-2])) #Polak\u2013Ribi\u00e8re\n",
" g= max(0,g) #provide steepest descent reset with polak-ribiere\n",
" return g\n",
" while(True):\n",
" i=i+1\n",
" dnew = -df(x[-1]) + gamma(x)*dnew\n",
" dnew /= mag(dnew)\n",
" a = golden(x[-1],dnew,f,alpha=alpha)\n",
" x = vstack((x,array([x[-1]+a*dnew])))\n",
" if abs(a)max_iter: break\n",
" return x\n",
"def plot_ros_cg(pts,n_frames = 20, figsize=(8, 6), dpi=80):\n",
" n_pts = shape(pts)[0]\n",
" fig = plt.figure(figsize=figsize, dpi=dpi)\n",
" ax = fig.add_subplot(111, projection='3d')\n",
" ax.axis('off')\n",
" ax.set_xlim((-1, 1))\n",
" ax.set_ylim((-1, 1))\n",
" N = 100;\n",
" x = linspace(-1.5,2,N); y = linspace(-1.5,2,N);\n",
" X,Y = meshgrid(x,y)\n",
" Z = rosenbrock(dstack((X,Y)))\n",
" x_plot, = ax.plot([], [], marker='.',linestyle='-', c='r')\n",
" def init():\n",
" ax.plot_wireframe(X, Y, Z, rstride=10, cstride=10)\n",
" x_plot.set_data([], [])\n",
" x_plot.set_3d_properties([])\n",
" def animate(i):\n",
" ax.azim = 360/n_frames*i\n",
" ax.elev = 30\n",
" xi = ceil(i*n_pts/n_frames)\n",
" x_plot.set_data(pts[:xi,0],pts[:xi,1]);\n",
" x_plot.set_3d_properties(rosenbrock(pts[:xi]))\n",
" return animation.FuncAnimation(fig, animate, init_func=init,\n",
" frames=n_frames, interval=20, blit=True)\n",
"def plot_ros_cg_2d(pts, f,n_frames = 20, figsize=(8, 6), dpi=80):\n",
" n_pts = shape(pts)[0]\n",
" fig = plt.figure(figsize=figsize, dpi=dpi)\n",
" ax.set_xlim((-1, 1))\n",
" ax.set_ylim((-1, 1))\n",
" N = 100;\n",
" x = linspace(-1.5,2,N); y = linspace(-1.5,2,N);\n",
" X,Y = meshgrid(x,y)\n",
" Z = f(dstack((X,Y)))\n",
" plt.contour(X, Y, Z, levels=[.3,1.,2.,10.,50])\n",
" x_plot, = plot([], [], c='g',marker='.')\n",
" def init():\n",
" x_plot.set_data([], [])\n",
" def animate(i):\n",
" xi = ceil(i*n_pts/n_frames)\n",
" x_plot.set_data(pts[:xi,0],pts[:xi,1])\n",
" return animation.FuncAnimation(fig, animate, init_func=init,\n",
" frames=n_frames, interval=20, blit=True)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 271
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"a = [1,10]\n",
"def quad(v):\n",
" if len(shape(v))==1:\n",
" return a[0]*v[0]**2 + a[1]*v[1]**2\n",
" else:\n",
" return a[0]*v[...,0]**2 + a[1]*v[...,1]**2\n",
"def d_quad(v): \n",
" if len(shape(v))==1:\n",
" return array([2*a[0]*v[0],2*a[1]*v[1]])\n",
" else:\n",
" return array([2*a[0]*v[...,0],2*a[1]*v[...,1]])"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 225
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"x = conj_grad(rosenbrock,d_rosenbrock,[-1.,-1],max_iter=200)\n",
"anim= plot_ros_cg_2d(x,rosenbrock,n_frames=20)\n",
"#anim = plot_ros_cg(x,n_frames=50,dpi=150,figsize=(10,8))\n",
"display_animation(anim)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
""
],
"metadata": {},
"output_type": "pyout",
"prompt_number": 272,
"text": [
""
]
}
],
"prompt_number": 272
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"x0 = [-1,-1]\n",
"d0 = d_rosenbrock(x0)\n",
"d0 /= mag(d0)\n",
"#linemin(x0,-d_rosenbrock(x0),rosenbrock,just_result=True)\n",
"print d0\n",
"def plot_linemin(xi,d,f):\n",
" fs = lambda a: f(xi+a*d) #scalar\n",
" t = linspace(-1,1,100).reshape(-1,1)\n",
" print shape(t),shape(fs(t))\n",
" plot(t,fs(t))\n",
" a = linemin(xi,d,f,just_result=True) \n",
" print 'a=',a\n",
" x = array([0.,a])\n",
" print x\n",
" print fs(x)\n",
" plot(x,[fs(xx) for xx in x],linestyle='',marker='o',color='g')\n",
"plot_linemin(x0,-d0,rosenbrock)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"[-0.89531628 -0.44543098]\n",
"(100, 1) (100,)\n",
"a="
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
" 1.31451932566\n",
"[ 0. 1.31451933]\n",
"141.308556937\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEACAYAAABbMHZzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8FPW9//HXYgJYlYtcNrDLaWouhIQQuQWLoqsQkFDS\nFDUaqEQi9RR6BK0Xqu3vUfAcSaieVlBTWxswastFbRMUjAE1ahWCgkdbYk2MYJPNRUtICVIMkPn9\nMWUlksRkk+xsdt/Px2MeJN+dnf3sZPnMd7/zvdgMwzAQEZGg0sfqAERExPeU/EVEgpCSv4hIEFLy\nFxEJQkr+IiJBSMlfRCQItZv8MzMzsdvtxMfHtyh/+OGHGTNmDGPHjmXFihWe8qysLKKiooiJiaGo\nqMhTvnfvXuLj44mKimL58uXd/BZERKSz2k3+ixYtorCwsEXZq6++ytatW3n//ff561//yp133glA\naWkpmzdvprS0lMLCQpYuXcrpIQRLliwhNzeX8vJyysvLzzqmiIj4VrvJf9q0aQwePLhF2a9//Wvu\nueceQkNDARg2bBgABQUFpKenExoaSnh4OJGRkZSUlFBTU0NjYyOJiYkALFy4kPz8/J54LyIi0kGd\nbvMvLy/n9ddf55JLLsHlcvHOO+8AUF1djdPp9OzndDpxu91nlTscDtxudzeELiIi3grp7BNOnjzJ\n4cOH2b17N2+//TZpaWl8/PHHPRGbiIj0kE4nf6fTybx58wCYPHkyffr04R//+AcOh4PKykrPflVV\nVTidThwOB1VVVS3KHQ5Hq8ceMCCSxsaKzoYkIhLUIiIi+Oijjzr1nE43+6SmpvLKK68AUFZWRlNT\nE0OHDiUlJYVNmzbR1NTEgQMHKC8vJzExkbCwMAYMGEBJSQmGYfDUU0+Rmpra6rEbGyvYv9/AMLT9\n/Oc/tzwGf9l0LnQudC7a3yoqOl9pbrfmn56ezmuvvcahQ4cYNWoU9913H5mZmWRmZhIfH0/fvn15\n8sknAYiNjSUtLY3Y2FhCQkLIycnBZrMBkJOTw0033cS//vUvkpOTufrqq9t8zRUr4PnnO/0+RESk\nE9pN/hs3bmy1/Kmnnmq1/N577+Xee+89q3zixIn85S9/6VBA+/fDK6/AVVd1aHcREfGC343wzc6G\nO++E5marI7GWy+WyOgS/oXPxJZ2LL+lcdI3NMAy/WczFZrPR3GwwdSosWQILF1odkYiI/7PZbHQ2\nlftd8jcMg127IC0N/vY3OO88q6MSEfFv3iR/v2v2Afj2t+HSS+HBB62OREQkMPllzR/g4EGYOBHe\nfx/aGBYgIiIEULPPaT/5CdTWwhNPWBeTiIi/C7jkf+QIjB4NL7xgfgsQEZGzBUyb/2kDBsB998Ht\nt4P/XKJERHo/v07+AJmZ5jeAZ5+1OhIRkcDh180+p736qnkR+OAD6N/fgsBERPxYwDX7nHbllTB+\nPPzyl1ZHIiISGHpFzR/g448hMRHee09dP0VEzhRwvX2+6t57oaoK/j2RqIiIEATJv7ERYmLguefg\nkkt8GJiIiB8L2Db/0y64wJz1c9kyzfopItIVvSr5AyxYAH36qOlHRKQrelWzz2lvvw3f/a7Z9XPg\nQB8EJiLixwK+zf9MixebI4DV/VNEgl1QJf9PP4W4OHjtNYiN7eHARET8WLff8M3MzMRutxMfH3/W\nY//7v/9Lnz59qK+v95RlZWURFRVFTEwMRUVFnvK9e/cSHx9PVFQUy5cv71SAbRk+HH72M1i+XPP+\niIh0VrvJf9GiRRQWFp5VXllZyY4dO/jmN7/pKSstLWXz5s2UlpZSWFjI0qVLPVeiJUuWkJubS3l5\nOeXl5a0e0xs/+pE55fNzz3XL4UREgka7yX/atGkMHjz4rPIf//jH/OIXv2hRVlBQQHp6OqGhoYSH\nhxMZGUlJSQk1NTU0NjaSmJgIwMKFC8nPz++W4ENC4JFH4Mc/hs8/75ZDiogEhU539SwoKMDpdDJu\n3LgW5dXV1TidTs/vTqcTt9t9VrnD4cDtdnch5JauuAIuuwxWr+62Q4qIBLyQzux87NgxVq9ezY4d\nOzxl/nC/+MEHYdw4uOkmiIqyOhoREf/XqeRfUVHBwYMHSUhIAKCqqoqJEydSUlKCw+GgsrLSs29V\nVRVOpxOHw0FVVVWLckc7M7OtXLnS87PL5cLlcn1tXCNHwj33wK23wosvgs3WmXclItK7FBcXU1xc\n3LWDGF/jwIEDxtixY1t9LDw83Dh06JBhGIaxf/9+IyEhwfjiiy+Mjz/+2LjooouM5uZmwzAMIzEx\n0di9e7fR3NxszJ4923jxxRdbPV4HwmlTU5NhxMUZxrPPen0IEZFeyZvc2W6bf3p6OlOnTqWsrIxR\no0axYcOGFo/bzqhix8bGkpaWRmxsLLNnzyYnJ8fzeE5ODosXLyYqKorIyEiuvvrqrl2xWhEaCjk5\ncNttcPRotx9eRCSg9NpBXm1ZuBDCwuArnZFERAJWUI3wbUtdHYwday79OHZsNwUmIuLHAn5K546w\n22HVKliyRNM+i4i0JeCSP8B//id88QXk5VkdiYiIfwq4Zp/T9u2D2bNh/34YOrRbDiki4pfU5v8V\ny5ebPX9yc7vtkCIifkfJ/yuOHDGne964EaZN67bDioj4Fd3w/YoBA2DtWrjlFvMegIiImAI6+QPM\nmweRkfDAA1ZHIiLiPwK62ee0Tz6BiRNh1y5N/CYigUfNPm345jfNid9++EOt+iUiAkGS/MHs+dPQ\noL7/IiIQJM0+p53u+/+Xv5hrAIuIBAJ19eyAu+6C6mr4/e979GVERHxGyb8DPv8c4uPh0UfNbwEi\nIr2dbvh2wHnnwW9+Y0781thodTQiItYIupr/aTfdZA4CW7fOJy8nItJj1OzTCfX1EBcHzz0HU6f6\n5CVFRHqEmn064cILzVr/zTdr6gcRCT5Bm/wBrr0WYmLgv//b6khERHyr3eSfmZmJ3W4nPj7eU3bX\nXXcxZswYEhISmDdvHv/85z89j2VlZREVFUVMTAxFRUWe8r179xIfH09UVBTLly/vgbfhHZvNXPT9\nt7+Fd9+1OhoREd9pN/kvWrSIwsLCFmUzZ85k//79vPfee0RHR5OVlQVAaWkpmzdvprS0lMLCQpYu\nXeppg1qyZAm5ubmUl5dTXl5+1jGtNGKEOelbZiacOGF1NCIivtFu8p82bRqDBw9uUZaUlESfPubT\npkyZQlVVFQAFBQWkp6cTGhpKeHg4kZGRlJSUUFNTQ2NjI4mJiQAsXLiQ/Pz8nngvXlu40LwIrFlj\ndSQiIr7RpTb/9evXk5ycDEB1dTVOp9PzmNPpxO12n1XucDhwu91dedluZ7OZff/XrjWnfhARCXQh\n3j7x/vvvp2/fvsyfP78742HlypWen10uFy6Xq1uP35ZRoyAry+z/v3s3hIb65GVFRDqtuLiY4uLi\nLh3Dq+T/xBNPsH37dl5++WVPmcPhoLKy0vN7VVUVTqcTh8PhaRo6Xe5wONo89pnJ39duvhmefRZ+\n8Qv46U8tC0NEpF1frRivWrWq08fodLNPYWEhDzzwAAUFBfTv399TnpKSwqZNm2hqauLAgQOUl5eT\nmJhIWFgYAwYMoKSkBMMweOqpp0hNTe10oL5gs8Hjj8NDD6n5R0QCW7vJPz09nalTp/Lhhx8yatQo\n1q9fz6233srRo0dJSkpi/PjxLF26FIDY2FjS0tKIjY1l9uzZ5OTkYLPZAMjJyWHx4sVERUURGRnJ\n1Vdf3fPvzEujRkF2NmRkqPePiASuoJ3eoT2GAXPmwJQp8POfWx2NiEj7NLdPN6quhosvhsJCmDDB\n6mhERNqmuX260ciR8KtfmWMANPePiAQa1fzbYRjm/D+RkRoAJiL+S80+PeCzzyAhAbZsgcsuszoa\nEZGzqdmnBwwbBo89Zvb+0cpfIhIoVPPvoJtvhnPOMWcAFRHxJ2r26UFHjpjNP+vWwdy5VkcjIvIl\nJf8e9sYbkJYG770Hw4dbHY2IiEnJ3wfuuQf274eCAnM6CBERq+mGrw+sWgVVVeYcQCIivZVq/l4o\nLYXLL4c334TRo62ORkSCnWr+PhIbC/fdBwsWQFOT1dGIiHSeav5eMgxISYGxY81FYERErKIbvj72\n6afm5G9/+AP4aMExEZGzqNnHx4YPhw0b4MYb4dAhq6MREek41fy7wR13QEUF/OlP6v4pIr6nmr9F\nVq+Gv//dnANIRKQ3UM2/m5SVwaWXwquvmjeBRUR8RTV/C0VHw4MPwvXXw7FjVkcjItK+dpN/ZmYm\ndrud+Ph4T1l9fT1JSUlER0czc+ZMGhoaPI9lZWURFRVFTEwMRUVFnvK9e/cSHx9PVFQUy5cv74G3\n4R8yMmDiRAjgtygiAaLd5L9o0SIKCwtblGVnZ5OUlERZWRnTp08nOzsbgNLSUjZv3kxpaSmFhYUs\nXbrU8zVkyZIl5ObmUl5eTnl5+VnHDCSPPgqvvQabNlkdiYhI29pN/tOmTWPw4MEtyrZu3UpGRgYA\nGRkZ5OfnA1BQUEB6ejqhoaGEh4cTGRlJSUkJNTU1NDY2kpiYCMDChQs9zwlEF1wAmzfDsmXw0UdW\nRyMi0rpOt/nX1dVht9sBsNvt1NXVAVBdXY3T6fTs53Q6cbvdZ5U7HA7cbndX4/Zr48fDypXm9M/H\nj1sdjYjI2UK68mSbzYatmzu2r1y50vOzy+XC1UuHzi5ZYvb8ufNOeOQRq6MRkUBSXFxMcXFxl47R\n6eRvt9upra0lLCyMmpoahv97VROHw0FlZaVnv6qqKpxOJw6Hg6qqqhblDoejzeOfmfx7M5sNfvc7\nmDABnnkGrrvO6ohEJFB8tWK8atWqTh+j080+KSkp5OXlAZCXl0dqaqqnfNOmTTQ1NXHgwAHKy8tJ\nTEwkLCyMAQMGUFJSgmEYPPXUU57nBLqBA83E/6MfmeMARET8htGOG264wRgxYoQRGhpqOJ1OY/36\n9cahQ4eM6dOnG1FRUUZSUpJx+PBhz/7333+/ERERYYwePdooLCz0lL/zzjvG2LFjjYiICOPWW29t\n8/W+Jpxe69e/Noxx4wzj2DGrIxGRQORN7tQIXx8wDHPu/3PPhdxcq6MRkUCjEb5+ymaD3/4W3noL\n1q+3OhoREc3t41OlpXDFFVBUZHYHFRHpDqr5+7nYWHj4Ybj2WjhjVgwREZ9Tzd8Cy5fDgQOQnw99\ndPkVkS5Szb+XeOABc+Wv1autjkREgpVq/haprobJk82BYLNnWx2NiPRmWsC9l/nzn+Gaa2DXLrjo\nIqujEZHeSs0+vcxll8HPfgbf+x58/rnV0YhIMFHN32KGAYsWmbN/btyoBeBFpPNU8++FbDZz4feK\nCnMZSBERX1DN309UVsKUKbBhA8yaZXU0ItKbqObfi40aZa4AtnAhlJdbHY2IBDolfz8ybRrcdx+k\npMCRI1ZHIyKBTM0+fmjpUrMZKD8fzjnH6mhExN+p2SdArF0LjY1w771WRyIigUrJ3w+FhsJzz8Gz\nz8KTT1odjYgEIjX7+LH9++HKK83mn6lTrY5GRPyVmn0CTFwcPPGEOQX0wYNWRyMigUTJ388lJ8OK\nFTB3rnoAiUj38Tr5Z2VlERcXR3x8PPPnz+eLL76gvr6epKQkoqOjmTlzJg1nrFiSlZVFVFQUMTEx\nFBUVdUvwwWLZMnMeoBtugJMnrY5GRAKBV23+Bw8e5KqrruKDDz6gX79+XH/99SQnJ7N//36GDh3K\n3XffzZo1azh8+DDZ2dmUlpYyf/583n77bdxuNzNmzKCsrIw+X1nJRG3+bTtxAubMgdGjzdXARERO\n81mb/4ABAwgNDeXYsWOcPHmSY8eOMXLkSLZu3UpGRgYAGRkZ5OfnA1BQUEB6ejqhoaGEh4cTGRnJ\nnj17vHnpoBUaCs88A6+8YnYFFRHpCq+S/4UXXsgdd9zBf/zHfzBy5EgGDRpEUlISdXV12O12AOx2\nO3V1dQBUV1fjdDo9z3c6nbjd7m4IP7gMHAjbtsGaNVBQYHU0ItKbhXjzpIqKCh566CEOHjzIwIED\nue6663j66adb7GOz2bC1Mz9xW4+tXLnS87PL5cLlcnkTYsAKDzcTf3IyOBwwaZLVEYmIrxUXF1Nc\nXNylY3iV/N955x2mTp3KkCFDAJg3bx67du0iLCyM2tpawsLCqKmpYfjw4QA4HA4qKys9z6+qqsLh\ncLR67DOTv7Ru8mR4/HFzDqA334RvfcvqiETEl75aMV61alWnj+FVs09MTAy7d+/mX//6F4ZhsHPn\nTmJjY5k7dy55eXkA5OXlkZqaCkBKSgqbNm2iqamJAwcOUF5eTmJiojcvLf+Wmgr33GN+A6ivtzoa\nEeltvKr5JyQksHDhQiZNmkSfPn2YMGECt9xyC42NjaSlpZGbm0t4eDhbtmwBIDY2lrS0NGJjYwkJ\nCSEnJ6fdJiHpmFtvhU8+MS8ERUXQv7/VEYlIb6HpHXq55maYPx9OnYJNmzQLqEgw0vQOQahPH8jL\ng3/8A26/3VwTWETk6yj5B4B+/eBPf4LiYvjFL6yORkR6A6/a/MX/DBoEL74Il14KdjvcdJPVEYmI\nP1PyDyAOBxQWgssFQ4fCd75jdUQi4q/U7BNgYmLMQWCLFpljAEREWqPkH4CmTIGnn4Z58+C996yO\nRkT8kZJ/gJo1y5z9MzkZPvrI6mhExN+ozT+ApaVBQwPMnAlvvGHeExARASX/gHfLLXD4MCQlwWuv\nwbBhVkckIv5AI3yDxE9/avYEeuUVc2poEQkc3uROJf8gYRjmcpDvvgsvvQTnnWd1RCLSXZT8pV3N\nzbB4sTkZ3AsvwLnnWh2RiHQHJX/5WqdOwfe/D0eOmFNC9O1rdUQi0lVK/tIhJ07A9debP2/ebK4P\nLCK9l2b1lA4JDTWnf25qggUL4ORJqyMSEV9T8g9SffvCs8+azT8LF5rNQSISPJT8g1j//ma7/2ef\n6QIgEmyU/IPcuefC1q3w6ae6AIgEEyV/OesCoHsAIoHP6+Tf0NDAtddey5gxY4iNjaWkpIT6+nqS\nkpKIjo5m5syZNDQ0ePbPysoiKiqKmJgYioqKuiV46T6nLwD/+IfZFfTECasjEpGe5HXyX758OcnJ\nyXzwwQe8//77xMTEkJ2dTVJSEmVlZUyfPp3s7GwASktL2bx5M6WlpRQWFrJ06VKam5u77U1I9zj3\nXHMtgCNHID1dFwCRQOZV8v/nP//JG2+8QWZmJgAhISEMHDiQrVu3kpGRAUBGRgb5+fkAFBQUkJ6e\nTmhoKOHh4URGRrJnz55uegvSnU7fBG5qgmuvhS++sDoiEekJXiX/AwcOMGzYMBYtWsSECRP4wQ9+\nwOeff05dXR12ux0Au91OXV0dANXV1TidTs/znU4nbre7G8KXntCvn9kNtG9fSEmBY8esjkhEuptX\nUzqfPHmSffv28cgjjzB58mRuu+02TxPPaTabDZvN1uYx2nps5cqVnp9dLhcul8ubEKWL+vaFjRvN\n5SCTk+H55+GCC6yOSkQAiouLKS4u7tIxvEr+TqcTp9PJ5MmTAbj22mvJysoiLCyM2tpawsLCqKmp\nYfjw4QA4HA4qKys9z6+qqsLRxsoiZyZ/sVZICDzxBCxZYq4HsH07XHih1VGJyFcrxqtWrer0Mbxq\n9gkLC2PUqFGUlZUBsHPnTuLi4pg7dy55eXkA5OXlkZqaCkBKSgqbNm2iqamJAwcOUF5eTmJiojcv\nLT52zjnwm9/AZZeBywW1tVZHJCLdweuVvB5++GEWLFhAU1MTERERbNiwgVOnTpGWlkZubi7h4eFs\n2bIFgNjYWNLS0oiNjSUkJIScnJx2m4TEv9hs8MAD5iIw06bBjh0QHm51VCLSFZrVUzpl3TrzQvDS\nSxAba3U0IgLe5U6t4SudsmwZDBkCV11ljgmYMsXqiETEG5reQTptwQLIzYXvfMdcF1hEeh8lf/HK\nnDlmzT8jA55+2upoRKSz1OwjXps6FV55BWbPNnsB3XGHeXNYRPyfbvhKl1VWmheA6dPhl780u4eK\niO9oDV+xTEMDpKbCsGHw5JPmJHEi4htaw1csM2iQ2f0zJARmzDCnhhYR/6XkL92mXz/4/e/h8svh\n29+G8nKrIxKRtij5S7fq0weysuCuu8zRwG+8YXVEItIaJX/pEbfcYrb9X3ONuoKK+CPd8JUetX8/\nzJ1rDgxbtcr8ZiAi3Uu9fcQvffopzJsHYWGQlwfnnWd1RCKBRb19xC8NHw4vvwznn2/eBzhjaQcR\nsYiSv/hEv36wYYO5MPyUKfDmm1ZHJBLc1OwjPvfii+acQKtXw+LFVkcj0vupzV96jQ8/hO9+15wa\n+qGHzDWDRcQ7avOXXmP0aCgpAbfbvABoeUgR31LyF8sMHAh/+pO5OPykSboPIOJLavYRv7B9Oyxa\nBD/7GfzXf2lqaJHO8Hmzz6lTpxg/fjxz584FoL6+nqSkJKKjo5k5cyYNDQ2efbOysoiKiiImJoai\noqKuvKwEoORk2LUL1q+H+fPh6FHYtmMbsxbNwnWTi1mLZrFtxzarwxQJGF1K/mvXriU2Nhbbv6tp\n2dnZJCUlUVZWxvTp08nOzgagtLSUzZs3U1paSmFhIUuXLqW5ubnr0UtAuegieOstczzAmIRtLPnl\ncorCi3jtW69RFF7E8keX6wIg0k28Tv5VVVVs376dxYsXe75ubN26lYyMDAAyMjLIz88HoKCggPT0\ndEJDQwkPDycyMpI9e/Z0Q/gSaM49Fx5/HAZGrqPykooWj1WMr+DhjQ9bFJlIYPE6+d9+++088MAD\n9Dljspa6ujrsdjsAdruduro6AKqrq3E6nZ79nE4nbrfb25eWIDB0xBetlh9vPu7jSEQCk1dr+L7w\nwgsMHz6c8ePHU1xc3Oo+NpvN0xzU1uOtWblypednl8uFy+XyJkTp5frZ+rVa3r9Pfx9HIuJ/iouL\n28y9HeVV8n/rrbfYunUr27dv5/jx4xw5coQbb7wRu91ObW0tYWFh1NTUMHz4cAAcDgeVZ0zoUlVV\nhcPhaPXYZyZ/CV7L5i+j4tEKKsZ/2fTT57kIvnXlrRiGegNJcPtqxXjVqlWdPkaXu3q+9tprPPjg\ngzz//PPcfffdDBkyhBUrVpCdnU1DQwPZ2dmUlpYyf/589uzZg9vtZsaMGXz00Udn1f7V1VPOtG3H\nNh7e+DDHm4/Tv09/5l1xK79ZNwenE3JzYehQqyMU8Q/e5E6vav6tvTDAT37yE9LS0sjNzSU8PJwt\nW7YAEBsbS1paGrGxsYSEhJCTk9Nuk5AIwJykOcxJmtOi7KZ0cyzAxRebE8UlJVkUnEgvp0Fe0ivt\n3Ak33QTXXWcuG9lftwIkiGluHwkaM2bAe++ZcwNNmgT/939WRyTSuyj5S681ZAhs3gwrVsDMmXD/\n/XDypNVRifQOavaRgFBZCZmZcOSIuVRkTIzVEYn4jpp9JGiNGgUvvQQLF8Jll8GDD8KpU1ZHJeK/\nVPOXgPPxx3DzzXD8uDlR3JgxVkck0rNU8xfBnCDu5ZfhxhvNBePvvx9OnLA6KhH/opq/BLRPPoEf\n/hBqaswJ4yZPtjoike6nmr/IV3zzm+ZCMXfeCXPnwu23m2sFiAQ7JX8JeDYbfP/78Ne/Qn09xMVB\nQYHVUYlYS80+EnRefRWWLDG7g65da347EOnN1Owj0gFXXmmODp44ESZMgNWr4YvWlw8QCVhK/hKU\n+vWD//f/4O23zbWDx42DwkKroxLxHTX7iAAvvAC33WbeD/jVr8zuoiK9hZp9RLz0ne/A/v1wySVm\nd9B77oHGRqujEuk5Sv4i/9avn5n0//IXqK6G0aPNRWM0TYQEIjX7iLRhzx644w5zsrgHH9TCMeK/\nvMmdSv4i7TAM+OMfzWmjIyNhzRpISLA6KpGW1OYv0s1sNrjmGigthTlzYNYsc+bQgwetjkyka7xK\n/pWVlVx55ZXExcUxduxY1q1bB0B9fT1JSUlER0czc+ZMGhoaPM/JysoiKiqKmJgYioqKuid6ER/p\n2xduvRXKyiA83BwjsGwZ1NVZHZmId7xq9qmtraW2tpaLL76Yo0ePMnHiRPLz89mwYQNDhw7l7rvv\nZs2aNRw+fJjs7GxKS0uZP38+b7/9Nm63mxkzZlBWVkafPi2vPWr2kd7i00/NwWFPPQW33AJ33QUX\nXmh1VBKsfNbsExYWxsUXXwzA+eefz5gxY3C73WzdupWMjAwAMjIyyM/PB6CgoID09HRCQ0MJDw8n\nMjKSPXv2ePPSIn5h+HB46CFz7eD6eoiOhp//HA4ftjoykY7pcpv/wYMHeffdd5kyZQp1dXXY7XYA\n7HY7df/+TlxdXY3T6fQ8x+l04na7u/rSIpYbNQp+8xuzZ1BlJURFwcqVugiI/+tS8j969CjXXHMN\na9eu5YILLmjxmM1mw2aztfnc9h4T6W0uushcNaykBP7+d7Nn0D33wGefWR2ZSOtCvH3iiRMnuOaa\na7jxxhtJTU0FzNp+bW0tYWFh1NTUMHz4cAAcDgeVlZWe51ZVVeFwOFo97sqVKz0/u1wuXC6XtyGK\n+FxEhHkROHgQsrPNgWI33miuJzBqlNXRSaAoLi6muLi4S8fw6oavYRhkZGQwZMgQfvWrX3nK7777\nboYMGcKKFSvIzs6moaGhxQ3fPXv2eG74fvTRR2fV/nXDVwJNdbU5V1BuLnz3u+ZFIC7O6qgk0Phs\nkNef//xnLr/8csaNG+dJ4FlZWSQmJpKWlsbf//53wsPD2bJlC4MGDQJg9erVrF+/npCQENauXcus\nWbO65Q2I9Ab19fDYY7BuHUyaZI4cdrnMcQQiXaURviJ+7vhxePJJ+OUv4RvfMJeVvP56cxyBiLeU\n/EV6ieZmePFFs0lo/35zkfkf/hD+3VlOpFM0vYNIL9GnjzldxM6dsGMHuN3mspILFpiLy6gOJD1N\nNX8RP3H4MGzYADk5cMEF5jeBBQvg/POtjkz8nZp9RAJAczO8/LJ5g/jVV+Haa+EHPzBvFOsGsbRG\nyV8kwFRXwxNPwO9+BwMGQGam+W1gyBCrIxN/ouQvEqCam6G42BxA9sILMH06ZGTA7NkQGmp1dGI1\nJX+RINDVUDI7AAAHf0lEQVTQAM88Y3YZ/fBDs6voggUwZYqahYKVkr9IkKmogD/8AZ5+Gk6ehBtu\ngPR0GDvW6sjEl5T8RYKUYcC778LGjbB5s9lb6LrrzE3TSQQ+JX8RobnZnF10yxZ49lmzq+i8eeY2\nYYKahgKRkr+ItNDcDO+8A889Zy5Ef/y4OcFcSgpccQX062d1hNIdlPxFpE2GAX/7G+Tnw/PPm9NK\nzJgByclmr6GRI62OULyl5C8iHfbZZ7B9uznH0I4d5noDs2aZ26WX6ltBb6LkLyJeOXnSXIrypZfM\nrbQUpk41xxNcdRVcfDGcc47VUUpblPxFpFs0NJiDynbuNKeYqK6Gyy837xNcfrl5MQjxeh1A6W5K\n/iLSI2przYvBG2/A66/DJ5+Yg8ouu8z8hjBlijn9hFhDyV9EfOLQIXPq6TfeMP/dtw++9S3zIpCY\nCJMnmwPNNPWEbyj5i4glTpyA994z7xuc3j75xLwATJwI48ebTUVjx8K551odbeBR8hcRv3H0qDnq\neN8+899334WyMvMbwrhxEB9vXgzi4swy3VD2nt8n/8LCQm677TZOnTrF4sWLWbFiRctglPxFAlpT\nE3zwAbz/Pvz1r19un30G0dHmamajR5tbVJS5DRpkddT+z6+T/6lTpxg9ejQ7d+7E4XAwefJkNm7c\nyJgxY74MRsnfo7i4GJfLZXUYfkHn4kuBei6OHjVnKP3b3778t7zc3Pr3h4gIuOgicwsPN7eTJ4uZ\nPdtlceT+wZvc6bPOWnv27CEyMpLw8HAAbrjhBgoKClokf/lSoP4n94bOxZcC9Vycf755b2DixJbl\nhgF1dfDxx+YMpgcPwu7dsGkTxMf7d/LftmMb6/6wji+ML+hn68ey+cuYkzTH6rA8fJb83W43o0aN\n8vzudDopKSnx1cuLSC9ks0FYmLlNndrysZUrLQmpQ7bt2MbyR5dTMb7CU1bxqPmzv1wA+vjqhWya\nSlBEgsS6P6xrkfgBKsZX8PDGhy2KqBWGj+zatcuYNWuW5/fVq1cb2dnZLfaJiIgwAG3atGnT1okt\nIiKi0znZZzd8T548yejRo3n55ZcZOXIkiYmJZ93wFRER3/BZm39ISAiPPPIIs2bN4tSpU9x8881K\n/CIiFvGrQV4iIuIbPrvh25pnnnmGuLg4zjnnHPbt29fmfoWFhcTExBAVFcWaNWt8GKHv1NfXk5SU\nRHR0NDNnzqShoaHV/cLDwxk3bhzjx48nMTHRx1H2nI78jZctW0ZUVBQJCQm8++67Po7Qd77uXBQX\nFzNw4EDGjx/P+PHj+Z//+R8LovSNzMxM7HY78fHxbe4TLJ+LrzsXnf5cdOkubhd98MEHxocffmi4\nXC5j7969re5z8uRJIyIiwjhw4IDR1NRkJCQkGKWlpT6OtOfdddddxpo1awzDMIzs7GxjxYoVre4X\nHh5uHDp0yJeh9biO/I23bdtmzJ492zAMw9i9e7cxZcoUK0LtcR05F6+++qoxd+5ciyL0rddff93Y\nt2+fMXbs2FYfD5bPhWF8/bno7OfC0pp/TEwM0dHR7e5z5uCw0NBQz+CwQLN161YyMjIAyMjIID8/\nv819jQBrqevI3/jM8zNlyhQaGhqoq6uzItwe1dHPe6B9Btoybdo0Bg8e3ObjwfK5gK8/F9C5z4Wl\nyb8jWhsc5na7LYyoZ9TV1WG32wGw2+1tfoBtNhszZsxg0qRJPP74474Mscd05G/c2j5VVVU+i9FX\nOnIubDYbb731FgkJCSQnJ1NaWurrMP1GsHwuOqKzn4se7+2TlJREbW3tWeWrV69m7ty5X/v8QBoc\n1ta5uP/++1v8brPZ2nzfb775JiNGjOCzzz4jKSmJmJgYpk2b1iPx+kpH/8ZfrdUE0mfjtI68pwkT\nJlBZWck3vvENXnzxRVJTUykrK/NBdP4pGD4XHdHZz0WPJ/8dO3Z06fkOh4PKykrP75WVlTidzq6G\nZYn2zoXdbqe2tpawsDBqamoYPnx4q/uNGDECgGHDhvG9732PPXv29Prk35G/8Vf3qaqqwuFw+CxG\nX+nIubjgggs8P8+ePZulS5dSX1/PhRde6LM4/UWwfC46orOfC79p9mmrrWrSpEmUl5dz8OBBmpqa\n2Lx5MykpKT6OruelpKSQl5cHQF5eHqmpqWftc+zYMRobGwH4/PPPKSoqarcXRG/Rkb9xSkoKTz75\nJAC7d+9m0KBBnmayQNKRc1FXV+f5/7Jnzx4MwwjKxA/B87noiE5/Lrp2/7lr/vjHPxpOp9Po37+/\nYbfbjauvvtowDMNwu91GcnKyZ7/t27cb0dHRRkREhLF69Wqrwu1Rhw4dMqZPn25ERUUZSUlJxuHD\nhw3DaHkuKioqjISEBCMhIcGIi4sLqHPR2t/4scceMx577DHPPj/60Y+MiIgIY9y4cW32DgsEX3cu\nHnnkESMuLs5ISEgwvv3tbxu7du2yMtwedcMNNxgjRowwQkNDDafTaeTm5gbt5+LrzkVnPxca5CUi\nEoT8ptlHRER8R8lfRCQIKfmLiAQhJX8RkSCk5C8iEoSU/EVEgpCSv4hIEFLyFxEJQv8fVJyX/+lV\nvDoAAAAASUVORK5CYII=\n",
"text": [
""
]
}
],
"prompt_number": 156
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"x0 = [-1.1,-1]\n",
"d0 = d_quad(x0)\n",
"d0 /= mag(d0)\n",
"#linemin(x0,-d_rosenbrock(x0),rosenbrock,just_result=True)\n",
"plot_linemin(x0,-d0,quad)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"(100, 1) (100,)\n",
"a="
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
" 1.01659614527\n",
"[ 0. 1.01659615]\n",
"1.2111027093\n"
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEACAYAAABBDJb9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VFWexvFvBVBgQNkrSHSCkcgeilVtlopQCYtBEKFB\nxDRLO62o3WLL4owSdKSjMjKA2NCKbdpxAWxAWjCySOEDgtiyiNKIIiBgEpYQZRGR5M4fR4JIQipL\n1U3dej/Pcx+Syq1bLzfXH8dzzz3HZVmWhYiIhLUouwOIiEj5qZiLiDiAirmIiAOomIuIOICKuYiI\nA6iYi4g4QEDFPD8/H4/HQ0pKCgBpaWnExMTg8XjweDxkZmYGNaSIiFxa1UB2mjFjBi1btuT48eMA\nuFwuxo0bx7hx44IaTkREAlNiy/zAgQMsX76cMWPGcO75Isuy0LNGIiKVR4nF/MEHH+SZZ54hKur8\nri6Xi1mzZpGQkMDo0aPJy8sLakgREbm0Sxbzt99+m0aNGuHxeC5oid9zzz3s2bOHrVu30rhxYx56\n6KGgBxURkUuwLmHSpElWTEyMFRsba0VHR1s1a9a0RowYccE+e/bssVq3bl3k++Pi4ixAmzZt2rSV\nYouLi7tUaS7SJYv5z/n9fuuWW26xLMuyvvnmm8LXn332WWvYsGFFH5yAD+94kydPtjtCpaFzcZ7O\nxXk6F+eVpXYGNJoFc2RcLhcA48ePZ9u2bbhcLpo2bcrcuXMDPYyIiARBwMXc6/Xi9XoBeOWVV4KV\nR0REykBPgIbIuX8IRefi53QuztO5KB/XT/0zwTm4y6Xx6CIipVSW2qmWuYiIA6iYi4g4gIq5iIgD\nqJiLiDiAirmIiAOomIuIOICKuYiIA6iYi4g4gIq5iIgDqJiLiDiAirmIiAOomIuIOICKuYiIAwS9\nmG/bFuxPEBGRgIp5fn4+Ho+HlJQUAHJzc/H5fMTHx5OUlEReXl6x773rLjhzpmLCiohI0QIq5jNm\nzKBly5aFy8alp6fj8/nYtWsXPXv2JD09vdj3xsbC449XSFYRESlGicX8wIEDLF++nDFjxhROlr50\n6VJSU1MBSE1NZcmSJcW+f+5cePFF2LixghKLiMhFSizmDz74IM888wxRUed3zcnJwe12A+B2u8nJ\nySn2/dHRMHu26W45daoCEouIyEUuuaDz22+/TaNGjfB4PPj9/iL3cblchd0vRUlLSwOgenUYPtzL\n4sXesmYVEXEkv99fbI0N1CXXAH3kkUd45ZVXqFq1KqdPn+a7777jtttu46OPPsLv9xMdHU1WVhaJ\niYns3Lnz4oP/bB27vDxo0wbmzYOkpHJlFhFxtLKsARrwgs5r165l2rRp/OMf/2D8+PHUr1+fCRMm\nkJ6eTl5eXpE3QX8ZaOVKGDUKPvkE6tYtVU4RkYgR9AWdz3WnTJw4kZUrVxIfH897773HxIkTA3q/\nzwcDB8LYsaXKKCIiJQi4ZV6mgxfxr8upU9ChA0yeDEOHBuuTRUTCV1C7WcqiuED//Cf07QubN0NM\nTLA+XUQkPAW9m6WidOwI998PI0dCQYEdCUREnMW2ibYmTYITJ2DmTLsSiIg4hy3dLOfs3g033ABr\n1kDr1sFKISISXsKmm+WcuDhIT4c774QffrAziYhIeLO1ZQ5gWXDbbXDddfDMM8FKIiISPsJmNMsv\nHTkCCQnwt79Bz57BSiMiEh7CrpvlnAYN4K9/hd/8BnJz7U4jIhJ+KkXL/JwHH4T9+2HhQrjE3F0i\nIo4Wti3zc/70J/jiC3jpJbuTiIiEl0rVMgf47DPwemHdOrj++uDkEhGpzMK+ZQ7QqhVMmQLDh2vt\nUBGRQFW6ljmY4Yq33grNm8PTTwchmIhIJRa2QxOLcvgweDxmlIvPV8HBREQqMUd0s5zTsCFkZJjh\niocO2Z1GRKRyq7Qt83MmToTt2+HttzVcUUQiQ1Ba5qdPn6ZLly60a9eOli1bMmnSJMAs1BwTE4PH\n48Hj8ZCZmVm21CV44gnzhOiMGUE5vIiIIwTUMj916hQ1a9bk7NmzdO3alWnTprF69Wpq167NuHHj\nij94BbTMAb76Crp0gXffhfbty304EZFKLWh95jVr1gTgzJkz5OfnU/en1ZiD2ENzgWuvhVmzzDJz\nx4+H5CNFRMJKQMW8oKCAdu3a4Xa7SUxMpFWrVgDMmjWLhIQERo8eTV5eXlCDDh0K3bvDffcF9WNE\nRMJSqW6AfvvttyQnJ5Oenk7Lli1p2LAhAI8++ihZWVnMmzfvwoO7XEyePLnwe6/Xi9frLXPYkyeh\nUydzU/Suu8p8GBGRSsXv9+P3+wu/nzJlSvDHmT/xxBPUqFGDP/7xj4Wv7d27l5SUFLZv337hwSuo\nz/zntm+Hm2/W4/4i4lxB6TM/cuRIYRfK999/z8qVK/F4PGRnZxfus3jxYtq0aVPKuGXTpg08+SQM\nGQKnT4fkI0VEKr0SW+bbt28nNTWVgoICCgoKGDFiBA8//DB33XUXW7duxeVy0bRpU+bOnYvb7b7w\n4EFomYN53H/YMKhXD55/vsIPLyJiK0c9zl+Sb7+FDh1g6lTTShcRcYqIKuYAmzdD797wwQdmDVER\nESdw1NwsgWjfHh57TP3nIiJh3TIH038+ZAg0agSzZwf1o0REQiLiWuZgJt968UXzqP/8+XanERGx\nR9i3zM/ZvBmSk2H9eoiPD8lHiogERUS2zM9p397MsDh4MHz/vd1pRERCyzEtczD958OHQ40a8IuZ\nBUREwkZEt8zB9J//5S+wYQO89JLdaUREQsdRLfNzduyAHj1g1SpISAj5x4uIlEvEt8zPadnSrEx0\n++3mSVEREadzZMv8nLFj4ZtvYNEirR8qIuFDLfNfePZZU8ynTbM7iYhIcDm6ZQ7w9dfQuTO88QaU\nY10MEZGQUcu8CNdcA6+8AnfcAQcP2p1GRCQ4HF/MAXw+s3bo4MFw5ozdaUREKp7ju1nOKSiAgQNN\nS33WLLvTiIgUr8K7WU6fPk2XLl1o164dLVu2ZNKkSQDk5ubi8/mIj48nKSmpcFm5yiwqCjIyzIRc\nr7xidxoRkYpVYsv81KlT1KxZk7Nnz9K1a1emTZvG0qVLadCgAePHj+epp57i2LFjpKenX3zwStQy\nP+fTTyExEVasAI/H7jQiIhcLyg3QmjVrAnDmzBny8/OpW7cuS5cuJTU1FYDU1FSWLFlShrj2aN3a\nzHt+221w9KjdaUREKkaJxbygoIB27drhdrtJTEykVatW5OTkFC7e7Ha7ycnJCXrQijRkiLkZOmwY\nnD1rdxoRkfKrWtIOUVFRbN26lW+//Zbk5GTWrFlzwc9dLheuSzxemZaWVvi11+vFW0kGe0+dCn36\nwCOPwNNP251GRCKZ3+/H7/eX6xilGs3yxBNPUKNGDV588UX8fj/R0dFkZWWRmJjIzp07Lz54Jewz\n/7mjR6FTJ1PYhw61O42IiFHhfeZHjhwpHKny/fffs3LlSjweD/379ycjIwOAjIwMBgwYUMbI9qpf\nHxYvhvvvh23b7E4jIlJ2l2yZb9++ndTUVAoKCigoKGDEiBE8/PDD5ObmMmTIEL7++mtiY2NZsGAB\nderUufjglbxlfs4bb5julk2boEEDu9OISKQrS+2MmIeGSjJxoinm774L1arZnUZEIpmKeTnk50NK\nClx3HcycaXcaEYlkmmirHKpUgddeMy1zrR8qIuFGLfNf2LkTunc3N0Z/9Su704hIJFLLvAI0b27m\ncBk82MyFLiISDlTMi9CnD4wbB7feCidP2p1GRKRk6mYphmXBb35jivmCBWbWRRGRUFA3SwVyuWDu\nXLOG6JQpdqcREbk0FfNLqF7d3AjNyID58+1OIyJSPHWzBGDbNujVC5YtM4tDi4gEk7pZgiQhwYw9\nHzgQ9u+3O42IyMVKnAJXjP79Ydcu85TounVQq5bdiUREzlM3SylYFvz2t3D4MCxaZJ4aFRGpaOpm\nCTKXC55/Ho4fh/Hj7U4jInKeinkpXXYZ/P3v5mbonDl2pxERMdRnXgZ165pi3rUrNG0Kycl2JxKR\nSFdiy3z//v2FCzm3bt2amT/ND5uWlkZMTAwejwePx0NmZmbQw1YmcXHw5pswYgR88ondaUQk0pV4\nAzQ7O5vs7GzatWvHiRMn6NChA0uWLGHBggXUrl2bcePGFX9wh90ALcr8+ab/fMMGuOoqu9OIiBOU\npXaW2M0SHR1NdHQ0ALVq1aJFixYcPHgQwPGFOhC//jV89RXccgu8/76GLIqIPUp1A3Tv3r1s2bKF\nG264AYBZs2aRkJDA6NGjCxd+jkQTJ0KHDqawnz1rdxoRiUQBF/MTJ05w++23M2PGDGrVqsU999zD\nnj172Lp1K40bN+ahhx4KZs5K7dyQxfx8GDvWjEcXEQmlgEaz/PjjjwwaNIg777yTAQMGANCoUaPC\nn48ZM4aUlJQi35uWllb4tdfrxev1lj1tJVatGixcaFYpSk+HSZPsTiQi4cLv9+P3+8t1jBJvgFqW\nRWpqKvXr12f69OmFr2dlZdG4cWMApk+fzkcffcRrr7124cEj4AboL33zDdx0EzzxhBnpIiJSWmWp\nnSUW83Xr1tG9e3fatm2Ly+UCYOrUqbz++uts3boVl8tF06ZNmTt3Lm63u9yBnGDHDkhMhFdfNbMt\nioiURlCKeXlEajEHM7Ll9tth5Uoz66KISKA0N0sl0r07zJ4N/frB3r12pxERp9Pj/EE0eDBkZ0Pv\n3rB+PdSvb3ciEXEqdbOEwKRJ4PfDqlXwb/9mdxoRqezUZ15JWRaMHAmHDsFbb5lhjCIixVGfeSXl\ncsELL0BUFIweDQUFdicSEadRMQ+RatVgwQL48kszMZf+h0VEKpKKeQjVrAlvvw2ZmfD003anEREn\n0WiWEKtXD9591yxsUb8+jBljdyIRcQIVcxs0aQIrVkCPHmbVokGD7E4kIuFOxdwmzZqZpeeSk+GK\nK8DnszuRiIQz9ZnbyOOBRYvgjjvMSkUiImWlYm6zrl3hlVdgwADYts3uNCISrlTMK4HeveG556BP\nH/j8c7vTiEg4Up95JTF4MJw4YfrO338fYmPtTiQi4UTFvBIZOdIU9F69TEG/6iq7E4lIuFAxr2Tu\nvx9OnoSePWHtWvjZ6nwiIsUqsc98//79JCYm0qpVK1q3bs3MmTMByM3NxefzER8fT1JSEnl5eUEP\nGykmTjTdLklJkJtrdxoRCQclzpqYnZ1NdnY27dq148SJE3To0IElS5bw17/+lQYNGjB+/Hieeuop\njh07Rnp6+oUH16yJZWZZ8PDDpnW+ahVceaXdiUQkVIIya2J0dDTt2rUDoFatWrRo0YKDBw+ydOlS\nUlNTAUhNTWXJkiVliCzFcbngmWfghhvMaJfjx+1OJCKVWanmM9+7dy89evTg008/5ZprruHYsWMA\nWJZFvXr1Cr8vPLha5uVWUAC/+x3s3AnvvKPFLUQiQVDnMz9x4gSDBg1ixowZ1K5d+6IPdrlcpfpg\nCUxUFMyZA3FxcMstcOqU3YlEpDIKaDTLjz/+yKBBgxgxYgQDBgwAwO12k52dTXR0NFlZWTQqZthF\nWlpa4dderxev11vu0JEmKgpefNEMXUxJgX/8w0ynKyLO4Pf78fv95TpGid0slmWRmppK/fr1mT59\neuHr48ePp379+kyYMIH09HTy8vJ0AzTI8vMhNRVycmDpUqhRw+5EIhIMQVkDdN26dXTv3p22bdsW\ndqX86U9/onPnzgwZMoSvv/6a2NhYFixYQJ06dcodSC4tPx/uuuv8eqJqoYs4jxZ0jhA/b6GroIs4\njxZ0jhBVqkBGBkRHmz70kyftTiQidlMxD1NVqsDLL8PVV0O/fmZOFxGJXCrmYaxKFXjpJbjuOjN9\nrh4sEolcKuZhLioK/vIXaNXKzOWiKXJEIpOKuQNERcGf/wxdupjZFo8etTuRiISairlDuFwwfbpp\nnXu9kJ1tdyIRCSXNZ+4gLhdMnWoeJureHVavNjdIRcT5VMwdxuWCxx6DWrVMQV+50twgFRFnUzF3\nqHHjTEH3eiEzE1q3tjuRiASTirmD3X031K5t1hRduhQ6d7Y7kYgEi26AOtywYWbGxVtugffeszuN\niASLinkEuOUWWLgQhg6FxYvtTiMiwaBulgjRo4fpO+/XD44dg1Gj7E4kIhVJxTyCtG9vFohOSjJT\n6E6YYEa/iEj40xS4EejgQbNIdK9e8D//Y54gFZHKQ/OZS8COHYP+/c1DRS+/DJddZnciETknKPOZ\njxo1CrfbTZs2bQpfS0tLIyYmBo/Hg8fjITMzs/RpxVZ168KKFfD999C3L3z3nd2JRKQ8SizmI0eO\nvKhYu1wuxo0bx5YtW9iyZQu9e/cOWkAJnho14M03IT7ePC2alWV3IhEpqxKLebdu3ahbt+5Fr6v7\nxBmqVIHZs2HwYLjxRtixw+5EIlIWZb71NWvWLBISEhg9ejR5mkQ7rLlc8J//CY8/DomJZsSLiISX\ngG6A7t27l5SUFLZv3w7AoUOHaNiwIQCPPvooWVlZzJs37+KDu1xMnjy58Huv14vX662g6BIMq1eb\np0ZnzDB/ikjw+f1+/H5/4fdTpkwJzmiWXxbzQH+m0Szhaft289To3XfDI49oLLpIqAVlNEtRsn52\np2zx4sUXjHSR8NemDWzYAIsWwejRcOaM3YlEpCQltsyHDRvG2rVrOXLkCG63mylTpuD3+9m6dSsu\nl4umTZsyd+5c3G73xQdXyzysnTgBd9xh/nzzTahXz+5EIpFBDw1JhcvPh/Hj4e23zdasmd2JRJwv\nZN0sEjmqVDGP/D/0EHTrBmvW2J1IRIqiYi4BuftueO01M43u3Ll2pxGRX1I3i5TKF19ASgr4fPDs\ns1Ctmt2JRJxH3SwSdM2awcaN8OWXZubFo0ftTiQioGIuZVCnjrkZ2r69WVe0iEcMRCTEVMylTKpU\ngWeegSlT4Oab4e9/tzuRSGRTn7mU28cfw223wZ13mvldqlSxO5FIeNM4c7HNoUPw619D9erw6qt6\nwEikPHQDVGzTqBGsXAktW0KnTrB1q92JRCKLirlUmKpVzQNGTz5phi6+/LLdiUQih7pZJCg++wwG\nDTIrGM2cabpfRCQw6jOXSuX4cTPr4pdfwsKFsPOrZcx8bSY/WD9wuetyHrjjAfr5+tkdU6TSKUvt\nrBqkLCLUrg3z58Nzz0H7G5dRo+3vyem2u/Dnu2ebr1XQRcpPLXMJiRtuS+bDhBUXvZ68L5nMlzKL\neIdI5NJoFqm0ql/xQ5Gvny44HeIkIs6kYi4hcbnr8iJfrx6lO6MiFaHEYj5q1CjcbvcFS8Pl5ubi\n8/mIj48nKSmJvLy8oIaU8PfAHQ8QtyXugteqLYkjP+d+TpywKZSIg5RYzEeOHElm5oV9munp6fh8\nPnbt2kXPnj1JT08PWkBxhn6+fswYO4Pkfcn02NOD5H3JvDZlBle7+9G+Pfzzn3YnFAlvAd0A3bt3\nLykpKWz/aXq85s2bs3btWtxuN9nZ2Xi9Xnbu3HnxwXUDVAKwYAHcd59ZzeiPf9TcLiIhuwGak5NT\nuICz2+0mJyenLIcRAWDIEPjoI1i2zMzAuG+f3YlEwk+5x5m7XC5cLlexP09LSyv82uv14vV6y/uR\n4kD//u9mfdFp06BjR7OK0Z13wiUuLRHH8Pv9+P3+ch2jzN0sfr+f6OhosrKySExMVDeLVJgtW2DE\nCGjeHP78Z2jY0O5EIqEVsm6W/v37k5GRAUBGRgYDBgwoy2FEiuTxmBui114LCQmwdKndiUQqvxJb\n5sOGDWPt2rUcOXIEt9vN448/zq233sqQIUP4+uuviY2NZcGCBdSpU+fig6tlLuW0bh2MHAk33ggz\nZkDdunYnEgk+TbQljnTyJEyaBIsWwfPPQ//+dicSCS4Vc3G0tWvNLIxduphWeoMGdicSCQ7NzSKO\n1qMHfPIJREdDmzbw+uugtoKIoZa5hKUPP4Tf/hauvtqMeLnmGrsTiVQctcwlYnTpYka83HQTtG8P\n06fD2bN2pxKxj1rmEvZ27YJ77oG8PJgzxywoLRLO1DKXiBQfD6tWwR/+YEa63HuvKewikUTFXBzB\n5TJPje7YYW6KtmgBL78MBQV2JxMJDXWziCN99BGMHQtVq8Ls2eapUpFwoW4WkZ906gQbN8KoUdCn\nD/zHf8Dhw3anEgkeFXNxrKgoGDMGdu6EmjWhZUsz6uXMGbuTiVQ8FXNxvDp1TBF//31zo7R1a3jr\nLT1wJM6iPnOJOO++a1Y1atDg/PzpIpWJ+sxFApCcDFu3wvDhZijj8OHw1Vd2pxIpHxVziUhVq5rp\nAHbtguuvNzdMf/97OHTI7mQiZaNiLhGtVi147DH417/Oj0//r//SQ0cSflTMRYBGjWDmTPj4Y/jm\nG2jWDJ58Eo4ftzuZSGDKdQM0NjaWK664gipVqlCtWjU2bdp04cF1A1TC1Oefw5QpsHq1uVl6772m\nFS8SCiFfnKJp06Z8/PHH1KtXr8ICiVQmn30Gjz8Ofr+Z+2XsWLjiCrtTidPZMppFxVqcrFUrmD8f\n1qyBTz+FuDiYPBmOHrU7mciFylXMXS4XvXr1omPHjrzwwgsVlUmk0mnZEl59FT74wPSpx8fDuHGw\nf7/dyUSMquV58/r162ncuDGHDx/G5/PRvHlzunXrdsE+aWlphV97vV68Xm95PlLEVs2awQsvmNb5\n//4vtGsHKSmmsLdta3c6CVd+vx+/31+uY1TYE6BTpkyhVq1aPPTQQ+cPrj5zcbhjx8yCGLNmmWkC\nxo2DpCQzL4xIWYW0z/zUqVMc/2nc1smTJ1mxYgVt2rQp6+FEwlLdujBpEuzZY54knTjR9LM//zyc\nOGF3OokkZW6Z79mzh4EDBwJw9uxZhg8fzqRJky48uFrmEmEsC9auhRkzzMRed95phjVef73dySSc\nhHxoYokHVzGXCLZvH8ydCy++aPrTf/c7uPVWqFbN7mRS2amYi1RCP/wAixaZwr5zJ6SmwujRZkSM\nSFE0a6JIJXT55TBsmHnwyO83XTHdu0PXrjBvHnz3nd0JxQnUMhexwY8/wjvvmEWn33sP+vY1/es+\nn7phRN0sImHp8GFYsAD+7//MvOq3325a8jfdpCGOkUrFXCTM7d4Nb7wBr78O334LgwbB4MFw440q\n7JFExVzEQXbsgDffhIUL4cgRMxJm4EDwek0/vDiXirmIQ33xBSxZAosXm5kce/Uy0wj07g3R0Xan\nk4qmYi4SAQ4dguXLYdkyWLUKrr3WrGualGS6Y9RqD38q5iIR5scfYeNGePddWLnSLH93001w882Q\nmAgej1nvVMKLirlIhMvNNdMJrFljtn37TGu9a1f41a+gc2etmBQOVMxF5AJHj8K6dbB+vdm2bjXT\n+HbuDJ06QYcOZrbHyy6zO6n8nIq5iFzSDz/AJ5/Ahx/Cpk2webMZ296ihZk/JiEB2rQxMz+63eBy\n2Z04MqmYi0ipnToF27fDtm1m+/RTM2LGssxsj82bm3lk4uLguuvMDdcrr6z4HMtWLmPmazP5wfqB\ny12X88AdD9DP16/iPygMlKV26taISISrWRO6dDHbOZZlRs18/vn57cMPzUNNe/ZAlSoQGwtXXw0x\nMdCkiRkiGR0NjRpBgwbQsKHpnw+kdb9s5TJ+P/v37PbsLnxt92zzdaQW9NJSy1xESsWyzI3Wffvg\nwAGzHTwI2dlmy8kxffWHD5tunTp1zFa7ttlq1YLq1aFGDdNXf9llsHxbMvt7r7jos5L3JZP5UqYN\nf0t7hbxlnpmZyR/+8Afy8/MZM2YMEyZMKM/hRCQMuFxQv77Z2re/9L5nzkBentmOHzerL504Ad9/\nb7YzZ8zwylVf/lDk+08XnA7C38CZylzM8/Pzue+++1i1ahVNmjShU6dO9O/fnxYtWlRkPsfw+/1a\nzPonOhfnOf1cXHaZ6XZp1OjS+y3+8HJ27wGaXvh69ajqQcvmNGWeumfTpk1cd911xMbGUq1aNYYO\nHcpbb71VkdkcpbwrbzuJzsV5OhfGA3c8QN2NdS94LW5zHPcPu9+mROGnzC3zgwcPcvXVVxd+HxMT\nw4cfflghoUQksvTz9aN3p97k7svldMFpqkdV5/777tfNz1IoczF3aQCqiFSg+Lh40tLS7I4Rvqwy\n2rBhg5WcnFz4/dSpU6309PQL9omLi7MAbdq0adNWii0uLq7UNbnMQxPPnj3L9ddfz+rVq7nqqqvo\n3Lkzr7/+um6AiojYoMzdLFWrVuW5554jOTmZ/Px8Ro8erUIuImKToD40JCIioVGhqwouXLiQVq1a\nUaVKFTZv3lzsfpmZmTRv3pxmzZrx1FNPVWSESiM3Nxefz0d8fDxJSUnk5eUVuV9sbCxt27bF4/HQ\nuXPnEKcMnkB+xw888ADNmjUjISGBLVu2hDhh6JR0Lvx+P1deeSUejwePx8N///d/25AyNEaNGoXb\n7aZNmzbF7hMp10VJ56LU10VZb4AW5V//+pf1+eefW16v1/r444+L3Ofs2bNWXFyctWfPHuvMmTNW\nQkKCtWPHjoqMUSk8/PDD1lNPPWVZlmWlp6dbEyZMKHK/2NhY6+jRo6GMFnSB/I6XLVtm9enTx7Is\ny9q4caPVpUsXO6IGXSDnYs2aNVZKSopNCUPr/ffftzZv3my1bt26yJ9HynVhWSWfi9JeFxXaMm/e\nvDnx8fGX3CdSHjZaunQpqampAKSmprJkyZJi97Uc1tMVyO/45+enS5cu5OXlkZOTY0fcoAr0enfa\nNVCcbt26Ubdu3WJ/HinXBZR8LqB010WFFvNAFPWw0cGDB0MdI+hycnJwu90AuN3uYi9Il8tFr169\n6NixIy+88EIoIwZNIL/jovY5cOBAyDKGSiDnwuVy8cEHH5CQkEDfvn3ZsWNHqGNWGpFyXQSitNdF\nqUez+Hw+srOzL3p96tSppKSkBBTQKYo7F08++eQF37tcrmL/3uvXr6dx48YcPnwYn89H8+bN6dat\nW1Dyhkqgv+NftjqcdG2cE8jfqX379uzfv5+aNWvyzjvvMGDAAHbt2hWCdJVTJFwXgSjtdVHqYr5y\n5cpyBWzZKW7QAAAB10lEQVTSpAn79+8v/H7//v3ExMSU65h2udS5cLvdZGdnEx0dTVZWFo2KmWmo\ncePGADRs2JCBAweyadOmsC/mgfyOf7nPgQMHaNKkScgyhkog56J27dqFX/fp04d7772X3Nxc6tWr\nF7KclUWkXBeBKO11EbRuluL6ejp27MgXX3zB3r17OXPmDPPnz6d///7BimGb/v37k5GRAUBGRgYD\nBgy4aJ9Tp05x/PhxAE6ePMmKFSsueZc/XATyO+7fvz9/+9vfANi4cSN16tQp7JZykkDORU5OTuF/\nL5s2bcKyrIgs5BA510UgSn1dlO9+7IUWLVpkxcTEWNWrV7fcbrfVu3dvy7Is6+DBg1bfvn0L91u+\nfLkVHx9vxcXFWVOnTq3ICJXG0aNHrZ49e1rNmjWzfD6fdezYMcuyLjwXu3fvthISEqyEhASrVatW\njjoXRf2O58yZY82ZM6dwn7Fjx1pxcXFW27Ztix395AQlnYvnnnvOatWqlZWQkGDdeOON1oYNG+yM\nG1RDhw61GjdubFWrVs2KiYmx5s2bF7HXRUnnorTXhR4aEhFxgJCPZhERkYqnYi4i4gAq5iIiDqBi\nLiLiACrmIiIOoGIuIuIAKuYiIg6gYi4i4gD/D5EXl1b8RQCFAAAAAElFTkSuQmCC\n",
"text": [
""
]
}
],
"prompt_number": 172
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"N=100\n",
"x = linspace(-1.5,2,N); y = linspace(-1.5,2,N);\n",
"X,Y = meshgrid(x,y)\n",
"plt.contour(X, Y, quad(dstack((X,Y))), levels=[.3,1.,2.,10.,50])\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 171,
"text": [
""
]
},
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4m9WdL/CvbMuLLMmyZVte5C3el8RxEjCBhDiAgSRg\nQuGWcEsnw9LJ0IXhPtOBTtu5JA+Qwtx2uDBMSYGWC7e9hbKnJKSBggkkcRyclCxe4iWOZdmWLUuW\nrX0794+jxU6cxLYcS7J+n+c5z3klv5aOX7/v9xy9Oq8kYIwxEEIIiRhRwW4AIYSQhUXBTwghEYaC\nnxBCIgwFPyGERBgKfkIIiTAU/IQQEmECCn6VSoX169ejsrISVVVVeOGFF6Zd75FHHkFxcTGqq6tx\n/PjxQJ6SEEJIgGIC+WWhUIjnnnsOy5cvh9FoxMqVK1FfX4/y8nLfOnv37kVXVxc6Oztx5MgRPPzw\nw2hqagq44YQQQuYmoBF/RkYGli9fDgAQi8UoLy/HwMDAlHV2796NrVu3AgBqa2sxNjYGjUYTyNMS\nQggJwLyd4+/t7cXx48dRW1s75X61Wo2cnBzfbaVSif7+/vl6WkIIIbM0L8FvNBpx99134/nnn4dY\nLL7g5+d/KoRAIJiPpyWEEDIHAZ3jBwCHw4G77roL9913HzZv3nzBz7Ozs6FSqXy3+/v7kZ2dfcF6\nRUVF6O7uDrQ5hBASUQoLC9HV1TWr3wloxM8Yw4MPPoiKigo8+uij067T0NCAN954AwDQ1NQEmUwG\nhUJxwXrd3d1gjIVteeKJJ4LehkhsO7U/+IXaH9wylwFzQCP+gwcP4ve//z2WLVuGmpoaAMDOnTvR\n19cHANi2bRs2btyIvXv3oqioCImJiXjttdcCeUpCCCEBCij416xZA7fbfdn1XnzxxUCehhBCyDyi\nK3fnSV1dXbCbMGfh3HaA2h9s1P7wI2CMhcQXsQgEAoRIUwghJGzMJTtpxE8IIRGGgp8QQiIMBT8h\nhEQYCn5CCIkwFPyEEBJhKPgJISTCUPATQkiEoeAnhJAIQ8FPCCERhoKfEEIiDAU/IYREGAp+QgiJ\nMBT8hBASYSj4CSEkwlDwE0JIhKHgJ4SQCEPBTwghEYaCnxBCIgwFPyGERBgKfkIIiTAU/IQQEmEC\nDv4HHngACoUCS5cunfbnjY2NSEpKQk1NDWpqavDUU08F+pSEEEICEBPoA9x///340Y9+hL/7u7+7\n6Drr1q3D7t27A30qQggh8yDgEf/atWuRnJx8yXUYY4E+DSGEkHlyxc/xCwQCHDp0CNXV1di4cSNa\nW1uv9FMSQgi5hIBP9VzOihUroFKpIBKJ8PHHH2Pz5s04c+bMtOtu377dt1xXV4e6uror3TxCCAkr\njY2NaGxsDOgxBGwezsP09vbi9ttvx8mTJy+7bkFBAVpaWpCSkjK1IQIBnRIihJBZmkt2XvFTPRqN\nxteo5uZmMMYuCH1CCCELJ+BTPffeey+++OILaLVa5OTkYMeOHXA4HACAbdu24Z133sFLL72EmJgY\niEQivPnmmwE3mhBCyNzNy6me+UCnegghZPZC8lQPIYSQ0ELBTwghEYaCnxBCIgwFPyGERBgKfkII\niTAU/IQQEmEo+AkhJMJQ8BNCSISh4CeEkAhDwU8IIRGGgp8QQiIMBT8hhEQYCn5CCIkwFPyEEBJh\nKPgJISTCUPATQkiEoeAnhJAIQ8FPCCERhoKfEEIiDAU/IYREGAp+QgiJMBT8hBASYSj4CSEkwgQc\n/A888AAUCgWWLl160XUeeeQRFBcXo7q6GsePHw/0KQkhhAQg4OC///77sW/fvov+fO/evejq6kJn\nZydefvllPPzww4E+JSGEkAAEHPxr165FcnLyRX++e/dubN26FQBQW1uLsbExaDSaQJ+WEELIHMVc\n6SdQq9XIycnx3VYqlejv74dCobjST01CGGOAyw24XLx2e267Pcvn18xbJv3+dAQCQDB5WQBETaqj\novx1tKecv0zIYnfFgx8A2HlHqUAgmHa97du3+5br6upQV1d3BVsVGRgDbA7AbAXMNl4sk2qLHbBe\nrDj471rtvLZ7btscgN3Jb3trh8t/2+EEnG5P7eLFe5/T5Q/7i4WvQOC577zg9hXP33b+bnR+x+At\n7kmdhtvt72Bc7kmdjqdNAgEQEz1NieK1MIaXWE8tjAZihfx2bAwQF+uphf46PtZfT1cSvHUcXxbF\n8WVRHCCK5/dFRy/gTkNCWmNjIxobGwN6jCse/NnZ2VCpVL7b/f39yM7OnnbdycEfiRjjAT1u5mXC\nW1uAcROvJzy3jd5iBUxWvmzyhLvJyh/He1sYAyTGTw0Vb7B4gydh0vLksJIk8NpXzgu2WKE/CH1h\n6AlEb2gKYy4MUm/AhxpvB+B08c7M5fIvOzyd2vnL3o7P5gBsdt75eTtHb6dp83SkxnH/ssXGf26x\n+5cnd8hmG/8/Wux8eybG845AFMeXE+MBcYKnnrQsEfFlqYj//yQiviwVAdJE//3UmYSn8wfFO3bs\nmPVjXPHgb2howIsvvogtW7agqakJMpls0Z7mYYwHs94I6Cd4PeYpvmWT/z6DyVPMvB438SCVJABJ\n3gNU5K+9B7E4HshTeA7yhAtDQBTvv49Gi7MTHc1LrDDYLfGb/KrNNKl4O3vvstHqHxAMj/kHCeOm\nSQMIM79ttPIOxLufJSUCMrGnTvTflomB5El1ssRfx9B+FbYCDv57770XX3zxBbRaLXJycrBjxw44\nHA4AwLZt27Bx40bs3bsXRUVFSExMxGuvvRZwoxeC1Q5oDZ4yzuvRcUA3cWHtLfoJPlJOkUw9QGST\nDqLCTP8Bdn6RikIrcEhoEAj8p4VSpPPzmG437yAMnk7BOwjxDkjGTHx/7tVMHbxMHtAkxvP9O2VS\nkUv9tbekSoHUJF6SEkPzlV6kEbDzT8AHiUAguOC9gPnkdAEjY4BGz0dD55cRg6f2LNscfIdNk/Fa\nPk1JOW+nl4kpuElkcLv5KwjdeYMf74Do/DJi4MVq58dOWhKQLuPFt5wMKGRTlxMTgv2Xhr65ZGdY\nBz9jfEcb0gGDOl4P6S+sNXo+QkmRAIpkXi7Y6WQ85NOSeC0V0ciEkPlms/NX0CPnDbjOLxrPcRsd\nxY/XjGQgI8W/nCn335fpuV+4IFNVQs+iCn6bHVCPAv0jgFrrXx4YnVoS4qbuAJkp/h3Eu0Mokvmo\nnc51ExI+vO+ZafSAZswzkNNNWtYDg6N80Dc8xk87ZcuBLDmQneopnmVlGqBM5a/KF9uALuyD//af\nM/RrgX4tH6Fnyaf+47zLWZPqhLhgt5wQEmwuFw//gVH/QNG73D/CM0Wt5VOKlam85KTxXMnxlNx0\nICedvw8RTsI++N/7kiEnjYe6IpkupiGEzK8Js78jUI3wZZWn9A3zEh3FO4E8BZCv8Nf5GUBBBn+P\nIpReNYR98IdIUwghEYoxfrbhnAY4N8zrXg2vzw7xZYeTdwQFno5gSaa/FGTw6dQLiYKfEEKuMIMJ\nODvIO4KzQ0DPIC/ejsE7bXtJJlCYBRRn+4tMPP/toeAnhJAgcrv5ewvdA0D3IK871f6SEAeUZAMl\nSn8pVQJF2fw6jbmg4CeEkBDFGJ+N1KkGzvQDZzx1h4q/WsiU806gVAmU5QJlOUBFHp9qfqn3FCj4\nCSEkDDld/PRRh6cjaPeU1nP85xV5vFR6Sz6fAMM/uJCCnxBCFg3G+MVurX3A6V7eEZzqBU57OoQn\nvgs8cicFPyGELHqM8esW3G4gK5WCnxBCIspcspMukSKEkAhDwU8IIRGGgp8QQiIMBT8hhEQYCn5C\nCIkwFPyEEBJhKPgJISTCUPATQkiEoeAnhJAIQ8FPCCERhoKfEEIiTMDBv2/fPpSVlaG4uBjPPvvs\nBT9vbGxEUlISampqUFNTg6eeeirQpySEEBKAmEB+2eVy4Yc//CE+/fRTZGdn46qrrkJDQwPKy8un\nrLdu3Trs3r07oIYSQgiZHwGN+Jubm1FUVIT8/HwIhUJs2bIFH3744QXr0aduEkJI6Ago+NVqNXJy\ncny3lUol1Gr1lHUEAgEOHTqE6upqbNy4Ea2trYE8JSGEkAAFdKpHcKkvgvRYsWIFVCoVRCIRPv74\nY2zevBlnzpyZdt3t27f7luvq6lBXVxdI8wghZNFpbGxEY2NjQI8R0BexNDU1Yfv27di3bx8A4Be/\n+AWioqLw+OOPX/R3CgoK0NLSgpSUlKkNoS9iIYSQWVvwL2JZtWoVOjs70dvbC7vdjrfeegsNDQ1T\n1tFoNL5GNTc3gzF2QegTQghZOAGd6omJicGLL76IW265BS6XCw8++CDKy8vxm9/8BgCwbds2vPPO\nO3jppZcQExMDkUiEN998c14aTkKD0+nGxIQDExNOGI1OGI0OGI1OmEz8ttnshNns8tVWqwsWC6+t\nVhfsdjdsNl7b7W44HLw4ncxT3HC5GNxu5qn5ZAHG+PeOnk8gAKKiBBAIgOhoAaKiBIiOFiAmRoCY\nmCjfcmxsFITCKMTGRiEuLhpxcbyOj49CQkKMrxaJoiES8ToxMcZXxOIYSCRCT4lBfHz0jE59EhIK\n6Dt3CaxWF0ZHbdBqbdDpbBgdtUOns0Gns0Ov52VsjBeDwYGxMTvGxx0YH3fAanVDIonxBaBYLIRY\n7A9Ib2CKRDFISIhGQkI04uN58YattxYKoyAUCjx1FGJivKEdhagoAaKi/KHurSfzdgaM8Q7C21l4\ni9PJOxRv5+JwuH0djtXqgs3m9nVIFgsv3g7LZOKdmMnEi7ej452eAy4Xg1QqRFJSLKRSIWQyIWSy\nWE8RIiUlDsnJsUhJ4UUuj4NcHofU1DhIpULqNMiczSU7KfgXKbPZCbXajIEBCwYGLBgaskCjsUKj\nsWJ42IqRERtGRqzQam2w2dxITY2DXB6L1NR4XzglJ/uLTMbrpKRYJCXxUJNKhRCJaKQLAHa7C+Pj\nDhgM3sI7SW/H6S2TO1at1obRURusVjfk8likpcUjLS0O6enxSE+Ph0IRj4yMBGRmxiMrS4SsrASk\npsbR9iZTUPBHCKvVhb4+k6eYoVKZoFKZ0d/Pi1ptgdXqQlZWgq9kZiZAoYiHQpGAtLQ4X8ikpcVD\nIomhMAkim80FrZZ3BMPD/o5Zo7FgaMiKgQELBgctUKvNMBqdyMpKgFIp8pWcHBFychKRmytCXl4i\n5HLqHCIJBf8i4XS6ce6cCd3dRvT0TKC314SzZ43o7TWht9cEg8EOpVKE3NxEz0HPD/ycHH8YJCfH\n0sG/CFmtLgwM8A7e29n39fGO/9w5PhCw213Iy0tEfr4Y+fmJKCgQo6BAjCVLxCgqkkAqFQb7zyDz\niII/zOj1NrS3j6OtbRzt7QZ0dEygo2Mcvb1GZGQkoKhIgoICfuB6D+C8vERkZCQgKopCnUxvfNyB\nc+f8A4WzZ43o6eGlq2sCUqkQpaVSlJZKUVbGS3l5EvLyEmm/CkMU/CFqdNSG06fHcPq0AadPG9Da\nyovZ7ERZWZLv4PMejEVFEsTHRwe72WQRcrsZ1GozOjrGPWUC7e0GtLWNY3TUhtJSKSoqklBZmYSK\niiRUVclQUJCI6Gj6IN9QRcEfZIwx9PebceyYDi0tOhw7psOxY3qYTE5UVSV5DigZKiqSUF4uhVIp\notMxJGRMTDjQ3j7uG5icPj2GU6cMGBmxYulSGVasSMGKFSlYuTIFlZVJiI2lwUkooOBfYFqtFc3N\nozh6lJevv9bB7WZYtUqOFSuSfQdKXl4iBTx4x2izuWAyOWAyOWA2O2Gx+IvN5vRMq3TBbnd5plvy\n2jsd0zunn0/ZvHB/4dM+eYmJifIUAYTCaMTGRiE2Nhqxsd7ppNGIj49BQkKMZ85+DEQiIcRiIUQi\nIZ328Bgfd+Cbb/S+AU1Liw69vUZUVclw1VVyXHWVHFdfLUdpqZS2WRBQ8F9BjDF0dIzjiy+GceDA\nMJqatNBqbVi1KgVXX813/lWr5MjJWfyjeJPJgeFhM0ZGzBgZsUCrtWB01IrRUV7r9VYYDDaMjdkw\nPm7HxIQdExMOGI12REdHeeb483D1hi6/aCraM68/xhfS/jn9/OKr6OgoX7ALBJgyl987j5/P3/d3\nFvyCMO+cfd6x2GwuWK1Oz5x93vGYzQ5fp2SxOJGQEAOJJBYSSSykUl5ksjgkJ8cjOTkeqakJkMvj\nIZfHIy1NhLS0BKSni5CcHL/oA9BodOD4cT2+/poPepqbRz3Hgxxr1qShrk6Ba65JpVOWC4CCfx4x\nxtDVNYHPP9fgs880aGzUIC4uCuvWKbBuXTquvTZt0Y1wJibs6O+fgFpthFptxMCACQMDRgwOmjA4\naMLQkAkajRkulxvp6SJf2KWmJvhCMCWFl6SkOCQlxUEqjfWEpxASSSyEwvAIArebwWTindX4uL8Y\nDDbo9VbodFZfZ6fVWqDVWn2dodHoQGpqAjIyEqFQiJCVlYjMTDGyshKRnS2GUilBdrYY6emiRbX/\naLVWHDkyigMHhvHFFxqcOmXAypUpWL9egRtuUKC2NhVxceHx/w8nFPwB0utt2LdvEPv3D+Kvfx2C\ny8Vwww0KrF+fgfXrFSgoEAe1fYGyWBzo6hpDZ+cYurvH0Ns7jnPnxtHXNwGVagI2mwtKpT+YsrPF\nyMxM9BVvkEkkNFX0Uux2F4aHzdBozJ5Ok3eeAwMmX6fa3z/hmZYrRk6OBLm5EuTnS1FQkITi4mQU\nF8uQnh7erx4nJhw4eHAEn3+uweefa9DWZsA116TippsycNtt2aioSArrvy9UUPDPwblzRnz4YT8+\n/LAfR4+OYt06BW69NRM33ZSJkhJJWO6YIyNmnD49irY2HdrbdWhr06GjQweNxowlS3iwFBYmoaAg\nCXl5Ul/wpKTEh+XfG66sVif6+yfQ1zeBc+d4J9zdbUBX1xjOnNHD5WIoKZGhrCwF5eVylJenoKJC\njiVLkhATE36zbMbG7PjiCw327x/ERx+pERMThTvuUGLzZiWuuy6NZg7NEQX/DLW3G/Duuyq8954K\nfX0m3HZbNjZvVqK+PhMiUUCfW7egLBYHTpzQ4sSJEZw8qcWpU6M4dUoLh8ONigoeFN5SWpqCvDxp\nWAZGpBodtaCjQ+/pvHlH3tqqw9CQCSUlyaislGPp0lQsW5aKmpp0ZGaGzytSxhi++UbvG3SpVGY0\nNChx1105uOmmDJoxNAsU/JcwNmbHH//Yi9de64FabcZdd+XirrtycN11aWERhowx9PQYcPjwAJqa\nBtHUNITW1lGUlaWgujoNS5emoqpKjqqqVGRm0iyixcxkcqCtbdTX2Z84MYJjx4YhEsXgmmsysXp1\nFq65JhMrVqQjPj48BjK9vUa8/74K77zTh87OCfz3/56PBx4oxLJlycFuWsij4D8PYwyNjRr89rfd\n+OgjNW6+ORMPPFCI+vqMkH9Z6XS6cfz4ML78Uo2vvlLj0KEBREcLcO21Wb4Du6YmDQkJdPk94ft6\nd/eYZ1AwiMOHB9HersPy5elYsyYLa9ZkY82abCQnxwe7qZfV3T2B11/vwWuv9SAjIx4PPliIe+/N\nR1JSbLCbFpIo+D3sdhf+8Ide/Pu/tyI6WoCHHirCffflIzU1tHf6oSETPvqoBx991IPPP1chJ0eC\ntWv5AXvddVnIy5PSSJ7MmNFox5EjQzh4UI0vv1SjqWkQJSXJ2LixAA0NhVi5UhHSs4pcLjc+/XQI\nr77ajU8/HcJ3v1uAH/+4HLm5icFuWkiJ+OC3Wl347W+78O//3oriYil++tNKrF+vCNmwZIyhrU2H\n99/vwu7d3ThzRo9bbsnHbbcV4JZb8pGWJgp2ExcUYwwWixMTEzbPXHoHzGYHrFan76Iuu93luZCL\nz9OfzHvRlvcz/OPiYhAXxy/YSkgQ+q4dSEzkU0xD/VXffHM4XDh8eBB79vTgz3/ugV5vw+23L8Ed\ndxTixhtzQ/q00OCgBc8914ZXX+3G5s1K/PSnVSgqkgS7WSEhYoPf7WZ45ZUu7NhxEqtWpeBnP6tC\nbW3qPLdw/pw7N45XXz2Jt97qgNXqwp13FqGhoRDXX58dNvPcZ2piwgaVahwqlQH9/eOeqY0TGBoy\nQaezQK+3YGzMirExK4xGO4TCaEgksUhMjPVd4BUfH+MLcaEwesrFXJP7dP6FK7xDcDhcky7WmnyB\nlh1Gox0mkwMJCTGQSuMgk8UjOTkBMlk80tP5nHs+716CnJwk5OYmIT198X2AWVeXHn/+cw8++KAL\n33yjxa235uPBB6tw4425Ifu36nQ2PP98B/7rv85g06YsPPNMDTIzE4LdrKCKyOD/5hs9tm07gpiY\nKLzwwiqsWBGa3+frcrmxb18vdu06gUOHBvCd75Rh69ZKrFiRHrKvSGaCMQaNxoQzZ0bR2TmKM2dG\n0dMzhrNn9ejp0cNqdSI3Nwk5OUlQKqXIyhIjI4MXuVzkCd14JCXFL+gFXm43g9nsgMHAOx293gq9\n3oLhYZOvcxoYMEKlMqCvz4DxcRtyc5OwZEkylixJRmFhMkpK5CgpkaOgIDnsZ6EMD5vxzjtn8PLL\nJ2EyObBt2zLcf38l5PLQDFWDwY5nnjmNV1/txo4dy7BtW1HEvYLziqjgt9lc+PnPv8Hrr/dg587l\neOCBwpAcpVgsDuzadQL/+Z9/Q0pKPL7//Wrcc08pEhPD701Zu92FlpYBNDX149SpYbS1adHaOoKY\nmCiUlqaipESO4uIUXzgWFMiQmhreFyF5WSwOnDtnQE+PHt3dOnR16dDZqcOZM6Po7x9Hfr4MlZXp\nqKhIRU1NJtasyUV6evidi2aMoalpEC+99A3+/OcebN5chMcfvwplZaE5oGptNWDbtiOw29343e+u\nQWWlLNhNWnARE/xarRV33nkAqalx+M1vapGeHnpv2jLG8OabHfjJT77EypUKPP74VaitzQx2s2bF\nbnehuVmNzz47i88/78XRo2oUF8tx7bVKLFumQEVFGsrL05CaGlnvRZzPZnOis1OH1tYRnD49jKNH\nB3DwoApZWRJcf30ubrihADfcUIC0tPDqCLRaC15++QSee+4YtmwpxRNPrEZqaui9AnC7GV59tQs/\n+9k3eO21a3DbbcpgN2lBRUTwt7UZcNttjdiyJQ9PPlkdkqP848eH8cgjn8NkcuD55+uwdm347Ig9\nPXrs3duJffu6cODAORQXy3HDDflYv74A112Xg6Sk0OtkQ5HL5caJExo0Nvbis896ceDAORQUyHDL\nLYW49dYirF2bFxbXjwC8A9i+/RD+9KczeOKJ1fjHf1wWkqdVmpq0+Na3DuCxxyrwT/9Uuiheac7E\nog9+tdqMq6/eh6efrsbf/33hArVsdl5//TT+5V8O4Omnr8MDD1SF5AFyPsYY/vKXbvzHfxzG3/42\nhE2bSrBhQxFuvLEAcnlkj+bni8PhwtGjA/jLX7rw0Ued0GrN+P73V2HbtlWQycKjMz15cgTf//5n\nkMni8Oabm0LydGVvrxG33vo5/umfSvHwwyXBbs6CWNTBb7O5UFf3KRoasvGv/1q1gC2bGcYYnn76\nCH73u9PYs2czysvlwW7SZdlsTvzhDyfxH/9xGFFRAvzzP6/Gli1ViIsL3Wl9i0VLywD+9/8+gj17\nzmDr1mo8+ug1yMsL/fPTDocL27Z9ihMntNizZzMUitA7fdXVNYHrrtuP995bi+uuSw92c664OU2M\nYQH6+OOPWWlpKSsqKmLPPPPMtOv86Ec/YkVFRWzZsmXs2LFj065zuaY88cQ37I47Gpnb7Q60yVfE\n//pfR1l19RtsYGAi2E2Zkfb2EbZkyfOsvv4Ntn9/V8hu18Wur2+M/fjHf2EpKc+yX/3qULCbMyNu\nt5v9z/95kFVU/B9mszmD3Zxp7dnTz3Jy3gvZ9s2nucR4QMHvdDpZYWEhO3v2LLPb7ay6upq1trZO\nWWfPnj1sw4YNjDHGmpqaWG1t7fQNuUTj7XYXy8h4h7W1jQXS3CtGrZ5gcvl/se5ufbCbMiMtLQMs\nI+OX7He/m74TJgvv3LkxVlLyn+zf/u2zsOmEN216j/3yl0eD3YyLWrduP3v77XPBbsYVN5fgD+gE\ndHNzM4qKipCfnw+hUIgtW7bgww8/nLLO7t27sXXrVgBAbW0txsbGoNFoZvU8H32kRmmpFGVlSYE0\n94rZseMwHnpoKZYsCf2X6m1tI7j11t/j17/eiPvvrwl2c4hHbm4Svvzyfvz5z2fw5JMHgt2cGfnV\nr9bhmWeOwmCwBbsp0/qHfyjCq692BbsZISmg4Fer1cjJyfHdViqVUKvVl12nv79/Vs/T02PEypWh\nOY8YAE6c0OKOO0LzzebzNTX1Y8OGYtx5Z3mwm0LOk56eiOefvxX793cHuykzUlqaAoVCBJVqIthN\nmdbKlXL09BiD3YyQFNC7eDOdLsXOe+PhYr+3fft233JdXR3q6uoA8Dd2Q/kr28xmfvl/ONDpLJDJ\n4oLdDHIRcnkCdDpLsJsxYyJRDMxmR7CbMa24uCjYbK5gN2PeNTY2orGxMaDHCCitsrOzoVKpfLdV\nKhWUSuUl1+nv70d2dva0jzc5+CfLzEzARx+pp/1ZKCgoSMLBgwNYvjz0ZxBUV2fgl788jB/84GqU\nlIT+zKNI4nK58eSTB7BiRXhc6Dc6akF3twHZ2aH5BTAdHeOL8nN8Jg+KAWDHjh2zfoyATvWsWrUK\nnZ2d6O3thd1ux1tvvYWGhoYp6zQ0NOCNN94AADQ1NUEmk0GhUMzqeTZvVuLTT4dgMNgDae4V87Of\n1eKZZ5phszmD3ZTLuummJXjqqfW4+eb/C5XKEOzmEA/GGH74w73QaEx49dWGy/9CCHjuuWP41reK\nkJ0dmp+S+cc/nsM99+QFuxmhKdB3lPfu3ctKSkpYYWEh27lzJ2OMsV27drFdu3b51vnBD37ACgsL\n2bJly1hLS8u0j3O5ptxzz5fspz89Hmhzr5i77trN7rtvL3M4XMFuyow899xhlpr67+zJJ79gY2OW\nYDcnYrndbvbRRx3s6qtfYatXv8oMBmuwmzQjH3/cw9LSfs16ekJzpl17u4HJ5W8ztdoU7KZccXOJ\n8YCDf77/6X8HAAAeGklEQVRcrvGDg2aWmfku++yzwQVq0eyYTHZ2883vsLvv3h02c4fPnNGy7373\nPSaXP8t+9rO/spGRxX+QhAqXy83efbeV1dTsYkuX/pq99dYp5nKFxzTOd989w9LSfs0OHVIHuynT\nslicrLp6D3vppY5gN2VBzCX4w+bKXQD4y18GcP/9Tdi3b31IfhenzebEli17MDJiwcsv16OiIjzO\noff06PHss1/h7bdbceONS7BpUzE2biwOy0+XDGUOhwtffdWHjz46gw8+6EBKSgL+7d+ux223lYTk\nZ06dz2JxYOfOZrz66ins2bMZK1bM7pTtQrDbXfjudw8BAN58c01EfF7Pov7IBq833+zFI498jbfe\nWoP16zMWoGWz43K5sWvXCWzffhgPPVSFn//8mpD8TJPpDA+bsGfPGezZ04lPP+1BSYkc9fVLcPPN\nhVi9OifsP3M+GHp69PjLX7qwf38PPv/8LEpK5LjtthLcdlsJamoywiaYPv74LH74w8+wcqUCzz23\nLiTP6xsMdtx115eQSGLw//7fdWEz0y5QERH8ANDYqMG3v/0lnnyyGv/wD0UhefAMDhrx4x8fwJdf\nqvHzn9fiu98tD6svRrfbXTh0SIVPPunGJ5/0oL1di9Wrc3D11VlYvjwDS5cqUFiYHBYfQrdQRkZM\nOHVqGCdOaNDSMogDB87BZnPh5psLcfPNS1BfXxh2r6IOHx7Azp3NaG/X4YUX1mPDhoJgN2lara0G\n3HvvV1izJg0vvLAqovbLiAl+gP+jt249BIlEiJdfrg3Z79/86is1nn32KI4cGcT3vrcUDz9cDaUy\nNNt6KTqdBV9+eQ4tLYP429+GcOrUMIaGjCgpkaOyMh2lpXJPSUVhYTIkksV5rYDT6UZfnwGdnaPo\n6BhFe7sW7e1anD49ArvdhaqqdCxdmo6amgysXZuH0lJ5SA5MLsVud+Gdd87ghReOY3jYgkcfrcG2\nbctC8sP77HYXnnmmFS+80IGnnlqGbduKw257Byqigh/gB+ELL3Rg587TeOyxcjzySBni40PzdMSZ\nM3r8538exx/+0IYbbsjFtm3LUFenDOvv2DWZ7L5v4ero0KKjg4dhT48eiYlCz7dwJSM3V4rc3CRk\nZkqQnp6I1FQR0tL41y6GysiMMQaTyYHRUTNGRswYGTFBozFBpTLg3DmD79u3+vvHkZEhRlFRCsrK\n5CgrS0VpaSqqqtKRmSkO69Dp6tLj979vw29+cxIVFSn40Y9qcPvtS0Lmf3S+r74axsMPNyMvLxEv\nvXQ1cnLC69XUfIm44Pfq6ZnA//gfx9DSMoqf/KQSDz1UFLIdwMSEHW+80YrXXz+Nnh4DGhoK8a1v\nFaO+PjckR1RzwRjD0JARPT16nD075gvPoSGjL1RHRswYH7dBIolFcnICkpPjIZHEQSqN83zZuhAi\nES+Tv2w9Li7G90Xr0dGCKUHrdjM4nW5fsdmcU75s3WTiX7JuMjkwPm7zFb3eAp3OgpiYKMjlvFNK\nT09EWloicnKkyMtL8n3fbn6+bFH9n9radHjvvU68804nhoZMuPvuYjz8cDUqK1OD3byLOnx4BNu3\nn0RHxzh+8Yvl2LIlL6w73EBFbPB7tbSMYseOkzh2TIfHHqvAQw8VQSQK3YO0r28c773XhXff7cTJ\nk1rcfHMe7rijEBs3FiA5OTy+nCMQLpcbBgMP3rExKyYm7L4wNpsdMJsdMJnsvvD21i4Xg8vlhst1\n4UeBCIW8Q4iJiZrSWSQkxCAxkXcoiYmxkErjkJQUB4kkDsnJ8UhOTkB8fOjuK/PF6XTjyJFBfPhh\nNz74oAsWixN33lmEu+8uwXXXZYXs6J4xhoMHR/Dkk6fQ0TGOn/2sClu3FtCEA1Dw+7S0jOKpp07h\nyy9HcN99+di2rRjl5aH5yZ5ew8NmfPRRDz74oAuff65CUZEM69YpsW6dEmvWZCMtjb4Ji8ye3e7C\nsWMaHDigRmOjCgcPDiA/X4o77ijE5s1FqKlJD+nR8sSEA3/8Yy927erE+LgDP/lJJf7u7yjwJ6Pg\nP09vrxGvvtqF3/62G4WFEnzve0X4b/8tN6RfBQD8YG1p0eCLL/rR2KjC4cODyMoSY82aLKxZk43V\nq7NQXCwL6QOWBIdeb0VT0yAOHRrAV1+pcfSoBkVFMlx/fTbq6nKwdm3oDyIYY2huHsUrr3Th3XdV\nqKtLxz/+YzHq6zPD4nqHhUbBfxEOhxt796rxyitdOHRIi9tvz8bdd+eivj4zZN8LmMzlcuPkSS2+\n+kqNr74awOHDAzCZHFi5UoHq6jQsW5aK6uo0lJam0EgoQjDGMDBgxIkTWnzzzYivqFQTuOqqDKxe\nnekZJGSGxXf6MsZw+rQB77zTh7ff7oPN5sJDDxXh7/9+CTIyFt8Hrc0nCv4ZUKvNeO89Fd59tw9/\n+5seGzZk4c47c7BhQxYkkvCZZz8wYMSxY8NTDvq+vnEUFclQVZWKqqpUlJenoKJCjsLCpLCePRTJ\nGGMYHjajtXUUra06tLaO4uRJLU6d0iI6OgrLlqVi+fJ0VFenobo6DZWVcsTEhOZ5+vMxxvD11zq8\n/74K772ngtnsxN135+Kuu3KwenUaje5niIJ/loaGLPjww368/74Khw6N4Prr07FhQxbq6zNRXCwJ\nu1MpFosD7e16nDrFg6GtTYe2Nh1UqgkUFCShtDQZxcUyLFkiQ36+FPn5UuTlSSEShU+Htxi5XG4M\nDprQ2zuOs2cN6O0dR3f3GDo69Ghv10EgACorU1FRwTvyqqpULF2aivT00D5lMx2DwY7GRg327x/E\nn/+shkgUgzvvVOLOO3Nw1VXhd81DKKDgD8DYmB0ffzyA/fsH8ckng4iOFuCmmzJxww0K3HBDRlh/\nrrfV6kRnpx7t7Xp0dY2ht5eHS2/vOM6dG4dEEovcXAlycyXIyZEgO1sMpZLXWVliZGUlQiyODfaf\nEZacTjc0GhMGBkxQq41Qq43o759Af78RfX3j6OubwMCACSkp8cjPl6KgQIr8/CQsWcI76rKyFKSm\nJoRtINpsLhw+rMVf/zqEv/51CCdPjmH16lTU12di06ZsVFSE9qSLcEDBP08YY+joGPfsrBo0NmqQ\nkRGP9esVWLdOgWuvTYNSGX6jrem43QwjI2b09U2gr28c/f3+YOrvN2Jw0IiBAROiowXIyEhERkYi\nFAoRFAoR0tISkJ4uQmpqgmdKJJ8aKZPxqZLhcsphphhjMBodGBuzYWzMCr3eBr3eCp3OiuFhM4aH\nLRgeNkOjMUOjMWFoyAydzoq0tARkZYmRmZkIpZJ3qjk5EiiVYuTlSaFUihfNtQFmsxNffz2KAweG\n0dg4jCNHtKioSMKNN2bgxhszcO21qRHzGToLhYL/CnG53Pjb3/RobBzGF19ocPiwFnFxUaitTUVt\nrRy1talYuTIFYvHiPGXCGIPBYPMEmhlDQyYMD5sxMsKDTqu1+EJQr+ehOD5uR1xcNKTSOEilsZBI\nYiGRCCEWe+fSx3gu0OJz7BMSYjwXakUjLi4asbHeEoWYGF74hVsCREX5y+SRMGMMbjcvfK6//4Iu\nh8MNu90Fu91/YZfV6oLVyi/uMpudnmsHnDAa7TCZeD0+bvdcX8Dr+PgYyGRxkMnifJ1dSko80tNF\nSE/nnaFCIfJ1kGlpokXXAXq53Qzt7eNobtbiyJFRHDmiRXv7OJYuleG669Kwfr0C11+fjqQkerV4\nJVHwLxDGGM6eNaKpie/wzc2jOHFCj4ICMa66So5Vq1KwcqUcy5bJQn7q6JXi/QgEg8GGiQkHJiZ4\ncHqvnDUaHb6gtVh44WHMa7vdBYfDDZvNBafT7Qtxl4sHO2O8Qz5fdHQUoqKAqCjBlA5jckciFEZN\n6WhEIiESEmIgEvHOSCz2Ft5Z8U4rdlG+ipkpt5vv8y0tOnz99Si+/lqHlhYd0tLicNVVfPBzzTWp\nWL48OSxmyi0mFPxBZLe7cOqUAUePjvoOjra2ceTmilBVJcPSpTJUVclQWZmEoiIJhMLIDBAS2hhj\n0GisOH3agNOnx3DyJC+nTxuQnByLlStTsGoVH9ysWiWHXL44P4wvnFDwhxiHw40zZ8Zx6hQ/eE6d\nMuD0aQP6+80oKhKjvDwJ5eVJqKiQoqwsCcXFkoh9hUAWltvN0NdnQkfHONraxtHaakBbmwGtrfx7\nmL2DFP+gJQnJyRTyoYiCP0xYLE60t4+jrc3gO+ja28fR02NEenocioulWLJEjMJCMZYsEaOggJeU\nlNiwnd1BFp7V6sK5cyacPWvE2bNGdHcb0dNjRHf3BLq6JpCcHIvSUqlnACJFRQUfiCgU8bSfhREK\n/jDncrnR22tCZ+fElAP17FkjentNcLncyMtLnFJyc3nJyREhMzMhYs9BRxrGGMbG7FCpzOjrM/nq\n3l4Tzp3jRau1ISdH5Bs4FBaKUVgowZIlYhQXS8LqgkVycRT8i9zYmB29vUbfgd3b6z/g+/pMGB21\nQ6GIh1IpQnZ2AjIz/UWhiEd6eryvjoujN+BCkdvNoNPZoNFYp5TBQQsGBy0YGLCgv9+M/n4zoqMF\nUCpFyM0V+QYAkwcF2dkJIftpm2T+UPBHOIfDjcFBfzAMDFgwNMQDY3jYHyIjIzYkJEQjLS0OaWnx\nSE2N8xW5PA4pKbG+OjnZW+IgFsfQZfSzYLW6MDZmh05n90x35cujozZfPTJig1Zrg1ZrxfCwDTqd\nDRKJEApF/KSSgMzMeGRkJCA7W4ScHBGys0WQSmnETij4yQzxefkOjIxYMTxsxeio3Rc+o6N26HQ2\nX63X232BZbG4IJUKIZMJkZQUC6k0BlKpEFKp0DPtUQixOMZXEhN5SUiI9szVj0Z8vLdEITY2GnFx\nUYiNjfJNs1yoc8uM8Xn+fMqo2zOFlC/z+f0uWCze4oTF4vLM7fcWByYmnJiYcMBodMJgcGB83AGD\nwQ6DwYGxMTvcbkAmEyIlJc7Xgcrl3k41DnJ57JSO17tMM77IbCxo8Ot0Otxzzz04d+4c8vPz8ac/\n/QkymeyC9fLz8yGVShEdHQ2hUIjm5uZ5azxZWE6n2xdq4+MOXzEYHFOCcHJATg5Oi8U1JVS9Qeu9\nsMrpZIiO5l+mIhR65+Hz2js3n1+0BV994QVcky/kgu9CLv/FXDzsHQ6354tbBJ6LxnjnExcXPaWD\n4p0WL3xuP+/MxOIYT2fHa6lUiKQkbx0LmUyI+PhoepOUXHELGvyPPfYYUlNT8dhjj+HZZ5+FXq/H\nM888c8F6BQUFaGlpQUpKyqUbQsEf8RhjcDoZ7Hb/BVsOh/eiLf7mtzfYGePnw883uXPwXuXLv5Fr\nakciFAro/DdZFBY0+MvKyvDFF19AoVBgaGgIdXV1aG9vv2C9goICfP3115DL5ZduCAU/IYTM2lyy\nc85DHo1GA4VCAQBQKBTQaDQXbdRNN92EVatW4ZVXXpnr0xFCCJknl7xMtL6+HkNDQxfc//TTT0+5\nLRAILnou8+DBg8jMzMTIyAjq6+tRVlaGtWvXTrvu9u3bfct1dXWoq6u7TPMJISSyNDY2orGxMaDH\nCOhUT2NjIzIyMjA4OIj169dPe6pnsh07dkAsFuOf//mfL2wIneohhJBZW9BTPQ0NDXj99dcBAK+/\n/jo2b958wTpmsxkTExMAAJPJhP3792Pp0qVzfUpCCCHzIKDpnN/+9rfR19c3ZTrnwMAAvve972HP\nnj3o6enBt771LQCA0+nEd77zHfzrv/7r9A2hET8hhMwaXcBFCCERZkFP9RBCCAlPFPyEEBJhKPgJ\nISTCUPATQkiEoeAnhJAIQ8FPCCERhoKfEEIiDAU/IYREGAp+QgiJMBT8hBASYSj4CSEkwlDwE0JI\nhLnkF7EQQggJPrcb6NUAp3uB0+eAU718+a7pv9PqsujTOQkhJERY7UCXGujoB9r6gHaVv5ZLgIo8\noDIPqMoHqgqAilxALKKPZSaEkJDGGDAwCnSogDP9POQ7VLxWa4H8DKBUCZTnAmU5vFTkAdLE6R+P\nPo+fEEJCgNsNDOqAnkE+gu8aADrV/iIRAcXZPOBLc4CSbF4vyQSEszwBT8FPCCELxGzl5917BoGz\nQ0D3AF/u9tyWioDCTKAom9fF2Z6iBJIuMnqfCwp+QgiZJ1Y7oBoGzg0DvUM8zCeXMSOQp+Cj9IIM\nHu5LJhWJaGHaGfbBb7YyJMQFuyWEkMXO5QI0ekA14i99w/66bxjQGwFlKpCbzoM9X8Hrgky+nCUH\nokJgQnzYB3/cBgZxApCTBijT+EbPTr2wlogAgSDYLSaEhCK7AxjS8zdK+0cA9Siv+7X+elAHJIt5\npuSkA7lpvM5JA/LSedhnpIRGsF9O2Ae/y8WgNfBeV63ldb+WL6u1/mUAyEwBMuW8zvLUGSn+OiMZ\nkEvD4x9HCLk8owUY0vFQH9LxMjipDIzyeswIKJJ5LignDRi9A8rsVCBbDsTFBvsvmh9hH/wzbcqE\nmf+TB0an/vMnL2v0wLgZSJXyjiBdxncGRTKgkPlvp3uWU5OAWOEV/iMJIT6MAfoJYHgMGDHwY3Z4\njNcaPaCZvKwH3IwP6HwDvGT/4M97X3YqP+ajo4P91y2ciAn+mbI7+I40pJu0Q435l731iAHQGgBx\nApCWxEuqp8ilfEeST1OSxdRZEALwEDeYAN0EMDp+YdEaAK2nHvGU0XEgMZ4PvNKSpg7EfIM0z0At\nI4Ufn3SK90ILGvxvv/02tm/fjvb2dhw9ehQrVqyYdr19+/bh0UcfhcvlwkMPPYTHH398+oYEeVaP\n281fIo4YeIeg9eyY3p1Va5i0U0/wZd04kBDn7wRSJECyxFOL+XKyGJCJAVkiv+1dTkpcPC81yeLg\ncgETFn4cjBl5kOs9y3ojH53rJy3rPMV723sspEj4VabeAVJqkn/wlOoZWKXJ+H10DARuQYO/vb0d\nUVFR2LZtG371q19NG/wulwulpaX49NNPkZ2djauuugp//OMfUV5ePi+NDza3mx8oo+P+g2J0fOqB\nMmby/8x7MI2Z+HKUgHcAMjGvpSJ/LRXxK/UkCfzN7Mm1eHKJ56OmWCGNhiKRywWYrLwYLYDRU0+Y\n+b55fj1u5vvguMlTm/37pMnK9yeZ2L9PJov9gxfvQMY7qPEOeLyDndleeETmx1yyc87/qrKyssuu\n09zcjKKiIuTn5wMAtmzZgg8//HDa4A9HUVH84EhKBJA5u99ljM8T9nYGkw9Cg8lzkJr4K5CeQX7b\naJlUew5k7wEP8A4gMR4QxU1annRfQpy/Toj1L8fH8tve5Tghr73L3tuxMXw51nNfTDR1NpO5XIDd\nCdgcvNgd/mWbg/+/vbXVDljsgMU2ddls47XFzi8QMnvu8y57Q95b7E7+fxQn8P/z5IHB+YOFtCSg\nKIvf591vpZOWJQmRdW48kl3RPlqtViMnJ8d3W6lU4siRI1fyKcOGQOAJ4Dj+BlWg7I4LQ8EbGCbr\nhUFisfFTVhbt1CCaEky2qaE1OdQcTsDp4p1AbAwf7XlrYbSn9izHnFeio/z15OWoKH8dJfCUKL6t\nBJ5t5i0Av28yBt6hwlMz5r/P7ebLbjd/k9Dl5kHtmua20zWpeG47nJ5y3rLd4a/dzN85xk3qKCd3\noAlxQPyk+7ydcLyn4xXH84AWxfFO29s5ezvvySGfGM9/Rp0vma1LBn99fT2GhoYuuH/nzp24/fbb\nL/vgglnukdu3b/ct19XVoa6ubla/H8liPSPxZMnCPafbzTuDyeF3fjA6nDxMvR3F5GB1uf23Xe5J\noezyh7TLPTXApwQ7Lgz/KR3DpA7D24lEeX4eHcVHt1ECfwcU7emAvJ1VdPSk5aipnZnQE+re23Gx\nfB0KYXKlNTY2orGxMaDHuGTwf/LJJwE9eHZ2NlQqle+2SqWCUqm86PqTg5+Evqgo/ykhQsjCOH9Q\nvGPHjlk/xrxc3nSxNxZWrVqFzs5O9Pb2wm6346233kJDQ8N8PCUhhJA5mnPwv//++8jJyUFTUxM2\nbdqEDRs2AAAGBgawadMmAEBMTAxefPFF3HLLLaioqMA999yzaN7YJYSQcLWoL+AihJDFbi7ZSZ9k\nQwghEYaCnxBCIgwFPyGERBgKfkIIiTAU/IQQEmEo+AkhJMJQ8BNCSISh4CeEkAhDwU8IIRGGgp8Q\nQiIMBT8hhEQYCn5CCIkwFPyEEBJhKPgJISTCUPATQkiEoeAnhJAIQ8FPCCERhoKfEEIiDAU/IYRE\nGAp+QgiJMBT8hBASYSj4CSEkwsw5+N9++21UVlYiOjoax44du+h6+fn5WLZsGWpqanD11VfP9ekI\nIYTMkzkH/9KlS/H+++/j+uuvv+R6AoEAjY2NOH78OJqbm+f6dCGvsbEx2E2Ys3BuO0DtDzZqf/iZ\nc/CXlZWhpKRkRusyxub6NGEjnHeecG47QO0PNmp/+Lni5/gFAgFuuukmrFq1Cq+88sqVfjpCCCGX\nEXOpH9bX12NoaOiC+3fu3Inbb799Rk9w8OBBZGZmYmRkBPX19SgrK8PatWvn1lpCCCGBYwGqq6tj\nLS0tM1p3+/bt7Je//OW0PyssLGQAqFChQoXKLEphYeGsc/uSI/6ZYhc5h282m+FyuSCRSGAymbB/\n/3488cQT067b1dU1H00hhBByGXM+x//+++8jJycHTU1N2LRpEzZs2AAAGBgYwKZNmwAAQ0NDWLt2\nLZYvX47a2lrcdtttuPnmm+en5YQQQuZEwC42XCeEELIoBe3K3XC+AGymbd+3bx/KyspQXFyMZ599\ndgFbeGk6nQ719fUoKSnBzTffjLGxsWnXC7VtP5Pt+cgjj6C4uBjV1dU4fvz4Arfw0i7X/sbGRiQl\nJaGmpgY1NTV46qmngtDK6T3wwANQKBRYunTpRdcJ5W1/ufaH8rZXqVRYv349KisrUVVVhRdeeGHa\n9Wa1/Wf9rsA8aWtrYx0dHZd9czg/P5+Njo4uYMsubyZtdzqdrLCwkJ09e5bZ7XZWXV3NWltbF7il\n0/uXf/kX9uyzzzLGGHvmmWfY448/Pu16obTtZ7I99+zZwzZs2MAYY6ypqYnV1tYGo6nTmkn7P//8\nc3b77bcHqYWXduDAAXbs2DFWVVU17c9Dedszdvn2h/K2HxwcZMePH2eMMTYxMcFKSkoC3veDNuIP\n5wvAZtL25uZmFBUVIT8/H0KhEFu2bMGHH364QC28tN27d2Pr1q0AgK1bt+KDDz646Lqhsu1nsj0n\n/121tbUYGxuDRqMJRnMvMNP9IVS29/nWrl2L5OTki/48lLc9cPn2A6G77TMyMrB8+XIAgFgsRnl5\nOQYGBqasM9vtH/If0hauF4Cp1Wrk5OT4biuVSqjV6iC2yE+j0UChUAAAFArFRXeQUNr2M9me063T\n39+/YG28lJm0XyAQ4NChQ6iursbGjRvR2tq60M2cs1De9jMRLtu+t7cXx48fR21t7ZT7Z7v952U6\n58WE8wVggbZdIBBciWbN2MXa//TTT0+5LRAILtrWULr4bqbb8/xRW7D/D14zaceKFSugUqkgEonw\n8ccfY/PmzThz5swCtG5+hOq2n4lw2PZGoxF33303nn/+eYjF4gt+Ppvtf0WD/5NPPgn4MTIzMwEA\naWlpuPPOO9Hc3Lwg4RNo27Ozs6FSqXy3VSoVlEploM2asUu1X6FQYGhoCBkZGRgcHER6evq06wVr\n209nJtvz/HX6+/uRnZ29YG28lJm0XyKR+JY3bNiA73//+9DpdEhJSVmwds5VKG/7mQj1be9wOHDX\nXXfhvvvuw+bNmy/4+Wy3f0ic6rnYuTWz2YyJiQkA8F0AdqlZBcFwsbavWrUKnZ2d6O3thd1ux1tv\nvYWGhoYFbt30Ghoa8PrrrwMAXn/99Wl3pFDb9jPZng0NDXjjjTcAAE1NTZDJZL5TWsE2k/ZrNBrf\n/tTc3AzGWMgEz+WE8rafiVDe9owxPPjgg6ioqMCjjz467Tqz3v7z9c7zbL333ntMqVSy+Ph4plAo\n2K233soYY0ytVrONGzcyxhjr7u5m1dXVrLq6mlVWVrKdO3cGq7lTzKTtjDG2d+9eVlJSwgoLC0Om\n7YwxNjo6ym688UZWXFzM6uvrmV6vZ4yF/rafbnvu2rWL7dq1y7fOD37wA1ZYWMiWLVs2448SWSiX\na/+LL77IKisrWXV1NVu9ejU7fPhwMJs7xZYtW1hmZiYTCoVMqVSy3/72t2G17S/X/lDe9l9++SUT\nCASsurqaLV++nC1fvpzt3bs3oO1PF3ARQkiECYlTPYQQQhYOBT8hhEQYCn5CCIkwFPyEEBJhKPgJ\nISTCUPATQkiEoeAnhJAIQ8FPCCER5v8DC532Kmk7KrAAAAAASUVORK5CYII=\n",
"text": [
""
]
}
],
"prompt_number": 171
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": []
}
],
"metadata": {}
}
]
}