Python linalg.logm方法代码示例

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import logm [as 别名]
def test_round_trip_random_float(self):
        np.random.seed(1234)
        for n in range(1, 6):
            M_unscaled = np.random.randn(n, n)
            for scale in np.logspace(-4, 4, 9):
                M = M_unscaled * scale

                # Eigenvalues are related to the branch cut.
                W = np.linalg.eigvals(M)
                err_msg = 'M:{0} eivals:{1}'.format(M, W)

                # Check sqrtm round trip because it is used within logm.
                M_sqrtm, info = sqrtm(M, disp=False)
                M_sqrtm_round_trip = M_sqrtm.dot(M_sqrtm)
                assert_allclose(M_sqrtm_round_trip, M)

                # Check logm round trip.
                M_logm, info = logm(M, disp=False)
                M_logm_round_trip = expm(M_logm)
                assert_allclose(M_logm_round_trip, M, err_msg=err_msg) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import logm [as 别名]
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) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import logm [as 别名]
def test_random_matrices_and_powers(self):
        # Each independent iteration of this fuzz test picks random parameters.
        # It tries to hit some edge cases.
        np.random.seed(1234)
        nsamples = 20
        for i in range(nsamples):
            # Sample a matrix size and a random real power.
            n = random.randrange(1, 5)
            p = np.random.randn()

            # Sample a random real or complex matrix.
            matrix_scale = np.exp(random.randrange(-4, 5))
            A = np.random.randn(n, n)
            if random.choice((True, False)):
                A = A + 1j * np.random.randn(n, n)
            A = A * matrix_scale

            # Check a couple of analytically equivalent ways
            # to compute the fractional matrix power.
            # These can be compared because they both use the principal branch.
            A_power = fractional_matrix_power(A, p)
            A_logm, info = logm(A, disp=False)
            A_power_expm_logm = expm(A_logm * p)
            assert_allclose(A_power, A_power_expm_logm) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import logm [as 别名]
def compute_velocity_from_msg(self, P0, P1):
        p0, q0, t0 = self.p_q_t_from_msg(P0)
        p1, q1, t1 = self.p_q_t_from_msg(P1)

        # There's something wrong with the current function to go from quat to matrix.
        # Using the TF version instead.
        q0_ros = [q0.x, q0.y, q0.z, q0.w]
        q1_ros = [q1.x, q1.y, q1.z, q1.w]
        
        import tf
        H0 = tf.transformations.quaternion_matrix(q0_ros)
        H0[:3, 3] = p0

        H1 = tf.transformations.quaternion_matrix(q1_ros)
        H1[:3, 3] = p1

        # Let the homogeneous matrix handle the inversion etc. Guaranteed correctness.
        H01 = np.dot(np.linalg.inv(H0), H1)
        dt = t1 - t0

        V = H01[:3, 3] / dt
        w_hat = logm(H01[:3, :3]) / dt
        Omega = np.array([w_hat[2,1], w_hat[0,2], w_hat[1,0]])

        return V, Omega, dt 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import logm [as 别名]
def compute_velocity_from_msg(self, P0, P1):
        p0, q0, t0 = self.p_q_t(P0)
        p1, q1, t1 = self.p_q_t(P1)

        # There's something wrong with the current function to go from quat to matrix.
        # Using the TF version instead.
        q0_ros = [q0.x, q0.y, q0.z, q0.w]
        q1_ros = [q1.x, q1.y, q1.z, q1.w]
        
        import tf
        H0 = tf.transformations.quaternion_matrix(q0_ros)
        H0[:3, 3] = p0

        H1 = tf.transformations.quaternion_matrix(q1_ros)
        H1[:3, 3] = p1

        # Let the homogeneous matrix handle the inversion etc. Guaranteed correctness.
        H01 = np.dot(np.linalg.inv(H0), H1)
        dt = t1 - t0

        V = H01[:3, 3] / dt
        w_hat = logm(H01[:3, :3]) / dt
        Omega = np.array([w_hat[2,1], w_hat[0,2], w_hat[1,0]])

        return V, Omega, dt 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import logm [as 别名]
def exact_matrix_logarithm_trace(Fx, x):
    """
    Computes slow-ass Tr(Ln(d(Fx)/dx))
    :param Fx: output of f(x)
    :param x: input
    :return: Tr(Ln(I + df/dx))
    """
    bs = Fx.size(0)
    outVector = torch.sum(Fx, 0).view(-1)
    outdim = outVector.size()[0]
    indim = x.view(bs, -1).size()
    jac = torch.empty([bs, outdim, indim[1]], dtype=torch.float)
    # for each output Fx[i] compute d(Fx[i])/d(x)
    for i in range(outdim):
        zero_gradients(x)
        jac[:, i, :] = torch.autograd.grad(outVector[i], x,
                                           retain_graph=True)[0].view(bs, outdim)
    jac = jac.cpu().numpy()
    iden = np.eye(jac.shape[1])
    log_jac = np.stack([logm(jac[i] + iden) for i in range(bs)])
    trace_jac = np.diagonal(log_jac, axis1=1, axis2=2).sum(1)
    return trace_jac 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import logm [as 别名]
def unitary_superoperator_matrix_log(M, mxBasis):
    """
    Construct the logarithm of superoperator matrix `M`
    that acts as a unitary on density-matrix space,
    (`M: rho -> U rho Udagger`) so that log(M) can be
    written as the action by Hamiltonian `H`:
    `log(M): rho -> -i[H,rho]`.


    Parameters
    ----------
    M : numpy array
        The superoperator matrix whose logarithm is taken

    mxBasis : {'std', 'gm', 'pp', 'qt'} or Basis object
        The source and destination basis, respectively.  Allowed
        values are Matrix-unit (std), Gell-Mann (gm), Pauli-product (pp),
        and Qutrit (qt) (or a custom basis object).

    Returns
    -------
    numpy array
        A matrix `logM`, of the same shape as `M`, such that `M = exp(logM)`
        and `logM` can be written as the action `rho -> -i[H,rho]`.
    """
    from . import lindbladtools as _lt  # (would create circular imports if at top)
    from . import optools as _gt  # (would create circular imports if at top)

    M_std = change_basis(M, mxBasis, "std")
    evals = _np.linalg.eigvals(M_std)
    assert(_np.allclose(_np.abs(evals), 1.0))  # simple but technically incomplete check for a unitary superop
    # (e.g. could be anti-unitary: diag(1, -1, -1, -1))
    U = _gt.process_mx_to_unitary(M_std)
    H = _spl.logm(U) / -1j  # U = exp(-iH)
    logM_std = _lt.hamiltonian_to_lindbladian(H)  # rho --> -i[H, rho] * sqrt(d)/2
    logM = change_basis(logM_std * (2.0 / _np.sqrt(H.shape[0])), "std", mxBasis)
    assert(_np.linalg.norm(_spl.expm(logM) - M) < 1e-8)  # expensive b/c of expm - could comment for performance
    return logM 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import logm [as 别名]
def near_identity_matrix_log(M, TOL=1e-8):
    """
    Construct the logarithm of superoperator matrix `M` that is
    near the identity.  If `M` is real, the resulting logarithm will be real.

    Parameters
    ----------
    M : numpy array
        The superoperator matrix whose logarithm is taken

    TOL : float, optional
        The tolerance used when testing for zero imaginary parts.

    Returns
    -------
    numpy array
        An matrix `logM`, of the same shape as `M`, such that `M = exp(logM)`
        and `logM` is real when `M` is real.
    """
    # A near-identity matrix should have a unique logarithm, and it should be
    # real if the original matrix is real
    M_is_real = bool(_np.linalg.norm(M.imag) < TOL)
    logM = _spl.logm(M)
    if M_is_real:
        assert(_np.linalg.norm(logM.imag) < TOL), \
            "Failed to construct a real logarithm! " \
            + "This is probably because M is not near the identity.\n" \
            + "Its eigenvalues are: " + str(_np.linalg.eigvals(M))
        logM = logM.real
    return logM 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import logm [as 别名]
def operation_from_error_generator(error_gen, target_op, typ="logG-logT"):
    """
    Construct a gate from an error generator and a target gate.

    Inverts the computation fone in :func:`error_generator` and
    returns the value of the gate given by
    gate = target_op * exp(error_gen).

    Parameters
    ----------
    error_gen : ndarray
      The error generator matrix

    target_op : ndarray
      The target operation matrix

    typ : {"logG-logT", "logTiG"}
      The type of error generator to compute.  Allowed values are:

      - "logG-logT" : errgen = log(gate) - log(target_op)
      - "logTiG" : errgen = log( dot(inv(target_op), gate) )


    Returns
    -------
    ndarray
      The operation matrix.
    """
    if typ == "logG-logT":
        return _spl.expm(error_gen + _spl.logm(target_op))
    elif typ == "logTiG":
        return _np.dot(target_op, _spl.expm(error_gen))
    elif typ == "logGTi":
        return _np.dot(_spl.expm(error_gen), target_op)
    else:
        raise ValueError("Invalid error-generator type: %s" % typ) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import logm [as 别名]
def compute_inplane_from_rotation(R, vertex):
    rot = compute_rotation_from_vertex(vertex)
    angle_axis = logm(np.matmul(R, rot.T))
    return angle_axis[1, 0] 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import logm [as 别名]
def rot_mag(R):
    angle = (1.0 / math.sqrt(2)) * \
        norm(logm(R), 'fro') * 180 / (math.pi)
    return angle 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import logm [as 别名]
def AccViewCls(output, target, numBins, specificView): 
  #unified
  binSize = 360. / numBins
  if specificView:
    acc = 0
    for t in range(target.shape[0]):
      idx = np.where(target[t] != numBins)
      ps = idx[0][0] / 3 * 3
      _, pred = output[t].view(-1, numBins)[ps: ps + 3].topk(1, 1, True, True)
      pred = pred.view(3).float() * binSize / 180. * pi 
      gt = target[t][ps: ps + 3].float() * binSize / 180. * pi
      R_pred = angle2dcm(pred)
      R_gt = angle2dcm(gt)
      err = ((logm(np.dot(np.transpose(R_pred), R_gt)) ** 2).sum()) ** 0.5 / sqrt2
      acc += 1 if err < pi / 6. else 0
    return 1.0 * acc / target.shape[0]
  else:
    _, pred = output.view(target.shape[0] * 3, numBins).topk(1, 1, True, True)
    pred = pred.view(target.shape[0], 3).float() * binSize / 180. * pi
    target = target.float() * binSize / 180. * pi
    acc = 0
    for t in range(target.shape[0]):
      R_pred = angle2dcm(pred[t])
      R_gt = angle2dcm(target[t])
      err = ((logm(np.dot(np.transpose(R_pred), R_gt)) ** 2).sum()) ** 0.5 / sqrt2
      acc += 1 if err < pi / 6. else 0
    return 1.0 * acc / target.shape[0] 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import logm [as 别名]
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) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import logm [as 别名]
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)
        #XXX: what would be the correct result? 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import logm [as 别名]
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)