Python scipy.sparse.kron() Examples

def get2DMatrix(N, h, bc_hor, bc_ver, order):
    assert np.size(N) == 2, 'N needs to be an array with two entries: N[0]=Nx and N[1]=Nz'
    assert np.size(h) == 2, 'h needs to be an array with two entries: h[0]=dx and h[1]=dz'

    Ax = getMatrix(N[0], h[0], bc_hor[0], bc_hor[1], order)
    Az = getMatrix(N[1], h[1], bc_ver[0], bc_ver[1], order)

    Dx = sp.kron(Ax, sp.eye(N[1]), format="csr")
    Dz = sp.kron(sp.eye(N[0]), Az, format="csr")

    return Dx, Dz


#
# NOTE: So far only constant dirichlet values can be prescribed, i.e. one fixed value for a whole segment
# 

def test_diagonal(self):
        # Does the matrix's .diagonal() method work?
        mats = []
        mats.append([[1,0,2]])
        mats.append([[1],[0],[2]])
        mats.append([[0,1],[0,2],[0,3]])
        mats.append([[0,0,1],[0,0,2],[0,3,0]])

        mats.append(kron(mats[0],[[1,2]]))
        mats.append(kron(mats[0],[[1],[2]]))
        mats.append(kron(mats[1],[[1,2],[3,4]]))
        mats.append(kron(mats[2],[[1,2],[3,4]]))
        mats.append(kron(mats[3],[[1,2],[3,4]]))
        mats.append(kron(mats[3],[[1,2,3,4]]))

        for m in mats:
            assert_equal(self.spmatrix(m).diagonal(),diag(m)) 

def test_sparse_format_conversions(self):
        A = sparse.kron([[1,0,2],[0,3,4],[5,0,0]], [[1,2],[0,3]])
        D = A.todense()
        A = self.spmatrix(A)

        for format in ['bsr','coo','csc','csr','dia','dok','lil']:
            a = A.asformat(format)
            assert_equal(a.format,format)
            assert_array_equal(a.todense(), D)

            b = self.spmatrix(D+3j).asformat(format)
            assert_equal(b.format,format)
            assert_array_equal(b.todense(), D+3j)

            c = eval(format + '_matrix')(A)
            assert_equal(c.format,format)
            assert_array_equal(c.todense(), D) 

def test_constructor1(self):
        # check native BSR format constructor
        indptr = array([0,2,2,4])
        indices = array([0,2,2,3])
        data = zeros((4,2,3))

        data[0] = array([[0, 1, 2],
                         [3, 0, 5]])
        data[1] = array([[0, 2, 4],
                         [6, 0, 10]])
        data[2] = array([[0, 4, 8],
                         [12, 0, 20]])
        data[3] = array([[0, 5, 10],
                         [15, 0, 25]])

        A = kron([[1,0,2,0],[0,0,0,0],[0,0,4,5]], [[0,1,2],[3,0,5]])
        Asp = bsr_matrix((data,indices,indptr),shape=(6,12))
        assert_equal(Asp.todense(),A)

        # infer shape from arrays
        Asp = bsr_matrix((data,indices,indptr))
        assert_equal(Asp.todense(),A) 

def whiten(self, X):
        """
        SUR whiten method.

        Parameters
        -----------
        X : list of arrays
            Data to be whitened.

        Returns
        -------
        If X is the exogenous RHS of the system.
        ``np.dot(np.kron(cholsigmainv,np.eye(M)),np.diag(X))``

        If X is the endogenous LHS of the system.

        """
        nobs = self.nobs
        if X is self.endog: # definitely not a robust check
            return np.dot(np.kron(self.cholsigmainv,np.eye(nobs)),
                X.reshape(-1,1))
        elif X is self.sp_exog:
            return (sparse.kron(self.cholsigmainv,
                sparse.eye(nobs,nobs))*X).toarray()#*=dot until cast to array 

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 _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 sparse_d_potential_d_flux_mat(self):
        r"""Calculates a of the potential energy w.r.t flux, at the current value of flux,
        as stored in the object.

        The flux is assumed to be given in the units of the ratio \Phi_{ext}/\Phi_0.
        So if \frac{\partial U}{ \partial \Phi_{\rm ext}}, is needed, the expression returned
        by this function, needs to be multiplied by 1/\Phi_0.

        Returns
        -------
        scipy.sparse.csc_matrix
            matrix representing the derivative of the potential energy
        """
        op_1 = sparse.kron(self._sin_phi_operator(x=- 2.0 * np.pi * self.flux / 2.0),
                           self._cos_theta_operator(), format='csc')
        op_2 = sparse.kron(self._cos_phi_operator(x=- 2.0 * np.pi * self.flux / 2.0),
                           self._sin_theta_operator(), format='csc')
        return - 2.0 * np.pi * self.EJ * op_1 - np.pi * self.EJ * self.dEJ * op_2 

def aveCC2F(self):
        "Construct the averaging operator on cell centers to faces."
        if getattr(self, '_aveCC2F', None) is None:
            if self.dim == 1:
                self._aveCC2F = av_extrap(self.nCx)
            elif self.dim == 2:
                self._aveCC2F = sp.vstack((
                    sp.kron(speye(self.nCy), av_extrap(self.nCx)),
                    sp.kron(av_extrap(self.nCy), speye(self.nCx))
                ), format="csr")
            elif self.dim == 3:
                self._aveCC2F = sp.vstack((
                    kron3(
                        speye(self.nCz), speye(self.nCy), av_extrap(self.nCx)
                    ),
                    kron3(
                        speye(self.nCz), av_extrap(self.nCy), speye(self.nCx)
                    ),
                    kron3(
                        av_extrap(self.nCz), speye(self.nCy), speye(self.nCx)
                    )
                ), format="csr")
        return self._aveCC2F 

def test_diagonal(self):
        # Does the matrix's .diagonal() method work?
        mats = []
        mats.append([[1,0,2]])
        mats.append([[1],[0],[2]])
        mats.append([[0,1],[0,2],[0,3]])
        mats.append([[0,0,1],[0,0,2],[0,3,0]])

        mats.append(kron(mats[0],[[1,2]]))
        mats.append(kron(mats[0],[[1],[2]]))
        mats.append(kron(mats[1],[[1,2],[3,4]]))
        mats.append(kron(mats[2],[[1,2],[3,4]]))
        mats.append(kron(mats[3],[[1,2],[3,4]]))
        mats.append(kron(mats[3],[[1,2,3,4]]))

        for m in mats:
            rows, cols = array(m).shape
            sparse_mat = self.spmatrix(m)
            for k in range(-rows + 1, cols):
                assert_equal(sparse_mat.diagonal(k=k), diag(m, k=k))
            assert_raises(ValueError, sparse_mat.diagonal, -rows)
            assert_raises(ValueError, sparse_mat.diagonal, cols)

        # Test all-zero matrix.
        assert_equal(self.spmatrix((40, 16130)).diagonal(), np.zeros(40)) 

def test_constructor1(self):
        # check native BSR format constructor
        indptr = array([0,2,2,4])
        indices = array([0,2,2,3])
        data = zeros((4,2,3))

        data[0] = array([[0, 1, 2],
                         [3, 0, 5]])
        data[1] = array([[0, 2, 4],
                         [6, 0, 10]])
        data[2] = array([[0, 4, 8],
                         [12, 0, 20]])
        data[3] = array([[0, 5, 10],
                         [15, 0, 25]])

        A = kron([[1,0,2,0],[0,0,0,0],[0,0,4,5]], [[0,1,2],[3,0,5]])
        Asp = bsr_matrix((data,indices,indptr),shape=(6,12))
        assert_equal(Asp.todense(),A)

        # infer shape from arrays
        Asp = bsr_matrix((data,indices,indptr))
        assert_equal(Asp.todense(),A) 

def _areaFzFull(self):
        """
        area of z-faces prior to deflation
        """
        if self.isSymmetric:
            return np.kron(
                    np.ones_like(self.vectorNz), pi*(
                        self.vectorNx**2 -
                        np.r_[0, self.vectorNx[:-1]]**2
                    )
            )
        return np.kron(
            np.ones(self._ntNz), np.kron(
                self.hy,
                0.5 * (self.vectorNx[1:]**2 - self.vectorNx[:-1]**2)
            )
        ) 

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 h_gridded(self):
        """
        Returns an (nC, dim) numpy array with the widths of all cells in order
        """

        if self.dim == 1:
            return np.reshape(self.h, (self.nC, 1))
        elif self.dim == 2:
            hx = np.kron(np.ones(self.nCy), self.h[0])
            hy = np.kron(self.h[1], np.ones(self.nCx))
            return np.c_[hx, hy]
        elif self.dim == 3:
            hx = np.kron(np.ones(self.nCy*self.nCz), self.h[0])
            hy = np.kron(np.ones(self.nCz), np.kron(self.h[1], np.ones(self.nCx)))
            hz = np.kron(self.h[2], np.ones(self.nCx*self.nCy))
            return np.c_[hx, hy, hz] 

def test_kron_spmatrix(a_spmatrix, b_spmatrix):
    sa = sparse.random((3, 4), density=0.5)
    a = sa.todense()
    sb = sparse.random((5, 6), density=0.5)
    b = sb.todense()

    if a_spmatrix:
        sa = sa.tocsr()

    if b_spmatrix:
        sb = sb.tocsr()

    sol = np.kron(a, b)
    assert_eq(sparse.kron(sa, sb), sol)
    assert_eq(sparse.kron(sa, b), sol)
    assert_eq(sparse.kron(a, sb), sol)

    with pytest.raises(ValueError):
        assert_eq(sparse.kron(a, b), sol) 

def test_kron(a_ndim, b_ndim):
    a_shape = (2, 3, 4)[:a_ndim]
    b_shape = (5, 6, 7)[:b_ndim]

    sa = sparse.random(a_shape, density=0.5)
    a = sa.todense()
    sb = sparse.random(b_shape, density=0.5)
    b = sb.todense()

    sol = np.kron(a, b)
    assert_eq(sparse.kron(sa, sb), sol)
    assert_eq(sparse.kron(sa, b), sol)
    assert_eq(sparse.kron(a, sb), sol)

    with pytest.raises(ValueError):
        assert_eq(sparse.kron(a, b), sol) 

def cellGradBC(self):
        """
        The cell centered Gradient boundary condition matrix
        """
        if getattr(self, '_cellGradBC', None) is None:
            BC = self.setCellGradBC(self._cellGradBC_list)
            n = self.vnC
            if self.dim == 1:
                G = ddxCellGradBC(n[0], BC[0])
            elif self.dim == 2:
                G1 = sp.kron(speye(n[1]), ddxCellGradBC(n[0], BC[0]))
                G2 = sp.kron(ddxCellGradBC(n[1], BC[1]), speye(n[0]))
                G = sp.block_diag((G1, G2), format="csr")
            elif self.dim == 3:
                G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGradBC(n[0], BC[0]))
                G2 = kron3(speye(n[2]), ddxCellGradBC(n[1], BC[1]), speye(n[0]))
                G3 = kron3(ddxCellGradBC(n[2], BC[2]), speye(n[1]), speye(n[0]))
                G = sp.block_diag((G1, G2, G3), format="csr")
            # Compute areas of cell faces & volumes
            S = self.area
            V = self.aveCC2F*self.vol  # Average volume between adjacent cells
            self._cellGradBC = sdiag(S/V)*G
        return self._cellGradBC 

def sparse_cH(terms, ldim=2):
    """Construct a sparse cyclic nearest-neighbour Hamiltonian

    :param terms: List of nearst-neighbour terms (square array or MPO,
        see return value of :func:`cXY_local_terms`)
    :param ldim: Local dimension

    :returns: The Hamiltonian as sparse matrix

    """
    H = 0
    N = len(terms)
    for pos, term in enumerate(terms[:-1]):
        if hasattr(term, 'lt'):
            # Convert MPO to regular matrix
            term = term.to_array_global().reshape((ldim**2, ldim**2))
        left = sp.eye(ldim**pos)
        right = sp.eye(ldim**(N - pos - 2))
        H += sp.kron(left, sp.kron(term, right))
    # The last term acts on the first and last site.
    cyc = terms[-1]
    middle = sp.eye(ldim**pos)
    for i in range(cyc.ranks[0]):
        H += sp.kron(cyc.lt[0][0, ..., i], sp.kron(middle, cyc.lt[1][i, ..., 0]))
    return H 

def test_constructor2(self):
        # construct from dense

        # test zero mats
        for shape in [(1,1), (5,1), (1,10), (10,4), (3,7), (2,1)]:
            A = zeros(shape)
            assert_equal(bsr_matrix(A).todense(),A)
        A = zeros((4,6))
        assert_equal(bsr_matrix(A,blocksize=(2,2)).todense(),A)
        assert_equal(bsr_matrix(A,blocksize=(2,3)).todense(),A)

        A = kron([[1,0,2,0],[0,0,0,0],[0,0,4,5]], [[0,1,2],[3,0,5]])
        assert_equal(bsr_matrix(A).todense(),A)
        assert_equal(bsr_matrix(A,shape=(6,12)).todense(),A)
        assert_equal(bsr_matrix(A,blocksize=(1,1)).todense(),A)
        assert_equal(bsr_matrix(A,blocksize=(2,3)).todense(),A)
        assert_equal(bsr_matrix(A,blocksize=(2,6)).todense(),A)
        assert_equal(bsr_matrix(A,blocksize=(2,12)).todense(),A)
        assert_equal(bsr_matrix(A,blocksize=(3,12)).todense(),A)
        assert_equal(bsr_matrix(A,blocksize=(6,12)).todense(),A)

        A = kron([[1,0,2,0],[0,1,0,0],[0,0,0,0]], [[0,1,2],[3,0,5]])
        assert_equal(bsr_matrix(A,blocksize=(2,3)).todense(),A) 

def test_sparse_format_conversions(self):
        A = sparse.kron([[1,0,2],[0,3,4],[5,0,0]], [[1,2],[0,3]])
        D = A.todense()
        A = self.spmatrix(A)

        for format in ['bsr','coo','csc','csr','dia','dok','lil']:
            a = A.asformat(format)
            assert_equal(a.format,format)
            assert_array_equal(a.todense(), D)

            b = self.spmatrix(D+3j).asformat(format)
            assert_equal(b.format,format)
            assert_array_equal(b.todense(), D+3j)

            c = eval(format + '_matrix')(A)
            assert_equal(c.format,format)
            assert_array_equal(c.todense(), D) 

def _cellGradxStencil(self):
        # TODO: remove this hard-coding
        BC = ['neumann', 'neumann']
        if self.dim == 1:
            G1 = ddxCellGrad(self.nCx, BC)
        elif self.dim == 2:
            G1 = sp.kron(speye(self.nCy), ddxCellGrad(self.nCx, BC))
        elif self.dim == 3:
            G1 = kron3(
                speye(self.nCz),
                speye(self.nCy),
                ddxCellGrad(self.nCx, BC)
            )
        return G1 

def _nodalGradStencilx(self):
        """
        Stencil for the nodal grad in the x-direction (nodes to x-edges)
        """
        if self.dim == 1:
            Gx = ddx(self.nCx)
        elif self.dim == 2:
            Gx = sp.kron(speye(self.nNy), ddx(self.nCx))
        elif self.dim == 3:
            Gx = kron3(speye(self.nNz), speye(self.nNy), ddx(self.nCx))
        return Gx 

def _nodalLaplaciany(self):
        Hy = sdiag(1./self.hy)
        if self.dim == 1:
            return None
        elif self.dim == 2:
            Hy = sp.kron(Hy, speye(self.nNx))
        elif self.dim == 3:
            Hy = kron3(speye(self.nNz), Hy, speye(self.nNx))
        return Hy.T * self._nodalGradStencily * Hy 

def _nodalLaplacianx(self):
        Hx = sdiag(1./self.hx)
        if self.dim == 2:
            Hx = sp.kron(speye(self.nNy), Hx)
        elif self.dim == 3:
            Hx = kron3(speye(self.nNz), speye(self.nNy), Hx)
        return Hx.T * self._nodalGradStencilx * Hx 

def _nodalLaplacianStencily(self):
        warnings.warn('Laplacian has not been tested rigorously.')

        if self.dim == 1:
            return None

        Dy = ddx(self.nCy)
        Ly = - Dy.T * Dy

        if self.dim == 2:
            Ly = sp.kron(Ly, speye(self.nNx))
        elif self.dim == 3:
            Ly = kron3(speye(self.nNz), Ly, speye(self.nNx))
        return Ly 

def __get_A(N, dx):
        """
        Helper function to assemble FD matrix A in sparse format

        Args:
            N (list): number of dofs
            dx (float): distance between two spatial nodes

        Returns:
            scipy.sparse.csc_matrix: matrix A in CSC format
        """

        stencil = [1, -2, 1]
        zero_pos = 2
        dstencil = np.concatenate((stencil, np.delete(stencil, zero_pos - 1)))
        offsets = np.concatenate(([N[0] - i - 1 for i in reversed(range(zero_pos - 1))],
                                  [i - zero_pos + 1 for i in range(zero_pos - 1, len(stencil))]))
        doffsets = np.concatenate((offsets, np.delete(offsets, zero_pos - 1) - N[0]))

        A = sp.diags(dstencil, doffsets, shape=(N[0], N[0]), format='csc')
        A = sp.kron(A, sp.eye(N[0])) + sp.kron(sp.eye(N[1]), A)
        A *= 1.0 / (dx ** 2)

        return A

    # noinspection PyTypeChecker 

def _areaFyFull(self):
        """
        Area of y-faces (Azimuthal faces), prior to deflation.
        """
        return np.kron(self.hz, np.kron(np.ones(self._ntNy), self.hx)) 

def _faceDivStencilx(self):
        """
        Face divergence operator in the x-direction (x-faces to cell centers)
        """
        if self.dim == 1:
            Dx = ddx(self.nCx)
        elif self.dim == 2:
            Dx = sp.kron(speye(self.nCy), ddx(self.nCx))
        elif self.dim == 3:
            Dx = kron3(speye(self.nCz), speye(self.nCy), ddx(self.nCx))
        return Dx 

def finite_gs_energy(L, J, g):
    """For comparison: obtain ground state energy from exact diagonalization.

    Exponentially expensive in L, only works for small enough `L` <~ 20.
    """
    if L >= 20:
        warnings.warn("Large L: Exact diagonalization might take a long time!")
    # get single site operaors
    sx = sparse.csr_matrix(np.array([[0., 1.], [1., 0.]]))
    sz = sparse.csr_matrix(np.array([[1., 0.], [0., -1.]]))
    id = sparse.csr_matrix(np.eye(2))
    sx_list = []  # sx_list[i] = kron([id, id, ..., id, sx, id, .... id])
    sz_list = []
    for i_site in range(L):
        x_ops = [id] * L
        z_ops = [id] * L
        x_ops[i_site] = sx
        z_ops[i_site] = sz
        X = x_ops[0]
        Z = z_ops[0]
        for j in range(1, L):
            X = sparse.kron(X, x_ops[j], 'csr')
            Z = sparse.kron(Z, z_ops[j], 'csr')
        sx_list.append(X)
        sz_list.append(Z)
    H_xx = sparse.csr_matrix((2**L, 2**L))
    H_z = sparse.csr_matrix((2**L, 2**L))
    for i in range(L - 1):
        H_xx = H_xx + sx_list[i] * sx_list[(i + 1) % L]
    for i in range(L):
        H_z = H_z + sz_list[i]
    H = -J * H_xx - g * H_z
    E, V = arp.eigsh(H, k=1, which='SA', return_eigenvectors=True, ncv=20)
    return E[0]