graphkit-learn/pygraph/kernels/randomWalkKernel.py

"""
@author: linlin
@references: S Vichy N Vishwanathan, Nicol N Schraudolph, Risi Kondor, and Karsten M Borgwardt. Graph kernels. Journal of Machine Learning Research, 11(Apr):1201–1242, 2010.
"""

import sys
sys.path.insert(0, "../")
import time
from functools import partial
from tqdm import tqdm

import networkx as nx
import numpy as np
from scipy.sparse import identity, kron
from scipy.sparse.linalg import cg
from scipy.optimize import fixed_point

from pygraph.utils.graphdataset import get_dataset_attributes
from pygraph.utils.parallel import parallel_gm

def randomwalkkernel(*args,
                     # params for all method.
                     compute_method=None,
                     weight=1, 
                     p=None, 
                     q=None,
                     edge_weight=None,
                     # params for conjugate and fp method.
                     node_kernels=None, 
                     edge_kernels=None,
                     node_label='atom',
                     edge_label='bond_type',
                     # params for spectral method.
                     sub_kernel=None,                                          
                     n_jobs=None,
                     verbose=True):
    """Calculate random walk graph kernels.
    Parameters
    ----------
    Gn : List of NetworkX graph
        List of graphs between which the kernels are calculated.
    /
    G1, G2 : NetworkX graphs
        Two graphs between which the kernel is calculated.
    compute_method : string
        Method used to compute kernel. The Following choices are 
        available:
        'sylvester' - Sylvester equation method.
        'conjugate' - conjugate gradient method.
        'fp' - fixed-point iterations.
        'spectral' - spectral decomposition.
    weight : float
        A constant weight set for random walks of length h.
    p : None
        Initial probability distribution on the unlabeled direct product graph 
        of two graphs. It is set to be uniform over all vertices in the direct 
        product graph.
    q : None
        Stopping probability distribution on the unlabeled direct product graph 
        of two graphs. It is set to be uniform over all vertices in the direct 
        product graph.
    edge_weight: float
        Edge attribute name corresponding to the edge weight.
        
    node_kernels: dict
        A dictionary of kernel functions for nodes, including 3 items: 'symb' 
        for symbolic node labels, 'nsymb' for non-symbolic node labels, 'mix' 
        for both labels. The first 2 functions take two node labels as 
        parameters, and the 'mix' function takes 4 parameters, a symbolic and a
        non-symbolic label for each the two nodes. Each label is in form of 2-D
        dimension array (n_samples, n_features). Each function returns a number
        as the kernel value. Ignored when nodes are unlabeled. This argument
        is designated to conjugate gradient method and fixed-point iterations.
    edge_kernels: dict
        A dictionary of kernel functions for edges, including 3 items: 'symb' 
        for symbolic edge labels, 'nsymb' for non-symbolic edge labels, 'mix' 
        for both labels. The first 2 functions take two edge labels as 
        parameters, and the 'mix' function takes 4 parameters, a symbolic and a
        non-symbolic label for each the two edges. Each label is in form of 2-D
        dimension array (n_samples, n_features). Each function returns a number
        as the kernel value. Ignored when edges are unlabeled. This argument
        is designated to conjugate gradient method and fixed-point iterations.
    node_label: string
        Node attribute used as label. The default node label is atom. This 
        argument is designated to conjugate gradient method and fixed-point 
        iterations.
    edge_label : string
        Edge attribute used as label. The default edge label is bond_type. This 
        argument is designated to conjugate gradient method and fixed-point 
        iterations.
        
    sub_kernel: string
        Method used to compute walk kernel. The Following choices are 
        available:
        'exp' : method based on exponential serials.
        'geo' : method based on geometric serials.
        
    n_jobs: int
        Number of jobs for parallelization. 

    Return
    ------
    Kmatrix : Numpy matrix
        Kernel matrix, each element of which is the path kernel up to d between 2 praphs.
    """
    compute_method = compute_method.lower()
    Gn = args[0] if len(args) == 1 else [args[0], args[1]]
    Gn = [g.copy() for g in Gn]

    eweight = None
    if edge_weight == None:
        if verbose:
            print('\n None edge weight specified. Set all weight to 1.\n')
    else:
        try:
            some_weight = list(
                nx.get_edge_attributes(Gn[0], edge_weight).values())[0]
            if isinstance(some_weight, float) or isinstance(some_weight, int):
                eweight = edge_weight
            else:
                if verbose:
                    print('\n Edge weight with name %s is not float or integer. Set all weight to 1.\n'
                          % edge_weight)
        except:
            if verbose:
                print('\n Edge weight with name "%s" is not found in the edge attributes. Set all weight to 1.\n'
                      % edge_weight)

    ds_attrs = get_dataset_attributes(
        Gn,
        attr_names=['node_labeled', 'node_attr_dim', 'edge_labeled',
                    'edge_attr_dim', 'is_directed'],
        node_label=node_label,
        edge_label=edge_label)
    
    # remove graphs with no edges, as no walk can be found in their structures, 
    # so the weight matrix between such a graph and itself might be zero.
    len_gn = len(Gn)
    Gn = [(idx, G) for idx, G in enumerate(Gn) if nx.number_of_edges(G) != 0]
    idx = [G[0] for G in Gn]
    Gn = [G[1] for G in Gn]
    if len(Gn) != len_gn:
        if verbose:
            print('\n %d graphs are removed as they don\'t contain edges.\n' %
                  (len_gn - len(Gn)))

    start_time = time.time()
    
#    # get vertex and edge concatenated labels for each graph
#    label_list, d = getLabels(Gn, node_label, edge_label, ds_attrs['is_directed'])
#    gmf = filterGramMatrix(A_wave_list[0], label_list[0], ('C', '0', 'O'), ds_attrs['is_directed'])

    if compute_method == 'sylvester':
        if verbose:
            import warnings
            warnings.warn('All labels are ignored.')
        Kmatrix = _sylvester_equation(Gn, weight, p, q, eweight, n_jobs, verbose=verbose)

    elif compute_method == 'conjugate':
        Kmatrix = _conjugate_gradient(Gn, weight, p, q, ds_attrs, node_kernels, 
                                      edge_kernels, node_label, edge_label, 
                                      eweight, n_jobs, verbose=verbose)
        
    elif compute_method == 'fp':
        Kmatrix = _fixed_point(Gn, weight, p, q, ds_attrs, node_kernels, 
                               edge_kernels, node_label, edge_label, 
                               eweight, n_jobs, verbose=verbose)

    elif compute_method == 'spectral':
        if verbose:
            import warnings
            warnings.warn('All labels are ignored. Only works for undirected graphs.')
        Kmatrix = _spectral_decomposition(Gn, weight, p, q, sub_kernel, 
                                          eweight, n_jobs, verbose=verbose)

    elif compute_method == 'kron':
        pass
        for i in range(0, len(Gn)):
            for j in range(i, len(Gn)):
                Kmatrix[i][j] = _randomwalkkernel_kron(Gn[i], Gn[j],
                                                       node_label, edge_label)
                Kmatrix[j][i] = Kmatrix[i][j]
    else:
        raise Exception(
            'compute method name incorrect. Available methods: "sylvester", "conjugate", "fp", "spectral" and "kron".'
        )

    run_time = time.time() - start_time
    if verbose:
        print("\n --- kernel matrix of random walk kernel of size %d built in %s seconds ---"
              % (len(Gn), run_time))

    return Kmatrix, run_time, idx


###############################################################################
def _sylvester_equation(Gn, lmda, p, q, eweight, n_jobs, verbose=True):
    """Calculate walk graph kernels up to n between 2 graphs using Sylvester method.

    Parameters
    ----------
    G1, G2 : NetworkX graph
        Graphs between which the kernel is calculated.
    node_label : string
        node attribute used as label.
    edge_label : string
        edge attribute used as label.

    Return
    ------
    kernel : float
        Kernel between 2 graphs.
    """
    Kmatrix = np.zeros((len(Gn), len(Gn)))

    if q == None:
        # don't normalize adjacency matrices if q is a uniform vector. Note
        # A_wave_list actually contains the transposes of the adjacency matrices.
        A_wave_list = [
            nx.adjacency_matrix(G, eweight).todense().transpose() for G in 
            (tqdm(Gn, desc='compute adjacency matrices', file=sys.stdout) if
             verbose else Gn)
        ]
#        # normalized adjacency matrices
#        A_wave_list = []
#        for G in tqdm(Gn, desc='compute adjacency matrices', file=sys.stdout):
#            A_tilde = nx.adjacency_matrix(G, eweight).todense().transpose()   
#            norm = A_tilde.sum(axis=0)
#            norm[norm == 0] = 1
#            A_wave_list.append(A_tilde / norm)
        if p == None: # p is uniform distribution as default.
            def init_worker(Awl_toshare):
                global G_Awl
                G_Awl = Awl_toshare
            do_partial = partial(wrapper_se_do, lmda)   
            parallel_gm(do_partial, Kmatrix, Gn, init_worker=init_worker, 
                        glbv=(A_wave_list,), n_jobs=n_jobs, verbose=verbose)
            
#            pbar = tqdm(
#                total=(1 + len(Gn)) * len(Gn) / 2,
#                desc='calculating kernels',
#                file=sys.stdout)
#            for i in range(0, len(Gn)):
#                for j in range(i, len(Gn)):
#                    S = lmda * A_wave_list[j]
#                    T_t = A_wave_list[i]
#                    # use uniform distribution if there is no prior knowledge.
#                    nb_pd = len(A_wave_list[i]) * len(A_wave_list[j])
#                    p_times_uni = 1 / nb_pd
#                    M0 = np.full((len(A_wave_list[j]), len(A_wave_list[i])), p_times_uni)
#                    X = dlyap(S, T_t, M0)
#                    X = np.reshape(X, (-1, 1), order='F')
#                    # use uniform distribution if there is no prior knowledge.
#                    q_times = np.full((1, nb_pd), p_times_uni)
#                    Kmatrix[i][j] = np.dot(q_times, X)
#                    Kmatrix[j][i] = Kmatrix[i][j]
#                    pbar.update(1)

    return Kmatrix


def wrapper_se_do(lmda, itr):
    i = itr[0]
    j = itr[1]
    return i, j, _se_do(G_Awl[i], G_Awl[j], lmda)


def _se_do(A_wave1, A_wave2, lmda):
    from control import dlyap
    S = lmda * A_wave2
    T_t = A_wave1
    # use uniform distribution if there is no prior knowledge.
    nb_pd = len(A_wave1) * len(A_wave2)
    p_times_uni = 1 / nb_pd
    M0 = np.full((len(A_wave2), len(A_wave1)), p_times_uni)
    X = dlyap(S, T_t, M0)
    X = np.reshape(X, (-1, 1), order='F')
    # use uniform distribution if there is no prior knowledge.
    q_times = np.full((1, nb_pd), p_times_uni)
    return np.dot(q_times, X)


###############################################################################
def _conjugate_gradient(Gn, lmda, p, q, ds_attrs, node_kernels, edge_kernels, 
                        node_label, edge_label, eweight, n_jobs, verbose=True):
    """Calculate walk graph kernels up to n between 2 graphs using conjugate method.

    Parameters
    ----------
    G1, G2 : NetworkX graph
        Graphs between which the kernel is calculated.
    node_label : string
        node attribute used as label.
    edge_label : string
        edge attribute used as label.

    Return
    ------
    kernel : float
        Kernel between 2 graphs.
    """
    Kmatrix = np.zeros((len(Gn), len(Gn)))
    
#    if not ds_attrs['node_labeled'] and ds_attrs['node_attr_dim'] < 1 and \
#        not ds_attrs['edge_labeled'] and ds_attrs['edge_attr_dim'] < 1:
#        # this is faster from unlabeled graphs. @todo: why?
#        if q == None:
#            # don't normalize adjacency matrices if q is a uniform vector. Note
#            # A_wave_list actually contains the transposes of the adjacency matrices.
#            A_wave_list = [
#                nx.adjacency_matrix(G, eweight).todense().transpose() for G in 
#                    tqdm(Gn, desc='compute adjacency matrices', file=sys.stdout)
#            ]
#            if p == None: # p is uniform distribution as default.
#                def init_worker(Awl_toshare):
#                    global G_Awl
#                    G_Awl = Awl_toshare
#                do_partial = partial(wrapper_cg_unlabled_do, lmda)   
#                parallel_gm(do_partial, Kmatrix, Gn, init_worker=init_worker, 
#                            glbv=(A_wave_list,), n_jobs=n_jobs)
#    else:  
    # reindex nodes using consecutive integers for convenience of kernel calculation.
    Gn = [nx.convert_node_labels_to_integers(
            g, first_label=0, label_attribute='label_orignal') for g in (tqdm(
                Gn, desc='reindex vertices', file=sys.stdout) if verbose else Gn)]
    
    if p == None and q == None: # p and q are uniform distributions as default.
        def init_worker(gn_toshare):
            global G_gn
            G_gn = gn_toshare
        do_partial = partial(wrapper_cg_labled_do, ds_attrs, node_kernels, 
                             node_label, edge_kernels, edge_label, lmda)   
        parallel_gm(do_partial, Kmatrix, Gn, init_worker=init_worker, 
                    glbv=(Gn,), n_jobs=n_jobs, verbose=verbose)  
            
#            pbar = tqdm(
#                total=(1 + len(Gn)) * len(Gn) / 2,
#                desc='calculating kernels',
#                file=sys.stdout)
#            for i in range(0, len(Gn)):
#                for j in range(i, len(Gn)):
#                    result = _cg_labled_do(Gn[i], Gn[j], ds_attrs, node_kernels,
#                                           node_label, edge_kernels, edge_label, lmda)
#                    Kmatrix[i][j] = result
#                    Kmatrix[j][i] = Kmatrix[i][j]
#                    pbar.update(1)
    return Kmatrix


def wrapper_cg_unlabled_do(lmda, itr):
    i = itr[0]
    j = itr[1]
    return i, j, _cg_unlabled_do(G_Awl[i], G_Awl[j], lmda)


def _cg_unlabled_do(A_wave1, A_wave2, lmda):
    nb_pd = len(A_wave1) * len(A_wave2)
    p_times_uni = 1 / nb_pd
    w_times = kron(A_wave1, A_wave2).todense()
    A = identity(w_times.shape[0]) - w_times * lmda
    b = np.full((nb_pd, 1), p_times_uni)
    x, _ = cg(A, b)
    # use uniform distribution if there is no prior knowledge.
    q_times = np.full((1, nb_pd), p_times_uni)
    return np.dot(q_times, x)


def wrapper_cg_labled_do(ds_attrs, node_kernels, node_label, edge_kernels, 
                         edge_label, lmda, itr):
    i = itr[0]
    j = itr[1]
    return i, j, _cg_labled_do(G_gn[i], G_gn[j], ds_attrs, node_kernels, 
                               node_label, edge_kernels, edge_label, lmda)


def _cg_labled_do(g1, g2, ds_attrs, node_kernels, node_label, 
                  edge_kernels, edge_label, lmda):
    # Frist, compute kernels between all pairs of nodes, method borrowed
    # from FCSP. It is faster than directly computing all edge kernels 
    # when $d_1d_2>2$, where $d_1$ and $d_2$ are vertex degrees of the
    # graphs compared, which is the most case we went though. For very 
    # sparse graphs, this would be slow.
    vk_dict = computeVK(g1, g2, ds_attrs, node_kernels, node_label)
                           
    # Compute weight matrix of the direct product graph.   
    w_times, w_dim = computeW(g1, g2, vk_dict, ds_attrs,
                              edge_kernels, edge_label)                                                            
    # use uniform distribution if there is no prior knowledge.
    p_times_uni = 1 / w_dim
    A = identity(w_times.shape[0]) - w_times * lmda
    b = np.full((w_dim, 1), p_times_uni)
    x, _ = cg(A, b)
    # use uniform distribution if there is no prior knowledge.
    q_times = np.full((1, w_dim), p_times_uni)
    return np.dot(q_times, x)


###############################################################################
def _fixed_point(Gn, lmda, p, q, ds_attrs, node_kernels, edge_kernels, 
                         node_label, edge_label, eweight, n_jobs, verbose=True):
    """Calculate walk graph kernels up to n between 2 graphs using Fixed-Point method.

    Parameters
    ----------
    G1, G2 : NetworkX graph
        Graphs between which the kernel is calculated.
    node_label : string
        node attribute used as label.
    edge_label : string
        edge attribute used as label.

    Return
    ------
    kernel : float
        Kernel between 2 graphs.
    """
    

    Kmatrix = np.zeros((len(Gn), len(Gn)))
    
#    if not ds_attrs['node_labeled'] and ds_attrs['node_attr_dim'] < 1 and \
#        not ds_attrs['edge_labeled'] and ds_attrs['edge_attr_dim'] > 1:
#        # this is faster from unlabeled graphs. @todo: why?
#        if q == None:
#            # don't normalize adjacency matrices if q is a uniform vector. Note
#            # A_wave_list actually contains the transposes of the adjacency matrices.
#            A_wave_list = [
#                nx.adjacency_matrix(G, eweight).todense().transpose() for G in 
#                    tqdm(Gn, desc='compute adjacency matrices', file=sys.stdout)
#            ]
#            if p == None: # p is uniform distribution as default.
#                pbar = tqdm(
#                    total=(1 + len(Gn)) * len(Gn) / 2,
#                    desc='calculating kernels',
#                    file=sys.stdout)
#                for i in range(0, len(Gn)):
#                    for j in range(i, len(Gn)):                   
#                        # use uniform distribution if there is no prior knowledge.
#                        nb_pd = len(A_wave_list[i]) * len(A_wave_list[j])
#                        p_times_uni = 1 / nb_pd
#                        w_times = kron(A_wave_list[i], A_wave_list[j]).todense()
#                        p_times = np.full((nb_pd, 1), p_times_uni)
#                        x = fixed_point(func_fp, p_times, args=(p_times, lmda, w_times))
#                        # use uniform distribution if there is no prior knowledge.
#                        q_times = np.full((1, nb_pd), p_times_uni)
#                        Kmatrix[i][j] = np.dot(q_times, x)
#                        Kmatrix[j][i] = Kmatrix[i][j]
#                        pbar.update(1)
#    else:  
    # reindex nodes using consecutive integers for convenience of kernel calculation.
    Gn = [nx.convert_node_labels_to_integers(
            g, first_label=0, label_attribute='label_orignal') for g in (tqdm(
                Gn, desc='reindex vertices', file=sys.stdout) if verbose else Gn)]
    
    if p == None and q == None: # p and q are uniform distributions as default.
        def init_worker(gn_toshare):
            global G_gn
            G_gn = gn_toshare
        do_partial = partial(wrapper_fp_labled_do, ds_attrs, node_kernels, 
                             node_label, edge_kernels, edge_label, lmda)   
        parallel_gm(do_partial, Kmatrix, Gn, init_worker=init_worker, 
                    glbv=(Gn,), n_jobs=n_jobs, verbose=verbose)
    return Kmatrix


def wrapper_fp_labled_do(ds_attrs, node_kernels, node_label, edge_kernels, 
                         edge_label, lmda, itr):
    i = itr[0]
    j = itr[1]
    return i, j, _fp_labled_do(G_gn[i], G_gn[j], ds_attrs, node_kernels, 
                               node_label, edge_kernels, edge_label, lmda)


def _fp_labled_do(g1, g2, ds_attrs, node_kernels, node_label, 
                  edge_kernels, edge_label, lmda):
    # Frist, compute kernels between all pairs of nodes, method borrowed
    # from FCSP. It is faster than directly computing all edge kernels 
    # when $d_1d_2>2$, where $d_1$ and $d_2$ are vertex degrees of the
    # graphs compared, which is the most case we went though. For very 
    # sparse graphs, this would be slow.
    vk_dict = computeVK(g1, g2, ds_attrs, node_kernels, node_label)
                           
    # Compute weight matrix of the direct product graph.   
    w_times, w_dim = computeW(g1, g2, vk_dict, ds_attrs,
                              edge_kernels, edge_label)                                                            
    # use uniform distribution if there is no prior knowledge.
    p_times_uni = 1 / w_dim
    p_times = np.full((w_dim, 1), p_times_uni)
    x = fixed_point(func_fp, p_times, args=(p_times, lmda, w_times),
                    xtol=1e-06, maxiter=1000)
    # use uniform distribution if there is no prior knowledge.
    q_times = np.full((1, w_dim), p_times_uni)
    return np.dot(q_times, x)


def func_fp(x, p_times, lmda, w_times):
    haha = w_times * x
    haha = lmda * haha
    haha = p_times + haha
    return p_times + lmda * np.dot(w_times, x)


###############################################################################
def _spectral_decomposition(Gn, weight, p, q, sub_kernel, eweight, n_jobs, verbose=True):
    """Calculate walk graph kernels up to n between 2 unlabeled graphs using 
    spectral decomposition method. Labels will be ignored.

    Parameters
    ----------
    G1, G2 : NetworkX graph
        Graphs between which the kernel is calculated.
    node_label : string
        node attribute used as label.
    edge_label : string
        edge attribute used as label.

    Return
    ------
    kernel : float
        Kernel between 2 graphs.
    """
    Kmatrix = np.zeros((len(Gn), len(Gn)))

    if q == None:
        # precompute the spectral decomposition of each graph.
        P_list = []
        D_list = []
        for G in (tqdm(Gn, desc='spectral decompose', file=sys.stdout) if 
                  verbose else Gn):
            # don't normalize adjacency matrices if q is a uniform vector. Note
            # A actually is the transpose of the adjacency matrix.
            A = nx.adjacency_matrix(G, eweight).todense().transpose()
            ew, ev = np.linalg.eig(A)
            D_list.append(ew)
            P_list.append(ev)
#        P_inv_list = [p.T for p in P_list] # @todo: also works for directed graphs?

        if p == None: # p is uniform distribution as default.            
            q_T_list = [np.full((1, nx.number_of_nodes(G)), 1 / nx.number_of_nodes(G)) for G in Gn]
#            q_T_list = [q.T for q in q_list]
            def init_worker(q_T_toshare, P_toshare, D_toshare):
                global G_q_T, G_P, G_D 
                G_q_T = q_T_toshare
                G_P = P_toshare
                G_D = D_toshare
            do_partial = partial(wrapper_sd_do, weight, sub_kernel)   
            parallel_gm(do_partial, Kmatrix, Gn, init_worker=init_worker, 
                        glbv=(q_T_list, P_list, D_list), n_jobs=n_jobs, 
                        verbose=verbose)
            
            
#            pbar = tqdm(
#                total=(1 + len(Gn)) * len(Gn) / 2,
#                desc='calculating kernels',
#                file=sys.stdout)
#            for i in range(0, len(Gn)):
#                for j in range(i, len(Gn)):
#                    result = _sd_do(q_T_list[i], q_T_list[j], P_list[i], P_list[j], 
#                                    D_list[i], D_list[j], weight, sub_kernel)
#                    Kmatrix[i][j] = result
#                    Kmatrix[j][i] = Kmatrix[i][j]
#                    pbar.update(1)
    return Kmatrix


def wrapper_sd_do(weight, sub_kernel, itr):
    i = itr[0]
    j = itr[1]
    return i, j, _sd_do(G_q_T[i], G_q_T[j], G_P[i], G_P[j], G_D[i], G_D[j], 
                        weight, sub_kernel)


def _sd_do(q_T1, q_T2, P1, P2, D1, D2, weight, sub_kernel):    
    # use uniform distribution if there is no prior knowledge.
    kl = kron(np.dot(q_T1, P1), np.dot(q_T2, P2)).todense()
    # @todo: this is not be needed when p = q (kr = kl.T) for undirected graphs
#    kr = kron(np.dot(P_inv_list[i], q_list[i]), np.dot(P_inv_list[j], q_list[j])).todense()
    if sub_kernel == 'exp':
        D_diag = np.array([d1 * d2 for d1 in D1 for d2 in D2])
        kmiddle = np.diag(np.exp(weight * D_diag))
    elif sub_kernel == 'geo':
        D_diag = np.array([d1 * d2 for d1 in D1 for d2 in D2])
        kmiddle = np.diag(weight * D_diag)
        kmiddle = np.identity(len(kmiddle)) - weight * kmiddle
        kmiddle = np.linalg.inv(kmiddle)
    return np.dot(np.dot(kl, kmiddle), kl.T)[0, 0]


###############################################################################
def _randomwalkkernel_kron(G1, G2, node_label, edge_label):
    """Calculate walk graph kernels up to n between 2 graphs using nearest Kronecker product approximation method.

    Parameters
    ----------
    G1, G2 : NetworkX graph
        Graphs between which the kernel is calculated.
    node_label : string
        node attribute used as label.
    edge_label : string
        edge attribute used as label.

    Return
    ------
    kernel : float
        Kernel between 2 graphs.
    """
    pass


###############################################################################
def getLabels(Gn, node_label, edge_label, directed):
    """Get symbolic labels of a graph dataset, where vertex labels are dealt
    with by concatenating them to the edge labels of adjacent edges.
    """
    label_list = []
    label_set = set()
    for g in Gn:
        label_g = {}
        for e in g.edges(data=True):
            nl1 = g.node[e[0]][node_label]
            nl2 = g.node[e[1]][node_label]
            if not directed and nl1 > nl2:
                nl1, nl2 = nl2, nl1
            label = (nl1, e[2][edge_label], nl2)
            label_g[(e[0], e[1])] = label
        label_list.append(label_g)  
    label_set = set([l for lg in label_list for l in lg.values()])
    return label_list, len(label_set)


def filterGramMatrix(gmt, label_dict, label, directed):
    """Compute (the transpose of) the Gram matrix filtered by a label.
    """
    gmf = np.zeros(gmt.shape)
    for (n1, n2), l in label_dict.items():
        if l == label:
            gmf[n2, n1] = gmt[n2, n1]
            if not directed:
                gmf[n1, n2] = gmt[n1, n2]
    return gmf


def computeVK(g1, g2, ds_attrs, node_kernels, node_label):
    '''Compute vertex kernels between vertices of two graphs.
    '''
    vk_dict = {}  # shortest path matrices dict
    if ds_attrs['node_labeled']:
        # node symb and non-synb labeled
        if ds_attrs['node_attr_dim'] > 0:
            kn = node_kernels['mix']
            for n1 in g1.nodes(data=True):
                for n2 in g2.nodes(data=True):
                    vk_dict[(n1[0], n2[0])] = kn(
                        n1[1][node_label], n2[1][node_label],
                        n1[1]['attributes'], n2[1]['attributes'])
        # node symb labeled
        else:
            kn = node_kernels['symb']
            for n1 in g1.nodes(data=True):
                for n2 in g2.nodes(data=True):
                    vk_dict[(n1[0], n2[0])] = kn(n1[1][node_label],
                                                 n2[1][node_label])
    else:
        # node non-synb labeled
        if ds_attrs['node_attr_dim'] > 0:
            kn = node_kernels['nsymb']
            for n1 in g1.nodes(data=True):
                for n2 in g2.nodes(data=True):
                    vk_dict[(n1[0], n2[0])] = kn(n1[1]['attributes'],
                                                 n2[1]['attributes'])
        # node unlabeled
        else:
            pass
    return vk_dict


def computeW(g1, g2, vk_dict, ds_attrs, edge_kernels, edge_label):
    '''Compute weight matrix of the direct product graph.
    '''
    w_dim = nx.number_of_nodes(g1) * nx.number_of_nodes(g2)
    w_times = np.zeros((w_dim, w_dim))
    if vk_dict: # node labeled
        if ds_attrs['is_directed']:
            if ds_attrs['edge_labeled']:
                # edge symb and non-synb labeled
                if ds_attrs['edge_attr_dim'] > 0:
                    ke = edge_kernels['mix']
                    for e1 in g1.edges(data=True):
                        for e2 in g2.edges(data=True):
                            ek_temp = ke(e1[2][edge_label], e2[2][edge_label],
                                         e1[2]['attributes'], e2[2]['attributes'])
                            w_idx = (e1[0] * nx.number_of_nodes(g2) + e2[0],
                                     e1[1] * nx.number_of_nodes(g2) + e2[1])
                            w_times[w_idx] = vk_dict[(e1[0], e2[0])] \
                                * ek_temp * vk_dict[(e1[1], e2[1])]
                # edge symb labeled
                else:
                    ke = edge_kernels['symb']
                    for e1 in g1.edges(data=True):
                        for e2 in g2.edges(data=True):
                            ek_temp = ke(e1[2][edge_label], e2[2][edge_label])
                            w_idx = (e1[0] * nx.number_of_nodes(g2) + e2[0],
                                     e1[1] * nx.number_of_nodes(g2) + e2[1])
                            w_times[w_idx] = vk_dict[(e1[0], e2[0])] \
                                * ek_temp * vk_dict[(e1[1], e2[1])]
            else:
                # edge non-synb labeled
                if ds_attrs['edge_attr_dim'] > 0:
                    ke = edge_kernels['nsymb']
                    for e1 in g1.edges(data=True):
                        for e2 in g2.edges(data=True):
                            ek_temp = ke(e1[2]['attributes'], e2[2]['attributes'])
                            w_idx = (e1[0] * nx.number_of_nodes(g2) + e2[0],
                                     e1[1] * nx.number_of_nodes(g2) + e2[1])
                            w_times[w_idx] = vk_dict[(e1[0], e2[0])] \
                                * ek_temp * vk_dict[(e1[1], e2[1])]
                # edge unlabeled
                else:
                    for e1 in g1.edges(data=True):
                        for e2 in g2.edges(data=True):
                            w_idx = (e1[0] * nx.number_of_nodes(g2) + e2[0],
                                     e1[1] * nx.number_of_nodes(g2) + e2[1])
                            w_times[w_idx] = vk_dict[(e1[0], e2[0])] \
                                * vk_dict[(e1[1], e2[1])]                                
        else: # undirected
            if ds_attrs['edge_labeled']:
                # edge symb and non-synb labeled
                if ds_attrs['edge_attr_dim'] > 0:
                    ke = edge_kernels['mix']
                    for e1 in g1.edges(data=True):
                        for e2 in g2.edges(data=True):
                            ek_temp = ke(e1[2][edge_label], e2[2][edge_label],
                                         e1[2]['attributes'], e2[2]['attributes'])
                            w_idx = (e1[0] * nx.number_of_nodes(g2) + e2[0],
                                     e1[1] * nx.number_of_nodes(g2) + e2[1])
                            w_times[w_idx] = vk_dict[(e1[0], e2[0])] \
                                * ek_temp * vk_dict[(e1[1], e2[1])] \
                                + vk_dict[(e1[0], e2[1])] \
                                * ek_temp * vk_dict[(e1[1], e2[0])]
                            w_times[w_idx[1], w_idx[0]] = w_times[w_idx[0], w_idx[1]]
                            w_idx2 = (e1[0] * nx.number_of_nodes(g2) + e2[1],
                                     e1[1] * nx.number_of_nodes(g2) + e2[0])
                            w_times[w_idx2[0], w_idx2[1]] = w_times[w_idx[0], w_idx[1]]
                            w_times[w_idx2[1], w_idx2[0]] = w_times[w_idx[0], w_idx[1]]
                # edge symb labeled
                else:
                    ke = edge_kernels['symb']
                    for e1 in g1.edges(data=True):
                        for e2 in g2.edges(data=True):
                            ek_temp = ke(e1[2][edge_label], e2[2][edge_label])
                            w_idx = (e1[0] * nx.number_of_nodes(g2) + e2[0],
                                     e1[1] * nx.number_of_nodes(g2) + e2[1])
                            w_times[w_idx] = vk_dict[(e1[0], e2[0])] \
                                * ek_temp * vk_dict[(e1[1], e2[1])] \
                                + vk_dict[(e1[0], e2[1])] \
                                * ek_temp * vk_dict[(e1[1], e2[0])]
                            w_times[w_idx[1], w_idx[0]] = w_times[w_idx[0], w_idx[1]]
                            w_idx2 = (e1[0] * nx.number_of_nodes(g2) + e2[1],
                                     e1[1] * nx.number_of_nodes(g2) + e2[0])
                            w_times[w_idx2[0], w_idx2[1]] = w_times[w_idx[0], w_idx[1]]
                            w_times[w_idx2[1], w_idx2[0]] = w_times[w_idx[0], w_idx[1]]
            else:
                # edge non-synb labeled
                if ds_attrs['edge_attr_dim'] > 0:
                    ke = edge_kernels['nsymb']
                    for e1 in g1.edges(data=True):
                        for e2 in g2.edges(data=True):
                            ek_temp = ke(e1[2]['attributes'], e2[2]['attributes'])
                            w_idx = (e1[0] * nx.number_of_nodes(g2) + e2[0],
                                     e1[1] * nx.number_of_nodes(g2) + e2[1])
                            w_times[w_idx] = vk_dict[(e1[0], e2[0])] \
                                * ek_temp * vk_dict[(e1[1], e2[1])] \
                                + vk_dict[(e1[0], e2[1])] \
                                * ek_temp * vk_dict[(e1[1], e2[0])]
                            w_times[w_idx[1], w_idx[0]] = w_times[w_idx[0], w_idx[1]]
                            w_idx2 = (e1[0] * nx.number_of_nodes(g2) + e2[1],
                                     e1[1] * nx.number_of_nodes(g2) + e2[0])
                            w_times[w_idx2[0], w_idx2[1]] = w_times[w_idx[0], w_idx[1]]
                            w_times[w_idx2[1], w_idx2[0]] = w_times[w_idx[0], w_idx[1]]
                # edge unlabeled
                else:
                    for e1 in g1.edges(data=True):
                        for e2 in g2.edges(data=True):
                            w_idx = (e1[0] * nx.number_of_nodes(g2) + e2[0],
                                     e1[1] * nx.number_of_nodes(g2) + e2[1])
                            w_times[w_idx] = vk_dict[(e1[0], e2[0])] \
                                * vk_dict[(e1[1], e2[1])] \
                                + vk_dict[(e1[0], e2[1])] \
                                * vk_dict[(e1[1], e2[0])]
                            w_times[w_idx[1], w_idx[0]] = w_times[w_idx[0], w_idx[1]]
                            w_idx2 = (e1[0] * nx.number_of_nodes(g2) + e2[1],
                                     e1[1] * nx.number_of_nodes(g2) + e2[0])
                            w_times[w_idx2[0], w_idx2[1]] = w_times[w_idx[0], w_idx[1]]
                            w_times[w_idx2[1], w_idx2[0]] = w_times[w_idx[0], w_idx[1]]
    else: # node unlabeled
        if ds_attrs['is_directed']:
            if ds_attrs['edge_labeled']:
                # edge symb and non-synb labeled
                if ds_attrs['edge_attr_dim'] > 0:
                    ke = edge_kernels['mix']
                    for e1 in g1.edges(data=True):
                        for e2 in g2.edges(data=True):
                            ek_temp = ke(e1[2][edge_label], e2[2][edge_label],
                                         e1[2]['attributes'], e2[2]['attributes'])
                            w_idx = (e1[0] * nx.number_of_nodes(g2) + e2[0],
                                     e1[1] * nx.number_of_nodes(g2) + e2[1])
                            w_times[w_idx] = ek_temp
                # edge symb labeled
                else:
                    ke = edge_kernels['symb']
                    for e1 in g1.edges(data=True):
                        for e2 in g2.edges(data=True):
                            ek_temp = ke(e1[2][edge_label], e2[2][edge_label])
                            w_idx = (e1[0] * nx.number_of_nodes(g2) + e2[0],
                                     e1[1] * nx.number_of_nodes(g2) + e2[1])
                            w_times[w_idx] = ek_temp
            else:
                # edge non-synb labeled
                if ds_attrs['edge_attr_dim'] > 0:
                    ke = edge_kernels['nsymb']
                    for e1 in g1.edges(data=True):
                        for e2 in g2.edges(data=True):
                            ek_temp = ke(e1[2]['attributes'], e2[2]['attributes'])
                            w_idx = (e1[0] * nx.number_of_nodes(g2) + e2[0],
                                     e1[1] * nx.number_of_nodes(g2) + e2[1])
                            w_times[w_idx] = ek_temp
                # edge unlabeled
                else:
                    for e1 in g1.edges(data=True):
                        for e2 in g2.edges(data=True):
                            w_idx = (e1[0] * nx.number_of_nodes(g2) + e2[0],
                                     e1[1] * nx.number_of_nodes(g2) + e2[1])
                            w_times[w_idx] = 1                                
        else: # undirected
            if ds_attrs['edge_labeled']:
                # edge symb and non-synb labeled
                if ds_attrs['edge_attr_dim'] > 0:
                    ke = edge_kernels['mix']
                    for e1 in g1.edges(data=True):
                        for e2 in g2.edges(data=True):
                            ek_temp = ke(e1[2][edge_label], e2[2][edge_label],
                                         e1[2]['attributes'], e2[2]['attributes'])
                            w_idx = (e1[0] * nx.number_of_nodes(g2) + e2[0],
                                     e1[1] * nx.number_of_nodes(g2) + e2[1])
                            w_times[w_idx] = ek_temp
                            w_times[w_idx[1], w_idx[0]] = w_times[w_idx[0], w_idx[1]]
                            w_idx2 = (e1[0] * nx.number_of_nodes(g2) + e2[1],
                                     e1[1] * nx.number_of_nodes(g2) + e2[0])
                            w_times[w_idx2[0], w_idx2[1]] = w_times[w_idx[0], w_idx[1]]
                            w_times[w_idx2[1], w_idx2[0]] = w_times[w_idx[0], w_idx[1]]
                # edge symb labeled
                else:
                    ke = edge_kernels['symb']
                    for e1 in g1.edges(data=True):
                        for e2 in g2.edges(data=True):
                            ek_temp = ke(e1[2][edge_label], e2[2][edge_label])
                            w_idx = (e1[0] * nx.number_of_nodes(g2) + e2[0],
                                     e1[1] * nx.number_of_nodes(g2) + e2[1])
                            w_times[w_idx] = ek_temp
                            w_times[w_idx[1], w_idx[0]] = w_times[w_idx[0], w_idx[1]]
                            w_idx2 = (e1[0] * nx.number_of_nodes(g2) + e2[1],
                                     e1[1] * nx.number_of_nodes(g2) + e2[0])
                            w_times[w_idx2[0], w_idx2[1]] = w_times[w_idx[0], w_idx[1]]
                            w_times[w_idx2[1], w_idx2[0]] = w_times[w_idx[0], w_idx[1]]
            else:
                # edge non-synb labeled
                if ds_attrs['edge_attr_dim'] > 0:
                    ke = edge_kernels['nsymb']
                    for e1 in g1.edges(data=True):
                        for e2 in g2.edges(data=True):
                            ek_temp = ke(e1[2]['attributes'], e2[2]['attributes'])
                            w_idx = (e1[0] * nx.number_of_nodes(g2) + e2[0],
                                     e1[1] * nx.number_of_nodes(g2) + e2[1])
                            w_times[w_idx] = ek_temp
                            w_times[w_idx[1], w_idx[0]] = w_times[w_idx[0], w_idx[1]]
                            w_idx2 = (e1[0] * nx.number_of_nodes(g2) + e2[1],
                                     e1[1] * nx.number_of_nodes(g2) + e2[0])
                            w_times[w_idx2[0], w_idx2[1]] = w_times[w_idx[0], w_idx[1]]
                            w_times[w_idx2[1], w_idx2[0]] = w_times[w_idx[0], w_idx[1]]
                # edge unlabeled
                else:
                    for e1 in g1.edges(data=True):
                        for e2 in g2.edges(data=True):
                            w_idx = (e1[0] * nx.number_of_nodes(g2) + e2[0],
                                     e1[1] * nx.number_of_nodes(g2) + e2[1])
                            w_times[w_idx] = 1
                            w_times[w_idx[1], w_idx[0]] = w_times[w_idx[0], w_idx[1]]
                            w_idx2 = (e1[0] * nx.number_of_nodes(g2) + e2[1],
                                     e1[1] * nx.number_of_nodes(g2) + e2[0])
                            w_times[w_idx2[0], w_idx2[1]] = w_times[w_idx[0], w_idx[1]]
                            w_times[w_idx2[1], w_idx2[0]] = w_times[w_idx[0], w_idx[1]]
    return w_times, w_dim