mp574_lecture_008.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lecture 8: Line search algorithms\n",
    "\n",
    "Goal: to describe 1D minimization algorithms. \n",
    "\n",
    "These are important for two main reasons: 1) they provide a preview of higher-dimensional optimization algorithms, and 2) they are often at the core of individual iterations within higher-dimensional optimization algorithms. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## MP 574 Lecture 8: Line search algorithms\n",
    "##\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import display, Image\n",
    "import matplotlib.image as mpimg\n",
    "from os.path import dirname, join as pjoin\n",
    "from scipy import signal\n",
    "import scipy.io as sio\n",
    "import scipy.optimize as opt\n",
    "import numpy.random as rnd\n",
    "import time\n",
    "\n",
    "font = {'weight' : 'normal',\n",
    "        'size'   : 20}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's define a function to be minimized\n",
    "def myfunction(x):\n",
    "    p = np.array([0.005, -0.42, 4, 400, 12])\n",
    "    p1 = np.polyder(p)\n",
    "    p2 = np.polyder(p1)\n",
    "    f = np.polyval(p,x)\n",
    "    f1 = np.polyval(p1,x)\n",
    "    f2 = np.polyval(p2,x)\n",
    "    return f, f1, f2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Minimizer xopt\n",
      "-13.75\n",
      "Value of f at xopt\n",
      "-3461.19091796875\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD4CAYAAAAO9oqkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXRV53nv8e+jESEQmoWQhMUgbCYzKYCNE+qJwRNOmrh205gmbsjgNMlqblqnXfe6cdIu9+amSZy2bp3YDU7TOM5UEyc2wXgewBYGGwMGyWBAIIRAQggNaHruH2fjHNsCCZC0dXR+n7XOOme/+91Hzw6Oftr7fffe5u6IiEh8Swi7ABERCZ/CQEREFAYiIqIwEBERFAYiIgIkhV3AucrNzfXS0tKwyxARiRmbNm064u55Pa2L2TAoLS2loqIi7DJERGKGme093TqdJhIREYWBiIgoDEREBIWBiIigMBARERQGIiKCwkBERIizMGjv7Obep9/iucq6sEsRERlS4ioMkhONHzy3mzVbDoZdiojIkBJXYWBmzC7JZMv+Y2GXIiIypMRVGADMKcmkqu4Ex9s6wi5FRGTIiLswmD0+E3d4fX9j2KWIiAwZcRcGs0oyMYPN+xrCLkVEZMiIuzDIGJHM5LxRGjcQEYkSd2EAMLskk837j+HuYZciIjIk9BoGZnahmW2Jeh03sy+bWbaZrTOzyuA9K+hvZnaPmVWZ2etmNjfqu1YG/SvNbGVU+zwz2xpsc4+Z2cDsbsSc8VnUN7ezr75lIH+MiEjM6DUM3H2nu89299nAPKAF+DVwB7De3cuA9cEywHKgLHitAu4FMLNs4E5gATAfuPNUgAR9VkVtt6xf9u40ZpdkAuhUkYhI4GxPE10JvOXue4EVwOqgfTVwY/B5BfCgR2wAMs2sEFgKrHP3endvANYBy4J1Ge7+kkfO2zwY9V0DYkrBKEamJLJ5n8JARATOPgxuBn4afC5w9xqA4D0/aC8C9kdtUx20nam9uof29zGzVWZWYWYVdXXnfkuJpMQEZhaNYbOODEREgLMIAzNLAW4Aft5b1x7a/Bza39/ofp+7l7t7eV5ej8907rM547PYfrCRto6u8/oeEZHh4GyODJYDr7p7bbBcG5ziIXg/HLRXAyVR2xUDB3tpL+6hfUDNGZ9JR5ezveb4QP8oEZEh72zC4Bb+cIoIYA1wakbQSuCRqPZbg1lFC4HG4DTSWmCJmWUFA8dLgLXBuiYzWxjMIro16rsGzJxgEFnjBiIikNSXTmY2Erga+ExU893Aw2Z2G7AP+FjQ/jvgGqCKyMyjTwK4e72ZfQN4Jeh3l7vXB58/B/wISAMeC14DKj9jBEWZacGVyBMG+seJiAxpfQoDd28Bct7TdpTI7KL39nXg9tN8zwPAAz20VwAz+lJLf9IdTEVEIuLyCuRT5ozPpLqhlbqmk2GXIiISqrgPA9BN60RE4joMpo8bQ0piApsUBiIS5+I6DEYkJzKjKIOKtxUGIhLf4joMAMpLs9larYvPRCS+KQwuyKK9q5s3DujJZyISv+I+DOZdELlxasVenSoSkfgV92GQMyqVibnpVLxd33tnEZFhKu7DACJHB5v2NujJZyIStxQGQHlpFg0tHbxV1xx2KSIioVAYEJlRBLBpr04ViUh8UhgAE3PTyU5P4RVdbyAicUphAJgZc8dHxg1EROKRwiBQXprFniPNHDmhm9aJSPxRGATKg+sNdHQgIvFIYRCYWTyGlKQEXW8gInGpT2FgZplm9gsze9PMdpjZJWaWbWbrzKwyeM8K+pqZ3WNmVWb2upnNjfqelUH/SjNbGdU+z8y2BtvcEzz+clClJiVycdEYXYksInGpr0cG3wMed/eLgFnADuAOYL27lwHrg2WA5UBZ8FoF3AtgZtnAncACYD5w56kACfqsitpu2fnt1rmZV5rFGwcaaW3XTetEJL70GgZmlgF8CLgfwN3b3f0YsAJYHXRbDdwYfF4BPOgRG4BMMysElgLr3L3e3RuAdcCyYF2Gu78UPDLzwajvGlQLJmTT0eV62I2IxJ2+HBlMBOqA/zSzzWb2QzNLBwrcvQYgeM8P+hcB+6O2rw7aztRe3UP7+5jZKjOrMLOKurq6PpR+dspLs0kw2LBH4wYiEl/6EgZJwFzgXnefAzTzh1NCPenpfL+fQ/v7G93vc/dydy/Py8s7c9XnIGNEMtPGZbBx99F+/24RkaGsL2FQDVS7+8Zg+RdEwqE2OMVD8H44qn9J1PbFwMFe2ot7aA/Fwgk5bN5/TA+7EZG40msYuPshYL+ZXRg0XQlsB9YAp2YErQQeCT6vAW4NZhUtBBqD00hrgSVmlhUMHC8B1gbrmsxsYTCL6Nao7xp0Cybm0N7ZzWv7j4VVgojIoEvqY7+/BH5iZinAbuCTRILkYTO7DdgHfCzo+zvgGqAKaAn64u71ZvYN4JWg313ufurk/OeAHwFpwGPBKxTzS7Mxg4176lkwMSesMkREBlWfwsDdtwDlPay6soe+Dtx+mu95AHigh/YKYEZfahloY0Ymc9HYDDbuOUpklquIyPCnK5B7sGBCNpv2NtDe2R12KSIig0Jh0IOFE7Np6+jm9WqNG4hIfFAY9GD+hMhYwUZdbyAicUJh0IPs9BSmFIxig643EJE4oTA4jQUTcti0t4GOLo0biMjwpzA4jYUTc2hp7+KNA41hlyIiMuAUBqcxf0I2oHEDEYkPCoPTyBudyuT8Ubz0lsYNRGT4UxicwaWTcnh5T72uNxCRYU9hcAaLJufS2tHFFt2nSESGOYXBGSycmEOCwQtVR8IuRURkQCkMzmBMWjIzi8YoDERk2FMY9GLR5Fy27D/GiZOdYZciIjJgFAa9WDQ5l85u5+U9mlUkIsOXwqAX8y7IIiUpgReqFAYiMnwpDHoxIjmRD5RmadxARIY1hUEfXDoplzcPNXHkxMmwSxERGRB9CgMze9vMtprZFjOrCNqyzWydmVUG71lBu5nZPWZWZWavm9ncqO9ZGfSvNLOVUe3zgu+vCra1/t7R83HZ5FwAXtTVyCIyTJ3NkcHl7j7b3U89/vIOYL27lwHrg2WA5USeF1kGrALuhUh4AHcCC4D5wJ2nAiTosypqu2XnvEcDYEbRGDJGJPGiThWJyDB1PqeJVgCrg8+rgRuj2h/0iA1AppkVAkuBde5e7+4NwDpgWbAuw91fCp6f/GDUdw0JiQnGwok5PK8wEJFhqq9h4MDvzWyTma0K2grcvQYgeM8P2ouA/VHbVgdtZ2qv7qH9fcxslZlVmFlFXV1dH0vvH5eV5VLd0Mq+oy2D+nNFRAZDX8NgkbvPJXIK6HYz+9AZ+vZ0vt/Pof39je73uXu5u5fn5eX1VnO/WhSMGzxbObghJCIyGPoUBu5+MHg/DPyayDn/2uAUD8H74aB7NVAStXkxcLCX9uIe2oeUibnpFGel8ewuhYGIDD+9hoGZpZvZ6FOfgSXAG8Aa4NSMoJXAI8HnNcCtwayihUBjcBppLbDEzLKCgeMlwNpgXZOZLQxmEd0a9V1DhpmxeEoeL751VLe0FpFhpy9HBgXA82b2GvAy8Ft3fxy4G7jazCqBq4NlgN8Bu4Eq4AfA5wHcvR74BvBK8LoraAP4HPDDYJu3gMfOf9f63+IpeZw42cmr+xrCLkVEpF8l9dbB3XcDs3poPwpc2UO7A7ef5rseAB7oob0CmNGHekN1yaQckhKMZ3bVsXBiTtjliIj0G12BfBZGj0hm3gVZPLNT4wYiMrwoDM7S4gvz2F5znMPH28IuRUSk3ygMztLiKZEprc9W6gI0ERk+FAZnaVphBnmjUzXFVESGFYXBWTIzPlSWx3OVdXR193htnIhIzFEYnIPFF+bR0NLB1gONYZciItIvFAbn4IOTczFDs4pEZNhQGJyDrPQUZhVn8tTOw713FhGJAQqDc3T5hfm8Vn1MTz8TkWFBYXCOrpyajzs89aaODkQk9ikMztH0cRkUjhnB+h0KAxGJfQqDc2RmXHFRPs9V1nGysyvsckREzovC4DxcNbWA5vYuNuyu772ziMgQpjA4D5dMyiEtOZH1O2rDLkVE5LwoDM7DiOREFk3OZf2Ow0Tu3C0iEpsUBufpqqn5HDjWys7aprBLERE5ZwqD83TFRfkAmlUkIjGtz2FgZolmttnMHg2WJ5jZRjOrNLOfmVlK0J4aLFcF60ujvuNrQftOM1sa1b4saKsyszv6b/cGXn7GCGYVj+EJjRuISAw7myODLwE7opb/CfiOu5cBDcBtQfttQIO7Twa+E/TDzKYBNwPTgWXAvwUBkwj8K7AcmAbcEvSNGVdOLWDLfl2NLCKxq09hYGbFwLVEHlqPmRlwBfCLoMtq4Mbg84pgmWD9lUH/FcBD7n7S3fcAVcD84FXl7rvdvR14KOgbM05djaxZRSISq/p6ZPBd4K+B7mA5Bzjm7p3BcjVQFHwuAvYDBOsbg/7vtL9nm9O1v4+ZrTKzCjOrqKsbOncMnVaYQUl2Gmu3KQxEJDb1GgZmdh1w2N03RTf30NV7WXe27e9vdL/P3cvdvTwvL+8MVQ8uM2PZ9LE8X3mEpraOsMsRETlrfTkyWATcYGZvEzmFcwWRI4VMM0sK+hQDB4PP1UAJQLB+DFAf3f6ebU7XHlOWzRhLe1c3T+rGdSISg3oNA3f/mrsXu3spkQHgJ93948BTwEeDbiuBR4LPa4JlgvVPeuSKrDXAzcFsowlAGfAy8ApQFsxOSgl+xpp+2btBNKcki7zRqazddijsUkREzlpS711O62+Ah8zsm8Bm4P6g/X7gx2ZWReSI4GYAd99mZg8D24FO4HZ37wIwsy8Aa4FE4AF333YedYUiIcFYOr2AX246QFtHFyOSE8MuSUSkzyxWb6NQXl7uFRUVYZfxLi9UHeHjP9zIf3xiHkunjw27HBGRdzGzTe5e3tM6XYHcj+ZPyCZzZDJr39CpIhGJLQqDfpScmMBVUwt4Ykct7Z3dvW8gIjJEKAz62bLpYzne1smG3UfDLkVEpM8UBv3ssrJc0lMSeUynikQkhigM+tmI5EQuvyiftdsO0dmlU0UiEhsUBgPg+lnjqG9u58W3dKpIRGKDwmAALJ6Sx+jUJH7zWsxdSC0icUphMABGJCeyZPpYHt92iJOdXWGXIyLSK4XBALl+ViFNbZ08s3Po3F1VROR0FAYDZNHkXLLTU/jN6zVhlyIi0iuFwQBJTkxg+YyxPLG9lpb2zt43EBEJkcJgAF0/axytHV2s36HbWovI0KYwGEAfKM2mICNVs4pEZMhTGAygxATj2pnjeHpnHY2tegKaiAxdCoMBdv2sQtq7unn8DQ0ki8jQpTAYYLNLMpmYm84vXz0QdikiIqfVaxiY2Qgze9nMXjOzbWb29aB9gpltNLNKM/tZ8MhKgsda/szMqoL1pVHf9bWgfaeZLY1qXxa0VZnZHf2/m+ExMz4yt4iX99Szv74l7HJERHrUlyODk8AV7j4LmA0sM7OFwD8B33H3MqABuC3ofxvQ4O6Tge8E/TCzaUQegTkdWAb8m5klmlki8K/AcmAacEvQd9j48NxizOBXOjoQkSGq1zDwiBPBYnLwcuAK4BdB+2rgxuDzimCZYP2VZmZB+0PuftLd9wBVwPzgVeXuu929HXgo6DtsFGWmccnEHH61uZpYfcyoiAxvfRozCP6C3wIcBtYBbwHH3P3U1VTVQFHwuQjYDxCsbwRyotvfs83p2nuqY5WZVZhZRV1dbN3m4SNzi9l7tIWKvQ1hlyIi8j59CgN373L32UAxkb/kp/bULXi306w72/ae6rjP3cvdvTwvL6/3woeQ5TPGMjIlkV+9Wh12KSIi73NWs4nc/RjwNLAQyDSzpGBVMXDqyqpqoAQgWD8GqI9uf882p2sfVtJTk1g2YyyPvlZDW4fuZCoiQ0tfZhPlmVlm8DkNuArYATwFfDTothJ4JPi8JlgmWP+kR06UrwFuDmYbTQDKgJeBV4CyYHZSCpFB5jX9sXNDzR/PLabpZCe/314bdikiIu+S1HsXCoHVwayfBOBhd3/UzLYDD5nZN4HNwP1B//uBH5tZFZEjgpsB3H2bmT0MbAc6gdvdvQvAzL4ArAUSgQfcfVu/7eEQcsnEHMaNGcHPK/Zzw6xxYZcjIvIOi9XZLeXl5V5RURF2GWftO+t2cc+TlTz71cspyR4ZdjkiEkfMbJO7l/e0TlcgD7KbPlCCAT97ZX+vfUVEBovCYJAVZaaxeEoeD1fsp7OrO+xyREQAhUEobpk/nsNNJ3nyTT3nQESGBoVBCK64KJ/80ak8pFNFIjJEKAxCkJSYwE3lJTy98zAHj7WGXY6IiMIgLH/ygRIceLhCRwciEj6FQUhKskdy2eRcfvaKBpJFJHwKgxB9fMEF1DS28cQODSSLSLgUBiG6amo+RZlp/OjFPWGXIiJxTmEQoqTEBD5xyQVs2F3Pm4eOh12OiMQxhUHI/qS8hNSkBFa/uDfsUkQkjikMQpaVnsKNs4v49eZqjrW0h12OiMQphcEQsPLSUto6ujXNVERCozAYAqaNy2D+hGwefGkvXd2xeRdZEYltCoMh4pOXllLd0Mo6PfhGREKgMBgirp5WQHFWGj94bnfYpYhIHFIYDBFJiQl8+oMT2bS3gYq368MuR0TiTF+egVxiZk+Z2Q4z22ZmXwras81snZlVBu9ZQbuZ2T1mVmVmr5vZ3KjvWhn0rzSzlVHt88xsa7DNPWZmA7GzQ93HyovJGpnMvz+jowMRGVx9OTLoBL7i7lOBhcDtZjYNuANY7+5lwPpgGWA5kYfdlwGrgHshEh7AncACYD5w56kACfqsitpu2fnvWuwZmZLErZeU8sSOWqoON4VdjojEkV7DwN1r3P3V4HMTsAMoAlYAq4Nuq4Ebg88rgAc9YgOQaWaFwFJgnbvXu3sDsA5YFqzLcPeXPPJA5gejvivu3HrJBYxITuC+Z3V0ICKD56zGDMysFJgDbAQK3L0GIoEB5AfdioDoCfPVQduZ2qt7aO/p568yswozq6irqzub0mNGzqhUbiov4debD1B7vC3sckQkTvQ5DMxsFPBL4MvufqYb6fR0vt/Pof39je73uXu5u5fn5eX1VnLM+ovLJtLV7TzwvG5gJyKDo09hYGbJRILgJ+7+q6C5NjjFQ/B+6j7M1UBJ1ObFwMFe2ot7aI9b43NGct3F4/jxhr3UN+sWFSIy8Poym8iA+4Ed7v7PUavWAKdmBK0EHolqvzWYVbQQaAxOI60FlphZVjBwvARYG6xrMrOFwc+6Neq74tZfXjGZ1o4uXXcgIoOiL0cGi4BPAFeY2ZbgdQ1wN3C1mVUCVwfLAL8DdgNVwA+AzwO4ez3wDeCV4HVX0AbwOeCHwTZvAY/1w77FtLKC0Vw7s5AHX3ybBh0diMgAs8gEnthTXl7uFRUVYZcxoHbVNrH0u8/y+T+axFeXXhR2OSIS48xsk7uX97ROVyAPYVMKRnPNzEJ+9IKODkRkYCkMhrgvXlFGc3sX92tmkYgMIIXBEHfh2MjYwX++sIcjJ06GXY6IDFMKgxjwV0um0NbZzb88WRV2KSIyTCkMYsCkvFHcVF7CTzbuZd/RlrDLEZFhSGEQI758VRmJCca31+0MuxQRGYYUBjGiIGMEn1o0gUe2HOSNA41hlyMiw4zCIIZ8ZvEkMkcm83/X6uhARPqXwiCGjElL5guXT+bZXXU8vfNw7xuIiPSRwiDG3HpJKRNy07nr0e20d3aHXY6IDBMKgxiTkpTA/7luGrvrmln94tthlyMiw4TCIAZdflE+l1+Yxz3rK6lr0oVoInL+FAYx6n9fN422zi6+tfbNsEsRkWFAYRCjJuaN4lOLJvBwRTWb9zWEXY6IxDiFQQz7yyvLGJsxgq/9aisdXRpMFpFzpzCIYaNSk7hrxXTePNSkJ6KJyHnpy2MvHzCzw2b2RlRbtpmtM7PK4D0raDczu8fMqszsdTObG7XNyqB/pZmtjGqfZ2Zbg23uCR59KX20ZPpYls8Yy/eeqOTtI81hlyMiMaovRwY/Apa9p+0OYL27lwHrg2WA5UBZ8FoF3AuR8ADuBBYA84E7TwVI0GdV1Hbv/VnSi7+/YTopiQn83f9sJVafXCci4eo1DNz9WaD+Pc0rgNXB59XAjVHtD3rEBiDTzAqBpcA6d6939wZgHbAsWJfh7i955LfYg1HfJX1UkDGCv1l+ES9UHeXnFdVhlyMiMehcxwwK3L0GIHjPD9qLgP1R/aqDtjO1V/fQ3iMzW2VmFWZWUVdXd46lD09/On88CyZkc9ej29lfr9tci8jZ6e8B5J7O9/s5tPfI3e9z93J3L8/LyzvHEoenhATj2zfNAuArP3+Nrm6dLhKRvjvXMKgNTvEQvJ+6a1o1UBLVrxg42Et7cQ/tcg6Ks0Zy5/XTeHlPPfc/r9lFIsNNU1vHgB35J53jdmuAlcDdwfsjUe1fMLOHiAwWN7p7jZmtBf4xatB4CfA1d683syYzWwhsBG4Fvn+ONQnw0XnFPLGjlv+3dhcfLMtjamFG2CWJDJjubqe9q5uOrm46upyOrm7aO9+93NOcip7mLJqBYSQkBO8WtJmRYIZB5P1d7X/oi/GufgmRLwy+5w99u91p7eiitb2Lto4u2jq6ae3ooqW9k+NtnRxraae+uZ2G5nYaWjqob27n0PE2DjW2ceJkJwUZqWz826v6/X/LXsPAzH4K/BGQa2bVRGYF3Q08bGa3AfuAjwXdfwdcA1QBLcAnAYJf+t8AXgn63eXupwalP0dkxlIa8FjwknNkZvzjh2ey9LvP8cWfbuaRLyxiZMq5Zr7IwHN3Gls7OHIi8kuwvvkkR5vbqT/RztHmdhpa2jnR1smJk500t3fSfLIr8vlkJy3tXWGXP2AyRiSRlZ5C1sgUJuWlc9nkXMZljmBcZtqA/DyL1amI5eXlXlFREXYZQ9ZzlXXc+sDLfHhOEd/+2Cx0+YaEpbvbqTnexp66Zg4ca+HAsTZqjrVS09jGwWOtHGxspa2j5yvoR6dGfiGOSk1iVGoS6amJpL/zOfIakZxASmICye+8jJSkyOekhOAv9Cg9/cZzd/zUu0O3g+OR93faot5P0xf/wzbdQV+I3jZypDAyJZHU5ETSTr1SEhmRnMCYtGQyR6aQmZZMUmL/XxNsZpvcvbyndfqTcZj6YFkeX7yijO+tr2ThhBxu+kBJ7xuJnIeW9k52Hmqi6vAJ9hxpftfrZNSzN8wgf3Qq4zLTmFqYwZVT8xk7Jo3cUSlkp0deOempZKUnk5qUGOIexReFwTD2xSvLqNhbz/9+5A1mFo/R+IH0C3fnwLFWdtQ08WbNcXYcOs6OmibePtr8zvn5pARjfM5IJuaO4kNT8piQm05pTjol2WkUZIwgeQD+6pXzo9NEw1xd00muvec50lIS+Z/PLyIrPSXskiTGHG/r4PX9jWze18CW/cfYsv8YR5vb31l/Qc5Ipo7NYGphBhcVjmZKwWhKstIG5DSHnB+dJopjeaNTuffP5nHLfRu4/b9fZfWn5uuvMjmj/fUtbNxTz8bdR9m8/xhv1Z145y/+yfmjuPyifGaVZDKtMIOLxo4mPVW/RoYD/SvGgXkXZHH3H8/krx5+ja//ZhvfvHFm2CXJEHLql/+G3UfZsPso1Q2tAGSNTGbO+CxWzBrH7PGZXFycyZi05JCrlYGiMIgTH5lbzK7aE/z7M29Rlj+alZeWhl2ShOR4WwcvVh3hmV11PFd55F2//BdOzOHTH5zIwok5lOWPIiFBs9DihcIgjnx16YVUHT7B3/9mG3mjU7lmZmHYJckg6O52ttcc55lddTyzs45N+xro6nZGpyZx6WT98pcIhUEcSUwwvn/LHP7s/o18+aEtZKYlc+nk3LDLkgHQ0t7Js7vq+P32Wp7dVceRE5EB3xlFGXxu8SQWX5jH7JJMjR/JOzSbKA41tnRw03+8RHVDCw+tuoSZxWPCLkn6wZETJ1m/o5bfb6vluaojtHd2kzkymcVT8lg8JY8PluWRNzo17DIlRGeaTaQwiFOHGtv443tfpKW9k5/8xUKmjdM1CLHo7SPN/H77IX6/rZZN+xpwh6LMNJZML2DJtLF8oDRLUzzlHQoD6dHeo83cfN8GWju6+K/bFjCjSEcIsWB33Qkefb2G375ew87aJgCmFWa8EwBTC0fr9iPSI4WBnNa+oy3c8oMNnDjZyX/dtkCnjIaovUebefT1Gh59vYYdNccxgw9ckM2yGWO5eloBJdkjwy5RYoDCQM5of30LN9+3geOtHfzHrfO4dJIGlYeC/fUt/HZr5Ahg64FGAOaOz+S6i8dxzcxCxo4ZEXKFEmsUBtKrg8da+fP/fJk9R5r59k2zuWHWuLBLiksHj7Xy29dreHRrDa/tPwbArJJMrptZyDUXF1I0QLcvlvig21FIr8ZlpvHzz1zKp39cwRd/upmDx1r5zIcm6tzzIKg93hYJgNcP8uq+SADMLBrDHcsv4tqZhToFJINCYSDvGDMymQc/NZ+v/Pw17n7sTXbUHOfuj1xMWopuI9zfGprbeeyNQ/zmtYNs2HMUd5hamMFXl17ItTMLKc1ND7tEiTMKA3mXEcmJfP/mOUwdO5pvr9vFrtoT3PeJefrrtB80tXWwbnsta147yPOVR+jsdibmpfOlK8u47uJxTM4fFXaJEseGzJiBmS0DvgckAj9097vP1F9jBgPvqZ2H+dJPNwPwDx+eyfUaRzhrbR1dPPnmYdZsOciTOw/T3tlNUWYa188ax/WzCplWmKFTcTJohvwAspklAruAq4FqIs9KvsXdt59uG4XB4Nh7tJkvPbSFLfuP8eE5RXx9xXQyRujOlWfS3tnNc5V1/Oa1g6zbXktzexd5o1O5dmYh188ax9zxmQoACUUsDCDPB6rcfTeAmT0ErABOGwYyOC7ISecXn72Ef3mqiu8/WcXG3Uf5+ooZXD2tIOzShpT2zm5eeOsIj289xOPbDtHY2sGYtGRumD2O6y8ex4KJOSTqJnAyhA2VMCgC9kctVwML3tvJzFYBqwDGjx8/OJUJSYkJfPmqKSyekscdv9zKpx+s4OppBfz9DdPjeqpjW0cXz+yq4/E3DuxEV3YAAAc9SURBVPHEjlqa2joZlZrEVVPzuWH2OC6bnEdKkm4FIbFhqIRBT38yve/8lbvfB9wHkdNEA12UvNuc8Vk8+sXLeOD5PXz3iUqu/PbTfHLRBD67eFLcPPTkxMlOnnrzMI+/cYindh6mpb2LMWnJLJ0+luUzxnJZWa4e4i4xaaiEQTVQErVcDBwMqRY5g+TEBD6zeBLXXlzIt9bu5N6n3+K/N+7js4sn8YlLLmDUMHwE4r6jLTz5Zi3r3zzMxt31tHd1kzsqhRvnFLF8xlgWTszRraAl5g2VAeQkIgPIVwIHiAwg/6m7bzvdNhpAHhq2HWzkW2t38vTOOkaPSOLjCy7gzy8tjelbJXR0dVPxdgNP7TzM+h21vFXXDMDEvHSuuDCfq6cVUF6arTEAiTlDfjYRgJldA3yXyNTSB9z9H87UX2EwtGze18APn9vDY2/UkJhgXD2tgI/OK+ZDZXlD/hbKXd3OjprjvPjWEV6oOsorb9fT0t5FcqKxcGIOl1+YzxUX5etCMIl5MREGZ0thMDTtO9rC6pfe5tebD1Df3P7OlMqrphYwf0L2kBhQbW3vYuuBRjbva2DT3gY27qmnsbUDgMn5o7h0Ug6XTsrlsrLcYXnaS+KXwkAGXXtnN0/tPMwvN1XzzK46TnZ2Mzo1iUWTc5k/IZvy0iymFWYM+FFDY2sHlbVN7KxtYkfNcTbvO8abh5ro6o78dz8+eyQLJ2Zz6aRcLpmUQ0FG7J7eEumNwkBC1drexfNVR1i/o5bnKo9w4FgrAGnJiUwpGMXk/NGUFYyiNCedgoxUCjJGkDc6tddB2e5u53hbB8daOjjafJLqhtaoVwtVh09Q09j2Tv/RqUlcXDKGOSVZzBmfyeySTHJG6TGQEj8UBjKk1DS2UvF2A5v3HWNXbRO7aps43HTyff1GJCcwMiWJtOREUpMS6HKns8vp6nbaOrtobO2gp/98s9NTKMpMY3L+KKYUjObCsZH3osw0XfkrcS0WrkCWOFI4Jo3rZ6W9615HjS0d7G9oofZ4G7XHT3K4qY3mk520dnTR0t5Fe2c3yYkJJCYYSQlGalICY0amkJmWTObIZLLSUyjOTKMoK42RKfrPWuRs6f81MiSMGZnMmJFj9BxmkZCEP7VDRERCpzAQERGFgYiIKAxERASFgYiIoDAQEREUBiIigsJARESI4dtRmFkdsDfsOs5SLnAk7CIGmfY5PmifY8MF7p7X04qYDYNYZGYVp7svyHClfY4P2ufYp9NEIiKiMBAREYXBYLsv7AJCoH2OD9rnGKcxAxER0ZGBiIgoDEREBIXBoDKz/2Vmbma5wbKZ2T1mVmVmr5vZ3LBr7C9m9i0zezPYr1+bWWbUuq8F+7zTzJaGWWd/MrNlwT5VmdkdYdczEMysxMyeMrMdZrbNzL4UtGeb2Tozqwzes8Kutb+ZWaKZbTazR4PlCWa2Mdjnn5lZStg1ng+FwSAxsxLgamBfVPNyoCx4rQLuDaG0gbIOmOHuFwO7gK8BmNk04GZgOrAM+DczSwytyn4S7MO/Evk3nQbcEuzrcNMJfMXdpwILgduD/bwDWO/uZcD6YHm4+RKwI2r5n4DvBPvcANwWSlX9RGEweL4D/DUQPWK/AnjQIzYAmWZWGEp1/czdf+/uncHiBqA4+LwCeMjdT7r7HqAKmB9Gjf1sPlDl7rvdvR14iMi+DivuXuPurwafm4j8ciwisq+rg26rgRvDqXBgmFkxcC3ww2DZgCuAXwRdYn6fFQaDwMxuAA64+2vvWVUE7I9arg7ahptPAY8Fn4frPg/X/TotMysF5gAbgQJ3r4FIYAD54VU2IL5L5I+57mA5BzgW9QdPzP97J4VdwHBhZk8AY3tY9XfA3wJLetqsh7aYmet7pn1290eCPn9H5NTCT05t1kP/mNnnMxiu+9UjMxsF/BL4srsfj/yhPDyZ2XXAYXffZGZ/dKq5h64x/e+tMOgn7n5VT+1mNhOYALwW/B+mGHjVzOYT+WuiJKp7MXBwgEvtN6fb51PMbCVwHXCl/+GClpje5zMYrvv1PmaWTCQIfuLuvwqaa82s0N1rglOdh8OrsN8tAm4ws2uAEUAGkSOFTDNLCo4OYv7fW6eJBpi7b3X3fHcvdfdSIr805rr7IWANcGswq2gh0HjqUDvWmdky4G+AG9y9JWrVGuBmM0s1swlEBs9fDqPGfvYKUBbMMEkhMki+JuSa+l1wrvx+YIe7/3PUqjXAyuDzSuCRwa5toLj719y9OPj/783Ak+7+ceAp4KNBt5jfZx0ZhOt3wDVEBlFbgE+GW06/+hcgFVgXHBFtcPfPuvs2M3sY2E7k9NHt7t4VYp39wt07zewLwFogEXjA3beFXNZAWAR8AthqZluCtr8F7gYeNrPbiMyY+1hI9Q2mvwEeMrNvApuJhGTM0u0oREREp4lERERhICIiKAxERASFgYiIoDAQEREUBiIigsJARESA/w87/ufJupBUfgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let's plot the function\n",
    "xmax = 50\n",
    "x = np.linspace(-xmax,xmax,2001)\n",
    "f,f1,f2 = myfunction(x)\n",
    "\n",
    "plt.plot(x,f)\n",
    "#plt.axis([-xmax, xmax, -20, 6000])\n",
    "\n",
    "imin = np.argmin(f)\n",
    "xopt = x[imin]\n",
    "print('Minimizer xopt')\n",
    "print(xopt)\n",
    "print('Value of f at xopt')\n",
    "print(f[imin])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of iterations = \n",
      "37\n",
      "Final answer xhat = \n",
      "-13.746092693010832\n",
      "Value of f at xopt\n",
      "-3461.191330045192\n",
      "Value of first derivative at xopt\n",
      "8.26144571419718e-07\n",
      "Value of second derivative at xopt\n",
      "53.97745744587804\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEWCAYAAAAzcgPFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxV9Z3/8dcnK2EJWQhrAgEBBUREIyLaWnfcnU7r6DhKW6d0sdP2N51pbTvzs7XLOO2vm+3U1qqjdGyt3UZqq4i4V0GD4IKIBBAIayBhDSHb5/fH+Uav4WaB5ObeJO/n43Ef957v+Z5zPt97zz2fs3zvuebuiIiIpJq0ZAcgIiISjxKUiIikJCUoERFJSUpQIiKSkpSgREQkJSlBiYhISuqzCcrMPmJmz7Uz/ikz+8eejKm7mNlXzOyuZMeRCGZ2r5l9s5N1x5rZATNLT3Rc7cRwnZk9loTl/oeZfb4T9d42s/MTGEfSP4PeysxKzczNLKOT9X9mZv+e6Lg6iGGVmX2gG+bzfTP7ZEf1UjpBmdk1ZrbMzA6a2c7w+tNmZsmOraeY2QfMrDK2zN2/7e69Mrl2J3ff5O6D3b0JEr/TEW+D4u73u/uFiVpmG3EUATcAP+/J5cbT+jPoSWZ2p5mtMbNmM/tIF+ZzVIkiWdz9k+7+DYi/Xehu8XYW3X2auz/VDbP/LvBVM8tqr1LKJigz+wLwI6KGjARGAJ8EzgTabZSkrlTeCPSio4CPAH9x90PJDiTJXgE+Dbyc7EB6m2R/D919G/AmcEVHFVPuAQwFDgJ/24l6C4AqYCPwb0BaGPcR4LmYuheEN2Qv8BPgaeAfY8Z/DFgN1ACLgHEx45woOa4N4/8LsDZimgWUA/uAHcD3Y8bNBp4H9hB9uT4QM64A+G9ga1jG/wKDgENAM3AgPEYDXwP+J2baK4BVYb5PAVNixr0N/Avwamj7b4ABbcQ+Mbwve4FdwG9ixp0ALAaqgTXA1THjLgVWhDZvBr4WM640vH83ApuAZ0L5WTHvxWbgI6H83vD+/hnYDywDjmsj3pZ5ZwDfApqAuvA+/aQTcd8L3AH8hWh9O7+DtmwKy2v5LM7gyPVsDvBSeA9fAubEjHsK+Abw19C2x4BhYdwA4H+A3eE9eQkY0Ua7nwD+IWZ4GPBwmK4aeJZ3vwftfv7Ax4GKMN1CYHQo/zrw4/A6M7w/3wnDOeF9zo/9DDpqYxh/A9F3dTfw7yG+89toZ5ufRat6zxHWn3a2FUf1ucaZPg24GVgXYn8QKIgZ/1tge3iPnwGmxYzLAb4X2r03xJsT897NCzHsAr7aThvuBb5J29uFNmOk7e9h3LiB+UADUB/m/6eY9en88Dob+CHRNmtreJ0dxn0AqAS+AOwEtgEfbdWerwL/3e7n1t7IZD2AuUBjy0rfTr0FwEPAkPABvAXcGMZ9hLDhIPoC7wM+RPRl+z9h/v8Yxl9F9CWdQrSx+zfg+ZjlONEGIA8YS5QQ57YR0wvA9eH1YGB2eD0mrDSXhBXpgjBcFMb/mWjjkR9iPDv2g261jK8REhQwmWjjcUGY7ouhLVkxK9SLYQUuIErCn2wj9l+HlSaNaIN5VigfRPSl/mh4f04h+jJNi4lxepjuJKLEfFWrL8aCMJ+c8B7uB64NMRcCJ8d8CauJEn0GcD/wQBvxtsw7duMYu9PRUdz3En0xz4xpc2fakhGzjI/w7npWQLRzcX1Y3rVhuDAmvnXhM8sJw7eFcZ8A/gQMBNKBU4HcNtpdBZwWM/wfwM/Ce5kJvI+wA9Xe5w+cG96PU4g2Nj/m3Q3XucBr4fWcEPeymHGvtPMZtNXGqUQbu7OIzoL8P6KNYFsJqs3PolW9ziSoo/pc40z/eWApUBzeq58Dv44Z/zGi7VDLRntlzLj/Cu/DmPDZzgn1Wpb7i/BezQAOE7OD2SqGe4FvtrNdaDNG4nwPOxH3O8uLKXubdxPUrWF5w4Eioh3Ob8TE1xjqZBJt92qB/Jh5fRB4ud3PraNkkYwH8A/A9lZlLXvbh4D3hw/6MDA1ps4ngKfibDhuAJbG1DOi7N6SoB4hJLYwnBbezHFh2Akb6zD8IHBzG7E/Q7T3OaxV+ZeAX7YqW0S09zSKaG8oP8784q2IX+PdBPXvwIOtYt9CODoLK1Ts3vZ3gJ+1EfsC4E6guFX53wHPtir7OXBLG/P5IfCDVl+MCTHjvwz8sZ0v4V0xw5cAb7ZRt2XebSWoduMOy1rQwboYry1tJajrgRdbTf8C7x4dPgX8W8y4TwOPhtcfI1rHT+rE96MBOCFm+FaiHbWJceq2+fkDdxOOisLw4DDvUt49Siok2iv/CtF3ZjDR+n17O59BW238v7x3oz6QaA89boJq77NoVd5hgjrazzVO/dXAeTHDo8J7dcQ0RDuyTnSGJ41omzWjnfW3OKbsReCadr4b7SWoNmMkzvewvbhbL6/V+tSSoNYBl8SMuwh4Oya+Q7z3u7KTsMMehi8A1rf3OaXqNajdwLBWF6PnuHteGJdGdFSURXTY3GIj0V5Ka6OJ9qRb5uWxw8A44EdmtsfMWk6TWKt5bY95XUv0RY3nRqK9xzfN7CUzuyxmGR9uWUZYzllEK1EJUO3uNW3Msz2jiXkP3L05tO1YYv8iUbtfDL11PhYT++mtYr+O6NogZna6mT1pZlVmtpfodOiwVvOOfb9LiFbutnQ23o60G3ecuDrblra857MIWq+TbbXtl0Q7LA+Y2VYz+46ZZbaxnBqivd4W3yU6an7MzNab2c2t6re1zNbrzgGi79cYj65vlQNnE+0QPk2UQM8MZU+3EVtHy4v9HtaG5cXVxc+iu+c1DvhjzHq0muiU8ggzSzez28xsnZntI9qIE+Y/jOjIvKfW97gxxtR55/3vIO7OaL2+bwxlLXa7e2PMcOu2DSE66GhTqiaoF4iOjq5sp84uor2DcTFlY4mOHlrbRrRRBCD0AiyJGb8Z+IS758U8ctz9+aMN3N3Xuvu1RIe9/wn8zsxaTjX9stUyBrn7bWFcgZnlxZtlB4vcSsx7ENO2eO9DR7Fvd/ePu/tooqPRn5rZxBDf061iH+zunwqT/oro+kWJuw8lOt3UuqdlbDs2A8cdbXydaUKr4Y7ijjdNe205qs8iaGudfG/g7g3u/nV3n0p0CugyoiP/eF4l2glqmXa/u3/B3ScAlwP/bGbndbTM1vGG9bQwJt6niU7nzSS6JvY00V7yLKIzBUdrG9Hpp5bl5YTltaUz61VndeVzhWhdurjVujTA3bcAf0+0rTqf6KipNExjRNupOrp/fY8Xc3sxxpuuvbjbWkas1uv72FDWWVOIrsW3KSUTlLvvITqN8FMz+5CZDTazNDM7mej8KR51a30Q+JaZDTGzccA/E11obu3PwDQz+2A4Kvss792L/hnwZTObBmBmQ83sw8cSu5n9g5kVhSOZlr2DphDX5WZ2UdhzGRC6ihZ71KPlkdDefDPLNLP3h2l3AIVmNrSNRT4IXGpm54U97i8QJfejTq5m9mEza9mA1BCtoE1E198mm9n1IbZMMzvNzKaEukOIjgDrzGwW0YrfnvuB883sajPLMLPC8Nl21Q5gQsxwR3HH015bqohOxU6IO2XU2WKymf19aNffEV13ebijwM3sHDObHnoS7iPa+Wqr6/ZfiI5iWqa9zMwmhp2TfWG6znT7/hXwUTM72cyygW8TXWd6O4x/mihJvuHu9YRTqMAGd6/qxPxb+x3Rd2BO6F78ddpPOO2uV2aWZWYDwjwyw3eqrW1aVz5XiLYR3wrbGcysyMxadqCHEH3ndhOdtvx2y0RhO3AP8H0zGx2++2eE97sr4m0X2osxnjbjjllGe+/Jr4F/C8sZRnQKN972ty1nE2332pSSCQrA3b9DlHC+SHTucgfR9YMv8e7G95+IOgisJzoP/SuilaH1vHYBHwZuI/owJhH1MmoZ/0eio50HwqHu68DFxxj6XGCVmR0g6iZ/jbvXuftmor2VrxB9ITYD/8q7n8H1RBulN0N7Px9ie5NoRVgfDt1jD6Fx9zVE1+x+TLS3djlwedigHK3TgGUh9oXA59x9g7vvBy4EriHaQ9pO9H61fMk+DdxqZvuJVtIH21uIu28iurb0BaLTqSuJLhB31Y+AD5lZjZnd3om442mzLeGU1LeAv4bPYnardu0mOvL5AtF69kXgsrD+dWQk0QZ8H9Gpmadp+8u+ALgkHIFAtD4/TtQB4QXgp96J36q4+xKia5i/Jzq6OY7ovWrxPNG1qJajpTeIjgaO5egJd19F9J19ICxvP9G6friNSTparx4jus4xh+jaacv16aOaV0efa/Ajou/EY2EeS4HTw7gFRKe3thC9R0tbTfsvwGtER6HVROtgl7a9bWwX2osxno7ivhuYGub/v3Gm/ybRaeBXidr3cijrkJmNItp5izffd+uFi1Ui0ouY2beBne7+w2THcqzMbDDRWYZJ7r4h2fFIzzGz7wHr3P2n7dZTghKRnmJmlwNLiE7LfY9oD/8U14ZI4kjZU3wi0iddybs/7JxEdApcyUni0hGUiIikJB1BiYhISkrZG3cmyrBhw7y0tDTZYYiI9BrLly/f5e5FPb3cfpegSktLKS8vT3YYIiK9hpm1vkNKj9ApPhERSUlKUCIikpKUoEREJCUpQYmISEpSghIRkZSkBCUiIilJCUpERFKSElQn1Dc2c8dT63h27bH8BY6IiBwLJahOyEw3fvHsehauPJo/ixQRka5QguoEM+PkkjxWbt7TcWUREekWSlCdNLMkj4qqA+yra0h2KCIi/YISVCedPDYPd3h1895khyIi0i8oQXXSjJI8zGDFpppkhyIi0i8oQXVS7oBMJhYN1nUoEZEeogR1FE4uyWPF5j3oX4hFRBIvYQnKzI43s5Uxj31m9nkzKzCzxWa2Njznh/pmZrebWYWZvWpmp8TMa16ov9bM5sWUn2pmr4VpbjczS1R7AGaOzaf6YD2bqmsTuRgRESGBCcrd17j7ye5+MnAqUAv8EbgZWOLuk4AlYRjgYmBSeMwH7gAwswLgFuB0YBZwS0tSC3Xmx0w3N1HtgegICtBpPhGRHtBTp/jOA9a5+0bgSuC+UH4fcFV4fSWwwCNLgTwzGwVcBCx292p3rwEWA3PDuFx3f8Gjc24LYuaVEJNHDGZgVjorNilBiYgkWk8lqGuAX4fXI9x9G0B4Hh7KxwCbY6apDGXtlVfGKT+Cmc03s3IzK6+qOvbbFWWkpzF9zFBW6AhKRCThEp6gzCwLuAL4bUdV45T5MZQfWeh+p7uXuXtZUVFRB2G0b+bYfN7Yupe6hqYuzUdERNrXE0dQFwMvu/uOMLwjnJ4jPO8M5ZVAScx0xcDWDsqL45Qn1MyxeTQ0OW9s25foRYmI9Gs9kaCu5d3TewALgZaeePOAh2LKbwi9+WYDe8MpwEXAhWaWHzpHXAgsCuP2m9ns0Hvvhph5JczM0FFC16FERBIrI5EzN7OBwAXAJ2KKbwMeNLMbgU3Ah0P5X4BLgAqiHn8fBXD3ajP7BvBSqHeru1eH158C7gVygEfCI6GG5w5gTF5OuKPE+EQvTkSk30pognL3WqCwVdluol59res6cFMb87kHuCdOeTlwYrcEexR0Z3MRkcTTnSSOwcyxeVTWHKJq/+FkhyIi0mcpQR2DmWNbrkPpxrEiIomiBHUMpo0eSlZ6GsuVoEREEkYJ6hgMyEznxDG5lL+tBCUikihKUMeorLSA1yr1g10RkURRgjpGZePyqW9q5vUt+oddEZFEUII6RqeOi26oXr5Rp/lERBJBCeoYFQ7OZsKwQZS/Xd1xZREROWpKUF1w6rh8lm+s0T/siogkgBJUF5SV5lNT28C6qoPJDkVEpM9RguqCstICAJZv1Gk+EZHupgTVBROGDaJgUBYv6fdQIiLdTgmqC8yMU8ZG16FERKR7KUF1UVlpPht2HWTXAd04VkSkOylBdVFZ+D2UjqJERLqXElQXTS8eSlZGmn4PJSLSzRKaoMwsz8x+Z2ZvmtlqMzvDzArMbLGZrQ3P+aGumdntZlZhZq+a2Skx85kX6q81s3kx5aea2WthmtvDX7/3qOyMdE4aM1R3lBAR6WaJPoL6EfCou58AzABWAzcDS9x9ErAkDANcDEwKj/nAHQBmVgDcApwOzAJuaUlqoc78mOnmJrg9cZ1ams/rW/ZyqF43jhUR6S4JS1Bmlgu8H7gbwN3r3X0PcCVwX6h2H3BVeH0lsMAjS4E8MxsFXAQsdvdqd68BFgNzw7hcd38h/F38gph59ajTxxfQ0OT6A0MRkW6UyCOoCUAV8N9mtsLM7jKzQcAId98GEJ6Hh/pjgM0x01eGsvbKK+OUH8HM5ptZuZmVV1VVdb1lrZSVFpBmsHSDrkOJiHSXRCaoDOAU4A53nwkc5N3TefHEu37kx1B+ZKH7ne5e5u5lRUVF7Ud9DHIHZDJ1dC7L1u/u9nmLiPRXiUxQlUCluy8Lw78jSlg7wuk5wvPOmPolMdMXA1s7KC+OU54Us8cXsmLzHv2BoYhIN0lYgnL37cBmMzs+FJ0HvAEsBFp64s0DHgqvFwI3hN58s4G94RTgIuBCM8sPnSMuBBaFcfvNbHbovXdDzLx63OkTCqlvbOaVzXuSFYKISJ+SkeD5/xNwv5llAeuBjxIlxQfN7EZgE/DhUPcvwCVABVAb6uLu1Wb2DeClUO9Wd2+52PMp4F4gB3gkPJJiVmkBZrBsQzWnTyhMVhgiIn1GQhOUu68EyuKMOi9OXQduamM+9wD3xCkvB07sYpjdYujATE4YmcuyDbuJeryLiEhX6E4S3ej08QUs31hDfWNzskMREen1lKC60ewJBdQ1NPNqpa5DiYh0lRJUN5o1Prr2tEy/hxIR6TIlqG5UMCiLySMGs1S/hxIR6TIlqG52+vhClm+soaFJ16FERLpCCaqbzZ5QSG19E69v2ZvsUEREejUlqG42a3wBoOtQIiJdpQTVzYqGZDNx+GBeWKfrUCIiXaEElQBzjivkxQ3V+j2UiEgXKEElwJkTh3GooYmVui+fiMgxU4JKgNkTCkkz+GvFrmSHIiLSaylBJcDQnEymjxmqBCUi0gVKUAly5sRhrNy8hwOHG5MdiohIr6QElSBnThxGY7Pz4gb15hMRORZKUAly6rh8sjLS+GuFEpSIyLFQgkqQAZnpnFaar+tQIiLHSAkqgeYcN4w3t+9n14HDyQ5FRKTXSWiCMrO3zew1M1tpZuWhrMDMFpvZ2vCcH8rNzG43swoze9XMTomZz7xQf62ZzYspPzXMvyJMa4lsz9E6a+IwAJ7XXSVERI5aTxxBnePuJ7t7y1+/3wwscfdJwJIwDHAx0X+lTwLmA3dAlNCAW4DTgVnALS1JLdSZHzPd3MQ3p/NOHDOU3AEZPK/TfCIiRy0Zp/iuBO4Lr+8DroopX+CRpUCemY0CLgIWu3u1u9cAi4G5YVyuu7/g7g4siJlXSkhPM2ZPKOQ5JSgRkaOW6ATlwGNmttzM5oeyEe6+DSA8Dw/lY4DNMdNWhrL2yivjlB/BzOabWbmZlVdVVXWxSUfnrEnDqKw5xKbdtT26XBGR3i7RCepMdz+F6PTdTWb2/nbqxrt+5MdQfmSh+53uXubuZUVFRR3F3K3ODNehnlnbs4lRRKS3S2iCcvet4Xkn8Eeia0g7wuk5wvPOUL0SKImZvBjY2kF5cZzylDJh2CCK83N45i0lKBGRo5GwBGVmg8xsSMtr4ELgdWAh0NITbx7wUHi9ELgh9OabDewNpwAXAReaWX7oHHEhsCiM229ms0PvvRti5pUyzIyzJxfx/Lrd+vsNEZGjkMgjqBHAc2b2CvAi8Gd3fxS4DbjAzNYCF4RhgL8A64EK4BfApwHcvRr4BvBSeNwaygA+BdwVplkHPJLA9hyzsycXceBwIy9vqkl2KCIivUZGombs7uuBGXHKdwPnxSl34KY25nUPcE+c8nLgxC4Hm2BnHFdIRprx9FtVzJ5QmOxwRER6Bd1JogcMGZDJqePyeXqNrkOJiHSWElQPOfv4It7Yto+d++qSHYqISK+gBNVDzp4cdW9/Zq1+tCsi0hlKUD1k6qhcioZkq7u5iEgnKUH1EDPj/ZOKeHZtFU3NcX9PLCIiMZSgetDZxxdRU9vAa1v2JjsUEZGUpwTVg943cRhmqDefiEgnKEH1oPxBWcwozuPJNTs7riwi0s8pQfWwc44fziuVe/QvuyIiHVCC6mHnTRmOOzz5po6iRETaowTVw6aNzmXU0AEsWa0EJSLSHiWoHmZmnHvCcJ5dW8XhxqZkhyMikrKUoJLg/CkjOFjfxNL11R1XFhHpp5SgkuCM4wrJyUxnyeodyQ5FRCRlKUElwYDMdM6cOIwlq3cS/cuIiIi0pgSVJOdPGc6WPYdYs2N/skMREUlJSlBJcu4JwwHUm09EpA0JT1Bmlm5mK8zs4TA83syWmdlaM/uNmWWF8uwwXBHGl8bM48uhfI2ZXRRTPjeUVZjZzYluS3canjuAGcVDeVzXoURE4uqJI6jPAatjhv8T+IG7TwJqgBtD+Y1AjbtPBH4Q6mFmU4FrgGnAXOCnIemlA/8FXAxMBa4NdXuN86aMYOVm3VVCRCSehCYoMysGLgXuCsMGnAv8LlS5D7gqvL4yDBPGnxfqXwk84O6H3X0DUAHMCo8Kd1/v7vXAA6Fur9FyVwn15hMROVKij6B+CHwRaA7DhcAed28Mw5XAmPB6DLAZIIzfG+q/U95qmrbKj2Bm882s3MzKq6pS507iU0flUlKQw6JVSlAiIq0lLEGZ2WXATndfHlscp6p3MO5oy48sdL/T3cvcvayoqKidqHuWmTF32kieW7uL/XUNyQ5HRCSlJPII6kzgCjN7m+j027lER1R5ZpYR6hQDW8PrSqAEIIwfClTHlreapq3yXmXuiSOpb2rmCd08VkTkPRKWoNz9y+5e7O6lRJ0cnnD364AngQ+FavOAh8LrhWGYMP4Jj37FuhC4JvTyGw9MAl4EXgImhV6BWWEZCxPVnkSZWZJP0ZBsFq3anuxQRERSSkbHVbrdl4AHzOybwArg7lB+N/BLM6sgOnK6BsDdV5nZg8AbQCNwk7s3AZjZZ4BFQDpwj7uv6tGWdIO0NOOiaSP4/fIt1DU0MSAzPdkhiYikBOtvt9opKyvz8vLyZIfxHn+t2MV1dy3j59efykXTRiY7HBGR9zCz5e5e1tPL1Z0kUsCs8QXkDcxk0es6zSci0kIJKgVkpqdx/pQRPL56B/WNzR1PICLSDyhBpYi500ayr66Rpet3JzsUEZGUoASVIs6aNIxBWek8otN8IiKAElTKGJCZzjknDGfRqu00Nuk0n4iIElQKuXzGaKoP1vP8Op3mExFpM0GZ2S/D8+d6Lpz+7ezJRQzJzuBPr/S6G2KIiHS79o6gTjWzccDHzCzfzApiHz0VYH8yIDOdC6eN5NFV2znc2JTscEREkqq9BPUz4FHgBGB5q0dq/dK1D7l8xij21zXy9JrUueu6iEgytJmg3P12d59CdAuhCe4+PuYxoQdj7FfOnDiMgkFZ/OnVbckORUQkqTrsJOHun+qJQCSSmZ7GxSeO5PE3dlBb39jxBCIifZR68aWgy2eM5lBDE0tW6y84RKT/UoJKQaeVFjAiN1u9+USkX1OCSkHpacal00fz1Joq9h7SP+2KSP+kBJWiLp8xivqmZh59XZ0lRKR/UoJKUSeX5DFh2CB+//KWZIciIpIUCUtQZjbAzF40s1fMbJWZfT2UjzezZWa21sx+E/6unfCX7r8xs4owvjRmXl8O5WvM7KKY8rmhrMLMbk5UW5LBzPjgKWN4cUM1m6trkx2OiEiPS+QR1GHgXHefAZwMzDWz2cB/Aj9w90lADXBjqH8jUOPuE4EfhHqY2VSiv3+fBswFfmpm6WaWDvwXcDEwFbg21O0z/uaUYszgDzqKEpF+KGEJyiMHwmBmeDhwLvC7UH4fcFV4fWUYJow/z8wslD/g7ofdfQNQAcwKjwp3X+/u9cADoW6fMSYvhzMmFPKHFZW4e7LDERHpUQm9BhWOdFYCO4HFwDpgj7u3/AK1EhgTXo8BNgOE8XuBwtjyVtO0VR4vjvlmVm5m5VVVvesWQh88pZiNu2sp31iT7FBERHpUQhOUuze5+8lAMdERz5R41cKztTHuaMvjxXGnu5e5e1lRUVHHgaeQi08cycCsdP7wcmWyQxER6VE90ovP3fcATwGzgTwzywijioGWX6NWAiUAYfxQoDq2vNU0bZX3KYOyM5h74kgefmUbdQ26w7mI9B+J7MVXZGZ54XUOcD6wGngS+FCoNg94KLxeGIYJ45/w6MLLQuCa0MtvPDAJeBF4CZgUegVmEXWkWJio9iTT355SzP7DjTz2xo5khyIi0mMyOq5yzEYB94XedmnAg+7+sJm9ATxgZt8EVgB3h/p3A780swqiI6drANx9lZk9CLwBNAI3uXsTgJl9BlgEpBPddX1VAtuTNGdMKGT00AH8tnwzV8wYnexwRER6hPW33mFlZWVeXt77/s7qB4vf4vYn1vLMv55DScHAZIcjIv2ImS1397KeXq7uJNFLXH1aCQb85qXNHdYVEekLlKB6iTF5OZw9uYgHyzfT2NSc7HBERBJOCaoXuXbWWHbuP8wTb+p/okSk71OC6kXOPWE4w4dk84BO84lIP6AE1YtkpKdxdVkJT63ZydY9h5IdjohIQilB9TJ/d1oJDjxYrqMoEenblKB6mZKCgZw1cRi/eUmdJUSkb1OC6oWuO30c2/bW8fhqdZYQkb5LCaoXOn/KcMbk5XDv8xuSHYqISMIoQfVCGelpXH/GOJaur+bN7fuSHY6ISEIoQfVSf1dWQnZGGvc9vzHZoYiIJIQSVC+VPyiLq04ewx9XVLKntj7Z4YiIdDslqF5s3pxS6hqa1eVcRPokJahebOroXGaNL2DBCxtpau5fd6UXkb5PCaqX++icUiprDrFYf2YoIn2MElQvd8HUERTn5/CLZ9cnOxQRkW6lBNXLZaSn8fH3TWD5xhrK365OdjgiIt0mYQnKzErM7EkzW21mq8zsc6G8wMwWm9na8Jwfys3MbjezCjN71cxOiSnkWccAABX+SURBVJnXvFB/rZnNiyk/1cxeC9PcbmaWqPaksg+XFZM/MJOfPa2jKBHpOxJ5BNUIfMHdpwCzgZvMbCpwM7DE3ScBS8IwwMXApPCYD9wBUUIDbgFOB2YBt7QktVBnfsx0cxPYnpQ1MCuDG84o5fHVO6jYuT/Z4YiIdIuEJSh33+buL4fX+4HVwBjgSuC+UO0+4Krw+kpggUeWAnlmNgq4CFjs7tXuXgMsBuaGcbnu/oK7O7AgZl79zg1njGNAZhp3PqOjKBHpG3rkGpSZlQIzgWXACHffBlESA4aHamOA2B/0VIay9sor45THW/58Mys3s/KqqqquNiclFQ7O5uqyEv64Ygs79tUlOxwRkS5LeIIys8HA74HPu3t7N46Ld/3Ij6H8yEL3O929zN3LioqKOgq51/rHsybQ1Ozc85xuIisivV9CE5SZZRIlp/vd/Q+heEc4PUd4bvnPiEqgJGbyYmBrB+XFccr7rbGFA7nspNH8culGqg/q9kci0rslshefAXcDq939+zGjFgItPfHmAQ/FlN8QevPNBvaGU4CLgAvNLD90jrgQWBTG7Tez2WFZN8TMq9/6p3MncqihSb+LEpFeL5FHUGcC1wPnmtnK8LgEuA24wMzWAheEYYC/AOuBCuAXwKcB3L0a+AbwUnjcGsoAPgXcFaZZBzySwPb0CpNGDOHS6aNY8Pzb1OgoSkR6MYs6wPUfZWVlXl5enuwwEuqtHfu56IfP8OkPHMe/XnRCssMRkV7OzJa7e1lPL1d3kuiDJo8YwiXTR3HvX3UUJSK9lxJUH/XZcydxsL6Ju9WjT0R6KSWoPur4kdG1qP/+6wZ2HTic7HBERI6aElQf9s8XTqausZmfPFGR7FBERI6aElQfdlzRYK4uK+H+ZRvZtLs22eGIiBwVJag+7vPnTyI9zfje4jXJDkVE5KgoQfVxI3IH8LEzx/PQyq28vmVvssMREek0Jah+4BNnH0fewEy+s0hHUSLSeyhB9QNDczL5zDkTeeatKp5as7PjCUREUoASVD9xwxmljB82iFsffoP6xuZkhyMi0iElqH4iKyON/3vZVNZXHeS+599OdjgiIh1SgupHzjlhOOccX8TtS9ZStV8/3hWR1KYE1c/8+2VTqWts4ruL3kx2KCIi7VKC6mcmFA3mY2eO58HySlZsqkl2OCIibVKC6of+6bxJjMwdwJf/8BoNTeowISKpSQmqHxqcncGtV07jze379c+7IpKyEvmX7/eY2U4zez2mrMDMFpvZ2vCcH8rNzG43swoze9XMTomZZl6ov9bM5sWUn2pmr4Vpbg9/+y6ddOG0kVx84kh+9Pha3t51MNnhiIgcIZFHUPcCc1uV3QwscfdJwJIwDHAxMCk85gN3QJTQgFuA04FZwC0tSS3UmR8zXetlSQe+dsU0stLT+Or/vkZ/+2dlEUl9CUtQ7v4MUN2q+ErgvvD6PuCqmPIFHlkK5JnZKOAiYLG7V7t7DbAYmBvG5br7Cx5tWRfEzEs6aUTuAL508Qn8tWI3vy2vTHY4IiLv0dPXoEa4+zaA8Dw8lI8BNsfUqwxl7ZVXximPy8zmm1m5mZVXVVV1uRF9yd/PGsvp4wu49eE32Fytv+QQkdSRKp0k4l0/8mMoj8vd73T3MncvKyoqOsYQ+6a0NON7V88A4Au/fYWmZp3qE5HU0NMJakc4PUd4brlzaSVQElOvGNjaQXlxnHI5BsX5A7nl8qm8uKGau59Trz6RvmZ/XUOvPEOS0cPLWwjMA24Lzw/FlH/GzB4g6hCx1923mdki4NsxHSMuBL7s7tVmtt/MZgPLgBuAH/dkQ/qaD51azOOrd/D/Fr3F+yYVMWVUbrJDEkmY5manvqmZhqZmGpqchqZm6hvfOxyv31C8vsJmYBhpaeHZQpkZaWYYRM/vKX+3LsZ76qVFMwzzebduszuHGpo4VN9EXUMTdQ3NHGpoora+kX11jeypraf6YD01B+upqW2g+mA92/fVsX1vHQcONzIiN5tlXzk/we9s90pYgjKzXwMfAIaZWSVRb7zbgAfN7EZgE/DhUP0vwCVABVALfBQgJKJvAC+Fere6e0vHi08R9RTMAR4JDzlGZsa3/2Y6F/3wWT776xU89JkzGZjV0/svIp3n7uw91MCuA9GGufrgYXYfrKf6QD27D9ZTU1vPgbpGDhxu5GB9IwcPN0WvDzdSW9+U7PATJndABvmDssgfmMVxRYM4a+IwRucNYHReTrJDO2rW37oXl5WVeXl5ebLDSFnPrq3ihnte5G9mjuF7H56Bfl4mydLc7GzbV8eGqoNs2VPLlj11bNtziG1769i65xBb9x6iriH+nVCGZEcb6cHZGQzOzmBQdjqD3nkdPQZkppGVnkbmOw8jKyN6nZEWjmRixNtSujve8uzQ7OB49PxOWcxzG3Xxd6dpDnUhdtroiGpgVjrZmenktDyy0hmQmcbQnEzyBmaRl5NJRnr3X7kxs+XuXtbtM+6AdpHlPd43qYjPnjuJHy1Zy+zxhVx9WknHE4l0QW19I2u276di5wE27Dr4nsfhmP8uM4PhQ7IZnZfDlFG5nDdlOCOH5jBscBYFg6JH4aBs8gdlkp2RnsQWSXdRgpIjfPa8SZRvrObfH3qd6cVDdT1KuoW7s2XPIVZv28+b2/axevs+Vm/bz9u7D75zvScjzRhbOJAJwwbz/slFjB82iNLCQZQU5DAidwCZCTg6kNSlU3wSV9X+w1x6+7PkZKXzv58+k/xBWckOSXqZfXUNvLp5Lys21bBy8x5Wbt7D7oP174wfVziQKSNzmTIqlxNGDWHyiCGU5Ock5BSVdI1O8UlKKRqSzR3/cCrX3rmUm371Mvd9bJb2XqVdm6trWbahmmXrd7Ni8x7WVR1458ho4vDBnHPCcGaU5DF1VC4njBzCoGxtfqR9WkOkTaeOy+e2v53OPz/4Cl//0yq+edX0ZIckKaQlIS1dv5ul63dTWXMIgPyBmcwcm8+VM0Zz8tg8TirOY2hOZpKjld5ICUra9cFTinlrxwF+9vQ6Jg0fwrw5pckOSZJkX10Dz1fs4um3qnh27a73JKTZEwr5+PsmMHtCIZOGDyYtTb0/peuUoKRD/3rR8VTsPMDX/rSKoiHZXDJ9VLJDkh7Q3Oy8sW0fT79VxdNrqli+qYamZmdIdgZzJiohSeIpQUmH0tOMH187k3+4exmff2AleTmZzJk4LNlhSQLU1jfyzFtVPPbGDp55q4pdB6JODSeOyeVTZx/H2ccXcXJJnq5HSo9QLz7ptL21DVz98xeorKnlgflnML14aLJDkm6w68BhlqzewWOrdvBsxS7qG5vJG5jJ2ZOLOHtyEe+bVETRkOxkhylJlKxefEpQclS2763jb+94ntr6Ru7/x9lMHa3fSPVGb+86yGNvbOexVTtYvqkGdxiTl8OF00Zw4dSRnFaar+7e8g4lqB6iBNV1G3cf5Jo7l3KooYn/ufF0ThyjI6neYH3VAR5+dRt/fnUba3bsB2DqqNx3ktKUUUN0ayuJSwmqhyhBdY9Nu2u59hdLOXC4kf+58XSd7ktRG3cf5OFXt/Hwq9tYvW0fZnDauALmnjiSC6aOoKRgYLJDlF5ACaqHKEF1n83VtVxz51L2HWrg5zecypzj1HEiFWyuruXPr0VHSq9t2QvAKWPzuOyk0VwyfRQjhw5IcoTS2yhB9RAlqO61dc8hPvLfL7Jh10G+d/XJXDFjdLJD6pe27jnEn1/dxsOvbeOVzXsAmFGSx2XTR3HJSaMY0wv/akFSh251JL3S6LwcfvuJOXz8l+V89tcr2LrnEJ94/wRdy+gBO/bVRUnp1a28vClKStPHDOXmi0/g0umjdPpOej0lKOmyoQMzWfCxWXzht69w2yNvsnrbPm774EnkZOkvD7pbzcF6Hnl9O396ZStLN+zGHaaMyuVfLzqeS6ePonTYoGSHKNJtlKCkWwzITOfH18xkysghfG/xW7y14wB3Xn+q9uK7wf66Bha/sYOFr2zlubW7aGx2JhQN4nPnTeKyk0YzcfjgZIcokhC9/hqUmc0FfgSkA3e5+23t1dc1qMR7cs1OPvfrFQB862+mc7muSx21uoYmnnhzJwtXbuWJNTupb2xmTF4Ol88YzeUzRjF1VK5Oo0qPUSeJY2Bm6cBbwAVAJfAScK27v9HWNEpQPWPj7oN87oGVrNy8h7+ZOYavXzmN3AG6o3V76hubeXZtFX96ZSuL39jBwfomioZkc+n0UVw+YzSnjM1TUpKkUCeJYzMLqHD39QBm9gBwJdBmgpKeMa5wEL/75Bn85MkKfvxEBcvW7+brV57IBVNHJDu0lFLf2Mxf1+3i0de28+iq7ew91MDQnEyuOHk0l580mtMnFJKuG7FKP9XbE9QYYHPMcCVweutKZjYfmA8wduzYnolMyEhP4/PnT+bsyUXc/PvX+PiCci6YOoKvXTGtX3d7rmto4um3qnj09e08vnoH++saGZydwflThnPFyaM5a2IRWRm6zZBIb09Q8XYtjzhn6e53AndCdIov0UHJe80cm8/Dnz2Le57bwA8fX8t533uKj545nk+efVy/+SO7A4cbefLNnTz6+naeXLOT2vomhuZkctG0kVx84kjOmjSM7Az1ehSJ1dsTVCVQEjNcDGxNUizSjsz0ND5x9nFcetIovrtoDXc8tY5fLdvEJ88+juvPGMfgvvD33/ffD1/9KmzaBGPHsumW/+CJybNZ8uZOlq2vpr6pmWGDs7hq5hguPnEksycU6m8rRNrR2ztJZBB1kjgP2ELUSeLv3X1VW9Ook0RqWLV1L99dtIan1lQxZEAG150+jo/MKe29t+G5/34aPvkpyvNLefK4MpZMnMW6wmjfaULRIM49fjgXTB1BWWmBrilJr6NefMfIzC4BfkjUzfwed/9We/WVoFLLik013PXsBh55fRvpacYFU0fwoVOLef+kopT/u4emZmf1tn08v24Xf73r97w0bAK1WTlkNjUwe9NrnLOunHPrtlL6+kvJDlWkS5SgeogSVGratLuW+154mz+u2EL1wfp3ulefP2UEs8YXxO800OqUGt/6Flx33dEvvJPzOVTfxGtb9rJiUw3LN9awbEM1ew81ADBx1ybmbHyVORtf4ayNKxlcfyiayAyam48+JpEUogTVQ5SgUlt9YzNPrtnJ75dX8vRbVRxubGZIdgZnThzGrPEFlJXmM3VULhkP/Brmz4fa2ncnzsyE3Fyoroa0NGhqOnIBrdf3++8/cj4DB7L3jl+w9v1zWbNjP6u37WPFpj28uX0/Tc3R9GMLBjJ7QgFzjhvGGccVMuKkE2DjxiOXN24cvP12198YkSRSguohSlC9x6H6Jp6r2MWS1Tt4du0utuyJjkpyMtOZvH0dE7dWMGnXZkprtjLiwG5GHKim6GANmc1xElOM5qZm9tU1sKe2gd0XXUblIady6PDwGEFFYQnbcoveqT8kO4OTSoYysySfmWPzOLkkj8LBrf4CvY1Ex513HttRnUgKUYLqIUpQvde2vYcof7uGFZv28Nb9f+CtwrHsHFJ4RL0BDXUMbDhMTkMd2Y0NNKWl0ZiWTlNaOnUZWewdmHvEgRRAQe1exuzdycTdm5m8axPH3/sTJo8Ywpi8nM7dwaG7TjmKpBglqB6iBNVHlJbCxo3szR7E5ryR7BhcwI7BBewcXMDBrBwOZWRTmzmA+oxMMpsaSfcmMpqbyW6sZ+gX/5m8nEzyBmaSf9N8ite+zph9OxnYcPjd+evUnMg7dKsjkaPxrW/B/PkMrT3I0B3rOHHHus5Pu/hn777+zHXRqbnY5DRwYDR/EUmq1O7HK9KW666Lru+MGxf1lCsshKysrs9n3DhdNxJJETrFJ31H7DWgzvbiE5EOJesUn46gpO+47rroulFzMzQ2Rsmo9UNEeg0lKBERSUlKUCIikpKUoEREJCUpQYmISEpSghIRkZTU77qZm1kVEOeuniltGLAr2UH0MLW5f1Cbe4dx7l7UcbXu1e8SVG9kZuXJ+A1CMqnN/YPaLO3RKT4REUlJSlAiIpKSlKB6hzuTHUASqM39g9osbdI1KBERSUk6ghIRkZSkBCUiIilJCaoXMLN/MTM3s2Fh2MzsdjOrMLNXzeyUZMfYXczsu2b2ZmjXH80sL2bcl0Ob15jZRcmMszuZ2dzQpgozuznZ8SSCmZWY2ZNmttrMVpnZ50J5gZktNrO14Tk/2bF2NzNLN7MVZvZwGB5vZstCm39jZsfwR2b9gxJUijOzEuACYFNM8cXApPCYD9yRhNASZTFworufBLwFfBnAzKYC1wDTgLnAT80sPWlRdpPQhv8i+kynAteGtvY1jcAX3H0KMBu4KbTzZmCJu08CloThvuZzwOqY4f8EfhDaXAPcmJSoegElqNT3A+CLQGxvliuBBR5ZCuSZ2aikRNfN3P0xd28Mg0uB4vD6SuABdz/s7huACmBWMmLsZrOACndf7+71wANEbe1T3H2bu78cXu8n2mCPIWrrfaHafcBVyYkwMcysGLgUuCsMG3Au8LtQpc+1uTspQaUwM7sC2OLur7QaNQbYHDNcGcr6mo8Bj4TXfbXNfbVdbTKzUmAmsAwY4e7bIEpiwPDkRZYQPyTawWwOw4XAnpidsD7/eXdFRrID6O/M7HFgZJxRXwW+AlwYb7I4Zb3m9wLttdndHwp1vkp0Wuj+lsni1O81bW5HX21XXGY2GPg98Hl33xcdUPRNZnYZsNPdl5vZB1qK41Tts593VylBJZm7nx+v3MymA+OBV8KXuBh42cxmEe11lcRULwa2JjjUbtNWm1uY2TzgMuA8f/eHer26ze3oq+06gpllEiWn+939D6F4h5mNcvdt4TT1zuRF2O3OBK4ws0uAAUAu0RFVnpllhKOoPvt5dwed4ktR7v6auw9391J3LyXakJ3i7tuBhcANoTffbGBvy2mS3s7M5gJfAq5w99qYUQuBa8ws28zGE3UQeTEZMXazl4BJoWdXFlFHkIVJjqnbhWsvdwOr3f37MaMWAvPC63nAQz0dW6K4+5fdvTh8f68BnnD364AngQ+Fan2qzd1NR1C901+AS4g6CtQCH01uON3qJ0A2sDgcOS5190+6+yozexB4g+jU303u3pTEOLuFuzea2WeARUA6cI+7r0pyWIlwJnA98JqZrQxlXwFuAx40sxuJeqp+OEnx9aQvAQ+Y2TeBFUSJW+LQrY5ERCQl6RSfiIikJCUoERFJSUpQIiKSkpSgREQkJSlBiYhISlKCEhGRlKQEJSIiKUkJSiSJzOy08N9XA8xsUPivpBOTHZdIKtAPdUWSLNxRYACQA1S6+38kOSSRlKAEJZJk4R58LwF1wJy+cAsnke6gU3wiyVcADAaGEB1JiQg6ghJJOjNbSPRPuuOBUe7+mSSHJJISdDdzkSQysxuARnf/lZmlA8+b2bnu/kSyYxNJNh1BiYhIStI1KBERSUlKUCIikpKUoEREJCUpQYmISEpSghIRkZSkBCUiIilJCUpERFLS/wfz7JYrYihFSwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Golden section search\n",
    "\n",
    "# Set algorithm parameters and initialize arrays\n",
    "MAX_ITER = 100\n",
    "MIN_XDIF = 1E-6\n",
    "a0 = -20;\n",
    "b0 = 20;\n",
    "rho = 0.382;\n",
    "alla = np.zeros((MAX_ITER,1))\n",
    "allb = np.zeros((MAX_ITER,1))\n",
    "allfa = np.zeros((MAX_ITER,1))\n",
    "allfb = np.zeros((MAX_ITER,1))\n",
    "\n",
    "## Iterate\n",
    "k=0;\n",
    "dx = MIN_XDIF + 1;\n",
    "while k<MAX_ITER and np.abs(dx) > MIN_XDIF:\n",
    "    a1 = a0 + rho*(b0-a0); # Note my implementation does not take advantage of the golden section cleverness, for simplicity\n",
    "    b1 = b0 - rho*(b0-a0);\n",
    "    fa,fa1,fa2 = myfunction(a1);\n",
    "    fb,fb1,fb2 = myfunction(b1);\n",
    "    alla[k] = a1\n",
    "    allb[k] = b1\n",
    "    allfa[k] = fa\n",
    "    allfb[k] = fb\n",
    "    \n",
    "    if fb<fa:\n",
    "        a0 = a1;\n",
    "    else:\n",
    "        b0 = b1;\n",
    "    k = k+1;\n",
    "    dx = b0-a0;\n",
    "\n",
    "## Plot iterations and print final xhat    \n",
    "plt.plot(alla[0:k],allfa[0:k],'or')\n",
    "plt.plot(x,f)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('f')\n",
    "plt.title('Golden section search iterations (showing a1 at each iteration)')\n",
    "print('Number of iterations = ')\n",
    "print(k)\n",
    "print('Final answer xhat = ')\n",
    "print((a1+b1)/2)   \n",
    "fhat,fhat1,fhat2 = myfunction((a1+b1)/2);\n",
    "print('Value of f at xopt')\n",
    "print(fhat)\n",
    "print('Value of first derivative at xopt')\n",
    "print(fhat1)\n",
    "print('Value of second derivative at xopt')\n",
    "print(fhat2)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of iterations = \n",
      "25\n",
      "Final answer xhat = \n",
      "[-13.74609179]\n",
      "Value of f at xopt\n",
      "[-3461.19133005]\n",
      "Value of first derivative at xopt\n",
      "[4.95599898e-05]\n",
      "Value of second derivative at xopt\n",
      "[53.97745368]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEWCAYAAACnlKo3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZhU5Zn///fdO3uvbN1Ag7TKJiAtIJrgLq44k8QxcYQxJiSOZrkm850xyczXiUnm60x+mSROMmZMNMGMozGLEY1KEPcFpBEVAYFmb5amoZu1aXq7f3/U06aEZiu6urq6Pq/rqqvOec5zTt2nlbrrWc455u6IiIjEIi3RAYiISPJSEhERkZgpiYiISMyUREREJGZKIiIiEjMlERERiZmSiIiIxExJRCQJmdlPzeyfExzDCjO7KJExSOIpiUhCmdlGM6s2s15RZZ8zs5c6+HNeMrPPdeQx2/mMb5jZv7ZT/jdm9tpx4nIzG39E+R9C+UXt7efuX3T3b4e6F5lZVQecwjGZ2S/N7DtHxDDG3V+K5+dK16ckIl1BBvCVRAfRAa4GnolhvzXArLYVMysApgI1HRTXcZlZRmd8jnRPSiLSFXwP+Hszy21vo5mdbWYLzKzWzFab2Y2hfLiZ7TGztLD+czPbGbXf/5jZV83su8DHgB+b2QEz+3HYPs3MlpjZ3vA+LWrfl8zs22b2upntN7M/mVnhsU7AzPKAM4E3Yzj/R4C/MrP0sP5p4Amg8Tif90sz+05owT0LDA7ndsDMBptZmpndZWbrzGy3mT1uZvlh39LQyrnNzDYDL4Ty35jZjvD3eMXMxoTyOcDNwD+E4z8Vyjea2WVhOdvMfmhm28Lrh2aWHbZdZGZVZvY1M9tpZtvN7Naoc7nazFaGv/NWM/v7GP6GkiBKItIVVAAvAUd9eYQvyQXA/wL9iXzB/peZjXH3DcA+YGKo/jHggJmNCusfB152928CrwJ3untvd78zfKH+EbgPKAD+A/hjaAW0+Qxwa/jcrPbii3IlsNDdW0715IFtwErgirA+C3j4ZHZ094PAVcC2cG693X0b8GXgBmA6MBioA35yxO7TgVEhdogkozIi5/s2keSGuz8Qlv89HP+6dkL5JpHW0wRgPDAZ+Keo7QOBfkAxcBvwk5B4AR4EvuDufYCxhKQmyUFJRLqK/wt8ycyKjii/Ftjo7r9w92Z3fxv4HfDJsP1lYLqZDQzrvw3rw4G+wLvH+LxrgLXu/qtw3EeBD4DoL8hfuPsadz8EPE7kC/JYriG2rqw2DwOzzOwsINfdY2nRRPsC8E13r3L3w8C/AJ88ouvqX9z9YDg/3P0hd98fVX+8mfU7yc+7GbjH3Xe6ew3wLeCWqO1NYXuTuz8DHADOito22sz6untd+G8sSUJJRLoEd38feBq464hNw4Apodtqj5ntIfKF1ZY0XgYuItLqeIVIi2Z6eL3q7q3H+MjBwKYjyjYR+aXcZkfUcj3Qu70Dhe60y4HnjvFZJ+P3wCXAl4BfncZx2gwDnoj6m60CWoABUXW2tC2YWbqZ3Ru6v/YBG8OmY3bhHeHIv+emUNZmt7s3R61H/z0/QWQ8aZOZvWxm55/kZ0oXoCQiXcndwOf56Bf5FiJdUrlRr97ufnvY/jKRbqyLwvJrwAVEksjLUcc58pkH24h80UYbCmyNIe7ziLSWYh4Id/d6It1Jt3PqSaS95zlsAa464u+W4+5bj7HfZ4CZwGVEup1KQ7kd5zOiHfn3HBrKThy8+xJ3n0mkG+0PRFp9kiSURKTLcPdK4NdE+vPbPA2caWa3mFlmeJ3XNu7h7muBQ8BfA6+4+z6gmsiv2+gkUg2MiFp/Jhz3M2aWYWZ/BYwOn3eqTqYry8wsJ/rVTp1vANPdfeMpfn41UHBE19NPge+a2bDw4UVmNvM4x+gDHAZ2Az2BI6cqH/n3O9KjwD+Fzykk0j35PycK3MyyzOxmM+vn7k1ExrhiGVeSBFESka7mHuDDa0bcfT+RAeebiPyy3QH8G5Adtc/LRLpLNketG7Asqs6PiIwJ1JnZfe6+m8h4y9eIfHH+A3Ctu++KIeaTmdo7jUiy+/B15NRad9/m7u1eT3I87v4BkS/x9aH7ajCR850H/MnM9gOLgCnHOczDRLqgthIZ5F90xPYHiYxb7DGzP7Sz/3eITJB4D1hOZGD+O+3Ua88twMbQjfZFIj8IJEmYnmwoEjszGwC8Awx2/WOSFKSWiMjp6Qf8nRKIpCq1REREJGZqiYiISMxS7p45hYWFXlpamugwRESSxtKlS3e5+5EXAgMpmERKS0upqKhIdBgiIknDzI68MPdD6s4SEZGYKYmIiEjMlERERCRmSiIiIhIzJREREYmZkoiIiMRMSURERGKmJHISGptbuf+ldby6NubHRYiIdEtKIichM9342avrmffOST1jR0QkZSiJnAQzY8KQXN7ZsifRoYiIdClKIidp4pBcKmsOsK+hKdGhiIh0GUoiJ2nC0Fzc4b0texMdiohIl6EkcpLGD8nFDJZtrkt0KCIiXYaSyEnqm5PJyKLeGhcREYmiJHIKJgzJZdmWPehpkCIiEXFLImZ2lpm9E/XaZ2ZfNbN8M1tgZmvDe16ob2Z2n5lVmtl7ZnZu1LFmh/przWx2VPkkM1se9rnPzCxe5wMwcWgetQcb2VxbH8+PERFJGnFLIu6+2t0nuPsEYBJQDzwB3AUsdPcyYGFYB7gKKAuvOcD9AGaWD9wNTAEmA3e3JZ5QZ07UfjPidT4QaYkA6tISEQk6qzvrUmCdu28CZgJzQ/lc4IawPBN42CMWAblmNgi4Eljg7rXuXgcsAGaEbX3d/U2P9C89HHWsuDhzQG96ZqWzbLOSiIgIdF4SuQl4NCwPcPftAOG9fygvBrZE7VMVyo5XXtVO+VHMbI6ZVZhZRU1N7LcuyUhPY1xxP5apJSIiAnRCEjGzLOB64DcnqtpOmcdQfnSh+wPuXu7u5UVF7T5r/qRNHJrHym17aWhqOa3jiIh0B53RErkKeNvdq8N6deiKIrzvDOVVwJCo/UqAbScoL2mnPK4mDs2lqcVZuX1fvD9KRKTL64wk8mn+3JUFMA9om2E1G3gyqnxWmKU1FdgburvmA1eYWV4YUL8CmB+27TezqWFW1qyoY8XNxDC4rnERERHIiOfBzawncDnwhajie4HHzew2YDPwqVD+DHA1UElkJtetAO5ea2bfBpaEeve4e21Yvh34JdADeDa84qp/3xyKc3uEK9eHx/vjRES6tLgmEXevBwqOKNtNZLbWkXUduOMYx3kIeKid8gpgbIcEewp0R18RkQhdsR6DiUNzqao7RM3+w4kORUQkoZREYjBxaNu4iG7GKCKpTUkkBmMG9yMrPY2lSiIikuKURGKQk5nO2OK+VGxUEhGR1KYkEqPy0nyWV+miQxFJbUoiMSoflkdjSyvvb9WTDkUkdSmJxGjSsMiNhCs2qUtLRFKXkkiMCnpnM6KwFxUba09cWUSkm1ISOQ2ThuWxdFOdnnQoIilLSeQ0lJfmUVffxLqag4kORUQkIZRETkN5aT4ASzepS0tEUpOSyGkYUdiL/F5ZLNH1IiKSopREToOZce7QyLiIiEgqUhI5TeWleWzYdZBdB3QzRhFJPUoip6k8XC+i1oiIpCIlkdM0rqQfWRlpul5ERFJSXJOImeWa2W/N7AMzW2Vm55tZvpktMLO14T0v1DUzu8/MKs3sPTM7N+o4s0P9tWY2O6p8kpktD/vcFx6T26myM9I5p7ifrlwXkZQU75bIj4Dn3P1sYDywCrgLWOjuZcDCsA5wFVAWXnOA+wHMLB+4G5gCTAbubks8oc6cqP1mxPl82jWpNI/3t+7lUKNuxigiqSVuScTM+gIfBx4EcPdGd98DzATmhmpzgRvC8kzgYY9YBOSa2SDgSmCBu9e6ex2wAJgRtvV19zfDo3UfjjpWp5oyPJ+mFtdDqkQk5cSzJTICqAF+YWbLzOznZtYLGODu2wHCe/9QvxjYErV/VSg7XnlVO+VHMbM5ZlZhZhU1NTWnf2ZHKC/NJ81g0QaNi4hIaolnEskAzgXud/eJwEH+3HXVnvbGMzyG8qML3R9w93J3Ly8qKjp+1DHom5PJ6MF9Wbx+d4cfW0SkK4tnEqkCqtx9cVj/LZGkUh26ogjvO6PqD4navwTYdoLyknbKE2Lq8AKWbdmjh1SJSEqJWxJx9x3AFjM7KxRdCqwE5gFtM6xmA0+G5XnArDBLayqwN3R3zQeuMLO8MKB+BTA/bNtvZlPDrKxZUcfqdFNGFNDY3Mq7W/YkKgQRkU6XEefjfwl4xMyygPXArUQS1+NmdhuwGfhUqPsMcDVQCdSHurh7rZl9G1gS6t3j7m2DD7cDvwR6AM+GV0JMLs3HDBZvqGXKiIJEhSEi0qnimkTc/R2gvJ1Nl7ZT14E7jnGch4CH2imvAMaeZpgdol/PTM4e2JfFG3YTmW0sItL96Yr1DjRleD5LN9XR2Nya6FBERDqFkkgHmjoin4amVt6r0riIiKQGJZEONHl4ZCxksa4XEZEUoSTSgfJ7ZXHmgN4s0vUiIpIilEQ62JThBSzdVEdTi8ZFRKT7UxLpYFNHFFDf2ML7W/cmOhQRkbhTEulgk4fnAxoXEZHUoCTSwYr6ZDOyf2/eXKdxERHp/pRE4mDaGQW8taFW14uISLenJBIHF4ws5FBTC+/oPloi0s0picTB1BEFpBm8Xrkr0aGIiMSVkkgc9OuRybjifkoiItLtKYnEyQUjC3lnyx4OHG5OdCgiInGjJBInF4wspLnVeWuDZmmJSPelJBInk4blkZWRxuuVSiIi0n0picRJTmY655XmaVxERLo1JZE4mnZGIR/s2M+uA4cTHYqISFzENYmY2UYzW25m75hZRSjLN7MFZrY2vOeFcjOz+8ys0szeM7Nzo44zO9Rfa2azo8onheNXhn0tnudzqi4cWQjAG7p6XUS6qc5oiVzs7hPcve0xuXcBC929DFgY1gGuIvJc2TJgDnA/RJIOcDcwBZgM3N2WeEKdOVH7zYj/6Zy8scX96JuTwRvq0hKRbioR3VkzgblheS5wQ1T5wx6xCMg1s0HAlcACd6919zpgATAjbOvr7m+G57M/HHWsLiE9zZg6ooDXlEREpJuKdxJx4E9mttTM5oSyAe6+HSC89w/lxcCWqH2rQtnxyqvaKT+Kmc0xswozq6ipqTnNUzo1F5YVUlV3iM276zv1c0VEOkO8k8gF7n4uka6qO8zs48ep2954hsdQfnSh+wPuXu7u5UVFRSeKuUNdEMZFXlnbuclLRKQzxDWJuPu28L4TeILImEZ16IoivO8M1auAIVG7lwDbTlBe0k55lzKisBcleT14ZY2SiIh0P3FLImbWy8z6tC0DVwDvA/OAthlWs4Enw/I8YFaYpTUV2Bu6u+YDV5hZXhhQvwKYH7btN7OpYVbWrKhjdRlmxvQzi3hj3W7dGl5Eup14tkQGAK+Z2bvAW8Af3f054F7gcjNbC1we1gGeAdYDlcDPgL8FcPda4NvAkvC6J5QB3A78POyzDng2jucTs+lnFnHgcDNvb65LdCgiIh0qI14Hdvf1wPh2yncDl7ZT7sAdxzjWQ8BD7ZRXAGNPO9g4O/+MAjLSjJfX1DB1REGiwxER6TC6Yr0T9MnJZNKwPF5erXEREelelEQ6yfSzili5fR879zUkOhQRkQ6jJNJJpp8ZmVr8ylpdeCgi3YeSSCcZPagvRX2yNdVXRLoVJZFOYmZ8vKyIV9fW0NLa7jWRIiJJR0mkE00/q4i6+iaWb92b6FBERDqEkkgn+tjIQszQLC0R6TaURDpRXq8sxpfk8uLqnSeuLCKSBJREOtnFZ/Xn3ao9etqhiHQLSiKd7NJR/XGHFz9Qa0REkp+SSCcbM7gvg/rlsHCVkoiIJD8lkU5mZlxydn9eXVvD4eaWRIcjInJalEQS4LJRAzjY2MKi9bUnriwi0oUpiSTA+WcU0CMznYWrqhMdiojIaVESSYCczHQuGFnIwlU7idwBX0QkOSmJJMhlo/qzdc8hVlfvT3QoIiIxUxJJkEvO7g+gWVoiktTinkTMLN3MlpnZ02F9uJktNrO1ZvZrM8sK5dlhvTJsL406xtdD+WozuzKqfEYoqzSzu+J9Lh2pf98cxpf043mNi4hIEuuMlshXgFVR6/8G/MDdy4A64LZQfhtQ5+4jgR+EepjZaOAmYAwwA/ivkJjSgZ8AVwGjgU+Huknj0lEDeGeLrl4XkeQV1yRiZiXANcDPw7oBlwC/DVXmAjeE5ZlhnbD90lB/JvCYux929w1AJTA5vCrdfb27NwKPhbpJo+3qdc3SEpFkFe+WyA+BfwBaw3oBsMfdm8N6FVAclouBLQBh+95Q/8PyI/Y5VvlRzGyOmVWYWUVNTde5g+7oQX0Zkt+D+SuUREQkOcUtiZjZtcBOd18aXdxOVT/BtlMtP7rQ/QF3L3f38qKiouNE3bnMjBljBvLa2l3sb2hKdDgiIqcsni2RC4DrzWwjka6mS4i0THLNLCPUKQG2heUqYAhA2N4PqI0uP2KfY5UnlRljB9LY0soLuiGjiCShuCURd/+6u5e4eymRgfEX3P1m4EXgk6HabODJsDwvrBO2v+CRK/HmATeF2VvDgTLgLWAJUBZme2WFz5gXr/OJl4lD8ijqk838FTsSHYqIyCnLOHGVDvePwGNm9h1gGfBgKH8Q+JWZVRJpgdwE4O4rzOxxYCXQDNzh7i0AZnYnMB9IBx5y9xWdeiYdIC3NuHLMAH63dCsNTS3kZKYnOiQRkZNmqXbbjfLycq+oqEh0GB/xeuUubv75Yv77lklcOWZgosMREfkIM1vq7uXtbdMV613A5OH55PbMZP776tISkeSiJNIFZKancdmoATy/qprG5tYT7yAi0kUoiXQRM8YMZF9DM4vW7050KCIiJ01JpIu4sKyQXlnpPKsuLRFJIkoiXUROZjoXn92f+St20NyiLi0RSQ5KIl3IdeMHU3uwkTfWqUtLRJLDMZOImf0qvH+l88JJbdPPLKJPdgZPvZt0F96LSIo6XktkkpkNAz5rZnlmlh/96qwAU0lOZjpXjBnIcyt2cLi5JdHhiIic0PGSyE+B54CzgaVHvLrW1XrdyHXjB7G/oZmXV3eduw2LiBzLMZOIu9/n7qOI3E5khLsPj3qN6MQYU8oFIwvJ75XFU+9tT3QoIiIndMKBdXe/vTMCkYjM9DSuGjuQ51dWU9/YfOIdREQSSLOzuqDrxg/mUFMLC1fp9vAi0rUpiXRB55XmM6BvtmZpiUiXpyTSBaWnGdeMG8xLq2vYe0hPPBSRrktJpIu6bvwgGltaee59DbCLSNelJNJFTRiSy4jCXvzu7a2JDkVE5JjilkTMLMfM3jKzd81shZl9K5QPN7PFZrbWzH4dHm1LePztr82sMmwvjTrW10P5ajO7Mqp8RiirNLO74nUuiWBm/OW5xby1oZYttfWJDkdEpF3xbIkcBi5x9/HABGCGmU0F/g34gbuXAXXAbaH+bUCdu48EfhDqYWajiTwqdwwwA/gvM0s3s3TgJ8BVwGjg06Fut/EX55ZgBr9Xa0REuqi4JRGPOBBWM8PLgUuA34byucANYXlmWCdsv9TMLJQ/5u6H3X0DUAlMDq9Kd1/v7o3AY6Fut1Gc24PzRxTw+2VVpNpjjEUkOcR1TCS0GN4BdgILgHXAHndvu4quCigOy8XAFoCwfS9QEF1+xD7HKm8vjjlmVmFmFTU1yXU7kb88t4RNu+up2FSX6FBERI4S1yTi7i3uPgEoIdJyGNVetfBux9h2quXtxfGAu5e7e3lRUdGJA+9Crho7kJ5Z6fz+7apEhyIicpROmZ3l7nuAl4CpQK6ZZYRNJUDbFXVVwBCAsL0fUBtdfsQ+xyrvVnplZzBj7ECefnc7DU26s6+IdC3xnJ1VZGa5YbkHcBmwCngR+GSoNht4MizPC+uE7S94ZCBgHnBTmL01HCgD3gKWAGVhtlcWkcH3efE6n0T6xLkl7D/czJ9WVic6FBGRj8g4cZWYDQLmhllUacDj7v60ma0EHjOz7wDLgAdD/QeBX5lZJZEWyE0A7r7CzB4HVgLNwB3u3gJgZncC84F0IncbXhHH80mY80cUMLhfDr+p2ML14wcnOhwRkQ9Zqs36KS8v94qK5Hscyg8WrOG+F9byyv+5mCH5PRMdjoikEDNb6u7l7W3TFetJ4sbzhmDAr5dsOWFdEZHOoiSSJIpzezD9zCIer9hCc0trosMREQGURJLKpycPZef+w7zwgZ4zIiJdg5JIErnk7P7075PNY+rSEpEuQkkkiWSkp3Fj+RBeWr2TbXsOJTocERElkWTzV+cNwYHHK9QaEZHEUxJJMkPye3LhyEJ+vUQD7CKSeEoiSejmKcPYvreB51dpgF1EEktJJAldNqo/xbk9+OUbGxIdioikOCWRJJSRnsYt5w9j0fpaPtixL9HhiEgKUxJJUn9VPoTsjDTmvrEp0aGISApTEklSeb2yuGFCMU8sq2JPfWOiwxGRFKUkksRmTyuloalV031FJGGURJLY6MF9mTw8n4ff3ERLa2rdjVlEugYlkSR367RSquoOsUAPrBKRBFASSXKXjx5ASV4Pfvbq+kSHIiIpSEkkyWWkp/H5j41g6aY6KjbWJjocEUkx8XzG+hAze9HMVpnZCjP7SijPN7MFZrY2vOeFcjOz+8ys0szeM7Nzo441O9Rfa2azo8onmdnysM99ZmbxOp+u7FPlJeT1zOSnL6s1IiKdK54tkWbga+4+CpgK3GFmo4G7gIXuXgYsDOsAVwFl4TUHuB8iSQe4G5gCTAbubks8oc6cqP1mxPF8uqyeWRnMOr+U51dVU7lzf6LDEZEUErck4u7b3f3tsLwfWAUUAzOBuaHaXOCGsDwTeNgjFgG5ZjYIuBJY4O617l4HLABmhG193f1Njzwo/uGoY6WcWecPIyczjQdeUWtERDpPp4yJmFkpMBFYDAxw9+0QSTRA/1CtGIi+4KEqlB2vvKqd8vY+f46ZVZhZRU1NzemeTpdU0DubG8uH8MSyrVTva0h0OCKSIuKeRMysN/A74KvufrwbPbU3nuExlB9d6P6Au5e7e3lRUdGJQk5an7twBC2tzkOv6caMItI54ppEzCyTSAJ5xN1/H4qrQ1cU4b3tfuZVwJCo3UuAbScoL2mnPGUNLejJtecM5leLNlF7ULdCEZH4i+fsLAMeBFa5+39EbZoHtM2wmg08GVU+K8zSmgrsDd1d84ErzCwvDKhfAcwP2/ab2dTwWbOijpWyvnTJSA41tei6ERHpFPFsiVwA3AJcYmbvhNfVwL3A5Wa2Frg8rAM8A6wHKoGfAX8L4O61wLeBJeF1TygDuB34edhnHfBsHM8nKZQN6MM14wbx8BsbqVNrRETizCITm1JHeXm5V1RUJDqMuFpTvZ8rf/gKf3vRGfyfK89OdDgikuTMbKm7l7e3TVesd0NnDujD1eMG8cvX1RoRkfhSEummvnxJGQcbW3hQM7VEJI6URLqpswZGxkZ+8foGdh04nOhwRKSbUhLpxv7uijNpaG7lxy9UJjoUEemmlES6sTOKenNj+RAeWbyJzbvrEx2OiHRDSiLd3FcvKyM9zfj+gtWJDkVEuiElkW5uQN8cPnvBcJ58Zxvvb92b6HBEpJtREkkBX5h+Brk9M/n3+WqNiEjHUhJJAf16ZHLnxSN5ZU0NL63eeeIdREROkpJIiph1finDC3txz9MraWxuTXQ4ItJNKImkiKyMNP7vtaNZX3OQuW9sTHQ4ItJNKImkkIvP7s/FZxVx38K11OzXBYgicvqURFLMP187mobmFr43/4NEhyIi3YCSSIoZUdSbz14wnMcrqli2uS7R4YhIklMSSUFfurSMgX1z+Prvl9PUokF2EYmdkkgK6p2dwT0zx/DBjv16AqKInJZ4Ph73ITPbaWbvR5Xlm9kCM1sb3vNCuZnZfWZWaWbvmdm5UfvMDvXXmtnsqPJJZrY87HNfeESunKQrxgzkqrED+dHza9m462CiwxGRJBXPlsgvgRlHlN0FLHT3MmBhWAe4CigLrznA/RBJOsDdwBRgMnB3W+IJdeZE7XfkZ8kJ/Mv1Y8hKT+Obf1hOqj3hUkQ6RtySiLu/AtQeUTwTmBuW5wI3RJU/7BGLgFwzGwRcCSxw91p3rwMWADPCtr7u/qZHvv0ejjqWnKQBfXP4x6vO5vXK3fymoirR4YhIEursMZEB7r4dILz3D+XFwJaoelWh7HjlVe2Ut8vM5phZhZlV1NTUnPZJdCefmTyUKcPzuefplWyp1e3iReTUdJWB9fbGMzyG8na5+wPuXu7u5UVFRTGG2D2lpRnfv3E8AF/7zbu0tKpbS0ROXmcnkerQFUV4b7sbYBUwJKpeCbDtBOUl7ZRLDEryenL3daN5a0MtD76m2Voi3c3+hqa49TRkxOWoxzYPmA3cG96fjCq/08weIzKIvtfdt5vZfOBfowbTrwC+7u61ZrbfzKYCi4FZwH925ol0N5+cVMLzq6r5/+av4WNlRYwa1DfRIYnETWur09jSSlNLK00tTlNLK43NH11vb65Je3NAzcAw0tLCu4UyM9LMMIi8f6T8z3UxPlIvLXLAcJw/121151BTC4caW2hoaqGhqZVDTS3UNzazr6GZPfWN1B5spO5gI3X1TdQebGTHvgZ27G3gwOFmBvTNZvE3Luvwv2XckoiZPQpcBBSaWRWRWVb3Ao+b2W3AZuBTofozwNVAJVAP3AoQksW3gSWh3j3u3jZYfzuRGWA9gGfDS2JkZvzrX4zjyh++ypcfXcaTd15Az6zO/o0hcvLcnb2Hmth1IPLlWXvwMLsPNlJ7oJHdBxupq2/kQEMzBw43c7CxmYOHWyLLh5upb2xJdPhx0zcng7xeWeT1zOKMol5cOLKQwbk5DM7tEZfPs1Sb2lleXu4VFRWJDqPLenVtDbMeeou/mFjM9z81Hl1+I3HxyCPwzW/C5s0wdCh897tw880fqdLa6mzf18CGmoNs3VPP1j0NbN9ziO17G9i25xDb9h6ioan9Oy70yY58kfbOzqB3dga9stPp9eFy5JWTmQUV4u4AAA9fSURBVEZWehqZH76MrIzIckZaaBFEae+b0t3xtneHVgfHI+8flkW9H6Mu/ud9WkNdiN430jLpmZVOdmY6PdpeWenkZKbRr0cmuT2zyO2RSUZ6x49SmNlSdy9vb5t+aspHfKysiC9fUsaPFq5l6vACbjxvyIl3EjkVjzwCc+ZAfaSPvn7bDlb/879TWZfFhpHj2LDr4Ievw1HPvjGD/n2yGZzbg1GD+nLpqP4M7NeDwt5Z5PeKvAp6ZZPXK5PsjPREnV3KURKRo3z50jIqNtXyz0++z7iSfhofkQ7h7mzdc4hVP3mMD8Zfx6r+w1nVfzgb8wbhlgZVkLFtPUMLejKisDcfP7OI4YW9KC3oxZD8Hgzom0NmHH5ly+lRd5a0q2b/Ya6571V6ZKXzh7+9gLxeWYkOSZLMvoYm3tuyl2Wb63hnyx7e2bKH3QcbP9w+rG4bo3ZuYNTODZxds5Ezd29myM4tcemOkdOj7iw5ZUV9srn/ryfx6QcWccf/vs3cz07Wr0A5ri219SzeUMvi9btZtmUP62oOfDjDaWT/3lx8dn/GD8ll9B1/w9nLF9GrqeGjBxg2DPT/WNJREpFjmjQsj3s/MY6/e/xdvvXUCr5zw7hEhySJcIxB8LaksWj9bhat301V3SEA8npmMnFoHjPHD2bC0FzOKcmlX4/MPx/v7z4Hc96CpqjP6NkzclxJOkoiclx/eW4Ja6oP8NOX11HWvw+zp5UmOiTpTFGD4PuyevJG9mBefrSCVz/oQ1VTZPA6r2cmU0cU8PmPjWDqiALK+vcmLe04s/raZmGdYHaWJAeNicgJtbQ6X/jVUhZ+UM1PPnMuV48blOiQpBO0tjorJ32cl3MG8fKISSwtHkVLWjp9Dh9k2s61TLvj5pNLGpL0jjcmoiQiJ+VQYwt//eBillft5Ze3nse0kYWJDknioL6xmVfW1PCnldW8sqaGXQciA+Fjd1Ry0foKpq9/mwnbV5PprdCqp2KmCg2sy2nrkZXOQ7PP48b/fpPPP1zBY3POZ1xJv0SHJR1g14HDLFxVzZ9WVPNq5S4am1vJ7ZnJ9DOLmP7T/8fHliygqH7PR3caNiwxwUqXoyQiJ61fz0zmfnYyn7j/DWY9tJhHPjeV0YN1DUky2rjrIH9auYM/rahm6eY63KE4twc3TxnKFaMHcl5pXmSqbctMWPLHj+6sQXCJou4sOWWbdh/kpgcWcaiphf+5bQpji9UiSQbraw7w9Hvb+eN721ldvR+A0YP6csWYAVwxeiCjBvVp/zY3J3GLEuneNCYSRUmkY2zeXc+nf7aIA4eb+Z/bpqhrq4vatPsgT7+3naff286q7fswg/OG5TNj7EAuHz2AIfk9Ex2iJAElkShKIh1nS209Nz2wiH2HmvjvWZOYdoYG27uCLbX1/HF5pMWxfOteAM4dmsu15wzm6nGDGNgvJ8ERSrJREomiJNKxtu05xN/84i027DrI92+cwPXjByc6pO7tGF1L2/Yc4o/vbefp5dt5d0tkEHz8kFyuHTeIq88ZRHGcbgMuqUGzsyRuBuf24DdfmMbnf1XBlx9dxrY9h/jCx0foFvLxcMTdb6t37+ePP3mCp9f15u36yD/lccX9uOuqs7lm3CB1VUmnUEtEOkRDUwtf+827/PG97cycMJh7//IcemTpdtwdqrSUuupanj1rGk+N+jiLho7DLY1RdVVce9OlXDNuEKWFvRIdpXRDaolI3OVkpvOfN01k1MA+fH/BGtZUH+CBWybp13AH2N/QxIKV1cw771ZeK51Ac3oGI3Zv4SuvP8q1q15lZN1W+G9d+CeJkfQtETObAfwISAd+7u73Hq++WiLx9+LqnXzl0WUAfPcvxnGdxklOWUNTCy98sJN572zjhdU7aWxupfjgbq5b/gLXrXqF0Ts38GGH4bBhsHFjAqOV7q7btkTMLB34CXA5UAUsMbN57r4ysZGltovP6s9TX7qQrzz2Dl96dBkvfLCTb80cQ9+czBPvnMIam1t5dW0NT727jQUrqznY2EJRn2w+M3ko140fzLmvPoM99JsPx0QAXfgnCZfUSQSYDFS6+3oAM3sMmAkoiSTYsIJe/PaL5/PjFyv5zxcqWbx+N9+aOZbLRw9IdGhdSmNzK6+v28Vzy3fw3Iod7D3URL8emVw/YTDXnTOYKSMKSG+7ueGwm8HQhX/SpSR1d5aZfRKY4e6fC+u3AFPc/c4j6s0B5gAMHTp00qZNmzo91lS2bHMdd/1uOaur93P56AH8y/VjUnrKaUNTCy+vqeG593fw/Kpq9jc00zs7g8tG9ef6CYO5cGQRWRl6OJN0Hd22Owtobx7pUVnR3R8AHoDImEi8g5KPmjg0j6e/fCEPvbaBHz6/lku//xK3XjCcL04/46MPK0p2x7k9yIHDzbz4wU6ee38HL67eSX1jC/16ZHLlmIFcNXYgF5YVkp2h2WySfJI9iVQBQ6LWS4BtCYpFjiMzPY0vTD+Da84ZxPfmr+b+l9bxv4s388XpZ3DL+cPonZ3k/ysecQ0Hmzax+e//iRd2ZbGw91AWr6+lsaWVwt5Z3DCxmKvGDmTqiAI9cliSXrJ3Z2UAa4BLga3AEuAz7r7iWPtodlbXsGLbXr43fzUvra6hT04GN08Zxt9MK03eW3KUltK0pYqK4tG8eEY5C0dOZl1B5PfNiKJeXHJWfy4fPYDy0vw/j3GIJIlufdsTM7sa+CGRKb4Puftxp6ooiXQtyzbX8fNXN/Ds+9tJTzMuHz2AT04q4eNlRZFbkXdhLa3Oqu37eGPdLl7/8SMsKRlNfVYPMluamLp5ORevq+CS9RWU1m5NdKgip6VbJ5FTpSTSNW3eXc/cNzfyxLKt1B5spKhPNteMG8RlowYweXh+5w40H2Ns41BjC8u37mXZ5jqWbqpj8YZa9h5qAmDk3u1Mq1zKtE3vcuGmd+jdeChyLF3DId2AkkgUJZGurbG5lRdX7+R3S6t4eU0Nh5tb6ZOdwQUjC5k8PJ/y0jxGD+p7cq2UY92/q6AAbrwRnnnm6EHwMLaxt8VYWziU1YXDWDW4jGWTL+WDxkxaWiP/Xobm92TqiHymnVHI+WcUMOCp3310TAQi13A88ICm4ErSUxKJoiSSPA41tvBa5S4Wrqrm1bW72Lon8uu+R2Y6Zw7ozcj+fSgb0JvSgl4M6JvNgL45FPXJjgxWH+cGkK0Y+3J6sSenD7t79aOq7wCqCoupuvYTVK1YR2WvIrb3Lfqwfp/DBzmnbgsTP30tE4fmMmFILgW9s48+sB7eJN2UkkgUJZHktX3vISo21rFs8x7WVO9nTfV+du4/fFS9nMw0eu6to0dTA9nNTbSkpdGclk5LWjoNGVnszemN29EtmfyG/RTv2cHIXVs4c9cmzqrZxJm7NlG8ryZyV+JW3Z9KUlN3vk5EUsigfj24bnyPj9yLa299E1vq6qne10D1vsPs3N/AwcPNHPrRU9Rn5tCYkUlmSzPp3kJGayvZzY30a9hP7qED5DbsJ+/QPkr27qR43056NjdGWhDtXYw6dGgnnqlI8lASkaTWr2cm/Xr2O/o579fef+oHGzYs0gXV3tiG7k8l0q6uPYdSpLO0JYqbb44Mhg8bFhlXGTZMg+Mix6EkIt3T8cb6Cgrg9tuPnShuvjkyLbe1NfKuBCJyTOrOku4rxSaNiCSCWiIiIhIzJREREYmZkoiIiMRMSURERGKmJCIiIjFLuduemFkNkGzPxy0EdiU6iE6mc04NOufkMMzdi9rbkHJJJBmZWcWx7lvTXemcU4POOfmpO0tERGKmJCIiIjFTEkkODyQ6gATQOacGnXOS05iIiIjETC0RERGJmZKIiIjETEkkCZjZ35uZm1lhWDczu8/MKs3sPTM7N9ExdhQz+56ZfRDO6wkzy43a9vVwzqvN7MpExtmRzGxGOKdKM7sr0fHEg5kNMbMXzWyVma0ws6+E8nwzW2Bma8N7XqJj7Whmlm5my8zs6bA+3MwWh3P+tZllJTrG06Ek0sWZ2RDgcmBzVPFVQFl4zQFieIxfl7UAGOvu5wBrgK8DmNlo4CZgDDAD+C8zS09YlB0knMNPiPw3HQ18Opxrd9MMfM3dRwFTgTvCed4FLHT3MmBhWO9uvgKsilr/N+AH4ZzrgNsSElUHURLp+n4A/AMQPQNiJvCwRywCcs1sUEKi62Du/id3bw6ri4CSsDwTeMzdD7v7BqASmJyIGDvYZKDS3de7eyPwGJFz7Vbcfbu7vx2W9xP5Ui0mcq5zQ7W5wA2JiTA+zKwEuAb4eVg34BLgt6FK0p+zkkgXZmbXA1vd/d0jNhUDW6LWq0JZd/NZ4Nmw3F3Pubue1zGZWSkwEVgMDHD37RBJNED/xEUWFz8k8iOwNawXAHuifigl/X9vPdkwwczseWBgO5u+CXwDuKK93dopS5q52sc7Z3d/MtT5JpEukEfadmunftKc83F01/Nql5n1Bn4HfNXd90V+mHdPZnYtsNPdl5rZRW3F7VRN6v/eSiIJ5u6XtVduZuOA4cC74R9aCfC2mU0m8utlSFT1EmBbnEPtMMc65zZmNhu4FrjU/3whU1Kf83F01/M6ipllEkkgj7j770NxtZkNcvftoUt2Z+Ii7HAXANeb2dVADtCXSMsk18wyQmsk6f97qzuri3L35e7e391L3b2UyJfNue6+A5gHzAqztKYCe9u6BJKdmc0A/hG43t3rozbNA24ys2wzG05kUsFbiYixgy0BysKMnSwikwfmJTimDhfGAh4EVrn7f0RtmgfMDsuzgSc7O7Z4cfevu3tJ+Pd7E/CCu98MvAh8MlRL+nNWSyQ5PQNcTWRwuR64NbHhdKgfA9nAgtACW+TuX3T3FWb2OLCSSDfXHe7eksA4O4S7N5vZncB8IB14yN1XJDiseLgAuAVYbmbvhLJvAPcCj5vZbURmIH4qQfF1pn8EHjOz7wDLiCTXpKXbnoiISMzUnSUiIjFTEhERkZgpiYiISMyUREREJGZKIiIiEjMlERERiZmSiIiIxExJRCSBzOy88OyUHDPrFZ61MTbRcYmcLF1sKJJg4crlHKAHUOXu/y/BIYmcNCURkQQL98xaAjQA07rD7Vwkdag7SyTx8oHeQB8iLRKRpKGWiEiCmdk8Ik80HA4Mcvc7ExySyEnTXXxFEsjMZgHN7v6/4Xnrb5jZJe7+QqJjEzkZaomIiEjMNCYiIiIxUxIREZGYKYmIiEjMlERERCRmSiIiIhIzJREREYmZkoiIiMTs/wf4HcwQQ8Y/EwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Newton's method\n",
    "\n",
    "## Set algorithm parameters\n",
    "MAX_ITER = 100 # Max number of iterations\n",
    "MIN_XDIF = 1E-6 # Min change in x allowable (otherwise stop)\n",
    "x0 = 0 # Initial guess\n",
    "mu0 = 0 # mu0 > 0 for LM instead of Newton's\n",
    "k=0 # Counter\n",
    "xk = np.zeros((MAX_ITER,1)) # Here we will store our x values as we iterate\n",
    "fk = np.zeros((MAX_ITER,1)) # Here we will store our f values\n",
    "xk[0] = x0 # Initial guess\n",
    "dx = MIN_XDIF + 1 # Just so it doesn't stop before starting the iterations\n",
    "\n",
    "## Iteratively try to find the optimum\n",
    "while k<(MAX_ITER-1) and np.abs(dx) > MIN_XDIF: # Keep iterating unless we reach NUM_ITER or x doesn't change much\n",
    "    f0,f1,f2 = myfunction(xk[k]) # Calculate f, f', f'' for quadratic approximation\n",
    "    fk[k] = f0\n",
    "    k = k+1\n",
    "    xk[k] = xk[k-1] - f1/(f2+mu0) # Update x based on current quadratic approximation\n",
    "    dx = xk[k] - xk[k-1]\n",
    "\n",
    "fhat,fhat1,fhat2 = myfunction(xk[k]) # Calculate the final f value, for plot\n",
    "fk[k] = fhat\n",
    "\n",
    "## Plot iterations and print final xhat\n",
    "plt.plot(xk[0:k+1],fk[0:k+1],'or')\n",
    "plt.plot(x,f)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('f')\n",
    "plt.title('Newton / LM iterations')\n",
    "print('Number of iterations = ')\n",
    "print(k)\n",
    "print('Final answer xhat = ')\n",
    "print(xk[k])\n",
    "print('Value of f at xopt')\n",
    "print(fhat)\n",
    "print('Value of first derivative at xopt')\n",
    "print(fhat1)\n",
    "print('Value of second derivative at xopt')\n",
    "print(fhat2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}