{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Rank and nullspace demo\n\nThis demo shows an example on how to use ``numpy`` in ``itom``.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nfrom numpy.linalg import svd\nfrom numpy.typing import ArrayLike"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Function to estimate the rank (i.e. the dimension of the nullspace) of a matrix.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def rank(A: ArrayLike, atol: float = 1e-13, rtol: int = 0) -> int:\n    \"\"\"Estimate the rank (i.e. the dimension of the nullspace) of a matrix.\n\n       The algorithm used by this function is based on the singular value\n       decomposition of `A`.\n       If both `atol` and `rtol` are positive, the combined tolerance is the\n       maximum of the two; that is::\n           tol = max(atol, rtol * smax)\n       Singular values smaller than `tol` are considered to be zero.\n\n    Args:\n        A (ArrayLike): A should be at most 2-D.  A 1-D array with length n will be treated\n            as a 2-D with shape (1, n)\n        atol (float, optional): The absolute tolerance for a zero singular value.  Singular values\n            smaller than `atol` are considered to be zero.\n        rtol (int, optional): The relative tolerance.  Singular values less than rtol*smax are\n            considered to be zero, where smax is the largest singular value.\n\n    Returns:\n        int: The estimated rank of the matrix.\n\n    See also\n    --------\n    numpy.linalg.matrix_rank\n        matrix_rank is basically the same as this function, but it does not\n        provide the option of the absolute tolerance.\"\"\"\n\n    A = np.atleast_2d(A)\n    s = svd(A, compute_uv=False)\n    tol = max(atol, rtol * s[0])\n    rank = int((s >= tol).sum())\n    return rank"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Function to compute an approximate basis for the nullspace of A.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def nullspace(A: ArrayLike, atol: float = 1e-13, rtol: int = 0) -> ArrayLike:\n    \"\"\"Compute an approximate basis for the nullspace of A.\n\n       The algorithm used by this function is based on the singular value\n       decomposition of `A`.\n\n       If both `atol` and `rtol` are positive, the combined tolerance is the\n       maximum of the two; that is::\n           tol = max(atol, rtol * smax)\n       Singular values smaller than `tol` are considered to be zero.\n\n    Args:\n        A (ArrayLike): A should be at most 2-D.  A 1-D array with length k will be treated\n            as a 2-D with shape (1, k)\n        atol (float, optional): The absolute tolerance for a zero singular value.  Singular values\n            smaller than `atol` are considered to be zero.\n        rtol (int, optional): The relative tolerance.  Singular values less than rtol*smax are\n            considered to be zero, where smax is the largest singular value.\n\n    Returns:\n        ArrayLike: If `A` is an array with shape (m, k), then `ns` will be an array\n            with shape (k, n), where n is the estimated dimension of the\n            nullspace of `A`.  The columns of `ns` are a basis for the\n            nullspace; each element in numpy.dot(A, ns) will be approximately\n            zero.\n    \"\"\"\n\n    A = np.atleast_2d(A)\n    u, s, vh = svd(A)\n    tol = max(atol, rtol * s[0])\n    nnz = (s >= tol).sum()\n    ns = vh[nnz:].conj().T\n    return ns"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Function to check rank and nullspace of the matrix. \n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def checkit(mat):\n    \"\"\"This method calculates the rank and nullspace of matrix mat. The results\n    are printed.\"\"\"\n    print(\"mat:\")\n    print(mat)\n    r = rank(mat)\n    print(\"rank is\", r)\n    ns = nullspace(mat)\n    print(\"nullspace:\")\n    print(ns)\n    if ns.size > 0:\n        res = np.abs(np.dot(mat, ns)).max()\n        print(\"max residual is\", res)\n\n\nprint(\"-\" * 25)\n\n# checks rank and nulllspace of a. The rank is limited since the 3rd row\n# is equal to 2x the 2nd row minus 1x the 1st row.\na = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]])\ncheckit(a)\n\nb = 2\n\nprint(\"-\" * 25)\n\n# full rank matrix\na = np.array([[0.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]])\ncheckit(a)\n\nprint(\"-\" * 25)\n\n# the rank is 2 since the matrix has only two rows\na = np.array([[0.0, 1.0, 2.0, 4.0], [1.0, 2.0, 3.0, 4.0]])\ncheckit(a)\n\nprint(\"-\" * 25)\n\n# check the rank and nullspace of a complex matrix\na = np.array(\n    [\n        [1.0, 1.0j, 2.0 + 2.0j],\n        [1.0j, -1.0, -2.0 + 2.0j],\n        [0.5, 0.5j, 1.0 + 1.0j],\n    ]\n)\ncheckit(a)\n\nprint(\"-\" * 25)"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.8.10"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}