Python numpy.diagflat() Examples

def kernel_qchem_inter_rf_pos_neg(self, **kw):
    """ This is constructing the E_m-E_n and E_n-E_m matrices """
    h_rpa = diagflat(concatenate((ravel(self.FmE),-ravel(self.FmE))))
    print(h_rpa.shape)

    nf = self.nfermi[0]
    nv = self.norbs-self.vstart[0]
    vs = self.vstart[0]
    neh = nf*nv
    x = self.mo_coeff[0,0,:,:,0]
    pab2v = self.pb.get_ac_vertex_array()
    self.pmn2v = pmn2v = einsum('nb,pmb->pmn', x[:nf,:], einsum('ma,pab->pmb', x[vs:,:], pab2v))
    pmn2c = einsum('qp,pmn->qmn', self.hkernel_den, pmn2v)
    meri = einsum('pmn,pik->mnik', pmn2c, pmn2v).reshape((nf*nv,nf*nv))
    #print(meri.shape)
    #meri.fill(0.0)
    h_rpa[:neh, :neh] = h_rpa[:neh, :neh]+meri
    h_rpa[:neh, neh:] = h_rpa[:neh, neh:]+meri
    h_rpa[neh:, :neh] = h_rpa[neh:, :neh]-meri
    h_rpa[neh:, neh:] = h_rpa[neh:, neh:]-meri
    edif, s2z = np.linalg.eig(h_rpa)
    print(abs(h_rpa-h_rpa.transpose()).sum())
    print('edif', edif.real*27.2114)
    
    return 

def learn_embedding(self):

        graph = self.g.G
        A = nx.to_numpy_matrix(graph)

        # self._beta = 0.0728

        # M_g = np.eye(graph.number_of_nodes()) - self._beta * A
        # M_l = self._beta * A

        M_g = np.eye(graph.number_of_nodes())
        M_l = np.dot(A, A)

        S = np.dot(np.linalg.inv(M_g), M_l)
        # s: \sigma_k
        u, s, vt = lg.svds(S, k=self._d // 2)
        sigma = np.diagflat(np.sqrt(s))
        X1 = np.dot(u, sigma)
        X2 = np.dot(vt.T, sigma)
        # self._X = X2
        self._X = np.concatenate((X1, X2), axis=1) 

def get_bounds_fit_matrix(xmin, xmax, ymin, ymax, zmin, zmax):
    """Create a 4x4 transform matrix which maps the default bounding box ([-0.5, 0.5] in all dimensions) into
    a custom bounding box ([xmin, xmax, ymin, ymax, zmin, zmax]).

    It is used for fitting fields (VectorField, Surface, MachingCubes) into custom domains."""

    for name, value in locals().items():
        try:
            float(value)
        except (TypeError, ValueError):
            raise TypeError('%s: expected float, %s given' % (name, type(value).__name__))

    matrix = np.diagflat(np.array((xmax - xmin, ymax - ymin, zmax - zmin, 1.0), np.float32, order='C'))
    matrix[0:3, 3] = ((xmax + xmin) / 2.0, (ymax + ymin) / 2.0, (zmax + zmin) / 2.0)

    return matrix 

def threeassetportfolio(plength=5000, SRlist=[1.0, 1.0, 1.0], annual_vol=.15, clist=[.0,.0,.0], index_start=pd.datetime(2000,1,1)):

    (c1, c2, c3)=clist
    dindex=arbitrary_timeindex(plength, index_start)

    daily_vol=annual_vol/16.0
    means=[x*annual_vol/250.0 for x in SRlist]
    stds = np.diagflat([daily_vol]*3)
    corr=np.array([[1.0, c1, c2], [c1, 1.0, c3], [c2, c3, 1.0]])
 
    covs=np.dot(stds, np.dot(corr, stds))
    plength=len(dindex)

    m = np.random.multivariate_normal(means, covs, plength).T

    portreturns=pd.DataFrame(dict(one=m[0], two=m[1], three=m[2]), dindex)
    portreturns=portreturns[['one', 'two', 'three']]
    
    return portreturns 

def diag_hamiltonian(self, subsystem, evals=None):
        """Returns a `qutip.Qobj` which has the eigenenergies of the object `subsystem` on the diagonal.

        Parameters
        ----------
        subsystem: object derived from `QuantumSystem`
            Subsystem for which the Hamiltonian is to be provided.
        evals: ndarray, optional
            Eigenenergies can be provided as `evals`; otherwise, they are calculated.

        Returns
        -------
        qutip.Qobj operator
        """
        evals_count = subsystem.truncated_dim
        if evals is None:
            evals = subsystem.eigenvals(evals_count=evals_count)
        diag_qt_op = qt.Qobj(inpt=np.diagflat(evals[0:evals_count]))
        return self.identity_wrap(diag_qt_op, subsystem) 

def number(dimension, prefactor=None):
    """Number operator matrix of size dimension x dimension in sparse matrix representation. An additional prefactor
    can be directly included in the generation of the matrix by supplying 'prefactor'.

    Parameters
    ----------
    dimension: int
    prefactor: float or complex, optional
        prefactor multiplying the number operator matrix


    Returns
    -------
    ndarray
        number operator matrix, size dimension x dimension
    """
    diag_elements = np.arange(dimension)
    if prefactor:
        diag_elements *= prefactor
    return np.diagflat(diag_elements) 

def _bwd(self, dLdy, X):
        """
        Actual computation of the gradient of the loss wrt. the input X.

        The Jacobian, J, of the softmax for input x = [x1, ..., xn] is:
            J[i, j] =
                softmax(x_i)  * (1 - softmax(x_j))  if i = j
                -softmax(x_i) * softmax(x_j)        if i != j
            where
                x_n is input example n (ie., the n'th row in X)
        """
        dX = []
        for dy, x in zip(dLdy, X):
            dxi = []
            for dyi, xi in zip(*np.atleast_2d(dy, x)):
                yi = self._fwd(xi.reshape(1, -1)).reshape(-1, 1)
                dyidxi = np.diagflat(yi) - yi @ yi.T  # jacobian wrt. input sample xi
                dxi.append(dyi @ dyidxi)
            dX.append(dxi)
        return np.array(dX).reshape(*X.shape) 

def learn_embedding(self):

        graph = self.g.G
        A = nx.to_numpy_matrix(graph)

        # self._beta = 0.0728

        # M_g = np.eye(graph.number_of_nodes()) - self._beta * A
        # M_l = self._beta * A

        M_g = np.eye(graph.number_of_nodes())
        M_l = np.dot(A, A)

        S = np.dot(np.linalg.inv(M_g), M_l)
        # s: \sigma_k
        u, s, vt = lg.svds(S, k=self._d // 2)
        sigma = np.diagflat(np.sqrt(s))
        X1 = np.dot(u, sigma)
        X2 = np.dot(vt.T, sigma)
        # self._X = X2
        self._X = np.concatenate((X1, X2), axis=1) 

def test_flip():
    X = np.diagflat([-2, -1, 0, 1, 2])
    # eigenvalues -2, -1, 0, 1, 2; eigenvectors are I

    Xflip = FlipPSD().fit_transform(X)
    assert np.allclose(Xflip, np.diagflat([2, 1, 0, 1, 2]))

    Xflip2 = FlipPSD().fit(X).transform(X)
    assert np.allclose(Xflip2, np.diagflat([2, 1, 0, 1, 2]))

    Xflip3 = FlipPSD().fit(X).transform(X[:3, :])
    assert np.allclose(Xflip3, [[2,0,0,0,0], [0,1,0,0,0], [0,0,0,0,0]])

    assert_raises(TypeError, lambda: FlipPSD().fit(X[:2, :]))
    assert_raises(TypeError, lambda: FlipPSD().fit_transform(X[:2, :]))
    assert_raises(TypeError, lambda: FlipPSD().fit(X).transform(X[:, :2])) 

def test_shift():
    X = np.diagflat([-2., -1, 0, 1, 2])
    # eigenvalues -2, -1, 0, 1, 2; eigenvectors are I

    Xshift = ShiftPSD().fit_transform(X)
    assert np.allclose(Xshift, np.diagflat([0, 1, 2, 3, 4]))

    Xshift2 = ShiftPSD().fit(X).transform(X)
    assert np.allclose(Xshift2, np.diagflat([0, 1, 2, 3, 4]))

    Xshift3 = ShiftPSD().fit(X).transform(X[:3, :])
    assert np.allclose(Xshift3, X[:3, :])

    Xshift4 = ShiftPSD(min_eig=2).fit_transform(X)
    assert np.allclose(Xshift4, np.diagflat([2, 3, 4, 5, 6]))

    assert_raises(TypeError, lambda: ShiftPSD().fit(X[:2, :]))
    assert_raises(TypeError, lambda: ShiftPSD().fit_transform(X[:2, :]))
    assert_raises(TypeError, lambda: ShiftPSD().fit(X).transform(X[:, :2])) 

def compute_shares_by_xi_jacobian(self, probabilities: Array, conditionals: Optional[Array]) -> Array:
        """Compute the Jacobian of shares with respect to xi (equivalently, to delta)."""
        diagonal_shares = np.diagflat(self.products.shares)
        weighted_probabilities = self.agents.weights * probabilities.T
        jacobian = diagonal_shares - probabilities @ weighted_probabilities

        if self.epsilon_scale != 1:
            jacobian /= self.epsilon_scale

        if self.H > 0:
            membership = self.get_membership_matrix()
            jacobian += self.rho / (1 - self.rho) * (
                diagonal_shares - membership * (conditionals @ weighted_probabilities)
            )

        return jacobian 

def get_matrix_elements():

    # In[5]: np.diagflat([1, 2, 3, 4])[:3, :]

    matrix = np.diagflat([1.0, 2.0, 3.0, 4.0])[:3, :]

    # Now matrix is:
    # array([[1, 0, 0, 0],
    #        [0, 2, 0, 0],
    #        [0, 0, 3, 0]])

    mc_energies = [1.0, 2.0, 3.0, 4.0, 5.0]

    ebounds = [1.0, 2.5, 4.5, 5.0]

    return matrix, mc_energies, ebounds 

def get_matrix_elements():

    # In[5]: np.diagflat([1, 2, 3, 4])[:3, :]

    matrix = np.diagflat([1.0, 2.0, 3.0, 4.0])[:3, :]

    # Now matrix is:
    # array([[1, 0, 0, 0],
    #        [0, 2, 0, 0],
    #        [0, 0, 3, 0]])

    mc_energies = [1.0, 2.0, 3.0, 4.0, 5.0]

    ebounds = [1.0, 2.5, 4.5, 5.0]

    return matrix, mc_energies, ebounds 

def learn_embedding(self):

        graph = self.g.G
        A = nx.to_numpy_matrix(graph)

        # self._beta = 0.0728

        # M_g = np.eye(graph.number_of_nodes()) - self._beta * A
        # M_l = self._beta * A

        M_g = np.eye(graph.number_of_nodes())
        M_l = np.dot(A, A)

        S = np.dot(np.linalg.inv(M_g), M_l)
        # s: \sigma_k
        u, s, vt = lg.svds(S, k=self._d // 2)
        sigma = np.diagflat(np.sqrt(s))
        X1 = np.dot(u, sigma)
        X2 = np.dot(vt.T, sigma)
        # self._X = X2
        self._X = np.concatenate((X1, X2), axis=1) 

def getLap(self):
        # degree_mat = np.diagflat(np.sum(self.adj_mat, axis=1))
        # print('np.diagflat(np.sum(self.adj_mat, axis=1))')
        # deg_trans = np.diagflat(np.reciprocal(np.sqrt(np.sum(self.adj_mat, axis=1))))
        # print('np.diagflat(np.reciprocal(np.sqrt(np.sum(self.adj_mat, axis=1))))')
        # deg_trans = np.nan_to_num(deg_trans)
        # L = degree_mat-self.adj_mat
        # print('begin norm_lap_mat')
        # # eye = np.eye(self.node_size)
        #
        # norm_lap_mat = np.matmul(np.matmul(deg_trans, L), deg_trans)
        G = self.g.G.to_undirected()
        print('begin norm_lap_mat')
        norm_lap_mat = nx.normalized_laplacian_matrix(G)
        print('finish norm_lap_mat')
        return norm_lap_mat 

def kernel_qchem_inter_rf(self, **kw):
    from pyscf.gw.gw import rpa_AB_matrices
    """ This is constructing A B matrices and diagonalizes the problem """

    nf = self.nfermi[0]
    nv = self.norbs-self.vstart[0]
    vs = self.vstart[0]
    
    x = self.mo_coeff[0,0,:,:,0]
    pab2v = self.pb.get_ac_vertex_array()
    self.pmn2v = pmn2v = einsum('nb,pmb->pmn', x[:nf,:], einsum('ma,pab->pmb', x[vs:,:], pab2v))
    pmn2c = einsum('qp,pmn->qmn', self.hkernel_den, pmn2v)
    meri = einsum('pmn,pik->mnik', pmn2c, pmn2v)
    #meri.fill(0.0)

    A = (diagflat( self.FmE ).reshape([nv,nf,nv,nf]) + meri).reshape([nv*nf,nv*nf])
    B = meri.reshape([nv*nf,nv*nf])

    assert np.allclose(A, A.transpose())
    assert np.allclose(B, B.transpose())

    AmB = A-B
    n = len(AmB)
    print(__name__)
    for i in range(n): print(i, AmB[i,i])
    
    ham_rpa = np.multiply(self.sqrt_FmE[:,None], np.multiply(A+B, self.sqrt_FmE))
    esq, self.s2z = np.linalg.eigh(ham_rpa)
    self.s2omega = np.sqrt(esq)
    print(self.s2omega*27.2114)

    self.s2z = self.s2z.T
    self.s2xpy = np.zeros_like(self.s2z)
    for s,(e2,z) in enumerate(zip(esq, self.s2z)):
      w = np.sqrt(e2)
      self.s2xpy[s] = np.multiply(self.sqrt_FmE, self.s2z[s])/np.sqrt(w)
      #print(e2, abs(np.dot(ham_rpa,z)-e2*z).sum())
    return self.s2omega,self.s2z 

def getLap(self):
        degree_mat = np.diagflat(np.sum(self.adj_mat, axis=1))
        deg_trans = np.diagflat(np.reciprocal(np.sqrt(np.sum(self.adj_mat, axis=1))))
        deg_trans = np.nan_to_num(deg_trans)
        L = degree_mat-self.adj_mat

        # eye = np.eye(self.node_size)
        norm_lap_mat = np.matmul(np.matmul(deg_trans, L), deg_trans)
        return norm_lap_mat 

def test_reduce_reorder(self):
        # gh 10370, 11029 Some compilers reorder the call to npy_getfloatstatus
        # and put it before the call to an intrisic function that causes
        # invalid status to be set. Also make sure warnings are not emitted
        for n in (2, 4, 8, 16, 32):
            for dt in (np.float32, np.float16, np.complex64):
                for r in np.diagflat(np.array([np.nan] * n, dtype=dt)):
                    assert_equal(np.min(r), np.nan) 

def testDiagflatExecution(self):
        a = diagflat([[1, 2], [3, 4]], chunk_size=1)

        res = self.executor.execute_tensor(a, concat=True)[0]
        expected = np.diagflat([[1, 2], [3, 4]])
        np.testing.assert_equal(res, expected)

        d = tensor([[1, 2], [3, 4]], chunk_size=1)
        a = diagflat(d)

        res = self.executor.execute_tensor(a, concat=True)[0]
        expected = np.diagflat([[1, 2], [3, 4]])
        np.testing.assert_equal(res, expected)

        a = diagflat([1, 2], 1, chunk_size=1)

        res = self.executor.execute_tensor(a, concat=True)[0]
        expected = np.diagflat([1, 2], 1)
        np.testing.assert_equal(res, expected)

        d = tensor([[1, 2]], chunk_size=1)
        a = diagflat(d, 1)

        res = self.executor.execute_tensor(a, concat=True)[0]
        expected = np.diagflat([1, 2], 1)
        np.testing.assert_equal(res, expected) 

def eye(n): return diagflat(ones(n)) 

def diagflat(a, k=0):
 if isinstance(a, garray): return a.diagflat(k)
 else: return numpy.diagflat(a,k) 

def diagflat(self, k=0):
  if self.ndim!=1: return self.ravel().diagflat(k)
  if k!=0: raise NotImplementedError('k!=0 for garray.diagflat')
  selfSize = self.size
  ret = zeros((selfSize, selfSize))
  ret.ravel()[:-1].reshape((selfSize-1, selfSize+1))[:, 0] = self[:-1]
  if selfSize!=0: ret.ravel()[-1] = self[-1]
  return ret 

def diagonal(self):
  if self.ndim==1: return self.diagflat()
  if self.ndim==2:
   if self.shape[0] > self.shape[1]: return self[:self.shape[1]].diagonal()
   if self.shape[1] > self.shape[0]: return self[:, :self.shape[0]].diagonal()
   return self.ravel()[::self.shape[0]+1]
  raise NotImplementedError('garray.diagonal for arrays with ndim other than 1 or 2.') 

def test_reduce_warns(self):
        # gh 10370, 11029 Some compilers reorder the call to npy_getfloatstatus
        # and put it before the call to an intrisic function that causes
        # invalid status to be set. Also make sure warnings are emitted
        for n in (2, 4, 8, 16, 32):
            with suppress_warnings() as sup:
                sup.record(RuntimeWarning)
                for r in np.diagflat([np.nan] * n):
                    assert_equal(np.min(r), np.nan)
                assert_equal(len(sup.log), n) 

def gen_roots_and_weights(n, an_func, sqrt_bn_func, mu):
    """[x,w] = gen_roots_and_weights(n,an_func,sqrt_bn_func,mu)

    Returns the roots (x) of an nth order orthogonal polynomial,
    and weights (w) to use in appropriate Gaussian quadrature with that
    orthogonal polynomial.

    The polynomials have the recurrence relation
          P_n+1(x) = (x - A_n) P_n(x) - B_n P_n-1(x)

    an_func(n)          should return A_n
    sqrt_bn_func(n)     should return sqrt(B_n)
    mu ( = h_0 )        is the integral of the weight over the orthogonal interval
    """
    nn = np.arange(1.0,n)
    sqrt_bn = sqrt_bn_func(nn)
    an = an_func(np.concatenate(([0], nn)))
    x, v = eig((np.diagflat(an) +
                np.diagflat(sqrt_bn,1) +
                np.diagflat(sqrt_bn,-1)))
    answer = []
    sortind = x.real.argsort()
    answer.append(x[sortind])
    answer.append((mu*v[0]**2)[sortind])
    return answer

# Jacobi Polynomials 1               P^(alpha,beta)_n(x) 

def _do_rescaled_decomposition(self, S):
        """
        Decomposing the similarity matrix.
        """
        U, sigmas, Vt = sps.linalg.svds(S, k=int(self.dimensions/2))
        sigmas = np.diagflat(np.sqrt(sigmas))
        self._left_embedding = np.dot(U, sigmas)
        self._right_embedding = np.dot(Vt.T, sigmas) 

def test_reduce_reorder(self):
        # gh 10370, 11029 Some compilers reorder the call to npy_getfloatstatus
        # and put it before the call to an intrisic function that causes
        # invalid status to be set. Also make sure warnings are not emitted
        for n in (2, 4, 8, 16, 32):
            for dt in (np.float32, np.float16, np.complex64):
                for r in np.diagflat(np.array([np.nan] * n, dtype=dt)):
                    assert_equal(np.min(r), np.nan) 

def testDiagFlat(self):
    array_transforms = [
        lambda x: x,  # Identity,
        tf.convert_to_tensor,
        np.array,
        lambda x: np.array(x, dtype=np.float32),
        lambda x: np.array(x, dtype=np.float64),
        array_ops.array,
        lambda x: array_ops.array(x, dtype=np.float32),
        lambda x: array_ops.array(x, dtype=np.float64)
    ]

    def run_test(arr):
      for fn in array_transforms:
        arr = fn(arr)
        self.match(
            array_ops.diagflat(arr),
            np.diagflat(arr),
            msg='diagflat({})'.format(arr))
        for k in range(-3, 3):
          self.match(
              array_ops.diagflat(arr, k),
              np.diagflat(arr, k),
              msg='diagflat({}, k={})'.format(arr, k))

    # 1-d arrays.
    run_test([])
    run_test([1])
    run_test([1, 2])
    # 2-d arrays.
    run_test([[]])
    run_test([[5]])
    run_test([[], []])
    run_test(np.arange(4).reshape((2, 2)).tolist())
    run_test(np.arange(2).reshape((2, 1)).tolist())
    run_test(np.arange(2).reshape((1, 2)).tolist())
    # 3-d arrays
    run_test(np.arange(8).reshape((2, 2, 2)).tolist()) 

def annihilation(dimension):
    """
    Returns a dense matrix of size dimension x dimension representing the annihilation operator in number basis.

    Parameters
    ----------
    dimension: int

    Returns
    -------
    ndarray
        annihilation operator matrix, size dimension x dimension
    """
    offdiag_elements = np.sqrt(range(1, dimension))
    return np.diagflat(offdiag_elements, 1) 

def test_reduce_warns(self):
        # gh 10370, 11029 Some compilers reorder the call to npy_getfloatstatus
        # and put it before the call to an intrisic function that causes
        # invalid status to be set. Also make sure warnings are emitted
        for n in (2, 4, 8, 16, 32):
            with suppress_warnings() as sup:
                sup.record(RuntimeWarning)
                for r in np.diagflat([np.nan] * n):
                    assert_equal(np.min(r), np.nan)
                assert_equal(len(sup.log), n)