def netStats(opt_fname,dummy_fname):
    print 'Evaluating statistics from %s,%s'%(opt_fname,dummy_fname)

    opt_tab=[a.strip().split('\tpp\t') for a in open(opt_fname,'r').readlines()]
    dummy_tab=[[a.strip().split('\t')[0],a.strip().split('\t')[2]] for a in open(dummy_fname,'r').readlines()]


    stats={}
    ##get optimal number of edges
    opt_g=networkx.from_edgelist([(a[0],a[1]) for a in opt_tab])
    dum_g=networkx.from_edgelist([(a[0],a[1]) for a in dummy_tab])

    stats['numEdges'] = len(opt_tab)
    ##get number of trees (dummy - opt)
    stats['numTrees'] = networkx.number_connected_components(opt_g) #len(dummy_tab) - len(opt_tab)

    conn_comps=networkx.connected_components(opt_g)

    ##number of nodes
    nodes=opt_g.nodes()
    stats['numNodes']=len(nodes)

    ##get tree sizes
    stats['treeSizes']=','.join([str(len(c)) for c in sorted(conn_comps, key=len, reverse=True)])

    ##find ubiquitin?
    if 'UBC' in nodes:
    	stats['hasUBC']='True'
    else:
	stats['hasUBC']='False'

    return stats

def split_nx(mesh, check_watertight=True, only_count=False):
    '''
    Given a mesh, will split it up into a list of meshes based on face connectivity
    If check_watertight is true, it will only return meshes where each face has
    exactly 3 adjacent faces, which is a simple metric for being watertight.
    '''
    def mesh_from_components(connected_faces):
        if check_watertight:
            subgraph   = nx.subgraph(face_adjacency, connected_faces)
            watertight = np.equal(list(subgraph.degree().values()), 3).all()
            if not watertight: return
        faces  = mesh.faces[[connected_faces]]
        unique = np.unique(faces.reshape(-1))
        replacement = dict()
        replacement.update(np.column_stack((unique, np.arange(len(unique)))))
        faces = replace_references(faces, replacement).reshape((-1,3))
        new_meshes.append(mesh.__class__(vertices     = mesh.vertices[[unique]],
                                         faces        = faces,
                                         face_normals = mesh.face_normals[[connected_faces]]))
    face_adjacency = nx.from_edgelist(mesh.face_adjacency())
    new_meshes     = deque()
    components     = list(nx.connected_components(face_adjacency))
    if only_count: return len(components)

    for component in components: mesh_from_components(component)
    log.info('split mesh into %i components.',
             len(new_meshes))
    return list(new_meshes)

 def _process(self, element, key=None):
     if self.p.layout and isinstance(self.p.layout, FunctionType):
         import networkx as nx
         edges = element.array([0, 1])
         graph = nx.from_edgelist(edges)
         if 'weight' in self.p.kwargs:
             weight = self.p.kwargs['weight']
             for (s, t), w in zip(edges, element[weight]):
                 graph.edges[s, t][weight] = w
         positions = self.p.layout(graph, **self.p.kwargs)
         nodes = [tuple(pos)+(idx,) for idx, pos in sorted(positions.items())]
     else:
         source = element.dimension_values(0, expanded=False)
         target = element.dimension_values(1, expanded=False)
         nodes = np.unique(np.concatenate([source, target]))
         if self.p.layout:
             import pandas as pd
             df = pd.DataFrame({'index': nodes})
             nodes = self.p.layout(df, element.dframe(), **self.p.kwargs)
             nodes = nodes[['x', 'y', 'index']]
         else:
             nodes = circular_layout(nodes)
     nodes = element.node_type(nodes)
     if element._nodes:
         for d in element.nodes.vdims:
             vals = element.nodes.dimension_values(d)
             nodes = nodes.add_dimension(d, len(nodes.vdims), vals, vdim=True)
     if self.p.only_nodes:
         return nodes
     return element.clone((element.data, nodes))

    def _addNonbondedForceToSystem(self, sys, verbose):
        '''Create the nonbonded force
        '''
        nb = mm.NonbondedForce()
        sys.addForce(nb)

        q = '''SELECT charge, sigma, epsilon
        FROM particle INNER JOIN nonbonded_param
        ON particle.nbtype=nonbonded_param.id'''
        for charge, sigma, epsilon in self._conn.execute(q):
            nb.addParticle(charge, sigma*angstrom, epsilon*kilocalorie_per_mole)

        if verbose:
            # Bond graph (for debugging)
            g = nx.from_edgelist(self._conn.execute('SELECT p0, p1 FROM stretch_harm_term').fetchall())
            nbnames = {1: '1-2', 2:'1-3', 3:'1-4'}

        q = '''SELECT p0, p1, aij, bij, qij
        FROM pair_12_6_es_term INNER JOIN pair_12_6_es_param
        ON pair_12_6_es_term.param=pair_12_6_es_param.id;'''
        for p0, p1, a_ij, b_ij, q_ij in self._conn.execute(q):
            if verbose:
                l = nx.algorithms.shortest_path_length(g, p0, p1)
                print 'Scaling interaction for a %d-%d (%s) interaction' % (p0, p1, nbnames[l])
            a_ij = (a_ij*kilocalorie_per_mole*(angstrom**12)).in_units_of(kilojoule_per_mole*(nanometer**12))
            b_ij = (b_ij*kilocalorie_per_mole*(angstrom**6)).in_units_of(kilojoule_per_mole*(nanometer**6))
            q_ij = q_ij*elementary_charge**2

            if (b_ij._value == 0.0) or (a_ij._value == 0.0):
                new_epsilon = 0
                new_sigma = 1
            else:
                new_epsilon =  b_ij**2/(4*a_ij)
                new_sigma = (a_ij / b_ij)**(1.0/6.0)
            nb.addException(p0, p1, q_ij, new_sigma, new_epsilon)

        n_total = self._conn.execute('''SELECT COUNT(*) FROM pair_12_6_es_term''').fetchone()
        n_in_exclusions= self._conn.execute('''SELECT COUNT(*)
        FROM exclusion INNER JOIN pair_12_6_es_term
        ON exclusion.p0==pair_12_6_es_term.p0 AND exclusion.p1==pair_12_6_es_term.p1''').fetchone()
        if not n_total == n_in_exclusions:
            raise NotImplementedError('All pair_12_6_es_terms must have a corresponding exclusion')

        # Desmond puts scaled 1-4 interactions in the pair_12_6_es
        # table, and then adds a corresponding exception here. We are
        # using the exception part of NonbondedForce, so we're just
        # adding the 1-4 interaction as an exception when its
        # registered, and then NOT registering it as an exception here.
        q = '''SELECT E.p0, E.p1
        FROM exclusion E LEFT OUTER JOIN pair_12_6_es_term P ON
        E.p0 = P.p0 and E.p1 = P.p1
        WHERE P.p0 is NULL'''
        # http://stackoverflow.com/questions/5464131/finding-pairs-that-do-not-exist-in-a-different-table
        for p0, p1 in self._conn.execute(q):
            if verbose:
                l = nx.algorithms.shortest_path_length(g, p0, p1)
                print 'Creating exception for a %d-%d (%s) interaction' % (p0, p1, nbnames[l])
            nb.addException(p0, p1, 0.0, 1.0, 0.0)

        return nb

def fix_face_winding(mesh):
    '''
    Traverse and change mesh faces in-place to make sure winding is coherent, 
    or that edges on adjacent faces are in opposite directions
    '''
    # we create the face adjacency graph: 
    # every node in g is an index of mesh.faces
    # every edge in g represents two faces which are connected
    graph_all = nx.from_edgelist(mesh.face_adjacency)
    flipped   = 0
    # we are going to traverse the graph using BFS, so we have to start
    # a traversal for every connected component
    for graph in nx.connected_component_subgraphs(graph_all):
        start = graph.nodes()[0]
        # we traverse every pair of faces in the graph
        # we modify mesh.faces and mesh.face_normals in place 
        for face_pair in nx.bfs_edges(graph, start):
            # for each pair of faces, we convert them into edges,
            # find the edge that both faces share, and then see if the edges
            # are reversed in order as you would expect in a well constructed mesh
            pair    = mesh.faces[[face_pair]]            
            edges   = faces_to_edges(pair)
            overlap = group_rows(np.sort(edges,axis=1), require_count=2)
            if len(overlap) == 0:
                # only happens on non-watertight meshes
                continue
            edge_pair = edges[[overlap[0]]]
            if edge_pair[0][0] == edge_pair[1][0]:
                # if the edges aren't reversed, invert the order of one of the faces
                flipped += 1
                mesh.faces[face_pair[1]] = mesh.faces[face_pair[1]][::-1]
    log.info('Flipped %d/%d edges', flipped, len(mesh.faces)*3)

def angle_connectivity(ibonds):
    """Given the bonds, get the indices of the atoms defining all the bond
    angles

    A 'bond angle' is defined as any set of 3 atoms, `i`, `j`, `k` such that
    atom `i` is bonded to `j` and `j` is bonded to `k`

    Parameters
    ----------
    ibonds : np.ndarray, shape=[n_bonds, 2], dtype=int
        Each row in `ibonds` is a pair of indicies `i`, `j`, indicating that
        atoms `i` and `j` are bonded

    Returns
    -------
    iangles : np.ndarray, shape[n_angles, 3], dtype=int
        n_angles x 3 array of indices, where each row is the index of three
        atoms m,n,o such that n is bonded to both m and o.
    """

    graph = nx.from_edgelist(ibonds)
    iangles = []

    for i in graph.nodes():
        for (m, n) in combinations(graph.neighbors(i), 2):
            # so now the there is a bond angle m-i-n
            iangles.append((m, i, n))

    return np.array(iangles)

def get_dihedral_connectivity(ibonds):
    """Given the bonds, get the indices of the atoms defining all the dihedral
    angles
    
    Parameters
    ----------
    ibonds : np.ndarray, shape=[n_bonds, 2], dtype=int
        n_bonds x 2 array of indices, where each row is the index of two
        atom who participate in a bond.
    
    Returns
    -------
    idihedrals : np.ndarray, shape[n_dihedrals, 4], dtype=int
        All sets of 4 atoms A,B,C,D such that A is bonded to B, B is bonded
        to C, and C is bonded to D
    """
    graph = nx.from_edgelist(ibonds)
    n_atoms = graph.number_of_nodes()
    idihedrals = []
    
    # TODO: CHECK FOR DIHEDRAL ANGLES THAT ARE 180 and recover
    # conf : msmbuilder.Trajectory
    #    An msmbuilder trajectory, only the first frame will be used. This
    #    is used purely to make the check for angle(ABC) != 180.

    for a in xrange(n_atoms):
        for b in graph.neighbors(a):
            for c in ifilter(lambda c: c not in [a, b], graph.neighbors(b)):
                for d in ifilter(lambda d: d not in [a, b, c], graph.neighbors(c)):
                    idihedrals.append((a, b, c, d))

    return np.array(idihedrals)

def get_angle_connectivity(ibonds):
    """Given the bonds, get the indices of the atoms defining all the bond
    angles
    
    Parameters
    ----------
    ibonds : np.ndarray, shape=[n_bonds, 2], dtype=int
        n_bonds x 2 array of indices, where each row is the index of two
        atom who participate in a bond.
    
    Returns
    -------
    iangles : np.ndarray, shape[n_angles, 3], dtype=int
        n_angles x 3 array of indices, where each row is the index of three
        atoms m,n,o such that n is bonded to both m and o.
    """

    graph = nx.from_edgelist(ibonds)
    n_atoms = graph.number_of_nodes()
    iangles = []

    for i in xrange(n_atoms):
        for (m, n) in combinations(graph.neighbors(i), 2):
            # so now the there is a bond angle m-i-n
            iangles.append((m, i, n))

    return np.array(iangles)

def network_analysis(gene_list,network_file,outdir):
    outfn = "%s/output" % outdir
    f = open(outfn,'w')
    f.write("gene\tdegrees\tbtw_centrality\n")
    network = networkx.read_adjlist(network_file)
    print "Number of edges in input graph: %s" % network.number_of_edges()
    print "Number of nodes in input graph: %s" % network.number_of_nodes()
    subnetwork = network.subgraph(gene_list)
    print "Number of edges in subgraph: %s" % subnetwork.number_of_edges()
    print "Number of nodes in subgraph: %s" % subnetwork.number_of_nodes()
    bwt_central = networkx.betweenness_centrality(subnetwork)
    degrees = subnetwork.degree(gene_list)
    for gene in gene_list:
        # Number of degrees
        if gene in degrees:
            num_degrees = degrees[gene]
        else:
            num_degress = "NA"
        # Betweenness centrality
        if gene in bwt_central:
            btw_gene = bwt_central[gene]
        else:
            btw_gene = "NA"
        # File with neighbor nodes
        if subnetwork.has_node(gene):
            neighbors = list(networkx.all_neighbors(subnetwork,gene))
            edges = [(unicode(gene),neighbor) for neighbor in neighbors]
            neighbor_networks = networkx.from_edgelist(edges)
            write_networks(neighbor_networks,gene,outdir)
        f.write("%s\t%s\t%s\n" % (gene,num_degrees,btw_gene))
    f.close()

def broken_faces(mesh, color=None):
    """
    Return the index of faces in the mesh which break the
    watertight status of the mesh.

    Parameters
    --------------
    mesh: Trimesh object
    color: (4,) uint8, will set broken faces to this color
           None,       will not alter mesh colors

    Returns
    ---------------
    broken: (n, ) int, indexes of mesh.faces
    """
    adjacency = nx.from_edgelist(mesh.face_adjacency)
    broken = [k for k, v in dict(adjacency.degree()).items()
              if v != 3]
    broken = np.array(broken)
    if color is not None:
        # if someone passed a broken color
        color = np.array(color)
        if not (color.shape == (4,) or color.shape == (3,)):
            color = [255, 0, 0, 255]
        mesh.visual.face_colors[broken] = color
    return broken

def fix_winding(mesh):
    """
    Traverse and change mesh faces in-place to make sure winding
    is correct, with edges on adjacent faces in
    opposite directions.

    Parameters
    -------------
    mesh: Trimesh object

    Alters
    -------------
    mesh.face: will reverse columns of certain faces
    """
    # anything we would fix is already done
    if mesh.is_winding_consistent:
        return

    graph_all = nx.from_edgelist(mesh.face_adjacency)
    flipped = 0

    faces = mesh.faces.view(np.ndarray).copy()

    # we are going to traverse the graph using BFS
    # start a traversal for every connected component
    for components in nx.connected_components(graph_all):
        # get a subgraph for this component
        g = graph_all.subgraph(components)
        # get the first node in the graph in a way that works on nx's
        # new API and their old API
        start = next(iter(g.nodes()))

        # we traverse every pair of faces in the graph
        # we modify mesh.faces and mesh.face_normals in place
        for face_pair in nx.bfs_edges(g, start):
            # for each pair of faces, we convert them into edges,
            # find the edge that both faces share and then see if edges
            # are reversed in order as you would expect
            # (2, ) int
            face_pair = np.ravel(face_pair)
            # (2, 3) int
            pair = faces[face_pair]
            # (6, 2) int
            edges = faces_to_edges(pair)
            overlap = group_rows(np.sort(edges, axis=1),
                                 require_count=2)
            if len(overlap) == 0:
                # only happens on non-watertight meshes
                continue
            edge_pair = edges[overlap[0]]
            if edge_pair[0][0] == edge_pair[1][0]:
                # if the edges aren't reversed, invert the order of one face
                flipped += 1
                faces[face_pair[1]] = faces[face_pair[1]][::-1]

    if flipped > 0:
        mesh.faces = faces

    log.debug('flipped %d/%d edges', flipped, len(mesh.faces) * 3)

def load_from_edge_list(filename):
    edgelist = []
    with open(filename, "r") as fh:
        for line in fh.readlines():
            source, target = line.split(",")
            edgelist.append((int(source), int(target)))

    return nx.from_edgelist(edgelist)

 def __init__(self, name, group_type, size, hours, school, data):
     self.name = name
     self.group_type = group_type
     self.size = size
     self.hours = hours
     self.school = school
     self.data = data
     self.G = nx.from_edgelist(self.data)

def drawGraph(inputFile, outputFile):
	listEdge = pd.read_csv(inputFile)
	arrayList = listEdge.ix[:,[2,3]].as_matrix()
	arrayList
	G = nx.from_edgelist(arrayList)
	limits=plt.axis('off')
	
	nx.draw_networkx(G,with_labels=True, font_size=5, node_size=500)
	plt.savefig(outputFile)

def fix_normals(mesh):
    '''
    Find and fix problems with mesh.face_normals and mesh.faces winding direction.
    
    For face normals ensure that vectors are consistently pointed outwards,
    and that mesh.faces is wound in the correct direction for all connected components.
    '''
    mesh.generate_face_normals()
    # we create the face adjacency graph: 
    # every node in g is an index of mesh.faces
    # every edge in g represents two faces which are connected
    graph = nx.from_edgelist(mesh.face_adjacency())
    
    # we are going to traverse the graph using BFS, so we have to start
    # a traversal for every connected component
    for connected in nx.connected_components(graph):
        # we traverse every pair of faces in the graph
        # we modify mesh.faces and mesh.face_normals in place 
        for face_pair in nx.bfs_edges(graph, connected[0]):
            # for each pair of faces, we convert them into edges,
            # find the edge that both faces share, and then see if the edges
            # are reversed in order as you would expect in a well constructed mesh
            pair      = mesh.faces[[face_pair]]
            edges     = faces_to_edges(pair, sort=False)
            overlap   = group_rows(np.sort(edges,axis=1), require_count=2)
            edge_pair = edges[[overlap[0]]]
            reversed  = edge_pair[0][0] != edge_pair[1][0]
            if reversed: continue
            # if the edges aren't reversed, invert the order of one of the faces
            # and negate its normal vector
            mesh.faces[face_pair[1]] = mesh.faces[face_pair[1]][::-1]
            mesh.face_normals[face_pair[1]] *= (reversed*2) - 1
            
        # the normals of every connected face now all pointed in 
        # the same direction, but there is no guarantee that they aren't all
        # pointed in the wrong direction
        faces           = mesh.faces[[connected]]
        faces_x         = np.min(mesh.vertices[:,0][[faces]], axis=1)
        left_order      = np.argsort(faces_x)
        left_values     = faces_x[left_order]
        left_candidates = np.abs(left_values - left_values[0]) < TOL_ZERO
        backwards       = None
        
        # note that we have to find a face which ISN'T perpendicular to the x axis 
        # thus we go through all the candidate faces that are at the extreme left
        # until we find one that has a nonzero dot product with the x axis
        for leftmost in left_order[left_candidates]:                
            face_dot = np.dot([-1.0,0,0], mesh.face_normals[leftmost]) 
            if abs(face_dot) > TOL_ZERO: 
                backwards = face_dot < 0.0
                break
        if backwards: mesh.face_normals[[connected]] *= -1.0
        
        winding_tri  = connected[0]
        winding_test = np.diff(mesh.vertices[[mesh.faces[winding_tri]]], axis=0)
        winding_dir  = np.dot(unitize(np.cross(*winding_test)), mesh.face_normals[winding_tri])
        if winding_dir < 0: mesh.faces[[connected]] = np.fliplr(mesh.faces[[connected]])

def construct_constrained_graph(adj, r, n):
    """
    given an adjacency matrix adj in the form of a condensed distance matrix
    (of the kind returned by pdist) for n observations, returns the similarity
    graph for all distances less than or equal to r.
    """
    ij = row_col_from_condensed_index(n, np.where(adj<=r)[0])
    g = nx.from_edgelist(zip(*ij))
    g.add_nodes_from(range(n))
    return g

 def facets_nx():
     graph_parallel = nx.from_edgelist(face_idx[parallel])
     facets_idx = np.array([list(i) for i in nx.connected_components(graph_parallel)])
     #### commented by weiwei
     # should also return the single triangles
     facets_idx_extra = copy.deepcopy(facets_idx.tolist())
     for item in range(mesh.faces.shape[0]):
         if item not in [i for subitem in facets_idx.tolist() for i in subitem]:
             facets_idx_extra.append([item])
     return np.array(facets_idx_extra)

def cube_neighbor_graph(file_name):
  # Builds and returns the 3x3x3 cube Game of Life neighbor graph given a file
  # of line adjacencies.
  with open(file_name, 'rb') as f:
    line, line_prev, edge_list = [], [], []
    for line in (map(int, line.strip().split()) for line in f):
      if line and line_prev:
        for edge in line_edges(line_prev, line): edge_list.append(edge)
      line_prev = line
  return nx.from_edgelist(edge_list, create_using=nx.Graph())

 def social_following(self):
     # FIXME: add proper site context to ploneintranet.network graph
     graph = queryUtility(INetworkGraph)
     result = []
     # FIXME: add proper API accessor to ploneintranet.network graph
     for user in graph._following["user"].keys():
         for following in graph.get_following("user", user):
             result.append((user, following))
     return nx.from_edgelist(result,
                             create_using=nx.DiGraph())

    def calculate(self):
        """Lazy initialization."""

        # TODO: ensure that this is fully unrestricted / runs as admin

        catalog = api.portal.get_tool('portal_catalog')
        content_tree = []
        content_authors = []
        content_tags = []
        for brain in catalog():
            context = brain.getObject()
            context_path = 'path:%s' % '/'.join(context.getPhysicalPath())

            for child_id in context.objectIds():
                child_path = '%s/%s' % (context_path, child_id)
                # containment is a bidirectional relationship
                content_tree.append((context_path, child_path))
                content_tree.append((child_path, context_path))

            # TODO: add reference links
            # TODO: add text links

            for author in context.creators:
                # authorship is a bidirectional relationship
                content_authors.append((context_path, 'user:%s' % author))
                content_authors.append(('user:%s' % author, context_path))

            # TODO: add sharing

            for tag in context.Subject():
                # tagging is bidirectional between context and tag
                content_tags.append((context_path, 'tag:%s' % tag))
                content_tags.append(('tag:%s' % tag, context_path))

        self._cache['content_tree'] = nx.from_edgelist(
            content_tree,
            create_using=nx.DiGraph())
        self._cache['content_authors'] = nx.from_edgelist(
            content_authors,
            create_using=nx.DiGraph())
        self._cache['content_tags'] = nx.from_edgelist(
            content_tags,
            create_using=nx.DiGraph())