Python networkx.cycle_basis() Examples

def decode_graph(adj, prefix):
    adj = np.asmatrix(adj)
    G = nx.from_numpy_matrix(adj)
    # G.remove_nodes_from(nx.isolates(G))
    print('num of nodes: {}'.format(G.number_of_nodes()))
    print('num of edges: {}'.format(G.number_of_edges()))
    G_deg = nx.degree_histogram(G)
    G_deg_sum = [a * b for a, b in zip(G_deg, range(0, len(G_deg)))]
    print('average degree: {}'.format(sum(G_deg_sum) / G.number_of_nodes()))
    if nx.is_connected(G):
        print('average path length: {}'.format(nx.average_shortest_path_length(G)))
        print('average diameter: {}'.format(nx.diameter(G)))
    G_cluster = sorted(list(nx.clustering(G).values()))
    print('average clustering coefficient: {}'.format(sum(G_cluster) / len(G_cluster)))
    cycle_len = []
    cycle_all = nx.cycle_basis(G, 0)
    for item in cycle_all:
        cycle_len.append(len(item))
    print('cycles', cycle_len)
    print('cycle count', len(cycle_len))
    draw_graph(G, prefix=prefix) 

def test_cycle_basis(self):
        G = self.G
        cy = networkx.cycle_basis(G, 0)
        sort_cy = sorted(sorted(c) for c in cy)
        assert_equal(sort_cy, [[0, 1, 2, 3], [0, 1, 6, 7, 8], [0, 3, 4, 5]])
        cy = networkx.cycle_basis(G, 1)
        sort_cy = sorted(sorted(c) for c in cy)
        assert_equal(sort_cy, [[0, 1, 2, 3], [0, 1, 6, 7, 8], [0, 3, 4, 5]])
        cy = networkx.cycle_basis(G, 9)
        sort_cy = sorted(sorted(c) for c in cy)
        assert_equal(sort_cy, [[0, 1, 2, 3], [0, 1, 6, 7, 8], [0, 3, 4, 5]])
        # test disconnected graphs
        nx.add_cycle(G, "ABC")
        cy = networkx.cycle_basis(G, 9)
        sort_cy = sorted(sorted(c) for c in cy[:-1]) + [sorted(cy[-1])]
        assert_equal(sort_cy, [[0, 1, 2, 3], [0, 1, 6, 7, 8], [0, 3, 4, 5],
                     ['A', 'B', 'C']]) 

def decode_graph(adj, prefix):
    adj = np.asmatrix(adj)
    G = nx.from_numpy_matrix(adj)
    # G.remove_nodes_from(nx.isolates(G))
    print('num of nodes: {}'.format(G.number_of_nodes()))
    print('num of edges: {}'.format(G.number_of_edges()))
    G_deg = nx.degree_histogram(G)
    G_deg_sum = [a * b for a, b in zip(G_deg, range(0, len(G_deg)))]
    print('average degree: {}'.format(sum(G_deg_sum) / G.number_of_nodes()))
    if nx.is_connected(G):
        print('average path length: {}'.format(nx.average_shortest_path_length(G)))
        print('average diameter: {}'.format(nx.diameter(G)))
    G_cluster = sorted(list(nx.clustering(G).values()))
    print('average clustering coefficient: {}'.format(sum(G_cluster) / len(G_cluster)))
    cycle_len = []
    cycle_all = nx.cycle_basis(G, 0)
    for item in cycle_all:
        cycle_len.append(len(item))
    print('cycles', cycle_len)
    print('cycle count', len(cycle_len))
    draw_graph(G, prefix=prefix) 

def stabilizers(fer_op):
    """ stabilizers """
    edge_list = bravyi_kitaev_fast_edge_list(fer_op)
    num_qubits = edge_list.shape[1]

    graph = networkx.Graph()
    graph.add_edges_from(tuple(edge_list.transpose()))
    stabs = np.asarray(networkx.cycle_basis(graph))
    stabilizer_ops = []
    for stab in stabs:
        a_op = WeightedPauliOperator(paulis=[[1.0, Pauli.from_label('I' * num_qubits)]])
        stab = np.asarray(stab)
        for i in range(np.size(stab)):
            a_op = a_op * edge_operator_aij(edge_list, stab[i], stab[(i + 1) % np.size(stab)]) * 1j
        stabilizer_ops.append(a_op)

    return stabilizer_ops 

def test_cycle_basis(self):
        G = self.G
        cy = networkx.cycle_basis(G, 0)
        sort_cy = sorted(sorted(c) for c in cy)
        assert_equal(sort_cy, [[0, 1, 2, 3], [0, 1, 6, 7, 8], [0, 3, 4, 5]])
        cy = networkx.cycle_basis(G, 1)
        sort_cy = sorted(sorted(c) for c in cy)
        assert_equal(sort_cy, [[0, 1, 2, 3], [0, 1, 6, 7, 8], [0, 3, 4, 5]])
        cy = networkx.cycle_basis(G, 9)
        sort_cy = sorted(sorted(c) for c in cy)
        assert_equal(sort_cy, [[0, 1, 2, 3], [0, 1, 6, 7, 8], [0, 3, 4, 5]])
        # test disconnected graphs
        nx.add_cycle(G, "ABC")
        cy = networkx.cycle_basis(G, 9)
        sort_cy = sorted(sorted(c) for c in cy[:-1]) + [sorted(cy[-1])]
        assert_equal(sort_cy, [[0, 1, 2, 3], [0, 1, 6, 7, 8], [0, 3, 4, 5],
                               ['A', 'B', 'C']]) 

def test_cycle_basis(self):
        G=self.G
        cy=networkx.cycle_basis(G,0)
        sort_cy= sorted( sorted(c) for c in cy )
        assert_equal(sort_cy, [[0,1,2,3],[0,1,6,7,8],[0,3,4,5]])
        cy=networkx.cycle_basis(G,1)
        sort_cy= sorted( sorted(c) for c in cy )
        assert_equal(sort_cy, [[0,1,2,3],[0,1,6,7,8],[0,3,4,5]])
        cy=networkx.cycle_basis(G,9)
        sort_cy= sorted( sorted(c) for c in cy )
        assert_equal(sort_cy, [[0,1,2,3],[0,1,6,7,8],[0,3,4,5]])
        # test disconnected graphs
        G.add_cycle(list("ABC"))
        cy=networkx.cycle_basis(G,9)
        sort_cy= sorted(sorted(c) for c in cy[:-1]) + [sorted(cy[-1])]
        assert_equal(sort_cy, [[0,1,2,3],[0,1,6,7,8],[0,3,4,5],['A','B','C']]) 

def loops_within_feeder(self, *args):
        """Returns the number of loops within a feeder."""
        if args:
            return len(nx.cycle_basis(args[0]))
        else:
            return len(nx.cycle_basis(self.G.graph)) 

def test_cycle_basis(self):
        G = nx.MultiGraph()
        cy = networkx.cycle_basis(G, 0) 

def test_cycle_basis(self):
        G = nx.DiGraph()
        cy = networkx.cycle_basis(G, 0) 

def minimum_cycle_basis(G, weight=None):
    """ Returns a minimum weight cycle basis for G

    Minimum weight means a cycle basis for which the total weight
    (length for unweighted graphs) of all the cycles is minimum.

    Parameters
    ----------
    G : NetworkX Graph
    weight: string
        name of the edge attribute to use for edge weights

    Returns
    -------
    A list of cycle lists.  Each cycle list is a list of nodes
    which forms a cycle (loop) in G. Note that the nodes are not
    necessarily returned in a order by which they appear in the cycle

    Examples
    --------
    >>> G=nx.Graph()
    >>> G.add_cycle([0,1,2,3])
    >>> G.add_cycle([0,3,4,5])
    >>> print(nx.minimum_cycle_basis(G))
    [[0, 1, 2, 3], [0, 3, 4, 5]]

    References:
        [1] Kavitha, Telikepalli, et al. "An O(m^2n) Algorithm for
        Minimum Cycle Basis of Graphs."
        http://link.springer.com/article/10.1007/s00453-007-9064-z
        [2] de Pina, J. 1995. Applications of shortest path methods.
        Ph.D. thesis, University of Amsterdam, Netherlands

    See Also
    --------
    simple_cycles, cycle_basis
    """
    # We first split the graph in commected subgraphs
    return sum((_min_cycle_basis(c, weight) for c in
                nx.connected_component_subgraphs(G)), []) 

def test_cycle_basis(self):
        G = nx.MultiGraph()
        cy = networkx.cycle_basis(G, 0) 

def test_cycle_basis(self):
        G = nx.DiGraph()
        cy = networkx.cycle_basis(G, 0) 

def penalized_logp(s):
    if s is None: return -100.0
    mol = Chem.MolFromSmiles(s)
    if mol is None: return -100.0

    logP_mean = 2.4570953396190123
    logP_std = 1.434324401111988
    SA_mean = -3.0525811293166134
    SA_std = 0.8335207024513095
    cycle_mean = -0.0485696876403053
    cycle_std = 0.2860212110245455

    log_p = Descriptors.MolLogP(mol)
    SA = -sascorer.calculateScore(mol)

    # cycle score
    cycle_list = nx.cycle_basis(nx.Graph(Chem.rdmolops.GetAdjacencyMatrix(mol)))
    if len(cycle_list) == 0:
        cycle_length = 0
    else:
        cycle_length = max([len(j) for j in cycle_list])
    if cycle_length <= 6:
        cycle_length = 0
    else:
        cycle_length = cycle_length - 6
    cycle_score = -cycle_length

    normalized_log_p = (log_p - logP_mean) / logP_std
    normalized_SA = (SA - SA_mean) / SA_std
    normalized_cycle = (cycle_score - cycle_mean) / cycle_std
    return normalized_log_p + normalized_SA + normalized_cycle 

def test_cycle_basis(self):
        G=nx.MultiGraph()
        cy=networkx.cycle_basis(G,0) 

def test_cycle_basis(self):
        G=nx.DiGraph()
        cy=networkx.cycle_basis(G,0) 

def get_graph_metrics(self):
        nxGraph = self.nxGraph.to_undirected()
        islands = list(nx.connected_component_subgraphs(nxGraph))
        print('Number of islands: ', len(islands))
        loops = nx.cycle_basis(nxGraph)
        print('Number of loops: ', len(loops))
        print('')
        for i, eachCycle in enumerate(loops):
            print('Loop ' + str(i) + ' : starts with node "' + eachCycle[0] + '"') 

def penalized_logp(s):
    if s is None: return -100.0
    mol = Chem.MolFromSmiles(s)
    if mol is None: return -100.0

    logP_mean = 2.4570953396190123
    logP_std = 1.434324401111988
    SA_mean = -3.0525811293166134
    SA_std = 0.8335207024513095
    cycle_mean = -0.0485696876403053
    cycle_std = 0.2860212110245455

    log_p = Descriptors.MolLogP(mol)
    SA = -sascorer.calculateScore(mol)

    # cycle score
    cycle_list = nx.cycle_basis(nx.Graph(Chem.rdmolops.GetAdjacencyMatrix(mol)))
    if len(cycle_list) == 0:
        cycle_length = 0
    else:
        cycle_length = max([len(j) for j in cycle_list])
    if cycle_length <= 6:
        cycle_length = 0
    else:
        cycle_length = cycle_length - 6
    cycle_score = -cycle_length

    normalized_log_p = (log_p - logP_mean) / logP_std
    normalized_SA = (SA - SA_mean) / SA_std
    normalized_cycle = (cycle_score - cycle_mean) / cycle_std
    return normalized_log_p + normalized_SA + normalized_cycle 

def correctLoops(code, loops_graph, charge_type):
    while nx.cycle_basis(loops_graph) != []:
        cycle = nx.cycle_basis(loops_graph)[0]
        loop = path.Path(cycle)
        for data in code.Primal.nodes():
            if loop.contains_points([data]) == [True]:
                charge = code.Primal.node[data]['charge'][charge_type]
                code.Primal.node[data]['charge'][charge_type] = (charge + 1)%2

        l = len(cycle)
        for i in range(l):
            n1, n2 = cycle[i], cycle[(i+1)%l]
            loops_graph.remove_edge(*(n1,n2))

    return code, loops_graph 

def cleanup_fixed_graph(pruned):
    # remove cycles
    while True:
        cycles = networkx.cycle_basis(pruned)
        if len(cycles) == 0:
            break

        to_delete = None
        worst_val = None
        for cycle in cycles:
            cur_worst_val = None
            cur_worst = None
            for a,b,data in pruned.edges(cycle, data=True):
                cur_val = data["p"]
                if cur_worst_val is None or cur_val < cur_worst_val:
                    cur_worst_val = cur_val
                    cur_worst = (a,b)
            if worst_val is None or cur_worst_val < worst_val:
                worst_val = cur_worst_val
                to_delete = cur_worst

        pruned.remove_edge(*to_delete)

    # remove all cis-edges at the ends of subgraphs
    degrees = pruned.degree()
    to_delete = []
    for node, degree in dict(degrees).items():
        if degree == 1:
            edge = list(pruned.edges([node], data=True))[0]
            if edge[2]["kind"]=="facing":
                to_delete.append(node)
    pruned.remove_nodes_from(to_delete)

    # remove unconnected nodes
    pruned.remove_nodes_from([node for (node, degree) in dict(pruned.degree()).items() if degree==0])

    return pruned 

def vacuum_operator(fer_op):
    """Use the stabilizers to find the vacuum state in bravyi_kitaev_fast.

    Args:
        fer_op (FermionicOperator): the fermionic operator in the second quantized form

    Returns:
        WeightedPauliOperator: the qubit operator
    """
    edge_list = bravyi_kitaev_fast_edge_list(fer_op)
    num_qubits = edge_list.shape[1]
    vac_operator = WeightedPauliOperator(paulis=[[1.0, Pauli.from_label('I' * num_qubits)]])

    graph = networkx.Graph()
    graph.add_edges_from(tuple(edge_list.transpose()))
    stabs = np.asarray(networkx.cycle_basis(graph))
    for stab in stabs:
        a_op = WeightedPauliOperator(paulis=[[1.0, Pauli.from_label('I' * num_qubits)]])
        stab = np.asarray(stab)
        for i in range(np.size(stab)):
            a_op = a_op * edge_operator_aij(edge_list, stab[i], stab[(i + 1) % np.size(stab)]) * 1j
            # a_op.scaling_coeff(1j)
        a_op += WeightedPauliOperator(paulis=[[1.0, Pauli.from_label('I' * num_qubits)]])
        vac_operator = vac_operator * a_op * np.sqrt(2)
        # vac_operator.scaling_coeff()

    return vac_operator 

def enumerate_cycle_basis(molecule):
        """Enumerate a closed cycle basis of bonds in molecule.

        This uses cycle_basis from NetworkX:
        https://networkx.github.io/documentation/networkx-1.10/reference/generated/networkx.algorithms.cycles.cycle_basis.html#networkx.algorithms.cycles.cycle_basis

        Parameters
        ----------
        molecule : OEMol
            The molecule for a closed cycle basis of Bonds is to be identified

        Returns
        -------
        bond_cycle_basis : list of list of OEBond
            bond_cycle_basis[cycle_index] is a list of OEBond objects that define a cycle in the basis
            You can think of these as the minimal spanning set of ring systems to check.
        """
        g = nx.Graph()
        for atom in molecule.GetAtoms():
            g.add_node(atom.GetIdx())
        for bond in molecule.GetBonds():
            g.add_edge(bond.GetBgnIdx(), bond.GetEndIdx(), bond=bond)
        bond_cycle_basis = list()
        for cycle in nx.cycle_basis(g):
            bond_cycle = list()
            for i in range(len(cycle)):
                atom_index_1 = cycle[i]
                atom_index_2 = cycle[(i+1)%len(cycle)]
                edge = g.edges[atom_index_1,atom_index_2]
                bond = edge['bond']
                bond_cycle.append(bond)
            bond_cycle_basis.append(bond_cycle)
        return bond_cycle_basis 

def vacuum_operator(edge_matrix_indices):
    """Use the stabilizers to find the vacuum state in bravyi_kitaev_fast.

    Args:
        edge_matrix_indices(numpy array): specifying the edges

    Return:
        An instance of QubitOperator

    """
    # Initialize qubit operator.
    g = networkx.Graph()
    g.add_edges_from(tuple(edge_matrix_indices.transpose()))
    stabs = numpy.array(networkx.cycle_basis(g))
    vac_operator = QubitOperator(())
    for stab in stabs:

        A = QubitOperator(())
        stab = numpy.array(stab)
        for i in range(numpy.size(stab)):
            if i == (numpy.size(stab) - 1):
                A = (1j)*A * edge_operator_aij(edge_matrix_indices,
                                               stab[i], stab[0])
            else:
                A = (1j)*A * edge_operator_aij(edge_matrix_indices,
                                               stab[i], stab[i+1])
        vac_operator = vac_operator*(QubitOperator(()) + A)/numpy.sqrt(2)

    return vac_operator 

def cycles_count(variables, relations):

    g = as_networkx_graph(variables, relations)
    cycles = nx.cycle_basis(g)
    return len(cycles) 

def find_cycles(sub_network, weight='x_pu'):
    """
    Find all cycles in the sub_network and record them in sub_network.C.

    networkx collects the cycles with more than 2 edges; then the 2-edge cycles
    from the MultiGraph must be collected separately (for cases where there
    are multiple lines between the same pairs of buses).

    Cycles with infinite impedance are skipped.
    """
    branches_bus0 = sub_network.branches()["bus0"]
    branches_i = branches_bus0.index

    #reduce to a non-multi-graph for cycles with > 2 edges
    mgraph = sub_network.graph(weight=weight, inf_weight=False)
    graph = nx.OrderedGraph(mgraph)

    cycles = nx.cycle_basis(graph)

    #number of 2-edge cycles
    num_multi = len(mgraph.edges()) - len(graph.edges())

    sub_network.C = dok_matrix((len(branches_bus0),len(cycles)+num_multi))

    for j,cycle in enumerate(cycles):

        for i in range(len(cycle)):
            branch = next(iterkeys(mgraph[cycle[i]][cycle[(i+1)%len(cycle)]]))
            branch_i = branches_i.get_loc(branch)
            sign = +1 if branches_bus0.iat[branch_i] == cycle[i] else -1
            sub_network.C[branch_i,j] += sign

    #counter for multis
    c = len(cycles)

    #add multi-graph 2-edge cycles for multiple branches between same pairs of buses
    for u,v in graph.edges():
        bs = list(mgraph[u][v].keys())
        if len(bs) > 1:
            first = bs[0]
            first_i = branches_i.get_loc(first)
            for b in bs[1:]:
                b_i = branches_i.get_loc(b)
                sign = -1 if branches_bus0.iat[b_i] == branches_bus0.iat[first_i] else +1
                sub_network.C[first_i,c] = 1
                sub_network.C[b_i,c] = sign
                c+=1 

def identify_possible_highly_variable(G,
                                      cycle_threshold_max=20,
                                      cycle_threshold_min=5,
                                      size_diff_threshold=0.5):

    # add family paralog attribute to nodes
    for node in G.nodes():
        G.nodes[node]['highVar'] = 0

    # find all the cycles shorter than cycle_threshold
    complete_basis = []
    for c in nx.connected_components(G):
        sub_G = G.subgraph(c)
        basis = nx.cycle_basis(sub_G, list(sub_G.nodes())[0])
        complete_basis += [
            set(b) for b in basis if len(b) <= cycle_threshold_max
        ]

    # remove cycles that are too short
    complete_basis = [b for b in complete_basis if len(b) >= 3]

    # merge cycles with more than one node in common (nested)
    if len(complete_basis) < 1:
        return G

    merged_basis = [[1, set(complete_basis[0])]]
    for b in complete_basis[1:]:
        b = set(b)
        merged = False
        for i, mb in enumerate(merged_basis):
            if len(mb[1].intersection(b)) > 1:
                merged = True
                merged_basis[i][0] += 1
                merged_basis[i][1] |= b
        if not merged:
            merged_basis.append([1, b])

    for b in merged_basis:
        if b[0] < cycle_threshold_min: continue
        max_size = max([G.nodes[node]['size'] for node in b[1]])
        for node in b[1]:
            if G.nodes[node]['size'] < (size_diff_threshold * max_size):
                G.nodes[node]['highVar'] = 1

    return G 

def find_loops(self, edge_node_df=None):
        """
        This function converts the input matrix into a networkx type graph and identifies all fundamental loops
        of the network. The group of fundamental loops is defined as the series of linear independent loops which
        can be combined to form all other loops.

        :param pd.DataFrame edge_node_df: DataFrame consisting of n rows (number of nodes) and e columns (number of edges)
                            and indicating the direction of flow of each edge e at node n: if e points to n,
                            value is 1; if e leaves node n, -1; else, 0. E.g. a plant will only have exiting flows,
                            so only negative values                               (n x e)

        :return: loops: list of all fundamental loops in the network
        :return: graph: networkx dictionary type graph of network

        :rtype: loops: list
        :rtype: graph: dictionary
        """
        if edge_node_df is None:
            edge_node_df = self.edge_node_df

        # check to see if we've already computed the loops
        edge_node_str = edge_node_df.to_string()
        if edge_node_str in ThermalNetwork.loops_cache:
            return ThermalNetwork.loops_cache[edge_node_str]

        edge_node_df_t = np.transpose(edge_node_df)  # transpose matrix to more intuitively setup graph

        graph = nx.Graph()  # set up networkx type graph

        for i in range(edge_node_df_t.shape[0]):
            new_edge = [0, 0]
            for j in range(0, edge_node_df_t.shape[1]):
                if edge_node_df_t.iloc[i][edge_node_df_t.columns[j]] == 1:
                    new_edge[0] = j
                elif edge_node_df_t.iloc[i][edge_node_df_t.columns[j]] == -1:
                    new_edge[1] = j
            graph.add_edge(new_edge[0], new_edge[1], edge_number=i)  # add edges to graph
            # edge number necessary to later identify which edges are in loop since graph is a dictionary

        loops = nx.cycle_basis(graph, 0)  # identifies all linear independent loops

        ThermalNetwork.loops_cache[edge_node_str] = (loops, graph)
        return loops, graph


# collect the results of each call to hourly_thermal_calculation in a record 

def get_normalized_values():
    fname = '/home/bowen/pycharm_deployment_directory/rl_graph_generation/gym-molecule/gym_molecule/dataset/250k_rndm_zinc_drugs_clean.smi'
    with open(fname) as f:
        smiles = f.readlines()

    for i in range(len(smiles)):
        smiles[i] = smiles[i].strip()
    smiles_rdkit = []

    for i in range(len(smiles)):
        smiles_rdkit.append(Chem.MolToSmiles(Chem.MolFromSmiles(smiles[i])))
    print(i)

    logP_values = []
    for i in range(len(smiles)):
        logP_values.append(MolLogP(Chem.MolFromSmiles(smiles_rdkit[i])))
    print(i)

    SA_scores = []
    for i in range(len(smiles)):
        SA_scores.append(
            -calculateScore(Chem.MolFromSmiles(smiles_rdkit[i])))
    print(i)

    cycle_scores = []
    for i in range(len(smiles)):
        cycle_list = nx.cycle_basis(nx.Graph(
            Chem.rdmolops.GetAdjacencyMatrix(Chem.MolFromSmiles(smiles_rdkit[
                                                                  i]))))
        if len(cycle_list) == 0:
            cycle_length = 0
        else:
            cycle_length = max([len(j) for j in cycle_list])
        if cycle_length <= 6:
            cycle_length = 0
        else:
            cycle_length = cycle_length - 6
        cycle_scores.append(-cycle_length)
    print(i)

    SA_scores_normalized = (np.array(SA_scores) - np.mean(SA_scores)) / np.std(
        SA_scores)
    logP_values_normalized = (np.array(logP_values) - np.mean(
        logP_values)) / np.std(logP_values)
    cycle_scores_normalized = (np.array(cycle_scores) - np.mean(
        cycle_scores)) / np.std(cycle_scores)

    return np.mean(SA_scores), np.std(SA_scores), np.mean(
        logP_values), np.std(logP_values), np.mean(
        cycle_scores), np.std(cycle_scores)




# smile = 'C'*38 

def reward_penalized_log_p(mol):
    """
    Reward that consists of log p penalized by SA and # long cycles,
    as described in (Kusner et al. 2017). Scores are normalized based on the
    statistics of 250k_rndm_zinc_drugs_clean.smi dataset
    :param mol: rdkit mol object
    :return: float
    """
    # normalization constants, statistics from 250k_rndm_zinc_drugs_clean.smi
    logP_mean = 2.4570953396190123
    logP_std = 1.434324401111988
    SA_mean = -3.0525811293166134
    SA_std = 0.8335207024513095
    cycle_mean = -0.0485696876403053
    cycle_std = 0.2860212110245455

    log_p = MolLogP(mol)
    SA = -calculateScore(mol)

    # cycle score
    cycle_list = nx.cycle_basis(nx.Graph(
        Chem.rdmolops.GetAdjacencyMatrix(mol)))
    if len(cycle_list) == 0:
        cycle_length = 0
    else:
        cycle_length = max([len(j) for j in cycle_list])
    if cycle_length <= 6:
        cycle_length = 0
    else:
        cycle_length = cycle_length - 6
    cycle_score = -cycle_length

    normalized_log_p = (log_p - logP_mean) / logP_std
    normalized_SA = (SA - SA_mean) / SA_std
    normalized_cycle = (cycle_score - cycle_mean) / cycle_std

    return normalized_log_p + normalized_SA + normalized_cycle


# # TEST compare with junction tree paper examples from Figure 7
# assert round(reward_penalized_log_p(Chem.MolFromSmiles('ClC1=CC=C2C(C=C(C('
#                                                        'C)=O)C(C(NC3=CC(NC('
#                                                        'NC4=CC(C5=C('
#                                                        'C)C=CC=C5)=CC=C4)=O)=CC=C3)=O)=C2)=C1')), 2) == 5.30
# assert round(reward_penalized_log_p(Chem.MolFromSmiles('CC(NC1=CC(C2=CC=CC('
#                                                        'NC(NC3=CC=CC(C4=CC('
#                                                        'F)=CC=C4)=C3)=O)=C2)=CC=C1)=O')), 2) == 4.49
# assert round(reward_penalized_log_p(Chem.MolFromSmiles('ClC(C('
#                                                        'Cl)=C1)=CC=C1NC2=CC=CC=C2C(NC(NC3=C(C(NC4=C(Cl)C=CC=C4)=S)C=CC=C3)=O)=O')), 2) == 4.93