Merge pull request #118 from mschauer/cpdag

Performance improvements

Author: mschauer

src/ges.jl

@@ -206,7 +206,8 @@ function tails_and_adj_neighbors(g, x, y)
     Nb[.~ a], Nb[a]
 end 
 function adj_neighbors(g, x, y)
-    intersect(inneighbors(g,y), outneighbors(g,y), all_neighbors(g,x))
+#    a = intersect(inneighbors(g,y), outneighbors(g,y), all_neighbors(g,x))
+    sorted_intersect_(neighbors_undirected(g,y), all_neighbors(g,x))
 end 
 
 

src/meek.jl

@@ -28,35 +28,38 @@ function meek_rules!(g; rule4=false)
 end
 
 """
-    meek_rule1(dg, v, w)
+    meek_rule1(g, v, w)
 
 Rule 1: Orient v-w into v->w whenever there is u->v
 such that u and w are not adjacent
 (otherwise a new v-structure is created.)
 """
-function meek_rule1(dg, v, w)
-    for u in inneighbors(dg, v)
-        has_edge(dg, v => u) && continue # not directed
-        isadjacent(dg, u, w) && continue
+function meek_rule1(g, v, w)
+    for u in inneighbors(g, v)
+        u == w && continue 
+        has_edge(g, v => u) && continue # not directed
+        isadjacent(g, u, w) && continue
         return true
     end
     return false
 end
 
 """
-    meek_rule2(dg, v, w)
+    meek_rule2(g, v, w)
 
 Rule 2: Orient v-w into v->w whenever there is a chain v->k->w
 (otherwise a directed cycle is created.)
 """
-function meek_rule2(dg, v, w)
+function meek_rule2(g, v, w)
     outs = Int[]
-    for k in outneighbors(dg, v)
-        !has_edge(dg, k => v) && push!(outs, k)
+    for k in outneighbors(g, v)
+        k == w && continue 
+        !has_edge(g, k => v) && push!(outs, k)
     end
     ins = Int[]
-    for k in inneighbors(dg, w)
-        !has_edge(dg, w => k) && push!(ins, k)
+    for k in inneighbors(g, w)
+        k == v && continue 
+        !has_edge(g, w => k) && push!(ins, k)
     end
     if !disjoint_sorted(ins, outs)
         return true
@@ -65,42 +68,42 @@ function meek_rule2(dg, v, w)
 end
 
 """
-    meek_rule3(dg, v, w)
+    meek_rule3(g, v, w)
 
 Rule 3 (Diagonal): Orient v-w into v->w whenever there are two chains
 v-k->w and v-l->w such that k and l are nonadjacent
 (otherwise a new v-structure or a directed cycle is created.)
 """
-function meek_rule3(dg, v, w)
-    fulls = [] # Find nodes k where v-k
-    for k in outneighbors(dg, v)
-        has_edge(dg, k => v) || continue 
+function meek_rule3(g, v, w)
+    fulls = Int[] # Find nodes k where v-k
+    for k in outneighbors(g, v)
+        has_edge(g, k => v) || continue 
         # Skip if not k->w (or if not l->w)
-        if has_edge(dg, w => k) || !has_edge(dg, k => w)
+        if has_edge(g, w => k) || !has_edge(g, k => w)
             continue
         end
         push!(fulls, k)
     end
     for (k, l) in combinations(fulls, 2) # FIXME: 
-        isadjacent(dg, k, l) && continue
+        isadjacent(g, k, l) && continue
         return true
     end
     return false
 end
 
 """
-    meek_rule4(dg, v, w)
+    meek_rule4(g, v, w)
 
 Rule 4: Orient v-w into v→w if v-k→l→w where adj(v,l) and not adj(k,w) [check].
 """
-function meek_rule4(dg, v, w)
-    for l in inneighbors(dg, w)
-        has_edge(dg, w => l) && continue # undirected
-        !isadjacent(dg, v, l) && continue # not adjacent to v      
-        for k in inneighbors(dg, l)
-            has_edge(dg, l => k) && continue # undirected
-            !has_both(dg, v, k) && continue # not undirected to v
-            isadjacent(dg, k, w) && continue # adjacent to w
+function meek_rule4(g, v, w)
+    for l in inneighbors(g, w)
+        has_edge(g, w => l) && continue # undirected
+        !isadjacent(g, v, l) && continue # not adjacent to v      
+        for k in inneighbors(g, l)
+            has_edge(g, l => k) && continue # undirected
+            !has_both(g, v, k) && continue # not undirected to v
+            isadjacent(g, k, w) && continue # adjacent to w
             return true
         end
     end

src/misc.jl

@@ -4,8 +4,7 @@ using Random
 
 Create `DiGraph` from edge-list.
 """
-function digraph(E)
-    d = maximum(flatten(E))
+function digraph(E, d = maximum(flatten(E)))
     g = DiGraph(d)
     for (i, j) in E
         add_edge!(g, i, j)
@@ -28,11 +27,12 @@ function graph(E)
 end
 
 """
-    vpairs(g)
+    arrows(g)
 
 Return the edge-list as `Pair`s.
 """
-vpairs(g) = map(Pair, collect(edges(g)))
+arrows(g::SimpleDiGraph{T}) where T = nv(g) > 0 ? map(Pair, edges(g)) : Pair{T,T}[]
+const vpairs = arrows
 
 """
     skel_oracle(g; stable=true)

src/pdag.jl

@@ -122,15 +122,71 @@ orientedge!(g, x, y) = @assert has_edge(g, x, y) && rem_edge!(g, y, x)
 All vertices in `g` connected to `x` by an undirected edge.
 Returns sorted array.
 """
-neighbors_undirected(g, x) = inneighbors(g, x) ∩ outneighbors(g, x)
+function neighbors_undirected(g, x) 
+    a = outneighbors(g, x) 
+    b = inneighbors(g, x)
+    z = sorted_intersect_(a, b)
+#    @assert z == a ∩ b
+    z
+end
 
 """
     neighbors_adjacent(g, x)
 
 All vertices in `g` connected to `x` by an any edge.
 Returns sorted array.
 """
-neighbors_adjacent(g, x) = sort(outneighbors(g, x) ∪ inneighbors(g, x))
+function neighbors_adjacent(g::SimpleDiGraph, x) 
+    a = outneighbors(g, x) 
+    b = inneighbors(g, x)
+    z = sorted_union_(a, b)
+    #@assert z == sort(a ∪ b)
+    z
+end
+function sorted_union_(x::Vector{T}, y::Vector{T}) where T
+    z = T[]
+    i = j = 1
+    while true
+        if j > length(y)
+            append!(z, @view x[i:end])
+            return z
+        end
+        if i > length(x)
+            append!(z, @view y[j:end])
+            return z
+        end        
+        if x[i] < y[j]
+            push!(z, x[i])
+            i += 1
+        elseif x[i] == y[j]
+            push!(z, x[i])
+            i += 1
+            j += 1
+        else #x[i] < y[]
+            push!(z, y[j])
+            j += 1
+        end
+    end
+end
+function sorted_intersect_(x::Vector{T}, y::Vector{T}) where T
+    z = T[]
+    i = j = 1
+    while true
+        if j > length(y) ||  i > length(x)
+            return z
+        end    
+        if x[i] < y[j]
+            i += 1
+        elseif x[i] == y[j]
+            push!(z, x[i])
+            i += 1
+            j += 1
+        else #x[i] < y[]
+            j += 1
+        end
+    end
+end
+
 
 """
     parents(g, x)

test/cpdag.jl

@@ -23,7 +23,7 @@ C3 = copy(E3)
 append!(C3, reverse.(E3))
 
 h, S = skel_oracle(digraph(E1a))
-@test vpairs(h) == E1a
+@test map(Pair, edges(h)) == E1a
 @test length(S) == 2
 @test S[Edge(2, 3)] == [1]
 @test S[Edge(1, 4)] == [2, 3]