Python scipy.sparse.identity() Examples

def _compute(self, data_1, data_2):
        data_1 = basic.graphs_to_adjacency_lists(data_1)
        data_2 = basic.graphs_to_adjacency_lists(data_2)
        res = np.zeros((len(data_1), len(data_2)))
        N = len(data_1) * len(data_2)
        for i, graph1 in enumerate(data_1):
            for j, graph2 in enumerate(data_2):
                # norm1, norm2 - normalized adjacency matrixes
                norm1 = _norm(graph1)
                norm2 = _norm(graph2)
                # if graph is unweighted, W_prod = kron(a_norm(g1)*a_norm(g2))
                w_prod = kron(lil_matrix(norm1), lil_matrix(norm2))
                starting_prob = np.ones(w_prod.shape[0]) / (w_prod.shape[0])
                stop_prob = starting_prob
                # first solve (I - lambda * W_prod) * x = p_prod
                A = identity(w_prod.shape[0]) - (w_prod * self._lmb)
                x = lsqr(A, starting_prob)
                res[i, j] = stop_prob.T.dot(x[0])
                # print float(len(data_2)*i + j)/float(N), "%"
        return res 

def suggest(self, method=s.RANDOM, eps=1e-8, *args, **kwargs):
        if method not in s.suggest_methods:
            raise Exception("Unknown suggest method: %s\n", method)
        if method == s.RANDOM:
            x = np.random.randn(self.n)
        elif method == s.SPECTRAL:
            if self.spectral_sol is None:
                self.spectral_sol, self.spectral_bound = solve_spectral(self.qcqp_form, *args, **kwargs)
                if self.maximize_flag:
                    self.spectral_bound *= -1
            x = self.spectral_sol
        elif method == s.SDR:
            if self.sdr_sol is None:
                self.sdr_sol, self.sdr_bound = solve_sdr(self.qcqp_form, *args, **kwargs)
                if self.maximize_flag:
                    self.sdr_bound *= -1
                self.mu = np.asarray(self.sdr_sol[:-1, -1]).flatten()
                self.Sigma = self.sdr_sol[:-1, :-1] - self.mu*self.mu.T + eps*sp.identity(self.n)
            x = np.random.multivariate_normal(self.mu, self.Sigma)

        assign_vars(self.prob.variables(), x)
        f0 = self.qcqp_form.f0.eval(x)
        if self.maximize_flag: f0 *= -1
        return (f0, max(self.qcqp_form.violations(x))) 

def get_measurements(model, domain_shape):
        # model is a set of contingency tables to calculate
        # each contingency table is a list of [(attribute, size)] 
        M = []
        for table in model:
            Q = [np.ones((1,size)) for size in domain_shape]
            for attribute, size in table:
                full_size = domain_shape[attribute]
                I = sparse.identity(size) 
                if size != full_size:
                    P = PrivBayesSelect.domain_transform(size, full_size)
                    Q[attribute] = I * P
                elif size == full_size:
                    Q[attribute] = I
                else:
                    print('bug here')
            M.append(reduce(sparse.kron, Q))
        return sparse.vstack(M) 

def lu_factor(sigma, A):
    """
    LU Factorisation wrapper of:

    .. math:: LU = (\sigma \mathbf{I} - \mathbf{A})

    In the case of ``A`` being a sparse matrix, the sparse methods in scipy are employed

    Args:
        sigma (float): Expansion frequency
        A (csc_matrix or np.ndarray): Dynamics matrix

    Returns:
        tuple or SuperLU: tuple (dense) or SuperLU (sparse) objects containing the LU factorisation
    """
    n = A.shape[0]
    if type(A) == libsp.csc_matrix:
        return scsp.linalg.splu(sigma * scsp.identity(n, dtype=complex, format='csc') - A)
    else:
        return sclalg.lu_factor(sigma * np.eye(n) - A) 

def __init__(self, graph, **kwargs):

        if graph.is_weighted():
            logger.warning('Your graph is weighted, and is considered '
                           'unweighted to build a binary line graph.')

        graph.compute_differential_operator()
        # incidence = np.abs(graph.D)  # weighted?
        incidence = (graph.D != 0)

        adjacency = incidence.T.dot(incidence).astype(np.int)
        adjacency -= sparse.identity(graph.n_edges, dtype=np.int)

        try:
            coords = incidence.T.dot(graph.coords) / 2
        except AttributeError:
            coords = None

        super(LineGraph, self).__init__(adjacency, coords=coords,
                plotting=graph.plotting, **kwargs) 

def sites2coupling_sparse(mij,spins):
  """Create the couplings of the HP-hamiltonian"""
  # SzSz terms
  # There are two matrices, one that renormalizes onsite energies
  # and another one that creates hoppings
  n = len(spins)
  ons = coo_matrix(([],([],[])),shape=(n,n),dtype=np.complex)
  hop = coo_matrix(([],([],[])),shape=(n,n),dtype=np.complex)
  iden = sparseiden(n,dtype=np.complex) # identity matrix
  raise # not finished
  
  for i in range(n): # loop over spins
    for j in range(n): # loop over spins
      cij = genij(i,j) # matrix with the couplings SiSj
      # the quantization axis is z
      ons[i,i] += cij[2,2] # Sz Sz component
      # there are missing terms!!!!!!
  # S+ S- terms
      sij = np.sqrt(spins[i]*spins[j]) # denominator
      hop[i,j] = -(cij[0,0] + cij[1,1])*sij/2.
  return (np.matrix(ons),np.matrix(hop)) # return matrices 

def test_sparse_lindblad_param(self):
        #Test sparse Lindblad gates

        print("\nGate Test:")
        SparseId = sps.identity(4**2,'d','csr')
        gate = LindbladDenseOp.from_operation_matrix( np.identity(4**2,'d') )
        print("gate Errgen type (should be dense):",type(gate.errorgen.todense()))
        self.assertIsInstance(gate.errorgen.todense(), np.ndarray)
        sparseOp = LindbladOp.from_operation_matrix( SparseId )
        print("spareGate Errgen type (should be sparse):",type(sparseOp.errorgen.tosparse()))
        self.assertIsInstance(sparseOp.errorgen.tosparse(), sps.csr_matrix)
        self.assertArraysAlmostEqual(gate.errorgen.todense(),sparseOp.errorgen.todense())

        perfectG = std2Q_XYICNOT.target_model().operations['Gix'].copy()
        noisyG = std2Q_XYICNOT.target_model().operations['Gix'].copy()
        noisyG.depolarize(0.9)
        Sparse_noisyG = sps.csr_matrix(noisyG,dtype='d')
        Sparse_perfectG = sps.csr_matrix(perfectG,dtype='d')
        op2 = LindbladDenseOp.from_operation_matrix( noisyG, perfectG )
        sparseGate2 = LindbladOp.from_operation_matrix( Sparse_noisyG, Sparse_perfectG )
        print("spareGate2 Errgen type (should be sparse):",type(sparseGate2.errorgen.tosparse()))
        self.assertIsInstance(sparseGate2.errorgen.tosparse(), sps.csr_matrix)
        #print("errgen = \n"); pygsti.tools.print_mx(op2.err_gen,width=4,prec=1)
        #print("sparse errgen = \n"); pygsti.tools.print_mx(sparseGate2.err_gen.toarray(),width=4,prec=1)
        self.assertArraysAlmostEqual(op2.errorgen.todense(),sparseGate2.errorgen.todense()) 

def get_interface(h,fun=None):
  """Return an operator that projects onte the interface"""
  dind = 1 # index to which divide the positions
  if h.has_spin:  dind *= 2 # duplicate for spin
  if h.has_eh:  dind *= 2  # duplicate for eh
  iden = csc(np.matrix(np.identity(dind,dtype=np.complex))) # identity matrix
  r = h.geometry.r # positions
  out = [[None for ri in r] for rj in r] # initialize
  if fun is None: # no input function
    cut = 2.0 # cutoff
    if h.dimensionality==1: index = 1
    elif h.dimensionality==2: index = 2
    else: raise
    def fun(ri): # define the function
      if np.abs(ri[index])<cut: return 1.0
      else: return 0.0
  for i in range(len(r)): # loop over positions
    out[i][i] = fun(r[i])*iden 
  return bmat(out) # return matrix 

def to_unitary(scaled_unitary):
    """
    Compute the scaling factor required to turn a scalar multiple of a unitary matrix
    to a unitary matrix.

    Parameters
    ----------
    scaled_unitary : ndarray
        A scaled unitary matrix

    Returns
    -------
    scale : float
    unitary : ndarray
        Such that `scale * unitary == scaled_unitary`.

    """
    scaled_identity = _np.dot(scaled_unitary, _np.conjugate(scaled_unitary.T))
    scale = _np.sqrt(scaled_identity[0, 0])
    assert(_np.allclose(scaled_identity / (scale**2), _np.identity(scaled_identity.shape[0], 'd'))), \
        "Given `scaled_unitary` does not appear to be a scaled unitary matrix!"
    return scale, (scaled_unitary / scale) 

def add_swave(delta=0.0,is_sparse=False,rs=None):
  """ Adds swave pairing """
  if rs is None: raise # raise error to signal that this is temporal
  n = len(rs) # number of sites
  if callable(delta): # delta is a function
    datar = [delta(ri) for ri in rs] # generate data for the different positions
    data = []
    # the coupling involves up and down
    for dr in datar: # loop over positions
      data.append(dr)  # up e dn h
      data.append(dr)  # dne up h
    iis = range(n*2) # indexes
    coupling = csc_matrix((data,(iis,iis)),dtype=np.complex)  # generate matrix
  else:
    coupling = sp.identity(n*2)*delta # delta matrix
  zero = coupling*0.
  return build_eh(zero,coupling=coupling,is_sparse=is_sparse) # return matrix 

def _fit(self):

        # Versions using sparse matrices
        # adj = nx.adjacency_matrix(self._G)
        # ident = sparse.identity(len(self._G.nodes)).tocsc()
        # sim = inv(ident - adj.multiply(self.beta).T) - ident
        # adj = nx.adjacency_matrix(self._G)
        # aux = adj.multiply(-self.beta).T
        # aux.setdiag(1+aux.diagonal(), k=0)
        # sim = inv(aux)
        # sim.setdiag(sim.diagonal()-1)
        # print(sim.nnz)
        # print(adj.nnz)

        # Version using dense matrices
        adj = nx.adjacency_matrix(self._G)
        aux = adj.T.multiply(-self.beta).todense()
        np.fill_diagonal(aux, 1+aux.diagonal())
        sim = np.linalg.inv(aux)
        np.fill_diagonal(sim, sim.diagonal()-1)
        return sim 

def laplacian_matrix(self, weighted=False, transposed=False):
        """
        Returns the transposed normalized Laplacian matrix corresponding to the network.

        Parameters
        ----------

        Returns
        -------

        """
        if weighted:
            A = self.transition_matrix().transpose()
            D = _sparse.identity(self.ncount())
        else:
            A = self.adjacency_matrix(weighted=False)
            D = _sparse.diags(_np.array([float(self.nodes[v]['degree']) for v in self.nodes]))
        L = D - A
        if transposed:
            return L.transpose()
        return L 

def laplacian_matrix(self, include_subpaths=True):
        """
        Returns the transposed Laplacian matrix corresponding to the higher-order network.

        Parameters
        ----------
        include_subpaths: bool
            Whether or not subpath statistics shall be included in the calculation of
            matrix weights

        Returns
        -------

        """
        transition_matrix = self.transition_matrix(include_subpaths)
        identity_matrix = _sparse.identity(self.ncount())

        return identity_matrix - transition_matrix 

def partial_trace_super(dim1: int, dim2: int) -> np.array:
    """
    Return the partial trace superoperator in the column-major basis.

    This returns the superoperator S_TrB such that:
        S_TrB * vec(rho_AB) = vec(rho_A)
    for rho_AB = kron(rho_A, rho_B)

    Args:
        dim1: the dimension of the system not being traced
        dim2: the dimension of the system being traced over

    Returns:
        A Numpy array of the partial trace superoperator S_TrB.
    """

    iden = sps.identity(dim1)
    ptr = sps.csr_matrix((dim1 * dim1, dim1 * dim2 * dim1 * dim2))

    for j in range(dim2):
        v_j = sps.coo_matrix(([1], ([0], [j])), shape=(1, dim2))
        tmp = sps.kron(iden, v_j.tocsr())
        ptr += sps.kron(tmp, tmp)

    return ptr 

def _zeropi_operator_in_product_basis(self, zeropi_operator, zeropi_evecs=None):
        """Helper method that converts a zeropi operator into one in the product basis.

        Returns
        -------
        scipy.sparse.csc_matrix
            operator written in the product basis
        """
        zeropi_dim = self.zeropi_cutoff
        zeta_dim = self.zeta_cutoff

        if zeropi_evecs is None:
            _, zeropi_evecs = self._zeropi.eigensys(evals_count=zeropi_dim)

        op_eigen_basis = sparse.dia_matrix((zeropi_dim, zeropi_dim),
                                           dtype=np.complex_)  # is this guaranteed to be zero?

        op_zeropi = spec_utils.get_matrixelement_table(zeropi_operator, zeropi_evecs)
        for n in range(zeropi_dim):
            for m in range(zeropi_dim):
                op_eigen_basis += op_zeropi[n, m] * op.hubbard_sparse(n, m, zeropi_dim)

        return sparse.kron(op_eigen_basis, sparse.identity(zeta_dim, format='csc', dtype=np.complex_), format='csc') 

def sparse_kinetic_mat(self):
        """
        Kinetic energy portion of the Hamiltonian.
        TODO: update this method to use single-variable operator methods

        Returns
        -------
        scipy.sparse.csc_matrix
            matrix representing the kinetic energy operator
        """
        pt_count = self.grid.pt_count
        dim_theta = 2 * self.ncut + 1
        identity_phi = sparse.identity(pt_count, format='csc', dtype=np.complex_)
        identity_theta = sparse.identity(dim_theta, format='csc', dtype=np.complex_)

        kinetic_matrix_phi = self.grid.second_derivative_matrix(prefactor=-2.0 * self.ECJ)

        diag_elements = 2.0 * self.ECS * np.square(np.arange(-self.ncut + self.ng, self.ncut + 1 + self.ng))
        kinetic_matrix_theta = sparse.dia_matrix((diag_elements, [0]), shape=(dim_theta, dim_theta)).tocsc()

        kinetic_matrix = (sparse.kron(kinetic_matrix_phi, identity_theta, format='csc')
                          + sparse.kron(identity_phi, kinetic_matrix_theta, format='csc'))

        kinetic_matrix -= 2.0 * self.ECS * self.dCJ * self.i_d_dphi_operator() * self.n_theta_operator()
        return kinetic_matrix 

def A_to_diffusion_kernel(A, k):
    """
    Computes [A**0, A**1, ..., A**k]

    :param A: 2d numpy array
    :param k: integer, degree of series
    :return: 3d numpy array [A**0, A**1, ..., A**k]
    """
    assert k >= 0

    Apow = [np.identity(A.shape[0])]

    if k > 0:
        d = A.sum(0)

        Apow.append(A / (d + 1.0))

        for i in range(2, k + 1):
            Apow.append(np.dot(A / (d + 1.0), Apow[-1]))

    return np.transpose(np.asarray(Apow, dtype='float32'), (1, 0, 2)) 

def sparse_A_to_diffusion_kernel(A, k):
    assert k >= 0

    num_nodes = A.shape[0]

    Apow = [sp.identity(num_nodes)]

    if k > 0:
        d = A.sum(0)

        Apow.append(A / (d + 1.0))

        for i in range(2, k + 1):
            Apow.append((A / (d + 1.0)).dot(Apow[-1]))

    return Apow 

def _compute_deepwalk_matrix(self, A, window, b):
        # directly compute deepwalk matrix
        n = A.shape[0]
        vol = float(A.sum())
        L, d_rt = sp.csgraph.laplacian(A, normed=True, return_diag=True)
        # X = D^{-1/2} A D^{-1/2}
        X = sp.identity(n) - L
        S = np.zeros_like(X)
        X_power = sp.identity(n)
        for i in range(window):
            print("Compute matrix %d-th power", i + 1)
            X_power = X_power.dot(X)
            S += X_power
        S *= vol / window / b
        D_rt_inv = sp.diags(d_rt ** -1)
        M = D_rt_inv.dot(D_rt_inv.dot(S).T).todense()
        M[M <= 1] = 1
        Y = np.log(M)
        return sp.csr_matrix(Y) 

def test_invPropertyTensor3D(self):
        M = discretize.TensorMesh([6, 6, 6])
        a1 = np.random.rand(M.nC)
        a2 = np.random.rand(M.nC)
        a3 = np.random.rand(M.nC)
        a4 = np.random.rand(M.nC)
        a5 = np.random.rand(M.nC)
        a6 = np.random.rand(M.nC)
        prop1 = a1
        prop2 = np.c_[a1, a2, a3]
        prop3 = np.c_[a1, a2, a3, a4, a5, a6]

        for prop in [4, prop1, prop2, prop3]:
            b = invPropertyTensor(M, prop)
            A = makePropertyTensor(M, prop)
            B1 = makePropertyTensor(M, b)
            B2 = invPropertyTensor(M, prop, returnMatrix=True)

            Z = B1*A - sp.identity(M.nC*3)
            self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
            Z = B2*A - sp.identity(M.nC*3)
            self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL) 

def test_invPropertyTensor2D(self):
        M = discretize.TensorMesh([6, 6])
        a1 = np.random.rand(M.nC)
        a2 = np.random.rand(M.nC)
        a3 = np.random.rand(M.nC)
        prop1 = a1
        prop2 = np.c_[a1, a2]
        prop3 = np.c_[a1, a2, a3]

        for prop in [4, prop1, prop2, prop3]:
            b = invPropertyTensor(M, prop)
            A = makePropertyTensor(M, prop)
            B1 = makePropertyTensor(M, b)
            B2 = invPropertyTensor(M, prop, returnMatrix=True)

            Z = B1*A - sp.identity(M.nC*2)
            self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
            Z = B2*A - sp.identity(M.nC*2)
            self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL) 

def test_invXXXBlockDiagonal(self):
        a = [np.random.rand(5, 1) for i in range(4)]

        B = inv2X2BlockDiagonal(*a)

        A = sp.vstack((sp.hstack((sdiag(a[0]), sdiag(a[1]))),
                       sp.hstack((sdiag(a[2]), sdiag(a[3])))))

        Z2 = B*A - sp.identity(10)
        self.assertTrue(np.linalg.norm(Z2.todense().ravel(), 2) < TOL)

        a = [np.random.rand(5, 1) for i in range(9)]
        B = inv3X3BlockDiagonal(*a)

        A = sp.vstack((sp.hstack((sdiag(a[0]), sdiag(a[1]),  sdiag(a[2]))),
                       sp.hstack((sdiag(a[3]), sdiag(a[4]),  sdiag(a[5]))),
                       sp.hstack((sdiag(a[6]), sdiag(a[7]),  sdiag(a[8])))))

        Z3 = B*A - sp.identity(15)

        self.assertTrue(np.linalg.norm(Z3.todense().ravel(), 2) < TOL) 

def get_movielens_100k(min_positive_score=4, negative_value=0):
    movielens_100k_dict = datasets.fetch_movielens(indicator_features=True, genre_features=True)

    def flip_ratings(ratings_matrix):
        ratings_matrix.data = np.array([1 if rating >= min_positive_score else negative_value
                                        for rating in ratings_matrix.data])
        return ratings_matrix

    test_interactions = flip_ratings(movielens_100k_dict['test'])
    train_interactions = flip_ratings(movielens_100k_dict['train'])

    # Create indicator features for all users
    num_users = train_interactions.shape[0]
    user_features = sp.identity(num_users)

    # Movie titles
    titles = movielens_100k_dict['item_labels']

    return train_interactions, test_interactions, user_features, movielens_100k_dict['item_features'], titles 

def test_overflow_predict():

    no_users, no_items = (1000, 1000)

    train = sp.rand(no_users, no_items, format="csr", random_state=42)

    model = LightFM(loss="warp")

    model.fit(train)

    with pytest.raises((ValueError, OverflowError)):
        print(
            model.predict(
                1231241241231241414,
                np.arange(no_items),
                user_features=sp.identity(no_users),
            )
        ) 

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 load_adj(pkl_filename, adjtype):
    sensor_ids, sensor_id_to_ind, adj_mx = load_pickle(pkl_filename)
    if adjtype == "scalap":
        adj = [calculate_scaled_laplacian(adj_mx)]
    elif adjtype == "normlap":
        adj = [calculate_normalized_laplacian(adj_mx).astype(np.float32).todense()]
    elif adjtype == "symnadj":
        adj = [sym_adj(adj_mx)]
    elif adjtype == "transition":
        adj = [asym_adj(adj_mx)]
    elif adjtype == "doubletransition":
        adj = [asym_adj(adj_mx), asym_adj(np.transpose(adj_mx))]
    elif adjtype == "identity":
        adj = [np.diag(np.ones(adj_mx.shape[0])).astype(np.float32)]
    else:
        error = 0
        assert error, "adj type not defined"
    return sensor_ids, sensor_id_to_ind, adj 

def direct_compute_deepwalk_matrix(A, window, b):
    n = A.shape[0]
    vol = float(A.sum())
    L, d_rt = csgraph.laplacian(A, normed=True, return_diag=True)
    # X = D^{-1/2} A D^{-1/2}
    X = sparse.identity(n) - L
    S = np.zeros_like(X)
    X_power = sparse.identity(n)
    for i in range(window):
        logger.info("Compute matrix %d-th power", i+1)
        X_power = X_power.dot(X)
        S += X_power
    S *= vol / window / b
    D_rt_inv = sparse.diags(d_rt ** -1)
    M = D_rt_inv.dot(D_rt_inv.dot(S).T)
    m = T.matrix()
    f = theano.function([m], T.log(T.maximum(m, 1)))
    Y = f(M.todense().astype(theano.config.floatX))
    return sparse.csr_matrix(Y) 

def permuteE(self):
        """Permutation matrix re-ordering of edges sorted by x, then y, then z"""
        # TODO: cache these?
        Px = np.lexsort(self.gridEx.T)
        Py = np.lexsort(self.gridEy.T) + self.nEx
        if self.dim == 2:
            P = np.r_[Px, Py]
        if self.dim == 3:
            Pz = np.lexsort(self.gridEz.T) + (self.nEx+self.nEy)
            P = np.r_[Px, Py, Pz]
        return sp.identity(self.nE).tocsr()[P] 

def get_user_representations(self, features=None):
        """
        Get the latent representations for users given model and features.

        Arguments
        ---------

        features: np.float32 csr_matrix of shape [n_users, n_user_features], optional
             Each row contains that user's weights over features.
             An identity matrix will be used if not supplied.

        Returns
        -------

        (user_biases, user_embeddings):
                (np.float32 array of shape n_users
                 np.float32 array of shape [n_users, num_components]
            Biases and latent representations for users.
        """

        self._check_initialized()

        if features is None:
            return self.user_biases, self.user_embeddings

        features = sp.csr_matrix(features, dtype=CYTHON_DTYPE)

        return features * self.user_biases, features * self.user_embeddings 

def permuteCC(self):
        """Permutation matrix re-ordering of cells sorted by x, then y, then z"""
        # TODO: cache these?
        P = np.lexsort(self.gridCC.T) # sort by x, then y, then z
        return sp.identity(self.nC).tocsr()[P]