{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "d22a20a1",
            "cell_type": "markdown",
            "source": "<div class='alert alert-warning'>\n\nSciPy's interactive examples with Jupyterlite are experimental and may not always work as expected. Execution of cells containing imports may result in large downloads (up to 60MB of content for the first import from SciPy). Load times when importing from SciPy may take roughly 10-20 seconds. If you notice any problems, feel free to open an [issue](https://github.com/scipy/scipy/issues/new/choose).\n\n</div>",
            "metadata": {}
        },
        {
            "id": "565f69d8",
            "cell_type": "markdown",
            "source": "Evaluate the derivative of ``np.exp`` at several points ``x``.\n",
            "metadata": {}
        },
        {
            "id": "c4824526",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import numpy as np\nfrom scipy.differentiate import derivative\nf = np.exp\ndf = np.exp  # true derivative\nx = np.linspace(1, 2, 5)\nres = derivative(f, x)\nres.df  # approximation of the derivative",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([2.71828183, 3.49034296, 4.48168907, 5.75460268, 7.3890561 ])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "25c6a544",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res.error  # estimate of the error",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([7.13740178e-12, 9.16600129e-12, 1.17594823e-11, 1.51061386e-11,\n       1.94262384e-11])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "1f7ca6ce",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "abs(res.df - df(x))  # true error",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([2.53130850e-14, 3.55271368e-14, 5.77315973e-14, 5.59552404e-14,\n       6.92779167e-14])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "edf90a9f",
            "cell_type": "markdown",
            "source": "Show the convergence of the approximation as the step size is reduced.\nEach iteration, the step size is reduced by `step_factor`, so for\nsufficiently small initial step, each iteration reduces the error by a\nfactor of ``1/step_factor**order`` until finite precision arithmetic\ninhibits further improvement.\n",
            "metadata": {}
        },
        {
            "id": "a01b50ff",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import matplotlib.pyplot as plt\niter = list(range(1, 12))  # maximum iterations\nhfac = 2  # step size reduction per iteration\nhdir = [-1, 0, 1]  # compare left-, central-, and right- steps\norder = 4  # order of differentiation formula\nx = 1\nref = df(x)\nerrors = []  # true error\nfor i in iter:\n    res = derivative(f, x, maxiter=i, step_factor=hfac,\n                     step_direction=hdir, order=order,\n                     # prevent early termination\n                     tolerances=dict(atol=0, rtol=0))\n    errors.append(abs(res.df - ref))\nerrors = np.array(errors)\nplt.semilogy(iter, errors[:, 0], label='left differences')\nplt.semilogy(iter, errors[:, 1], label='central differences')\nplt.semilogy(iter, errors[:, 2], label='right differences')\nplt.xlabel('iteration')\nplt.ylabel('error')\nplt.legend()\nplt.show()",
            "outputs": []
        },
        {
            "id": "b88ded16",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "(errors[1, 1] / errors[0, 1], 1 / hfac**order)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "(0.06215223140159822, 0.0625)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "4e78a74d",
            "cell_type": "markdown",
            "source": "The implementation is vectorized over `x`, `step_direction`, and `args`.\nThe function is evaluated once before the first iteration to perform input\nvalidation and standardization, and once per iteration thereafter.\n",
            "metadata": {}
        },
        {
            "id": "81d64612",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def f(x, p):\n    f.nit += 1\n    return x**p\nf.nit = 0\ndef df(x, p):\n    return p*x**(p-1)\nx = np.arange(1, 5)\np = np.arange(1, 6).reshape((-1, 1))\nhdir = np.arange(-1, 2).reshape((-1, 1, 1))\nres = derivative(f, x, args=(p,), step_direction=hdir, maxiter=1)\nnp.allclose(res.df, df(x, p))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "a2e092ac",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res.df.shape",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "(3, 5, 4)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "c84f8f9a",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "f.nit",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "2"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "1b560781",
            "cell_type": "markdown",
            "source": "By default, `preserve_shape` is False, and therefore the callable\n`f` may be called with arrays of any broadcastable shapes.\nFor example:\n",
            "metadata": {}
        },
        {
            "id": "369b065d",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "shapes = []\ndef f(x, c):\n   shape = np.broadcast_shapes(x.shape, c.shape)\n   shapes.append(shape)\n   return np.sin(c*x)\n\nc = [1, 5, 10, 20]\nres = derivative(f, 0, args=(c,))\nshapes",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[(4,), (4, 8), (4, 2), (3, 2), (2, 2), (1, 2)]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "b4152a39",
            "cell_type": "markdown",
            "source": "To understand where these shapes are coming from - and to better\nunderstand how `derivative` computes accurate results - note that\nhigher values of ``c`` correspond with higher frequency sinusoids.\nThe higher frequency sinusoids make the function's derivative change\nfaster, so more function evaluations are required to achieve the target\naccuracy:\n",
            "metadata": {}
        },
        {
            "id": "345d216c",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res.nfev",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([11, 13, 15, 17], dtype=int32)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "de5ebcdc",
            "cell_type": "markdown",
            "source": "The initial ``shape``, ``(4,)``, corresponds with evaluating the\nfunction at a single abscissa and all four frequencies; this is used\nfor input validation and to determine the size and dtype of the arrays\nthat store results. The next shape corresponds with evaluating the\nfunction at an initial grid of abscissae and all four frequencies.\nSuccessive calls to the function evaluate the function at two more\nabscissae, increasing the effective order of the approximation by two.\nHowever, in later function evaluations, the function is evaluated at\nfewer frequencies because the corresponding derivative has already\nconverged to the required tolerance. This saves function evaluations to\nimprove performance, but it requires the function to accept arguments of\nany shape.\n\n\"Vector-valued\" functions are unlikely to satisfy this requirement.\nFor example, consider\n",
            "metadata": {}
        },
        {
            "id": "42b8733a",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def f(x):\n   return [x, np.sin(3*x), x+np.sin(10*x), np.sin(20*x)*(x-1)**2]",
            "outputs": []
        },
        {
            "id": "39e88675",
            "cell_type": "markdown",
            "source": "This integrand is not compatible with `derivative` as written; for instance,\nthe shape of the output will not be the same as the shape of ``x``. Such a\nfunction *could* be converted to a compatible form with the introduction of\nadditional parameters, but this would be inconvenient. In such cases,\na simpler solution would be to use `preserve_shape`.\n",
            "metadata": {}
        },
        {
            "id": "491d3060",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "shapes = []\ndef f(x):\n    shapes.append(x.shape)\n    x0, x1, x2, x3 = x\n    return [x0, np.sin(3*x1), x2+np.sin(10*x2), np.sin(20*x3)*(x3-1)**2]\n\nx = np.zeros(4)\nres = derivative(f, x, preserve_shape=True)\nshapes",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "[(4,), (4, 8), (4, 2), (4, 2), (4, 2), (4, 2)]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "f98fdcf0",
            "cell_type": "markdown",
            "source": "Here, the shape of ``x`` is ``(4,)``. With ``preserve_shape=True``, the\nfunction may be called with argument ``x`` of shape ``(4,)`` or ``(4, n)``,\nand this is what we observe.\n",
            "metadata": {}
        }
    ]
}