Python scipy.sparse 模块，diags() 实例源码

def solve(self, x, w):
        # Check that w is a vector or a nx1 matrix
        if w.ndim == 2:
            assert(w.shape[1] == 1)
        elif w.dim == 1:
            w = w.reshape(w.shape[0], 1)
        A_smooth = (self.Dm - self.Dn.dot(self.grid.blur(self.Dn)))
        w_splat = self.grid.splat(w)
        A_data = diags(w_splat[:,0], 0)
        A = self.params["lam"] * A_smooth + A_data
        xw = x * w
        b = self.grid.splat(xw)
        # Use simple Jacobi preconditioner
        A_diag = np.maximum(A.diagonal(), self.params["A_diag_min"])
        M = diags(1 / A_diag, 0)
        # Flat initialization
        y0 = self.grid.splat(xw) / w_splat
        yhat = np.empty_like(y0)
        for d in xrange(x.shape[-1]):
            yhat[..., d], info = cg(A, b[..., d], x0=y0[..., d], M=M, maxiter=self.params["cg_maxiter"], tol=self.params["cg_tol"])
        xhat = self.grid.slice(yhat)
        return xhat

def generate_laplacian_score_scalar(X_ent, X_word, kNeighbors):
    # Generate cosine similarity graph
    n = X_ent.shape[0]
    cosX = cosine_similarity(X_word)
    graph = np.zeros((n, n))
    for i in range(n):
        for j in np.argpartition(cosX[i], -kNeighbors)[-kNeighbors:]:
            if j == i:
                continue
            graph[i, j] = cosX[i, j]
            graph[j, i] = cosX[i, j]

    D = sparse.diags([graph.sum(axis=0)], [0])
    L = D - graph
    f_tilde = X_ent - (float(X_ent.transpose() * D * np.ones((n, 1))) / D.sum().sum()) * np.ones((n, 1))
    score = float(f_tilde.transpose() * L * f_tilde) / float(f_tilde.transpose() * D * f_tilde + 1e-10)
    laplacian_score = score
    return laplacian_score

def iteration(self, user, fixed_vecs):
        num_solve = self.num_users if user else self.num_items
        num_fixed = fixed_vecs.shape[0]
        YTY = fixed_vecs.T.dot(fixed_vecs)
        eye = sparse.eye(num_fixed)
        lambda_eye = self.reg_param * sparse.eye(self.num_factors)
        solve_vecs = np.zeros((num_solve, self.num_factors))

        for i in xrange(num_solve):
            if user:
                matrix_i = self.matrix[i].toarray()
            else:
                matrix_i = self.matrix[:, i].T.toarray()
            CuI = sparse.diags(matrix_i, [0])
            pu = matrix_i.copy()
            pu[np.where(pu != 0)] = 1.0
            YTCuIY = fixed_vecs.T.dot(CuI).dot(fixed_vecs)
            YTCupu = fixed_vecs.T.dot(CuI + eye).dot(sparse.csr_matrix(pu).T)
            xu = spsolve(YTY + YTCuIY + lambda_eye, YTCupu)
            solve_vecs[i] = xu

        return solve_vecs

def _calculate_jacobian_mask(nx, ny, dof):
    from scipy.sparse import diags, kron

    N = nx * ny

    j_structure = diags(
        np.ones(shape=(7,)),
        [-nx + 1, -nx, -1, 0, 1, +nx, +nx - 1],
        shape=(N, N),
        format='csr',
    )

    j_structure = kron(
        j_structure,
        np.ones(shape=(dof, dof)),
        format='csr',
    )

    return j_structure

def hessian_wrt_mean(self):
        """ The Hessian of the multivariate Gaussian w.r.t. its mean, potentially including the
        linear projection.

        :return: The Hessian w.r.t. the mean
        :rtype: float|np.ndarray|sp.spmatrix
        """
        if self.hessian_mean_cache is None:
            hessian = self.hessian
            if self.mean_lin_op is not None:
                if np.isscalar(hessian) and isinstance(self.mean_lin_op, IndexOperator):
                    # index operator preserves diagonality
                    hessian = sp.diags(self.mean_lin_op.rmatvec(hessian * np.ones(self.dim)), 0)
                elif np.isscalar(hessian):
                    hessian = hessian * np.eye(self.dim)
                    hessian = rmatvec_nd(self.mean_lin_op, hessian)
                else:
                    hessian = rmatvec_nd(self.mean_lin_op, hessian)
            self.hessian_mean_cache = hessian
        return self.hessian_mean_cache

def make_linearOperator(shape, Xn, K):

    M,N = shape

    fx = K[0,0]
    fy = K[1,1]

    x_hat = Xn[0,:]
    y_hat = Xn[1,:]

    Kx,Ky = make_derivatives_2D_complete(shape)       # use one-sided differences with backward diff at image border
    Kx = Kx.tocsr()
    Ky = Ky.tocsr()

    spId = sparse.eye(M*N, M*N, format='csr')
    spXhat = sparse.diags(x_hat.flatten(), 0).tocsr()
    spYhat = sparse.diags(y_hat.flatten(), 0).tocsr()

    L = sparse.vstack([-Kx/fy, -Ky/fx, 
                       spXhat*Kx/fy + spYhat*Ky/fx +
                       2*spId/(fx*fy)
                        ])

    return L.tocsr()

def diags(self, format='dia'):
        """Return a regular sparse matrix of specified format

        kwargs:
            format  ('dia', 'csr')

        """
        if self._diags is None:
            self._diags = sp_diags(list(self.values()), list(self.keys()),
                                   shape=self.shape, format=format)

        if self._diags.format != format:
            self._diags = sp_diags(list(self.values()), list(self.keys()),
                                   shape=self.shape, format=format)

        return self._diags

def fit(self, X, y, L):
        """Fit the model according to the given training data.

        Prameters
        ---------
        X : array-like, shpae = [n_samples, n_features]
            Training data.

        y : array-like, shpae = [n_samples]
            Target values (unlabeled points are marked as 0).

        L : array-like, shpae = [n_samples, n_samples]
            Graph Laplacian.
        """

        labeled               = y != 0
        y_labeled             = y[labeled]
        n_samples, n_features = X.shape
        n_labeled_samples     = y_labeled.size
        I                     = sp.eye(n_samples)
        J                     = sp.diags(labeled.astype(np.float64))
        K                     = rbf_kernel(X, gamma=self.gamma_k)
        M                     = J @ K \
            + self.gamma_a * n_labeled_samples * I \
            + self.gamma_i * n_labeled_samples / n_samples**2 * L**self.p @ K

        # Train a classifer
        self.dual_coef_       = LA.solve(M, y)

        return self

def preprocess_graph(adj):
    adj = sp.coo_matrix(adj)
    adj_ = adj + sp.eye(adj.shape[0])
    rowsum = np.array(adj_.sum(1))
    degree_mat_inv_sqrt = sp.diags(np.power(rowsum, -0.5).flatten())
    adj_normalized = adj_.dot(degree_mat_inv_sqrt).transpose().dot(degree_mat_inv_sqrt).tocoo()
    return sparse_to_tuple(adj_normalized)

def getDiffusionMap(SSM, Kappa, t = -1, includeDiag = True, thresh = 5e-4, NEigs = 51):
    """
    :param SSM: Metric between all pairs of points
    :param Kappa: Number in (0, 1) indicating a fraction of nearest neighbors
                used to autotune neighborhood size
    :param t: Diffusion parameter.  If -1, do Autotuning
    :param includeDiag: If true, include recurrence to diagonal in the markov
        chain.  If false, zero out diagonal
    :param thresh: Threshold below which to zero out entries in markov chain in
        the sparse approximation
    :param NEigs: The number of eigenvectors to use in the approximation
    """
    N = SSM.shape[0]
    #Use the letters from the delaPorte paper
    K = getW(SSM, int(Kappa*N))
    if not includeDiag:
        np.fill_diagonal(K, np.zeros(N))
    RowSumSqrt = np.sqrt(np.sum(K, 1))
    DInvSqrt = sparse.diags([1/RowSumSqrt], [0])

    #Symmetric normalized Laplacian
    Pp = (K/RowSumSqrt[None, :])/RowSumSqrt[:, None]
    Pp[Pp < thresh] = 0
    Pp = sparse.csr_matrix(Pp)

    lam, X = sparse.linalg.eigsh(Pp, NEigs, which='LM')
    lam = lam/lam[-1] #In case of numerical instability

    #Check to see if autotuning
    if t > -1:
        lamt = lam**t
    else:
        #Autotuning diffusion time
        lamt = np.array(lam)
        lamt[0:-1] = lam[0:-1]/(1-lam[0:-1])

    #Do eigenvector version
    V = DInvSqrt.dot(X) #Right eigenvectors
    M = V*lamt[None, :]
    return M/RowSumSqrt[:, None] #Put back into orthogonal Euclidean coordinates

def getC(mask, c):
    scribble_mask = mask != 0
    return c * sparse.diags(c * scribble_mask.ravel())

def test_solve_range_range(dtype):
    b = np.arange(solver.N, dtype=dtype)
    mat = solver.ones(dtype=dtype) + sp.diags([b], offsets=[0], shape=(solver.N, solver.N),
                                              format='csr', dtype=dtype)
    obs = solver.solve(mat, b)
    exp = spla.spsolve(mat, b)
    assert np.allclose(exp, obs)

def get_tfidf_matrix(cnts):
    """Convert the word count matrix into tfidf one.

    tfidf = log(tf + 1) * log((N - Nt + 0.5) / (Nt + 0.5))
    * tf = term frequency in document
    * N = number of documents
    * Nt = number of occurences of term in all documents
    """
    Ns = get_doc_freqs(cnts)
    idfs = np.log((cnts.shape[1] - Ns + 0.5) / (Ns + 0.5))
    idfs[idfs < 0] = 0
    idfs = sp.diags(idfs, 0)
    tfs = cnts.log1p()
    tfidfs = idfs.dot(tfs)
    return tfidfs

def normalize_adj(adj):
    """Symmetrically normalize adjacency matrix."""
    adj = sp.coo_matrix(adj)
    rowsum = np.array(adj.sum(1))
    d_inv_sqrt = np.power(rowsum, -0.5).flatten()
    d_inv_sqrt[np.isinf(d_inv_sqrt)] = 0.
    d_mat_inv_sqrt = sp.diags(d_inv_sqrt)
    return adj.dot(d_mat_inv_sqrt).transpose().dot(d_mat_inv_sqrt).tocoo()

def normalize_adj(adj):
    """Symmetrically normalize adjacency matrix."""
    adj = sp.coo_matrix(adj)
    rowsum = np.array(adj.sum(1))
    d_inv_sqrt = np.power(rowsum, -0.5).flatten()
    d_inv_sqrt[np.isinf(d_inv_sqrt)] = 0.
    d_mat_inv_sqrt = sp.diags(d_inv_sqrt)
    return adj.dot(d_mat_inv_sqrt).transpose().dot(d_mat_inv_sqrt).tocoo()

def variational_kalman( observations, mask, state_mask, uncertainty, H_matrix, n_params,
            x_forecast, P_forecast, P_forecast_inv, the_metadata, approx_diagonal=True):
    """We can just use """
    if len(H_matrix) == 2:
        non_linear = True
        H0, H_matrix = H_matrix
    else:
        H0 = 0.
        non_linear = False
    R_mat = sp.diags(uncertainty.diagonal()[state_mask.flatten()])
    LOG.info("Creating linear problem")
    y = observations[state_mask]
    y = np.where(mask[state_mask], y, 0.)
    if non_linear:
        y = y + H_matrix.dot(x_forecast) - H0


    #Aa = matrix_squeeze (P_forecast_inv, mask=maska.ravel())
    A = H_matrix.T.dot(R_mat).dot(H_matrix) + P_forecast_inv
    b = H_matrix.T.dot(R_mat).dot(y) + P_forecast_inv.dot (x_forecast)
    # Here we can either do a spLU of A, and solve, or we can have a first go
    # by assuming P_forecast_inv is diagonal, and use the inverse of A_approx as
    # a preconditioner
    LOG.info("Solving")
    AI = sp.linalg.splu ( A )
    x_analysis = AI.solve (b)
    # So retval is the solution vector and A is the Hessian 
    # (->inv(A) is posterior cov)
    innovations = y - H_matrix.dot(x_analysis)

    LOG.info("Inflating analysis state")
    #x_analysis = reconstruct_array ( x_analysis_prime, x_forecast,
    #                                    mask.ravel(), n_params=n_params)

    return x_analysis, None, A, innovations

def on_changed(self, which):
            # pylint: disable=access-member-before-definition, attribute-defined-outside-init

            if not hasattr(self, 'normalized'):
                self.normalized = True

            if hasattr(self, 'v') and hasattr(self, 'f'):
                if 'f' not in which and hasattr(self, 'faces_by_vertex') and self.faces_by_vertex.shape[0]/3 == self.v.shape[0]:
                    self.tns.v = self.v
                else: # change in f or in size of v. shouldn't happen often.
                    f = self.f

                    IS = f.ravel()
                    JS = np.array([range(f.shape[0])]*3).T.ravel()
                    data = np.ones(len(JS))

                    IS = np.concatenate((IS*3, IS*3+1, IS*3+2))
                    JS = np.concatenate((JS*3, JS*3+1, JS*3+2))
                    data = np.concatenate((data, data, data))

                    sz = self.v.size
                    self.faces_by_vertex = sp.csc_matrix((data, (IS, JS)), shape=(sz, f.size))

                    self.tns = ch.Ch(lambda v: CrossProduct(TriEdges(f, 1, 0, v), TriEdges(f, 2, 0, v)))
                    self.tns.v = self.v

                    if self.normalized:
                        tmp = ch.MatVecMult(self.faces_by_vertex, self.tns)
                        self.normals = NormalizedNx3(tmp)
                    else:
                        test = self.faces_by_vertex.dot(np.ones(self.faces_by_vertex.shape[1]))
                        faces_by_vertex = sp.diags([1. / test], [0]).dot(self.faces_by_vertex).tocsc()
                        normals = ch.MatVecMult(faces_by_vertex, self.tns).reshape((-1, 3))
                        normals = normals / (ch.sum(normals ** 2, axis=1) ** .25).reshape((-1, 1))
                        self.normals = normals

def _degrees_inv(g):
    """
    Given graph g. Returns inverse of degree matrix
    where 0 entries are ignored.
    """
    n, degs = g.order(), g.degree().values()
    degs = list(map(lambda x : 0 if x == 0 else 1.0/float(x), degs))
    return sp.diags(degs, format = "csc")

def normalize_rows(self, matrix):
        """
        Normalizes the rows of the given matrix, such that every row sums up to 1

        :param matrix:
            Matrix of which the rows should be normalized
        :return:
            Row-normalized version of the given matrix
        """

        # We'll normalize the row to sum up to 1 by pre-multiplying by a matrix
        # that contains the scalars for every row on the diagonal
        scalars = [1.0 / matrix.getrow(i).sum() for i in range(matrix.shape[0])]

        return diags(scalars, 0) * matrix

def subtract_rowsum_on_diag(spmat):
    spmat = spmat.tocsr() - sparse.diags(numpy.array(spmat.sum(axis=1)).T, offsets=[0], format="csr")
    return spmat.tocsr()

def emptyGraph(hin):
        n = len(hin.Ids)
        return sparse.diags(np.zeros(n))

def generate_meta_graph(scope, scope_name, type_list, count):
    split_path = 'data/local/split/' + scope_name + '/'
    pred_path = 'data/local/metagraph/' + scope_name + '/'
    if not os.path.exists('data/local/metagraph/' + scope_name + '/'):
        os.makedirs('data/local/metagraph/' + scope_name + '/')

    with open('data/local/' + scope_name + '.dmp') as f:
        hin = pk.load(f)
    tf_param = {'word':True, 'entity':True, 'we_weight':0.1}

    for t in type_list:
        #print t
        X, newIds, entitynewIds = GraphGenerator.getTFVectorX(hin, tf_param, t)
        n = X.shape[0]
        e = X.shape[1]
        with open('data/local/laplacian/' + scope_name + '/' + str(t) + '_scores') as f:
            laplacian_score = pk.load(f)
        laplacian_score = 20 * np.exp(-laplacian_score * 0.01)
        D = sparse.diags(laplacian_score)
        X = X * D
        X = X.toarray()
        graph = np.zeros((n+e,n+e))
        graph[0:n,n:n+e] = X
        graph[n:n+e,0:n] = X.transpose()
        graph = sparse.csc_matrix(graph)

        newLabel = GraphGenerator.getNewLabels(hin)
        lp_param = {'alpha':0.98, 'normalization_factor':5}
        # 3-class classification
        lp_candidate = [5]
        for lp in lp_candidate:
            ssl = SSLClassifier(graph, newLabel, scope, lp_param, repeatTimes=50, trainNumbers=lp, classCount=count)
            if not os.path.exists(pred_path + str(t) + '/'):
                os.makedirs(pred_path + str(t) + '/')
            ssl.repeatedFixedExperimentwithNewIds(pathPrefix=split_path + 'lb' + str(lp).zfill(3) + '_',newIds=newIds, saveProb=True,savePathPrefix=pred_path + str(t) + '/' + 'lb' + str(lp).zfill(3))

def generate_laplacian_score(X_ent, X_word, kNeighbors):
    # Generate cosine similarity graph
    n = X_ent.shape[0]
    m = X_ent.shape[1]
    cosX = cosine_similarity(X_word)
    graph = np.zeros((n, n))
    t = cosX.sum().sum() / n/n
    for i in range(n):
        for j in np.argpartition(cosX[i], -kNeighbors)[-kNeighbors:]:
            if j == i:
                continue
            # diff = (X_word[i, :] - X_word[j, :]).toarray().flatten()

            # dist = np.exp(np.dot(diff, diff) / t)
            graph[i, j] = cosX[i, j] #np.exp(- (1 - cosX[i, j]) / 0.03) #
            graph[j, i] = cosX[i, j] #np.exp(- (1 - cosX[i, j]) / 0.03) #

    D = sparse.diags([graph.sum(axis=0)], [0])
    L = D - graph

    laplacian_score = np.zeros(m)
    for i in range(m):
        f_tilde = X_ent[:, i] - (float(X_ent[:, i].transpose() * D * np.ones((n, 1))) / D.sum().sum()) * np.ones(
            (n, 1))
        score = float(f_tilde.transpose() * L * f_tilde) / float(f_tilde.transpose() * D * f_tilde + 1e-10)
        laplacian_score[i] = score


    return (laplacian_score)

def generate_laplacian_score(X_ent, X_word, kNeighbors):
    # Generate cosine similarity graph
    n = X_ent.shape[0]
    m = X_ent.shape[1]
    cosX = cosine_similarity(X_word)
    graph = np.zeros((n, n))
    t = cosX.sum().sum() / n/n
    for i in range(n):
        for j in np.argpartition(cosX[i], -kNeighbors)[-kNeighbors:]:
            if j == i:
                continue
            # diff = (X_word[i, :] - X_word[j, :]).toarray().flatten()

            # dist = np.exp(np.dot(diff, diff) / t)
            graph[i, j] = cosX[i, j] #np.exp(- (1 - cosX[i, j]) / 0.03) #
            graph[j, i] = cosX[i, j] #np.exp(- (1 - cosX[i, j]) / 0.03) #

    D = sparse.diags([graph.sum(axis=0)], [0])
    L = D - graph

    laplacian_score = np.zeros(m)
    for i in range(m):
        f_tilde = X_ent[:, i] - (float(X_ent[:, i].transpose() * D * np.ones((n, 1))) / D.sum().sum()) * np.ones(
            (n, 1))
        score = float(f_tilde.transpose() * L * f_tilde) / float(f_tilde.transpose() * D * f_tilde + 1e-10)
        laplacian_score[i] = score


    return (laplacian_score)

def extract_DF_F(Y, A, C, i=None):
    """Extract DF/F values from spatial/temporal components and background

     Parameters
     -----------
     Y: np.ndarray
           input data (d x T)
     A: sparse matrix of np.ndarray
           Set of spatial including spatial background (d x K)
     C: matrix
           Set of temporal components including background (K x T)

     Returns
     -----------
     C_df: matrix
          temporal components in the DF/F domain
     Df:  np.ndarray
          vector with baseline values for each trace
    """
    A2 = A.copy()
    A2.data **= 2
    nA2 = np.squeeze(np.array(A2.sum(axis=0)))
    A = A * diags(1 / nA2, 0)
    C = diags(nA2, 0) * C

    # if i is None:
    #    i = np.argmin(np.max(A,axis=0))

    Y = np.matrix(Y)
    Yf = A.transpose() * (Y - A * C)  # + A[:,i]*C[i,:])
    Df = np.median(np.array(Yf), axis=1)
    C_df = diags(1 / Df, 0) * C

    return C_df, Df

def Diag(self, v, **kwargs):
        """ A := diag(v) """
        v = np.require(v, requirements='F')
        if v.ndim > 1:
            v = v.flatten(order='A')
        dtype = kwargs.get('dtype', np.dtype('complex64'))
        M = spp.diags( v, offsets=0 ).astype(dtype)
        return self.SpMatrix(M, **kwargs)

def TikhonovDiff(f, dx, lam, d = 1):
    """
    Tikhonov differentiation.

    return argmin_g \|Ag-f\|_2^2 + lam*\|Dg\|_2^2
    where A is trapezoidal integration and D is finite differences for first dervative

    It looks like it will work well and does for the ODE case but 
    tends to introduce too much bias to work well for PDEs.  If the data is noisy, try using
    polynomials instead.
    """

    # Initialize a few things    
    n = len(f)
    f = np.matrix(f - f[0]).reshape((n,1))

    # Get a trapezoidal approximation to an integral
    A = np.zeros((n,n))
    for i in range(1, n):
        A[i,i] = dx/2
        A[i,0] = dx/2
        for j in range(1,i): A[i,j] = dx

    e = np.ones(n-1)
    D = sparse.diags([e, -e], [1, 0], shape=(n-1, n)).todense() / dx

    # Invert to find derivative
    g = np.squeeze(np.asarray(np.linalg.lstsq(A.T.dot(A) + lam*D.T.dot(D),A.T.dot(f))[0]))

    if d == 1: return g

    # If looking for a higher order derivative, this one should be smooth so now we can use finite differences
    else: return FiniteDiff(g, dx, d-1)

def get_tfidf_matrix(cnts):
    """Convert the word count matrix into tfidf one.

    tfidf = log(tf + 1) * log((N - Nt + 0.5) / (Nt + 0.5))
    * tf = term frequency in document
    * N = number of documents
    * Nt = number of occurences of term in all documents
    """
    Ns = get_doc_freqs(cnts)
    idfs = np.log((cnts.shape[1] - Ns + 0.5) / (Ns + 0.5))
    idfs[idfs < 0] = 0
    idfs = sp.diags(idfs, 0)
    tfs = cnts.log1p()
    tfidfs = idfs.dot(tfs)
    return tfidfs

def normalize_adj(adj, symmetric=True):
    if symmetric:
        d = sp.diags(np.power(np.array(adj.sum(1)), -0.5).flatten(), 0)
        a_norm = adj.dot(d).transpose().dot(d).tocsr()
    else:
        d = sp.diags(np.power(np.array(adj.sum(1)), -1).flatten(), 0)
        a_norm = d.dot(adj).tocsr()
    return a_norm

def bistochastize(grid, maxiter=10):
    """Compute diagonal matrices to bistochastize a bilateral grid"""
    m = grid.splat(np.ones(grid.npixels))
    n = np.ones(grid.nvertices)
    for i in xrange(maxiter):
        n = np.sqrt(n * m / grid.blur(n))
    # Correct m to satisfy the assumption of bistochastization regardless
    # of how many iterations have been run.
    m = n * grid.blur(n)
    Dm = diags(m, 0)
    Dn = diags(n, 0)
    return Dn, Dm

def factor_aug(z, DPhival, G, A):    
    M, N = G.shape
    P, N = A.shape
    """Multiplier for inequality constraints"""
    l = z[N+P:N+P+M]

    """Slacks"""
    s = z[N+P+M:]

    """Sigma matrix"""
    SIG = diags(l/s, 0)

    """Condensed system"""
    if issparse(DPhival):
        if not issparse(A):
            A = csr_matrix(A)        
        H = DPhival + mydot(G.T, mydot(SIG, G))
        J = bmat([[H, A.T], [A, None]])
    else:
        if issparse(A):
            A = A.toarray()
        J = np.zeros((N+P, N+P))
        J[0:N, 0:N] = DPhival + mydot(G.T, mydot(SIG, G))            
        J[0:N, N:] = A.T
        J[N:, 0:N] = A

    LU = myfactor(J)    
    return LU

def factor_schur(z, DPhival, G, A):
    M, N = G.shape
    P, N = A.shape
    """Multiplier for inequality constraints"""
    l = z[N+P:N+P+M]

    """Slacks"""
    s = z[N+P+M:]

    """Sigma matrix"""
    SIG = diags(l/s, 0)

    """Augmented Jacobian"""
    H = DPhival + mydot(G.T, mydot(SIG, G))

    """Factor H"""
    LU_H = myfactor(H)

    """Compute H^{-1}A^{T}"""
    HinvAt = mysolve(LU_H, A.T)

    """Compute Schur complement AH^{-1}A^{T}"""
    S = mydot(A, HinvAt)

    """Factor Schur complement"""
    LU_S = myfactor(S)

    LU = (LU_S, LU_H)
    return LU

def contrib_prod_e2_x(docs, test_x, n_docs):
    ''' Compute contribution of the partition to the product of E2 and X

    Parameters
    -----------
    docs : m-by-vocab_size array or csr_matrix
        Current partition of word count vectors.
    test_x : vocab_size-by-k array
        Test matrix where k is the number of factors.
    n_docs : int
        Total number of documents, could be greater than m.

    Returns
    ----------
    out : vocab_size-by-k array
        Contribution of the partition to the product of E2 and X.
    '''
    partition_size, vocab_size = docs.shape
    _vocab_size, num_factors = test_x.shape
    assert partition_size >= 1 and vocab_size >= 1
    assert partition_size <= n_docs
    assert vocab_size == _vocab_size and num_factors >= 1

    docs = docs[np.squeeze(np.array(docs.sum(axis=1))) >= 3]

    document_lengths = np.squeeze(np.array(docs.sum(axis=1)))
    diag_l = spsp.diags(1.0 / document_lengths / (document_lengths - 1))

    scaled_docs = diag_l.dot(docs)
    prod_x = scaled_docs.T.dot(docs.dot(test_x))

    sum_scaled_docs = np.squeeze(np.array(scaled_docs.sum(axis=0)))
    prod_x_adj = spsp.diags(sum_scaled_docs).dot(test_x)

    return (prod_x - prod_x_adj) / n_docs

def test_get_posterior_hessian(self):
        """ Tests the computation of the log-posterior Hessian with a simple graph,
              X
             /|
            / |
           v  v
           Y  Z
        where all nodes contain 2D Gaussians and node X encodes the mean for nodes Y and Z
        """
        for k in range(NUM_TESTS):
            prec_x = np.diag(np.random.rand(2), 0)
            prec_y = sp.diags(np.random.rand(2), 0)  # throw in a sparse precision
            prec_z = np.diag(np.random.rand(2), 0)
            node_x = Node(name='x', data=np.random.randn(2), cpd=GaussianCPD(precision=prec_x))
            node_y = Node(name='y', data=np.random.randn(2), cpd=GaussianCPD(precision=prec_y),
                          param_nodes={GaussianCPD.MEAN_KEY: node_x})
            node_z = Node(name='z', data=np.random.randn(2), cpd=GaussianCPD(precision=prec_z),
                          param_nodes={GaussianCPD.MEAN_KEY: node_x}, held_out=True)
            learner = undertest.BayesNetLearner(nodes=[node_x, node_y, node_z])
            np.testing.assert_almost_equal(learner.get_posterior_hessian('x', use_held_out=True),
                                           -prec_x - prec_y - prec_z)
            np.testing.assert_almost_equal(learner.get_posterior_hessian('x', use_held_out=False),
                                           -prec_x - prec_y)
            np.testing.assert_almost_equal(learner.get_posterior_hessian('x'), -prec_x - prec_y)
            np.testing.assert_almost_equal(learner.get_posterior_hessian('x', np.random.randn(2)),
                                           -prec_x - prec_y)

def counts_to_ppmi(counts_csr, smoothing=0.75):
    """
    Converts a sparse matrix of co-occurrences into a sparse matrix of positive
    pointwise mutual information. Context distributional smoothing is applied
    to the resulting matrix.
    """
    # word_counts adds up the total amount of association for each term.
    word_counts = np.asarray(counts_csr.sum(axis=1)).flatten()

    # smooth_context_freqs represents the relative frequency of occurrence
    # of each term as a context (a column of the table).
    smooth_context_freqs = np.asarray(counts_csr.sum(axis=0)).flatten() ** smoothing
    smooth_context_freqs /= smooth_context_freqs.sum()

    # Divide each row of counts_csr by the word counts. We accomplish this by
    # multiplying on the left by the sparse diagonal matrix of 1 / word_counts.
    ppmi = sparse.diags(1 / word_counts).dot(counts_csr)

    # Then, similarly divide the columns by smooth_context_freqs, by the same
    # method except that we multiply on the right.
    ppmi = ppmi.dot(sparse.diags(1 / smooth_context_freqs))

    # Take the log of the resulting entries to give pointwise mutual
    # information. Discard those whose PMI is less than 0, to give positive
    # pointwise mutual information (PPMI).
    ppmi.data = np.maximum(np.log(ppmi.data), 0)
    ppmi.eliminate_zeros()
    return ppmi

def fit(self, X, y, L):
        """Fit the model according to the given training data.

        Prameters
        ---------
        X : array-like, shpae = [n_samples, n_features]
            Training data.

        y : array-like, shpae = [n_samples]
            Target values (unlabeled points are marked as 0).

        L : array-like, shpae = [n_samples, n_samples]
            Graph Laplacian.
        """

        labeled               = y != 0
        y_labeled             = y[labeled]
        n_samples, n_features = X.shape
        n_labeled_samples     = y_labeled.size
        I                     = sp.eye(n_samples)
        Y                     = sp.diags(y_labeled)
        J                     = sp.eye(n_labeled_samples, n_samples)
        K                     = rbf_kernel(X, gamma=self.gamma_k)
        M                     = 2 * self.gamma_a * I \
            + 2 * self.gamma_i / n_samples**2 * L**self.p @ K

        # Construct the QP, invoke solver
        solvers.options['show_progress'] = False
        sol                   = solvers.qp(
            P                 = matrix(Y @ J @ K @ LA.inv(M) @ J.T @ Y),
            q                 = matrix(-1 * np.ones(n_labeled_samples)),
            G                 = matrix(np.vstack((
                -1 * np.eye(n_labeled_samples),
                n_labeled_samples * np.eye(n_labeled_samples)
            ))),
            h                 = matrix(np.hstack((
                np.zeros(n_labeled_samples),
                np.ones(n_labeled_samples)
            ))),
            A                 = matrix(y_labeled, (1, n_labeled_samples), 'd'),
            b                 = matrix(0.0)
        )

        # Train a classifer
        self.dual_coef_       = LA.solve(M, J.T @ Y @ np.array(sol['x']).ravel())

        return self

def factor_aug(z, DPhival, G, A):
    r"""Set up augmented system and return.

    Parameters
    ----------
    z : (N+P+M+M,) ndarray
        Current iterate, z = (x, nu, l, s)
    DPhival : LinearOperator
        Jacobian of the variational inequality mapping
    G : (M, N) ndarray or sparse matrix
        Inequality constraints
    A : (P, N) ndarray or sparse matrix
        Equality constraints  

    Returns
    -------
    J : LinearOperator
        Augmented system

    """
    M, N = G.shape
    P, N = A.shape
    """Multiplier for inequality constraints"""
    l = z[N+P:N+P+M]

    """Slacks"""
    s = z[N+P+M:]

    """Sigma matrix"""
    SIG = diags(l/s, 0)
    # SIG = diags(l*s, 0)

    """Convert A"""
    if not issparse(A):
        A = csr_matrix(A)

    """Convert G"""
    if not issparse(G):
        G = csr_matrix(G)

    """Since we expect symmetric DPhival, we need to change A"""
    sign = np.zeros(N)
    sign[0:N/2] = 1.0
    sign[N/2:] = -1.0
    T = diags(sign, 0)

    A_new = A.dot(T)

    W = AugmentedSystem(DPhival, G, SIG, A_new)
    return W

def solve_factorized_aug(z, Fval, LU, G, A):
    M, N=G.shape
    P, N=A.shape

    """Total number of inequality constraints"""
    m = M    

    """Primal variable"""
    x = z[0:N]

    """Multiplier for equality constraints"""
    nu = z[N:N+P]

    """Multiplier for inequality constraints"""
    l = z[N+P:N+P+M]

    """Slacks"""
    s = z[N+P+M:]

    """Dual infeasibility"""
    rd = Fval[0:N]

    """Primal infeasibility"""
    rp1 = Fval[N:N+P]
    rp2 = Fval[N+P:N+P+M]

    """Centrality"""
    rc = Fval[N+P+M:]

    """Sigma matrix"""
    SIG = diags(l/s, 0)

    """RHS for condensed system"""
    b1 = -rd - mydot(G.T, mydot(SIG, rp2)) + mydot(G.T, rc/s)
    b2 = -rp1
    b = np.hstack((b1, b2))
    dxnu = mysolve(LU, b)
    dx = dxnu[0:N]
    dnu = dxnu[N:]

    """Obtain search directions for l and s"""
    ds = -rp2 - mydot(G, dx)
    dl = -mydot(SIG, ds) - rc/s

    dz = np.hstack((dx, dnu, dl, ds))
    return dz

###############################################################################
# Solve via normal equations (Schur complement)
###############################################################################

def solve_factorized_schur(z, Fval, LU, G, A):
    M, N=G.shape
    P, N=A.shape

    """Total number of inequality constraints"""
    m = M    

    """Primal variable"""
    x = z[0:N]

    """Multiplier for equality constraints"""
    nu = z[N:N+P]

    """Multiplier for inequality constraints"""
    l = z[N+P:N+P+M]

    """Slacks"""
    s = z[N+P+M:]

    """Dual infeasibility"""
    rd = Fval[0:N]

    """Primal infeasibility"""
    rp1 = Fval[N:N+P]
    rp2 = Fval[N+P:N+P+M]

    """Centrality"""
    rc = Fval[N+P+M:]

    """Sigma matrix"""
    SIG = diags(l/s, 0)

    """Assemble right hand side of augmented system"""
    r1 = rd + mydot(G.T, mydot(SIG, rp2)) - mydot(G.T, rc/s)
    r2 = rp1

    """Unpack LU-factors"""
    LU_S, LU_H = LU

    """Assemble right hand side for normal equation"""
    b = r2 - mydot(A, mysolve(LU_H, r1))  

    """Solve for dnu"""
    dnu = mysolve(LU_S, b)

    """Solve for dx"""
    dx = mysolve(LU_H, -(r1 + mydot(A.T, dnu)))    

    """Obtain search directions for l and s"""
    ds = -rp2 - mydot(G, dx)
    dl = -mydot(SIG, ds) - rc/s

    dz = np.hstack((dx, dnu, dl, ds))
    return dz

def test_rmatvec_nd(self):
        """ Test that given an n x k linear operator and n x n matrix rmatvec_nd yields a k x k
        matrix"""

        def rev_index(index_map, x, output_dim):
            intermediate = np.empty((output_dim, x.shape[1]))
            final = np.empty((output_dim, output_dim))
            for i in range(x.shape[1]):
                intermediate[:, i] = np.bincount(index_map, weights=x[:, i], minlength=output_dim)
            for i in range(output_dim):
                final[i, :] = np.bincount(index_map, weights=intermediate[i, :],
                                          minlength=output_dim)
            return final

        n = 10
        x = np.random.randn(n, n)
        k = np.random.randint(1, 5)
        index_map = np.random.randint(k, size=n)
        lin_op = undertest.IndexOperator(index_map=index_map, dim_x=k)
        actual = undertest.rmatvec_nd(lin_op, x)
        expected_backproj = rev_index(index_map, x, k)
        np.testing.assert_array_equal(actual, expected_backproj)

        # Sparse, non-diagonal
        x_sp = sp.csr_matrix(x)
        actual = undertest.rmatvec_nd(lin_op, x_sp)
        np.testing.assert_array_equal(actual, expected_backproj)

        # Sparse diagonal
        x_sp_diag = sp.diags(np.diag(x), 0)
        actual = undertest.rmatvec_nd(lin_op, x_sp_diag)
        self.assertEqual(actual.shape, (k, k))
        expected_backproj = np.diag(np.bincount(index_map, weights=np.diag(x), minlength=k))
        np.testing.assert_array_equal(actual, expected_backproj)

        # Non-sparse diagonal
        x_diag = np.diag(np.random.randn(n))
        actual = undertest.rmatvec_nd(lin_op, x_diag)
        self.assertEqual(actual.shape, (k, k))
        # The result should also be sparse and diagonal
        expected_backproj = np.diag(np.bincount(index_map, weights=np.diag(x_diag), minlength=k))
        np.testing.assert_array_equal(actual, expected_backproj)

        # scalar
        x = 1.3
        k = 5
        index_map = np.random.randint(k, size=1)
        lin_op = undertest.IndexOperator(index_map=index_map, dim_x=k)
        actual = undertest.rmatvec_nd(lin_op, x)
        expected_backproj = np.zeros(k)
        expected_backproj[index_map] = x
        np.testing.assert_array_equal(actual, expected_backproj)