{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0a574a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "\n",
    "# this is the example function that we try to minimize:\n",
    "# f(x) = 4 * (1/x^12 - 1/x^6)\n",
    "#\n",
    "# we consider only the domain x > 0 here, for negative x or zero we return positive infinity\n",
    "#\n",
    "def f(x):\n",
    "    if x>0:\n",
    "        invx6 = 1/x**6\n",
    "        return 4*invx6*(invx6 - 1)\n",
    "    else:\n",
    "        return np.inf  \n",
    "\n",
    "# this is f'(x), the first derivative of the function f(x)\n",
    "#\n",
    "def fprime(x):\n",
    "    if x>0:\n",
    "        invx6 = 1/x**6\n",
    "        return 24*invx6*(1 - 2*invx6)/x\n",
    "    else:\n",
    "        return math.nan  # this is NaN (\"not a number\"), reflecting that f'(x) is ill-defined here\n",
    "\n",
    "# this is f''(x), the second derivative of the function f(x)\n",
    "#\n",
    "def fsecond(x):\n",
    "    if x>0:\n",
    "        invx7 = 1/x**7\n",
    "        return 24*invx7*(26*invx7 - 7/x)\n",
    "    else:\n",
    "        return math.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07cf6ccf",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# this plot will show us how f(x) looks in the relevant region\n",
    "#\n",
    "import seaborn as sbn\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "absc = []\n",
    "ordn = []\n",
    "size = 1000\n",
    "init = 1\n",
    "end = 1.35\n",
    "\n",
    "for i in range(size+1):\n",
    "    x = init + i*(end-init)/size\n",
    "    absc.append(x)\n",
    "    ordn.append(f(x))\n",
    "    \n",
    "frame, axes = plt.subplots()\n",
    "\n",
    "sbn.lineplot(x=absc, y=ordn, color=\"b\", ax=axes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26f11bfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# hill climbing, first variant, where in each step\n",
    "# we aim at having y = f(x) decrease by the amount delta\n",
    "#\n",
    "# whenever this fails, we reduce delta by the factor alpha and try again\n",
    "#\n",
    "# this stops when delta has become very small (controlled by choice of epsilon)\n",
    "\n",
    "x, delta = 1, 0.25\n",
    "epsilon, alpha = 1.0e-06, 0.125\n",
    "\n",
    "y = f(x)\n",
    "i = 0\n",
    "print(0, x, y, \"initial\", sep=\"\\t\")\n",
    "\n",
    "x_accepted = [x]\n",
    "y_accepted = [y]\n",
    "x_rejected = []\n",
    "y_rejected = []\n",
    "\n",
    "while delta > epsilon:\n",
    "    i += 1\n",
    "    yprime = fprime(x)\n",
    "    if yprime == 0:\n",
    "        print(\"plateau reached; terminating\\n\")\n",
    "        break\n",
    "    xshift = -delta/yprime\n",
    "    xtest = x + xshift\n",
    "    ytest = f(xtest)\n",
    "    print(i, xtest, ytest, sep=\"\\t\", end=\"\\t\")\n",
    "    if ytest < y:\n",
    "        x, y = xtest, ytest\n",
    "        print(\"accepted\")\n",
    "        x_accepted.append(xtest)\n",
    "        y_accepted.append(ytest)\n",
    "    else:\n",
    "        delta *= alpha\n",
    "        print(\"rejected\")\n",
    "        x_rejected.append(xtest)\n",
    "        y_rejected.append(ytest)\n",
    "\n",
    "print(\"\\nFinal outcome: x =\", x, \", y =\", y)\n",
    "print(\"Compare the analytic solution: x =\", 2**(1/6), \" y = - 1\")\n",
    "print(\"\\nDeviation in x:\", abs(x - 2**(1/6)), \"\\nDeviation in y:\", y+1)\n",
    "\n",
    "# redraw the line plot for f(x)\n",
    "#\n",
    "frame, axes = plt.subplots()\n",
    "sbn.lineplot(x=absc, y=ordn, color=\"b\", ax=axes)\n",
    "\n",
    "# add accepted and rejected points to the diagram\n",
    "#\n",
    "sbn.scatterplot(x=x_accepted, y=y_accepted, color=\"g\", ax=axes)\n",
    "sbn.scatterplot(x=x_rejected, y=y_rejected, color=\"r\", ax=axes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e80eb3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# hill climbing, second variant, where in each step\n",
    "# we change the value of x by the amount +delta or -delta\n",
    "#\n",
    "# whenever this fails, we reduce delta by the factor alpha and try again\n",
    "#\n",
    "# this stops when delta has become very small (controlled by choice of epsilon)\n",
    "\n",
    "x, delta = 1, 0.25\n",
    "epsilon, alpha = 1e-04, 0.125\n",
    "\n",
    "y = f(x)\n",
    "i = 0\n",
    "print(0, x, y, \"initial\", sep=\"\\t\")\n",
    "\n",
    "x_accepted = [x]\n",
    "y_accepted = [y]\n",
    "x_rejected = []\n",
    "y_rejected = []\n",
    "\n",
    "while delta > epsilon:\n",
    "    i += 1\n",
    "    yprime = fprime(x)\n",
    "    if yprime < 0:\n",
    "        xtest = x + delta\n",
    "    elif yprime > 0:\n",
    "        xtest = x - delta\n",
    "    else:\n",
    "        print(\"plateau reached; terminating\\n\")\n",
    "        break\n",
    "    ytest = f(xtest)\n",
    "    print(i, xtest, ytest, sep=\"\\t\", end=\"\\t\")\n",
    "    if ytest < y:\n",
    "        x, y = xtest, ytest\n",
    "        print(\"accepted\")\n",
    "        x_accepted.append(xtest)\n",
    "        y_accepted.append(ytest)\n",
    "    else:\n",
    "        delta *= alpha\n",
    "        print(\"rejected\")\n",
    "        x_rejected.append(xtest)\n",
    "        y_rejected.append(ytest)\n",
    "\n",
    "print(\"\\nFinal outcome: x =\", x, \", y =\", y)\n",
    "print(\"Compare the analytic solution: x =\", 2**(1/6), \" y = - 1\")\n",
    "print(\"\\nDeviation in x:\", abs(x - 2**(1/6)), \"\\nDeviation in y:\", y+1)\n",
    "\n",
    "# redraw the line plot for f(x)\n",
    "#\n",
    "frame, axes = plt.subplots()\n",
    "sbn.lineplot(x=absc, y=ordn, color=\"b\", ax=axes)\n",
    "\n",
    "# add accepted and rejected points to the diagram\n",
    "#\n",
    "sbn.scatterplot(x=x_accepted, y=y_accepted, color=\"g\", ax=axes)\n",
    "sbn.scatterplot(x=x_rejected, y=y_rejected, color=\"r\", ax=axes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adc393b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Newton's method tries to find a zero of the first derivative\n",
    "# Accordingly, let us look at a plot of f'(x)\n",
    "\n",
    "absc_prime = []\n",
    "ordn_prime = []\n",
    "size = 1000\n",
    "init = 1\n",
    "end = 1.35\n",
    "\n",
    "for i in range(size+1):\n",
    "    x = init + i*(end-init)/size\n",
    "    absc_prime.append(x)\n",
    "    ordn_prime.append(fprime(x))\n",
    "    \n",
    "frame, axes = plt.subplots()\n",
    "\n",
    "sbn.lineplot(x=[init, end], y=[0, 0], color=\"black\", ax=axes)\n",
    "sbn.lineplot(x=absc_prime, y=ordn_prime, color=\"b\", ax=axes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a152197d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Newton's method, where in each step we aim to reach the point where f'(x) = 0,\n",
    "# or in other words, we aim at a change in f'(x) by the amount -f'(x)\n",
    "#\n",
    "# whenever this fails, we reduce a scaling factor delta\n",
    "# (initially 1) by the factor alpha and try again;\n",
    "# upon success, delta is set back to 1\n",
    "#\n",
    "# this stops when delta has become very small (controlled by choice of epsilon)\n",
    "\n",
    "x = 1\n",
    "epsilon, alpha = 1.0e-3, 0.125\n",
    "\n",
    "yprime = fprime(x)\n",
    "\n",
    "print(0, x, yprime, \"initial\", sep=\"\\t\")\n",
    "\n",
    "x_accepted = [x]\n",
    "yprime_accepted = [yprime]\n",
    "x_rejected = []\n",
    "yprime_rejected = []\n",
    "\n",
    "i, delta = 0, 1\n",
    "while abs(yprime) > epsilon:\n",
    "    i += 1\n",
    "    ysecond = fsecond(x)\n",
    "    if ysecond == 0:\n",
    "        print(\"plateau reached; terminating\\n\")\n",
    "        break\n",
    "    xshift = -delta*yprime/ysecond\n",
    "    xtest = x + xshift\n",
    "    yprimetest = fprime(xtest)\n",
    "    print(i, xtest, yprimetest, sep=\"\\t\", end=\"\\t\")\n",
    "    if abs(yprimetest) < abs(yprime):\n",
    "        x, yprime = xtest, yprimetest\n",
    "        delta = 1\n",
    "        print(\"accepted\")\n",
    "        x_accepted.append(xtest)\n",
    "        yprime_accepted.append(yprimetest)\n",
    "    else:\n",
    "        delta *= alpha\n",
    "        print(\"rejected\")\n",
    "        x_rejected.append(xtest)\n",
    "        yprime_rejected.append(yprimetest)\n",
    "\n",
    "print(\"\\nFinal outcome: x =\", x, \", y =\", f(x))\n",
    "print(\"Compare the analytic solution: x =\", 2**(1/6), \" y = - 1\")\n",
    "print(\"\\nDeviation in x:\", abs(x - 2**(1/6)), \"\\nDeviation in y:\", f(x)+1)\n",
    "\n",
    "# redraw the line plot for f(x)\n",
    "#\n",
    "frame, axes = plt.subplots()\n",
    "sbn.lineplot(x=[init, end], y=[0, 0], color=\"black\", ax=axes)\n",
    "sbn.lineplot(x=absc_prime, y=ordn_prime, color=\"b\", ax=axes)\n",
    "\n",
    "# add accepted and rejected points to the diagram\n",
    "#\n",
    "sbn.scatterplot(x=x_accepted, y=yprime_accepted, color=\"g\", ax=axes)\n",
    "sbn.scatterplot(x=x_rejected, y=yprime_rejected, color=\"r\", ax=axes)"
   ]
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}