Update tellurium.py
Added nullspace and rank helper functions
@@ -10,6 +10,8 @@ import warnings
import antimony
import matplotlib.pyplot as plt
+import numpy as np
+from numpy.linalg import svd
# check availability of property cycler (matplotlib 1.5ish)
if True: # create dummy scope
@@ -960,6 +962,87 @@ class ExtendedRoadRunner(roadrunner.RoadRunner):
return s
+def rank(A, atol=1e-13, rtol=0):
+ """Estimate the rank (i.e. the dimension of the nullspace) of a matrix.
+
+ The algorithm used by this function is based on the singular value
+ decomposition of `A`.
+
+ Parameters
+ ----------
+ A : ndarray
+ A should be at most 2-D. A 1-D array with length n will be treated
+ as a 2-D with shape (1, n)
+ atol : float
+ The absolute tolerance for a zero singular value. Singular values
+ smaller than `atol` are considered to be zero.
+ rtol : float
+ The relative tolerance. Singular values less than rtol*smax are
+ considered to be zero, where smax is the largest singular value.
+
+ If both `atol` and `rtol` are positive, the combined tolerance is the
+ maximum of the two; that is::
+ tol = max(atol, rtol * smax)
+ Singular values smaller than `tol` are considered to be zero.
+
+ Return value
+ ------------
+ r : int
+ The estimated rank of the matrix.
+
+ See also
+ --------
+ numpy.linalg.matrix_rank
+ matrix_rank is basically the same as this function, but it does not
+ provide the option of the absolute tolerance.
+ """
+
+ A = np.atleast_2d(A)
+ s = svd(A, compute_uv=False)
+ tol = max(atol, rtol * s[0])
+ rank = int((s >= tol).sum())
+ return rank
+
+def nullspace(A, atol=1e-13, rtol=0):
+ """Compute an approximate basis for the nullspace of A.
+
+ The algorithm used by this function is based on the singular value
+ decomposition of `A`.
+
+ Parameters
+ ----------
+ A : ndarray
+ A should be at most 2-D. A 1-D array with length k will be treated
+ as a 2-D with shape (1, k)
+ atol : float
+ The absolute tolerance for a zero singular value. Singular values
+ smaller than `atol` are considered to be zero.
+ rtol : float
+ The relative tolerance. Singular values less than rtol*smax are
+ considered to be zero, where smax is the largest singular value.
+
+ If both `atol` and `rtol` are positive, the combined tolerance is the
+ maximum of the two; that is::
+ tol = max(atol, rtol * smax)
+ Singular values smaller than `tol` are considered to be zero.
+
+ Return value
+ ------------
+ ns : ndarray
+ If `A` is an array with shape (m, k), then `ns` will be an array
+ with shape (k, n), where n is the estimated dimension of the
+ nullspace of `A`. The columns of `ns` are a basis for the
+ nullspace; each element in numpy.dot(A, ns) will be approximately
+ zero.
+ """
+
+ A = np.atleast_2d(A)
+ u, s, vh = svd(A)
+ tol = max(atol, rtol * s[0])
+ nnz = (s >= tol).sum()
+ ns = vh[nnz:].conj().T
+ return ns
+
# ---------------------------------------------------------------
# End of routines
# ---------------------------------------------------------------
