Building an allometric food web - help with matrices

New to Julia
1

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

2

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
3

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:
 [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
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
4

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
5

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 @c42f meant by scalar operations?

image

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?