{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "3be1438d",
            "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": "3919874d",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import numpy as np\nfrom scipy.sparse import csgraph",
            "outputs": []
        },
        {
            "id": "181030a5",
            "cell_type": "markdown",
            "source": "Our first illustration is the symmetric graph\n",
            "metadata": {}
        },
        {
            "id": "63e353bf",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "G = np.arange(4) * np.arange(4)[:, np.newaxis]\nG",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[0, 0, 0, 0],\n       [0, 1, 2, 3],\n       [0, 2, 4, 6],\n       [0, 3, 6, 9]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "b1c78c36",
            "cell_type": "markdown",
            "source": "and its symmetric Laplacian matrix\n",
            "metadata": {}
        },
        {
            "id": "1cf1be30",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "csgraph.laplacian(G)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[ 0,  0,  0,  0],\n       [ 0,  5, -2, -3],\n       [ 0, -2,  8, -6],\n       [ 0, -3, -6,  9]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "e242b8ed",
            "cell_type": "markdown",
            "source": "The non-symmetric graph\n",
            "metadata": {}
        },
        {
            "id": "f2c5eee6",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "G = np.arange(9).reshape(3, 3)\nG",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[0, 1, 2],\n       [3, 4, 5],\n       [6, 7, 8]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "780e94f3",
            "cell_type": "markdown",
            "source": "has different row- and column sums, resulting in two varieties\nof the Laplacian matrix, using an in-degree, which is the default\n",
            "metadata": {}
        },
        {
            "id": "60990a85",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "L_in_degree = csgraph.laplacian(G)\nL_in_degree",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[ 9, -1, -2],\n       [-3,  8, -5],\n       [-6, -7,  7]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "fc868372",
            "cell_type": "markdown",
            "source": "or alternatively an out-degree\n",
            "metadata": {}
        },
        {
            "id": "ec8de0f2",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "L_out_degree = csgraph.laplacian(G, use_out_degree=True)\nL_out_degree",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[ 3, -1, -2],\n       [-3,  8, -5],\n       [-6, -7, 13]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "98adb5c6",
            "cell_type": "markdown",
            "source": "Constructing a symmetric Laplacian matrix, one can add the two as\n",
            "metadata": {}
        },
        {
            "id": "1c95b3d6",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "L_in_degree + L_out_degree.T",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[ 12,  -4,  -8],\n        [ -4,  16, -12],\n        [ -8, -12,  20]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "70db4d08",
            "cell_type": "markdown",
            "source": "or use the ``symmetrized=True`` option\n",
            "metadata": {}
        },
        {
            "id": "7a88e842",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "csgraph.laplacian(G, symmetrized=True)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[ 12,  -4,  -8],\n       [ -4,  16, -12],\n       [ -8, -12,  20]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "3248cfff",
            "cell_type": "markdown",
            "source": "that is equivalent to symmetrizing the original graph\n",
            "metadata": {}
        },
        {
            "id": "f4585e73",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "csgraph.laplacian(G + G.T)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[ 12,  -4,  -8],\n       [ -4,  16, -12],\n       [ -8, -12,  20]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "b8e354ee",
            "cell_type": "markdown",
            "source": "The goal of normalization is to make the non-zero diagonal entries\nof the Laplacian matrix to be all unit, also scaling off-diagonal\nentries correspondingly. The normalization can be done manually, e.g.,\n",
            "metadata": {}
        },
        {
            "id": "677c2aa7",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "G = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]])\nL, d = csgraph.laplacian(G, return_diag=True)\nL",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[ 2, -1, -1],\n       [-1,  2, -1],\n       [-1, -1,  2]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "aaa6c2f6",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "d",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([2, 2, 2])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "817b8baa",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "scaling = np.sqrt(d)\nscaling",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([1.41421356, 1.41421356, 1.41421356])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "2ced1663",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "(1/scaling)*L*(1/scaling)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[ 1. , -0.5, -0.5],\n       [-0.5,  1. , -0.5],\n       [-0.5, -0.5,  1. ]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "0526027f",
            "cell_type": "markdown",
            "source": "Or using ``normed=True`` option\n",
            "metadata": {}
        },
        {
            "id": "b5aafd86",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "L, d = csgraph.laplacian(G, return_diag=True, normed=True)\nL",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[ 1. , -0.5, -0.5],\n       [-0.5,  1. , -0.5],\n       [-0.5, -0.5,  1. ]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "648c97ad",
            "cell_type": "markdown",
            "source": "which now instead of the diagonal returns the scaling coefficients\n",
            "metadata": {}
        },
        {
            "id": "f8e21115",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "d",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([1.41421356, 1.41421356, 1.41421356])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "5ebf2cef",
            "cell_type": "markdown",
            "source": "Zero scaling coefficients are substituted with 1s, where scaling\nhas thus no effect, e.g.,\n",
            "metadata": {}
        },
        {
            "id": "8595432c",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "G = np.array([[0, 0, 0], [0, 0, 1], [0, 1, 0]])\nG",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[0, 0, 0],\n       [0, 0, 1],\n       [0, 1, 0]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "98e2e324",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "L, d = csgraph.laplacian(G, return_diag=True, normed=True)\nL",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[ 0., -0., -0.],\n       [-0.,  1., -1.],\n       [-0., -1.,  1.]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "6a71e39a",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "d",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([1., 1., 1.])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "487e075b",
            "cell_type": "markdown",
            "source": "Only the symmetric normalization is implemented, resulting\nin a symmetric Laplacian matrix if and only if its graph is symmetric\nand has all non-negative degrees, like in the examples above.\n\nThe output Laplacian matrix is by default a dense array or a sparse\narray or matrix inferring its class, shape, format, and dtype from\nthe input graph matrix:\n",
            "metadata": {}
        },
        {
            "id": "381b047c",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "G = np.array([[0, 1, 1], [1, 0, 1], [1, 1, 0]]).astype(np.float32)\nG",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[0., 1., 1.],\n       [1., 0., 1.],\n       [1., 1., 0.]], dtype=float32)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "47fa4431",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "csgraph.laplacian(G)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[ 2., -1., -1.],\n       [-1.,  2., -1.],\n       [-1., -1.,  2.]], dtype=float32)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "c0822d31",
            "cell_type": "markdown",
            "source": "but can alternatively be generated matrix-free as a LinearOperator:\n",
            "metadata": {}
        },
        {
            "id": "ee433b4d",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "L = csgraph.laplacian(G, form=\"lo\")\nL",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "<3x3 _CustomLinearOperator with dtype=float32>"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "fe57e33f",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "L(np.eye(3))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[ 2., -1., -1.],\n       [-1.,  2., -1.],\n       [-1., -1.,  2.]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "a0ea000b",
            "cell_type": "markdown",
            "source": "or as a lambda-function:\n",
            "metadata": {}
        },
        {
            "id": "4212c99e",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "L = csgraph.laplacian(G, form=\"function\")\nL",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "<function _laplace.<locals>.<lambda> at 0x0000012AE6F5A598>"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "5732797e",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "L(np.eye(3))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[ 2., -1., -1.],\n       [-1.,  2., -1.],\n       [-1., -1.,  2.]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "b8c5ad4d",
            "cell_type": "markdown",
            "source": "The Laplacian matrix is used for\nspectral data clustering and embedding\nas well as for spectral graph partitioning.\nOur final example illustrates the latter\nfor a noisy directed linear graph.\n",
            "metadata": {}
        },
        {
            "id": "6814266e",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "from scipy.sparse import diags_array, random_array\nfrom scipy.sparse.linalg import lobpcg",
            "outputs": []
        },
        {
            "id": "4171b0ae",
            "cell_type": "markdown",
            "source": "Create a directed linear graph with ``N=35`` vertices\nusing a sparse adjacency matrix ``G``:\n",
            "metadata": {}
        },
        {
            "id": "ca529351",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "N = 35\nG = diags_array(np.ones(N - 1), offsets=1, format=\"csr\")",
            "outputs": []
        },
        {
            "id": "47d1362a",
            "cell_type": "markdown",
            "source": "Fix a random seed ``rng`` and add a random sparse noise to the graph ``G``:\n",
            "metadata": {}
        },
        {
            "id": "5b4085de",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "rng = np.random.default_rng()\nG += 1e-2 * random_array((N, N), density=0.1, rng=rng)",
            "outputs": []
        },
        {
            "id": "ed8587fa",
            "cell_type": "markdown",
            "source": "Set initial approximations for eigenvectors:\n",
            "metadata": {}
        },
        {
            "id": "4e73a9e8",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "X = rng.random((N, 2))",
            "outputs": []
        },
        {
            "id": "e1e08c6b",
            "cell_type": "markdown",
            "source": "The constant vector of ones is always a trivial eigenvector\nof the non-normalized Laplacian to be filtered out:\n",
            "metadata": {}
        },
        {
            "id": "204d1e50",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "Y = np.ones((N, 1))",
            "outputs": []
        },
        {
            "id": "a05c253f",
            "cell_type": "markdown",
            "source": "Alternating (1) the sign of the graph weights allows determining\nlabels for spectral max- and min- cuts in a single loop.\nSince the graph is undirected, the option ``symmetrized=True``\nmust be used in the construction of the Laplacian.\nThe option ``normed=True`` cannot be used in (2) for the negative weights\nhere as the symmetric normalization evaluates square roots.\nThe option ``form=\"lo\"`` in (2) is matrix-free, i.e., guarantees\na fixed memory footprint and read-only access to the graph.\nCalling the eigenvalue solver ``lobpcg`` (3) computes the Fiedler vector\nthat determines the labels as the signs of its components in (5).\nSince the sign in an eigenvector is not deterministic and can flip,\nwe fix the sign of the first component to be always +1 in (4).\n",
            "metadata": {}
        },
        {
            "id": "75316647",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "for cut in [\"max\", \"min\"]:\n    G = -G  # 1.\n    L = csgraph.laplacian(G, symmetrized=True, form=\"lo\")  # 2.\n    _, eves = lobpcg(L, X, Y=Y, largest=False, tol=1e-2)  # 3.\n    eves *= np.sign(eves[0, 0])  # 4.\n    print(cut + \"-cut labels:\\n\", 1 * (eves[:, 0]>0))  # 5.",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "max-cut labels:\n[1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1]\nmin-cut labels:\n[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "4ad5fc02",
            "cell_type": "markdown",
            "source": "As anticipated for a (slightly noisy) linear graph,\nthe max-cut strips all the edges of the graph coloring all\nodd vertices into one color and all even vertices into another one,\nwhile the balanced min-cut partitions the graph\nin the middle by deleting a single edge.\nBoth determined partitions are optimal.",
            "metadata": {}
        }
    ]
}