{
    "nbformat": 4,
    "nbformat_minor": 5,
    "metadata": {
        "kernelspec": {
            "display_name": "Python",
            "language": "python",
            "name": "python"
        },
        "language_info": {
            "name": "python"
        }
    },
    "cells": [
        {
            "id": "6f5ceace",
            "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": "55590aaa",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "import numpy as np\nfrom scipy.sparse.linalg import LaplacianNd\nfrom scipy.sparse import diags_array, csgraph\nfrom scipy.linalg import eigvalsh",
            "outputs": []
        },
        {
            "id": "9306e60b",
            "cell_type": "markdown",
            "source": "The one-dimensional Laplacian demonstrated below for pure Neumann boundary\nconditions on a regular grid with ``n=6`` grid points is exactly the\nnegative graph Laplacian for the undirected linear graph with ``n``\nvertices using the sparse adjacency matrix ``G`` represented by the\nfamous tri-diagonal matrix:\n",
            "metadata": {}
        },
        {
            "id": "a846664d",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "n = 6\nG = diags_array(np.ones(n - 1), offsets=1, format='csr')\nLf = csgraph.laplacian(G, symmetrized=True, form='function')\ngrid_shape = (n, )\nlap = LaplacianNd(grid_shape, boundary_conditions='neumann')\nnp.array_equal(lap.matmat(np.eye(n)), -Lf(np.eye(n)))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "2862455e",
            "cell_type": "markdown",
            "source": "Since all matrix entries of the Laplacian are integers, ``'int8'`` is\nthe default dtype for storing matrix representations.\n",
            "metadata": {}
        },
        {
            "id": "d839787d",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "lap.tosparse()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "<DIAgonal sparse array of dtype 'int8'\n    with 16 stored elements (3 diagonals) and shape (6, 6)>"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "09c32f3f",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "lap.toarray()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[-1,  1,  0,  0,  0,  0],\n       [ 1, -2,  1,  0,  0,  0],\n       [ 0,  1, -2,  1,  0,  0],\n       [ 0,  0,  1, -2,  1,  0],\n       [ 0,  0,  0,  1, -2,  1],\n       [ 0,  0,  0,  0,  1, -1]], dtype=int8)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "0a630ac6",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "np.array_equal(lap.matmat(np.eye(n)), lap.toarray())",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "1f2ee06b",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "np.array_equal(lap.tosparse().toarray(), lap.toarray())",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "836c54df",
            "cell_type": "markdown",
            "source": "Any number of extreme eigenvalues and/or eigenvectors can be computed.\n",
            "metadata": {}
        },
        {
            "id": "c6f44885",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "lap = LaplacianNd(grid_shape, boundary_conditions='periodic')\nlap.eigenvalues()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([-4., -3., -3., -1., -1.,  0.])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "38180e2a",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "lap.eigenvalues()[-2:]",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([-1.,  0.])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "36f2c9e3",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "lap.eigenvalues(2)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([-1.,  0.])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "7c11c963",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "lap.eigenvectors(1)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[0.40824829],\n       [0.40824829],\n       [0.40824829],\n       [0.40824829],\n       [0.40824829],\n       [0.40824829]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "e9331cdd",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "lap.eigenvectors(2)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[ 0.5       ,  0.40824829],\n       [ 0.        ,  0.40824829],\n       [-0.5       ,  0.40824829],\n       [-0.5       ,  0.40824829],\n       [ 0.        ,  0.40824829],\n       [ 0.5       ,  0.40824829]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "fa8ae3f8",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "lap.eigenvectors()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[ 0.40824829,  0.28867513,  0.28867513,  0.5       ,  0.5       ,\n         0.40824829],\n       [-0.40824829, -0.57735027, -0.57735027,  0.        ,  0.        ,\n         0.40824829],\n       [ 0.40824829,  0.28867513,  0.28867513, -0.5       , -0.5       ,\n         0.40824829],\n       [-0.40824829,  0.28867513,  0.28867513, -0.5       , -0.5       ,\n         0.40824829],\n       [ 0.40824829, -0.57735027, -0.57735027,  0.        ,  0.        ,\n         0.40824829],\n       [-0.40824829,  0.28867513,  0.28867513,  0.5       ,  0.5       ,\n         0.40824829]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "938400df",
            "cell_type": "markdown",
            "source": "The two-dimensional Laplacian is illustrated on a regular grid with\n``grid_shape = (2, 3)`` points in each dimension.\n",
            "metadata": {}
        },
        {
            "id": "34afe8b1",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "grid_shape = (2, 3)\nn = np.prod(grid_shape)",
            "outputs": []
        },
        {
            "id": "4b6917fc",
            "cell_type": "markdown",
            "source": "Numeration of grid points is as follows:\n",
            "metadata": {}
        },
        {
            "id": "852db3e3",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "np.arange(n).reshape(grid_shape + (-1,))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[[0],\n        [1],\n        [2]],\n<BLANKLINE>\n       [[3],\n        [4],\n        [5]]])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "61201c72",
            "cell_type": "markdown",
            "source": "Each of the boundary conditions ``'dirichlet'``, ``'periodic'``, and\n``'neumann'`` is illustrated separately; with ``'dirichlet'``\n",
            "metadata": {}
        },
        {
            "id": "199f680a",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "lap = LaplacianNd(grid_shape, boundary_conditions='dirichlet')\nlap.tosparse()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "<Compressed Sparse Row sparse array of dtype 'int8'\n    with 20 stored elements and shape (6, 6)>"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "a48ea945",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "lap.toarray()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[-4,  1,  0,  1,  0,  0],\n       [ 1, -4,  1,  0,  1,  0],\n       [ 0,  1, -4,  0,  0,  1],\n       [ 1,  0,  0, -4,  1,  0],\n       [ 0,  1,  0,  1, -4,  1],\n       [ 0,  0,  1,  0,  1, -4]], dtype=int8)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "2317c3bf",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "np.array_equal(lap.matmat(np.eye(n)), lap.toarray())",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "17355eee",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "np.array_equal(lap.tosparse().toarray(), lap.toarray())",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "40dda065",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "lap.eigenvalues()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([-6.41421356, -5.        , -4.41421356, -3.58578644, -3.        ,\n       -1.58578644])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "300209d0",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "eigvals = eigvalsh(lap.toarray().astype(np.float64))\nnp.allclose(lap.eigenvalues(), eigvals)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "54470585",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "np.allclose(lap.toarray() @ lap.eigenvectors(),\n            lap.eigenvectors() @ np.diag(lap.eigenvalues()))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "46242f20",
            "cell_type": "markdown",
            "source": "with ``'periodic'``\n",
            "metadata": {}
        },
        {
            "id": "f56dcc56",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "lap = LaplacianNd(grid_shape, boundary_conditions='periodic')\nlap.tosparse()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "<Compressed Sparse Row sparse array of dtype 'int8'\n    with 24 stored elements and shape (6, 6)>"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "841483d1",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "lap.toarray()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "    array([[-4,  1,  1,  2,  0,  0],\n           [ 1, -4,  1,  0,  2,  0],\n           [ 1,  1, -4,  0,  0,  2],\n           [ 2,  0,  0, -4,  1,  1],\n           [ 0,  2,  0,  1, -4,  1],\n           [ 0,  0,  2,  1,  1, -4]], dtype=int8)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "c628d9e0",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "np.array_equal(lap.matmat(np.eye(n)), lap.toarray())",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "4e21f801",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "np.array_equal(lap.tosparse().toarray(), lap.toarray())",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "92537ec7",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "lap.eigenvalues()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([-7., -7., -4., -3., -3.,  0.])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "2c5c056f",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "eigvals = eigvalsh(lap.toarray().astype(np.float64))\nnp.allclose(lap.eigenvalues(), eigvals)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "3a8f78df",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "np.allclose(lap.toarray() @ lap.eigenvectors(),\n            lap.eigenvectors() @ np.diag(lap.eigenvalues()))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "c34a888f",
            "cell_type": "markdown",
            "source": "and with ``'neumann'``\n",
            "metadata": {}
        },
        {
            "id": "094e3180",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "lap = LaplacianNd(grid_shape, boundary_conditions='neumann')\nlap.tosparse()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "<Compressed Sparse Row sparse array of dtype 'int8'\n    with 20 stored elements and shape (6, 6)>"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "bba79b69",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "lap.toarray()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([[-2,  1,  0,  1,  0,  0],\n       [ 1, -3,  1,  0,  1,  0],\n       [ 0,  1, -2,  0,  0,  1],\n       [ 1,  0,  0, -2,  1,  0],\n       [ 0,  1,  0,  1, -3,  1],\n       [ 0,  0,  1,  0,  1, -2]], dtype=int8)"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "aa4ce715",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "np.array_equal(lap.matmat(np.eye(n)), lap.toarray())",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "0c41fb21",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "np.array_equal(lap.tosparse().toarray(), lap.toarray())",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "449d3fa6",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "lap.eigenvalues()",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "array([-5., -3., -3., -2., -1.,  0.])"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "429baae1",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "eigvals = eigvalsh(lap.toarray().astype(np.float64))\nnp.allclose(lap.eigenvalues(), eigvals)",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        },
        {
            "id": "8975c6ae",
            "cell_type": "code",
            "metadata": {},
            "execution_count": null,
            "source": "np.allclose(lap.toarray() @ lap.eigenvectors(),\n            lap.eigenvectors() @ np.diag(lap.eigenvalues()))",
            "outputs": [
                {
                    "output_type": "execute_result",
                    "metadata": {},
                    "data": {
                        "text/plain": "True"
                    },
                    "execution_count": null
                }
            ]
        }
    ]
}