In this example, we will work on another plant and pollination network. As in the previous example, we will download it and turn it into an adjacency matrix.

Because there will be some randomizations involved, we will start Julia with julia -p n, where n is the number of CPUs to run on. If you have a 4-core machine, then julia -p 3 is fine. This will not give excellent scaling because there is some overhead when using a small number of cores, but it will still cut the computing time.

Brim has to be loaded on all workers, so it should be imported with

@everywhere using Brim

Have a look at the Parallel Computing page of Julia's documentation if you need more information.

download(
  "https://www.nceas.ucsb.edu/interactionweb/data/plant_pollinator/text/ino_matr_f.txt",
  "ino.txt"
)
@everywhere A = map((x) -> x>0?1:0, int(readdlm("ino.txt")))

The first step is to optimize its modularity. Since this is a small network, we will use a random starting partition:

@everywhere rndbrim = (A) -> A |> partition_random |> recursive_brim!
ino_partitions = pmap(rndbrim, [A for i in 1:1000])

For reference, on my machine (runing three workers), the pmap (parallel) version runs 1000 replicates within 2 seconds, and the map (serial) does the same thing in 4 seconds. The gain in time becomes greater when using more replicates (or more cores...).

We can then filter only the partitions giving the best modularity:

ino_q = map(Q, ino_partitions)
best_partition_id = filter((x) -> ino_q[x] == maximum(ino_q), 1:length(ino_q))[1]
ino_mod = ino_partitions[best_partition_id]

The modularity of this partition (Q(ino_mod)) should be around 0.39 -- which is low, and may not be significant.

We will test this using permutations of the original matrix.

To generate a single randomization keeping the rows and columns marginals, we can use:

A_rnd = null_preserve_marginals(A)

This can take a few seconds to run (around 8 sec. on my machine). We will need to do this in parallel, which is not difficult:

random_matrices = pmap(null_preserve_marginals, [A for i in 1:50])

Note that I used 50 replicates for the sake of not spending minutes waiting for the permutations to finish. It is of course unreasonably low.

Now, we can apply the same steps as above, to measure the distribution of modularities in the shuffled matrices. Because we have not assigned random_matrices everywhere, we will not use pmap:

expected_q = map(Q, map(rndbrim, random_matrices))

This array will have the expected distribution of Q, for this network, under the null hypothesis that interactions are randomly distributed (this is a way to test the impact of the degree distribution on the modularity). There are different ways to get estimates of the p-value, but an easy one is to measure the proportion of randomized matrices that are less modular than the original one:

approx_pval = maximum([sum(expected_q .>= Q(ino_mod))/length(expected_q) 1/length(expected_q)])

Note that the p-value cannot be lower than 1/m, where m is the number of shuffled matrices tested. In this example, this should return a value around 0.02. This indicates that the modular structure of this network, although not strong, is a significant deviation from the random expectation.