Python numpy.matmul() Examples

def _BtDB(self,s,r):
        """
        dot product of B^T, D, B
        params:
            s,r:natural position of evalue point.2-array.
        returns:
            3x3 matrix.
        """
        print(self._B(s,r).transpose(2,0,1).shape)
        print(
            np.matmul(
                np.dot(self._B(s,r).T,self._D),
                self._B(s,r).transpose(2,0,1)).shape
            )
        print(self._D.shape)
       
        
        return np.matmul(np.dot(self._B(s,r).T,self._D),self._B(s,r).transpose(2,0,1)).transpose(1,2,0) 

def DataProjection(X, r, type='NormalProj'):
    Xp = None
    D, N = X.shape
    if r == 0:
        Xp = X
    else:
        if type == 'PCA':
            isEcon = False
            if D > N:
                isEcon = True
            U, S, V = np.linalg.svd(X.T, full_matrices=isEcon)
            Xp = U[:, 0:r].T
        if type == 'NormalProj':
            normP = (1.0 / math.sqrt(r)) * np.random.randn(r * D, 1)
            PrN = normP.reshape(r, D, order='F')
            Xp = np.matmul(PrN, X)
        if type == 'BernoulliProj':
            bp = np.random.rand(r * D, 1)
            Bp = (1.0 / math.sqrt(r)) * (bp >= .5) - (1.0 / math.sqrt(r)) * (bp < .5)
            PrB = Bp.reshape(r, D, order='F')
            Xp = np.matmul(PrB, X)
    return Xp 

def SpectralClustering(CKSym, n):
    # This is direct port of JHU vision lab code. Could probably use sklearn SpectralClustering.
    CKSym = CKSym.astype(float)
    N, _ = CKSym.shape
    MAXiter = 1000  # Maximum number of iterations for KMeans
    REPlic = 20  # Number of replications for KMeans

    DN = np.diag(np.divide(1, np.sqrt(np.sum(CKSym, axis=0) + np.finfo(float).eps)))
    LapN = identity(N).toarray().astype(float) - np.matmul(np.matmul(DN, CKSym), DN)
    _, _, vN = np.linalg.svd(LapN)
    vN = vN.T
    kerN = vN[:, N - n:N]
    normN = np.sqrt(np.sum(np.square(kerN), axis=1))
    kerNS = np.divide(kerN, normN.reshape(len(normN), 1) + np.finfo(float).eps)
    km = KMeans(n_clusters=n, n_init=REPlic, max_iter=MAXiter, n_jobs=-1).fit(kerNS)
    return km.labels_ 

def _evaluate(self, X, out, *args, **kwargs):
        if self.w == 0:
            X1 = X[:, 0]
            X2 = X[:, 1]
        else:
            # If rotated, we rotate it back by applying the inverted rotation matrix to X
            Y = np.array([np.matmul(self.IRM, x) for x in X])
            X1 = Y[:, 0]
            X2 = Y[:, 1]

        a, b, c = self.a, self.b, self.c
        t1_hat = sign(X1) * ceil((abs(X1) - a - c / 2) / (2 * a + c))
        t2_hat = sign(X2) * ceil((abs(X2) - b / 2) / b)
        one = ones(len(X))
        t1 = sign(t1_hat) * min(np.vstack((abs(t1_hat), one)), axis=0)
        t2 = sign(t2_hat) * min(np.vstack((abs(t2_hat), one)), axis=0)

        p1 = X1 - t1 * c
        p2 = X2 - t2 * b

        f1 = (p1 + a) ** 2 + p2 ** 2
        f2 = (p1 - a) ** 2 + p2 ** 2
        out["F"] = np.vstack((f1, f2)).T 

def _calc_pareto_set(self, n_pareto_points=500):
        # The SYM-PART test problem has 9 equivalent Pareto subsets.
        h = int(n_pareto_points / 9)
        PS = zeros((h * 9, self.n_var))
        cnt = 0
        for row in [-1, 0, 1]:
            for col in [1, 0, -1]:
                X1 = np.linspace(row * self.c - self.a, row * self.c + self.a, h)
                X2 = np.tile(col * self.b, h)
                PS[cnt * h:cnt * h + h, :] = np.vstack((X1, X2)).T
                cnt = cnt + 1
        if self.w != 0:
            # If rotated, we apply the rotation matrix to PS
            # Calculate the rotation matrix
            RM = np.array([
                [cos(self.w), -sin(self.w)],
                [sin(self.w), cos(self.w)]
            ])
            PS = np.array([np.matmul(RM, x) for x in PS])
        return PS 

def overlap_ni(me, sp1,R1, sp2,R2, **kw):
    """
      Computes overlap for an atom pair. The atom pair is given by a pair of species indices
      and the coordinates of the atoms.
      Args: 
        sp1,sp2 : specie indices, and
        R1,R2 :   respective coordinates in Bohr, atomic units
      Result:
        matrix of orbital overlaps
      The procedure uses the numerical integration in coordinate space.
    """
    from pyscf.nao.m_ao_matelem import build_3dgrid
    from pyscf.nao.m_ao_eval_libnao import ao_eval_libnao as ao_eval #_libnao

    grids = build_3dgrid(me, sp1, R1, sp2, R2, **kw)

    ao1 = ao_eval(me.ao1, R1, sp1, grids.coords)
    ao1 = ao1 * grids.weights

    ao2 = ao_eval(me.ao2, R2, sp2, grids.coords)

    overlaps = np.einsum("ij,kj->ik", ao1, ao2) #      overlaps = np.matmul(ao1, ao2.T)

    return overlaps 

def get_3d_box_batch(box_size, heading_angle, center):
    ''' box_size: [x1,x2,...,xn,3]
        heading_angle: [x1,x2,...,xn]
        center: [x1,x2,...,xn,3]
    Return:
        [x1,x3,...,xn,8,3]
    '''
    input_shape = heading_angle.shape
    R = roty_batch(heading_angle)
    l = np.expand_dims(box_size[...,0], -1) # [x1,...,xn,1]
    w = np.expand_dims(box_size[...,1], -1)
    h = np.expand_dims(box_size[...,2], -1)
    corners_3d = np.zeros(tuple(list(input_shape)+[8,3]))
    corners_3d[...,:,0] = np.concatenate((l/2,l/2,-l/2,-l/2,l/2,l/2,-l/2,-l/2), -1)
    corners_3d[...,:,1] = np.concatenate((h/2,h/2,h/2,h/2,-h/2,-h/2,-h/2,-h/2), -1)
    corners_3d[...,:,2] = np.concatenate((w/2,-w/2,-w/2,w/2,w/2,-w/2,-w/2,w/2), -1)
    tlist = [i for i in range(len(input_shape))]
    tlist += [len(input_shape)+1, len(input_shape)]
    corners_3d = np.matmul(corners_3d, np.transpose(R, tuple(tlist)))
    corners_3d += np.expand_dims(center, -2)
    return corners_3d 

def rotate_molecule(path, target_path, count=10):
    # Load dataset
    mols = Chem.SDMolSupplier(path)
    rotated_mols = []

    print("Loaded {} Molecules from {}".format(len(mols), path))

    print("Rotating Molecules...")
    for mol in mols:
        for _ in range(count):
            for atom in mol.GetAtoms():
                atom_idx = atom.GetIdx()

                pos = list(mol.GetConformer().GetAtomPosition(atom_idx))
                pos_rotated = np.matmul(random_rotation_matrix(), pos)

                mol.GetConformer().SetAtomPosition(atom_idx, pos_rotated)

            rotated_mols.append(mol)

    w = Chem.SDWriter(target_path)
    for m in rotated_mols:
        if m is not None:
            w.write(m)
    print("Saved {} Molecules to {}".format(len(rotated_mols), target_path)) 

def backward(self, Z, O, E, W):
        outdim = W.shape[-1]
        time, batch, zdim = Z.shape
        indim = zdim - outdim

        bwO = self.actfn.backward(O)

        for t in range(time-1, -1, -1):
            E[t] *= bwO[t]
            deltaZ = np.dot(E[t], W.T)
            E[t-1] += deltaZ[:, indim:] if t else 0.

        dX = E[:, :, :indim]
        nablaW = np.matmul(Z.transpose(0, 2, 1), E).sum(axis=0)
        nablab = E.sum(axis=(0, 1))
        return dX, nablaW, nablab 

def valid(A, F):
        im, ic, iy, ix = A.shape
        nf, fc, fy, fx = F.shape
        # fx, fy, fc, nf = F.shape
        recfield_size = fx * fy * fc
        oy, ox = iy - fy + 1, ix - fx + 1
        rfields = np.zeros((im, oy*ox, recfield_size))
        Frsh = F.reshape(nf, recfield_size)

        if fc != ic:
            err = "Supplied filter (F) is incompatible with supplied input! (X)\n"
            err += "input depth: {} != {} :filter depth".format(ic, fc)
            raise ValueError(err)

        for i, sy, sx in ((idx, shy, shx) for shx in range(ox) for shy in range(oy) for idx in range(im)):
            rfields[i][sy*ox + sx] = A[i, :, sy:sy+fy, sx:sx+fx].ravel()

        # output = np.zeros((im, oy*ox, nf))
        # for m in range(im):
        #     output[m] = np.dot(rfields[m], Frsh.T)

        output = np.matmul(rfields, Frsh.T)
        output = output.transpose((0, 2, 1)).reshape((im, nf, oy, ox))
        return output 

def _parse_gufunc_signature(signature):
    """
    Parse string signatures for a generalized universal function.

    Arguments
    ---------
    signature : string
        Generalized universal function signature, e.g., ``(m,n),(n,p)->(m,p)``
        for ``np.matmul``.

    Returns
    -------
    Tuple of input and output core dimensions parsed from the signature, each
    of the form List[Tuple[str, ...]].
    """
    if not re.match(_SIGNATURE, signature):
        raise ValueError(
            'not a valid gufunc signature: {}'.format(signature))
    return tuple([tuple(re.findall(_DIMENSION_NAME, arg))
                  for arg in re.findall(_ARGUMENT, arg_list)]
                 for arg_list in signature.split('->')) 

def test_basic_property(self):
        # Check A = L L^H
        shapes = [(1, 1), (2, 2), (3, 3), (50, 50), (3, 10, 10)]
        dtypes = (np.float32, np.float64, np.complex64, np.complex128)

        for shape, dtype in itertools.product(shapes, dtypes):
            np.random.seed(1)
            a = np.random.randn(*shape)
            if np.issubdtype(dtype, np.complexfloating):
                a = a + 1j*np.random.randn(*shape)

            t = list(range(len(shape)))
            t[-2:] = -1, -2

            a = np.matmul(a.transpose(t).conj(), a)
            a = np.asarray(a, dtype=dtype)

            c = np.linalg.cholesky(a)

            b = np.matmul(c, c.transpose(t).conj())
            assert_allclose(b, a,
                            err_msg="{} {}\n{}\n{}".format(shape, dtype, a, c),
                            atol=500 * a.shape[0] * np.finfo(dtype).eps) 

def dictionary_update(self, Y: np.ndarray):
        # iterate rows
        D = self.dictionary.matrix
        n, K = D.shape
        for k in range(K):
            wk = np.nonzero(self.alphas[k, :])[0]
            if len(wk) == 0:
                continue
            D[:, k] = 0
            g = np.transpose(self.alphas)[wk, k]
            d = np.matmul(Y[:, wk], g) - np.matmul(D, self.alphas[:, wk]).dot(g)
            d = d / np.linalg.norm(d)
            g = np.matmul(Y[:, wk].T, d) - np.transpose(np.matmul(D, self.alphas[:, wk])).dot(d)
            D[:, k] = d
            self.alphas[k, wk] = g.T
        self.dictionary = Dictionary(D) 

def rgb2ycbcr(img, only_y=True):
    '''same as matlab rgb2ycbcr
    only_y: only return Y channel
    Input:
        uint8, [0, 255]
        float, [0, 1]
    '''
    in_img_type = img.dtype
    img.astype(np.float32)
    if in_img_type != np.uint8:
        img *= 255.
    # convert
    if only_y:
        rlt = np.dot(img, [65.481, 128.553, 24.966]) / 255.0 + 16.0
    else:
        rlt = np.matmul(img, [[65.481, -37.797, 112.0], [128.553, -74.203, -93.786],
                              [24.966, 112.0, -18.214]]) / 255.0 + [16, 128, 128]
    if in_img_type == np.uint8:
        rlt = rlt.round()
    else:
        rlt /= 255.
    return rlt.astype(in_img_type) 

def ycbcr2rgb(img):
    '''same as matlab ycbcr2rgb
    Input:
        uint8, [0, 255]
        float, [0, 1]
    '''
    in_img_type = img.dtype
    img.astype(np.float32)
    if in_img_type != np.uint8:
        img *= 255.
    # convert
    rlt = np.matmul(img, [[0.00456621, 0.00456621, 0.00456621], [0, -0.00153632, 0.00791071],
                          [0.00625893, -0.00318811, 0]]) * 255.0 + [-222.921, 135.576, -276.836]
    if in_img_type == np.uint8:
        rlt = rlt.round()
    else:
        rlt /= 255.
    return rlt.astype(in_img_type) 

def bgr2ycbcr(img, only_y=True):
    '''bgr version of rgb2ycbcr
    only_y: only return Y channel
    Input:
        uint8, [0, 255]
        float, [0, 1]
    '''
    in_img_type = img.dtype
    img.astype(np.float32)
    if in_img_type != np.uint8:
        img *= 255.
    # convert
    if only_y:
        rlt = np.dot(img, [24.966, 128.553, 65.481]) / 255.0 + 16.0
    else:
        rlt = np.matmul(img, [[24.966, 112.0, -18.214], [128.553, -74.203, -93.786],
                              [65.481, -37.797, 112.0]]) / 255.0 + [16, 128, 128]
    if in_img_type == np.uint8:
        rlt = rlt.round()
    else:
        rlt /= 255.
    return rlt.astype(in_img_type) 

def rotateBox(box_list,rotate_matrix,h,w):
    trans_box_list = []
    for bbx in box_list:
        bbx = [[bbx[0]-w//2,bbx[2]-w//2,bbx[4]-w//2,bbx[6]-w//2],
               [bbx[1]-h//2,bbx[3]-h//2,bbx[5]-h//2,bbx[7]-h//2]]
        trans_box_list.append(bbx)
    if len(trans_box_list) == 0:
        return []
    else:
        res_box_list = []
        for bbx in trans_box_list:
            bbx = np.matmul(rotate_matrix,np.array(bbx))
            bbx = bbx + np.array([
                [w//2,w//2,w//2,w//2],
                [h//2,h//2,h//2,h//2]
            ])
            x_mean = np.mean(bbx[0])
            y_mean = np.mean(bbx[1])
            if 0 < x_mean < w and 0 < y_mean < h:
                bbx = [bbx[0,0],bbx[1,0],bbx[0,1],bbx[1,1],bbx[0,2],bbx[1,2],bbx[0,3],bbx[1,3]]
                res_box_list.append(bbx)
        return res_box_list 

def _rmatvec(self, x):
        x = np.squeeze(x.reshape(self.nsl, self.nx, self.nz))
        if self.usematmul:
            if self.nz == 1:
                x = x[..., np.newaxis]
            if hasattr(self, 'GT'):
                y = np.matmul(self.GT, x)
            else:
                y = np.matmul(self.G.transpose((0, 2, 1)).conj(), x)
        else:
            y = np.squeeze(np.zeros((self.nsl, self.ny, self.nz),
                                    dtype=self.dtype))
            if hasattr(self, 'GT'):
                for isl in range(self.nsl):
                    y[isl] = np.dot(self.GT[isl], x[isl])
            else:
                for isl in range(self.nsl):
                    y[isl] = np.dot(self.G[isl].conj().T, x[isl])
        return y.ravel() 

def compute_conjunctions(self, X):
        """Compute conjunctions of features as specified in self.z.

        Args:
            X (DataFrame): Binarized features with MultiIndex column labels
        Returns:
            array: A -- Feature conjunction values, shape (X.shape[0], self.z.shape[1])
        """
        try:
            #A = 1 - (np.matmul(1 - X, self.z) > 0)
            A = 1 - (((1 - X) @ self.z) > 0)
        except AttributeError:
            print("Attribute 'z' does not exist, please fit model first.")
        # Scale conjunctions as done in training
        A = A * self.lambda0 / (self.lambda0 + self.lambda1 * self.z.sum().values)
        return A 

def make_geocentric(XYZ, sensor_height, camera_elevation_degree):
  """Transforms the point cloud into geocentric coordinate frame.
  Input:
    XYZ                     : ...x3
    sensor_height           : height of the sensor
    camera_elevation_degree : camera elevation to rectify.
  Output:
    XYZ : ...x3
  """
  R = ru.get_r_matrix([1.,0.,0.], angle=np.deg2rad(camera_elevation_degree))
  XYZ = np.matmul(XYZ.reshape(-1,3), R.T).reshape(XYZ.shape)
  XYZ[...,2] = XYZ[...,2] + sensor_height
  return XYZ 

def _step(self, a):
		self.robot.apply_action(a)
		self.scene.global_step()

		state = self.robot.calc_state()  # sets self.to_target_vec

		potential_old = self.potential
		self.potential = self.robot.calc_potential()

		joint_vel = np.array([
			self.robot.shoulder_pan_joint.get_velocity(),
			self.robot.shoulder_lift_joint.get_velocity(),
			self.robot.upper_arm_roll_joint.get_velocity(),
			self.robot.elbow_flex_joint.get_velocity(),
			self.robot.upper_arm_roll_joint.get_velocity(),
			self.robot.wrist_flex_joint.get_velocity(),
			self.robot.wrist_roll_joint.get_velocity()
		])

		action_product = np.matmul(np.abs(a), np.abs(joint_vel))
		action_sum = np.sum(a)

		electricity_cost = (
			-0.10 * action_product  # work torque*angular_velocity
			- 0.01 * action_sum  # stall torque require some energy
		)

		stuck_joint_cost = 0
		for j in self.robot.ordered_joints:
			if np.abs(j.current_relative_position()[0]) - 1 < 0.01:
				stuck_joint_cost += -0.1

		self.rewards = [float(self.potential - potential_old), float(electricity_cost), float(stuck_joint_cost)]
		self.HUD(state, a, False)
		return state, sum(self.rewards), False, {} 

def icp_test(S, M, tree_M, M_sampled, tree_M_sampled, S_given, gt_pose, MAX_LOOPS, ftol):
	from math import sin, cos
	import matplotlib.pyplot as plt
	from mpl_toolkits.mplot3d import Axes3D
 
	icp = ICP(M_sampled, S, tree_M_sampled)
	start = time.time()

	matrix, points, itr, neighbor_idx = icp.compute(MAX_LOOPS, ftol)		# Perform the iterations.

	end = time.time()

	final_pose = helper.find_final_pose((np.linalg.inv(matrix)).reshape((1,4,4)))

	title = "Actual T (Red->Green): "
	for i in range(len(gt_pose[0])):
		if i>2:
			title += str(round(gt_pose[0][i]*(180/np.pi),2))
		else:
			title += str(gt_pose[0][i])
		title += ', '
	title += "\nPredicted T (Red->Blue): "

	for i in range(len(final_pose[0])):
		if i>2:
			title += str(round(final_pose[0,i]*(180/np.pi),3))
		else:
			title += str(round(final_pose[0,i],3))
		title += ', '

	title += '\nElapsed Time: '+str(np.round((end-start)*1000,3))+' ms'+' & Iterations: '+str(itr)
	title += ' & Method: ICP'

	rot = matrix[0:3,0:3]
	tra = np.reshape(np.transpose(matrix[0:3,3]), (3,1))
	transformed = np.matmul(rot, np.transpose(S_given)) + tra
	return final_pose, M, S_given, transformed.T, title, end-start, itr 

def transform(point, center, scale, resolution, rotation=0, invert=False):
    _pt = np.ones(3)
    _pt[0] = point[0]
    _pt[1] = point[1]

    h = 200.0 * scale
    t = np.eye(3)
    t[0, 0] = resolution / h
    t[1, 1] = resolution / h
    t[0, 2] = resolution * (-center[0] / h + 0.5)
    t[1, 2] = resolution * (-center[1] / h + 0.5)

    if rotation != 0:
        rotation = -rotation
        r = np.eye(3)
        ang = rotation * math.pi / 180.0
        s = math.sin(ang)
        c = math.cos(ang)
        r[0][0] = c
        r[0][1] = -s
        r[1][0] = s
        r[1][1] = c

        t_ = np.eye(3)
        t_[0][2] = -resolution / 2.0
        t_[1][2] = -resolution / 2.0
        t_inv = torch.eye(3)
        t_inv[0][2] = resolution / 2.0
        t_inv[1][2] = resolution / 2.0
        t = reduce(np.matmul, [t_inv, r, t_, t])

    if invert:
        t = np.linalg.inv(t)
    new_point = (np.matmul(t, _pt))[0:2]

    return new_point.astype(float) 

def transform(point, center, scale, resolution, rotation=0, invert=False):
    _pt = np.ones(3)
    _pt[0] = point[0]
    _pt[1] = point[1]

    h = 200.0 * scale
    t = np.eye(3)
    t[0, 0] = resolution / h
    t[1, 1] = resolution / h
    t[0, 2] = resolution * (-center[0] / h + 0.5)
    t[1, 2] = resolution * (-center[1] / h + 0.5)

    if rotation != 0:
        rotation = -rotation
        r = np.eye(3)
        ang = rotation * math.pi / 180.0
        s = math.sin(ang)
        c = math.cos(ang)
        r[0][0] = c
        r[0][1] = -s
        r[1][0] = s
        r[1][1] = c

        t_ = np.eye(3)
        t_[0][2] = -resolution / 2.0
        t_[1][2] = -resolution / 2.0
        t_inv = torch.eye(3)
        t_inv[0][2] = resolution / 2.0
        t_inv[1][2] = resolution / 2.0
        t = reduce(np.matmul, [t_inv, r, t_, t])

    if invert:
        t = np.linalg.inv(t)
    new_point = (np.matmul(t, _pt))[0:2]

    return new_point.astype(int) 

def random_rotation_matrix():
    theta = np.random.rand()
    r_x = np.array([1, 0, 0, 0, np.cos(theta), -np.sin(theta), 0, np.sin(theta), np.cos(theta)]).reshape([3, 3])
    theta = np.random.rand()
    r_y = np.array([np.cos(theta), 0, np.sin(theta), 0, 1, 0, -np.sin(theta), 0, np.cos(theta)]).reshape([3, 3])
    theta = np.random.rand()
    r_z = np.array([np.cos(theta), -np.sin(theta), 0, np.sin(theta), np.cos(theta), 0, 0, 0, 1]).reshape([3, 3])

    return np.matmul(np.matmul(r_x, r_y), r_z) 

def tensorize(self, batch_x, batch_c):
        atom_tensor = np.zeros((len(batch_x), self.num_atoms, self.get_num_features()))
        adjm_tensor = np.zeros((len(batch_x), self.num_atoms, self.num_atoms))
        posn_tensor = np.zeros((len(batch_x), self.num_atoms, self.num_atoms, 3))

        for mol_idx, mol in enumerate(batch_x):
            Chem.RemoveHs(mol)
            mol_atoms = mol.GetNumAtoms()

            # Atom features
            atom_tensor[mol_idx, :mol_atoms, :] = self.get_atom_features(mol)

            # Adjacency matrix
            adjms = np.array(rdmolops.GetAdjacencyMatrix(mol), dtype="float")

            # Normalize adjacency matrix by D^(-1/2) * A_hat * D^(-1/2), Kipf et al. 2016
            adjms += np.eye(mol_atoms)
            degree = np.array(adjms.sum(1))
            deg_inv_sqrt = np.power(degree, -0.5)
            deg_inv_sqrt[np.isinf(deg_inv_sqrt)] = 0.
            deg_inv_sqrt = np.diag(deg_inv_sqrt)

            adjms = np.matmul(np.matmul(deg_inv_sqrt, adjms), deg_inv_sqrt)

            adjm_tensor[mol_idx, : mol_atoms, : mol_atoms] = adjms

            # Relative position matrix
            for atom_idx in range(mol_atoms):
                pos_c = batch_c[mol_idx][atom_idx]

                for neighbor_idx in range(mol_atoms):
                    pos_n = batch_c[mol_idx][neighbor_idx]

                    # Direction should be Neighbor -> Center
                    n_to_c = [pos_c[0] - pos_n[0], pos_c[1] - pos_n[1], pos_c[2] - pos_n[2]]
                    posn_tensor[mol_idx, atom_idx, neighbor_idx, :] = n_to_c

        return [atom_tensor, adjm_tensor, posn_tensor] 

def random_rotation_matrix():
    theta = np.random.rand() * 2 * np.pi
    r_x = np.array([1, 0, 0, 0, np.cos(theta), -np.sin(theta), 0, np.sin(theta), np.cos(theta)]).reshape([3, 3])
    theta = np.random.rand() * 2 * np.pi
    r_y = np.array([np.cos(theta), 0, np.sin(theta), 0, 1, 0, -np.sin(theta), 0, np.cos(theta)]).reshape([3, 3])
    theta = np.random.rand() * 2 * np.pi
    r_z = np.array([np.cos(theta), -np.sin(theta), 0, np.sin(theta), np.cos(theta), 0, 0, 0, 1]).reshape([3, 3])

    return np.matmul(np.matmul(r_x, r_y), r_z) 

def rotate_molecule(mol, axis, degree):
    rotation_matrix = degree_rotation_matrix(axis, float(degree))
    for atom in mol.GetAtoms():
        atom_idx = atom.GetIdx()
        pos = list(mol.GetConformer().GetAtomPosition(atom_idx))
        pos_rotated = np.matmul(rotation_matrix, pos)
        mol.GetConformer().SetAtomPosition(atom_idx, pos_rotated)
    return mol 

def transformPoint(point):
    """
    :param points: (numpy array)
    :return: (numpy array)
    """
    return np.matmul(complete_transform, point) 

def test_update_cholcovs_with_error(i, kalman_results):
    inp = kalman_results["cov_error"][i]["input"]
    calculated_chol = update_cholcov_with_measurement_error(**inp)
    expected_chol = kalman_results["cov_error"][i]["output"]

    calculated_cov = np.matmul(
        calculated_chol, np.transpose(calculated_chol, axes=(0, 2, 1))
    )
    expected_cov = np.matmul(expected_chol, np.transpose(expected_chol, axes=(0, 2, 1)))

    aaae(calculated_cov, expected_cov)