WIP: Algorithm for counting the number of consistent extensions (size of an MEC given a CPDAG)

Project.toml

@@ -11,6 +11,7 @@ Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Memoization = "6fafb56a-5788-4b4e-91ca-c0cea6611c73"
+LinkedLists = "70f5e60a-1556-5f34-a19e-a48b3e4aaee9"
 MetaGraphs = "626554b9-1ddb-594c-aa3c-2596fe9399a5"
 NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
@@ -22,6 +23,15 @@ Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 TabularDisplay = "3eeacb1d-13c2-54cc-9b18-30c86af3cadb"
 ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d"
 
+[weakdeps]
+GraphRecipes = "bd48cda9-67a9-57be-86fa-5b3c104eda73"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+TikzGraphs = "b4f28e30-c73f-5eaf-a395-8a9db949a742"
+
+[extensions]
+GraphRecipesExt = ["GraphRecipes", "Plots"]
+TikzGraphsExt = "TikzGraphs"
+
 [compat]
 Combinatorics = "1.0"
 DelimitedFiles = "1.6, 1.7, 1.8, 1.9"
@@ -43,10 +53,6 @@ ThreadsX = "0.1"
 TikzGraphs = "1.3, 1.4"
 julia = "1.6, 1.7, 1.8, 1.9"
 
-[extensions]
-GraphRecipesExt = ["GraphRecipes", "Plots"]
-TikzGraphsExt = "TikzGraphs"
-
 [extras]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 GraphRecipes = "bd48cda9-67a9-57be-86fa-5b3c104eda73"

src/CausalInference.jl

@@ -4,6 +4,7 @@ using Graphs
 using Graphs.SimpleGraphs
 using Combinatorics
 using Base.Iterators
+using LinkedLists
 using Memoization, LRUCache
 using ThreadsX
 
@@ -12,7 +13,7 @@ import Base: ==, show
 export ancestors, descendants, alt_test_dsep, test_covariate_adjustment, alt_test_backdoor, find_dsep, find_min_dsep, find_covariate_adjustment, find_backdoor_adjustment, find_frontdoor_adjustment, find_min_covariate_adjustment, find_min_backdoor_adjustment, find_min_frontdoor_adjustment, list_dseps, list_covariate_adjustment, list_backdoor_adjustment, list_frontdoor_adjustment
 export dsep, skeleton, gausscitest, dseporacle, partialcor
 export unshielded, pcalg, vskel, vskel!, alt_vskel
-export cpdag, alt_cpdag, meek_rules!
+export cpdag, alt_cpdag, meek_rules!, MECsize
 export digraph, vpairs, skel_oracle, pc_oracle, randdag
 export cmitest, kl_entropy, kl_renyi, kl_mutual_information
 export kl_cond_mi, kl_perm_mi_test, kl_perm_cond_mi_test
@@ -44,6 +45,7 @@ include("backdoor.jl")
 include("ges.jl")
 include("gensearch.jl")
 include("workloads.jl")
+include("chordal.jl")
 
 # Compatibility with the new "Package Extensions" (https://github.com/JuliaLang/julia/pull/47695)
 const EXTENSIONS_SUPPORTED = isdefined(Base, :get_extension)

src/chordal.jl

@@ -0,0 +1,165 @@
+import LinkedLists
+
+"""
+    ischordal(g)
+
+Return true if the given graph is chordal
+"""
+function ischordal(G)
+    mcsorder, invmcsorder, _ = mcs(G, Set())
+
+    n = length(mcsorder)
+
+    f = zeros(Int, n)
+    index = zeros(Int, n)
+    for i=n:-1:1
+        w = mcsorder[i]
+        f[w] = w
+        index[w] = i
+        for v in neighbors(G, w)
+            if invmcsorder[v] > i
+                index[v] = i
+                if f[v] == v
+                    f[v] = w
+                end
+            end
+        end
+        for v in neighbors(G, w)
+            if invmcsorder[v] > i
+                if index[f[v]] > i
+                    return false
+                end
+            end
+        end
+    end
+    return true
+end
+
+function cliquetreefrommcs(G, mcsorder, invmcsorder)
+    n = nv(G)
+    # data structures for the algorithm
+    K = Vector{Set}()
+    push!(K, Set())
+    s = 1
+    edgelist = Set{Edge}()
+    visited = falses(n)
+    clique = zeros(Int, n)
+
+    for i = 1:n
+        x = mcsorder[i]
+        S = Set{Int}()
+        for w in inneighbors(G, x)
+            if visited[w]
+                push!(S, w)
+            end
+        end
+
+        # if necessary create new maximal clique
+        if K[s] != S
+            s += 1
+            push!(K, S)
+            k, _ = findmax(map(x -> invmcsorder[x], collect(S)))
+            p = clique[mcsorder[k]]
+            push!(edgelist, Edge(p, s))
+        end
+
+        union!(K[s], x)
+        clique[x] = s 
+        visited[x] = true;
+    end
+
+    T = SimpleGraphFromIterator(edgelist)
+    # ensure graph is not empty
+    nv(T) == 0 && add_vertices!(T,1)
+    return K, T
+end
+
+@inline function vispush!(l::LinkedList, pointers, x, vis)
+    if vis
+        pointers[x] = push!(l,x)
+    else
+        pointers[x] = pushfirst!(l,x)
+    end
+end
+
+# TODO: separate mcs and mcs plus cgk??
+# Returns the visit order of the vertices, its inverse and the subgraphs C_G(K) (see Def. 1 in [1,2]). If K is empty a normal MCS is performed.
+function mcs(G, K)
+    n = nv(G)
+    copy_K = copy(K)
+
+    # data structures for MCS
+    sets = [LinkedList{Int}() for _ = 1:n+1]
+    pointers = Vector(undef,n)
+    size = Vector{Int}(undef, n)
+    visited = falses(n)
+
+    # output data structures
+    mcsorder = Vector{Int}(undef, n)
+    invmcsorder = Vector{Int}(undef, n)
+    subgraphs = Array[]
+
+    # init
+    visited[collect(copy_K)] .= true
+    for v in vertices(G)
+        size[v] = 1
+        vispush!(sets[1], pointers, v, visited[v])
+    end
+    maxcard = 1
+
+    for i = 1:n
+        # first, the vertices in K are chosen
+        # they are always in the set of maximum cardinality vertices
+        if !isempty(copy_K)
+            v = pop!(copy_K)
+        # afterwards, the algorithm chooses any vertex from maxcard
+        else
+            v = first(sets[maxcard])
+        end
+        # v is the ith vertex in the mcsorder
+        mcsorder[i] = v
+        invmcsorder[v] = i
+        size[v] = -1
+
+        # immediately append possible subproblems to the output
+        if !visited[v]
+            vertexset = Vector{Int}()
+            for x in sets[maxcard]
+                visited[x] && break
+                visited[x] = true
+                push!(vertexset, x)
+            end
+            sg = induced_subgraph(G, vertexset)
+            subgraphs = vcat(subgraphs, (map(x -> sg[2][x], connected_components(sg[1]))))
+        end
+
+        deleteat!(sets[maxcard], pointers[v])
+
+        # update the neighbors
+        for w in inneighbors(G, v)
+            if size[w] >= 1
+                deleteat!(sets[size[w]], pointers[w])
+                size[w] += 1
+                vispush!(sets[size[w]], pointers, w, visited[w])
+            end
+        end
+        maxcard += 1
+        while maxcard >= 1 && isempty(sets[maxcard])
+            maxcard -= 1
+        end
+    end
+
+    return mcsorder, invmcsorder, subgraphs
+end
+
+"""
+    cliquetree(G)
+
+Computes a clique tree of a graph G. A vector K of maximal cliques and a tree T on 1,2,...,|K| is returned.
+
+"""
+function cliquetree(G)
+    mcsorder, invmcsorder, _ = mcs(G, Set())
+    K, T = cliquetreefrommcs(G, mcsorder, invmcsorder)
+    return K, T
+end