You are here: Home Code

Code

Code produced in the lab

Allesina, S., Alonso, D., Pascual, M., 2008 - A general model for food web structure - Science 320(5876):658-661.

Here's a program that computes the likelihood of a matrix produced by the niche model. It takes as an input an adjacency matrix and the position of the species.

For example, say that you want to compute the probability of producing the following network using the niche model:

0 0 0 0 0 0 0 0 0 0
0 0 1 1 1 1 0 0 1 0
0 0 1 0 1 1 1 0 0 0
0 0 0 0 1 1 1 0 0 0
0 0 0 0 1 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 1

Using the positons:

0.00652295188046992 0.194482436170802 0.233010006370023 0.257559132762253 0.327586651546881 0.357965625124052 0.38448193622753 0.587343979859725 0.68662214349024 0.777188333682716

Here you can find the matrix and here the positions

Download the C code and compile using:

gcc NewNicheLikelihood.c -o LNICHE -lgsl -lgslcblas -Wall -DHAVE_INLINE -O3

you will need the GSL library.

Then you can run the code:

./LNICHE 10 ExampleMatrix-10-0.25.txt ExamplePositions-10-0.25.txt

This will create the output Output-ExampleMatrix-10-0.25.txt:

Compatible Adjacency Matrix
0 0 0 0 0 0 0 0 0 0
0 0 1 1 1 1 0 0 1 0
0 0 1 0 1 1 1 0 0 0
0 0 0 0 1 1 1 0 0 0
0 0 0 0 1 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 1
Species Positions
0.007 0.194 0.233 0.258 0.328 0.358 0.384 0.587 0.687 0.777
Beta: 1.500000
First    Sp 0, LogLikelihood Diet 0.0 (1.0)
Basal    Sp 1, LogLikelihood Diet -0.400166 (0.670209)
Consumer Sp 2, LogLikelihood Diet -2.290796 (0.101186)
Consumer Sp 3, LogLikelihood Diet -1.881361 (0.152383)
Consumer Sp 4, LogLikelihood Diet -2.625511 (0.072403)
Consumer Sp 5, LogLikelihood Diet -1.967941 (0.139744)
Consumer Sp 6, LogLikelihood Diet -3.605495 (0.027174)
Basal    Sp 7, LogLikelihood Diet -1.689461 (0.184619)
Consumer Sp 8, LogLikelihood Diet -3.775930 (0.022916)
Consumer Sp 9, LogLikelihood Diet -2.257073 (0.104656)
LogLikelihood: -20.493734

 

Allesina, S., Pascual, M., 2009 - Food web Models: a Plea for Groups - Ecology Letters12:652-662.

PDF

Here is an R script to visualize the results (GroupBasedRand.R).

Here is the data (AllData).

Usage in R:
Copy all the files in the same directory and start R. The code requires the "Social Network Analysis" package (sna)

>require(sna)

load the code

>source("GroupBasedRand.R")

load the data

>load("AllData")

now you can visualize the results. For example, calling:

>GroupBasedRandDigraph(skip,part_skip)

returns:

 [1] "L_ij matrix"

  1 2 3 4  5  6 7

1 0 0 0 0  3  0 0

2 0 3 4 1  1  0 0

3 0 0 0 0 13  2 0

4 9 0 0 2  0  0 9

5 0 0 0 0  1  2 0

6 0 0 0 0  0  2 0

7 0 0 0 0  1 12 3

[1] "Size_ij matrix"

 1  2  3  4  5  6  7

1  9 15 12  9 12 24 12

2 15 25 20 15 20 40 20

3 12 20 16 12 16 32 16

4  9 15 12  9 12 24 12

5 12 20 16 12 16 32 16

6 24 40 32 24 32 64 32

7 12 20 16 12 16 32 16

[1] "P_ij matrix"

1 0 0.00 0.0 0.00000000 0.2500 0.00000 0.0000

2 0 0.12 0.2 0.06666667 0.0500 0.00000 0.0000

3 0 0.00 0.0 0.00000000 0.8125 0.06250 0.0000

4 1 0.00 0.0 0.22222222 0.0000 0.00000 0.7500

5 0 0.00 0.0 0.00000000 0.0625 0.06250 0.0000

6 0 0.00 0.0 0.00000000 0.0000 0.03125 0.0000

7 0 0.00 0.0 0.00000000 0.0625 0.37500 0.1875

[1] "LogLikelihood"

[1] -113.0452

[1] "AIC"

[1] 386.0904

and produces the figure:

Skipwoth Groups

the list of networks and their partitions can be accessed through

>ls()

[1] "benguela" "bridge" "broom"

[4] "chesapeake" "coachella" "grass"

[7] "GroupBasedRandDigraph" "part_benguela" "part_bridge"

[10] "part_broom" "part_chesapeake" "part_coachella"

[13] "part_grass" "part_reef" "part_skip"

[16] "part_stmarks" "part_stmart" "plot.mygroupmatrix"

[19] "reef" "skip" "stmarks"

[22] "stmart"

Finding groups in networks: HowTo

Here you can find the code I used to find the groups in the networks, and here is a sample network that will be used in the examples below. The code has been tested on Ubuntu Linux, but works on Mac OSX and should be easy to port into Windows. If you make the code run on your Mac or Win machine, please drop me an email with a detailed description and I'll be happy to post it here.

The code is provided as is, and is distributed under GPL license (free software). Use the code at your own risk.

If you use this code for your research, please do cite the paper and, if possible, drop me an email. I'll be delighted to know that you find this approach useful.

If you have any question on the code, the algorithm etc. do not hesitate to write me, but be patient, as it could take me some days to get back to you.

Compiling the code:

You will need a gcc compiler and the GNU Scientific Library. Both are available for free.

To compile, open a Terminal and type:

gcc -lgsl -lgslcblas GA-Groups.c -o Groups -Wall -O3 -DHAVE_INLINE

This will create the executable "Groups".

Launching the code:

In the terminal, type something like:

GSL_RNG_SEED=44 ./Groups 31 chesapeake 5 500 250
  • GSL_RNG_SEED=44: choose 44 as the seed for the number generator (this allows you to repeat the exact same computation multiple times).
  • ./Groups: the name of the executable.
  • 31: number of nodes in the network.
  • chesapeake: file containing the adjacency matrix for the network. For the format, see the example file.
  • 5: desired number of groups.
  • 500: population size for the Genetic Algorithm.
  • 250: number of generations for the Genetic Algorithm.

Output:

While running the code, you should see in the terminal something like:

N=31

File=chesapeake

Groups=5

Pop=500

Gen=250

GSL_RNG_SEED=34

-Gen0MinAIC 393.241284 MinAIC-Gen 393.241284 Groups 5
----------Gen 10MinAIC 393.241284 MinAIC-Gen 393.241284 Groups 5
----------Gen 20MinAIC 393.241284 MinAIC-Gen 393.241284 Groups 5
----------Gen 30MinAIC 393.241284 MinAIC-Gen 393.241284 Groups 5
----------Gen 40MinAIC 393.241284 MinAIC-Gen 393.241284 Groups 5

....

Once the computation has finished, you should also have a file named "Out-[filename]-[groups]-[AIC].txt" (e.g. Out-chesapeake-5-393.241284.txt) that contains the vector specifying the way of grouping the nodes.

Finding the best configuration:

In general, start with only 1 group, and then progress until the AIC starts increasing: in this way you can find the minimum AIC.

The larger (or more complex) the network, the slower the computation.

Increasing the population size and the number of generations will probably provide better results. If this is the case, increase the population size until no better solution is found. Generally speaking, using more groups will reuquire larger populations, as the number of possible configurations increases dramatically.

Producing figures as the ones in the paper:

You will need to have R and GraphViz installed. Both are freely available.

Copy this script in the same folder where you have the data and the output from the Groups code. Open R and type:

>source("DrawGroupsGraphViz.R")
>DrawGroups("chesapeake","Out-chesapeake-5-393.241284.txt")

Where the first argument is the file containing the adjacency matrix, and the second is the file containing the solution found by the Groups program. You can pass a third argument: a vector contaning the names of the nodes (ordered as in the adjacency matrix).

Running the script specified above creates a ".dot" file that can now be interpreted by GraphViz.

Exit R and type in the terminal:

fdp -Tpng -O Out-chesapeake-5-393.241284.txt.dot

Where

fdp is the layout program;
-Tpng is the desired output format (e.g. -Tpdf, -Tsvg, ...);
-O means that you want to conserve the file name (otherwise use -o myoutputfilename);
Out-chesapeake-5-393.241284.txt.dot is the file created by R.

You now have your figure. Enjoy!

Document Actions