Merge pull request #143 from ShikharJ/DenseMatrix

Additions to DenseMatrixBase and MutableDenseMatrix

Author: isuruf

symengine/lib/symengine.pxd

@@ -769,6 +769,12 @@ cdef extern from "<symengine/matrix.h>" namespace "SymEngine":
         DenseMatrix(unsigned i, unsigned j) nogil
         DenseMatrix(unsigned i, unsigned j, const vec_basic &v) nogil
         void resize(unsigned i, unsigned j) nogil
+        void row_join(const DenseMatrix &B) nogil
+        void col_join(const DenseMatrix &B) nogil
+        void row_insert(const DenseMatrix &B, unsigned pos) nogil
+        void col_insert(const DenseMatrix &B, unsigned pos) nogil
+        void row_del(unsigned k) nogil
+        void col_del(unsigned k) nogil
 
     bool is_a_DenseMatrix "SymEngine::is_a<SymEngine::DenseMatrix>"(const MatrixBase &b) nogil
     DenseMatrix* static_cast_DenseMatrix "static_cast<SymEngine::DenseMatrix*>"(const MatrixBase *a)
@@ -792,6 +798,12 @@ cdef extern from "<symengine/matrix.h>" namespace "SymEngine":
     void diag(DenseMatrix &A, vec_basic &v, int k) nogil
     void ones(DenseMatrix &A) nogil
     void zeros(DenseMatrix &A) nogil
+    void row_exchange_dense(DenseMatrix &A, unsigned i, unsigned j) nogil
+    void row_mul_scalar_dense(DenseMatrix &A, unsigned i, RCP[const Basic] &c) nogil
+    void row_add_row_dense(DenseMatrix &A, unsigned i, unsigned j, RCP[const Basic] &c) nogil
+    void column_exchange_dense(DenseMatrix &A, unsigned i, unsigned j) nogil
+    void dot(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &C) nogil
+    void cross(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &C) nogil
 
 cdef extern from "<symengine/ntheory.h>" namespace "SymEngine":
     int probab_prime_p(const Integer &a, int reps)

symengine/lib/symengine_wrapper.pyx

@@ -2505,32 +2505,62 @@ cdef class DenseMatrixBase(MatrixBase):
         cdef DenseMatrixBase o = sympify(rhs)
         if self.rows != o.rows:
             raise ShapeError("`self` and `rhs` must have the same number of rows.")
-        cdef DenseMatrixBase result = self.__class__(self.rows, self.cols + o.cols)
-        cdef Basic e_
-        for i in range(self.rows):
-            for j in range(self.cols):
-                e_ = self._get(i, j)
-                deref(result.thisptr).set(i, j, e_.thisptr)
-        for i in range(o.rows):
-            for j in range(o.cols):
-                e_ = sympify(o._get(i, j))
-                deref(result.thisptr).set(i, j + self.cols, e_.thisptr)
-        return result
+        cdef DenseMatrixBase d = self.__class__(self)
+        deref(symengine.static_cast_DenseMatrix(d.thisptr)).row_join(deref(symengine.static_cast_DenseMatrix(o.thisptr)))
+        return d
 
     def col_join(self, bott):
         cdef DenseMatrixBase o = sympify(bott)
         if self.cols != o.cols:
             raise ShapeError("`self` and `rhs` must have the same number of columns.")
-        cdef DenseMatrixBase result = self.__class__(self.rows + o.rows, self.cols)
-        cdef Basic e_
-        for i in range(self.rows):
-            for j in range(self.cols):
-                e_ = self._get(i, j)
-                deref(result.thisptr).set(i, j, e_.thisptr)
-        for i in range(o.rows):
-            for j in range(o.cols):
-                e_ = sympify(o._get(i, j))
-                deref(result.thisptr).set(i + self.rows, j, e_.thisptr)
+        cdef DenseMatrixBase d = self.__class__(self)
+        deref(symengine.static_cast_DenseMatrix(d.thisptr)).col_join(deref(symengine.static_cast_DenseMatrix(o.thisptr)))
+        return d
+
+    def row_insert(self, pos, bott):
+        cdef DenseMatrixBase o = sympify(bott)
+        if pos < 0:
+            pos = self.rows + pos
+        if pos < 0:
+            pos = 0
+        elif pos > self.rows:
+            pos = self.rows
+        if self.cols != o.cols:
+            raise ShapeError("`self` and `other` must have the same number of columns.")
+        cdef DenseMatrixBase d = self.__class__(self)
+        deref(symengine.static_cast_DenseMatrix(d.thisptr)).row_insert(deref(symengine.static_cast_DenseMatrix(o.thisptr)), pos)
+        return d
+
+    def col_insert(self, pos, bott):
+        cdef DenseMatrixBase o = sympify(bott)
+        if pos < 0:
+            pos = self.cols + pos
+        if pos < 0:
+            pos = 0
+        elif pos > self.cols:
+            pos = self.cols
+        if self.rows != o.rows:
+            raise ShapeError("`self` and `other` must have the same number of rows.")
+        cdef DenseMatrixBase d = self.__class__(self)
+        deref(symengine.static_cast_DenseMatrix(d.thisptr)).col_insert(deref(symengine.static_cast_DenseMatrix(o.thisptr)), pos)
+        return d
+
+    def dot(self, b):
+        cdef DenseMatrixBase o = sympify(b)
+        cdef DenseMatrixBase result = self.__class__(self.rows, self.cols)
+        symengine.dot(deref(symengine.static_cast_DenseMatrix(self.thisptr)), deref(symengine.static_cast_DenseMatrix(o.thisptr)), deref(symengine.static_cast_DenseMatrix(result.thisptr)))
+        if len(result) == 1:
+            return result[0, 0]
+        else:
+            return result
+
+    def cross(self, b):
+        cdef DenseMatrixBase o = sympify(b)
+        if self.cols * self.rows != 3 or o.cols * o.rows != 3:
+            raise ShapeError("Dimensions incorrect for cross product: %s x %s" %
+                             ((self.rows, self.cols), (b.rows, b.cols)))
+        cdef DenseMatrixBase result = self.__class__(self.rows, self.cols)
+        symengine.cross(deref(symengine.static_cast_DenseMatrix(self.thisptr)), deref(symengine.static_cast_DenseMatrix(o.thisptr)), deref(symengine.static_cast_DenseMatrix(result.thisptr)))
         return result
 
     @property
@@ -2541,6 +2571,10 @@ cdef class DenseMatrixBase(MatrixBase):
     def cols(self):
         return self.ncols()
 
+    @property
+    def is_square(self):
+        return self.rows == self.cols
+
     def nrows(self):
         return deref(self.thisptr).nrows()
 
@@ -2589,6 +2623,12 @@ cdef class DenseMatrixBase(MatrixBase):
         # No error checking is done
         return c2py(deref(self.thisptr).get(i, j))
 
+    def col(self, j):
+        return self[:, j]
+
+    def row(self, i):
+        return self[i, :]
+
     def set(self, i, j, e):
         i, j = self._get_index(i, j)
         return self._set(i, j, e)
@@ -2665,6 +2705,13 @@ cdef class DenseMatrixBase(MatrixBase):
                     deref(out.thisptr).set(i, j, e_.thisptr)
         return out
 
+    def _applyfunc(self, f):
+        cdef int nr = self.nrows()
+        cdef int nc = self.ncols()
+        for i in range(nr):
+            for j in range(nc):
+                self._set(i, j, f(self._get(i, j)))
+
     def msubs(self, *args):
         cdef _DictBasic D = get_dict(*args)
         return self.applyfunc(lambda x: x.msubs(D))
@@ -2826,6 +2873,9 @@ cdef class DenseMatrixBase(MatrixBase):
     def tolist(self):
         return self[:]
 
+    def _mat(self):
+        return self
+
     def atoms(self, *types):
         if types:
             s = set()
@@ -2864,24 +2914,43 @@ class DenseMatrixBaseIter(object):
 cdef class MutableDenseMatrix(DenseMatrixBase):
 
     def col_swap(self, i, j):
-        for k in range(0, self.rows):
-            self[k, i], self[k, j] = self[k, j], self[k, i]
+        symengine.column_exchange_dense(deref(symengine.static_cast_DenseMatrix(self.thisptr)), i, j)
 
     def fill(self, value):
         for i in range(self.rows):
             for j in range(self.cols):
                 self[i, j] = value
 
     def row_swap(self, i, j):
-        for k in range(0, self.cols):
-            self[i, k], self[j, k] = self[j, k], self[i, k]
+        symengine.row_exchange_dense(deref(symengine.static_cast_DenseMatrix(self.thisptr)), i, j)
 
-    def _applyfunc(self, f):
-        cdef int nr = self.nrows()
-        cdef int nc = self.ncols()
-        for i in range(nr):
-            for j in range(nc):
-                self._set(i, j, f(self._get(i, j)))
+    def rowmul(self, i, c, *args):
+        cdef Basic _c = sympify(c)
+        symengine.row_mul_scalar_dense(deref(symengine.static_cast_DenseMatrix(self.thisptr)), i, _c.thisptr)
+        return self
+
+    def rowadd(self, i, j, c, *args):
+        cdef Basic _c = sympify(c)
+        symengine.row_add_row_dense(deref(symengine.static_cast_DenseMatrix(self.thisptr)), i, j, _c.thisptr)
+        return self
+
+    def row_del(self, i):
+        if i < -self.rows or i >= self.rows:
+            raise IndexError("Index out of range: 'i = %s', valid -%s <= i"
+                             " < %s" % (i, self.rows, self.rows))
+        if i < 0:
+            i += self.rows
+        deref(symengine.static_cast_DenseMatrix(self.thisptr)).row_del(i)
+        return self
+
+    def col_del(self, i):
+        if i < -self.cols or i >= self.cols:
+            raise IndexError("Index out of range: 'i=%s', valid -%s <= i < %s"
+                             % (i, self.cols, self.cols))
+        if i < 0:
+            i += self.cols
+        deref(symengine.static_cast_DenseMatrix(self.thisptr)).col_del(i)
+        return self
 
 Matrix = DenseMatrix = MutableDenseMatrix
 
@@ -2890,12 +2959,6 @@ cdef class ImmutableDenseMatrix(DenseMatrixBase):
     def __setitem__(self, key, value):
         raise TypeError("Cannot set values of {}".format(self.__class__))
 
-    def set(self, i, j, e):
-        raise TypeError("Cannot set values of {}".format(self.__class__))
-
-    def _set(self, i, j, e):
-        raise TypeError("Cannot set values of {}".format(self.__class__))
-
 ImmutableMatrix = ImmutableDenseMatrix
 
 cdef matrix_to_vec(DenseMatrixBase d, symengine.vec_basic& v):

symengine/tests/test_matrices.py

@@ -369,6 +369,95 @@ def test_row_swap():
     assert A == B
 
 
+def test_row_col_del():
+    e = DenseMatrix(3, 3, [1, 2, 3, 4, 5, 6, 7, 8, 9])
+    raises(IndexError, lambda: e.row_del(5))
+    raises(IndexError, lambda: e.row_del(-5))
+    raises(IndexError, lambda: e.col_del(5))
+    raises(IndexError, lambda: e.col_del(-5))
+
+    assert e.row_del(-1) == DenseMatrix([[1, 2, 3], [4, 5, 6]])
+    assert e.col_del(-1) == DenseMatrix([[1, 2], [4, 5]])
+
+    e = DenseMatrix(3, 3, [1, 2, 3, 4, 5, 6, 7, 8, 9])
+    assert e.row_del(1) == DenseMatrix([[1, 2, 3], [7, 8, 9]])
+    assert e.col_del(1) == DenseMatrix([[1, 3], [7, 9]])
+
+
+def test_row_join():
+    assert eye(3).row_join(DenseMatrix([7, 7, 7])) == \
+           DenseMatrix([[1, 0, 0, 7],
+                        [0, 1, 0, 7],
+                        [0, 0, 1, 7]])
+
+
+def test_col_join():
+    assert eye(3).col_join(DenseMatrix([[7, 7, 7]])) == \
+           DenseMatrix([[1, 0, 0],
+                        [0, 1, 0],
+                        [0, 0, 1],
+                        [7, 7, 7]])
+
+
+def test_row_insert():
+    M = zeros(3)
+    V = ones(1, 3)
+    assert M.row_insert(1, V) == DenseMatrix([[0, 0, 0],
+                                              [1, 1, 1],
+                                              [0, 0, 0],
+                                              [0, 0, 0]])
+
+
+def test_col_insert():
+    M = zeros(3)
+    V = ones(3, 1)
+    assert M.col_insert(1, V) == DenseMatrix([[0, 1, 0, 0],
+                                              [0, 1, 0, 0],
+                                              [0, 1, 0, 0]])
+
+
+def test_rowmul():
+    M = ones(3)
+    assert M.rowmul(2, 2) == DenseMatrix([[1, 1, 1],
+                                          [1, 1, 1], 
+                                          [2, 2, 2]])
+
+
+def test_rowadd():
+    M = ones(3)
+    assert M.rowadd(2, 1, 1) == DenseMatrix([[1, 1, 1],
+                                             [1, 1, 1], 
+                                             [2, 2, 2]])
+
+
+def test_row_col():
+    m = DenseMatrix(3, 3, [1, 2, 3, 4, 5, 6, 7, 8, 9])
+    assert m.row(0) == DenseMatrix(1, 3, [1, 2, 3])
+    assert m.col(0) == DenseMatrix(3, 1, [1, 4, 7])
+
+
+def test_is_square():
+    m = DenseMatrix([[1],[1]])
+    m2 = DenseMatrix([[2, 2], [2, 2]])
+    assert not m.is_square
+    assert m2.is_square
+
+
+def test_dot():
+    A = DenseMatrix(2, 3, [1, 2, 3, 4, 5, 6])
+    B = DenseMatrix(2, 1, [7, 8])
+    assert A.dot(B) == DenseMatrix(1, 3, [39, 54, 69])
+    assert ones(1, 3).dot(ones(3, 1)) == 3
+
+
+def test_cross():
+    M = DenseMatrix(1, 3, [1, 2, 3])
+    V = DenseMatrix(1, 3, [3, 4, 5])
+    assert M.cross(V) == DenseMatrix(1, 3, [-2, 4, -2])
+    raises(ShapeError, lambda:
+        DenseMatrix(1, 2, [1, 1]).cross(DenseMatrix(1, 2, [1, 1])))
+
+
 def test_immutablematrix():
     A = ImmutableMatrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
     assert A.shape == (3, 3)

symengine_version.txt

@@ -1 +1 @@
-d13fec95c651bbce195988a6d9a146e9b726b2c2
+978dfd1656d5dbc722f5dc448bbe96e97b2d7be9