developping vector

symla/dot.py

@@ -0,0 +1,90 @@
+from operator  import mul, add
+from functools import reduce
+
+from sympy.core               import Basic
+from sympy.core               import Mul, Add
+
+from sympde.calculus.core import BasicOperator
+from symla.kronecker      import Kron, LinearOperator, Vector
+
+#==============================================================================
+class Dot(BasicOperator):
+    """
+    """
+    def __new__(cls, *args, **options):
+        # (Try to) sympify args first
+
+        if options.pop('evaluate', True):
+            r = cls.eval(*args)
+        else:
+            r = None
+
+        if r is None:
+            return Basic.__new__(cls, *args, **options)
+        else:
+            return r
+
+    @classmethod
+    def eval(cls, *_args):
+
+        if not _args:
+            return
+
+        if not len(_args) == 2:
+            raise ValueError('Expecting one argument')
+
+        args = _args
+        args = list(args)
+        originals = args.copy()
+
+        # ...
+
+        # ... treate the case where there is and Mul node
+        mul_args = [arg for arg in args if isinstance(arg, Mul) and 
+                    not all([isinstance(i, LinearOperator) for i in arg.args])]
+        if mul_args:
+            arg = mul_args[0]
+            index = originals.index(arg)
+
+            linops = [i for i in arg.args if isinstance(i, (LinearOperator, Kron, Vector, Add))]
+            coeffs = [i for i in arg.args if i not in linops]
+
+            linop = reduce(mul, linops)
+            if len(coeffs) > 0:
+                coeff = reduce(mul, coeffs)
+            else:
+                coeff = 1
+
+            newargs = originals.copy()
+            newargs[index] = linop
+            return coeff*cls(*newargs)
+        # ...
+
+        add_args = [arg for arg in args if isinstance(arg, Add)]
+        if add_args:
+            arg = add_args[0]
+            index = originals.index(arg)
+
+            _args = []
+            for i in arg.args:
+                newargs = originals.copy()
+                newargs[index] = i
+                newexpr = cls(*newargs)
+                _args.append(newexpr)
+
+            return reduce(add, _args)
+
+        return cls(*args, evaluate=False)
+
+    def _hashable_content(self):
+        return tuple(sorted(self._args, key= lambda x:x.__hash__()))
+
+    def __hash__(self):
+        h = self._mhash
+        if h is None:
+            h = hash((type(self).__name__,)) 
+        return h
+
+    def _sympystr(self, printer):
+        sstr = printer.doprint
+        return f'Dot({self.args[0]}, {self.args[1]})'

symla/inverse.py

@@ -1,16 +1,9 @@
-from operator  import mul, add
+from operator  import mul
 from functools import reduce
 
-from sympy                    import Indexed, sympify, Symbol
 from sympy                    import Matrix as sp_Matrix
-from sympy                    import ImmutableDenseMatrix
-from sympy                    import cacheit
 from sympy.core               import Basic
-from sympy.core               import Add, Mul, Pow, Expr
-from sympy.core.containers    import Tuple
-from sympy.core.singleton     import S
-from sympy.core.decorators    import call_highest_priority
-from sympy.core.compatibility import is_sequence
+from sympy.core               import Mul
 
 from sympde.calculus.core import BasicOperator
 from symla.kronecker      import Kron, LinearOperator
@@ -23,6 +16,8 @@ def __new__(cls, *args, **options):
         # (Try to) sympify args first
 
         if options.pop('evaluate', True):
+            args = [arg.expand() if hasattr(arg, 'expand') else arg
+                    for arg in args]
             r = cls.eval(*args)
         else:
             r = None

symla/kronecker.py

@@ -1,16 +1,12 @@
 from operator  import mul, add
 from functools import reduce
 
-from sympy                    import Indexed, sympify, Symbol
+from sympy                    import Symbol
 from sympy                    import Matrix as sp_Matrix
 from sympy                    import ImmutableDenseMatrix
 from sympy                    import cacheit
 from sympy.core               import Basic
-from sympy.core               import Add, Mul, Pow, Expr
-from sympy.core.containers    import Tuple
-from sympy.core.singleton     import S
-from sympy.core.decorators    import call_highest_priority
-from sympy.core.compatibility import is_sequence
+from sympy.core               import Add, Mul, Expr
 
 from sympde.calculus.core import BasicOperator
 
@@ -33,11 +29,12 @@ def is_zero(x):
 class FiniteVectorSpace(Basic):
     """
     """
-    def __new__(cls, name, shape=None):
+    def __new__(cls, name, dimension=None, field=None):
 
         obj = Basic.__new__(cls)
         obj._name  = name
-        obj._shape = shape
+        obj._dimension = dimension
+        obj._field = field
 
         return obj
 
@@ -46,8 +43,17 @@ def name(self):
         return self._name
 
     @property
-    def shape(self):
-        return self._shape
+    def dimension(self):
+        return self._dimension
+
+    @property
+    def field(self):
+        return self._field
+
+    @property
+    def dtype(self):
+        # TODO assign a data type for each field
+        pass
 
     def _sympystr(self, printer):
         sstr = printer.doprint
@@ -57,12 +63,31 @@ def __mul__(self, other):
         raise NotImplementedError('TODO')
 
     def __hash__(self):
-        return hash((self.name, self.shape))
+        return hash((self.name, self._dimension, self._field))
 
 #==============================================================================
 # TODO improve
 class Vector(Symbol):
-    pass
+    def __new__(cls, name, space=None, **assumptions):
+        obj = Symbol.__new__(cls, name)
+        obj._space = space
+        obj._discretizable = False
+        if assumptions.get('values'):
+            obj._values = assumptions.get('values')
+            # TODO validate values type and shape
+            obj._discretizable = True
+
+        return obj
+
+    @property
+    def shape(self):
+        return None if self._space == None else self._space.dimension()
+
+    @property
+    def dtype(self):
+        return None if self._space == None else self._space.dtype()
+
+
 
 #==============================================================================
 class LinearOperator(Symbol):
@@ -111,8 +136,9 @@ class Kron(BasicOperator):
 
     def __new__(cls, *args, **options):
         # (Try to) sympify args first
-
         if options.pop('evaluate', True):
+            args = [arg.expand() if hasattr(arg, 'expand') else arg
+                    for arg in args]
             r = cls.eval(*args)
         else:
             r = None
@@ -169,12 +195,13 @@ def eval(cls, *_args):
         # ...
 
         # ... treate the case where there is and Mul node
-        mul_args = [arg for arg in args if isinstance(arg, Mul)]
+        mul_args = [arg for arg in args if isinstance(arg, Mul) and 
+                    not all([isinstance(i, LinearOperator) for i in arg.args])]
         if mul_args:
             arg = mul_args[0]
             index = originals.index(arg)
 
-            linops = [i for i in arg.args if isinstance(i, LinearOperator)]
+            linops = [i for i in arg.args if isinstance(i, (LinearOperator, Kron))]
             coeffs = [i for i in arg.args if i not in linops]
 
             linop = reduce(mul, linops)

symla/transpose.py

@@ -0,0 +1,69 @@
+from operator  import mul
+from functools import reduce
+
+from sympy.core               import Basic
+from sympy.core               import Add, Mul, Pow
+
+from sympde.calculus.core import BasicOperator
+from symla.kronecker      import Kron, LinearOperator
+
+#==============================================================================
+class Transpose(BasicOperator):
+    """
+    """
+    def __new__(cls, *args, **options):
+        # (Try to) sympify args first
+
+        if options.pop('evaluate', True):
+            args = [arg.expand() if hasattr(arg, 'expand') else arg
+                    for arg in args]
+            r = cls.eval(*args)
+        else:
+            r = None
+
+        if r is None:
+            return Basic.__new__(cls, *args, **options)
+        else:
+            return r
+
+    @classmethod
+    def eval(cls, *_args):
+
+        if not _args:
+            return
+        if not len(_args) == 1:
+            raise ValueError('Expecting one argument')
+
+        expr = _args[0]
+
+        if isinstance(expr, Transpose):
+            return _args[0].args[0]
+
+        if isinstance(expr, Kron):
+            args = [cls(a, evaluate=False) for a in expr.args]
+            return Kron(*args)
+
+        elif isinstance(expr, Pow):
+            return cls(expr.base) ** expr.exp
+
+        elif isinstance(expr, Add):
+            return reduce(Add, [cls(arg) for arg in expr.args])
+
+        elif isinstance(expr, Mul):
+            linops = [i for i in expr.args if isinstance(i, (LinearOperator, Kron, Transpose))]
+            coeffs = [i for i in expr.args if i not in linops]
+
+            linops = [cls(i) for i in linops]
+            linop = reduce(mul, linops[::-1])
+            if len(coeffs) > 0:
+                coeff = reduce(mul, coeffs)
+            else:
+                coeff = 1
+
+            return coeff * linop
+
+        return cls(expr, evaluate=False)
+
+    def _sympystr(self, printer):
+        sstr = printer.doprint
+        return 'Trans({arg})'.format(arg=sstr(self.args[0]))

test_dot.py

@@ -0,0 +1,26 @@
+from sympy import Symbol
+
+from symla.kronecker import FiniteVectorSpace
+from symla.kronecker import Matrix
+from symla.kronecker import Vector
+from symla.dot   import Dot
+from symla.transpose   import Transpose
+
+# ====================================================================
+def test_dot_1():
+    V = FiniteVectorSpace('V')
+    u, v, w = [Vector(c,V) for c in ['u','v', 'w']]
+    alpha, beta, gamma = [Symbol(i) for i in ['alpha', 'beta', 'gamma']]
+    # .....................................
+    assert(Dot(alpha*v,beta*u) == alpha*beta*Dot(v,u))
+    assert(Dot(alpha*(v+w),beta*u) == alpha*beta*(Dot(v,u) + Dot(w,u)))
+    assert(Dot(u,v) == Dot(v,u))
+    assert(Dot(u,v) - Dot(v,u) == 0)
+    assert(Dot(u,v) + Dot(v,u) == 2*Dot(v, u))
+    assert(Dot(u,u) +2*Dot(u,v) + Dot(v,v) == Dot(u+v,u+v))
+
+
+#######################################
+if __name__ == '__main__':
+
+    test_dot_1()

test_inverse.py

@@ -1,4 +1,5 @@
 from sympy import Symbol
+from sympy import Expr
 
 from symla.kronecker import FiniteVectorSpace
 from symla.kronecker import Matrix
@@ -26,23 +27,18 @@ def test_inverse_1():
     assert(Inverse(Kron(alpha*a,b)) == Kron(Inverse(a), Inverse(b))/alpha)
     assert(Inverse(Kron(alpha*a,beta*b)) == Kron(Inverse(a), Inverse(b))/(alpha*beta))
     assert(Inverse(Kron(alpha*a,beta*b,gamma*c)) == Kron(Inverse(a), Inverse(b), Inverse(c))/(alpha*beta*gamma))
+
+    # Inverse of Inverse
+    assert(Inverse(Inverse(a)) == a)
+    assert(Inverse(Inverse(Inverse(a))) == Inverse(a))
+    assert(Inverse(a*b*Inverse(c)*d) == Inverse(d)*c*Inverse(b)*Inverse(a))
+
+    # Inverse of Power
+    assert(Inverse(a**alpha) == Inverse(a) ** alpha)
+
 
-    expr = Inverse(a+b)
-    print(expr)
 
 #######################################
 if __name__ == '__main__':
 
-    test_inverse_1()
-
-    # .....................................
-#    # TODO not working
-#    expr = Kron(a+d*a, beta*b*c)
-#    print(expr)
-#    print('')
-
-#    # TODO not working
-#    expr = Kron(alpha*a,Kron(b, gamma*c) + Kron(beta*b,d))
-#    print(expr)
-#    print('')
-    # .....................................
+    test_inverse_1()