mp573_lecture_optim_007.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimization Lecture 7: 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": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## MP 573 Part 2, Lecture 7: 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": 2,
   "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": 3,
   "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": 4,
   "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": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of iterations = \n",
      "8\n",
      "Final answer xhat = \n",
      "[-13.74609271]\n",
      "Value of f at xopt\n",
      "[-3461.19133005]\n",
      "Value of first derivative at xopt\n",
      "[5.68434189e-14]\n",
      "Value of second derivative at xopt\n",
      "[53.97745751]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEWCAYAAACnlKo3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxV9Z3/8dcnewKEhBDCTkBAWRSRFHBpaV1xxW6Olp8wrR3ajnZ5jDMdW3+/n62tU2f667R12rFD1RY7Tq1dpbaKiHsVNAiKgEBkDWsgYQ0h2+f3x/2GXiBsITcnN/f9fDzu457zPd9z7udEuZ/7Xc455u6IiIi0RVrUAYiISPJSEhERkTZTEhERkTZTEhERkTZTEhERkTZTEhERkTZTEhERkTZTEhFJQmb2EzP7PxHHsNzMPhxlDBI9JRGJlJmtN7PtZtYtruyzZvZiO3/Oi2b22fY8Ziuf8XUz+5dWyv/WzF49QVxuZuOOKv9DKP9wa/u5++fd/Vuh7ofNrLIdTuG4zOznZvbto2IY4+4vJvJzpfNTEpHOIAP4ctRBtINrgD+3Yb/VwIyWFTMrAiYDVe0U1wmZWUZHfI50TUoi0hl8F/hHMytobaOZnWNm882s2sxWmdlNoXyome02s7Sw/pCZ7Yjb77/N7Ctmdh/wQeBHZrbfzH4Utl9kZm+a2Z7wflHcvi+a2bfM7C9mts/MnjWz3sc7ATMrBEYCr7fh/B8D/sbM0sP6LcDvgfoTfN7PzezboQX3NNA/nNt+M+tvZmlmdpeZvW9mu8zsCTPrFfYtDa2c28xsI/B8KP+1mW0Lf4+XzWxMKJ8FTAe+Go7/x1C+3swuD8vZZvYDM9sSXj8ws+yw7cNmVmlmd5rZDjPbamafjjuXa8xsRfg7bzazf2zD31AioiQinUE58CJwzJdH+JKcD/wP0IfYF+x/mtkYd18H7AXGh+ofBPab2aiw/iHgJXe/G3gFuMPdu7v7HeEL9U/AA0AR8O/An0IroMWngE+Hz81qLb44VwEL3L3pdE8e2AKsAK4M6zOAR09lR3c/AFwNbAnn1t3dtwBfAm4EpgD9gRrgx0ftPgUYFWKHWDIaQex83yKW3HD32WH538Lxr28llLuJtZ7OB8YBE4H/Hbe9L9ATGADcBvw4JF6Ah4HPuXsPYCwhqUlyUBKRzuL/Al80s+Kjyq8D1rv7z9y90d3fAn4LfCJsfwmYYmZ9w/pvwvpQIB94+zifdy2wxt1/EY77S+A9IP4L8mfuvtrdDwJPEPuCPJ5raVtXVotHgRlmdjZQ4O5tadHE+xxwt7tXuvsh4BvAJ47quvqGux8I54e7P+Lu++LqjzOznqf4edOBe919h7tXAd8Ebo3b3hC2N7j7n4H9wNlx20abWb6714T/xpIklESkU3D3d4GngLuO2jQEmBS6rXab2W5iX1gtSeMl4MPEWh0vE2vRTAmvV9y9+Tgf2R/YcFTZBmK/lFtsi1uuBbq3dqDQnXYF8MxxPutU/A64FPgi8IszOE6LIcDv4/5mK4EmoCSuzqaWBTNLN7P7Q/fXXmB92HTcLryjHP333BDKWuxy98a49fi/58eJjSdtMLOXzOzCU/xM6QSURKQzuQf4O478It9ErEuqIO7V3d2/ELa/RKwb68Nh+VXgYmJJ5KW44xz9zIMtxL5o4w0GNrch7g8Qay21eSDc3WuJdSd9gdNPIq09z2ETcPVRf7ccd998nP0+BUwDLifW7VQayu0EnxHv6L/n4FB28uDd33T3acS60f5ArNUnSUJJRDoNd68AfkWsP7/FU8BIM7vVzDLD6wMt4x7uvgY4CPwv4GV33wtsJ/brNj6JbAeGxa3/ORz3U2aWYWZ/A4wOn3e6TqUry8wsJ/7VSp2vA1Pcff1pfv52oOiorqefAPeZ2ZDw4cVmNu0Ex+gBHAJ2AXnA0VOVj/77He2XwP8On9ObWPfkf58scDPLMrPpZtbT3RuIjXG1ZVxJIqIkIp3NvcDha0bcfR+xAeebif2y3Qb8K5Adt89LxLpLNsatG7Akrs4PiY0J1JjZA+6+i9h4y53Evji/Clzn7jvbEPOpTO29iFiyO/w6emqtu29x91avJzkRd3+P2Jf42tB91Z/Y+c4FnjWzfcBCYNIJDvMosS6ozcQG+Rcetf1hYuMWu83sD63s/21iEyTeAZYRG5j/div1WnMrsD50o32e2A8CSRKmJxuKtJ2ZlQBLgf6uf0ySgtQSETkzPYF/UAKRVKWWiIiItJlaIiIi0mYpd8+c3r17e2lpadRhiIgkjcWLF+9096MvBAZSMImUlpZSXl4edRgiIknDzI6+MPcwdWeJiEibKYmIiEibKYmIiEibKYmIiEibKYmIiEibKYmcisceg9JSSEuLvT/2WNQRiYh0Cik3xfe0PfYYzJoFtbWx9Q0bYusA06dHF5eISCeglsjJ3H039XX1PDjp47xSGh5sV1sLd98dbVwiIp2AksjJbNxIZnMjP534MeaOmnJEuYhIqlMSOZnBgzHg/C2rWNp/5BHlIiKpTknkZO67D/LyGL9lFRVFg9iblQd5ebFyEZEUpyRyMtOnw+zZnN+0G7c03jn/Epg9W4PqIiIoiZya6dMZ99IfMYMl//IjJRARkUBJ5BTl52QyvLg7SzftjjoUEZFOQ0nkNJw/qIAlm3ajp0GKiMQkLImY2dlmtjTutdfMvmJmvcxsvpmtCe+Fob6Z2QNmVmFm75jZBXHHmhnqrzGzmXHlE8xsWdjnATOzRJ0PwPjBhVQfqGdjdW0iP0ZEJGkkLIm4+yp3P9/dzwcmALXA74G7gAXuPgJYENYBrgZGhNcs4EEAM+sF3ANMAiYC97QknlBnVtx+UxN1PhBriQDq0hIRCTqqO+sy4H133wBMA+aE8jnAjWF5GvCoxywECsysH3AVMN/dq929BpgPTA3b8t39dY/1Lz0ad6yEGFnSnbysdJZsVBIREYGOSyI3A78MyyXuvhUgvPcJ5QOATXH7VIayE5VXtlJ+DDObZWblZlZeVVXV5pPISE/j3AE9WaKWiIgI0AFJxMyygBuAX5+saitl3obyYwvdZ7t7mbuXFRe3+qz5UzZ+cCErtuyhrqHpjI4jItIVdERL5GrgLXffHta3h64owvuOUF4JDIrbbyCw5STlA1spT6jxgwtoaHJWbN2b6I8SEen0OiKJ3MJfu7IA5gItM6xmAk/Glc8Is7QmA3tCd9c84EozKwwD6lcC88K2fWY2OczKmhF3rIQZHwbXNS4iIpLg54mYWR5wBfC5uOL7gSfM7DZgI/DJUP5n4BqggthMrk8DuHu1mX0LeDPUu9fdq8PyF4CfA7nA0+GVUH3ycxhQkMuSjTXA0ER/nIhIp5bQJOLutUDRUWW7iM3WOrquA7cf5ziPAI+0Ul4OjG2XYE/D+YMKNM1XRARdsd4m4wcXUFlzkKp9h6IORUQkUkoibTB+cMu4SE3EkYiIREtJpA3G9O9JVnoai5VERCTFKYm0QU5mOmMH5FO+XklERFKbkkgblZX2YlmlLjoUkdSmJNJGZUMKqW9q5t3Ne6IORUQkMkoibTRhSOxGwuUb1KUlIqlLSaSNirpnM6x3N8rXV5+8sohIF6UkcgYmDClk8YYaPelQRFKWksgZKCstpKa2gferDkQdiohIJJREzkBZaS8AFm9Ql5aIpCYlkTMwrHc3enXL4k1dLyIiKUpJ5AyYGRcMjo2LiIikIiWRM1RWWsi6nQfYuV83YxSR1KMkcobKwvUiao2ISCpSEjlD5w7sSVZGmq4XEZGUlNAkYmYFZvYbM3vPzFaa2YVm1svM5pvZmvBeGOqamT1gZhVm9o6ZXRB3nJmh/hozmxlXPsHMloV9HgiPye1Q2RnpnDegp65cF5GUlOiWyA+BZ9z9HGAcsBK4C1jg7iOABWEd4GpgRHjNAh4EMLNewD3AJGAicE9L4gl1ZsXtNzXB59OqCaWFvLt5DwfrdTNGEUktCUsiZpYPfAh4GMDd6919NzANmBOqzQFuDMvTgEc9ZiFQYGb9gKuA+e5e7e41wHxgatiW7+6vh0frPhp3rA41aWgvGppcD6kSkZSTyJbIMKAK+JmZLTGzh8ysG1Di7lsBwnufUH8AsClu/8pQdqLyylbKj2Fms8ys3MzKq6qqzvzMjlJW2os0g4XrNC4iIqklkUkkA7gAeNDdxwMH+GvXVWtaG8/wNpQfW+g+293L3L2suLj4xFG3QX5OJqP757No7a52P7aISGeWyCRSCVS6+6Kw/htiSWV76IoivO+Iqz8obv+BwJaTlA9spTwSk4cWsWTTbj2kSkRSSsKSiLtvAzaZ2dmh6DJgBTAXaJlhNRN4MizPBWaEWVqTgT2hu2secKWZFYYB9SuBeWHbPjObHGZlzYg7VoebNKyI+sZm3t60O6oQREQ6XEaCj/9F4DEzywLWAp8mlrieMLPbgI3AJ0PdPwPXABVAbaiLu1eb2beAN0O9e929ZfDhC8DPgVzg6fCKxMTSXpjBonXVTBpWFFUYIiIdKqFJxN2XAmWtbLqslboO3H6c4zwCPNJKeTkw9gzDbBc98zI5p28+i9btIjbbWESk69MV6+1o0tBeLN5QQ31jc9ShiIh0CCWRdjR5WC/qGpp5p1LjIiKSGpRE2tHEobGxkEW6XkREUoSSSDvq1S2LkSXdWajrRUQkRSiJtLNJQ4tYvKGGhiaNi4hI16ck0s4mDyuitr6JdzfviToUEZGEUxJpZxOH9gI0LiIiqUFJpJ0V98hmeJ/uvP6+xkVEpOtTEkmAi84q4o111bpeRES6PCWRBLh4eG8ONjSxVPfREpEuTkkkASYPKyLN4C8VO6MORUQkoZREEqBnbibnDuipJCIiXZ6SSIJcPLw3SzftZv+hxqhDERFJGCWRBLl4eG8am5031mmWloh0XUoiCTJhSCFZGWn8pUJJRES6LiWRBMnJTOcDpYUaFxGRLk1JJIEuOqs3723bx879h6IORUQkIRKaRMxsvZktM7OlZlYeynqZ2XwzWxPeC0O5mdkDZlZhZu+Y2QVxx5kZ6q8xs5lx5RPC8SvCvpbI8zldlwzvDcBrunpdRLqojmiJfMTdz3f3lsfk3gUscPcRwIKwDnA1sefKjgBmAQ9CLOkA9wCTgInAPS2JJ9SZFbff1MSfzqkbO6An+TkZvKYuLRHpoqLozpoGzAnLc4Ab48of9ZiFQIGZ9QOuAua7e7W71wDzgalhW767vx6ez/5o3LE6hfQ0Y/KwIl5VEhGRLirRScSBZ81ssZnNCmUl7r4VILz3CeUDgE1x+1aGshOVV7ZSfgwzm2Vm5WZWXlVVdYandHouGdGbypqDbNxV26GfKyLSERKdRC529wuIdVXdbmYfOkHd1sYzvA3lxxa6z3b3MncvKy4uPlnM7eriMC7y8pqOTV4iIh0hoUnE3beE9x3A74mNaWwPXVGE9x2heiUwKG73gcCWk5QPbKW8UxnWuxsDC3N5ebWSiIh0PQlLImbWzcx6tCwDVwLvAnOBlhlWM4Enw/JcYEaYpTUZ2BO6u+YBV5pZYRhQvxKYF7btM7PJYVbWjLhjdRpmxpSRxbz2/i7dGl5EupxEtkRKgFfN7G3gDeBP7v4McD9whZmtAa4I6wB/BtYCFcBPgb8HcPdq4FvAm+F1bygD+ALwUNjnfeDpBJ5Pm00ZWcz+Q428tbEm6lBERNpVRqIO7O5rgXGtlO8CLmul3IHbj3OsR4BHWikvB8aecbAJduFZRWSkGS+trmLysKKowxERaTe6Yr0D9MjJZMKQQl5apXEREelalEQ6yJSzi1mxdS879tZFHYqISLtREukgU0bGpha/vEYXHopI16Ek0kFG98unuEe2pvqKSJeiJNJBzIwPjSjmlTVVNDW3ek2kiEjSURLpQFPOLqamtoFlm/dEHYqISLtQEulAHxzeGzM0S0tEugwlkQ5U2C2LcQMLeGHVjpNXFhFJAkoiHewjZ/fh7crdetqhiHQJSiId7LJRfXCHF95Ta0REkp+SSAcb0z+ffj1zWLBSSUREkp+SSAczMy49pw+vrKniUGNT1OGIiJwRJZEIXD6qhAP1TSxcW33yyiIinZiSSAQuPKuI3Mx0FqzcHnUoIiJnREkkAjmZ6Vw8vDcLVu4gdgd8EZHkpCQSkctH9WHz7oOs2r4v6lBERNpMSSQil57TB0CztEQkqSU8iZhZupktMbOnwvpQM1tkZmvM7FdmlhXKs8N6RdheGneMr4XyVWZ2VVz51FBWYWZ3Jfpc2lOf/BzGDezJcxoXEZEk1hEtkS8DK+PW/xX4vruPAGqA20L5bUCNuw8Hvh/qYWajgZuBMcBU4D9DYkoHfgxcDYwGbgl1k8Zlo0pYuklXr4tI8kpoEjGzgcC1wENh3YBLgd+EKnOAG8PytLBO2H5ZqD8NeNzdD7n7OqACmBheFe6+1t3rgcdD3aTRcvW6ZmmJSLJKdEvkB8BXgeawXgTsdvfGsF4JDAjLA4BNAGH7nlD/cPlR+xyv/BhmNsvMys2svKqq89xBd3S/fAb1ymXeciUREUlOCUsiZnYdsMPdF8cXt1LVT7LtdMuPLXSf7e5l7l5WXFx8gqg7lpkxdUxfXl2zk311DVGHIyJy2hLZErkYuMHM1hPrarqUWMukwMwyQp2BwJawXAkMAgjbewLV8eVH7XO88qQydWxf6puaeV43ZBSRJJSwJOLuX3P3ge5eSmxg/Hl3nw68AHwiVJsJPBmW54Z1wvbnPXYl3lzg5jB7aygwAngDeBMYEWZ7ZYXPmJuo80mU8YMKKe6Rzbzl26IORUTktGWcvEq7+2fgcTP7NrAEeDiUPwz8wswqiLVAbgZw9+Vm9gSwAmgEbnf3JgAzuwOYB6QDj7j78g49k3aQlmZcNaaE3y7eTF1DEzmZ6VGHJCJyyizVbrtRVlbm5eXlUYdxhL9U7GT6Q4v4r1sncNWYvlGHIyJyBDNb7O5lrW3TFeudwMShvSjIy2Teu+rSEpHkoiTSCWSmp3H5qBKeW7md+sbmk+8gItJJKIl0ElPH9GVvXSML1+6KOhQRkVOmJNJJXDKiN92y0nlaXVoikkSURDqJnMx0PnJOH+Yt30Zjk7q0RCQ5KIl0IteP60/1gXpee19dWiKSHI6bRMzsF+H9yx0XTmqbMrKYHtkZ/PHtpLvwXkRS1IlaIhPMbAjwGTMrNLNe8a+OCjCV5GSmc+WYvjyzfBuHGpuiDkdE5KROlER+AjwDnAMsPurVua7W60KuH9ePfXWNvLSq89xtWETkeI6bRNz9AXcfRex2IsPcfWjca1gHxphSLh7em17dsvjjO1ujDkVE5KROOrDu7l/oiEAkJjM9javH9uW5FduprW88+Q4iIhHS7KxO6Ppx/TnY0MSClbo9vIh0bkoindAHSntRkp+tWVoi0ukpiXRC6WnGtef258VVVew5qCceikjnpSTSSV0/rh/1Tc08864G2EWk81IS6aTOH1TAsN7d+O1bm6MORUTkuBKWRMwsx8zeMLO3zWy5mX0zlA81s0VmtsbMfhUebUt4/O2vzKwibC+NO9bXQvkqM7sqrnxqKKsws7sSdS5RMDM+dsEA3lhXzabq2qjDERFpVSJbIoeAS919HHA+MNXMJgP/Cnzf3UcANcBtof5tQI27Dwe+H+phZqOJPSp3DDAV+E8zSzezdODHwNXAaOCWULfL+OgFAzGD36k1IiKdVMKSiMfsD6uZ4eXApcBvQvkc4MawPC2sE7ZfZmYWyh9390Puvg6oACaGV4W7r3X3euDxULfLGFCQy4XDivjdkkpS7THGIpIcEjomEloMS4EdwHzgfWC3u7dcRVcJDAjLA4BNAGH7HqAovvyofY5X3locs8ys3MzKq6qS63YiH7tgIBt21VK+oSbqUEREjpHQJOLuTe5+PjCQWMthVGvVwrsdZ9vplrcWx2x3L3P3suLi4pMH3olcPbYveVnp/O6tyqhDERE5RofMznL33cCLwGSgwMwywqaBQMsVdZXAIICwvSdQHV9+1D7HK+9SumVnMHVsX556eyt1Dbqzr4h0LomcnVVsZgVhORe4HFgJvAB8IlSbCTwZlueGdcL25z02EDAXuDnM3hoKjADeAN4ERoTZXlnEBt/nJup8ovTxCway71Ajz67YHnUoIiJHyDh5lTbrB8wJs6jSgCfc/SkzWwE8bmbfBpYAD4f6DwO/MLMKYi2QmwHcfbmZPQGsABqB2929CcDM7gDmAenE7ja8PIHnE5kLhxXRv2cOvy7fxA3j+kcdjojIYZZqs37Kysq8vDz5Hofy/fmreeD5Nbz8Tx9hUK+8qMMRkRRiZovdvay1bbpiPUnc9IFBGPCrNzedtK6ISEdREkkSAwpymTKymCfKN9HY1Bx1OCIigJJIUrll4mB27DvE8+/pOSMi0jkoiSSRS8/pQ58e2TyuLi0R6SSURJJIRnoaN5UN4sVVO9iy+2DU4YiIKIkkm7/5wCAceKJcrRERiZ6SSJIZ1CuPS4b35ldvaoBdRKKnJJKEpk8awtY9dTy3UgPsIhItJZEkdPmoPgwoyOXnr62LOhQRSXFKIkkoIz2NWy8cwsK11by3bW/U4YhIClMSSVJ/UzaI7Iw05ry2IepQRCSFKYkkqcJuWdx4/gB+v6SS3bX1UYcjIilKSSSJzbyolLqGZk33FZHIKIkksdH985k4tBePvr6BpubUuhuziHQOSiJJ7tMXlVJZc5D5emCViERASSTJXTG6hIGFufz0lbVRhyIiKUhJJMllpKfxdx8cxuINNZSvr446HBFJMYl8xvogM3vBzFaa2XIz+3Io72Vm881sTXgvDOVmZg+YWYWZvWNmF8Qda2aov8bMZsaVTzCzZWGfB8zMEnU+ndknywZSmJfJT15Sa0REOlYiWyKNwJ3uPgqYDNxuZqOBu4AF7j4CWBDWAa4GRoTXLOBBiCUd4B5gEjARuKcl8YQ6s+L2m5rA8+m08rIymHFhKc+t3E7Fjn1RhyMiKSRhScTdt7r7W2F5H7ASGABMA+aEanOAG8PyNOBRj1kIFJhZP+AqYL67V7t7DTAfmBq25bv76x57UPyjccdKOTMuHEJOZhqzX1ZrREQ6ToeMiZhZKTAeWASUuPtWiCUaoE+oNgCIv+ChMpSdqLyylfLWPn+WmZWbWXlVVdWZnk6nVNQ9m5vKBvH7JZvZvrcu6nBEJEUkPImYWXfgt8BX3P1EN3pqbTzD21B+bKH7bHcvc/ey4uLik4WctD57yTCamp1HXtWNGUWkYyQ0iZhZJrEE8pi7/y4Ubw9dUYT3lvuZVwKD4nYfCGw5SfnAVspT1uCiPK47rz+/WLiB6gO6FYqIJF4iZ2cZ8DCw0t3/PW7TXKBlhtVM4Mm48hlhltZkYE/o7poHXGlmhWFA/UpgXti2z8wmh8+aEXeslPXFS4dzsKFJ142ISIdIZEvkYuBW4FIzWxpe1wD3A1eY2RrgirAO8GdgLVAB/BT4ewB3rwa+BbwZXveGMoAvAA+Ffd4Hnk7g+SSFESU9uPbcfjz62npq1BoRkQSz2MSm1FFWVubl5eVRh5FQq7fv46ofvMzff/gs/umqc6IOR0SSnJktdvey1rbpivUuaGRJD645tx8//4taIyKSWEoiXdSXLh3BgfomHtZMLRFJICWRLursvrGxkZ/9ZR079x+KOhwR6aKURLqwf7hyJHWNzfzo+YqoQxGRLkpJpAs7q7g7N5UN4rFFG9i4qzbqcESkC1IS6eK+cvkI0tOM781fFXUoItIFKYl0cSX5OXzm4qE8uXQL727eE3U4ItLFKImkgM9NOYuCvEz+bZ5aIyLSvpREUkDP3Ezu+MhwXl5dxYurdpx8BxGRU6QkkiJmXFjK0N7duPepFdQ3Nkcdjoh0EUoiKSIrI43/e91o1lYdYM5r66MOR0S6CCWRFPKRc/rwkbOLeWDBGqr26QJEETlzSiIp5v9cN5q6xia+O++9qEMRkS5ASSTFDCvuzmcuHsoT5ZUs2VgTdTgikuSURFLQFy8bQd/8HL72u2U0NGmQXUTaTkkkBXXPzuDeaWN4b9s+PQFRRM5IIh+P+4iZ7TCzd+PKepnZfDNbE94LQ7mZ2QNmVmFm75jZBXH7zAz115jZzLjyCWa2LOzzQHhErpyiK8f05eqxffnhc2tYv/NA1OGISJJKZEvk58DUo8ruAha4+whgQVgHuBoYEV6zgAchlnSAe4BJwETgnpbEE+rMitvv6M+Sk/jGDWPISk/j7j8sI9WecCki7SNhScTdXwaqjyqeBswJy3OAG+PKH/WYhUCBmfUDrgLmu3u1u9cA84GpYVu+u7/usW+/R+OOJaeoJD+Hf776HP5SsYtfl1dGHY6IJKGOHhMpcfetAOG9TygfAGyKq1cZyk5UXtlKeavMbJaZlZtZeVVV1RmfRFfyqYmDmTS0F/c+tYJN1bpdvIicns4ysN7aeIa3obxV7j7b3cvcvay4uLiNIXZNaWnG924aB8Cdv36bpmZ1a4nIqevoJLI9dEUR3lvuBlgJDIqrNxDYcpLyga2USxsMLMzjnutH88a6ah5+VbO1RLqafXUNCetpyEjIUY9vLjATuD+8PxlXfoeZPU5sEH2Pu281s3nAv8QNpl8JfM3dq81sn5lNBhYBM4D/6MgT6Wo+MWEgz63czv+bt5oPjihmVL/8qEMSSZjmZqe+qZmGpmYampyGpmbqG49cb22uSWtzQM3AMNLSwruFMjPSzDCIvR9R/te6GEfUS4sdMBznr3Wb3TnY0MTB+ibqGpqoa2jmYEMTtfWN7K1rZHdtPdUH6qk5UE9NbQPVB+rZtreObXvq2H+okZL8bBZ9/fJ2/1smLImY2S+BDwO9zayS2Cyr+4EnzOw2YCPwyVD9z8A1QAVQC3waICSLbwFvhnr3unvLYP0XiM0AywWeDi9pIzPjXz56Llf94BW+9MslPHnHxeRldfRvDJFT5+7sOdjAzv2xL8/qA4fYdaCe6v317DpQT01tPfvrGtl/qJED9Y0cONQUWz7USG19U9ThJ0x+TgaF3bIozMvirOJuXDK8N/0LcuhfkJuQz7NUm9pZVlbm5eXlUYfRab2ypooZDy/io+sW8b1f30th1AsAAA9TSURBVIcNHgz33QfTp0cdmqSY5mZn69461lUdYPPuWjbvrmPr7oNs3VPHlt0H2bLnIHUNrd9xoUd27Iu0e3YG3bMz6JadTrfDy7FXTmYaWelpZB5+GVkZseWMtNAiiNPaN6W74y3vDs0OjsfeD5fFvR+nLv7XfZpDXYjfN9YyyctKJzszndyWV1Y6OZlp9MzNpCAvi4LcTDLS23+UwswWu3tZa9v0U1OO8ME3nuVLi+bzw0mfZPLYy7lp2XyYNSu2UYlEEqC2vpFV2/ZRsWM/63YeOOJ1KO7ZN2bQp0c2/QtyGdUvn8tG9aFvz1x6d8+iV7fYq6hbNoXdMsnOSI/wjFKLWiJypNJSmjZuYsZN91I+YBR/+MWdjKpaD0OGwPr1UUcnSczd2bz7ICu37uO9rXtZuW0vK7fuY/2uA4fHHzLSjMFFeQzr3Z1hxd0Y2rsbpUXdGNQrl5L8HDIT8CtbTu5ELRElETlSWhq4U5VXwLV/+0NyGw/xh0fvpPDQfmjWzRrl1O2ta+CdTXtYsrGGpZt2s3TTbnYdqD+8fUhRHqP65jOqXz7n9OvByJIeDCrMTUh3jJwZdWfJqRs8GDZsoLh2Nw/+4Tvccst3uH3aXcxZ9BCZUccmndqm6loWratm0dpdLNm0m/er9h9uYQzv052PnNOHcYMKGN0vn3P69qBbtr5+ugL9V5Qj3XdfbAyktpYJW97j/mce4B+uu5NvXvD/+HbUsUmn0pI0Fq7dxcK1u6isOQhAYV4m4wcXMm1cf84fXMB5AwvomaufIF2VkogcqWXw/O67YeNGPrZ/LauL6/hJVQ4jXlvPzItKIw1PorO3roHXKnby0uoqXlmz84ikMXlYEX/3wWFMHlbEiD7dSUvTTbVThZKIHGv69CNmYv1Ts1Pxi8V844/LKe6RzTXn9oswOOkozc3Oiq17eWl1FS+tqmLxxhqamp0e2RlcNFxJQ2KUROSk0tOM/7hlPP/r4UV85fGlFORmctHw3lGHJQlQW9/Iy6ureHbFdl5eXcXO/bGB8LED8vnClLOYcnYx5w8q0CwpOUyzs+SU7alt4Kb/ep3Kmloen3Uh5w7sGXVI0g527j/EgpXbeXb5dl6p2El9YzMFeZlMGVnMlJHFfHBEMcU9sqMOUyKkKb5xlETOzLY9dXz8wdeorW/ksc9OZnR/3WMrGa3feYBnV2zj2eXbWbyxBncYUJDLlWNKuHJ0Xz5QWqiptnKYkkgcJZEzt2HXAW6evZCDDU38922TGDtALZJksLZqP0+9s5U/vbOVVdv3ATC6X/7hxDGqXw/0lGlpjZJIHCWR9rFxVy23/HQh+w818t+3TVLXVie1YdcBnnpnK0+9s5WVW/diBh8Y0oupY/tyxegSBvXKizpESQJKInGURNrPpupabp69kL0HG/ivGRO46CwNtncGm6pr+dOyWItj2eY9AFwwuIDrzuvPNef2o2/PnIgjlGSjJBJHSaR9bdl9kL/92Rus23mA7910PjeM6x91SClpy+6D/OmdrTy1bCtvb9oNwLhBBVx3bj+uOa8fAxJ0G3BJDbrtiSRM/4Jcfv25i/i7X5TzpV8uYcvug3zuQ8PUt94Btu+tiyWOd7bw1sZY4jh3QE/uuvocrj23n7qqpEMoicgZ65mXyaOfmcidv36b+59+j5Vb93L/x84jN0u3425vNQfqefrdbfzx7S0sXLcLdxjVL59/uupsrj23H6W9u0UdoqQYJRFpFzmZ6fzHzeMZ1bcH35u/mtXb9zP71gn6NdwO9tU1MH/Fdua+vYVX1+yksdkZVtyNL182guvO68/wPt2jDlFSWNKPiZjZVOCHQDrwkLvff6L6GhNJvBdW7eDLv1wCwH0fPZfrNU5yYo89dvheZYQnSdbddDPPv7eDuUu38PyqHdQ3NjOgIJfrx/Xn+nH9GN0vX12G0mG67MC6maUDq4ErgEpiz2K/xd1XHG8fJZGOsWHXAb78+FKWbtrNR8cP4JvTxpCfozu5HuOxxw7fNbk+LYNXho7nj2MvZf7oSzjQbBT3yObac/tx/bj+XDC4QIlDItGVk8iFwDfc/aqw/jUAd//O8fZREuk4jU3N/OiFCv7j+QpKemTzzWljuWJ0SdRhdSr1w87iL1bIMyMv4pmRF7Entwc9D+7jmi1vc/137mTSsCLSdXNDiVhXTiKfAKa6+2fD+q3AJHe/46h6s4BZAIMHD56wYcOGDo81lS3ZWMNdv13Gqu37uGJ0Cd+4YUxKTzmta2jipdVVPPPuNp57fRX7crrT/VAtl1cs4oYVL3HJ+qVkeZOeJCmdRlee4tvaT7RjsqK7zwZmQ6wlkuig5EjjBxfy1Jcu4ZFX1/GD59Zw2fde5NMXD+XzU85KmYcV7T/UyAvv7eCZd7fxwqod1NY30TM3k6u2LOPqt57lkvVLyG5q/OsOQ4ZEF6zIaUj2JFIJDIpbHwhsiSgWOYHM9DQ+N+Usrj2vH9+dt4oHX3yf/1m0kc9POYtbLxxC9y74qNSNu2p5/r3tLHhvB4vWVlPf1Ezv7lncOH4AV4/ty+RhRWQ+XgV/+j7EJ5C8vNgTJkWSQLJ3Z2UQG1i/DNhMbGD9U+6+/Hj7aEykc1i+ZQ/fnbeKF1dV0SMng+mThvC3F5Um9S05GpqaKV9fwwurdrBg5XberzoAwLDiblx6dh+uGF1CWWmvY8c4WpmdFf9QMJGoddkxEQAzuwb4AbEpvo+4+wl/wimJdC5LNtbw0CvrePrdraSnGVeMLuETEwbyoRHFZDz+y0795drU7KzcupfX3t/JXyp28eb6amrrm8hMNyYPK+IjZ/fh0nP66AJASXpdOomcLiWRzmnjrlrmvL6e3y/ZTPWBeoozmrl28TwuX/kqEzctJ6u5MdbNM3t2ZInkYH0TyzbvYcnGGhZvqGHRumr2HGwAYHif7lx0VhEXndWbS0b07pLdc5K6lETiKIl0bvWNzbywage//c7DvFQyikOZ2fQ4dICL1y9l4qbllDXVMLr8xRM/MKmle+joWXijR8Py4/Z0HmHPwQbWbN/Hqu37WLl1L0s27ua9bftoao79exncK4/Jw3px0Vm9ufCsIkryk7cbTuRklETiKIkkibQ0DqZn8WrpOBYMn8QrpeezuWfsGpPczHRGlnRneJ8ejCjpTmlRN0rysynJz6H4qd+R+bnYxXutGj2a5mXvsreugd21Dew6cIjKmoNxr1oqduxn6566w7v0yM7gvEE9GT+okPGDCzh/UAFF3fW4WEkdSiJxlESSRGnpMS2JrT2KKB/3IZZ89dus3r6P1dv3sWPfoWN2zWmoI6/hELkNdWQ3NtCUlkZjWjpNaenUZWSxJy+f1v6379UtiwEFuQzv052RJT04u2/sfUBBrq4Ul5TWla8Tka7qvvsO3w6kRb+mg1z/+Y9z/fWjD5ftqW1gU00t2/fWsX3vIXbceRcHsnI5mJFNbWYO9RmZZDY1ku5NZDQ3k91YT8+v/gMFuZkU5GVS2C2LgQW5DCjMJS9L/xxETpf+1Ujn1DJ4fpLZWT3zMumZ1/Ovz3mvfO3YsZCjzf9JAgIWSU0nGJ0Uidj06bB+fez2H+vXn9qsrPvui83iOp7Ro4+/TUROm5KIdC3Tp8emAbd225DTmJ0lIqdGSUS6npYWjPuRLyUQkXanJCIiIm2mJCIiIm2mJCIiIm2mJCIiIm2mJCIiIm2Wcrc9MbMqINmej9sb2Bl1EB1M55wadM7JYYi7F7e2IeWSSDIys/Lj3bemq9I5pwadc/JTd5aIiLSZkoiIiLSZkkhymB11ABHQOacGnXOS05iIiIi0mVoiIiLSZkoiIiLSZkoiScDM/tHM3Mx6h3UzswfMrMLM3jGzC6KOsb2Y2XfN7L1wXr83s4K4bV8L57zKzK6KMs72ZGZTwzlVmNldUceTCGY2yMxeMLOVZrbczL4cynuZ2XwzWxPeC6OOtb2ZWbqZLTGzp8L6UDNbFM75V2aWFXWMZ0JJpJMzs0HAFcDGuOKrgRHhNQt4MILQEmU+MNbdzwNWA18DMLPRwM3AGGAq8J9mlh5ZlO0knMOPif03HQ3cEs61q2kE7nT3UcBk4PZwnncBC9x9BLAgrHc1XwZWxq3/K/D9cM41wG2RRNVOlEQ6v+8DXwXiZ0BMAx71mIVAgZn1iyS6dubuz7p7Y1hdCAwMy9OAx939kLuvAyqAiVHE2M4mAhXuvtbd64HHiZ1rl+LuW939rbC8j9iX6gBi5zonVJsD3BhNhIlhZgOBa4GHwroBlwK/CVWS/pyVRDoxM7sB2Ozubx+1aQCwKW69MpR1NZ8Bng7LXfWcu+p5HZeZlQLjgUVAibtvhViiAfpEF1lC/IDYj8DmsF4E7I77oZT0/70zog4g1ZnZc0DfVjbdDXwduLK13VopS5q52ic6Z3d/MtS5m1gXyGMtu7VSP2nO+QS66nm1ysy6A78FvuLue2M/zLsmM7sO2OHui83swy3FrVRN6v/eSiIRc/fLWys3s3OBocDb4R/aQOAtM5tI7NfLoLjqA4EtCQ613RzvnFuY2UzgOuAy/+uFTEl9zifQVc/rGGaWSSyBPObuvwvF282sn7tvDV2yO6KLsN1dDNxgZtcAOUA+sZZJgZllhNZI0v/3VndWJ+Xuy9y9j7uXunspsS+bC9x9GzAXmBFmaU0G9rR0CSQ7M5sK/DNwg7vXxm2aC9xsZtlmNpTYpII3ooixnb0JjAgzdrKITR6YG3FM7S6MBTwMrHT3f4/bNBeYGZZnAk92dGyJ4u5fc/eB4d/vzcDz7j4deAH4RKiW9Oeslkhy+jNwDbHB5Vrg09GG065+BGQD80MLbKG7f97dl5vZE8AKYt1ct7t7U4Rxtgt3bzSzO4B5QDrwiLsvjzisRLgYuBVYZmZLQ9nXgfuBJ8zsNmIzED8ZUXwd6Z+Bx83s28ASYsk1aem2JyIi0mbqzhIRkTZTEhERkTZTEhERkTZTEhERkTZTEhERkTZTEhERkTZTEhERkTZTEhGJkJl9IDw7JcfMuoVnbYyNOi6RU6WLDUUiFq5czgFygUp3/07EIYmcMiURkYiFe2a9CdQBF3WF27lI6lB3lkj0egHdgR7EWiQiSUMtEZGImdlcYk80HAr0c/c7Ig5J5JTpLr4iETKzGUCju/9PeN76a2Z2qbs/H3VsIqdCLREREWkzjYmIiEibKYmIiEibKYmIiEibKYmIiEibKYmIiEibKYmIiEibKYmIiEib/X+weNLjMW5znQAAAABJRU5ErkJggg==\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": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}