def setUp(self):
     self.G_basic = tg.node_generate(
         [[10, 11], [20], [30, 31], [40, 41], [50]],
         itertools.count(step=100))
     self.G_basic.add_path([10, 20, 30, 40, 50])
     self.G_basic.add_path([11, 20, 31])
     self.G_basic.add_path([30, 41, 50])
     nx.freeze(self.G_basic)

def MakeAdjacencyGraph(n, type, freq=0.5, wexc=0.1):
  '''Temporary implementation that only generates complete graph.'''
  G = nx.complete_graph(n, nx.DiGraph())
  # Should profile this.
  for edge in G.edges_iter():
    nx.set_edge_attributes(G, 'weight', {edge: wexc})
  nx.freeze(G) # To prevent unintended modification of G's structure
  return G

def build_meetups_graph(groups, topics, members, events):
    graph = nx.Graph()
    for group in groups:
        key = group['id']
        graph.add_node(key, name=group['name'], type=Type.Groups, members=group['members'])
        link_members(graph, group,members)
        link_topics(graph, group, topics)
    link_members_to_topics(graph, members, topics)
    link_members_to_events(graph, members, events)
    nx.freeze(graph)

    logging.info('graph built')
    return graph

 def setUp(self):
     G=nx.Graph(name="test")
     e=[('a','b'),('b','c'),('c','d'),('d','e'),('e','f'),('a','f')]
     G.add_edges_from(e,width=10)
     G.add_node('g',color='green')
     G.graph['number']=1
     DG=nx.DiGraph(G)
     MG=nx.MultiGraph(G)
     MG.add_edge('a', 'a')
     MDG=nx.MultiDiGraph(G)
     MDG.add_edge('a', 'a')
     fG = G.copy()
     fDG = DG.copy()
     fMG = MG.copy()
     fMDG = MDG.copy()
     nx.freeze(fG)
     nx.freeze(fDG)
     nx.freeze(fMG)
     nx.freeze(fMDG)
     self.G=G
     self.DG=DG
     self.MG=MG
     self.MDG=MDG
     self.fG=fG
     self.fDG=fDG
     self.fMG=fMG
     self.fMDG=fMDG

def main():
  seed(0) #set seed
  #get graph info
  G = nx.read_gpickle("input/graphMTC_CentroidsLength6.gpickle") #noCentroidsLength15.gpickle") #does not have centroidal links. There is also the choice of a proper multidigraph: nx.read_gpickle("input/graphMTC_CentroidsLength5.gpickle")
  G = nx.freeze(G) #prevents edges or nodes to be added or deleted
  #get od info. This is in format of a dict keyed by od, like demand[sd1][sd2] = 200000.
  demand = bd.build_demand('input/BATS2000_34SuperD_TripTableData.csv', 'input/superdistricts_centroids.csv')
  #get earthquake info
  q = QuakeMaps('input/20130210_mtc_total_lnsas3.pkl', 'input/20130210_mtc_magnitudes3.pkl', 'input/20130210_mtc_faults3.pkl', 'input/20130210_mtc_weights3.pkl', 'input/20130210_mtc_scenarios3.pkl') #input/20130107_mtc_total_lnsas1.pkl', 'input/20130107_mtc_magnitudes1.pkl','input/20130107_mtc_faults1.pkl', 'input/20130107_mtc_weights1.pkl', 'input/20130107_mtc_scenarios1.pkl') #'input/20130210_mtc_total_lnsas3.pkl', 'input/20130210_mtc_magnitudes3.pkl', 'input/20130210_mtc_faults3.pkl', 'input/20130210_mtc_weights3.pkl', 'input/20130210_mtc_scenarios3.pkl') #('input/20130107_mtc_total_lnsas1.pkl', 'input/20130107_mtc_magnitudes1.pkl',  #totalfilename=None, magfilename=None, faultfilename=None, weightsfilename=None, scenariofilename=None):
  print 'weights: ', q.weights
  q.num_sites = len(q.lnsas[0])
  #determine which scenarios you want to run
  good_indices = pick_scenarios(q.lnsas, q.weights)
  
  travel_index_times = []
  index = 0
  #loop over scenarios
  print 'size of lnsas: ', len(q.lnsas)
  for scenario in q.lnsas: #each 'scenario' has 1557 values of lnsa, i.e. one per site
    if index in good_indices:
      print 'index: ', index
      (bridges, flow, path, path2) = run_simple_iteration(G, scenario, demand, False)
      travel_index_times.append((index, bridges, flow, path, path2))
#      print 'new travel times: ', travel_index_times
      if index%1000 ==0:
        util.write_2dlist(time.strftime("%Y%m%d")+'_bridges_flow_paths4.txt',travel_index_times)
    index += 1 #IMPORTANT
  util.write_2dlist(time.strftime("%Y%m%d")+'_bridges_flow_paths4.txt',travel_index_times)
  print 'the number of scenarios I considered doing: ', index
  print 'the number of scenarios I actually did: ', len(travel_index_times)

def generic_graph_view(G, create_using=None):
    if create_using is None:
        newG = G.__class__()
    else:
        newG = nx.empty_graph(0, create_using)
    if G.is_multigraph() != newG.is_multigraph():
        raise NetworkXError("Multigraph for G must agree with create_using")
    newG = nx.freeze(newG)

    # create view by assigning attributes from G
    newG._graph = G
    newG.graph = G.graph

    newG._node = G._node
    if newG.is_directed():
        if G.is_directed():
            newG._succ = G._succ
            newG._pred = G._pred
            newG._adj = G._succ
        else:
            newG._succ = G._adj
            newG._pred = G._adj
            newG._adj = G._adj
    elif G.is_directed():
        if G.is_multigraph():
            newG._adj = UnionMultiAdjacency(G._succ, G._pred)
        else:
            newG._adj = UnionAdjacency(G._succ, G._pred)
    else:
        newG._adj = G._adj
    return newG

def main():
  seed(0) #set seed
  #get graph info
  G = nx.read_gpickle("input/graphMTC_CentroidsLength5.gpickle") #noCentroidsLength15.gpickle") #does not have centroidal links
  print '|V| = ', len(G.nodes())
  print '|E| = ', len(G.edges())
  G = nx.freeze(G) #prevents edges or nodes to be added or deleted
  #get od info. This is in format of a dict keyed by od, like demand[sd1][sd2] = 200000.
  demand = bd.build_demand('input/BATS2000_34SuperD_TripTableData.csv', 'input/superdistricts_centroids.csv') #bd.build_demand('input/BATS2000_34SuperD_TripTableData.csv', 'input/superdistricts_centroids.csv')
  #get earthquake info
  q = QuakeMaps('input/20130210_mtc_total_lnsas3.pkl', 'input/20130210_mtc_magnitudes3.pkl', 'input/20130210_mtc_faults3.pkl', 'input/20130210_mtc_weights3.pkl', 'input/20130210_mtc_scenarios3.pkl') #(input/20130107_mtc_total_lnsas1.pkl', 'input/20130107_mtc_magnitudes1.pkl', 'input/20130107_mtc_faults1.pkl', 'input/20130107_mtc_weights1.pkl', 'input/20130107_mtc_scenarios1.pkl') #totalfilename=None, magfilename=None, faultfilename=None, weightsfilename=None, scenariofilename=None): 'input/20130210_mtc_total_lnsas3.pkl', 'input/20130210_mtc_magnitudes3.pkl', 'input/20130210_mtc_faults3.pkl', 'input/20130210_mtc_weights3.pkl', 'input/20130210_mtc_scenarios3.pkl') #(


  q.num_sites = len(q.lnsas[0])
  #determine which scenarios you want to run
  good_indices = pick_scenarios(q.lnsas, q.weights)
  
  travel_index_times = []
  index = 0
  #loop over scenarios
  for scenario in q.lnsas: #each 'scenario' has 1557 values of lnsa, i.e. one per site
    if index in good_indices:
      print 'index: ', index
      (travel_time, vmt) = run_iteration(G, scenario, demand)
      travel_index_times.append((index, travel_time, vmt))
#      print 'new travel times: ', travel_index_times
      if index%100 ==0:
        util.write_2dlist(time.strftime("%Y%m%d")+'_travel_time.txt',travel_index_times)
    index += 1 #IMPORTANT
  util.write_2dlist(time.strftime("%Y%m%d")+'_travel_time.txt',travel_index_times)

def subgraph_view(G, filter_node=no_filter, filter_edge=no_filter):
    newG = nx.freeze(G.__class__())
    newG._NODE_OK = filter_node
    newG._EDGE_OK = filter_edge

    # create view by assigning attributes from G
    newG._graph = G
    newG.graph = G.graph

    newG._node = FilterAtlas(G._node, filter_node)
    if G.is_multigraph():
        Adj = FilterMultiAdjacency

        def reverse_edge(u, v, k): return filter_edge(v, u, k)
    else:
        Adj = FilterAdjacency

        def reverse_edge(u, v): return filter_edge(v, u)
    if G.is_directed():
        newG._succ = Adj(G._succ, filter_node, filter_edge)
        newG._pred = Adj(G._pred, filter_node, reverse_edge)
        newG._adj = newG._succ
    else:
        newG._adj = Adj(G._adj, filter_node, filter_edge)
    return newG

 def _swap(self, replacement_graph):
     """Validates the replacement graph and then swaps the underlying graph
     with a frozen version of the replacement graph (this maintains the
     invariant that the underlying graph is immutable).
     """
     self._validate(replacement_graph)
     self._graph = nx.freeze(replacement_graph)

    def __init__(self, graph:str or networkx.Graph):
        if isinstance(graph, networkx.Graph):
            nxgraph = networkx.Graph(graph)
        elif isinstance(graph, str):
            nxgraph = phasme.build_graph.graph_from_file(graph)
        else:
            raise ValueError("Unexpected {}".format(graph))
        # internal graph representation
        self.__edges = map(frozenset, nxgraph.edges)
        if const.TEST_INTEGRITY:
            self.__edges = tuple(self.__edges)
            for args in self.__edges:
                if len(args) > 2:
                    print('WARNING: Weird edge: {}. It will be filtered.'.format(args))
                else:
                    assert len(args) == 2, args
        self.__edges = set(edge for edge in self.__edges if len(edge) == 2)

        # data
        self.__nodes = set(nxgraph.nodes)
        self.__uid = str(min(self.__nodes, key=str))
        self.__nb_node = len(self.__nodes)
        self.__nb_cc = networkx.number_connected_components(nxgraph)
        self._nxgraph = nxgraph if constants.KEEP_NX_GRAPH else networkx.freeze(nxgraph)
        self.__initial_number_of_edge = len(self.__edges)
        self.__active_recipe = None

        self.__hierarchy = set()  # inclusions between powernodes
        self.__powernodes = defaultdict(set)  # (step, set) -> {node in powernode}
        self.__poweredges = defaultdict(set)  # (step, set) -> (step, set)

def flatten(item, freeze=True):
    """Flattens a item (a task or flow) into a single execution graph."""
    graph = _post_flatten(_flatten(item, set()))
    if freeze:
        # Frozen graph can't be modified...
        return nx.freeze(graph)
    return graph

def get_graph():
  import networkx
  '''loads full mtc highway graph with dummy links and then adds a few fake centroidal nodes for max flow and traffic assignment'''
  G = networkx.read_gpickle("input/graphMTC_CentroidsLength3int.gpickle")
  G = add_superdistrict_centroids(G)
   # Directed! only one edge between nodes
  G = nx.freeze(G) #prevents edges or nodes to be added or deleted
  return G

 def graph(self):
     if self._subgraph is not None:
         return self._subgraph
     if self._target is None:
         return self._graph
     nodes = [self._target]
     nodes.extend(dst for _src, dst in
                  traversal.dfs_edges(self._graph.reverse(), self._target))
     self._subgraph = nx.freeze(self._graph.subgraph(nodes))
     return self._subgraph

 def flatten(self):
     """Flattens a item (a task or flow) into a single execution graph."""
     if self._graph is not None:
         return self._graph
     self._pre_flatten()
     graph = self._flatten(self._root)
     self._post_flatten(graph)
     if self._freeze:
         self._graph = nx.freeze(graph)
     else:
         self._graph = graph
     return self._graph

 def test_freeze(self):
     G=networkx.freeze(self.G)
     assert_equal(G.frozen,True)
     assert_raises(networkx.NetworkXError, G.add_node, 1)
     assert_raises(networkx.NetworkXError, G.add_nodes_from, [1])
     assert_raises(networkx.NetworkXError, G.remove_node, 1)
     assert_raises(networkx.NetworkXError, G.remove_nodes_from, [1])
     assert_raises(networkx.NetworkXError, G.add_edge, 1,2)
     assert_raises(networkx.NetworkXError, G.add_edges_from, [(1,2)])
     assert_raises(networkx.NetworkXError, G.remove_edge, 1,2)
     assert_raises(networkx.NetworkXError, G.remove_edges_from, [(1,2)])
     assert_raises(networkx.NetworkXError, G.clear)

def main():
  '''can change the number of epsilons below'''
  seed(0) #set seed
  simple = False  #simple is just %bridges out, which is computationally efficient
  number_of_highway_bridges = 1743
  numeps = 3 #the number of epsilons
  tol = 0.00001 #the minimum annual rate that you care about in the original event set (the weight now is the original annual rate / number of epsilons per event)
  demand = bd.build_demand('input/BATS2000_34SuperD_TripTableData.csv', 'input/superdistricts_centroids.csv') #we just take a percentage in ita.py, namely  #to get morning flows, take 5.3% of daily driver values. 11.5/(4.5*6+11.5*10+14*4+4.5*4) from Figure S10 of http://www.nature.com/srep/2012/121220/srep01001/extref/srep01001-s1.pdf
  #figure out ground motions
  lnsas, weights = ground_motions(numeps, tol, '/Users/mahalia/Documents/matlab/Research/Herbst2011/output_data/SF2_mtc_total_3909scenarios_1743bridgesPlusBART_3eps.txt')
  bart_dict = transit_to_damage.make_bart_dict()
  muni_dict = transit_to_damage.make_muni_dict()
  set_main_path('/Users/mahaliamiller/Desktop/trn/transit_lines/', None) #TODO: need to change THREE file paths (these plus bart)

  print 'the number of ground motion events we are considering: ', len(lnsas)
  index = 0
  bridge_array = []
  travel_index_times = []

  # G = nx.read_gpickle("input/graphMTC_noCentroidsLength15.gpickle")
  G = nx.read_gpickle("input/graphMTC_CentroidsLength6.gpickle")
   # Directed! only one edge between nodes
  G = nx.freeze(G) #prevents edges or nodes to be added or deleted
  print 'am I a multi graph? ', G.is_multigraph()
  no_damage_travel_time, no_damage_vmt = compute_tt_vmt(G, demand)
  if not os.path.isdir(time.strftime("%Y%m%d")+'_filesForCube/'):
    os.mkdir(time.strftime("%Y%m%d")+'_filesForCube/')
  if not os.path.isdir(time.strftime("%Y%m%d")+'_filesForCube/transit/'):
    os.mkdir(time.strftime("%Y%m%d")+'_filesForCube/transit/')
  if not os.path.isdir(time.strftime("%Y%m%d")+'_filesForCube/modCapacities/'):
    os.mkdir(time.strftime("%Y%m%d")+'_filesForCube/modCapacities/')

  for scenario in lnsas:
    print index
    #figure out bridge damage for each scenario
    damaged_bridges, num_bridges_out = damage_bridges(scenario) #e.g., [1, 89, 598] #num_bridges_out is highway bridges only
    bridge_array.append(damaged_bridges)

    #figure out network damage and output Cube files to this effect
    G = damage_network(damaged_bridges, G, time.strftime("%Y%m%d")+'_filesForCube/', index)

    #figure out impact (performance metrics)
    flow, shortest_paths, travel_time, vmt = measure_performance(G, damaged_bridges, demand, no_damage_travel_time, no_damage_vmt)
    travel_index_times.append((index, num_bridges_out, flow, shortest_paths, travel_time, vmt, num_bridges_out/float(number_of_highway_bridges)))
    G = util.clean_up_graph(G)
    index +=1

    # if index%3909 == 0:
    if index%100 == 0:
      save_results(bridge_array, travel_index_times, int(index/float(3909)))

  test(numeps, lnsas, damaged_bridges, damaged_graph, num_bridges_out, flow, shortest_paths, travel_time, vmt)

def contour(m, levels, nest=False, degrees=False, simplify=True):
    import pkg_resources
    try:
        pkg_resources.require('healpy >= 1.9.0')
    except:
        raise RuntimeError('This function requires healpy >= 1.9.0.')
    try:
        import networkx as nx
    except:
        raise RuntimeError('This function requires the networkx package.')

    # Determine HEALPix resolution
    npix = len(m)
    nside = hp.npix2nside(npix)
    min_area = 0.4 * hp.nside2pixarea(nside)

    # Compute faces, vertices, and neighbors.
    # vertices is an N X 3 array of the distinct vertices of the HEALPix faces.
    # faces is an npix X 4 array mapping HEALPix faces to their vertices.
    # neighbors is an npix X 4 array mapping faces to their nearest neighbors.
    faces = np.ascontiguousarray(
            np.rollaxis(hp.boundaries(nside, np.arange(npix), nest=nest), 2, 1))
    dtype = faces.dtype
    faces = faces.view(np.dtype((np.void, dtype.itemsize * 3)))
    vertices, faces = np.unique(faces.ravel(), return_inverse=True)
    faces = faces.reshape(-1, 4)
    vertices = vertices.view(dtype).reshape(-1, 3)
    neighbors = hp.get_all_neighbours(nside, np.arange(npix), nest=nest)[::2].T

    # Loop over the requested contours.
    paths = []
    for level in levels:

        # Find credible region
        indicator = (m >= level)

        # Construct a graph of the eges of the contour.
        graph = nx.Graph()
        face_pairs = set()
        for ipix1, ipix2 in enumerate(neighbors):
            for ipix2 in ipix2:
                # Determine if we have already considered this pair of faces.
                new_face_pair = frozenset((ipix1, ipix2))
                if new_face_pair in face_pairs: continue
                face_pairs.add(new_face_pair)

                # Determine if this pair of faces are on a boundary of the
                # credible level.
                if indicator[ipix1] == indicator[ipix2]: continue

                # Add all common edges of this pair of faces.
                i1 = np.concatenate((faces[ipix1], [faces[ipix1][0]]))
                i2 = np.concatenate((faces[ipix2], [faces[ipix2][0]]))
                edges1 = frozenset(frozenset(_) for _ in zip(i1[:-1], i1[1:]))
                edges2 = frozenset(frozenset(_) for _ in zip(i2[:-1], i2[1:]))
                for edge in edges1 & edges2:
                    graph.add_edge(*edge)
        graph = nx.freeze(graph)

        # Record a closed path for each cycle in the graph.
        cycles = [
            np.take(vertices, cycle, axis=0) for cycle in nx.cycle_basis(graph)]

        # Simplify paths if requested
        if simplify:
            cycles = [_simplify(cycle, min_area) for cycle in cycles]
            cycles = [cycle for cycle in cycles if len(cycle) > 2]

        # Convert to lists
        cycles = [
            _vec2radec(cycle, degrees=degrees).tolist() for cycle in cycles]

        # Add to output paths
        paths.append([cycle + [cycle[0]] for cycle in cycles])

    return paths

def flatten(item, freeze=True):
    graph = _flatten(item, set())
    if freeze:
        # Frozen graph can't be modified...
        return nx.freeze(graph)
    return graph

 def freeze(self):
     nx.freeze(self)

 def __init__(self, name):
     super(Flow, self).__init__(name)
     self._graph = nx.freeze(nx.DiGraph())