diff --git a/.travis.yml b/.travis.yml
index aefddee..86015b5 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -11,11 +11,8 @@ notifications:
 git:
   depth: 99999999
 
-## uncomment the following lines to allow failures on nightly julia
-## (tests will run but not make your overall status red)
-#matrix:
-#  allow_failures:
-#  - julia: nightly
+ allow_failures:
+ - julia: nightly
 
 ## uncomment and modify the following lines to manually install system packages
 #addons:
diff --git a/src/LightGraphsFlows.jl b/src/LightGraphsFlows.jl
index 476beb4..fe5ac3c 100644
--- a/src/LightGraphsFlows.jl
+++ b/src/LightGraphsFlows.jl
@@ -1,5 +1,24 @@
 module LightGraphsFlows
 
-# package code goes here
+import LightGraphs
+const lg = LightGraphs
+using SimpleTraits: @traitfn, @traitimpl
+import SimpleTraits
 
-end # module
+import Base: getindex, size, transpose, ctranspose
+
+include("maximum_flow.jl")
+include("edmonds_karp.jl")
+include("dinic.jl")
+include("boykov_kolmogorov.jl")
+include("push_relabel.jl")
+include("multiroute_flow.jl")
+include("kishimoto.jl")
+include("ext_multiroute_flow.jl")
+
+# flow
+export
+maximum_flow, EdmondsKarpAlgorithm, DinicAlgorithm, BoykovKolmogorovAlgorithm, PushRelabelAlgorithm,
+multiroute_flow, KishimotoAlgorithm, ExtendedMultirouteFlowAlgorithm
+
+end
diff --git a/src/boykov_kolmogorov.jl b/src/boykov_kolmogorov.jl
new file mode 100644
index 0000000..f904454
--- /dev/null
+++ b/src/boykov_kolmogorov.jl
@@ -0,0 +1,200 @@
+"""
+    boykov_kolmogorov_impl(residual_graph, source, target, capacity_matrix)
+
+Compute the max-flow/min-cut between `source` and `target` for `residual_graph`
+using the Boykov-Kolmogorov algorithm.
+
+Return the maximum flow in the network, the flow matrix and the partition
+`{S,T}` in the form of a vector of 0's, 1's and 2's.
+
+### References
+- BOYKOV, Y.; KOLMOGOROV, V., 2004. An Experimental Comparison of
+Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision.
+
+### Author
+- Júlio Hoffimann Mendes (juliohm@stanford.edu)
+"""
+function boykov_kolmogorov_impl end
+# see https://github.com/mauro3/SimpleTraits.jl/issues/47#issuecomment-327880153 for syntax
+@traitfn function boykov_kolmogorov_impl{T, U, AG<:lg.AbstractGraph{U}}(
+    residual_graph::AG::lg.IsDirected,          # the input graph
+    source::Integer,                           # the source vertex
+    target::Integer,                           # the target vertex
+    capacity_matrix::AbstractMatrix{T}    # edge flow capacities
+    )
+    n = lg.nv(residual_graph)
+
+    flow = 0
+    flow_matrix = zeros(T, n, n)
+
+    TREE = zeros(U, n)
+    TREE[source] = U(1)
+    TREE[target] = U(2)
+
+    PARENT = zeros(U, n)
+
+    A = [source, target]
+    O = Vector{U}()
+
+    while true
+        # growth stage
+        path = find_path!(residual_graph, source, target, flow_matrix, capacity_matrix, PARENT, TREE, A)
+
+        isempty(path) && break
+
+        # augmentation stage
+        flow += augment!(path, flow_matrix, capacity_matrix, PARENT, TREE, O)
+
+        # adoption stage
+        adopt!(residual_graph, source, target, flow_matrix, capacity_matrix, PARENT, TREE, A, O)
+    end
+
+    return flow, flow_matrix, TREE
+end
+
+# see https://github.com/mauro3/SimpleTraits.jl/issues/47#issuecomment-327880153 for syntax
+@traitfn function find_path!{T, AG<:lg.AbstractGraph{T}}(
+        residual_graph::AG::lg.IsDirected,            # the input graph
+        source::Integer,                       # the source vertex
+        target::Integer,                       # the target vertex
+        flow_matrix::AbstractMatrix,       # the current flow matrix
+        capacity_matrix::AbstractMatrix,   # edge flow capacities
+        PARENT::Vector,                     # parent table
+        TREE::Vector,                       # tree table
+        A::Vector                           # active set
+    )
+    tree_cap(p, q) = TREE[p] == one(T) ? capacity_matrix[p, q] - flow_matrix[p, q] :
+    capacity_matrix[q, p] - flow_matrix[q, p]
+    while !isempty(A)
+        p = last(A)
+        for q in lg.neighbors(residual_graph, p)
+            if tree_cap(p, q) > 0
+                if TREE[q] == zero(T)
+                    TREE[q] = TREE[p]
+                    PARENT[q] = p
+                    unshift!(A, q)
+                end
+                if TREE[q] ≠ zero(T) && TREE[q] ≠ TREE[p]
+                    # p -> source
+                    path_to_source = [p]
+                    while PARENT[p] ≠ zero(T)
+                        p = PARENT[p]
+                        push!(path_to_source, p)
+                    end
+
+                    # q -> target
+                    path_to_target = [q]
+                    while PARENT[q] ≠ zero(T)
+                        q = PARENT[q]
+                        push!(path_to_target, q)
+                    end
+
+                    # source -> target
+                    path = [reverse!(path_to_source); path_to_target]
+
+                    if path[1] == source && path[end] == target
+                        return path
+                    elseif path[1] == target && path[end] == source
+                        return reverse!(path)
+                    end
+                end
+            end
+        end
+        pop!(A)
+    end
+
+    return Vector{T}()
+end
+
+function augment!(
+    path::AbstractVector,                  # path from source to target
+    flow_matrix::AbstractMatrix,       # the current flow matrix
+    capacity_matrix::AbstractMatrix,   # edge flow capacities
+    PARENT::Vector,                     # parent table
+    TREE::Vector,                       # tree table
+    O::Vector                           # orphan set
+    )
+
+    T = eltype(path)
+    # bottleneck capacity
+    Δ = Inf
+    for i = 1:(length(path) - 1)
+        p, q = path[i:(i + 1)]
+        cap = capacity_matrix[p, q] - flow_matrix[p, q]
+        cap < Δ && (Δ = cap)
+    end
+
+    # update residual graph
+    for i = 1:(length(path) - 1)
+        p, q = path[i:(i + 1)]
+        flow_matrix[p, q] += Δ
+        flow_matrix[q, p] -= Δ
+
+        # create orphans
+        if flow_matrix[p, q] == capacity_matrix[p, q]
+            if TREE[p] == TREE[q] == one(T)
+                PARENT[q] = zero(T)
+                unshift!(O, q)
+            end
+            if TREE[p] == TREE[q] == 2
+                PARENT[p] = zero(T)
+                unshift!(O, p)
+            end
+        end
+    end
+
+    return Δ
+end
+
+@traitfn function adopt!{T, AG<:lg.AbstractGraph{T}}(
+        residual_graph::AG::lg.IsDirected,               # the input graph
+        source::Integer,                             # the source vertex
+        target::Integer,                             # the target vertex
+        flow_matrix::AbstractMatrix,       # the current flow matrix
+        capacity_matrix::AbstractMatrix,   # edge flow capacities
+        PARENT::Vector,                     # parent table
+        TREE::Vector,                       # tree table
+        A::Vector,                             # active set
+        O::Vector                           # orphan set
+    )
+    tree_cap(p, q) = TREE[p] == 1 ? capacity_matrix[p, q] - flow_matrix[p, q] :
+    capacity_matrix[q, p] - flow_matrix[q, p]
+    while !isempty(O)
+        p = pop!(O)
+        # try to find parent that is not an orphan
+        parent_found = false
+        for q in lg.neighbors(residual_graph, p)
+            if TREE[q] == TREE[p] && tree_cap(q, p) > 0
+                # check if "origin" is either source or target
+                o = q
+                while PARENT[o] ≠ 0
+                    o = PARENT[o]
+                end
+                if o == source || o == target
+                    parent_found = true
+                    PARENT[p] = q
+                    break
+                end
+            end
+        end
+
+        if !parent_found
+            # scan all neighbors and make the orphan a free node
+            for q in lg.neighbors(residual_graph, p)
+                if TREE[q] == TREE[p]
+                    if tree_cap(q, p) > 0
+                        unshift!(A, q)
+                    end
+                    if PARENT[q] == p
+                        PARENT[q] = zero(T)
+                        unshift!(O, q)
+                    end
+                end
+            end
+
+            TREE[p] = zero(T)
+            B = setdiff(A, p)
+            resize!(A, length(B))[:] = B
+        end
+    end
+end
diff --git a/src/dinic.jl b/src/dinic.jl
new file mode 100644
index 0000000..514bd08
--- /dev/null
+++ b/src/dinic.jl
@@ -0,0 +1,118 @@
+"""
+    function dinic_impl(residual_graph, source, target, capacity_matrix)
+
+Compute the maximum flow between the `source` and `target` for `residual_graph`
+with edge flow capacities in `capacity_matrix` using
+[Dinic\'s Algorithm](https://en.wikipedia.org/wiki/Dinic%27s_algorithm).
+Return the value of the maximum flow as well as the final flow matrix.
+"""
+function dinic_impl end
+@traitfn function dinic_impl{T}(
+    residual_graph::::lg.IsDirected,               # the input graph
+    source::Integer,                       # the source vertex
+    target::Integer,                       # the target vertex
+    capacity_matrix::AbstractMatrix{T}    # edge flow capacities
+    )
+    n = lg.nv(residual_graph)                     # number of vertexes
+    flow_matrix = zeros(T, n, n)           # initialize flow matrix
+    P = zeros(Int, n)                      # Sharable parent vector
+
+    flow = 0
+
+    while true
+        augment = blocking_flow!(residual_graph, source, target, capacity_matrix, flow_matrix, P)
+        augment == 0 && break
+        flow += augment
+    end
+    return flow, flow_matrix
+end
+
+
+
+
+"""
+    blocking_flow!(residual_graph, source, target, capacity_matrix, flow-matrix, P)
+
+Like `blocking_flow`, but requires a preallocated parent vector `P`.
+"""
+function blocking_flow! end
+@traitfn function blocking_flow!{T}(
+    residual_graph::::lg.IsDirected,               # the input graph
+    source::Integer,                           # the source vertex
+    target::Integer,                           # the target vertex
+    capacity_matrix::AbstractMatrix{T},   # edge flow capacities
+    flow_matrix::AbstractMatrix,       # the current flow matrix
+    P::AbstractVector{Int}                 # Parent vector to store Level Graph
+    )
+    n = lg.nv(residual_graph)                     # number of vertexes
+    fill!(P, -1)
+    P[source] = -2
+
+    Q = [source]
+    sizehint!(Q, n)
+
+    while length(Q) > 0                   # Construct the Level Graph using BFS
+        u = pop!(Q)
+        for v in lg.out_neighbors(residual_graph, u)
+            if P[v] == -1 && capacity_matrix[u, v] > flow_matrix[u, v]
+                P[v] = u
+                unshift!(Q, v)
+            end
+        end
+    end
+
+    P[target] == -1 && return 0                    # BFS couldn't reach the target
+
+    total_flow = 0
+
+    for bv in lg.in_neighbors(residual_graph, target)    # Trace all possible routes to source
+        flow = typemax(T)
+        v = target
+        u = bv
+        while v != source
+            if u == -1                    # Vertex unreachable from source
+                flow = 0
+                break
+            else
+                flow = min(flow, capacity_matrix[u, v] - flow_matrix[u, v])
+                v = u
+                u = P[u]
+            end
+        end
+
+        flow == 0 && continue                      # Flow cannot be augmented along path
+
+        v = target
+        u = bv
+        while v != source             # Augment flow along path
+            flow_matrix[u, v] += flow
+            flow_matrix[v, u] -= flow
+            v = u
+            u = P[u]
+        end
+
+        total_flow += flow
+    end
+    return total_flow
+end
+
+"""
+    blocking_flow(residual_graph, source, target, capacity_matrix, flow-matrix)
+
+Use BFS to identify a blocking flow in the `residual_graph` with current flow
+matrix `flow_matrix`and then backtrack from `target` to `source`,
+augmenting flow along all possible paths.
+"""
+blocking_flow(
+    residual_graph::lg.AbstractGraph,               # the input graph
+    source::Integer,                       # the source vertex
+    target::Integer,                       # the target vertex
+    capacity_matrix::AbstractMatrix,   # edge flow capacities
+    flow_matrix::AbstractMatrix,       # the current flow matrix
+    ) = blocking_flow!(
+            residual_graph,
+            source,
+            target,
+            capacity_matrix,
+            flow_matrix,
+            zeros(Int, lg.nv(residual_graph)))
diff --git a/src/edmonds_karp.jl b/src/edmonds_karp.jl
new file mode 100644
index 0000000..9a99852
--- /dev/null
+++ b/src/edmonds_karp.jl
@@ -0,0 +1,162 @@
+"""
+    edmonds_karp_impl(residual_graph, source, target, capacity_matrix)
+
+Compute the maximum flow in flow graph `residual_graph` between `source` and
+`target` and capacities defined in `capacity_matrix` using the
+[Edmonds-Karp algorithm](https://en.wikipedia.org/wiki/Edmondss%E2%80%93Karp_algorithm).
+Return the value of the maximum flow as well as the final flow matrix.
+"""
+function edmonds_karp_impl end
+@traitfn function edmonds_karp_impl{T}(
+    residual_graph::::lg.IsDirected,               # the input graph
+    source::Integer,                       # the source vertex
+    target::Integer,                       # the target vertex
+    capacity_matrix::AbstractMatrix{T}    # edge flow capacities
+    )
+    n = lg.nv(residual_graph)                     # number of vertexes
+    flow = 0
+    flow_matrix = zeros(T, n, n)           # initialize flow matrix
+    P = zeros(Int, n)
+    S = zeros(Int, n)
+    while true
+        fill!(P, -1)
+        fill!(S, -1)
+        v, P, S, flag = fetch_path!(residual_graph, source, target, flow_matrix, capacity_matrix, P, S)
+
+        if flag != 0                       # no more valid paths
+            break
+        else
+            path = [Int(v)]                     # initialize path
+            sizehint!(path, n)
+
+            u = v
+            while u != source                # trace path from v to source
+                u = P[u]
+                push!(path, u)
+            end
+            reverse!(path)
+
+            u = v                          # trace path from v to target
+            while u != target
+                u = S[u]
+                push!(path, Int(u))
+            end
+            # augment flow along path
+            flow += augment_path!(path, flow_matrix, capacity_matrix)
+        end
+    end
+
+    return flow, flow_matrix
+end
+
+"""
+    augment_path!(path, flow_matrix, capacity_matrix)
+
+Calculate the amount by which flow can be augmented in the given path.
+Augment the flow and returns the augment value.
+"""
+function augment_path!(
+    path::Vector{Int},                     # input path
+    flow_matrix::AbstractMatrix{T},       # the current flow matrix
+    capacity_matrix::AbstractMatrix    # edge flow capacities
+    ) where T
+    augment = typemax(T)                   # initialize augment
+    for i in 1:(length(path) - 1)              # calculate min capacity along path
+        u = path[i]
+        v = path[i + 1]
+        augment = min(augment, capacity_matrix[u, v] - flow_matrix[u, v])
+    end
+
+    for i in 1:(length(path) - 1)              # augment flow along path
+        u = path[i]
+        v = path[i + 1]
+        flow_matrix[u, v] += augment
+        flow_matrix[v, u] -= augment
+    end
+
+    return augment
+end
+
+"""
+    fetch_path!(residual_graph, source, target, flow_matrix, capacity_matrix, P, S)
+
+Like `fetch_path`, but requires preallocated parent vector `P` and successor
+vector `S`.
+"""
+function fetch_path! end
+@traitfn function fetch_path!(
+    residual_graph::::lg.IsDirected,               # the input graph
+    source::Integer,                       # the source vertex
+    target::Integer,                       # the target vertex
+    flow_matrix::AbstractMatrix,       # the current flow matrix
+    capacity_matrix::AbstractMatrix,   # edge flow capacities
+    P::Vector{Int},                        # parent table of path init to -1s
+    S::Vector{Int}                         # successor table of path init to -1s
+    )
+    n = lg.nv(residual_graph)
+
+    P[source] = -2
+    S[target] = -2
+
+    Q_f = [source]                         # forward queue
+    sizehint!(Q_f, n)
+
+    Q_r = [target]                         # reverse queue
+    sizehint!(Q_r, n)
+
+    while true
+        if length(Q_f) <= length(Q_r)
+            u = pop!(Q_f)
+            for v in lg.out_neighbors(residual_graph, u)
+                if capacity_matrix[u, v] - flow_matrix[u, v] > 0 && P[v] == -1
+                    P[v] = u
+                    if S[v] == -1
+                        unshift!(Q_f, v)
+                    else
+                        return v, P, S, 0 # 0 indicates success
+                    end
+                end
+            end
+
+            length(Q_f) == 0 && return 0, P, S, 1 # No paths to target
+        else
+            v = pop!(Q_r)
+            for u in lg.in_neighbors(residual_graph, v)
+                if capacity_matrix[u, v] - flow_matrix[u, v] > 0 && S[u] == -1
+                    S[u] = v
+                    P[u] != -1 && return u, P, S, 0 # 0 indicates success
+
+                    unshift!(Q_r, u)
+                end
+
+            end
+
+            length(Q_r) == 0 && return 0, P, S, 2 # No paths to source
+        end
+    end
+end
+
+
+"""
+    fetch_path(residual_graph, source, target, flow_matrix, capacity_matrix)
+
+
+Use bidirectional BFS to look for augmentable paths from `source` to `target` in
+`residual_graph`. Return the vertex where the two BFS searches intersect,
+the parent table of the path, the successor table of the path found, and a
+flag indicating success (0 => success; 1 => no path to target, 2 => no path
+to source).
+"""
+function fetch_path end
+@traitfn function fetch_path(
+    residual_graph::::lg.IsDirected,               # the input graph
+    source::Integer,                           # the source vertex
+    target::Integer,                           # the target vertex
+    flow_matrix::AbstractMatrix,       # the current flow matrix
+    capacity_matrix::AbstractMatrix    # edge flow capacities
+    )
+    n = lg.nv(residual_graph)
+    P = fill(-1, n)
+    S = fill(-1, n)
+    return fetch_path!(residual_graph, source, target, flow_matrix, capacity_matrix, P, S)
+end
diff --git a/src/ext_multiroute_flow.jl b/src/ext_multiroute_flow.jl
new file mode 100644
index 0000000..095dc1f
--- /dev/null
+++ b/src/ext_multiroute_flow.jl
@@ -0,0 +1,252 @@
+"""
+    emrf(flow_graph, source, target, capacity_matrix, flow_algorithm, routes=0)
+
+Compute the maximum multiroute flow (for any number of `route`s)
+between `source` and `target` in `flow_graph` via flow algorithm `flow_algorithm`.
+
+If a number of routes is given, return the value of the multiroute flow as
+well as the final flow matrix, along with a multiroute cut if the
+Boykov-Kolmogorov max-flow algorithm is used as a subroutine.
+Otherwise, return the vector of breaking points of the parametric
+multiroute flow function.
+
+### References
+- [Extended Multiroute Flow algorithm](http://dx.doi.org/10.1016/j.disopt.2016.05.002)
+"""
+# EMRF (Extended Multiroute Flow) algorithms
+function emrf(
+        flow_graph::lg.AbstractGraph,          # the input graph
+        source::Integer,                       # the source vertex
+        target::Integer,                       # the target vertex
+        capacity_matrix::AbstractMatrix,       # edge flow capacities
+        flow_algorithm::AbstractFlowAlgorithm, # keyword argument for algorithm
+        routes::Real = 0
+    )
+    breakingpoints = breakingPoints(flow_graph, source, target, capacity_matrix)
+    if routes > 0
+        x, f = intersection(breakingpoints, routes)
+        return maximum_flow(flow_graph, source, target, capacity_matrix,
+        algorithm = flow_algorithm, restriction = x)
+    end
+    return breakingpoints
+end
+
+@doc_str """
+    auxiliaryPoints(flow_graph, source, target, capacity_matrix)
+
+Output a set of (point, slope) that compose the restricted max-flow function
+of `flow_graph` from `source to `target` using capacities in `capacity_matrix`.
+
+### Performance
+One point by possible slope is enough (hence ``\\mathcal{O}(λ×maximum_flow)`` complexity).
+"""
+function auxiliaryPoints end
+@traitfn function auxiliaryPoints(
+    flow_graph::::lg.IsDirected,                   # the input graph
+    source::Integer,                           # the source vertex
+    target::Integer,                           # the target vertex
+    capacity_matrix::AbstractMatrix   # edge flow capacities
+    )
+    # Problem descriptors
+    λ = maximum_flow(flow_graph, source, target)[1] # Connectivity
+    n = lg.nv(flow_graph) # number of vertices
+    r1, r2 = minmaxCapacity(capacity_matrix) # restriction left (1) and right (2)
+    auxpoints = fill((0., 0.), λ + 1)
+
+    # Initialisation of left side (1)
+    f1, F1, cut1 = maximum_flow(flow_graph, source, target, capacity_matrix,
+    algorithm = BoykovKolmogorovAlgorithm(), restriction = r1)
+    s1 = slope(flow_graph, capacity_matrix, cut1, r1) # left slope
+    auxpoints[λ + 1 - s1] = (r1, f1) # Add left initial auxiliary point
+
+    # Initialisation of right side (2)
+    f2, F2, cut2 = maximum_flow(flow_graph, source, target, capacity_matrix,
+    algorithm = BoykovKolmogorovAlgorithm(), restriction = r2)
+    s2 = slope(flow_graph, capacity_matrix, cut2, r2) # right slope
+    auxpoints[λ + 1 - s2] = (r2, f2) # Add right initial auxiliary point
+
+    # Loop if the slopes are distinct by at least 2
+    if s1 > s2 + 1
+        queue = [((f1, s1, r1), (f2, s2, r2))]
+
+        while !isempty(queue)
+            # Computes an intersection (middle) with a new associated slope
+            (f1, s1, r1), (f2, s2, r2) = pop!(queue)
+            r, expectedflow = intersection(r1, f1, s1, r2, f2, s2)
+            f, F, cut = maximum_flow(flow_graph, source, target, capacity_matrix,
+            algorithm = BoykovKolmogorovAlgorithm(), restriction = r)
+            s = slope(flow_graph, capacity_matrix, max.(cut, 1), r) # current slope
+            auxpoints[λ + 1 - s] = (r, f)
+            # If the flow at the intersection (middle) is as expected
+            if expectedflow ≉ f # approximatively not equal (enforced by floating precision)
+                # if the slope difference between (middle) and left is at least 2
+                # push (left),(middle)
+                if s1 > s + 1 && !approximately_equal((r2, f2), (r, f))
+                    q = (f1, s1, r1), (f, s, r)
+                    push!(queue, q)
+                end
+                # if the slope difference between (middle) and right is at least 2
+                # push (middle),(right)
+                if s > s2 + 1 && !approximately_equal((r1, f1), (r, f))
+                    q = (f, s, r), (f2, s2, r2)
+                    push!(queue, q)
+                end
+            end
+        end
+    end
+    return auxpoints
+end
+
+"""
+    breakingPoints(flow_graph::::IsDirected, source, target, capacity_matrix)
+
+Calculates the breaking of the restricted max-flow from a set of auxiliary points
+for `flow_graph` from `source to `target` using capacities in `capacity_matrix`.
+"""
+function breakingPoints end
+@traitfn function breakingPoints{T}(
+        flow_graph::::lg.IsDirected,                   # the input graph
+        source::Integer,                           # the source vertex
+        target::Integer,                           # the target vertex
+        capacity_matrix::AbstractMatrix{T}   # edge flow capacities
+    )
+    auxpoints = auxiliaryPoints(flow_graph, source, target, capacity_matrix)
+    λ = length(auxpoints) - 1
+    left_index = 1
+    breakingpoints = Vector{Tuple{T,T,Int}}()
+
+    for (id, point) in enumerate(auxpoints)
+        if id == 1
+            push!(breakingpoints, (0., 0., λ))
+        else
+            pleft = breakingpoints[left_index]
+            if point[1] != 0
+                x, y = intersection(pleft[1], pleft[2], pleft[3],
+                point[1], point[2], λ + 1 - id)
+                push!(breakingpoints, (x, y, λ + 1 - id))
+                left_index += 1
+            end
+        end
+    end
+    return breakingpoints
+end
+
+"""
+    minmaxCapacity(capacity_matrix)
+
+Return the nonzero min and max function of `capacity_matrix`.
+
+Note: this is more efficient than maximum() / minimum() / extrema()
+since we have to ignore zero values.
+"""
+
+# Function to get the nonzero min and max function of a Matrix
+# note: this is more efficient than maximum() / minimum() / extrema()
+#       since we have to ignore zero values.
+function minmaxCapacity(
+    capacity_matrix::AbstractMatrix{T}    # edge flow capacities
+    ) where T
+    cmin, cmax = typemax(T), typemin(T)
+    for c in capacity_matrix
+        if c > zero(T)
+            cmin = min(cmin, c)
+        end
+        cmax = max(cmax, c)
+    end
+    return cmin, cmax
+end
+
+"""
+    slope(flow_graph, capacity_matrix, cut, restriction)
+
+Return the slope of `flow_graph` using capacities in `capacity_matrix` and
+a cut vector `cut`. The slope is initialized at 0 and is incremented for
+each edge whose capacity does not exceed `restriction`.
+"""
+function slope end
+# Function to get the slope of the restricted flow
+@traitfn function slope(
+    flow_graph::::lg.IsDirected,                   # the input graph
+    capacity_matrix::AbstractMatrix,  # edge flow capacities
+    cut::Vector,                      # cut information for vertices
+    restriction::Number                 # value of the restriction
+    )
+    slope = 0
+    for e in lg.edges(flow_graph)
+        ## Chain comparison to wether an edge cross the cut from the source side of
+        # the cut to the target side of the cut. Then the edge is selected iff the
+        # capacity of the edge is larger then the restriction argument.
+        # cut[dst(e)] == 2 > cut[src(e)] is equivalent to
+        # cut[dst(e)] == 2 && 2 > cut[src(e)]
+        # Description of chain comparisons can be found at https://goo.gl/IJpCqe
+        if cut[lg.dst(e)] == 2 > cut[lg.src(e)] &&
+            capacity_matrix[lg.src(e), lg.dst(e)] > restriction
+            slope += 1
+        end
+    end
+    return slope
+end
+
+"""
+        intersection(x1, y1, a1, x2, y2, a2)
+
+Return the intersection of two lines defined by `x` and `y` with slopes `a`.
+2) A set of segments and a linear function of slope k passing by the origin.
+Requires argument:
+1) - x1, y1, a1, x2, y2, a2::T<:AbstractFloat # Coordinates/slopes
+2) - points::Vector{Tuple{T, T, Int}}         # vector of points with T<:AbstractFloat
+- k::R<:Real                             # number of routes (slope of the line)
+"""
+function intersection(
+    x1::T,          # x coordinate of point 1
+    y1::T,          # y coordinate of point 1
+    a1::Integer,        # slope passing by point 1
+    x2::T,          # x coordinate of point 2
+    y2::T,          # y coordinate of point 2
+    a2::R           # slope passing by point 2
+    ) where T<:AbstractFloat where R<:Real
+
+    (a1 == a2) && return -1., -1. # result will be ignored in other intersection method
+    b1 = y1 - a1 * x1
+    b2 = y2 - a2 * x2
+    x = (b2 - b1) / (a1 - a2)
+    y = a1 * x + b1
+    return x, y
+end
+
+
+"""
+    intersection(points, k)
+
+Return the intersection of a set of line segments and a line of slope `k`
+passing by the origin. Segments are defined as a triple (x, y, slope).
+"""
+function intersection(
+    points::Vector{Tuple{T,T,I}},  # vector of breaking points
+    k::R                               # number of routes (slope of the line)
+    ) where T<:AbstractFloat where I<:Integer where R<:Real
+    λ = points[1][1] # Connectivity
+
+    # Loop over the segments (pair of breaking points)
+    for (id, p) in enumerate(points[1:(end - 1)])
+        if id == 1
+            k ≈ λ && return points[2]
+        else
+            x, y = intersection(p[1], p[2], p[3], 0., 0., k)
+            (p[1] ≤ x ≤ points[id + 1][1]) && return x, y
+        end
+    end
+    p = points[end]
+    return intersection(p[1], p[2], p[3], 0., 0., k)
+end
+
+"""
+    approximately_equal(a, b)
+
+Return true if each element in the tuple is approximately equal to its counterpart.
+
+### Implementation Notes:
+This is a separate function because we don't want to hijack isapprox for tuples.
+"""
+approximately_equal(a::Tuple{T,T}, b::Tuple{T,T}) where T <: AbstractFloat =
+    a[1] ≈ b[1] && a[2] ≈ b[2]
diff --git a/src/kishimoto.jl b/src/kishimoto.jl
new file mode 100644
index 0000000..2a2d8ff
--- /dev/null
+++ b/src/kishimoto.jl
@@ -0,0 +1,70 @@
+# Method when using Boykov-Kolmogorov as a subroutine
+# Kishimoto algorithm
+
+@traitfn function kishimoto(
+    flow_graph::::lg.IsDirected,                       # the input graph
+    source::Integer,                           # the source vertex
+    target::Integer,                           # the target vertex
+    capacity_matrix::AbstractMatrix,      # edge flow capacities
+    flow_algorithm::BoykovKolmogorovAlgorithm, # keyword argument for algorithm
+    routes::Int                                # keyword argument for routes
+    )
+    # Initialisation
+    flow, F, labels = maximum_flow(flow_graph, source, target,
+    capacity_matrix, algorithm = flow_algorithm)
+    restriction = flow / routes
+    flow, F, labels = maximum_flow(flow_graph, source, target, capacity_matrix,
+    algorithm = flow_algorithm, restriction = restriction)
+
+    # Loop condition : approximatively not equal is enforced by floating precision
+    i = 1
+    while flow < routes * restriction && flow ≉ routes * restriction
+        restriction = (flow - i * restriction) / (routes - i)
+        i += 1
+        flow, F, labels = maximum_flow(flow_graph, source, target, capacity_matrix,
+        algorithm = flow_algorithm, restriction = restriction)
+    end
+
+    # End
+    return flow, F, labels
+end
+
+
+"""
+    kishimoto(flow_graph, source, target, capacity_matrix, flow_algorithm, routes)
+
+Compute the maximum multiroute flow (for an integer number of `route`s)
+between `source` and `target` in `flow_graph` with capacities in `capacity_matrix`
+using the [Kishimoto algorithm](http://dx.doi.org/10.1109/ICCS.1992.255031).
+Return the value of the multiroute flow as well as the final flow matrix,
+along with a multiroute cut if Boykov-Kolmogorov is used as a subroutine.
+"""
+function kishimoto end
+@traitfn function kishimoto(
+    flow_graph::::lg.IsDirected,                   # the input graph
+    source::Integer,                       # the source vertex
+    target::Integer,                       # the target vertex
+    capacity_matrix::AbstractMatrix,  # edge flow capacities
+    flow_algorithm::AbstractFlowAlgorithm, # keyword argument for algorithm
+    routes::Int                            # keyword argument for routes
+    )
+    # Initialisation
+    flow, F = maximum_flow(flow_graph, source, target,
+    capacity_matrix, algorithm = flow_algorithm)
+    restriction = flow / routes
+
+    flow, F = maximum_flow(flow_graph, source, target, capacity_matrix,
+    algorithm = flow_algorithm, restriction = restriction)
+
+    # Loop condition : approximatively not equal is enforced by floating precision
+    i = 1
+    while flow < routes * restriction && flow ≉ routes * restriction
+        restriction = (flow - i * restriction) / (routes - i)
+        i += 1
+        flow, F = maximum_flow(flow_graph, source, target, capacity_matrix,
+        algorithm = flow_algorithm, restriction = restriction)
+    end
+
+    # End
+    return flow, F
+end
diff --git a/src/maximum_flow.jl b/src/maximum_flow.jl
new file mode 100644
index 0000000..56a2fab
--- /dev/null
+++ b/src/maximum_flow.jl
@@ -0,0 +1,180 @@
+"""
+    AbstractFlowAlgorithm
+
+Abstract type that allows users to pass in their preferred algorithm
+"""
+abstract type AbstractFlowAlgorithm end
+
+"""
+    EdmondsKarpAlgorithm <: AbstractFlowAlgorithm
+
+Forces the maximum_flow function to use the Edmonds–Karp algorithm.
+"""
+struct EdmondsKarpAlgorithm <: AbstractFlowAlgorithm end
+
+"""
+    DinicAlgorithm <: AbstractFlowAlgorithm
+
+Forces the maximum_flow function to use Dinic's algorithm.
+"""
+struct DinicAlgorithm <: AbstractFlowAlgorithm end
+
+"""
+    BoykovKolmogorovAlgorithm <: AbstractFlowAlgorithm
+
+Forces the maximum_flow function to use the Boykov-Kolmogorov algorithm.
+"""
+struct BoykovKolmogorovAlgorithm <: AbstractFlowAlgorithm end
+
+"""
+Forces the maximum_flow function to use the Push-Relabel algorithm.
+"""
+struct PushRelabelAlgorithm <: AbstractFlowAlgorithm end
+
+"""
+    DefaultCapacity{T}
+
+Structure that returns `1` if a forward edge exists in `flow_graph`, and `0` otherwise.
+"""
+struct DefaultCapacity{T<:Integer} <: AbstractMatrix{T}
+    flow_graph::lg.DiGraph
+    nv::T
+end
+
+@traitfn DefaultCapacity(flow_graph::::lg.IsDirected) =
+    DefaultCapacity(lg.DiGraph(flow_graph), lg.nv(flow_graph))
+
+getindex(d::DefaultCapacity{T}, s::Integer, t::Integer) where T = if lg.has_edge(d.flow_graph, s, t) one(T) else zero(T) end
+# isassigned{T<:Integer}(d::DefaultCapacity{T}, u::T, v::T) = (u in 1:d.nv) && (v in 1:d.nv)
+size(d::DefaultCapacity) = (Int(d.nv), Int(d.nv))
+transpose(d::DefaultCapacity) = DefaultCapacity(reverse(d.flow_graph))
+ctranspose(d::DefaultCapacity) = DefaultCapacity(reverse(d.flow_graph))
+
+"""
+    residual(flow_graph)
+
+Return a directed residual graph for a directed `flow_graph`.
+
+The residual graph comprises the same node list as the orginal flow graph, but
+ensures that for each edge (u,v), (v,u) also exists in the graph. This allows
+flow in the reverse direction.
+
+If only the forward edge exists, a reverse edge is created with capacity 0.
+If both forward and reverse edges exist, their capacities are left unchanged.
+Since the capacities in [`LightGraphs.DefaultDistance`](@ref) cannot be changed, an array of ones
+is created.
+"""
+function residual end
+@traitfn residual(flow_graph::::lg.IsDirected) = lg.DiGraph(lg.Graph(flow_graph))
+
+# Method for Edmonds–Karp algorithm
+
+@traitfn function maximum_flow(
+    flow_graph::::lg.IsDirected,                   # the input graph
+    source::Integer,                       # the source vertex
+    target::Integer,                       # the target vertex
+    capacity_matrix::AbstractMatrix,   # edge flow capacities
+    algorithm::EdmondsKarpAlgorithm        # keyword argument for algorithm
+    )
+    residual_graph = residual(flow_graph)
+    return edmonds_karp_impl(residual_graph, source, target, capacity_matrix)
+end
+
+# Method for Dinic's algorithm
+
+@traitfn function maximum_flow(
+    flow_graph::::lg.IsDirected,                   # the input graph
+    source::Integer,                       # the source vertex
+    target::Integer,                       # the target vertex
+    capacity_matrix::AbstractMatrix,   # edge flow capacities
+    algorithm::DinicAlgorithm              # keyword argument for algorithm
+    )
+    residual_graph = residual(flow_graph)
+    return dinic_impl(residual_graph, source, target, capacity_matrix)
+end
+
+# Method for Boykov-Kolmogorov algorithm
+
+@traitfn function maximum_flow(
+    flow_graph::::lg.IsDirected,                   # the input graph
+    source::Integer,                       # the source vertex
+    target::Integer,                       # the target vertex
+    capacity_matrix::AbstractMatrix,   # edge flow capacities
+    algorithm::BoykovKolmogorovAlgorithm   # keyword argument for algorithm
+    )
+    residual_graph = residual(flow_graph)
+    return boykov_kolmogorov_impl(residual_graph, source, target, capacity_matrix)
+end
+
+# Method for Push-relabel algorithm
+
+@traitfn function maximum_flow(
+    flow_graph::::lg.IsDirected,                   # the input graph
+    source::Integer,                       # the source vertex
+    target::Integer,                       # the target vertex
+    capacity_matrix::AbstractMatrix,   # edge flow capacities
+    algorithm::PushRelabelAlgorithm        # keyword argument for algorithm
+    )
+    residual_graph = residual(flow_graph)
+    return push_relabel(residual_graph, source, target, capacity_matrix)
+end
+
+"""
+    maximum_flow(flow_graph, source, target[, capacity_matrix][, algorithm][, restriction])
+
+Generic maximum_flow function for `flow_graph` from `source` to `target` with
+capacities in `capacity_matrix`.
+Uses flow algorithm `algorithm` and cutoff restriction `restriction`.
+
+- If `capacity_matrix` is not specified, `DefaultCapacity(flow_graph)` will be used.
+- If `algorithm` is not specified, it will default to [`PushRelabelAlgorithm`](@ref).
+- If `restriction` is not specified, it will default to `0`.
+
+Return a tuple of (maximum flow, flow matrix). For the Boykov-Kolmogorov
+algorithm, the associated mincut is returned as a third output.
+
+### Usage Example:
+
+```jldoctest
+julia> flow_graph = lg.DiGraph(8) # Create a flow-graph
+julia> flow_edges = [
+(1,2,10),(1,3,5),(1,4,15),(2,3,4),(2,5,9),
+(2,6,15),(3,4,4),(3,6,8),(4,7,16),(5,6,15),
+(5,8,10),(6,7,15),(6,8,10),(7,3,6),(7,8,10)
+]
+
+julia> capacity_matrix = zeros(Int, 8, 8)  # Create a capacity matrix
+
+julia> for e in flow_edges
+    u, v, f = e
+    lg.add_edge!(flow_graph, u, v)
+    capacity_matrix[u,v] = f
+end
+
+julia> f, F = maximum_flow(flow_graph, 1, 8) # Run default maximum_flow (push-relabel) without the capacity_matrix
+
+julia> f, F = maximum_flow(flow_graph, 1, 8, capacity_matrix) # Run default maximum_flow with the capacity_matrix
+
+julia> f, F = maximum_flow(flow_graph, 1, 8, capacity_matrix, algorithm=EdmondsKarpAlgorithm()) # Run Edmonds-Karp algorithm
+
+julia> f, F = maximum_flow(flow_graph, 1, 8, capacity_matrix, algorithm=DinicAlgorithm()) # Run Dinic's algorithm
+
+julia> f, F, labels = maximum_flow(flow_graph, 1, 8, capacity_matrix, algorithm=BoykovKolmogorovAlgorithm()) # Run Boykov-Kolmogorov algorithm
+
+```
+"""
+function maximum_flow(
+    flow_graph::lg.AbstractGraph,                   # the input graph
+    source::Integer,                       # the source vertex
+    target::Integer,                       # the target vertex
+    capacity_matrix::AbstractMatrix =  # edge flow capacities
+    DefaultCapacity(flow_graph);
+    algorithm::AbstractFlowAlgorithm  =    # keyword argument for algorithm
+    PushRelabelAlgorithm(),
+    restriction::Real = 0               # keyword argument for restriction max-flow
+    )
+    if restriction > 0
+        return maximum_flow(flow_graph, source, target, min.(restriction, capacity_matrix), algorithm)
+    end
+    return maximum_flow(flow_graph, source, target, capacity_matrix, algorithm)
+end
diff --git a/src/multiroute_flow.jl b/src/multiroute_flow.jl
new file mode 100644
index 0000000..b60cd98
--- /dev/null
+++ b/src/multiroute_flow.jl
@@ -0,0 +1,226 @@
+"""
+    AbstractMultirouteFlowAlgorithm
+
+Abstract type that allows users to pass in their preferred algorithm.
+"""
+abstract type AbstractMultirouteFlowAlgorithm end
+
+"""
+    KishimotoAlgorithm
+
+Used to specify the Kishimoto algorithm.
+"""
+struct KishimotoAlgorithm <: AbstractMultirouteFlowAlgorithm end
+
+"""
+    ExtendedMultirouteFlowAlgorithm
+
+Used to specify the Extended Multiroute Flow algorithm.
+"""
+struct ExtendedMultirouteFlowAlgorithm <: AbstractMultirouteFlowAlgorithm end
+
+# Methods when the number of routes is more than the connectivity
+# 1) When using Boykov-Kolmogorov as a flow subroutine
+# 2) Other flow algorithm
+function empty_flow(
+    capacity_matrix::AbstractMatrix{T},     # edge flow capacities
+    flow_algorithm::BoykovKolmogorovAlgorithm # keyword argument for algorithm
+    ) where T<:Real
+    n = size(capacity_matrix, 1)
+    return zero(T), zeros(T, n, n), zeros(T, n)
+end
+# 2) Other flow algorithm
+function empty_flow(
+    capacity_matrix::AbstractMatrix{T},     # edge flow capacities
+    flow_algorithm::AbstractFlowAlgorithm     # keyword argument for algorithm
+    ) where T<:Real
+    n = size(capacity_matrix, 1)
+    return zero(T), zeros(T, n, n)
+end
+
+# Method for Kishimoto algorithm
+@traitfn function multiroute_flow(
+    flow_graph::::lg.IsDirected,                   # the input graph
+    source::Integer,                       # the source vertex
+    target::Integer,                       # the target vertex
+    capacity_matrix::AbstractMatrix,  # edge flow capacities
+    flow_algorithm::AbstractFlowAlgorithm, # keyword argument for algorithm
+    mrf_algorithm::KishimotoAlgorithm,     # keyword argument for algorithm
+    routes::Int                            # keyword argument for routes
+    )
+    return kishimoto(flow_graph, source, target, capacity_matrix, flow_algorithm, routes)
+end
+
+## Methods for Extended Multiroute Flow Algorithm
+#1 When the breaking points are not already known
+@traitfn function multiroute_flow(
+    flow_graph::::lg.IsDirected,                       # the input graph
+    source::Integer,                                # the source vertex
+    target::Integer,                                # the target vertex
+    capacity_matrix::AbstractMatrix,           # edge flow capacities
+    flow_algorithm::AbstractFlowAlgorithm,          # keyword argument for algorithm
+    mrf_algorithm::ExtendedMultirouteFlowAlgorithm, # keyword argument for algorithm
+    routes::Real                                       # keyword argument for routes
+    )
+    return emrf(flow_graph, source, target, capacity_matrix, flow_algorithm, routes)
+end
+#2 When the breaking points are already known
+#2-a Output: flow value (paired with the associated restriction)
+multiroute_flow(
+    breakingpoints::Vector{Tuple{T,T,Int}},       # vector of breaking points
+    routes::R                                       # keyword argument for routes
+    ) where T<:Real where R<:Real =
+    intersection(breakingpoints, routes)
+
+#2-b Output: flow value, flows(, labels)
+function multiroute_flow(
+    breakingpoints::AbstractVector{Tuple{T1,T1,Int}}, # vector of breaking points
+    routes::R,                                # keyword argument for routes
+    flow_graph::lg.AbstractGraph,                      # the input graph
+    source::Integer,                          # the source vertex
+    target::Integer,                          # the target vertex
+    capacity_matrix::AbstractMatrix{T2} =   # edge flow capacities
+    DefaultCapacity(flow_graph);
+    flow_algorithm::AbstractFlowAlgorithm  =  # keyword argument for algorithm
+    PushRelabelAlgorithm()
+    ) where T2 where T1<:Real where R<:Real
+    x, f = intersection(breakingpoints, routes)
+    # For other cases, capacities need to be Floats
+    if !(T2<:AbstractFloat)
+        capacity_matrix = convert(AbstractMatrix{Float64}, capacity_matrix)
+    end
+
+    return maximum_flow(flow_graph, source, target, capacity_matrix,
+    algorithm = flow_algorithm, restriction = x)
+end
+
+### TODO: CLEAN UP THIS FUNCTION AND DOCUMENTATION. THERE SHOULD BE NO NEED TO
+### HAVE A TYPE-UNSTABLE FUNCTION HERE. (sbromberger 2017-03-26)
+"""
+    multiroute_flow(flow_graph, source, target[, DefaultCapacity][, flow_algorithm][, mrf_algorithm][, routes])
+
+The generic multiroute_flow function.
+
+The output will vary depending on the input:
+
+- When the number of `route`s is `0`, return the set of breaking points of
+the multiroute flow.
+- When the number of `route`s is `1`, return a flow with a set of 1-disjoint paths
+(this is the classical max-flow implementation).
+- When the input is limited to a set of breaking points and a route value `k`,
+return only the k-route flow.
+- Otherwise, a tuple with 1) the maximum flow and 2) the flow matrix. When the
+max-flow subroutine is the Boykov-Kolmogorov algorithm, the associated mincut is
+returned as a third output.
+
+When the input is a network, it requires the following arguments:
+
+- `flow_graph`: the input graph
+- `source`: the source vertex
+- `target`: the target vertex
+- `capacity_matrix`: matrix of edge flow capacities
+- `flow_algorithm`: keyword argument for flow algorithm
+- `mrf_algorithm`: keyword argument for multiroute flow algorithm
+- `routes`: keyword argument for the number of routes
+
+When the input is only the set of (breaking) points and the number of route,
+it requires the following arguments:
+
+- `breakingpoints`: vector of breaking points
+- `routes`: number of routes
+
+When the input is the set of (breaking) points, the number of routes,
+and the network descriptors, it requires the following arguments:
+
+- `breakingpoints`: vector of breaking points
+- `routes`: number of routes
+- `flow_graph`: the input graph
+- `source`: the source vertex
+- `target`: the target vertex
+- `capacity_matrix`: matrix of edge flow capacities
+- `flow_algorithm`: keyword argument for flow algorithm
+
+The function defaults to the Push-relabel (classical flow) and Kishimoto
+(multiroute) algorithms. Alternatively, the algorithms to be used can also
+be specified through  keyword arguments. A default capacity of `1` is assumed
+for each link if no capacity matrix is provided.
+
+The `mrf_algorithm` keyword is inforced to Extended Multiroute Flow
+in the following cases:
+
+- The number of routes is non-integer
+- The number of routes is 0 or non-specified
+
+### Usage Example :
+(please consult the  [`maximum_flow`](@ref) section for options about flow_algorithm
+and capacity_matrix)
+
+```jldoctest
+julia> flow_graph = lg.DiGraph(8) # Create a flow graph
+
+julia> flow_edges = [
+(1, 2, 10), (1, 3, 5),  (1, 4, 15), (2, 3, 4),  (2, 5, 9),
+(2, 6, 15), (3, 4, 4),  (3, 6, 8),  (4, 7, 16), (5, 6, 15),
+(5, 8, 10), (6, 7, 15), (6, 8, 10), (7, 3, 6),  (7, 8, 10)
+]
+
+julia> capacity_matrix = zeros(Int, 8, 8) # Create a capacity matrix
+
+julia> for e in flow_edges
+    u, v, f = e
+    add_edge!(flow_graph, u, v)
+    capacity_matrix[u, v] = f
+end
+
+julia> f, F = multiroute_flow(flow_graph, 1, 8, capacity_matrix, routes = 2) # Run default multiroute_flow with an integer number of routes = 2
+
+julia> f, F = multiroute_flow(flow_graph, 1, 8, capacity_matrix, routes = 1.5) # Run default multiroute_flow with a noninteger number of routes = 1.5
+
+julia> points = multiroute_flow(flow_graph, 1, 8, capacity_matrix) # Run default multiroute_flow for all the breaking points values
+
+julia> f, F = multiroute_flow(points, 1.5) # Then run multiroute flow algorithm for any positive number of routes
+
+julia> f = multiroute_flow(points, 1.5, valueonly = true)
+
+julia> f, F, labels = multiroute_flow(flow_graph, 1, 8, capacity_matrix, algorithm = BoykovKolmogorovAlgorithm(), routes = 2) # Run multiroute flow algorithm using Boykov-Kolmogorov algorithm as maximum_flow routine
+
+```
+"""
+function multiroute_flow(
+    flow_graph::lg.AbstractGraph,                    # the input graph
+    source::Integer,                        # the source vertex
+    target::Integer,                        # the target vertex
+    capacity_matrix::AbstractMatrix{T} =  # edge flow capacities
+        DefaultCapacity(flow_graph);
+    flow_algorithm::AbstractFlowAlgorithm  =    # keyword argument for algorithm
+    PushRelabelAlgorithm(),
+    mrf_algorithm::AbstractMultirouteFlowAlgorithm  =    # keyword argument for algorithm
+    KishimotoAlgorithm(),
+    routes::R = 0              # keyword argument for number of routes (0 = all values)
+    ) where T where R <: Real
+
+    # a flow with a set of 1-disjoint paths is a classical max-flow
+    (routes == 1) &&
+    return maximum_flow(flow_graph, source, target, capacity_matrix, flow_algorithm)
+
+    # routes > λ (connectivity) → f = 0
+    λ = maximum_flow(flow_graph, source, target, DefaultCapacity(flow_graph),
+    algorithm = flow_algorithm)[1]
+    (routes > λ) && return empty_flow(capacity_matrix, flow_algorithm)
+
+    # For other cases, capacities need to be Floats
+    if !(T<:AbstractFloat)
+        capacity_matrix = convert(AbstractMatrix{Float64}, capacity_matrix)
+    end
+
+    # Ask for all possible values (breaking points)
+    (routes == 0) &&
+    return emrf(flow_graph, source, target, capacity_matrix, flow_algorithm)
+    # The number of routes is a float → EMRF
+    (R <: AbstractFloat) &&
+    return emrf(flow_graph, source, target, capacity_matrix, flow_algorithm, routes)
+
+    # Other calls
+    return multiroute_flow(flow_graph, source, target, capacity_matrix,
+    flow_algorithm, mrf_algorithm, routes)
+end
diff --git a/src/push_relabel.jl b/src/push_relabel.jl
new file mode 100644
index 0000000..db80bbc
--- /dev/null
+++ b/src/push_relabel.jl
@@ -0,0 +1,205 @@
+@doc_str """
+    push_relabel(residual_graph, source, target, capacity_matrix)
+
+Return the maximum flow of `residual_graph` from `source` to `target` using the
+FIFO push relabel algorithm with gap heuristic.
+
+### Performance
+Takes approximately ``\\mathcal{O}(|V|^{3})`` time.
+"""
+function push_relabel end
+@traitfn function push_relabel{T}(
+    residual_graph::::lg.IsDirected,       # the input graph
+    source::Integer,                       # the source vertex
+    target::Integer,                       # the target vertex
+    capacity_matrix::AbstractMatrix{T}     # edge flow capacities
+    )
+
+    n = lg.nv(residual_graph)
+    flow_matrix = zeros(T, n, n)
+
+    height = zeros(Int, n)
+    height[source] = n
+
+    count = zeros(Int, 2 * n + 1)
+    count[0 + 1] = n - 1
+    count[n + 1] = 1
+
+    excess = zeros(T, n)
+    excess[source] = typemax(T)
+
+    active = falses(n)
+    active[source] = true
+    active[target] = true
+
+    Q = Array{Int,1}()
+    sizehint!(Q, n)
+
+
+    for v in lg.out_neighbors(residual_graph, source)
+        push_flow!(residual_graph, source, v, capacity_matrix, flow_matrix, excess, height, active, Q)
+    end
+
+    while length(Q) > 0
+        v = pop!(Q)
+        active[v] = false
+        discharge!(residual_graph, v, capacity_matrix, flow_matrix, excess, height, active, count, Q)
+    end
+
+    return sum([flow_matrix[v, target] for v in lg.in_neighbors(residual_graph, target)]), flow_matrix
+end
+
+"""
+    enqueue_vertex!(Q, v, active, excess)
+
+Push inactive node `v` into queue `Q` and activates it. Requires preallocated
+`active` and `excess` vectors.
+"""
+
+function enqueue_vertex!(
+    Q::AbstractVector,
+    v::Integer,                                # input vertex
+    active::AbstractVector{Bool},
+    excess::AbstractVector
+    )
+    if !active[v] && excess[v] > 0
+        active[v] = true
+        unshift!(Q, v)
+    end
+    return nothing
+end
+
+"""
+    push_flow!(residual_graph, u, v, capacity_matrix, flow_matrix, excess, height, active, Q)
+
+Using `residual_graph` with capacities in `capacity_matrix`, push as much flow
+as possible through the given edge(`u`, `v`). Requires preallocated `flow_matrix`
+matrix, and `excess`, `height, `active`, and `Q` vectors.
+"""
+function push_flow! end
+@traitfn function push_flow!(
+    residual_graph::::lg.IsDirected,             # the input graph
+    u::Integer,                              # input from-vertex
+    v::Integer,                              # input to-vetex
+    capacity_matrix::AbstractMatrix,
+    flow_matrix::AbstractMatrix,
+    excess::AbstractVector,
+    height::AbstractVector{Int},
+    active::AbstractVector{Bool},
+    Q::AbstractVector
+    )
+    flow = min(excess[u], capacity_matrix[u, v] - flow_matrix[u, v])
+
+    flow == 0 && return nothing
+    height[u] <= height[v] && return nothing
+
+    flow_matrix[u, v] += flow
+    flow_matrix[v, u] -= flow
+
+    excess[u] -= flow
+    excess[v] += flow
+
+    enqueue_vertex!(Q, v, active, excess)
+    nothing
+end
+
+"""
+    gap!(residual_graph, h, excess, height, active, count, Q)
+
+Implement the push-relabel gap heuristic. Relabel all vertices above a cutoff height.
+Reduce the number of relabels required.
+
+Requires arguments:
+
+- residual_graph::DiGraph                # the input graph
+- h::Int                                 # cutoff height
+- excess::AbstractVector
+- height::AbstractVector{Int}
+- active::AbstractVector{Bool}
+- count::AbstractVector{Int}
+- Q::AbstractVector
+"""
+function gap! end
+@traitfn function gap!(
+    residual_graph::::lg.IsDirected,               # the input graph
+    h::Int,                                # cutoff height
+    excess::AbstractVector,
+    height::AbstractVector{Int},
+    active::AbstractVector{Bool},
+    count::AbstractVector{Int},
+    Q::AbstractVector                # FIFO queue
+    )
+    n = lg.nv(residual_graph)
+    for v in lg.vertices(residual_graph)
+        height[v] < h && continue
+        count[height[v] + 1] -= 1
+        height[v] = max(height[v], n + 1)
+        count[height[v] + 1] += 1
+        enqueue_vertex!(Q, v, active, excess)
+    end
+    nothing
+end
+
+"""
+    relabel!(residual_graph, v, capacity_matrix, flow_matrix, excess, height, active, count, Q)
+
+Relabel a node `v` with respect to its neighbors to produce an admissable edge.
+"""
+function relabel! end
+@traitfn function relabel!(
+    residual_graph::::lg.IsDirected,                # the input graph
+    v::Integer,                                 # input vertex to be relabeled
+    capacity_matrix::AbstractMatrix,
+    flow_matrix::AbstractMatrix,
+    excess::AbstractVector,
+    height::AbstractVector{Int},
+    active::AbstractVector{Bool},
+    count::AbstractVector{Int},
+    Q::AbstractVector
+    )
+    n = lg.nv(residual_graph)
+    count[height[v] + 1] -= 1
+    height[v] = 2 * n
+    for to in lg.out_neighbors(residual_graph, v)
+        if capacity_matrix[v, to] > flow_matrix[v, to]
+            height[v] = min(height[v], height[to] + 1)
+        end
+    end
+    count[height[v] + 1] += 1
+    enqueue_vertex!(Q, v, active, excess)
+    nothing
+end
+
+
+"""
+    discharge!(residual_graph, v, capacity_matrix, flow_matrix, excess, height, active, count, Q)
+
+Drain the excess flow out of node `v`. Run the gap heuristic or relabel the
+vertex if the excess remains non-zero.
+"""
+function discharge! end
+@traitfn function discharge!(
+    residual_graph::::lg.IsDirected,                # the input graph
+    v::Integer,                                 # vertex to be discharged
+    capacity_matrix::AbstractMatrix,
+    flow_matrix::AbstractMatrix,
+    excess::AbstractVector,
+    height::AbstractVector{Int},
+    active::AbstractVector{Bool},
+    count::AbstractVector{Int},
+    Q::AbstractVector                 # FIFO queue
+    )
+    for to in lg.out_neighbors(residual_graph, v)
+        excess[v] == 0 && break
+        push_flow!(residual_graph, v, to, capacity_matrix, flow_matrix, excess, height, active, Q)
+    end
+
+    if excess[v] > 0
+        if count[height[v] + 1] == 1
+            gap!(residual_graph, height[v], excess, height, active, count, Q)
+        else
+            relabel!(residual_graph, v, capacity_matrix, flow_matrix, excess, height, active, count, Q)
+        end
+    end
+    nothing
+end
diff --git a/test/boykov_kolmogorov.jl b/test/boykov_kolmogorov.jl
new file mode 100644
index 0000000..45813b2
--- /dev/null
+++ b/test/boykov_kolmogorov.jl
@@ -0,0 +1,30 @@
+@testset "Boykov Kolmogorov" begin
+    # construct graph
+    gg = lg.DiGraph(3)
+    lg.add_edge!(gg, 1, 2)
+    lg.add_edge!(gg, 2, 3)
+
+    # source and sink terminals
+    source, target = 1, 3
+
+
+    for g in testdigraphs(gg)
+    # default capacity
+      capacity_matrix = lg.DefaultCapacity(g)
+      residual_graph = @inferred(lg.residual(g))
+      T = eltype(g)
+      flow_matrix = zeros(T, 3, 3)
+      TREE = zeros(T, 3)
+      TREE[source] = T(1)
+      TREE[target] = T(2)
+      PARENT = zeros(T, 3)
+      A = [T(source), T(target)]
+# see https://github.com/JuliaLang/julia/issues/21077
+# @show("testing $g with eltype $T, residual_graph type is $(eltype(residual_graph)), flow_matrix type is $(eltype(flow_matrix)), capacity_matrix type is $(eltype(capacity_matrix))")
+      path = lg.find_path!(
+        residual_graph, source, target, flow_matrix,
+        capacity_matrix, PARENT, TREE, A)
+
+      @test path == [1, 2, 3]
+  end
+end
diff --git a/test/dinic.jl b/test/dinic.jl
new file mode 100644
index 0000000..cb6a04b
--- /dev/null
+++ b/test/dinic.jl
@@ -0,0 +1,62 @@
+@testset "Dinic" begin
+    # Construct DiGraph
+    flow_graph = lg.DiGraph(8)
+
+    # Load custom dataset
+    flow_edges = [
+        (1, 2, 10), (1, 3, 5), (1, 4, 15), (2, 3, 4), (2, 5, 9),
+        (2, 6, 15), (3, 4, 4), (3, 6, 8), (4, 7, 16), (5, 6, 15),
+        (5, 8, 10), (6, 7, 15), (6, 8, 10), (7, 3, 6), (7, 8, 10)
+    ]
+
+    capacity_matrix = zeros(Int, lg.nv(flow_graph), lg.nv(flow_graph))
+
+    for e in flow_edges
+        u, v, f = e
+        lg.add_edge!(flow_graph, u, v)
+        capacity_matrix[u, v] = f
+    end
+
+    # Construct the residual graph
+    for fg in (flow_graph, lg.DiGraph{UInt8}(flow_graph), lg.DiGraph{Int16}(flow_graph))
+      residual_graph = @inferred(lg.residual(fg))
+
+      # Test with default distances
+      @test @inferred(LightGraphsFlows.dinic_impl(residual_graph, 1, 8, LightGraphsFlows.DefaultCapacity(residual_graph)))[1] == 3
+
+      # Test with capacity matrix
+      @test @inferred(LightGraphsFlows.dinic_impl(residual_graph, 1, 8, capacity_matrix))[1] == 28
+
+      # Test on disconnected graphs
+      function test_blocking_flow(residual_graph, source, target, capacity_matrix, flow_matrix)
+          #disconnect source
+          h = copy(residual_graph)
+          for dst in collect(lg.neighbors(residual_graph, source))
+              lg.rem_edge!(h, source, dst)
+          end
+          @test @inferred(LightGraphsFlows.blocking_flow(h, source, target, capacity_matrix, flow_matrix)) == 0
+
+          #disconnect target and add unreachable vertex
+          h = copy(residual_graph)
+          for src in collect(lg.in_neighbors(residual_graph, target))
+              lg.rem_edge!(h, src, target)
+          end
+          @test @inferred(LightGraphsFlows.blocking_flow(h, source, target, capacity_matrix, flow_matrix)) == 0
+
+          # unreachable vertex (covers the case where a vertex isn't reachable from the source)
+          h = copy(residual_graph)
+          lg.add_vertex!(h)
+          lg.add_edge!(h, lg.nv(residual_graph) + 1, target)
+          capacity_matrix_ = vcat(hcat(capacity_matrix, zeros(Int, lg.nv(residual_graph))), zeros(Int, 1, lg.nv(residual_graph) + 1))
+          flow_graph_  = vcat(hcat(flow_matrix, zeros(Int, lg.nv(residual_graph))), zeros(Int, 1, lg.nv(residual_graph) + 1))
+
+          @test @inferred(LightGraphsFlows.blocking_flow(h, source, target, capacity_matrix_, flow_graph_)) > 0
+
+          #test with connected graph
+          @test @inferred(LightGraphsFlows.blocking_flow(residual_graph, source, target, capacity_matrix, flow_matrix)) > 0
+      end
+
+      flow_matrix = zeros(Int, lg.nv(residual_graph), lg.nv(residual_graph))
+      test_blocking_flow(residual_graph, 1, 8, capacity_matrix, flow_matrix)
+    end
+end
diff --git a/test/edmonds_karp.jl b/test/edmonds_karp.jl
new file mode 100644
index 0000000..bca4148
--- /dev/null
+++ b/test/edmonds_karp.jl
@@ -0,0 +1,59 @@
+@testset "Edmonds Karp" begin
+    # Construct DiGraph
+    flow_graph = lg.DiGraph(8)
+
+    # Load custom dataset
+    flow_edges = [
+        (1, 2, 10), (1, 3, 5), (1, 4, 15), (2, 3, 4), (2, 5, 9),
+        (2, 6, 15), (3, 4, 4), (3, 6, 8), (4, 7, 16), (5, 6, 15),
+        (5, 8, 10), (6, 7, 15), (6, 8, 10), (7, 3, 6), (7, 8, 10)
+    ]
+
+    for fg in (flow_graph, lg.DiGraph{UInt8}(flow_graph), lg.DiGraph{Int16}(flow_graph))
+        capacity_matrix = zeros(Int, 8, 8)
+        for e in flow_edges
+            u, v, f = e
+            lg.add_edge!(fg, u, v)
+            capacity_matrix[u, v] = f
+        end
+        residual_graph = @inferred(LightGraphsFlows.residual(fg))
+
+        # Test with default distances
+        @test @inferred(LightGraphsFlows.edmonds_karp_impl(residual_graph, 1, 8, LightGraphsFlows.DefaultCapacity(residual_graph)))[1] == 3
+
+        # Test with capacity matrix
+        @test @inferred(LightGraphsFlows.edmonds_karp_impl(residual_graph, 1, 8, capacity_matrix))[1] == 28
+
+        # Test the types of the values returned by fetch_path
+        function test_find_path_types(residual_graph, s, t, flow_matrix, capacity_matrix)
+            v, P, S, flag = LightGraphsFlows.fetch_path(residual_graph, s, t, flow_matrix, capacity_matrix)
+            @test typeof(P) == Vector{Int}
+            @test typeof(S) == Vector{Int}
+            @test typeof(flag) == Int
+            @test typeof(v) == eltype(residual_graph)
+        end
+
+        # Test the value of the flags returned.
+        function test_find_path_disconnected(residual_graph, s, t, flow_matrix, capacity_matrix)
+            h = copy(residual_graph)
+            for dst in collect(lg.neighbors(residual_graph, s))
+                lg.rem_edge!(residual_graph, s, dst)
+            end
+            v, P, S, flag = LightGraphsFlows.fetch_path(residual_graph, s, t, flow_matrix, capacity_matrix)
+            @test flag == 1
+            for dst in collect(lg.neighbors(h, t))
+                lg.rem_edge!(h, t, dst)
+            end
+            v, P, S, flag = LightGraphsFlows.fetch_path(h, s, t, flow_matrix, capacity_matrix)
+            @test flag == 0
+            for i in collect(lg.in_neighbors(h, t))
+                lg.rem_edge!(h, i, t)
+            end
+            v, P, S, flag = LightGraphsFlows.fetch_path(h, s, t, flow_matrix, capacity_matrix)
+            @test flag == 2
+        end
+        flow_matrix = zeros(Int, lg.nv(residual_graph), lg.nv(residual_graph))
+        test_find_path_types(residual_graph, 1, 8, flow_matrix, capacity_matrix)
+        test_find_path_disconnected(residual_graph, 1, 8, flow_matrix, capacity_matrix)
+    end
+end
diff --git a/test/maximum_flow.jl b/test/maximum_flow.jl
new file mode 100644
index 0000000..ced907f
--- /dev/null
+++ b/test/maximum_flow.jl
@@ -0,0 +1,57 @@
+@testset "Maximum flow" begin
+    #### Graphs for testing
+    graphs = [
+    # Graph with 8 vertices
+    (8,
+    [
+    (1, 2, 10), (1, 3, 5), (1, 4, 15), (2, 3, 4), (2, 5, 9),
+    (2, 6, 15), (3, 4, 4), (3, 6, 8), (4, 7, 16), (5, 6, 15),
+    (5, 8, 10), (6, 7, 15), (6, 8, 10), (7, 3, 6), (7, 8, 10)
+    ],
+    1, 8, # source/target
+    3,    # answer for default capacity
+    28,   # answer for custom capacity
+    15, 5 # answer for restricted capacity/restriction
+    ),
+
+    # Graph with 6 vertices
+    (6,
+    [
+    (1, 2, 9), (1, 3, 9), (2, 3, 10), (2, 4, 8), (3, 4, 1),
+    (3, 5, 3), (5, 4, 8), (4, 6, 10), (5, 6, 7)
+    ],
+    1, 6, # source/target
+    2,   # answer for default capacity
+    12,  # answer for custom capacity
+    8, 5  # answer for restricted capacity/restriction
+    )
+    ]
+
+    for (nvertices, flow_edges, s, t, fdefault, fcustom, frestrict, caprestrict) in graphs
+        flow_graph = lg.DiGraph(nvertices)
+        for g in testdigraphs(flow_graph)
+            capacity_matrix = zeros(Int, nvertices, nvertices)
+            for e in flow_edges
+                u, v, f = e
+                lg.add_edge!(g, u, v)
+                capacity_matrix[u, v] = f
+            end
+
+            # Test DefaultCapacity
+            d = @inferred(LightGraphsFlows.DefaultCapacity(g))
+            T = eltype(d)
+            @test typeof(d) <: AbstractMatrix{T}
+            @test d[s, t] == 0
+            @test size(d) == (nvertices, nvertices)
+            @test typeof(transpose(d)) <: LightGraphsFlows.DefaultCapacity
+            @test typeof(ctranspose(d)) <: LightGraphsFlows.DefaultCapacity
+
+            # Test all algorithms - type instability in PushRelabel #553
+            for ALGO in [EdmondsKarpAlgorithm, DinicAlgorithm, BoykovKolmogorovAlgorithm, PushRelabelAlgorithm]
+                @test maximum_flow(g, s, t, algorithm=ALGO())[1] == fdefault
+                @test maximum_flow(g, s, t, capacity_matrix, algorithm=ALGO())[1] == fcustom
+                @test maximum_flow(g, s, t, capacity_matrix, algorithm=ALGO(), restriction=caprestrict)[1] == frestrict
+            end
+        end
+    end
+end
diff --git a/test/multiroute_flow.jl b/test/multiroute_flow.jl
new file mode 100644
index 0000000..b4a7e0b
--- /dev/null
+++ b/test/multiroute_flow.jl
@@ -0,0 +1,86 @@
+@testset "Multiroute flow" begin
+    #### Graphs for testing
+    graphs = [
+    # Graph with 8 vertices
+    (8,
+    [
+    (1, 2, 10), (1, 3, 5),  (1, 4, 15), (2, 3, 4),  (2, 5, 9),
+    (2, 6, 15), (3, 4, 4),  (3, 6, 8),  (4, 7, 16), (5, 6, 15),
+    (5, 8, 10), (6, 7, 15), (6, 8, 10), (7, 3, 6),  (7, 8, 10),
+    (8, 1, 8)                        # Reverse edge to test the slope in EMRF
+    ],
+    1, 8,                              # source/target
+    [28, 28, 15, 0],                   # answer for 1 to 4 route(s) flows
+    [(0., 0., 3), (5., 15., 2),        # breaking points
+    (10., 25., 1), (13., 28., 0)],
+    28.                                # value for 1.5 routes
+    ),
+
+    # Graph with 6 vertices
+    (6,
+    [
+    (1, 2, 9), (1, 3, 9), (2, 3, 10), (2, 4, 8), (3, 4, 1),
+    (3, 5, 3), (5, 4, 8), (4, 6, 10), (5, 6, 7)
+    ],
+    1, 6,                                      # source/target
+    [12, 6, 0],                                # answer for 1 to 3 route(s) flows
+    [(0., 0., 2), (3., 6., 1), (9., 12., 0)],  # breaking points
+    9.                                         # value for 1.5 routes
+    ),
+
+    # Graph with 7 vertices
+    (7,
+    [
+    (1, 2, 1), (1, 3, 2), (1, 4, 3), (1, 5, 4), (1, 6, 5),
+    (2, 7, 1), (3, 7, 2), (4, 7, 3), (5, 7, 4), (6, 7, 5)
+    ],
+    1, 7,                                      # source/target
+    [15, 15, 15, 12, 5, 0],                    # answer for 1 to 6 route(s) flows
+    [(0., 0., 5), (1., 5., 4), (2., 9., 3),    # breaking points
+    (3., 12., 2), (4., 14., 1), (5., 15., 0)],
+    15.                                        # value for 1.5 routes
+    ),
+
+    # Graph with 6 vertices
+    (6,
+    [
+    (1, 2, 1), (1, 3, 1), (1, 4, 2), (1, 5, 2),
+    (2, 6, 1), (3, 6, 1), (4, 6, 2), (5, 6, 2),
+    ],
+    1, 6,                                      # source/target
+    [6, 6, 6, 4, 0],                           # answer for 1 to 5 route(s) flows
+    [(0., 0., 4), (1., 4., 2), (2., 6., 0)],   # breaking points
+    6.                                         # value for 1.5 routes
+    )
+    ]
+
+    for (nvertices, flow_edges, s, t, froutes, breakpts, ffloat) in graphs
+        flow_graph = lg.DiGraph(nvertices)
+        for g in testdigraphs(flow_graph)
+            capacity_matrix = zeros(Int, nvertices, nvertices)
+            for e in flow_edges
+                u, v, f = e
+                lg.add_edge!(g, u, v)
+                capacity_matrix[u, v] = f
+            end
+            # Test ExtendedMultirouteFlowAlgorithm when the number of routes is either
+            # Noninteger or 0 (the algorithm returns the breaking points)
+            @test multiroute_flow(g, s, t, capacity_matrix) == breakpts
+            @test multiroute_flow(g, s, t, capacity_matrix, routes=1.5)[1] ≈ ffloat
+            @test multiroute_flow(breakpts, 1.5)[2] ≈ ffloat
+
+            # Test all other algorithms - PR is unstable - see #553
+            for AlgoFlow in [EdmondsKarpAlgorithm, DinicAlgorithm, BoykovKolmogorovAlgorithm, PushRelabelAlgorithm]
+                # When the input are breaking points and routes number
+                @test multiroute_flow(breakpts, 1.5, g, s, t, capacity_matrix)[1] ≈ ffloat
+                for AlgoMrf in [ExtendedMultirouteFlowAlgorithm, KishimotoAlgorithm]
+                    for (k, val) in enumerate(froutes)
+                        @test multiroute_flow(g, s, t, capacity_matrix,
+                            flow_algorithm = AlgoFlow(), mrf_algorithm = AlgoMrf(),
+                            routes = k)[1] ≈ val
+                    end
+                end
+            end
+        end
+    end
+end
diff --git a/test/push_relabel.jl b/test/push_relabel.jl
new file mode 100644
index 0000000..5ab8f14
--- /dev/null
+++ b/test/push_relabel.jl
@@ -0,0 +1,94 @@
+@testset "Push relabel" begin
+    # Construct DiGraph
+    flow_graph = lg.DiGraph(8)
+
+    # Load custom dataset
+    flow_edges = [
+        (1, 2, 10), (1, 3, 5), (1, 4, 15), (2, 3, 4), (2, 5, 9),
+        (2, 6, 15), (3, 4, 4), (3, 6, 8), (4, 7, 16), (5, 6, 15),
+        (5, 8, 10), (6, 7, 15), (6, 8, 10), (7, 3, 6), (7, 8, 10)
+    ]
+
+    capacity_matrix = zeros(Int, 8, 8)
+
+    for e in flow_edges
+        u, v, f = e
+        lg.add_edge!(flow_graph, u, v)
+        capacity_matrix[u, v] = f
+    end
+    for g in testdigraphs(flow_graph)
+      residual_graph = @inferred(LightGraphsFlows.residual(g))
+
+      # Test enqueue_vertex
+      Q = Array{Int,1}()
+      excess = [0, 1, 0, 1]
+      active = [false, false, true, true]
+      @test @inferred(LightGraphsFlows.enqueue_vertex!(Q, 1, active, excess)) == nothing
+      @test @inferred(LightGraphsFlows.enqueue_vertex!(Q, 3, active, excess)) == nothing
+      @test @inferred(LightGraphsFlows.enqueue_vertex!(Q, 4, active, excess)) == nothing
+      @test length(Q) == 0
+      @test @inferred(LightGraphsFlows.enqueue_vertex!(Q, 2, active, excess)) == nothing
+      @test length(Q) == 1
+
+      # Test push_flow
+      Q = Array{Int,1}()
+      excess = [15, 1, 1, 0, 0, 0, 0, 0]
+      height = [8, 0, 0, 0, 0, 0, 0, 0]
+      active = [true, false, false, false, false, false, false, true]
+      flow_matrix = zeros(Int, 8, 8)
+      @test @inferred(LightGraphsFlows.push_flow!(residual_graph, 1, 2, capacity_matrix, flow_matrix, excess, height, active, Q)) == nothing
+      @test length(Q) == 1
+      @test flow_matrix[1, 2] == 10
+      @test @inferred(LightGraphsFlows.push_flow!(residual_graph, 2, 3, capacity_matrix, flow_matrix, excess, height, active, Q)) == nothing
+      @test length(Q) == 1
+      @test flow_matrix[2, 3] == 0
+
+      # Test gap
+      Q = Array{Int,1}()
+      excess = [15, 1, 1, 0, 0, 0, 0, 0]
+      height = [8, 2, 2, 1, 3, 3, 4, 5]
+      active = [true, false, false, false, false, false, false, true]
+      count  = [0, 1, 2, 2, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0]
+      flow_matrix = zeros(Int, 8, 8)
+
+      @test @inferred(LightGraphsFlows.gap!(residual_graph, 1, excess, height, active, count, Q)) == nothing
+      @test length(Q) == 2
+
+      # Test relabel
+      Q = Array{Int,1}()
+      excess = [15, 1, 1, 0, 0, 0, 0, 0]
+      height = [8, 1, 1, 1, 1, 1, 1, 0]
+      active = [true, false, false, false, false, false, false, true]
+      count  = [1, 6, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0]
+      flow_matrix = zeros(Int, 8, 8)
+
+      @test @inferred(LightGraphsFlows.relabel!(residual_graph, 2, capacity_matrix, flow_matrix, excess, height, active, count, Q)) == nothing
+      @test length(Q) == 1
+
+      # Test discharge
+      Q = Array{Int,1}()
+      excess = [50, 1, 1, 0, 0, 0, 0, 0]
+      height = [8, 0, 0, 0, 0, 0, 0, 0]
+      active = [true, false, false, false, false, false, false, true]
+      count  = [7, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0]
+      flow_matrix = zeros(Int, 8, 8)
+
+      @test @inferred(LightGraphsFlows.discharge!(residual_graph, 1, capacity_matrix, flow_matrix, excess, height, active, count, Q)) == nothing
+      @test length(Q) == 3
+
+      # Test with default distances
+      @test LightGraphsFlows.push_relabel(residual_graph, 1, 8, LightGraphsFlows.DefaultCapacity(residual_graph))[1] == 3
+
+      # Test with capacity matrix
+      @test LightGraphsFlows.push_relabel(residual_graph, 1, 8, capacity_matrix)[1] == 28
+    end
+    # Non regression test added for #448
+    M448 = [0 1 0 0 1 1
+            1 0 0 0 1 0
+            0 0 0 1 0 0
+            0 0 0 0 0 0
+            1 0 1 0 0 1
+            0 0 0 0 1 0]
+    g448 = lg.DiGraph(M448)
+    @test maximum_flow(g448, 1, 2, M448, algorithm=PushRelabelAlgorithm())[1] == 1
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 466e658..70faad0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,21 @@
 using LightGraphsFlows
 using Base.Test
+import LightGraphs
+import SimpleTraits
+const lg = LightGraphs
 
-# write your own tests here
-@test 1 == 2
+const testdir = dirname(@__FILE__)
+
+testgraphs(g) = [g, lg.Graph{UInt8}(g), lg.Graph{Int16}(g)]
+testdigraphs(g) = [g, lg.DiGraph{UInt8}(g), lg.DiGraph{Int16}(g)]
+
+for t in [
+        "edmonds_karp",
+        "dinic",
+        "boykov_kolmogorov",
+        "push_relabel",
+        "maximum_flow",
+        "multiroute_flow"]
+    tp = joinpath(testdir, "$(t).jl")
+    include(tp)
+end