ipynb • Lines: 1438{
"cells": [
{
"cell_type": "markdown",
"id": "reverse-interview",
"metadata": {},
"source": [
"# Differentiation in Python: Symbolic, Numerical and Automatic"
]
},
{
"cell_type": "markdown",
"id": "parental-conclusion",
"metadata": {},
"source": [
"In this lab you explore which tools and libraries are available in Python to compute derivatives. You will perform symbolic differentiation with `SymPy` library, numerical with `NumPy` and automatic with `JAX` (based on `Autograd`). Comparing the speed of calculations, you will investigate the computational efficiency of those three methods."
]
},
{
"cell_type": "markdown",
"id": "looking-barcelona",
"metadata": {},
"source": [
"# Table of Contents\n",
"- [ 1 - Functions in Python](#1)\n",
"- [ 2 - Symbolic Differentiation](#2)\n",
" - [ 2.1 - Introduction to Symbolic Computation with `SymPy`](#2.1)\n",
" - [ 2.2 - Symbolic Differentiation with `SymPy`](#2.2)\n",
" - [ 2.3 - Limitations of Symbolic Differentiation](#2.3)\n",
"- [ 3 - Numerical Differentiation](#3)\n",
" - [ 3.1 - Numerical Differentiation with `NumPy`](#3.1)\n",
" - [ 3.2 - Limitations of Numerical Differentiation](#3.2)\n",
"- [ 4 - Automatic Differentiation](#4)\n",
" - [ 4.1 - Introduction to `JAX`](#4.1)\n",
" - [ 4.2 - Automatic Differentiation with `JAX` ](#4.2)\n",
"- [ 5 - Computational Efficiency of Symbolic, Numerical and Automatic Differentiation](#5)"
]
},
{
"cell_type": "markdown",
"id": "101116ab",
"metadata": {},
"source": [
"<a name='1'></a>\n",
"## 1 - Functions in Python"
]
},
{
"cell_type": "markdown",
"id": "bc6140d8",
"metadata": {},
"source": [
"This is just a reminder how to define functions in Python. A simple function $f\\left(x\\right) = x^2$, it can be set up as:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "d07a15ef",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"9\n"
]
}
],
"source": [
"def f(x):\n",
" return x**2\n",
"\n",
"print(f(3))"
]
},
{
"cell_type": "markdown",
"id": "06330bd1",
"metadata": {},
"source": [
"You can easily find the derivative of this function analytically. You can set it up as a separate function:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "1ff4ffb5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"6\n"
]
}
],
"source": [
"def dfdx(x):\n",
" return 2*x\n",
"\n",
"print(dfdx(3))"
]
},
{
"cell_type": "markdown",
"id": "8301af3f",
"metadata": {},
"source": [
"Since you have been working with the `NumPy` arrays, you can apply the function to each element of an array:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "9f5831d8",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x: \n",
" [1 2 3]\n",
"f(x) = x**2: \n",
" [1 4 9]\n",
"f'(x) = 2x: \n",
" [2 4 6]\n"
]
}
],
"source": [
"import numpy as np\n",
"\n",
"x_array = np.array([1, 2, 3])\n",
"\n",
"print(\"x: \\n\", x_array)\n",
"print(\"f(x) = x**2: \\n\", f(x_array))\n",
"print(\"f'(x) = 2x: \\n\", dfdx(x_array))"
]
},
{
"cell_type": "markdown",
"id": "8428f910",
"metadata": {},
"source": [
"Now you can apply those functions `f` and `dfdx` to an array of a larger size. The following code will plot function and its derivative (you don't have to understand the details of the `plot_f1_and_f2` function at this stage):"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "5c255f4e",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAGFCAYAAABg2vAPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAA9hAAAPYQGoP6dpAABQ0UlEQVR4nO3deZyNdf/H8dcxGIYx1jFkX0LWUiQtwo1yd0fpLlFGsjWUpU0bpZLlVlmylK0QbUiLsqssSakoQvaM3YyxzIyZ6/fH5zdOk20wM9dZ3s/H4zy65pozZz5OM9f1nu/qcRzHQURERIJWDrcLEBEREXcpDIiIiAQ5hQEREZEgpzAgIiIS5BQGREREgpzCgIiISJBTGBAREQlyCgMiQc5xHOLj49GSIyLBS2FAJMgdPXqUiIgIjh496nYpIuIShQEREZEgpzAgIiIS5BQGREREgpzCgIiISJBTGBAREQlyCgMiIiJBTmFAREQkyCkMiIiIBDmFARERkSCnMCAiIhLkFAZERESCnMKAiIhIkFMYEPFhgwYN4rrrriM8PJzIyEhatWrFxo0b0z2nUaNGeDyedI9u3bq5VLGI+COFAREftnTpUmJiYli5ciXz588nOTmZZs2acezYsXTP69y5M3v27Dn9GDJkiEsVi4g/yul2AZw8CTlz2kNE0pk3b166jydPnkxkZCRr1qzh5ptvPn0+LCyMqKio7C5PRDLDqVP2yJPHtRLcbRkYPx7KlYP333e1DBF/ERcXB0DhwoXTnZ82bRpFixalRo0a9OvXj+PHj5/zNRITE4mPj0/3EBEXTZ8O5cvbPdEl7v45fugQ7N0LgwZBu3aQQ70WIueSmppKr169aNiwITVq1Dh9/v7776ds2bKULFmSX375haeeeoqNGzfyySefnPV1Bg0axIsvvphdZYvI+aSmwmuvQWwsHD7sWhkex3Ec1757XByUKQPx8TBrFrRq5VopIr6ue/fufPnll3z77beUKlXqnM9btGgRTZo0YfPmzVSsWPGMzycmJpKYmHj64/j4eEqXLk1cXBwFChTIktpF5BxmzYK77oKICNixA1z6HXT3T/GICOjRw45ffRVczCUivqxHjx589tlnLF68+LxBAKB+/foAbN68+ayfDw0NpUCBAukeIuICx7F7H9i90MXfRffb5R97zAZNrF4Nixa5XY2IT3Echx49ejBr1iwWLVpE+fLlL/g1a9euBaBEiRJZXJ2IXJaFC+GHHyBvXrsXusj9MBAZCZ0723FaQhIRAGJiYpg6dSrTp08nPDyc2NhYYmNjOXHiBABbtmxh4MCBrFmzhm3btvHpp5/y4IMPcvPNN1OrVi2XqxeR80q753XuDMWKuVqKu2MG0uzYARUr2tSKlSvh/5s5RYKdx+M56/lJkyYRHR3Nzp07ad++PevWrePYsWOULl2a1q1b89xzz2W4+T8+Pp6IiAiNGRDJTitXQoMGNq1+yxYbP+ci3wgDAB07wuTJcOedMHu229WIBA2FAREX3HknfPqp3fsmTnS7Gh8KAxs2wFVX2YCKdeugenW3KxIJCgoDItls3TqoWRM8Hvj9d6hSxe2KfGDMQJqqVW16BdicSxERkUCUdo+7+26fCALgSy0DAGvWwLXXQkgI/PEHVKjgdkUiAU8tAyLZ6M8/4corISXF7nnXXON2RYAvtQwA1K0LzZvbm6SNVkREJNAMHmz3uObNfSYIgK+1DAB8+y3cdBPkzm0jLC+wwIqIXB61DIhkk127rMU7OdnudQ0bul3Rab7VMgBw441wyy2QlATDhrldjYiISOYYOtSCQKNGPhUEwBdbBgAWLIB//ctWZdq6FYoXd7sikYCllgGRbLB3r+3Se/Kk3eOaNHG7onR8r2UA7E2qXx9OnIDXX3e7GhERkcszfLgFgeuvh8aN3a7mDL4ZBjweeO45Ox492rY6FhER8UcHD8Jbb9nxc8/ZPc7H+GYYAGjZEmrXhoQEGDHC7WpEREQuzYgRdi+rUwduv93tas7Kd8PA31sH3nwT4uPdrUdERORixcV5/6D10VYB8OUwALYiYbVqcOSIdReIiIj4k7fesntYtWrQurXb1ZyTb4eBHDngmWfsePhwa2YRERHxBwkJ8L//2fEzz9g9zUf5bmVp7rsPKlWCAwdgzBi3qxEREcmYt96ywYOVK9u9zIf5fhjImROefdaOhw6FY8fcrUdERORCjh3zLpz37LN2L/Nhvh8GANq1syUc9++HcePcrkZEROT8xo61e1aFCnYP83H+EQZy5fKOHRgyBI4fd7ceERGRczl+3LvZnh+0CoC/hAGABx+0pRz37oXx492uRkRE5OzGjYN9+6B8eXjgAberyRD/CQN/bx0YPNiWKhYREfElJ054WwWeecbuXX7Af8IAQIcOUKYMxMbCO++4XY2IiEh6b79t96gyZaxF20/4VxjInRv69bPj116zTR9ERER8wcmT1nIN1iqQO7e79VwE/woDAB07QqlS8Ndfah0QERHf8fbbdm8qXRqio92u5qL4XxgIDfWOHRg0SK0DIiLivhMn7J4Edo8KDXW3novkf2EA4KGHLHn99ZdmFoiIiPvGj4c9e2yswEMPuV3NRfPPMBAa6t3RcNAgzSwQERH3HD/ubRV47jm/GiuQxj/DAFh/TNmyNmpz7Fi3qxERkWA1dqytgVOunN+NFUjjv2Egd254/nk7fu017VkgIiLZ79gxuweB3ZP8ZF2Bf/LfMAA2h7NCBVvpSTsaiohIdnvrLe8eBH6y2uDZ+HcYyJXLO3Zg8GDbO1pERCQ7JCR4Vxv041YB8PcwAJbEKlaEAwdg1Ci3qxERkWAxapTdeypVgvbt3a7msvh/GMiZE154wY6HDoX4eHfrERGRwBcX520VeOEFv9iZ8Hz8PwwA3H8/VKkChw7BG2+4XY2IiAS6N96Aw4ehalW7B/m5wAgDOXPCSy/Z8f/+Z6FAREQkKxw8aPcasHtPSIi79WSCwAgDAG3aQK1a1k0wdKjb1YiISKAaOhSOHoXateHuu92uJlMEThjIkQMGDrTjESNsAQgREZHMFBtr9xiwe06OwLiNBsa/Is0dd0C9erY0ZNoiECIiIpnltddsCfx69eDf/3a7mkzjcRzHcbuITDV/PjRrZvsXbN5s2x2LyDnFx8cTERFBXFwcBQoUcLscEd+1c6dNI0xKsntN06ZuV5RpAqtlAOx/zs03Q2IivPyy29WIXJZBgwZx3XXXER4eTmRkJK1atWLjxo3pnnPy5EliYmIoUqQI+fPn5+6772avuslEMt/LL1sQuOUWaNLE7WoyVeCFAY/HO3ZgwgT480936xG5DEuXLiUmJoaVK1cyf/58kpOTadasGcf+thdH7969mTt3Lh9++CFLly7lr7/+4q677nKxapEAtGULTJxoxwMH2r0mgAReN0Ga5s3h669tVaj33nO7GpFMsX//fiIjI1m6dCk333wzcXFxFCtWjOnTp9OmTRsANmzYQLVq1VixYgXXX3/9BV9T3QQiGdC+PUybZveWefPcribTBV7LQJpXX7X/TpsGv/7qbi0imSQuLg6AwoULA7BmzRqSk5Np+re+y6pVq1KmTBlWrFhx1tdITEwkPj4+3UNEzuPXX2H6dDtOu7cEmMANA3Xr2toDjuPdzEjEj6WmptKrVy8aNmxIjRo1AIiNjSV37twULFgw3XOLFy9ObGzsWV9n0KBBREREnH6ULl06q0sX8W/PPmv3knvugWuucbuaLBG4YQC8c0A//RSWL3e7GpHLEhMTw7p165gxY8ZlvU6/fv2Ii4s7/di5c2cmVSgSgL77DubOtVUG08ajBaDADgNVq0LHjnb8zDOW7ET8UI8ePfjss89YvHgxpf42XTYqKoqkpCSOHDmS7vl79+4lKirqrK8VGhpKgQIF0j1E5Cwcx+4dYPeSKlXcrScLBXYYAOjf39YcWLrUBhSK+BHHcejRowezZs1i0aJFlC9fPt3n69atS65cuVi4cOHpcxs3bmTHjh00aNAgu8sVCSxffQXLltk9pH9/t6vJUoEfBkqXhkceseN+/SA11d16RC5CTEwMU6dOZfr06YSHhxMbG0tsbCwnTpwAICIigk6dOtGnTx8WL17MmjVr6NixIw0aNMjQTAIROYfUVLtnAMTEBPwCdoE7tfDv9u+HChUgIQFmzoT//tftikQyxHOOucyTJk0iOjoasEWH+vbty/vvv09iYiLNmzfnrbfeOmc3wT9paqHIWcycCffdB+Hhtl5N0aJuV5SlgiMMAAwYAC++CJUrw/r1kCuX2xWJ+ASFAZF/SE6G6tVh0ya7dwR4FwEEQzdBmj59oFgx+587YYLb1YiIiK965x27VxQrZveOIBA8YaBAAXj+eTseMMC6DERERP4uIcFakQFeeMG6CYJA8IQBgK5dbezA3r3w+utuVyMiIr5m+HC7R1SsCF26uF1NtgmuMJA7N7zyih0PGWIDC0VERAD27YOhQ+34lVfsnhEkgisMgM0kqFvXmoK0xbGIiKR5+WW7N1x7rS09HESCZzbB3y1cCE2b2oyCDRus60AkSGk2gQi2RXG1ajaTYOFCaNzY7YqyVfC1DAA0aWLbUCYnaxMjERGxe0FyMrRoEXRBAIK1ZQBg7Vrbfcpx4IcfrOtAJAipZUCC3po11jXg8cBPP0Ht2m5XlO2Cs2UAoE4daNfOjp94QpsYiYgEI8eBxx+343btgjIIQDCHAbDBIqGhsHgxfPGF29WIiEh2+/xzWLLE7gVps82CUHCHgbJl4bHH7PiJJ+DUKXfrERGR7HPqlF37AXr1gjJlXC3HTcEdBsB2pSpSBH7/HSZOdLsaERHJLhMm2IyyIkW8OxQGKYWBggW9m1C88AIcPepqOSIikg2OHrVrPtgS9RERrpbjNoUBsGWKK1WyJSjTVp8SEZHANWSIrThYubLdA4KcwgDYkpODB9vxsGGwe7e79YiISNbZvRv+9z87HjxYW9qjMODVujU0bAgnTnh3NxQRkcDz/PN2rb/xRmjVyu1qfELwLjp0NitXQoMGtvDEjz/aWgQiAU6LDklQ+eknW2TOceyaX7++2xX5BLUM/N3118O999oPSZ8+WohIRCSQ/P3aft99CgJ/ozDwT4MHexci+vRTt6sREZHMMmeOLTCUJw+89prb1fgUhYF/KlsW+va148cfh6Qkd+sREZHLl5joXXa4b1+71stpCgNn8/TTEBUFmzfDqFFuVyMiIpdr1Cjbpjgqyq7xko7CwNmEh9u+BQAvvQQHDrhbj4iIXLr9+2HgQDt+5RXIn9/denyQwsC5REfbbIK4OFudSkRE/NOAAXYtr1MHOnRwuxqfpDBwLiEhMHy4HY8dC7/95m49IiJy8davh3Hj7Pj11+3aLmdQGDifW2+FO++ElBTo3VtTDUVE/EnaVMKUFFtcqFEjtyvyWQoDFzJsmC1X/PXX8NlnblcjIiIZNXeuXbtz59a+MxegMHAhlSpZqwDYfxMT3a1HREQuLDHRWgXA/lupkrv1+DiFgYx49lkoUcKmpbz5ptvViIjIhbzxhl2zS5SAZ55xuxqfpzCQEeHh3tWqBg6EPXvcrUdERM5tzx7v9PDXXrNruJyXwkBGtW8P9epBQgL06+d2NSIici5PP23X6vr17dotF6QwkFE5csCIEXY8ZQp8/7279YiIyJlWrYJ337XjESPs2i0XpHfpYtSv712w4tFHITXV3XpERMQrNdWuzWALx9Wr52o5/kRh4GINGmRLWf49fYqIiPvSWm3Dw+1aLRmmMHCxSpSA/v3t+Mkn4cgRV8sRERHg8GF46ik77t/fNiSSDFMYuBSPPQbVqtnmFy+84HY1IiLywgt2Ta5WzdtVIBmmMHApcuWCkSPtePRo+Plnd+sREQlma9fCW2/Z8ahRdo2Wi6IwcKmaNIF77rEBKz16aN8CERE3OI5dg1NT4b//hcaN3a7ILykMXI7//Q/CwuDbb2HaNLerkQC0bNky7rjjDkqWLInH42H27NnpPh8dHY3H40n3aNGihTvFirhh6lT47ju7Fg8b5nY1fkth4HKULg3PPWfHTzwB8fHu1iMB59ixY9SuXZvRo0ef8zktWrRgz549px/vv/9+NlYo4qK4OLv2Ajz/vF2T5ZLkdLsAv9enD0yaBJs2wYABMHy42xVJALntttu47bbbzvuc0NBQojRyWoLRgAGwdy9UruzdUE4uiVoGLldoqHdlwhEj4Jdf3K1Hgs6SJUuIjIykSpUqdO/enYMHD573+YmJicTHx6d7iPidn3/2DuQeOdKuxXLJFAYyQ4sWcPfdkJIC3bppZULJNi1atODdd99l4cKFDB48mKVLl3LbbbeRkpJyzq8ZNGgQERERpx+l1bQq/iY1Fbp3t2tumzbQvLnbFfk9j+NoGHym2LXL5rcmJMA770CnTm5XJAHG4/Ewa9YsWrVqdc7n/Pnnn1SsWJEFCxbQpEmTsz4nMTGRxMTE0x/Hx8dTunRp4uLiKFCgQGaXLZL53nkHOne21WA3bIArrnC7Ir+nloHMUqoUvPiiHT/5JBw44G49EpQqVKhA0aJF2bx58zmfExoaSoECBdI9RPzGgQPelQZfeklBIJMoDGSmnj2hZk04dMi20BTJZrt27eLgwYOUKFHC7VJEssZTT9k1tlYtu+ZKplAYyEy5csGYMXY8YYLNfRW5DAkJCaxdu5a1a9cCsHXrVtauXcuOHTtISEjgiSeeYOXKlWzbto2FCxdy5513UqlSJZqrD1UC0XffwcSJdjxmDOTUhLjMojEDWeHhhy0M1KwJa9ZoaUy5ZEuWLOHWW28943yHDh0YM2YMrVq14qeffuLIkSOULFmSZs2aMXDgQIoXL57h7xEfH09ERITGDIhvS06Ga66BdevsGvv2225XFFAUBrLCwYNQpYr9d8gQ76IYIj5IYUD8wpAh1kVQpAhs3Gj/lUyjboKsUKQIDB1qx/37w9at7tYjIuLPtm61BYbAlhxWEMh0CgNZJToaGjWCEyfgkUe0kZGIyKVwHFtT4MQJuPVW6NDB7YoCksJAVvF4YOxYyJ0b5s2DDz5wuyIREf8zcyZ89ZWtMDh2rF1bJdMpDGSlKlXg2Wft+LHH4PBhd+sREfEnhw/btRPsWnrlle7WE8AUBrLaU09B1aq2mYbWHhARybinnoJ9+2x117SFhiRLKAxktdBQGD/ejsePh2+/dbceERF/8M033umD48dbl6tkGYWB7HDTTTYvFqBLF/jbuvAiIvIPiYnQtasdd+4MN97obj1BQGEguwwZAsWLw++/w6uvul2NiIjveuUVu1ZGRcHgwW5XExQUBrJLoUIwapQdv/oq/Pqru/WIiPiiX36BQYPseNQou3ZKllMYyE533w2tWsGpU7bF8Xn2nBcRCTpp18ZTp6B1a7tmSrZQGMhOHg+MHg0REbB6Nbz5ptsViYj4jjffhB9+sGtkWkuqZAuFgexWsqR3qeLnnoM//3S3HhERX7BlCzz/vB0PG2bXSsk2CgNuePhhW1bzxAkbKaulikUkmDmOzbRKW3K4Uye3Kwo6CgNu8Hhs3myePLBoEbzzjtsViYi455137FqYN6+tLaAlh7OdwoBbKlWCgQPtuG9f2LHD3XpERNywY4ddA8GuiRUrultPkFIYcFPv3nD99XD0qDWRqbtARIKJ41hX6dGjdi3s1cvtioKWwoCbQkJg4kRbsvirr2DSJLcrEhHJPhMnwtdf2zVw0iS7JoorFAbcVq0avPSSHffuDbt2uVuPiEh22LkT+vSx44EDbUM3cY3CgC/o0wfq1YP4eHUXiEjgS5s9EB8P9et7Q4G4RmHAF+TMaU1kuXPDl1/ClCluVyQiknUmT4Z589Q94EMUBnzFVVd5uwt69bImNBGRQLNzp3eg4EsvWVepuE5hwJf07WtNZnFxtuiGugtEJJCkpsJDD6l7wAcpDPiSnDmtiyBvXpg/H8aMcbsiEZHMM2YMLFhg17h337VrnvgEhQFfU6UKvPaaHT/xBGze7G49IiKZYdMmu6YBDB4MV17pbj2SjsKAL+rRw9bnPn4cOnTQVsci4t9SUuxaduIENG4MMTFuVyT/oDDgi3LksBG24eGwfLnt4CUi4q+GDoUVK6BAAbu25dCtx9fo/4ivKlvW9vYGeOEF+OUXd+sREbkUv/xi1zCwa1qZMu7WI2elMODLoqPhP/+BpCRo3x5OnnS7IhGRjDt5Etq1g+Rku5Z16OB2RXIOCgO+zOOx7TwjI+HXX+HZZ92uSEQk4555Btats2uYtib2aQoDvi4yEiZMsOPhw21ajoiIr5s/H15/3Y4nTbJrmfgshQF/8O9/Q/fudtyhAxw65G49IiLnc/CgdXMCPPII3H67q+XIhSkM+Ithw2wNgr/+gq5dtTqhiPgmx7Fr1F9/2U6EQ4e6XZFkgMKAvwgLg2nTbMWujz6y1btERHzNlCnw8cd2rZo2za5d4vMUBvxJ3brezYx69NDqhCLiWzZvhp497XjgQLjmGnfrkQxTGPA3Tz4Jt9wCCQlw//027VBExG1JSdC2rV2bbrnFu/Sw+AWFAX8TEgJTp0LhwrB6NTz/vNsViYjAc8/BDz/YtWnqVLtWid9QGPBHpUp5pxsOGQJff+1uPSIS3L76yjtQcOJEu0aJX1EY8FetWtmUHYAHH4S9e10tR7LGsmXLuOOOOyhZsiQej4fZs2en+7zjOLzwwguUKFGCvHnz0rRpUzZt2uROsRKc9u61axDYBkR33uluPXJJFAb82bBhUKOG/TJGR0NqqtsVSSY7duwYtWvXZvTo0Wf9/JAhQxgxYgRjx45l1apV5MuXj+bNm3NSS1dLdkhNtbVP9u2DmjU1jdCPeRxHE9b92vr1cO21tgb40KHw+ONuVyRZxOPxMGvWLFq1agVYq0DJkiXp27cvj/////e4uDiKFy/O5MmTue+++zL0uvHx8URERBAXF0eBAgWyqnwJRMOG2UDBvHltvMBVV7ldkd86cQJCQ93b0FEtA/6uenV44w077tcPVq50tRzJPlu3biU2NpamTZuePhcREUH9+vVZsWLFOb8uMTGR+Pj4dA+Ri7ZihV1zwJYdVhC4JDt2wNNPQ+nSMG+ee3UoDASCLl3g3nvh1Cn7r5YrDgqxsbEAFC9ePN354sWLn/7c2QwaNIiIiIjTj9KlS2dpnRKADh3yXnPuu8+uQZJhjgOLF8Ndd0H58jB4sK3gPHOmezUpDAQCjwfGj4dKlSxmRkdruWI5p379+hEXF3f6sXPnTrdLEn/iODZOYOdOqFwZxo3TboQZdOyYvV21akHjxjBrlg27SDueONG92hQGAkWBAvDBB5A7N8yd690tTAJWVFQUAHv/MZNk7969pz93NqGhoRQoUCDdQyTDhg+Hzz6zDu4PPrBrj5zXn39C374247JbN9vVOSzMe7xwoU0Qc3NpBoWBQHL11d7xA089pfEDAa58+fJERUWxcOHC0+fi4+NZtWoVDRo0cLEyCVgrVlgHN9i1pk4dN6vxaY5juzj/5z/WaDt8OBw5AhUq2PHu3TBmjA378gU53S5AMlm3brBkiSX2e++FH3+EIkXcrkouUUJCApv/tgfF1q1bWbt2LYULF6ZMmTL06tWLl19+mcqVK1O+fHmef/55SpYseXrGgUimOXjQxgekjU3q2tXtinzS0aO2V9OoUbBxo/d88+a2bcNtt7k3Y+C8HAk8cXGOU6mS44DjtGjhOKdOuV2RXKLFixc7wBmPDh06OI7jOKmpqc7zzz/vFC9e3AkNDXWaNGnibNy48aK+R1xcnAM4cXFxWfAvkIBw6pTjNG9u15RKlewaI+ls3Og4PXs6Tni4vU1gxz17Os6GDW5Xd2FaZyBQ/fILXH+9TV4dMAD693e7IvFRWmdALmjAAHjxRVtPYOVKGwEnpKbadMCRI9NPC6xSxTaWffBB/xlSoTAQyN57z34aPR744gto0cLtisQHKQzIeX35JbRsaX/svvcetG/vdkWui4uDSZNg9GjvTvIej71NPXtC06Y+2hVwHgoDga57dxg71nYSW7MGypVzuyLxMQoDck7btsE118Dhw3Yteesttyty1W+/WSvAe+/ZNEGAiAjo1Mm2iqlY0d36LofCQKBLTISbbrLtjuvWhW+/hTx53K5KfIjCgJzVyZNw4432R0S9erBsmU0nDDIpKTZbe+RIWLTIe756dWsFaN8e8uVzr77M4mcNGXLRQkPho49sRsGaNdaRpfwnIufjOHatWLPGrh0ffhh0QeDQIdvupWJFaN3agkCOHN7jX3+1CRWBEARAUwuDQ5kyMH26jRmYMME2NurWze2qRMRXjR1r14ocOezaUaaM2xVlm19+sVaAadNs/DVYL2vnztZTUrasu/VlFXUTBJPBg23BkFy5LNreeKPbFYkPUDeBpPPNN7Y+7qlTds148km3K8pyyckwe7aFgG++8Z6vU8e6Atq2tYkUgUxhIJg4ji0a8sEHULy4bTlaqpTbVYnLFAbktF27bGzRvn22sND77wf0vgP798Pbb9u4yN277VxICNx9t4WAhg0D+p+fjsJAsDl2DBo0sA6vevVg6VINKAxyCgMC2IDBm2+2wca1asHy5YHTIf4Pa9ZYK8CMGTbGGqBYMRsD0K0bXHGFu/W5QQMIg02+fNYeVqgQfP89xMRoQKFIsHMcmxu3erV1kM+aFXBBICnJGjpuuMGGTU2ZYkHguuvg3XdtE8aBA4MzCIAGEAanChVs4+wWLWzPzNq14dFH3a5KRNwyYoStopMjh10bKlRwu6JMs2ePbRs8bhzExtq5XLngv/+1roD69d2tz1eomyCY/e9/8PjjdgH44gvbSUOCjroJgty8ebZ0XmoqDBtme+36OceBVausK+DDD22AIECJEtYN0KULnGeX76CkMBDMHMeWzpo0yZbRWrkSqlZ1uyrJZgoDQWzDBvvTOD4eOna06YR+PGIuMdEaNkaOtPHRaW64wVoB7roLcud2rz5fpjAQ7BITbSHtb7+1TbdXrbI+QwkaCgNB6tAhCwKbN9s04wUL/HZhoV27bGmE8eNthgDYP6VtWwsB11zjbn3+QGFAbBpRvXqwfTs0aWIbk+TK5XZVkk0UBoJQcrKNGVq0yFbRWb3ahtP7Ecexv2FGjoRPPrFlg8FmSz/yCDz8sN/9k1yl2QQCkZHw6ac2enjhQnjsMc0wEAlUjmMDhhctgvz5beF9P7prnjhhvRlXX20zIT/80ILALbfYyutbt0K/fn71T/IJmk0gplYtW3a0VSsYM8a6DPr0cbsqEclsw4dbm7rHY2vu1qzpdkUZsn27LQ70zjvWwwG2KmC7dtYVUKuWu/X5O3UTSHrDh9toYo8HPv7YduWQgKZugiDyySfQpo21DgwfDr17u13ReTkOLF5sXQGffmoTHsB2Yo+JgYce0hCnzKIwIOml7Vb21lsWu5cssfEEErAUBoLE999Do0bWzh4TY3dYH505kJAAU6fCqFGwfr33fNOm1grQsqUtGyyZR2FAznTqFPznPzaQMDLSZhiUK+d2VZJFFAaCwLZtNnNg3z64/XaYMwdy+l4v8ZYtMHq0rYUWF2fn8uWDDh3sb5Rq1dytL5ApDMjZHT0KN90EP/8MV10F330HBQu6XZVkAYWBAHfkiO2489tvtg3fsmUQHu52VaelpsL8+dZQ8cUX3rHLlSpZA0bHjrYMimQthQE5t1277K+Jv/6yobrz5mlTowCkMBDATp60lUWXLbNF91et8pnF9+PjbX+AUaPgjz+852+7zboCmje3xVEle+itlnMrVcqieni47W744IPeETwi4ttSUuCBBywIFCgAn3/uE0Fg40a72V9xhc1w/OMPK++xx+z4iy8sECgIZC/f6zQS31K7tu1y2KKFTegtUQLeeMNnBx6JCNbW3ru3TbzPlct2Iaxd27VyUlPtJj9yJHz9tfd81aoWDB54wKd6LoKSspdcWOPGtscn2O5mw4a5W4+InN/QoXbnBfvdbdzYlTKOHLEZjFdeCXfcYUHA47HxyV9/bcMYHnlEQcAXqGVAMua++2zsQN++8OST1kLQvr3bVYnIP733Hjz1lB0PH26/u9ls3TobC/Dee3D8uJ0rVMj2RXvkEShfPttLkgtQGJCM69PHBhW+/roN8S1UyCb8iohv+Owz+90E+33NxkWFTp2ylY1HjrSFgtLUrGldAe3aQVhYtpUjF0mzCeTipKbaQMJp02xmwddf2xRE8VuaTRAgli2zIfgnT9qd9913s2UU3sGDtkTwW2/Bjh12LkcOW7y0Z0/bP0BDjHyfwoBcvORk+03//HMbBrx4sfYI9WMKAwHgxx/h1lttvt6//23LDmfxzqNr11orwPTplj8AihSBLl2gWzcoUyZLv71kMoUBuTQnTtgMg2XLbHuwb76BKlXcrkougcKAn9uwwVrnDhyw9UC+/NKWEs8Cyck2MWHkSNs+OM0111grwH33aSkSf6UwIJcuLs5GKf/4I5QubVcH/TngdxQG/NiOHXDjjbBzp92RFy+21rpMtm8fjB9vmx3u3m3ncua0PY969oQGDdQV4O8UBuTy7N9vf5Vs3Gjrhy5dCiVLul2VXASFAT/111/WIb9li7XKffONtdJlotWrrRVg5kxISrJzkZHWDdC1q37VA4lmE8jlKVbMFha/+WbYvBmaNLFAEBnpdmUigWvvXvtd27LF5unNn59pQSApydYXGznSVi9OU7++tQK0aQOhoZnyrcSHqGVAMsfWrRYIdu2yuUSLF9toIvF5ahnwMwcO2GDBdeuse27ZskzZVXTPHusGGDfOsgZA7txw770WAq677rK/hfgwhQHJPJs2WSCIjbX+y4ULtdOhH1AY8CNHjtg4nZ9+soW/li6FypUv+eUcB1assFaAjz6ytQLAmv+7d4fOnaF48cwpXXybwoBkrt9+g0aNbCxB/frw1Vfaf9THKQz4ibg4W0dg1SrrEli6FKpVu6SXOnkSZsywEPDjj97zDRtaK8Bdd2X5zETxMQoDkvl++cWaMQ8dgnr1LBCohcBnKQz4gSNHLAh8/z0ULmzdcLVqXfTL7NwJY8bA229bbwNY///991sIuPrqzC1b/IfCgGSNtWuhaVNbnuzaa22lwkKF3K5KzkJhwMcdPgzNmsEPP9g4nIULL2oHQsexYQUjR9oGpCkpdr5MGdsnoFMnKFo0a0oX/6FdCyVr1KkDixbZVeaHHywYHDrkdlUBacCAAXg8nnSPqlWrul2WZIZDh+x354cf7Hdp0aIMB4Hjx60FoHZt67n7+GMLArfeagsUbtli+xkpCAhoaqFkpVq17OLVpIl1TDZpAgsWaJZBFqhevToLFiw4/XHOnPrV9nsHD1oQWLvWxggsWgQ1alzwy7ZutX0CJkywRgWwDYLat4cePWyyj8g/6YohWSttmmHjxnZRa9TIugxKlHC7soCSM2dOoqKi3C5DMsuePfCvf8H69bZmx6JFUL36OZ/uONZ7MGKEbVyY1vlbvjzExMBDD6mXTs5P3QSS9apXhyVLLACsW2fTD7dvd7uqgLJp0yZKlixJhQoVaNeuHTvSto87i8TEROLj49M9xIds326req5fb3P8liw5ZxBISLBWgOrVLTvMnWtB4F//gk8/tdm+ffsqCMiFaQChZJ8//7Sugm3bbLGUBQvgyivdrsrvffnllyQkJFClShX27NnDiy++yO7du1m3bh3h4eFnPH/AgAG8+OKLZ5zXAEIfsHGjdQ3s2mULCS1cCBUqnPG0TZtg9GiYNMk2KgTInx86dLCuAA0ZkYulMCDZa9cu+7NlwwZr/pw//5KmSMm5HTlyhLJlyzJ8+HA6dep0xucTExNJTEw8/XF8fDylS5dWGHDbzz/brIF9++xuvmABXHHF6U+nptos3ZEjbWPCNJUrWwCIjs6SPYokSGjMgGSvUqVssZTmzW0MwS23WNvmjTe6XVnAKFiwIFdeeSWbN28+6+dDQ0MJ1eLyvuXbb+GOO2w9gTp1bFzN/+81EBcHkydbS8CmTfZ0jwduu83WBmjWDHKow1cuk36EJPtFRtqgwhtusItfWgenZIqEhAS2bNlCCQ3S9A9z5tjvwJEj9juxeDEUK8bvv9vgvyuugF69LAhERNjxH3/A559DixYKApI59GMk7ihY0LoI/v1vWxu1dWt45x23q/JLjz/+OEuXLmXbtm0sX76c1q1bExISQtu2bd0uTS7k7bdt7d+TJ+GOO0iZN59PlxXkX/+Cq66ywYHHjtnxmDHWy/b667ZbuEhmUjeBuCcsDGbNso3RJ060XVH27IHnnrN2UMmQXbt20bZtWw4ePEixYsW48cYbWblyJcUyeW97yUSOAwMHQv/+ABxq/ygTaw5ndK0Qtm2zp+TIYT0HPXvazFz9SkhW0gBCcZ/jwPPPwyuv2Mddu8KoUaCFc7KFliPOZqdOWfv/+PH8Sg1GXj2RqRuu5cQJu9sXKgQPP2xLBWfCzsQiGaIwIL5j1Ch49FELBy1awMyZGh6dDRQGslF8PKfa3Mec+XkZSU+W0uj0p2rVslaA+++3RjOR7KQwIL5l9my7Gp44YasXfv65rUkgWUZhIHsc+Hk3bzf/iDF7W7OTMgCEhNhwmZ49bZ0hdQWIWxQGxPf88IN1lsbG2qqFc+dC3bpuVxWwFAay1o8/wsgBB3h/bn4SyQNA0YLJdHkkF926KeuKb1AYEN+0fTu0bGlLsoaFwbvvwt13u11VQFIYyHzJybZL4MiRsHy593zdPOvp+UoU9z5ShDx53KtP5J80tVB8U9my8N13tqLK8ePQpg0MGGDLsIn4qNhYeOkl+/Ft29aCQE6Sact0ltfrxeo9pejQR0FAfI/CgPiuiAgbM9Crl3384otwzz22O4uID1m1yrYILlPGZgvu2QNRoYfpzwB2UIbpvX+gwXfD8BSMcLtUkbNSN4H4h0mToFs3SEqyYdezZ9v+rHLZ1E1waRIT4YMPrCtg9Wrv+euvPknPAwNos3M4uXN7YNw42zhAxIcpDIj/WLHChl7v3WuTsadNswXa5bIoDFyc3bth7FgYP972FALInRvuuw96XvMd1754Bxw+DFFR8Mkn0KCBuwWLZIC6CcR/NGhgMw2uu84uti1bWptsSorblUmAcxzbS+jee20hoJdftiBwxRV2vHNbClPKvsC1vW60n8169ay5QEFA/IRaBsT/JCZC7962WDvYIMNp06BoUXfr8lNqGTi3Eyfg/fetK2DtWu/5m26ytQFatYJccQdsbYz58+2TjzwCw4eDdoYUP6IwIP7rvfds6eITJ2yy9owZtuubXBSFgTPt2GGbBL3zDhw8aOfy5IF27aBHD9tlGLDpAvfdBzt3Qt681nfQvr1bZYtcMoUB8W+//mrrD2zaZMu5vfQSPPWUHUuGKAwYx4ElS6wVYM4c7yzWsmXtj/1OnaBIkf9/ckoKvPaat5uqcmVbWKBmTbfKF7ksCgPi/+LjoXt3mD7dPm7c2FoNSpZ0ty4/Eexh4NgxmDrVtsZYt857vnFj6wq4445/ZMu//rK//hcvto/vv9+6rILwvZPAoTAggcFxYMoU2w3u+HEbPzB5sg0ylPMK1jDw558werTtnn3kiJ0LC4MHH7SugOrVz/JFn39u0wQPHLAnjx4NHTpoUwHxewoDElg2bLA+3J9/to+7dIH//Q/y53e3Lh8WTGHAcWyc38iRdl9Pu/pVrGg5smNHKFjwLF+YkAB9+9qYALBBAzNmQJUq2VS5SNbS1EIJLFWrwsqV3lULx4+H2rVtaWMJWkePWjdAtWrQvDl89pkFgbTjP/6wCSpnDQLffms/Q2lBoHdvW/NCQUACiFoGJHAtXmxNujt2WDPuk0/aksaa8pVOILcM/PGHhYDJky0QAISH249FTMwF7ucnT9oAwaFDLTmUKWMvdOutWV+4SDZTGJDAFhdnrQSTJ9vHVavafLGGDd2syqcEWhhITYUvv7SugK++8p6vUsXGAnToYIHgvL77Dh5+2LqdwNLDG2/YfhkiAUhhQILD7Nm2t8HevdZKEBMDr76agbtC4AuUMHDkiG1hMXo0bNli5zwe+Pe/bVZA06YZGOd39Cj062eLDDgOFC9uaw+3apXF1Yu4S2FAgsehQ/D443bHAFuoaOxYuP12d+tymb+Hgd9+s1aA996zaYJgf8B36mSZr0KFDL7QF19YYNy50z5+6CEYNsz2wRAJcAoDEnwWLLBZBlu32setW8Prr9vqMkHIH8NASgrMnWvjARYu9J6vXt1aAdq3h3z5Mvhi27dbV9Ls2fZx+fI2WLBp00yuWsR3aTaBBJ+mTW3lwr59bTWZWbNsmPmrr9q+B+KzDh2CIUNsKmDr1hYEcuSw40WL7H9r164ZDAKJifDKK/b/fvZs+1no29deREFAgoxaBiS4rVtno8qWLrWPK1e2VoLbbw+ahWT8oWXg55+tK2DaNBvkD1C4MHTubItPXlSjjuPYIgN9+tgy1gC33GLNDDVqZHrtIv5AYUDEcWwp48cfh9hYO9ekiS1WVLu2u7VlA18NA6dO2R/sI0fCsmXe83XqWFdA27a2N9BFWbvW/vpftMg+joqy/89t2wZN+BM5G4UBkTRxcdZV8MYbkJRkN4foaNuwPoD3OfC1MLB/v3XZjx0Lu3bZuZAQ24+qZ0+bFXrR9+2//oLnnrMppo4DuXPbOIFnntF0QREUBkTOtHWrTS+bOdM+zpvXuhKefNL2PAgwvhIG1qyBESNsld+kJDsXGWljPbt1gyuuuIQXPXDABhmMGmVbXYMtVz1oEJQrl1mli/g9hQGRc1m50pqUly+3j/Pnt6Vo+/Q5x7q1/snNMJCUBB99ZF0BK1d6z193HTz6KNxzzyUuGHnkCAwfbuM/EhLsXMOG1iVQv35mlC4SUBQGRM7HcWw5u+eeg59+snOFClkTc48eNorNz7kRBmJjYdw46wpIG6aRKxf897/WFXDJ9+tDh6wV4PXXvVsRXn21dfXcdpvGBYicg8KASEakptoUxBdesFVuwFoKuna1lgI/HlOQXWHAceyv/5EjrTUgOdnOlyhh3QBduth4vkvy11/WEjBunLcl4KqrYOBAm3eoECByXgoDIhcjJcXuZIMGebdJzp3bFrx/7DFb9cbPZHUYOHnShl+MHGnjAtLccIO1Atx1l72Fl2T9enjzTZgyxTvQoHZtG/PRpo2NPBSRC1IYELkUad0HgwbZFrdpmjSxzu6WLf3mRpRVYWDXLhgzBt5+22YIgPX/t21rIeCaay7xhVNSbJ2AESPSLz944402O6BFC7UEiFwkhQGRy/Xtt9ZHPXu2dSeALWnbtau1GFxy23f2yMww4Dj2dowcCZ98YvdtgFKl4JFHbCPAYsUu8cVjY60FYNw471LSacsP9uplYUBELonCgEhm2b7ddrt7+204fNjOhYTYtnmdOtkAtpw53a3xLDIjDBw/bus2jRrl7T0BW9ivZ0+4885L/KefOmUbCE2YYK0BaemiUCEbZHDRyw+KyNkoDIhktuPH4f337Qa2YoX3fFQU3HuvzXOvX99nmrIvJwxs22b5Z8IEG8gPtixD+/Y22aJWrUsoyHFg1SpbcGDmTO90A4AGDSxYtW0LYWGX8OIicjYKAyJZ6bff7E757ru2AE6acuUsGLRpY53nOdzbM+xiw4Dj2Gq+I0fazoFpPSPlytmWwQ89dAkzLlNT4ccfbXDmjBnWypKmWDF48EF74auuusgXFpGMUBgQyQ5JSfD113ajmz0bjh3zfq5kSRtweMcdNgAxm//izWgYSEiAqVOtK2D9eu/5pk2tK+Cix0weP24DAOfOhc8+gz17vJ/Llw9atbJWlGbNLmO6gYhkhMKASHY7ftz6v2fOhHnz0geDPHlsIFzjxvaoWzfLxxlcKAxs2QKjR8PEibZ9A9i9ukMH6wqoVi2D3+jUKZtbuGiRPb791rsFYdqLtmhhLSYtW6obQCQbKQyIuOnkSViyxP46njsXdu5M//nwcFtG9/rrbZxBvXpnbYMfPXo0Q4cOJTY2ltq1azNy5Ejq1auXoRLOFgZSU2H+fOsK+OIL6xoAqFTJAkB0dAb29zl0CL7/3lYaWrUKvvsOjh5N/5zSpeE//7FWkUaNLnHtYRG5XAoDIr7CcWyMweLF9pfzkiXeWQl/V7myjcyrUQNq1OCLHTto8/TTjB43jvr16/PGG2/w4YcfsnHjRiIjIy/4bf8eBqAAU6ZYV8Aff3ifc9tt1hXQvPlZhjecOgWbN8O6dd7HL7/Apk1nfrNCheym37gx3HqrjQHwkYGUIsFMYUDEV6Wk2Dy9FSu8f12f7QYLpHg8hJQtC+XL45Qvz5APP6Rms2bcHh1tW/9FRkKRItb0/o+O/bQw0KVLHO+/X+D0H+8FCjh0bJdMzD37qJzvL9i3zx5//WXz/NMeO3d6p/z905VXWotG/fo2E6B2bb9ZjEkkmGQoDDiOw9F/Nu+JSPY7dAjWroXff4f160n97TdO/PQT+S7iJZzQUJywMFLy5GN+0i28dbQ13yS1AXYCBbgyxya6eN7mvpSphHPsQi9n8ua1wQNXXWWPatWgTp2A2MhJJBCEh4fjOU8rXIbCQNpfDiIiIuJ/LjRbSC0DlyA+Pp7SpUuzc+fObN//PZjofb6wPXv2ULVqVebPn59uwODzzz/Pd999x6JFi+yE4/D72kTGjXOYMSuUEyet4z8iXzJ3NvyNd7++mY3vv09UZKSN6g8Ls10Zw8JshoP69S+bfp6zh97ns7tQy0CG5ix5PB69qWdRoEABvS/ZQO/zueXJk4eQkBASEhLSvUdHjhzhiiuuICysAHPn2oDAtFwAULOmDQi8/344fLg875aGsBtvpECpUi78K4KLfp6zh97ni+Pesmcictly585N3bp1Wfi33ftSU1P5+us1JCX1omJF2yJ40SIbt3f33TZJ4eefoXNnawQQEfG9XVNE5KL06dOHDh06cO211xIefhPPPPMXsbGrmTcvL2CTCNL29Cld2uViRcQnKQxcgtDQUPr370+oFkjJUnqfM+auu+7lq68i6Nq1MImJFYAKAFx9NTz6qK3omyfPub8+7f3V+5y19POcPfQ+XxqtMyDip/bts92Sx4yB3bvtXM6c1hXQsyfccEPGxv1lxhbGIuLf1DIg4mdWr7ZlgmfOtP2PwNYU6tYNuna1fY9ERC6GwoCIH0hKgg8/tBCwapX3fL161gpwzz1a1l9ELp3CgIgP27MHxo2DsWNh7147lyuXbezXs6eFARGRy6UwIOJjHMe2Ixg5Ej76yPYBAmv+79bNZgYUL+5ujSISWBQGRHzEyZMwY4aFgB9/9J5v2NBaAe66y1oFREQymxYdyiSJiYnUqVMHj8fD2rVr3S4noGzbto1OnTpRvnx58ubNS8WKFenfvz9JaaPn/NzOnfDMM7YGQMeOFgRCQ73H335r3QJZEQRGjx5NjRo1AGjcuDHff/995n+TIDdo0CCuu+46wsPDiYyMpFWrVmzcuNHtsgLea6+9hsfjoVevXm6X4hcUBjLJk08+SUkN484SGzZsIDU1lXHjxrF+/Xpef/11xo4dyzPPPON2aZfMcWDpUmjTBsqXh0GD4MABCwSDBsGuXTBxoq0VkFVmzpxJnz59ePrppwGoUaMGzZs3Z9++fVn3TYPQ0qVLiYmJYeXKlcyfP5/k5GSaNWvGsWMZ3BFSLtrq1asZN24ctWrVcrsU/+HIZfviiy+cqlWrOuvXr3cA56effnK7pIA3ZMgQp3z58m6XcdGOHXOc8eMdp2ZNx7FIYI9GjRzn448dJzk5+2qpV6+eExMT48TFxTmAc/jwYadkyZLOoEGDsq+IILRv3z4HcJYuXep2KQHp6NGjTuXKlZ358+c7t9xyi/PYY4+5XZJfUMvAZdq7dy+dO3fmvffeIywszO1ygkZcXByFCxd2u4wM27YNnngCSpWyAYC//mobAnbpAr/8AosX25iAnNk0iicpKYk1a9bQtGnT0+dy5MhB06ZNWbFiRfYUEaTi4uIA/Orn15/ExMTQsmXLdD/bcmEaQHgZHMchOjqabt26ce2117Jt2za3SwoKmzdvZuTIkQwbNsztUs7LcWDhQhsQOHeufQzWLRATAw89BIUKuVPbgQMHSElJofg/piUUL16cDRs2uFNUEEhNTaVXr140bNjw9FgNyTwzZszgxx9/ZPXq1W6X4nfUMnAWTz/9NB6P57yPDRs2MHLkSI4ePUq/fv3cLtkvZfR9/rvdu3fTokUL7rnnHjp37uxS5eeXkABvvQXVq8O//gWffmpBIO140ybo29e9ICDuiYmJYd26dcyYMcPtUgLOzp07eeyxx5g2bRp5zrcZh5yV9iY4i/3793Pw4MHzPqdChQr897//Ze7cuXj+tgB8SkoKISEhtGvXjilTpmR1qX4to+9z7ty5Afjrr79o1KgR119/PZMnTyZHDt/Ksps2wejRMGkSxMfbufz5oUMH6NEDqlZ1t76/S0pKIiwsjI8++ojGjRuf3pugZ8+eHDlyhDlz5rhdYsDp0aMHc+bMYdmyZZQvX97tcgLO7Nmzad26NSEhIafPpaSk4PF4yJEjB4mJiek+J+kpDFyGHTt2EJ921cduVs2bN+ejjz6ifv36lCpVysXqAsvu3bu59dZbqVu3LlOnTvWZX+rUVPjqK+sK+PJL7/nKlS0AREeDr+79U79+ferVq8crr7xCREQEhw8fpkaNGvTo0eP0DAO5fI7j0LNnT2bNmsWSJUuoXLmy2yUFpKNHj7J9+/Z05zp27EjVqlV56qmn1C1zARozcBnKlCmT7uP8+fMDULFiRQWBTLR7924aNWpE2bJlGTZsGPv37z/9uaioKFdqiouDyZOtJWDTJu/522+3BYKaNQMfa7g4Q58+fejQoQPVq1cHoHfv3hw7doyOHTu6XFlgiYmJYfr06cyZM4fw8HBiY2MBiIiIIG/evC5XFzjCw8PPuOHny5ePIkWKKAhkgMKA+Lz58+ezefNmNm/efEbIyu6Grd9/h1GjYMoUSJsmXqCADQaMiYFKlbK1nMty7733sn//fl555RUAfv31V+bNm3fGoEK5PGPGjAGgUaNG6c5PmjSJ6Ojo7C9I5CzUTSByASkp8Pnn1hWwYIH3/FVXWVfAAw/Y2AB/FR8ff3rMQAFf7dMQkSyllgGRczh8GCZMsJkBW7fauRw54I47rCugcWP429hRERG/pTAg8g+//mqtAFOnwokTdq5QIXj4YXjkEShXztXyREQyncKACLZN8Jw5FgKWLvWer1XLWgHuv99WDBQRCUQKAxLUDhyAt9+GMWNs90CAkBBo3dpCwE03qStARAKfwoAEpR9/tFkB06dDYqKdK1rU9gro1s12DxQRCRYKAxI0kpPh44+tK2D5cu/5unWtFeDee0GrmIpIMFIYkIC3dy+MGwdjx8KePXYuZ0645x6bGtiggboCRCS4KQxIwFq1yloBPvjAWgUAoqKga1d7lCjhbn0iIr5CYUACSmKi3fxHjYLvv/eev/566wpo0wb+f98jERH5fwoDEhB277ZugPHjYd8+O5c7N9x3n4WAa691tz4REV+mMCB+y3Hgu++sK+CTT2ytAIArroDu3aFzZ4iMdLdGERF/oDAgfufECXj/fQsBa9d6z990k7UCtGoFuXK5VZ2IiP9RGBC/sWOHLQ709ttw8KCdy5MH2rWzWQF16rhanoiI31IYEJ/mOLBkibUCzJkDqal2vkwZ2zK4UycoUsTVEkVE/J7CgPikY8dso6BRo2DdOu/5xo2tK+Df/7a1AkRE5PLpcio+5c8/bcvgCRPgyBE7FxYGDz5oXQHVq7tanohIQFIYENc5DixYYF0Bn31mHwNUrGhdAR07QsGCrpYoIhLQFAbENUePwpQp1hWwcaP3fPPm1hVw222QI4d79YmIBAuFAcl2f/xhAWDyZAsEAOHhEB1tLQFVqrhZnYhI8FEYkGyRmgrz5llXwLx53vNVqthYgA4dLBCIiEj2UxiQLBUXB5MmwejRsHmznfN4oGVL6wpo2lRdASIibtNlWLLEb7/ZksBXXAG9e1sQiIiAPn1g0yaYOxeaNVMQuFzlypXD4/Gke7z22mtulyUifkYtA5JpUlJsNsCIEbBokfd89erWCtC+PeTL5159geqll16ic+fOpz8OV3+LiFwkhQG5bIcO2boAb70F27bZuRw54M47LQQ0amRdA5I1wsPDiYqKcrsMEfFjHsdJm9UtcnF++cUGBE6bZpsHARQubLsFdu8OZcu6W18wKFeuHCdPniQ5OZkyZcpw//3307t3b3KeZ3nGxMREEhMTT38cHx9P6dKliYuLo0CBAtlRtoj4GLUMyEU5dQpmz7augG++8Z6vU8daAdq2hbx53aou+Dz66KNcc801FC5cmOXLl9OvXz/27NnD8OHDz/k1gwYN4sUXX8zGKkXE16llQDJk/37bLXDMGNi1y86FhMDdd1sIaNhQXQGZ5emnn2bw4MHnfc7vv/9O1apVzzg/ceJEunbtSkJCAqGhoWf9WrUMiMg/KQzIea1ZY10BM2ZA2v2jWDHo2hW6dbPZApK59u/fz8G0PZrPoUKFCuTOnfuM8+vXr6dGjRps2LCBKhlcvSk+Pp6IiAiFAZEgpm4COUNSEnz8sYWAFSu856+7zloB/vtfOMcfnZIJihUrRrFixS7pa9euXUuOHDmIjIzM5KpEJJApDMhpsbEwbhyMHWvHALly2c2/Z0+oX9/d+iS9FStWsGrVKm699VbCw8NZsWIFvXv3pn379hQqVMjt8kTEjygMBDnHgVWrrBXgww8hOdnOlyhh3QBduoBmrfmm0NBQZsyYwYABA0hMTKR8+fL07t2bPn36uF2aiPgZjRkIUomJMHOmhYAffvCev+EG2yvg7rvhLF3SEoA0ZkBE1DIQZHbtsm6A8eNthgBY/3/bttYVcM017tYnIiLZT2EgCDgOfPuttQJ88oktGwxQqhQ88gg8/LDNEBARkeCkMBDATpyA6dMtBPz8s/f8LbdYK8Cdd8J5FqoTEZEgoVtBANq+3fYJeOcd2zcAbFXA9u1tPECtWu7WJyIivkVhIEA4DixebK0An34Kqal2vlw5iImBhx6yfQNERPxBSkoKyWnTm+SccuXKRUhIyGW/jsKAn0tIgKlTYdQoWL/ee75pU+sKaNnSlg0WEfEHjuMQGxvLkSNH3C7FbxQsWJCoqCg8l7EmvMKAn9qyBUaPhokTIS7OzuXLBx06WFdAtWru1icicinSgkBkZCRhYWGXdYMLdI7jcPz4cfbt2wdAiRIlLvm1FAb8SGoqzJ9vXQFffGFdAwCVKlkAiI6GiAhXSxQRuWQpKSmng0CRIkXcLscv5P3/bWL37dtHZGTkJXcZKAz4gfh4mDLFugL++MN7/rbbrCugeXPIkcO9+kREMkPaGIGwsDCXK/Evae9XcnKywkAg2rjRAsDkyTY2ACA8HDp2tEGBV17pankiIllCXQMXJzPeL4UBH5OSAl9+aV0BX3/tPV+tmnUFPPCABQIREZHMojDgI44cscGAo0fDn3/aOY8H7rjDugKaNLGPRUREMpt6ml22bp3tDnjFFdC3rwWBggXh8cdtxsCcOTZNUEFARMR3OY5Dly5dKFy4MB6Ph7Vr13Lw4EEiIyPZtm1bhl4jKSmJcuXK8cPfd4/LJmoZcMGpUzB3rnUFLF7sPV+jhrUCtGtn0wRFRMQ/zJs3j8mTJ7NkyRIqVKhA0aJFefLJJ7nzzjspV65chl4jd+7cPP744zz11FMsXLgwawv+B4WBbHTwoC0R/NZbsGOHncuRA1q3thBw881qARAR8UdbtmyhRIkS3HDDDQAcP36cCRMm8NVXX13U67Rr146+ffuyfv16qlevnhWlnpXCQDZYu9ZaAaZPh5Mn7VyRItClC3TvDqVLu1qeiIjvchw4fjz7v29YWIb/OouOjmbKlCmAjewvW7Ysw4YNIzQ0lOuvv/7081566SXGjh3Lr7/+enodhZYtW3L8+HEWLlxIjhw5KFSoEA0bNmTGjBkMHDgw8/9d56AwkEWSk2HWLAsB337rPX/11fDoo3DffZAnj3v1iYj4hePHIX/+7P++CQkZ7q998803qVixIuPHj2f16tWEhITw8ssvU7du3XTPe/bZZ5k3bx4PP/wws2bNYvTo0Sxfvpyff/6ZHH9bLKZevXp88803mfrPuRCFgUy2bx+MHw9jx8Lu3XYuZ05o08a6Aho0UFeAiEggiYiIIDw8nJCQEKKiogDYvn07JUuWTPe8kJAQpk6dSp06dXj66acZMWIE77zzDmXKlEn3vJIlS7J9+/Zsqx8UBjLNDz/AiBEwcyYkJdm5yEjo2tVmC/zjZ0JERDIiLMy76lp2f9/LcOLECfKcpfm3QoUKDBs2jK5du3Lvvfdy//33n/GcvHnzcjybu0YUBi5DUhJ8+KF1Baxa5T1fr561AtxzD4SGulefiIjf83j8cnpV0aJFOXz48Fk/t2zZMkJCQti2bRunTp0iZ870t+JDhw5RrFix7CjzNK0zcAn27IH+/aFMGWjf3oJA7tze41Wr7FhBQEQkOF199dX89ttvZ5yfOXMmn3zyCUuWLGHHjh1nHSS4bt06rr766uwo8zSFgQxyHFi+HNq2tRDw0kuwd681/w8caFMF33vPWgVERCS4NW/enPXr16drHdi1axfdu3dn8ODB3HjjjUyaNIlXX32VlStXpvvab775hmbNmmVrvQoDF3DypG0UdO210LAhzJhhiwY1bGjjA7Ztg+eeg+LF3a5URER8Rc2aNbnmmmv44IMPAFuhMDo6mnr16tGjRw/AAkP37t1p3749Cf8/LmLFihXExcXRpk2bbK3X4ziOk63f0U/s3AljxsDbb8OBA3YuNBTuv9/GA2RzC45IlomPjyciIoK4uDgKFCjgdjkSxE6ePMnWrVspX778WQff+ZvPP/+cJ554gnXr1qWbOng+9957L7Vr1+aZZ57J8PfJjPdNAwj/xnFg2TIbEDh7tu0gCLYo0COPwMMPQ9GirpYoIiJ+omXLlmzatIndu3dTOgOryyUlJVGzZk169+6dDdWlp5YBbE2LadMsBPz6q/d8o0bWCvCf/9haASKBSC0D4isCrWUgu6hl4DJt3Wr7BEyYAGljPPLmhQcegB49oGZNd+sTERHJDkEXBhwHFiyAUaNs58C0dpHy5SEmBh56CAoVcrdGERGR7BQ0YSAhAd5910LA7797z//rX9YVcPvtEBLiXn0iIiJuCfgwsGkTjB4NkyZBfLydy58fOnSwroCqVd2tT0RExG0BGQZSU+Hrr22vgC+/9J6vXNkCQHQ0aJyUiIiICagwEBdnCwSNHm0tAmluv926Apo1gwxO9RQREQkaAREGfv/dxgK8+653c6sCBWwwYEwMVKrkbn0iIiK+zG//Tk5JgU8/tQGAV11lUwQTEqBaNTvevRtef11BQEREsp7jOHTp0oXChQvj8XhYu3btGc/ZuHEjUVFRHD16NEOveeDAASIjI9m1a1cmV3smvwsDhw7BsGF2k7/zTpsmmCOH93j9euje3QYJioiIZId58+YxefJkPvvsM/bs2cOwYcMYMGBAuuf069ePnj17Eh4enqHXLFq0KA8++CD9+/fPgorT85tugl9/tRUCp06FEyfsXKFCtkTwI49AuXKuliciIkFsy5YtlChRghtuuAGAnP9YtnbHjh189tlnjBw58qJet2PHjtStW5ehQ4dSuHDhTKv3n3w6DJw6BXPmWAhYutR7vlYtGxB4//0QFuZefSIikrUcx5aMz25hYeDxZOy50dHRTJkyBQCPx0PZsmVp1KhRuud88MEH1K5dmyuuuOL0uYceeogffviB1atXExoaSlJSEvXr16dmzZq8++67AFSvXp2SJUsya9YsOnXqlCn/trPxyW6CAwdg0CCoUAHatLEgEBLiPV671loEFAQkkL3yyivccMMNhIWFUbBgwbM+Z8eOHbRs2ZKwsDAiIyN54oknOHXqVPYWKpKFjh+3bt/sflxMAHnzzTd56aWXKFWqFHv27GH16tVnPOebb77h2muvTXduxIgRHDt2jKeffhqAZ599liNHjjBq1Kh0z6tXrx7ffPPNxb95F8GnWgZ+/NFaAd5/HxIT7VzRotClC3TrZrsHigSLpKQk7rnnHho0aMCECRPO+HxKSgotW7YkKiqK5cuXs2fPHh588EFy5crFq6++6kLFIsEpIiKC8PBwQkJCiIqKAmDy5MnpnrN9+/YzwkD+/PmZOnUqt9xyC+Hh4bzxxhssXrz4jA3DSpYsyU8//ZSl/wbXw0ByMnz8sYWA5cu95+vWta6Ae+8FbV4lwejFF18EzryopPn666/57bffWLBgAcWLF6dOnToMHDiQp556igEDBpA7d+5srFYka4SFeaeMZ/f3zUwnTpw4646CDRo04PHHHz/9u3vjjTee8Zy8efNyPIv7SlwNA2PGwMCBsGfP/xeTE+65x0LA9ddnvL9GJBitWLGCmjVrUrx48dPnmjdvTvfu3Vm/fj1XX331Wb8uMTGRxLSmN2wLYxFf5fFAvnxuV3H5ihYtyuG07XH/JjU1le+++46QkBA2b9581q89dOgQxYoVy9L6XB0zcPSoBYGoKOjfH3bsgOnToUEDBQGRC4mNjU0XBIDTH8fGxp7z6wYNGkRERMTpR2n1v4lkuauvvprffvvtjPNDhw5lw4YNLF26lHnz5jFp0qQznrNu3bpzhvvM4moY6NQJpk2D7dthwAAoUcLNakSy3tNPP43H4znvY8OGDVlaQ79+/YiLizv92LlzZ5Z+PxGxVrsVK1aQkpJy+txPP/3ECy+8wDvvvEPDhg0ZPnw4jz32GH/++efp5xw/fpw1a9bQrFmzLK3P1W6CIkVseqBIsOjbty/R0dHnfU6FChUy9FpRUVF8//336c7t3bv39OfOJTQ0lNDQ0Ax9DxHJHLfddhs5c+ZkwYIFNG/enJMnT9K+fXuio6O54447AOjSpQuff/45DzzwAMuWLSMkJIQ5c+ZQpkwZbrrppiytz/UBhCLBpFixYpnW99egQQNeeeUV9u3bR2RkJADz58+nQIECXHXVVZnyPUQkY3r16kWvXr3O+fmcOXPyzDPPMHz4cJo3b06ePHlYv379Gc+bM2dOuo/ffPNNXnjhhcwu98z6svw7iMgl2bFjB4cOHWLHjh2kpKScXuu8UqVK5M+fn2bNmnHVVVfxwAMPMGTIEGJjY3nuueeIiYnRX/4iPqhr164cOXKEo0ePZmhJ4gMHDnDXXXfRtm3bLK/N4ziOk+XfRUQu2t9XNfu7xYsXn17dbPv27XTv3p0lS5aQL18+OnTowGuvvXbGUqjnEx8fT0REBHFxcWfMbxbJTidPnmTr1q2UL1/+rNPw5Owy431TGBAJcgoD4isUBi5NZrxvPrkcsYiIBC/9jXpxMuP9UhgQERGfkCtXLoAsX20v0KS9X2nv36XQAEIREfEJISEhFCxYkH379gEQFhaGRyvQnZPjOBw/fpx9+/ZRsGBBQkJCLvm1FAZERMRnpK2RkRYI5MIKFix43rVFMkJhQEREfIbH46FEiRJERkaSnJzsdjk+L1euXJfVIpBGYUBERHxOSEhIptzkJGM0gFBERCTIKQyIiIgEOYUBERGRIKcwICIiEuQUBkRERIKc9iYQCXKO45zeRU0LvIgEJ4UBERGRIKduAhERkSCnMCAiIhLkFAZERESCnMKAiIhIkFMYEBERCXIKAyIiIkFOYUBERCTI/R8NnCOO7lHTiwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"# Output of plotting commands is displayed inline within the Jupyter notebook.\n",
"%matplotlib inline\n",
"\n",
"def plot_f1_and_f2(f1, f2=None, x_min=-5, x_max=5, label1=\"f(x)\", label2=\"f'(x)\"):\n",
" x = np.linspace(x_min, x_max,100)\n",
"\n",
" # Setting the axes at the centre.\n",
" fig = plt.figure()\n",
" ax = fig.add_subplot(1, 1, 1)\n",
" ax.spines['left'].set_position('center')\n",
" ax.spines['bottom'].set_position('zero')\n",
" ax.spines['right'].set_color('none')\n",
" ax.spines['top'].set_color('none')\n",
" ax.xaxis.set_ticks_position('bottom')\n",
" ax.yaxis.set_ticks_position('left')\n",
"\n",
" plt.plot(x, f1(x), 'r', label=label1)\n",
" if not f2 is None:\n",
" # If f2 is an array, it is passed as it is to be plotted as unlinked points.\n",
" # If f2 is a function, f2(x) needs to be passed to plot it. \n",
" if isinstance(f2, np.ndarray):\n",
" plt.plot(x, f2, 'bo', markersize=3, label=label2,)\n",
" else:\n",
" plt.plot(x, f2(x), 'b', label=label2)\n",
" plt.legend()\n",
"\n",
" plt.show()\n",
" \n",
"plot_f1_and_f2(f, dfdx)"
]
},
{
"cell_type": "markdown",
"id": "ffb711b8",
"metadata": {},
"source": [
"In real life the functions are more complicated and it is not possible to calculate the derivatives analytically every time. Let's explore which tools and libraries are available in Python for the computation of derivatives without manual derivation."
]
},
{
"cell_type": "markdown",
"id": "8d17f76f",
"metadata": {},
"source": [
"<a name='2'></a>\n",
"## 2 - Symbolic Differentiation"
]
},
{
"cell_type": "markdown",
"id": "severe-studio",
"metadata": {},
"source": [
"**Symbolic computation** deals with the computation of mathematical objects that are represented exactly, not approximately (e.g. $\\sqrt{2}$ will be written as it is, not as $1.41421356237$). For differentiation it would mean that the output will be somehow similar to if you were computing derivatives by hand using rules (analytically). Thus, symbolic differentiation can produce exact derivatives."
]
},
{
"cell_type": "markdown",
"id": "5b51ca07",
"metadata": {},
"source": [
"<a name='2.1'></a>\n",
"### 2.1 - Introduction to Symbolic Computation with `SymPy`\n",
"\n",
"Let's explore symbolic differentiation in Python with commonly used `SymPy` library.\n",
"\n",
"If you want to compute the approximate decimal value of $\\sqrt{18}$, you could normally do it in the following way:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "52d0b0a6",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4.242640687119285"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import math\n",
"\n",
"math.sqrt(18)"
]
},
{
"cell_type": "markdown",
"id": "a5227c24",
"metadata": {},
"source": [
"The output $4.242640687119285$ is an approximate result. You may recall that $\\sqrt{18} = \\sqrt{9 \\cdot 2} = 3\\sqrt{2}$ and see that it is pretty much impossible to deduct it from the approximate result. But with the symbolic computation systems the roots are not approximated with a decimal number but rather only simplified, so the output is exact:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "d5f1747d",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 3 \\sqrt{2}$"
],
"text/plain": [
"3*sqrt(2)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# This format of module import allows to use the sympy functions without sympy. prefix.\n",
"from sympy import *\n",
"\n",
"# This is actually sympy.sqrt function, but sympy. prefix is omitted.\n",
"sqrt(18)"
]
},
{
"cell_type": "markdown",
"id": "f5495348",
"metadata": {},
"source": [
"Numerical evaluation of the result is available, and you can set number of the digits to show in the approximated output:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "925375b9",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 4.2426407$"
],
"text/plain": [
"4.2426407"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N(sqrt(18),8)"
]
},
{
"cell_type": "markdown",
"id": "fb625f44",
"metadata": {},
"source": [
"In `SymPy` variables are defined using **symbols**. In this particular library they need to be predefined (a list of them should be provided). Have a look in the cell below, how the symbolic expression, correspoinding to the mathematical expression $2x^2 - xy$, is defined:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "b52721c0",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 2 x^{2} - x y$"
],
"text/plain": [
"2*x**2 - x*y"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# List of symbols.\n",
"x, y = symbols('x y')\n",
"# Definition of the expression.\n",
"expr = 2 * x**2 - x * y\n",
"expr"
]
},
{
"cell_type": "markdown",
"id": "edffde06",
"metadata": {},
"source": [
"Now you can perform various manipulations with this expression: add or subtract some terms, multiply by other expressions etc., just like if you were doing it by hands:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "ac575c24",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle x \\left(x^{3} + 2 x^{2}\\right)$"
],
"text/plain": [
"x*(x**3 + 2*x**2)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"expr_manip = x * (expr + x * y + x**3)\n",
"expr_manip"
]
},
{
"cell_type": "markdown",
"id": "db151265",
"metadata": {},
"source": [
"You can also expand the expression:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "afb04770",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle x^{4} + 2 x^{3}$"
],
"text/plain": [
"x**4 + 2*x**3"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"expand(expr_manip)"
]
},
{
"cell_type": "markdown",
"id": "77b3750c",
"metadata": {},
"source": [
"Or factorise it:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "0fc456cb",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle x^{3} \\left(x + 2\\right)$"
],
"text/plain": [
"x**3*(x + 2)"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"factor(expr_manip)"
]
},
{
"cell_type": "markdown",
"id": "c885c385",
"metadata": {},
"source": [
"To substitute particular values for the variables in the expression, you can use the following code:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "3c7d239e",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 4.0$"
],
"text/plain": [
"4.00000000000000"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"expr.evalf(subs={x:-1, y:2})"
]
},
{
"cell_type": "markdown",
"id": "4acf4536",
"metadata": {},
"source": [
"This can be used to evaluate a function $f\\left(x\\right) = x^2$:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "622e4335",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 9.0$"
],
"text/plain": [
"9.00000000000000"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f_symb = x ** 2\n",
"f_symb.evalf(subs={x:3})"
]
},
{
"cell_type": "markdown",
"id": "1346fb9f",
"metadata": {},
"source": [
"You might be wondering now, is it possible to evaluate the symbolic functions for each element of the array? At the beginning of the lab you have defined a `NumPy` array `x_array`:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "af6562d3",
"metadata": {},
"outputs": [],
"source": [
"print(x_array)"
]
},
{
"cell_type": "markdown",
"id": "f13ae4e0",
"metadata": {},
"source": [
"Now try to evaluate function `f_symb` for each element of the array. You will get an error:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "11df38ae",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"'Pow' object is not callable\n"
]
}
],
"source": [
"try:\n",
" f_symb(x_array)\n",
"except TypeError as err:\n",
" print(err)"
]
},
{
"cell_type": "markdown",
"id": "4f7800fe",
"metadata": {},
"source": [
"It is possible to evaluate the symbolic functions for each element of the array, but you need to make a function `NumPy`-friendly first:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "3d7d8a58",
"metadata": {},
"outputs": [],
"source": [
"from sympy.utilities.lambdify import lambdify\n",
"\n",
"f_symb_numpy = lambdify(x, f_symb, 'numpy')"
]
},
{
"cell_type": "markdown",
"id": "631e3e71",
"metadata": {},
"source": [
"The following code should work now:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "1ab8f0da",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x: \n",
" [1 2 3]\n",
"f(x) = x**2: \n",
" [1 4 9]\n"
]
}
],
"source": [
"print(\"x: \\n\", x_array)\n",
"print(\"f(x) = x**2: \\n\", f_symb_numpy(x_array))"
]
},
{
"cell_type": "markdown",
"id": "681062e1",
"metadata": {},
"source": [
"`SymPy` has lots of great functions to manipulate expressions and perform various operations from calculus. More information about them can be found in the official documentation [here](https://docs.sympy.org/)."
]
},
{
"cell_type": "markdown",
"id": "6c6c9d29",
"metadata": {},
"source": [
"<a name='2.2'></a>\n",
"### 2.2 - Symbolic Differentiation with `SymPy`\n",
"\n",
"Let's try to find a derivative of a simple power function using `SymPy`:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "abb93aeb",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 3 x^{2}$"
],
"text/plain": [
"3*x**2"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"diff(x**3,x)"
]
},
{
"cell_type": "markdown",
"id": "bd8a81b1",
"metadata": {},
"source": [
"Some standard functions can be used in the expression, and `SymPy` will apply required rules (sum, product, chain) to calculate the derivative:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "63ce2cd2",
"metadata": {},
"outputs": [],
"source": [
"dfdx_composed = diff(exp(-2*x) + 3*sin(3*x), x)\n",
"dfdx_composed"
]
},
{
"cell_type": "markdown",
"id": "13407c2b",
"metadata": {},
"source": [
"Now calculate the derivative of the function `f_symb` defined in [2.1](#2.1) and make it `NumPy`-friendly:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "597d5879",
"metadata": {},
"outputs": [],
"source": [
"dfdx_symb = diff(f_symb, x)\n",
"dfdx_symb_numpy = lambdify(x, dfdx_symb, 'numpy')"
]
},
{
"cell_type": "markdown",
"id": "0f9bee28",
"metadata": {},
"source": [
"Evaluate function `dfdx_symb_numpy` for each element of the `x_array`:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b74b4d04",
"metadata": {},
"outputs": [],
"source": [
"print(\"x: \\n\", x_array)\n",
"print(\"f'(x) = 2x: \\n\", dfdx_symb_numpy(x_array))"
]
},
{
"cell_type": "markdown",
"id": "ada41a99",
"metadata": {},
"source": [
"You can apply symbolically defined functions to the arrays of larger size. The following code will plot function and its derivative, you can see that it works:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "031a757c",
"metadata": {},
"outputs": [],
"source": [
"plot_f1_and_f2(f_symb_numpy, dfdx_symb_numpy)"
]
},
{
"cell_type": "markdown",
"id": "de01feee",
"metadata": {},
"source": [
"<a name='2.3'></a>\n",
"### 2.3 - Limitations of Symbolic Differentiation\n",
"\n",
"Symbolic Differentiation seems to be a great tool. But it also has some limitations. Sometimes the output expressions are too complicated and even not possible to evaluate. For example, find the derivative of the function $$\\left|x\\right| = \\begin{cases} x, \\ \\text{if}\\ x > 0\\\\ -x, \\ \\text{if}\\ x < 0 \\\\ 0, \\ \\text{if}\\ x = 0\\end{cases}$$ Analytically, its derivative is:\n",
"$$\\frac{d}{dx}\\left(\\left|x\\right|\\right) = \\begin{cases} 1, \\ \\text{if}\\ x > 0\\\\ -1, \\ \\text{if}\\ x < 0\\\\\\ \\text{does not exist}, \\ \\text{if}\\ x = 0\\end{cases}$$\n",
"\n",
"Have a look the output from the symbolic differentiation:"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "collect-needle",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{\\left(\\operatorname{re}{\\left(x\\right)} \\frac{d}{d x} \\operatorname{re}{\\left(x\\right)} + \\operatorname{im}{\\left(x\\right)} \\frac{d}{d x} \\operatorname{im}{\\left(x\\right)}\\right) \\operatorname{sign}{\\left(x \\right)}}{x}$"
],
"text/plain": [
"(re(x)*Derivative(re(x), x) + im(x)*Derivative(im(x), x))*sign(x)/x"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dfdx_abs = diff(abs(x),x)\n",
"dfdx_abs"
]
},
{
"cell_type": "markdown",
"id": "f9c3d8e5",
"metadata": {},
"source": [
"Looks complicated, but it would not be a problem if it was possible to evaluate. But check, that for $x=-2$ instead of the derivative value $-1$ it outputs some unevaluated expression:"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "d53e7c64",
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\frac{\\left(\\operatorname{re}{\\left(x\\right)} \\frac{d}{d x} \\operatorname{re}{\\left(x\\right)} + \\operatorname{im}{\\left(x\\right)} \\frac{d}{d x} \\operatorname{im}{\\left(x\\right)}\\right) \\operatorname{sign}{\\left(x \\right)}}{x}$"
],
"text/plain": [
"(re(x)*Derivative(re(x), x) + im(x)*Derivative(im(x), x))*sign(x)/x"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dfdx_abs.evalf(subs={x:-2})"
]
},
{
"cell_type": "markdown",
"id": "f4a9140c",
"metadata": {},
"source": [
"And in the `NumPy` friendly version it also will give an error:"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "644f00ce",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"name 'Derivative' is not defined\n"
]
}
],
"source": [
"dfdx_abs_numpy = lambdify(x, dfdx_abs,'numpy')\n",
"\n",
"try:\n",
" dfdx_abs_numpy(np.array([1, -2, 0]))\n",
"except NameError as err:\n",
" print(err)"
]
},
{
"cell_type": "markdown",
"id": "3e1f3c94",
"metadata": {},
"source": [
"In fact, there are problems with the evaluation of the symbolic expressions wherever there is a \"jump\" in the derivative (e.g. function expressions are different for different intervals of $x$), like it happens with $\\frac{d}{dx}\\left(\\left|x\\right|\\right)$. \n",
"\n",
"Also, you can see in this example, that you can get a very complicated function as an output of symbolic computation. This is called **expression swell**, which results in unefficiently slow computations. You will see the example of that below after learning other differentiation libraries in Python."
]
},
{
"cell_type": "markdown",
"id": "1cb85963",
"metadata": {},
"source": [
"<a name='3'></a>\n",
"## 3 - Numerical Differentiation\n",
"\n",
"This method does not take into account the function expression. The only important thing is that the function can be evaluated in the nearby points $x$ and $x+\\Delta x$, where $\\Delta x$ is sufficiently small. Then $\\frac{df}{dx}\\approx\\frac{f\\left(x + \\Delta x\\right) - f\\left(x\\right)}{\\Delta x}$, which can be called a **numerical approximation** of the derivative. \n",
"\n",
"Based on that idea there are different approaches for the numerical approximations, which somehow vary in the computation speed and accuracy. However, for all of the methods the results are not accurate - there is a round off error. At this stage there is no need to go into details of various methods, it is enough to investigate one of the numerial differentiation functions, available in `NumPy` package."
]
},
{
"cell_type": "markdown",
"id": "1cc2a87e",
"metadata": {},
"source": [
"<a name='3.1'></a>\n",
"### 3.1 - Numerical Differentiation with `NumPy`"
]
},
{
"cell_type": "markdown",
"id": "c469b76c",
"metadata": {},
"source": [
"You can call function `np.gradient` to find the derivative of function $f\\left(x\\right) = x^2$ defined above. The first argument is an array of function values, the second defines the spacing $\\Delta x$ for the evaluation. Here pass it as an array of $x$ values, the differences will be calculated automatically. You can find the documentation [here](https://numpy.org/doc/stable/reference/generated/numpy.gradient.html)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b275519f",
"metadata": {},
"outputs": [],
"source": [
"x_array_2 = np.linspace(-5, 5, 100)\n",
"dfdx_numerical = np.gradient(f(x_array_2), x_array_2)\n",
"\n",
"plot_f1_and_f2(dfdx_symb_numpy, dfdx_numerical, label1=\"f'(x) exact\", label2=\"f'(x) approximate\")"
]
},
{
"cell_type": "markdown",
"id": "6a1d5843",
"metadata": {},
"source": [
"Try to do numerical differentiation for more complicated function:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9fa0d7cf",
"metadata": {},
"outputs": [],
"source": [
"def f_composed(x):\n",
" return np.exp(-2*x) + 3*np.sin(3*x)\n",
"\n",
"plot_f1_and_f2(lambdify(x, dfdx_composed, 'numpy'), np.gradient(f_composed(x_array_2), x_array_2),\n",
" label1=\"f'(x) exact\", label2=\"f'(x) approximate\")"
]
},
{
"cell_type": "markdown",
"id": "826da796",
"metadata": {},
"source": [
"The results are pretty impressive, keeping in mind that it does not matter at all how the function was calculated - only the final values of it!"
]
},
{
"cell_type": "markdown",
"id": "bc60825b",
"metadata": {},
"source": [
"<a name='3.2'></a>\n",
"### 3.2 - Limitations of Numerical Differentiation"
]
},
{
"cell_type": "markdown",
"id": "8dbf76a0",
"metadata": {},
"source": [
"Obviously, the first downside of the numerical differentiation is that it is not exact. However, the accuracy of it is normally enough for machine learning applications. At this stage there is no need to evaluate errors of the numerical differentiation.\n",
"\n",
"Another problem is similar to the one which appeared in the symbolic differentiation: it is inaccurate at the points where there are \"jumps\" of the derivative. Let's compare the exact derivative of the absolute value function and with numerical approximation:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "28bb6a5f",
"metadata": {},
"outputs": [],
"source": [
"def dfdx_abs(x):\n",
" if x > 0:\n",
" return 1\n",
" else:\n",
" if x < 0:\n",
" return -1\n",
" else:\n",
" return None\n",
"\n",
"plot_f1_and_f2(np.vectorize(dfdx_abs), np.gradient(abs(x_array_2), x_array_2))"
]
},
{
"cell_type": "markdown",
"id": "229fdab5",
"metadata": {},
"source": [
"\n"
]
},
{
"cell_type": "markdown",
"id": "00ccb653",
"metadata": {},
"source": [
"You can see that the results near the \"jump\" are $0.5$ and $-0.5$, while they should be $1$ and $-1$. These cases can give significant errors in the computations.\n",
"\n",
"But the biggest problem with the numerical differentiation is slow speed. It requires function evalutation every time. In machine learning models there are hundreds of parameters and there are hundreds of derivatives to be calculated, performing full function evaluation every time slows down the computation process. You will see the example of it below."
]
},
{
"cell_type": "markdown",
"id": "2caeb33f",
"metadata": {},
"source": [
"<a name='4'></a>\n",
"## 4 - Automatic Differentiation"
]
},
{
"cell_type": "markdown",
"id": "eba8f444",
"metadata": {},
"source": [
"**Automatic differentiation** (autodiff) method breaks down the function into common functions ($sin$, $cos$, $log$, power functions, etc.), and constructs the computational graph consisting of the basic functions. Then the chain rule is used to compute the derivative at any node of the graph. It is the most commonly used approach in machine learning applications and neural networks, as the computational graph for the function and its derivatives can be built during the construction of the neural network, saving in future computations.\n",
"\n",
"The main disadvantage of it is implementational difficulty. However, nowadays there are libraries that are convenient to use, such as [MyGrad](https://mygrad.readthedocs.io/en/latest/index.html), [Autograd](https://autograd.readthedocs.io/en/latest/) and [JAX](https://jax.readthedocs.io/en/latest/). `Autograd` and `JAX` are the most commonly used in the frameworks to build neural networks. `JAX` brings together `Autograd` functionality for optimization problems, and `XLA` (Accelerated Linear Algebra) compiler for parallel computing.\n",
"\n",
"The syntax of `Autograd` and `JAX` are slightly different. It would be overwhelming to cover both at this stage. In this notebook you will be performing automatic differentiation using one of them: `JAX`."
]
},
{
"cell_type": "markdown",
"id": "20071067",
"metadata": {},
"source": [
"<a name='4.1'></a>\n",
"### 4.1 - Introduction to `JAX`"
]
},
{
"cell_type": "markdown",
"id": "f444d827",
"metadata": {},
"source": [
"To begin with, load the required libraries. From `jax` package you need to load just a couple of functions for now (`grad` and `vmap`). Package `jax.numpy` is a wrapped `NumPy`, which pretty much replaces `NumPy` when `JAX` is used. It can be loaded as `np` as if it was an original `NumPy` in most of the cases. However, in this notebook you'll upload it as `jnp` to distinguish them for now."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "85d818ea",
"metadata": {},
"outputs": [],
"source": [
"from jax import grad, vmap\n",
"import jax.numpy as jnp"
]
},
{
"cell_type": "markdown",
"id": "d42ede8e",
"metadata": {},
"source": [
"Create a new `jnp` array and check its type."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8856647c",
"metadata": {},
"outputs": [],
"source": [
"x_array_jnp = jnp.array([1., 2., 3.])\n",
"\n",
"print(\"Type of NumPy array:\", type(x_array))\n",
"print(\"Type of JAX NumPy array:\", type(x_array_jnp))\n",
"# Please ignore the warning message if it appears."
]
},
{
"cell_type": "markdown",
"id": "730a2dd3",
"metadata": {},
"source": [
"The same array can be created just converting previously defined `x_array = np.array([1, 2, 3])`, although in some cases `JAX` does not operate with integers, thus the values need to be converted to floats. You will see an example of it below."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3008671b",
"metadata": {},
"outputs": [],
"source": [
"x_array_jnp = jnp.array(x_array.astype('float32'))\n",
"print(\"JAX NumPy array:\", x_array_jnp)\n",
"print(\"Type of JAX NumPy array:\", type(x_array_jnp))"
]
},
{
"cell_type": "markdown",
"id": "f81ce077",
"metadata": {},
"source": [
"Note, that `jnp` array has a specific type `jaxlib.xla_extension.DeviceArray`. In most of the cases the same operators and functions are applicable to them as in the original `NumPy`, for example:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "742003ec",
"metadata": {},
"outputs": [],
"source": [
"print(x_array_jnp * 2)\n",
"print(x_array_jnp[2])"
]
},
{
"cell_type": "markdown",
"id": "3c7ef8a4",
"metadata": {},
"source": [
"But sometimes working with `jnp` arrays the approach needs to be changed. In the following code, trying to assign a new value to one of the elements, you will get an error:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3fc00cab",
"metadata": {},
"outputs": [],
"source": [
"try:\n",
" x_array_jnp[2] = 4.0\n",
"except TypeError as err:\n",
" print(err)"
]
},
{
"cell_type": "markdown",
"id": "cf9e29fe",
"metadata": {},
"source": [
"To assign a new value to an element in the `jnp` array you need to apply functions `.at[i]`, stating which element to update, and `.set(value)` to set a new value. These functions also operate **out-of-place**, the updated array is returned as a new array and the original array is not modified by the update."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ffc53ad2",
"metadata": {},
"outputs": [],
"source": [
"y_array_jnp = x_array_jnp.at[2].set(4.0)\n",
"print(y_array_jnp)"
]
},
{
"cell_type": "markdown",
"id": "05a07ce0",
"metadata": {},
"source": [
"Although, some of the `JAX` functions will work with arrays defined with `np` and `jnp`. In the following code you will get the same result in both lines:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5b80429d",
"metadata": {},
"outputs": [],
"source": [
"print(jnp.log(x_array))\n",
"print(jnp.log(x_array_jnp))"
]
},
{
"cell_type": "markdown",
"id": "89397092",
"metadata": {},
"source": [
"This is probably confusing - which `NumPy` to use then? Usually when `JAX` is used, only `jax.numpy` gets imported as `np`, and used instead of the original one."
]
},
{
"cell_type": "markdown",
"id": "20f12b94",
"metadata": {},
"source": [
" <a name='4.2'></a>\n",
"### 4.2 - Automatic Differentiation with `JAX` "
]
},
{
"cell_type": "markdown",
"id": "9cd26792",
"metadata": {},
"source": [
"Time to do automatic differentiation with `JAX`. The following code will calculate the derivative of the previously defined function $f\\left(x\\right) = x^2$ at the point $x = 3$:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "070e417a",
"metadata": {},
"outputs": [],
"source": [
"print(\"Function value at x = 3:\", f(3.0))\n",
"print(\"Derivative value at x = 3:\",grad(f)(3.0))"
]
},
{
"cell_type": "markdown",
"id": "3514bda9",
"metadata": {},
"source": [
"Very easy, right? Keep in mind, please, that this cannot be done using integers. The following code will output an error:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a50295a3",
"metadata": {},
"outputs": [],
"source": [
"try:\n",
" grad(f)(3)\n",
"except TypeError as err:\n",
" print(err)"
]
},
{
"cell_type": "markdown",
"id": "872bbbc6",
"metadata": {},
"source": [
"Try to apply the `grad` function to an array, calculating the derivative for each of its elements: "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "caf0e431",
"metadata": {},
"outputs": [],
"source": [
"try:\n",
" grad(f)(x_array_jnp)\n",
"except TypeError as err:\n",
" print(err)"
]
},
{
"cell_type": "markdown",
"id": "9452ebc2",
"metadata": {},
"source": [
"There is some broadcasting issue there. You don't need to get into more details of this at this stage, function `vmap` can be used here to solve the problem.\n",
"\n",
"*Note*: Broadcasting is covered in the Course 1 of this Specialization \"Linear Algebra\". You can also review it in the documentation [here](https://numpy.org/doc/stable/user/basics.broadcasting.html#:~:text=The%20term%20broadcasting%20describes%20how,that%20they%20have%20compatible%20shapes.)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f9b28641",
"metadata": {},
"outputs": [],
"source": [
"dfdx_jax_vmap = vmap(grad(f))(x_array_jnp)\n",
"print(dfdx_jax_vmap)"
]
},
{
"cell_type": "markdown",
"id": "933e382f",
"metadata": {},
"source": [
"Great, now `vmap(grad(f))` can be used to calculate the derivative of function `f` for arrays of larger size and you can plot the output:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "da0a1262",
"metadata": {},
"outputs": [],
"source": [
"plot_f1_and_f2(f, vmap(grad(f)))"
]
},
{
"cell_type": "markdown",
"id": "4162d5e5",
"metadata": {},
"source": [
"In the following code you can comment/uncomment lines to visualize the common derivatives. All of them are found using `JAX` automatic differentiation. The results look pretty good!"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f68b4c0e",
"metadata": {},
"outputs": [],
"source": [
"def g(x):\n",
"# return x**3\n",
"# return 2*x**3 - 3*x**2 + 5\n",
"# return 1/x\n",
"# return jnp.exp(x)\n",
"# return jnp.log(x)\n",
"# return jnp.sin(x)\n",
"# return jnp.cos(x)\n",
" return jnp.abs(x)\n",
"# return jnp.abs(x)+jnp.sin(x)*jnp.cos(x)\n",
"\n",
"plot_f1_and_f2(g, vmap(grad(g)))"
]
},
{
"cell_type": "markdown",
"id": "a58ee858",
"metadata": {},
"source": [
"<a name='5'></a>\n",
"## 5 - Computational Efficiency of Symbolic, Numerical and Automatic Differentiation"
]
},
{
"cell_type": "markdown",
"id": "2211158e",
"metadata": {},
"source": [
"In sections [2.3](#2.3) and [3.2](#3.2) low computational efficiency of symbolic and numerical differentiation was discussed. Now it is time to compare speed of calculations for each of three approaches. Try to find the derivative of the same simple function $f\\left(x\\right) = x^2$ multiple times, evaluating it for an array of a larger size, compare the results and time used:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "36c42dac",
"metadata": {},
"outputs": [],
"source": [
"import timeit, time\n",
"\n",
"x_array_large = np.linspace(-5, 5, 1000000)\n",
"\n",
"tic_symb = time.time()\n",
"res_symb = lambdify(x, diff(f(x),x),'numpy')(x_array_large)\n",
"toc_symb = time.time()\n",
"time_symb = 1000 * (toc_symb - tic_symb) # Time in ms.\n",
"\n",
"tic_numerical = time.time()\n",
"res_numerical = np.gradient(f(x_array_large),x_array_large)\n",
"toc_numerical = time.time()\n",
"time_numerical = 1000 * (toc_numerical - tic_numerical)\n",
"\n",
"tic_jax = time.time()\n",
"res_jax = vmap(grad(f))(jnp.array(x_array_large.astype('float32')))\n",
"toc_jax = time.time()\n",
"time_jax = 1000 * (toc_jax - tic_jax)\n",
"\n",
"print(f\"Results\\nSymbolic Differentiation:\\n{res_symb}\\n\" + \n",
" f\"Numerical Differentiation:\\n{res_numerical}\\n\" + \n",
" f\"Automatic Differentiation:\\n{res_jax}\")\n",
"\n",
"print(f\"\\n\\nTime\\nSymbolic Differentiation:\\n{time_symb} ms\\n\" + \n",
" f\"Numerical Differentiation:\\n{time_numerical} ms\\n\" + \n",
" f\"Automatic Differentiation:\\n{time_jax} ms\")"
]
},
{
"cell_type": "markdown",
"id": "493e5457",
"metadata": {},
"source": [
"The results are pretty much the same, but the time used is different. Numerical approach is obviously inefficient when differentiation needs to be performed many times, which happens a lot training machine learning models. Symbolic and automatic approach seem to be performing similarly for this simple example. But if the function becomes a little bit more complicated, symbolic computation will experiance significant expression swell and the calculations will slow down.\n",
"\n",
"*Note*: Sometimes the execution time results may vary slightly, especially for automatic differentiation. You can run the code above a few time to see different outputs. That does not influence the conclusion that numerical differentiation is slower. `timeit` module can be used more efficiently to evaluate execution time of the codes, but that would unnecessary overcomplicate the codes here.\n",
"\n",
"Try to define some polynomial function, which should not be that hard to differentiate, and compare the computation time for its differentiation symbolically and automatically:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "13047a93",
"metadata": {},
"outputs": [],
"source": [
"def f_polynomial_simple(x):\n",
" return 2*x**3 - 3*x**2 + 5\n",
"\n",
"def f_polynomial(x):\n",
" for i in range(3):\n",
" x = f_polynomial_simple(x)\n",
" return x\n",
"\n",
"tic_polynomial_symb = time.time()\n",
"res_polynomial_symb = lambdify(x, diff(f_polynomial(x),x),'numpy')(x_array_large)\n",
"toc_polynomial_symb = time.time()\n",
"time_polynomial_symb = 1000 * (toc_polynomial_symb - tic_polynomial_symb)\n",
"\n",
"tic_polynomial_jax = time.time()\n",
"res_polynomial_jax = vmap(grad(f_polynomial))(jnp.array(x_array_large.astype('float32')))\n",
"toc_polynomial_jax = time.time()\n",
"time_polynomial_jax = 1000 * (toc_polynomial_jax - tic_polynomial_jax)\n",
"\n",
"print(f\"Results\\nSymbolic Differentiation:\\n{res_polynomial_symb}\\n\" + \n",
" f\"Automatic Differentiation:\\n{res_polynomial_jax}\")\n",
"\n",
"print(f\"\\n\\nTime\\nSymbolic Differentiation:\\n{time_polynomial_symb} ms\\n\" + \n",
" f\"Automatic Differentiation:\\n{time_polynomial_jax} ms\")"
]
},
{
"cell_type": "markdown",
"id": "231c9da0",
"metadata": {},
"source": [
"Again, the results are similar, but automatic differentiation is times faster. \n",
"\n",
"With the increase of function computation graph, the efficiency of automatic differentiation compared to other methods raises, because autodiff method uses chain rule!\n",
"\n",
"Congratulations! Now you are equiped with Python tools to perform differentiation."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "93dfd229",
"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.8.8"
}
},
"nbformat": 4,
"nbformat_minor": 5
}