From c1cd68462b803bfbd97ff9cb6ef2fc9015607954 Mon Sep 17 00:00:00 2001
From: Ken Alan Berkpinar <ken.berkpinar@gmail.com>
Date: Sun, 2 Aug 2026 15:54:24 +0200
Subject: [PATCH 1/9] feat: Add PATH-Algorithm feature for graph matching

---
 Project.toml              |   4 +-
 src/GraphsOptim.jl        |   5 +-
 src/pathGraphMatching.jl  | 444 ++++++++++++++++++++++++++++++++++++++
 test/pathGraphMatching.jl |  23 ++
 test/runtests.jl          |   4 +
 5 files changed, 478 insertions(+), 2 deletions(-)
 create mode 100644 src/pathGraphMatching.jl
 create mode 100644 test/pathGraphMatching.jl

diff --git a/Project.toml b/Project.toml
index 8e9a517..cef71db 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,10 +1,11 @@
 name = "GraphsOptim"
 uuid = "e79ef3ae-79c8-4b97-b165-63f338db30c2"
-authors = ["Guillaume Dalle, Aurora Rossi and contributors"]
 version = "0.1.0"
+authors = ["Guillaume Dalle, Aurora Rossi and contributors"]
 
 [deps]
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+FrankWolfe = "f55ce6ea-fdc5-4628-88c5-0087fe54bd30"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
@@ -17,6 +18,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Aqua = "0.6"
 Documenter = "0.27"
 FillArrays = "0.13, 1.15"
+FrankWolfe = "0.6.4"
 Graphs = "1.7"
 HiGHS = "1"
 JuMP = "1"
diff --git a/src/GraphsOptim.jl b/src/GraphsOptim.jl
index 7156d42..ca47bdd 100644
--- a/src/GraphsOptim.jl
+++ b/src/GraphsOptim.jl
@@ -15,10 +15,11 @@ 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 LinearAlgebra: norm, tr, dot
+using LinearAlgebra: norm, tr, dot, I
 using MathOptInterface: OPTIMAL
 using SparseArrays: sparse
 using OptimalTransport: sinkhorn
+using FrankWolfe
 
 export min_cost_flow
 export min_cost_assignment
@@ -28,6 +29,7 @@ export maximum_weight_independent_set
 export fractional_chromatic_number, fractional_clique_number
 export shortest_path
 export maximum_weight_clique
+export pathAlgorithm
 
 include("utils.jl")
 include("flow.jl")
@@ -38,5 +40,6 @@ include("fractional_coloring.jl")
 include("shortest_path.jl")
 include("maximum_clique.jl")
 include("independent_set.jl")
+include("pathGraphMatching.jl")
 
 end
diff --git a/src/pathGraphMatching.jl b/src/pathGraphMatching.jl
new file mode 100644
index 0000000..365f0b8
--- /dev/null
+++ b/src/pathGraphMatching.jl
@@ -0,0 +1,444 @@
+
+"""
+pathAlgorithm(G::Matrix{Float64}, H::Matrix{Float64}, ϵ_λ_f::Float64=0.1, ϵ_λ_p::Float64=0.1; kwargs...)
+
+Solve the graph matching problem (or Quadratic Assignment Problem / QAP) using a path-following algorithm.
+
+The algorithm tracks a convex combination (via parameter λ ∈ [0, 1]) between an easily solvable convex relaxation and a concave function that is as hard as the problem. If input matrices differ in size, the smaller matrix is automatically padded with zero rows and columns.
+
+# Arguments
+- `G::Matrix{Float64}`: Adjacency or cost matrix of the first graph (n × n).
+- `H::Matrix{Float64}`: Adjacency or cost matrix of the second graph (m × m).
+- `ϵ_λ_f::Float64=0.1`: Threshold for the maximum normalized change in function value during dynamic step-size control of λ.
+- `ϵ_λ_p::Float64=0.1`: Threshold for the normalized permutation matrix change (||P_{new} - P_{opt}|| / √{2n}) during step-size control.
+
+# Keywords
+- `dλ_min::Float64=1.0e-5`: Minimum step size for incrementing the path parameter λ.
+- `solveQAP::Bool=false`: If `true`, adjusts sign logic to solve a general QAP.
+- `return_log::Bool=false`: If `true`, returns a formatted summary log string (runtime, costs, and iterations).
+- `return_dataPoints::Bool=false`: If `true`, returns a `NamedTuple` containing trace histories for λ, f_0, f_1, and f_λ.
+- `verbose::Bool=false`: Enable detailed console output tracking path-following progress.
+- `verbose_FW::Bool=false`: Enable console output for individual Frank-Wolfe optimization steps.
+
+# Returns
+- `p_vec::Vector{Int}`: Resulting permutation vector indicating node assignments.
+- `log_string::Union{String, Nothing}`: Formatted summary string if `return_log=true`, otherwise `nothing`.
+- `dataPoints::Union{NamedTuple, Nothing}`: `NamedTuple` containing `λ_list`, `f0_list`, `f1_list`, and `fλ_list` if `return_dataPoints=true`, otherwise `nothing`.
+"""
+function pathAlgorithm(G::Matrix{Float64}, H::Matrix{Float64}, ϵ_λ_f::Float64=0.1, ϵ_λ_p::Float64=0.1;
+dλ_min::Float64=1.0e-5,
+solveQAP::Bool=false,
+return_log::Bool=false,
+return_dataPoints::Bool=false,
+verbose::Bool=false,
+verbose_FW::Bool=false
+)
+# extend the smaller matrix by zero rows and columns
+diffSize = size(G,1)-size(H,1)
+
+if diffSize > 0
+    # G is larger
+    verbose && println("G is larger than H by ", diffSize, " rows and columns. Adding zeros to H.")
+    H = cat(H,zeros(diffSize,diffSize); dims=(1,2))
+elseif diffSize < 0
+    # H is larger
+    verbose && println("H is larger than G by ", abs(diffSize), " rows and columns. Adding zeros to G.")
+    diffSize = abs(diffSize)
+    G = cat(G,zeros(diffSize,diffSize); dims=(1,2))
+end
+m_size = size(G,1)
+verbose && println("Size of G and H: ", m_size, " x ", m_size)
+
+t1 = time()
+
+# allocate fixed space for the gradient matrices so that they don't allocate new space in each calculation
+storage0 = Matrix{Float64}(undef, m_size, m_size)
+storage1 = Matrix{Float64}(undef, m_size, m_size)
+
+# Start with P as the identity matrix
+p_start = Matrix(1.0I, m_size, m_size)
+lmo = FrankWolfe.BirkhoffPolytopeLMO() #via Hungarian algorithm
+
+verbose && println("Starting path-following algorithm with λ = 0.0")
+
+# find initial minimum of F0 ( -F1 for QAP)
+# TODO use Newton instead of FrankWolfe for initialization
+verbose && println("Finding initial minimum of F0 with FrankWolfe")
+if !solveQAP
+    init_f   = FλForP(0.0, G, H)
+    init_∇!  = ∇FλForP!(storage0, storage1, 0.0, G, H)
+else
+    init_f   = FλForP_QAP(0.0, G, H)
+    init_∇!  = ∇FλForP_QAP!(storage0, storage1, 0.0, G, H)
+end
+
+p_opt, _ = FrankWolfe.frank_wolfe(
+init_f, init_∇!, lmo, p_start;
+epsilon = 1e-8,
+max_iteration = 10_000,
+verbose=verbose_FW
+)
+
+# change in λ is dynamically adjusted; starts at minimum
+dλ = dλ_min
+# begin with λ=0; iteratively increase up until 1
+λ = 0.0
+
+# redefine f0, f1 and fλ depending on whether the QAP should be solved or not, s.t. f0 is always convex and f1 is always concave.
+if !solveQAP
+    fλNormalizedFinal = fλNormalized
+    fλNormalizedFinal = fλ_QAP
+end
+
+count_iter = 0
+
+λ_list = Float64[]
+f0_list = Float64[]
+f1_list = Float64[]
+fλ_list = Float64[]
+if return_dataPoints
+    push!(λ_list, λ)
+    push!(f0_list, f0(p_opt,G,H))
+    push!(f1_list, f1(p_opt,G,H))
+    push!(fλ_list, fλ(p_opt,λ,G,H))
+end
+
+verbose && println("λ = ", λ)
+verbose && println()
+while(λ < 1.0)
+    count_iter += 1
+    # set first possible value for λ_new
+    λ_new = λ + dλ
+
+    # calculate local optimum w.r.t. initial λ_new
+    verbose && println("   dλ = ", dλ)
+    if !solveQAP
+        fλ_new_minimize = FλForP(λ_new, G, H)
+        ∇fλ_new_minimize = ∇FλForP!(storage0, storage1, λ_new, G, H)
+    else
+        fλ_new_minimize = FλForP_QAP(λ_new, G, H)
+        ∇fλ_new_minimize = ∇FλForP_QAP!(storage0, storage1, λ_new, G, H)
+    end
+    p_new, _ = frank_wolfe(
+        fλ_new_minimize, ∇fλ_new_minimize, lmo, p_opt; 
+        epsilon = 1e-8,
+        max_iteration = 10_000,
+        verbose=verbose_FW
+    )
+    p_change_normalized = norm(p_new - p_opt) / sqrt(2 * m_size)
+
+    p_last::Union{Nothing, Matrix{Float64}} = nothing
+
+    # update dλ until criterion is met
+    # TODO implemented new stopping criterion. Need to still find out ϵ_f and ϵ_p values from FrankWolfe implementation and calculate ϵ_λ_f and ϵ_λ_p with added input M.
+    # d_λ is doubled until one value is larger than it's threshold (or new λ is already 1)
+    while abs(fλNormalizedFinal(p_new,λ_new,G,H)-fλNormalizedFinal(p_opt,λ,G,H)) < ϵ_λ_f   &&   p_change_normalized < ϵ_λ_p   &&   λ_new < one(Float64)
+        dλ = 2*dλ
+        λ_new = min(λ + dλ, one(Float64))
+
+        verbose && println("   dλ = ", dλ)
+        if !solveQAP
+            fλ_new_minimize = FλForP(λ_new, G, H)
+            ∇fλ_new_minimize = ∇FλForP!(storage0, storage1, λ_new, G, H)
+        else
+            fλ_new_minimize = FλForP_QAP(λ_new, G, H)
+            ∇fλ_new_minimize = ∇FλForP_QAP!(storage0, storage1, λ_new, G, H)
+        end
+        p_last = p_new
+        p_new, _ = frank_wolfe(
+            fλ_new_minimize, ∇fλ_new_minimize, lmo, p_opt; 
+            epsilon = 1e-8,
+            max_iteration = 10_000,
+            verbose = verbose_FW
+        )
+        p_change_normalized = norm(p_new - p_opt) / sqrt(2 * m_size)
+    end
+    
+    # if the last while loop's condition is not met (anymore), dλ is one step too large and can be halved once directly
+    dλ = max(dλ/2,dλ_min)
+    λ_new = λ + dλ
+    verbose && println("   dλ = ", dλ)
+    if !isnothing(p_last)
+        p_new = p_last
+    else
+        if !solveQAP
+            fλ_new_minimize = FλForP(λ_new, G, H)
+            ∇fλ_new_minimize = ∇FλForP!(storage0, storage1, λ_new, G, H)
+        else
+            fλ_new_minimize = FλForP_QAP(λ_new, G, H)
+            ∇fλ_new_minimize = ∇FλForP_QAP!(storage0, storage1, λ_new, G, H)
+        end
+        p_new, _ = frank_wolfe(
+            fλ_new_minimize, ∇fλ_new_minimize, lmo, p_opt; 
+            epsilon = 1e-8,
+            max_iteration = 10_000,
+            verbose = verbose_FW
+        )
+    end
+    p_change_normalized = norm(p_new - p_opt) / sqrt(2 * m_size)
+
+    # d_λ is halved until both values are smaller than their thresholds (or dλ is already at minimum)
+    while (abs(fλNormalizedFinal(p_new,λ_new,G,H)-fλNormalizedFinal(p_opt,λ,G,H)) > ϵ_λ_f   ||   p_change_normalized > ϵ_λ_p)   &&   dλ > dλ_min
+        dλ = max(dλ/2,dλ_min)
+        λ_new = min(λ + dλ, one(Float64))
+        verbose && println("   dλ = ", dλ)
+
+        if !solveQAP
+            fλ_new_minimize = FλForP(λ_new, G, H)
+            ∇fλ_new_minimize = ∇FλForP!(storage0, storage1, λ_new, G, H)
+        else
+            fλ_new_minimize = FλForP_QAP(λ_new, G, H)
+            ∇fλ_new_minimize = ∇FλForP_QAP!(storage0, storage1, λ_new, G, H)
+        end
+        p_new, _ = frank_wolfe(
+            fλ_new_minimize, ∇fλ_new_minimize, lmo, p_opt; 
+            epsilon = 1e-8,
+            max_iteration = 10_000,
+            verbose = verbose_FW
+        )
+        p_change_normalized = norm(p_new - p_opt) / sqrt(2 * m_size)
+    end
+    λ = λ_new
+    verbose && println("λ = ", λ)
+    verbose && println()
+    # criterion is met, λ is set correctly and p_new contans the local optimum w.r.t. the new λ. Set p_opt to p_new for next iteration.
+
+    p_opt = p_new
+
+    if return_dataPoints
+        push!(λ_list, λ)
+        push!(f0_list, f0(p_opt,G,H))
+        push!(f1_list, f1(p_opt,G,H))
+        push!(fλ_list, fλ(p_opt,λ,G,H))
+    end
+
+    # stop immediately if FrankWolfe arrives at a Permutationmatrix as this is a feasible minimum
+    if isPerm(p_opt)
+        verbose && println("Found a Permutationmatrix as local optimum, stopping path-following algorithm")
+        verbose && println("P:")
+        break
+    end
+end
+p_vec = permMtV(p_opt)
+verbose && display(p_vec)
+
+elapsed_time = time() - t1
+verbose && println("Elapsed time: ", elapsed_time, " seconds")
+
+
+log_stream = IOBuffer()
+if return_log
+    function write_log(msg)
+        println(log_stream, msg) # Schreibt in den Buffer
+    end
+
+    write_log("="^60)
+    write_log("Results for Graph Matching/QAP")
+    write_log("="^60)
+    write_log("")
+    write_log("ϵ_λ_f: $(ϵ_λ_f)")
+    write_log("ϵ_λ_p: $(ϵ_λ_p)")
+    write_log("solveQAP: $(solveQAP)")
+    write_log("")
+    write_log("Runtime: $(elapsed_time) seconds")
+    write_log("λ Iterations: $(count_iter)")
+    write_log("")
+    write_log("Cost:")
+    if !solveQAP
+        write_log("F0: $(f0(p_opt, G, H))")
+        write_log("F1: $(f1(p_opt, G, H))")
+    else
+        write_log("$(qapVal(p_opt, G, H))")
+    end
+    write_log("")
+    write_log("-"^60)
+    write_log("Resulting Matrix P")
+    write_log("-"^60)
+    write_log(p_vec)
+end
+log_string = return_log ? String(take!(log_stream)) : nothing
+
+dataPoints = return_dataPoints ? (;
+        λ_list = λ_list,
+        f0_list = f0_list,
+        f1_list = f1_list,
+        fλ_list = fλ_list
+    ) : nothing
+
+return p_vec, log_string, dataPoints
+
+end
+
+# returns true if P contains only zeros and ones and false if not
+function isPerm(P)
+return all(x -> x == 0.0 || x == 1.0, P)
+end
+
+# returns the permutation vector of a permutation matrix P
+function permMtV(P)
+return [argmax(row) for row in eachrow(P)]
+end
+
+# returns the permutation matrix of a permutation vector P
+function permVtM(P)
+return Matrix{Float64}(I(length(P))[P, :])
+end
+
+# returns the squared frobenius norm of matrix A
+function sqd_frob(A)
+val = norm(A,2)
+return val^2
+end
+
+# returns the diagonal degree matrix of G
+# column by column is quicker to go through in julia
+function diagonal_degree(G)
+D = zeros(size(G))
+for j = 1:size(D,1)
+    sum = 0.0
+    for i = 1:size(D,1)
+        sum += G[i, j]
+    end
+    D[j, j] = sum
+end
+return D
+end
+
+# returns the matrix Δ as stated in the paper
+function Δ(G,H)
+D_G = diagonal_degree(G)
+D_H = diagonal_degree(H)
+
+Δ_G_H = zeros(size(G))
+for i = 1:size(G,1)
+    for j = 1:size(G,1)
+        Δ_G_H[i, j] = D_H[j, j] - D_G[i, i]
+    end
+end
+return Δ_G_H .^ 2
+end
+
+# returns the laplacian matrix of G
+function laplacian(G)
+return diagonal_degree(G) .- G
+end
+
+# convex function F0
+# algorithm uses only normalized version. This is just for plotting and displaying the correct data.
+function f0(P,G,H)
+return sqd_frob(G*P .- P*H)
+end
+
+# F0 normalized for values between 0 and 1.
+function f0Normalized(P,G,H)
+value = f0(P,G,H)
+return value ./ (sqd_frob(G) + sqd_frob(H))
+end
+
+# gradient of F0 but normalized for values between 0 and 1
+# save solution value in variable "storage" for space economy
+function ∇f0Normalized!(storage, P, G, H)
+    value = 2.0 .* ((G^2) * P .- 2.0 .* G * P * H .+ P * (H^2))
+    storage .= value ./ (sqd_frob(G) + sqd_frob(H))
+end
+
+# concave function F1.
+# algorithm uses only normalized version. This is just for plotting and displaying the correct data.
+function f1(P, G, H)
+constantTerm = tr(laplacian(G)^2)+tr(laplacian(H)^2)
+return .- tr(Δ(G,H)'*P) .- 2.0 .* (vec(P)' * vec(laplacian(G) * P * laplacian(H))) + constantTerm
+end
+
+# F1 normalized for values between 0 and 1.
+function f1Normalized(P, G, H)
+value = f1(P, G, H)
+return value ./ (sqd_frob(G) + sqd_frob(H))
+end
+
+# gradient of F1 abut normalized for values between 0 and 1
+# save solution value in variable "storage" for space economy
+function ∇f1Normalized!(storage, P, G, H)
+    # the PATH-Algorithm paper has 2.0 in front of the second term, but 4.0 should be correct.
+value = .- Δ(G,H)' .- 4.0 .* laplacian(G) * P * laplacian(H)
+storage .= value ./ (sqd_frob(G) + sqd_frob(H))
+end
+
+# Fλ is convex combination of F0 and F1.
+# algorithm uses only normalized version. This is just for plotting and displaying the correct data.
+function fλ(P, λ, G, H)
+return (1-λ) * f0(P, G, H)  +  λ * f1(P, G, H)
+end
+
+# Fλ normalized for values between 0 and 1.
+function fλNormalized(P, λ, G, H)
+return (1-λ) * f0Normalized(P, G, H)  +  λ * f1Normalized(P, G, H)
+end
+
+struct FλForP
+    λ::Float64
+    G::Matrix{Float64}
+    H::Matrix{Float64}
+end
+# for the FW-algorithm we need a function that takes only P as input.
+function(fλ_struct::FλForP)(P) 
+    return fλNormalized(P, fλ_struct.λ, fλ_struct.G, fλ_struct.H)
+end
+
+# function flipped for maximization of the initial function and thus solving QAP
+function fλ_QAP(P, λ, G, H)
+return (1-λ) * (-f1Normalized(P, G, H))  +  λ * (-f0Normalized(P, G, H))
+end
+
+struct FλForP_QAP
+    λ::Float64
+    G::Matrix{Float64}
+    H::Matrix{Float64}
+end
+# for the FW-algorithm we need a function that takes only P as input.
+function(fλ_struct::FλForP_QAP)(P) 
+    return fλ_QAP(P, fλ_struct.λ, fλ_struct.G, fλ_struct.H)
+end
+
+# gradient of FλNormalized
+# save solution value in variable "storage" for space economy
+function ∇fλ!(storageλ, storage0, storage1, P, λ, G, H)
+∇f0Normalized!(storage0, P, G, H)
+∇f1Normalized!(storage1, P, G, H)
+storageλ .= (1.0-λ) .* storage0 .+ λ .* storage1
+end
+struct ∇FλForP!
+    storage0::Matrix{Float64}
+    storage1::Matrix{Float64}
+    λ::Float64
+    G::Matrix{Float64}
+    H::Matrix{Float64}
+end
+# for the FW-algorithm we need a function that takes only P as input.
+function(∇fλ_struct::∇FλForP!)(storageλ, P)
+    ∇fλ!(storageλ, ∇fλ_struct.storage0, ∇fλ_struct.storage1, P, ∇fλ_struct.λ, ∇fλ_struct.G, ∇fλ_struct.H)
+end
+
+# gradient flipped for maximization and solving QAP
+# save solution value in variable "storage" for space economy
+function ∇fλ_QAP!(storageλ, storage0, storage1, P, λ, G, H)
+∇f0Normalized!(storage0, P, G, H)
+∇f1Normalized!(storage1, P, G, H)
+storageλ .= (1.0-λ) .* (-storage1) .+ λ .* (-storage0)
+end
+
+struct ∇FλForP_QAP!
+    storage0::Matrix{Float64}
+    storage1::Matrix{Float64}
+    λ::Float64
+    G::Matrix{Float64}
+    H::Matrix{Float64}
+end
+# for the FW-algorithm we need a function that takes only P as input.
+function(∇fλ_struct::∇FλForP_QAP!)(storageλ, P)
+    ∇fλ_QAP!(storageλ, ∇fλ_struct.storage0, ∇fλ_struct.storage1, P, ∇fλ_struct.λ, ∇fλ_struct.G, ∇fλ_struct.H)
+end
+
+# returns the value of the QAP objective function for a given permutation matrix P and adjacency matrices G and H
+function qapVal(P,G,H)
+return tr(G*P*H'*P')
+end
\ No newline at end of file
diff --git a/test/pathGraphMatching.jl b/test/pathGraphMatching.jl
new file mode 100644
index 0000000..9064432
--- /dev/null
+++ b/test/pathGraphMatching.jl
@@ -0,0 +1,23 @@
+using GraphsOptim
+using LinearAlgebra
+using Test
+
+P, _, _ = GraphsOptim.pathAlgorithm([
+    1.0 2.0;
+    3.0 4.0
+    ],[
+    1.0 2.0;
+    3.0 4.0
+    ],
+    0.1, 0.1)
+@test P == [1, 2]
+
+P, _, _ = GraphsOptim.pathAlgorithm([
+    1.0 2.0;
+    3.0 4.0
+    ],[
+    4.0 3.0;
+    2.0 1.0
+    ],
+    0.1, 0.1)
+@test P == [2, 1]
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 6a46e91..88eb3a7 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 verbose = true "Path graph matching" begin
+        include("pathGraphMatching.jl")
+    end
 end;

From 91f9b912b510d62c8da8df7d0e9d0a64dff7f227 Mon Sep 17 00:00:00 2001
From: Ken Alan Berkpinar <ken.berkpinar@gmail.com>
Date: Mon, 3 Aug 2026 12:30:03 +0200
Subject: [PATCH 2/9] Add docstring to docs

---
 docs/src/algorithms.md | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/docs/src/algorithms.md b/docs/src/algorithms.md
index 7131a94..5b0ca7e 100644
--- a/docs/src/algorithms.md
+++ b/docs/src/algorithms.md
@@ -116,3 +116,9 @@ GraphsOptim.is_permutation_matrix
 GraphsOptim.flat_doubly_stochastic
 GraphsOptim.indvec
 ```
+
+## Path Graph Matching
+
+```@docs
+GraphsOptim.pathAlgorithm
+```

From 53278b453e501dedf68d73c4348fed56453b08c8 Mon Sep 17 00:00:00 2001
From: Ken Alan Berkpinar <ken.berkpinar@gmail.com>
Date: Mon, 3 Aug 2026 13:11:15 +0200
Subject: [PATCH 3/9] Format Code

---
 src/independent_set.jl    |   2 +-
 src/maximum_clique.jl     |   2 +-
 src/pathGraphMatching.jl  | 554 +++++++++++++++++++++-----------------
 test/pathGraphMatching.jl |  38 +--
 4 files changed, 325 insertions(+), 271 deletions(-)

diff --git a/src/independent_set.jl b/src/independent_set.jl
index 8f7762a..a45587c 100644
--- a/src/independent_set.jl
+++ b/src/independent_set.jl
@@ -16,7 +16,7 @@ function maximum_weight_independent_set!(
     model[Symbol(var_name)] = f
     @constraint(
         model,
-        covering_constraint[i=1:nv(g), j=1:nv(g); i ≠ j && has_edge(g, i, j)],
+        covering_constraint[i = 1:nv(g), j = 1:nv(g); i ≠ j && has_edge(g, i, j)],
         f[i] + f[j] <= 1,
     )
     obj = objective_function(model)
diff --git a/src/maximum_clique.jl b/src/maximum_clique.jl
index e0aad86..8568afc 100644
--- a/src/maximum_clique.jl
+++ b/src/maximum_clique.jl
@@ -16,7 +16,7 @@ function maximum_weight_clique!(
     model[Symbol(var_name)] = f
     @constraint(
         model,
-        packing_constraint[i=1:nv(g), j=1:nv(g); i ≠ j && !has_edge(g, i, j)],
+        packing_constraint[i = 1:nv(g), j = 1:nv(g); i ≠ j && !has_edge(g, i, j)],
         f[i] + f[j] <= 1,
     )
     obj = objective_function(model)
diff --git a/src/pathGraphMatching.jl b/src/pathGraphMatching.jl
index 365f0b8..7f706b1 100644
--- a/src/pathGraphMatching.jl
+++ b/src/pathGraphMatching.jl
@@ -25,117 +25,105 @@ The algorithm tracks a convex combination (via parameter λ ∈ [0, 1]) between
 - `log_string::Union{String, Nothing}`: Formatted summary string if `return_log=true`, otherwise `nothing`.
 - `dataPoints::Union{NamedTuple, Nothing}`: `NamedTuple` containing `λ_list`, `f0_list`, `f1_list`, and `fλ_list` if `return_dataPoints=true`, otherwise `nothing`.
 """
-function pathAlgorithm(G::Matrix{Float64}, H::Matrix{Float64}, ϵ_λ_f::Float64=0.1, ϵ_λ_p::Float64=0.1;
-dλ_min::Float64=1.0e-5,
-solveQAP::Bool=false,
-return_log::Bool=false,
-return_dataPoints::Bool=false,
-verbose::Bool=false,
-verbose_FW::Bool=false
+function pathAlgorithm(
+    G::Matrix{Float64},
+    H::Matrix{Float64},
+    ϵ_λ_f::Float64=0.1,
+    ϵ_λ_p::Float64=0.1;
+    dλ_min::Float64=1.0e-5,
+    solveQAP::Bool=false,
+    return_log::Bool=false,
+    return_dataPoints::Bool=false,
+    verbose::Bool=false,
+    verbose_FW::Bool=false,
 )
-# extend the smaller matrix by zero rows and columns
-diffSize = size(G,1)-size(H,1)
-
-if diffSize > 0
-    # G is larger
-    verbose && println("G is larger than H by ", diffSize, " rows and columns. Adding zeros to H.")
-    H = cat(H,zeros(diffSize,diffSize); dims=(1,2))
-elseif diffSize < 0
-    # H is larger
-    verbose && println("H is larger than G by ", abs(diffSize), " rows and columns. Adding zeros to G.")
-    diffSize = abs(diffSize)
-    G = cat(G,zeros(diffSize,diffSize); dims=(1,2))
-end
-m_size = size(G,1)
-verbose && println("Size of G and H: ", m_size, " x ", m_size)
-
-t1 = time()
-
-# allocate fixed space for the gradient matrices so that they don't allocate new space in each calculation
-storage0 = Matrix{Float64}(undef, m_size, m_size)
-storage1 = Matrix{Float64}(undef, m_size, m_size)
-
-# Start with P as the identity matrix
-p_start = Matrix(1.0I, m_size, m_size)
-lmo = FrankWolfe.BirkhoffPolytopeLMO() #via Hungarian algorithm
-
-verbose && println("Starting path-following algorithm with λ = 0.0")
-
-# find initial minimum of F0 ( -F1 for QAP)
-# TODO use Newton instead of FrankWolfe for initialization
-verbose && println("Finding initial minimum of F0 with FrankWolfe")
-if !solveQAP
-    init_f   = FλForP(0.0, G, H)
-    init_∇!  = ∇FλForP!(storage0, storage1, 0.0, G, H)
-else
-    init_f   = FλForP_QAP(0.0, G, H)
-    init_∇!  = ∇FλForP_QAP!(storage0, storage1, 0.0, G, H)
-end
+    # extend the smaller matrix by zero rows and columns
+    diffSize = size(G, 1)-size(H, 1)
 
-p_opt, _ = FrankWolfe.frank_wolfe(
-init_f, init_∇!, lmo, p_start;
-epsilon = 1e-8,
-max_iteration = 10_000,
-verbose=verbose_FW
-)
+    if diffSize > 0
+        # G is larger
+        verbose && println(
+            "G is larger than H by ", diffSize, " rows and columns. Adding zeros to H."
+        )
+        H = cat(H, zeros(diffSize, diffSize); dims=(1, 2))
+    elseif diffSize < 0
+        # H is larger
+        verbose && println(
+            "H is larger than G by ",
+            abs(diffSize),
+            " rows and columns. Adding zeros to G.",
+        )
+        diffSize = abs(diffSize)
+        G = cat(G, zeros(diffSize, diffSize); dims=(1, 2))
+    end
+    m_size = size(G, 1)
+    verbose && println("Size of G and H: ", m_size, " x ", m_size)
 
-# change in λ is dynamically adjusted; starts at minimum
-dλ = dλ_min
-# begin with λ=0; iteratively increase up until 1
-λ = 0.0
+    t1 = time()
 
-# redefine f0, f1 and fλ depending on whether the QAP should be solved or not, s.t. f0 is always convex and f1 is always concave.
-if !solveQAP
-    fλNormalizedFinal = fλNormalized
-    fλNormalizedFinal = fλ_QAP
-end
+    # allocate fixed space for the gradient matrices so that they don't allocate new space in each calculation
+    storage0 = Matrix{Float64}(undef, m_size, m_size)
+    storage1 = Matrix{Float64}(undef, m_size, m_size)
 
-count_iter = 0
-
-λ_list = Float64[]
-f0_list = Float64[]
-f1_list = Float64[]
-fλ_list = Float64[]
-if return_dataPoints
-    push!(λ_list, λ)
-    push!(f0_list, f0(p_opt,G,H))
-    push!(f1_list, f1(p_opt,G,H))
-    push!(fλ_list, fλ(p_opt,λ,G,H))
-end
+    # Start with P as the identity matrix
+    p_start = Matrix(1.0I, m_size, m_size)
+    lmo = FrankWolfe.BirkhoffPolytopeLMO() #via Hungarian algorithm
 
-verbose && println("λ = ", λ)
-verbose && println()
-while(λ < 1.0)
-    count_iter += 1
-    # set first possible value for λ_new
-    λ_new = λ + dλ
+    verbose && println("Starting path-following algorithm with λ = 0.0")
 
-    # calculate local optimum w.r.t. initial λ_new
-    verbose && println("   dλ = ", dλ)
+    # find initial minimum of F0 ( -F1 for QAP)
+    # TODO use Newton instead of FrankWolfe for initialization
+    verbose && println("Finding initial minimum of F0 with FrankWolfe")
     if !solveQAP
-        fλ_new_minimize = FλForP(λ_new, G, H)
-        ∇fλ_new_minimize = ∇FλForP!(storage0, storage1, λ_new, G, H)
+        init_f = FλForP(0.0, G, H)
+        init_∇! = ∇FλForP!(storage0, storage1, 0.0, G, H)
     else
-        fλ_new_minimize = FλForP_QAP(λ_new, G, H)
-        ∇fλ_new_minimize = ∇FλForP_QAP!(storage0, storage1, λ_new, G, H)
+        init_f = FλForP_QAP(0.0, G, H)
+        init_∇! = ∇FλForP_QAP!(storage0, storage1, 0.0, G, H)
     end
-    p_new, _ = frank_wolfe(
-        fλ_new_minimize, ∇fλ_new_minimize, lmo, p_opt; 
-        epsilon = 1e-8,
-        max_iteration = 10_000,
-        verbose=verbose_FW
+
+    p_opt, _ = FrankWolfe.frank_wolfe(
+        init_f,
+        init_∇!,
+        lmo,
+        p_start;
+        epsilon=1e-8,
+        max_iteration=10_000,
+        verbose=verbose_FW,
     )
-    p_change_normalized = norm(p_new - p_opt) / sqrt(2 * m_size)
 
-    p_last::Union{Nothing, Matrix{Float64}} = nothing
+    # change in λ is dynamically adjusted; starts at minimum
+    dλ = dλ_min
+    # begin with λ=0; iteratively increase up until 1
+    λ = 0.0
 
-    # update dλ until criterion is met
-    # TODO implemented new stopping criterion. Need to still find out ϵ_f and ϵ_p values from FrankWolfe implementation and calculate ϵ_λ_f and ϵ_λ_p with added input M.
-    # d_λ is doubled until one value is larger than it's threshold (or new λ is already 1)
-    while abs(fλNormalizedFinal(p_new,λ_new,G,H)-fλNormalizedFinal(p_opt,λ,G,H)) < ϵ_λ_f   &&   p_change_normalized < ϵ_λ_p   &&   λ_new < one(Float64)
-        dλ = 2*dλ
-        λ_new = min(λ + dλ, one(Float64))
+    # redefine f0, f1 and fλ depending on whether the QAP should be solved or not, s.t. f0 is always convex and f1 is always concave.
+    if !solveQAP
+        fλNormalizedFinal = fλNormalized
+        fλNormalizedFinal = fλ_QAP
+    end
+
+    count_iter = 0
+
+    λ_list = Float64[]
+    f0_list = Float64[]
+    f1_list = Float64[]
+    fλ_list = Float64[]
+    if return_dataPoints
+        push!(λ_list, λ)
+        push!(f0_list, f0(p_opt, G, H))
+        push!(f1_list, f1(p_opt, G, H))
+        push!(fλ_list, fλ(p_opt, λ, G, H))
+    end
 
+    verbose && println("λ = ", λ)
+    verbose && println()
+    while (λ < 1.0)
+        count_iter += 1
+        # set first possible value for λ_new
+        λ_new = λ + dλ
+
+        # calculate local optimum w.r.t. initial λ_new
         verbose && println("   dλ = ", dλ)
         if !solveQAP
             fλ_new_minimize = FλForP(λ_new, G, H)
@@ -144,234 +132,276 @@ while(λ < 1.0)
             fλ_new_minimize = FλForP_QAP(λ_new, G, H)
             ∇fλ_new_minimize = ∇FλForP_QAP!(storage0, storage1, λ_new, G, H)
         end
-        p_last = p_new
         p_new, _ = frank_wolfe(
-            fλ_new_minimize, ∇fλ_new_minimize, lmo, p_opt; 
-            epsilon = 1e-8,
-            max_iteration = 10_000,
-            verbose = verbose_FW
+            fλ_new_minimize,
+            ∇fλ_new_minimize,
+            lmo,
+            p_opt;
+            epsilon=1e-8,
+            max_iteration=10_000,
+            verbose=verbose_FW,
         )
         p_change_normalized = norm(p_new - p_opt) / sqrt(2 * m_size)
-    end
-    
-    # if the last while loop's condition is not met (anymore), dλ is one step too large and can be halved once directly
-    dλ = max(dλ/2,dλ_min)
-    λ_new = λ + dλ
-    verbose && println("   dλ = ", dλ)
-    if !isnothing(p_last)
-        p_new = p_last
-    else
-        if !solveQAP
-            fλ_new_minimize = FλForP(λ_new, G, H)
-            ∇fλ_new_minimize = ∇FλForP!(storage0, storage1, λ_new, G, H)
-        else
-            fλ_new_minimize = FλForP_QAP(λ_new, G, H)
-            ∇fλ_new_minimize = ∇FλForP_QAP!(storage0, storage1, λ_new, G, H)
+
+        p_last::Union{Nothing,Matrix{Float64}} = nothing
+
+        # update dλ until criterion is met
+        # TODO implemented new stopping criterion. Need to still find out ϵ_f and ϵ_p values from FrankWolfe implementation and calculate ϵ_λ_f and ϵ_λ_p with added input M.
+        # d_λ is doubled until one value is larger than it's threshold (or new λ is already 1)
+        while abs(fλNormalizedFinal(p_new, λ_new, G, H)-fλNormalizedFinal(p_opt, λ, G, H)) <
+              ϵ_λ_f &&
+              p_change_normalized < ϵ_λ_p &&
+              λ_new < one(Float64)
+            dλ = 2*dλ
+            λ_new = min(λ + dλ, one(Float64))
+
+            verbose && println("   dλ = ", dλ)
+            if !solveQAP
+                fλ_new_minimize = FλForP(λ_new, G, H)
+                ∇fλ_new_minimize = ∇FλForP!(storage0, storage1, λ_new, G, H)
+            else
+                fλ_new_minimize = FλForP_QAP(λ_new, G, H)
+                ∇fλ_new_minimize = ∇FλForP_QAP!(storage0, storage1, λ_new, G, H)
+            end
+            p_last = p_new
+            p_new, _ = frank_wolfe(
+                fλ_new_minimize,
+                ∇fλ_new_minimize,
+                lmo,
+                p_opt;
+                epsilon=1e-8,
+                max_iteration=10_000,
+                verbose=verbose_FW,
+            )
+            p_change_normalized = norm(p_new - p_opt) / sqrt(2 * m_size)
         end
-        p_new, _ = frank_wolfe(
-            fλ_new_minimize, ∇fλ_new_minimize, lmo, p_opt; 
-            epsilon = 1e-8,
-            max_iteration = 10_000,
-            verbose = verbose_FW
-        )
-    end
-    p_change_normalized = norm(p_new - p_opt) / sqrt(2 * m_size)
 
-    # d_λ is halved until both values are smaller than their thresholds (or dλ is already at minimum)
-    while (abs(fλNormalizedFinal(p_new,λ_new,G,H)-fλNormalizedFinal(p_opt,λ,G,H)) > ϵ_λ_f   ||   p_change_normalized > ϵ_λ_p)   &&   dλ > dλ_min
-        dλ = max(dλ/2,dλ_min)
-        λ_new = min(λ + dλ, one(Float64))
+        # if the last while loop's condition is not met (anymore), dλ is one step too large and can be halved once directly
+        dλ = max(dλ/2, dλ_min)
+        λ_new = λ + dλ
         verbose && println("   dλ = ", dλ)
-
-        if !solveQAP
-            fλ_new_minimize = FλForP(λ_new, G, H)
-            ∇fλ_new_minimize = ∇FλForP!(storage0, storage1, λ_new, G, H)
+        if !isnothing(p_last)
+            p_new = p_last
         else
-            fλ_new_minimize = FλForP_QAP(λ_new, G, H)
-            ∇fλ_new_minimize = ∇FλForP_QAP!(storage0, storage1, λ_new, G, H)
+            if !solveQAP
+                fλ_new_minimize = FλForP(λ_new, G, H)
+                ∇fλ_new_minimize = ∇FλForP!(storage0, storage1, λ_new, G, H)
+            else
+                fλ_new_minimize = FλForP_QAP(λ_new, G, H)
+                ∇fλ_new_minimize = ∇FλForP_QAP!(storage0, storage1, λ_new, G, H)
+            end
+            p_new, _ = frank_wolfe(
+                fλ_new_minimize,
+                ∇fλ_new_minimize,
+                lmo,
+                p_opt;
+                epsilon=1e-8,
+                max_iteration=10_000,
+                verbose=verbose_FW,
+            )
         end
-        p_new, _ = frank_wolfe(
-            fλ_new_minimize, ∇fλ_new_minimize, lmo, p_opt; 
-            epsilon = 1e-8,
-            max_iteration = 10_000,
-            verbose = verbose_FW
-        )
         p_change_normalized = norm(p_new - p_opt) / sqrt(2 * m_size)
-    end
-    λ = λ_new
-    verbose && println("λ = ", λ)
-    verbose && println()
-    # criterion is met, λ is set correctly and p_new contans the local optimum w.r.t. the new λ. Set p_opt to p_new for next iteration.
-
-    p_opt = p_new
 
-    if return_dataPoints
-        push!(λ_list, λ)
-        push!(f0_list, f0(p_opt,G,H))
-        push!(f1_list, f1(p_opt,G,H))
-        push!(fλ_list, fλ(p_opt,λ,G,H))
-    end
+        # d_λ is halved until both values are smaller than their thresholds (or dλ is already at minimum)
+        while (
+            abs(fλNormalizedFinal(p_new, λ_new, G, H)-fλNormalizedFinal(p_opt, λ, G, H)) >
+            ϵ_λ_f || p_change_normalized > ϵ_λ_p
+        ) && dλ > dλ_min
+            dλ = max(dλ/2, dλ_min)
+            λ_new = min(λ + dλ, one(Float64))
+            verbose && println("   dλ = ", dλ)
+
+            if !solveQAP
+                fλ_new_minimize = FλForP(λ_new, G, H)
+                ∇fλ_new_minimize = ∇FλForP!(storage0, storage1, λ_new, G, H)
+            else
+                fλ_new_minimize = FλForP_QAP(λ_new, G, H)
+                ∇fλ_new_minimize = ∇FλForP_QAP!(storage0, storage1, λ_new, G, H)
+            end
+            p_new, _ = frank_wolfe(
+                fλ_new_minimize,
+                ∇fλ_new_minimize,
+                lmo,
+                p_opt;
+                epsilon=1e-8,
+                max_iteration=10_000,
+                verbose=verbose_FW,
+            )
+            p_change_normalized = norm(p_new - p_opt) / sqrt(2 * m_size)
+        end
+        λ = λ_new
+        verbose && println("λ = ", λ)
+        verbose && println()
+        # criterion is met, λ is set correctly and p_new contans the local optimum w.r.t. the new λ. Set p_opt to p_new for next iteration.
+
+        p_opt = p_new
+
+        if return_dataPoints
+            push!(λ_list, λ)
+            push!(f0_list, f0(p_opt, G, H))
+            push!(f1_list, f1(p_opt, G, H))
+            push!(fλ_list, fλ(p_opt, λ, G, H))
+        end
 
-    # stop immediately if FrankWolfe arrives at a Permutationmatrix as this is a feasible minimum
-    if isPerm(p_opt)
-        verbose && println("Found a Permutationmatrix as local optimum, stopping path-following algorithm")
-        verbose && println("P:")
-        break
+        # stop immediately if FrankWolfe arrives at a Permutationmatrix as this is a feasible minimum
+        if isPerm(p_opt)
+            verbose && println(
+                "Found a Permutationmatrix as local optimum, stopping path-following algorithm",
+            )
+            verbose && println("P:")
+            break
+        end
     end
-end
-p_vec = permMtV(p_opt)
-verbose && display(p_vec)
+    p_vec = permMtV(p_opt)
+    verbose && display(p_vec)
 
-elapsed_time = time() - t1
-verbose && println("Elapsed time: ", elapsed_time, " seconds")
+    elapsed_time = time() - t1
+    verbose && println("Elapsed time: ", elapsed_time, " seconds")
 
+    log_stream = IOBuffer()
+    if return_log
+        function write_log(msg)
+            return println(log_stream, msg) # Schreibt in den Buffer
+        end
 
-log_stream = IOBuffer()
-if return_log
-    function write_log(msg)
-        println(log_stream, msg) # Schreibt in den Buffer
+        write_log("="^60)
+        write_log("Results for Graph Matching/QAP")
+        write_log("="^60)
+        write_log("")
+        write_log("ϵ_λ_f: $(ϵ_λ_f)")
+        write_log("ϵ_λ_p: $(ϵ_λ_p)")
+        write_log("solveQAP: $(solveQAP)")
+        write_log("")
+        write_log("Runtime: $(elapsed_time) seconds")
+        write_log("λ Iterations: $(count_iter)")
+        write_log("")
+        write_log("Cost:")
+        if !solveQAP
+            write_log("F0: $(f0(p_opt, G, H))")
+            write_log("F1: $(f1(p_opt, G, H))")
+        else
+            write_log("$(qapVal(p_opt, G, H))")
+        end
+        write_log("")
+        write_log("-"^60)
+        write_log("Resulting Matrix P")
+        write_log("-"^60)
+        write_log(p_vec)
     end
+    log_string = return_log ? String(take!(log_stream)) : nothing
 
-    write_log("="^60)
-    write_log("Results for Graph Matching/QAP")
-    write_log("="^60)
-    write_log("")
-    write_log("ϵ_λ_f: $(ϵ_λ_f)")
-    write_log("ϵ_λ_p: $(ϵ_λ_p)")
-    write_log("solveQAP: $(solveQAP)")
-    write_log("")
-    write_log("Runtime: $(elapsed_time) seconds")
-    write_log("λ Iterations: $(count_iter)")
-    write_log("")
-    write_log("Cost:")
-    if !solveQAP
-        write_log("F0: $(f0(p_opt, G, H))")
-        write_log("F1: $(f1(p_opt, G, H))")
+    dataPoints = if return_dataPoints
+        (; λ_list=λ_list, f0_list=f0_list, f1_list=f1_list, fλ_list=fλ_list)
     else
-        write_log("$(qapVal(p_opt, G, H))")
+        nothing
     end
-    write_log("")
-    write_log("-"^60)
-    write_log("Resulting Matrix P")
-    write_log("-"^60)
-    write_log(p_vec)
-end
-log_string = return_log ? String(take!(log_stream)) : nothing
-
-dataPoints = return_dataPoints ? (;
-        λ_list = λ_list,
-        f0_list = f0_list,
-        f1_list = f1_list,
-        fλ_list = fλ_list
-    ) : nothing
-
-return p_vec, log_string, dataPoints
 
+    return p_vec, log_string, dataPoints
 end
 
 # returns true if P contains only zeros and ones and false if not
 function isPerm(P)
-return all(x -> x == 0.0 || x == 1.0, P)
+    return all(x -> x == 0.0 || x == 1.0, P)
 end
 
 # returns the permutation vector of a permutation matrix P
 function permMtV(P)
-return [argmax(row) for row in eachrow(P)]
+    return [argmax(row) for row in eachrow(P)]
 end
 
 # returns the permutation matrix of a permutation vector P
 function permVtM(P)
-return Matrix{Float64}(I(length(P))[P, :])
+    return Matrix{Float64}(I(length(P))[P, :])
 end
 
 # returns the squared frobenius norm of matrix A
 function sqd_frob(A)
-val = norm(A,2)
-return val^2
+    val = norm(A, 2)
+    return val^2
 end
 
 # returns the diagonal degree matrix of G
 # column by column is quicker to go through in julia
 function diagonal_degree(G)
-D = zeros(size(G))
-for j = 1:size(D,1)
-    sum = 0.0
-    for i = 1:size(D,1)
-        sum += G[i, j]
+    D = zeros(size(G))
+    for j in 1:size(D, 1)
+        sum = 0.0
+        for i in 1:size(D, 1)
+            sum += G[i, j]
+        end
+        D[j, j] = sum
     end
-    D[j, j] = sum
-end
-return D
+    return D
 end
 
 # returns the matrix Δ as stated in the paper
-function Δ(G,H)
-D_G = diagonal_degree(G)
-D_H = diagonal_degree(H)
-
-Δ_G_H = zeros(size(G))
-for i = 1:size(G,1)
-    for j = 1:size(G,1)
-        Δ_G_H[i, j] = D_H[j, j] - D_G[i, i]
+function Δ(G, H)
+    D_G = diagonal_degree(G)
+    D_H = diagonal_degree(H)
+
+    Δ_G_H = zeros(size(G))
+    for i in 1:size(G, 1)
+        for j in 1:size(G, 1)
+            Δ_G_H[i, j] = D_H[j, j] - D_G[i, i]
+        end
     end
-end
-return Δ_G_H .^ 2
+    return Δ_G_H .^ 2
 end
 
 # returns the laplacian matrix of G
 function laplacian(G)
-return diagonal_degree(G) .- G
+    return diagonal_degree(G) .- G
 end
 
 # convex function F0
 # algorithm uses only normalized version. This is just for plotting and displaying the correct data.
-function f0(P,G,H)
-return sqd_frob(G*P .- P*H)
+function f0(P, G, H)
+    return sqd_frob(G*P .- P*H)
 end
 
 # F0 normalized for values between 0 and 1.
-function f0Normalized(P,G,H)
-value = f0(P,G,H)
-return value ./ (sqd_frob(G) + sqd_frob(H))
+function f0Normalized(P, G, H)
+    value = f0(P, G, H)
+    return value ./ (sqd_frob(G) + sqd_frob(H))
 end
 
 # gradient of F0 but normalized for values between 0 and 1
 # save solution value in variable "storage" for space economy
 function ∇f0Normalized!(storage, P, G, H)
     value = 2.0 .* ((G^2) * P .- 2.0 .* G * P * H .+ P * (H^2))
-    storage .= value ./ (sqd_frob(G) + sqd_frob(H))
+    return storage .= value ./ (sqd_frob(G) + sqd_frob(H))
 end
 
 # concave function F1.
 # algorithm uses only normalized version. This is just for plotting and displaying the correct data.
 function f1(P, G, H)
-constantTerm = tr(laplacian(G)^2)+tr(laplacian(H)^2)
-return .- tr(Δ(G,H)'*P) .- 2.0 .* (vec(P)' * vec(laplacian(G) * P * laplacian(H))) + constantTerm
+    constantTerm = tr(laplacian(G)^2)+tr(laplacian(H)^2)
+    return .- tr(Δ(G, H)'*P) .- 2.0 .* (vec(P)' * vec(laplacian(G) * P * laplacian(H))) + constantTerm
 end
 
 # F1 normalized for values between 0 and 1.
 function f1Normalized(P, G, H)
-value = f1(P, G, H)
-return value ./ (sqd_frob(G) + sqd_frob(H))
+    value = f1(P, G, H)
+    return value ./ (sqd_frob(G) + sqd_frob(H))
 end
 
 # gradient of F1 abut normalized for values between 0 and 1
 # save solution value in variable "storage" for space economy
 function ∇f1Normalized!(storage, P, G, H)
     # the PATH-Algorithm paper has 2.0 in front of the second term, but 4.0 should be correct.
-value = .- Δ(G,H)' .- 4.0 .* laplacian(G) * P * laplacian(H)
-storage .= value ./ (sqd_frob(G) + sqd_frob(H))
+    value = .- Δ(G, H)' .- 4.0 .* laplacian(G) * P * laplacian(H)
+    return storage .= value ./ (sqd_frob(G) + sqd_frob(H))
 end
 
 # Fλ is convex combination of F0 and F1.
 # algorithm uses only normalized version. This is just for plotting and displaying the correct data.
 function fλ(P, λ, G, H)
-return (1-λ) * f0(P, G, H)  +  λ * f1(P, G, H)
+    return (1-λ) * f0(P, G, H) + λ * f1(P, G, H)
 end
 
 # Fλ normalized for values between 0 and 1.
 function fλNormalized(P, λ, G, H)
-return (1-λ) * f0Normalized(P, G, H)  +  λ * f1Normalized(P, G, H)
+    return (1-λ) * f0Normalized(P, G, H) + λ * f1Normalized(P, G, H)
 end
 
 struct FλForP
@@ -380,13 +410,13 @@ struct FλForP
     H::Matrix{Float64}
 end
 # for the FW-algorithm we need a function that takes only P as input.
-function(fλ_struct::FλForP)(P) 
+function (fλ_struct::FλForP)(P)
     return fλNormalized(P, fλ_struct.λ, fλ_struct.G, fλ_struct.H)
 end
 
 # function flipped for maximization of the initial function and thus solving QAP
 function fλ_QAP(P, λ, G, H)
-return (1-λ) * (-f1Normalized(P, G, H))  +  λ * (-f0Normalized(P, G, H))
+    return (1-λ) * (-f1Normalized(P, G, H)) + λ * (-f0Normalized(P, G, H))
 end
 
 struct FλForP_QAP
@@ -395,16 +425,16 @@ struct FλForP_QAP
     H::Matrix{Float64}
 end
 # for the FW-algorithm we need a function that takes only P as input.
-function(fλ_struct::FλForP_QAP)(P) 
+function (fλ_struct::FλForP_QAP)(P)
     return fλ_QAP(P, fλ_struct.λ, fλ_struct.G, fλ_struct.H)
 end
 
 # gradient of FλNormalized
 # save solution value in variable "storage" for space economy
 function ∇fλ!(storageλ, storage0, storage1, P, λ, G, H)
-∇f0Normalized!(storage0, P, G, H)
-∇f1Normalized!(storage1, P, G, H)
-storageλ .= (1.0-λ) .* storage0 .+ λ .* storage1
+    ∇f0Normalized!(storage0, P, G, H)
+    ∇f1Normalized!(storage1, P, G, H)
+    return storageλ .= (1.0-λ) .* storage0 .+ λ .* storage1
 end
 struct ∇FλForP!
     storage0::Matrix{Float64}
@@ -414,16 +444,24 @@ struct ∇FλForP!
     H::Matrix{Float64}
 end
 # for the FW-algorithm we need a function that takes only P as input.
-function(∇fλ_struct::∇FλForP!)(storageλ, P)
-    ∇fλ!(storageλ, ∇fλ_struct.storage0, ∇fλ_struct.storage1, P, ∇fλ_struct.λ, ∇fλ_struct.G, ∇fλ_struct.H)
+function (∇fλ_struct::∇FλForP!)(storageλ, P)
+    return ∇fλ!(
+        storageλ,
+        ∇fλ_struct.storage0,
+        ∇fλ_struct.storage1,
+        P,
+        ∇fλ_struct.λ,
+        ∇fλ_struct.G,
+        ∇fλ_struct.H,
+    )
 end
 
 # gradient flipped for maximization and solving QAP
 # save solution value in variable "storage" for space economy
 function ∇fλ_QAP!(storageλ, storage0, storage1, P, λ, G, H)
-∇f0Normalized!(storage0, P, G, H)
-∇f1Normalized!(storage1, P, G, H)
-storageλ .= (1.0-λ) .* (-storage1) .+ λ .* (-storage0)
+    ∇f0Normalized!(storage0, P, G, H)
+    ∇f1Normalized!(storage1, P, G, H)
+    return storageλ .= (1.0-λ) .* (-storage1) .+ λ .* (-storage0)
 end
 
 struct ∇FλForP_QAP!
@@ -434,11 +472,19 @@ struct ∇FλForP_QAP!
     H::Matrix{Float64}
 end
 # for the FW-algorithm we need a function that takes only P as input.
-function(∇fλ_struct::∇FλForP_QAP!)(storageλ, P)
-    ∇fλ_QAP!(storageλ, ∇fλ_struct.storage0, ∇fλ_struct.storage1, P, ∇fλ_struct.λ, ∇fλ_struct.G, ∇fλ_struct.H)
+function (∇fλ_struct::∇FλForP_QAP!)(storageλ, P)
+    return ∇fλ_QAP!(
+        storageλ,
+        ∇fλ_struct.storage0,
+        ∇fλ_struct.storage1,
+        P,
+        ∇fλ_struct.λ,
+        ∇fλ_struct.G,
+        ∇fλ_struct.H,
+    )
 end
 
 # returns the value of the QAP objective function for a given permutation matrix P and adjacency matrices G and H
-function qapVal(P,G,H)
-return tr(G*P*H'*P')
-end
\ No newline at end of file
+function qapVal(P, G, H)
+    return tr(G*P*H'*P')
+end
diff --git a/test/pathGraphMatching.jl b/test/pathGraphMatching.jl
index 9064432..1f70dfc 100644
--- a/test/pathGraphMatching.jl
+++ b/test/pathGraphMatching.jl
@@ -2,22 +2,30 @@ using GraphsOptim
 using LinearAlgebra
 using Test
 
-P, _, _ = GraphsOptim.pathAlgorithm([
-    1.0 2.0;
-    3.0 4.0
-    ],[
-    1.0 2.0;
-    3.0 4.0
+P, _, _ = GraphsOptim.pathAlgorithm(
+    [
+        1.0 2.0;
+        3.0 4.0
     ],
-    0.1, 0.1)
+    [
+        1.0 2.0;
+        3.0 4.0
+    ],
+    0.1,
+    0.1,
+)
 @test P == [1, 2]
 
-P, _, _ = GraphsOptim.pathAlgorithm([
-    1.0 2.0;
-    3.0 4.0
-    ],[
-    4.0 3.0;
-    2.0 1.0
+P, _, _ = GraphsOptim.pathAlgorithm(
+    [
+        1.0 2.0;
+        3.0 4.0
+    ],
+    [
+        4.0 3.0;
+        2.0 1.0
     ],
-    0.1, 0.1)
-@test P == [2, 1]
\ No newline at end of file
+    0.1,
+    0.1,
+)
+@test P == [2, 1]

From 0532c1022e3e79108460a009bb3e2032b0d1aebb Mon Sep 17 00:00:00 2001
From: Ken Alan Berkpinar <ken.berkpinar@gmail.com>
Date: Mon, 3 Aug 2026 14:15:40 +0200
Subject: [PATCH 4/9] Format Code (JuliaFormatter v1)

---
 src/independent_set.jl    |  2 +-
 src/maximum_clique.jl     |  2 +-
 src/pathGraphMatching.jl  | 41 +++++++++++++++++++++------------------
 test/pathGraphMatching.jl |  8 ++++----
 4 files changed, 28 insertions(+), 25 deletions(-)

diff --git a/src/independent_set.jl b/src/independent_set.jl
index a45587c..8f7762a 100644
--- a/src/independent_set.jl
+++ b/src/independent_set.jl
@@ -16,7 +16,7 @@ function maximum_weight_independent_set!(
     model[Symbol(var_name)] = f
     @constraint(
         model,
-        covering_constraint[i = 1:nv(g), j = 1:nv(g); i ≠ j && has_edge(g, i, j)],
+        covering_constraint[i=1:nv(g), j=1:nv(g); i ≠ j && has_edge(g, i, j)],
         f[i] + f[j] <= 1,
     )
     obj = objective_function(model)
diff --git a/src/maximum_clique.jl b/src/maximum_clique.jl
index 8568afc..e0aad86 100644
--- a/src/maximum_clique.jl
+++ b/src/maximum_clique.jl
@@ -16,7 +16,7 @@ function maximum_weight_clique!(
     model[Symbol(var_name)] = f
     @constraint(
         model,
-        packing_constraint[i = 1:nv(g), j = 1:nv(g); i ≠ j && !has_edge(g, i, j)],
+        packing_constraint[i=1:nv(g), j=1:nv(g); i ≠ j && !has_edge(g, i, j)],
         f[i] + f[j] <= 1,
     )
     obj = objective_function(model)
diff --git a/src/pathGraphMatching.jl b/src/pathGraphMatching.jl
index 7f706b1..67ed188 100644
--- a/src/pathGraphMatching.jl
+++ b/src/pathGraphMatching.jl
@@ -38,7 +38,7 @@ function pathAlgorithm(
     verbose_FW::Bool=false,
 )
     # extend the smaller matrix by zero rows and columns
-    diffSize = size(G, 1)-size(H, 1)
+    diffSize = size(G, 1) - size(H, 1)
 
     if diffSize > 0
         # G is larger
@@ -148,11 +148,13 @@ function pathAlgorithm(
         # update dλ until criterion is met
         # TODO implemented new stopping criterion. Need to still find out ϵ_f and ϵ_p values from FrankWolfe implementation and calculate ϵ_λ_f and ϵ_λ_p with added input M.
         # d_λ is doubled until one value is larger than it's threshold (or new λ is already 1)
-        while abs(fλNormalizedFinal(p_new, λ_new, G, H)-fλNormalizedFinal(p_opt, λ, G, H)) <
-              ϵ_λ_f &&
-              p_change_normalized < ϵ_λ_p &&
-              λ_new < one(Float64)
-            dλ = 2*dλ
+        while abs(
+                      fλNormalizedFinal(p_new, λ_new, G, H) -
+                      fλNormalizedFinal(p_opt, λ, G, H),
+                  ) < ϵ_λ_f &&
+                  p_change_normalized < ϵ_λ_p &&
+                  λ_new < one(Float64)
+            dλ = 2 * dλ
             λ_new = min(λ + dλ, one(Float64))
 
             verbose && println("   dλ = ", dλ)
@@ -177,7 +179,7 @@ function pathAlgorithm(
         end
 
         # if the last while loop's condition is not met (anymore), dλ is one step too large and can be halved once directly
-        dλ = max(dλ/2, dλ_min)
+        dλ = max(dλ / 2, dλ_min)
         λ_new = λ + dλ
         verbose && println("   dλ = ", dλ)
         if !isnothing(p_last)
@@ -204,10 +206,10 @@ function pathAlgorithm(
 
         # d_λ is halved until both values are smaller than their thresholds (or dλ is already at minimum)
         while (
-            abs(fλNormalizedFinal(p_new, λ_new, G, H)-fλNormalizedFinal(p_opt, λ, G, H)) >
+            abs(fλNormalizedFinal(p_new, λ_new, G, H) - fλNormalizedFinal(p_opt, λ, G, H)) >
             ϵ_λ_f || p_change_normalized > ϵ_λ_p
         ) && dλ > dλ_min
-            dλ = max(dλ/2, dλ_min)
+            dλ = max(dλ / 2, dλ_min)
             λ_new = min(λ + dλ, one(Float64))
             verbose && println("   dλ = ", dλ)
 
@@ -356,7 +358,7 @@ end
 # convex function F0
 # algorithm uses only normalized version. This is just for plotting and displaying the correct data.
 function f0(P, G, H)
-    return sqd_frob(G*P .- P*H)
+    return sqd_frob(G * P .- P * H)
 end
 
 # F0 normalized for values between 0 and 1.
@@ -375,8 +377,9 @@ end
 # concave function F1.
 # algorithm uses only normalized version. This is just for plotting and displaying the correct data.
 function f1(P, G, H)
-    constantTerm = tr(laplacian(G)^2)+tr(laplacian(H)^2)
-    return .- tr(Δ(G, H)'*P) .- 2.0 .* (vec(P)' * vec(laplacian(G) * P * laplacian(H))) + constantTerm
+    constantTerm = tr(laplacian(G)^2) + tr(laplacian(H)^2)
+    return .-tr(Δ(G, H)' * P) .- 2.0 .* (vec(P)' * vec(laplacian(G) * P * laplacian(H))) +
+           constantTerm
 end
 
 # F1 normalized for values between 0 and 1.
@@ -389,19 +392,19 @@ end
 # save solution value in variable "storage" for space economy
 function ∇f1Normalized!(storage, P, G, H)
     # the PATH-Algorithm paper has 2.0 in front of the second term, but 4.0 should be correct.
-    value = .- Δ(G, H)' .- 4.0 .* laplacian(G) * P * laplacian(H)
+    value = .-Δ(G, H)' .- 4.0 .* laplacian(G) * P * laplacian(H)
     return storage .= value ./ (sqd_frob(G) + sqd_frob(H))
 end
 
 # Fλ is convex combination of F0 and F1.
 # algorithm uses only normalized version. This is just for plotting and displaying the correct data.
 function fλ(P, λ, G, H)
-    return (1-λ) * f0(P, G, H) + λ * f1(P, G, H)
+    return (1 - λ) * f0(P, G, H) + λ * f1(P, G, H)
 end
 
 # Fλ normalized for values between 0 and 1.
 function fλNormalized(P, λ, G, H)
-    return (1-λ) * f0Normalized(P, G, H) + λ * f1Normalized(P, G, H)
+    return (1 - λ) * f0Normalized(P, G, H) + λ * f1Normalized(P, G, H)
 end
 
 struct FλForP
@@ -416,7 +419,7 @@ end
 
 # function flipped for maximization of the initial function and thus solving QAP
 function fλ_QAP(P, λ, G, H)
-    return (1-λ) * (-f1Normalized(P, G, H)) + λ * (-f0Normalized(P, G, H))
+    return (1 - λ) * (-f1Normalized(P, G, H)) + λ * (-f0Normalized(P, G, H))
 end
 
 struct FλForP_QAP
@@ -434,7 +437,7 @@ end
 function ∇fλ!(storageλ, storage0, storage1, P, λ, G, H)
     ∇f0Normalized!(storage0, P, G, H)
     ∇f1Normalized!(storage1, P, G, H)
-    return storageλ .= (1.0-λ) .* storage0 .+ λ .* storage1
+    return storageλ .= (1.0 - λ) .* storage0 .+ λ .* storage1
 end
 struct ∇FλForP!
     storage0::Matrix{Float64}
@@ -461,7 +464,7 @@ end
 function ∇fλ_QAP!(storageλ, storage0, storage1, P, λ, G, H)
     ∇f0Normalized!(storage0, P, G, H)
     ∇f1Normalized!(storage1, P, G, H)
-    return storageλ .= (1.0-λ) .* (-storage1) .+ λ .* (-storage0)
+    return storageλ .= (1.0 - λ) .* (-storage1) .+ λ .* (-storage0)
 end
 
 struct ∇FλForP_QAP!
@@ -486,5 +489,5 @@ end
 
 # returns the value of the QAP objective function for a given permutation matrix P and adjacency matrices G and H
 function qapVal(P, G, H)
-    return tr(G*P*H'*P')
+    return tr(G * P * H' * P')
 end
diff --git a/test/pathGraphMatching.jl b/test/pathGraphMatching.jl
index 1f70dfc..8b37acc 100644
--- a/test/pathGraphMatching.jl
+++ b/test/pathGraphMatching.jl
@@ -4,11 +4,11 @@ using Test
 
 P, _, _ = GraphsOptim.pathAlgorithm(
     [
-        1.0 2.0;
+        1.0 2.0
         3.0 4.0
     ],
     [
-        1.0 2.0;
+        1.0 2.0
         3.0 4.0
     ],
     0.1,
@@ -18,11 +18,11 @@ P, _, _ = GraphsOptim.pathAlgorithm(
 
 P, _, _ = GraphsOptim.pathAlgorithm(
     [
-        1.0 2.0;
+        1.0 2.0
         3.0 4.0
     ],
     [
-        4.0 3.0;
+        4.0 3.0
         2.0 1.0
     ],
     0.1,

From b5b136d580dd2bb073d83b36a5fc9c4d07fa4c34 Mon Sep 17 00:00:00 2001
From: Ken Alan Berkpinar <ken.berkpinar@gmail.com>
Date: Mon, 3 Aug 2026 14:48:23 +0200
Subject: [PATCH 5/9] Reformat Project.TOML

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index cef71db..7b921ea 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "GraphsOptim"
 uuid = "e79ef3ae-79c8-4b97-b165-63f338db30c2"
-version = "0.1.0"
 authors = ["Guillaume Dalle, Aurora Rossi and contributors"]
+version = "0.1.0"
 
 [deps]
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"

From 979ab9ee8337554407c3e2a1af35deb7d800a04e Mon Sep 17 00:00:00 2001
From: Ken Alan Berkpinar <ken.berkpinar@gmail.com>
Date: Mon, 3 Aug 2026 15:10:43 +0200
Subject: [PATCH 6/9] fix wrong case distinction

---
 src/pathGraphMatching.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/pathGraphMatching.jl b/src/pathGraphMatching.jl
index 67ed188..56f5035 100644
--- a/src/pathGraphMatching.jl
+++ b/src/pathGraphMatching.jl
@@ -98,8 +98,8 @@ function pathAlgorithm(
     λ = 0.0
 
     # redefine f0, f1 and fλ depending on whether the QAP should be solved or not, s.t. f0 is always convex and f1 is always concave.
-    if !solveQAP
-        fλNormalizedFinal = fλNormalized
+    fλNormalizedFinal = fλNormalized
+    if solveQAP
         fλNormalizedFinal = fλ_QAP
     end
 

From e8572998812f1f3c398922eede60214e40674b86 Mon Sep 17 00:00:00 2001
From: Ken Alan Berkpinar <ken.berkpinar@gmail.com>
Date: Tue, 4 Aug 2026 10:43:13 +0200
Subject: [PATCH 7/9] Change Code linting test case (temporarily) to check if
 git accepts

---
 test/runtests.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 88eb3a7..890dc0e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -19,7 +19,7 @@ using Test
     end
 
     @testset "Code linting" begin
-        JET.test_package(GraphsOptim; target_defined_modules=true)
+        JET.test_package(GraphsOptim; target_modules=(GraphsOptim,))
     end
 
     @testset "Doctests" begin

From e0c3e462ae93e647b95f17e83a89834d4e9a928b Mon Sep 17 00:00:00 2001
From: Ken Alan Berkpinar <ken.berkpinar@gmail.com>
Date: Tue, 4 Aug 2026 14:11:53 +0200
Subject: [PATCH 8/9] change back code linting test

---
 test/runtests.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 890dc0e..88eb3a7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -19,7 +19,7 @@ using Test
     end
 
     @testset "Code linting" begin
-        JET.test_package(GraphsOptim; target_modules=(GraphsOptim,))
+        JET.test_package(GraphsOptim; target_defined_modules=true)
     end
 
     @testset "Doctests" begin

From 642c03b473548c3dc6af389019693b47de7b3bec Mon Sep 17 00:00:00 2001
From: Ken Alan Berkpinar <ken.berkpinar@gmail.com>
Date: Thu, 6 Aug 2026 20:21:16 +0200
Subject: [PATCH 9/9] JuliaFormatter on Julia 1.10

---
 src/independent_set.jl   |  2 +-
 src/maximum_clique.jl    |  2 +-
 src/pathGraphMatching.jl | 12 +++++-------
 3 files changed, 7 insertions(+), 9 deletions(-)

diff --git a/src/independent_set.jl b/src/independent_set.jl
index 8f7762a..a45587c 100644
--- a/src/independent_set.jl
+++ b/src/independent_set.jl
@@ -16,7 +16,7 @@ function maximum_weight_independent_set!(
     model[Symbol(var_name)] = f
     @constraint(
         model,
-        covering_constraint[i=1:nv(g), j=1:nv(g); i ≠ j && has_edge(g, i, j)],
+        covering_constraint[i = 1:nv(g), j = 1:nv(g); i ≠ j && has_edge(g, i, j)],
         f[i] + f[j] <= 1,
     )
     obj = objective_function(model)
diff --git a/src/maximum_clique.jl b/src/maximum_clique.jl
index e0aad86..8568afc 100644
--- a/src/maximum_clique.jl
+++ b/src/maximum_clique.jl
@@ -16,7 +16,7 @@ function maximum_weight_clique!(
     model[Symbol(var_name)] = f
     @constraint(
         model,
-        packing_constraint[i=1:nv(g), j=1:nv(g); i ≠ j && !has_edge(g, i, j)],
+        packing_constraint[i = 1:nv(g), j = 1:nv(g); i ≠ j && !has_edge(g, i, j)],
         f[i] + f[j] <= 1,
     )
     obj = objective_function(model)
diff --git a/src/pathGraphMatching.jl b/src/pathGraphMatching.jl
index 56f5035..a2208e1 100644
--- a/src/pathGraphMatching.jl
+++ b/src/pathGraphMatching.jl
@@ -149,11 +149,10 @@ function pathAlgorithm(
         # TODO implemented new stopping criterion. Need to still find out ϵ_f and ϵ_p values from FrankWolfe implementation and calculate ϵ_λ_f and ϵ_λ_p with added input M.
         # d_λ is doubled until one value is larger than it's threshold (or new λ is already 1)
         while abs(
-                      fλNormalizedFinal(p_new, λ_new, G, H) -
-                      fλNormalizedFinal(p_opt, λ, G, H),
-                  ) < ϵ_λ_f &&
-                  p_change_normalized < ϵ_λ_p &&
-                  λ_new < one(Float64)
+                  fλNormalizedFinal(p_new, λ_new, G, H) - fλNormalizedFinal(p_opt, λ, G, H)
+              ) < ϵ_λ_f &&
+              p_change_normalized < ϵ_λ_p &&
+              λ_new < one(Float64)
             dλ = 2 * dλ
             λ_new = min(λ + dλ, one(Float64))
 
@@ -378,8 +377,7 @@ end
 # algorithm uses only normalized version. This is just for plotting and displaying the correct data.
 function f1(P, G, H)
     constantTerm = tr(laplacian(G)^2) + tr(laplacian(H)^2)
-    return .-tr(Δ(G, H)' * P) .- 2.0 .* (vec(P)' * vec(laplacian(G) * P * laplacian(H))) +
-           constantTerm
+    return .-tr(Δ(G, H)' * P) .- 2.0 .* (vec(P)' * vec(laplacian(G) * P * laplacian(H))) + constantTerm
 end
 
 # F1 normalized for values between 0 and 1.