Python numpy.trace() Examples
The following are 30
code examples of numpy.trace().
You can vote up the ones you like or vote down the ones you don't like,
and go to the original project or source file by following the links above each example.
You may also want to check out all available functions/classes of the module
numpy
, or try the search function
.
Example #1
Source File: transformations.py From rai-python with MIT License | 6 votes |
def reflection_matrix(point, normal): """Return matrix to mirror at plane defined by point and normal vector. >>> v0 = numpy.random.random(4) - 0.5 >>> v0[3] = 1.0 >>> v1 = numpy.random.random(3) - 0.5 >>> R = reflection_matrix(v0, v1) >>> numpy.allclose(2., numpy.trace(R)) True >>> numpy.allclose(v0, numpy.dot(R, v0)) True >>> v2 = v0.copy() >>> v2[:3] += v1 >>> v3 = v0.copy() >>> v2[:3] -= v1 >>> numpy.allclose(v2, numpy.dot(R, v3)) True """ normal = unit_vector(normal[:3]) M = numpy.identity(4) M[:3, :3] -= 2.0 * numpy.outer(normal, normal) M[:3, 3] = (2.0 * numpy.dot(point[:3], normal)) * normal return M
Example #2
Source File: fciqmc.py From pyscf with Apache License 2.0 | 6 votes |
def read_neci_1dms(fciqmcci, norb, nelec, filename='OneRDM.1', directory='.'): ''' Read spinned rdms, as they are in the neci output ''' f = open(os.path.join(directory, filename),'r') dm1a = numpy.zeros((norb,norb)) dm1b = numpy.zeros((norb,norb)) for line in f.readlines(): linesp = line.split() i, j = int(linesp[0]), int(linesp[1]) assert((i % 2) == (j % 2)) if i % 2 == 1: # alpha assert(all(x<norb for x in (i/2,j/2))) dm1a[i/2,j/2] = float(linesp[2]) dm1a[j/2,i/2] = float(linesp[2]) else: assert(all(x<norb for x in (i/2 - 1,j/2 - 1))) dm1b[i/2 - 1,j/2 - 1] = float(linesp[2]) dm1b[j/2 - 1,i/2 - 1] = float(linesp[2]) f.close() assert(numpy.allclose(dm1a.trace()+dm1b.trace(),sum(nelec))) return dm1a, dm1b
Example #3
Source File: test_nlinalg.py From D-VAE with MIT License | 6 votes |
def test_trace(): rng = numpy.random.RandomState(utt.fetch_seed()) x = theano.tensor.matrix() g = trace(x) f = theano.function([x], g) for shp in [(2, 3), (3, 2), (3, 3)]: m = rng.rand(*shp).astype(config.floatX) v = numpy.trace(m) assert v == f(m) xx = theano.tensor.vector() ok = False try: trace(xx) except TypeError: ok = True assert ok
Example #4
Source File: test_pose.py From SfmLearner-Pytorch with MIT License | 6 votes |
def compute_pose_error(gt, pred): RE = 0 snippet_length = gt.shape[0] scale_factor = np.sum(gt[:,:,-1] * pred[:,:,-1])/np.sum(pred[:,:,-1] ** 2) ATE = np.linalg.norm((gt[:,:,-1] - scale_factor * pred[:,:,-1]).reshape(-1)) for gt_pose, pred_pose in zip(gt, pred): # Residual matrix to which we compute angle's sin and cos R = gt_pose[:,:3] @ np.linalg.inv(pred_pose[:,:3]) s = np.linalg.norm([R[0,1]-R[1,0], R[1,2]-R[2,1], R[0,2]-R[2,0]]) c = np.trace(R) - 1 # Note: we actually compute double of cos and sin, but arctan2 is invariant to scale RE += np.arctan2(s,c) return ATE/snippet_length, RE/snippet_length
Example #5
Source File: transform.py From bop_toolkit with MIT License | 6 votes |
def reflection_matrix(point, normal): """Return matrix to mirror at plane defined by point and normal vector. >>> v0 = numpy.random.random(4) - 0.5 >>> v0[3] = 1. >>> v1 = numpy.random.random(3) - 0.5 >>> R = reflection_matrix(v0, v1) >>> numpy.allclose(2, numpy.trace(R)) True >>> numpy.allclose(v0, numpy.dot(R, v0)) True >>> v2 = v0.copy() >>> v2[:3] += v1 >>> v3 = v0.copy() >>> v2[:3] -= v1 >>> numpy.allclose(v2, numpy.dot(R, v3)) True """ normal = unit_vector(normal[:3]) M = numpy.identity(4) M[:3, :3] -= 2.0 * numpy.outer(normal, normal) M[:3, 3] = (2.0 * numpy.dot(point[:3], normal)) * normal return M
Example #6
Source File: basis.py From pyGSTi with Apache License 2.0 | 6 votes |
def is_normalized(self): ''' Check if a basis is normalized, meaning that Tr(Bi Bi) = 1.0. Available only to bases whose elements are *matrices* for now. Returns ------- bool ''' if self.elndim == 2: for i, mx in enumerate(self.elements): t = _np.trace(_np.dot(mx, mx)) t = _np.real(t) if not _np.isclose(t, 1.0): return False return True elif self.elndim == 1: raise NotImplementedError("TODO: add code so this works for *vector*-valued bases too!") else: raise ValueError("I don't know what normalized means for elndim == %d!" % self.elndim)
Example #7
Source File: transformations.py From cozmo_driver with Apache License 2.0 | 6 votes |
def reflection_matrix(point, normal): """Return matrix to mirror at plane defined by point and normal vector. >>> v0 = numpy.random.random(4) - 0.5 >>> v0[3] = 1.0 >>> v1 = numpy.random.random(3) - 0.5 >>> R = reflection_matrix(v0, v1) >>> numpy.allclose(2., numpy.trace(R)) True >>> numpy.allclose(v0, numpy.dot(R, v0)) True >>> v2 = v0.copy() >>> v2[:3] += v1 >>> v3 = v0.copy() >>> v2[:3] -= v1 >>> numpy.allclose(v2, numpy.dot(R, v3)) True """ normal = unit_vector(normal[:3]) M = numpy.identity(4) M[:3, :3] -= 2.0 * numpy.outer(normal, normal) M[:3, 3] = (2.0 * numpy.dot(point[:3], normal)) * normal return M
Example #8
Source File: optools.py From pyGSTi with Apache License 2.0 | 6 votes |
def tracenorm(A): """ Compute the trace norm of matrix A given by: Tr( sqrt{ A^dagger * A } ) Parameters ---------- A : numpy array The matrix to compute the trace norm of. """ if _np.linalg.norm(A - _np.conjugate(A.T)) < 1e-8: #Hermitian, so just sum eigenvalue magnitudes return _np.sum(_np.abs(_np.linalg.eigvals(A))) else: #Sum of singular values (positive by construction) return _np.sum(_np.linalg.svd(A, compute_uv=False))
Example #9
Source File: optools.py From pyGSTi with Apache License 2.0 | 6 votes |
def jtracedist(A, B, mxBasis='pp'): # Jamiolkowski trace distance: Tr(|J(A)-J(B)|) """ Compute the Jamiolkowski trace distance between operation matrices A and B, given by: D = 0.5 * Tr( sqrt{ (J(A)-J(B))^2 } ) where J(.) is the Jamiolkowski isomorphism map that maps a operation matrix to it's corresponding Choi Matrix. Parameters ---------- A, B : numpy array The matrices to compute the distance between. 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). """ JA = _jam.fast_jamiolkowski_iso_std(A, mxBasis) JB = _jam.fast_jamiolkowski_iso_std(B, mxBasis) return tracedist(JA, JB)
Example #10
Source File: optools.py From pyGSTi with Apache License 2.0 | 6 votes |
def povm_jtracedist(model, targetModel, povmlbl): """ Computes the Jamiolkowski trace distance between POVM maps using :func:`jtracedist`. Parameters ---------- model, targetModel : Model Models containing the two POVMs to compare. povmlbl : str The POVM label Returns ------- float """ povm_mx = get_povm_map(model, povmlbl) target_povm_mx = get_povm_map(targetModel, povmlbl) return jtracedist(povm_mx, target_povm_mx, targetModel.basis)
Example #11
Source File: geometry.py From connecting_the_dots with MIT License | 6 votes |
def axisangle_from_rotm(R): # logarithm of rotation matrix # R = R.reshape(-1,3,3) # tr = np.trace(R, axis1=1, axis2=2) # phi = np.arccos(np.clip((tr - 1) / 2, -1, 1)) # scale = np.zeros_like(phi) # div = 2 * np.sin(phi) # np.divide(phi, div, out=scale, where=np.abs(div) > 1e-6) # A = (R - R.transpose(0,2,1)) * scale.reshape(-1,1,1) # aa = np.stack((A[:,2,1], A[:,0,2], A[:,1,0]), axis=1) # return aa.squeeze() R = R.reshape(-1,3,3) omega = np.empty((R.shape[0], 3), dtype=R.dtype) omega[:,0] = R[:,2,1] - R[:,1,2] omega[:,1] = R[:,0,2] - R[:,2,0] omega[:,2] = R[:,1,0] - R[:,0,1] r = np.linalg.norm(omega, axis=1).reshape(-1,1) t = np.trace(R, axis1=1, axis2=2).reshape(-1,1) omega = np.arctan2(r, t-1) * omega aa = np.zeros_like(omega) np.divide(omega, r, out=aa, where=r != 0) return aa.squeeze()
Example #12
Source File: optimal_givens_decomposition_test.py From OpenFermion-Cirq with Apache License 2.0 | 6 votes |
def test_circuit_generation_and_accuracy(): for dim in range(2, 10): qubits = cirq.LineQubit.range(dim) u_generator = numpy.random.random( (dim, dim)) + 1j * numpy.random.random((dim, dim)) u_generator = u_generator - numpy.conj(u_generator).T assert numpy.allclose(-1 * u_generator, numpy.conj(u_generator).T) unitary = scipy.linalg.expm(u_generator) circuit = cirq.Circuit() circuit.append(optimal_givens_decomposition(qubits, unitary)) fermion_generator = QubitOperator(()) * 0.0 for i, j in product(range(dim), repeat=2): fermion_generator += jordan_wigner( FermionOperator(((i, 1), (j, 0)), u_generator[i, j])) true_unitary = scipy.linalg.expm( get_sparse_operator(fermion_generator).toarray()) assert numpy.allclose(true_unitary.conj().T.dot(true_unitary), numpy.eye(2 ** dim, dtype=complex)) test_unitary = cirq.unitary(circuit) assert numpy.isclose( abs(numpy.trace(true_unitary.conj().T.dot(test_unitary))), 2 ** dim)
Example #13
Source File: utils.py From lsm with MIT License | 5 votes |
def rot2quat(M): if M.shape[0] < 4 or M.shape[1] < 4: newM = np.zeros((4, 4)) newM[:3, :3] = M[:3, :3] newM[3, 3] = 1 M = newM q = np.empty((4, )) t = np.trace(M) if t > M[3, 3]: q[0] = t q[3] = M[1, 0] - M[0, 1] q[2] = M[0, 2] - M[2, 0] q[1] = M[2, 1] - M[1, 2] else: i, j, k = 0, 1, 2 if M[1, 1] > M[0, 0]: i, j, k = 1, 2, 0 if M[2, 2] > M[i, i]: i, j, k = 2, 0, 1 t = M[i, i] - (M[j, j] + M[k, k]) + M[3, 3] q[i] = t q[j] = M[i, j] + M[j, i] q[k] = M[k, i] + M[i, k] q[3] = M[k, j] - M[j, k] q = q[[3, 0, 1, 2]] q *= 0.5 / math.sqrt(t * M[3, 3]) return q
Example #14
Source File: transform.py From bop_toolkit with MIT License | 5 votes |
def rotation_from_matrix(matrix): """Return rotation angle and axis from rotation matrix. >>> angle = (random.random() - 0.5) * (2*math.pi) >>> direc = numpy.random.random(3) - 0.5 >>> point = numpy.random.random(3) - 0.5 >>> R0 = rotation_matrix(angle, direc, point) >>> angle, direc, point = rotation_from_matrix(R0) >>> R1 = rotation_matrix(angle, direc, point) >>> is_same_transform(R0, R1) True """ R = numpy.array(matrix, dtype=numpy.float64, copy=False) R33 = R[:3, :3] # direction: unit eigenvector of R33 corresponding to eigenvalue of 1 w, W = numpy.linalg.eig(R33.T) i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0] if not len(i): raise ValueError("no unit eigenvector corresponding to eigenvalue 1") direction = numpy.real(W[:, i[-1]]).squeeze() # point: unit eigenvector of R33 corresponding to eigenvalue of 1 w, Q = numpy.linalg.eig(R) i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0] if not len(i): raise ValueError("no unit eigenvector corresponding to eigenvalue 1") point = numpy.real(Q[:, i[-1]]).squeeze() point /= point[3] # rotation angle depending on direction cosa = (numpy.trace(R33) - 1.0) / 2.0 if abs(direction[2]) > 1e-8: sina = (R[1, 0] + (cosa - 1.0) * direction[0] * direction[1]) / direction[2] elif abs(direction[1]) > 1e-8: sina = (R[0, 2] + (cosa - 1.0) * direction[0] * direction[2]) / direction[1] else: sina = (R[2, 1] + (cosa - 1.0) * direction[1] * direction[2]) / direction[0] angle = math.atan2(sina, cosa) return angle, direction, point
Example #15
Source File: util.py From RelativePose with BSD 3-Clause "New" or "Revised" License | 5 votes |
def angular_distance_np(R_hat, R): # measure the angular distance between two rotation matrice # R1,R2: [n, 3, 3] if R_hat.shape == (3,3): R_hat = R_hat[np.newaxis,:] if R.shape == (3,3): R = R[np.newaxis,:] n = R.shape[0] trace_idx = [0,4,8] trace = np.matmul(R_hat, R.transpose(0,2,1)).reshape(n,-1)[:,trace_idx].sum(1) metric = np.arccos(((trace - 1)/2).clip(-1,1)) / np.pi * 180.0 return metric
Example #16
Source File: transform.py From bop_toolkit with MIT License | 5 votes |
def identity_matrix(): """Return 4x4 identity/unit matrix. >>> I = identity_matrix() >>> numpy.allclose(I, numpy.dot(I, I)) True >>> numpy.sum(I), numpy.trace(I) (4.0, 4.0) >>> numpy.allclose(I, numpy.identity(4)) True """ return numpy.identity(4)
Example #17
Source File: transform.py From bop_toolkit with MIT License | 5 votes |
def scale_from_matrix(matrix): """Return scaling factor, origin and direction from scaling matrix. >>> factor = random.random() * 10 - 5 >>> origin = numpy.random.random(3) - 0.5 >>> direct = numpy.random.random(3) - 0.5 >>> S0 = scale_matrix(factor, origin) >>> factor, origin, direction = scale_from_matrix(S0) >>> S1 = scale_matrix(factor, origin, direction) >>> is_same_transform(S0, S1) True >>> S0 = scale_matrix(factor, origin, direct) >>> factor, origin, direction = scale_from_matrix(S0) >>> S1 = scale_matrix(factor, origin, direction) >>> is_same_transform(S0, S1) True """ M = numpy.array(matrix, dtype=numpy.float64, copy=False) M33 = M[:3, :3] factor = numpy.trace(M33) - 2.0 try: # direction: unit eigenvector corresponding to eigenvalue factor w, V = numpy.linalg.eig(M33) i = numpy.where(abs(numpy.real(w) - factor) < 1e-8)[0][0] direction = numpy.real(V[:, i]).squeeze() direction /= vector_norm(direction) except IndexError: # uniform scaling factor = (factor + 2.0) / 3.0 direction = None # origin: any eigenvector corresponding to eigenvalue 1 w, V = numpy.linalg.eig(M) i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0] if not len(i): raise ValueError("no eigenvector corresponding to eigenvalue 1") origin = numpy.real(V[:, i[-1]]).squeeze() origin /= origin[3] return factor, origin, direction
Example #18
Source File: test_core.py From auto-alt-text-lambda-api with MIT License | 5 votes |
def test_trace(self): # Tests trace on MaskedArrays. (x, X, XX, m, mx, mX, mXX, m2x, m2X, m2XX) = self.d mXdiag = mX.diagonal() assert_equal(mX.trace(), mX.diagonal().compressed().sum()) assert_almost_equal(mX.trace(), X.trace() - sum(mXdiag.mask * X.diagonal(), axis=0)) assert_equal(np.trace(mX), mX.trace())
Example #19
Source File: transformations.py From cozmo_driver with Apache License 2.0 | 5 votes |
def identity_matrix(): """Return 4x4 identity/unit matrix. >>> I = identity_matrix() >>> numpy.allclose(I, numpy.dot(I, I)) True >>> numpy.sum(I), numpy.trace(I) (4.0, 4.0) >>> numpy.allclose(I, numpy.identity(4, dtype=numpy.float64)) True """ return numpy.identity(4, dtype=numpy.float64)
Example #20
Source File: _expm_multiply.py From lambda-packs with MIT License | 5 votes |
def _trace(A): # A compatibility function which should eventually disappear. if scipy.sparse.isspmatrix(A): return A.diagonal().sum() else: return np.trace(A)
Example #21
Source File: test_numeric.py From auto-alt-text-lambda-api with MIT License | 5 votes |
def test_trace(self): c = [[1, 2], [3, 4], [5, 6]] assert_equal(np.trace(c), 5)
Example #22
Source File: test_numeric.py From recruit with Apache License 2.0 | 5 votes |
def test_trace(self): c = [[1, 2], [3, 4], [5, 6]] assert_equal(np.trace(c), 5)
Example #23
Source File: test_core.py From recruit with Apache License 2.0 | 5 votes |
def test_trace(self): # Tests trace on MaskedArrays. (x, X, XX, m, mx, mX, mXX, m2x, m2X, m2XX) = self.d mXdiag = mX.diagonal() assert_equal(mX.trace(), mX.diagonal().compressed().sum()) assert_almost_equal(mX.trace(), X.trace() - sum(mXdiag.mask * X.diagonal(), axis=0)) assert_equal(np.trace(mX), mX.trace()) # gh-5560 arr = np.arange(2*4*4).reshape(2,4,4) m_arr = np.ma.masked_array(arr, False) assert_equal(arr.trace(axis1=1, axis2=2), m_arr.trace(axis1=1, axis2=2))
Example #24
Source File: fid.py From chainer-stylegan with MIT License | 5 votes |
def calc_FID(m0, c0, m1, c1): ret = 0 ret += np.sum((m0-m1)**2) ret += np.trace(c0 + c1 - 2.0 * scipy.linalg.sqrtm(np.dot(c0, c1))) return np.real(ret)
Example #25
Source File: test_core.py From vnpy_crypto with MIT License | 5 votes |
def test_trace(self): # Tests trace on MaskedArrays. (x, X, XX, m, mx, mX, mXX, m2x, m2X, m2XX) = self.d mXdiag = mX.diagonal() assert_equal(mX.trace(), mX.diagonal().compressed().sum()) assert_almost_equal(mX.trace(), X.trace() - sum(mXdiag.mask * X.diagonal(), axis=0)) assert_equal(np.trace(mX), mX.trace()) # gh-5560 arr = np.arange(2*4*4).reshape(2,4,4) m_arr = np.ma.masked_array(arr, False) assert_equal(arr.trace(axis1=1, axis2=2), m_arr.trace(axis1=1, axis2=2))
Example #26
Source File: test_numeric.py From vnpy_crypto with MIT License | 5 votes |
def test_trace(self): c = [[1, 2], [3, 4], [5, 6]] assert_equal(np.trace(c), 5)
Example #27
Source File: transformations.py From rai-python with MIT License | 5 votes |
def rotation_from_matrix(matrix): """Return rotation angle and axis from rotation matrix. >>> angle = (random.random() - 0.5) * (2*math.pi) >>> direc = numpy.random.random(3) - 0.5 >>> point = numpy.random.random(3) - 0.5 >>> R0 = rotation_matrix(angle, direc, point) >>> angle, direc, point = rotation_from_matrix(R0) >>> R1 = rotation_matrix(angle, direc, point) >>> is_same_transform(R0, R1) True """ R = numpy.array(matrix, dtype=numpy.float64, copy=False) R33 = R[:3, :3] # direction: unit eigenvector of R33 corresponding to eigenvalue of 1 l, W = numpy.linalg.eig(R33.T) i = numpy.where(abs(numpy.real(l) - 1.0) < 1e-8)[0] if not len(i): raise ValueError("no unit eigenvector corresponding to eigenvalue 1") direction = numpy.real(W[:, i[-1]]).squeeze() # point: unit eigenvector of R33 corresponding to eigenvalue of 1 l, Q = numpy.linalg.eig(R) i = numpy.where(abs(numpy.real(l) - 1.0) < 1e-8)[0] if not len(i): raise ValueError("no unit eigenvector corresponding to eigenvalue 1") point = numpy.real(Q[:, i[-1]]).squeeze() point /= point[3] # rotation angle depending on direction cosa = (numpy.trace(R33) - 1.0) / 2.0 if abs(direction[2]) > 1e-8: sina = (R[1, 0] + (cosa-1.0)*direction[0]*direction[1]) / direction[2] elif abs(direction[1]) > 1e-8: sina = (R[0, 2] + (cosa-1.0)*direction[0]*direction[2]) / direction[1] else: sina = (R[2, 1] + (cosa-1.0)*direction[1]*direction[2]) / direction[0] angle = math.atan2(sina, cosa) return angle, direction, point
Example #28
Source File: transformations.py From rai-python with MIT License | 5 votes |
def rotation_matrix(angle, direction, point=None): """Return matrix to rotate about axis defined by point and direction. >>> angle = (random.random() - 0.5) * (2*math.pi) >>> direc = numpy.random.random(3) - 0.5 >>> point = numpy.random.random(3) - 0.5 >>> R0 = rotation_matrix(angle, direc, point) >>> R1 = rotation_matrix(angle-2*math.pi, direc, point) >>> is_same_transform(R0, R1) True >>> R0 = rotation_matrix(angle, direc, point) >>> R1 = rotation_matrix(-angle, -direc, point) >>> is_same_transform(R0, R1) True >>> I = numpy.identity(4, numpy.float64) >>> numpy.allclose(I, rotation_matrix(math.pi*2, direc)) True >>> numpy.allclose(2., numpy.trace(rotation_matrix(math.pi/2, ... direc, point))) True """ sina = math.sin(angle) cosa = math.cos(angle) direction = unit_vector(direction[:3]) # rotation matrix around unit vector R = numpy.array(((cosa, 0.0, 0.0), (0.0, cosa, 0.0), (0.0, 0.0, cosa)), dtype=numpy.float64) R += numpy.outer(direction, direction) * (1.0 - cosa) direction *= sina R += numpy.array((( 0.0, -direction[2], direction[1]), ( direction[2], 0.0, -direction[0]), (-direction[1], direction[0], 0.0)), dtype=numpy.float64) M = numpy.identity(4) M[:3, :3] = R if point is not None: # rotation not around origin point = numpy.array(point[:3], dtype=numpy.float64, copy=False) M[:3, 3] = point - numpy.dot(R, point) return M
Example #29
Source File: svar_model.py From vnpy_crypto with MIT License | 5 votes |
def loglike(self, params): """ Loglikelihood for SVAR model Notes ----- This method assumes that the autoregressive parameters are first estimated, then likelihood with structural parameters is estimated """ #TODO: this doesn't look robust if A or B is None A = self.A B = self.B A_mask = self.A_mask B_mask = self.B_mask A_len = len(A[A_mask]) B_len = len(B[B_mask]) if A is not None: A[A_mask] = params[:A_len] if B is not None: B[B_mask] = params[A_len:A_len+B_len] nobs = self.nobs neqs = self.neqs sigma_u = self.sigma_u W = np.dot(npl.inv(B),A) trc_in = np.dot(np.dot(W.T,W),sigma_u) sign, b_logdet = slogdet(B**2) #numpy 1.4 compat b_slogdet = sign * b_logdet likl = -nobs/2. * (neqs * np.log(2 * np.pi) - \ np.log(npl.det(A)**2) + b_slogdet + \ np.trace(trc_in)) return likl
Example #30
Source File: transformations.py From cozmo_driver with Apache License 2.0 | 5 votes |
def rotation_from_matrix(matrix): """Return rotation angle and axis from rotation matrix. >>> angle = (random.random() - 0.5) * (2*math.pi) >>> direc = numpy.random.random(3) - 0.5 >>> point = numpy.random.random(3) - 0.5 >>> R0 = rotation_matrix(angle, direc, point) >>> angle, direc, point = rotation_from_matrix(R0) >>> R1 = rotation_matrix(angle, direc, point) >>> is_same_transform(R0, R1) True """ R = numpy.array(matrix, dtype=numpy.float64, copy=False) R33 = R[:3, :3] # direction: unit eigenvector of R33 corresponding to eigenvalue of 1 l, W = numpy.linalg.eig(R33.T) i = numpy.where(abs(numpy.real(l) - 1.0) < 1e-8)[0] if not len(i): raise ValueError("no unit eigenvector corresponding to eigenvalue 1") direction = numpy.real(W[:, i[-1]]).squeeze() # point: unit eigenvector of R33 corresponding to eigenvalue of 1 l, Q = numpy.linalg.eig(R) i = numpy.where(abs(numpy.real(l) - 1.0) < 1e-8)[0] if not len(i): raise ValueError("no unit eigenvector corresponding to eigenvalue 1") point = numpy.real(Q[:, i[-1]]).squeeze() point /= point[3] # rotation angle depending on direction cosa = (numpy.trace(R33) - 1.0) / 2.0 if abs(direction[2]) > 1e-8: sina = (R[1, 0] + (cosa-1.0)*direction[0]*direction[1]) / direction[2] elif abs(direction[1]) > 1e-8: sina = (R[0, 2] + (cosa-1.0)*direction[0]*direction[2]) / direction[1] else: sina = (R[2, 1] + (cosa-1.0)*direction[1]*direction[2]) / direction[0] angle = math.atan2(sina, cosa) return angle, direction, point