sympycore/matrices/functions.py

   1 """ Provides matrix utility functions.
   2 """
   3 #
   4 # Author: Pearu Peterson
   5 # Created: March 2008
   6 #
   7
   8 __docformat__ = "restructuredtext"
   9 __all__ = ['eye', 'concatenate', 'jacobian']
  10
  11 from ..utils import MATRIX, MATRIX_DICT
  12 from .algebra import MatrixDict, Matrix, MatrixBase
  13
  14 def jacobian(expr_list, var_list):
  15     """ Return a jacobian matrix of functions in expr_list with
  16     respect to variables in var_list.
  17     """
  18     m, n = len(expr_list), len(var_list)
  19     jac = Matrix(m, n)
  20     for i, e in enumerate(expr_list):
  21         for j, v in enumerate(var_list):
  22             jac[i,j] = e.diff(v)
  23     return jac
  24
  25 def eye(m, n=None, k=0):
  26     """ Return n x m matrix where the k-th diagonal is all ones,
  27     everything else is zeros.
  28     """
  29     if n is None:
  30         n = m
  31     d = {}
  32     if k<0:
  33         for i in range(min(m,n)+k):
  34             d[i-k,i] = 1
  35     else:
  36         for i in range(min(m,n)-k):
  37             d[i,i+k] = 1
  38     return MatrixDict(MATRIX(m, n, MATRIX_DICT), d)
  39
  40 def concatenate(*args, **kws):
  41     """ Join matrices together.
  42
  43     concatenate(m1, m2, ..) - join matrices m1, m2, .. together along rows
  44     concatenate(m1, m2, .., axis=1) - join matrices m1, m2, .. together along columns
  45     concatenate(m1, m2, .., diagonal=True) - join matrices m1, m2, .. together along diagonal
  46     """
  47     assert args,`args`
  48     assert set(kws).issubset(['axis','diagonal']) and len(kws)<=1,`kws`
  49     axis = kws.get('axis',0)
  50     diagonal = kws.get('diagonal',False)
  51     d = {}
  52     rows, cols = 0, 0
  53     if axis:
  54         for a in args:
  55             head, data = Matrix(a)[:].pair
  56             assert not (head.is_transpose or head.is_diagonal)
  57             m, n = head.shape
  58             for (i,j),x in data.items():
  59                 d[rows + i,j] = x
  60             rows += m
  61             cols = max(cols, n)
  62     elif diagonal:
  63         for a in args:
  64             head, data = Matrix(a)[:].pair
  65             assert not (head.is_transpose or head.is_diagonal)
  66             m, n = head.shape
  67             for (i,j),x in data.items():
  68                 d[rows + i, cols + j] = x
  69             rows += m
  70             cols += n
  71     else:
  72         for a in args:
  73             head, data = Matrix(a)[:].pair
  74             assert not (head.is_transpose or head.is_diagonal)
  75             m, n = head.shape
  76             for (i,j),x in data.items():
  77                 d[i,cols+j] = x
  78             cols += n
  79             rows = max(rows, m)
  80     return MatrixDict(MATRIX(rows, cols, MATRIX_DICT), d)