Python numpy.diagonal() Examples

The following are 30 code examples of numpy.diagonal(). 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: utils.py    From vnpy_crypto with MIT License 7 votes vote down vote up
def log_multivariate_normal_density(X, means, covars, min_covar=1.e-7):
    """Log probability for full covariance matrices. """
    if hasattr(linalg, 'solve_triangular'):
        # only in scipy since 0.9
        solve_triangular = linalg.solve_triangular
    else:
        # slower, but works
        solve_triangular = linalg.solve
    n_samples, n_dim = X.shape
    nmix = len(means)
    log_prob = np.empty((n_samples, nmix))
    for c, (mu, cv) in enumerate(zip(means, covars)):
        try:
            cv_chol = linalg.cholesky(cv, lower=True)
        except linalg.LinAlgError:
            # The model is most probabily stuck in a component with too
            # few observations, we need to reinitialize this components
            cv_chol = linalg.cholesky(cv + min_covar * np.eye(n_dim),
                                      lower=True)
        cv_log_det = 2 * np.sum(np.log(np.diagonal(cv_chol)))
        cv_sol = solve_triangular(cv_chol, (X - mu).T, lower=True).T
        log_prob[:, c] = - .5 * (np.sum(cv_sol ** 2, axis=1) + \
                                     n_dim * np.log(2 * np.pi) + cv_log_det)

    return log_prob 
Example #2
Source File: gmm.py    From bhmm with GNU Lesser General Public License v3.0 7 votes vote down vote up
def _log_multivariate_normal_density_full(X, means, covars, min_covar=1.e-7):
    """Log probability for full covariance matrices.
    """
    n_samples, n_dim = X.shape
    nmix = len(means)
    log_prob = np.empty((n_samples, nmix))
    for c, (mu, cv) in enumerate(zip(means, covars)):
        try:
            cv_chol = linalg.cholesky(cv, lower=True)
        except linalg.LinAlgError:
            # The model is most probably stuck in a component with too
            # few observations, we need to reinitialize this components
            cv_chol = linalg.cholesky(cv + min_covar * np.eye(n_dim),
                                      lower=True)
        cv_log_det = 2 * np.sum(np.log(np.diagonal(cv_chol)))
        cv_sol = linalg.solve_triangular(cv_chol, (X - mu).T, lower=True).T
        log_prob[:, c] = - .5 * (np.sum(cv_sol ** 2, axis=1) +
                                 n_dim * np.log(2 * np.pi) + cv_log_det)

    return log_prob 
Example #3
Source File: array_ops.py    From trax with Apache License 2.0 6 votes vote down vote up
def diagflat(v, k=0):
  """Returns a 2-d array with flattened `v` as diagonal.

  Args:
    v: array_like of any rank. Gets flattened when setting as diagonal. Could be
      an ndarray, a Tensor or any object that can be converted to a Tensor using
      `tf.convert_to_tensor`.
    k: Position of the diagonal. Defaults to 0, the main diagonal. Positive
      values refer to diagonals shifted right, negative values refer to
      diagonals shifted left.

  Returns:
    2-d ndarray.
  """
  v = asarray(v)
  return diag(tf.reshape(v.data, [-1]), k) 
Example #4
Source File: plot_barycenter_fgw.py    From POT with MIT License 6 votes vote down vote up
def sp_to_adjency(C, threshinf=0.2, threshsup=1.8):
    """ Thresholds the structure matrix in order to compute an adjency matrix.
    All values between threshinf and threshsup are considered representing connected nodes and set to 1. Else are set to 0
    Parameters
    ----------
    C : ndarray, shape (n_nodes,n_nodes)
        The structure matrix to threshold
    threshinf : float
        The minimum value of distance from which the new value is set to 1
    threshsup : float
        The maximum value of distance from which the new value is set to 1
    Returns
    -------
    C : ndarray, shape (n_nodes,n_nodes)
        The threshold matrix. Each element is in {0,1}
    """
    H = np.zeros_like(C)
    np.fill_diagonal(H, np.diagonal(C))
    C = C - H
    C = np.minimum(np.maximum(C, threshinf), threshsup)
    C[C == threshsup] = 0
    C[C != 0] = 1

    return C 
Example #5
Source File: random_matrix.py    From tenpy with GNU General Public License v3.0 6 votes vote down vote up
def O_close_1(size, a=0.01):
    r"""return an random orthogonal matrix 'close' to the Identity.

    Parameters
    ----------
    size : tuple
        ``(n, n)``, where `n` is the dimension of the output matrix.
    a : float
        Parameter determining how close the result is on `O`;
        :math:`\lim_{a \rightarrow 0} <|O-E|>_a = 0`` (where `E` is the identity).

    Returns
    -------
    O : ndarray
        Orthogonal matrix close to the identiy (for small `a`).
    """
    n, m = size
    assert n == m
    A = GOE(size) / (2. * n)**0.5  # scale such that eigenvalues are in [-1, 1]
    E = np.eye(size[0])
    Q, R = np.linalg.qr(E + a * A)
    L = np.diagonal(R)  # make QR decomposition unique & ensure Q is close to one for small `a`
    Q *= np.sign(L)
    return Q 
Example #6
Source File: random_matrix.py    From tenpy with GNU General Public License v3.0 6 votes vote down vote up
def CUE(size):
    r"""Circular unitary ensemble (CUE).

    Parameters
    ----------
    size : tuple
        ``(n, n)``, where `n` is the dimension of the output matrix.

    Returns
    -------
    U : ndarray
        Unitary matrix drawn from the CUE (=Haar measure on U(n)).
    """
    # almost same code as for CRE
    n, m = size
    assert n == m  # ensure that `mode` in qr doesn't matter.
    A = standard_normal_complex(size)
    Q, R = np.linalg.qr(A)
    # Q-R is not unique; to make it unique ensure that the diagonal of R is positive
    # Q' = Q*L; R' = L^{-1} *R, where L = diag(phase(diagonal(R)))
    L = np.diagonal(R).copy()
    L[np.abs(L) < 1.e-15] = 1.
    Q *= L / np.abs(L)
    return Q 
Example #7
Source File: random_matrix.py    From tenpy with GNU General Public License v3.0 6 votes vote down vote up
def CRE(size):
    r"""Circular real ensemble (CRE).

    Parameters
    ----------
    size : tuple
        ``(n, n)``, where `n` is the dimension of the output matrix.

    Returns
    -------
    U : ndarray
        Orthogonal matrix drawn from the CRE (=Haar measure on O(n)).
    """
    # almost same code as for CUE
    n, m = size
    assert n == m  # ensure that `mode` in qr doesn't matter.
    A = np.random.standard_normal(size)
    Q, R = np.linalg.qr(A)
    # Q-R is not unique; to make it unique ensure that the diagonal of R is positive
    # Q' = Q*L; R' = L^{-1} *R, where L = diag(phase(diagonal(R)))
    L = np.diagonal(R)
    Q *= np.sign(L)
    return Q 
Example #8
Source File: test_multiarray.py    From Computable with MIT License 6 votes vote down vote up
def test_diagonal(self):
        a = np.arange(12).reshape((3, 4))
        assert_equal(a.diagonal(), [0, 5, 10])
        assert_equal(a.diagonal(0), [0, 5, 10])
        assert_equal(a.diagonal(1), [1, 6, 11])
        assert_equal(a.diagonal(-1), [4, 9])

        b = np.arange(8).reshape((2, 2, 2))
        assert_equal(b.diagonal(), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(0), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(1), [[2], [3]])
        assert_equal(b.diagonal(-1), [[4], [5]])
        assert_raises(ValueError, b.diagonal, axis1=0, axis2=0)
        assert_equal(b.diagonal(0, 1, 2), [[0, 3], [4, 7]])
        assert_equal(b.diagonal(0, 0, 1), [[0, 6], [1, 7]])
        assert_equal(b.diagonal(offset=1, axis1=0, axis2=2), [[1], [3]])
        # Order of axis argument doesn't matter:
        assert_equal(b.diagonal(0, 2, 1), [[0, 3], [4, 7]]) 
Example #9
Source File: gp_algebra.py    From flare with MIT License 6 votes vote down vote up
def get_like_from_mats(ky_mat, l_mat, alpha, name):
    """ compute the likelihood from the covariance matrix

    :param ky_mat: the covariance matrix

    :return: float, likelihood
    """
    # catch linear algebra errors
    labels = _global_training_labels[name]

    # calculate likelihood
    like = (-0.5 * np.matmul(labels, alpha) -
            np.sum(np.log(np.diagonal(l_mat))) -
            math.log(2 * np.pi) * ky_mat.shape[1] / 2)

    return like


#######################################
##### KY MATRIX FUNCTIONS and gradients
####################################### 
Example #10
Source File: dct.py    From pyVSR with GNU General Public License v3.0 6 votes vote down vote up
def zz(matrix, nb):
    r"""Zig-zag traversal of the input matrix
    :param matrix: input matrix
    :param nb: number of coefficients to keep
    :return: an array of nb coefficients
    """
    flipped = np.fliplr(matrix)
    rows, cols = flipped.shape  # nb of columns

    coefficient_list = []

    for loop, i in enumerate(range(cols - 1, -rows, -1)):
        anti_diagonal = np.diagonal(flipped, i)

        # reversing even diagonals prioritizes the X resolution
        # reversing odd diagonals prioritizes the Y resolution
        # for square matrices, the information content is the same only when nb covers half of the matrix
        #  e.g. [ nb = n*(n+1)/2 ]
        if loop % 2 == 0:
            anti_diagonal = anti_diagonal[::-1]  # reverse anti_diagonal

        coefficient_list.extend([x for x in anti_diagonal])

    # flattened = [val for sublist in coefficient_list for val in sublist]
    return coefficient_list[:nb] 
Example #11
Source File: states.py    From strawberryfields with Apache License 2.0 6 votes vote down vote up
def cov(self):
        r"""The covariance matrix describing the Gaussian state.

        The diagonal elements of the covariance matrix correspond to the
        variance in the position and momentum quadratures:

        .. math::
            \mathbf{V}_{ii} = \begin{cases}
                (\Delta x_i)^2, & 0\leq i\leq N-1\\
                (\Delta p_{i-N})^2, & N\leq i\leq 2(N-1)
            \end{cases}

        where :math:`\Delta x_i` and :math:`\Delta p_i` refer to the
        position and momentum quadrature variance of mode :math:`i` respectively.

        Note that if the covariance matrix is purely diagonal, then this
        corresponds to squeezing :math:`z=re^{i\phi}` where :math:`\phi=0`,
        and :math:`\Delta x_i = e^{-2r}`, :math:`\Delta p_i = e^{2r}`.

        Returns:
          array: the :math:`2N\times 2N` covariance matrix.
        """
        return self._cov 
Example #12
Source File: patch_prior.py    From eht-imaging with GNU General Public License v3.0 5 votes vote down vote up
def loggausspdf2(X, sigma):
#log pdf of Gaussian with zero mena
#Based on code written by Mo Chen (mochen@ie.cuhk.edu.hk). March 2009.
    d = X.shape[0]

    R = np.linalg.cholesky(sigma).T;
    # todo check that sigma is psd

    q = np.sum( ( np.dot( np.linalg.inv(np.transpose(R)) , X ) )**2 , 0);  # quadratic term (M distance)
    c = d*np.log(2*np.pi)+2*np.sum(np.log( np.diagonal(R) ), 0);   # normalization constant
    y = -(c+q)/2.0;

    return y 
Example #13
Source File: mtag.py    From mtag with GNU General Public License v3.0 5 votes vote down vote up
def _posDef_adjustment(mat, scaling_factor=0.99,max_it=1000):
    '''
    Checks whether the provided is pos semidefinite. If it is not, then it performs the the adjustment procedure descried in 1.2.2 of the Supplementary Note

    scaling_factor: the multiplicative factor that all off-diagonal elements of the matrix are scaled by in the second step of the procedure.
    max_it: max number of iterations set so that
    '''
    logging.info('Checking for positive definiteness ..')
    assert mat.ndim == 2
    assert mat.shape[0] == mat.shape[1]
    is_pos_semidef = lambda m: np.all(np.linalg.eigvals(m) >= 0)
    if is_pos_semidef(mat):
        return mat
    else:
        logging.info('matrix is not positive definite, performing adjustment..')
        P = mat.shape[0]
        for i in range(P):
            for j in range(i,P):
                if np.abs(mat[i,j]) > np.sqrt(mat[i,i] * mat[j,j]):
                    mat[i,j] = scaling_factor*np.sign(mat[i,j])*np.sqrt(mat[i,i] * mat[j,j])
                    mat[j,i] = mat[i,j]
        n=0
        while not is_pos_semidef(mat) and n < max_it:
            dg = np.diag(mat)
            mat = scaling_factor * mat
            mat[np.diag_indices(P)] = dg
            n += 1
        if n == max_it:
            logging.info('Warning: max number of iterations reached in adjustment procedure. Sigma matrix used is still non-positive-definite.')
        else:
            logging.info('Completed in {} iterations'.format(n))
        return mat 
Example #14
Source File: test_numeric.py    From GraphicDesignPatternByPython with MIT License 5 votes vote down vote up
def test_diagonal(self):
        a = [[0, 1, 2, 3],
             [4, 5, 6, 7],
             [8, 9, 10, 11]]
        out = np.diagonal(a)
        tgt = [0, 5, 10]

        assert_equal(out, tgt) 
Example #15
Source File: test_numeric.py    From GraphicDesignPatternByPython with MIT License 5 votes vote down vote up
def test_diagonal():
    b1 = np.matrix([[1,2],[3,4]])
    diag_b1 = np.matrix([[1, 4]])
    array_b1 = np.array([1, 4])

    assert_equal(b1.diagonal(), diag_b1)
    assert_equal(np.diagonal(b1), array_b1)
    assert_equal(np.diag(b1), array_b1) 
Example #16
Source File: gp_algebra.py    From flare with MIT License 5 votes vote down vote up
def get_like_grad_from_mats(ky_mat, hyp_mat, name):
    """compute the gradient of likelihood to hyper-parameters
    from covariance matrix and its gradient

    :param ky_mat: covariance matrix
    :type ky_mat: np.array
    :param hyp_mat: dky/d(hyper parameter) matrix
    :type hyp_mat: np.array
    :param name: name of the gp instance.

    :return: float, list. the likelihood and its gradients
    """

    number_of_hyps = hyp_mat.shape[0]

    # catch linear algebra errors
    try:
        ky_mat_inv = np.linalg.inv(ky_mat)
        l_mat = np.linalg.cholesky(ky_mat)
    except:
        return -1e8, np.zeros(number_of_hyps)

    labels = _global_training_labels[name]

    alpha = np.matmul(ky_mat_inv, labels)
    alpha_mat = np.matmul(alpha.reshape(-1, 1),
                          alpha.reshape(1, -1))
    like_mat = alpha_mat - ky_mat_inv

    # calculate likelihood
    like = (-0.5 * np.matmul(labels, alpha) -
            np.sum(np.log(np.diagonal(l_mat))) -
            math.log(2 * np.pi) * ky_mat.shape[1] / 2)

    # calculate likelihood gradient
    like_grad = np.zeros(number_of_hyps)
    for n in range(number_of_hyps):
        like_grad[n] = 0.5 * \
                       np.trace(np.matmul(like_mat, hyp_mat[n, :, :]))

    return like, like_grad 
Example #17
Source File: mtag.py    From mtag with GNU General Public License v3.0 5 votes vote down vote up
def flatten_out_omega(omega_est):
    # stacks the lower part of the cholesky decomposition ROW_WISE [(0,0) (1,0) (1,1) (2,0) (2,1) (2,2) ...]
    P_c = len(omega_est)
    x_chol = np.linalg.cholesky(omega_est)

    # transform components of cholesky decomposition for better optimization
    lowTr_ind = np.tril_indices(P_c)
    x_chol_trf = np.zeros((P_c,P_c))
    for i in range(P_c):
        for j in range(i): # fill in lower triangular components not on diagonal
            x_chol_trf[i,j] = x_chol[i,j]/np.sqrt(x_chol[i,i]*x_chol[j,j])
    x_chol_trf[np.diag_indices(P_c)] = np.log(np.diag(x_chol))  # replace with log transformation on the diagonal
    return tuple(x_chol_trf[lowTr_ind]) 
Example #18
Source File: mtag.py    From mtag with GNU General Public License v3.0 5 votes vote down vote up
def jointEffect_probability(Z_score, omega_hat, sigma_hat,N_mats, S=None):
    ''' For each SNP m in each state s , computes the evaluates the multivariate normal distribution at the observed row of Z-scores
    Calculate the distribution of (Z_m | s ) for all s in S, m in M. --> M  x|S| matrix
    The output is a M x n_S matrix of joint probabilities
    '''

    DTYPE = np.float64
    (M,P) = Z_score.shape
    if S is None: # 2D dimensional form
        assert omega_hat.ndim == 2
        omega_hat = omega_hat.reshape(1,P,P)
        S = np.ones((1,P),dtype=bool)

    (n_S,_) = S.shape
    jointProbs = np.empty((M,n_S))

    xRinvs = np.zeros([M,n_S,P], dtype=DTYPE)
    logSqrtDetSigmas = np.zeros([M,n_S], dtype=DTYPE)
    Ls = np.zeros([M,n_S,P,P], dtype=DTYPE)
    cov_s = np.zeros([M,n_S,P,P], dtype=DTYPE)

    Zs_rep = np.einsum('mp,s->msp',Z_score,np.ones(n_S))  # functionally equivalent to repmat
    cov_s = np.einsum('mpq,spq->mspq',N_mats,omega_hat) + sigma_hat

    Ls = np.linalg.cholesky(cov_s)
    Rs = np.transpose(Ls, axes=(0,1,3,2))

    xRinvs = np.linalg.solve(Ls, Zs_rep)

    logSqrtDetSigmas = np.sum(np.log(np.diagonal(Rs,axis1=2,axis2=3)),axis=2).reshape(M,n_S)

    quadforms = np.sum(xRinvs**2,axis=2).reshape(M,n_S)
    jointProbs = np.exp(-0.5 * quadforms - logSqrtDetSigmas - P * np.log(2 * np.pi) / 2)

    if n_S == 1:
        jointProbs = jointProbs.flatten()

    return jointProbs 
Example #19
Source File: evaluate_contact_prediction_metrics.py    From tape-neurips2019 with MIT License 5 votes vote down vote up
def partition_contacts(full_contact_map: np.ndarray):
    """ Returns lists of short, medium, and long-range contacts.
    """
    length = full_contact_map.shape[0]
    short_contacts = [np.diagonal(full_contact_map, i) for i in range(6,12)] 
    medium_contacts = [np.diagonal(full_contact_map, i) for i in range(12,min(length, 24))]
    if length >= 24:
        long_contacts = [np.diagonal(full_contact_map, i) for i in range(24, length)] 
    else:
        return np.concatenate(short_contacts), np.concatenate(medium_contacts), []
    return np.concatenate(short_contacts), np.concatenate(medium_contacts), np.concatenate(long_contacts) 
Example #20
Source File: states.py    From strawberryfields with Apache License 2.0 5 votes vote down vote up
def mean_photon(self, mode, **kwargs):
        # pylint: disable=unused-argument
        n = np.arange(self._cutoff)
        probs = np.diagonal(self.reduced_dm(mode))
        mean = np.sum(n * probs).real
        var = np.sum(n ** 2 * probs).real - mean ** 2
        return mean, var 
Example #21
Source File: tfdeploy.py    From tfdeploy with MIT License 5 votes vote down vote up
def MatrixDiagPart(a):
    """
    Batched diag op that returns only the diagonal elements.
    """
    r = np.zeros(a.shape[:-2] + (min(a.shape[-2:]),))
    for coord in np.ndindex(a.shape[:-2]):
        pos = coord + (Ellipsis,)
        r[pos] = np.diagonal(a[pos])
    return r, 
Example #22
Source File: tfdeploy.py    From tfdeploy with MIT License 5 votes vote down vote up
def DiagPart(a):
    """
    Diag op that returns only the diagonal elements.
    """
    return np.diagonal(a), 
Example #23
Source File: DataProcessor.py    From RaptorX-Contact with GNU General Public License v3.0 5 votes vote down vote up
def PriorDistancePotential(sequence=None, paramfile=None):
    	##add pairwise distance potential here
    	## an example paramfile is data4contact/pdb25-pair-dist.pkl

    	if not os.path.isfile(paramfile):
		print 'cannot find the parameter file: ', paramfile
		exit(-1)

    	fh = open(paramfile,'rb')
    	potential = cPickle.load(fh)[0].astype(np.float32)
    	fh.close()
	assert (len(potential.shape) == 4)

    	potentialFeature = np.zeros((len(sequence), len(sequence), potential.shape[-1]), dtype=theano.config.floatX)

	##convert AAs to integers
	ids = [ ord(AA) - ord('A') for AA in sequence ]

	##the below procedure is not very effective. What we can do is to generate a full matrix of only long-range potential using OuterConcatenate and np.choose
	##and then using the np.diagonal() function to replace near-, short- and medium-range potential in the full matrix
    	for i, id0 in zip(xrange(len(ids)), ids):
		for j, id1 in zip(xrange(i+1, len(ids)), ids[i+1:]):
	    		if j-i<6:
	        		sepIndex = 0
	    		elif j-i < 12:
	        		sepIndex = 1
	    		elif j-i < 24:
	        		sepIndex = 2
            		else:
	        		sepIndex = 3
		
			if id0 <=id1:
	    			potentialFeature[i][j]=potential[sepIndex][id0][id1]
			else:
	    			potentialFeature[i][j]=potential[sepIndex][id1][id0]
	    		potentialFeature[j][i]=potentialFeature[i][j]

    	return potentialFeature

##d is a dictionary for a protein 
Example #24
Source File: gmm.py    From bhmm with GNU Lesser General Public License v3.0 5 votes vote down vote up
def _covar_mstep_diag(gmm, X, responsibilities, weighted_X_sum, norm,
                      min_covar):
    """Performing the covariance M step for diagonal cases"""
    avg_X2 = np.dot(responsibilities.T, X * X) * norm
    avg_means2 = gmm.means_ ** 2
    avg_X_means = gmm.means_ * weighted_X_sum * norm
    return avg_X2 - 2 * avg_X_means + avg_means2 + min_covar 
Example #25
Source File: states.py    From strawberryfields with Apache License 2.0 5 votes vote down vote up
def diagonal_expectation(self, modes, values):
        """Calculates the expectation value of an operator that is diagonal in the number basis"""
        if len(modes) != len(set(modes)):
            raise ValueError("There can be no duplicates in the modes specified.")

        cutoff = self._cutoff  # Fock space cutoff.
        num_modes = self._modes  # number of modes in the state.

        traced_modes = tuple(item for item in range(num_modes) if item not in modes)
        if self.is_pure:
            # state is a tensor of probability amplitudes
            ps = np.abs(self.ket()) ** 2
            ps = ps.sum(axis=traced_modes)
            for _ in modes:
                ps = np.tensordot(values, ps, axes=1)
            return float(ps)

        # state is a tensor of density matrix elements in the SF convention
        ps = np.real(self.dm())
        traced_modes = list(traced_modes)
        traced_modes.sort(reverse=True)
        for mode in traced_modes:
            ps = np.tensordot(np.identity(cutoff), ps, axes=((0, 1), (2 * mode, 2 * mode + 1)))
        for _ in range(len(modes)):
            ps = np.tensordot(np.diag(values), ps, axes=((0, 1), (0, 1)))
        return float(ps) 
Example #26
Source File: gmm.py    From bhmm with GNU Lesser General Public License v3.0 5 votes vote down vote up
def _log_multivariate_normal_density_diag(X, means, covars):
    """Compute Gaussian log-density at X for a diagonal model"""
    n_samples, n_dim = X.shape
    lpr = -0.5 * (n_dim * np.log(2 * np.pi) + np.sum(np.log(covars), 1)
                  + np.sum((means ** 2) / covars, 1)
                  - 2 * np.dot(X, (means / covars).T)
                  + np.dot(X ** 2, (1.0 / covars).T))
    return lpr 
Example #27
Source File: lax_numpy_test.py    From trax with Apache License 2.0 5 votes vote down vote up
def testDiagonal(self, shape, dtype, offset, axis1, axis2, rng_factory):
    rng = rng_factory()
    onp_fun = lambda arg: onp.diagonal(arg, offset, axis1, axis2)
    lnp_fun = lambda arg: lnp.diagonal(arg, offset, axis1, axis2)
    args_maker = lambda: [rng(shape, dtype)]
    self._CheckAgainstNumpy(onp_fun, lnp_fun, args_maker, check_dtypes=True)
    self._CompileAndCheck(
        lnp_fun, args_maker, check_dtypes=True, check_incomplete_shape=True) 
Example #28
Source File: test_multiarray.py    From Computable with MIT License 5 votes vote down vote up
def test_diagonal_memleak(self):
        # Regression test for a bug that crept in at one point
        a = np.zeros((100, 100))
        assert_(sys.getrefcount(a) < 50)
        for i in range(100):
            a.diagonal()
        assert_(sys.getrefcount(a) < 50) 
Example #29
Source File: ops.py    From strawberryfields with Apache License 2.0 5 votes vote down vote up
def __init__(self, S, vacuum=False, tol=1e-10):
        super().__init__([S])
        self.ns = S.shape[0] // 2
        self.vacuum = (
            vacuum  #: bool: if True, ignore the first unitary matrix when applying the gate
        )
        N = self.ns  # shorthand

        # check if input symplectic is passive (orthogonal)
        diffn = np.linalg.norm(S @ S.T - np.identity(2 * N))
        self.active = (
            np.abs(diffn) > _decomposition_tol
        )  #: bool: S is an active symplectic transformation

        if not self.active:
            # The transformation is passive, do Clements
            X1 = S[:N, :N]
            P1 = S[N:, :N]
            self.U1 = X1 + 1j * P1
        else:
            # transformation is active, do Bloch-Messiah
            O1, smat, O2 = dec.bloch_messiah(S, tol=tol)
            X1 = O1[:N, :N]
            P1 = O1[N:, :N]
            X2 = O2[:N, :N]
            P2 = O2[N:, :N]

            self.U1 = X1 + 1j * P1  #: array[complex]: unitary matrix corresponding to O_1
            self.U2 = X2 + 1j * P2  #: array[complex]: unitary matrix corresponding to O_2
            self.Sq = np.diagonal(smat)[
                :N
            ]  #: array[complex]: diagonal vector of the squeezing matrix R 
Example #30
Source File: ops.py    From funsor with Apache License 2.0 5 votes vote down vote up
def _diagonal(x, dim1, dim2):
    return np.diagonal(x, axis1=dim1, axis2=dim2)