diff --git a/Project.toml b/Project.toml
index 8e9a517..5827dbb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -19,7 +19,7 @@ Documenter = "0.27"
 FillArrays = "0.13, 1.15"
 Graphs = "1.7"
 HiGHS = "1"
-JuMP = "1"
+JuMP = "1.3"
 JuliaFormatter = "1"
 MathOptInterface = "1.18"
 OptimalTransport = "0.3"
diff --git a/docs/src/algorithms.md b/docs/src/algorithms.md
index 7131a94..9c51633 100644
--- a/docs/src/algorithms.md
+++ b/docs/src/algorithms.md
@@ -98,6 +98,63 @@ graph_matching
 GraphsOptim.graph_matching_step_size
 ```
 
+## Graph edit distance
+```@docs
+edit_distance
+GraphsOptim.edit_distance!
+```
+
+### Mathematical Details
+
+The formulations implemented in the code are described in `D'ascenzo, Andrea, et al. "Enhancing Graph Edit Distance Computation: Stronger and Orientation-based ILP Formulations." Proceedings of the VLDB Endowment 18.11 (2025): 4737-4749`. The paper also containts a more detailed description of the graph edit distance problem.
+
+In the graph edit distance problem, we search for a minimum cost *edit path* between
+two graphs $G$ and $H$, that is a sequence of edit operations transforming $G$ to $H$. Valid
+edit operations are inserting or deleting a node or edge, or relabeling an existing node or
+edge.
+
+This formulation of the problem is easy to visualize, but not very suitable for implementing
+with integer programming. We thus use another, equivalent definition instead: a *node map*
+is a relation $\pi \subset V_{G+ \epsilon} \times V_{H + \epsilon}$ on the vertex sets
+augmented by $\epsilon$, in which 
+
+* each node $v \in V_G$ is mapped to exactly one element of $V_{H + \epsilon}$
+* each node $w \in V_H$ has exactly one preimage in $V_{G + \epsilon}$
+
+Being mapped to $\epsilon$ is equivalent to being deleted or created. We will not give a
+proof here, but for metric cost functions this is equivalent to the edit path formulation
+(and any cost function can be transformed into an equivalent metric one).
+
+This version of the problem is much better suited for integer programming, and is the
+foundation for all formulations implemented here. The formulations all use the same basic
+principle, defining binary variables to encode the mapping between nodes and edges,
+constraints that ensure the variables encode valid node and edge maps, and finally *topology
+constraints* which link node and edge variables based on the topology of the input graphs.
+The formulations differ in the exact variable layout and more importantly the topology
+constraints used. As seen in the paper mentioned above, the different topology constraints
+not only provide massive speedups in practice but their relaxations provide different
+bounds.
+
+The best formulation both in theory and practive orients the graphs to provide stricter
+topology constraints. The graph $G$ is oriented in a canonical way, while in $H$ there are
+forward and backwards edges for each undirected edge in the input. We can then additionally
+demand that edges be mapped consistently with their orientation.
+
+### Function documentation for advanced use
+```@docs
+GraphsOptim.EditCosts
+GraphsOptim.get_default_edit_costs
+GraphsOptim.validate_cost_function
+GraphsOptim.Formulation
+GraphsOptim.F1
+GraphsOptim.F1prime
+GraphsOptim.F1plus
+GraphsOptim.F2minus
+GraphsOptim.F2
+GraphsOptim.F2plus
+GraphsOptim.FORI
+```
+
 ## Coloring
 
 ```@docs
diff --git a/src/GraphsOptim.jl b/src/GraphsOptim.jl
index 7156d42..584f97a 100644
--- a/src/GraphsOptim.jl
+++ b/src/GraphsOptim.jl
@@ -8,6 +8,7 @@ module GraphsOptim
 using Graphs: AbstractGraph, is_directed
 using Graphs: vertices, edges, nv, ne, src, dst, inneighbors, outneighbors, has_edge
 using Graphs: complement, maximal_cliques
+using Graphs: SimpleGraph
 using FillArrays: Zeros, Ones, Fill
 using HiGHS: HiGHS
 using JuMP: Model, AffExpr
@@ -15,6 +16,7 @@ using JuMP: objective_function, add_to_expression!
 using JuMP: set_silent, optimize!, termination_status, value
 using JuMP: set_optimizer, objective_value
 using JuMP: @variable, @constraint, @objective
+using JuMP: VariableRef
 using LinearAlgebra: norm, tr, dot
 using MathOptInterface: OPTIMAL
 using SparseArrays: sparse
@@ -29,6 +31,8 @@ export fractional_chromatic_number, fractional_clique_number
 export shortest_path
 export maximum_weight_clique
 
+export edit_distance, F1, F1prime, F1plus, F2, F2minus, F2plus, FORI
+
 include("utils.jl")
 include("flow.jl")
 include("assignment.jl")
@@ -38,5 +42,6 @@ include("fractional_coloring.jl")
 include("shortest_path.jl")
 include("maximum_clique.jl")
 include("independent_set.jl")
+include("graph_edit_distance.jl")
 
 end
diff --git a/src/graph_edit_distance.jl b/src/graph_edit_distance.jl
new file mode 100644
index 0000000..5df8cdc
--- /dev/null
+++ b/src/graph_edit_distance.jl
@@ -0,0 +1,570 @@
+"""
+Control flow types to communicate the ILP formulation to the edit distance code. The
+subtypes are the different formulations, with `FORI` being the best performing one. 
+"""
+abstract type Formulation end
+"""
+Control flow type signifying the F1 ILP formulation. See also [`Formulation`](@ref).
+"""
+struct F1 <: Formulation end
+"""
+Control flow type signifying the F1' ILP formulation. See also [`Formulation`](@ref).
+"""
+struct F1prime <: Formulation end
+"""
+Control flow type signifying the F1+ ILP formulation. See also [`Formulation`](@ref).
+"""
+struct F1plus <: Formulation end
+"""
+Control flow type signifying the F2- ILP formulation. See also [`Formulation`](@ref).
+"""
+struct F2minus <: Formulation end
+"""
+Control flow type signifying the F2 ILP formulation. See also [`Formulation`](@ref).
+"""
+struct F2 <: Formulation end
+"""
+Control flow type signifying the F2+ ILP formulation. See also [`Formulation`](@ref).
+"""
+struct F2plus <: Formulation end
+"""
+Control flow type signifying the FORI ILP formulation. See also [`Formulation`](@ref).
+"""
+struct FORI <: Formulation end
+
+"""
+Struct maintaining the variable references for formulations [`F2minus`](@ref), [`F2`](@ref)
+and [`F2plus`](@ref). Compared to [`FullVariables`](@ref), the deletion variables are
+inferred.
+"""
+struct ReducedVariables
+    x::Matrix{VariableRef}
+    y
+end
+
+"""
+Struct maintaining the variable references for formulations [`F1`](@ref), [`F1prime`](@ref)
+and [`F1plus`](@ref).
+"""
+struct FullVariables
+    x::Matrix{VariableRef}
+    y
+    nodeDelG::Vector{VariableRef}
+    nodeDelH::Vector{VariableRef}
+    edgeDelG
+    edgeDelH
+end
+
+"""
+Struct maintaining the variable references for formulation [`FORI`](@ref). This formulation
+implicitly adds two oriented edges for each edge the second graph to improve topology
+constraints, which is reflected in the edge map variables.
+"""
+struct OrientedVariables
+    x::Matrix{VariableRef}
+    z
+end
+
+"""
+Type used for functions applicable to both F1 and F2 derived formulations.
+"""
+Variables = Union{ReducedVariables,FullVariables}
+
+"""
+    create_model_vars_reduced!(model, G, H, bidirectional = false)
+
+Adds necessary variables for edit distance computation between `G` and `H` to `model`. If
+`bidirectional` is false, for formulations [`F2minus`](@ref), [`F2`](@ref),
+[`F2plus`](@ref). If `bidirectional` is true, `H` is bidirected and the variables are used
+for the [`FORI`](@ref) formulation.
+
+The variables model a full mapping between nodes and edges respectively.
+
+# Returns
+- [`OrientedVariables`](@ref) or [`ReducedVariables`](@ref) containing the variables.
+"""
+function create_model_vars_reduced!(
+    model::Model, G::AbstractGraph, H::AbstractGraph, bidirectional::Bool=false
+)
+    @variable(model, x[1:nv(G), 1:nv(H)], Bin)
+
+    # Use a sparse indexed edge variable set, to later do partial sums y[i, :, k, l].
+    # Edge variables are always oriented to avoid ambiguity.
+    @variable(
+        model,
+        y[
+            i in vertices(G),
+            j in vertices(G),
+            k in vertices(H),
+            l in vertices(H);
+            has_edge(G, i, j) && has_edge(H, k, l) && i < j && (k < l || bidirectional),
+        ],
+        Bin
+    )
+
+    if bidirectional
+        return OrientedVariables(x, y)
+    else
+        return ReducedVariables(x, y)
+    end
+end
+
+"""
+Convenience function bind. See [`create_model_vars_reduced`](@ref).
+"""
+function create_model_vars_bidirectional!(model::Model, G::AbstractGraph, H::AbstractGraph)
+    return create_model_vars_reduced!(model, G, H, true)
+end
+
+"""
+    create_model_vars_full!(model, G, H, bidirectional = false)
+
+Adds necessary variables for edit distance computation between `G` and `H` to `model` for
+formulations [`F1`](@ref), [`F1prime`](@ref) and [`F1plus`](@ref).
+
+In addition to the normal node and edge map variables (see
+[`create_model_vars_reduced`](@ref)), these formulations explicitly have variables to model
+nodes and edges being deleted or added.
+"""
+function create_model_vars_full!(model::Model, G::AbstractGraph, H::AbstractGraph)
+    vars = create_model_vars_reduced!(model, G, H)
+    @variable(model, nodeDelG[1:nv(G)], Bin)
+    @variable(model, nodeDelH[1:nv(H)], Bin)
+
+    @variable(model, edgeDelG[i in vertices(G), j in vertices(G); has_edge(G, i, j)], Bin)
+    @variable(model, edgeDelH[k in vertices(H), l in vertices(H); has_edge(H, k, l)], Bin)
+
+    return FullVariables(vars.x, vars.y, nodeDelG, nodeDelH, edgeDelG, edgeDelH)
+end
+
+"""
+Models arbitrary edit costs between two graphs. The implementation is taken from
+[FORI-GED](https://github.com/meffertj/FORI-GED).
+
+For details on valid variable dimensions see [`validate_cost_function`](@ref), for
+constructing the canonical cost function see [`get_default_edit_costs`](@ref).
+"""
+struct EditCosts
+    c_ik::Array{Number,2}
+    c_iε::Vector{Number}
+    c_εk::Vector{Number}
+    # note that for sparse graphs we might want to use a sparse array here
+    c_ijkl::Array{Number,4}
+    c_ijε::Array{Number,2}
+    c_εkl::Array{Number,2}
+end
+
+"""
+    validate_cost_function(c, G, H)
+
+Validates that cost function `c` has appropriate dimensions for graphs `G` and `H`. This 
+means
+
+- `c.c_ik` has a cost for each mapping `i ∈ V(G)` to `k ∈ V(H)`
+- `c.c_iε` has a cost for deleting each `i ∈ V(G)`
+- `c.c_εk` has a cost for creating each `k ∈ V(H)`
+- `c.c_ijkl` has a cost for each mapping from `ij ∈ E(G)` to  `kl ∈ E(H)`
+- `c.c_ijε` has a cost for deleting each `ij ∈ E(G)`
+- `c.c_εkl` has a cost for creating each `kl ∈ E(H)`
+"""
+function validate_cost_function(c::EditCosts, G::AbstractGraph, H::AbstractGraph)
+    @assert size(c.c_ik) == (nv(G), nv(H))
+    @assert size(c.c_iε) == (nv(G),)
+    @assert size(c.c_εk) == (nv(H),)
+    @assert size(c.c_ijkl) == (nv(G), nv(G), nv(H), nv(H))
+    @assert size(c.c_ijε) == (nv(G), nv(G))
+    @assert size(c.c_εkl) == (nv(H), nv(H))
+    return nothing
+end
+
+"""
+Returns the cost function where mappings are free, deleting or creating things costs
+uniformly 1. This corresponds to "number of edits needed to go from `G` to `H`".
+"""
+function get_default_edit_costs(G::AbstractGraph, H::AbstractGraph)
+    return EditCosts(
+        zeros(Int, nv(G), nv(H)),
+        ones(Int, nv(G)),
+        ones(Int, nv(H)),
+        zeros(Int, nv(G), nv(G), nv(H), nv(H)),
+        ones(Int, nv(G), nv(G)),
+        ones(Int, nv(H), nv(H)),
+    )
+end
+
+"""
+Adds objective function for F1 style formulations to `model`.
+"""
+function add_F1_objective!(
+    model, c::EditCosts, vars::FullVariables, G::AbstractGraph, H::AbstractGraph
+)
+    @objective(
+        model,
+        Min,
+        sum(vars.x[i, k] * c.c_ik[i, k] for i in 1:nv(G) for k in 1:nv(H)) +
+            sum(vars.nodeDelG[i] * c.c_iε[i] for i in 1:nv(G)) +
+            sum(vars.nodeDelH[k] * c.c_εk[k] for k in 1:nv(H)) +
+            sum(
+                vars.y[i, j, k, l] * c.c_ijkl[i, j, k, l] for i in 1:nv(G) for j in 1:nv(G)
+                for k in 1:nv(H) for
+                l in 1:nv(H) if has_edge(G, i, j) && has_edge(H, k, l) && i < j && k < l
+            ) +
+            sum(
+                vars.edgeDelG[i, j] * c.c_ijε[i, j] for i in 1:nv(G) for
+                j in 1:nv(G) if has_edge(G, i, j) && i < j
+            ) +
+            sum(
+                vars.edgeDelH[k, l] * c.c_εkl[k, l] for k in 1:nv(H) for
+                l in 1:nv(H) if has_edge(H, k, l) && k < l
+            )
+    )
+    return nothing
+end
+
+"""
+Adds objective function for F2 style formulations to `model`. `bidirectional` is used to
+support input for [`FORI`](@ref) variables.
+"""
+function add_F2_objective!(
+    model,
+    c::EditCosts,
+    vars::ReducedVariables,
+    G::AbstractGraph,
+    H::AbstractGraph,
+    bidirectional::Bool=false,
+)
+    # Note that init=0 is necessary to allow for empty graph edge cases.
+    K =
+        sum(c.c_iε[i] for i in 1:nv(G); init=0) +
+        sum(c.c_εk[k] for k in 1:nv(H); init=0) +
+        sum(
+            c.c_ijε[i, j] for i in 1:nv(G) for j in 1:nv(G) if has_edge(G, i, j) && i < j;
+            init=0,
+        ) +
+        sum(
+            c.c_εkl[k, l] for k in 1:nv(H) for l in 1:nv(H) if has_edge(H, k, l) && k < l;
+            init=0,
+        )
+    @objective(
+        model,
+        Min,
+        sum(
+                vars.x[i, k] * (c.c_ik[i, k] - c.c_iε[i] - c.c_εk[k]) for i in 1:nv(G) for
+                k in 1:nv(H);
+                init=0,
+            ) +
+            sum(
+                vars.y[i, j, k, l] * (c.c_ijkl[i, j, k, l] - c.c_ijε[i, j] - c.c_εkl[k, l])
+                for i in 1:nv(G) for j in 1:nv(G) for k in 1:nv(H) for l in 1:nv(H) if
+                has_edge(G, i, j) && has_edge(H, k, l) && i < j && (k < l || bidirectional);
+                init=0,
+            ) +
+            K
+    )
+    return nothing
+end
+
+"""
+Adds objective function for [`FORI`](@ref) formulation variables. The implementation is in
+[`add_F2_objective`](@ref).
+"""
+function add_FORI_objective!(
+    model, c::EditCosts, vars::OrientedVariables, G::AbstractGraph, H::AbstractGraph
+)
+    return add_F2_objective!(model, c, ReducedVariables(vars.x, vars.z), G, H, true)
+end
+
+"""
+Add node map constraints to model for F1 style formulations, i.e. each node is only mapped
+to exactly one other node or deleted/created.
+"""
+function add_node_map_constraints!(model::Model, vars::FullVariables, G, H)
+    @constraint(model, [i in 1:nv(G)], sum(vars.x[i, :]) + vars.nodeDelG[i] == 1)
+    @constraint(model, [j in 1:nv(H)], sum(vars.x[:, j]) + vars.nodeDelH[j] == 1)
+    return nothing
+end
+
+"""
+Add node map constraints to model for F2 style and FORI formulations, i.e. each node is only
+mapped to at most one other node (not being mapped implies being deleted or created).
+"""
+function add_node_map_constraints!(
+    model::Model, vars::Union{ReducedVariables,OrientedVariables}, G, H
+)
+    @constraint(model, [i in 1:nv(G)], sum(vars.x[i, :]) <= 1)
+    @constraint(model, [j in 1:nv(H)], sum(vars.x[:, j]) <= 1)
+    return nothing
+end
+
+"""
+Adds constraints on the edge map variables so each edge is mapped to exactly one other or be
+deleted/created. Note that these constraints are only necessary in F1 style formulations, as
+they are implied in F2 styled and FORI.
+"""
+function add_edge_map_constraints!(model::Model, vars::FullVariables, G, H)
+    # Note that init=0 is necessary to allow for empty graph edge cases.
+    @constraint(
+        model,
+        [i in vertices(G), j in vertices(G); has_edge(G, i, j) && i < j],
+        sum(vars.y[i, j, :, :]; init=0) + vars.edgeDelG[i, j] == 1
+    )
+    @constraint(
+        model,
+        [k in vertices(H), l in vertices(H); has_edge(H, k, l) && k < l],
+        sum(vars.y[:, :, k, l]; init=0) + vars.edgeDelH[k, l] == 1
+    )
+    return nothing
+end
+
+"""
+Add simplest form of topology constraints: If edge `ij` is mapped to `kl`, then `i` or `j`
+must map to `k`, and `i` or `j` must map to `l`.
+"""
+function add_simple_topology_constraints!(model::Model, vars::Variables, G, H)
+    @constraint(
+        model,
+        [
+            i in vertices(G),
+            j in vertices(G),
+            k in vertices(H),
+            l in vertices(H);
+            has_edge(G, i, j) && has_edge(H, k, l) && i < j && k < l,
+        ],
+        vars.y[i, j, k, l] <= vars.x[i, k] + vars.x[j, k]
+    )
+    @constraint(
+        model,
+        [
+            i in vertices(G),
+            j in vertices(G),
+            k in vertices(H),
+            l in vertices(H);
+            has_edge(G, i, j) && has_edge(H, k, l) && i < j && k < l,
+        ],
+        vars.y[i, j, k, l] <= vars.x[i, l] + vars.x[j, l]
+    )
+    return nothing
+end
+
+"""
+Add improved topology constraints, replacing [`add_simple_topology_constraints`](@ref): 
+If edge `ij` is mapped to any edge incident to `k`, then `i` or `j` must be mapped to `k`.
+"""
+function add_improved_topology_constraints_G_to_H!(model::Model, vars::Variables, G, H)
+    @constraint(
+        model,
+        [i in vertices(G), j in vertices(G), k in vertices(H); has_edge(G, i, j) && i < j],
+        sum(vars.y[i, j, k, :]; init=0) + sum(vars.y[i, j, :, k]; init=0) <=
+            vars.x[i, k] + vars.x[j, k]
+    )
+    return nothing
+end
+
+"""
+Add topology constraints mirroring [`add_improved_topology_constraints_G_to_H`](@ref), but
+backwards: If any edge incident to `i` is mapped to `kl`, then `i` must be mapped to `k` or
+`l`.
+"""
+function add_improved_topology_constraints_H_to_G!(model::Model, vars::Variables, G, H)
+    @constraint(
+        model,
+        [k in vertices(H), l in vertices(H), i in vertices(G); has_edge(H, k, l) && k < l],
+        sum(vars.y[i, :, k, l]; init=0) + sum(vars.y[:, i, k, l]; init=0) <=
+            vars.x[i, k] + vars.x[i, l]
+    )
+    return nothing
+end
+
+"""
+Add topology constraints for FORI. Uses the same implictions as used in
+[`add_improved_topology_constraints_G_to_H`](@ref) and
+[`add_improved_topology_constraints_H_to_G`](@ref), but since edges in `G` are oriented and
+each edge in `H` has two possible orientations, the implications are more precise.
+"""
+function add_oriented_topology_constraints!(model::Model, vars::OrientedVariables, G, H)
+    @constraint(
+        model,
+        [k in vertices(H), i in vertices(G), j in vertices(G); has_edge(G, i, j) && i < j],
+        sum(vars.z[i, j, k, :]; init=0) <= vars.x[i, k]
+    )
+    @constraint(
+        model,
+        [k in vertices(H), i in vertices(G), j in vertices(G); has_edge(G, i, j) && i < j],
+        sum(vars.z[i, j, :, k]; init=0) <= vars.x[j, k]
+    )
+    @constraint(
+        model,
+        [k in vertices(H), l in vertices(H), i in vertices(G); has_edge(H, k, l)],
+        sum(vars.z[i, :, k, l]; init=0) + sum(vars.z[:, i, l, k]; init=0) <= vars.x[i, k]
+    )
+    return nothing
+end
+
+"""
+    construct_formulation!(::Type{<:Formulation}, model, G, H;
+        c)
+
+Modify `model` to use a specific formulation to solve the graph edit distance problem. Each
+formulation has its own method implementing the required variables, constraints and
+objective.
+"""
+function construct_formulation! end
+
+function construct_formulation!(
+    ::Type{F1}, model, G, H, c::EditCosts=get_default_edit_costs(G, H)
+)
+    vars = create_model_vars_full!(model, G, H)
+
+    add_node_map_constraints!(model, vars, G, H)
+    add_edge_map_constraints!(model, vars, G, H)
+
+    add_simple_topology_constraints!(model, vars, G, H)
+
+    add_F1_objective!(model, c, vars, G, H)
+    return nothing
+end
+
+function construct_formulation!(
+    ::Type{F2minus}, model, G, H, c::EditCosts=get_default_edit_costs(G, H)
+)
+    vars = create_model_vars_reduced!(model, G, H)
+
+    add_node_map_constraints!(model, vars, G, H)
+    # Edgemap constraints are implied, so we can skip them.
+
+    add_simple_topology_constraints!(model, vars, G, H)
+
+    add_F2_objective!(model, c, vars, G, H)
+    return nothing
+end
+
+function construct_formulation!(
+    ::Type{F2}, model, G, H, c::EditCosts=get_default_edit_costs(G, H)
+)
+    vars = create_model_vars_reduced!(model, G, H)
+
+    add_node_map_constraints!(model, vars, G, H)
+    # Edgemap constraints are implied, so we can skip them.
+
+    add_improved_topology_constraints_G_to_H!(model, vars, G, H)
+
+    add_F2_objective!(model, c, vars, G, H)
+    return nothing
+end
+
+function construct_formulation!(
+    ::Type{F2plus}, model, G, H, c::EditCosts=get_default_edit_costs(G, H)
+)
+    vars = create_model_vars_reduced!(model, G, H)
+
+    add_node_map_constraints!(model, vars, G, H)
+    # Edgemap constraints are implied, so we can skip them.
+
+    add_improved_topology_constraints_G_to_H!(model, vars, G, H)
+    add_improved_topology_constraints_H_to_G!(model, vars, G, H)
+
+    add_F2_objective!(model, c, vars, G, H)
+    return nothing
+end
+
+function construct_formulation!(
+    ::Type{F1prime}, model, G, H, c::EditCosts=get_default_edit_costs(G, H)
+)
+    vars = create_model_vars_full!(model, G, H)
+
+    add_node_map_constraints!(model, vars, G, H)
+    add_edge_map_constraints!(model, vars, G, H)
+
+    add_improved_topology_constraints_G_to_H!(model, vars, G, H)
+
+    add_F1_objective!(model, c, vars, G, H)
+    return nothing
+end
+
+function construct_formulation!(
+    ::Type{F1plus}, model, G, H, c::EditCosts=get_default_edit_costs(G, H)
+)
+    vars = create_model_vars_full!(model, G, H)
+
+    add_node_map_constraints!(model, vars, G, H)
+    add_edge_map_constraints!(model, vars, G, H)
+
+    add_improved_topology_constraints_G_to_H!(model, vars, G, H)
+    add_improved_topology_constraints_H_to_G!(model, vars, G, H)
+
+    add_F1_objective!(model, c, vars, G, H)
+    return nothing
+end
+
+function construct_formulation!(
+    ::Type{FORI}, model, G, H, c::EditCosts=get_default_edit_costs(G, H)
+)
+    vars = create_model_vars_bidirectional!(model, G, H)
+    add_node_map_constraints!(model, vars, G, H)
+    # Edgemap constraints are implied, so we can skip them.
+
+    add_oriented_topology_constraints!(model, vars, G, H)
+
+    add_FORI_objective!(model, c, vars, G, H)
+    return nothing
+end
+
+"""
+    edit_distance(model, G, H;
+        c, formulation
+    )
+
+Modify a JuMP model to compute the graph edit distance between undirected graphs `G` and `H`
+given cost function `c` using formulation `formulation`. See [`edit_distance`](@ref) for
+more details.
+"""
+function edit_distance!(
+    model,
+    G::SimpleGraph,
+    H::SimpleGraph;
+    c::EditCosts=get_default_edit_costs(G, H),
+    formulation::Type{<:Formulation}=FORI,
+)
+    construct_formulation!(formulation, model, G, H, c)
+    return nothing
+end
+
+"""
+    edit_distance(G, H;
+        c, formulation, optimizer
+    )
+
+Compute the graph edit distance between undirected graphs `G` and `H` given edit costs `c`.
+
+# Arguments
+- `G::Graphs.SimpleGraph`
+- `H::Graphs.SimpleGraph`
+
+# Keywords
+- `c::EditCosts`: pairwise edit cost struct for the graphs. See
+    [`GraphsOptim.EditCosts`](@ref). Defaults to the canonical edit costs
+    [`GraphsOptim.get_default_edit_costs`](@ref).
+- `formulation::Type{<:Formulation}`: the ILP formulation to use. Defaults to the most 
+    powerful `FORI`.
+- `optimizer`: JuMP-compatible solver (default is `HiGHS.Optimizer`)
+
+# Returns
+- `Matrix{Int}`: the node map matrix, which encodes the optimal edit path
+"""
+function edit_distance(
+    G::SimpleGraph,
+    H::SimpleGraph;
+    c::EditCosts=get_default_edit_costs(G, H),
+    formulation::Type{<:Formulation}=FORI,
+    optimizer=HiGHS.Optimizer,
+)
+    model = Model(optimizer)
+    set_silent(model)
+    edit_distance!(model, G, H; c=c, formulation=formulation)
+    optimize!(model)
+    if termination_status(model) != OPTIMAL
+        error("Graph edit distance was not solved optimally.")
+    end
+    node_matching = convert(Matrix{Int}, round.(value.(model[:x])))
+    return (; objective_value=objective_value(model), node_matching)
+end
diff --git a/test/graph_edit_distance.jl b/test/graph_edit_distance.jl
new file mode 100644
index 0000000..f235292
--- /dev/null
+++ b/test/graph_edit_distance.jl
@@ -0,0 +1,177 @@
+using GraphsOptim, Graphs, JuMP, HiGHS
+using Test
+
+all_formulations = [F1, F1prime, F1plus, F2, F2minus, F2plus, FORI]
+
+@testset "Small Graph" begin
+    G = Graph(3)
+    add_edge!(G, 1, 2)
+
+    H = Graph(3)
+    add_edge!(H, 2, 3)
+    add_edge!(H, 1, 3)
+
+    for formulation in all_formulations
+        model = Model(HiGHS.Optimizer)
+        GraphsOptim.edit_distance!(model, G, H; formulation=formulation)
+        set_silent(model)
+        optimize!(model)
+        @test objective_value(model) == 1
+        # testing the normal entrypoint
+        @test GraphsOptim.edit_distance(G, H).objective_value == 1
+    end
+end
+
+@testset "Cost function validations" begin
+    G = Graph(3)
+    add_edge!(G, 1, 2)
+
+    H = Graph(3)
+    add_edge!(H, 2, 3)
+    add_edge!(H, 1, 3)
+
+    c = GraphsOptim.get_default_edit_costs(G, H)
+    @test isnothing(GraphsOptim.validate_cost_function(c, G, H))
+
+    broken_c = GraphsOptim.EditCosts(
+        zeros(Int, 1, 2), c.c_iε, c.c_εk, c.c_ijkl, c.c_ijε, c.c_εkl
+    )
+    @test_throws AssertionError GraphsOptim.validate_cost_function(broken_c, G, H)
+
+    Hprime = Graph(4)
+    @test_throws AssertionError GraphsOptim.validate_cost_function(c, G, Hprime)
+end
+
+@testset "Cost Functions Simple" begin
+    G = Graph(3)
+    add_edge!(G, 1, 2)
+
+    H = Graph(3)
+    add_edge!(H, 2, 3)
+    add_edge!(H, 1, 3)
+    c_origin = GraphsOptim.get_default_edit_costs(G, H)
+
+    # increating costs for node subsitutions 
+    c = GraphsOptim.EditCosts(
+        ones(Int, nv(G), nv(H)),
+        5 * c_origin.c_iε,
+        5 * c_origin.c_εk,
+        c_origin.c_ijkl,
+        c_origin.c_ijε,
+        c_origin.c_εkl,
+    )
+    for formulation in all_formulations
+        model = Model(HiGHS.Optimizer)
+        GraphsOptim.edit_distance!(model, G, H; c=c, formulation=formulation)
+        set_silent(model)
+        optimize!(model)
+        @test objective_value(model) == 4
+    end
+    # increating costs for edge subsitutions 
+    c = GraphsOptim.EditCosts(
+        c_origin.c_ik,
+        c_origin.c_iε,
+        c_origin.c_εk,
+        c_origin.c_ijkl,
+        100 * c_origin.c_ijε,
+        10 * c_origin.c_εkl,
+    )
+    for formulation in all_formulations
+        model = Model(HiGHS.Optimizer)
+        GraphsOptim.edit_distance!(model, G, H; c=c, formulation=formulation)
+        set_silent(model)
+        optimize!(model)
+        @test objective_value(model) == 10
+    end
+
+    # force a specific suboptimal node map with edit costs 
+    G = Graph(3)
+    add_edge!(G, 2, 3)
+
+    H = Graph(4)
+    add_edge!(H, 1, 3)
+    add_edge!(H, 3, 4)
+    c_origin = GraphsOptim.get_default_edit_costs(G, H)
+    custom_substitution_cost = 100 .+ c_origin.c_ik
+    custom_substitution_cost[1, 1] = 0
+    custom_substitution_cost[2, 2] = 0
+    custom_substitution_cost[3, 3] = 0
+    c = GraphsOptim.EditCosts(
+        custom_substitution_cost,
+        c_origin.c_iε,
+        c_origin.c_εk,
+        c_origin.c_ijkl,
+        c_origin.c_ijε,
+        c_origin.c_εkl,
+    )
+    for formulation in all_formulations
+        model = Model(HiGHS.Optimizer)
+        vars = GraphsOptim.edit_distance!(model, G, H; c=c, formulation=formulation)
+        set_silent(model)
+        optimize!(model)
+        @test objective_value(model) == 4
+        @test value(model[:x][1, 1]) == 1
+        @test value(model[:x][2, 2]) == 1
+        @test value(model[:x][3, 3]) == 1
+    end
+end
+
+@testset "Edgecases" begin
+    G = Graph(3)
+
+    H = Graph(3)
+    add_edge!(H, 1, 2)
+
+    # one graph empty
+    for formulation in all_formulations
+        model = Model(HiGHS.Optimizer)
+        GraphsOptim.edit_distance!(model, G, H; formulation=formulation)
+        set_silent(model)
+        optimize!(model)
+        @test objective_value(model) == 1
+
+        # reverse order
+        model = Model(HiGHS.Optimizer)
+        GraphsOptim.edit_distance!(model, H, G; formulation=formulation)
+        set_silent(model)
+        optimize!(model)
+        @test objective_value(model) == 1
+    end
+
+    # both graphs empty
+    H = Graph(5)
+    for formulation in all_formulations
+        model = Model(HiGHS.Optimizer)
+        GraphsOptim.edit_distance!(model, G, H; formulation=formulation)
+        set_silent(model)
+        optimize!(model)
+        @test objective_value(model) == 2
+    end
+
+    # no nodes in one graph
+    H = Graph(0)
+    for formulation in [F1, F1prime, F1plus]
+        model = Model(HiGHS.Optimizer)
+        GraphsOptim.edit_distance!(model, G, H; formulation=formulation)
+        set_silent(model)
+        optimize!(model)
+        @test objective_value(model) == 3
+    end
+    # for F2 type formulations, the formulation contains no variables
+    # we thus add a dummy variable to make the solver accept the model
+    for formulation in [F2, F2minus, F2plus, FORI]
+        model = Model(HiGHS.Optimizer)
+        GraphsOptim.edit_distance!(model, G, H; formulation=formulation)
+        @test num_variables(model) == 0
+        @variable(model, dummy)
+        set_silent(model)
+        optimize!(model)
+        @test objective_value(model) == 3
+    end
+end
+
+@testset "Invalid Input" begin
+    G = Graph(3)
+    H = DiGraph(4)
+    @test_throws MethodError GraphsOptim.edit_distance(G, H)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 6a46e91..3001d87 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -57,4 +57,8 @@ using Test
     @testset verbose = true "Shortest path" begin
         include("shortest_path.jl")
     end
+
+    @testset "Graph Edit Distance" begin
+        include("graph_edit_distance.jl")
+    end
 end;