def _create_sparse_precision(X, graph, n_features, n_features_per_vertex,
                             mode='concatenation', dtype=np.float32,
                             n_components=None, bias=0,
                             return_covariances=False, verbose=False):
    # check mode argument
    if mode not in ['concatenation', 'subtraction']:
        raise ValueError("mode must be either ''concatenation'' "
                         "or ''subtraction''; {} is given.".format(mode))

    # Initialize arrays
    all_blocks = np.zeros((graph.n_edges * 4,
                           n_features_per_vertex, n_features_per_vertex),
                          dtype=dtype)
    if return_covariances:
        if mode == 'concatenation':
            cov_shape = (graph.n_edges,
                         2 * n_features_per_vertex, 2 * n_features_per_vertex)
        else:
            cov_shape = (graph.n_edges,
                         n_features_per_vertex, n_features_per_vertex)
        all_covariances = np.zeros(cov_shape, dtype=dtype)
    columns = np.zeros(graph.n_edges - 4)
    rows = np.zeros(graph.n_edges * 4)

    # Print information if asked
    if verbose:
        edges = print_progress(range(graph.n_edges), n_items=graph.n_edges,
                               prefix='Precision per edge',
                               end_with_newline=False)
    else:
        edges = range(graph.n_edges)

    # Compute covariance matrix for each edge, invert it and store it
    count = -1
    for e in edges:
        # edge vertices
        v1 = graph.edges[e, 0]
        v2 = graph.edges[e, 1]

        # find indices in data matrix
        v1_from = v1 * n_features_per_vertex
        v1_to = (v1 + 1) * n_features_per_vertex
        v2_from = v2 * n_features_per_vertex
        v2_to = (v2 + 1) * n_features_per_vertex

        # data concatenation
        if mode == 'concatenation':
            edge_data = X[:, list(range(v1_from, v1_to)) +
                             list(range(v2_from, v2_to))]
        else:
            edge_data = X[:, v1_from:v1_to] - X[:, v2_from:v2_to]

        # compute covariance matrix
        covmat = np.cov(edge_data, rowvar=0, bias=bias)
        if return_covariances:
            all_covariances[e] = covmat

        # invert it
        covmat = _covariance_matrix_inverse(covmat, n_components)

        # store it
        if mode == 'concatenation':
            # v1, v1
            count += 1
            all_blocks[count] = \
                covmat[:n_features_per_vertex, :n_features_per_vertex]
            rows[count] = v1
            columns[count] = v1
            # v2, v2
            count += 1
            all_blocks[count] = \
                covmat[n_features_per_vertex::, n_features_per_vertex::]
            rows[count] = v2
            columns[count] = v2
            # v1, v2
            count += 1
            all_blocks[count] = \
                covmat[:n_features_per_vertex, n_features_per_vertex::]
            rows[count] = v1
            columns[count] = v2
            # v2, v1
            count += 1
            all_blocks[count] = \
                covmat[n_features_per_vertex::, :n_features_per_vertex]
            rows[count] = v2
            columns[count] = v1
        else:
            # v1, v1
            count += 1
            all_blocks[count] = covmat
            rows[count] = v1
            columns[count] = v1
            # v2, v2
            count += 1
            all_blocks[count] = covmat
            rows[count] = v2
            columns[count] = v2
            # v1, v2
            count += 1
            all_blocks[count] = -covmat
            rows[count] = v1
            columns[count] = v2
            # v2, v1
            count += 1
            all_blocks[count] = -covmat
            rows[count] = v2
            columns[count] = v1

    # sort rows, columns and all_blocks
    rows_arg_sort = rows.argsort()
    columns = columns[rows_arg_sort]
    all_blocks = all_blocks[rows_arg_sort]
    rows = rows[rows_arg_sort]

    # create indptr
    n_rows = graph.n_vertices
    indptr = np.zeros(n_rows + 1)
    for i in range(n_rows):
        inds, = np.where(rows == i)
        if inds.size == 0:
            indptr[i + 1] = indptr[i]
        else:
            indptr[i] = inds[0]
            indptr[i + 1] = inds[-1] + 1

    # create block sparse matrix
    if return_covariances:
        return (bsr_matrix((all_blocks, columns, indptr),
                           shape=(n_features, n_features), dtype=dtype),
                all_covariances)
    else:
        return bsr_matrix((all_blocks, columns, indptr),
                          shape=(n_features, n_features), dtype=dtype)
dataset/ETHPy150Open menpo/menpo/menpo/model/gmrf.py/_create_sparse_precision/operator `*` -> `-`