{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "395b3fe6",
            "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": "16925287",
            "cell_type": "markdown",
            "source": "Import all necessary modules.\n",
            "metadata": {}
        },
        {
            "id": "8d1f74a7",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import numpy as np\nfrom scipy.special import pseudo_huber, huber\nimport matplotlib.pyplot as plt",
            "outputs": []
        },
        {
            "id": "c22b2cd2",
            "cell_type": "markdown",
            "source": "Calculate the function for ``delta=1`` at ``r=2``.\n",
            "metadata": {}
        },
        {
            "id": "f188272b",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "pseudo_huber(1., 2.)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "1.2360679774997898"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "6d279cf4",
            "cell_type": "markdown",
            "source": "Calculate the function at ``r=2`` for different `delta` by providing\na list or NumPy array for `delta`.\n",
            "metadata": {}
        },
        {
            "id": "5aee7524",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "pseudo_huber([1., 2., 4.], 3.)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([2.16227766, 3.21110255, 4.        ])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "f781ff57",
            "cell_type": "markdown",
            "source": "Calculate the function for ``delta=1`` at several points by providing\na list or NumPy array for `r`.\n",
            "metadata": {}
        },
        {
            "id": "4ee30179",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "pseudo_huber(2., np.array([1., 1.5, 3., 4.]))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([0.47213595, 1.        , 3.21110255, 4.94427191])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "787117d6",
            "cell_type": "markdown",
            "source": "The function can be calculated for different `delta` and `r` by\nproviding arrays for both with compatible shapes for broadcasting.\n",
            "metadata": {}
        },
        {
            "id": "7527f784",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "r = np.array([1., 2.5, 8., 10.])\ndeltas = np.array([[1.], [5.], [9.]])\nprint(r.shape, deltas.shape)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "(4,) (3, 1)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "15323bb2",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "pseudo_huber(deltas, r)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[ 0.41421356,  1.6925824 ,  7.06225775,  9.04987562],\n       [ 0.49509757,  2.95084972, 22.16990566, 30.90169944],\n       [ 0.49846624,  3.06693762, 27.37435121, 40.08261642]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "4fc87b83",
            "cell_type": "markdown",
            "source": "Plot the function for different `delta`.\n",
            "metadata": {}
        },
        {
            "id": "4de9bb7d",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "x = np.linspace(-4, 4, 500)\ndeltas = [1, 2, 3]\nlinestyles = [\"dashed\", \"dotted\", \"dashdot\"]\nfig, ax = plt.subplots()\ncombined_plot_parameters = list(zip(deltas, linestyles))\nfor delta, style in combined_plot_parameters:\n    ax.plot(x, pseudo_huber(delta, x), label=rf\"$\\delta={delta}$\",\n            ls=style)\nax.legend(loc=\"upper center\")\nax.set_xlabel(\"$x$\")\nax.set_title(r\"Pseudo-Huber loss function $h_{\\delta}(x)$\")\nax.set_xlim(-4, 4)\nax.set_ylim(0, 8)\nplt.show()",
            "outputs": []
        },
        {
            "id": "de71da8e",
            "cell_type": "markdown",
            "source": "Finally, illustrate the difference between `huber` and `pseudo_huber` by\nplotting them and their gradients with respect to `r`. The plot shows\nthat `pseudo_huber` is continuously differentiable while `huber` is not\nat the points $\\pm\\delta$.\n",
            "metadata": {}
        },
        {
            "id": "4712b066",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def huber_grad(delta, x):\n    grad = np.copy(x)\n    linear_area = np.argwhere(np.abs(x) > delta)\n    grad[linear_area]=delta*np.sign(x[linear_area])\n    return grad\ndef pseudo_huber_grad(delta, x):\n    return x* (1+(x/delta)**2)**(-0.5)\nx=np.linspace(-3, 3, 500)\ndelta = 1.\nfig, ax = plt.subplots(figsize=(7, 7))\nax.plot(x, huber(delta, x), label=\"Huber\", ls=\"dashed\")\nax.plot(x, huber_grad(delta, x), label=\"Huber Gradient\", ls=\"dashdot\")\nax.plot(x, pseudo_huber(delta, x), label=\"Pseudo-Huber\", ls=\"dotted\")\nax.plot(x, pseudo_huber_grad(delta, x), label=\"Pseudo-Huber Gradient\",\n        ls=\"solid\")\nax.legend(loc=\"upper center\")\nplt.show()",
            "outputs": []
        }
    ]
}