compute_gram_matrices() (in module pygraph.utils.model_selection_precomputed)

model_selection_for_precomputed_kernel() (in module pygraph.utils.model_selection_precomputed)

parallel_trial_do() (in module pygraph.utils.model_selection_precomputed)

printResultsInTable() (in module pygraph.utils.model_selection_precomputed)

read_gram_matrices_from_file() (in module pygraph.utils.model_selection_precomputed)

trial_do() (in module pygraph.utils.model_selection_precomputed)

""" Obtain all kinds of attributes of a graph dataset.
"""

def get_dataset_attributes(Gn,
                           target=None,
                           attr_names=[],
                           node_label=None,
                           edge_label=None):
    """Returns the structure and property information of the graph dataset Gn.

    Parameters
    ----------
    Gn : List of NetworkX graph
        List of graphs whose information will be returned.
    target : list
        The list of classification targets corresponding to Gn. Only works for 
        classification problems.
    attr_names : list
        List of strings which indicate which informations will be returned. The
        possible choices includes:
        'substructures': sub-structures Gn contains, including 'linear', 'non 
            linear' and 'cyclic'.
        'node_labeled': whether vertices have symbolic labels.
        'edge_labeled': whether egdes have symbolic labels.
        'is_directed': whether graphs in Gn are directed.
        'dataset_size': number of graphs in Gn.
        'ave_node_num': average number of vertices of graphs in Gn.
        'min_node_num': minimum number of vertices of graphs in Gn.
        'max_node_num': maximum number of vertices of graphs in Gn.
        'ave_edge_num': average number of edges of graphs in Gn.
        'min_edge_num': minimum number of edges of graphs in Gn.
        'max_edge_num': maximum number of edges of graphs in Gn.
        'ave_node_degree': average vertex degree of graphs in Gn.
        'min_node_degree': minimum vertex degree of graphs in Gn.
        'max_node_degree': maximum vertex degree of graphs in Gn.
        'ave_fill_factor': average fill factor (number_of_edges / 
            (number_of_nodes ** 2)) of graphs in Gn.
        'min_fill_factor': minimum fill factor of graphs in Gn.
        'max_fill_factor': maximum fill factor of graphs in Gn.
        'node_label_num': number of symbolic vertex labels.
        'edge_label_num': number of symbolic edge labels.
        'node_attr_dim': number of dimensions of non-symbolic vertex labels. 
            Extracted from the 'attributes' attribute of graph nodes.
        'edge_attr_dim': number of dimensions of non-symbolic edge labels.
            Extracted from the 'attributes' attribute of graph edges.
        'class_number': number of classes. Only available for classification 
            problems.
    node_label : string
        Node attribute used as label. The default node label is atom. Mandatory
        when 'node_labeled' or 'node_label_num' is required.
    edge_label : string
        Edge attribute used as label. The default edge label is bond_type. 
        Mandatory when 'edge_labeled' or 'edge_label_num' is required.

    Return
    ------
    attrs : dict
        Value for each property.
    """
    import networkx as nx
    import numpy as np

    attrs = {}

    def get_dataset_size(Gn):
        return len(Gn)

    def get_all_node_num(Gn):
        return [nx.number_of_nodes(G) for G in Gn]

    def get_ave_node_num(all_node_num):
        return np.mean(all_node_num)

    def get_min_node_num(all_node_num):
        return np.amin(all_node_num)

    def get_max_node_num(all_node_num):
        return np.amax(all_node_num)

    def get_all_edge_num(Gn):
        return [nx.number_of_edges(G) for G in Gn]

    def get_ave_edge_num(all_edge_num):
        return np.mean(all_edge_num)

    def get_min_edge_num(all_edge_num):
        return np.amin(all_edge_num)

    def get_max_edge_num(all_edge_num):
        return np.amax(all_edge_num)

    def is_node_labeled(Gn):
        return False if node_label is None else True

    def get_node_label_num(Gn):
        nl = set()
        for G in Gn:
            nl = nl | set(nx.get_node_attributes(G, node_label).values())
        return len(nl)

    def is_edge_labeled(Gn):
        return False if edge_label is None else True

    def get_edge_label_num(Gn):
        el = set()
        for G in Gn:
            el = el | set(nx.get_edge_attributes(G, edge_label).values())
        return len(el)

    def is_directed(Gn):
        return nx.is_directed(Gn[0])

    def get_ave_node_degree(Gn):
        return np.mean([np.mean(list(dict(G.degree()).values())) for G in Gn])

    def get_max_node_degree(Gn):
        return np.amax([np.mean(list(dict(G.degree()).values())) for G in Gn])

    def get_min_node_degree(Gn):
        return np.amin([np.mean(list(dict(G.degree()).values())) for G in Gn])
    
    # get fill factor, the number of non-zero entries in the adjacency matrix.
    def get_ave_fill_factor(Gn):
        return np.mean([nx.number_of_edges(G) / (nx.number_of_nodes(G) 
                                           * nx.number_of_nodes(G)) for G in Gn])

    def get_max_fill_factor(Gn):
        return np.amax([nx.number_of_edges(G) / (nx.number_of_nodes(G) 

    def get_min_fill_factor(Gn):
        return np.amin([nx.number_of_edges(G) / (nx.number_of_nodes(G) 

    def get_substructures(Gn):
        subs = set()
        for G in Gn:
            degrees = list(dict(G.degree()).values())
            if any(i == 2 for i in degrees):
                subs.add('linear')
            if np.amax(degrees) >= 3:
                subs.add('non linear')
            if 'linear' in subs and 'non linear' in subs:
                break

        if is_directed(Gn):
            for G in Gn:
                if len(list(nx.find_cycle(G))) > 0:
                    subs.add('cyclic')
                    break
        # else:
        #     # @todo: this method does not work for big graph with large amount of edges like D&D, try a better way.
        #     upper = np.amin([nx.number_of_edges(G) for G in Gn]) * 2 + 10
        #     for G in Gn:
        #         if (nx.number_of_edges(G) < upper):
        #             cyc = list(nx.simple_cycles(G.to_directed()))
        #             if any(len(i) > 2 for i in cyc):
        #                 subs.add('cyclic')
        #                 break
        #     if 'cyclic' not in subs:
        #         for G in Gn:
        #                 break

        return subs

    def get_class_num(target):
        return len(set(target))

    def get_node_attr_dim(Gn):
        for G in Gn:
            for n in G.nodes(data=True):
                if 'attributes' in n[1]:
                    return len(n[1]['attributes'])
        return 0

    def get_edge_attr_dim(Gn):
        for G in Gn:
            if nx.number_of_edges(G) > 0:
                for e in G.edges(data=True):
                    if 'attributes' in e[2]:
                        return len(e[2]['attributes'])
        return 0

    if attr_names == []:
        attr_names = [
            'substructures',
            'node_labeled',
            'edge_labeled',
            'is_directed',
            'dataset_size',
            'ave_node_num',
            'min_node_num',
            'max_node_num',
            'ave_edge_num',
            'min_edge_num',
            'max_edge_num',
            'ave_node_degree',
            'min_node_degree',
            'max_node_degree',
            'ave_fill_factor',
            'min_fill_factor',
            'max_fill_factor',
            'node_label_num',
            'edge_label_num',
            'node_attr_dim',
            'edge_attr_dim',
            'class_number',
        ]

    # dataset size
    if 'dataset_size' in attr_names:

        attrs.update({'dataset_size': get_dataset_size(Gn)})

    # graph node number
    if any(i in attr_names
           for i in ['ave_node_num', 'min_node_num', 'max_node_num']):

        all_node_num = get_all_node_num(Gn)

    if 'ave_node_num' in attr_names:

        attrs.update({'ave_node_num': get_ave_node_num(all_node_num)})

    if 'min_node_num' in attr_names:

        attrs.update({'min_node_num': get_min_node_num(all_node_num)})

    if 'max_node_num' in attr_names:

        attrs.update({'max_node_num': get_max_node_num(all_node_num)})

    # graph edge number
    if any(i in attr_names for i in
           ['ave_edge_num', 'min_edge_num', 'max_edge_num']):

        all_edge_num = get_all_edge_num(Gn)

    if 'ave_edge_num' in attr_names:

        attrs.update({'ave_edge_num': get_ave_edge_num(all_edge_num)})

    if 'max_edge_num' in attr_names:

        attrs.update({'max_edge_num': get_max_edge_num(all_edge_num)})

    if 'min_edge_num' in attr_names:

        attrs.update({'min_edge_num': get_min_edge_num(all_edge_num)})

    # label number
    if any(i in attr_names for i in ['node_labeled', 'node_label_num']):
        is_nl = is_node_labeled(Gn)
        node_label_num = get_node_label_num(Gn)

    if 'node_labeled' in attr_names:
        # graphs are considered node unlabeled if all nodes have the same label.
        attrs.update({'node_labeled': is_nl if node_label_num > 1 else False})

    if 'node_label_num' in attr_names:
        attrs.update({'node_label_num': node_label_num})

    if any(i in attr_names for i in ['edge_labeled', 'edge_label_num']):
        is_el = is_edge_labeled(Gn)
        edge_label_num = get_edge_label_num(Gn)

    if 'edge_labeled' in attr_names:
        # graphs are considered edge unlabeled if all edges have the same label.
        attrs.update({'edge_labeled': is_el if edge_label_num > 1 else False})

    if 'edge_label_num' in attr_names:
        attrs.update({'edge_label_num': edge_label_num})

    if 'is_directed' in attr_names:
        attrs.update({'is_directed': is_directed(Gn)})

    if 'ave_node_degree' in attr_names:
        attrs.update({'ave_node_degree': get_ave_node_degree(Gn)})

    if 'max_node_degree' in attr_names:
        attrs.update({'max_node_degree': get_max_node_degree(Gn)})

    if 'min_node_degree' in attr_names:
        attrs.update({'min_node_degree': get_min_node_degree(Gn)})
        
    if 'ave_fill_factor' in attr_names:
        attrs.update({'ave_fill_factor': get_ave_fill_factor(Gn)})

    if 'max_fill_factor' in attr_names:
        attrs.update({'max_fill_factor': get_max_fill_factor(Gn)})

    if 'min_fill_factor' in attr_names:
        attrs.update({'min_fill_factor': get_min_fill_factor(Gn)})

    if 'substructures' in attr_names:
        attrs.update({'substructures': get_substructures(Gn)})

    if 'class_number' in attr_names:
        attrs.update({'class_number': get_class_num(target)})

    if 'node_attr_dim' in attr_names:
        attrs['node_attr_dim'] = get_node_attr_dim(Gn)

    if 'edge_attr_dim' in attr_names:
        attrs['edge_attr_dim'] = get_edge_attr_dim(Gn)

    from collections import OrderedDict
    return OrderedDict(
        sorted(attrs.items(), key=lambda i: attr_names.index(i[0])))

""" Utilities function to manage graph files
"""
from os.path import dirname, splitext

def loadCT(filename):
    """load data from a Chemical Table (.ct) file.

    Notes
    ------
    a typical example of data in .ct is like this:

     3 2  <- number of nodes and edges
        0.0000    0.0000    0.0000 C <- each line describes a node (x,y,z + label)
        0.0000    0.0000    0.0000 C
        0.0000    0.0000    0.0000 O
      1  3  1  1 <- each line describes an edge : to, from, bond type, bond stereo
      2  3  1  1
      
    Check https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=10&ved=2ahUKEwivhaSdjsTlAhVhx4UKHczHA8gQFjAJegQIARAC&url=https%3A%2F%2Fwww.daylight.com%2Fmeetings%2Fmug05%2FKappler%2Fctfile.pdf&usg=AOvVaw1cDNrrmMClkFPqodlF2inS
    for detailed format discription.
    """
    from os.path import basename
    g = nx.Graph()
    with open(filename) as f:
        content = f.read().splitlines()
        g = nx.Graph(
            name = str(content[0]),
            filename = basename(filename))  # set name of the graph
        tmp = content[1].split(" ")
        if tmp[0] == '':
            nb_nodes = int(tmp[1])  # number of the nodes
            nb_edges = int(tmp[2])  # number of the edges
        else:
            nb_nodes = int(tmp[0])
            nb_edges = int(tmp[1])
            # patch for compatibility : label will be removed later
        for i in range(0, nb_nodes):
            tmp = content[i + 2].split(" ")
            tmp = [x for x in tmp if x != '']
            g.add_node(i, atom=tmp[3].strip(), 
                       label=[item.strip() for item in tmp[3:]], 
                       attributes=[item.strip() for item in tmp[0:3]])
        for i in range(0, nb_edges):
            tmp = content[i + g.number_of_nodes() + 2].split(" ")
            g.add_edge(int(tmp[0]) - 1, int(tmp[1]) - 1,
                       bond_type=tmp[2].strip(),
                       label=[item.strip() for item in tmp[2:]])
    return g

def loadGXL(filename):
    import xml.etree.ElementTree as ET

    tree = ET.parse(filename)
    root = tree.getroot()
    index = 0
    g = nx.Graph(filename=basename(filename), name=root[0].attrib['id'])
    dic = {}  # used to retrieve incident nodes of edges
    for node in root.iter('node'):
        dic[node.attrib['id']] = index
        labels = {}
        for attr in node.iter('attr'):
            labels[attr.attrib['name']] = attr[0].text
        if 'chem' in labels:
            labels['label'] = labels['chem']
            labels['atom'] = labels['chem']
        g.add_node(index, **labels)
        index += 1

    for edge in root.iter('edge'):
        labels = {}
        for attr in edge.iter('attr'):
        if 'valence' in labels:
            labels['label'] = labels['valence']
            labels['bond_type'] = labels['valence']
        g.add_edge(dic[edge.attrib['from']], dic[edge.attrib['to']], **labels)
    return g

def saveGXL(graph, filename, method='benoit'):
    if method == 'benoit':
        import xml.etree.ElementTree as ET
        root_node = ET.Element('gxl')
        attr = dict()
        attr['id'] = str(graph.graph['name'])
        attr['edgeids'] = 'true'
        attr['edgemode'] = 'undirected'
        graph_node = ET.SubElement(root_node, 'graph', attrib=attr)
    
        for v in graph:
            current_node = ET.SubElement(graph_node, 'node', attrib={'id': str(v)})
            for attr in graph.nodes[v].keys():
                cur_attr = ET.SubElement(
                    current_node, 'attr', attrib={'name': attr})
                cur_value = ET.SubElement(cur_attr,
                                          graph.nodes[v][attr].__class__.__name__)
                cur_value.text = graph.nodes[v][attr]
    
        for v1 in graph:
            for v2 in graph[v1]:
                if (v1 < v2):  # Non oriented graphs
                    cur_edge = ET.SubElement(
                        graph_node,
                        'edge',
                        attrib={
                            'from': str(v1),
                            'to': str(v2)
                        })
                    for attr in graph[v1][v2].keys():
                        cur_attr = ET.SubElement(
                            cur_edge, 'attr', attrib={'name': attr})
                        cur_value = ET.SubElement(
                            cur_attr, graph[v1][v2][attr].__class__.__name__)
                        cur_value.text = str(graph[v1][v2][attr])
    
        tree = ET.ElementTree(root_node)
        tree.write(filename)
    elif method == 'gedlib':
        # reference: https://github.com/dbblumenthal/gedlib/blob/master/data/generate_molecules.py#L22
#        pass
        gxl_file = open(filename, 'w')
        gxl_file.write("<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n")
        gxl_file.write("<!DOCTYPE gxl SYSTEM \"http://www.gupro.de/GXL/gxl-1.0.dtd\">\n")
        gxl_file.write("<gxl xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n")
        gxl_file.write("<graph id=\"" + str(graph.graph['name']) + "\" edgeids=\"true\" edgemode=\"undirected\">\n")
        for v, attrs in graph.nodes(data=True):
            gxl_file.write("<node id=\"_" + str(v) + "\">")
            gxl_file.write("<attr name=\"" + "chem" + "\"><int>" + str(attrs['chem']) + "</int></attr>")
            gxl_file.write("</node>\n")
        for v1, v2, attrs in graph.edges(data=True):
            gxl_file.write("<edge from=\"_" + str(v1) + "\" to=\"_" + str(v2) + "\">")
            gxl_file.write("<attr name=\"valence\"><int>" + str(attrs['valence']) + "</int></attr>")
#            gxl_file.write("<attr name=\"valence\"><int>" + "1" + "</int></attr>")
            gxl_file.write("</edge>\n")
        gxl_file.write("</graph>\n")
        gxl_file.write("</gxl>")
        gxl_file.close()
    elif method == 'gedlib-letter':
        # and https://github.com/dbblumenthal/gedlib/blob/master/data/datasets/Letter/HIGH/AP1_0000.gxl
        gxl_file.write("<graph id=\"" + str(graph.graph['name']) + "\" edgeids=\"false\" edgemode=\"undirected\">\n")
            gxl_file.write("<attr name=\"x\"><float>" + str(attrs['attributes'][0]) + "</float></attr>")
            gxl_file.write("<attr name=\"y\"><float>" + str(attrs['attributes'][1]) + "</float></attr>")
            gxl_file.write("<edge from=\"_" + str(v1) + "\" to=\"_" + str(v2) + "\"/>\n")

def loadSDF(filename):
    """load data from structured data file (.sdf file).

    Notes
    ------
    A SDF file contains a group of molecules, represented in the similar way as in MOL format.
    Check http://www.nonlinear.com/progenesis/sdf-studio/v0.9/faq/sdf-file-format-guidance.aspx, 2018 for detailed structure.
    """
    from tqdm import tqdm
    import sys
    data = []
        index = 0
        pbar = tqdm(total=len(content) + 1, desc='load SDF', file=sys.stdout)
        while index < len(content):
            index_old = index

            g = nx.Graph(name=content[index].strip())  # set name of the graph

            tmp = content[index + 3]
            nb_nodes = int(tmp[:3])  # number of the nodes
            nb_edges = int(tmp[3:6])  # number of the edges

            for i in range(0, nb_nodes):
                tmp = content[i + index + 4]
                g.add_node(i, atom=tmp[31:34].strip())

            for i in range(0, nb_edges):
                tmp = content[i + index + g.number_of_nodes() + 4]
                tmp = [tmp[i:i + 3] for i in range(0, len(tmp), 3)]
                g.add_edge(
                    int(tmp[0]) - 1, int(tmp[1]) - 1, bond_type=tmp[2].strip())

            data.append(g)

            index += 4 + g.number_of_nodes() + g.number_of_edges()
            while content[index].strip() != '$$$$':  # seperator
                index += 1
            index += 1

            pbar.update(index - index_old)
        pbar.update(1)
        pbar.close()

    return data

def loadMAT(filename, extra_params):
    """Load graph data from a MATLAB (up to version 7.1) .mat file.

    Notes
    ------
    A MAT file contains a struct array containing graphs, and a column vector lx containing a class label for each graph.
    Check README in downloadable file in http://mlcb.is.tuebingen.mpg.de/Mitarbeiter/Nino/WL/, 2018 for detailed structure.
    """
    from scipy.io import loadmat
    data = []
    content = loadmat(filename)
    order = extra_params['am_sp_al_nl_el']
    # print(content)
    # print('----')
    for key, value in content.items():
        if key[0] == 'l':  # class label
            y = np.transpose(value)[0].tolist()
            # print(y)
        elif key[0] != '_':
            # print(value[0][0][0])
            # print()
            # print(value[0][0][1])
            # print()
            # print(value[0][0][2])
            # print()
            # if len(value[0][0]) > 3:
            #     print(value[0][0][3])
            # print('----')
            # if adjacency matrix is not compressed / edge label exists
            if order[1] == 0:
                for i, item in enumerate(value[0]):
                    # print(item)
                    # print('------')
                    g = nx.Graph(name=i)  # set name of the graph
                    nl = np.transpose(item[order[3]][0][0][0])  # node label
                    # print(item[order[3]])
                    # print()
                    for index, label in enumerate(nl[0]):
                        g.add_node(index, atom=str(label))
                    el = item[order[4]][0][0][0]  # edge label
                    for edge in el:
                        g.add_edge(
                            edge[0] - 1, edge[1] - 1, bond_type=str(edge[2]))
                    data.append(g)
            else:
                from scipy.sparse import csc_matrix
                    # print(nl)
                    # print()
                    sam = item[order[0]]  # sparse adjacency matrix
                    index_no0 = sam.nonzero()
                    for col, row in zip(index_no0[0], index_no0[1]):
                        # print(col)
                        # print(row)
                        g.add_edge(col, row)
                    # print(g.edges(data=True))
    return data, y

def loadTXT(dirname_dataset):
    """Load graph data from a .txt file.

    Notes
    ------
    The graph data is loaded from separate files.
    Check README in downloadable file http://tiny.cc/PK_MLJ_data, 2018 for detailed structure.
    """
    from os import listdir
    from os.path import dirname

    # load data file names
    for name in listdir(dirname_dataset):
        if '_A' in name:
            fam = dirname_dataset + '/' + name
        elif '_graph_indicator' in name:
            fgi = dirname_dataset + '/' + name
        elif '_graph_labels' in name:
            fgl = dirname_dataset + '/' + name
        elif '_node_labels' in name:
            fnl = dirname_dataset + '/' + name
        elif '_edge_labels' in name:
            fel = dirname_dataset + '/' + name
        elif '_edge_attributes' in name:
            fea = dirname_dataset + '/' + name
        elif '_node_attributes' in name:
            fna = dirname_dataset + '/' + name
        elif '_graph_attributes' in name:
            fga = dirname_dataset + '/' + name
        # this is supposed to be the node attrs, make sure to put this as the last 'elif'
        elif '_attributes' in name:

    content_gi = open(fgi).read().splitlines()  # graph indicator
    content_am = open(fam).read().splitlines()  # adjacency matrix
    content_gl = open(fgl).read().splitlines()  # lass labels

    # create graphs and add nodes
    data = [nx.Graph(name=i) for i in range(0, len(content_gl))]
    if 'fnl' in locals():
        content_nl = open(fnl).read().splitlines()  # node labels
        for i, line in enumerate(content_gi):
            # transfer to int first in case of unexpected blanks
            data[int(line) - 1].add_node(i, atom=str(int(content_nl[i])))
    else:
            data[int(line) - 1].add_node(i)
            
    # add edges
    for line in content_am:
        tmp = line.split(',')
        n1 = int(tmp[0]) - 1
        n2 = int(tmp[1]) - 1
        # ignore edge weight here.
        g = int(content_gi[n1]) - 1
        data[g].add_edge(n1, n2)

    # add edge labels
    if 'fel' in locals():
        content_el = open(fel).read().splitlines()
        for index, line in enumerate(content_el):
            label = line.strip()
            n = [int(i) - 1 for i in content_am[index].split(',')]
            g = int(content_gi[n[0]]) - 1
            data[g].edges[n[0], n[1]]['bond_type'] = label

    # add node attributes
    if 'fna' in locals():
        content_na = open(fna).read().splitlines()
        for i, line in enumerate(content_na):
            attrs = [i.strip() for i in line.split(',')]
            g = int(content_gi[i]) - 1
            data[g].nodes[i]['attributes'] = attrs

    # add edge attributes
    if 'fea' in locals():
        content_ea = open(fea).read().splitlines()
        for index, line in enumerate(content_ea):
            data[g].edges[n[0], n[1]]['attributes'] = attrs

    # load y
    y = [int(i) for i in content_gl]

def loadDataset(filename, filename_y=None, extra_params=None):
    """Read graph data from filename and load them as NetworkX graphs.

    filename : string
        The name of the file from where the dataset is read.
    filename_y : string
        The name of file of the targets corresponding to graphs.
    extra_params : dict
        Extra parameters only designated to '.mat' format.

    Return
    ------
    data : List of NetworkX graph.
    y : List
        Targets corresponding to graphs.
        
    Notes
    -----
    This function supports following graph dataset formats:
    'ds': load data from .ds file. See comments of function loadFromDS for a example.
    'cxl': load data from Graph eXchange Language file (.cxl file). See 
        http://www.gupro.de/GXL/Introduction/background.html, 2019 for detail.
    'sdf': load data from structured data file (.sdf file). See 
        http://www.nonlinear.com/progenesis/sdf-studio/v0.9/faq/sdf-file-format-guidance.aspx, 
        2018 for details.
    'mat': Load graph data from a MATLAB (up to version 7.1) .mat file. See
        README in downloadable file in http://mlcb.is.tuebingen.mpg.de/Mitarbeiter/Nino/WL/, 
    'txt': Load graph data from a special .txt file. See
        https://ls11-www.cs.tu-dortmund.de/staff/morris/graphkerneldatasets,
        2019 for details. Note here filename is the name of either .txt file in
        the dataset directory.
    """
    extension = splitext(filename)[1][1:]
    if extension == "ds":
        data, y = loadFromDS(filename, filename_y)
    elif extension == "cxl":
        
        dirname_dataset = dirname(filename)
        tree = ET.parse(filename)
        root = tree.getroot()
        data = []
        y = []
        for graph in root.iter('graph'):
            mol_filename = graph.attrib['file']
            mol_class = graph.attrib['class']
            data.append(loadGXL(dirname_dataset + '/' + mol_filename))
            y.append(mol_class)
    elif extension == 'xml':
        data, y = loadFromXML(filename, extra_params)
    elif extension == "sdf":
        import numpy as np
        from tqdm import tqdm
        import sys

        data = loadSDF(filename)

        y_raw = open(filename_y).read().splitlines()
        y_raw.pop(0)
        tmp0 = []
        tmp1 = []
        for i in range(0, len(y_raw)):
            tmp = y_raw[i].split(',')
            tmp0.append(tmp[0])
            tmp1.append(tmp[1].strip())

        y = []
        for i in tqdm(range(0, len(data)), desc='ajust data', file=sys.stdout):
            try:
                y.append(tmp1[tmp0.index(data[i].name)].strip())
            except ValueError:  # if data[i].name not in tmp0
                data[i] = []
        data = list(filter(lambda a: a != [], data))
    elif extension == "mat":
        data, y = loadMAT(filename, extra_params)
    elif extension == 'txt':
        data, y = loadTXT(dirname_dataset)
        # print(len(y))
        # print(y)
        # print(data[0].nodes(data=True))
        # print('----')
        # print(data[0].edges(data=True))
        # for g in data:
        #     print(g.nodes(data=True))
        #     print('----')
        #     print(g.edges(data=True))

def loadFromXML(filename, extra_params):
    
    if extra_params:
        dirname_dataset = extra_params
    else:
    data = []
    y = []
    for graph in root.iter('graph'):
        mol_filename = graph.attrib['file']
        mol_class = graph.attrib['class']
        data.append(loadGXL(dirname_dataset + '/' + mol_filename))
        y.append(mol_class)
        
def loadFromDS(filename, filename_y):
    """Load data from .ds file.
    Possible graph formats include:
        '.ct': see function loadCT for detail.
        '.gxl': see dunction loadGXL for detail.
    Note these graph formats are checked automatically by the extensions of 
    graph files.
    """
    dirname_dataset = dirname(filename)
    data = []
    y = []
    content = open(filename).read().splitlines()
    extension = splitext(content[0].split(' ')[0])[1][1:]
    if filename_y is None or filename_y == '':
        if extension == 'ct': 
            for i in range(0, len(content)):
                tmp = content[i].split(' ')
                # remove the '#'s in file names
                data.append(
                    loadCT(dirname_dataset + '/' + tmp[0].replace('#', '', 1)))
                y.append(float(tmp[1]))
        elif extension == 'gxl':
                    loadGXL(dirname_dataset + '/' + tmp[0].replace('#', '', 1)))
    else:  # y in a seperate file
        if extension == 'ct':
                tmp = content[i]
                    loadCT(dirname_dataset + '/' + tmp.replace('#', '', 1)))
                    loadGXL(dirname_dataset + '/' + tmp.replace('#', '', 1)))                                                       
                                                               
        content_y = open(filename_y).read().splitlines()
        # assume entries in filename and filename_y have the same order.
        for item in content_y:
            tmp = item.split(' ')
            # assume the 3rd entry in a line is y (for Alkane dataset)
            y.append(float(tmp[2]))
            
def saveDataset(Gn, y, gformat='gxl', group=None, filename='gfile', xparams=None):
    """Save list of graphs.
    """
    import os
    dirname_ds = os.path.dirname(filename)
    if dirname_ds != '':
        dirname_ds += '/'
        if not os.path.exists(dirname_ds) :
            os.makedirs(dirname_ds)
                
    if 'graph_dir' in xparams:
        graph_dir = xparams['graph_dir'] + '/'
        if not os.path.exists(graph_dir):
            os.makedirs(graph_dir)
    else:
        graph_dir = dirname_ds 
            
    if group == 'xml' and gformat == 'gxl':
        with open(filename + '.xml', 'w') as fgroup:
            fgroup.write("<?xml version=\"1.0\"?>")
            fgroup.write("\n<!DOCTYPE GraphCollection SYSTEM \"http://www.inf.unibz.it/~blumenthal/dtd/GraphCollection.dtd\">")
            fgroup.write("\n<GraphCollection>")
            for idx, g in enumerate(Gn):
                fname_tmp = "graph" + str(idx) + ".gxl"
                saveGXL(g, graph_dir + fname_tmp, method=xparams['method'])
                fgroup.write("\n\t<graph file=\"" + fname_tmp + "\" class=\"" + str(y[idx]) + "\"/>")
            fgroup.write("\n</GraphCollection>")
            fgroup.close()
            
if __name__ == '__main__':    
#    ### Load dataset from .ds file.
#    # .ct files.
#    ds = {'name': 'Alkane', 'dataset': '../../datasets/Alkane/dataset.ds',
#        'dataset_y': '../../datasets/Alkane/dataset_boiling_point_names.txt'}
#    Gn, y = loadDataset(ds['dataset'], filename_y=ds['dataset_y'])
##    ds = {'name': 'Acyclic', 'dataset': '../../datasets/acyclic/dataset_bps.ds'}  # node symb
##    Gn, y = loadDataset(ds['dataset'])
##    ds = {'name': 'MAO', 'dataset': '../../datasets/MAO/dataset.ds'} # node/edge symb
##    ds = {'name': 'PAH', 'dataset': '../../datasets/PAH/dataset.ds'} # unlabeled
#    print(Gn[1].nodes(data=True))
#    print(Gn[1].edges(data=True))
#    print(y[1])
    
#    # .gxl file.
#    ds = {'name': 'monoterpenoides', 
#          'dataset': '../../datasets/monoterpenoides/dataset_10+.ds'}  # node/edge symb
#    Gn, y = loadDataset(ds['dataset'])
    
    ### Convert graph from one format to another.
    # .gxl file.
    ds = {'name': 'monoterpenoides', 
          'dataset': '../../datasets/monoterpenoides/dataset_10+.ds'}  # node/edge symb
    Gn, y = loadDataset(ds['dataset'])
    y = [int(i) for i in y]
    print(Gn[1].nodes(data=True))
    print(Gn[1].edges(data=True))
    print(y[1])
    # Convert a graph to the proper NetworkX format that can be recognized by library gedlib.
    Gn_new = []
    for G in Gn:
        G_new = nx.Graph()
        for nd, attrs in G.nodes(data=True):
            G_new.add_node(str(nd), chem=attrs['atom'])
        for nd1, nd2, attrs in G.edges(data=True):
            G_new.add_edge(str(nd1), str(nd2), valence=attrs['bond_type'])
#            G_new.add_edge(str(nd1), str(nd2))
        Gn_new.append(G_new)
    print(Gn_new[1].nodes(data=True))
    print(Gn_new[1].edges(data=True))
    print(Gn_new[1])
    filename = '/media/ljia/DATA/research-repo/codes/others/gedlib/tests_linlin/generated_datsets/monoterpenoides/gxl/monoterpenoides'
    xparams = {'method': 'gedlib'}
    saveDataset(Gn, y, gformat='gxl', group='xml', filename=filename, xparams=xparams)
    
#    ds = {'name': 'MUTAG', 'dataset': '../../datasets/MUTAG/MUTAG.mat',
#          'extra_params': {'am_sp_al_nl_el': [0, 0, 3, 1, 2]}}  # node/edge symb
#    Gn, y = loadDataset(ds['dataset'], extra_params=ds['extra_params'])
#    saveDataset(Gn, y, group='xml', filename='temp/temp')

""" Functions for python system.
"""

def isNotebook():
    """check if code is executed in the IPython notebook.
    """
    try:
        shell = get_ipython().__class__.__name__
        if shell == 'ZMQInteractiveShell':
            return True   # Jupyter notebook or qtconsole
        elif shell == 'TerminalInteractiveShell':
            return False  # Terminal running IPython
        else:
            return False  # Other type (?)
    except NameError:
        return False      # Probably standard Python interpreter

"""Those who are not graph kernels. We can be kernels for nodes or edges!
These kernels are defined between pairs of vectors.
"""
import numpy as np

def deltakernel(x, y):
    """Delta kernel. Return 1 if x == y, 0 otherwise.

    x, y : any
        Two parts to compare.

    Return
    ------
    kernel : integer
        Delta kernel.

    References
    [1] H. Kashima, K. Tsuda, and A. Inokuchi. Marginalized kernels between 
    labeled graphs. In Proceedings of the 20th International Conference on 
    Machine Learning, Washington, DC, United States, 2003.
    """
    return x == y  #(1 if condition else 0)

def gaussiankernel(x, y, gamma=None):
    """Gaussian kernel.
    Compute the rbf (gaussian) kernel between x and y:

        K(x, y) = exp(-gamma ||x-y||^2).

    Read more in the :ref:`User Guide <rbf_kernel>`.

    x, y : array

    gamma : float, default None
        If None, defaults to 1.0 / n_features

    Returns
    -------
    kernel : float
    """
    if gamma is None:
        gamma = 1.0 / len(x)

    xt = np.array([float(itm) for itm in x])
    yt = np.array([float(itm) for itm in y])
    kernel = xt - yt
    kernel = kernel ** 2
    kernel = np.sum(kernel)
    kernel *= -gamma
    kernel = np.exp(kernel)
    return kernel

def polynomialkernel(x, y, d=1, c=0):
    """Polynomial kernel.
    Compute the polynomial kernel between x and y:

        K(x, y) = <x, y> ^d + c.

    x, y : array

    d : integer, default 1
    
    c : float, default 0

    Returns
    -------
    """
    return np.dot(x, y) ** d + c

def linearkernel(x, y):

        K(x, y) = <x, y>.

    x, y : array

    Returns
    -------
    """
    return np.dot(x, y)

def kernelsum(k1, k2, d11, d12, d21=None, d22=None, lamda1=1, lamda2=1):
    """Sum of a pair of kernels.

    k = lamda1 * k1(d11, d12) + lamda2 * k2(d21, d22)

    k1, k2 : function
        A pair of kernel functions.
    d11, d12:
        Inputs of k1. If d21 or d22 is None, apply d11, d12 to both k1 and k2.
    d21, d22:
        Inputs of k2.
    lamda1, lamda2: float
        Coefficients of the product.

    Return
    ------

    """
    if d21 == None or d22 == None:
        kernel = lamda1 * k1(d11, d12) + lamda2 * k2(d11, d12)
    else:
        kernel = lamda1 * k1(d11, d12) + lamda2 * k2(d21, d22)

def kernelproduct(k1, k2, d11, d12, d21=None, d22=None, lamda=1):
    """Product of a pair of kernels.

    k = lamda * k1(d11, d12) * k2(d21, d22)

    d11, d12:
    d21, d22:
    lamda: float
        Coefficient of the product.

    Return
    ------
    """
        kernel = lamda * k1(d11, d12) * k2(d11, d12)
    else:
        kernel = lamda * k1(d11, d12) * k2(d21, d22)

if __name__ == '__main__':
    o = polynomialkernel([1, 2], [3, 4], 2, 3)

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Nov  8 14:21:25 2019

@author: ljia
"""

import sys
import time

class Logger(object):
    def __init__(self):
        self.terminal = sys.stdout
        self.log = open("log." + str(time.time()) + ".log", "a")

    def write(self, message):
        self.terminal.write(message)
        self.log.write(message)  

    def flush(self):
        #this flush method is needed for python 3 compatibility.
        #this handles the flush command by doing nothing.
        #you might want to specify some extra behavior here.
        pass    

sys.stdout = Logger()

import matplotlib
matplotlib.use('Agg')
from matplotlib import pyplot as plt
from sklearn.kernel_ridge import KernelRidge
from sklearn.svm import SVC
from sklearn.metrics import accuracy_score, mean_squared_error
from sklearn.model_selection import KFold, train_test_split, ParameterGrid

#from joblib import Parallel, delayed
from multiprocessing import Pool, Array
from functools import partial
import sys
sys.path.insert(0, "../")
import os
import datetime
#from os.path import basename, splitext
from pygraph.utils.graphfiles import loadDataset
from tqdm import tqdm

#from memory_profiler import profile

#@profile
def model_selection_for_precomputed_kernel(datafile,
                                           estimator,
                                           param_grid_precomputed,
                                           param_grid,
                                           model_type,
                                           NUM_TRIALS=30,
                                           datafile_y=None,
                                           extra_params=None,
                                           ds_name='ds-unknown',
                                           n_jobs=1,
                                           read_gm_from_file=False,
                                           verbose=True):
    """Perform model selection, fitting and testing for precomputed kernels 
    using nested CV. Print out neccessary data during the process then finally 
    the results.

    datafile : string
        Path of dataset file.
    estimator : function
        kernel function used to estimate. This function needs to return a gram matrix.
    param_grid_precomputed : dictionary
        Dictionary with names (string) of parameters used to calculate gram 
        matrices as keys and lists of parameter settings to try as values. This 
        enables searching over any sequence of parameter settings. Params with 
        length 1 will be omitted.
    param_grid : dictionary
        Dictionary with names (string) of parameters used as penelties as keys 
        and lists of parameter settings to try as values. This enables 
        searching over any sequence of parameter settings. Params with length 1
        will be omitted.
    model_type : string
        Type of the problem, can be 'regression' or 'classification'.
    NUM_TRIALS : integer
        Number of random trials of outer cv loop. The default is 30.
    datafile_y : string
        Path of file storing y data. This parameter is optional depending on 
        the given dataset file.
        Extra parameters for loading dataset. See function pygraph.utils.
        graphfiles.loadDataset for detail.
    ds_name : string
        Name of the dataset.
    n_jobs : int
        Number of jobs for parallelization.
    read_gm_from_file : boolean
        Whether gram matrices are loaded from a file.

    Examples
    --------
    >>> import numpy as np
    >>> import sys
    >>> sys.path.insert(0, "../")
    >>> from pygraph.utils.model_selection_precomputed import model_selection_for_precomputed_kernel
    >>> from pygraph.kernels.untilHPathKernel import untilhpathkernel
    >>>
    >>> datafile = '../datasets/MUTAG/MUTAG_A.txt'
    >>> estimator = untilhpathkernel
    >>> param_grid_precomputed = {’depth’:  np.linspace(1, 10, 10), ’k_func’:
            [’MinMax’, ’tanimoto’], ’compute_method’:  [’trie’]}
    >>> # ’C’ for classification problems and ’alpha’ for regression problems.
    >>> param_grid = [{’C’: np.logspace(-10, 10, num=41, base=10)}, {’alpha’:
            np.logspace(-10, 10, num=41, base=10)}]
    >>>
    >>> model_selection_for_precomputed_kernel(datafile, estimator, 
            param_grid_precomputed, param_grid[0], 'classification', ds_name=’MUTAG’)
    """
    tqdm.monitor_interval = 0

    results_dir = '../notebooks/results/' + estimator.__name__
    if not os.path.exists(results_dir):
        os.makedirs(results_dir)
    # a string to save all the results.
    str_fw = '###################### log time: ' + datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S") + '. ######################\n\n'
    str_fw += '# This file contains results of ' + estimator.__name__ + ' on dataset ' + ds_name + ',\n# including gram matrices, serial numbers for gram matrix figures and performance.\n\n'

    # setup the model type
    model_type = model_type.lower()
    if model_type != 'regression' and model_type != 'classification':
        raise Exception(
            'The model type is incorrect! Please choose from regression or classification.'
        )
    if verbose:
        print()
        print('--- This is a %s problem ---' % model_type)
    str_fw += 'This is a %s problem.\n' % model_type
    
    # calculate gram matrices rather than read them from file.
    if read_gm_from_file == False:
        # Load the dataset
        if verbose:
            print()
            print('\n1. Loading dataset from file...')
        if isinstance(datafile, str):
            dataset, y_all = loadDataset(
                    datafile, filename_y=datafile_y, extra_params=extra_params)
        else: # load data directly from variable.
            dataset = datafile
            y_all = datafile_y                

        #     import matplotlib.pyplot as plt
        #     import networkx as nx
        #     nx.draw_networkx(dataset[30])
        #     plt.show()
    
        # Grid of parameters with a discrete number of values for each.
        param_list_precomputed = list(ParameterGrid(param_grid_precomputed))
        param_list = list(ParameterGrid(param_grid))
    
        gram_matrices = [
        ]  # a list to store gram matrices for all param_grid_precomputed
        gram_matrix_time = [
        ]  # a list to store time to calculate gram matrices
        param_list_pre_revised = [
        ]  # list to store param grids precomputed ignoring the useless ones
    
        # calculate all gram matrices
            print()
            print('2. Calculating gram matrices. This could take a while...')
        str_fw += '\nII. Gram matrices.\n\n'
        tts = time.time()  # start training time
        nb_gm_ignore = 0  # the number of gram matrices those should not be considered, as they may contain elements that are not numbers (NaN)
        for idx, params_out in enumerate(param_list_precomputed):
            y = y_all[:]
            params_out['n_jobs'] = n_jobs
            params_out['verbose'] = verbose
#            print(dataset)
#            import networkx as nx
#            nx.draw_networkx(dataset[1])
#            plt.show()
            rtn_data = estimator(dataset[:], **params_out)
            Kmatrix = rtn_data[0]
            current_run_time = rtn_data[1]
            # for some kernels, some graphs in datasets may not meet the 
            # kernels' requirements for graph structure. These graphs are trimmed. 
            if len(rtn_data) == 3:
                idx_trim = rtn_data[2]  # the index of trimmed graph list
                y = [y[idxt] for idxt in idx_trim] # trim y accordingly
#            Kmatrix = np.random.rand(2250, 2250)
#            current_run_time = 0.1
            
            # remove graphs whose kernels with themselves are zeros 
            # @todo: y not changed accordingly?
            Kmatrix_diag = Kmatrix.diagonal().copy()
            nb_g_ignore = 0
            for idxk, diag in enumerate(Kmatrix_diag):
                if diag == 0:
                    Kmatrix = np.delete(Kmatrix, (idxk - nb_g_ignore), axis=0)
                    Kmatrix = np.delete(Kmatrix, (idxk - nb_g_ignore), axis=1)
                    nb_g_ignore += 1
            # normalization
            # @todo: works only for undirected graph?
            for i in range(len(Kmatrix)):
                for j in range(i, len(Kmatrix)):
                    Kmatrix[i][j] /= np.sqrt(Kmatrix_diag[i] * Kmatrix_diag[j])
                    Kmatrix[j][i] = Kmatrix[i][j]
            if verbose:
                print()
            if params_out == {}:
                if verbose:
                    print('the gram matrix is: ')
                str_fw += 'the gram matrix is:\n\n'
            else:
                    print('the gram matrix with parameters', params_out, 'is: \n\n')
                str_fw += 'the gram matrix with parameters %s is:\n\n' % params_out
            if len(Kmatrix) < 2:
                nb_gm_ignore += 1
                    print('ignored, as at most only one of all its diagonal value is non-zero.')
                str_fw += 'ignored, as at most only one of all its diagonal value is non-zero.\n\n'
            else:                
                if np.isnan(Kmatrix).any(
                ):  # if the matrix contains elements that are not numbers
                    nb_gm_ignore += 1
                    if verbose:
                        print('ignored, as it contains elements that are not numbers.')
                    str_fw += 'ignored, as it contains elements that are not numbers.\n\n'
                else:
#                    print(Kmatrix)
                    str_fw += np.array2string(
                            Kmatrix,
                            separator=',') + '\n\n'
#                            separator=',',
#                            threshold=np.inf,
#                            floatmode='unique') + '\n\n'

                    fig_file_name = results_dir + '/GM[ds]' + ds_name
                    if params_out != {}:
                        fig_file_name += '[params]' + str(idx)
                    plt.imshow(Kmatrix)
                    plt.colorbar()
                    plt.savefig(fig_file_name + '.eps', format='eps', dpi=300)
#                    plt.show()
                    plt.clf()
                    gram_matrices.append(Kmatrix)
                    gram_matrix_time.append(current_run_time)
                    param_list_pre_revised.append(params_out)
                    if nb_g_ignore > 0:
                        if verbose:
                            print(', where %d graphs are ignored as their graph kernels with themselves are zeros.' % nb_g_ignore)
                        str_fw += ', where %d graphs are ignored as their graph kernels with themselves are zeros.' % nb_g_ignore
            print()
            print(
            '{} gram matrices are calculated, {} of which are ignored.'.format(
                len(param_list_precomputed), nb_gm_ignore))
        str_fw += '{} gram matrices are calculated, {} of which are ignored.\n\n'.format(len(param_list_precomputed), nb_gm_ignore)
        str_fw += 'serial numbers of gram matrix figures and their corresponding parameters settings:\n\n'
        str_fw += ''.join([
            '{}: {}\n'.format(idx, params_out)
            for idx, params_out in enumerate(param_list_precomputed)
        ])

            print()
        if len(gram_matrices) == 0:
                print('all gram matrices are ignored, no results obtained.')
            str_fw += '\nall gram matrices are ignored, no results obtained.\n\n'
        else:
            # save gram matrices to file.
#            np.savez(results_dir + '/' + ds_name + '.gm', 
#                     gms=gram_matrices, params=param_list_pre_revised, y=y, 
#                     gmtime=gram_matrix_time)
                print(
                '3. Fitting and predicting using nested cross validation. This could really take a while...'
                )
            
            # ---- use pool.imap_unordered to parallel and track progress. ----
#            train_pref = []
#            val_pref = []
#            test_pref = []
#            def func_assign(result, var_to_assign):
#                for idx, itm in enumerate(var_to_assign):
#                    itm.append(result[idx])                
#            trial_do_partial = partial(trial_do, param_list_pre_revised, param_list, y, model_type)
#                      
#            parallel_me(trial_do_partial, range(NUM_TRIALS), func_assign, 
#                        [train_pref, val_pref, test_pref], glbv=gram_matrices,
#                        method='imap_unordered', n_jobs=n_jobs, chunksize=1,
#                        itr_desc='cross validation')
            
            def init_worker(gms_toshare):
                global G_gms
                G_gms = gms_toshare
            
#            gram_matrices = np.array(gram_matrices)
#            gms_shape = gram_matrices.shape
#            gms_array = Array('d', np.reshape(gram_matrices.copy(), -1, order='C'))
#            pool = Pool(processes=n_jobs, initializer=init_worker, initargs=(gms_array, gms_shape))
            pool = Pool(processes=n_jobs, initializer=init_worker, initargs=(gram_matrices,))
            trial_do_partial = partial(parallel_trial_do, param_list_pre_revised, param_list, y, model_type)
            train_pref = []
            val_pref = []
            test_pref = []
#            if NUM_TRIALS < 1000 * n_jobs:
#                chunksize = int(NUM_TRIALS / n_jobs) + 1
#            else:
#                chunksize = 1000
            chunksize = 1
                iterator = tqdm(pool.imap_unordered(trial_do_partial, 
                        range(NUM_TRIALS), chunksize), desc='cross validation', file=sys.stdout)
            else:
                iterator = pool.imap_unordered(trial_do_partial, range(NUM_TRIALS), chunksize)
            for o1, o2, o3 in iterator:
                train_pref.append(o1)
                val_pref.append(o2)
                test_pref.append(o3)
            pool.close()
            pool.join()
    
#            # ---- use pool.map to parallel. ----
#            pool =  Pool(n_jobs)
#            trial_do_partial = partial(trial_do, param_list_pre_revised, param_list, gram_matrices, y[0:250], model_type)
#            result_perf = pool.map(trial_do_partial, range(NUM_TRIALS))
#            train_pref = [item[0] for item in result_perf]
#            val_pref = [item[1] for item in result_perf]
#            test_pref = [item[2] for item in result_perf]
    
#            # ---- direct running, normally use a single CPU core. ----
#            for i in tqdm(range(NUM_TRIALS), desc='cross validation', file=sys.stdout):
#                o1, o2, o3 = trial_do(param_list_pre_revised, param_list, gram_matrices, y, model_type, i)
#                train_pref.append(o1)
#                val_pref.append(o2)
#                test_pref.append(o3)
#            print()
    
                print()
                print('4. Getting final performance...')
            str_fw += '\nIII. Performance.\n\n'
            # averages and confidences of performances on outer trials for each combination of parameters
            average_train_scores = np.mean(train_pref, axis=0)
#            print('val_pref: ', val_pref[0][0])
            average_val_scores = np.mean(val_pref, axis=0)
#            print('test_pref: ', test_pref[0][0])
            average_perf_scores = np.mean(test_pref, axis=0)
            # sample std is used here
            std_train_scores = np.std(train_pref, axis=0, ddof=1)
            std_val_scores = np.std(val_pref, axis=0, ddof=1)
            std_perf_scores = np.std(test_pref, axis=0, ddof=1)
    
            if model_type == 'regression':
                best_val_perf = np.amin(average_val_scores)
            else:
                best_val_perf = np.amax(average_val_scores)
#            print('average_val_scores: ', average_val_scores)
#            print('best_val_perf: ', best_val_perf)
#            print()
            best_params_index = np.where(average_val_scores == best_val_perf)
            # find smallest val std with best val perf.
            best_val_stds = [
                std_val_scores[value][best_params_index[1][idx]]
                for idx, value in enumerate(best_params_index[0])
            ]
            min_val_std = np.amin(best_val_stds)
            best_params_index = np.where(std_val_scores == min_val_std)
            best_params_out = [
                param_list_pre_revised[i] for i in best_params_index[0]
            ]
            best_params_in = [param_list[i] for i in best_params_index[1]]
                print('best_params_out: ', best_params_out)
                print('best_params_in: ', best_params_in)
                print()
                print('best_val_perf: ', best_val_perf)
                print('best_val_std: ', min_val_std)
            str_fw += 'best settings of hyper-params to build gram matrix: %s\n' % best_params_out
            str_fw += 'best settings of other hyper-params: %s\n\n' % best_params_in
            str_fw += 'best_val_perf: %s\n' % best_val_perf
            str_fw += 'best_val_std: %s\n' % min_val_std
    
#            print(best_params_index)
#            print(best_params_index[0])
#            print(average_perf_scores)
            final_performance = [
                average_perf_scores[value][best_params_index[1][idx]]
            ]
            final_confidence = [
                std_perf_scores[value][best_params_index[1][idx]]
            ]
                print('final_performance: ', final_performance)
                print('final_confidence: ', final_confidence)
            str_fw += 'final_performance: %s\n' % final_performance
            str_fw += 'final_confidence: %s\n' % final_confidence
            train_performance = [
                average_train_scores[value][best_params_index[1][idx]]
            ]
            train_std = [
                std_train_scores[value][best_params_index[1][idx]]
            ]
                print('train_performance: %s' % train_performance)
                print('train_std: ', train_std)
            str_fw += 'train_performance: %s\n' % train_performance
            str_fw += 'train_std: %s\n\n' % train_std

                print()
            tt_total = time.time() - tts  # training time for all hyper-parameters
            average_gram_matrix_time = np.mean(gram_matrix_time)
            std_gram_matrix_time = np.std(gram_matrix_time, ddof=1) if len(gram_matrix_time) > 1 else 0
            best_gram_matrix_time = [
                gram_matrix_time[i] for i in best_params_index[0]
            ]
            ave_bgmt = np.mean(best_gram_matrix_time)
            std_bgmt = np.std(best_gram_matrix_time, ddof=1) if len(best_gram_matrix_time) > 1 else 0
                print('time to calculate gram matrix with different hyper-params: {:.2f}±{:.2f}s'
                      .format(average_gram_matrix_time, std_gram_matrix_time))
                print('time to calculate best gram matrix: {:.2f}±{:.2f}s'.format(
                        ave_bgmt, std_bgmt))
                print('total training time with all hyper-param choices: {:.2f}s'.format(
                        tt_total))
            str_fw += 'time to calculate gram matrix with different hyper-params: {:.2f}±{:.2f}s\n'.format(average_gram_matrix_time, std_gram_matrix_time)
            str_fw += 'time to calculate best gram matrix: {:.2f}±{:.2f}s\n'.format(ave_bgmt, std_bgmt)
            str_fw += 'total training time with all hyper-param choices: {:.2f}s\n\n'.format(tt_total)
    
            # # save results to file
            # np.savetxt(results_name_pre + 'average_train_scores.dt',
            #            average_train_scores)
            # np.savetxt(results_name_pre + 'average_val_scores', average_val_scores)
            # np.savetxt(results_name_pre + 'average_perf_scores.dt',
            #            average_perf_scores)
            # np.savetxt(results_name_pre + 'std_train_scores.dt', std_train_scores)
            # np.savetxt(results_name_pre + 'std_val_scores.dt', std_val_scores)
            # np.savetxt(results_name_pre + 'std_perf_scores.dt', std_perf_scores)
    
            # np.save(results_name_pre + 'best_params_index', best_params_index)
            # np.save(results_name_pre + 'best_params_pre.dt', best_params_out)
            # np.save(results_name_pre + 'best_params_in.dt', best_params_in)
            # np.save(results_name_pre + 'best_val_perf.dt', best_val_perf)
            # np.save(results_name_pre + 'best_val_std.dt', best_val_std)
            # np.save(results_name_pre + 'final_performance.dt', final_performance)
            # np.save(results_name_pre + 'final_confidence.dt', final_confidence)
            # np.save(results_name_pre + 'train_performance.dt', train_performance)
            # np.save(results_name_pre + 'train_std.dt', train_std)
    
            # np.save(results_name_pre + 'gram_matrix_time.dt', gram_matrix_time)
            # np.save(results_name_pre + 'average_gram_matrix_time.dt',
            #         average_gram_matrix_time)
            # np.save(results_name_pre + 'std_gram_matrix_time.dt',
            #         std_gram_matrix_time)
            # np.save(results_name_pre + 'best_gram_matrix_time.dt',
            #         best_gram_matrix_time)
    
    # read gram matrices from file.
    else:    
#        param_list_precomputed = list(ParameterGrid(param_grid_precomputed))
    
        # read gram matrices from file.
            print()
            print('2. Reading gram matrices from file...')
        str_fw += '\nII. Gram matrices.\n\nGram matrices are read from file, see last log for detail.\n'
        gmfile = np.load(results_dir + '/' + ds_name + '.gm.npz')
        gram_matrices = gmfile['gms'] # a list to store gram matrices for all param_grid_precomputed
        gram_matrix_time = gmfile['gmtime'] # time used to compute the gram matrices
        param_list_pre_revised = gmfile['params'] # list to store param grids precomputed ignoring the useless ones
        y = gmfile['y'].tolist()
        
#        nb_gm_ignore = 0  # the number of gram matrices those should not be considered, as they may contain elements that are not numbers (NaN)            
            print(
                    '3. Fitting and predicting using nested cross validation. This could really take a while...'
                    )
 
        # ---- use pool.imap_unordered to parallel and track progress. ----
        def init_worker(gms_toshare):
            global G_gms
            G_gms = gms_toshare

        pool = Pool(processes=n_jobs, initializer=init_worker, initargs=(gram_matrices,))
        trial_do_partial = partial(parallel_trial_do, param_list_pre_revised, param_list, y, model_type)
        train_pref = []
        val_pref = []
        test_pref = []
        chunksize = 1
            iterator = tqdm(pool.imap_unordered(trial_do_partial, 
                    range(NUM_TRIALS), chunksize), desc='cross validation', file=sys.stdout)
        else:
            iterator = pool.imap_unordered(trial_do_partial, range(NUM_TRIALS), chunksize)
        for o1, o2, o3 in iterator:
            train_pref.append(o1)
            val_pref.append(o2)
            test_pref.append(o3)
        pool.close()
        pool.join()
        
        # # ---- use pool.map to parallel. ----
        # result_perf = pool.map(trial_do_partial, range(NUM_TRIALS))
        # train_pref = [item[0] for item in result_perf]
        # val_pref = [item[1] for item in result_perf]
        # test_pref = [item[2] for item in result_perf]

        # # ---- use joblib.Parallel to parallel and track progress. ----
        # trial_do_partial = partial(trial_do, param_list_pre_revised, param_list, gram_matrices, y, model_type)
        # result_perf = Parallel(n_jobs=n_jobs, verbose=10)(delayed(trial_do_partial)(trial) for trial in range(NUM_TRIALS))

#        # ---- direct running, normally use a single CPU core. ----
#        train_pref = []
#        val_pref = []
#        test_pref = []
#        for i in tqdm(range(NUM_TRIALS), desc='cross validation', file=sys.stdout):
#            o1, o2, o3 = trial_do(param_list_pre_revised, param_list, gram_matrices, y, model_type, i)
#            train_pref.append(o1)
#            val_pref.append(o2)
#            test_pref.append(o3)

            print()
            print('4. Getting final performance...')
        str_fw += '\nIII. Performance.\n\n'
        # averages and confidences of performances on outer trials for each combination of parameters
        average_train_scores = np.mean(train_pref, axis=0)
        average_val_scores = np.mean(val_pref, axis=0)
        average_perf_scores = np.mean(test_pref, axis=0)
        # sample std is used here
        std_train_scores = np.std(train_pref, axis=0, ddof=1)
        std_val_scores = np.std(val_pref, axis=0, ddof=1)
        std_perf_scores = np.std(test_pref, axis=0, ddof=1)

        if model_type == 'regression':
            best_val_perf = np.amin(average_val_scores)
        else:
            best_val_perf = np.amax(average_val_scores)
        best_params_index = np.where(average_val_scores == best_val_perf)
        # find smallest val std with best val perf.
        best_val_stds = [
            std_val_scores[value][best_params_index[1][idx]]
            for idx, value in enumerate(best_params_index[0])
        ]
        min_val_std = np.amin(best_val_stds)
        best_params_index = np.where(std_val_scores == min_val_std)
        best_params_out = [
            param_list_pre_revised[i] for i in best_params_index[0]
        ]
        best_params_in = [param_list[i] for i in best_params_index[1]]
            print('best_params_out: ', best_params_out)
            print('best_params_in: ', best_params_in)
            print()
            print('best_val_perf: ', best_val_perf)
            print('best_val_std: ', min_val_std)
        str_fw += 'best settings of hyper-params to build gram matrix: %s\n' % best_params_out
        str_fw += 'best settings of other hyper-params: %s\n\n' % best_params_in
        str_fw += 'best_val_perf: %s\n' % best_val_perf
        str_fw += 'best_val_std: %s\n' % min_val_std

        final_performance = [
            average_perf_scores[value][best_params_index[1][idx]]
        ]
        final_confidence = [
            std_perf_scores[value][best_params_index[1][idx]]
        ]
            print('final_performance: ', final_performance)
            print('final_confidence: ', final_confidence)
        str_fw += 'final_performance: %s\n' % final_performance
        str_fw += 'final_confidence: %s\n' % final_confidence
        train_performance = [
            average_train_scores[value][best_params_index[1][idx]]
        ]
        train_std = [
            std_train_scores[value][best_params_index[1][idx]]
        ]
            print('train_performance: %s' % train_performance)
            print('train_std: ', train_std)
        str_fw += 'train_performance: %s\n' % train_performance
        str_fw += 'train_std: %s\n\n' % train_std

            print()
        average_gram_matrix_time = np.mean(gram_matrix_time)
        std_gram_matrix_time = np.std(gram_matrix_time, ddof=1) if len(gram_matrix_time) > 1 else 0
        best_gram_matrix_time = [
            gram_matrix_time[i] for i in best_params_index[0]
        ]
        ave_bgmt = np.mean(best_gram_matrix_time)
        std_bgmt = np.std(best_gram_matrix_time, ddof=1) if len(best_gram_matrix_time) > 1 else 0
        if verbose:        
            print(
                    'time to calculate gram matrix with different hyper-params: {:.2f}±{:.2f}s'
                    .format(average_gram_matrix_time, std_gram_matrix_time))
            print('time to calculate best gram matrix: {:.2f}±{:.2f}s'.format(
                    ave_bgmt, std_bgmt))
        tt_poster = time.time() - tts  # training time with hyper-param choices who did not participate in calculation of gram matrices
            print(
                    'training time with hyper-param choices who did not participate in calculation of gram matrices: {:.2f}s'.format(
                            tt_poster))
            print('total training time with all hyper-param choices: {:.2f}s'.format(
                    tt_poster + np.sum(gram_matrix_time)))
#        str_fw += 'time to calculate gram matrix with different hyper-params: {:.2f}±{:.2f}s\n'.format(average_gram_matrix_time, std_gram_matrix_time)
#        str_fw += 'time to calculate best gram matrix: {:.2f}±{:.2f}s\n'.format(ave_bgmt, std_bgmt)
        str_fw += 'training time with hyper-param choices who did not participate in calculation of gram matrices: {:.2f}s\n\n'.format(tt_poster)

        # open file to save all results for this dataset.
        if not os.path.exists(results_dir):
            os.makedirs(results_dir)
            
    # print out results as table.
    str_fw += printResultsInTable(param_list, param_list_pre_revised, average_val_scores,
              std_val_scores, average_perf_scores, std_perf_scores,
              average_train_scores, std_train_scores, gram_matrix_time,
              model_type, verbose)
            
    # open file to save all results for this dataset.
    if not os.path.exists(results_dir + '/' + ds_name + '.output.txt'):
        with open(results_dir + '/' + ds_name + '.output.txt', 'w') as f:
            f.write(str_fw)
    else:
        with open(results_dir + '/' + ds_name + '.output.txt', 'r+') as f:
            content = f.read()
            f.seek(0, 0)
            f.write(str_fw + '\n\n\n' + content)

def trial_do(param_list_pre_revised, param_list, gram_matrices, y, model_type, trial): # Test set level

#    # get gram matrices from global variables.
#    gram_matrices = np.reshape(G_gms.copy(), G_gms_shape, order='C')
    
    # Arrays to store scores
    train_pref = np.zeros((len(param_list_pre_revised), len(param_list)))
    val_pref = np.zeros((len(param_list_pre_revised), len(param_list)))
    test_pref = np.zeros((len(param_list_pre_revised), len(param_list)))

    # randomness added to seeds of split function below. "high" is "size" times
    # 10 so that at least 10 different random output will be yielded. Remove
    # these lines if identical outputs is required.
    rdm_out = np.random.RandomState(seed=None)
    rdm_seed_out_l = rdm_out.uniform(high=len(param_list_pre_revised) * 10, 
                                   size=len(param_list_pre_revised))
#    print(trial, rdm_seed_out_l)
#    print()
    # loop for each outer param tuple
    for index_out, params_out in enumerate(param_list_pre_revised):
        # get gram matrices from global variables.
#        gm_now = G_gms[index_out * G_gms_shape[1] * G_gms_shape[2]:(index_out + 1) * G_gms_shape[1] * G_gms_shape[2]]
#        gm_now = np.reshape(gm_now.copy(), (G_gms_shape[1], G_gms_shape[2]), order='C')
        gm_now = gram_matrices[index_out].copy()
    
        # split gram matrix and y to app and test sets.
        indices = range(len(y))
        # The argument "random_state" in function "train_test_split" can not be
        # set to None, because it will use RandomState instance used by 
        # np.random, which is possible for multiple subprocesses to inherit the
        # same seed if they forked at the same time, leading to identical 
        # random variates for different subprocesses. Instead, we use "trial" 
        # and "index_out" parameters to generate different seeds for different 
        # trials/subprocesses and outer loops. "rdm_seed_out_l" is used to add 
        # randomness into seeds, so that it yields a different output every 
        # time the program is run. To yield identical outputs every time,
        # remove the second line below. Same method is used to the "KFold"
        # function in the inner loop.
        rdm_seed_out = (trial + 1) * (index_out + 1)
        rdm_seed_out = (rdm_seed_out + int(rdm_seed_out_l[index_out])) % (2 ** 32 - 1)
#        print(trial, rdm_seed_out)
        X_app, X_test, y_app, y_test, idx_app, idx_test = train_test_split(
            gm_now, y, indices, test_size=0.1, 
            random_state=rdm_seed_out, shuffle=True)
#        print(trial, idx_app, idx_test)
#        print()
        X_app = X_app[:, idx_app]
        X_test = X_test[:, idx_app]
        y_app = np.array(y_app)
        y_test = np.array(y_test)

        rdm_seed_in_l = rdm_out.uniform(high=len(param_list) * 10, 
                                   size=len(param_list))
        # loop for each inner param tuple
        for index_in, params_in in enumerate(param_list):
#            if trial == 0:
#                print(index_out, index_in)
#                print('params_in: ', params_in)
#            st = time.time()
            rdm_seed_in = (trial + 1) * (index_out + 1) * (index_in + 1)
#            print("rdm_seed_in1: ", trial, index_in, rdm_seed_in)
            rdm_seed_in = (rdm_seed_in + int(rdm_seed_in_l[index_in])) % (2 ** 32 - 1)
#            print("rdm_seed_in2: ", trial, index_in, rdm_seed_in)
            inner_cv = KFold(n_splits=10, shuffle=True, random_state=rdm_seed_in)
            current_train_perf = []
            current_valid_perf = []
            current_test_perf = [] 

            # For regression use the Kernel Ridge method
#            try:
                kr = KernelRidge(kernel='precomputed', **params_in)
                # loop for each split on validation set level
                # validation set level
                for train_index, valid_index in inner_cv.split(X_app):
#                    print("train_index, valid_index: ", trial, index_in, train_index, valid_index)
#                    if trial == 0:
#                        print('train_index: ', train_index)
#                        print('valid_index: ', valid_index)
#                        print('idx_test: ', idx_test)
#                        print('y_app[train_index]: ', y_app[train_index])
#                        print('X_app[train_index, :][:, train_index]: ', X_app[train_index, :][:, train_index])
#                        print('X_app[valid_index, :][:, train_index]: ', X_app[valid_index, :][:, train_index])
                    kr.fit(X_app[train_index, :][:, train_index],
                           y_app[train_index])

                    # predict on the train, validation and test set
                    y_pred_train = kr.predict(
                        X_app[train_index, :][:, train_index])
                    y_pred_valid = kr.predict(
                        X_app[valid_index, :][:, train_index])
#                    if trial == 0:     
#                        print('y_pred_valid: ', y_pred_valid)
#                        print()
                    y_pred_test = kr.predict(
                        X_test[:, train_index])

                    # root mean squared errors
                    current_train_perf.append(
                        np.sqrt(
                            mean_squared_error(
                                y_app[train_index], y_pred_train)))
                    current_valid_perf.append(
                        np.sqrt(
                                y_app[valid_index], y_pred_valid)))
#                        print(mean_squared_error(
#                                y_app[valid_index], y_pred_valid))
                    current_test_perf.append(
                        np.sqrt(
                                y_test, y_pred_test)))
            # For clcassification use SVM
            else:
                svc = SVC(kernel='precomputed', cache_size=200, 
                          verbose=False, **params_in)
#                        np.savez("bug.npy",X_app[train_index, :][:, train_index],y_app[train_index])
                    svc.fit(X_app[train_index, :][:, train_index],
                    
                    y_pred_train = svc.predict(
                    y_pred_valid = svc.predict(
                    y_pred_test = svc.predict(

                        accuracy_score(y_app[train_index],
                                       y_pred_train))
                        accuracy_score(y_app[valid_index],
                                       y_pred_valid))
                        accuracy_score(y_test, y_pred_test))
#            except ValueError:
#                print(sys.exc_info()[0])
#                print(params_out, params_in)

            # average performance on inner splits
            train_pref[index_out][index_in] = np.mean(
                current_train_perf)
            val_pref[index_out][index_in] = np.mean(
                current_valid_perf)
            test_pref[index_out][index_in] = np.mean(
                current_test_perf)
#            print(time.time() - st)
#    if trial == 0:
#        print('val_pref: ', val_pref)
#        print('test_pref: ', test_pref)

    return train_pref, val_pref, test_pref

def parallel_trial_do(param_list_pre_revised, param_list, y, model_type, trial):
    train_pref, val_pref, test_pref = trial_do(param_list_pre_revised, 
                                               param_list, G_gms, y, 
                                               model_type, trial)

def compute_gram_matrices(dataset, y, estimator, param_list_precomputed, 
                          results_dir, ds_name,
                          n_jobs=1, str_fw='', verbose=True):
    gram_matrices = [
    gram_matrix_time = [
    param_list_pre_revised = [
    
    nb_gm_ignore = 0  # the number of gram matrices those should not be considered, as they may contain elements that are not numbers (NaN)
    for idx, params_out in enumerate(param_list_precomputed):
        params_out['n_jobs'] = n_jobs
        rtn_data = estimator(dataset[:], **params_out)
        Kmatrix = rtn_data[0]
        current_run_time = rtn_data[1]
        # for some kernels, some graphs in datasets may not meet the 
        # kernels' requirements for graph structure. These graphs are trimmed. 
        if len(rtn_data) == 3:
            idx_trim = rtn_data[2]  # the index of trimmed graph list
            y = [y[idxt] for idxt in idx_trim] # trim y accordingly

        Kmatrix_diag = Kmatrix.diagonal().copy()
        # remove graphs whose kernels with themselves are zeros
        nb_g_ignore = 0
        for idxk, diag in enumerate(Kmatrix_diag):
            if diag == 0:
                Kmatrix = np.delete(Kmatrix, (idxk - nb_g_ignore), axis=0)
                Kmatrix = np.delete(Kmatrix, (idxk - nb_g_ignore), axis=1)
                nb_g_ignore += 1
        # normalization
        for i in range(len(Kmatrix)):
            for j in range(i, len(Kmatrix)):
                Kmatrix[i][j] /= np.sqrt(Kmatrix_diag[i] * Kmatrix_diag[j])
                Kmatrix[j][i] = Kmatrix[i][j]

            print()
        if params_out == {}:
                print('the gram matrix is: ')
            str_fw += 'the gram matrix is:\n\n'
        else:
                print('the gram matrix with parameters', params_out, 'is: ')
            str_fw += 'the gram matrix with parameters %s is:\n\n' % params_out
        if len(Kmatrix) < 2:
            nb_gm_ignore += 1
                print('ignored, as at most only one of all its diagonal value is non-zero.')
            str_fw += 'ignored, as at most only one of all its diagonal value is non-zero.\n\n'
        else:                
            if np.isnan(Kmatrix).any(
            ):  # if the matrix contains elements that are not numbers
                    print('ignored, as it contains elements that are not numbers.')
                str_fw += 'ignored, as it contains elements that are not numbers.\n\n'
            else:
                str_fw += np.array2string(
                        Kmatrix,
                        separator=',') + '\n\n'

                fig_file_name = results_dir + '/GM[ds]' + ds_name
                if params_out != {}:
                    fig_file_name += '[params]' + str(idx)
                plt.imshow(Kmatrix)
                plt.colorbar()
                plt.savefig(fig_file_name + '.eps', format='eps', dpi=300)
                plt.clf()
                gram_matrices.append(Kmatrix)
                gram_matrix_time.append(current_run_time)
                param_list_pre_revised.append(params_out)
                if nb_g_ignore > 0:
                        print(', where %d graphs are ignored as their graph kernels with themselves are zeros.' % nb_g_ignore)
                    str_fw += ', where %d graphs are ignored as their graph kernels with themselves are zeros.' % nb_g_ignore
        print()
        print(
    str_fw += '{} gram matrices are calculated, {} of which are ignored.\n\n'.format(len(param_list_precomputed), nb_gm_ignore)
    str_fw += 'serial numbers of gram matrix figures and their corresponding parameters settings:\n\n'
    str_fw += ''.join([
        '{}: {}\n'.format(idx, params_out)
        for idx, params_out in enumerate(param_list_precomputed)
    ])
            
    return gram_matrices, gram_matrix_time, param_list_pre_revised, y, str_fw

def read_gram_matrices_from_file(results_dir, ds_name):
    gmfile = np.load(results_dir + '/' + ds_name + '.gm.npz')
    gram_matrices = gmfile['gms'] # a list to store gram matrices for all param_grid_precomputed
    param_list_pre_revised = gmfile['params'] # list to store param grids precomputed ignoring the useless ones
    y = gmfile['y'].tolist()
    return gram_matrices, param_list_pre_revised, y

def printResultsInTable(param_list, param_list_pre_revised, average_val_scores,
                        std_val_scores, average_perf_scores, std_perf_scores,
                        average_train_scores, std_train_scores, gram_matrix_time,
                        model_type, verbose):
    from tabulate import tabulate
    table_dict = {}
    if model_type == 'regression':
        for param_in in param_list:
            param_in['alpha'] = '{:.2e}'.format(param_in['alpha'])
    else:
            param_in['C'] = '{:.2e}'.format(param_in['C'])
    table_dict['params'] = [{**param_out, **param_in}
                            for param_in in param_list for param_out in param_list_pre_revised]
    table_dict['gram_matrix_time'] = [
        '{:.2f}'.format(gram_matrix_time[index_out])
        for param_in in param_list
        for index_out, _ in enumerate(param_list_pre_revised)
    ]
    table_dict['valid_perf'] = [
        '{:.2f}±{:.2f}'.format(average_val_scores[index_out][index_in],
                               std_val_scores[index_out][index_in])
        for index_in, _ in enumerate(param_list)
    ]
    table_dict['test_perf'] = [
        '{:.2f}±{:.2f}'.format(average_perf_scores[index_out][index_in],
                               std_perf_scores[index_out][index_in])
    ]
    table_dict['train_perf'] = [
        '{:.2f}±{:.2f}'.format(average_train_scores[index_out][index_in],
                               std_train_scores[index_out][index_in])
    ]
    
    keyorder = [
        'params', 'train_perf', 'valid_perf', 'test_perf',
        'gram_matrix_time'
    ]
        print()
    tb_print = tabulate(OrderedDict(sorted(table_dict.items(), 
                        key=lambda i: keyorder.index(i[0]))), headers='keys')
#            print(tb_print)
    return 'table of performance v.s. hyper-params:\n\n%s\n\n' % tb_print

"""
Created on Tue Dec 11 11:39:46 2018
Parallel aid functions.
"""
import multiprocessing
from multiprocessing import Pool
import sys

def parallel_me(func, func_assign, var_to_assign, itr, len_itr=None, init_worker=None, 
                glbv=None, method=None, n_jobs=None, chunksize=None, itr_desc='',
                verbose=True):
    '''
    '''
    if method == 'imap_unordered':
        if glbv: # global varibles required.
#            def init_worker(v_share):
#                global G_var
#                G_var = v_share
            if n_jobs == None:
                n_jobs = multiprocessing.cpu_count()
            with Pool(processes=n_jobs, initializer=init_worker, 
                      initargs=glbv) as pool:                
                if chunksize == None:
                    if len_itr < 100 * n_jobs:
                        chunksize = int(len_itr / n_jobs) + 1
                    else:
                        chunksize = 100
                for result in (tqdm(pool.imap_unordered(func, itr, chunksize),
                    desc=itr_desc, file=sys.stdout) if verbose else 
                    pool.imap_unordered(func, itr, chunksize)):
                    func_assign(result, var_to_assign)
        else:
            with Pool(processes=n_jobs) as pool:
                    else:
            
def parallel_gm(func, Kmatrix, Gn, init_worker=None, glbv=None, 
                method='imap_unordered', n_jobs=None, chunksize=None, 
    from itertools import combinations_with_replacement
    def func_assign(result, var_to_assign):
        var_to_assign[result[0]][result[1]] = result[2]
        var_to_assign[result[1]][result[0]] = result[2]
    itr = combinations_with_replacement(range(0, len(Gn)), 2)
    len_itr = int(len(Gn) * (len(Gn) + 1) / 2)
    parallel_me(func, func_assign, Kmatrix, itr, len_itr=len_itr,
        init_worker=init_worker, glbv=glbv, method=method, n_jobs=n_jobs, 
        chunksize=chunksize, itr_desc='calculating kernels', verbose=verbose)

"""
Created on Wed Jan 30 10:48:49 2019

Trie (prefix tree)
@references:
        https://viblo.asia/p/nlp-build-a-trie-data-structure-from-scratch-with-python-3P0lPzroKox, 2019.1
"""

import pickle
import json

""" Trie class
"""
class Trie:
    # init Trie class
        self.root = self.getNode()

    def getNode(self):
        return {"isEndOfWord": False, "children": {}}

    def insertWord(self, word):
        current = self.root
        for ch in word:

            if ch in current["children"]:
                node = current["children"][ch]
            else:
                node = self.getNode()
                current["children"][ch] = node

            current = node
        current["isEndOfWord"] = True
        if 'count' in current:
            current['count'] += 1
        else:
            current['count'] = 1

    def searchWord(self, word):
            if ch not in current["children"]:
                return 0
            node = current["children"][ch]

            return current["count"]
        else:
            return 0

    def searchWordPrefix(self, word):
            if not current["children"].has_key(ch):
                return False

        # return True if children contain keys and values
        return bool(current["children"])

    def deleteWord(self, word):
        self._delete(self.root, word, 0)

    def _delete(self, current, word, index):
        if(index == len(word)):
            if not current["isEndOfWord"]:
            current["isEndOfWord"] = False
            return len(current["children"].keys()) == 0

        ch = word[index]
        if not current["children"].has_key(ch):
            return False
        node = current["children"][ch]

        should_delete_current_node = self._delete(node, word, index + 1)

        if should_delete_current_node:
            current["children"].pop(ch)

        return False

    def save_to_pickle(self, file_name):
        f = open(file_name + ".pkl", "wb")
        pickle.dump(self.root, f)
        f.close()

    def load_from_pickle(self, file_name):
        f = open(file_name + ".pkl", "rb")
        self.root = pickle.load(f)
        f.close()
        
    def to_json(self):
        return json.dump(self.root)

    def save_to_json(self, file_name):
        json_data = json.dumps(self.root)
        f = open(file_name + ".json", "w")
        f.write(json_data)
        f.close()

    def load_from_json(self, file_name):
        json_file = open(file_name + ".json", "r")
        self.root = json.load(json_file)
        json_file.close()

import networkx as nx
from copy import deepcopy
#from itertools import product

# from tqdm import tqdm

def getSPLengths(G1):
    sp = nx.shortest_path(G1)
    distances = np.zeros((G1.number_of_nodes(), G1.number_of_nodes()))
    for i in sp.keys():
        for j in sp[i].keys():
            distances[i, j] = len(sp[i][j]) - 1
    return distances

def getSPGraph(G, edge_weight=None):
    """Transform graph G to its corresponding shortest-paths graph.

    G : NetworkX graph
        The graph to be tramsformed.
    edge_weight : string
        edge attribute corresponding to the edge weight.

    Return
    ------
    S : NetworkX graph
        The shortest-paths graph corresponding to G.

    Notes
    ------
    For an input graph G, its corresponding shortest-paths graph S contains the same set of nodes as G, while there exists an edge between all nodes in S which are connected by a walk in G. Every edge in S between two nodes is labeled by the shortest distance between these two nodes.

    [1] Borgwardt KM, Kriegel HP. Shortest-path kernels on graphs. InData Mining, Fifth IEEE International Conference on 2005 Nov 27 (pp. 8-pp). IEEE.
    """
    return floydTransformation(G, edge_weight=edge_weight)

def floydTransformation(G, edge_weight=None):
    """Transform graph G to its corresponding shortest-paths graph using Floyd-transformation.

        edge attribute corresponding to the edge weight. The default edge weight is bond_type.

    Return
    ------

    """
    spMatrix = nx.floyd_warshall_numpy(G, weight=edge_weight)
    S = nx.Graph()
    S.add_nodes_from(G.nodes(data=True))
    ns = list(G.nodes())
    for i in range(0, G.number_of_nodes()):
        for j in range(i + 1, G.number_of_nodes()):
            if spMatrix[i, j] != np.inf:
                S.add_edge(ns[i], ns[j], cost=spMatrix[i, j])
    return S

def untotterTransformation(G, node_label, edge_label):
    """Transform graph G according to Mahé et al.'s work to filter out tottering patterns of marginalized kernel and tree pattern kernel.

        node attribute used as label. The default node label is 'atom'.
        edge attribute used as label. The default edge label is 'bond_type'.

    Return
    ------
    gt : NetworkX graph
        The transformed graph corresponding to G.

    [1] Pierre Mahé, Nobuhisa Ueda, Tatsuya Akutsu, Jean-Luc Perret, and Jean-Philippe Vert. Extensions of marginalized graph kernels. In Proceedings of the twenty-first international conference on Machine learning, page 70. ACM, 2004.
    """
    # arrange all graphs in a list
    G = G.to_directed()
    gt = nx.Graph()
    gt.graph = G.graph
    gt.add_nodes_from(G.nodes(data=True))
    for edge in G.edges():
        gt.add_node(edge)
        gt.node[edge].update({node_label: G.node[edge[1]][node_label]})
        gt.add_edge(edge[0], edge)
        gt.edges[edge[0], edge].update({
            edge_label:
            G[edge[0]][edge[1]][edge_label]
        })
        for neighbor in G[edge[1]]:
            if neighbor != edge[0]:
                gt.add_edge(edge, (edge[1], neighbor))
                gt.edges[edge, (edge[1], neighbor)].update({
                    edge_label:
                    G[edge[1]][neighbor][edge_label]
                })
    # nx.draw_networkx(gt)
    # plt.show()

    # relabel nodes using consecutive integers for convenience of kernel calculation.
    gt = nx.convert_node_labels_to_integers(
        gt, first_label=0, label_attribute='label_orignal')
    return gt

def direct_product(G1, G2, node_label, edge_label):
    """Return the direct/tensor product of directed graphs G1 and G2.

    G1, G2 : NetworkX graph
        The original graphs.

    Return
    ------
        The direct product graph of G1 and G2.

    Notes
    -----
    This method differs from networkx.tensor_product in that this method only adds nodes and edges in G1 and G2 that have the same labels to the direct product graph.

    [1] Thomas Gärtner, Peter Flach, and Stefan Wrobel. On graph kernels: Hardness results and efficient alternatives. Learning Theory and Kernel Machines, pages 129–143, 2003.
    """
    from itertools import product
    # G = G.to_directed()
    gt = nx.DiGraph()
    # add nodes
    for u, v in product(G1, G2):
        if G1.nodes[u][node_label] == G2.nodes[v][node_label]:
            gt.add_node((u, v))
            gt.nodes[(u, v)].update({node_label: G1.nodes[u][node_label]})
    # add edges, faster for sparse graphs (no so many edges), which is the most case for now.
    for (u1, v1), (u2, v2) in product(G1.edges, G2.edges):
        if (u1, u2) in gt and (
                v1, v2
        ) in gt and G1.edges[u1, v1][edge_label] == G2.edges[u2,
                                                             v2][edge_label]:
            gt.add_edge((u1, u2), (v1, v2))
            gt.edges[(u1, u2), (v1, v2)].update({
                edge_label:
                G1.edges[u1, v1][edge_label]
            })

    # # add edges, faster for dense graphs (a lot of edges, complete graph would be super).
    # for u, v in product(gt, gt):
    #     if (u[0], v[0]) in G1.edges and (
    #             u[1], v[1]
    #     ) in G2.edges and G1.edges[u[0],
    #                                v[0]][edge_label] == G2.edges[u[1],
    #                                                              v[1]][edge_label]:
    #         gt.add_edge((u[0], u[1]), (v[0], v[1]))
    #         gt.edges[(u[0], u[1]), (v[0], v[1])].update({
    #             edge_label:
    #             G1.edges[u[0], v[0]][edge_label]
    #         })

    # gt = nx.convert_node_labels_to_integers(
    #     gt, first_label=0, label_attribute='label_orignal')
    return gt

def graph_deepcopy(G):
    """Deep copy a graph, including deep copy of all nodes, edges and 
    attributes of the graph, nodes and edges.
    
    Note
    ----
    It is the same as the NetworkX function graph.copy(), as far as I know.
    """
    # add graph attributes.
    labels = {}
    for k, v in G.graph.items():
        labels[k] = deepcopy(v)
    if G.is_directed():
        G_copy = nx.DiGraph(**labels)
    else:
        G_copy = nx.Graph(**labels)
        
    # add nodes    
    for nd, attrs in G.nodes(data=True):
        labels = {}
        for k, v in attrs.items():
            labels[k] = deepcopy(v)
        G_copy.add_node(nd, **labels)
        
    # add edges.
    for nd1, nd2, attrs in G.edges(data=True):
        labels = {}
        G_copy.add_edge(nd1, nd2, **labels)
    
    return G_copy

def graph_isIdentical(G1, G2):
    """Check if two graphs are identical, including: same nodes, edges, node
    labels/attributes, edge labels/attributes.
    
    Notes
    ----
    1. The type of graphs has to be the same.
    2. Global/Graph attributes are neglected as they may contain names for graphs.
    """
    # check nodes.
    nlist1 = [n for n in G1.nodes(data=True)]
    nlist2 = [n for n in G2.nodes(data=True)]
    if not nlist1 == nlist2:
    # check edges.
    elist1 = [n for n in G1.edges(data=True)]
    elist2 = [n for n in G2.edges(data=True)]
    if not elist1 == elist2:
    # check graph attributes.
    
    return True

def get_node_labels(Gn, node_label):
    """Get node labels of dataset Gn.
    """
    nl = set()
    for G in Gn:
        nl = nl | set(nx.get_node_attributes(G, node_label).values())
    return nl

def get_edge_labels(Gn, edge_label):
    """Get edge labels of dataset Gn.
    """
    el = set()
    for G in Gn:
        el = el | set(nx.get_edge_attributes(G, edge_label).values())
    return el