Python scipy.linalg.eigvals() Examples

def testMinrealSS(self):
        """Test a minreal model reduction"""
        #A = [-2, 0.5, 0; 0.5, -0.3, 0; 0, 0, -0.1]
        A = [[-2, 0.5, 0], [0.5, -0.3, 0], [0, 0, -0.1]]
        #B = [0.3, -1.3; 0.1, 0; 1, 0]
        B = [[0.3, -1.3], [0.1, 0.], [1.0, 0.0]]
        #C = [0, 0.1, 0; -0.3, -0.2, 0]
        C = [[0., 0.1, 0.0], [-0.3, -0.2, 0.0]]
        #D = [0 -0.8; -0.3 0]
        D = [[0., -0.8], [-0.3, 0.]]
        # sys = ss(A, B, C, D)

        sys = StateSpace(A, B, C, D)
        sysr = sys.minreal()
        self.assertEqual(sysr.states, 2)
        self.assertEqual(sysr.inputs, sys.inputs)
        self.assertEqual(sysr.outputs, sys.outputs)
        np.testing.assert_array_almost_equal(
            eigvals(sysr.A), [-2.136154, -0.1638459]) 

def check_stability(A, dt=True):
    """
    Checks the stability of the system.

    Args:
        A (np.ndarray): System plant matrix
        dt (bool): Discrete time system

    Returns:
        bool: True if the system is stable
    """
    eigvals = scalg.eigvals(A)
    if dt:
        criteria = np.abs(eigvals) > 1.
    else:
        criteria = np.real(eigvals) > 0.0

    if np.sum(criteria) >= 1.0:
        return True
    else:
        return False 

def test_dare_g(self):
        A = array([[-0.6, 0],[-0.1, -0.4]])
        Q = array([[2, 1],[1, 3]])
        B = array([[1, 5],[2, 4]])
        R = array([[1, 0],[0, 1]])
        S = array([[1, 0],[2, 0]])
        E = array([[2, 1],[1, 2]])

        X,L,G = dare(A,B,Q,R,S,E)
        # print("The solution obtained is", X)
        Gref = solve(B.T.dot(X).dot(B) + R, B.T.dot(X).dot(A) + S.T)
        assert_array_almost_equal(Gref,G)
        assert_array_almost_equal(
            A.T.dot(X).dot(A) - E.T.dot(X).dot(E)
            - (A.T.dot(X).dot(B) + S).dot(Gref) + Q,
            zeros((2,2)) )
        # check for stable closed loop
        lam = eigvals(A - B.dot(G), E)
        assert_array_less(abs(lam), 1.0)

        A = array([[-0.6, 0],[-0.1, -0.4]])
        Q = array([[2, 1],[1, 3]])
        B = array([[1],[2]])
        R = 1
        S = array([[1],[2]])
        E = array([[2, 1],[1, 2]])

        X,L,G = dare(A,B,Q,R,S,E)
        # print("The solution obtained is", X)
        assert_array_almost_equal(
            A.T.dot(X).dot(A) - E.T.dot(X).dot(E) -
            (A.T.dot(X).dot(B) + S).dot(solve(B.T.dot(X).dot(B) + R, B.T.dot(X).dot(A) + S.T)) + Q,
            zeros((2,2)) )
        assert_array_almost_equal((B.T.dot(X).dot(A) + S.T) / (B.T.dot(X).dot(B) + R), G)
        # check for stable closed loop
        lam = eigvals(A - B.dot(G), E)
        assert_array_less(abs(lam), 1.0) 

def modularity_spectrum(G):
    """Return eigenvalues of the modularity matrix of G.

    Parameters
    ----------
    G : Graph
       A NetworkX Graph or DiGraph

    Returns
    -------
    evals : NumPy array
      Eigenvalues

    See Also
    --------
    modularity_matrix

    References
    ----------
    .. [1] M. E. J. Newman, "Modularity and community structure in networks",
       Proc. Natl. Acad. Sci. USA, vol. 103, pp. 8577-8582, 2006.
    """
    from scipy.linalg import eigvals
    if G.is_directed():
        return eigvals(nx.directed_modularity_matrix(G))
    else:
        return eigvals(nx.modularity_matrix(G))

# fixture for nose tests 

def adjacency_spectrum(G, weight='weight'):
    """Return eigenvalues of the adjacency matrix of G.

    Parameters
    ----------
    G : graph
       A NetworkX graph

    weight : string or None, optional (default='weight')
       The edge data key used to compute each value in the matrix.
       If None, then each edge has weight 1.

    Returns
    -------
    evals : NumPy array
      Eigenvalues

    Notes
    -----
    For MultiGraph/MultiDiGraph, the edges weights are summed.
    See to_numpy_matrix for other options.

    See Also
    --------
    adjacency_matrix
    """
    from scipy.linalg import eigvals
    return eigvals(nx.adjacency_matrix(G, weight=weight).todense()) 

def modularity_spectrum(G):
    """Returns eigenvalues of the modularity matrix of G.

    Parameters
    ----------
    G : Graph
       A NetworkX Graph or DiGraph

    Returns
    -------
    evals : NumPy array
      Eigenvalues

    See Also
    --------
    modularity_matrix

    References
    ----------
    .. [1] M. E. J. Newman, "Modularity and community structure in networks",
       Proc. Natl. Acad. Sci. USA, vol. 103, pp. 8577-8582, 2006.
    """
    from scipy.linalg import eigvals
    if G.is_directed():
        return eigvals(nx.directed_modularity_matrix(G))
    else:
        return eigvals(nx.modularity_matrix(G))

# fixture for nose tests 

def adjacency_spectrum(G, weight='weight'):
    """Returns eigenvalues of the adjacency matrix of G.

    Parameters
    ----------
    G : graph
       A NetworkX graph

    weight : string or None, optional (default='weight')
       The edge data key used to compute each value in the matrix.
       If None, then each edge has weight 1.

    Returns
    -------
    evals : NumPy array
      Eigenvalues

    Notes
    -----
    For MultiGraph/MultiDiGraph, the edges weights are summed.
    See to_numpy_matrix for other options.

    See Also
    --------
    adjacency_matrix
    """
    from scipy.linalg import eigvals
    return eigvals(nx.adjacency_matrix(G, weight=weight).todense()) 

def z2_vanderbilt(h,nk=30,nt=100,nocc=None,full=False):
  """ Calculate Z2 invariant according to Vanderbilt algorithm"""
  out = [] # output list
  path = np.linspace(0.,1.,nk) # set of kpoints
  fo = open("WANNIER_CENTERS.OUT","w")
  if full:  ts = np.linspace(0.,1.0,nt,endpoint=False)
  else:  ts = np.linspace(0.,0.5,nt,endpoint=False)
  wfall = [[occ_states2d(h,np.array([k,t,0.,])) for k in path] for t in ts] 
  # select a continuos gauge for the first wave
  for it in range(len(ts)-1): # loop over ts
    wfall[it+1][0] = smooth_gauge(wfall[it][0],wfall[it+1][0]) 
  for it in range(len(ts)): # loop over t points
    row = [] # empty list for this row
    t = ts[it] # select the t point
    wfs = wfall[it] # get set of waves 
    for i in range(len(wfs)-1):
      wfs[i+1] = smooth_gauge(wfs[i],wfs[i+1]) # transform into a smooth gauge
#      m = uij(wfs[i],wfs[i+1]) # matrix of wavefunctions
    m = uij(wfs[0],wfs[len(wfs)-1]) # matrix of wavefunctions
    evals = lg.eigvals(m) # eigenvalues of the rotation 
    x = np.angle(evals) # phase of the eigenvalues
    fo.write(str(t)+"    ") # write pumping variable
    row.append(t) # store
    for ix in x: # loop over phases
      fo.write(str(ix)+"  ")
      row.append(ix) # store
    fo.write("\n")
    out.append(row) # store
  fo.close()
  return np.array(out).transpose() # transpose the map 

def det_linear_buckling(system):
    """
    Determine linear buckling by solving the generalized eigenvalue problem (k -λkg)x = 0.

    geometrical stiffness matrix at buckling point: Kg = f(N_max)
    1st order forces: N0
    Nmax = λN0
    Kg(Nmax) = λ(Kg(N0) = λKg0

    2nd order analysis is solved by:
    (K + λKg0)U = F

    We are interested in the point that there is nog additional load F and displacement U is possible.
    (K + λKg0)ΔU = ΔF = 0
    (K + λKg0) = 0

    Is the generalized eigenvalue problem:
    (A - λB)x = 0

    :param system: (SystemElements)
    :return: (flt) The factor the loads can be increased until the structure fails due to buckling.
    """
    system.solve()

    # buckling
    k0 = system.reduced_system_matrix * 1.0  # copy

    for el in system.element_map.values():
        el.compile_geometric_non_linear_stiffness_matrix()
        el.reset()

    system.solve()
    kg = system.reduced_system_matrix - k0
    # solve (k -λkg)x = 0

    eigenvalues = np.abs(linalg.eigvals(k0, kg))
    return np.min(eigenvalues) 

def modularity_spectrum(G):
    """Return eigenvalues of the modularity matrix of G.

    Parameters
    ----------
    G : Graph
       A NetworkX Graph or DiGraph

    Returns
    -------
    evals : NumPy array
      Eigenvalues

    See Also
    --------
    modularity_matrix

    References
    ----------
    .. [1] M. E. J. Newman, "Modularity and community structure in networks",
       Proc. Natl. Acad. Sci. USA, vol. 103, pp. 8577-8582, 2006.
    """
    from scipy.linalg import eigvals
    if G.is_directed():
        return eigvals(nx.directed_modularity_matrix(G))
    else:
        return eigvals(nx.modularity_matrix(G))

# fixture for nose tests 

def adjacency_spectrum(G, weight='weight'):
    """Return eigenvalues of the adjacency matrix of G.

    Parameters
    ----------
    G : graph
       A NetworkX graph

    weight : string or None, optional (default='weight')
       The edge data key used to compute each value in the matrix.
       If None, then each edge has weight 1.

    Returns
    -------
    evals : NumPy array
      Eigenvalues

    Notes
    -----
    For MultiGraph/MultiDiGraph, the edges weights are summed.
    See to_numpy_matrix for other options.

    See Also
    --------
    adjacency_matrix
    """
    from scipy.linalg import eigvals
    return eigvals(nx.adjacency_matrix(G,weight=weight).todense()) 

def run_hamiltonian(hessian, verbose=True):
    c = ClassicalHamiltonian()

    xopt = optimize.fmin(c.potential, c.initialposition(), xtol=1e-10)

    hessian.fun = c.potential
    hessian.full_output = True

    h, info = hessian(xopt)
    true_h = np.array([[5.23748385e-12, -2.61873829e-12],
                       [-2.61873829e-12, 5.23748385e-12]])
    eigenvalues = linalg.eigvals(h)
    normal_modes = c.normal_modes(eigenvalues)

    if verbose:
        print(c.potential([-0.5, 0.5]))
        print(c.potential([-0.5, 0.0]))
        print(c.potential([0.0, 0.0]))
        print(xopt)
        print('h', h)
        print('h-true_h', np.abs(h - true_h))
        print('error_estimate', info.error_estimate)

        print('eigenvalues', eigenvalues)
        print('normal_modes', normal_modes)
    return h, info.error_estimate, true_h 

def _default_response_times(A, n):
    """Compute a reasonable set of time samples for the response time.

    This function is used by `impulse`, `impulse2`, `step` and `step2`
    to compute the response time when the `T` argument to the function
    is None.

    Parameters
    ----------
    A : array_like
        The system matrix, which is square.
    n : int
        The number of time samples to generate.

    Returns
    -------
    t : ndarray
        The 1-D array of length `n` of time samples at which the response
        is to be computed.
    """
    # Create a reasonable time interval.
    # TODO: This could use some more work.
    # For example, what is expected when the system is unstable?
    vals = linalg.eigvals(A)
    r = min(abs(real(vals)))
    if r == 0.0:
        r = 1.0
    tc = 1.0 / r
    t = linspace(0.0, 7 * tc, n)
    return t 

def entropy(state, base=2):
    r"""Calculate the von-Neumann entropy of a quantum state.

    The entropy :math:`S` is given by

    .. math:

        S(\rho) = - Tr[\rho \log(\rho)]

    Args:
        state (Statevector or DensityMatrix): a quantum state.
        base (int): the base of the logarithm [Default: 2].

    Returns:
        float: The von-Neumann entropy S(rho).

    Raises:
        QiskitError: if the input state is not a valid QuantumState.
    """
    # pylint: disable=assignment-from-no-return
    state = _format_state(state, validate=True)
    if isinstance(state, Statevector):
        return 0
    # Density matrix case
    evals = np.maximum(np.real(la.eigvals(state.data)), 0.)
    return shannon_entropy(evals, base=base) 

def _default_response_times(A, n):
    """Compute a reasonable set of time samples for the response time.

    This function is used by `impulse`, `impulse2`, `step` and `step2`
    to compute the response time when the `T` argument to the function
    is None.

    Parameters
    ----------
    A : array_like
        The system matrix, which is square.
    n : int
        The number of time samples to generate.

    Returns
    -------
    t : ndarray
        The 1-D array of length `n` of time samples at which the response
        is to be computed.
    """
    # Create a reasonable time interval.
    # TODO: This could use some more work.
    # For example, what is expected when the system is unstable?
    vals = linalg.eigvals(A)
    r = min(abs(real(vals)))
    if r == 0.0:
        r = 1.0
    tc = 1.0 / r
    t = linspace(0.0, 7 * tc, n)
    return t 

def test_ackermann_controllable():
    #
    A = haroldcompanion([1, 6, 5, 1])
    B = eye(3)[:, [-1]]
    p = [-10, -9, -8]
    K = ackermann((A, B), p)
    pa = eigvals(A - B@K)
    assert_array_almost_equal(array(p, dtype=complex), sort(pa)) 

def bench_eigvals():
    numpy_eigvals = nl.eigvals
    scipy_eigvals = sl.eigvals
    print()
    print('           Finding matrix eigenvalues')
    print('      ==================================')
    print('      |    contiguous     |   non-contiguous ')
    print('----------------------------------------------')
    print(' size |  scipy  | numpy   |  scipy  | numpy ')

    for size,repeat in [(20,150),(100,7),(200,2)]:
        repeat *= 1
        print('%5s' % size, end=' ')
        sys.stdout.flush()

        a = random([size,size])

        print('| %6.2f ' % measure('scipy_eigvals(a)',repeat), end=' ')
        sys.stdout.flush()

        print('| %6.2f ' % measure('numpy_eigvals(a)',repeat), end=' ')
        sys.stdout.flush()

        a = a[-1::-1,-1::-1]  # turn into a non-contiguous array
        assert_(not a.flags['CONTIGUOUS'])

        print('| %6.2f ' % measure('scipy_eigvals(a)',repeat), end=' ')
        sys.stdout.flush()

        print('| %6.2f ' % measure('numpy_eigvals(a)',repeat), end=' ')
        sys.stdout.flush()

        print('   (secs for %s calls)' % (repeat)) 

def _default_response_times(A, n):
    """Compute a reasonable set of time samples for the response time.

    This function is used by `impulse`, `impulse2`, `step` and `step2`
    to compute the response time when the `T` argument to the function
    is None.

    Parameters
    ----------
    A : ndarray
        The system matrix, which is square.
    n : int
        The number of time samples to generate.

    Returns
    -------
    t : ndarray
        The 1-D array of length `n` of time samples at which the response
        is to be computed.
    """
    # Create a reasonable time interval.
    # TODO: This could use some more work.
    # For example, what is expected when the system is unstable?
    vals = linalg.eigvals(A)
    r = min(abs(real(vals)))
    if r == 0.0:
        r = 1.0
    tc = 1.0 / r
    t = linspace(0.0, 7 * tc, n)
    return t 

def _default_response_times(A, n):
    """Compute a reasonable set of time samples for the response time.

    This function is used by `impulse`, `impulse2`, `step` and `step2`
    to compute the response time when the `T` argument to the function
    is None.

    Parameters
    ----------
    A : array_like
        The system matrix, which is square.
    n : int
        The number of time samples to generate.

    Returns
    -------
    t : ndarray
        The 1-D array of length `n` of time samples at which the response
        is to be computed.
    """
    # Create a reasonable time interval.
    # TODO: This could use some more work.
    # For example, what is expected when the system is unstable?
    vals = linalg.eigvals(A)
    r = min(abs(real(vals)))
    if r == 0.0:
        r = 1.0
    tc = 1.0 / r
    t = linspace(0.0, 7 * tc, n)
    return t 

def dare(A, B, Q, R, S=None, E=None, stabilizing=True):
    """ (X,L,G) = dare(A,B,Q,R) solves the discrete-time algebraic Riccati
    equation

        :math:`A^T X A - X - A^T X B (B^T X B + R)^{-1} B^T X A + Q = 0`

    where A and Q are square matrices of the same dimension. Further, Q
    is a symmetric matrix. The function returns the solution X, the gain
    matrix G = (B^T X B + R)^-1 B^T X A and the closed loop eigenvalues L,
    i.e., the eigenvalues of A - B G.

    (X,L,G) = dare(A,B,Q,R,S,E) solves the generalized discrete-time algebraic
    Riccati equation

        :math:`A^T X A - E^T X E - (A^T X B + S) (B^T X B + R)^{-1} (B^T X A + S^T) + Q = 0`

    where A, Q and E are square matrices of the same dimension. Further, Q and
    R are symmetric matrices. The function returns the solution X, the gain
    matrix :math:`G = (B^T X B + R)^{-1} (B^T X A + S^T)` and the closed loop
    eigenvalues L, i.e., the eigenvalues of A - B G , E.
    """
    if S is not None or E is not None or not stabilizing:
        return dare_old(A, B, Q, R, S, E, stabilizing)
    else:
        Rmat = _ssmatrix(R)
        Qmat = _ssmatrix(Q)
        X = solve_discrete_are(A, B, Qmat, Rmat)
        G = solve(B.T.dot(X).dot(B) + Rmat, B.T.dot(X).dot(A))
        L = eigvals(A - B.dot(G))
        return _ssmatrix(X), L, _ssmatrix(G) 

def testMinrealBrute(self):
        for n, m, p in permutations(range(1,6), 3):
            s = matlab.rss(n, p, m)
            sr = s.minreal()
            if s.states > sr.states:
                self.nreductions += 1
            else:
                # Check to make sure that poles and zeros match

                # For poles, just look at eigenvalues of A
                np.testing.assert_array_almost_equal(
                    np.sort(eigvals(s.A)), np.sort(eigvals(sr.A)))

                # For zeros, need to extract SISO systems
                for i in range(m):
                    for j in range(p):
                        # Extract SISO dynamixs from input i to output j
                        s1 = matlab.ss(s.A, s.B[:,i], s.C[j,:], s.D[j,i])
                        s2 = matlab.ss(sr.A, sr.B[:,i], sr.C[j,:], sr.D[j,i])

                        # Check that the zeros match
                        # Note: sorting doesn't work => have to do the hard way
                        z1 = matlab.zero(s1)
                        z2 = matlab.zero(s2)

                        # Start by making sure we have the same # of zeros
                        self.assertEqual(len(z1), len(z2))

                        # Make sure all zeros in s1 are in s2
                        for zero in z1:
                            # Find the closest zero
                            self.assertAlmostEqual(min(abs(z2 - zero)), 0.)

                        # Make sure all zeros in s2 are in s1
                        for zero in z2:
                            # Find the closest zero
                            self.assertAlmostEqual(min(abs(z1 - zero)), 0.)

        # Make sure that the number of systems reduced is as expected
        # (Need to update this number if you change the seed at top of file)
        self.assertEqual(self.nreductions, 2) 

def test_dare(self):
        A = array([[-0.6, 0],[-0.1, -0.4]])
        Q = array([[2, 1],[1, 0]])
        B = array([[2, 1],[0, 1]])
        R = array([[1, 0],[0, 1]])

        X,L,G = dare(A,B,Q,R)
        # print("The solution obtained is", X)
        Gref = solve(B.T.dot(X).dot(B) + R, B.T.dot(X).dot(A))
        assert_array_almost_equal(Gref, G)
        assert_array_almost_equal(
            A.T.dot(X).dot(A) - X -
            A.T.dot(X).dot(B).dot(Gref) + Q,
            zeros((2,2)))
        # check for stable closed loop
        lam = eigvals(A - B.dot(G))
        assert_array_less(abs(lam), 1.0)

        A = array([[1, 0],[-1, 1]])
        Q = array([[0, 1],[1, 1]])
        B = array([[1],[0]])
        R = 2

        X,L,G = dare(A,B,Q,R)
        # print("The solution obtained is", X)
        assert_array_almost_equal(
            A.T.dot(X).dot(A) - X -
            A.T.dot(X).dot(B) * solve(B.T.dot(X).dot(B) + R, B.T.dot(X).dot(A)) + Q, zeros((2,2)))
        assert_array_almost_equal(B.T.dot(X).dot(A) / (B.T.dot(X).dot(B) + R), G)
        # check for stable closed loop
        lam = eigvals(A - B.dot(G))
        assert_array_less(abs(lam), 1.0) 

def roots(self):
        """
        Utilises Boyd's O(n^2) recursive subdivision algorithm. The chebfun
        is recursively subsampled until it is successfully represented to
        machine precision by a sequence of piecewise interpolants of degree
        100 or less. A colleague matrix eigenvalue solve is then applied to
        each of these pieces and the results are concatenated.
        See:
        J. P. Boyd, Computing zeros on a real interval through Chebyshev
        expansion and polynomial rootfinding, SIAM J. Numer. Anal., 40
        (2002), pp. 1666–1682.
        """
        if self.size() == 1:
            return np.array([])

        elif self.size() <= 100:
            ak = self.coefficients()
            v = np.zeros_like(ak[:-1])
            v[1] = 0.5
            C1 = linalg.toeplitz(v)
            C2 = np.zeros_like(C1)
            C1[0,1] = 1.
            C2[-1,:] = ak[:-1]
            C = C1 - .5/ak[-1] * C2
            eigenvalues = linalg.eigvals(C)
            roots = [eig.real for eig in eigenvalues
                    if np.allclose(eig.imag,0,atol=1e-10)
                        and np.abs(eig.real) <=1]
            scaled_roots = self._ui_to_ab(np.array(roots))
            return scaled_roots
        else:
            try:
                # divide at a close-to-zero split-point
                split_point = self._ui_to_ab(0.0123456789)
                return np.concatenate(
                    (self.restrict([self._domain[0],split_point]).roots(),
                     self.restrict([split_point,self._domain[1]]).roots()))
            except:
                # Seems to have many fake roots for high degree fits
                coeffs = self.coefficients()
                domain = self._domain
                possibilities =  Chebyshev(coeffs, domain).roots()
                return np.array([float(i.real) for i in possibilities if i.imag == 0.0])

    # ----------------------------------------------------------------
    # Interpolation and evaluation (go from values to coefficients)
    # ---------------------------------------------------------------- 

def concurrence(state):
    r"""Calculate the concurrence of a quantum state.

    The concurrence of a bipartite
    :class:`~qiskit.quantum_info.Statevector` :math:`|\psi\rangle` is
    given by

    .. math:

        C(|\psi\rangle) = \sqrt{2(1 - Tr[\rho_0^2])}

    where :math:`\rho_0 = Tr_1[|\psi\rangle\!\langle\psi|]` is the
    reduced state from by taking the
    :math:`~qiskit.quantum_info.partial_trace` of the input state.

    For density matrices the concurrence is only defined for
    2-qubit states, it is given by:

    .. math:

        C(\rho) = \max(0, \lambda_1 - \lambda_2 - \lamda_3 - \lambda_4)

    where  :math:`\lambda _1 \ge \lambda _2 \ge \lambda _3 \ge \lambda _4`
    are the ordered eigenvalues of the matrix
    :math:`R=\sqrt{\sqrt{\rho }(Y\otimes Y)\overline{\rho}(Y\otimes Y)\sqrt{\rho}}}`.

    Args:
        state (Statevector or DensityMatrix): a 2-qubit quantum state.

    Returns:
        float: The concurrence.

    Raises:
        QiskitError: if the input state is not a valid QuantumState.
        QiskitError: if input is not a bipartite QuantumState.
        QiskitError: if density matrix input is not a 2-qubit state.
    """
    # Concurrence computation requires the state to be valid
    state = _format_state(state, validate=True)
    if isinstance(state, Statevector):
        # Pure state concurrence
        dims = state.dims()
        if len(dims) != 2:
            raise QiskitError('Input is not a bipartite quantum state.')
        qargs = [0] if dims[0] > dims[1] else [1]
        rho = partial_trace(state, qargs)
        return float(np.sqrt(2 * (1 - np.real(purity(rho)))))
    # If input is a density matrix it must be a 2-qubit state
    if state.dim != 4:
        raise QiskitError('Input density matrix must be a 2-qubit state.')
    rho = DensityMatrix(state).data
    yy_mat = np.fliplr(np.diag([-1, 1, 1, -1]))
    sigma = rho.dot(yy_mat).dot(rho.conj()).dot(yy_mat)
    w = np.sort(np.real(la.eigvals(sigma)))
    w = np.sqrt(np.maximum(w, 0.))
    return max(0.0, w[-1] - np.sum(w[0:-1])) 

def spherical_integral(C,rho):
    """
    Calculate the integral of a function over a unit sphere.
    """
    # phi - azimuthal angle (angle in xy-plane)
    # theta - polar angle (angle between z and xy-plane)
    #       (  y  , x )
    def func(theta,phi,C,rho):  # Test function. Can I get 4*pi^2????
        x = sp.cos(phi)*sp.sin(theta)
        y = sp.sin(phi)*sp.sin(theta)
        z = sp.cos(theta)
        #dir = sp.array((x,y,z))
        #dc = dir_cosines(dir)
        dc = sp.array((x,y,z))  # Turns out these are direction cosines!
        Gamma = make_gamma(dc,C)
        rho_c_square = linalg.eigvals(Gamma).real  # GPa
        rho_c_square = rho_c_square*1e9  # Pa
        sound_vel = sp.sqrt(rho_c_square/rho)  # m/s
        integrand = 1/(sound_vel[0]**3) + 1/(sound_vel[1]**3) + 1/(sound_vel[2]**3)
        return integrand*sp.sin(theta)
    #        (  y  , x      )
    #def sfunc(theta,phi,args=()):
    #    return func(theta,phi,args)*sp.sin(theta)

    integral,error = dblquad(func,0,2*sp.pi,lambda g: 0,lambda h:
            sp.pi,args=(C,rho))
    return integral


#direction = sp.array((1.0,1.0,1.0))
#dc = dir_cosines(direction)
#C = read_file.read_file(sys.argv[1])
#C.pop(0)
#C = sp.array(C,float)
#Gamma = make_gamma(dc,C)
#density = 7500 #kg/m**3
#density = float(read_file.read_file(sys.argv[2])[0][0])
#rho_c_square = linalg.eigvals(Gamma) #GPa
#rho_c_square = rho_c_square*1e9 #Pa
#sound_vel = sp.sqrt(rho_c_square/density).real
#avg_vel = sp.average(sound_vel)
#print Gamma
#print direction
#print C
#print rho_c_square
#print rho_c_square.real
#print sound_vel," in m/s"
#print avg_vel
#print spherical_integral(C,density)