def fidelity(statevec, dm, u_dm, ustatevec=np.array([0,0])):
    '''
    returns the fidelity (and its uncertainty) of the measured density
    matrix with a given state vector.
    '''
   
    f = error.Formula()
    beta,a,x,b,alpha = sympy.symbols('beta,a,x,b,alpha')
    v = Matrix([alpha, beta])
    rho = Matrix([[x,a+1j*b],[a-1j*b, 1-x]])    
    f.formula = (v.conjugate().transpose() * rho * v)[0]
    
    f.values[alpha] = statevec[0]
    f.values[beta] = statevec[1]
    f.values[a]=float(np.real(dm[0,1]))
    f.values[x]=float(dm[0,0])
    f.values[b]=float(np.imag(dm[0,1]))

    f.uncertainties[alpha]=ustatevec[0]
    f.uncertainties[beta]=ustatevec[1]
    f.uncertainties[x]=u_dm[0,0]
    f.uncertainties[a]=float(np.real(u_dm[0,1]))
    f.uncertainties[b]=float(np.imag(u_dm[0,1]))
    
    _fid,_ufid = f.num_eval()
    fid = float(_fid.as_real_imag()[0])
    ufid = float(_ufid.as_real_imag()[0])
   
    return (fid,ufid)

def demo_symetricA_Msquare(A):
    print "\n1)Any symmetric matrix A can be written as A = M^2:\n"\
            "<= any symetric matrix A can be diagonalized as A = Q * D * Q.inv()\n"\
            "<= D is diagonal matrix with eigen values\n"\
            "   Q is eigen vectors with norm 1\n"\
            "   Q.inv() == Q.T\n"\
            "first take a look as eigen vectors:\n"
    pprint(A.eigenvects())
    print "\nthen sympy diagonalized result:\n"
    pprint(A.diagonalize())
    d = A.diagonalize()[1]
    q = A.diagonalize()[0]
    q = Matrix(q.rows, q.cols, lambda i, j: q[:,j].normalized()[i])
    print "\nthen normalized Q:\n"
    pprint(q)
    print "\nthen the transpose of Q:\n"
    pprint(q.T)
    print "\nthen the inverse of Q:\n"
    pprint(q.inv())
    print "\nthen Q * D * Q.inv():\n"
    pprint(q*d*q.inv())
    print "\nif we define a new diagonal matrix Droot:\n"
    d = d.applyfunc(lambda x: x**.5)
    pprint(d)
    print "\nthen let M = Q * Droot * Q.inv():\n"
    m = q*d*q.inv()
    pprint(m)
    print "\nthen A = M*M:\n"
    pprint(q*d*d*q.inv())

    def distance(self, o):
        """Distance beteen the plane and another geometric entity.

        Parameters
        ==========

        Point3D, LinearEntity3D, Plane.

        Returns
        =======

        distance

        Notes
        =====

        This method accepts only 3D entities as it's parameter, but if you want
        to calculate the distance between a 2D entity and a plane you should
        first convert to a 3D entity by projecting onto a desired plane and
        then proceed to calculate the distance.

        Examples
        ========

        >>> from sympy import Point, Point3D, Line, Line3D, Plane
        >>> a = Plane(Point3D(1, 1, 1), normal_vector=(1, 1, 1))
        >>> b = Point3D(1, 2, 3)
        >>> a.distance(b)
        sqrt(3)
        >>> c = Line3D(Point3D(2, 3, 1), Point3D(1, 2, 2))
        >>> a.distance(c)
        0

        """
        from sympy.geometry.line3d import LinearEntity3D
        x, y, z = map(Dummy, 'xyz')
        if self.intersection(o) != []:
            return S.Zero

        if isinstance(o, Point3D):
           x, y, z = map(Dummy, 'xyz')
           k = self.equation(x, y, z)
           a, b, c = [k.coeff(i) for i in (x, y, z)]
           d = k.xreplace({x: o.args[0], y: o.args[1], z: o.args[2]})
           t = abs(d/sqrt(a**2 + b**2 + c**2))
           return t
        if isinstance(o, LinearEntity3D):
            a, b = o.p1, self.p1
            c = Matrix(a.direction_ratio(b))
            d = Matrix(self.normal_vector)
            e = c.dot(d)
            f = sqrt(sum([i**2 for i in self.normal_vector]))
            return abs(e / f)
        if isinstance(o, Plane):
            a, b = o.p1, self.p1
            c = Matrix(a.direction_ratio(b))
            d = Matrix(self.normal_vector)
            e = c.dot(d)
            f = sqrt(sum([i**2 for i in self.normal_vector]))
            return abs(e / f)

def demo_corr_bound(n):
    def corr(i, j):
        if i == j:
            return 1
        else:
            return rho
    Rho = Matrix(n, n, corr)
    print "what's the bound of rho to make below a correlation matrix?\n"
    pprint(Rho)
    print "\ncan be viewed as the sum of 2 matrices Rho = A + B:\n"
    A = (1-rho)*eye(n)
    B = rho*ones(n,n)
    pprint(A)
    pprint(B)
    print "\nthe eigen value and its dimention of first matrix A is:"
    pprint(A.eigenvects())
    print "\nas for the seconde matrix B\n"\
            "it's product of any vector v:"
    v = IndexedBase('v')
    vec = Matrix(n, 1, lambda i, j: v[i+1])
    pprint(vec)
    pprint(B*vec)
    print "\nin order for it's to equal to a linear transform of v\n"\
            "we can see its eigen values, vectors and dimentions are:"
    pprint(B.eigenvects())
    print "\nfor any eigen vector of B v, we have Rho*v :\n"
    pprint(Rho*vec)
    print "\nwhich means v is also an eigen vector of Rho,\n"\
            "the eigen values, vectors and dimentions of Rho are:\n"
    pprint(Rho.eigenvects())
    print "\nsince have no negative eigen values <=> positive semidefinite:\n"\
            "the boundaries for rho are: [1/(%d-1),1]" %n

    def perpendicular_plane(self, pt):
        """
        Plane perpendicular to the given plane and passing through the point pt.

        Parameters
        ==========

        pt: Point3D

        Returns
        =======

        Plane

        Examples
        ========

        >>> from sympy import Plane, Point3D
        >>> a = Plane(Point3D(1,4,6), normal_vector=[2, 4, 6])
        >>> a.perpendicular_plane(Point3D(2, 3, 5))
        Plane(Point3D(2, 3, 5), [2, 8, -6])

        """
        a = Matrix(self.normal_vector)
        b, c = pt, self.p1
        d = Matrix(b.direction_ratio(c))
        e = list(d.cross(a))
        return Plane(pt, normal_vector=e)

 def _eval_transpose(self):
     # Flip all the individual matrices
     matrices = [Transpose(matrix) for matrix in self.blocks]
     # Make a copy
     M = Matrix(self.blockshape[0], self.blockshape[1], matrices)
     # Transpose the block structure
     M = M.transpose()
     return BlockMatrix(M)

def evalMat( Mlist,  blist  ):
	M = Matrix(Mlist)
	b = Matrix(blist)
	c = M.LUsolve( b )
	print 
	print " ================= Case: ",b.transpose()
	print 
	print c
	return c

def test_eigenvals():
    M = EigenOnlyMatrix([[0, 1, 1],
                [1, 0, 0],
                [1, 1, 1]])
    assert M.eigenvals() == {2*S.One: 1, -S.One: 1, S.Zero: 1}

    # if we cannot factor the char poly, we raise an error
    m = Matrix([[3, 0, 0, 0, -3], [0, -3, -3, 0, 3], [0, 3, 0, 3, 0], [0, 0, 3, 0, 3], [3, 0, 0, 3, 0]])
    raises(MatrixError, lambda: m.eigenvals())

	def cleanup(self):
		symbolDictionary = dict(zip(self.state + self.velocity, self.cleanState + self.cleanVelocity))
		N                = len(self._fullCoordinate)
		cleanMatrix      = Matrix([simplify(element).subs(symbolDictionary) for element in  self.mass])
		self._cleanMass  = list([cleanMatrix.reshape(N, N)])
		self._cleanForce = [simplify(element).subs(symbolDictionary) for element in self.force]
		self.SpecialFunctions(symbolDictionary)
		self.cleanConstraint      = simplify(self.constraint_equation.subs(symbolDictionary))
		self.cleanConstraintForce = [simplify(element).subs(symbolDictionary) for element in self.virtualForce]
		self.cleanNullForce       = self.nullForce

def euler101():
    #actual function that generates polynomials
    def u(n):
        return n ** 3

    order = 3  # largest polynomial order of u
    seq = [u(x) for x in range(1, order + 2)]
    A = Matrix([[1 ** 1, 1], [8 ** 1, 1]][::-1])
    b = Matrix([[x] for x in seq[:2]])
    return A.inv() * b

 def is_scalar_multiple(p1, p2):
     """Returns whether `p1` and `p2` are scalar multiples
     of eachother.
     """
     # if the vectors p1 and p2 are linearly dependent, then they must
     # be scalar multiples of eachother
     m = Matrix([p1.args, p2.args])
     # XXX: issue #9480 we need `simplify=True` otherwise the
     # rank may be computed incorrectly
     return m.rank(simplify=True) < 2

    def is_concyclic(self, *args):
        """Do `self` and the given sequence of points lie in a circle?

        Returns True if the set of points are concyclic and
        False otherwise. A trivial value of True is returned
        if there are fewer than 2 other points.

        Parameters
        ==========

        args : sequence of Points

        Returns
        =======

        is_concyclic : boolean


        Examples
        ========

        >>> from sympy import Point

        Define 4 points that are on the unit circle:

        >>> p1, p2, p3, p4 = Point(1, 0), (0, 1), (-1, 0), (0, -1)

        >>> p1.is_concyclic() == p1.is_concyclic(p2, p3, p4) == True
        True

        Define a point not on that circle:

        >>> p = Point(1, 1)

        >>> p.is_concyclic(p1, p2, p3)
        False

        """
        points = (self,) + args
        points = Point._normalize_dimension(*[Point(i) for i in points])
        points = list(uniq(points))
        if not Point.affine_rank(*points) <= 2:
            return False
        origin = points[0]
        points = [p - origin for p in points]
        # points are concyclic if they are coplanar and
        # there is a point c so that ||p_i-c|| == ||p_j-c|| for all
        # i and j.  Rearranging this equation gives us the following
        # condition: the matrix `mat` must not a pivot in the last
        # column.
        mat = Matrix([list(i) + [i.dot(i)] for i in points])
        rref, pivots = mat.rref()
        if len(origin) not in pivots:
            return True
        return False

def test_1xN_vecs():
    gl = glsl_code
    for i in range(1,10):
        A = Matrix(range(i))
        assert gl(A.transpose()) == gl(A)
        assert gl(A,mat_transpose=True) == gl(A)
        if i > 1:
            if i <= 4:
                assert gl(A) == 'vec%s(%s)' % (i,', '.join(str(s) for s in range(i)))
            else:
                assert gl(A) == 'float[%s](%s)' % (i,', '.join(str(s) for s in range(i)))

 def as_density_matrix(self):
     """
     NOT USED ANYMORE (but works fine)
     Represent this pure state as density matrix.
     """
     coeffs = self.get_coefficients_as_general_state()
     matrix = Matrix(coeffs)
     matrix = matrix * matrix.conjugate().transpose()
     coeff_comm = self.coeff_comm
     matrix = matrix * coeff_comm * coeff_comm
     return matrix

def msolve(args, f, x0, tol=None, maxsteps=None, verbose=False, norm=None,
           modules=['mpmath', 'sympy']):
    """
    Solves a nonlinear equation system numerically.

    f is a vector function of symbolic expressions representing the system.
    args are the variables.
    x0 is a starting vector close to a solution.

    Be careful with x0, not using floats might give unexpected results.

    Use modules to specify which modules should be used to evaluate the
    function and the Jacobian matrix. Make sure to use a module that supports
    matrices. For more information on the syntax, please see the docstring
    of lambdify.

    Currently only fully determined systems are supported.

    >>> from sympy import Symbol, Matrix
    >>> x1 = Symbol('x1')
    >>> x2 = Symbol('x2')
    >>> f1 = 3 * x1**2 - 2 * x2**2 - 1
    >>> f2 = x1**2 - 2 * x1 + x2**2 + 2 * x2 - 8
    >>> msolve((x1, x2), (f1, f2), (-1., 1.))
    [-1.19287309935246]
    [ 1.27844411169911]
    """
    if isinstance(f,  (list,  tuple)):
        f = Matrix(f).T
    if len(args) != f.cols:
        raise NotImplementedError('need exactly as many variables as equations')
    if verbose:
        print 'f(x):'
        print f
    # derive Jacobian
    J = f.jacobian(args)
    if verbose:
        print 'J(x):'
        print J
    # create functions
    f = lambdify(args, f.T, modules)
    J = lambdify(args, J, modules)
    # solve system using Newton's method
    kwargs = {}
    if tol:
        kwargs['tol'] = tol
    if maxsteps:
        kwargs['maxsteps'] = maxsteps
    kwargs['verbose'] = verbose
    if norm:
        kwargs['norm'] = norm
    x = newton(f, x0, J, **kwargs)
    return x

def vcv2corr(vcv):
    n = vcv.rows
    m = vcv.cols
    if n != m:
        raise ValueError('variance mattrix has to be square matrix')
    def var(i, j):
        if i == j:
            return vcv[i, j]**.5
        else:
            return 0
    D = Matrix(n, n, var)
    corr = D.inv() * vcv * D.inv()
    return corr, D

    def are_coplanar(*points):
        """

        This function tests whether passed points are coplanar or not.
        It uses the fact that the triple scalar product of three vectors
        vanishes if the vectors are coplanar. Which means that the volume
        of the solid described by them will have to be zero for coplanarity.

        Parameters
        ==========

        A set of points 3D points

        Returns
        =======

        boolean

        Examples
        ========

        >>> from sympy import Point3D
        >>> p1 = Point3D(1, 2, 2)
        >>> p2 = Point3D(2, 7, 2)
        >>> p3 = Point3D(0, 0, 2)
        >>> p4 = Point3D(1, 1, 2)
        >>> p5 = Point3D(1, 2, 2)
        >>> p1.are_coplanar(p2, p3, p4, p5)
        True
        >>> p6 = Point3D(0, 1, 3)
        >>> p1.are_coplanar(p2, p3, p4, p5, p6)
        False

        """
        if not all(isinstance(p, Point3D) for p in points):
            raise TypeError('Must pass only 3D Point objects')
        if(len(points) < 4):
            return True # These cases are always True
        points = list(set(points))
        for i in range(len(points) - 3):
            pv1 = [j - k for j, k in zip(points[i].args,   \
                points[i + 1].args)]
            pv2 = [j - k for j, k in zip(points[i + 1].args,
                points[i + 2].args)]
            pv3 = [j - k for j, k in zip(points[i + 2].args,
                points[i + 3].args)]
            pv1, pv2, pv3 = Matrix(pv1), Matrix(pv2), Matrix(pv3)
            stp = pv1.dot(pv2.cross(pv3))
            if stp != 0:
                return False
        return True

def test_converting_functions():
    arr_list = [1, 2, 3, 4]
    arr_matrix = Matrix(((1, 2), (3, 4)))

    # list
    arr_ndim_array = MutableDenseNDimArray(arr_list, (2, 2))
    assert (isinstance(arr_ndim_array, MutableDenseNDimArray))
    assert arr_matrix.tolist() == arr_ndim_array.tolist()

    # Matrix
    arr_ndim_array = MutableDenseNDimArray(arr_matrix)
    assert (isinstance(arr_ndim_array, MutableDenseNDimArray))
    assert arr_matrix.tolist() == arr_ndim_array.tolist()
    assert arr_matrix.shape == arr_ndim_array.shape

def verification_of_sufficient_second_order_conditions(K, df, x, l):
    if not K:
        print "The sufficient optimality condition of second kind is made."
    else:
        ddf = []
        for xi in x:
            ddf.append([ddfi.diff(xi) for ddfi in df])
        
        ddf = Matrix(ddf)
        ll = Matrix(l)
        if Matrix(ddf.dot(ll.T)).dot(ll).subs(K) > 0:
            print "The sufficient optimality condition of second kind is made."
        else:
            print "The sufficient optimality condition of second kind is not made."
            return

def demo_eigen_number():
    m = Matrix(3,3,[2,1,0,0,2,1,0,0,2])
    pprint(m)
    print "\nany n*n matrix has n complex eigen values, counted with multiplications\n"\
            "since the characteritic polynomial is n-dim\n"\
            "and Fundamental Theorem of Algebra gives us exactly n roots:"
    pprint(m.charpoly().as_expr())
    print"\neach distinct eigen value as at least 1 at most multiplication eigen vectors\n"\
            "but it's possible to have less than multiplication number of eigen vectors\n"\
            "the given matrix is an example with only 1 eigen vector:\n"
    pprint(m.eigenvects())
    print "to see this for any vector v we can compute M * v:\n"
    v = IndexedBase('v')
    vec = Matrix(m.rows, 1, lambda i, j: v[i+1])
    pprint(vec)
    pprint(m*vec)