Python scipy.linalg.det() Examples

def f_value_wilks_lambda(ER, EF, dfnum, dfden, a, b):
    """Calculation of Wilks lambda F-statistic for multivarite data, per
    Maxwell & Delaney p.657.
    """
    if isinstance(ER, (int, float)):
        ER = array([[ER]])
    if isinstance(EF, (int, float)):
        EF = array([[EF]])
    lmbda = linalg.det(EF) / linalg.det(ER)
    if (a-1)**2 + (b-1)**2 == 5:
        q = 1
    else:
        q = np.sqrt(((a-1)**2*(b-1)**2 - 2) / ((a-1)**2 + (b-1)**2 - 5))
    n_um = (1 - lmbda**(1.0/q))*(a-1)*(b-1)
    d_en = lmbda**(1.0/q) / (n_um*q - 0.5*(a-1)*(b-1) + 1)
    return n_um / d_en 

def check_one_qubit_euler_angles(self, operator, basis='U3',
                                     tolerance=1e-12,
                                     phase_equal=False):
        """Check euler_angles_1q works for the given unitary"""
        decomposer = OneQubitEulerDecomposer(basis)
        with self.subTest(operator=operator):
            target_unitary = operator.data
            decomp_unitary = Operator(decomposer(operator)).data
            if not phase_equal:
                target_unitary *= la.det(target_unitary)**(-0.5)
                decomp_unitary *= la.det(decomp_unitary)**(-0.5)
            maxdist = np.max(np.abs(target_unitary - decomp_unitary))
            if not phase_equal and maxdist > 0.1:
                maxdist = np.max(np.abs(target_unitary + decomp_unitary))
            self.assertTrue(np.abs(maxdist) < tolerance, "Worst distance {}".format(maxdist))

    # U3 basis 

def _params_zyz(mat):
        """Return the euler angles and phase for the ZYZ basis."""
        # We rescale the input matrix to be special unitary (det(U) = 1)
        # This ensures that the quaternion representation is real
        coeff = la.det(mat)**(-0.5)
        phase = -np.angle(coeff)
        su_mat = coeff * mat  # U in SU(2)
        # OpenQASM SU(2) parameterization:
        # U[0, 0] = exp(-i(phi+lambda)/2) * cos(theta/2)
        # U[0, 1] = -exp(-i(phi-lambda)/2) * sin(theta/2)
        # U[1, 0] = exp(i(phi-lambda)/2) * sin(theta/2)
        # U[1, 1] = exp(i(phi+lambda)/2) * cos(theta/2)
        theta = 2 * math.atan2(abs(su_mat[1, 0]), abs(su_mat[0, 0]))
        phiplambda = 2 * np.angle(su_mat[1, 1])
        phimlambda = 2 * np.angle(su_mat[1, 0])
        phi = (phiplambda + phimlambda) / 2.0
        lam = (phiplambda - phimlambda) / 2.0
        return theta, phi, lam, phase 

def __init__(self, sigma, mu, seed=42):
        self.sigma = sigma
        self.mu = mu
        if not (len(mu.shape) == 1):
            raise Exception('mu has shape ' + str(mu.shape) +
                            ' (it should be a vector)')

        self.sigma_inv = solve(self.sigma, N.identity(mu.shape[0]),
                               sym_pos=True)
        self.L = cholesky(self.sigma)

        self.s_rng = make_theano_rng(seed, which_method='normal')

        #Compute logZ
        #log Z = log 1/( (2pi)^(-k/2) |sigma|^-1/2 )
        # = log 1 - log (2pi^)(-k/2) |sigma|^-1/2
        # = 0 - log (2pi)^(-k/2) - log |sigma|^-1/2
        # = (k/2) * log(2pi) + (1/2) * log |sigma|
        k = float(self.mu.shape[0])
        self.logZ = 0.5 * (k * N.log(2. * N.pi) + N.log(det(sigma))) 

def test_compute_log_det_cholesky():
    n_features = 2
    rand_data = RandomData(np.random.RandomState(0))

    for covar_type in COVARIANCE_TYPE:
        covariance = rand_data.covariances[covar_type]

        if covar_type == 'full':
            predected_det = np.array([linalg.det(cov) for cov in covariance])
        elif covar_type == 'tied':
            predected_det = linalg.det(covariance)
        elif covar_type == 'diag':
            predected_det = np.array([np.prod(cov) for cov in covariance])
        elif covar_type == 'spherical':
            predected_det = covariance ** n_features

        # We compute the cholesky decomposition of the covariance matrix
        expected_det = _compute_log_det_cholesky(_compute_precision_cholesky(
            covariance, covar_type), covar_type, n_features=n_features)
        assert_array_almost_equal(expected_det, - .5 * np.log(predected_det)) 

def f_value_wilks_lambda(ER, EF, dfnum, dfden, a, b):
    """Calculation of Wilks lambda F-statistic for multivarite data, per
    Maxwell & Delaney p.657.
    """
    if isinstance(ER, (int, float)):
        ER = array([[ER]])
    if isinstance(EF, (int, float)):
        EF = array([[EF]])
    lmbda = linalg.det(EF) / linalg.det(ER)
    if (a-1)**2 + (b-1)**2 == 5:
        q = 1
    else:
        q = np.sqrt(((a-1)**2*(b-1)**2 - 2) / ((a-1)**2 + (b-1)**2 - 5))

    n_um = (1 - lmbda**(1.0/q))*(a-1)*(b-1)
    d_en = lmbda**(1.0/q) / (n_um*q - 0.5*(a-1)*(b-1) + 1)
    return n_um / d_en 

def berry_phase(h,nk=20,kpath=None):
  """ Calculates the Berry phase of a Hamiltonian"""
  if h.dimensionality==0: raise
  elif h.dimensionality == 1:
    ks = np.linspace(0.,1.,nk,endpoint=False) # list of kpoints
  elif h.dimensionality > 1: # you must provide a kpath
      if kpath is None: 
          print("You must provide a k-path")
          raise # error
      ks = kpath # continue
      nk = len(kpath) # redefine
  else: raise # otherwise
  hkgen = h.get_hk_gen() # get Hamiltonian generator
  wf0 = occupied_states(hkgen,ks[0]) # get occupied states, first k-point
  wfold = wf0.copy() # copy
  m = np.matrix(np.identity(len(wf0))) # initialize as the identity matrix
  for ik in range(1,len(ks)): # loop over k-points, except first one
    wf = occupied_states(hkgen,ks[ik])  # get waves
    m = m@uij(wfold,wf)   # get the uij   and multiply
    wfold = wf.copy() # this is the new old
  m = m@uij(wfold,wf0)   # last one
  d = lg.det(m) # calculate determinant
  phi = np.arctan2(d.imag,d.real)
  open("BERRY_PHASE.OUT","w").write(str(phi/np.pi)+"\n")
  return phi # return Berry phase 

def berry_curvature(h,k,dk=0.01,window=None,max_waves=None):
  """ Calculates the Berry curvature of a 2d hamiltonian"""
  if h.dimensionality != 2: raise # only for 2d
  k = np.array([k[0],k[1]]) 
  dx = np.array([dk,0.])
  dy = np.array([0.,dk])
# get the function that returns the occ states
  occf = occ_states_generator(h,k,window=window,max_waves=max_waves)  
  # get the waves
#  print("Doing k-point",k)
  wf1 = occf(k-dx-dy) 
  wf2 = occf(k+dx-dy) 
  wf3 = occf(k+dx+dy) 
  wf4 = occf(k-dx+dy) 
  dims = [len(wf1),len(wf2),len(wf3),len(wf4)] # number of vectors
  if max(dims)!=min(dims): # check that the dimensions are fine 
#    print("WARNING, skipping this k-point",k)
    return 0.0 # if different number of vectors
  # get the uij  
  m = uij(wf1,wf2)@uij(wf2,wf3)@uij(wf3,wf4)@uij(wf4,wf1)
  d = lg.det(m) # calculate determinant
  phi = np.arctan2(d.imag,d.real)/(4.*dk*dk)
  return phi 

def is_rotation_matrix(mat, show_diff=False):
    from scipy.linalg import det, norm

    dim = mat.shape[0]
    if dim != mat.shape[1]:
        return False

    determ = det(mat)
    right_handed = (np.abs(determ - 1.) < 1E-10)
    orthonorm_diff = mat * mat.T - np.eye(dim)
    diff_norm = norm(orthonorm_diff, 2)
    orthonormal = (diff_norm < 1E-10)
    if not right_handed or not orthonormal:
        if show_diff:
            print('matrix S:\n', mat)
            print('det(S): ', determ)
            print('S*S.T - eye:\n', orthonorm_diff)
            print('2-norm of difference: ', diff_norm)
        return False
    return True 

def test_compute_log_det_cholesky():
    n_features = 2
    rand_data = RandomData(np.random.RandomState(0))

    for covar_type in COVARIANCE_TYPE:
        covariance = rand_data.covariances[covar_type]

        if covar_type == 'full':
            predected_det = np.array([linalg.det(cov) for cov in covariance])
        elif covar_type == 'tied':
            predected_det = linalg.det(covariance)
        elif covar_type == 'diag':
            predected_det = np.array([np.prod(cov) for cov in covariance])
        elif covar_type == 'spherical':
            predected_det = covariance ** n_features

        # We compute the cholesky decomposition of the covariance matrix
        expected_det = _compute_log_det_cholesky(_compute_precision_cholesky(
            covariance, covar_type), covar_type, n_features=n_features)
        assert_array_almost_equal(expected_det, - .5 * np.log(predected_det)) 

def _likelihood(ca, D, y_demeaned, n, model):
    R = model(D, ca)
    R_inv = la.inv(R)
    L = np.log(la.det(R)) + y_demeaned.dot(R_inv).dot(y_demeaned)
    return L


# Kernel reg of empirical data 

def f_value_multivariate(ER, EF, dfnum, dfden):
    """
    Returns a multivariate F-statistic.

    Parameters
    ----------
    ER : ndarray
        Error associated with the null hypothesis (the Restricted model).
        From a multivariate F calculation.
    EF : ndarray
        Error associated with the alternate hypothesis (the Full model)
        From a multivariate F calculation.
    dfnum : int
        Degrees of freedom the Restricted model.
    dfden : int
        Degrees of freedom associated with the Restricted model.

    Returns
    -------
    fstat : float
        The computed F-statistic.

    """
    if isinstance(ER, (int, float)):
        ER = array([[ER]])
    if isinstance(EF, (int, float)):
        EF = array([[EF]])
    n_um = (linalg.det(ER) - linalg.det(EF)) / float(dfnum)
    d_en = linalg.det(EF) / float(dfden)
    return n_um / d_en


#####################################
#         SUPPORT FUNCTIONS         #
##################################### 

def logsqrtdet(self):
    return (1/2)*np.log(la.det(self._Sigma)) 

def mlogdet(self):
        return np.log(linalg.det(self.m)) 

def compute_phi_prior(self):
        phi = 0
        for d in range(self.Dout):
            # s, a = npalg.slogdet(self.Kuu[d])
            a = np.log(npalg.det(self.Kuu[d]))
            phi += 0.5 * a
        return phi 

def mdet(self):
        return linalg.det(self.m) 

def detomega(self):
        r"""
        Return determinant of white noise covariance with degrees of freedom
        correction:

        .. math::

            \hat \Omega = \frac{T}{T - Kp - 1} \hat \Omega_{\mathrm{MLE}}
        """
        return L.det(self.sigma_u) 

def is_stable(self, verbose=False):
        """Determine stability based on model coefficients

        Parameters
        ----------
        verbose : bool
            Print eigenvalues of the VAR(1) companion

        Notes
        -----
        Checks if det(I - Az) = 0 for any mod(z) <= 1, so all the eigenvalues of
        the companion matrix must lie outside the unit circle
        """
        return is_stable(self.coefs, verbose=verbose) 

def compute_phi_posterior(self):
        phi = 0
        for d in range(self.Dout):
            mud_val = self.mu[d].get_value()
            Sud_val = self.Su[d].get_value()
            # s, a = npalg.slogdet(Sud_val)
            a = np.log(npalg.det(Sud_val))
            phi += 0.5 * np.dot(mud_val.T, np.dot(npalg.inv(Sud_val), mud_val))[0, 0]
            phi += 0.5 * a
        return phi 

def test_matdet_grad_1():
    npr.seed(1)

    for ii in xrange(NUM_TRIALS):
        
        np_A = npr.randn(12,6)
        A    = np.dot(np_A.T, np_A) + 1e-6*np.eye(6)
        B    = kayak.Parameter(A)
        D    = kayak.MatDet(B)

        assert_less((D.value - spla.det(A))**2, 1e-6)

        assert_equal(D.grad(B).shape, B.shape)
        assert_less(kayak.util.checkgrad(B, D), MAX_GRAD_DIFF) 

def compute_phi_cavity(self):
        phi = 0
        for d in range(self.Dout):
            muhatd_val = self.muhat[d].get_value()
            Suhatd_val = self.Suhat[d].get_value()
            # s, a = npalg.slogdet(Sud_val)
            a = np.log(npalg.det(Suhatd_val))
            phi += 0.5 * np.dot(muhatd_val.T, np.dot(npalg.inv(Suhatd_val), muhatd_val))[0, 0]
            phi += 0.5 * a
        return phi 

def test_matdet_values_1():
    npr.seed(1)

    for ii in xrange(NUM_TRIALS):
        
        np_A = npr.randn(12,6)
        A    = np.dot(np_A.T, np_A) + 1e-6*np.eye(6)
        B    = kayak.Parameter(A)
        D    = kayak.MatDet(B)

        assert_less((D.value - spla.det(A))**2, 1e-6) 

def test_det(self):
        """Quaternion det = 1"""
        assert_allclose(la.det(self.mat2), 1) 

def check_one_qubit_euler_angles(self, operator, tolerance=1e-14):
        """Check euler_angles_1q works for the given unitary"""
        with self.subTest(operator=operator):
            target_unitary = operator.data
            angles = euler_angles_1q(target_unitary)
            decomp_unitary = U3Gate(*angles).to_matrix()
            target_unitary *= la.det(target_unitary)**(-0.5)
            decomp_unitary *= la.det(decomp_unitary)**(-0.5)
            maxdist = np.max(np.abs(target_unitary - decomp_unitary))
            if maxdist > 0.1:
                maxdist = np.max(np.abs(target_unitary + decomp_unitary))
            self.assertTrue(np.abs(maxdist) < tolerance, "Worst distance {}".format(maxdist)) 

def logsqrtdet(self):
    # Sigma = L^T L -- but be careful: only the lower triangle and
    # diagonal of L are actually initialized; the upper triangle is
    # garbage.
    L, _lower = self._cholesky

    # det Sigma = det L^T L = det L^T det L = (det L)^2.  Since L is
    # triangular, its determinant is the product of its diagonal.  To
    # compute log sqrt(det Sigma) = log det L, we sum the logs of its
    # diagonal.
    return np.sum(np.log(np.diag(L))) 

def logpdf(X, Mu, Sigma):
  """Multivariate normal log pdf."""
  # This is the multivariate normal log pdf for an array X of n
  # outputs, an array Mu of n means, and an n-by-n positive-definite
  # covariance matrix Sigma.  The direct-space density is:
  #
  #     P(X | Mu, Sigma)
  #       = ((2 pi)^n det Sigma)^(-1/2)
  #         exp((-1/2) (X - Mu)^T Sigma^-1 (X - Mu)),
  #
  # We want this in log-space, so we have
  #
  #     log P(X | Mu, Sigma)
  #     = (-1/2) (X - Mu)^T Sigma^-1 (X - Mu) - log ((2 pi)^n det Sigma)^(1/2)
  #     = (-1/2) (X - Mu)^T Sigma^-1 (X - Mu)
  #         - (n/2) log (2 pi) - (1/2) log det Sigma.
  #
  n = len(X)
  assert X.shape == (n,)
  assert Mu.shape == (n,)
  assert Sigma.shape == (n, n)
  assert np.all(np.isfinite(X))
  assert np.all(np.isfinite(Mu))
  assert np.all(np.isfinite(Sigma))

  X_ = X - Mu
  covf = _covariance_factor(Sigma)

  logp = -np.dot(X_.T, covf.solve(X_)/2.)
  logp -= (n/2.)*np.log(2*np.pi)
  logp -= covf.logsqrtdet()

  # Convert 1x1 matrix to float.
  return float(logp) 

def _compute_covariance(self):
        """Computes the covariance matrix for each Gaussian kernel using
        covariance_factor().
        """
        self.factor = self.covariance_factor()
        # Cache covariance and inverse covariance of the data
        if not hasattr(self, '_data_inv_cov'):
            self._data_covariance = atleast_2d(np.cov(self.dataset, rowvar=1,
                                               bias=False))
            self._data_inv_cov = linalg.inv(self._data_covariance)

        self.covariance = self._data_covariance * self.factor**2
        self.inv_cov = self._data_inv_cov / self.factor**2
        self._norm_factor = sqrt(linalg.det(2*pi*self.covariance)) * self.n 

def mdet(self):
        return linalg.det(self.m) 

def mlogdet(self):
        return np.log(linalg.det(self.m)) 

def calcWeight( jac ):

  n = jac.shape

  if n[0] == n[1]:
    return abs(det(jac))
  elif n[0] == 2 and n[1] == 1:
    return sqrt(sum(sum(jac*jac)))

#------------------------------------------------------------------------------
#
#------------------------------------------------------------------------------