Python numpy.diagonal() Examples

The following are 30 code examples for showing how to use numpy.diagonal(). These examples are extracted from open source projects. 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 check out the related API usage on the sidebar.

You may also want to check out all available functions/classes of the module numpy , or try the search function .

Example 1
Project: vnpy_crypto   Author: birforce   File: utils.py    License: 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
Project: bhmm   Author: bhmm   File: gmm.py    License: 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
Project: tenpy   Author: tenpy   File: random_matrix.py    License: 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 4
Project: tenpy   Author: tenpy   File: random_matrix.py    License: 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 5
Project: tenpy   Author: tenpy   File: random_matrix.py    License: 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
Project: strawberryfields   Author: XanaduAI   File: states.py    License: 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 7
Project: Computable   Author: ktraunmueller   File: test_multiarray.py    License: 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 8
Project: POT   Author: PythonOT   File: plot_barycenter_fgw.py    License: 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 9
Project: trax   Author: google   File: array_ops.py    License: 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 10
Project: flare   Author: mir-group   File: gp_algebra.py    License: 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 11
Project: pyVSR   Author: georgesterpu   File: dct.py    License: 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 12
Project: fishervector   Author: jonasrothfuss   File: FisherVector.py    License: MIT License 5 votes vote down vote up
def _fit(self, X, model_dump_path=None, verbose=True):
    """
    :param X: shape (n_videos, n_frames, n_descriptors_per_image, n_dim_descriptor)
    :param model_dump_path: (optional) path where the fitted model shall be dumped
    :param verbose - boolean that controls the verbosity
    :return: fitted Fisher vector object
    """
    assert X.ndim == 4
    self.feature_dim = X.shape[-1]

    X = X.reshape(-1, X.shape[-1])

    # fit GMM and store params of fitted model
    self.gmm = gmm = GaussianMixture(n_components=self.n_kernels, covariance_type=self.covariance_type, max_iter=1000).fit(X)
    self.covars = gmm.covariances_
    self.means = gmm.means_
    self.weights = gmm.weights_

    # if cov_type is diagonal - make sure that covars holds a diagonal matrix
    if self.covariance_type == 'diag':
      cov_matrices = np.empty(shape=(self.n_kernels, self.covars.shape[1], self.covars.shape[1]))
      for i in range(self.n_kernels):
        cov_matrices[i, :, :] = np.diag(self.covars[i, :])
      self.covars = cov_matrices

    assert self.covars.ndim == 3
    self.fitted = True
    if verbose:
      print('fitted GMM with %i kernels'%self.n_kernels)

    if model_dump_path:
      with open(model_dump_path, 'wb') as f:
        pickle.dump(self,f, protocol=4)
      if verbose:
        print('Dumped fitted model to', model_dump_path)

    return self 
Example 13
Project: pyscf   Author: pyscf   File: uccsd_slow.py    License: Apache License 2.0 5 votes vote down vote up
def init_amps(self, eris):
        time0 = time.clock(), time.time()
        mo_e = eris.fock.diagonal()
        nocc = self.nocc
        eia = mo_e[:nocc,None] - mo_e[None,nocc:]
        eijab = lib.direct_sum('ia,jb->ijab',eia,eia)
        t1 = eris.fock[:nocc,nocc:] / eia
        eris_oovv = np.array(eris.oovv)
        t2 = eris_oovv/eijab
        self.emp2 = 0.25*einsum('ijab,ijab',t2,eris_oovv.conj()).real
        logger.info(self, 'Init t2, MP2 energy = %.15g', self.emp2)
        logger.timer(self, 'init mp2', *time0)
        return self.emp2, t1, t2 
Example 14
Project: respy   Author: OpenSourceEconomics   File: conditional_draws.py    License: MIT License 5 votes vote down vote up
def make_cholesky_unique(chol):
    """Make a lower triangular cholesky factor unique.

    Cholesky factors are only unique with the additional requirement that all diagonal
    elements are positive. This is done automatically by np.linalg.cholesky.
    Since we calucate cholesky factors by QR decompositions we have to do it manually.

    It is obvious from that this is admissible because:

    chol sign_swither sign_switcher.T chol.T = chol chol.T

    """
    sign_switcher = np.sign(np.diag(np.diagonal(chol)))
    return chol @ sign_switcher 
Example 15
Project: recruit   Author: Frank-qlu   File: test_numeric.py    License: Apache License 2.0 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
Project: recruit   Author: Frank-qlu   File: test_numeric.py    License: Apache License 2.0 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 17
Project: ibllib   Author: int-brain-lab   File: cca.py    License: MIT License 5 votes vote down vote up
def get_correlations(cca, data_0, data_1):
    """

    :param cca:
    :param data_0:
    :param data_1:
    :return:
    """
    x_scores, y_scores = get_cca_projection(cca, data_0, data_1)
    corrs_tmp = np.corrcoef(x_scores.T, y_scores.T)
    corrs = np.diagonal(corrs_tmp, offset=data_0.shape[1])
    return corrs 
Example 18
Project: pase   Author: santi-pdp   File: prep_voxforge.py    License: MIT License 5 votes vote down vote up
def compute_nrg(xframes):
	n_frames = xframes.shape[1]
	return np.diagonal(np.dot(xframes,xframes.T))/float(n_frames) 
Example 19
Project: auto-alt-text-lambda-api   Author: abhisuri97   File: test_numeric.py    License: 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 20
Project: auto-alt-text-lambda-api   Author: abhisuri97   File: test_numeric.py    License: 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 21
Project: KPConv   Author: HuguesTHOMAS   File: metrics.py    License: MIT License 5 votes vote down vote up
def metrics(confusions, ignore_unclassified=False):
    """
    Computes different metrics from confusion matrices.
    :param confusions: ([..., n_c, n_c] np.int32). Can be any dimension, the confusion matrices should be described by
    the last axes. n_c = number of classes
    :param ignore_unclassified: (bool). True if the the first class should be ignored in the results
    :return: ([..., n_c] np.float32) precision, recall, F1 score, IoU score
    """

    # If the first class (often "unclassified") should be ignored, erase it from the confusion.
    if (ignore_unclassified):
        confusions[..., 0, :] = 0
        confusions[..., :, 0] = 0

    # Compute TP, FP, FN. This assume that the second to last axis counts the truths (like the first axis of a
    # confusion matrix), and that the last axis counts the predictions (like the second axis of a confusion matrix)
    TP = np.diagonal(confusions, axis1=-2, axis2=-1)
    TP_plus_FP = np.sum(confusions, axis=-1)
    TP_plus_FN = np.sum(confusions, axis=-2)

    # Compute precision and recall. This assume that the second to last axis counts the truths (like the first axis of
    # a confusion matrix), and that the last axis counts the predictions (like the second axis of a confusion matrix)
    PRE = TP / (TP_plus_FN + 1e-6)
    REC = TP / (TP_plus_FP + 1e-6)

    # Compute Accuracy
    ACC = np.sum(TP, axis=-1) / (np.sum(confusions, axis=(-2, -1)) + 1e-6)

    # Compute F1 score
    F1 = 2 * TP / (TP_plus_FP + TP_plus_FN + 1e-6)

    # Compute IoU
    IoU = F1 / (2 - F1)

    return PRE, REC, F1, IoU, ACC 
Example 22
Project: KPConv   Author: HuguesTHOMAS   File: metrics.py    License: MIT License 5 votes vote down vote up
def IoU_from_confusions(confusions):
    """
    Computes IoU from confusion matrices.
    :param confusions: ([..., n_c, n_c] np.int32). Can be any dimension, the confusion matrices should be described by
    the last axes. n_c = number of classes
    :param ignore_unclassified: (bool). True if the the first class should be ignored in the results
    :return: ([..., n_c] np.float32) IoU score
    """

    # Compute TP, FP, FN. This assume that the second to last axis counts the truths (like the first axis of a
    # confusion matrix), and that the last axis counts the predictions (like the second axis of a confusion matrix)
    TP = np.diagonal(confusions, axis1=-2, axis2=-1)
    TP_plus_FN = np.sum(confusions, axis=-1)
    TP_plus_FP = np.sum(confusions, axis=-2)

    # Compute IoU
    IoU = TP / (TP_plus_FP + TP_plus_FN - TP + 1e-6)

    # Compute mIoU with only the actual classes
    mask = TP_plus_FN < 1e-3
    counts = np.sum(1 - mask, axis=-1, keepdims=True)
    mIoU = np.sum(IoU, axis=-1, keepdims=True) / (counts + 1e-6)

    # If class is absent, place mIoU in place of 0 IoU to get the actual mean later
    IoU += mask * mIoU

    return IoU 
Example 23
Project: vnpy_crypto   Author: birforce   File: test_numeric.py    License: 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 24
Project: vnpy_crypto   Author: birforce   File: test_numeric.py    License: 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 25
Project: vnpy_crypto   Author: birforce   File: markov_switching.py    License: MIT License 5 votes vote down vote up
def expected_durations(self):
        """
        (array) Expected duration of a regime, possibly time-varying.
        """
        return 1. / (1 - np.diagonal(self.regime_transition).squeeze()) 
Example 26
Project: vnpy_crypto   Author: birforce   File: linalg_covmat.py    License: MIT License 5 votes vote down vote up
def logpdf_obs(self, x):
        x = x - self.mean
        x_whitened = self.whiten(x)

        #sigmainv = linalg.cholesky(sigma)
        logdetsigma = np.log(np.linalg.det(sigma))

        sigma2 = 1. # error variance is included in sigma

        llike  =  0.5 * (np.log(sigma2)
                         - 2.* np.log(np.diagonal(self.cholsigmainv))
                         + (x_whitened**2)/sigma2
                         +  np.log(2*np.pi))

        return llike 
Example 27
Project: vnpy_crypto   Author: birforce   File: linalg_covmat.py    License: MIT License 5 votes vote down vote up
def mvn_nloglike_obs(x, sigma):
    '''loglike multivariate normal

    assumes x is 1d, (nobs,) and sigma is 2d (nobs, nobs)

    brute force from formula
    no checking of correct inputs
    use of inv and log-det should be replace with something more efficient
    '''
    #see numpy thread
    #Sturla: sqmahal = (cx*cho_solve(cho_factor(S),cx.T).T).sum(axis=1)

    #Still wasteful to calculate pinv first
    sigmainv = linalg.inv(sigma)
    cholsigmainv = linalg.cholesky(sigmainv)
    #2 * np.sum(np.log(np.diagonal(np.linalg.cholesky(A)))) #Dag mailinglist
    # logdet not needed ???
    #logdetsigma = 2 * np.sum(np.log(np.diagonal(cholsigmainv)))
    x_whitened = np.dot(cholsigmainv, x)

    #sigmainv = linalg.cholesky(sigma)
    logdetsigma = np.log(np.linalg.det(sigma))

    sigma2 = 1. # error variance is included in sigma

    llike  =  0.5 * (np.log(sigma2) - 2.* np.log(np.diagonal(cholsigmainv))
                          + (x_whitened**2)/sigma2
                          +  np.log(2*np.pi))

    return llike, (x_whitened**2) 
Example 28
Project: vnpy_crypto   Author: birforce   File: try_mlecov.py    License: MIT License 5 votes vote down vote up
def mvn_loglike_chol(x, sigma):
    '''loglike multivariate normal

    assumes x is 1d, (nobs,) and sigma is 2d (nobs, nobs)

    brute force from formula
    no checking of correct inputs
    use of inv and log-det should be replace with something more efficient
    '''
    #see numpy thread
    #Sturla: sqmahal = (cx*cho_solve(cho_factor(S),cx.T).T).sum(axis=1)
    sigmainv = np.linalg.inv(sigma)
    cholsigmainv = np.linalg.cholesky(sigmainv).T
    x_whitened = np.dot(cholsigmainv, x)

    logdetsigma = np.log(np.linalg.det(sigma))
    nobs = len(x)
    from scipy import stats
    print('scipy.stats')
    print(np.log(stats.norm.pdf(x_whitened)).sum())

    llf = - np.dot(x_whitened.T, x_whitened)
    llf -= nobs * np.log(2 * np.pi)
    llf -= logdetsigma
    llf *= 0.5
    return llf, logdetsigma, 2 * np.sum(np.log(np.diagonal(cholsigmainv)))
#0.5 * np.dot(x_whitened.T, x_whitened) + nobs * np.log(2 * np.pi) + logdetsigma) 
Example 29
Project: vnpy_crypto   Author: birforce   File: try_mlecov.py    License: MIT License 5 votes vote down vote up
def mvn_nloglike_obs(x, sigma):
    '''loglike multivariate normal

    assumes x is 1d, (nobs,) and sigma is 2d (nobs, nobs)

    brute force from formula
    no checking of correct inputs
    use of inv and log-det should be replace with something more efficient
    '''
    #see numpy thread
    #Sturla: sqmahal = (cx*cho_solve(cho_factor(S),cx.T).T).sum(axis=1)

    #Still wasteful to calculate pinv first
    sigmainv = np.linalg.inv(sigma)
    cholsigmainv = np.linalg.cholesky(sigmainv).T
    #2 * np.sum(np.log(np.diagonal(np.linalg.cholesky(A)))) #Dag mailinglist
    # logdet not needed ???
    #logdetsigma = 2 * np.sum(np.log(np.diagonal(cholsigmainv)))
    x_whitened = np.dot(cholsigmainv, x)

    #sigmainv = linalg.cholesky(sigma)
    logdetsigma = np.log(np.linalg.det(sigma))

    sigma2 = 1. # error variance is included in sigma

    llike  =  0.5 * (np.log(sigma2) - 2.* np.log(np.diagonal(cholsigmainv))
                          + (x_whitened**2)/sigma2
                          +  np.log(2*np.pi))

    return llike 
Example 30
Project: strawberryfields   Author: XanaduAI   File: states.py    License: 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