<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 NetworkX is a Python package for creation, manipulation, and study of the structure, dynamics, and functions of complex networks.
A complex network is a graph (network) elements of which are connected neither purely regular nor purely random. Some cool examples
Step1:  There are two major components in a network - nodes (or vertices) and edges connecting nodes. They are not specified as networkx objects, leaving you to use something meaningful as nodes and edges. Nodes can be anything
Step2:  Depending on the nature of your task, you may decide to use different types of graphs. There are
Step3:  Data structure
networkx uses a “dictionary of dictionaries of dictionaries” as the basic network data structure. This allows fast lookup with reasonable storage for large sparse networks. To access nodes's or edges's attributes
Step4:  To add an attribute to nodes or edges
Step5:  Analysing network
The structure of the graph can be analyzed using various graph-theoretic functions. For example, to find out how many nodes are connected together in a cluster and how many clusters there are use nx.connected_components().
Step6:  Or find the shortest path between point A and point Dusing Dijkstra’s algorithm
Step7:  Drawing network
networkx is not primarily a graph drawing package but basic drawing with matplotlib as well as an interface to use the open source graphviz software package are included. The basic drawing functions essentially place the nodes on a scatterplot using the positions you provide via a dictionary or the positions are computed with a layout function. The edges are lines between those dots.
Step8:  dot draws hierarchical layouts of directed graphs while neato creates spring model layouts (http
Step9:  Task 1
Step10:  So, unsurprisingly, the hydroxyl radical (OH) has the maximum number of interactions in the network (actually in both, the MCM and the CRI). Ozone is on 11th place in the MCM, but only 40th in the CRI.
Task 2. Find species reactions between which directly lead to ozone formation.
Step11:  Task 3. Find species that are produced directly in ozone destroying reactions.
Step12:  We were able to find the last two answers only because at the beginning we used MultiDiGraph(), which is a directed graph.
If we load the CRI mechanism and rerun the same cells, the layout of the CRI network is
<ASSISTANT_TASK:>
Python Code:
import networkx as nx
import matplotlib.pyplot as plt
%matplotlib inline
Explanation: NetworkX is a Python package for creation, manipulation, and study of the structure, dynamics, and functions of complex networks.
A complex network is a graph (network) elements of which are connected neither purely regular nor purely random. Some cool examples: 
* Lord of the Rings
* Star Wars
* Game of Thrones
All these networks have different sizes, densities, degrees of centrality, clustarisation and so on, and one of the ways to study these properties is to use networkx.
End of explanation
G = nx.Graph() # create an empty graph
G.add_node('Luke') # add one node
G.add_nodes_from(['Leia', 'Han']) # add multiple nodes
G.add_edge('Luke', 'Leia') # add one edge
G.add_edges_from([('Luke', 'Han'), ('Leia', 'Han'), ('Leia', 'Lando')]) # add multiple edges
G.remove_edge('Leia', 'Lando') # remove one edge
G.remove_node('Lando') # remove one node
nx.draw(G, with_labels=True, node_size=3000) # draw network
Explanation: There are two major components in a network - nodes (or vertices) and edges connecting nodes. They are not specified as networkx objects, leaving you to use something meaningful as nodes and edges. Nodes can be anything: strings, floats, files, etc. Edges are the relationships between nodes, so they can be directed or undirected, and can hold arbitrary data (e.g., weights).
Constructing network
The graph can be grown in several ways. The simplest is to add one or more nodes and edges at time, without using any graph generator functions. Note that networkx automatically adds a node if there is an edge containing this node.
End of explanation
G = nx.DiGraph() # create an empty directed graph
G.add_edges_from([('Luke', 'Han'), ('Luke', 'Leia'), ('Leia', 'Han')])
nx.draw(G, with_labels=True, node_size=3000)
Explanation: Depending on the nature of your task, you may decide to use different types of graphs. There are:
* Graph() - undirected graph, ignores multiple edges between two nodes
* DiGraph() - graph with directed edges
* MultiGraph() - allows multiple undirected edges between pairs of nodes
* MultiDiGraph() - directed version of a MultiGraph()
End of explanation
G['Luke']
G['Luke']['Leia']
Explanation: Data structure
networkx uses a “dictionary of dictionaries of dictionaries” as the basic network data structure. This allows fast lookup with reasonable storage for large sparse networks. To access nodes's or edges's attributes:
End of explanation
G.node['Luke']['Force'] = 1
G.node['Leia']['Force'] = 1
G.node['Han']['Force'] = 0
G.nodes(data=True)
G['Luke']['Han']['weight'] = 10
G['Luke']['Leia']['weight'] = 10
G['Leia']['Han']['weight'] = 10
G.edges(data=True)
Explanation: To add an attribute to nodes or edges:
End of explanation
G = nx.Graph()
G.add_edges_from([('Luke', 'Han'), ('Luke', 'Leia'), ('Leia', 'Han')])
G.add_node('Lando')
list(nx.connected_components(G)) # generate a sorted list of connected components, largest first
nx.draw(G, with_labels=True, node_size=3000)
Explanation: Analysing network
The structure of the graph can be analyzed using various graph-theoretic functions. For example, to find out how many nodes are connected together in a cluster and how many clusters there are use nx.connected_components().
End of explanation
G = nx.Graph()
e = [('A', 'B'), ('B', 'C'), ('A', 'C'), ('C', 'D')]
G.add_edges_from(e)
print(nx.dijkstra_path(G, 'A', 'D'))
nx.draw(G, with_labels=True)
Explanation: Or find the shortest path between point A and point Dusing Dijkstra’s algorithm:
End of explanation
G = nx.complete_graph(10) # return a graph where all nodes are connected to all edges
fig, ax = plt.subplots(1, 3, figsize=(15,4))
nx.draw_random(G, ax=ax[0])
ax[0].set_title('random')
nx.draw_circular(G, ax=ax[1])
ax[1].set_title('circular')
nx.draw_spectral(G, ax=ax[2])
ax[2].set_title('spectral');
r = 2 # branching factor of the tree
h = 4 # height of the tree
btree = nx.balanced_tree(r,h) # return the perfectly balanced r-tree of height h
pos_neato = nx.nx_pydot.graphviz_layout(btree, prog='neato') # set initial positions for nodes using neato layout
pos_dot = nx.nx_pydot.graphviz_layout(btree, prog='dot') # set initial positions for nodes using dot layout
fig, ax = plt.subplots(1, 2, figsize=(15,4))
nx.draw(btree, pos_neato, ax=ax[0])
ax[0].set_title('neato')
nx.draw(btree, pos_dot, ax=ax[1])
ax[1].set_title('dot');
Explanation: Drawing network
networkx is not primarily a graph drawing package but basic drawing with matplotlib as well as an interface to use the open source graphviz software package are included. The basic drawing functions essentially place the nodes on a scatterplot using the positions you provide via a dictionary or the positions are computed with a layout function. The edges are lines between those dots.
End of explanation
import os
import json
import numpy as np
import collections
# Choose chemical mechanism
model_name = 'MCM_C1C5' # or 'CRI_C1C5
# Load reaction descriptions
eqs_json_path = '../data/'
with open(os.path.join(eqs_json_path, model_name+'.json'), 'r') as f:
    eqs = json.load(f)
eqs[0:3]
# Convert equations's info to nodes, edges and edge labels
nodes = []
edges = []
edge_labels = collections.OrderedDict()
major_reactants = ['CL', 'H2', 'HO2', 'NO', 'NO2', 'NO3', 'OH', 'SO2', 'SO3']
for eq in eqs:
    if len(eq['reac']) == 1:
        # photolysis
        if 'J' in eq['coef'] and len(eq['prod']) == 1:
            edge = (eq['reac'][0], eq['prod'][0])
            edges.append(edge)
            edge_labels[edge] = 'hv'
        elif 'J' in eq['coef'] and len(eq['prod']) != 1:
            for prod in eq['prod']:
                edge = (eq['reac'][0], prod)
                edges.append(edge)
                edge_labels[edge] = 'hv'
        # thermal decomposition
        elif 'J' not in eq['coef'] and len(eq['prod']) == 1:
            edge = (eq['reac'][0], eq['prod'][0])
            edges.append(edge)
            edge_labels[edge] = '' 
        elif 'J' not in eq['coef'] and len(eq['prod']) != 1:
            for prod in eq['prod']:
                edge = (eq['reac'][0], prod)
                edges.append(edge)
                edge_labels[edge] = ''
    # bimolecular
    else: # len(eq['reac']) == 2:
        reac1, reac2 = eq['reac']
        if reac1 in major_reactants and reac2 not in major_reactants:
            reac1, reac2 = reac2, reac1 # for consistency in edge labels (if possible, always use a major reactant as an edge label)
        if len(eq['prod']) == 1:
            for reac in eq['reac']:
                edge = (reac, eq['prod'][0])
                edges.append(edge)
                edge_labels[edge] = reac2
        elif len(eq['prod']) != 1:
            for reac in eq['reac']:
                for prod in eq['prod']:
                    edge = (reac, prod)
                    edges.append(edge)
                    edge_labels[edge] = reac2
    # make sure that all reactants and products have a node
    for reac in eq['reac']:
        if reac not in nodes:
            nodes.append(reac)
    for prod in eq['prod']:
        if prod not in nodes:
            nodes.append(prod)
# Create network layout
scheme = nx.MultiDiGraph()
scheme.add_edges_from(edges)
scheme.add_nodes_from(nodes)
pos = nx.nx_pydot.graphviz_layout(scheme)
# Choose x, y limits for the graph
if model_name.split('_')[0] == 'MCM':
    pass
if model_name.split('_')[0] == 'CRI':
    x1, x2 = -400, 400
    y1, y2 = -400, 400
# Draw network
fig, ax = plt.subplots(figsize=(15,15))
nx.draw_networkx_nodes(scheme, pos, ax=ax, node_color='#669999')
nx.draw_networkx_edges(scheme, pos, ax=ax, edge_color='k', width=0.1, arrows=True)
# Add node and edge labels, but it makes the graph too busy
# nx.draw_networkx_labels(scheme, pos, {k: k for k in nodes}, ax=ax, font_size=9)
# nx.draw_networkx_edge_labels(scheme, pos, edge_labels=edge_labels)
_ = ax.axis('off')
ax.collections[0].set_edgecolor('k') # change color of the nodes's outline
if model_name.split('_')[0] == 'CRI':
    ax.set_xlim(x1, x2)
    ax.set_ylim(y1, y2)
Explanation: dot draws hierarchical layouts of directed graphs while neato creates spring model layouts (http://www.adp-gmbh.ch/misc/tools/graphviz/index.html).
Using networkx for analysing the Master Chemical Mechanism (MCM) and the Common Representative Intermediates (CRI) Mechanism
The idea here is to see visually by how much the CRI mechanism condenses the same subset of tropospheric chemistry (namely inorganic+C${1}$-C${5}$ alkane oxidation) in comparison with the MCM.
The input data was downloaded from the corresponding websites and pre-processed into json files following the approach described here.
End of explanation
d = nx.degree_centrality(scheme) # calculate degree of centrality
{k: d[k] for k in sorted(d.keys())[:10]} # show the first 10 entries of a dictionary with degrees of centrality
max([val for key, val in d.items()]) # max degree of centrality in the whole network
# Sort a dictionary with degrees of centrality by value and print the first 11 key, value pairs
sorted_degree_centrality = [(k, d[k]) for k in sorted(d, key=d.get, reverse=True)]
for n, i in enumerate(sorted_degree_centrality[0:11]):
    print(n+1, i)
Explanation: Task 1: Find the species that has the maximum number of connections.
End of explanation
list(scheme.predecessors('O3'))
Explanation: So, unsurprisingly, the hydroxyl radical (OH) has the maximum number of interactions in the network (actually in both, the MCM and the CRI). Ozone is on 11th place in the MCM, but only 40th in the CRI.
Task 2. Find species reactions between which directly lead to ozone formation.
End of explanation
list(scheme.successors('O3'))
Explanation: Task 3. Find species that are produced directly in ozone destroying reactions.
End of explanation
HTML(html)
Explanation: We were able to find the last two answers only because at the beginning we used MultiDiGraph(), which is a directed graph.
If we load the CRI mechanism and rerun the same cells, the layout of the CRI network is:
End of explanation
<END_TASK>