Building an allometric food web - help with matrices


It’s like I have Lego pieces or better yet, I know what Lego pieces I need but I don’t know where to find them. I cannot find any examples to follow.

Please help me put the pieces together, or please give me an example structure that I could follow.

I need a matrix that calculates the result of

Lij = (mi/mj*Ropt * e^(1-mi/mj*Ropt))^gamma
Ropt = 100
gamma = 2

where mi and mj correspond to body masses that are of random uniform distribution, with values between 10^2 and 10^12. These are sorted from smallest to largest. Lets say for 50 species.

BM= rand(50,10^2,10^12)
BM = sort(BM)

If L[i,j] <=0.01, then L[i,j] = 0

Then, taking the results of matrix L and comparing them to a Food web matrix of random uniform distributed values between 0 and 1.

if L[i,j] >p, then Food Web matrix = 1

Then I need to know how to create a parameter matrix, and a parameter vector, as I will need to input the results into a formula.

I assume it should look something like:

parameter.matrix = function(food.web, m, I, b, c, E)
k = 8.6173324e-5 
T0 = 273.15
Temp = 293.15
exp(I)*m[i]^b*m[j]^c*exp((E*(T0 - Temp))/(k*T0*Temp)) * food.web.matrix[i,j]

where the number of rows and columns corresponds to the number of species (in this case 50)

And the vector

param.vector = function(food.web, m, I, b, E){ 
param.vector = exp(I)*m^b*exp((E*(T0 - Temp))/(k*T0*Temp))

Then using the matrix would something look like:

matrix.a = param.matrix(food.web.matrix, BM, -13.1, 0.25, -0.8, -0.38)

and vector

vector.K = param.vector(food.web.matrix, BM, 10, 0.28, 0.71) 

Any help is greatly appreciated

My matrix is actually a vector?

To get random body masses m in the cleanest way, you probably want to use the Distributions package (type ]add Distributions in the REPL):

using Distributions
m = sort(rand(Uniform(10^2, 10^12), 50))

Though do you really want this uniform distribution, or are you looking for a uniform distribution in the log of the mass? That would be

m = sort(10 .^ rand(Uniform(2,12), 50))

Now that you’ve got m one way or another, there’s various ways to express L. For example, using broadcast, operating on all elements with a single expression:

Ropt = 100
gamma = 2
# The @. macro here broadcasts the following scalar expression across `m` and `m'`
L = @. (m/m'*Ropt * exp(1-m/m'*Ropt))^gamma
# Use logical indexing to set values of L below threshold to 0:
L[L .< 0.01] .= 0

Or, if you prefer to keep all the scalar operations on elements of L in one place, you could write it as

function Lij(mi, mj)
    x = (mi/mj*Ropt * exp(1-mi/mj*Ropt))^gamma
    x < 0.01 ? 0 : x

L = Lij.(m, m')  # The extra . here broadcasts the scalar function `Lij` over the elements of `m` and `m'`

To create a random uniformly distributed (food web?) matrix F of the same size as L, with uniform values between 0 and 1, and to clamp it based on the values in L (I assume your p is F[i,j]):

F = rand(Uniform(0,1), size(L)...)
F[L .> F] .= 1    # Sets F[i,j] = 0  wherever  L[i,j] > F[i,j]

At this point, I get confused by what you want to do in the rest of your description, but I hope this gives you enough to get started with array operations.


Thank you for your reply. Based on the structure of your function, my attempt to create another matrix looks like

function ParameterMatrix(m, I, b, c, E)
    x = exp(I)*mi^b*mj^c*exp((E*(T0 - Temp))/(k*T0*Temp))

Attack = ParameterMatrix(m, -13.1, 0.25, -0.8, -0.38)

Which generates

ERROR: MethodError: no method matching ^(::Array{Float64,1}, ::Float64)
Closest candidates are:
  ^(::Missing, ::Number) at missing.jl:93
  ^(::Float64, ::Float64) at math.jl:780
  ^(::Irrational{:ℯ}, ::Number) at mathconstants.jl:91
 [1] ParameterMatrix(::Array{Float64,1}, ::Float64, ::Float64, ::Float64, ::Float64) at C:\Users\Helga\Documents\Julia\Allo Food Web.jl:22
 [2] top-level scope at none:0

I think it may have to do with mi and mj not being specified (defined) in Attack = ParameterMatrix(m, -13.1, 0.25, -0.8, -0.38) but am not sure. Errors are generated if I use m and m' or mi and mj as well.

Code for ease of reference

using Distributions
m = sort(rand(Uniform(10^2, 10^12), 50))
Ropt = 100
γ= 2
k = 8.6173324e-5 #Boltzman eV/K
T0 = 273.15 #K
Temp = 293.15 # example of 20 C degrees

function Lij(mi, mj)
    x = (mi/(mj*Ropt) * exp(1-mi/(mj*Ropt)))^γ
    x < 0.01 ? 0 : x

L = Lij.(m, m')  # The extra . here broadcasts the scalar function `Lij` over the elements of `m` and `m'`

FoodWeb = rand(Uniform(0,1), size(L)...)
FoodWeb[L .> FoodWeb] .= 1    # Sets F[i,j] = 0  wherever  L[i,j] < F[i,j]


Running your example above doesn’t work as mi and mj are not defined, should they somehow be inferred from the input m (which is otherwise unused in your function)?

The error you are getting stems from attempting to exponentiate an Array{Float64, 1} which doesn’t work. Consider:

[1., 2.]^2. # gives your error 
[1., 2.].^2. # this works, as the exponentiation is broadcasted
[1. 2.; 3. 4.]^2. # this works as well, as there is a method to exponentiate square matrices


I should have clarified that it is important for m to have differentiated i and j values (corresponding to consumer and prey, respectively). That is, m[ ,2] and m[2, ] are both the same species, but I need to be able to specify if its from the i or j side for certain formulas. Is this what @Chris_Foster meant by scalar operations?


So while it may be possible for me to write

function ParameterMatrix(m, I, b, c, E)
    x = exp(I)*m^b*m'^c*exp((E*(T0 - Temp))/(k*T0*Temp))

#attack rate a (m^2/s)
a = @. ParameterMatrix(m, -13.1, 0.25, -0.8, -0.38)

(I assume @.will take the same values of m and ‘paste’ them into this matrix)

Will I have trouble later when I need to call a [ ,i], or is it possible to use a' for this?