diff --git a/doc/index.md b/doc/index.md
index 8df24d4..6eac994 100644
--- a/doc/index.md
+++ b/doc/index.md
@@ -5,3 +5,4 @@ require external packages to function. Currently, the following functionality is
 - Spectral functions (Combinatorial Adjacencies, based on the [GraphMatrices](https://github.com/jpfairbanks/GraphMatrices.jl) package
 - Datasets
 - Community detection
+- Network Interdiction (flow) using [JuMP](https://github.com/JuliaOpt/JuMP.jl)
diff --git a/doc/interdiction.md b/doc/interdiction.md
new file mode 100644
index 0000000..0ef1d45
--- /dev/null
+++ b/doc/interdiction.md
@@ -0,0 +1,130 @@
+# Network Interdiction
+
+## Flow
+
+Network Interdiction is a family of problems where some elements of graph
+(vertices or links usually) are forbidden, destroyed or have failed. The goal is
+then to compute some objective function, as a maximum flow or a shortest path.
+
+Currently, `LightGraphsExtras.jl` implement methods for three kinds of Interdiction
+Flow's problems. The Interdiction flow problems can be seen as a game between two
+players, a defender (that maximizes the flow) and an attacker (that destroy some
+links). Called using respectively `NetworkInterdictionProblem()`,
+`AdaptiveFlowArcProblem()` and `AdaptiveFlowPathProblem()`, the three following
+variants, as described by
+[Bertsimas et al.](http://dx.doi.org/10.1016/j.orl.2015.11.005), are considered :
+
+- The Network Interdiction Flow where the attacker strikes first, then a maximum-flow
+  is computed on the remaining network. The goal is for the defender to predict the
+  best worst case. The problem is called by using
+- Adaptive Network Flow whre the defender computes a flow first. There are two
+  variants.
+
+    - Arc version : the flow is expressed as combination of flow on the links. The
+      destruction of an edge only implies the flow on that to be destroyed (it could
+      results with some flow on other egdes not reaching the target node).
+    - Path version : the flow is now combination of paths from the souce to the sink
+      where each edge that is destroyed also implies the destruction of all the paths
+      going through this edge.
+
+
+This Interdiction Flow (including all the variants) for general graphs is an NP-hard problem and several exact and approximative algorithms have been designed in the academic litterature. Called using respectively `MultilinkAttackAlgorithm()` and `BilevelMixedIntegerLinearProgram()`, `LightGraphsExtras.jl` provides the following methods (any help to implement other algorithms is welcome):
+
+- A deterministic pollynomial algorithm called Multilink Attack relying on the
+  Extended Multiroute Flow algorithm (available in `LightGraphs.jl`) introduced by
+  [Baffier et al.](http://dx.doi.org/10.1016/j.disopt.2016.05.002). For a certain
+  category of graph, it solves Network Interdiction Flow and Adaptive Flow (arc)
+  exactly. When it fails, it provides upper and lower bounds. It provides a
+  2-approximation for Adaptive Flow (path), see
+  [Suppakitpaisarn et al.](http://dx.doi.org/10.1109/HPSR.2015.7483079)
+  for more details.
+- A Bilevel Mixed Integer Linear Program (BMILP) framework using `JuMP.jl` that is
+  guaranteed to converge. (Only the adaptive flow (arc) vairant is covered yet, the
+  others use dummy functions). The results of this framework through
+  `LightGraphsExtras.jl` will be appear in the following month in the litterature.
+
+The `interdiction_flow{T<:AbstractFloat}` function takes the following arguments:
+
+- flow_graph::DiGraph                           # the input graph
+- source::Int                                   # the source vertex
+- target::Int                                   # the target vertex
+- capacity_matrix::AbstractArray{T, 2}          # edge flow capacities
+- attacks::Int                                  # argument for attacks
+- algorithm::AbstractInterdictionFlowAlgorithm  # argument for algorithm,
+- problem::NetworkInterdictionProblem           # argument for problem
+- solver::AbstractMathProgSolver                # argument for solver (BMILP only)
+- rtol::T                                       # relative tolerance (BMILP only)
+- atol::T                                       # absolute tolerance (BMILP only)
+- time_limit::Float64                           # time limit (BMILP only)
+
+If no number of attacks is given, it will be set to -1 (that means for any possible
+number of attacks between 0 and the source-sink connectivity).
+The function defaults to `MultilinkAttackAlgorithm()` and `NetworkInterdictionProblem()`
+for algorithm and problem. If `BilevelMixedIntegerLinearProgram()` is used,
+the solver, rtol, atol, and time_limit are also used and default as respectively
+`UnsetSolver()`, `sqrt(eps())`, `0.`, `Inf`. The relative and absolute tolerance are
+used to compare the upper and lower bounds for termination.
+Please consult the [`JuMP.jl` documentation](http://http://jump.readthedocs.io) for the
+use of the solver keyword, as it is similar to the `Model()` method there.
+
+All algorithms return a tuple with 1) a lower bound and 2) an upper bound.
+For the Multilink Attack algorithm, it also returns the restriction values (to use with
+the function maximum_flow in LightGraphs.jl) associated with 3-a) the lower bound and
+4-a) the upper bound. When the BMILP is used, the third element returned is 3-b) the
+time used by the algorithm.
+
+When the number of attacks is set to -1, an array with the results for any possible number of attacks will be output. Each result will be output as above.
+
+### Approximative comparison
+
+Approximative comparison to deal with Float precision. LaTeX command : `\precapproc`
+
+```julia
+function ⪷{T<:AbstractFloat}(
+  x::T,
+  y::T
+  )
+  return x ≈ y || x < y
+end
+```
+
+### Usage Example :
+
+```julia
+
+# Create a flow-graph and a capacity matrix
+flow_graph = DiGraph(8)
+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(Float64, 8, 8)
+for e in flow_edges
+    u, v, f = e
+    add_edge!(flow_graph, u, v)
+    capacity_matrix[u, v] = f
+end
+
+# Run default values
+lower_bound, upper_bound , lower_restriction, upper_restriction =
+  interdiction_flow(flow_graph, 1, 8, capacity_matrix)
+
+# Run default values but with a specific number of attacks
+lower_bound, upper_bound , lower_restriction, upper_restriction =
+  interdiction_flow(flow_graph, 1, 8, capacity_matrix, 1)
+
+# Run with BMILP for 1 attack and adaptive flow (arc) problem
+lower_bound, upper_bound , time_used =
+  interdiction_flow(flow_graph, 1, 8, capacity_matrix, 1,
+    algorithm = BilevelMixedIntegerLinearProgram(),
+    problem = AdaptiveFlowArcProblem())
+
+if lower_bound ⪷ upper_bound
+  print("Lower bound = $lower_bound and Upper Bound = $upper_bound")
+  println(" in $time_used seconds.")
+else
+  error("There is a bound error in the BMILP")
+end
+
+```
diff --git a/mkdocs.yml b/mkdocs.yml
index acfef0c..3e08c6d 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -6,6 +6,7 @@ pages:
   - 'Getting Started': 'index.md'
   - 'Matching': 'matching.md'
   - 'Spectral': 'spectral.md'
+  - 'Network Interdiction': 'interdiction.md'
 docs_dir: 'doc'
 markdown_extensions:
   - extra
diff --git a/src/LightGraphsExtras.jl b/src/LightGraphsExtras.jl
index a6d7000..0d7ac00 100644
--- a/src/LightGraphsExtras.jl
+++ b/src/LightGraphsExtras.jl
@@ -4,10 +4,12 @@ export Matching
 export CombinatorialAdjacencies
 export Datasets
 export Community
+export Interdiction
 
 include("matching/matching.jl")
 include("spectral/spectral.jl")
 include("datasets/datasets.jl")
 include("community/detection.jl")
+include("interdiction/interdiction.jl")
 
 end # module
diff --git a/src/interdiction/adaptive_arc.jl b/src/interdiction/adaptive_arc.jl
new file mode 100644
index 0000000..47952cb
--- /dev/null
+++ b/src/interdiction/adaptive_arc.jl
@@ -0,0 +1,32 @@
+# Method when the algorithm used is the Multilink Attack algorithm
+function adaptive_arc{T<:AbstractFloat}(
+  flow_graph::DiGraph,                           # the input graph
+  source::Int,                                   # the source vertex
+  target::Int,                                   # the target vertex
+  capacity_matrix::AbstractArray{T, 2},          # edge flow capacities
+  attacks::Int,                                  # argument for attacks
+  algorithm::MultilinkAttackAlgorithm,           # argument for algorithm
+  solver::AbstractMathProgSolver,                # argument for solver (unused)
+  rtol::T,                                       # absolute tolerance (unused)
+  atol::T,                                       # relative tolerance (unused)
+  time_limit::Float64                            # time limit (seconds) (unused)
+  )
+  return multilink_attack(flow_graph, source, target, capacity_matrix, attacks)
+end
+
+# Method when the algorithm used is a Bilevel Mixed Integer Linear Program
+function adaptive_arc{T<:AbstractFloat}(
+  flow_graph::DiGraph,                           # the input graph
+  source::Int,                                   # the source vertex
+  target::Int,                                   # the target vertex
+  capacity_matrix::AbstractArray{T, 2},          # edge flow capacities
+  attacks::Int,                                  # argument for attacks
+  algorithm::BilevelMixedIntegerLinearProgram,   # argument for algorithm
+  solver::AbstractMathProgSolver,                # argument for solver
+  rtol::T,                                       # absolute tolerance
+  atol::T,                                       # relative tolerance
+  time_limit::Float64                            # time limit (seconds)
+  )
+  return bilevel_adaptive_arc(flow_graph, source, target, capacity_matrix,
+                                  attacks, solver, rtol, atol, time_limit)
+end
diff --git a/src/interdiction/adaptive_path.jl b/src/interdiction/adaptive_path.jl
new file mode 100644
index 0000000..c033836
--- /dev/null
+++ b/src/interdiction/adaptive_path.jl
@@ -0,0 +1,36 @@
+# Method when the algorithm used is the Multilink Attack algorithm
+function adaptive_path{T<:AbstractFloat}(
+  flow_graph::DiGraph,                           # the input graph
+  source::Int,                                   # the source vertex
+  target::Int,                                   # the target vertex
+  capacity_matrix::AbstractArray{T, 2},          # edge flow capacities
+  attacks::Int,                                  # argument for attacks
+  algorithm::MultilinkAttackAlgorithm,           # argument for algorithm
+  solver::AbstractMathProgSolver,                # argument for solver (unused)
+  rtol::T,                                       # absolute tolerance (unused)
+  atol::T,                                       # relative tolerance (unused)
+  time_limit::Float64                            # time limit (seconds) (unused)
+  )
+  lower_bound, upper_bound , lower_restriction, upper_restriction =
+    multilink_attack(flow_graph, source, target, capacity_matrix, attacks)
+
+  # MultilinkAttackAlgorithm is a 2-approximation for AdaptiveFlowPathProblem
+  return lower_bound / 2, upper_bound, lower_restriction, upper_restriction
+end
+
+# Method when the algorithm used is a Bilevel Mixed Integer Linear Program
+function adaptive_path{T<:AbstractFloat}(
+  flow_graph::DiGraph,                           # the input graph
+  source::Int,                                   # the source vertex
+  target::Int,                                   # the target vertex
+  capacity_matrix::AbstractArray{T, 2},          # edge flow capacities
+  attacks::Int,                                  # argument for attacks
+  algorithm::BilevelMixedIntegerLinearProgram,   # argument for algorithm
+  solver::AbstractMathProgSolver,                # argument for solver (unused)
+  rtol::T,                                       # absolute tolerance (unused)
+  atol::T,                                       # relative tolerance (unused)
+  time_limit::Float64                            # time limit (seconds) (unused)
+  )
+  return bilevel_adaptive_path(flow_graph, source, target, capacity_matrix,
+                                   attacks, solver, rtol, atol, time_limit)
+end
diff --git a/src/interdiction/bilevel_adaptive_arc.jl b/src/interdiction/bilevel_adaptive_arc.jl
new file mode 100644
index 0000000..4c5ea6b
--- /dev/null
+++ b/src/interdiction/bilevel_adaptive_arc.jl
@@ -0,0 +1,140 @@
+function init_first_level_arc{T<:AbstractFloat}(
+  flow_graph::DiGraph,                          # the input graph
+  source::Int,                                  # the source vertex
+  target::Int,                                  # the target vertex
+  capacity_matrix::AbstractArray{T, 2},         # edge flow capacities
+  solver::AbstractMathProgSolver                # keyword for solver
+  )
+  n = nv(flow_graph)           # size of the network
+
+	first_level = Model(solver = solver)
+
+	# x : first level binary variables
+	@variable(first_level, 0 ≤ x[i = 1:n, j = 1:n] ≤ capacity_matrix[i, j])
+  # z : objective function
+	@variable(first_level, z)
+	@objective(first_level, Max, z)
+
+	# flow conservation constraints
+	for v in vertices(flow_graph)
+		if (v ≠ source && v ≠ target)
+      @constraint(first_level, sum{x[i, j], i = 1:n, j = 1:n;
+        capacity_matrix[i ,j] > 0} - sum{x[j, i], i = 1:n, j = 1:n;
+        capacity_matrix[j, i] > 0} == 0)
+	   	end
+	end
+
+  # objective function upper bound
+	# @constraint(first_level, z ≤ prevfloat(typemax(T)))
+	@constraint(first_level, z ≤ 10000000)
+
+  return first_level, x, z
+end
+
+function init_second_level_arc{T<:AbstractFloat}(
+  flow_graph::DiGraph,                          # the input graph
+  source::Int,                                  # the source vertex
+  target::Int,                                  # the target vertex
+  capacity_matrix::AbstractArray{T, 2},         # edge flow capacities
+  attacks::Int,                                 # argument for attacks
+  solver::AbstractMathProgSolver,               # keyword for solver
+  x::Matrix{JuMP.Variable}                      # flows from first_level
+  )
+  n = nv(flow_graph)           # size of the network
+  x_value = getvalue(x)        # extract the value of variable x
+
+  second_level = Model(solver = solver)
+
+  # variable that model if an arc is attacked or not
+  @variable(second_level, 0 ≤ μ[i = 1:n, j = 1:n] ≤ 1, Int)
+  # first cut variables
+  @variable(second_level, δ[i = 1:n,j = 1:n] ≥ 0)
+  # second cut variables
+  @variable(second_level, σ[i = 1:n] ≥ 0)
+  # linearization variables: ν = δμ
+  @variable(second_level, ν[i = 1:n, j = 1:n] ≥ 0)
+
+  # set objective
+  @objective(second_level, Min, sum{
+    x_value[i, j] * δ[i, j] - x_value[i, j] * ν[i, j], i = 1:n, j = 1:n;
+    capacity_matrix[i, j] > 0})
+
+  # constraints over the edges
+  for e in edges(flow_graph)
+    i, j = src(e), dst(e)
+    # first linearization constraint for v
+    @constraint(second_level, ν[i, j] ≤ μ[i, j])
+    # second linearization constraint for v
+    @constraint(second_level, ν[i, j] ≤ δ[i, j])
+    if i == source
+      # cut constraint for edge (i,j) when i is the source
+      if j != target
+        @constraint(second_level, δ[i, j] + σ[j] ≥ 1 )
+      else
+        @constraint(second_level, δ[i, j] ≥ 1 )
+      end
+    elseif j == target
+      # cut constraint for edge (i,j) when j is the destination
+      if i != source
+         @constraint(second_level, δ[i, j] - σ[i] ≥ 0 )
+      end
+    else
+      # cut constraint for edge (i,j) in the remaining cases
+      @constraint(second_level, δ[i, j] + σ[j] - σ[i] ≥ 0)
+    end
+  end
+
+  # constraint on the upper bound for the number of attacks
+  @constraint(second_level, sum{μ[i, j], i = 1:n, j = 1:n} ≤ attacks)
+
+  return second_level, δ, ν
+end
+
+function bilevel_adaptive_arc{T<:AbstractFloat}(
+  flow_graph::DiGraph,                          # the input graph
+  source::Int,                                  # the source vertex
+  target::Int,                                  # the target vertex
+  capacity_matrix::AbstractArray{T, 2},         # edge flow capacities
+  attacks::Int,                                 # argument for attacks
+  solver::AbstractMathProgSolver,               # keyword for solver
+  rtol::T,                                      # relative tolerance
+  atol::T,                                      # absolute tolerance
+  time_limit::Float64                           # time limit (seconds)
+  )
+	start_time = time()               # time stamp (seconds)
+  n = nv(flow_graph)                # size of the network
+  lower_bound = 0.
+  upper_bound = 0.
+
+  # Initialization : first level
+  first_level, x, z =
+    init_first_level_arc(flow_graph, source, target, capacity_matrix, solver)
+
+  # Loop over the first level while adding cuts from the second level
+  while solve(first_level) ≠ :Infeasible && time() - start_time < time_limit
+    # Store first level results
+    upper_bound, x_value = getobjectivevalue(first_level), getvalue(x)
+
+    # Initialization : second level
+    second_level, δ, ν =
+      init_second_level_arc(flow_graph, source, target, capacity_matrix,
+                                                     attacks, solver, x)
+
+    # Solve second level and update lower_bound
+		solve(second_level)
+		lower_bound = getobjectivevalue(second_level)
+
+    # If the bounds are tight, exit the loop
+    (isapprox(upper_bound, lower_bound, rtol = rtol, atol = atol)) && break
+
+    # Otherwise add the new cut to the first level
+		δ_value = getvalue(δ)
+		ν_value = getvalue(ν)
+    @constraint(first_level, z ≤ sum{
+      δ_value[i, j] * x[i, j] - ν_value[i, j] * x[i, j], i = 1:n, j = 1:n;
+      capacity_matrix[i, j] > 0})
+  end
+
+  # Return objective value and elapsed time
+  return lower_bound, upper_bound, time() - start_time
+end
diff --git a/src/interdiction/bilevel_adaptive_path.jl b/src/interdiction/bilevel_adaptive_path.jl
new file mode 100644
index 0000000..9c57209
--- /dev/null
+++ b/src/interdiction/bilevel_adaptive_path.jl
@@ -0,0 +1,20 @@
+function bilevel_adaptive_path{T<:AbstractFloat}(
+  flow_graph::DiGraph,                          # the input graph
+  source::Int,                                  # the source vertex
+  target::Int,                                  # the target vertex
+  capacity_matrix::AbstractArray{T, 2},         # edge flow capacities
+  attacks::Int,                                 # argument for attacks
+  solver::AbstractMathProgSolver,               # keyword for solver
+  rtol::T,                                      # absolute tolerance
+  atol::T,                                      # relative tolerance
+  time_limit::Float64                           # time limit
+  )
+	start_time = time()                           # time stamp
+  n = nv(flow_graph)                            # size of the network
+  lower_bound = 0.
+
+  # TODO: Write the proper algorithm
+
+  # Return objective value and elapsed time
+  return lower_bound, time() - start_time
+end
diff --git a/src/interdiction/bilevel_network_interdiction.jl b/src/interdiction/bilevel_network_interdiction.jl
new file mode 100644
index 0000000..864f5d9
--- /dev/null
+++ b/src/interdiction/bilevel_network_interdiction.jl
@@ -0,0 +1,20 @@
+function bilevel_network_interdiction{T<:AbstractFloat}(
+  flow_graph::DiGraph,                          # the input graph
+  source::Int,                                  # the source vertex
+  target::Int,                                  # the target vertex
+  capacity_matrix::AbstractArray{T, 2},         # edge flow capacities
+  attacks::Int,                                 # argument for attacks
+  solver::AbstractMathProgSolver,               # keyword for solver
+  rtol::T,                                      # absolute tolerance
+  atol::T,                                      # relative tolerance
+  time_limit::Float64                           # time limit
+  )
+	start_time = time()                           # time stamp
+  n = nv(flow_graph)                            # size of the network
+  lower_bound = 0.
+
+  # TODO: Write the proper algorithm
+
+  # Return objective value and elapsed time
+  return lower_bound, time() - start_time
+end
diff --git a/src/interdiction/interdiction.jl b/src/interdiction/interdiction.jl
new file mode 100644
index 0000000..fea05df
--- /dev/null
+++ b/src/interdiction/interdiction.jl
@@ -0,0 +1,250 @@
+__precompile__(true)
+module Interdiction
+
+using LightGraphs
+using JuMP
+using MathProgBase
+
+import MathProgBase.SolverInterface.AbstractMathProgSolver,
+       JuMP.UnsetSolver,
+       LightGraphs.DefaultCapacity
+
+export interdiction_flow, ⪷,
+       NetworkInterdictionProblem, AdaptiveFlowArcProblem, AdaptiveFlowPathProblem,
+       MultilinkAttackAlgorithm, BilevelMixedIntegerLinearProgram
+
+"""
+Abstract type that allows users to indicate their Problem
+"""
+abstract AbstractInterdictionFlowProblem
+
+"""
+Set the problem version to Network Interdiction (attacker strikes first)
+"""
+type NetworkInterdictionProblem <: AbstractInterdictionFlowProblem
+end
+
+"""
+Set the problem version to Adaptive Network Flow [arc version] (flow is computed first)
+"""
+type AdaptiveFlowArcProblem <: AbstractInterdictionFlowProblem
+end
+
+"""
+Set the problem version to Adaptive Network Flow [path version] (flow is computed first)
+"""
+type AdaptiveFlowPathProblem <: AbstractInterdictionFlowProblem
+end
+
+"""
+Abstract type that allows users to pass in their preferred Algorithm
+"""
+abstract AbstractInterdictionFlowAlgorithm
+
+"""
+Forces the network_interdiction, adaptive_arc, or adaptive_path functions to use the Multilink Attack (MLA) algorithm.
+"""
+type MultilinkAttackAlgorithm <: AbstractInterdictionFlowAlgorithm
+end
+
+"""
+Forces the adaptive_arc function to use a Bilevel Mixed-Integer Linear Program (BMILP)
+"""
+type BilevelMixedIntegerLinearProgram <: AbstractInterdictionFlowAlgorithm
+end
+
+# Includes : algorithms
+include("multilink_attack.jl")
+include("bilevel_adaptive_arc.jl")
+include("bilevel_adaptive_path.jl")
+include("bilevel_network_interdiction.jl")
+# Includes : Problems
+include("adaptive_arc.jl")
+include("network_interdiction.jl")
+include("adaptive_path.jl")
+
+"""
+Approximative comparison to deal with Float precision
+
+```julia
+function ⪷{T<:AbstractFloat}(
+  x::T,
+  y::T
+  )
+  return x ≈ y || x < y
+end
+```
+"""
+function ⪷{T<:AbstractFloat}(
+  x::T,
+  y::T
+  )
+  return x ≈ y || x < y
+end
+
+# Method when the problem considered is Network Interdiction
+function interdiction_flow{T<:AbstractFloat}(
+  flow_graph::DiGraph,                           # the input graph
+  source::Int,                                   # the source vertex
+  target::Int,                                   # the target vertex
+  capacity_matrix::AbstractArray{T, 2},          # edge flow capacities
+  attacks::Int,                                  # argument for attacks
+  algorithm::AbstractInterdictionFlowAlgorithm,  # argument for algorithm,
+  problem::NetworkInterdictionProblem,           # argument for problem
+  solver::AbstractMathProgSolver,                # argument for solver
+  rtol::T,                                       # relative tolerance
+  atol::T,                                       # absolute tolerance
+  time_limit::Float64                            # time limit (seconds)
+  )
+  return network_interdiction(flow_graph, source, target, capacity_matrix,
+                       attacks, algorithm, solver, rtol, atol, time_limit)
+end
+
+# Method when the problem considered is Adaptive Flow (arc version)
+function interdiction_flow{T<:AbstractFloat}(
+  flow_graph::DiGraph,                           # the input graph
+  source::Int,                                   # the source vertex
+  target::Int,                                   # the target vertex
+  capacity_matrix::AbstractArray{T, 2},          # edge flow capacities
+  attacks::Int,                                  # argument for attacks
+  algorithm::AbstractInterdictionFlowAlgorithm,  # argument for algorithm,
+  problem::AdaptiveFlowArcProblem,               # argument for problem
+  solver::AbstractMathProgSolver,                # argument for solver
+  rtol::T,                                       # relative tolerance
+  atol::T,                                       # absolute tolerance
+  time_limit::Float64                            # time limit (seconds)
+  )
+  return adaptive_arc(flow_graph, source, target, capacity_matrix, attacks,
+                                 algorithm, solver, rtol, atol, time_limit)
+end
+
+# Method when the problem considered is Adaptive Flow (path version)
+function interdiction_flow{T<:AbstractFloat}(
+  flow_graph::DiGraph,                           # the input graph
+  source::Int,                                   # the source vertex
+  target::Int,                                   # the target vertex
+  capacity_matrix::AbstractArray{T, 2},          # edge flow capacities
+  attacks::Int,                                  # argument for attacks
+  algorithm::AbstractInterdictionFlowAlgorithm,  # argument for algorithm,
+  problem::AdaptiveFlowPathProblem,              # argument for problem
+  solver::AbstractMathProgSolver,                # argument for solver
+  rtol::T,                                       # relative tolerance
+  atol::T,                                       # absolute tolerance
+  time_limit::Float64                            # time limit (seconds)
+  )
+  return adaptive_path(flow_graph, source, target, capacity_matrix, attacks,
+                                  algorithm, solver, rtol, atol, time_limit)
+end
+
+"""
+# Network Interdiction
+
+## Flow
+
+Network Interdiction is a family of problems where some elements of graph
+(vertices or links usually) are forbidden, destroyed or have failed. The goal is
+then to compute some objective function, as a maximum flow or a shortest path.
+
+Currently, `LightGraphsExtras.jl` implement methods for three kinds of Interdiction
+Flow's problems. The Interdiction flow problems can be seen as a game between two
+players, a defender (that maximizes the flow) and an attacker (that destroy some
+links). Called using respectively `NetworkInterdictionProblem()`,
+`AdaptiveFlowArcProblem()` and `AdaptiveFlowPathProblem()`, the three following
+variants, as described by
+[Bertsimas et al.](http://dx.doi.org/10.1016/j.orl.2015.11.005), are considered :
+
+- The Network Interdiction Flow where the attacker strikes first, then a maximum-flow
+  is computed on the remaining network. The goal is for the defender to predict the
+  best worst case. The problem is called by using
+- Adaptive Network Flow whre the defender computes a flow first. There are two
+  variants.
+
+    - Arc version : the flow is expressed as combination of flow on the links. The
+      destruction of an edge only implies the flow on that to be destroyed (it could
+      results with some flow on other egdes not reaching the target node).
+    - Path version : the flow is now combination of paths from the souce to the sink
+      where each edge that is destroyed also implies the destruction of all the paths
+      going through this edge.
+
+
+This Interdiction Flow (including all the variants) for general graphs is an NP-hard problem and several exact and approximative algorithms have been designed in the academic litterature. Called using respectively `MultilinkAttackAlgorithm()` and `BilevelMixedIntegerLinearProgram()`, `LightGraphsExtras.jl` provides the following methods (any help to implement other algorithms is welcome):
+
+- A deterministic pollynomial algorithm called Multilink Attack relying on the
+  Extended Multiroute Flow algorithm (available in `LightGraphs.jl`) introduced by
+  [Baffier et al.](http://dx.doi.org/10.1016/j.disopt.2016.05.002). For a certain
+  category of graph, it solves Network Interdiction Flow and Adaptive Flow (arc)
+  exactly. When it fails, it provides upper and lower bounds. It provides a
+  2-approximation for Adaptive Flow (path), see
+  [Suppakitpaisarn et al.](http://dx.doi.org/10.1109/HPSR.2015.7483079)
+  for more details.
+- A Bilevel Mixed Integer Linear Program (BMILP) framework using `JuMP.jl` that is
+  guaranteed to converge. (Only the adaptive flow (arc) vairant is covered yet, the
+  others use dummy functions). The results of this framework through
+  `LightGraphsExtras.jl` will be appear in the following month in the litterature.
+
+The `interdiction_flow{T<:AbstractFloat}` function takes the following arguments:
+
+- flow_graph::DiGraph                           # the input graph
+- source::Int                                   # the source vertex
+- target::Int                                   # the target vertex
+- capacity_matrix::AbstractArray{T, 2}          # edge flow capacities
+- attacks::Int                                  # argument for attacks
+- algorithm::AbstractInterdictionFlowAlgorithm  # argument for algorithm,
+- problem::NetworkInterdictionProblem           # argument for problem
+- solver::AbstractMathProgSolver                # argument for solver (BMILP only)
+- rtol::T                                       # relative tolerance (BMILP only)
+- atol::T                                       # absolute tolerance (BMILP only)
+- time_limit::Float64                           # time limit (BMILP only)
+
+If no number of attacks is given, it will be set to -1 (that means for any possible
+number of attacks between 0 and the source-sink connectivity).
+The function defaults to `MultilinkAttackAlgorithm()` and `NetworkInterdictionProblem()`
+for algorithm and problem. If `BilevelMixedIntegerLinearProgram()` is used,
+the solver, rtol, atol, and time_limit are also used and default as respectively
+`UnsetSolver()`, `sqrt(eps())`, `0.`, `Inf`. The relative and absolute tolerance are
+used to compare the upper and lower bounds for termination.
+Please consult the [`JuMP.jl` documentation](http://http://jump.readthedocs.io) for the
+use of the solver keyword, as it is similar to the `Model()` method there.
+
+All algorithms return a tuple with 1) a lower bound and 2) an upper bound.
+For the Multilink Attack algorithm, it also returns the restriction values (to use with
+the function maximum_flow in LightGraphs.jl) associated with 3-a) the lower bound and
+4-a) the upper bound. When the BMILP is used, the third element returned is 3-b) the
+time used by the algorithm.
+
+When the number of attacks is set to -1, an array with the results for any possible number of attacks will be output. Each result will be output as above.
+
+"""
+
+function interdiction_flow{T<:AbstractFloat}(
+  flow_graph::DiGraph,                           # the input graph
+  source::Int,                                   # the source vertex
+  target::Int,                                   # the target vertex
+  capacity_matrix::AbstractArray{T, 2},          # edge flow capacities
+  attacks::Int = -1;                             # argument for attacks
+  algorithm::AbstractInterdictionFlowAlgorithm = # keyword argument for algorithm
+    MultilinkAttackAlgorithm(),
+  problem::AbstractInterdictionFlowProblem =     # keyword argument for problem
+    NetworkInterdictionProblem(),
+  solver::AbstractMathProgSolver =               # keyword for solver
+    UnsetSolver(),
+  rtol::T = sqrt(eps()),                         # relative tolerance
+  atol::T = 0.,                                  # absolute tolerance
+  time_limit::Float64 = Inf                      # time limit (seconds)
+  )
+  # attacks ≥ λ (connectivity) → f = 0
+  λ = maximum_flow(flow_graph, source, target, DefaultCapacity(flow_graph), algorithm = EdmondsKarpAlgorithm())[1]
+  (attacks ≥ λ) && return 0., 0., 0., 0.
+
+  # compute the result when 0 ≤ attacks < λ
+  (attacks ≥ 0) && return interdiction_flow(flow_graph, source, target, capacity_matrix,
+                       attacks, algorithm, problem, solver, rtol, atol, time_limit)
+
+  # Iteration over all the relevant number of attacks k ∈ 0:λ
+  return [
+    interdiction_flow(flow_graph, source, target, capacity_matrix, k, algorithm,
+                                              problem, solver, rtol, atol, time_limit)
+    for k in 0:λ]
+end
+
+end #module
diff --git a/src/interdiction/multilink_attack.jl b/src/interdiction/multilink_attack.jl
new file mode 100644
index 0000000..625419e
--- /dev/null
+++ b/src/interdiction/multilink_attack.jl
@@ -0,0 +1,54 @@
+# Compute the projection of (x,y) on the y-axis along a given slope
+function projection{T<:AbstractFloat}(
+  x::T, y::T,       # Coordinates
+  slope::Int                    # Associated slope
+  )
+  return y - slope * x
+end
+
+# Compute the upper and lower bounds for a given number of attacks
+function bounds_from_points{T<:AbstractFloat}(
+  breaking_points::Vector{Tuple{T,T,Int}},       # Set of breaking points
+  attacks::Int                                   # argument for attacks
+  )
+  # Indices for upper and lower bounds
+  upper_index = 0
+  lower_index = 0
+  # Loop to find the correct indices
+  for (id, point) in enumerate(breaking_points)
+    if attacks ≤ point[3]
+      lower_index = id
+    end
+    if attacks ≥ point[3]
+      upper_index = id
+      break
+    end
+  end
+
+  # Assign the points corresponding to the upper and lower bounds
+  lower_point = breaking_points[lower_index]
+  upper_point = breaking_points[upper_index]
+
+  # A bound is described by its value and an associated restriction
+  # The restriction is useful to reconstruct a bound's flows and cuts
+  lower_value = projection(lower_point[1], lower_point[2], lower_point[3])
+  upper_value = projection(upper_point[1], upper_point[2], upper_point[3])
+  lower_restriction, upper_restriction = lower_point[3], upper_point[3]
+
+  return lower_value, upper_value, lower_restriction, upper_restriction
+end
+
+
+function multilink_attack{T<:AbstractFloat}(
+  flow_graph::DiGraph,                           # the input graph
+  source::Int,                                   # the source vertex
+  target::Int,                                   # the target vertex
+  capacity_matrix::AbstractArray{T, 2},          # edge flow capacities
+  attacks::Int                                   # argument for attacks
+  )
+  # Get all the breaking points from extended multiroute flow algorithm
+  breaking_points = multiroute_flow(flow_graph, source, target, capacity_matrix, flow_algorithm = EdmondsKarpAlgorithm())
+
+  # Return the lower and upper bounds : value, restriction
+  return bounds_from_points(breaking_points, attacks)
+end
diff --git a/src/interdiction/network_interdiction.jl b/src/interdiction/network_interdiction.jl
new file mode 100644
index 0000000..30f47bf
--- /dev/null
+++ b/src/interdiction/network_interdiction.jl
@@ -0,0 +1,33 @@
+# Method when the algorithm used is the Multilink Attack algorithm
+function network_interdiction{T<:AbstractFloat}(
+  flow_graph::DiGraph,                           # the input graph
+  source::Int,                                   # the source vertex
+  target::Int,                                   # the target vertex
+  capacity_matrix::AbstractArray{T, 2},          # edge flow capacities
+  attacks::Int,                                  # argument for attacks
+  algorithm::MultilinkAttackAlgorithm,           # argument for algorithm
+  solver::AbstractMathProgSolver,                # argument for solver (unused)
+  rtol::T,                                       # absolute tolerance (unused)
+  atol::T,                                       # relative tolerance (unused)
+  time_limit::Float64                            # time limit (seconds) (unused)
+  )
+  return multilink_attack(flow_graph, source, target, capacity_matrix, attacks)
+end
+
+
+# Method when the algorithm used is a Bilevel Mixed Integer Linear Program
+function network_interdiction{T<:AbstractFloat}(
+  flow_graph::DiGraph,                           # the input graph
+  source::Int,                                   # the source vertex
+  target::Int,                                   # the target vertex
+  capacity_matrix::AbstractArray{T, 2},          # edge flow capacities
+  attacks::Int,                                  # argument for attacks
+  algorithm::BilevelMixedIntegerLinearProgram,   # argument for algorithm
+  solver::AbstractMathProgSolver,                # argument for solver
+  rtol::T,                                       # absolute tolerance
+  atol::T,                                       # relative tolerance
+  time_limit::Float64                            # time limit (seconds)
+  )
+  return bilevel_network_interdiction(flow_graph, source, target, capacity_matrix,
+                                          attacks, solver, rtol, atol, time_limit)
+end
diff --git a/test/interdiction/runtests.jl b/test/interdiction/runtests.jl
new file mode 100644
index 0000000..479df28
--- /dev/null
+++ b/test/interdiction/runtests.jl
@@ -0,0 +1,85 @@
+using Base.Test
+using LightGraphs
+using LightGraphsExtras.Interdiction
+
+#### 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
+  ),
+
+  # 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
+  ),
+
+  # 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
+  ),
+
+  # Graph with 6 vertices
+  (6,
+   [
+     (1, 2, 1), (1, 6, 1), (1, 3, 1), (1, 4, 2), (1, 5, 2),
+     (2, 6, 1), (3, 6, 1), (4, 6, 2), (5, 6, 2), (6, 1, 1),
+   ],
+   1, 6                                      # source/target
+  )
+]
+
+for (nvertices, flow_edges, s, t) in graphs
+  flow_graph = DiGraph(nvertices)
+  capacity_matrix = zeros(nvertices, nvertices)
+  for e in flow_edges
+    u, v, f = e
+    add_edge!(flow_graph, u, v)
+    capacity_matrix[u, v] = f
+  end
+
+  # Connectivity λ
+  λ = maximum_flow(flow_graph, s, t)[1]
+  for k in 0:λ
+    # Compare MultilinkAttackAlgorithm and BilevelMixedIntegerLinearProgram
+    # a) TODO: When the problem is NetworkInterdictionProblem
+    mla = interdiction_flow(flow_graph, s, t, capacity_matrix, k)
+    @test mla[1] ⪷ mla[2]
+    bmilp = interdiction_flow(flow_graph, s, t, capacity_matrix, k,
+                              problem = NetworkInterdictionProblem(),
+                              algorithm = BilevelMixedIntegerLinearProgram())
+    @test bmilp[1] ≈ 0.
+    # b) When the problem is AdaptiveFlowArcProblem
+    mla = interdiction_flow(flow_graph, s, t, capacity_matrix, k,
+                            problem = AdaptiveFlowArcProblem())
+    bmilp = interdiction_flow(flow_graph, s, t, capacity_matrix, k,
+                              problem = AdaptiveFlowArcProblem(),
+                              algorithm = BilevelMixedIntegerLinearProgram())
+    @test mla[1] ⪷ bmilp[1] ⪷ bmilp[2] ⪷ mla[2]
+    # c) TODO: When the problem is AdaptiveFlowPathProblem
+    mla = interdiction_flow(flow_graph, s, t, capacity_matrix, k,
+                            problem = AdaptiveFlowPathProblem())
+    @test mla[1] * 2. ⪷ mla[2]
+    bmilp = interdiction_flow(flow_graph, s, t, capacity_matrix, k,
+                              problem = AdaptiveFlowPathProblem(),
+                              algorithm = BilevelMixedIntegerLinearProgram())
+    @test bmilp[1] ≈ 0.
+  end
+  # Test when attacks = -1
+  @test interdiction_flow(flow_graph, s, t, capacity_matrix)[2] ==
+        interdiction_flow(flow_graph, s, t, capacity_matrix, 1)
+
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 7077148..760ec5c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,6 +4,7 @@ using LightGraphsExtras
 using LightGraphsExtras.Matching
 using LightGraphsExtras.Datasets
 using LightGraphsExtras.Community
+using LightGraphsExtras.Interdiction
 using Base.Test
 
 testdir = dirname(@__FILE__)
@@ -12,7 +13,8 @@ tests = [
     "matching/runtests",
     "spectral/runtests",
     "datasets/runtests",
-    "community/detection"
+    "community/detection",
+    "interdiction/runtests"
 ]
 
 for t in tests