def test_digon(self):
        """Check if digons are handled properly. Taken from ticket
        #618 by arv."""
        nodes = [(1, {}),
                 (2, {'demand': -4}),
                 (3, {'demand': 4}),
                 ]
        edges = [(1, 2, {'capacity': 3, 'weight': 600000}),
                 (2, 1, {'capacity': 2, 'weight': 0}),
                 (2, 3, {'capacity': 5, 'weight': 714285}),
                 (3, 2, {'capacity': 2, 'weight': 0}),
                 ]
        G = nx.DiGraph(edges)
        G.add_nodes_from(nodes)
        flowCost, H = nx.network_simplex(G)
        soln = {1: {2: 0},
                2: {1: 0, 3: 4},
                3: {2: 0}}
        assert_equal(flowCost, 2857140)
        assert_equal(nx.min_cost_flow_cost(G), 2857140)
        assert_equal(H, soln)
        assert_equal(nx.min_cost_flow(G), soln)
        assert_equal(nx.cost_of_flow(G, H), 2857140)

        flowCost, H = nx.capacity_scaling(G)
        assert_equal(flowCost, 2857140)
        assert_equal(H, soln)
        assert_equal(nx.cost_of_flow(G, H), 2857140)

    def test_digraph1(self):
        # From Bradley, S. P., Hax, A. C. and Magnanti, T. L. Applied
        # Mathematical Programming. Addison-Wesley, 1977.
        G = nx.DiGraph()
        G.add_node(1, demand=-20)
        G.add_node(4, demand=5)
        G.add_node(5, demand=15)
        G.add_edges_from([(1, 2, {'capacity': 15, 'weight': 4}),
                          (1, 3, {'capacity': 8, 'weight': 4}),
                          (2, 3, {'weight': 2}),
                          (2, 4, {'capacity': 4, 'weight': 2}),
                          (2, 5, {'capacity': 10, 'weight': 6}),
                          (3, 4, {'capacity': 15, 'weight': 1}),
                          (3, 5, {'capacity': 5, 'weight': 3}),
                          (4, 5, {'weight': 2}),
                          (5, 3, {'capacity': 4, 'weight': 1})])
        flowCost, H = nx.network_simplex(G)
        soln = {1: {2: 12, 3: 8},
                2: {3: 8, 4: 4, 5: 0},
                3: {4: 11, 5: 5},
                4: {5: 10},
                5: {3: 0}}
        assert_equal(flowCost, 150)
        assert_equal(nx.min_cost_flow_cost(G), 150)
        assert_equal(H, soln)
        assert_equal(nx.min_cost_flow(G), soln)
        assert_equal(nx.cost_of_flow(G, H), 150)

        flowCost, H = nx.capacity_scaling(G)
        assert_equal(flowCost, 150)
        assert_equal(H, soln)
        assert_equal(nx.cost_of_flow(G, H), 150)

    def test_zero_capacity_edges(self):
        """Address issue raised in ticket #617 by arv."""
        G = nx.DiGraph()
        G.add_edges_from([(1, 2, {'capacity': 1, 'weight': 1}),
                          (1, 5, {'capacity': 1, 'weight': 1}),
                          (2, 3, {'capacity': 0, 'weight': 1}),
                          (2, 5, {'capacity': 1, 'weight': 1}),
                          (5, 3, {'capacity': 2, 'weight': 1}),
                          (5, 4, {'capacity': 0, 'weight': 1}),
                          (3, 4, {'capacity': 2, 'weight': 1})])
        G.nodes[1]['demand'] = -1
        G.nodes[2]['demand'] = -1
        G.nodes[4]['demand'] = 2

        flowCost, H = nx.network_simplex(G)
        soln = {1: {2: 0, 5: 1},
                2: {3: 0, 5: 1},
                3: {4: 2},
                4: {},
                5: {3: 2, 4: 0}}
        assert_equal(flowCost, 6)
        assert_equal(nx.min_cost_flow_cost(G), 6)
        assert_equal(H, soln)
        assert_equal(nx.min_cost_flow(G), soln)
        assert_equal(nx.cost_of_flow(G, H), 6)

        flowCost, H = nx.capacity_scaling(G)
        assert_equal(flowCost, 6)
        assert_equal(H, soln)
        assert_equal(nx.cost_of_flow(G, H), 6)

def constrained_kmeans(data, demand, maxiter=None, fixedprec=1e9):
    data = np.array(data)

    min_ = np.min(data, axis=0)
    max_ = np.max(data, axis=0)

    c = min_ + np.random.random((len(demand), data.shape[1])) * (max_ - min_)
    m = np.array([-1] * len(data), dtype=np.int)

    itercnt = 0
    while True:
        itercnt += 1

        # memberships
        g = nx.DiGraph()
        g.add_nodes_from(xrange(0, data.shape[0]), demand=-1)  # points
        for i in xrange(0, len(c)):
            g.add_node(len(data) + i, demand=demand[i])

        # Calculating cost...
        cost = np.array([np.linalg.norm(
            np.tile(data.T, len(c)).T - np.tile(c, len(data)).reshape(len(c) * len(data), c.shape[1]), axis=1)])
        # Preparing data_to_c_edges...
        data_to_c_edges = np.concatenate((np.tile([xrange(0, data.shape[0])], len(c)).T,
                                          np.tile(np.array([xrange(data.shape[0], data.shape[0] + c.shape[0])]).T,
                                                  len(data)).reshape(len(c) * len(data), 1), cost.T * fixedprec),
                                         axis=1).astype(np.uint64)
        # Adding to graph
        g.add_weighted_edges_from(data_to_c_edges)

        a = len(data) + len(c)
        g.add_node(a, demand=len(data) - np.sum(demand))
        c_to_a_edges = np.concatenate((np.array([xrange(len(data), len(data) + len(c))]).T, np.tile([[a]], len(c)).T),
                                      axis=1)
        g.add_edges_from(c_to_a_edges)

        # Calculating min cost flow...
        f = nx.min_cost_flow(g)

        # assign
        m_new = np.ones(len(data), dtype=np.int) * -1
        for i in xrange(len(data)):
            p = sorted(f[i].iteritems(), key=lambda x: x[1])[-1][0]
            m_new[i] = p - len(data)

        # stop condition
        if np.all(m_new == m):
            # Stop
            return c, m, f

        m = m_new

        # compute new centers
        for i in xrange(len(c)):
            c[i, :] = np.mean(data[m == i, :], axis=0)

        if maxiter is not None and itercnt >= maxiter:
            # Max iterations reached
            return c, m, f

def assign(photos, pois):
    """Assign photos to POIs with minimum cost"""
    #REF: en.wikipedia.org/wiki/Minimum-cost_flow_problem
    #assert(len(photos) == len(pois))

    dists = np.zeros((len(photos), len(pois)), dtype=np.float64)
    for i, d in enumerate(photos):
        for j, p in enumerate(pois):
            dists[i, j] = round(math.sqrt( (d[0] - p[0])**2 + (d[1] - p[1])**2 ))
    #print(dists)

    G = nx.DiGraph()
 
    # complete bipartite graph: photo -> POI
    # infinity capacity
    for i in range(len(photos)):
        for j in range(len(pois)):
            u = 'd' + str(i)
            v = 'p' + str(j)
            G.add_edge(u, v, weight=dists[i, j])

    # source -> photo
    # capacity = 1
    for i in range(len(photos)):
        u = 's'
        v = 'd' + str(i)
        G.add_edge(u, v, capacity=1, weight=0)

    # POI -> sink
    # infinity capacity
    for j in range(len(pois)):
        u = 'p' + str(j)
        v = 't'
        G.add_edge(u, v, weight=0)

    # demand for source and sink
    G.add_node('s', demand=-len(photos))
    G.add_node('t', demand=len(photos))

    #print(G.nodes())
    #print(G.edges())

    flowDict = nx.min_cost_flow(G)

    assignment = dict()
    for e in G.edges():
        u = e[0]
        v = e[1]
        if u != 's' and v != 't' and flowDict[u][v] > 0:
            #print(e, flowDict[u][v])
            assignment[u] = v
    return assignment

    def test_transshipment(self):
        G = nx.DiGraph()
        G.add_node('a', demand=1)
        G.add_node('b', demand=-2)
        G.add_node('c', demand=-2)
        G.add_node('d', demand=3)
        G.add_node('e', demand=-4)
        G.add_node('f', demand=-4)
        G.add_node('g', demand=3)
        G.add_node('h', demand=2)
        G.add_node('r', demand=3)
        G.add_edge('a', 'c', weight=3)
        G.add_edge('r', 'a', weight=2)
        G.add_edge('b', 'a', weight=9)
        G.add_edge('r', 'c', weight=0)
        G.add_edge('b', 'r', weight=-6)
        G.add_edge('c', 'd', weight=5)
        G.add_edge('e', 'r', weight=4)
        G.add_edge('e', 'f', weight=3)
        G.add_edge('h', 'b', weight=4)
        G.add_edge('f', 'd', weight=7)
        G.add_edge('f', 'h', weight=12)
        G.add_edge('g', 'd', weight=12)
        G.add_edge('f', 'g', weight=-1)
        G.add_edge('h', 'g', weight=-10)
        flowCost, H = nx.network_simplex(G)
        soln = {'a': {'c': 0},
                'b': {'a': 0, 'r': 2},
                'c': {'d': 3},
                'd': {},
                'e': {'r': 3, 'f': 1},
                'f': {'d': 0, 'g': 3, 'h': 2},
                'g': {'d': 0},
                'h': {'b': 0, 'g': 0},
                'r': {'a': 1, 'c': 1}}
        assert_equal(flowCost, 41)
        assert_equal(nx.min_cost_flow_cost(G), 41)
        assert_equal(H, soln)
        assert_equal(nx.min_cost_flow(G), soln)
        assert_equal(nx.cost_of_flow(G, H), 41)

        flowCost, H = nx.capacity_scaling(G)
        assert_equal(flowCost, 41)
        assert_equal(nx.cost_of_flow(G, H), 41)
        assert_equal(H, soln)

 def test_simple_digraph(self):
     G = nx.DiGraph()
     G.add_node('a', demand = -5)
     G.add_node('d', demand = 5)
     G.add_edge('a', 'b', weight = 3, capacity = 4)
     G.add_edge('a', 'c', weight = 6, capacity = 10)
     G.add_edge('b', 'd', weight = 1, capacity = 9)
     G.add_edge('c', 'd', weight = 2, capacity = 5)
     flowCost, H = nx.network_simplex(G)
     soln = {'a': {'b': 4, 'c': 1},
             'b': {'d': 4},
             'c': {'d': 1},
             'd': {}}
     assert_equal(flowCost, 24)
     assert_equal(nx.min_cost_flow_cost(G), 24)
     assert_equal(H, soln)
     assert_equal(nx.min_cost_flow(G), soln)
     assert_equal(nx.cost_of_flow(G, H), 24)

		total += int(exits) - int(entries)

	for n in G.nodes():
		if "demand" not in G.node[n]:
			G.remove_node(n)

	turnstile_stations = [record.strip().split(',')[2] for record in noondata]
	gtfs_stations = G.nodes()

	#print set(turnstile_stations) - set(gtfs_stations)
	extra_nodes = set(gtfs_stations) - set(turnstile_stations)

	nx.draw_spring(G, with_labels=True, node_color='w', node_size=350, font_size=7)
	#plt.show()

	flow = nx.min_cost_flow(G)

	inflowstation=[]
	inflow=[]
	
	for x in G.nodes():
		total = 0
		inflowstation.append(x)
		for p in G.predecessors(x):
			total = total+ flow[p][x]
		inflow.append(total)

    
	rows = zip(inflowstation,inflow)

	length= range(0,len(inflowstation))

        y = [ func(xx) for xx in x ]
        plt.figure()
        plt.plot(x,y)
        plt.show()
    
    
    def FLOWCOST( flow, cost ) :
        res = [ cc( flow.get( e, 0. ) ) for e, cc in cost.iteritems() ]    # big difference!
        #res = [ cost.get( e, line(0.) )( flow[e] ) for e in flow ]
        return sum( res )
    
    def FEAS( flow, capacity, network ) :
        for e in network.edges() :
            if flow.get(e, 0. ) > capacity.get(e, np.Inf ) : return False
        return True
    
    flow = MinConvexCostFlow( g, u, supply, cf, epsilon=.001 )
    print ( flow, FLOWCOST( flow, cf ), FEAS( flow, u, g ) )
    
    flowstar = { 'b' : 10., 'd' : 10. }
    print ( flowstar, FLOWCOST( flowstar, cf ), FEAS( flowstar, u, g ) )
    
    
    digraph = mincostflow_nx( g, u, supply, c )
    compare = nx.min_cost_flow( digraph )
    
    
    
    
    

# Let's now solve for the min-cost flow of the example we had in class. We have already added capacities and demands to our network. We finally need to introduce edge costs (weights).

# In[10]:

G.edge['s']['u']['weight'] = 1
G.edge['s']['v']['weight'] = 1
G.edge['u']['v']['weight'] = 1
G.edge['u']['t']['weight'] = 1
G.edge['v']['t']['weight'] = 2


# Now we can just call the min-cost flow function. Note that we could have used this to solve the flow-with-demands problem, just by setting all the costs (weights) to 0.

# In[11]:

minimum_cost_flow = nx.min_cost_flow(G)
print_flow(minimum_cost_flow)


# This is the same as the maximum flow found in class. For this example, the maximum flow is unique.

# ## Exercises
# 
# 1. Construct the bipartite graph from the example application in the slides and find a matching using max flow.
# 
# 2. Write a function to compute a conformal decomposition of a flow with demands, and run it on the flows found in the "Min-Cost Flow" section above.

# ### Solution 1
# 
# For the first problem, we form the network by hand then solve max flow using new source and sink vertices, imposing a capacity of 1 on all edges.

def constrcut_flow_graph(reach_graph,sorted_bbox,flow_graph=None,scores=None):
    def to_left_node(index):
        return 2*index+1

    def to_right_node(index):
        return 2*index+2

    def scale_to_int(score):
        MAX_VALUE = 10000
        return int(score*MAX_VALUE)

    def compute_box_center(box):
        center_y = (box[3]+box[1])/2
        center_x = (box[2]+box[0])/2
        return center_x,center_y

    def get_data_cost(score):
        return scale_to_int(score)

    def get_smoothness_cost(box1, box2):
        alpha=0.0
        h1=box1[3]-box1[1]+1
        h2=box2[3]-box2[1]+1
        x1, y1=compute_box_center(box1)
        x2, y2=compute_box_center(box2)
        d=((x1-x2)**2+(y1-y2)**2)**0.5/min(h1,h2)
        s=float(abs(h1-h2))/min(h1, h2)
        return scale_to_int(alpha*d+(1-alpha)*s)

    def get_entry_cost(index,reach_graph,scores):
        reach_node_costs = [scores[left] for left, right in reach_graph.edges() if right == index]
        sorted_costs = sorted(reach_node_costs)
        if len(sorted_costs) > 0:
            cost = scale_to_int(-sorted_costs[-1])
        else:
            cost = 0
        return cost

    def get_exit_cost(index,reach_graph,scores):
        reach_node_costs = [scores[right] for left, right in reach_graph.edges() if left == index]
        sorted_costs = sorted(reach_node_costs)
        if len(sorted_costs) > 0:
            cost = scale_to_int(-sorted_costs[-1])
        else:
            cost = 0
        return cost

    def overlaps(box1,box2):
        h1=box1[3]-box1[1]+1
        h2=box2[3]-box2[1]+1
        overlaps=min(box2[3], box1[3])-max(box1[1], box2[1])
        overlaps=overlaps if overlaps>0 else 0
        return float(overlaps)/h2

    length = len(sorted_bbox)
    scores = [-1 for i in range(length)]

    if flow_graph == None:
        #box
        print 'construct graph'
        flow_graph = nx.DiGraph()
        ENTRY_NODE = -1
        flow_graph.add_node(ENTRY_NODE,demand = -1)
        EXIT_NODE = -2
        flow_graph.add_node(EXIT_NODE,demand = 1)
        # data cost
        for index in range(length):
            left_node = to_left_node(index)
            right_node = to_right_node(index)
            flow_graph.add_node(left_node)
            flow_graph.add_node(right_node)
            data_cost = get_data_cost(scores[index])
            flow_graph.add_edge(left_node,right_node,weight=data_cost,capacity = 1)
        # smoothness cost
        for left,right in reach_graph.edges():
            left_node = to_right_node(left)
            right_node = to_left_node(right)
            smoothness_cost = get_smoothness_cost(sorted_bbox[left],sorted_bbox[right])
            flow_graph.add_edge(left_node,right_node,weight=smoothness_cost,capacity = 1)
        # entry cost
        for index in range(length):
            entry_cost = get_entry_cost(index,reach_graph,scores)
            left_node = to_left_node(index)
            flow_graph.add_edge(ENTRY_NODE,left_node,weight=entry_cost,capacity = 1)
        # exit cost
        for index in range(length):
            exit_cost = get_exit_cost(index,reach_graph,scores)
            right_node = to_right_node(index)
            flow_graph.add_edge(right_node,EXIT_NODE,weight=exit_cost,capacity = 1)
    box_inds=[]
    flowDict = nx.min_cost_flow(flow_graph)
    flowCost = nx.min_cost_flow_cost(flow_graph)
    for index in range(length):
        left_node=to_left_node(index)
        right_node=to_right_node(index)
        try:
            find = [node for node, value in flowDict[left_node].iteritems() if value == 1]
        except:
            find = []
        if len(find)>0:
#.........这里部分代码省略.........