Isomorphism

[scheinerman]

A natural addition is a function to find isomorphisms between graphs. I've written some code that does that. (I'm sure it works for simple graphs and simple digraphs. Not sure about weighted.) I hope you find this useful.

using Graphs
using JuMP
using HiGHS
using ChooseOptimizer

set_solver(HiGHS)
set_solver_verbose(false)

_graphs_not_iso() = error("The graphs are not isomorphic")

"""
    graph_iso(g::AbstractGraph, h::AbstractGraph)::Dict{Int,Int}

Return an isomorphism `f` from `g` to `h`, or throw an error if the 
graphs are not isomorphic. Here `f` is a `Dict` with the property that
`(u,v)` is an edge of `g` if and only if `(f[u],f[v])` is an edge of `h`.
"""
function graph_iso(g::AbstractGraph, h::AbstractGraph)::Dict{Int,Int}
    n = nv(g)
    if n != nv(h)
        _graphs_not_iso()
    end

    # other simple tests can go here, for example:
    if ne(g) != ne(h)
        _graphs_not_iso()
    end

    A = adjacency_matrix(g)
    B = adjacency_matrix(h)
    MOD = Model(get_solver())
    @variable(MOD, X[1:n, 1:n], Bin)

    # success when A*X = X*B and X is permutation 

    # X must be doubly stochastic ==> permutation
    for i = 1:n
        @constraint(MOD, sum(X[i, j] for j = 1:n) == 1)
        @constraint(MOD, sum(X[j, i] for j = 1:n) == 1)
    end

    # A*X = X*B 
    for i = 1:n
        for k = 1:n
            @constraint(
                MOD,
                sum(A[i, j] * X[j, k] for j = 1:n) == sum(X[i, j] * B[j, k] for j = 1:n)
            )
        end
    end

    optimize!(MOD)
    status = Int(termination_status(MOD))

    # status == 1 means success
    if status != 1
        _graphs_not_iso()
    end

    # get the matrix
    P = Int.(value.(X))

    # convert to a dictionary
    result = Dict{Int,Int}()
    for v=1:n 
        for w=1:n 
            if P[v,w] > 0
                result[v] = w 
            end
        end
    end

    return result
end