# 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.

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

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
end

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.

2 Likes

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))
end

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
...
Stacktrace:
 ParameterMatrix(::Array{Float64,1}, ::Float64, ::Float64, ::Float64, ::Float64) at C:\Users\Helga\Documents\Julia\Allo Food Web.jl:22
 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
end

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

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

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))
end

#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?