<commit_before>import warnings
from sympy import Basic, Symbol
from sympy.core import sympify

from sympy.core.basic import S, C
from sympy.polys import Poly, roots
from sympy.simplify import simplify

# from sympy.printing import StrPrinter /cyclic/

import random

class NonSquareMatrixException(Exception):
    pass

class ShapeError(ValueError):
    """Wrong matrix shape"""
    pass

class MatrixError(Exception):
    pass

def _dims_to_nm( dims ):
    """Converts dimensions tuple (or any object with length 1 or 2) or scalar
    in dims to matrix dimensions n and m."""

    try:
        l = len( dims )
    except TypeError:
        dims = (dims,)
        l = 1

    # This will work for nd-array too when they are added to sympy.
    try:
        for dim in dims:
            assert (dim > 0) and isinstance( dim, int )
    except AssertionError:
        raise ValueError("Matrix dimensions should positive integers!")

    if l == 2:
        n, m = dims
    elif l == 1:
        n, m = dims[0], dims[0]
    else:
        raise ValueError("Matrix dimensions should be a two-element tuple of ints or a single int!")

    return n, m

def _iszero(x):
    return x == 0

class DeferredVector(object):
    def __init__(self,name):
        self.name=name
    def __getitem__(self,i):
        component_name = '%s[%d]'%(self.name,i)
        return Symbol(component_name)
    def __str__(self):
        return StrPrinter.doprint(self)
    def __repr__(self):
        return StrPrinter.doprint(self)

class Matrix(object):

    # Added just for numpy compatibility
    # TODO: investigate about __array_priority__
    __array_priority__ = 10.0

    def __init__(self, *args):
        """
        Matrix can be constructed with values or a rule.

        >>> from sympy import *
        >>> Matrix( ((1,2+I), (3,4)) ) #doctest:+NORMALIZE_WHITESPACE
        [1, 2 + I]
        [3,     4]
        >>> Matrix(2, 2, lambda i,j: (i+1)*j ) #doctest:+NORMALIZE_WHITESPACE
        [0, 1]
        [0, 2]
        """
        if len(args) == 3 and callable(args[2]):
            operation = args[2]
            assert isinstance(args[0], int) and isinstance(args[1], int)
            self.lines = args[0]
            self.cols = args[1]
            self.mat = []
            for i in range(self.lines):
                for j in range(self.cols):
                    self.mat.append(sympify(operation(i, j)))
        elif len(args)==3 and isinstance(args[0],int) and \
                isinstance(args[1],int) and isinstance(args[2], (list, tuple)):
            self.lines=args[0]
            self.cols=args[1]
            mat = args[2]
            if len(mat) != self.lines*self.cols:
                raise MatrixError('List length should be equal to rows*columns')
            self.mat = map(lambda i: sympify(i), mat)
        elif len(args) == 1:
            mat = args[0]
            if isinstance(mat, Matrix):
                self.lines = mat.lines
                self.cols = mat.cols
                self.mat = mat[:]
                return
            elif hasattr(mat, "__array__"):
                # NumPy array or matrix or some other object that implements
                # __array__. So let's first use this method to get a
                # numpy.array() and then make a python list out of it.
                arr = mat.__array__()
                if len(arr.shape) == 2:
                    self.lines, self.cols = arr.shape[0], arr.shape[1]
                    self.mat = map(lambda i: sympify(i), arr.ravel())
                    return
                elif len(arr.shape) == 1:
                    self.lines, self.cols = 1, arr.shape[0]
                    self.mat = [0]*self.cols
                    for i in xrange(len(arr)):
                        self.mat[i] = sympify(arr[i])
                    return
                else:
                    raise NotImplementedError("Sympy supports just 1D and 2D matrices")
            elif not isinstance(mat, (list, tuple)):
                raise TypeError("Matrix constructor doesn't accept %s as input" % str(type(mat)))
            self.lines = len(mat)
            if len(mat) != 0:
                if not isinstance(mat[0], (list, tuple)):
                    self.cols = 1
                    self.mat = map(lambda i: sympify(i), mat)
                    return
                self.cols = len(mat[0])
            else:
                self.cols = 0
            self.mat = []
            for j in xrange(self.lines):
                assert len(mat[j])==self.cols
                for i in xrange(self.cols):
                    self.mat.append(sympify(mat[j][i]))
        elif len(args) == 0:
            # Empty Matrix
            self.lines = self.cols = 0
            self.mat = []
        else:
            # TODO: on 0.7.0 delete this and uncomment the last line
            mat = args
            if not isinstance(mat[0], (list, tuple)):
                # make each element a singleton
                mat = [ [element] for element in mat ]
            warnings.warn("Deprecated constructor, use brackets: Matrix(%s)" % str(mat))
            self.lines=len(mat)
            self.cols=len(mat[0])
            self.mat=[]
            for j in xrange(self.lines):
                assert len(mat[j])==self.cols
                for i in xrange(self.cols):
                    self.mat.append(sympify(mat[j][i]))
            #raise TypeError("Data type not understood")

    def key2ij(self,key):
        """Converts key=(4,6) to 4,6 and ensures the key is correct."""
        if not (isinstance(key,(list, tuple)) and len(key) == 2):
            raise TypeError("wrong syntax: a[%s]. Use a[i,j] or a[(i,j)]"
                    %repr(key))
        i,j=key
        if not (i>=0 and i<self.lines and j>=0 and j < self.cols):
            print self.lines, " ", self.cols
            raise IndexError("Index out of range: a[%s]"%repr(key))
        return i,j

    def transpose(self):
        """
        Matrix transposition.

        >>> from sympy import *
        >>> m=Matrix(((1,2+I),(3,4)))
        >>> m  #doctest: +NORMALIZE_WHITESPACE
        [1, 2 + I]
        [3,     4]
        >>> m.transpose() #doctest: +NORMALIZE_WHITESPACE
        [    1, 3]
        [2 + I, 4]
        >>> m.T == m.transpose()
        True
        """
        a = [0]*self.cols*self.lines
        for i in xrange(self.cols):
            a[i*self.lines:(i+1)*self.lines] = self.mat[i::self.cols]
        return Matrix(self.cols,self.lines,a)

    T = property(transpose,None,None,"Matrix transposition.")

    def conjugate(self):
        """By-element conjugation."""
        out = Matrix(self.lines,self.cols,
                lambda i,j: self[i,j].conjugate())
        return out

    C = property(conjugate,None,None,"By-element conjugation.")

    @property
    def H(self):
        """
        Hermite conjugation.

        >>> from sympy import *
        >>> m=Matrix(((1,2+I),(3,4)))
        >>> m  #doctest: +NORMALIZE_WHITESPACE
        [1, 2 + I]
        [3,     4]
        >>> m.H #doctest: +NORMALIZE_WHITESPACE
        [    1, 3]
        [2 - I, 4]
        """
        out = self.T.C
        return out

    @property
    def D(self):
        """Dirac conjugation."""
        from sympy.physics.matrices import mgamma
        out = self.H * mgamma(0)
        return out

    def __getitem__(self,key):
        """
        >>> from sympy import *
        >>> m=Matrix(((1,2+I),(3,4)))
        >>> m  #doctest: +NORMALIZE_WHITESPACE
        [1, 2 + I]
        [3,     4]
        >>> m[1,0]
        3
        >>> m.H[1,0]
        2 - I

        """
        if type(key) is tuple:
            i, j = key
            if type(i) is slice or type(j) is slice:
                return self.submatrix(key)

            else:
                # a2idx inlined
                try:
                    i = i.__int__()
                except AttributeError:
                    try:
                        i = i.__index__()
                    except AttributeError:
                        raise IndexError("Invalid index a[%r]" % (key,))

                # a2idx inlined
                try:
                    j = j.__int__()
                except AttributeError:
                    try:
                        j = j.__index__()
                    except AttributeError:
                        raise IndexError("Invalid index a[%r]" % (key,))


                if not (i>=0 and i<self.lines and j>=0 and j < self.cols):
                    raise IndexError("Index out of range: a[%s]" % (key,))
                else:
                    return self.mat[i*self.cols + j]


        else:
            # row-wise decomposition of matrix
            if type(key) is slice:
                return self.mat[key]
            else:
                k = a2idx(key)
                if k is not None:
                    return self.mat[k]
        raise IndexError("Invalid index: a[%s]" % repr(key))

    def __setitem__(self, key, value):
        """
        >>> from sympy import *
        >>> m=Matrix(((1,2+I),(3,4)))
        >>> m  #doctest: +NORMALIZE_WHITESPACE
        [1, 2 + I]
        [3,     4]
        >>> m[1,0]=9
        >>> m  #doctest: +NORMALIZE_WHITESPACE
        [1, 2 + I]
        [9,     4]

        """
        if type(key) is tuple:
            i, j = key
            if type(i) is slice or type(j) is slice:
                if isinstance(value, Matrix):
                    self.copyin_matrix(key, value)
                    return
                if isinstance(value, (list, tuple)):
                    self.copyin_list(key, value)
                    return
            else:
                # a2idx inlined
                try:
                    i = i.__int__()
                except AttributeError:
                    try:
                        i = i.__index__()
                    except AttributeError:
                        raise IndexError("Invalid index a[%r]" % (key,))

                # a2idx inlined
                try:
                    j = j.__int__()
                except AttributeError:
                    try:
                        j = j.__index__()
                    except AttributeError:
                        raise IndexError("Invalid index a[%r]" % (key,))


                if not (i>=0 and i<self.lines and j>=0 and j < self.cols):
                    raise IndexError("Index out of range: a[%s]" % (key,))
                else:
                    self.mat[i*self.cols + j] = sympify(value)
                    return

        else:
            # row-wise decomposition of matrix
            if type(key) is slice:
                raise IndexError("Vector slices not implemented yet.")
            else:
                k = a2idx(key)
                if k is not None:
                    self.mat[k] = sympify(value)
                    return
        raise IndexError("Invalid index: a[%s]"%repr(key))

    def __array__(self):
        return matrix2numpy(self)

    def tolist(self):
        """
        Return the Matrix converted in a python list.

        >>> from sympy import *
        >>> m=Matrix(3, 3, range(9))
        >>> m
        [0, 1, 2]
        [3, 4, 5]
        [6, 7, 8]
        >>> m.tolist()
        [[0, 1, 2], [3, 4, 5], [6, 7, 8]]
        """
        ret = [0]*self.lines
        for i in xrange(self.lines):
            ret[i] = self.mat[i*self.cols:(i+1)*self.cols]
        return ret

    def copyin_matrix(self, key, value):
        rlo, rhi = self.slice2bounds(key[0], self.lines)
        clo, chi = self.slice2bounds(key[1], self.cols)
        assert value.lines == rhi - rlo and value.cols == chi - clo
        for i in range(value.lines):
            for j in range(value.cols):
                self[i+rlo, j+clo] = sympify(value[i,j])

    def copyin_list(self, key, value):
        assert isinstance(value, (list, tuple))
        self.copyin_matrix(key, Matrix(value))

    def hash(self):
        """Compute a hash every time, because the matrix elements
        could change."""
        return hash(self.__str__() )

    @property
    def shape(self):
        return (self.lines, self.cols)

    def __rmul__(self,a):
        if hasattr(a, "__array__"):
            return matrix_multiply(a,self)
        out = Matrix(self.lines,self.cols,map(lambda i: a*i,self.mat))
        return out

    def expand(self):
        out = Matrix(self.lines,self.cols,map(lambda i: i.expand(), self.mat))
        return out

    def combine(self):
        out = Matrix(self.lines,self.cols,map(lambda i: i.combine(),self.mat))
        return out

    def subs(self, *args):
        out = Matrix(self.lines,self.cols,map(lambda i: i.subs(*args),self.mat))
        return out

    def __sub__(self,a):
        return self + (-a)

    def __mul__(self,a):
        if hasattr(a, "__array__"):
            return matrix_multiply(self,a)
        out = Matrix(self.lines,self.cols,map(lambda i: i*a,self.mat))
        return out

    def __pow__(self, num):
        if not self.is_square:
            raise NonSquareMatrixException()
        if isinstance(num, int) or isinstance(num, Integer):
            n = int(num)
            if n < 0:
                return self.inv() ** -n   # A**-2 = (A**-1)**2
            a = eye(self.cols)
            while n:
                if n % 2:
                    a = a * self
                    n -= 1
                self = self * self
                n = n // 2
            return a
        raise NotImplementedError('Can only rise to the power of an integer for now')

    def __add__(self,a):
        return matrix_add(self,a)

    def __radd__(self,a):
        return matrix_add(a,self)

    def __div__(self,a):
        return self * (S.One/a)

    def __truediv__(self,a):
        return self.__div__(a)

    def multiply(self,b):
        """Returns self*b """
        return matrix_multiply(self,b)

    def add(self,b):
        """Return self+b """
        return matrix_add(self,b)

    def __neg__(self):
        return -1*self

    def __eq__(self, a):
        if not isinstance(a, (Matrix, Basic)):
            a = sympify(a)
        if isinstance(a, Matrix):
            return self.hash() == a.hash()
        else:
            return False

    def __ne__(self,a):
        if not isinstance(a, (Matrix, Basic)):
            a = sympify(a)
        if isinstance(a, Matrix):
            return self.hash() != a.hash()
        else:
            return True

    def _format_str(self, strfunc, rowsep='\n'):
        # Build table of string representations of the elements
        res = []
        # Track per-column max lengths for pretty alignment
        maxlen = [0] * self.cols
        for i in range(self.lines):
            res.append([])
            for j in range(self.cols):
                string = strfunc(self[i,j])
                res[-1].append(string)
                maxlen[j] = max(len(string), maxlen[j])
        # Patch strings together
        for i, row in enumerate(res):
            for j, elem in enumerate(row):
                # Pad each element up to maxlen so the columns line up
                row[j] = elem.rjust(maxlen[j])
            res[i] = "[" + ", ".join(row) + "]"
        return rowsep.join(res)

    def __str__(self):
        return StrPrinter.doprint(self)

    def __repr__(self):
        return StrPrinter.doprint(self)

    def inv(self, method="GE", iszerofunc=_iszero):
        """
        Calculates the matrix inverse.

        According to the "method" parameter, it calls the appropriate method:

          GE .... inverse_GE()
          LU .... inverse_LU()
          ADJ ... inverse_ADJ()

        Note, the GE and LU methods may require the matrix to be simplified
        before it is inverted in order to properly detect zeros during
        pivoting. In difficult cases a custom zero detection function can
        be provided by setting the iszerosfunc argument to a function that
        should return True if its argument is zero.

        """
        assert self.cols==self.lines
        if method == "GE":
            return self.inverse_GE(iszerofunc=iszerofunc)
        elif method == "LU":
            return self.inverse_LU(iszerofunc=iszerofunc)
        elif method == "ADJ":
            return self.inverse_ADJ()
        else:
            raise Exception("Inversion method unrecognized")


    def __mathml__(self):
        mml = ""
        for i in range(self.lines):
            mml += "<matrixrow>"
            for j in range(self.cols):
                mml += self[i,j].__mathml__()
            mml += "</matrixrow>"
        return "<matrix>" + mml + "</matrix>"

    def row(self, i, f):
        """Elementary row operation using functor"""
        for j in range(0, self.cols):
            self[i, j] = f(self[i, j], j)

    def col(self, j, f):
        """Elementary column operation using functor"""
        for i in range(0, self.lines):
            self[i, j] = f(self[i, j], i)

    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]

    def col_swap(self, i, j):
        for k in range(0, self.lines):
            self[k, i], self[k, j] = self[k, j], self[k, i]

    def row_del(self, i):
        self.mat = self.mat[:i*self.cols] + self.mat[(i+1)*self.cols:]
        self.lines -= 1

    def col_del(self, i):
        """
        >>> import sympy
        >>> M = sympy.matrices.eye(3)
        >>> M.col_del(1)
        >>> M   #doctest: +NORMALIZE_WHITESPACE
        [1, 0]
        [0, 0]
        [0, 1]
        """
        for j in range(self.lines-1, -1, -1):
            del self.mat[i+j*self.cols]
        self.cols -= 1

    def row_join(self, rhs):
        """
        Concatenates two matrices along self's last and rhs's first column

        >>> from sympy import *
        >>> M = Matrix(3,3,lambda i,j: i+j)
        >>> V = Matrix(3,1,lambda i,j: 3+i+j)
        >>> M.row_join(V)
        [0, 1, 2, 3]
        [1, 2, 3, 4]
        [2, 3, 4, 5]
        """
        assert self.lines == rhs.lines
        newmat = self.zeros((self.lines, self.cols + rhs.cols))
        newmat[:,:self.cols] = self[:,:]
        newmat[:,self.cols:] = rhs
        return newmat

    def col_join(self, bott):
        """
        Concatenates two matrices along self's last and bott's first row

        >>> from sympy import *
        >>> M = Matrix(3,3,lambda i,j: i+j)
        >>> V = Matrix(1,3,lambda i,j: 3+i+j)
        >>> M.col_join(V)
        [0, 1, 2]
        [1, 2, 3]
        [2, 3, 4]
        [3, 4, 5]
        """
        assert self.cols == bott.cols
        newmat = self.zeros((self.lines+bott.lines, self.cols))
        newmat[:self.lines,:] = self[:,:]
        newmat[self.lines:,:] = bott
        return newmat

    def row_insert(self, pos, mti):
        """
        >>> from sympy import *
        >>> M = Matrix(3,3,lambda i,j: i+j)
        >>> M
        [0, 1, 2]
        [1, 2, 3]
        [2, 3, 4]
        >>> V = zeros((1, 3))
        >>> V
        [0, 0, 0]
        >>> M.row_insert(1,V)
        [0, 1, 2]
        [0, 0, 0]
        [1, 2, 3]
        [2, 3, 4]
        """
        if pos is 0:
            return mti.col_join(self)
        assert self.cols == mti.cols
        newmat = self.zeros((self.lines + mti.lines, self.cols))
        newmat[:pos,:] = self[:pos,:]
        newmat[pos:pos+mti.lines,:] = mti[:,:]
        newmat[pos+mti.lines:,:] = self[pos:,:]
        return newmat

    def col_insert(self, pos, mti):
        """
        >>> from sympy import *
        >>> M = Matrix(3,3,lambda i,j: i+j)
        >>> M
        [0, 1, 2]
        [1, 2, 3]
        [2, 3, 4]
        >>> V = zeros((3, 1))
        >>> V
        [0]
        [0]
        [0]
        >>> M.col_insert(1,V)
        [0, 0, 1, 2]
        [1, 0, 2, 3]
        [2, 0, 3, 4]
        """
        if pos is 0:
            return mti.row_join(self)
        assert self.lines == mti.lines
        newmat = self.zeros((self.lines, self.cols + mti.cols))
        newmat[:,:pos] = self[:,:pos]
        newmat[:,pos:pos+mti.cols] = mti[:,:]
        newmat[:,pos+mti.cols:] = self[:,pos:]
        return newmat

    def trace(self):
        assert self.cols == self.lines
        trace = 0
        for i in range(self.cols):
            trace += self[i,i]
        return trace

    def submatrix(self, keys):
        """
        >>> from sympy import *
        >>> m = Matrix(4,4,lambda i,j: i+j)
        >>> m   #doctest: +NORMALIZE_WHITESPACE
        [0, 1, 2, 3]
        [1, 2, 3, 4]
        [2, 3, 4, 5]
        [3, 4, 5, 6]
        >>> m[0:1, 1]   #doctest: +NORMALIZE_WHITESPACE
        [1]
        >>> m[0:2, 0:1] #doctest: +NORMALIZE_WHITESPACE
        [0]
        [1]
        >>> m[2:4, 2:4] #doctest: +NORMALIZE_WHITESPACE
        [4, 5]
        [5, 6]
        """
        assert isinstance(keys[0], slice) or isinstance(keys[1], slice)
        rlo, rhi = self.slice2bounds(keys[0], self.lines)
        clo, chi = self.slice2bounds(keys[1], self.cols)
        if not ( 0<=rlo<=rhi and 0<=clo<=chi ):
            raise IndexError("Slice indices out of range: a[%s]"%repr(keys))
        outLines, outCols = rhi-rlo, chi-clo
        outMat = [0]*outLines*outCols
        for i in xrange(outLines):
            outMat[i*outCols:(i+1)*outCols] = self.mat[(i+rlo)*self.cols+clo:(i+rlo)*self.cols+chi]
        return Matrix(outLines,outCols,outMat)

    def slice2bounds(self, key, defmax):
        """
            Takes slice or number and returns (min,max) for iteration
            Takes a default maxval to deal with the slice ':' which is (none, none)
        """
        if isinstance(key, slice):
            lo, hi = 0, defmax
            if key.start != None:
                if key.start >= 0:
                    lo = key.start
                else:
                    lo = defmax+key.start
            if key.stop != None:
                if key.stop >= 0:
                    hi = key.stop
                else:
                    hi = defmax+key.stop
            return lo, hi
        elif isinstance(key, int):
            if key >= 0:
                return key, key+1
            else:
                return defmax+key, defmax+key+1
        else:
            raise IndexError("Improper index type")

    def applyfunc(self, f):
        """
        >>> from sympy import *
        >>> m = Matrix(2,2,lambda i,j: i*2+j)
        >>> m   #doctest: +NORMALIZE_WHITESPACE
        [0, 1]
        [2, 3]
        >>> m.applyfunc(lambda i: 2*i)  #doctest: +NORMALIZE_WHITESPACE
        [0, 2]
        [4, 6]
        """
        assert callable(f)
        out = Matrix(self.lines,self.cols,map(f,self.mat))
        return out

    def evalf(self, prec=None, **options):
        if prec is None:
            return self.applyfunc(lambda i: i.evalf(**options))
        else:
            return self.applyfunc(lambda i: i.evalf(prec, **options))

    def reshape(self, _rows, _cols):
        """
        >>> from sympy import *
        >>> m = Matrix(2,3,lambda i,j: 1)
        >>> m   #doctest: +NORMALIZE_WHITESPACE
        [1, 1, 1]
        [1, 1, 1]
        >>> m.reshape(1,6)  #doctest: +NORMALIZE_WHITESPACE
        [1, 1, 1, 1, 1, 1]
        >>> m.reshape(3,2)  #doctest: +NORMALIZE_WHITESPACE
        [1, 1]
        [1, 1]
        [1, 1]
        """
        if self.lines*self.cols != _rows*_cols:
            print "Invalid reshape parameters %d %d" % (_rows, _cols)
        return Matrix(_rows, _cols, lambda i,j: self.mat[i*_cols + j])

    def print_nonzero (self, symb="X"):
        """
        Shows location of non-zero entries for fast shape lookup
        >>> from sympy import *
        >>> m = Matrix(2,3,lambda i,j: i*3+j)
        >>> m           #doctest: +NORMALIZE_WHITESPACE
        [0, 1, 2]
        [3, 4, 5]
        >>> m.print_nonzero()   #doctest: +NORMALIZE_WHITESPACE
        [ XX]
        [XXX]
        >>> m = matrices.eye(4)
        >>> m.print_nonzero("x")    #doctest: +NORMALIZE_WHITESPACE
        [x   ]
        [ x  ]
        [  x ]
        [   x]
        """
        s="";
        for i in range(self.lines):
            s+="["
            for j in range(self.cols):
                if self[i,j] == 0:
                    s+=" "
                else:
                    s+= symb+""
            s+="]\n"
        print s

    def LUsolve(self, rhs, iszerofunc=_iszero):
        """
        Solve the linear system Ax = b.
        self is the coefficient matrix A and rhs is the right side b.
        """
        assert rhs.lines == self.lines
        A, perm = self.LUdecomposition_Simple(iszerofunc=_iszero)
        n = self.lines
        b = rhs.permuteFwd(perm)
        # forward substitution, all diag entries are scaled to 1
        for i in range(n):
            for j in range(i):
                b.row(i, lambda x,k: x - b[j,k]*A[i,j])
        # backward substitution
        for i in range(n-1,-1,-1):
            for j in range(i+1, n):
                b.row(i, lambda x,k: x - b[j,k]*A[i,j])
            b.row(i, lambda x,k: x / A[i,i])
        return b

    def LUdecomposition(self, iszerofunc=_iszero):
        """
        Returns the decompositon LU and the row swaps p.
        """
        combined, p = self.LUdecomposition_Simple(iszerofunc=_iszero)
        L = self.zeros(self.lines)
        U = self.zeros(self.lines)
        for i in range(self.lines):
            for j in range(self.lines):
                if i > j:
                    L[i,j] = combined[i,j]
                else:
                    if i == j:
                        L[i,i] = 1
                    U[i,j] = combined[i,j]
        return L, U, p

    def LUdecomposition_Simple(self, iszerofunc=_iszero):
        """
        Returns A compused of L,U (L's diag entries are 1) and
        p which is the list of the row swaps (in order).
        """
        assert self.lines == self.cols
        n = self.lines
        A = self[:,:]
        p = []
        # factorization
        for j in range(n):
            for i in range(j):
                for k in range(i):
                    A[i,j] = A[i,j] - A[i,k]*A[k,j]
            pivot = -1
            for i in range(j,n):
                for k in range(j):
                    A[i,j] = A[i,j] - A[i,k]*A[k,j]
                # find the first non-zero pivot, includes any expression
                if pivot == -1 and not iszerofunc(A[i,j]):
                    pivot = i
            if pivot < 0:
                raise "Error: non-invertible matrix passed to LUdecomposition_Simple()"
            if pivot != j: # row must be swapped
                A.row_swap(pivot,j)
                p.append([pivot,j])
            assert not iszerofunc(A[j,j])
            scale = 1 / A[j,j]
            for i in range(j+1,n):
                A[i,j] = A[i,j] * scale
        return A, p


    def LUdecompositionFF(self):
        """
        Returns 4 matrices P, L, D, U such that PA = L D**-1 U.

        From the paper "fraction-free matrix factors..." by Zhou and Jeffrey
        """
        n, m = self.lines, self.cols
        U, L, P = self[:,:], eye(n), eye(n)
        DD = zeros(n) # store it smarter since it's just diagonal
        oldpivot = 1

        for k in range(n-1):
            if U[k,k] == 0:
                kpivot = k+1
                Notfound = True
                while kpivot < n and Notfound:
                    if U[kpivot, k] != 0:
                        Notfound = False
                    else:
                        kpivot = kpivot + 1
                if kpivot == n+1:
                    raise "Matrix is not full rank"
                else:
                    swap = U[k, k:]
                    U[k,k:] = U[kpivot,k:]
                    U[kpivot, k:] = swap
                    swap = P[k, k:]
                    P[k, k:] = P[kpivot, k:]
                    P[kpivot, k:] = swap
            assert U[k, k] != 0
            L[k,k] = U[k,k]
            DD[k,k] = oldpivot * U[k,k]
            assert DD[k,k] != 0
            Ukk = U[k,k]
            for i in range(k+1, n):
                L[i,k] = U[i,k]
                Uik = U[i,k]
                for j in range(k+1, m):
                    U[i,j] = (Ukk * U[i,j] - U[k,j]*Uik) / oldpivot
                U[i,k] = 0
            oldpivot = U[k,k]
        DD[n-1,n-1] = oldpivot
        return P, L, DD, U

    def cofactorMatrix(self, method="berkowitz"):
        out = Matrix(self.lines, self.cols, lambda i,j:
                self.cofactor(i, j, method))
        return out

    def minorEntry(self, i, j, method="berkowitz"):
        assert 0 <= i < self.lines and 0 <= j < self.cols
        return self.minorMatrix(i,j).det(method)

    def minorMatrix(self, i, j):
        assert 0 <= i < self.lines and 0 <= j < self.cols
        return self.delRowCol(i,j)

    def cofactor(self, i, j, method="berkowitz"):
        if (i+j) % 2 == 0:
            return self.minorEntry(i, j, method)
        else:
            return -1 * self.minorEntry(i, j, method)

    def jacobian(self, varlist):
        """
        Calculates the Jacobian matrix (derivative of a vectorial function).

        self is a vector of expression representing functions f_i(x_1, ...,
        x_n).  varlist is the set of x_i's in order.
        """
        assert self.lines == 1
        m = self.cols
        if isinstance(varlist, Matrix):
            assert varlist.lines == 1
            n = varlist.cols
        elif isinstance(varlist, (list, tuple)):
            n = len(varlist)
        assert n > 0 # need to diff by something
        J = self.zeros((m, n)) # maintain subclass type
        for i in range(m):
            if isinstance(self[i], (float, int)):
                continue    # constant function, jacobian row is zero
            try:
                tmp = self[i].diff(varlist[0])   # check differentiability
                J[i,0] = tmp
            except AttributeError:
                raise "Function %d is not differentiable" % i
            for j in range(1,n):
                J[i,j] = self[i].diff(varlist[j])
        return J
<commit_msg>raising string exceptions is deprecated<commit_after>import warnings
from sympy import Basic, Symbol
from sympy.core import sympify

from sympy.core.basic import S, C
from sympy.polys import Poly, roots
from sympy.simplify import simplify

# from sympy.printing import StrPrinter /cyclic/

import random

class NonSquareMatrixException(Exception):
    pass

class ShapeError(ValueError):
    """Wrong matrix shape"""
    pass

class MatrixError(Exception):
    pass

def _dims_to_nm( dims ):
    """Converts dimensions tuple (or any object with length 1 or 2) or scalar
    in dims to matrix dimensions n and m."""

    try:
        l = len( dims )
    except TypeError:
        dims = (dims,)
        l = 1

    # This will work for nd-array too when they are added to sympy.
    try:
        for dim in dims:
            assert (dim > 0) and isinstance( dim, int )
    except AssertionError:
        raise ValueError("Matrix dimensions should positive integers!")

    if l == 2:
        n, m = dims
    elif l == 1:
        n, m = dims[0], dims[0]
    else:
        raise ValueError("Matrix dimensions should be a two-element tuple of ints or a single int!")

    return n, m

def _iszero(x):
    return x == 0

class DeferredVector(object):
    def __init__(self,name):
        self.name=name
    def __getitem__(self,i):
        component_name = '%s[%d]'%(self.name,i)
        return Symbol(component_name)
    def __str__(self):
        return StrPrinter.doprint(self)
    def __repr__(self):
        return StrPrinter.doprint(self)

class Matrix(object):

    # Added just for numpy compatibility
    # TODO: investigate about __array_priority__
    __array_priority__ = 10.0

    def __init__(self, *args):
        """
        Matrix can be constructed with values or a rule.

        >>> from sympy import *
        >>> Matrix( ((1,2+I), (3,4)) ) #doctest:+NORMALIZE_WHITESPACE
        [1, 2 + I]
        [3,     4]
        >>> Matrix(2, 2, lambda i,j: (i+1)*j ) #doctest:+NORMALIZE_WHITESPACE
        [0, 1]
        [0, 2]
        """
        if len(args) == 3 and callable(args[2]):
            operation = args[2]
            assert isinstance(args[0], int) and isinstance(args[1], int)
            self.lines = args[0]
            self.cols = args[1]
            self.mat = []
            for i in range(self.lines):
                for j in range(self.cols):
                    self.mat.append(sympify(operation(i, j)))
        elif len(args)==3 and isinstance(args[0],int) and \
                isinstance(args[1],int) and isinstance(args[2], (list, tuple)):
            self.lines=args[0]
            self.cols=args[1]
            mat = args[2]
            if len(mat) != self.lines*self.cols:
                raise MatrixError('List length should be equal to rows*columns')
            self.mat = map(lambda i: sympify(i), mat)
        elif len(args) == 1:
            mat = args[0]
            if isinstance(mat, Matrix):
                self.lines = mat.lines
                self.cols = mat.cols
                self.mat = mat[:]
                return
            elif hasattr(mat, "__array__"):
                # NumPy array or matrix or some other object that implements
                # __array__. So let's first use this method to get a
                # numpy.array() and then make a python list out of it.
                arr = mat.__array__()
                if len(arr.shape) == 2:
                    self.lines, self.cols = arr.shape[0], arr.shape[1]
                    self.mat = map(lambda i: sympify(i), arr.ravel())
                    return
                elif len(arr.shape) == 1:
                    self.lines, self.cols = 1, arr.shape[0]
                    self.mat = [0]*self.cols
                    for i in xrange(len(arr)):
                        self.mat[i] = sympify(arr[i])
                    return
                else:
                    raise NotImplementedError("Sympy supports just 1D and 2D matrices")
            elif not isinstance(mat, (list, tuple)):
                raise TypeError("Matrix constructor doesn't accept %s as input" % str(type(mat)))
            self.lines = len(mat)
            if len(mat) != 0:
                if not isinstance(mat[0], (list, tuple)):
                    self.cols = 1
                    self.mat = map(lambda i: sympify(i), mat)
                    return
                self.cols = len(mat[0])
            else:
                self.cols = 0
            self.mat = []
            for j in xrange(self.lines):
                assert len(mat[j])==self.cols
                for i in xrange(self.cols):
                    self.mat.append(sympify(mat[j][i]))
        elif len(args) == 0:
            # Empty Matrix
            self.lines = self.cols = 0
            self.mat = []
        else:
            # TODO: on 0.7.0 delete this and uncomment the last line
            mat = args
            if not isinstance(mat[0], (list, tuple)):
                # make each element a singleton
                mat = [ [element] for element in mat ]
            warnings.warn("Deprecated constructor, use brackets: Matrix(%s)" % str(mat))
            self.lines=len(mat)
            self.cols=len(mat[0])
            self.mat=[]
            for j in xrange(self.lines):
                assert len(mat[j])==self.cols
                for i in xrange(self.cols):
                    self.mat.append(sympify(mat[j][i]))
            #raise TypeError("Data type not understood")

    def key2ij(self,key):
        """Converts key=(4,6) to 4,6 and ensures the key is correct."""
        if not (isinstance(key,(list, tuple)) and len(key) == 2):
            raise TypeError("wrong syntax: a[%s]. Use a[i,j] or a[(i,j)]"
                    %repr(key))
        i,j=key
        if not (i>=0 and i<self.lines and j>=0 and j < self.cols):
            print self.lines, " ", self.cols
            raise IndexError("Index out of range: a[%s]"%repr(key))
        return i,j

    def transpose(self):
        """
        Matrix transposition.

        >>> from sympy import *
        >>> m=Matrix(((1,2+I),(3,4)))
        >>> m  #doctest: +NORMALIZE_WHITESPACE
        [1, 2 + I]
        [3,     4]
        >>> m.transpose() #doctest: +NORMALIZE_WHITESPACE
        [    1, 3]
        [2 + I, 4]
        >>> m.T == m.transpose()
        True
        """
        a = [0]*self.cols*self.lines
        for i in xrange(self.cols):
            a[i*self.lines:(i+1)*self.lines] = self.mat[i::self.cols]
        return Matrix(self.cols,self.lines,a)

    T = property(transpose,None,None,"Matrix transposition.")

    def conjugate(self):
        """By-element conjugation."""
        out = Matrix(self.lines,self.cols,
                lambda i,j: self[i,j].conjugate())
        return out

    C = property(conjugate,None,None,"By-element conjugation.")

    @property
    def H(self):
        """
        Hermite conjugation.

        >>> from sympy import *
        >>> m=Matrix(((1,2+I),(3,4)))
        >>> m  #doctest: +NORMALIZE_WHITESPACE
        [1, 2 + I]
        [3,     4]
        >>> m.H #doctest: +NORMALIZE_WHITESPACE
        [    1, 3]
        [2 - I, 4]
        """
        out = self.T.C
        return out

    @property
    def D(self):
        """Dirac conjugation."""
        from sympy.physics.matrices import mgamma
        out = self.H * mgamma(0)
        return out

    def __getitem__(self,key):
        """
        >>> from sympy import *
        >>> m=Matrix(((1,2+I),(3,4)))
        >>> m  #doctest: +NORMALIZE_WHITESPACE
        [1, 2 + I]
        [3,     4]
        >>> m[1,0]
        3
        >>> m.H[1,0]
        2 - I

        """
        if type(key) is tuple:
            i, j = key
            if type(i) is slice or type(j) is slice:
                return self.submatrix(key)

            else:
                # a2idx inlined
                try:
                    i = i.__int__()
                except AttributeError:
                    try:
                        i = i.__index__()
                    except AttributeError:
                        raise IndexError("Invalid index a[%r]" % (key,))

                # a2idx inlined
                try:
                    j = j.__int__()
                except AttributeError:
                    try:
                        j = j.__index__()
                    except AttributeError:
                        raise IndexError("Invalid index a[%r]" % (key,))


                if not (i>=0 and i<self.lines and j>=0 and j < self.cols):
                    raise IndexError("Index out of range: a[%s]" % (key,))
                else:
                    return self.mat[i*self.cols + j]


        else:
            # row-wise decomposition of matrix
            if type(key) is slice:
                return self.mat[key]
            else:
                k = a2idx(key)
                if k is not None:
                    return self.mat[k]
        raise IndexError("Invalid index: a[%s]" % repr(key))

    def __setitem__(self, key, value):
        """
        >>> from sympy import *
        >>> m=Matrix(((1,2+I),(3,4)))
        >>> m  #doctest: +NORMALIZE_WHITESPACE
        [1, 2 + I]
        [3,     4]
        >>> m[1,0]=9
        >>> m  #doctest: +NORMALIZE_WHITESPACE
        [1, 2 + I]
        [9,     4]

        """
        if type(key) is tuple:
            i, j = key
            if type(i) is slice or type(j) is slice:
                if isinstance(value, Matrix):
                    self.copyin_matrix(key, value)
                    return
                if isinstance(value, (list, tuple)):
                    self.copyin_list(key, value)
                    return
            else:
                # a2idx inlined
                try:
                    i = i.__int__()
                except AttributeError:
                    try:
                        i = i.__index__()
                    except AttributeError:
                        raise IndexError("Invalid index a[%r]" % (key,))

                # a2idx inlined
                try:
                    j = j.__int__()
                except AttributeError:
                    try:
                        j = j.__index__()
                    except AttributeError:
                        raise IndexError("Invalid index a[%r]" % (key,))


                if not (i>=0 and i<self.lines and j>=0 and j < self.cols):
                    raise IndexError("Index out of range: a[%s]" % (key,))
                else:
                    self.mat[i*self.cols + j] = sympify(value)
                    return

        else:
            # row-wise decomposition of matrix
            if type(key) is slice:
                raise IndexError("Vector slices not implemented yet.")
            else:
                k = a2idx(key)
                if k is not None:
                    self.mat[k] = sympify(value)
                    return
        raise IndexError("Invalid index: a[%s]"%repr(key))

    def __array__(self):
        return matrix2numpy(self)

    def tolist(self):
        """
        Return the Matrix converted in a python list.

        >>> from sympy import *
        >>> m=Matrix(3, 3, range(9))
        >>> m
        [0, 1, 2]
        [3, 4, 5]
        [6, 7, 8]
        >>> m.tolist()
        [[0, 1, 2], [3, 4, 5], [6, 7, 8]]
        """
        ret = [0]*self.lines
        for i in xrange(self.lines):
            ret[i] = self.mat[i*self.cols:(i+1)*self.cols]
        return ret

    def copyin_matrix(self, key, value):
        rlo, rhi = self.slice2bounds(key[0], self.lines)
        clo, chi = self.slice2bounds(key[1], self.cols)
        assert value.lines == rhi - rlo and value.cols == chi - clo
        for i in range(value.lines):
            for j in range(value.cols):
                self[i+rlo, j+clo] = sympify(value[i,j])

    def copyin_list(self, key, value):
        assert isinstance(value, (list, tuple))
        self.copyin_matrix(key, Matrix(value))

    def hash(self):
        """Compute a hash every time, because the matrix elements
        could change."""
        return hash(self.__str__() )

    @property
    def shape(self):
        return (self.lines, self.cols)

    def __rmul__(self,a):
        if hasattr(a, "__array__"):
            return matrix_multiply(a,self)
        out = Matrix(self.lines,self.cols,map(lambda i: a*i,self.mat))
        return out

    def expand(self):
        out = Matrix(self.lines,self.cols,map(lambda i: i.expand(), self.mat))
        return out

    def combine(self):
        out = Matrix(self.lines,self.cols,map(lambda i: i.combine(),self.mat))
        return out

    def subs(self, *args):
        out = Matrix(self.lines,self.cols,map(lambda i: i.subs(*args),self.mat))
        return out

    def __sub__(self,a):
        return self + (-a)

    def __mul__(self,a):
        if hasattr(a, "__array__"):
            return matrix_multiply(self,a)
        out = Matrix(self.lines,self.cols,map(lambda i: i*a,self.mat))
        return out

    def __pow__(self, num):
        if not self.is_square:
            raise NonSquareMatrixException()
        if isinstance(num, int) or isinstance(num, Integer):
            n = int(num)
            if n < 0:
                return self.inv() ** -n   # A**-2 = (A**-1)**2
            a = eye(self.cols)
            while n:
                if n % 2:
                    a = a * self
                    n -= 1
                self = self * self
                n = n // 2
            return a
        raise NotImplementedError('Can only rise to the power of an integer for now')

    def __add__(self,a):
        return matrix_add(self,a)

    def __radd__(self,a):
        return matrix_add(a,self)

    def __div__(self,a):
        return self * (S.One/a)

    def __truediv__(self,a):
        return self.__div__(a)

    def multiply(self,b):
        """Returns self*b """
        return matrix_multiply(self,b)

    def add(self,b):
        """Return self+b """
        return matrix_add(self,b)

    def __neg__(self):
        return -1*self

    def __eq__(self, a):
        if not isinstance(a, (Matrix, Basic)):
            a = sympify(a)
        if isinstance(a, Matrix):
            return self.hash() == a.hash()
        else:
            return False

    def __ne__(self,a):
        if not isinstance(a, (Matrix, Basic)):
            a = sympify(a)
        if isinstance(a, Matrix):
            return self.hash() != a.hash()
        else:
            return True

    def _format_str(self, strfunc, rowsep='\n'):
        # Build table of string representations of the elements
        res = []
        # Track per-column max lengths for pretty alignment
        maxlen = [0] * self.cols
        for i in range(self.lines):
            res.append([])
            for j in range(self.cols):
                string = strfunc(self[i,j])
                res[-1].append(string)
                maxlen[j] = max(len(string), maxlen[j])
        # Patch strings together
        for i, row in enumerate(res):
            for j, elem in enumerate(row):
                # Pad each element up to maxlen so the columns line up
                row[j] = elem.rjust(maxlen[j])
            res[i] = "[" + ", ".join(row) + "]"
        return rowsep.join(res)

    def __str__(self):
        return StrPrinter.doprint(self)

    def __repr__(self):
        return StrPrinter.doprint(self)

    def inv(self, method="GE", iszerofunc=_iszero):
        """
        Calculates the matrix inverse.

        According to the "method" parameter, it calls the appropriate method:

          GE .... inverse_GE()
          LU .... inverse_LU()
          ADJ ... inverse_ADJ()

        Note, the GE and LU methods may require the matrix to be simplified
        before it is inverted in order to properly detect zeros during
        pivoting. In difficult cases a custom zero detection function can
        be provided by setting the iszerosfunc argument to a function that
        should return True if its argument is zero.

        """
        assert self.cols==self.lines
        if method == "GE":
            return self.inverse_GE(iszerofunc=iszerofunc)
        elif method == "LU":
            return self.inverse_LU(iszerofunc=iszerofunc)
        elif method == "ADJ":
            return self.inverse_ADJ()
        else:
            raise Exception("Inversion method unrecognized")


    def __mathml__(self):
        mml = ""
        for i in range(self.lines):
            mml += "<matrixrow>"
            for j in range(self.cols):
                mml += self[i,j].__mathml__()
            mml += "</matrixrow>"
        return "<matrix>" + mml + "</matrix>"

    def row(self, i, f):
        """Elementary row operation using functor"""
        for j in range(0, self.cols):
            self[i, j] = f(self[i, j], j)

    def col(self, j, f):
        """Elementary column operation using functor"""
        for i in range(0, self.lines):
            self[i, j] = f(self[i, j], i)

    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]

    def col_swap(self, i, j):
        for k in range(0, self.lines):
            self[k, i], self[k, j] = self[k, j], self[k, i]

    def row_del(self, i):
        self.mat = self.mat[:i*self.cols] + self.mat[(i+1)*self.cols:]
        self.lines -= 1

    def col_del(self, i):
        """
        >>> import sympy
        >>> M = sympy.matrices.eye(3)
        >>> M.col_del(1)
        >>> M   #doctest: +NORMALIZE_WHITESPACE
        [1, 0]
        [0, 0]
        [0, 1]
        """
        for j in range(self.lines-1, -1, -1):
            del self.mat[i+j*self.cols]
        self.cols -= 1

    def row_join(self, rhs):
        """
        Concatenates two matrices along self's last and rhs's first column

        >>> from sympy import *
        >>> M = Matrix(3,3,lambda i,j: i+j)
        >>> V = Matrix(3,1,lambda i,j: 3+i+j)
        >>> M.row_join(V)
        [0, 1, 2, 3]
        [1, 2, 3, 4]
        [2, 3, 4, 5]
        """
        assert self.lines == rhs.lines
        newmat = self.zeros((self.lines, self.cols + rhs.cols))
        newmat[:,:self.cols] = self[:,:]
        newmat[:,self.cols:] = rhs
        return newmat

    def col_join(self, bott):
        """
        Concatenates two matrices along self's last and bott's first row

        >>> from sympy import *
        >>> M = Matrix(3,3,lambda i,j: i+j)
        >>> V = Matrix(1,3,lambda i,j: 3+i+j)
        >>> M.col_join(V)
        [0, 1, 2]
        [1, 2, 3]
        [2, 3, 4]
        [3, 4, 5]
        """
        assert self.cols == bott.cols
        newmat = self.zeros((self.lines+bott.lines, self.cols))
        newmat[:self.lines,:] = self[:,:]
        newmat[self.lines:,:] = bott
        return newmat

    def row_insert(self, pos, mti):
        """
        >>> from sympy import *
        >>> M = Matrix(3,3,lambda i,j: i+j)
        >>> M
        [0, 1, 2]
        [1, 2, 3]
        [2, 3, 4]
        >>> V = zeros((1, 3))
        >>> V
        [0, 0, 0]
        >>> M.row_insert(1,V)
        [0, 1, 2]
        [0, 0, 0]
        [1, 2, 3]
        [2, 3, 4]
        """
        if pos is 0:
            return mti.col_join(self)
        assert self.cols == mti.cols
        newmat = self.zeros((self.lines + mti.lines, self.cols))
        newmat[:pos,:] = self[:pos,:]
        newmat[pos:pos+mti.lines,:] = mti[:,:]
        newmat[pos+mti.lines:,:] = self[pos:,:]
        return newmat

    def col_insert(self, pos, mti):
        """
        >>> from sympy import *
        >>> M = Matrix(3,3,lambda i,j: i+j)
        >>> M
        [0, 1, 2]
        [1, 2, 3]
        [2, 3, 4]
        >>> V = zeros((3, 1))
        >>> V
        [0]
        [0]
        [0]
        >>> M.col_insert(1,V)
        [0, 0, 1, 2]
        [1, 0, 2, 3]
        [2, 0, 3, 4]
        """
        if pos is 0:
            return mti.row_join(self)
        assert self.lines == mti.lines
        newmat = self.zeros((self.lines, self.cols + mti.cols))
        newmat[:,:pos] = self[:,:pos]
        newmat[:,pos:pos+mti.cols] = mti[:,:]
        newmat[:,pos+mti.cols:] = self[:,pos:]
        return newmat

    def trace(self):
        assert self.cols == self.lines
        trace = 0
        for i in range(self.cols):
            trace += self[i,i]
        return trace

    def submatrix(self, keys):
        """
        >>> from sympy import *
        >>> m = Matrix(4,4,lambda i,j: i+j)
        >>> m   #doctest: +NORMALIZE_WHITESPACE
        [0, 1, 2, 3]
        [1, 2, 3, 4]
        [2, 3, 4, 5]
        [3, 4, 5, 6]
        >>> m[0:1, 1]   #doctest: +NORMALIZE_WHITESPACE
        [1]
        >>> m[0:2, 0:1] #doctest: +NORMALIZE_WHITESPACE
        [0]
        [1]
        >>> m[2:4, 2:4] #doctest: +NORMALIZE_WHITESPACE
        [4, 5]
        [5, 6]
        """
        assert isinstance(keys[0], slice) or isinstance(keys[1], slice)
        rlo, rhi = self.slice2bounds(keys[0], self.lines)
        clo, chi = self.slice2bounds(keys[1], self.cols)
        if not ( 0<=rlo<=rhi and 0<=clo<=chi ):
            raise IndexError("Slice indices out of range: a[%s]"%repr(keys))
        outLines, outCols = rhi-rlo, chi-clo
        outMat = [0]*outLines*outCols
        for i in xrange(outLines):
            outMat[i*outCols:(i+1)*outCols] = self.mat[(i+rlo)*self.cols+clo:(i+rlo)*self.cols+chi]
        return Matrix(outLines,outCols,outMat)

    def slice2bounds(self, key, defmax):
        """
            Takes slice or number and returns (min,max) for iteration
            Takes a default maxval to deal with the slice ':' which is (none, none)
        """
        if isinstance(key, slice):
            lo, hi = 0, defmax
            if key.start != None:
                if key.start >= 0:
                    lo = key.start
                else:
                    lo = defmax+key.start
            if key.stop != None:
                if key.stop >= 0:
                    hi = key.stop
                else:
                    hi = defmax+key.stop
            return lo, hi
        elif isinstance(key, int):
            if key >= 0:
                return key, key+1
            else:
                return defmax+key, defmax+key+1
        else:
            raise IndexError("Improper index type")

    def applyfunc(self, f):
        """
        >>> from sympy import *
        >>> m = Matrix(2,2,lambda i,j: i*2+j)
        >>> m   #doctest: +NORMALIZE_WHITESPACE
        [0, 1]
        [2, 3]
        >>> m.applyfunc(lambda i: 2*i)  #doctest: +NORMALIZE_WHITESPACE
        [0, 2]
        [4, 6]
        """
        assert callable(f)
        out = Matrix(self.lines,self.cols,map(f,self.mat))
        return out

    def evalf(self, prec=None, **options):
        if prec is None:
            return self.applyfunc(lambda i: i.evalf(**options))
        else:
            return self.applyfunc(lambda i: i.evalf(prec, **options))

    def reshape(self, _rows, _cols):
        """
        >>> from sympy import *
        >>> m = Matrix(2,3,lambda i,j: 1)
        >>> m   #doctest: +NORMALIZE_WHITESPACE
        [1, 1, 1]
        [1, 1, 1]
        >>> m.reshape(1,6)  #doctest: +NORMALIZE_WHITESPACE
        [1, 1, 1, 1, 1, 1]
        >>> m.reshape(3,2)  #doctest: +NORMALIZE_WHITESPACE
        [1, 1]
        [1, 1]
        [1, 1]
        """
        if self.lines*self.cols != _rows*_cols:
            print "Invalid reshape parameters %d %d" % (_rows, _cols)
        return Matrix(_rows, _cols, lambda i,j: self.mat[i*_cols + j])

    def print_nonzero (self, symb="X"):
        """
        Shows location of non-zero entries for fast shape lookup
        >>> from sympy import *
        >>> m = Matrix(2,3,lambda i,j: i*3+j)
        >>> m           #doctest: +NORMALIZE_WHITESPACE
        [0, 1, 2]
        [3, 4, 5]
        >>> m.print_nonzero()   #doctest: +NORMALIZE_WHITESPACE
        [ XX]
        [XXX]
        >>> m = matrices.eye(4)
        >>> m.print_nonzero("x")    #doctest: +NORMALIZE_WHITESPACE
        [x   ]
        [ x  ]
        [  x ]
        [   x]
        """
        s="";
        for i in range(self.lines):
            s+="["
            for j in range(self.cols):
                if self[i,j] == 0:
                    s+=" "
                else:
                    s+= symb+""
            s+="]\n"
        print s

    def LUsolve(self, rhs, iszerofunc=_iszero):
        """
        Solve the linear system Ax = b.
        self is the coefficient matrix A and rhs is the right side b.
        """
        assert rhs.lines == self.lines
        A, perm = self.LUdecomposition_Simple(iszerofunc=_iszero)
        n = self.lines
        b = rhs.permuteFwd(perm)
        # forward substitution, all diag entries are scaled to 1
        for i in range(n):
            for j in range(i):
                b.row(i, lambda x,k: x - b[j,k]*A[i,j])
        # backward substitution
        for i in range(n-1,-1,-1):
            for j in range(i+1, n):
                b.row(i, lambda x,k: x - b[j,k]*A[i,j])
            b.row(i, lambda x,k: x / A[i,i])
        return b

    def LUdecomposition(self, iszerofunc=_iszero):
        """
        Returns the decompositon LU and the row swaps p.
        """
        combined, p = self.LUdecomposition_Simple(iszerofunc=_iszero)
        L = self.zeros(self.lines)
        U = self.zeros(self.lines)
        for i in range(self.lines):
            for j in range(self.lines):
                if i > j:
                    L[i,j] = combined[i,j]
                else:
                    if i == j:
                        L[i,i] = 1
                    U[i,j] = combined[i,j]
        return L, U, p

    def LUdecomposition_Simple(self, iszerofunc=_iszero):
        """
        Returns A compused of L,U (L's diag entries are 1) and
        p which is the list of the row swaps (in order).
        """
        assert self.lines == self.cols
        n = self.lines
        A = self[:,:]
        p = []
        # factorization
        for j in range(n):
            for i in range(j):
                for k in range(i):
                    A[i,j] = A[i,j] - A[i,k]*A[k,j]
            pivot = -1
            for i in range(j,n):
                for k in range(j):
                    A[i,j] = A[i,j] - A[i,k]*A[k,j]
                # find the first non-zero pivot, includes any expression
                if pivot == -1 and not iszerofunc(A[i,j]):
                    pivot = i
            if pivot < 0:
                raise "Error: non-invertible matrix passed to LUdecomposition_Simple()"
            if pivot != j: # row must be swapped
                A.row_swap(pivot,j)
                p.append([pivot,j])
            assert not iszerofunc(A[j,j])
            scale = 1 / A[j,j]
            for i in range(j+1,n):
                A[i,j] = A[i,j] * scale
        return A, p


    def LUdecompositionFF(self):
        """
        Returns 4 matrices P, L, D, U such that PA = L D**-1 U.

        From the paper "fraction-free matrix factors..." by Zhou and Jeffrey
        """
        n, m = self.lines, self.cols
        U, L, P = self[:,:], eye(n), eye(n)
        DD = zeros(n) # store it smarter since it's just diagonal
        oldpivot = 1

        for k in range(n-1):
            if U[k,k] == 0:
                kpivot = k+1
                Notfound = True
                while kpivot < n and Notfound:
                    if U[kpivot, k] != 0:
                        Notfound = False
                    else:
                        kpivot = kpivot + 1
                if kpivot == n+1:
                    raise "Matrix is not full rank"
                else:
                    swap = U[k, k:]
                    U[k,k:] = U[kpivot,k:]
                    U[kpivot, k:] = swap
                    swap = P[k, k:]
                    P[k, k:] = P[kpivot, k:]
                    P[kpivot, k:] = swap
            assert U[k, k] != 0
            L[k,k] = U[k,k]
            DD[k,k] = oldpivot * U[k,k]
            assert DD[k,k] != 0
            Ukk = U[k,k]
            for i in range(k+1, n):
                L[i,k] = U[i,k]
                Uik = U[i,k]
                for j in range(k+1, m):
                    U[i,j] = (Ukk * U[i,j] - U[k,j]*Uik) / oldpivot
                U[i,k] = 0
            oldpivot = U[k,k]
        DD[n-1,n-1] = oldpivot
        return P, L, DD, U

    def cofactorMatrix(self, method="berkowitz"):
        out = Matrix(self.lines, self.cols, lambda i,j:
                self.cofactor(i, j, method))
        return out

    def minorEntry(self, i, j, method="berkowitz"):
        assert 0 <= i < self.lines and 0 <= j < self.cols
        return self.minorMatrix(i,j).det(method)

    def minorMatrix(self, i, j):
        assert 0 <= i < self.lines and 0 <= j < self.cols
        return self.delRowCol(i,j)

    def cofactor(self, i, j, method="berkowitz"):
        if (i+j) % 2 == 0:
            return self.minorEntry(i, j, method)
        else:
            return -1 * self.minorEntry(i, j, method)

    def jacobian(self, varlist):
        """
        Calculates the Jacobian matrix (derivative of a vectorial function).

        self is a vector of expression representing functions f_i(x_1, ...,
        x_n).  varlist is the set of x_i's in order.
        """
        assert self.lines == 1
        m = self.cols
        if isinstance(varlist, Matrix):
            assert varlist.lines == 1
            n = varlist.cols
        elif isinstance(varlist, (list, tuple)):
            n = len(varlist)
        assert n > 0 # need to diff by something
        J = self.zeros((m, n)) # maintain subclass type
        for i in range(m):
            if isinstance(self[i], (float, int)):
                continue    # constant function, jacobian row is zero
            try:
                tmp = self[i].diff(varlist[0])   # check differentiability
                J[i,0] = tmp
            except AttributeError:
                raise ValueError("Function %d is not differentiable" % i)
            for j in range(1,n):
                J[i,j] = self[i].diff(varlist[j])
        return J

{'fraction_non_alphanumeric': 0.13069807679510215, 'fraction_numerical': 0.014374126572170094, 'mean_word_length': 1.4812796697626418, 'pattern_counts': {'":': 6, '<': 53, '<?xml version=': 0, '>': 461, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2103598', 'n_tokens_mistral': 19976, 'n_tokens_neox': 18753, 'n_words': 5176}