def test_logm_type_preservation_and_conversion(self):
        # The logm matrix function should preserve the type of a matrix
        # whose eigenvalues are positive with zero imaginary part.
        # Test this preservation for variously structured matrices.
        complex_dtype_chars = ('F', 'D', 'G')
        for matrix_as_list in (
                [[1, 0], [0, 1]],
                [[1, 0], [1, 1]],
                [[2, 1], [1, 1]],
                [[2, 3], [1, 2]]):

            # check that the spectrum has the expected properties
            W = scipy.linalg.eigvals(matrix_as_list)
            assert_(not any(w.imag or w.real < 0 for w in W))

            # check float type preservation
            A = np.array(matrix_as_list, dtype=float)
            A_logm, info = logm(A, disp=False)
            assert_(A_logm.dtype.char not in complex_dtype_chars)

            # check complex type preservation
            A = np.array(matrix_as_list, dtype=complex)
            A_logm, info = logm(A, disp=False)
            assert_(A_logm.dtype.char in complex_dtype_chars)

            # check float->complex type conversion for the matrix negation
            A = -np.array(matrix_as_list, dtype=float)
            A_logm, info = logm(A, disp=False)
            assert_(A_logm.dtype.char in complex_dtype_chars)

 def relentropy(self,idealrho=None):
   # attention - log not working on pure matrices with 0-eigenvalues
   self.check_idealrho(idealrho)
   a = np.trace(np.dot(self.idealrho,splnlg.logm(self.idealrho)))
   b = np.trace(np.dot(self.idealrho,splnlg.logm(self.exprho)))
   self.relativeentropy=np.real(a-b)
   return self.relativeentropy

 def test_nils(self):
     a = array([[-2., 25., 0., 0., 0., 0., 0.],
                [0., -3., 10., 3., 3., 3., 0.],
                [0., 0., 2., 15., 3., 3., 0.],
                [0., 0., 0., 0., 15., 3., 0.],
                [0., 0., 0., 0., 3., 10., 0.],
                [0., 0., 0., 0., 0., -2., 25.],
                [0., 0., 0., 0., 0., 0., -3.]])
     m = (identity(7)*3.1+0j)-a
     logm(m, disp=False)

 def test_nils(self):
     a = array(
         [
             [-2.0, 25.0, 0.0, 0.0, 0.0, 0.0, 0.0],
             [0.0, -3.0, 10.0, 3.0, 3.0, 3.0, 0.0],
             [0.0, 0.0, 2.0, 15.0, 3.0, 3.0, 0.0],
             [0.0, 0.0, 0.0, 0.0, 15.0, 3.0, 0.0],
             [0.0, 0.0, 0.0, 0.0, 3.0, 10.0, 0.0],
             [0.0, 0.0, 0.0, 0.0, 0.0, -2.0, 25.0],
             [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -3.0],
         ]
     )
     m = (identity(7) * 3.1 + 0j) - a
     logm(m, disp=False)

 def test_opposite_sign_complex_eigenvalues(self):
     # See gh-6113
     E = [[0, 1], [-1, 0]]
     L = [[0, np.pi*0.5], [-np.pi*0.5, 0]]
     assert_allclose(expm(L), E, atol=1e-14)
     assert_allclose(logm(E), L, atol=1e-14)
     E = [[1j, 4], [0, -1j]]
     L = [[1j*np.pi*0.5, 2*np.pi], [0, -1j*np.pi*0.5]]
     assert_allclose(expm(L), E, atol=1e-14)
     assert_allclose(logm(E), L, atol=1e-14)
     E = [[1j, 0], [0, -1j]]
     L = [[1j*np.pi*0.5, 0], [0, -1j*np.pi*0.5]]
     assert_allclose(expm(L), E, atol=1e-14)
     assert_allclose(logm(E), L, atol=1e-14)

def vn_entropy_b(psi_t, label, nos, nol_b, nop):
    """
    Calculates Von-Neumann entropy as S = - tr(rho * ln(rho)).
    Also calculates trace of square of density matrix (measure of
    entanglement).
    Uses a filter to suppress 'WARNING: The logm input matrix may be nearly
    singular'. Wraps loop in tqdm for progress bar.

    :param psi_t: Psi(t)
    :param label: Relabelled states
    :param nos: No. of states
    :param nol_b: No. of lattice sites in B
    :param nop: No. of particles
    :return: Real Von-Neumann entropy
    :return: Trace of density matrix of B
    """
    vn_entropy = np.zeros(len(psi_t), dtype=complex)
    tr_sqr = np.zeros(len(psi_t), dtype=float)

    warnings.filterwarnings("ignore")

    for idx, vec in enumerate(psi_t):
        d_matrix = rho_b_pbasis(label, vec, nos, nol_b, nop)
        vn_entropy[idx] = -np.trace(np.dot(d_matrix, la.logm(d_matrix)))
        tr_sqr[idx] = trace_squared(d_matrix)

    return vn_entropy.real, tr_sqr

def d2c(Ad, Bd, C, D, dt):
    """returns A, B, C, D
    Converts a set of digital state space system matrices to their continuous counterpart.
    """
    A = la.logm(Ad)/dt
    B = la.solve((Ad - sp.eye(A.shape[0])), A) @ Bd
    return A, B, C, D

def test_expm():
    """Testing expm function"""
    m = np.random.rand(15, 15)
    m = m + m.T
    m_exp = my_mfd.expm(m)
    assert_array_almost_equal(m_exp.T, m_exp)
    assert_array_almost_equal(linalg.logm(m_exp), m)

    def __init__(self):
        app.Canvas.__init__(self, position=(50, 50), keys='interactive')

        self.program = gloo.Program(VERT_SHADER, FRAG_SHADER)

        self.program['a_position'] = gloo.VertexBuffer(v_position)
        self.program['a_color'] = gloo.VertexBuffer(v_color)
        self.program['a_size'] = gloo.VertexBuffer(v_size)

        self.program['u_pan'] = (0., 0.)
        self.program['u_scale'] = (1., 1.)

        self.program['u_vec1'] = (1., 0., 0., 0.)
        self.program['u_vec2'] = (0., 1., 0., 0.)

        # Circulant matrix.
        circ = np.diagflat(np.ones(ndim-1), 1)
        circ[-1, 0] = -1 if ndim % 2 == 0 else 1
        self.logcirc = logm(circ)
        # We will solve the equation dX/dt = log(circ) * X in real time
        # to compute the matrix exponential expm(t*log(circ)).
        self.mat = np.eye(ndim)
        self.dt = .001
        gloo.set_state(clear_color=(1, 1, 1, 1), blend=True,
                       blend_func=('src_alpha', 'one_minus_src_alpha'))
        self._timer = app.Timer('auto', connect=self.on_timer, start=True)

def karcher_mean(x, tol=0.01):
	'''
	Determined the Karcher mean of rotations
	Implementation from Algorithm 1, Rotation Averaging, Hartley et al, IJCV 2013
	'''
	R = x[0]
	N = x.shape[0]
	normDeltaR = np.inf
	itr = 0
	while True:
		#Estimate the delta rotation between the current center and all points
		deltaR  = np.zeros((3,3))
		oldNorm = normDeltaR
		for i in range(N):
			deltaR += linalg.logm(np.dot(np.transpose(R),x[i]))
		deltaR     = deltaR / N
		normDeltaR = linalg.norm(deltaR, ord='fro')/np.sqrt(2)

		if oldNorm - normDeltaR < tol:
			break
	
		R = np.dot(R, linalg.expm(deltaR)) 
		#print itr
		itr += 1		
	
	return R

 def test_al_mohy_higham_2012_experiment_1_logm(self):
     # The logm completes the round trip successfully.
     # Note that the expm leg of the round trip is badly conditioned.
     A = _get_al_mohy_higham_2012_experiment_1()
     A_logm, info = logm(A, disp=False)
     A_round_trip = expm(A_logm)
     assert_allclose(A_round_trip, A, rtol=1e-5, atol=1e-14)

def compute_Rt(P):
    """Gen a P matrix, calculate the R*t matrix."""
    try:
        Rt = matrix(sl.logm(P, disp = False)[0])
    except Exception as E:
        raise MatCalcExcep(str(E))

    return makeReal(Rt)

    def test_logm_type_conversion_mixed_sign_or_complex_spectrum(self):
        complex_dtype_chars = ("F", "D", "G")
        for matrix_as_list in ([[1, 0], [0, -1]], [[0, 1], [1, 0]], [[0, 1, 0], [0, 0, 1], [1, 0, 0]]):

            # check that the spectrum has the expected properties
            W = scipy.linalg.eigvals(matrix_as_list)
            assert_(any(w.imag or w.real < 0 for w in W))

            # check complex->complex
            A = np.array(matrix_as_list, dtype=complex)
            A_logm, info = logm(A, disp=False)
            assert_(A_logm.dtype.char in complex_dtype_chars)

            # check float->complex
            A = np.array(matrix_as_list, dtype=float)
            A_logm, info = logm(A, disp=False)
            assert_(A_logm.dtype.char in complex_dtype_chars)

 def test_multilog(self):
     A = np.zeros((self.k, self.m, self.m))
     l = np.zeros((self.k, self.m, self.m))
     for i in range(self.k):
         a = np.diag(rnd.rand(self.m))
         q, r = la.qr(rnd.randn(self.m, self.m))
         A[i] = q.dot(a.dot(q.T))
         l[i] = logm(A[i])
     np_testing.assert_allclose(multilog(A, pos_def=True), l)

 def test_round_trip_random_complex(self):
     np.random.seed(1234)
     for n in range(1, 6):
         M_unscaled = np.random.randn(n, n) + 1j * np.random.randn(n, n)
         for scale in np.logspace(-4, 4, 9):
             M = M_unscaled * scale
             M_logm, info = logm(M, disp=False)
             M_round_trip = expm(M_logm)
             assert_allclose(M_round_trip, M)

 def log(parameter):
     n_matrices = len(parameter) / 9
     matrices = numpy.empty((n_matrices, 4, 4))
     for i in xrange(n_matrices):
         p = parameter[i * 9: (i + 1) * 9].copy()
         matrices[i, :, :] = logm(
             affine_from_parameters(p)
         )
     return matrices

 def test_real_mixed_sign_spectrum(self):
     # These matrices have real eigenvalues with mixed signs.
     # The output logm dtype is complex, regardless of input dtype.
     for M in (
             [[1, 0], [0, -1]],
             [[0, 1], [1, 0]]):
         for dt in float, complex:
             A = np.array(M, dtype=dt)
             A_logm, info = logm(A, disp=False)
             assert_(np.issubdtype(A_logm.dtype, np.complexfloating))

 def test_logm_nearly_singular(self):
     M = np.array([[1e-100]])
     expected_warning = _matfuncs_inv_ssq.LogmNearlySingularWarning
     with warnings.catch_warnings(record=True) as w:
         warnings.simplefilter('always')
         L, info = logm(M, disp=False)
         assert_equal(len(w), 1)
         assert_(issubclass(w[-1].category, expected_warning))
         E = expm(L)
         assert_allclose(E, M, atol=1e-14)

 def test_logm_consistency(self):
     random.seed(1234)
     for dtype in [np.float64, np.complex128]:
         for n in range(1, 10):
             for scale in [1e-4, 1e-3, 1e-2, 1e-1, 1, 1e1, 1e2]:
                 # make logm(A) be of a given scale
                 A = (eye(n) + random.rand(n, n) * scale).astype(dtype)
                 if np.iscomplexobj(A):
                     A = A + 1j * random.rand(n, n) * scale
                 assert_array_almost_equal(expm(logm(A)), A)

 def test_complex_spectrum_real_logm(self):
     # This matrix has complex eigenvalues and real logm.
     # Its output dtype depends on its input dtype.
     M = [[1, 1, 2], [2, 1, 1], [1, 2, 1]]
     for dt in float, complex:
         X = np.array(M, dtype=dt)
         w = scipy.linalg.eigvals(X)
         assert_(1e-2 < np.absolute(w.imag).sum())
         Y, info = logm(X, disp=False)
         assert_(np.issubdtype(Y.dtype, np.inexact))
         assert_allclose(expm(Y), X)