mschauer / CausalInference.jl

Causal inference, graphical models and structure learning in Julia
https://mschauer.github.io/CausalInference.jl/latest/
Other
190 stars 24 forks source link

Stable threaded skeleton search #101

Closed mschauer closed 1 year ago

mschauer commented 1 year ago

Timings versus pc-alg on 4 threads:

> system.time(skel.fit <- skeleton(suffStat = list(C = cor(dat), n = n),
+                                  indepTest = gaussCItest, ## (partial correlations)
+                                  alpha = 0.001, p = p,  method = "stable.fast"))#
   user  system elapsed 
 13.693   0.033  13.759 
> skel.fit
Object of class 'pcAlgo', from Call:
skeleton(suffStat = list(C = cor(dat), n = n), indepTest = gaussCItest, 
    alpha = 0.001, p = p, method = "stable.fast")
Number of undirected edges:  254 
Number of directed edges:    0 
Total number of edges:       254 
> 

versus

julia> h, s = @time CausalInference.skeleton(d, gausscitest, (cor(X), n), -quantile(Normal(), alpha/2)); h
 13.371163 seconds (42.97 M allocations: 10.909 GiB, 19.89% gc time)
{64, 254} undirected simple Int64 graph
mwien commented 1 year ago

Nice! Before we were slower?

mschauer commented 1 year ago

Yes, a bit slower as the C implementation "stable.fast" contained in pcalg, by a factor of 2 which is okay (and much faster than the native implementation "stable" of course which is the default in pcalg )

mwien commented 1 year ago

What if you pass numCores = 4 in the pcalg call?

E.g.: system.time(skel.fit <- skeleton(suffStat = list(C = cor(dat), n = n), indepTest = gaussCItest, alpha = 0.001, p = p, method = "stable.fast", numCores = 4))#

mschauer commented 1 year ago

Hm, you are right, but not much happens:

> system.time(skel.fit <- skeleton(suffStat = list(C = cor(dat), n = n),
+                                  indepTest = gaussCItest, ## (partial correlations)
+                                  alpha = 0.001, p = p,  method = "stable.fast", numCores = 4))#, labels = V
   user  system elapsed 
 13.838   0.063  14.031 
> 
mwien commented 1 year ago

Hmm weird... But not our problem I guess...

mschauer commented 1 year ago

That is more tricky on a second note. Not sure if I even have OpenMP support. See https://mac.r-project.org/openmp/#do

Agreed that this is above our pay grade.