def main():
    a=Symbol("a", real=True)
    b=Symbol("b", real=True)
    c=Symbol("c", real=True)

    p = (a,b,c)

    assert u(p, 1).D * u(p, 2) == Matrix(1, 1, [0])
    assert u(p, 2).D * u(p, 1) == Matrix(1, 1, [0])

    p1,p2,p3 =[Symbol(x, real=True) for x in ["p1","p2","p3"]]
    pp1,pp2,pp3 =[Symbol(x, real=True) for x in ["pp1","pp2","pp3"]]
    k1,k2,k3 =[Symbol(x, real=True) for x in ["k1","k2","k3"]]
    kp1,kp2,kp3 =[Symbol(x, real=True) for x in ["kp1","kp2","kp3"]]

    p = (p1,p2,p3)
    pp = (pp1,pp2,pp3)

    k = (k1,k2,k3)
    kp = (kp1,kp2,kp3)

    mu = Symbol("mu")

    e = (pslash(p)+m*ones(4))*(pslash(k)-m*ones(4))
    f = pslash(p)+m*ones(4)
    g = pslash(p)-m*ones(4)


    #pprint(e)
    xprint( 'Tr(f*g)', Tr(f*g) )
    #print Tr(pslash(p) * pslash(k)).expand()

    M0 = [ ( v(pp, 1).D * mgamma(mu) * u(p, 1) ) * ( u(k, 1).D * mgamma(mu,True) * \
             v(kp, 1) ) for mu in range(4)]
    M = M0[0]+M0[1]+M0[2]+M0[3]
    M = M[0]
    assert isinstance(M, Basic)
    #print M
    #print simplify(M)

    d=Symbol("d", real=True) #d=E+m

    xprint('M', M)
    print "-"*40
    M = ((M.subs(E,d-m)).expand() * d**2 ).expand()
    xprint('M2', 1/(E+m)**2 * M)
    print "-"*40
    x,y= M.as_real_imag()
    xprint('Re(M)', x)
    xprint('Im(M)', y)
    e = x**2+y**2
    xprint('abs(M)**2', e)
    print "-"*40
    xprint('Expand(abs(M)**2)', e.expand())

def test_creation():
    """
    Check that matrix dimensions can be specified using any reasonable type
    (see issue 1515).
    """
    raises(ValueError, 'zeros((3, 0))')
    raises(ValueError, 'zeros((1,2,3,4))')
    assert zeros(3L) == zeros(3)
    assert zeros(Integer(3)) == zeros(3)
    assert zeros(3.) == zeros(3)
    assert eye(3L) == eye(3)
    assert eye(Integer(3)) == eye(3)
    assert eye(3.) == eye(3)
    assert ones((3L, Integer(4))) == ones((3, 4))

def make_state_space_sym(num_signals, num_states, is_homo):

    mu_rho_dict_sym = {}
    state_state_sym_dict = {}

    mu_names = ['mu_'+str(i) for i in range(num_states)]
    rho_names = ['rho_'+str(i) for i in range(num_states)]
    mu_rho_names = mu_names + rho_names

    mu_names_sym = [sympy.Symbol(x) for x in mu_names]
    rho_names_sym = [sympy.Symbol(x) for x in rho_names]
    mu_rho_names_sym = [sympy.Symbol(x) for x in mu_rho_names]
    mu_rho_dict_sym.update(dict(zip(mu_rho_names, mu_rho_names_sym)))

    A_identity_sym = sympy.eye(num_states)
    A_rhoblock_sym = sympy.diag(*rho_names_sym)
    A_sym = sympy.diag(A_identity_sym, A_rhoblock_sym)

    D_sym_mu_part = sympy.ones(num_signals, num_states)
    D_sym_zeta_part = sympy.ones(num_signals, num_states)
    D_sym = D_sym_mu_part.row_join(D_sym_zeta_part)

    if is_homo:
        sigmas_signal_names = [
            'sigma_signal_'+str(i) for i in range(num_signals)]
        sigmas_state_names = ['sigma_state_'+str(i) for i in range(num_states)]
        sigmas_signal_sym = [sympy.Symbol(x) for x in sigmas_signal_names]
        sigmas_state_sym = [sympy.Symbol(x) for x in sigmas_state_names]

        C_nonsingularblock_sym = sympy.diag(*sigmas_state_sym)
        G_nonsingularblock_sym = sympy.diag(*sigmas_signal_sym)

        C_singularblock_sym = sympy.zeros(num_states, num_states)
        G_singularblock_sym = sympy.zeros(num_signals, num_states)

        G_sym = G_singularblock_sym.row_join(G_nonsingularblock_sym)
        C_sym = sympy.diag(C_singularblock_sym, C_nonsingularblock_sym)

    main_matrices_sym = {
        'A_z': A_sym, 'C_z': C_sym, 'D_s': D_sym, 'G_s': G_sym}
    sub_matrices_sym = {'A_z_stable': A_rhoblock_sym,
                        'C_z_nonsingular': C_nonsingularblock_sym,
                        'G_s_nonsingular': G_nonsingularblock_sym}

    state_state_sym_dict.update(main_matrices_sym)
    state_state_sym_dict.update(sub_matrices_sym)

    return state_state_sym_dict

def cholesky(A):
    """
    # A is positive definite mxm
    """
    assert A.shape[0] == A.shape[1]
    # assert all(A.eigenvals() > 0)
    m = A.shape[0]
    N = deepcopy(A)
    D = ones(*A.shape)
    for i in xrange(m - 1):
        for j in xrange(i + 1, m):
            N[j, i] = N[i, j]
            D[j, i] = D[i, j]
            n, d = ratior(N[i, j], D[i, j], N[i, i], D[i, i])
            N[i, j], D[i, j] = n, d
            if verbose_chol:
                print "i={}, j={}".format(i + 1, j + 1)
                print "N:"
                printnp(N)
                print "D:"
                printnp(D)
        for k in xrange(i + 1, m):
            for l in xrange(k, m):
                n, d = multr(N[k, i], D[k, i], N[i, l], D[i, l])
                N[k, l], D[k, l] = subr(N[k, l], D[k, l], n, d)
                if verbose_chol:
                    print "k={}, l={}".format(k + 1, l + 1)
                    print "N:"
                    printnp(N)
                    print "D:"
                    printnp(D)
    return N, D

    def test_matrix_tensor_product():
        l1 = zeros(4)
        for i in range(16):
            l1[i] = 2**i
        l2 = zeros(4)
        for i in range(16):
            l2[i] = i
        l3 = zeros(2)
        for i in range(4):
            l3[i] = i
        vec = Matrix([1,2,3])

        #test for Matrix known 4x4 matricies
        numpyl1 = np.matrix(l1.tolist())
        numpyl2 = np.matrix(l2.tolist())
        numpy_product = np.kron(numpyl1,numpyl2)
        args = [l1, l2]
        sympy_product = matrix_tensor_product(*args)
        assert numpy_product.tolist() == sympy_product.tolist()
        numpy_product = np.kron(numpyl2,numpyl1)
        args = [l2, l1]
        sympy_product = matrix_tensor_product(*args)
        assert numpy_product.tolist() == sympy_product.tolist()

        #test for other known matrix of different dimensions
        numpyl2 = np.matrix(l3.tolist())
        numpy_product = np.kron(numpyl1,numpyl2)
        args = [l1, l3]
        sympy_product = matrix_tensor_product(*args)
        assert numpy_product.tolist() == sympy_product.tolist()
        numpy_product = np.kron(numpyl2,numpyl1)
        args = [l3, l1]
        sympy_product = matrix_tensor_product(*args)
        assert numpy_product.tolist() == sympy_product.tolist()

        #test for non square matrix
        numpyl2 = np.matrix(vec.tolist())
        numpy_product = np.kron(numpyl1,numpyl2)
        args = [l1, vec]
        sympy_product = matrix_tensor_product(*args)
        assert numpy_product.tolist() == sympy_product.tolist()
        numpy_product = np.kron(numpyl2,numpyl1)
        args = [vec, l1]
        sympy_product = matrix_tensor_product(*args)
        assert numpy_product.tolist() == sympy_product.tolist()

        #test for random matrix with random values that are floats
        random_matrix1 = np.random.rand(np.random.rand()*5+1,np.random.rand()*5+1)
        random_matrix2 = np.random.rand(np.random.rand()*5+1,np.random.rand()*5+1)
        numpy_product = np.kron(random_matrix1,random_matrix2)
        args = [Matrix(random_matrix1.tolist()),Matrix(random_matrix2.tolist())]
        sympy_product = matrix_tensor_product(*args)
        assert not (sympy_product - Matrix(numpy_product.tolist())).tolist() > \
        (ones((sympy_product.rows,sympy_product.cols))*epsilon).tolist()

        #test for three matrix kronecker
        sympy_product = matrix_tensor_product(l1,vec,l2)

        numpy_product = np.kron(l1,np.kron(vec,l2))
        assert numpy_product.tolist() == sympy_product.tolist()

def met_zeyd(A, b):
    C, d = iteration_view(A, b)
    H = sympy.zeros(3)
    F = sympy.zeros(3)
    c = sympy.zeros(3,1)
    for i in xrange(3):
        c[i] = d[i]
        for j in xrange(3):
            if i > j:
                H[i, j] = C[i, j]
            else:
                F[i, j] = C[i, j]
    print "\nx = Cx + d\n"
    print "C = \n", C, "\n", "\nd = \n", d
    print "\nConvergence: ", convergence_mzeyd(C, d)
    if convergence_mzeyd(C, d):
        E = sympy.eye(3)
        x0 = sympy.ones(3, 1)
        x1 = (E-H).inv()*F*x0 + (E-H).inv()*c
        while ((x1-x0)[0] > 0.00001 or (x1-x0)[1] > 0.00001 or\
              (x1-x0)[2] > 0.00001 or (x0-x1)[0] > 0.00001 or\
              (x0-x1)[1] > 0.00001 or (x0-x1)[2] > 0.00001):
              x0 = x1
              x1 = (E-H).inv()*F*x0 + (E-H).inv()*c
        print "\nSolution:" 

    return [element for element in x1]

def get_matrix_of_converted_atoms(Nu, positions, pending_conversion, natural_influence, Omicron, D):
    """
    :param Nu: A matrix with a shape=(<number of Matters in Universe>, <number of Atoms in Universe>) where
    each Nu[i,j] stands for how many atoms of type j in matter of type i.
    :type Nu: Matrix
    :param ps: positions of matters
    :type ps: [Matrix]
    :return:
    """

    x, y = symbols('x y')

    number_of_matters, number_of_atoms = Nu.shape

    M = zeros(0, number_of_atoms)

    if number_of_matters != len(positions):
        raise Exception("Parameters shapes mismatch.")

    for (i, position) in enumerate(positions):
        (a, b) = tuple(position)
        K = get_conversion_ratio_matrix(pending_conversion, Nu[i, :])
        M = M.col_join(((diag(*(ones(1, number_of_atoms)*diag(*K)*Omicron.transpose()))*D).transpose() *
                        natural_influence).transpose().subs({x: a, y: b}))

    return M.evalf()

def initialise_working_matrices(G):
    """  G is a nonzero matrix with at least two rows.  """
    B = eye(G.shape[0])
    # Lower triang matrix
    L = zeros(G.shape[0], G.shape[0])
    D = ones(G.shape[0] + 1, 1)
    A = Matrix(G)
    return A, B, L, D

def elementary_weight(tree,s,arrays,method):
    """
        Constructs elementary weights for a Volterra Runge-Kutta method,
        supposing the row sum condition.
        The output needs to be multiplied by b^T and equated to LHS
        to obtain the order condition.
        
        It is used by gen_order_conditions in vrk_methods
        
        INPUT:
            
        - tree   -- input tree, must be a RootedTree.
        - s      -- number of stages
        - arrays -- it depends on the method input, is a list containing two
                    or three arrays, the order should be: c,A,e,D (if VRK)
                    c,A,d (if BVRK) or c,A if (if PVRK)
        - method -- select which type of method to use: VRK, BVRK or if a PVRK
                    method is wanted, the three must be created with the xa='a' flag.
        
        OUTPUT: it is a column vector belonging to 
                <type 'numpy.ndarray'>
                e.g. like array([[],...,[]], dtype=object)

    """
    from sympy import eye, ones
    if tree=='': return ''    
    
    u=np.array(ones((s,1)))   #np.ones((s, 1), dtype=np.int)
    if tree=='a': return u # Matrix(u)

    I=np.array(eye(s))    # np.eye(s, dtype=np.int)
    ew=u.copy()
    
    c=arrays[0]
    A=arrays[1]
    if method=='VRK':
        d=arrays[2]
        D=arrays[3]
    elif method=='BVRK':
        d=arrays[2]
    else: #method=='PVRK'
        d=arrays[0]
    
    
    
    nx,na,subtrees=tree._parse_subtrees()
    ew*=c**na #na and nx can also be zero, then we have a vector of ones.
    ew*=d**nx
    
    # Two curly bracket contain at least a symbol, 'a' or 'x', i.e {} is not a leaf
    if len(subtrees)>0:
        for subtree in subtrees:
            if method=='VRK':
                ew=ew*_elem_weight_sub3(subtree,method,c,D,A,I,u)
            else:
                ew=ew*_elem_weight_sub3(subtree,method,c,d,A,I,u)
    
    return ew #returns a column

def _pressure_tensor(grid):
    press = [Symbol('flux[%d]' % i) for i in range(grid.dim * (grid.dim + 1) / 2)]
    P = sympy.ones(grid.dim)  # P_ab - rho cs^2 \delta_ab
    k = 0
    for i in range(grid.dim):
        for j in range(i, grid.dim):
            P[i, j] = press[k]
            P[j, i] = press[k]
            k += 1
    return P

def construct_beta(Kbeta):
    # symbole associé au petit t (c'est à dire l'instant dans la période de suivi)
    t = sy.Symbol("t")
    # smbole associé au grand T (c'est à dire la durée de suivie)
    s = sy.Symbol("s")
    syPhi = sy.ones(Kbeta, 1)
    syb = sy.ones(1, Kbeta)
    b = [[] for k in range(Kbeta)]
    v = [np.arange(np.sqrt(Kbeta)), np.arange(np.sqrt(Kbeta))]
    expo = cg.expandnp(v)
    for x in range(len(expo[:, 0])):
        syPhi[x] = (t ** expo[x, 0]) * (s ** expo[x, 1])
        syb[x] = sy.Symbol("b" + str(x))
        b[x] = sy.Symbol("b" + str(x))
    syBeta = syb * syPhi
    syBeta = syBeta[0, 0]
    arg = [t, s] + b
    Beta_fonc_est = sy.lambdify(tuple(arg), syBeta, "numpy")
    return Beta_fonc_est

def z(A):
    """Return the (symbolic) vector of rates toward absorbing state.

    Args:
        A (SymPy dxd-matrix): compartment matrix

    Returns:
        SymPy dx1-matrix: :math:`\\bf{z} = -B^T\\,\\bf{1}`
    """
    o = ones(A.rows, 1)
    return -A.transpose()*o

 def lagrange(self, nodes):
     """
     Lagrange polynomial 
     """
     length = len(nodes)
     r = sympy.Symbol('r')
     phi = sympy.ones(1, length)
     for k in range(length):
         for l in range(length):
             if (k != l):
                 phi[k] *= (r - nodes[l])/(nodes[k] - nodes[l])
     return phi 

def test_zeros_ones_fill():
    n, m = 3, 5

    a = zeros( (n, m) )
    a.fill( 5 )

    b = 5 * ones( (n, m) )

    assert a == b
    assert a.rows == b.rows == 3
    assert a.cols == b.cols == 5
    assert a.shape == b.shape == (3, 5)

def Pre_Comp_YX(L, T, Xdata, Y, Kbeta, J):
    t = sy.Symbol("t")
    s = sy.Symbol("s")
    # Récupération des variables et paramètres
    N = len(L)
    D = len(L[0])
    # ----------------- Construction de la base fonctionnelle
    syPhi = sy.ones(Kbeta ** 2, 1)
    syb = sy.ones(1, Kbeta ** 2)
    v = [np.arange(Kbeta), np.arange(Kbeta)]
    expo = cg.expandnp(v)
    Phi_fonc = [[] for j in range(Kbeta ** 2)]
    for x in range(len(expo[:, 0])):
        syPhi[x] = (t ** expo[x, 0]) * (s ** expo[x, 1])
        Phi_fonc[x] = sy.lambdify((t, s), syPhi[x], "numpy")
        syb[x] = sy.Symbol("b" + str(x))
    syBeta = syb * syPhi
    Phi_mat = Comp_Phi(Phi_fonc, T, J)
    I_pen = J22_fast(syPhi, np.max(T), 50)[3]
    # ----------------- Construction des noyaux et leurs dérivées
    # Construction de la forme du noyau
    el1 = sy.Symbol("el1")
    per1 = sy.Symbol("per1")
    sig1 = sy.Symbol("sig1")
    args1 = [el1, per1, sig1]
    el2 = sy.Symbol("el2")
    sig2 = sy.Symbol("sig2")
    args2 = [el2, sig2]
    syk = cg.sy_Periodic((s, t), *args1) + cg.sy_RBF((s, t), *args2)
    args = [t, s] + args1 + args2
    # Dérivation et construction des fonctions vectorielles associées
    k_fonc = sy.lambdify(tuple(args), syk, "numpy")
    n_par = len(args) - 2
    k_der = [[] for i in range(n_par)]
    for i in range(n_par):
        func = syk.diff(args[i + 2])
        k_der[i] = sy.lambdify(tuple(args), func, "numpy")
    return (Phi_mat, k_fonc, k_der, I_pen)

def cum_dist_func(beta, A, Qt):
    """Return the (symbolic) cumulative distribution function of phase-type.

    Args:
        beta (SymPy dx1-matrix): initial distribution vector
        A (SymPy dxd-matrix): transition rate matrix
        Qt (SymPy dxd-matrix): Qt = :math:`e^{t\\,\\bf{A}}`

    Returns:
        SymPy expression: cumulative distribution function of PH(:math:`\\bf{\\beta}`, :math:`\\bf{A}`)
            :math:`F_T(t) = 1 - \\bf{1}^T\\,e^{t\\,\\bf{A}}\\,\\bf{\\beta}`
    """
    o = ones(1, A.cols)
    return 1 - (o * (Qt * beta))[0]

def nth_moment(beta, A, n):
    """Return the (symbolic) nth moment of the phase-type distribution.

    Args:
        beta (SymPy dx1-matrix): initial distribution vector
        A (SymPy dxd-matrix): transition rate matrix
        n (positive int): order of the moment
    
    Returns:
        SymPy expression: nth moment of PH(:math:`\\bf{\\beta}`, :math:`\\bf{A}`)
            :math:`\mathbb{E}[T^n] = (-1)^n\\,n!\\,\\bf{1}^T\\,\\bf{A}^{-1}\\,\\bf{\\beta}`
    """
    o = ones(1, A.cols)
    return ((-1)**n*factorial(n)*o*(A**-n)*beta)[0]

def cholesky(A):
    """
    # A is positive definite mxm
    """
    assert A.shape[0] == A.shape[1]
    # assert all(A.eigenvals() > 0)
    m = A.shape[0]
    N = deepcopy(A)
    D = ones(*A.shape)
    for i in xrange(m - 1):
        for j in xrange(i + 1, m):
            N[j, i] = N[i, j]
            D[j, i] = D[i, j]
            n, d = ratior(N[i, j], D[i, j], N[i, i], D[i, i])
            N[i, j], D[i, j] = n, d
        for k in xrange(i + 1, m):
            for l in xrange(k, m):
                n, d = multr(N[k, i], D[k, i], N[i, l], D[i, l])
                N[k, l], D[k, l] = subr(N[k, l], D[k, l], n, d)
    return N, D

    def relaxation_local(self, m, with_rel_velocity=False):
        """
        Return symbolic expression which computes the relaxation operator.

        Parameters
        ----------

        m : SymPy Matrix
            indexed objects for the moments

        with_rel_velocity : boolean
            check if the scheme uses relative velocity.
            (default is False)

        """
        if with_rel_velocity:
            eq = (self.Tu*self.eq).subs(list(zip(self.mv, m)))
        else:
            eq = self.eq.subs(list(zip(self.mv, m)))
        relax = (sp.ones(*self.s.shape) - self.s).multiply_elementwise(sp.Matrix(m)) + self.s.multiply_elementwise(eq)
        alltogether(relax)
        return Eq(m, relax)

 def lagrangeDeri(self, nodes):
     """
     Lagrange matrix at the nodes is just an Identity
     We'll come back to interpolation at points other than nodes
     at a later time
     Here we create derivative operator at the nodes
     Lagrange polynomial is
     phi = Product(l, l.neq.k) (r - r_l)/(r_k - r_l)
     r_i are the nodes
     """
     length = len(nodes)
     r = sympy.Symbol('r')
     phi = sympy.ones(1, length)
     dPhi = sympy.zeros(length, length)
     for k in range(length):
         for l in range(length):
             if (k != l):
                 phi[k] *= (r - nodes[l])/(nodes[k] - nodes[l])
     for k in range(length):
         for l in range(length):
             dPhi[k, l] = sympy.diff(phi[l]).evalf(subs = {r: nodes[k]})
     return dPhi