Python numpy.matmul() Examples

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 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 _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 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 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 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 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 _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 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 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 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 _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 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 _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 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 rotation_matrix_along_axis(x, angle):
    cx = np.cos(angle)
    sx = np.sin(angle)
    x_cpm = np.array([
        [0, -x[2], x[1]],
        [x[2], 0, -x[0]],
        [-x[1], x[0], 0]
    ], dtype='float')
    x = x.reshape(-1, 1)
    mat33_x = cx * np.eye(3) + sx * x_cpm + (1.0 - cx) * np.matmul(x, x.T)
    return mat33_x 

def compute_metrical_clip(self, pose, diameter):

        width = self.cam[0, 0] * diameter / pose[2, 3]  # X coordinate == shape[1]
        height = self.cam[1, 1] * diameter / pose[2, 3]  # Y coordinate == shape[0]
        proj = np.matmul(self.cam, pose[0:3, 3])
        proj /= proj[2]
        cut = np.asarray([proj[1] - height//2, proj[0] - width//2, proj[1] + height//2, proj[0] + width//2], dtype=int)

        # Can lead to offsetted extractions, not really nice...
        cut[0] = np.clip(cut[0], 0, self.shape[0])
        cut[2] = np.clip(cut[2], 0, self.shape[0])
        cut[1] = np.clip(cut[1], 0, self.shape[1])
        cut[3] = np.clip(cut[3], 0, self.shape[1])
        return cut 

def render_obj(self,ply_model,R=np.eye(3),t=np.array([0,0,1])):
        pts = (np.matmul(R,ply_model['pts'].T) + np.expand_dims(t,axis=1)).T #Nx3
        face_set = ply_model['faces'].astype(np.int)
        return self.cuda_render(pts,face_set) 

def _matvec(self, x):
        x = np.squeeze(x.reshape(self.nsl, self.ny, self.nz))
        if self.usematmul:
            if self.nz == 1:
                x = x[..., np.newaxis]
            y = np.matmul(self.G, x)
        else:
            y = np.squeeze(np.zeros((self.nsl, self.nx, self.nz),
                                    dtype=self.dtype))
            for isl in range(self.nsl):
                y[isl] = np.dot(self.G[isl], x[isl])
        return y.ravel() 

def test_power_is_minus_one(self, dt):
        def tz(mat):
            invmat = matrix_power(mat, -1)
            mmul = matmul if mat.dtype != object else dot
            assert_almost_equal(
                mmul(invmat, mat), identity_like_generalized(mat))

        for mat in self.rshft_all:
            if dt not in self.dtnoinv:
                tz(mat.astype(dt)) 

def test_power_is_two(self, dt):
        def tz(mat):
            mz = matrix_power(mat, 2)
            mmul = matmul if mat.dtype != object else dot
            assert_equal(mz, mmul(mat, mat))
            assert_equal(mz.dtype, mat.dtype)

        for mat in self.rshft_all:
            tz(mat.astype(dt))
            if dt != object:
                tz(self.stacked.astype(dt)) 

def test_matmul(self):
        array = np.array([1, 2], dtype=np.float64)
        array_like = ArrayLike(array)
        expected = ArrayLike(np.float64(5))
        _assert_equal_type_and_value(expected, np.matmul(array_like, array))
        if not PY2:
            _assert_equal_type_and_value(
                expected, operator.matmul(array_like, array))
            _assert_equal_type_and_value(
                expected, operator.matmul(array, array_like)) 

def precompute(self):
        # for quicker norm, sim
        self.qdims = np.arange(self.qdim)

        # codes  # qdim x centroids x -1
        # indexes  # vectors x qdim
        norm = np.sum(self.codes ** 2, axis=-1)  # qdim x centroids
        norm = np.sum(norm[self.qdims, self.indexes], axis=-1)  # vectors x 1
        self.norm = np.sqrt(norm)

        self.ab = np.matmul(
            self.codes,  # qdim x centroids x -1
            np.transpose(self.codes, axes=[0, 2, 1])  # qdim x -1 x centroids
        )  # qdim x centroids x centroids 

def rotate_coordinates(local3d, angles):
    """
    Rotate local rectangular coordinates from given view_angles.

    :param local3d: numpy array. Unit vectors for local rectangular coordinates's , shape (3, 3).
    :param angles: tuple of length 3. Rotation angles around each axis.
    :return:
    """
    cx, cy, cz = np.cos(angles)
    sx, sy, sz = np.sin(angles)

    x = local3d[0]
    x_cpm = np.array([
        [0, -x[2], x[1]],
        [x[2], 0, -x[0]],
        [-x[1], x[0], 0]
    ], dtype='float')
    x = x.reshape(-1, 1)
    mat33_x = cx * np.eye(3) + sx * x_cpm + (1.0 - cx) * np.matmul(x, x.T)

    mat33_z = np.array([
        [cz, sz, 0],
        [-sz, cz, 0],
        [0, 0, 1]
    ], dtype='float')

    local3d = local3d @ mat33_x.T @ mat33_z
    return local3d 

def runOptimiser(K, u, preOptw, initialValue, maxWeight=10000):
    """
    Args:
        K (double 2d array): Similarity/distance matrix
        u (double array): Mean similarity of each prototype
        preOptw (double): Weight vector
        initialValue (double): Initialize run
        maxWeight (double): Upper bound on weight

    Returns:
        Prototypes, weights and objective values
    """
    d = u.shape[0]
    lb = np.zeros((d, 1))
    ub = maxWeight * np.ones((d, 1))
    x0 = np.append( preOptw, initialValue/K[d-1, d-1] )

    G = np.vstack((np.identity(d), -1*np.identity(d)))
    h = np.vstack((ub, -1*lb))

    #     Solve a QP defined as follows:
    #         minimize
    #             (1/2) * x.T * P * x + q.T * x
    #         subject to
    #             G * x <= h
    #             A * x == b
    sol = solve_qp(K, -u, G, h, A=None, b=None, solver='cvxopt', initvals=x0)

    # compute objective function value
    x = sol.reshape(sol.shape[0], 1)
    P = K
    q = - u.reshape(u.shape[0], 1)
    obj_value = 1/2 * np.matmul(np.matmul(x.T, P), x) + np.matmul(q.T, x)
    return(sol, obj_value[0,0]) 

def execute(cls, ctx, op):
        (a, b), device_id, xp = as_same_device(
            [ctx[c.key] for c in op.inputs], device=op.device, ret_extra=True)

        with device(device_id):
            if not op.sparse and is_sparse_module(xp):
                # tell sparse to do calculation on numpy or cupy matmul
                ctx[op.outputs[0].key] = xp.matmul(a, b, sparse=False)
            else:
                try:
                    # `np.matmul` support `order` argument in version 1.16
                    ctx[op.outputs[0].key] = xp.matmul(a, b, casting=op.casting, order=op.order)
                except TypeError:  # pragma: no cover
                    ctx[op.outputs[0].key] = xp.matmul(a, b).astype(dtype=op.dtype,
                                                                    casting=op.casting, order=op.order)