Merge pull request #166 from mschauer/do

Add do operator, make new version, improve API of pcalg.

Author: Moritz Schauer

Project.toml

@@ -26,54 +26,56 @@ Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 TabularDisplay = "3eeacb1d-13c2-54cc-9b18-30c86af3cadb"
 ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d"
 
+[weakdeps]
+GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2"
+GraphRecipes = "bd48cda9-67a9-57be-86fa-5b3c104eda73"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+TikzGraphs = "b4f28e30-c73f-5eaf-a395-8a9db949a742"
+
+[extensions]
+GraphMakieExt = "GraphMakie"
+GraphRecipesExt = ["GraphRecipes", "Plots"]
+TikzGraphsExt = "TikzGraphs"
+
 [compat]
 Combinatorics = "1.0"
 DelimitedFiles = "1.6, 1.7, 1.8, 1.9"
 Distances = "0.8, 0.9, 0.10"
 Distributions = "0.22, 0.23, 0.24, 0.25"
 GraphMakie = "0.5"
 GraphRecipes = "0.5"
-Graphs = "1.5, 1.6, 1.7, 1.8, 1.9"
-LRUCache = "1.4, 1.5"
+Graphs = "1.5, 1.6, 1.7, 1.8, 1.9, 1.10, 1.11, 1.12"
+LRUCache = "1.4, 1.5, 1.6"
 LinearAlgebra = "1"
 LinkedLists = "0.1"
 LogarithmicNumbers = "1.4"
 Memoization = "0.2"
-MetaGraphs = "0.7"
+MetaGraphs = "0.7, 0.8"
 NearestNeighbors = "0.4"
-OffsetArrays = "1.12"
+OffsetArrays = "1.12, 1.13, 1.14"
+OrderedCollections = "1.7.0"
 Plots = "1"
 PrecompileTools = "1.1, 1.2"
-ProgressMeter = "1.9"
+ProgressMeter = "1.9, 1.10"
 Random = "1"
 Requires = "1.3"
-SpecialFunctions = "1.8, 2.0, 2.1, 2.2, 2.3"
+SpecialFunctions = "1.8, 2.0, 2.1, 2.2, 2.3, 2.4"
 Statistics = "1"
-Tables = "1.6, 1.7, 1.8, 1.9, 1.10, 1.11"
+Tables = "1.6, 1.7, 1.8, 1.9, 1.10, 1.11, 1.12"
 TabularDisplay = "1.2"
 ThreadsX = "0.1"
 TikzGraphs = "1.3, 1.4"
 julia = "1.6, 1.7, 1.8, 1.9, 1.10"
 
-[extensions]
-GraphMakieExt = "GraphMakie"
-GraphRecipesExt = ["GraphRecipes", "Plots"]
-TikzGraphsExt = "TikzGraphs"
-
 [extras]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2"
 GraphRecipes = "bd48cda9-67a9-57be-86fa-5b3c104eda73"
+OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TikzGraphs = "b4f28e30-c73f-5eaf-a395-8a9db949a742"
 
 [targets]
-test = ["Test", "StatsBase", "DelimitedFiles"]
-
-[weakdeps]
-GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2"
-GraphRecipes = "bd48cda9-67a9-57be-86fa-5b3c104eda73"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-TikzGraphs = "b4f28e30-c73f-5eaf-a395-8a9db949a742"
+test = ["Test", "StatsBase", "DelimitedFiles", "OrderedCollections"]

docs/src/library.md

@@ -33,6 +33,7 @@ CausalInference.meek_rule4
 pdag2dag!
 pdag_to_dag_dortarsi!
 pdag_to_dag_meek!
+CausalInference.do!
 ```
 
 ## PC algorithm

src/CausalInference.jl

@@ -41,6 +41,7 @@ include("combinations_without.jl")
 include("klentropy.jl")
 include("skeleton.jl")
 include("pdag.jl")
+include("do.jl")
 include("dsep.jl")
 include("meek.jl")
 include("pc.jl")

src/do.jl

@@ -0,0 +1,14 @@
+
+""""
+    do!(g, v)
+
+Graphical do operator, removes all 
+incoming edges to vertex `v` in a DiGraph `g`.
+Returns the modified graph.
+"""
+function do!(g::DiGraph{T}, v::T) where {T}
+    for u in collect(inneighbors(g, v))
+        rem_edge!(g, u, v)
+    end
+    return g
+end

src/pc.jl

@@ -114,9 +114,30 @@ function _vskel(n::V, I, par...) where {V}
     dg
 end
 
+"""
+    pcalg(g, I; stable=true)
+
+Perform the PC algorithm for a set of 1:n variables using the tests
+
+    I(u, v, [s1, ..., sn])
+
+Use `IClosure(I, args)` to wrap a function f with signature
+
+    f(u, v, [s1, ..., sn], par...)
+
+Returns the CPDAG as DiGraph. By default uses a stable and threaded versions
+of the skeleton algorithm. (This is the most recent interface)
+"""
+function pcalg(g::Graph, I; stable=true)
+    g, S = stable ? skeleton_stable(g, I) : skeleton(g, I)
+    dg = DiGraph(g) # use g to keep track of unoriented edges
+    _, dg = orient_unshielded(g, dg, S)
+    meek_rules!(dg)
+    dg
+end
+
 """
     pcalg(n::V, I, par...; stable=true)
-    pcalg(g, I, par...; stable=true)
 
 Perform the PC algorithm for a set of 1:n variables using the tests
 

src/skeleton.jl

@@ -83,7 +83,7 @@ function skeleton(g::SimpleGraph, I)
             @label nextedge
         end
         d = d + 1
-        if isdone
+        if isdone || d > n-2 
             return g, S
         end
     end
@@ -136,7 +136,7 @@ function skeleton_stable(g, I)
         end
         empty!(remove)
         d = d + 1
-        if isdone
+        if isdone || d > n-2 
             return g, S
         end
     end

test/dag_sampler.jl

@@ -1,6 +1,7 @@
 using Random, CausalInference, Statistics, Test, Graphs, LinearAlgebra
+using OrderedCollections
+sort2(d; kwargs...) = sort!(OrderedDict(d); kwargs...)
 @testset "Dag-Zig-Zag" begin
-
     Random.seed!(1)
 
     N = 100 # number of data points
@@ -28,12 +29,12 @@ using Random, CausalInference, Statistics, Test, Graphs, LinearAlgebra
             graphs, graph_pairs, hs, τs, ws, ts, scores = CausalInference.unzipgs(gs)
             @test score_dag(graphs[end], score) ≈ scores[end] + score_dag(DiGraph(n), score)
             @test !any(Graphs.is_cyclic.(graphs))
-            posterior = sort(keyedreduce(+, graph_pairs, ws); byvalue=true, rev=true)
+            posterior = sort2(keyedreduce(+, graph_pairs, ws); byvalue=true, rev=true)
             logΠ = map(g->score_dag(digraph(g, n), score), collect(keys(posterior)))
             Π = normalize(exp.(logΠ .- maximum(logΠ) ), 1)
             @test norm(collect(values(posterior)) - Π, 1) < 30/sqrt(iterations)
 
-            posterior = sort(keyedreduce(+, vpairs.(cpdag.(graphs)), ws); byvalue=true, rev=true)
+            posterior = sort2(keyedreduce(+, vpairs.(cpdag.(graphs)), ws); byvalue=true, rev=true)
             @test first(posterior).first == [1=>2, 1=>3, 2=>1, 2=>4, 3=>1, 3=>4, 4=>5] 
         end
     end
@@ -50,7 +51,7 @@ end #testset
     n = 3
     gs = dagzigzag(n; iterations);
     graphs, graph_pairs, hs, τs, ws, ts, scores = CausalInference.unzipgs(gs);
-    posterior = sort(keyedreduce(+, graph_pairs, ws); byvalue=true, rev=true)
+    posterior = sort2(keyedreduce(+, graph_pairs, ws); byvalue=true, rev=true)
     @test length(posterior) == 25
     @test maximum(values(posterior)) < 1/25 + 0.01
     @test minimum(values(posterior)) > 1/25 - 0.01

test/equations.jl

@@ -17,24 +17,21 @@ df = (x=x, v=v, w=w, z=z, s=s)
 
 est_g, score = ges(df; penalty=1.0, parallel=true)
 
-
 est_dag= pdag2dag!(est_g)
 
 scm= estimate_equations(df,est_dag)
 
-display(scm)
+#display(scm)
 
 #println(CI.SCM)
 
-df_generated= generate_data(scm, 2000)
+df_generated = generate_data(scm, 2000)
 
 println("df: ")
 
-display(df)
+#display(df)
 
 println("df_generated: ")
 
-
-
-display(df_generated)
+#display(df_generated)
 

test/runtests.jl

@@ -1,7 +1,7 @@
 using CausalInference
 using Graphs
 using Test
-
+include("zintervention.jl")
 include("exact.jl")
 include("operators.jl")
 include("ges.jl")

test/sampler.jl

@@ -1,4 +1,6 @@
 using Random, CausalInference, Statistics, Test, Graphs, LinearAlgebra
+using OrderedCollections
+sort2(d; kwargs...) = sort!(OrderedDict(d); kwargs...)
 @testset "Zig-Zag" begin
     Random.seed!(1)
 
@@ -21,7 +23,7 @@ using Random, CausalInference, Statistics, Test, Graphs, LinearAlgebra
     score = GaussianScore(C, N, penalty)
     gs = causalzigzag(n; score, κ, iterations)
     graphs, graph_pairs, hs, τs, ws, ts, scores = CausalInference.unzipgs(gs)
-    posterior = sort(keyedreduce(+, graph_pairs, ws); byvalue=true, rev=true)
+    posterior = sort2(keyedreduce(+, graph_pairs, ws); byvalue=true, rev=true)
 
     # maximum aposteriori estimate
     @test first(posterior).first == [1=>2, 1=>3, 2=>1, 2=>4, 3=>1, 3=>4, 4=>5] 
@@ -46,14 +48,14 @@ end #testset
 
     gs = causalzigzag(n; iterations);
     graphs, graph_pairs, hs, τs, ws, ts, scores = CausalInference.unzipgs(gs);
-    posterior = sort(keyedreduce(+, graph_pairs, ws); byvalue=true, rev=true)
+    posterior = sort2(keyedreduce(+, graph_pairs, ws); byvalue=true, rev=true)
     @test length(posterior) == m
     @test norm(values(posterior) .- 1/m, 1)/2 < 0.05
     T_skew = sum(τs)
 
     gs = causalzigzag(n; iterations, σ=1.0, ρ=0.0);
     graphs, graph_pairs, hs, τs, ws, ts, scores = CausalInference.unzipgs(gs);
-    posterior = sort(keyedreduce(+, graph_pairs, ws); byvalue=true, rev=true)
+    posterior = sort2(keyedreduce(+, graph_pairs, ws); byvalue=true, rev=true)
     @test length(posterior) == m
     @test norm(values(posterior) .- 1/m, 1)/2 < 0.05
     T = sum(τs)