{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "526dc313",
            "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": "55b9c4e3",
            "cell_type": "markdown",
            "source": "In this example we find a minimum of the Rosenbrock function without bounds\non independent variables.\n",
            "metadata": {}
        },
        {
            "id": "54d7b463",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import numpy as np\ndef fun_rosenbrock(x):\n    return np.array([10 * (x[1] - x[0]**2), (1 - x[0])])",
            "outputs": []
        },
        {
            "id": "0d66157c",
            "cell_type": "markdown",
            "source": "Notice that we only provide the vector of the residuals. The algorithm\nconstructs the cost function as a sum of squares of the residuals, which\ngives the Rosenbrock function. The exact minimum is at ``x = [1.0, 1.0]``.\n",
            "metadata": {}
        },
        {
            "id": "bb51f098",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.optimize import least_squares\nx0_rosenbrock = np.array([2, 2])\nres_1 = least_squares(fun_rosenbrock, x0_rosenbrock)\nres_1.x",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([ 1.,  1.])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "88541608",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res_1.cost",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "9.8669242910846867e-30"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "17657e26",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res_1.optimality",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "8.8928864934219529e-14"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "bd28cf2b",
            "cell_type": "markdown",
            "source": "We now constrain the variables, in such a way that the previous solution\nbecomes infeasible. Specifically, we require that ``x[1] >= 1.5``, and\n``x[0]`` left unconstrained. To this end, we specify the `bounds` parameter\nto `least_squares` in the form ``bounds=([-np.inf, 1.5], np.inf)``.\n\nWe also provide the analytic Jacobian:\n",
            "metadata": {}
        },
        {
            "id": "dc98ff94",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def jac_rosenbrock(x):\n    return np.array([\n        [-20 * x[0], 10],\n        [-1, 0]])",
            "outputs": []
        },
        {
            "id": "71415ea7",
            "cell_type": "markdown",
            "source": "Putting this all together, we see that the new solution lies on the bound:\n",
            "metadata": {}
        },
        {
            "id": "8f66f828",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res_2 = least_squares(fun_rosenbrock, x0_rosenbrock, jac_rosenbrock,\n                      bounds=([-np.inf, 1.5], np.inf))\nres_2.x",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([ 1.22437075,  1.5       ])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "9fd56fd0",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res_2.cost",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "0.025213093946805685"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "d060e48e",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res_2.optimality",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "1.5885401433157753e-07"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "6dfc7242",
            "cell_type": "markdown",
            "source": "Now we solve a system of equations (i.e., the cost function should be zero\nat a minimum) for a Broyden tridiagonal vector-valued function of 100000\nvariables:\n",
            "metadata": {}
        },
        {
            "id": "1a50a0a6",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def fun_broyden(x):\n    f = (3 - x) * x + 1\n    f[1:] -= x[:-1]\n    f[:-1] -= 2 * x[1:]\n    return f",
            "outputs": []
        },
        {
            "id": "d41cd6c2",
            "cell_type": "markdown",
            "source": "The corresponding Jacobian matrix is sparse. We tell the algorithm to\nestimate it by finite differences and provide the sparsity structure of\nJacobian to significantly speed up this process.\n",
            "metadata": {}
        },
        {
            "id": "cb079909",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.sparse import lil_array\ndef sparsity_broyden(n):\n    sparsity = lil_array((n, n), dtype=int)\n    i = np.arange(n)\n    sparsity[i, i] = 1\n    i = np.arange(1, n)\n    sparsity[i, i - 1] = 1\n    i = np.arange(n - 1)\n    sparsity[i, i + 1] = 1\n    return sparsity\n\nn = 100000\nx0_broyden = -np.ones(n)\n\nres_3 = least_squares(fun_broyden, x0_broyden,\n                      jac_sparsity=sparsity_broyden(n))\nres_3.cost",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "4.5687069299604613e-23"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "af82934d",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res_3.optimality",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "1.1650454296851518e-11"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "8805125d",
            "cell_type": "markdown",
            "source": "Let's also solve a curve fitting problem using robust loss function to\ntake care of outliers in the data. Define the model function as\n``y = a + b * exp(c * t)``, where t is a predictor variable, y is an\nobservation and a, b, c are parameters to estimate.\n\nFirst, define the function which generates the data with noise and\noutliers, define the model parameters, and generate data:\n",
            "metadata": {}
        },
        {
            "id": "3e11b8b2",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from numpy.random import default_rng\nrng = default_rng()\ndef gen_data(t, a, b, c, noise=0., n_outliers=0, seed=None):\n    rng = default_rng(seed)\n\n    y = a + b * np.exp(t * c)\n\n    error = noise * rng.standard_normal(t.size)\n    outliers = rng.integers(0, t.size, n_outliers)\n    error[outliers] *= 10\n\n    return y + error\n\na = 0.5\nb = 2.0\nc = -1\nt_min = 0\nt_max = 10\nn_points = 15\n\nt_train = np.linspace(t_min, t_max, n_points)\ny_train = gen_data(t_train, a, b, c, noise=0.1, n_outliers=3)",
            "outputs": []
        },
        {
            "id": "c981b168",
            "cell_type": "markdown",
            "source": "Define function for computing residuals and initial estimate of\nparameters.\n",
            "metadata": {}
        },
        {
            "id": "27dd7b3e",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def fun(x, t, y):\n    return x[0] + x[1] * np.exp(x[2] * t) - y\n\nx0 = np.array([1.0, 1.0, 0.0])",
            "outputs": []
        },
        {
            "id": "322659a3",
            "cell_type": "markdown",
            "source": "Compute a standard least-squares solution:\n",
            "metadata": {}
        },
        {
            "id": "093d8d8a",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res_lsq = least_squares(fun, x0, args=(t_train, y_train))",
            "outputs": []
        },
        {
            "id": "b0e0557f",
            "cell_type": "markdown",
            "source": "Now compute two solutions with two different robust loss functions. The\nparameter `f_scale` is set to 0.1, meaning that inlier residuals should\nnot significantly exceed 0.1 (the noise level used).\n",
            "metadata": {}
        },
        {
            "id": "0e753eec",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "res_soft_l1 = least_squares(fun, x0, loss='soft_l1', f_scale=0.1,\n                            args=(t_train, y_train))\nres_log = least_squares(fun, x0, loss='cauchy', f_scale=0.1,\n                        args=(t_train, y_train))",
            "outputs": []
        },
        {
            "id": "daa87a36",
            "cell_type": "markdown",
            "source": "And, finally, plot all the curves. We see that by selecting an appropriate\n`loss`  we can get estimates close to optimal even in the presence of\nstrong outliers. But keep in mind that generally it is recommended to try\n'soft_l1' or 'huber' losses first (if at all necessary) as the other two\noptions may cause difficulties in optimization process.\n",
            "metadata": {}
        },
        {
            "id": "fe810ac1",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "t_test = np.linspace(t_min, t_max, n_points * 10)\ny_true = gen_data(t_test, a, b, c)\ny_lsq = gen_data(t_test, *res_lsq.x)\ny_soft_l1 = gen_data(t_test, *res_soft_l1.x)\ny_log = gen_data(t_test, *res_log.x)\n\nimport matplotlib.pyplot as plt\nplt.plot(t_train, y_train, 'o')\nplt.plot(t_test, y_true, 'k', linewidth=2, label='true')\nplt.plot(t_test, y_lsq, label='linear loss')\nplt.plot(t_test, y_soft_l1, label='soft_l1 loss')\nplt.plot(t_test, y_log, label='cauchy loss')\nplt.xlabel(\"t\")\nplt.ylabel(\"y\")\nplt.legend()\nplt.show()",
            "outputs": []
        },
        {
            "id": "514f1926",
            "cell_type": "markdown",
            "source": "In the next example, we show how complex-valued residual functions of\ncomplex variables can be optimized with ``least_squares()``. Consider the\nfollowing function:\n",
            "metadata": {}
        },
        {
            "id": "4fc2a802",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def f(z):\n    return z - (0.5 + 0.5j)",
            "outputs": []
        },
        {
            "id": "2ab52e55",
            "cell_type": "markdown",
            "source": "We wrap it into a function of real variables that returns real residuals\nby simply handling the real and imaginary parts as independent variables:\n",
            "metadata": {}
        },
        {
            "id": "0e19d45f",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "def f_wrap(x):\n    fx = f(x[0] + 1j*x[1])\n    return np.array([fx.real, fx.imag])",
            "outputs": []
        },
        {
            "id": "6a85f075",
            "cell_type": "markdown",
            "source": "Thus, instead of the original m-D complex function of n complex\nvariables we optimize a 2m-D real function of 2n real variables:\n",
            "metadata": {}
        },
        {
            "id": "b1e98740",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.optimize import least_squares\nres_wrapped = least_squares(f_wrap, (0.1, 0.1), bounds=([0, 0], [1, 1]))\nz = res_wrapped.x[0] + res_wrapped.x[1]*1j\nz",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "(0.49999999999925893+0.49999999999925893j)"
                    },
                    "execution_count": null
                }
            ]
        }
    ]
}