from collections import defaultdict
from ortools.linear_solver.pywraplp import Solver

flow_network = {
    's': [('o', 3), ('p', 2)],
    'o': [('p', 2), ('q', 3)],
    'p': [('r', 2)],
    'q': [('r', 4), ('t', 2)],
    'r': [('t', 3)],
    # Feedback edge (t, s)
    't': [('s', float('inf'))],
}

solver = Solver('maxflows', Solver.GLOP_LINEAR_PROGRAMMING)

# f[u, v] = flow through edge u, v
f = {}
for u in flow_network.keys():
    for (v, capacity_uv) in flow_network[u]:
        f[u, v] = solver.NumVar(0, capacity_uv, f'f({u},{v})')

# flow_in[u] = flow into vertex u
# flow_out[u] = flow out of vertex u
flow_in = defaultdict(int)
flow_out = defaultdict(int)
for u in flow_network.keys():
    for (v, _ ) in flow_network[u]:
        flow_in[v] += f[u, v]
        flow_out[u] += f[u, v]

for u in flow_network.keys():
    solver.Add(flow_in[u] == flow_out[u])

solver.Maximize(f['t','s'])

if solver.Solve() == Solver.OPTIMAL:
    for fe in f.values():
        print(fe, '=', fe.solution_value())